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

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