source/unres/src_MD/md-diff/mts/sddir_verlet1.f

   1 c---------------------------------------------------------------------
   2       subroutine sddir_verlet1
   3 c Applying velocity Verlet algorithm - step 1 to velocities
   4       implicit real*8 (a-h,o-z)
   5       include 'DIMENSIONS'
   6       include 'COMMON.CONTROL'
   7       include 'COMMON.VAR'
   8       include 'COMMON.MD'
   9 #ifndef LANG0
  10       include 'COMMON.LANGEVIN'
  11 #else
  12       include 'COMMON.LANGEVIN.lang0'
  13 #endif
  14       include 'COMMON.CHAIN'
  15       include 'COMMON.DERIV'
  16       include 'COMMON.GEO'
  17       include 'COMMON.LOCAL'
  18       include 'COMMON.INTERACT'
  19       include 'COMMON.IOUNITS'
  20       include 'COMMON.NAMES'
  21 c Revised 3/31/05 AL: correlation between random contributions to
  22 c position and velocity increments included.
  23       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
  24       double precision adt,adt2
  25 c
  26 c Add the contribution from BOTH friction and stochastic force to the
  27 c coordinates, but ONLY the contribution from the friction forces to velocities
  28 c
  29       do j=1,3
  30         adt=(d_a_old(j,0)+d_af_work(j))*d_time
  31         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
  32         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
  33         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
  34         d_t(j,0)=d_t_old(j,0)+adt
  35       enddo
  36       ind=3
  37       do i=nnt,nct-1
  38         do j=1,3
  39           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
  40           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
  41           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
  42           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
  43           d_t(j,i)=d_t_old(j,i)+adt
  44         enddo
  45         ind=ind+3
  46       enddo
  47       do i=nnt,nct
  48         if (itype(i).ne.10) then
  49           inres=i+nres
  50           do j=1,3
  51             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
  52             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
  53             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
  54             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
  55             d_t(j,inres)=d_t_old(j,inres)+adt
  56           enddo
  57           ind=ind+3
  58         endif
  59       enddo
  60       return
  61       end