c setup_conf.f
         

c-------------------------------------------------------------------
c setup initial conformation
c initial-1 flat circle
c         2 from init_conf 'no euler'
c         3 from Init_conf with euler
c setup information in x(i,j) as well as r(i,j)and (f,v,u)
c for linear or nicked DNA, no (f,v)
      	subroutine setup_conf(n,initial)
      	implicit double precision(a-h, o-z)
       	parameter(maxa=500,maxn=maxa*3,pi2=2.d0*3.1415926535d0)
       	integer initial,n
	common/matrix/x(maxa,3),bangle(maxa)
        common/matrix2/ r(maxa,maxa)
        common/ddx/d_x(maxa,3)
      	common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3)
        common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma,
     +amass,am(maxn),coef,dtime2
	common/tor_param/dr,cg,phi0,dlk,phii,subseg
          common/ischeme/icircular,iele,ihydro,iclose,iknot
     +,neuler
      common/xoldtw/olddlk,xold(maxa,3),fold(maxa,3),vold(maxa,3)
     +,uold(maxa,3)

	do i=1,n
	 r(i,i)=0
	enddo
        write(*,*) 'icircular=',icircular,'initial=',initial
c flat circle or straight line
      if (initial.eq.1) then
	if (icircular.eq.1) then
          call init_circular(n)
        else
	  call init_linear(n)
      endif  
      endif

c old configuration but assign new euler angles
      if (initial.eq.2) then
       open(unit=139, file='Init_conf',status='old')
       do k=1,n
        read(139,*) x(k,1),x(k,2),x(k,3)
       enddo
       call update_bond(n)
       call update_rij(n)

c d_x and r(i,i+1) are updated at update_bond(n) 
c calculate the initial Wr and deduct it from the initial phi0
  
      if ((icircular.eq.1).and.(iclose.eq.1)) then
      do i = 1, 3
         x(n+1,i) = x(1,i)
      enddo
      do i = 1, n
         do j=1,3
            d_x(i,j)=x(i+1,j)-x(i,j)
         enddo
      enddo  
       call kwr(n,topl,jwrold,ierr,0)
       call bwr(n,2,n,beeta)
      writhe = jwrold + beeta     

      dlk=dlk-writhe
      call init_phi0(n)
      call update_bond(n)
      
       do i = 1, n
           angle = (i-1)*phii
           cos1 = cos(angle)
           sin1 = sin(angle)
           f(i,1) = -sin1*u(i,2)
           f(i,2) = sin1*u(i,1)
           f(i,3) = -cos1
           v(i,1) = -cos1*u(i,2)
           v(i,2) = cos1*u(i,1)
           v(i,3) = sin1
        enddo
        call update_rij(n)
        call update_euler(n)  
c      adjust the dlk if it is not as same as the initial olddlk
        call update_euler(n)
        call cal_twist(n,twist)
        write(*,*) 'before adjusted twist=',twist,dlk
      endif
      endif

c old 'Init_conf'
      if (initial.eq.3) then
       open(unit=139, file='Init_conf',status='old')
       read(139,*) ((x(mi,mii),mii=1,3),mi=1,n)
       if ((icircular.eq.1).and.(iclose.eq.1)) then
         read(139,*) ((f(mi,mii),mii=1,3),mi=1,n)
         read(139,*) ((v(mi,mii),mii=1,3),mi=1,n)
         read(139,*) ((u(mi,mii),mii=1,3),mi=1,n)
        call update_euler(n)
           endif
          call update_rij(n)
          call update_bond(n)  
      endif
      return
      end