source/unres/src_MD/md-diff/md/sd_verlet_ciccotti_setup.f

   1 c-----------------------------------------------------------
   2       subroutine sd_verlet_ciccotti_setup
   3 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
   4 c version
   5       implicit real*8 (a-h,o-z)
   6       include 'DIMENSIONS'
   7       include 'COMMON.CONTROL'
   8       include 'COMMON.VAR'
   9       include 'COMMON.MD'
  10       include 'COMMON.LANGEVIN'
  11       include 'COMMON.CHAIN'
  12       include 'COMMON.DERIV'
  13       include 'COMMON.GEO'
  14       include 'COMMON.LOCAL'
  15       include 'COMMON.INTERACT'
  16       include 'COMMON.IOUNITS'
  17       include 'COMMON.NAMES'
  18       include 'COMMON.TIME1'
  19       double precision emgdt(MAXRES6),
  20      & pterm,vterm,rho,rhoc,vsig,
  21      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
  22      & afric_vec(MAXRES6),prand_vec(MAXRES6),
  23      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
  24       logical lprn /.false./
  25       double precision zero /1.0d-8/, gdt_radius /0.05d0/
  26       double precision ktm
  27       tt0 = tcpu()
  28 c
  29 c AL 8/17/04 Code adapted from tinker
  30 c
  31 c Get the frictional and random terms for stochastic dynamics in the
  32 c eigenspace of mass-scaled UNRES friction matrix
  33 c
  34       do i = 1, dimen
  35             write (iout,*) "i",i," fricgam",fricgam(i)
  36             gdt = fricgam(i) * d_time
  37 c
  38 c Stochastic dynamics reduces to simple MD for zero friction
  39 c
  40             if (gdt .le. zero) then
  41                pfric_vec(i) = 1.0d0
  42                vfric_vec(i) = d_time
  43                afric_vec(i) = 0.5d0*d_time*d_time
  44                prand_vec(i) = afric_vec(i)
  45                vrand_vec2(i) = vfric_vec(i)
  46 c
  47 c Analytical expressions when friction coefficient is large
  48 c
  49             else
  50                egdt = dexp(-gdt)
  51                pfric_vec(i) = egdt
  52                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
  53                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
  54                prand_vec(i) = afric_vec(i)
  55                vrand_vec2(i) = vfric_vec(i)
  56 c
  57 c Compute the scaling factors of random terms for the nonzero friction case
  58 c
  59 c               ktm = 0.5d0*d_time/fricgam(i)
  60 c               psig = dsqrt(ktm*pterm) / fricgam(i)
  61 c               vsig = dsqrt(ktm*vterm)
  62 c               prand_vec(i) = psig*afric_vec(i)
  63 c               vrand_vec2(i) = vsig*vfric_vec(i)
  64             end if
  65       end do
  66       if (lprn) then
  67       write (iout,*)
  68      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
  69      &  " vrand_vec2"
  70       do i=1,dimen
  71         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
  72      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
  73       enddo
  74       endif
  75 c
  76 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
  77 c
  78       call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
  79       call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
  80       call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
  81       call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
  82       call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
  83       t_sdsetup=t_sdsetup+tcpu()-tt0
  84       return
  85       end