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

   1 c-----------------------------------------------------------
   2       subroutine sd_verlet_p_setup
   3 c Sets up the parameters of stochastic Verlet algorithm
   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       include 'COMMON.TIME1'
  18       double precision emgdt(MAXRES6),
  19      & pterm,vterm,rho,rhoc,vsig,
  20      & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
  21      & afric_vec(MAXRES6),prand_vec(MAXRES6),
  22      & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
  23       logical lprn /.false./
  24       double precision zero /1.0d-8/, gdt_radius /0.05d0/
  25       double precision ktm
  26       tt0 = tcpu()
  27 c
  28 c AL 8/17/04 Code adapted from tinker
  29 c
  30 c Get the frictional and random terms for stochastic dynamics in the
  31 c eigenspace of mass-scaled UNRES friction matrix
  32 c
  33       do i = 1, dimen
  34             gdt = fricgam(i) * d_time
  35 c
  36 c Stochastic dynamics reduces to simple MD for zero friction
  37 c
  38             if (gdt .le. zero) then
  39                pfric_vec(i) = 1.0d0
  40                vfric_vec(i) = d_time
  41                afric_vec(i) = 0.5d0 * d_time * d_time
  42                prand_vec(i) = 0.0d0
  43                vrand_vec1(i) = 0.0d0
  44                vrand_vec2(i) = 0.0d0
  45 c
  46 c Analytical expressions when friction coefficient is large
  47 c
  48             else
  49                if (gdt .ge. gdt_radius) then
  50                   egdt = dexp(-gdt)
  51                   pfric_vec(i) = egdt
  52                   vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
  53                   afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
  54                   pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
  55                   vterm = 1.0d0 - egdt**2
  56                   rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
  57 c
  58 c Use series expansions when friction coefficient is small
  59 c
  60                else
  61                   gdt2 = gdt * gdt
  62                   gdt3 = gdt * gdt2
  63                   gdt4 = gdt2 * gdt2
  64                   gdt5 = gdt2 * gdt3
  65                   gdt6 = gdt3 * gdt3
  66                   gdt7 = gdt3 * gdt4
  67                   gdt8 = gdt4 * gdt4
  68                   gdt9 = gdt4 * gdt5
  69                   afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
  70      &                          - gdt5/120.0d0 + gdt6/720.0d0
  71      &                          - gdt7/5040.0d0 + gdt8/40320.0d0
  72      &                          - gdt9/362880.0d0) / fricgam(i)**2
  73                   vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
  74                   pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
  75                   pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
  76      &                       + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
  77      &                       + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
  78      &                       + 127.0d0*gdt9/90720.0d0
  79                   vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
  80      &                       - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
  81      &                       - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
  82      &                       - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
  83                   rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
  84      &                       - 17.0d0*gdt2/1280.0d0
  85      &                       + 17.0d0*gdt3/6144.0d0
  86      &                       + 40967.0d0*gdt4/34406400.0d0
  87      &                       - 57203.0d0*gdt5/275251200.0d0
  88      &                       - 1429487.0d0*gdt6/13212057600.0d0)
  89                end if
  90 c
  91 c Compute the scaling factors of random terms for the nonzero friction case
  92 c
  93                ktm = 0.5d0*d_time/fricgam(i)
  94                psig = dsqrt(ktm*pterm) / fricgam(i)
  95                vsig = dsqrt(ktm*vterm)
  96                rhoc = dsqrt(1.0d0 - rho*rho)
  97                prand_vec(i) = psig
  98                vrand_vec1(i) = vsig * rho
  99                vrand_vec2(i) = vsig * rhoc
 100             end if
 101       end do
 102       if (lprn) then
 103       write (iout,*)
 104      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
 105      &  " vrand_vec2"
 106       do i=1,dimen
 107         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
 108      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
 109       enddo
 110       endif
 111 c
 112 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
 113 c
 114       call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
 115       call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
 116       call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
 117       call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
 118       call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1)
 119       call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
 120 c      call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat)
 121 c      call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat)
 122 c      call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat)
 123 c      call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat)
 124 c      call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1)
 125 c      call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2)
 126       t_sdsetup=t_sdsetup+tcpu()-tt0
 127       return
 128       end