source/unres/src_MD/md-diff/md/init_md.f

   1 c---------------------------------------------------------
   2       subroutine init_MD
   3 c  Set up the initial conditions of a MD simulation
   4       implicit real*8 (a-h,o-z)
   5       include 'DIMENSIONS'
   6 #ifdef MP
   7       include 'mpif.h'
   8       include 'COMMON.INFO'
   9       include 'COMMON.SETUP'
  10       character*4 liczba
  11       character*16 form
  12 #endif
  13       include 'COMMON.CONTROL'
  14       include 'COMMON.VAR'
  15       include 'COMMON.MD'
  16       include 'COMMON.LANGEVIN'
  17       include 'COMMON.CHAIN'
  18       include 'COMMON.DERIV'
  19       include 'COMMON.GEO'
  20       include 'COMMON.LOCAL'
  21       include 'COMMON.INTERACT'
  22       include 'COMMON.IOUNITS'
  23       include 'COMMON.NAMES'
  24       real*8 energia(0:n_ene),energia_long(0:n_ene),
  25      &  energia_short(0:n_ene),vcm(3),incr(3)
  26       double precision cm(3),L(3),xv,sigv,lowb,highb
  27       character*256 qstr
  28       integer ilen
  29       external ilen
  30       character*50 tytul
  31       common /gucio/ cm
  32       d_time0=d_time
  33       write(iout,*) "d_time", d_time
  34 c Compute the standard deviations of stochastic forces for Langevin dynamics
  35 c if the friction coefficients do not depend on surface area
  36       if (lang.gt.0 .and. .not.surfarea) then
  37         do i=nnt,nct-1
  38           stdforcp(i)=stdfp*dsqrt(gamp)
  39         enddo
  40         do i=nnt,nct
  41           stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
  42         enddo
  43       endif
  44 c Open the pdb file for snapshotshots
  45 #ifdef MPI
  46       if (nprocs.eq.1) then
  47         npos=3
  48       else
  49         npos = dlog10(dfloat(nprocs-1))+1
  50       endif
  51       if (npos.lt.3) npos=3
  52       write (liczba,'(i1)') npos
  53       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
  54      &  //')'
  55       write (liczba,form) myrank
  56       if(mdpdb) then
  57         open(ipdb,
  58      &  file=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
  59      &  //".pdb")
  60       else
  61         cartname=prefix(:ilen(prefix))//"_MD"//liczba(:ilen(liczba))
  62      &  //".cx"
  63       endif
  64 #else
  65       if(mdpdb) then
  66          open(ipdb,file=prefix(:ilen(prefix))//"_MD.pdb")
  67       else
  68          cartname=prefix(:ilen(prefix))//"_MD.cx"
  69       endif
  70 #endif
  71       if (usampl) then
  72         write (qstr,'(256(1h ))')
  73         ipos=1
  74         do i=1,nfrag
  75           iq = qinfrag(i)*10
  76           iw = wfrag(i)/100
  77           if (iw.gt.0) then
  78             write (iout,*) "Frag",qinfrag(i),wfrag(i),iq,iw
  79             write (qstr(ipos:ipos+6),'(2h_f,i1,1h_,i1,1h_,i1)') i,iq,iw
  80             ipos=ipos+7
  81           endif
  82         enddo
  83         do i=1,npair
  84           iq = qinpair(i)*10
  85           iw = wpair(i)/100
  86           if (iw.gt.0) then
  87             write (iout,*) "Pair",i,qinpair(i),wpair(i),iq,iw
  88             write (qstr(ipos:ipos+6),'(2h_p,i1,1h_,i1,1h_,i1)') i,iq,iw
  89             ipos=ipos+7
  90           endif
  91         enddo
  92         pdbname=pdbname(:ilen(pdbname)-4)//qstr(:ipos-1)//'.pdb'
  93         cartname=cartname(:ilen(cartname)-2)//qstr(:ipos-1)//'.x'
  94         statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
  95       endif
  96       icg=1
  97       if (rest) then
  98         write(iout,*) "Initial state will be read from file ",
  99      &    rest2name(:ilen(rest2name))
 100         call readrst
 101       else
 102 c Generate initial velocities
 103         write(iout,*) "Initial velocities randomly generated"
 104         call random_vel
 105         totT=0.0d0
 106       endif
 107 c      rest2name = prefix(:ilen(prefix))//'.rst'
 108       write(iout,*) "Initial backbone velocities"
 109       do i=nnt,nct-1
 110         write(iout,*) (d_t(j,i),j=1,3)
 111       enddo
 112       write(iout,*) "Initial side-chain velocities"
 113       do i=nnt,nct
 114         write(iout,*) (d_t(j,i+nres),j=1,3)
 115       enddo
 116 c  Zeroing the total angular momentum of the system
 117       write(iout,*) "Calling the zero-angular
 118      & momentum subroutine"
 119       call inertia_tensor
 120 c  Getting the potential energy and forces and velocities and accelerations
 121       call vcm_vel(vcm)
 122 c      write (iout,*) "velocity of the center of the mass:"
 123 c      write (iout,*) (vcm(j),j=1,3)
 124       do j=1,3
 125         d_t(j,0)=d_t(j,0)-vcm(j)
 126       enddo
 127 c Removing the velocity of the center of mass
 128       call vcm_vel(vcm)
 129       write (iout,*) "vcm right after adjustment:"
 130       write (iout,*) (vcm(j),j=1,3)
 131       if (.not.rest) then
 132          call chainbuild
 133       endif
 134       call chainbuild_cart
 135       call kinetic(EK)
 136       if (tbf) then
 137         call verlet_bath(EK)
 138       endif
 139       kinetic_T=2.0d0/(dimen*Rb)*EK
 140       call cartprint
 141       call intout
 142       call zerograd
 143       call etotal(energia)
 144       potE=energia(0)
 145       call cartgrad
 146       call lagrangian
 147       call max_accel
 148       if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
 149       write(iout,*) "Potential energy"
 150       write(iout,*) (energia(i),i=0,n_ene)
 151       potE=energia(0)-energia(20)
 152       totE=EK+potE
 153       call statout(itime)
 154       write (iout,*) "Initial:",
 155      &   " Kinetic energy",EK," potential energy",potE,
 156      &   " total energy",totE," maximum acceleration ",
 157      &   amax
 158       if (large) then
 159         write (iout,*) "Initial velocities"
 160         do i=0,nres
 161           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 162      &    (d_t(j,i+nres),j=1,3)
 163         enddo
 164         write (iout,*) "Initial accelerations"
 165         do i=0,nres
 166           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 167      &    (d_a(j,i+nres),j=1,3)
 168         enddo
 169       endif
 170       do i=0,2*nres
 171         do j=1,3
 172           dc_old(j,i)=dc(j,i)
 173           d_t_old(j,i)=d_t(j,i)
 174           d_a_old(j,i)=d_a(j,i)
 175         enddo
 176 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
 177       enddo
 178       if (RESPA) then
 179         call zerograd
 180         call etotal_long(energia_long)
 181         call cartgrad
 182         call lagrangian
 183 c        call etotal_short(energia_short)
 184 c        write (iout,*) "energia_long",energia_long(0),
 185 c     &   " energia_short",energia_short(0),
 186 c     &   " total",energia_long(0)+energia_short(0)
 187       endif
 188       return
 189       end