Skip to content

some bug fixes #13

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
dmargul merged 1 commit into from
Jun 28, 2015
Merged

some bug fixes #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 6 additions & 5 deletions source/bath.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
c
module bath
implicit none
use sizes
integer maxnose
parameter (maxnose=4)
integer len_nhc,isok_L,isok_M
Expand All @@ -58,11 +59,11 @@ module bath
real*8 compress,collide
real*8 eta,volmove
real*8 vbar,qbar,gbar
real*8 q_iso1(isok_chain,3,10000)
real*8 q_iso2(isok_chain,3,10000)
real*8 v_iso1(isok_chain,3,10000)
real*8 v_iso2(isok_chain,3,10000)
real*8 f_iso1(isok_chain,3,10000)
real*8 q_iso1(isok_chain,3,maxatm)
real*8 q_iso2(isok_chain,3,maxatm)
real*8 v_iso1(isok_chain,3,maxatm)
real*8 v_iso2(isok_chain,3,maxatm)
real*8 f_iso1(isok_chain,3,maxatm)
real*8 v_sinr1(isok_chain)
real*8 v_sinr2(isok_chain)
real*8 q_sinr1(isok_chain)
Expand Down
180 changes: 178 additions & 2 deletions source/ehal1.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1493,7 +1493,7 @@ subroutine ehalSR
c
c check for an interaction distance less than the cutoff
c
if (rik2 .le. off2) then
if (rik2 .le. Rcut2) then
rik = sqrt(rik2)
rv = radmin(kt,it)
eps = epsilon(kt,it)
Expand Down Expand Up @@ -1546,6 +1546,18 @@ subroutine ehalSR
de = e*dtaper + de*taper
e = e * taper
end if
c use RESPA quintic switching in switching domain
c
if (rik2.gt.REtaper2) then
res_u = (rik - REvdwcut + res_lam) / res_lam
ru2 = res_u * res_u
ru3 = ru2 * res_u
ru4 = ru3 * res_u
ru5 = ru4 * res_u
q_switch = 1.0d0 + (15.0d0*ru4) - (6.0d0*ru5)
& - (10.0d0*ru3)
de = de*q_switch
end if
c
c scale the interaction based on its group membership
c
Expand Down Expand Up @@ -1748,7 +1760,6 @@ subroutine ehalLR
devLR(3,i) = 0.0d0
end do

icount = 0
c
c perform dynamic allocation of some local arrays
c
Expand Down Expand Up @@ -1842,6 +1853,171 @@ subroutine ehalLR
c
c decide whether to compute the current interaction
c

c
c decide whether to compute the current interaction
c
do kk = 1, REnvlst(ii)
k = ivdw(REvlst(kk,ii))
kv = ired(k)
mutk = mut(k)
proceed = .true.
if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0)
if (proceed) proceed = (usei .or. use(k) .or. use(kv))
c
c compute the energy contribution for this interaction
c
if (proceed) then
kt = jvdw(k)
xr = xi - xred(k)
yr = yi - yred(k)
zr = zi - zred(k)
call image (xr,yr,zr)
rik2 = xr*xr + yr*yr + zr*zr
c
c check for an interaction distance less than the cutoff
c
if (rik2 .le. REcut2) then
rik = sqrt(rik2)
rv = radmin(kt,it)
eps = epsilon(kt,it)
if (iv14(k) .eq. i) then
rv = radmin4(kt,it)
eps = epsilon4(kt,it)
end if
eps = eps * vscale(k)
c
c get the energy and gradient, via soft core if necessary
c
if ((muti .and. .not.mutk) .or.
& (mutk .and. .not.muti)) then
rho = rik / rv
rho6 = rho**6
rho7 = rho6 * rho
eps = eps * vlambda**scexp
scal = scalpha * (1.0d0-vlambda)**2
s1 = 1.0d0 / (scal+(rho+dhal)**7)
s2 = 1.0d0 / (scal+rho7+ghal)
t1 = (1.0d0+dhal)**7 * s1
t2 = (1.0d0+ghal) * s2
dt1drho = -7.0d0*(rho+dhal)**6 * t1 * s1
dt2drho = -7.0d0*rho6 * t2 * s2
e = 0.0d0
de = eps * (dt1drho*(t2-2.0d0)+t1*dt2drho) / rv
else
rv7 = rv**7
rik6 = rik2**3
rik7 = rik6 * rik
rho = rik7 + ghal*rv7
tau = (dhal+1.0d0) / (rik + dhal*rv)
tau7 = tau**7
dtau = tau / (dhal+1.0d0)
gtau = eps*tau7*rik6*(ghal+1.0d0)*(rv7/rho)**2
e = 0.0d0
de = -7.0d0 * (dtau*e+gtau)
end if
c
c use energy switching if near the cutoff distance
c
if (rik2 .gt. cut2) then
rik3 = rik2 * rik
rik4 = rik2 * rik2
rik5 = rik2 * rik3
taper = c5*rik5 + c4*rik4 + c3*rik3
& + c2*rik2 + c1*rik + c0
dtaper = 5.0d0*c5*rik4 + 4.0d0*c4*rik3
& + 3.0d0*c3*rik2 + 2.0d0*c2*rik + c1
de = e*dtaper + de*taper
e = e * taper
end if
c
c use RESPA quintic switching in switching domain
c
if (rik2.gt.REtaper2) then
res_u = (rik - REvdwcut + res_lam) / res_lam
ru2 = res_u * res_u
ru3 = ru2 * res_u
ru4 = ru3 * res_u
ru5 = ru4 * res_u
q_switch = 1.0d0 + (15.0d0*ru4) - (6.0d0*ru5)
& - (10.0d0*ru3)
de = de*(1.0d0-q_switch)
else
e = 0.0d0
de = 0.0d0
end if
c
c scale the interaction based on its group membership
c
if (use_group) then
e = e * fgrp
de = de * fgrp
end if
c
c find the chain rule terms for derivative components
c
de = de / rik
dedx = de * xr
dedy = de * yr
dedz = de * zr
c
c increment the total van der Waals energy and derivatives
c
evo = evo + e
if (i .eq. iv) then
devo(1,i) = devo(1,i) + dedx
devo(2,i) = devo(2,i) + dedy
devo(3,i) = devo(3,i) + dedz
else
devo(1,i) = devo(1,i) + dedx*redi
devo(2,i) = devo(2,i) + dedy*redi
devo(3,i) = devo(3,i) + dedz*redi
devo(1,iv) = devo(1,iv) + dedx*rediv
devo(2,iv) = devo(2,iv) + dedy*rediv
devo(3,iv) = devo(3,iv) + dedz*rediv
end if
if (k .eq. kv) then
devo(1,k) = devo(1,k) - dedx
devo(2,k) = devo(2,k) - dedy
devo(3,k) = devo(3,k) - dedz
else
redk = kred(k)
redkv = 1.0d0 - redk
devo(1,k) = devo(1,k) - dedx*redk
devo(2,k) = devo(2,k) - dedy*redk
devo(3,k) = devo(3,k) - dedz*redk
devo(1,kv) = devo(1,kv) - dedx*redkv
devo(2,kv) = devo(2,kv) - dedy*redkv
devo(3,kv) = devo(3,kv) - dedz*redkv
end if
c
c increment the internal virial tensor components
c
vxx = xr * dedx
vyx = yr * dedx
vzx = zr * dedx
vyy = yr * dedy
vzy = zr * dedy
vzz = zr * dedz
viro(1,1) = viro(1,1) + vxx
viro(2,1) = viro(2,1) + vyx
viro(3,1) = viro(3,1) + vzx
viro(1,2) = viro(1,2) + vyx
viro(2,2) = viro(2,2) + vyy
viro(3,2) = viro(3,2) + vzy
viro(1,3) = viro(1,3) + vzx
viro(2,3) = viro(2,3) + vzy
viro(3,3) = viro(3,3) + vzz
c
c increment the total intermolecular energy
c
if (molcule(i) .ne. molcule(k)) then
eintero = eintero + e
end if
end if
end if
end do

