source/unres/src_MD/src/md-diff/mts/respa_step.f

   1 c-------------------------------------------------------------------------------
   2       subroutine RESPA_step(itime)
   3 c-------------------------------------------------------------------------------
   4 c  Perform a single RESPA step.
   5 c-------------------------------------------------------------------------------
   6       implicit real*8 (a-h,o-z)
   7       include 'DIMENSIONS'
   8 #ifdef MPI
   9       include 'mpif.h'
  10       integer IERROR,ERRCODE
  11 #endif
  12       include 'COMMON.SETUP'
  13       include 'COMMON.CONTROL'
  14       include 'COMMON.VAR'
  15       include 'COMMON.MD'
  16 #ifndef LANG0
  17       include 'COMMON.LANGEVIN'
  18 #else
  19       include 'COMMON.LANGEVIN.lang0'
  20 #endif
  21       include 'COMMON.CHAIN'
  22       include 'COMMON.DERIV'
  23       include 'COMMON.GEO'
  24       include 'COMMON.LOCAL'
  25       include 'COMMON.INTERACT'
  26       include 'COMMON.IOUNITS'
  27       include 'COMMON.NAMES'
  28       include 'COMMON.TIME1'
  29       double precision energia_short(0:n_ene),
  30      & energia_long(0:n_ene)
  31       double precision cm(3),L(3),vcm(3),incr(3)
  32       double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
  33      & d_a_old0(3,0:maxres2)
  34       integer ilen,count,rstcount
  35       external ilen
  36       character*50 tytul
  37       integer maxcount_scale /10/
  38       common /gucio/ cm
  39       double precision stochforcvec(MAXRES6)
  40       common /stochcalc/ stochforcvec
  41       integer itime
  42       logical scale
  43       common /cipiszcze/ itt
  44       itt=itime
  45       if (ntwe.ne.0) then
  46       if (large.and. mod(itime,ntwe).eq.0) then
  47         write (iout,*) "***************** RESPA itime",itime
  48         write (iout,*) "Cartesian and internal coordinates: step 0"
  49 c        call cartprint
  50         call pdbout(0.0d0,"cipiszcze",iout)
  51         call intout
  52         write (iout,*) "Accelerations from long-range forces"
  53         do i=0,nres
  54           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
  55      &      (d_a(j,i+nres),j=1,3)
  56         enddo
  57         write (iout,*) "Velocities, step 0"
  58         do i=0,nres
  59           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
  60      &      (d_t(j,i+nres),j=1,3)
  61         enddo
  62       endif
  63       endif
  64 c
  65 c Perform the initial RESPA step (increment velocities)
  66 c      write (iout,*) "*********************** RESPA ini"
  67       call RESPA_vel
  68       if (ntwe.ne.0) then
  69       if (mod(itime,ntwe).eq.0 .and. large) then
  70         write (iout,*) "Velocities, end"
  71         do i=0,nres
  72           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
  73      &      (d_t(j,i+nres),j=1,3)
  74         enddo
  75       endif
  76       endif
  77 c Compute the short-range forces
  78 #ifdef MPI
  79       tt0 =MPI_Wtime()
  80 #else
  81       tt0 = tcpu()
  82 #endif
  83 C 7/2/2009 commented out
  84 c      call zerograd
  85 c      call etotal_short(energia_short)
  86 c      call cartgrad
  87 c      call lagrangian
  88 C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
  89         do i=0,2*nres
  90           do j=1,3
  91             d_a(j,i)=d_a_short(j,i)
  92           enddo
  93         enddo
  94       if (ntwe.ne.0) then
  95       if (large.and. mod(itime,ntwe).eq.0) then
  96         write (iout,*) "energia_short",energia_short(0)
  97         write (iout,*) "Accelerations from short-range forces"
  98         do i=0,nres
  99           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 100      &      (d_a(j,i+nres),j=1,3)
 101         enddo
 102       endif
 103       endif
 104 #ifdef MPI
 105         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 106 #else
 107         t_enegrad=t_enegrad+tcpu()-tt0
 108 #endif
 109       do i=0,2*nres
 110         do j=1,3
 111           dc_old(j,i)=dc(j,i)
 112           d_t_old(j,i)=d_t(j,i)
 113           d_a_old(j,i)=d_a(j,i)
 114         enddo
 115       enddo
 116 c 6/30/08 A-MTS: attempt at increasing the split number
 117       do i=0,2*nres
 118         do j=1,3
 119           dc_old0(j,i)=dc_old(j,i)
 120           d_t_old0(j,i)=d_t_old(j,i)
 121           d_a_old0(j,i)=d_a_old(j,i)
 122         enddo
 123       enddo
 124       if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
 125       if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
 126 c
 127       scale=.true.
 128       d_time0=d_time
 129       do while (scale)
 130
 131       scale=.false.
 132 c      write (iout,*) "itime",itime," ntime_split",ntime_split
 133 c Split the time step
 134       d_time=d_time0/ntime_split
 135 c Perform the short-range RESPA steps (velocity Verlet increments of
 136 c positions and velocities using short-range forces)
 137 c      write (iout,*) "*********************** RESPA split"
 138       do itsplit=1,ntime_split
 139         if (lang.eq.1) then
 140           call sddir_precalc
 141         else if (lang.eq.2 .or. lang.eq.3) then
 142 #ifndef LANG0
 143           call stochastic_force(stochforcvec)
 144 #else
 145           write (iout,*)
 146      &      "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
 147 #ifdef MPI
 148           call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 149 #endif
 150           stop
 151 #endif
 152         endif
 153 c First step of the velocity Verlet algorithm
 154         if (lang.eq.2) then
 155 #ifndef LANG0
 156           call sd_verlet1
 157 #endif
 158         else if (lang.eq.3) then
 159 #ifndef LANG0
 160           call sd_verlet1_ciccotti
 161 #endif
 162         else if (lang.eq.1) then
 163           call sddir_verlet1
 164         else
 165           call verlet1
 166         endif
 167 c Build the chain from the newly calculated coordinates
 168         call chainbuild_cart
 169         if (rattle) call rattle1
 170         if (ntwe.ne.0) then
 171         if (large.and. mod(itime,ntwe).eq.0) then
 172           write (iout,*) "***** ITSPLIT",itsplit
 173           write (iout,*) "Cartesian and internal coordinates: step 1"
 174           call pdbout(0.0d0,"cipiszcze",iout)
 175 c          call cartprint
 176           call intout
 177           write (iout,*) "Velocities, step 1"
 178           do i=0,nres
 179             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 180      &        (d_t(j,i+nres),j=1,3)
 181           enddo
 182         endif
 183         endif
 184 #ifdef MPI
 185         tt0 = MPI_Wtime()
 186 #else
 187         tt0 = tcpu()
 188 #endif
 189 c Calculate energy and forces
 190         call zerograd
 191         call etotal_short(energia_short)
 192 #ifdef TIMING_ENE
 193 #ifdef MPI
 194         t_eshort=t_eshort+MPI_Wtime()-tt0
 195 #else
 196         t_eshort=t_eshort+tcpu()-tt0
 197 #endif
 198 #endif
 199         call cartgrad
 200 c Get the new accelerations
 201         call lagrangian
 202 C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
 203         do i=0,2*nres
 204           do j=1,3
 205             d_a_short(j,i)=d_a(j,i)
 206           enddo
 207         enddo
 208         if (ntwe.ne.0) then
 209         if (large.and. mod(itime,ntwe).eq.0) then
 210           write (iout,*)"energia_short",energia_short(0)
 211           write (iout,*) "Accelerations from short-range forces"
 212           do i=0,nres
 213             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 214      &        (d_a(j,i+nres),j=1,3)
 215           enddo
 216         endif
 217         endif
 218 c 6/30/08 A-MTS
 219 c Determine maximum acceleration and scale down the timestep if needed
 220         call max_accel
 221         amax=amax/ntime_split**2
 222         call predict_edrift(epdrift)
 223         if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
 224      &   write (iout,*) "amax",amax," damax",damax,
 225      &   " epdrift",epdrift," epdriftmax",epdriftmax
 226 c Exit loop and try with increased split number if the change of
 227 c acceleration is too big
 228         if (amax.gt.damax .or. epdrift.gt.edriftmax) then
 229           if (ntime_split.lt.maxtime_split) then
 230             scale=.true.
 231             ntime_split=ntime_split*2
 232             do i=0,2*nres
 233               do j=1,3
 234                 dc_old(j,i)=dc_old0(j,i)
 235                 d_t_old(j,i)=d_t_old0(j,i)
 236                 d_a_old(j,i)=d_a_old0(j,i)
 237               enddo
 238             enddo
 239             write (iout,*) "acceleration/energy drift too large",amax,
 240      &      epdrift," split increased to ",ntime_split," itime",itime,
 241      &       " itsplit",itsplit
 242             exit
 243           else
 244             write (iout,*)
 245      &      "Uh-hu. Bumpy landscape. Maximum splitting number",
 246      &       maxtime_split,
 247      &      " already reached!!! Trying to carry on!"
 248           endif
 249         endif
 250 #ifdef MPI
 251         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 252 #else
 253         t_enegrad=t_enegrad+tcpu()-tt0
 254 #endif
 255 c Second step of the velocity Verlet algorithm
 256         if (lang.eq.2) then
 257 #ifndef LANG0
 258           call sd_verlet2
 259 #endif
 260         else if (lang.eq.3) then
 261 #ifndef LANG0
 262           call sd_verlet2_ciccotti
 263 #endif
 264         else if (lang.eq.1) then
 265           call sddir_verlet2
 266         else
 267           call verlet2
 268         endif
 269         if (rattle) call rattle2
 270 c Backup the coordinates, velocities, and accelerations
 271         do i=0,2*nres
 272           do j=1,3
 273             dc_old(j,i)=dc(j,i)
 274             d_t_old(j,i)=d_t(j,i)
 275             d_a_old(j,i)=d_a(j,i)
 276           enddo
 277         enddo
 278       enddo
 279
 280       enddo ! while scale
 281
 282 c Restore the time step
 283       d_time=d_time0
 284 c Compute long-range forces
 285 #ifdef MPI
 286       tt0 =MPI_Wtime()
 287 #else
 288       tt0 = tcpu()
 289 #endif
 290       call zerograd
 291       call etotal_long(energia_long)
 292 #ifdef TIMING_ENE
 293 #ifdef MPI
 294         t_elong=t_elong+MPI_Wtime()-tt0
 295 #else
 296         t_elong=t_elong+tcpu()-tt0
 297 #endif
 298 #endif
 299       call cartgrad
 300       call lagrangian
 301 #ifdef MPI
 302         t_enegrad=t_enegrad+MPI_Wtime()-tt0
 303 #else
 304         t_enegrad=t_enegrad+tcpu()-tt0
 305 #endif
 306 c Compute accelerations from long-range forces
 307       if (ntwe.ne.0) then
 308       if (large.and. mod(itime,ntwe).eq.0) then
 309         write (iout,*) "energia_long",energia_long(0)
 310         write (iout,*) "Cartesian and internal coordinates: step 2"
 311 c        call cartprint
 312         call pdbout(0.0d0,"cipiszcze",iout)
 313         call intout
 314         write (iout,*) "Accelerations from long-range forces"
 315         do i=0,nres
 316           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 317      &      (d_a(j,i+nres),j=1,3)
 318         enddo
 319         write (iout,*) "Velocities, step 2"
 320         do i=0,nres
 321           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 322      &      (d_t(j,i+nres),j=1,3)
 323         enddo
 324       endif
 325       endif
 326 c Compute the final RESPA step (increment velocities)
 327 c      write (iout,*) "*********************** RESPA fin"
 328       call RESPA_vel
 329 c Compute the complete potential energy
 330       do i=0,n_ene
 331         potEcomp(i)=energia_short(i)+energia_long(i)
 332       enddo
 333       potE=potEcomp(0)-potEcomp(20)
 334 c      potE=energia_short(0)+energia_long(0)
 335       totT=totT+d_time
 336 c Calculate the kinetic and the total energy and the kinetic temperature
 337       call kinetic(EK)
 338       totE=EK+potE
 339 c Couple the system to Berendsen bath if needed
 340       if (tbf .and. lang.eq.0) then
 341         call verlet_bath
 342       endif
 343       kinetic_T=2.0d0/(dimen3*Rb)*EK
 344 c Backup the coordinates, velocities, and accelerations
 345       if (ntwe.ne.0) then
 346       if (mod(itime,ntwe).eq.0 .and. large) then
 347         write (iout,*) "Velocities, end"
 348         do i=0,nres
 349           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 350      &      (d_t(j,i+nres),j=1,3)
 351         enddo
 352       endif
 353       endif
 354       return
 355       end