source/unres/src_MD/src/md-diff/mts/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 #ifdef MPI
  10       include 'mpif.h'
  11       integer ierror,ierrcode
  12 #endif
  13       include 'COMMON.SETUP'
  14       include 'COMMON.CONTROL'
  15       include 'COMMON.VAR'
  16       include 'COMMON.MD'
  17 #ifndef LANG0
  18       include 'COMMON.LANGEVIN'
  19 #else
  20       include 'COMMON.LANGEVIN.lang0'
  21 #endif
  22       include 'COMMON.CHAIN'
  23       include 'COMMON.DERIV'
  24       include 'COMMON.GEO'
  25       include 'COMMON.LOCAL'
  26       include 'COMMON.INTERACT'
  27       include 'COMMON.IOUNITS'
  28       include 'COMMON.NAMES'
  29       include 'COMMON.TIME1'
  30       include 'COMMON.MUCA'
  31       double precision vcm(3),incr(3)
  32       double precision cm(3),L(3)
  33       integer ilen,count,rstcount
  34       external ilen
  35       character*50 tytul
  36       integer maxcount_scale /20/
  37       common /gucio/ cm
  38       double precision stochforcvec(MAXRES6)
  39       common /stochcalc/ stochforcvec
  40       integer itime
  41       logical scale
  42 c
  43       scale=.true.
  44       icount_scale=0
  45       if (lang.eq.1) then
  46         call sddir_precalc
  47       else if (lang.eq.2 .or. lang.eq.3) then
  48 #ifndef LANG0
  49         call stochastic_force(stochforcvec)
  50 #else
  51         write (iout,*)
  52      &   "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
  53 #ifdef MPI
  54         call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
  55 #endif
  56         stop
  57 #endif
  58       endif
  59       itime_scal=0
  60       do while (scale)
  61         icount_scale=icount_scale+1
  62         if (icount_scale.gt.maxcount_scale) then
  63           write (iout,*)
  64      &      "ERROR: too many attempts at scaling down the time step. ",
  65      &      "amax=",amax,"epdrift=",epdrift,
  66      &      "damax=",damax,"edriftmax=",edriftmax,
  67      &      "d_time=",d_time
  68           call flush(iout)
  69 #ifdef MPI
  70           call MPI_Abort(MPI_COMM_WORLD,IERROR,IERRCODE)
  71 #endif
  72           stop
  73         endif
  74 c First step of the velocity Verlet algorithm
  75         if (lang.eq.2) then
  76 #ifndef LANG0
  77           call sd_verlet1
  78 #endif
  79         else if (lang.eq.3) then
  80 #ifndef LANG0
  81           call sd_verlet1_ciccotti
  82 #endif
  83         else if (lang.eq.1) then
  84           call sddir_verlet1
  85         else
  86           call verlet1
  87         endif
  88 c Build the chain from the newly calculated coordinates
  89         call chainbuild_cart
  90         if (rattle) call rattle1
  91         if (ntwe.ne.0) then
  92         if (large.and. mod(itime,ntwe).eq.0) then
  93           write (iout,*) "Cartesian and internal coordinates: step 1"
  94           call cartprint
  95           call intout
  96           write (iout,*) "dC"
  97           do i=0,nres
  98             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(dc(j,i),j=1,3),
  99      &      (dc(j,i+nres),j=1,3)
 100           enddo
 101           write (iout,*) "Accelerations"
 102           do i=0,nres
 103             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 104      &      (d_a(j,i+nres),j=1,3)
 105           enddo
 106           write (iout,*) "Velocities, step 1"
 107           do i=0,nres
 108             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 109      &      (d_t(j,i+nres),j=1,3)
 110           enddo
 111         endif
 112         endif
 113 #ifdef MPI
 114         tt0 = MPI_Wtime()
 115 #else
 116         tt0 = tcpu()
 117 #endif
 118 c Calculate energy and forces
 119         call zerograd
 120         call etotal(potEcomp)
 121 #ifdef TIMING_ENE
 122 #ifdef MPI
 123         t_etotal=t_etotal+MPI_Wtime()-tt0
 124 #else
 125         t_etotal=t_etotal+tcpu()-tt0
 126 #endif
 127 #endif
 128         potE=potEcomp(0)-potEcomp(20)
 129         call cartgrad
 130 c Get the new accelerations
 131         call lagrangian
 132 #ifdef MPI
 133         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 134 #else
 135         t_enegrad=t_enegrad+tcpu()-tt0
 136 #endif
 137 c Determine maximum acceleration and scale down the timestep if needed
 138         call max_accel
 139         amax=amax/(itime_scal+1)**2
 140         call predict_edrift(epdrift)
 141         if (amax/(itime_scal+1).gt.damax .or. epdrift.gt.edriftmax) then
 142 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
 143           scale=.true.
 144           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
 145      &      /dlog(2.0d0)+1
 146           itime_scal=itime_scal+ifac_time
 147 c          fac_time=dmin1(damax/amax,0.5d0)
 148           fac_time=0.5d0**ifac_time
 149           d_time=d_time*fac_time
 150           if (lang.eq.2 .or. lang.eq.3) then
 151 #ifndef LANG0
 152 c            write (iout,*) "Calling sd_verlet_setup: 1"
 153 c Rescale the stochastic forces and recalculate or restore
 154 c the matrices of tinker integrator
 155             if (itime_scal.gt.maxflag_stoch) then
 156               if (large) write (iout,'(a,i5,a)')
 157      &         "Calculate matrices for stochastic step;",
 158      &         " itime_scal ",itime_scal
 159               if (lang.eq.2) then
 160                 call sd_verlet_p_setup
 161               else
 162                 call sd_verlet_ciccotti_setup
 163               endif
 164               write (iout,'(2a,i3,a,i3,1h.)')
 165      &         "Warning: cannot store matrices for stochastic",
 166      &         " integration because the index",itime_scal,
 167      &         " is greater than",maxflag_stoch
 168               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
 169      &         " integration Langevin algorithm for better efficiency."
 170             else if (flag_stoch(itime_scal)) then
 171               if (large) write (iout,'(a,i5,a,l1)')
 172      &         "Restore matrices for stochastic step; itime_scal ",
 173      &         itime_scal," flag ",flag_stoch(itime_scal)
 174               do i=1,dimen
 175                 do j=1,dimen
 176                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
 177                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
 178                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
 179                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
 180                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
 181                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
 182                 enddo
 183               enddo
 184             else
 185               if (large) write (iout,'(2a,i5,a,l1)')
 186      &         "Calculate & store matrices for stochastic step;",
 187      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
 188               if (lang.eq.2) then
 189                 call sd_verlet_p_setup
 190               else
 191                 call sd_verlet_ciccotti_setup
 192               endif
 193               flag_stoch(ifac_time)=.true.
 194               do i=1,dimen
 195                 do j=1,dimen
 196                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
 197                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
 198                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
 199                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
 200                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
 201                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
 202                 enddo
 203               enddo
 204             endif
 205             fac_time=1.0d0/dsqrt(fac_time)
 206             do i=1,dimen
 207               stochforcvec(i)=fac_time*stochforcvec(i)
 208             enddo
 209 #endif
 210           else if (lang.eq.1) then
 211 c Rescale the accelerations due to stochastic forces
 212             fac_time=1.0d0/dsqrt(fac_time)
 213             do i=1,dimen
 214               d_as_work(i)=d_as_work(i)*fac_time
 215             enddo
 216           endif
 217           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
 218      &      "itime",itime," Timestep scaled down to ",
 219      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
 220         else
 221 c Second step of the velocity Verlet algorithm
 222           if (lang.eq.2) then
 223 #ifndef LANG0
 224             call sd_verlet2
 225 #endif
 226           else if (lang.eq.3) then
 227 #ifndef LANG0
 228             call sd_verlet2_ciccotti
 229 #endif
 230           else if (lang.eq.1) then
 231             call sddir_verlet2
 232           else
 233             call verlet2
 234           endif
 235           if (rattle) call rattle2
 236           totT=totT+d_time
 237           if (d_time.ne.d_time0) then
 238             d_time=d_time0
 239 #ifndef   LANG0
 240             if (lang.eq.2 .or. lang.eq.3) then
 241               if (large) write (iout,'(a)')
 242      &         "Restore original matrices for stochastic step"
 243 c              write (iout,*) "Calling sd_verlet_setup: 2"
 244 c Restore the matrices of tinker integrator if the time step has been restored
 245               do i=1,dimen
 246                 do j=1,dimen
 247                   pfric_mat(i,j)=pfric0_mat(i,j,0)
 248                   afric_mat(i,j)=afric0_mat(i,j,0)
 249                   vfric_mat(i,j)=vfric0_mat(i,j,0)
 250                   prand_mat(i,j)=prand0_mat(i,j,0)
 251                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
 252                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
 253                 enddo
 254               enddo
 255             endif
 256 #endif
 257           endif
 258           scale=.false.
 259         endif
 260       enddo
 261 c Calculate the kinetic and the total energy and the kinetic temperature
 262       call kinetic(EK)
 263       totE=EK+potE
 264 c diagnostics
 265 c      call kinetic1(EK1)
 266 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
 267 c end diagnostics
 268 c Couple the system to Berendsen bath if needed
 269       if (tbf .and. lang.eq.0) then
 270         call verlet_bath
 271       endif
 272       kinetic_T=2.0d0/(dimen3*Rb)*EK
 273 c Backup the coordinates, velocities, and accelerations
 274       do i=0,2*nres
 275         do j=1,3
 276           dc_old(j,i)=dc(j,i)
 277           d_t_old(j,i)=d_t(j,i)
 278           d_a_old(j,i)=d_a(j,i)
 279         enddo
 280       enddo
 281       if (ntwe.ne.0) then
 282       if (mod(itime,ntwe).eq.0 .and. large) then
 283         write (iout,*) "Velocities, step 2"
 284         do i=0,nres
 285           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 286      &      (d_t(j,i+nres),j=1,3)
 287         enddo
 288       endif
 289       endif
 290       return
 291       end