source/unres/src_MD/md-diff/md/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       include 'COMMON.CONTROL'
   9       include 'COMMON.VAR'
  10       include 'COMMON.MD'
  11       include 'COMMON.LANGEVIN'
  12       include 'COMMON.CHAIN'
  13       include 'COMMON.DERIV'
  14       include 'COMMON.GEO'
  15       include 'COMMON.LOCAL'
  16       include 'COMMON.INTERACT'
  17       include 'COMMON.IOUNITS'
  18       include 'COMMON.NAMES'
  19       include 'COMMON.TIME1'
  20       double precision energia(0:n_ene),energia_short(0:n_ene),
  21      & energia_long(0:n_ene)
  22       double precision cm(3),L(3),vcm(3),incr(3)
  23       integer ilen,count,rstcount
  24       external ilen
  25       character*50 tytul
  26       integer maxcount_scale /10/
  27       common /gucio/ cm
  28       double precision stochforcvec(MAXRES6)
  29       common /stochcalc/ stochforcvec
  30       integer itime
  31       logical scale
  32       if (large.and. mod(itime,ntwe).eq.0) then
  33         write (iout,*) "***************** RESPA itime",itime
  34         write (iout,*) "Cartesian and internal coordinates: step 0"
  35         call cartprint
  36         call intout
  37         write (iout,*) "Accelerations"
  38         do i=0,nres
  39           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
  40      &      (d_a(j,i+nres),j=1,3)
  41         enddo
  42         write (iout,*) "Velocities, step 0"
  43         do i=0,nres
  44           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
  45      &      (d_t(j,i+nres),j=1,3)
  46         enddo
  47       endif
  48 c
  49 c Perform the initial RESPA step (increment velocities)
  50 c      write (iout,*) "*********************** RESPA ini"
  51       call RESPA_vel
  52       if (mod(itime,ntwe).eq.0 .and. large) then
  53         write (iout,*) "Velocities, end"
  54         do i=0,nres
  55           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
  56      &      (d_t(j,i+nres),j=1,3)
  57         enddo
  58       endif
  59 c Compute the short-range forces
  60       call zerograd
  61       call etotal_short(energia_short)
  62       call cartgrad
  63       call lagrangian
  64       do i=0,2*nres
  65         do j=1,3
  66           dc_old(j,i)=dc(j,i)
  67           d_t_old(j,i)=d_t(j,i)
  68           d_a_old(j,i)=d_a(j,i)
  69         enddo
  70       enddo
  71       d_time0=d_time
  72 c Split the time step
  73       d_time=d_time/ntime_split
  74 c Perform the short-range RESPSA steps (velocity Verlet increments of
  75 c positions and velocities using short-range forces)
  76 c      write (iout,*) "*********************** RESPA split"
  77       do itsplit=1,ntime_split
  78         if (lang.eq.1) then
  79           call sddir_precalc
  80         else if (lang.eq.2 .or. lang.eq.3) then
  81           call stochastic_force(stochforcvec)
  82         endif
  83 c First step of the velocity Verlet algorithm
  84         if (lang.eq.2) then
  85           call sd_verlet1
  86         else if (lang.eq.3) then
  87           call sd_verlet1_ciccotti
  88         else if (lang.eq.1) then
  89           call sddir_verlet1
  90         else
  91           call verlet1
  92         endif
  93 c Build the chain from the newly calculated coordinates
  94         call chainbuild_cart
  95         if (rattle) call rattle1
  96         if (large.and. mod(itime,ntwe).eq.0) then
  97           write (iout,*) "Cartesian and internal coordinates: step 1"
  98           call cartprint
  99           call intout
 100           write (iout,*) "Accelerations"
 101           do i=0,nres
 102             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 103      &        (d_a(j,i+nres),j=1,3)
 104           enddo
 105           write (iout,*) "Velocities, step 1"
 106           do i=0,nres
 107             write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 108      &        (d_t(j,i+nres),j=1,3)
 109           enddo
 110         endif
 111         tt0 = tcpu()
 112 c Calculate energy and forces
 113         call zerograd
 114         call etotal_short(energia_short)
 115         call cartgrad
 116 c Get the new accelerations
 117         call lagrangian
 118         t_enegrad=t_enegrad+tcpu()-tt0
 119 c Second step of the velocity Verlet algorithm
 120         if (lang.eq.2) then
 121           call sd_verlet2
 122         else if (lang.eq.3) then
 123           call sd_verlet2_ciccotti
 124         else if (lang.eq.1) then
 125           call sddir_verlet2
 126         else
 127           call verlet2
 128         endif
 129         if (rattle) call rattle2
 130 c Backup the coordinates, velocities, and accelerations
 131         do i=0,2*nres
 132           do j=1,3
 133             dc_old(j,i)=dc(j,i)
 134             d_t_old(j,i)=d_t(j,i)
 135             d_a_old(j,i)=d_a(j,i)
 136           enddo
 137         enddo
 138       enddo
 139 c Restore the time step
 140       d_time=d_time0
 141 c Compute long-range forces
 142       call zerograd
 143       call etotal_long(energia_long)
 144       call cartgrad
 145 c Compute accelerations from long-range forces
 146       call lagrangian
 147       if (large.and. mod(itime,ntwe).eq.0) then
 148         write (iout,*) "Cartesian and internal coordinates: step 2"
 149         call cartprint
 150         call intout
 151         write (iout,*) "Accelerations"
 152         do i=0,nres
 153           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
 154      &      (d_a(j,i+nres),j=1,3)
 155         enddo
 156         write (iout,*) "Velocities, step 2"
 157         do i=0,nres
 158           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 159      &      (d_t(j,i+nres),j=1,3)
 160         enddo
 161       endif
 162 c Compute the final RESPA step (increment velocities)
 163 c      write (iout,*) "*********************** RESPA fin"
 164       call RESPA_vel
 165 c Compute the complete potential energy
 166       potE=energia_short(0)+energia_long(0)
 167       totT=totT+d_time
 168 c Calculate the kinetic and the total energy and the kinetic temperature
 169       call kinetic(EK)
 170       totE=EK+potE
 171 c Couple the system to Berendsen bath if needed
 172       if (tbf .and. lang.eq.0) then
 173         call verlet_bath
 174       endif
 175       kinetic_T=2.0d0/(dimen*Rb)*EK
 176 c Backup the coordinates, velocities, and accelerations
 177       if (mod(itime,ntwe).eq.0 .and. large) then
 178         write (iout,*) "Velocities, end"
 179         do i=0,nres
 180           write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
 181      &      (d_t(j,i+nres),j=1,3)
 182         enddo
 183       endif
 184       return
 185       end