1 c-----------------------------------------------------------
2 subroutine sd_verlet_ciccotti_setup
3 c Sets up the parameters of stochastic velocity Verlet algorithmi; Ciccotti's
5 implicit real*8 (a-h,o-z)
7 include 'COMMON.CONTROL'
10 include 'COMMON.LANGEVIN'
11 include 'COMMON.CHAIN'
12 include 'COMMON.DERIV'
14 include 'COMMON.LOCAL'
15 include 'COMMON.INTERACT'
16 include 'COMMON.IOUNITS'
17 include 'COMMON.NAMES'
18 include 'COMMON.TIME1'
19 double precision emgdt(MAXRES6),
20 & pterm,vterm,rho,rhoc,vsig,
21 & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
22 & afric_vec(MAXRES6),prand_vec(MAXRES6),
23 & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
24 logical lprn /.false./
25 double precision zero /1.0d-8/, gdt_radius /0.05d0/
29 c AL 8/17/04 Code adapted from tinker
31 c Get the frictional and random terms for stochastic dynamics in the
32 c eigenspace of mass-scaled UNRES friction matrix
35 write (iout,*) "i",i," fricgam",fricgam(i)
36 gdt = fricgam(i) * d_time
38 c Stochastic dynamics reduces to simple MD for zero friction
40 if (gdt .le. zero) then
43 afric_vec(i) = 0.5d0*d_time*d_time
44 prand_vec(i) = afric_vec(i)
45 vrand_vec2(i) = vfric_vec(i)
47 c Analytical expressions when friction coefficient is large
52 vfric_vec(i) = dexp(-0.5d0*gdt)*d_time
53 afric_vec(i) = 0.5d0*dexp(-0.25d0*gdt)*d_time*d_time
54 prand_vec(i) = afric_vec(i)
55 vrand_vec2(i) = vfric_vec(i)
57 c Compute the scaling factors of random terms for the nonzero friction case
59 c ktm = 0.5d0*d_time/fricgam(i)
60 c psig = dsqrt(ktm*pterm) / fricgam(i)
61 c vsig = dsqrt(ktm*vterm)
62 c prand_vec(i) = psig*afric_vec(i)
63 c vrand_vec2(i) = vsig*vfric_vec(i)
68 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
71 write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
72 & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
76 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
78 call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
79 call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
80 call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
81 call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
82 call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
83 t_sdsetup=t_sdsetup+tcpu()-tt0