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)
6 include 'COMMON.CONTROL'
9 include 'COMMON.LANGEVIN'
10 include 'COMMON.CHAIN'
11 include 'COMMON.DERIV'
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/
28 c AL 8/17/04 Code adapted from tinker
30 c Get the frictional and random terms for stochastic dynamics in the
31 c eigenspace of mass-scaled UNRES friction matrix
34 gdt = fricgam(i) * d_time
36 c Stochastic dynamics reduces to simple MD for zero friction
38 if (gdt .le. zero) then
41 afric_vec(i) = 0.5d0 * d_time * d_time
46 c Analytical expressions when friction coefficient is large
49 if (gdt .ge. gdt_radius) then
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)
58 c Use series expansions when friction coefficient is small
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)
91 c Compute the scaling factors of random terms for the nonzero friction case
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)
98 vrand_vec1(i) = vsig * rho
99 vrand_vec2(i) = vsig * rhoc
104 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
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)
112 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
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