source/unres/src_MD/src/md-diff/np/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       double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6)
  34       double precision vtnp_(maxres6),vtnp_a(maxres6)
  35 c
  36       scale=.true.
  37       icount_scale=0
  38       if (lang.eq.1) then
  39         call sddir_precalc
  40       else if (lang.eq.2 .or. lang.eq.3) then
  41         call stochastic_force(stochforcvec)
  42       endif
  43       itime_scal=0
  44       do while (scale)
  45         icount_scale=icount_scale+1
  46         if (icount_scale.gt.maxcount_scale) then
  47           write (iout,*)
  48      &      "ERROR: too many attempts at scaling down the time step. ",
  49      &      "amax=",amax,"epdrift=",epdrift,
  50      &      "damax=",damax,"edriftmax=",edriftmax,
  51      &      "d_time=",d_time
  52             stop
  53         endif
  54 c First step of the velocity Verlet algorithm
  55         if (lang.eq.2) then
  56           call sd_verlet1
  57         else if (lang.eq.3) then
  58           call sd_verlet1_ciccotti
  59         else if (lang.eq.1) then
  60           call sddir_verlet1
  61         else if (tnp1) then
  62           call tnp1_step1
  63         else if (tnp) then
  64           call tnp_step1
  65         else
  66           if (tnh) then
  67
  68             call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
  69
  70             do i=0,2*nres
  71              do j=1,3
  72               d_t_old(j,i)=d_t_old(j,i)*scale_nh
  73              enddo
  74             enddo
  75           endif
  76           call verlet1
  77         endif
  78 c Build the chain from the newly calculated coordinates
  79         call chainbuild_cart
  80         if (rattle) call rattle1
  81         if (large.and. mod(itime,ntwe).eq.0) then
  82           write (iout,*) "Cartesian and internal coordinates: step 1"
  83           call cartprint
  84           call intout
  85           write (iout,*) "Accelerations"
  86           do i=0,nres
  87             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
  88      &      (d_a(j,i+nres),j=1,3)
  89           enddo
  90           write (iout,*) "Velocities, step 1"
  91           do i=0,nres
  92             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
  93      &      (d_t(j,i+nres),j=1,3)
  94           enddo
  95         endif
  96         tt0 = tcpu()
  97 c Calculate energy and forces
  98         call zerograd
  99         E_old=potE
 100         call etotal(energia)
 101         potE=energia(0)-energia(20)
 102         call cartgrad
 103 c Get the new accelerations
 104         call lagrangian
 105         t_enegrad=t_enegrad+tcpu()-tt0
 106 c Determine maximum acceleration and scale down the timestep if needed
 107         call max_accel
 108         call predict_edrift(epdrift)
 109         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
 110 c Maximum acceleration or maximum predicted energy drift exceeded, rescale the time step
 111           scale=.true.
 112           ifac_time=dmax1(dlog(amax/damax),dlog(epdrift/edriftmax))
 113      &      /dlog(2.0d0)+1
 114           itime_scal=itime_scal+ifac_time
 115 c          fac_time=dmin1(damax/amax,0.5d0)
 116           fac_time=0.5d0**ifac_time
 117           d_time=d_time*fac_time
 118           if (lang.eq.2 .or. lang.eq.3) then
 119 c            write (iout,*) "Calling sd_verlet_setup: 1"
 120 c Rescale the stochastic forces and recalculate or restore
 121 c the matrices of tinker integrator
 122             if (itime_scal.gt.maxflag_stoch) then
 123               if (large) write (iout,'(a,i5,a)')
 124      &         "Calculate matrices for stochastic step;",
 125      &         " itime_scal ",itime_scal
 126               if (lang.eq.2) then
 127                 call sd_verlet_p_setup
 128               else
 129                 call sd_verlet_ciccotti_setup
 130               endif
 131               write (iout,'(2a,i3,a,i3,1h.)')
 132      &         "Warning: cannot store matrices for stochastic",
 133      &         " integration because the index",itime_scal,
 134      &         " is greater than",maxflag_stoch
 135               write (iout,'(2a)')"Increase MAXFLAG_STOCH or use direct",
 136      &         " integration Langevin algorithm for better efficiency."
 137             else if (flag_stoch(itime_scal)) then
 138               if (large) write (iout,'(a,i5,a,l1)')
 139      &         "Restore matrices for stochastic step; itime_scal ",
 140      &         itime_scal," flag ",flag_stoch(itime_scal)
 141               do i=1,dimen
 142                 do j=1,dimen
 143                   pfric_mat(i,j)=pfric0_mat(i,j,itime_scal)
 144                   afric_mat(i,j)=afric0_mat(i,j,itime_scal)
 145                   vfric_mat(i,j)=vfric0_mat(i,j,itime_scal)
 146                   prand_mat(i,j)=prand0_mat(i,j,itime_scal)
 147                   vrand_mat1(i,j)=vrand0_mat1(i,j,itime_scal)
 148                   vrand_mat2(i,j)=vrand0_mat2(i,j,itime_scal)
 149                 enddo
 150               enddo
 151             else
 152               if (large) write (iout,'(2a,i5,a,l1)')
 153      &         "Calculate & store matrices for stochastic step;",
 154      &         " itime_scal ",itime_scal," flag ",flag_stoch(itime_scal)
 155               if (lang.eq.2) then
 156                 call sd_verlet_p_setup
 157               else
 158                 call sd_verlet_ciccotti_setup
 159               endif
 160               flag_stoch(ifac_time)=.true.
 161               do i=1,dimen
 162                 do j=1,dimen
 163                   pfric0_mat(i,j,itime_scal)=pfric_mat(i,j)
 164                   afric0_mat(i,j,itime_scal)=afric_mat(i,j)
 165                   vfric0_mat(i,j,itime_scal)=vfric_mat(i,j)
 166                   prand0_mat(i,j,itime_scal)=prand_mat(i,j)
 167                   vrand0_mat1(i,j,itime_scal)=vrand_mat1(i,j)
 168                   vrand0_mat2(i,j,itime_scal)=vrand_mat2(i,j)
 169                 enddo
 170               enddo
 171             endif
 172             fac_time=1.0d0/dsqrt(fac_time)
 173             do i=1,dimen
 174               stochforcvec(i)=fac_time*stochforcvec(i)
 175             enddo
 176           else if (lang.eq.1) then
 177 c Rescale the accelerations due to stochastic forces
 178             fac_time=1.0d0/dsqrt(fac_time)
 179             do i=1,dimen
 180               d_as_work(i)=d_as_work(i)*fac_time
 181             enddo
 182           endif
 183           if (large) write (iout,'(a,i10,a,f8.6,a,i3,a,i3)')
 184      &      "itime",itime," Timestep scaled down to ",
 185      &      d_time," ifac_time",ifac_time," itime_scal",itime_scal
 186         else
 187 c Second step of the velocity Verlet algorithm
 188           if (lang.eq.2) then
 189             call sd_verlet2
 190           else if (lang.eq.3) then
 191             call sd_verlet2_ciccotti
 192           else if (lang.eq.1) then
 193             call sddir_verlet2
 194           else if (tnp1) then
 195             call tnp1_step2
 196           else if (tnp) then
 197             call tnp_step2
 198           else
 199             call verlet2
 200             if (tnh) then
 201               call kinetic(EK)
 202               call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8)
 203               do i=0,2*nres
 204                do j=1,3
 205                 d_t(j,i)=d_t(j,i)*scale_nh
 206                enddo
 207               enddo
 208             endif
 209           endif
 210           if (rattle) call rattle2
 211           totT=totT+d_time
 212           if (d_time.ne.d_time0) then
 213             d_time=d_time0
 214             if (lang.eq.2 .or. lang.eq.3) then
 215               if (large) write (iout,'(a)')
 216      &         "Restore original matrices for stochastic step"
 217 c              write (iout,*) "Calling sd_verlet_setup: 2"
 218 c Restore the matrices of tinker integrator if the time step has been restored
 219               do i=1,dimen
 220                 do j=1,dimen
 221                   pfric_mat(i,j)=pfric0_mat(i,j,0)
 222                   afric_mat(i,j)=afric0_mat(i,j,0)
 223                   vfric_mat(i,j)=vfric0_mat(i,j,0)
 224                   prand_mat(i,j)=prand0_mat(i,j,0)
 225                   vrand_mat1(i,j)=vrand0_mat1(i,j,0)
 226                   vrand_mat2(i,j)=vrand0_mat2(i,j,0)
 227                 enddo
 228               enddo
 229             endif
 230           endif
 231           scale=.false.
 232         endif
 233       enddo
 234 c Calculate the kinetic and the total energy and the kinetic temperature
 235       if (tnp .or. tnp1) then
 236        do i=0,2*nres
 237         do j=1,3
 238           d_t_old(j,i)=d_t(j,i)
 239           d_t(j,i)=d_t(j,i)/s_np
 240         enddo
 241        enddo
 242       endif
 243       call kinetic(EK)
 244       totE=EK+potE
 245 c diagnostics
 246 c      call kinetic1(EK1)
 247 c      write (iout,*) "step",itime," EK",EK," EK1",EK1
 248 c end diagnostics
 249 c Couple the system to Berendsen bath if needed
 250       if (tbf .and. lang.eq.0) then
 251         call verlet_bath
 252       endif
 253       kinetic_T=2.0d0/(dimen*Rb)*EK
 254 c Backup the coordinates, velocities, and accelerations
 255       do i=0,2*nres
 256         do j=1,3
 257           dc_old(j,i)=dc(j,i)
 258           if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i)
 259           d_a_old(j,i)=d_a(j,i)
 260         enddo
 261       enddo
 262 c      if (mod(itime,ntwe).eq.0 .and. large) then
 263       if (mod(itime,ntwe).eq.0) then
 264
 265        if(tnp .or. tnp1) then
 266         HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
 267         H=(HNose1-H0)*s_np
 268 cd        write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0
 269 cd     &   ,EK+potE+pi_np**2/(2*Q_np)+dimen*0.001986d0*t_bath*log(s_np)
 270         write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
 271        endif
 272
 273        if(tnh) then
 274         HNose1=Hnose_nh(EK,potE)
 275         H=HNose1-H0
 276         write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
 277        endif
 278
 279        if (large) then
 280         itnp=0
 281         do j=1,3
 282          itnp=itnp+1
 283          vtnp(itnp)=d_t(j,0)
 284         enddo
 285         do i=nnt,nct-1
 286          do j=1,3
 287           itnp=itnp+1
 288           vtnp(itnp)=d_t(j,i)
 289          enddo
 290         enddo
 291         do i=nnt,nct
 292          if (itype(i).ne.10) then
 293           inres=i+nres
 294           do j=1,3
 295            itnp=itnp+1
 296            vtnp(itnp)=d_t(j,inres)
 297           enddo
 298          endif
 299         enddo
 300
 301 c Transform velocities from UNRES coordinate space to cartesian and Gvec
 302 c eigenvector space
 303
 304         do i=1,dimen
 305           vtnp_(i)=0.0d0
 306           vtnp_a(i)=0.0d0
 307           do j=1,dimen
 308             vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j)
 309             vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j)
 310           enddo
 311           vtnp_(i)=vtnp_(i)*dsqrt(geigen(i))
 312         enddo
 313
 314         do i=1,dimen
 315          write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
 316         enddo
 317
 318        endif
 319       endif
 320       return
 321       end