From 72cb95e812d15e6f3f9f7a44df250d146dcf0aab Mon Sep 17 00:00:00 2001 From: "Daniel T. Margul" Date: Sun, 28 Jun 2015 02:49:49 -0400 Subject: [PATCH] some bug fixes --- source/bath.f | 11 +-- source/ehal1.f | 180 +++++++++++++++++++++++++++++++++++++++++- source/stoch_respa2.f | 48 +++++++++-- 3 files changed, 226 insertions(+), 13 deletions(-) diff --git a/source/bath.f b/source/bath.f index 9e8a8be67..7acc7062a 100644 --- a/source/bath.f +++ b/source/bath.f @@ -44,6 +44,7 @@ c module bath implicit none + use sizes integer maxnose parameter (maxnose=4) integer len_nhc,isok_L,isok_M @@ -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) diff --git a/source/ehal1.f b/source/ehal1.f index 632244cbc..535b24ab5 100644 --- a/source/ehal1.f +++ b/source/ehal1.f @@ -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) @@ -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 @@ -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 @@ -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)) diff --git a/source/stoch_respa2.f b/source/stoch_respa2.f index a4bca2861..f529efa5e 100644 --- a/source/stoch_respa2.f +++ b/source/stoch_respa2.f @@ -262,7 +262,7 @@ 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) @@ -270,17 +270,53 @@ subroutine applyisok (dt) 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