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

   1 c-------------------------------------------------------------------------------
   2       subroutine velverlet_step(itime)
   3 c-------------------------------------------------------------------------------
   4 c  Perform a single velocity Verlet step; the time step can be rescaled if
   5 c  increments in accelerations exceed the threshold
   6 c-------------------------------------------------------------------------------
   7       implicit real*8 (a-h,o-z)
   8       include 'DIMENSIONS'
   9       include 'COMMON.CONTROL'
  10       include 'COMMON.VAR'
  11       include 'COMMON.MD'
  12       include 'COMMON.LANGEVIN'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.DERIV'
  15       include 'COMMON.GEO'
  16       include 'COMMON.LOCAL'
  17       include 'COMMON.INTERACT'
  18       include 'COMMON.IOUNITS'
  19       include 'COMMON.NAMES'
  20       include 'COMMON.TIME1'
  21       include 'COMMON.MUCA'
  22       double precision energia(0:n_ene),vcm(3),incr(3)
  23       double precision cm(3),L(3)
  24       integer ilen,count,rstcount
  25       external ilen
  26       character*50 tytul
  27       integer maxcount_scale /20/
  28       common /gucio/ cm, energia
  29       double precision stochforcvec(MAXRES6)
  30       common /stochcalc/ stochforcvec
  31       integer itime
  32       logical scale
  33 c
  34       scale=.true.
  35       icount_scale=0
  36       if (lang.eq.1) then
  37         call sddir_precalc
  38       else if (lang.eq.2 .or. lang.eq.3) then
  39         call stochastic_force(stochforcvec)
  40       endif
  41       itime_scal=0
  42       do while (scale)
  43         icount_scale=icount_scale+1
  44         if (icount_scale.gt.maxcount_scale) then
  45           write (iout,*)
  46      &      "ERROR: too many attempts at scaling down the time step. ",
  47      &      "amax=",amax,"epdrift=",epdrift,
  48      &      "damax=",damax,"edriftmax=",edriftmax,
  49      &      "d_time=",d_time
  50             stop
  51         endif
  52 c First step of the velocity Verlet algorithm
  53         if (lang.eq.2) then
  54           call sd_verlet1
  55         else if (lang.eq.3) then
  56           call sd_verlet1_ciccotti
  57         else if (lang.eq.1) then
  58           call sddir_verlet1
  59         else
  60           call verlet1
  61         endif
  62 c Build the chain from the newly calculated coordinates
  63         call chainbuild_cart
  64         if (rattle) call rattle1
  65         if (large.and. mod(itime,ntwe).eq.0) then
  66           write (iout,*) "Cartesian and internal coordinates: step 1"
  67           call cartprint
  68           call intout
  69           write (iout,*) "Accelerations"
  70           do i=0,nres
  71             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
  72      &      (d_a(j,i+nres),j=1,3)
  73           enddo
  74           write (iout,*) "Velocities, step 1"
  75           do i=0,nres
  76             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
  77      &      (d_t(j,i+nres),j=1,3)
  78           enddo
  79         endif
  80         tt0 = tcpu()
  81 c Calculate energy and forces
  82         call zerograd
  83         call etotal(energia)
  84         potE=energia(0)-energia(20)
  85         call cartgrad
  86 c Get the new accelerations
  87         call lagrangian
  88         t_enegrad=t_enegrad+tcpu()-tt0
  89 c Determine maximum acceleration and scale down the timestep if needed
  90         call max_accel
  91         call predict_edrift(epdrift)
  92         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
  93 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
  94           scale=.true.
  95           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
  96      &      /dlog(2.0d0)+1
  97           itime_scal=itime_scal+ifac_time
  98 c          fac_time=dmin1(damax/amax,0.5d0)
  99           fac_time=0.5d0**ifac_time
 100           d_time=d_time*fac_time
 101           if (lang.eq.2 .or. lang.eq.3) then
 102 c            write (iout,*) "Calling sd_verlet_setup: 1"
 103 c Rescale the stochastic forces and recalculate or restore
 104 c the matrices of tinker integrator
 105             if (itime_scal.gt.maxflag_stoch) then
 106               if (large) write (iout,'(a,i5,a)')
 107      &         "Calculate matrices for stochastic step;",
 108      &         " itime_scal ",itime_scal
 109               if (lang.eq.2) then
 110                 call sd_verlet_p_setup
 111               else
 112                 call sd_verlet_ciccotti_setup
 113               endif
 114               write (iout,'(2a,i3,a,i3,1h.)')
 115      &         "Warning: cannot store matrices for stochastic",
 116      &         " integration because the index",itime_scal,
 117      &         " is greater than",maxflag_stoch
 118               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
 119      &         " integration Langevin algorithm for better efficiency."
 120             else if (flag_stoch(itime_scal)) then
 121               if (large) write (iout,'(a,i5,a,l1)')
 122      &         "Restore matrices for stochastic step; itime_scal ",
 123      &         itime_scal," flag ",flag_stoch(itime_scal)
 124               do i=1,dimen
 125                 do j=1,dimen
 126                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
 127                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
 128                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
 129                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
 130                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
 131                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
 132                 enddo
 133               enddo
 134             else
 135               if (large) write (iout,'(2a,i5,a,l1)')
 136      &         "Calculate & store matrices for stochastic step;",
 137      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
 138               if (lang.eq.2) then
 139                 call sd_verlet_p_setup
 140               else
 141                 call sd_verlet_ciccotti_setup
 142               endif
 143               flag_stoch(ifac_time)=.true.
 144               do i=1,dimen
 145                 do j=1,dimen
 146                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
 147                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
 148                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
 149                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
 150                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
 151                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
 152                 enddo
 153               enddo
 154             endif
 155             fac_time=1.0d0/dsqrt(fac_time)
 156             do i=1,dimen
 157               stochforcvec(i)=fac_time*stochforcvec(i)
 158             enddo
 159           else if (lang.eq.1) then
 160 c Rescale the accelerations due to stochastic forces
 161             fac_time=1.0d0/dsqrt(fac_time)
 162             do i=1,dimen
 163               d_as_work(i)=d_as_work(i)*fac_time
 164             enddo
 165           endif
 166           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
 167      &      "itime",itime," Timestep scaled down to ",
 168      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
 169         else
 170 c Second step of the velocity Verlet algorithm
 171           if (lang.eq.2) then
 172             call sd_verlet2
 173           else if (lang.eq.3) then
 174             call sd_verlet2_ciccotti
 175           else if (lang.eq.1) then
 176             call sddir_verlet2
 177           else
 178             call verlet2
 179           endif
 180           if (rattle) call rattle2
 181           totT=totT+d_time
 182           if (d_time.ne.d_time0) then
 183             d_time=d_time0
 184             if (lang.eq.2 .or. lang.eq.3) then
 185               if (large) write (iout,'(a)')
 186      &         "Restore original matrices for stochastic step"
 187 c              write (iout,*) "Calling sd_verlet_setup: 2"
 188 c Restore the matrices of tinker integrator if the time step has been restored
 189               do i=1,dimen
 190                 do j=1,dimen
 191                   pfric_mat(i,j)=pfric0_mat(i,j,0)
 192                   afric_mat(i,j)=afric0_mat(i,j,0)
 193                   vfric_mat(i,j)=vfric0_mat(i,j,0)
 194                   prand_mat(i,j)=prand0_mat(i,j,0)
 195                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
 196                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
 197                 enddo
 198               enddo
 199             endif
 200           endif
 201           scale=.false.
 202         endif
 203       enddo
 204 c Calculate the kinetic and the total energy and the kinetic temperature
 205       call kinetic(EK)
 206       totE=EK+potE
 207 c diagnostics
 208 c      call kinetic1(EK1)
 209 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
 210 c end diagnostics
 211 c Couple the system to Berendsen bath if needed
 212       if (tbf .and. lang.eq.0) then
 213         call verlet_bath
 214       endif
 215       kinetic_T=2.0d0/(dimen*Rb)*EK
 216 c Backup the coordinates, velocities, and accelerations
 217       do i=0,2*nres
 218         do j=1,3
 219           dc_old(j,i)=dc(j,i)
 220           d_t_old(j,i)=d_t(j,i)
 221           d_a_old(j,i)=d_a(j,i)
 222         enddo
 223       enddo
 224       if (mod(itime,ntwe).eq.0 .and. large) then
 225         write (iout,*) "Velocities, step 2"
 226         do i=0,nres
 227           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 228      &      (d_t(j,i+nres),j=1,3)
 229         enddo
 230       endif
 231       return
 232       end