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)
10 include 'COMMON.CONTROL'
14 include 'COMMON.LANGEVIN'
16 include 'COMMON.LANGEVIN.lang0'
18 include 'COMMON.CHAIN'
19 include 'COMMON.DERIV'
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/
40 c AL 8/17/04 Code adapted from tinker
42 c Get the frictional and random terms for stochastic dynamics in the
43 c eigenspace of mass-scaled UNRES friction matrix
46 gdt = fricgam(i) * d_time
48 c Stochastic dynamics reduces to simple MD for zero friction
50 if (gdt .le. zero) then
53 afric_vec(i) = 0.5d0 * d_time * d_time
58 c Analytical expressions when friction coefficient is large
61 if (gdt .ge. gdt_radius) then
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)
70 c Use series expansions when friction coefficient is small
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)
103 c Compute the scaling factors of random terms for the nonzero friction case
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)
110 vrand_vec1(i) = vsig * rho
111 vrand_vec2(i) = vsig * rhoc
116 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
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)
124 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
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)
135 t_sdsetup=t_sdsetup+MPI_Wtime()
137 t_sdsetup=t_sdsetup+tcpu()-tt0