source/unres/src_MD/src/md-diff/np/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),E_short
  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
 146       if(tnp .or. tnp1) then
 147        s_np=1.0
 148        pi_np=0.0
 149        HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
 150        H0=Hnose1
 151        write(iout,*) 'H0= ',H0
 152       endif
 153
 154        if(tnh) then
 155         HNose1=Hnose_nh(EK,potE)
 156         H0=HNose1
 157         write (iout,*) 'H0= ',H0
 158        endif
 159
 160
 161       call cartgrad
 162       call lagrangian
 163       call max_accel
 164       if (amax*d_time .gt. dvmax) d_time=d_time*dvmax/amax
 165       write(iout,*) "Potential energy"
 166       write(iout,*) (energia(i),i=0,n_ene)
 167       potE=energia(0)-energia(20)
 168       totE=EK+potE
 169       itime=0
 170       call statout(itime)
 171       write (iout,*) "Initial:",
 172      &   " Kinetic energy",EK," potential energy",potE,
 173      &   " total energy",totE," maximum acceleration ",
 174      &   amax
 175       if (large) then
 176         write (iout,*) "Initial velocities"
 177         do i=0,nres
 178           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 179      &    (d_t(j,i+nres),j=1,3)
 180         enddo
 181         write (iout,*) "Initial accelerations"
 182         do i=0,nres
 183           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 184      &    (d_a(j,i+nres),j=1,3)
 185         enddo
 186       endif
 187       do i=0,2*nres
 188         do j=1,3
 189           dc_old(j,i)=dc(j,i)
 190           d_t_old(j,i)=d_t(j,i)
 191           d_a_old(j,i)=d_a(j,i)
 192         enddo
 193 c        write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3)
 194       enddo
 195       if (RESPA) then
 196         call zerograd
 197         call etotal_long(energia_long)
 198         E_long=energia_long(0)
 199         if(tnp .or. tnp1) then
 200          call etotal_short(energia_short)
 201          E_short=energia_short(0)
 202          HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen)
 203          Csplit=Hnose1
 204 c         Csplit =110
 205 c_new_var_csplit          Csplit=H0-E_long
 206 c          Csplit = H0-energia_short(0)
 207           write(iout,*) 'Csplit= ',Csplit
 208         endif
 209         call cartgrad
 210         call lagrangian
 211 c        call etotal_short(energia_short)
 212 c        write (iout,*) "energia_long",energia_long(0),
 213 c     &   " energia_short",energia_short(0),
 214 c     &   " total",energia_long(0)+energia_short(0)
 215       endif
 216       return
 217       end