source/unres/src_MD/md-diff/mts/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 #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             write (iout,*) "i",i," fricgam",fricgam(i)
  47             gdt = fricgam(i) * d_time
  48 c
  49 c Stochastic dynamics reduces to simple MD for zero friction
  50 c
  51             if (gdt .le. zero) then
  52                pfric_vec(i) = 1.0d0
  53                vfric_vec(i) = d_time
  54                afric_vec(i) = 0.5d0*d_time*d_time
  55                prand_vec(i) = afric_vec(i)
  56                vrand_vec2(i) = vfric_vec(i)
  57 c
  58 c Analytical expressions when friction coefficient is large
  59 c
  60             else
  61                egdt = dexp(-gdt)
  62                pfric_vec(i) = egdt
  63                vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
  64                afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
  65                prand_vec(i) = afric_vec(i)
  66                vrand_vec2(i) = vfric_vec(i)
  67 c
  68 c Compute the scaling factors of random terms for the nonzero friction case
  69 c
  70 c               ktm = 0.5d0*d_time/fricgam(i)
  71 c               psig = dsqrt(ktm*pterm) / fricgam(i)
  72 c               vsig = dsqrt(ktm*vterm)
  73 c               prand_vec(i) = psig*afric_vec(i)
  74 c               vrand_vec2(i) = vsig*vfric_vec(i)
  75             end if
  76       end do
  77       if (lprn) then
  78       write (iout,*)
  79      &  "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
  80      &  " vrand_vec2"
  81       do i=1,dimen
  82         write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
  83      &      afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
  84       enddo
  85       endif
  86 c
  87 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
  88 c
  89       call eigtransf(dimen,maxres2,mt3,mt2,pfric_vec,pfric_mat)
  90       call eigtransf(dimen,maxres2,mt3,mt2,vfric_vec,vfric_mat)
  91       call eigtransf(dimen,maxres2,mt3,mt2,afric_vec,afric_mat)
  92       call eigtransf(dimen,maxres2,mt3,mt1,prand_vec,prand_mat)
  93       call eigtransf(dimen,maxres2,mt3,mt1,vrand_vec2,vrand_mat2)
  94 #ifdef MPI
  95       t_sdsetup=t_sdsetup+MPI_Wtime()
  96 #else
  97       t_sdsetup=t_sdsetup+tcpu()-tt0
  98 #endif
  99       return
 100       end