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

   1 c---------------------------------------------------------------------
   2       subroutine sddir_verlet2
   3 c  Calculating the adjusted velocities for accelerations
   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       double precision stochforcvec(MAXRES6),d_as_work1(MAXRES6)
  22       double precision cos60 /0.5d0/, sin60 /0.86602540378443864676d0/
  23 c Revised 3/31/05 AL: correlation between random contributions to
  24 c position and velocity increments included.
  25 c The correlation coefficients are calculated at low-friction limit.
  26 c Also, friction forces are now not calculated with new velocities.
  27
  28 c      call friction_force
  29       call stochastic_force(stochforcvec)
  30 c
  31 c Compute the acceleration due to friction forces (d_af_work) and stochastic
  32 c forces (d_as_work)
  33 c
  34       call ginv_mult(stochforcvec, d_as_work1)
  35
  36 c
  37 c Update velocities
  38 c
  39       do j=1,3
  40         d_t(j,0)=d_t_new(j,0)+(0.5d0*(d_a(j,0)+d_af_work(j))
  41      &    +sin60*d_as_work(j)+cos60*d_as_work1(j))*d_time
  42       enddo
  43       ind=3
  44       do i=nnt,nct-1
  45         do j=1,3
  46           d_t(j,i)=d_t_new(j,i)+(0.5d0*(d_a(j,i)+d_af_work(ind+j))
  47      &     +sin60*d_as_work(ind+j)+cos60*d_as_work1(ind+j))*d_time
  48         enddo
  49         ind=ind+3
  50       enddo
  51       do i=nnt,nct
  52         if (itype(i).ne.10) then
  53           inres=i+nres
  54           do j=1,3
  55             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
  56      &       +d_af_work(ind+j))+sin60*d_as_work(ind+j)
  57      &       +cos60*d_as_work1(ind+j))*d_time
  58           enddo
  59           ind=ind+3
  60         endif
  61       enddo
  62       return
  63       end