source/unres/src_MD/src/md-diff/md/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       include 'COMMON.LANGEVIN'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.DERIV'
  12       include 'COMMON.GEO'
  13       include 'COMMON.LOCAL'
  14       include 'COMMON.INTERACT'
  15       include 'COMMON.IOUNITS'
  16       include 'COMMON.NAMES'
  17 c Revised 3/31/05 AL: correlation between random contributions to
  18 c position and velocity increments included.
  19       double precision sqrt13 /0.57735026918962576451d0/ ! 1/sqrt(3)
  20       double precision adt,adt2
  21 c
  22 c Add the contribution from BOTH friction and stochastic force to the
  23 c coordinates, but ONLY the contribution from the friction forces to velocities
  24 c
  25       do j=1,3
  26         adt=(d_a_old(j,0)+d_af_work(j))*d_time
  27         adt2=0.5d0*adt+sqrt13*d_as_work(j)*d_time
  28         dc(j,0)=dc_old(j,0)+(d_t_old(j,0)+adt2)*d_time
  29         d_t_new(j,0)=d_t_old(j,0)+0.5d0*adt
  30         d_t(j,0)=d_t_old(j,0)+adt
  31       enddo
  32       ind=3
  33       do i=nnt,nct-1
  34         do j=1,3
  35           adt=(d_a_old(j,i)+d_af_work(ind+j))*d_time
  36           adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
  37           dc(j,i)=dc_old(j,i)+(d_t_old(j,i)+adt2)*d_time
  38           d_t_new(j,i)=d_t_old(j,i)+0.5d0*adt
  39           d_t(j,i)=d_t_old(j,i)+adt
  40         enddo
  41         ind=ind+3
  42       enddo
  43       do i=nnt,nct
  44         if (itype(i).ne.10) then
  45           inres=i+nres
  46           do j=1,3
  47             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
  48             adt2=0.5d0*adt+sqrt13*d_as_work(ind+j)*d_time
  49             dc(j,inres)=dc_old(j,inres)+(d_t_old(j,inres)+adt2)*d_time
  50             d_t_new(j,inres)=d_t_old(j,inres)+0.5d0*adt
  51             d_t(j,inres)=d_t_old(j,inres)+adt
  52           enddo
  53           ind=ind+3
  54         endif
  55       enddo
  56       return
  57       end