do kk = 1, nvlst(ii)
icount = icount + 1
k = ivdw(vlst(kk,ii))
Expand Down
48 changes: 42 additions & 6 deletions source/stoch_respa2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -262,25 +262,61 @@ subroutine applyisok (dt)
real*8 dt,dtc,dts
real*8 dt2,dt4,dt8
real*8 tempstor,kBT,LkT
real*8 w(3)
real*8 w(respa_therm_nsy)
real*8 len_fac

len_fac = dble(len_nhc)/dble(len_nhc+1)

kBT = boltzmann*kelvin
LkT = (dble(len_nhc))*kBT

ns = 3
dtc = dt / dble(respa_therm_nc)
w(1) = 1.0d0 / (2.0d0-2.0d0**(1.0d0/3.0d0))
w(2) = 1.0d0 - 2.0d0*w(1)
w(3) = w(1)
if (respa_therm_nsy .eq. 3) then
w(1) = 1.0d0 / (2.0d0-2.0d0**(1.0d0/3.0d0))
w(2) = 1.0d0 - 2.0d0*w(1)
w(3) = w(1)
else if (respa_therm_nsy .eq. 5) then
w(1) = 1.0d0 / (4.0d0-4.0d0**(1.0d0/3.0d0))
w(2) = w(1)
w(4) = w(1)
w(5) = w(1)
w(3) = 1 - (4.0d0 * w(1))
else if (respa_therm_nsy .eq. 7) then
w(1) = 0.784513610477560
w(2) = 0.235573213359357
w(3) = -1.17767998417887
w(5) = w(3)
w(6) = w(2)
w(7) = w(1)
w(4) = 1.0d0 - 2.0d0*(w(1)+w(2)+w(3))
else if (respa_therm_nsy.eq.15) then
w(1) = 0.102799849391985
w(2) = -1.96061023297549
w(3) = 1.93813913762276
w(4) = -0.158240635368243
w(5) = -1.44485223686048
w(6) = 0.253693336566229
w(7) = 0.914844246229740
w(8) = 1.0d0 - 2.0d0*(w(1)+w(2)+w(3)+w(4)+w(5)+w(6)+w(7))
w(9) = w(1)
w(10) = w(2)
w(11) = w(3)
w(12) = w(4)
w(13) = w(5)
w(14) = w(6)
w(15) = w(7)
else
write (iout,10)
10 format (/,' Suzuki-Yoshida decomposition of selected order
& not supported. Try RESPA-THERM-SY = 3,5,7, or 15')
call fatal
end if

c
c use multiple time steps and Suzuki-Yoshida decomposition
c
do k = 1, respa_therm_nc
do j = 1, ns
do j = 1, respa_therm_nsy

dts = w(j) * dtc
dt2 = 0.5d0 * dts
Expand Down