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)
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 write (iout,*) "i",i," fricgam",fricgam(i)
47 gdt = fricgam(i) * d_time
49 c Stochastic dynamics reduces to simple MD for zero friction
51 if (gdt .le. zero) then
54 afric_vec(i) = 0.5d0*d_time*d_time
55 prand_vec(i) = afric_vec(i)
56 vrand_vec2(i) = vfric_vec(i)
58 c Analytical expressions when friction coefficient is large
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)
68 c Compute the scaling factors of random terms for the nonzero friction case
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)
79 & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
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)
87 c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
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)
95 t_sdsetup=t_sdsetup+MPI_Wtime()
97 t_sdsetup=t_sdsetup+tcpu()-tt0