+++ /dev/null
-c-----------------------------------------------------------
- subroutine sd_verlet_p_setup
-c Sets up the parameters of stochastic Verlet algorithm
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
- include 'COMMON.LANGEVIN'
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.NAMES'
- include 'COMMON.TIME1'
- double precision emgdt(MAXRES6),
- & pterm,vterm,rho,rhoc,vsig,
- & pfric_vec(MAXRES6),vfric_vec(MAXRES6),
- & afric_vec(MAXRES6),prand_vec(MAXRES6),
- & vrand_vec1(MAXRES6),vrand_vec2(MAXRES6)
- logical lprn /.false./
- double precision zero /1.0d-8/, gdt_radius /0.05d0/
- double precision ktm
- tt0 = tcpu()
-c
-c AL 8/17/04 Code adapted from tinker
-c
-c Get the frictional and random terms for stochastic dynamics in the
-c eigenspace of mass-scaled UNRES friction matrix
-c
- do i = 1, dimen
- gdt = fricgam(i) * d_time
-c
-c Stochastic dynamics reduces to simple MD for zero friction
-c
- if (gdt .le. zero) then
- pfric_vec(i) = 1.0d0
- vfric_vec(i) = d_time
- afric_vec(i) = 0.5d0 * d_time * d_time
- prand_vec(i) = 0.0d0
- vrand_vec1(i) = 0.0d0
- vrand_vec2(i) = 0.0d0
-c
-c Analytical expressions when friction coefficient is large
-c
- else
- if (gdt .ge. gdt_radius) then
- egdt = dexp(-gdt)
- pfric_vec(i) = egdt
- vfric_vec(i) = (1.0d0-egdt) / fricgam(i)
- afric_vec(i) = (d_time-vfric_vec(i)) / fricgam(i)
- pterm = 2.0d0*gdt - 3.0d0 + (4.0d0-egdt)*egdt
- vterm = 1.0d0 - egdt**2
- rho = (1.0d0-egdt)**2 / sqrt(pterm*vterm)
-c
-c Use series expansions when friction coefficient is small
-c
- else
- gdt2 = gdt * gdt
- gdt3 = gdt * gdt2
- gdt4 = gdt2 * gdt2
- gdt5 = gdt2 * gdt3
- gdt6 = gdt3 * gdt3
- gdt7 = gdt3 * gdt4
- gdt8 = gdt4 * gdt4
- gdt9 = gdt4 * gdt5
- afric_vec(i) = (gdt2/2.0d0 - gdt3/6.0d0 + gdt4/24.0d0
- & - gdt5/120.0d0 + gdt6/720.0d0
- & - gdt7/5040.0d0 + gdt8/40320.0d0
- & - gdt9/362880.0d0) / fricgam(i)**2
- vfric_vec(i) = d_time - fricgam(i)*afric_vec(i)
- pfric_vec(i) = 1.0d0 - fricgam(i)*vfric_vec(i)
- pterm = 2.0d0*gdt3/3.0d0 - gdt4/2.0d0
- & + 7.0d0*gdt5/30.0d0 - gdt6/12.0d0
- & + 31.0d0*gdt7/1260.0d0 - gdt8/160.0d0
- & + 127.0d0*gdt9/90720.0d0
- vterm = 2.0d0*gdt - 2.0d0*gdt2 + 4.0d0*gdt3/3.0d0
- & - 2.0d0*gdt4/3.0d0 + 4.0d0*gdt5/15.0d0
- & - 4.0d0*gdt6/45.0d0 + 8.0d0*gdt7/315.0d0
- & - 2.0d0*gdt8/315.0d0 + 4.0d0*gdt9/2835.0d0
- rho = sqrt(3.0d0) * (0.5d0 - 3.0d0*gdt/16.0d0
- & - 17.0d0*gdt2/1280.0d0
- & + 17.0d0*gdt3/6144.0d0
- & + 40967.0d0*gdt4/34406400.0d0
- & - 57203.0d0*gdt5/275251200.0d0
- & - 1429487.0d0*gdt6/13212057600.0d0)
- end if
-c
-c Compute the scaling factors of random terms for the nonzero friction case
-c
- ktm = 0.5d0*d_time/fricgam(i)
- psig = dsqrt(ktm*pterm) / fricgam(i)
- vsig = dsqrt(ktm*vterm)
- rhoc = dsqrt(1.0d0 - rho*rho)
- prand_vec(i) = psig
- vrand_vec1(i) = vsig * rho
- vrand_vec2(i) = vsig * rhoc
- end if
- end do
- if (lprn) then
- write (iout,*)
- & "pfric_vec, vfric_vec, afric_vec, prand_vec, vrand_vec1,",
- & " vrand_vec2"
- do i=1,dimen
- write (iout,'(i5,6e15.5)') i,pfric_vec(i),vfric_vec(i),
- & afric_vec(i),prand_vec(i),vrand_vec1(i),vrand_vec2(i)
- enddo
- endif
-c
-c Transform from the eigenspace of mass-scaled friction matrix to UNRES variables
-c
- call eigtransf(dimen,maxres6,mt3,mt2,pfric_vec,pfric_mat)
- call eigtransf(dimen,maxres6,mt3,mt2,vfric_vec,vfric_mat)
- call eigtransf(dimen,maxres6,mt3,mt2,afric_vec,afric_mat)
- call eigtransf(dimen,maxres6,mt3,mt1,prand_vec,prand_mat)
- call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec1,vrand_mat1)
- call eigtransf(dimen,maxres6,mt3,mt1,vrand_vec2,vrand_mat2)
-c call eigtransf1(dimen,maxres6,mt3mt2,pfric_vec,pfric_mat)
-c call eigtransf1(dimen,maxres6,mt3mt2,vfric_vec,vfric_mat)
-c call eigtransf1(dimen,maxres6,mt3mt2,afric_vec,afric_mat)
-c call eigtransf1(dimen,maxres6,mt3mt1,prand_vec,prand_mat)
-c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec1,vrand_mat1)
-c call eigtransf1(dimen,maxres6,mt3mt1,vrand_vec2,vrand_mat2)
- t_sdsetup=t_sdsetup+tcpu()-tt0
- return
- end
projects / unres.git / blobdiff