c c c ################################################### c ## COPYRIGHT (C) 1993 by Jay William Ponder ## c ## All Rights Reserved ## c ################################################### c c ############################################################### c ## ## c ## subroutine ehal -- buffered 14-7 van der Waals energy ## c ## ## c ############################################################### c c c "ehal" calculates the buffered 14-7 van der Waals energy c c modified to optionally compute WCA repulsion between the c solute and solvent molecules c c subroutine ehal implicit none include 'sizes.i' include 'cutoff.i' include 'iounit.i' include 'mutant.i' real*8 esol c c c zero out the solute-solvent Weeks-Chandler-Andersen energy c esol = 0.0d0 c c choose the method for summing over pairwise interactions c if (use_lights) then call ehal0b (esol) else if (use_vlist) then call ehal0c (esol) else call ehal0a (esol) end if c c print the solute-solvent Weeks-Chandler-Anderson energy c write (iout,10) esol,lambda 10 format (' Solute-Solvent WCA Energy :',9x,f12.4,3x,f12.4) return end c c c ################################################################ c ## ## c ## subroutine ehal0a -- buffered 14-7 vdw via double loop ## c ## ## c ################################################################ c c c "ehal0a" calculates the buffered 14-7 van der Waals energy c using a pairwise double loop c c subroutine ehal0a (esol) implicit none include 'sizes.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'cell.i' include 'couple.i' include 'energi.i' include 'group.i' include 'molcul.i' include 'mutant.i' include 'shunt.i' include 'usage.i' include 'vdw.i' include 'vdwpot.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer im,km,ikmax,ikmin integer iv14(maxatm) real*8 e,rv,rv7 real*8 eps,rdn,fgrp real*8 xi,yi,zi real*8 xr,yr,zr real*8 rho,tau,taper real*8 rik,rik2,rik3 real*8 rik4,rik5,rik7 real*8 esol,erep real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) logical proceed,usei c c c zero out the van der Waals energy contribution c ev = 0.0d0 c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c apply any reduction factor to the atomic coordinates c do k = 1, nvdw i = ivdw(k) iv = ired(i) rdn = kred(i) xred(i) = rdn*(x(i)-x(iv)) + x(iv) yred(i) = rdn*(y(i)-y(iv)) + y(iv) zred(i) = rdn*(z(i)-z(iv)) + z(iv) end do c c find the van der Waals energy via double loop search c do ii = 1, nvdw-1 i = ivdw(ii) iv = ired(i) it = class(i) im = molcule(i) xi = xred(i) yi = yred(i) zi = zred(i) usei = (use(i) .or. use(iv)) do j = ii+1, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = ii+1, nvdw k = ivdw(kk) kv = ired(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 = class(k) km = molcule(k) xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) if (use_image) call image (xr,yr,zr,0) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) 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) rv7 = rv**7 rik7 = rik**7 rho = rik7 + ghal*rv7 tau = (dhal+1.0d0) / (rik + dhal*rv) e = eps * rv7 * tau**7 & * ((ghal+1.0d0)*rv7/rho-2.0d0) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 e = e * taper end if c c scale the interaction based on its group membership c if (use_group) e = e * fgrp c c increment the overall van der Waals energy components c ev = ev + e end if end if end do end do c c for periodic boundary conditions with large cutoffs c neighbors must be found by the replicates method c if (.not. use_replica) return c c calculate interaction energy with other unit cells c do ii = 1, nvdw i = ivdw(ii) iv = ired(i) it = class(i) xi = xred(i) yi = yred(i) zi = zred(i) usei = (use(i) .or. use(iv)) do j = ii, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = ii, nvdw k = ivdw(kk) kv = ired(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 = class(k) do j = 2, ncell xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) call image (xr,yr,zr,j) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then rik = sqrt(rik2) rv = radmin(kt,it) eps = epsilon(kt,it) if (use_polymer) then if (rik2 .le. polycut2) then if (iv14(k) .eq. i) then rv = radmin4(kt,it) eps = epsilon4(kt,it) end if eps = eps * vscale(k) end if end if rv7 = rv**7 rik7 = rik**7 rho = rik7 + ghal*rv7 tau = (dhal+1.0d0) / (rik + dhal*rv) e = eps * rv7 * tau**7 & * ((ghal+1.0d0)*rv7/rho-2.0d0) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 e = e * taper end if c c scale the interaction based on its group membership c if (use_group) e = e * fgrp c c increment the overall van der Waals energy component; c interaction of an atom with its own image counts half c if (i .eq. k) e = 0.5d0 * e ev = ev + e end if end do end if end do end do return end c c c ################################################################## c ## ## c ## subroutine ehal0b -- buffered 14-7 vdw energy via lights ## c ## ## c ################################################################## c c c "ehal0b" calculates the buffered 14-7 van der Waals energy c using the method of lights c c subroutine ehal0b (esol) implicit none include 'sizes.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'boxes.i' include 'cell.i' include 'couple.i' include 'energi.i' include 'group.i' include 'light.i' include 'molcul.i' include 'mutant.i' include 'shunt.i' include 'usage.i' include 'vdw.i' include 'vdwpot.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer kgy,kgz integer start,stop integer im,km,ikmax,ikmin integer iv14(maxatm) real*8 e,rv,rv7 real*8 eps,rdn,fgrp real*8 xi,yi,zi real*8 xr,yr,zr real*8 rho,tau,taper real*8 rik,rik2,rik3 real*8 rik4,rik5,rik7 real*8 esol,erep real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) real*8 xsort(maxlight) real*8 ysort(maxlight) real*8 zsort(maxlight) logical proceed,usei logical prime,repeat c c c zero out the van der Waals energy contribution c ev = 0.0d0 c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c apply any reduction factor to the atomic coordinates c do j = 1, nvdw i = ivdw(j) iv = ired(i) rdn = kred(i) xred(j) = rdn*(x(i)-x(iv)) + x(iv) yred(j) = rdn*(y(i)-y(iv)) + y(iv) zred(j) = rdn*(z(i)-z(iv)) + z(iv) end do c c transfer the interaction site coordinates to sorting arrays c do i = 1, nvdw xsort(i) = xred(i) ysort(i) = yred(i) zsort(i) = zred(i) end do c c use the method of lights to generate neighbors c call lights (off,nvdw,xsort,ysort,zsort) c c loop over all atoms computing the interactions c do ii = 1, nvdw i = ivdw(ii) iv = ired(i) it = class(i) xi = xsort(rgx(ii)) yi = ysort(rgy(ii)) zi = zsort(rgz(ii)) usei = (use(i) .or. use(iv)) do j = 1, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do if (kbx(ii) .le. kex(ii)) then repeat = .false. start = kbx(ii) + 1 stop = kex(ii) else repeat = .true. start = 1 stop = kex(ii) end if 10 continue do j = start, stop kk = locx(j) kgy = rgy(kk) if (kby(ii) .le. key(ii)) then if (kgy.lt.kby(ii) .or. kgy.gt.key(ii)) goto 20 else if (kgy.lt.kby(ii) .and. kgy.gt.key(ii)) goto 20 end if kgz = rgz(kk) if (kbz(ii) .le. kez(ii)) then if (kgz.lt.kbz(ii) .or. kgz.gt.kez(ii)) goto 20 else if (kgz.lt.kbz(ii) .and. kgz.gt.kez(ii)) goto 20 end if k = ivdw(kk-((kk-1)/nvdw)*nvdw) kv = ired(k) prime = (kk .le. nvdw) c c decide whether to compute the current interaction c 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 = class(k) xr = xi - xsort(j) yr = yi - ysort(kgy) zr = zi - zsort(kgz) if (use_image) then if (abs(xr) .gt. xcell2) xr = xr - sign(xcell,xr) if (abs(yr) .gt. ycell2) yr = yr - sign(ycell,yr) if (abs(zr) .gt. zcell2) zr = zr - sign(zcell,zr) if (monoclinic) then xr = xr + zr*beta_cos zr = zr * beta_sin else if (triclinic) then xr = xr + yr*gamma_cos + zr*beta_cos yr = yr*gamma_sin + zr*beta_term zr = zr * gamma_term end if end if rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then rik = sqrt(rik2) rv = radmin(kt,it) eps = epsilon(kt,it) if (prime) then if (iv14(k) .eq. i) then rv = radmin4(kt,it) eps = epsilon4(kt,it) end if eps = eps * vscale(k) end if rv7 = rv**7 rik7 = rik**7 rho = rik7 + ghal*rv7 tau = (dhal+1.0d0) / (rik + dhal*rv) e = eps * rv7 * tau**7 & * ((ghal+1.0d0)*rv7/rho-2.0d0) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 e = e * taper end if c c scale the interaction based on its group membership c if (use_group) e = e * fgrp c c increment the overall van der Waals energy components c ev = ev + e end if end if 20 continue end do if (repeat) then repeat = .false. start = kbx(ii) + 1 stop = nlight goto 10 end if end do return end c c c ################################################################ c ## ## c ## subroutine ehal0c -- buffered 14-7 vdw energy via list ## c ## ## c ################################################################ c c c "ehal0c" calculates the buffered 14-7 van der Waals energy c using a pairwise neighbor list c c subroutine ehal0c (esol) implicit none include 'sizes.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'couple.i' include 'energi.i' include 'group.i' include 'molcul.i' include 'mutant.i' include 'neigh.i' include 'shunt.i' include 'usage.i' include 'vdw.i' include 'vdwpot.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer im,km,ikmax,ikmin integer iv14(maxatm) real*8 e,eps,rdn real*8 fgrp,rv,rv7 real*8 xi,yi,zi real*8 xr,yr,zr real*8 rho,tau,taper real*8 rik,rik2,rik3 real*8 rik4,rik5,rik7 real*8 esol,erep real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) logical proceed,usei c c c zero out the van der Waals energy contribution c ev = 0.0d0 c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c apply any reduction factor to the atomic coordinates c do k = 1, nvdw i = ivdw(k) iv = ired(i) rdn = kred(i) xred(i) = rdn*(x(i)-x(iv)) + x(iv) yred(i) = rdn*(y(i)-y(iv)) + y(iv) zred(i) = rdn*(z(i)-z(iv)) + z(iv) end do c c find the van der Waals energy via neighbor list search c do ii = 1, nvdw i = ivdw(ii) iv = ired(i) it = class(i) xi = xred(i) yi = yred(i) zi = zred(i) usei = (use(i) .or. use(iv)) do j = 1, nvlst(ii) vscale(ivdw(vlst(j,ii))) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = 1, nvlst(ii) k = ivdw(vlst(kk,ii)) kv = ired(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 = class(k) xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) if (use_image) call image (xr,yr,zr,0) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) 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) rv7 = rv**7 rik7 = rik**7 rho = rik7 + ghal*rv7 tau = (dhal+1.0d0) / (rik + dhal*rv) e = eps * rv7 * tau**7 & * ((ghal+1.0d0)*rv7/rho-2.0d0) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 e = e * taper end if c c scale the interaction based on its group membership c if (use_group) e = e * fgrp c c increment the overall van der Waals energy components c ev = ev + e end if end if end do end do return end c c c ################################################### c ## COPYRIGHT (C) 1993 by Jay William Ponder ## c ## All Rights Reserved ## c ################################################### c c ################################################################ c ## ## c ## subroutine ehal1 -- buffered 14-7 energy & derivatives ## c ## ## c ################################################################ c c c "ehal1" calculates the buffered 14-7 van der Waals energy and c its first derivatives with respect to Cartesian coordinates c c modified to optionally compute WCA repulsion between the c solute and solvent molecules c c subroutine ehal1 implicit none include 'sizes.i' include 'cutoff.i' include 'iounit.i' include 'mutant.i' real*8 esol c c c zero out the solute-solvent Weeks-Chandler-Andersen energy c esol = 0.0d0 c c choose the method for summing over pairwise interactions c if (use_lights) then call ehal1b (esol) else if (use_vlist) then call ehal1c (esol) else call ehal1a (esol) end if c c print the solute-solvent Weeks-Chandler-Anderson energy c write (iout,10) esol,lambda 10 format (' Solute-Solvent WCA Energy :',9x,f12.4,3x,f12.4) return end c c c ################################################################# c ## ## c ## subroutine ehal1a -- double loop buffer 14-7 vdw derivs ## c ## ## c ################################################################# c c c "ehal1a" calculates the buffered 14-7 van der Waals energy and c its first derivatives with respect to Cartesian coordinates c using a pairwise double loop c c subroutine ehal1a (esol) implicit none include 'sizes.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'cell.i' include 'couple.i' include 'deriv.i' include 'energi.i' include 'group.i' include 'inter.i' include 'molcul.i' include 'mutant.i' include 'shunt.i' include 'usage.i' include 'vdw.i' include 'vdwpot.i' include 'virial.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer im,km,ikmax,ikmin integer iv14(maxatm) real*8 e,de,eps,rdn real*8 fgrp,rv,rv7 real*8 xi,yi,zi real*8 xr,yr,zr real*8 redi,rediv real*8 redk,redkv real*8 dedx,dedy,dedz real*8 rho,tau,tau7 real*8 dtau,gtau real*8 taper,dtaper real*8 rik,rik2,rik3 real*8 rik4,rik5 real*8 rik6,rik7 real*8 vxx,vyy,vzz real*8 vyx,vzx,vzy real*8 esol,erep real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) logical proceed,usei c c c zero out the van der Waals energy and first derivatives c ev = 0.0d0 do i = 1, n dev(1,i) = 0.0d0 dev(2,i) = 0.0d0 dev(3,i) = 0.0d0 end do c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c apply any reduction factor to the atomic coordinates c do k = 1, nvdw i = ivdw(k) iv = ired(i) rdn = kred(i) xred(i) = rdn*(x(i)-x(iv)) + x(iv) yred(i) = rdn*(y(i)-y(iv)) + y(iv) zred(i) = rdn*(z(i)-z(iv)) + z(iv) end do c c find van der Waals energy and derivatives via double loop c do ii = 1, nvdw-1 i = ivdw(ii) iv = ired(i) redi = kred(i) rediv = 1.0d0 - redi it = class(i) im = molcule(i) xi = xred(i) yi = yred(i) zi = zred(i) usei = (use(i) .or. use(iv)) do j = ii+1, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = ii+1, nvdw k = ivdw(kk) kv = ired(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 = class(k) km = molcule(k) xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) if (use_image) call image (xr,yr,zr,0) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then 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) rik = sqrt(rik2) 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 = eps*tau7*rv7*((ghal+1.0d0)*rv7/rho-2.0d0) de = -7.0d0 * (dtau*e+gtau) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 de = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 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 ev = ev + e if (i .eq. iv) then dev(1,i) = dev(1,i) + dedx dev(2,i) = dev(2,i) + dedy dev(3,i) = dev(3,i) + dedz else dev(1,i) = dev(1,i) + dedx*redi dev(2,i) = dev(2,i) + dedy*redi dev(3,i) = dev(3,i) + dedz*redi dev(1,iv) = dev(1,iv) + dedx*rediv dev(2,iv) = dev(2,iv) + dedy*rediv dev(3,iv) = dev(3,iv) + dedz*rediv end if if (k .eq. kv) then dev(1,k) = dev(1,k) - dedx dev(2,k) = dev(2,k) - dedy dev(3,k) = dev(3,k) - dedz else redk = kred(k) redkv = 1.0d0 - redk dev(1,k) = dev(1,k) - dedx*redk dev(2,k) = dev(2,k) - dedy*redk dev(3,k) = dev(3,k) - dedz*redk dev(1,kv) = dev(1,kv) - dedx*redkv dev(2,kv) = dev(2,kv) - dedy*redkv dev(3,kv) = dev(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 vir(1,1) = vir(1,1) + vxx vir(2,1) = vir(2,1) + vyx vir(3,1) = vir(3,1) + vzx vir(1,2) = vir(1,2) + vyx vir(2,2) = vir(2,2) + vyy vir(3,2) = vir(3,2) + vzy vir(1,3) = vir(1,3) + vzx vir(2,3) = vir(2,3) + vzy vir(3,3) = vir(3,3) + vzz c c increment the total intermolecular energy c if (molcule(i) .ne. molcule(k)) then einter = einter + e end if end if end if end do end do c c for periodic boundary conditions with large cutoffs c neighbors must be found by the replicates method c if (.not. use_replica) return c c calculate interaction energy with other unit cells c do ii = 1, nvdw i = ivdw(ii) iv = ired(i) redi = kred(i) rediv = 1.0d0 - redi it = class(i) im = molcule(i) xi = xred(i) yi = yred(i) zi = zred(i) usei = (use(i) .or. use(iv)) do j = ii, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = ii, nvdw k = ivdw(kk) kv = ired(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 = class(k) km = molcule(k) do j = 2, ncell xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) call image (xr,yr,zr,j) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then rv = radmin(kt,it) eps = epsilon(kt,it) if (use_polymer) then if (rik2 .le. polycut2) then if (iv14(k) .eq. i) then rv = radmin4(kt,it) eps = epsilon4(kt,it) end if eps = eps * vscale(k) end if end if rik = sqrt(rik2) 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 = eps*tau7*rv7*((ghal+1.0d0)*rv7/rho-2.0d0) de = -7.0d0 * (dtau*e+gtau) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 de = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 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 if (i .eq. k) e = 0.5d0 * e ev = ev + e if (i .eq. iv) then dev(1,i) = dev(1,i) + dedx dev(2,i) = dev(2,i) + dedy dev(3,i) = dev(3,i) + dedz else dev(1,i) = dev(1,i) + dedx*redi dev(2,i) = dev(2,i) + dedy*redi dev(3,i) = dev(3,i) + dedz*redi dev(1,iv) = dev(1,iv) + dedx*rediv dev(2,iv) = dev(2,iv) + dedy*rediv dev(3,iv) = dev(3,iv) + dedz*rediv end if if (i .ne. k) then if (k .eq. kv) then dev(1,k) = dev(1,k) - dedx dev(2,k) = dev(2,k) - dedy dev(3,k) = dev(3,k) - dedz else redk = kred(k) redkv = 1.0d0 - redk dev(1,k) = dev(1,k) - dedx*redk dev(2,k) = dev(2,k) - dedy*redk dev(3,k) = dev(3,k) - dedz*redk dev(1,kv) = dev(1,kv) - dedx*redkv dev(2,kv) = dev(2,kv) - dedy*redkv dev(3,kv) = dev(3,kv) - dedz*redkv end if 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 vir(1,1) = vir(1,1) + vxx vir(2,1) = vir(2,1) + vyx vir(3,1) = vir(3,1) + vzx vir(1,2) = vir(1,2) + vyx vir(2,2) = vir(2,2) + vyy vir(3,2) = vir(3,2) + vzy vir(1,3) = vir(1,3) + vzx vir(2,3) = vir(2,3) + vzy vir(3,3) = vir(3,3) + vzz c c increment the total intermolecular energy c einter = einter + e end if end do end if end do end do return end c c c ################################################################## c ## ## c ## subroutine ehal1b -- buffered 14-7 vdw derivs via lights ## c ## ## c ################################################################## c c c "ehal1b" calculates the buffered 14-7 van der Waals energy and c its first derivatives with respect to Cartesian coordinates c using the method of lights c c subroutine ehal1b (esol) implicit none include 'sizes.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'boxes.i' include 'cell.i' include 'couple.i' include 'deriv.i' include 'energi.i' include 'group.i' include 'inter.i' include 'iounit.i' include 'light.i' include 'molcul.i' include 'mutant.i' include 'shunt.i' include 'usage.i' include 'vdw.i' include 'vdwpot.i' include 'virial.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer kgy,kgz integer start,stop integer im,km,ikmax,ikmin integer iv14(maxatm) real*8 e,de,eps,rdn real*8 fgrp,rv,rv7 real*8 xi,yi,zi real*8 xr,yr,zr real*8 redi,rediv real*8 redk,redkv real*8 dedx,dedy,dedz real*8 rho,tau,tau7 real*8 dtau,gtau real*8 taper,dtaper real*8 rik,rik2,rik3 real*8 rik4,rik5 real*8 rik6,rik7 real*8 vxx,vyy,vzz real*8 vyx,vzx,vzy real*8 esol,erep real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) real*8 xsort(maxlight) real*8 ysort(maxlight) real*8 zsort(maxlight) logical proceed,usei logical prime,repeat c c c zero out the van der Waals energy and first derivatives c ev = 0.0d0 do i = 1, n dev(1,i) = 0.0d0 dev(2,i) = 0.0d0 dev(3,i) = 0.0d0 end do c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c apply any reduction factor to the atomic coordinates c do j = 1, nvdw i = ivdw(j) iv = ired(i) rdn = kred(i) xred(j) = rdn*(x(i)-x(iv)) + x(iv) yred(j) = rdn*(y(i)-y(iv)) + y(iv) zred(j) = rdn*(z(i)-z(iv)) + z(iv) end do c c transfer the interaction site coordinates to sorting arrays c do i = 1, nvdw xsort(i) = xred(i) ysort(i) = yred(i) zsort(i) = zred(i) end do c c use the method of lights to generate neighbors c call lights (off,nvdw,xsort,ysort,zsort) c c loop over all atoms computing the interactions c do ii = 1, nvdw i = ivdw(ii) iv = ired(i) redi = kred(i) rediv = 1.0d0 - redi it = class(i) im = molcule(i) xi = xsort(rgx(ii)) yi = ysort(rgy(ii)) zi = zsort(rgz(ii)) usei = (use(i) .or. use(iv)) do j = 1, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do if (kbx(ii) .le. kex(ii)) then repeat = .false. start = kbx(ii) + 1 stop = kex(ii) else repeat = .true. start = 1 stop = kex(ii) end if 10 continue do j = start, stop kk = locx(j) kgy = rgy(kk) if (kby(ii) .le. key(ii)) then if (kgy.lt.kby(ii) .or. kgy.gt.key(ii)) goto 20 else if (kgy.lt.kby(ii) .and. kgy.gt.key(ii)) goto 20 end if kgz = rgz(kk) if (kbz(ii) .le. kez(ii)) then if (kgz.lt.kbz(ii) .or. kgz.gt.kez(ii)) goto 20 else if (kgz.lt.kbz(ii) .and. kgz.gt.kez(ii)) goto 20 end if k = ivdw(kk-((kk-1)/nvdw)*nvdw) kv = ired(k) prime = (kk .le. nvdw) c c decide whether to compute the current interaction c 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 = class(k) km = molcule(k) xr = xi - xsort(j) yr = yi - ysort(kgy) zr = zi - zsort(kgz) if (use_image) then if (abs(xr) .gt. xcell2) xr = xr - sign(xcell,xr) if (abs(yr) .gt. ycell2) yr = yr - sign(ycell,yr) if (abs(zr) .gt. zcell2) zr = zr - sign(zcell,zr) if (monoclinic) then xr = xr + zr*beta_cos zr = zr * beta_sin else if (triclinic) then xr = xr + yr*gamma_cos + zr*beta_cos yr = yr*gamma_sin + zr*beta_term zr = zr * gamma_term end if end if rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then rv = radmin(kt,it) eps = epsilon(kt,it) if (prime) then if (iv14(k) .eq. i) then rv = radmin4(kt,it) eps = epsilon4(kt,it) end if eps = eps * vscale(k) end if rik = sqrt(rik2) 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 = eps*tau7*rv7*((ghal+1.0d0)*rv7/rho-2.0d0) de = -7.0d0 * (dtau*e+gtau) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 de = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 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 ev = ev + e if (i .eq. iv) then dev(1,i) = dev(1,i) + dedx dev(2,i) = dev(2,i) + dedy dev(3,i) = dev(3,i) + dedz else dev(1,i) = dev(1,i) + dedx*redi dev(2,i) = dev(2,i) + dedy*redi dev(3,i) = dev(3,i) + dedz*redi dev(1,iv) = dev(1,iv) + dedx*rediv dev(2,iv) = dev(2,iv) + dedy*rediv dev(3,iv) = dev(3,iv) + dedz*rediv end if if (k .eq. kv) then dev(1,k) = dev(1,k) - dedx dev(2,k) = dev(2,k) - dedy dev(3,k) = dev(3,k) - dedz else redk = kred(k) redkv = 1.0d0 - redk dev(1,k) = dev(1,k) - dedx*redk dev(2,k) = dev(2,k) - dedy*redk dev(3,k) = dev(3,k) - dedz*redk dev(1,kv) = dev(1,kv) - dedx*redkv dev(2,kv) = dev(2,kv) - dedy*redkv dev(3,kv) = dev(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 vir(1,1) = vir(1,1) + vxx vir(2,1) = vir(2,1) + vyx vir(3,1) = vir(3,1) + vzx vir(1,2) = vir(1,2) + vyx vir(2,2) = vir(2,2) + vyy vir(3,2) = vir(3,2) + vzy vir(1,3) = vir(1,3) + vzx vir(2,3) = vir(2,3) + vzy vir(3,3) = vir(3,3) + vzz c c increment the total intermolecular energy c if (.not.prime .or. molcule(i).ne.molcule(k)) then einter = einter + e end if end if end if 20 continue end do if (repeat) then repeat = .false. start = kbx(ii) + 1 stop = nlight goto 10 end if end do return end c c c ################################################################ c ## ## c ## subroutine ehal1c -- buffered 14-7 vdw derivs via list ## c ## ## c ################################################################ c c c "ehal1c" calculates the buffered 14-7 van der Waals energy and c its first derivatives with respect to Cartesian coordinates c using a pairwise neighbor list c c subroutine ehal1c (esol) implicit none include 'sizes.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'couple.i' include 'deriv.i' include 'energi.i' include 'group.i' include 'inter.i' include 'molcul.i' include 'mutant.i' include 'neigh.i' include 'shunt.i' include 'usage.i' include 'vdw.i' include 'vdwpot.i' include 'virial.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer im,km,ikmax,ikmin integer iv14(maxatm) real*8 e,de,eps,rdn real*8 fgrp,rv,rv7 real*8 xi,yi,zi real*8 xr,yr,zr real*8 redi,rediv real*8 redk,redkv real*8 dedx,dedy,dedz real*8 rho,tau,tau7 real*8 dtau,gtau real*8 taper,dtaper real*8 rik,rik2,rik3 real*8 rik4,rik5 real*8 rik6,rik7 real*8 vxx,vyy,vzz real*8 vyx,vzx,vzy real*8 esol,erep real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) logical proceed,usei c c c zero out the van der Waals energy and first derivatives c ev = 0.0d0 do i = 1, n dev(1,i) = 0.0d0 dev(2,i) = 0.0d0 dev(3,i) = 0.0d0 end do c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c apply any reduction factor to the atomic coordinates c do k = 1, nvdw i = ivdw(k) iv = ired(i) rdn = kred(i) xred(i) = rdn*(x(i)-x(iv)) + x(iv) yred(i) = rdn*(y(i)-y(iv)) + y(iv) zred(i) = rdn*(z(i)-z(iv)) + z(iv) end do c c find van der Waals energy and derivatives via neighbor list c do ii = 1, nvdw i = ivdw(ii) iv = ired(i) redi = kred(i) rediv = 1.0d0 - redi it = class(i) im = molcule(i) xi = xred(i) yi = yred(i) zi = zred(i) usei = (use(i) .or. use(iv)) do j = 1, nvlst(ii) vscale(ivdw(vlst(j,ii))) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = 1, nvlst(ii) k = ivdw(vlst(kk,ii)) kv = ired(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 = class(k) km = molcule(k) xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) if (use_image) call image (xr,yr,zr,0) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then 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) rik = sqrt(rik2) 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 = eps*tau7*rv7*((ghal+1.0d0)*rv7/rho-2.0d0) de = -7.0d0 * (dtau*e+gtau) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 de = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 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 ev = ev + e if (i .eq. iv) then dev(1,i) = dev(1,i) + dedx dev(2,i) = dev(2,i) + dedy dev(3,i) = dev(3,i) + dedz else dev(1,i) = dev(1,i) + dedx*redi dev(2,i) = dev(2,i) + dedy*redi dev(3,i) = dev(3,i) + dedz*redi dev(1,iv) = dev(1,iv) + dedx*rediv dev(2,iv) = dev(2,iv) + dedy*rediv dev(3,iv) = dev(3,iv) + dedz*rediv end if if (k .eq. kv) then dev(1,k) = dev(1,k) - dedx dev(2,k) = dev(2,k) - dedy dev(3,k) = dev(3,k) - dedz else redk = kred(k) redkv = 1.0d0 - redk dev(1,k) = dev(1,k) - dedx*redk dev(2,k) = dev(2,k) - dedy*redk dev(3,k) = dev(3,k) - dedz*redk dev(1,kv) = dev(1,kv) - dedx*redkv dev(2,kv) = dev(2,kv) - dedy*redkv dev(3,kv) = dev(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 vir(1,1) = vir(1,1) + vxx vir(2,1) = vir(2,1) + vyx vir(3,1) = vir(3,1) + vzx vir(1,2) = vir(1,2) + vyx vir(2,2) = vir(2,2) + vyy vir(3,2) = vir(3,2) + vzy vir(1,3) = vir(1,3) + vzx vir(2,3) = vir(2,3) + vzy vir(3,3) = vir(3,3) + vzz c c increment the total intermolecular energy c if (molcule(i) .ne. molcule(k)) then einter = einter + e end if end if end if end do end do return end c c c ################################################### c ## COPYRIGHT (C) 1993 by Jay William Ponder ## c ## All Rights Reserved ## c ################################################### c c ################################################################ c ## ## c ## subroutine ehal2 -- atom-by-atom buffered 14-7 Hessian ## c ## ## c ################################################################ c c c "ehal2" calculates the buffered 14-7 van der Waals second c derivatives for a single atom at a time c c subroutine ehal2 (iatom,xred,yred,zred) implicit none include 'sizes.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'cell.i' include 'couple.i' include 'group.i' include 'hessn.i' include 'molcul.i' include 'mutant.i' include 'shunt.i' include 'vdw.i' include 'vdwpot.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer iatom,jcell integer im,km,ikmax,ikmin integer nlist,list(5) integer iv14(maxatm) real*8 e,de,d2e,fgrp real*8 eps,rv,rv7 real*8 xi,yi,zi real*8 xr,yr,zr real*8 redi,rediv real*8 redk,redkv real*8 redi2,rediv2,rediiv real*8 redik,redivk real*8 redikv,redivkv real*8 rho,tau,tau7 real*8 dtau,gtau real*8 taper,dtaper,d2taper real*8 rik,rik2,rik3 real*8 rik4,rik5 real*8 rik6,rik7 real*8 d2edx,d2edy,d2edz real*8 esol,erep real*8 term(3,3) real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) logical proceed c c c zero out the solute-solvent Weeks-Chandler-Andersen energy c esol = 0.0d0 c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c check to see if the atom of interest is a vdw site c nlist = 0 do k = 1, nvdw if (ivdw(k) .eq. iatom) then nlist = nlist + 1 list(nlist) = iatom goto 10 end if end do return 10 continue c c determine the atoms involved via reduction factors c nlist = 1 list(nlist) = iatom do k = 1, n12(iatom) i = i12(k,iatom) if (ired(i) .eq. iatom) then nlist = nlist + 1 list(nlist) = i end if end do c c find van der Waals Hessian elements for involved atoms c do ii = 1, nlist i = list(ii) iv = ired(i) redi = kred(i) if (i .ne. iv) then rediv = 1.0d0 - redi redi2 = redi * redi rediv2 = rediv * rediv rediiv = redi * rediv end if it = class(i) im = molcule(i) xi = xred(i) yi = yred(i) zi = zred(i) do j = 1, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = 1, nvdw k = ivdw(kk) kv = ired(k) proceed = .true. if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) if (proceed) proceed = (k .ne. i) c c compute the Hessian elements for this interaction c if (proceed) then kt = class(k) km = molcule(k) xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) if (use_image) call image (xr,yr,zr,0) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then 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) rik = sqrt(rik2) 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 = eps*tau7*rv7*((ghal+1.0d0)*rv7/rho-2.0d0) de = -7.0d0 * (dtau*e+gtau) d2e = 56.0d0*dtau*dtau*e - 42.0d0*gtau/rik & + 98.0d0*gtau*(dtau+rik6/rho) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 de = 0.0d0 d2e = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 = rik3 * rik2 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 d2taper = 20.0d0*c5*rik3 + 12.0d0*c4*rik2 & + 6.0d0*c3*rik + 2.0d0*c2 d2e = e*d2taper + 2.0d0*de*dtaper + d2e*taper de = e*dtaper + de*taper end if c c scale the interaction based on its group membership c if (use_group) then de = de * fgrp d2e = d2e * fgrp end if c c get chain rule terms for van der Waals Hessian elements c de = de / rik d2e = (d2e-de) / rik2 d2edx = d2e * xr d2edy = d2e * yr d2edz = d2e * zr term(1,1) = d2edx*xr + de term(1,2) = d2edx*yr term(1,3) = d2edx*zr term(2,1) = term(1,2) term(2,2) = d2edy*yr + de term(2,3) = d2edy*zr term(3,1) = term(1,3) term(3,2) = term(2,3) term(3,3) = d2edz*zr + de c c increment diagonal and non-diagonal Hessian elements c if (i .eq. iatom) then if (i.eq.iv .and. k.eq.kv) then do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j) hessy(j,i) = hessy(j,i) + term(2,j) hessz(j,i) = hessz(j,i) + term(3,j) hessx(j,k) = hessx(j,k) - term(1,j) hessy(j,k) = hessy(j,k) - term(2,j) hessz(j,k) = hessz(j,k) - term(3,j) end do else if (k .eq. kv) then do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j)*redi2 hessy(j,i) = hessy(j,i) + term(2,j)*redi2 hessz(j,i) = hessz(j,i) + term(3,j)*redi2 hessx(j,k) = hessx(j,k) - term(1,j)*redi hessy(j,k) = hessy(j,k) - term(2,j)*redi hessz(j,k) = hessz(j,k) - term(3,j)*redi hessx(j,iv) = hessx(j,iv) + term(1,j)*rediiv hessy(j,iv) = hessy(j,iv) + term(2,j)*rediiv hessz(j,iv) = hessz(j,iv) + term(3,j)*rediiv end do else if (i .eq. iv) then redk = kred(k) redkv = 1.0d0 - redk do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j) hessy(j,i) = hessy(j,i) + term(2,j) hessz(j,i) = hessz(j,i) + term(3,j) hessx(j,k) = hessx(j,k) - term(1,j)*redk hessy(j,k) = hessy(j,k) - term(2,j)*redk hessz(j,k) = hessz(j,k) - term(3,j)*redk hessx(j,kv) = hessx(j,kv) - term(1,j)*redkv hessy(j,kv) = hessy(j,kv) - term(2,j)*redkv hessz(j,kv) = hessz(j,kv) - term(3,j)*redkv end do else redk = kred(k) redkv = 1.0d0 - redk redik = redi * redk redikv = redi * redkv do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j)*redi2 hessy(j,i) = hessy(j,i) + term(2,j)*redi2 hessz(j,i) = hessz(j,i) + term(3,j)*redi2 hessx(j,k) = hessx(j,k) - term(1,j)*redik hessy(j,k) = hessy(j,k) - term(2,j)*redik hessz(j,k) = hessz(j,k) - term(3,j)*redik hessx(j,iv) = hessx(j,iv) + term(1,j)*rediiv hessy(j,iv) = hessy(j,iv) + term(2,j)*rediiv hessz(j,iv) = hessz(j,iv) + term(3,j)*rediiv hessx(j,kv) = hessx(j,kv) - term(1,j)*redikv hessy(j,kv) = hessy(j,kv) - term(2,j)*redikv hessz(j,kv) = hessz(j,kv) - term(3,j)*redikv end do end if else if (iv .eq. iatom) then if (k .eq. kv) then do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j)*rediiv hessy(j,i) = hessy(j,i) + term(2,j)*rediiv hessz(j,i) = hessz(j,i) + term(3,j)*rediiv hessx(j,k) = hessx(j,k) - term(1,j)*rediv hessy(j,k) = hessy(j,k) - term(2,j)*rediv hessz(j,k) = hessz(j,k) - term(3,j)*rediv hessx(j,iv) = hessx(j,iv) + term(1,j)*rediv2 hessy(j,iv) = hessy(j,iv) + term(2,j)*rediv2 hessz(j,iv) = hessz(j,iv) + term(3,j)*rediv2 end do else redk = kred(k) redkv = 1.0d0 - redk redivk = rediv * redk redivkv = rediv * redkv do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j)*rediiv hessy(j,i) = hessy(j,i) + term(2,j)*rediiv hessz(j,i) = hessz(j,i) + term(3,j)*rediiv hessx(j,k) = hessx(j,k) - term(1,j)*redivk hessy(j,k) = hessy(j,k) - term(2,j)*redivk hessz(j,k) = hessz(j,k) - term(3,j)*redivk hessx(j,iv) = hessx(j,iv) + term(1,j)*rediv2 hessy(j,iv) = hessy(j,iv) + term(2,j)*rediv2 hessz(j,iv) = hessz(j,iv) + term(3,j)*rediv2 hessx(j,kv) = hessx(j,kv) - term(1,j)*redivkv hessy(j,kv) = hessy(j,kv) - term(2,j)*redivkv hessz(j,kv) = hessz(j,kv) - term(3,j)*redivkv end do end if end if end if end if end do end do c c for periodic boundary conditions with large cutoffs c neighbors must be found by the replicates method c if (.not. use_replica) return c c calculate interaction energy with other unit cells c do ii = 1, nlist i = list(ii) iv = ired(i) redi = kred(i) if (i .ne. iv) then rediv = 1.0d0 - redi redi2 = redi * redi rediv2 = rediv * rediv rediiv = redi * rediv end if it = class(i) im = molcule(i) xi = xred(i) yi = yred(i) zi = zred(i) do j = 1, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = 1, nvdw k = ivdw(kk) kv = ired(k) proceed = .true. if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) c c compute the Hessian elements for this interaction c if (proceed) then kt = class(k) km = molcule(k) do jcell = 2, ncell xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) call image (xr,yr,zr,jcell) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then rv = radmin(kt,it) eps = epsilon(kt,it) if (use_polymer) then if (rik2 .le. polycut2) then if (iv14(k) .eq. i) then rv = radmin4(kt,it) eps = epsilon4(kt,it) end if eps = eps * vscale(k) end if end if rik = sqrt(rik2) 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 = eps*rv7*tau7*((ghal+1.0d0)*rv7/rho-2.0d0) de = -7.0d0 * (dtau*e+gtau) d2e = 56.0d0*dtau*dtau*e - 42.0d0*gtau/rik & + 98.0d0*gtau*(dtau+rik6/rho) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 de = 0.0d0 d2e = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 = rik3 * rik2 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 d2taper = 20.0d0*c5*rik3 + 12.0d0*c4*rik2 & + 6.0d0*c3*rik + 2.0d0*c2 d2e = e*d2taper + 2.0d0*de*dtaper + d2e*taper de = e*dtaper + de*taper end if c c scale the interaction based on its group membership c if (use_group) then de = de * fgrp d2e = d2e * fgrp end if c c get chain rule terms for van der Waals Hessian elements c de = de / rik d2e = (d2e-de) / rik2 d2edx = d2e * xr d2edy = d2e * yr d2edz = d2e * zr term(1,1) = d2edx*xr + de term(1,2) = d2edx*yr term(1,3) = d2edx*zr term(2,1) = term(1,2) term(2,2) = d2edy*yr + de term(2,3) = d2edy*zr term(3,1) = term(1,3) term(3,2) = term(2,3) term(3,3) = d2edz*zr + de c c increment diagonal and non-diagonal Hessian elements c if (i .eq. iatom) then if (i.eq.iv .and. k.eq.kv) then do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j) hessy(j,i) = hessy(j,i) + term(2,j) hessz(j,i) = hessz(j,i) + term(3,j) hessx(j,k) = hessx(j,k) - term(1,j) hessy(j,k) = hessy(j,k) - term(2,j) hessz(j,k) = hessz(j,k) - term(3,j) end do else if (k .eq. kv) then do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j)*redi2 hessy(j,i) = hessy(j,i) + term(2,j)*redi2 hessz(j,i) = hessz(j,i) + term(3,j)*redi2 hessx(j,k) = hessx(j,k) - term(1,j)*redi hessy(j,k) = hessy(j,k) - term(2,j)*redi hessz(j,k) = hessz(j,k) - term(3,j)*redi hessx(j,iv) = hessx(j,iv) & + term(1,j)*rediiv hessy(j,iv) = hessy(j,iv) & + term(2,j)*rediiv hessz(j,iv) = hessz(j,iv) & + term(3,j)*rediiv end do else if (i .eq. iv) then redk = kred(k) redkv = 1.0d0 - redk do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j) hessy(j,i) = hessy(j,i) + term(2,j) hessz(j,i) = hessz(j,i) + term(3,j) hessx(j,k) = hessx(j,k) - term(1,j)*redk hessy(j,k) = hessy(j,k) - term(2,j)*redk hessz(j,k) = hessz(j,k) - term(3,j)*redk hessx(j,kv) = hessx(j,kv) & - term(1,j)*redkv hessy(j,kv) = hessy(j,kv) & - term(2,j)*redkv hessz(j,kv) = hessz(j,kv) & - term(3,j)*redkv end do else redk = kred(k) redkv = 1.0d0 - redk redik = redi * redk redikv = redi * redkv do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j)*redi2 hessy(j,i) = hessy(j,i) + term(2,j)*redi2 hessz(j,i) = hessz(j,i) + term(3,j)*redi2 hessx(j,k) = hessx(j,k) - term(1,j)*redik hessy(j,k) = hessy(j,k) - term(2,j)*redik hessz(j,k) = hessz(j,k) - term(3,j)*redik hessx(j,iv) = hessx(j,iv) & + term(1,j)*rediiv hessy(j,iv) = hessy(j,iv) & + term(2,j)*rediiv hessz(j,iv) = hessz(j,iv) & + term(3,j)*rediiv hessx(j,kv) = hessx(j,kv) & - term(1,j)*redikv hessy(j,kv) = hessy(j,kv) & - term(2,j)*redikv hessz(j,kv) = hessz(j,kv) & - term(3,j)*redikv end do end if else if (iv .eq. iatom) then if (k .eq. kv) then do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j)*rediiv hessy(j,i) = hessy(j,i) + term(2,j)*rediiv hessz(j,i) = hessz(j,i) + term(3,j)*rediiv hessx(j,k) = hessx(j,k) - term(1,j)*rediv hessy(j,k) = hessy(j,k) - term(2,j)*rediv hessz(j,k) = hessz(j,k) - term(3,j)*rediv hessx(j,iv) = hessx(j,iv) & + term(1,j)*rediv2 hessy(j,iv) = hessy(j,iv) & + term(2,j)*rediv2 hessz(j,iv) = hessz(j,iv) & + term(3,j)*rediv2 end do else redk = kred(k) redkv = 1.0d0 - redk redivk = rediv * redk redivkv = rediv * redkv do j = 1, 3 hessx(j,i) = hessx(j,i) + term(1,j)*rediiv hessy(j,i) = hessy(j,i) + term(2,j)*rediiv hessz(j,i) = hessz(j,i) + term(3,j)*rediiv hessx(j,k) = hessx(j,k) - term(1,j)*redivk hessy(j,k) = hessy(j,k) - term(2,j)*redivk hessz(j,k) = hessz(j,k) - term(3,j)*redivk hessx(j,iv) = hessx(j,iv) & + term(1,j)*rediv2 hessy(j,iv) = hessy(j,iv) & + term(2,j)*rediv2 hessz(j,iv) = hessz(j,iv) & + term(3,j)*rediv2 hessx(j,kv) = hessx(j,kv) & - term(1,j)*redivkv hessy(j,kv) = hessy(j,kv) & - term(2,j)*redivkv hessz(j,kv) = hessz(j,kv) & - term(3,j)*redivkv end do end if end if end if end do end if end do end do return end c c c ################################################### c ## COPYRIGHT (C) 1993 by Jay William Ponder ## c ## All Rights Reserved ## c ################################################### c c ################################################################# c ## ## c ## subroutine ehal3 -- buffered 14-7 vdw energy & analysis ## c ## ## c ################################################################# c c c "ehal3" calculates the buffered 14-7 van der Waals energy c and partitions the energy among the atoms c c modified to optionally compute WCA repulsion between the c solute and solvent molecules c c subroutine ehal3 implicit none include 'sizes.i' include 'cutoff.i' include 'iounit.i' include 'mutant.i' real*8 esol c c c zero out the solute-solvent Weeks-Chandler-Andersen energy c esol = 0.0d0 c c choose the method for summing over pairwise interactions c if (use_lights) then call ehal3b (esol) else if (use_vlist) then call ehal3c (esol) else call ehal3a (esol) end if c c print the solute-solvent Weeks-Chandler-Anderson energy c write (iout,10) esol,lambda 10 format (' Solute-Solvent WCA Energy :',9x,f12.4,3x,f12.4) return end c c c ################################################################# c ## ## c ## subroutine ehal3a -- double loop buffered 14-7 analysis ## c ## ## c ################################################################# c c c "ehal3a" calculates the buffered 14-7 van der Waals energy c and partitions the energy among the atoms using a pairwise c double loop c c subroutine ehal3a (esol) implicit none include 'sizes.i' include 'action.i' include 'analyz.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'cell.i' include 'couple.i' include 'energi.i' include 'group.i' include 'inform.i' include 'inter.i' include 'iounit.i' include 'molcul.i' include 'mutant.i' include 'shunt.i' include 'usage.i' include 'vdw.i' include 'vdwpot.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer im,km,ikmax,ikmin integer iv14(maxatm) real*8 e,rv,rv7 real*8 eps,rdn,fgrp real*8 xi,yi,zi real*8 xr,yr,zr real*8 rho,tau,taper real*8 rik,rik2,rik3 real*8 rik4,rik5,rik7 real*8 esol,erep real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) logical proceed,usei logical header,huge c c c zero out the van der Waals energy and partitioning terms c nev = 0 ev = 0.0d0 do i = 1, n aev(i) = 0.0d0 end do header = .true. c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c apply any reduction factor to the atomic coordinates c do k = 1, nvdw i = ivdw(k) iv = ired(i) rdn = kred(i) xred(i) = rdn*(x(i)-x(iv)) + x(iv) yred(i) = rdn*(y(i)-y(iv)) + y(iv) zred(i) = rdn*(z(i)-z(iv)) + z(iv) end do c c find the van der Waals energy via double loop search c do ii = 1, nvdw-1 i = ivdw(ii) iv = ired(i) it = class(i) im = molcule(i) xi = xred(i) yi = yred(i) zi = zred(i) usei = (use(i) .or. use(iv)) do j = ii+1, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = ii+1, nvdw k = ivdw(kk) kv = ired(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)) if (proceed) proceed = (vscale(k) .ne. 0.0d0) c c compute the energy contribution for this interaction c if (proceed) then kt = class(k) km = molcule(k) xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) if (use_image) call image (xr,yr,zr,0) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) 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) rv7 = rv**7 rik7 = rik**7 rho = rik7 + ghal*rv7 tau = (dhal+1.0d0) / (rik + dhal*rv) e = eps * rv7 * tau**7 & * ((ghal+1.0d0)*rv7/rho-2.0d0) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e end if c c use energy switching if near the cutoff distance c if (rik2 .gt. cut2) then rik = sqrt(rik2) rik3 = rik2 * rik rik4 = rik2 * rik2 rik5 = rik2 * rik3 taper = c5*rik5 + c4*rik4 + c3*rik3 & + c2*rik2 + c1*rik + c0 e = e * taper end if c c scale the interaction based on its group membership c if (use_group) e = e * fgrp c c increment the overall van der Waals energy components c if (e .ne. 0.0d0) then nev = nev + 1 ev = ev + e aev(i) = aev(i) + 0.5d0*e aev(k) = aev(k) + 0.5d0*e end if c c increment the total intermolecular energy c if (molcule(i) .ne. molcule(k)) then einter = einter + e end if c c print a message if the energy of this interaction is large c huge = (e .gt. 10.0d0) if ((debug.and.e.ne.0.0d0) & .or. (verbose.and.huge)) then if (header) then header = .false. write (iout,10) 10 format (/,' Individual van der Waals', & ' Interactions :', & //,' Type',13x,'Atom Names', & 18x,'Minimum',4x,'Actual', & 6x,'Energy',/) end if write (iout,20) i,name(i),k,name(k), & rv,sqrt(rik2),e 20 format (' VDW-Hal',4x,i5,'-',a3,1x,i5,'-', & a3,12x,2f10.4,f12.4) end if end if end if end do end do c c for periodic boundary conditions with large cutoffs c neighbors must be found by the replicates method c if (.not. use_replica) return c c calculate interaction energy with other unit cells c do ii = 1, nvdw i = ivdw(ii) iv = ired(i) it = class(i) im = molcule(i) xi = xred(i) yi = yred(i) zi = zred(i) usei = (use(i) .or. use(iv)) do j = ii, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = ii, nvdw k = ivdw(kk) kv = ired(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 = class(k) km = molcule(k) do j = 2, ncell xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) call image (xr,yr,zr,j) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then rik = sqrt(rik2) rv = radmin(kt,it) eps = epsilon(kt,it) if (use_polymer) then if (rik2 .le. polycut2) then if (iv14(k) .eq. i) then rv = radmin4(kt,it) eps = epsilon4(kt,it) end if eps = eps * vscale(k) end if end if rv7 = rv**7 rik7 = rik**7 rho = rik7 + ghal*rv7 tau = (dhal+1.0d0) / (rik + dhal*rv) e = eps * rv7 * tau**7 & * ((ghal+1.0d0)*rv7/rho-2.0d0) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e end if c c use energy switching if near the cutoff distance c if (rik2 .gt. cut2) then rik = sqrt(rik2) rik3 = rik2 * rik rik4 = rik2 * rik2 rik5 = rik2 * rik3 taper = c5*rik5 + c4*rik4 + c3*rik3 & + c2*rik2 + c1*rik + c0 e = e * taper end if c c scale the interaction based on its group membership c if (use_group) e = e * fgrp c c increment the overall van der Waals energy component c if (e .ne. 0.0d0) then nev = nev + 1 if (i .eq. k) then ev = ev + 0.5d0*e aev(i) = aev(i) + 0.5d0*e else ev = ev + e aev(i) = aev(i) + 0.5d0*e aev(k) = aev(k) + 0.5d0*e end if end if c c increment the total intermolecular energy c einter = einter + e c c print a message if the energy of this interaction is large c huge = (e .gt. 10.0d0) if ((debug.and.e.ne.0.0d0) & .or. (verbose.and.huge)) then if (header) then header = .false. write (iout,30) 30 format (/,' Individual van der Waals', & ' Interactions :', & //,' Type',13x,'Atom Names', & 18x,'Minimum',4x,'Actual', & 6x,'Energy',/) end if write (iout,40) i,name(i),k,name(k), & rv,sqrt(rik2),e 40 format (' VDW-Hal',4x,i5,'-',a3,1x,i5,'-', & a3,3x,'(XTAL)',3x,2f10.4,f12.4) end if end if end do end if end do end do return end c c c ################################################################ c ## ## c ## subroutine ehal3b -- buffered 14-7 analysis via lights ## c ## ## c ################################################################ c c c "ehal3b" calculates the buffered 14-7 van der Waals energy c and also partitions the energy among the atoms using the c method of lights c c subroutine ehal3b (esol) implicit none include 'sizes.i' include 'action.i' include 'analyz.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'boxes.i' include 'cell.i' include 'couple.i' include 'energi.i' include 'group.i' include 'inform.i' include 'inter.i' include 'iounit.i' include 'light.i' include 'molcul.i' include 'mutant.i' include 'shunt.i' include 'usage.i' include 'vdw.i' include 'vdwpot.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer kgy,kgz integer start,stop integer im,km,ikmax,ikmin integer iv14(maxatm) real*8 e,rv,rv7 real*8 eps,rdn,fgrp real*8 xi,yi,zi real*8 xr,yr,zr real*8 rho,tau,taper real*8 rik,rik2,rik3 real*8 rik4,rik5,rik7 real*8 esol,erep real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) real*8 xsort(maxlight) real*8 ysort(maxlight) real*8 zsort(maxlight) logical proceed,usei logical prime,repeat logical header,huge c c c zero out the van der Waals energy and partitioning terms c nev = 0 ev = 0.0d0 do i = 1, n aev(i) = 0.0d0 end do header = .true. c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c apply any reduction factor to the atomic coordinates c do j = 1, nvdw i = ivdw(j) iv = ired(i) rdn = kred(i) xred(j) = rdn*(x(i)-x(iv)) + x(iv) yred(j) = rdn*(y(i)-y(iv)) + y(iv) zred(j) = rdn*(z(i)-z(iv)) + z(iv) end do c c transfer the interaction site coordinates to sorting arrays c do i = 1, nvdw xsort(i) = xred(i) ysort(i) = yred(i) zsort(i) = zred(i) end do c c use the method of lights to generate neighbors c call lights (off,nvdw,xsort,ysort,zsort) c c loop over all atoms computing the interactions c do ii = 1, nvdw i = ivdw(ii) iv = ired(i) it = class(i) im = molcule(i) xi = xsort(rgx(ii)) yi = ysort(rgy(ii)) zi = zsort(rgz(ii)) usei = (use(i) .or. use(iv)) do j = 1, nvdw vscale(ivdw(j)) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do if (kbx(ii) .le. kex(ii)) then repeat = .false. start = kbx(ii) + 1 stop = kex(ii) else repeat = .true. start = 1 stop = kex(ii) end if 10 continue do j = start, stop kk = locx(j) kgy = rgy(kk) if (kby(ii) .le. key(ii)) then if (kgy.lt.kby(ii) .or. kgy.gt.key(ii)) goto 50 else if (kgy.lt.kby(ii) .and. kgy.gt.key(ii)) goto 50 end if kgz = rgz(kk) if (kbz(ii) .le. kez(ii)) then if (kgz.lt.kbz(ii) .or. kgz.gt.kez(ii)) goto 50 else if (kgz.lt.kbz(ii) .and. kgz.gt.kez(ii)) goto 50 end if k = ivdw(kk-((kk-1)/nvdw)*nvdw) kv = ired(k) prime = (kk .le. nvdw) c c decide whether to compute the current interaction c 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 = class(k) km = molcule(k) xr = xi - xsort(j) yr = yi - ysort(kgy) zr = zi - zsort(kgz) if (use_image) then if (abs(xr) .gt. xcell2) xr = xr - sign(xcell,xr) if (abs(yr) .gt. ycell2) yr = yr - sign(ycell,yr) if (abs(zr) .gt. zcell2) zr = zr - sign(zcell,zr) if (monoclinic) then xr = xr + zr*beta_cos zr = zr * beta_sin else if (triclinic) then xr = xr + yr*gamma_cos + zr*beta_cos yr = yr*gamma_sin + zr*beta_term zr = zr * gamma_term end if end if rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) then rik = sqrt(rik2) rv = radmin(kt,it) eps = epsilon(kt,it) if (prime) then if (iv14(k) .eq. i) then rv = radmin4(kt,it) eps = epsilon4(kt,it) end if eps = eps * vscale(k) end if rv7 = rv**7 rik7 = rik**7 rho = rik7 + ghal*rv7 tau = (dhal+1.0d0) / (rik + dhal*rv) e = eps * rv7 * tau**7 & * ((ghal+1.0d0)*rv7/rho-2.0d0) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 e = e * taper end if c c scale the interaction based on its group membership c if (use_group) e = e * fgrp c c increment the overall van der Waals energy components c if (e .ne. 0.0d0) then nev = nev + 1 ev = ev + e aev(i) = aev(i) + 0.5d0*e aev(k) = aev(k) + 0.5d0*e end if c c increment the total intermolecular energy c if (.not.prime .or. molcule(i).ne.molcule(k)) then einter = einter + e end if c c print a message if the energy of this interaction is large c huge = (e .gt. 10.0d0) if ((debug.and.e.ne.0.0d0) & .or. (verbose.and.huge)) then if (header) then header = .false. write (iout,20) 20 format (/,' Individual van der Waals', & ' Interactions :', & //,' Type',13x,'Atom Names', & 18x,'Minimum',4x,'Actual', & 6x,'Energy',/) end if if (prime) then write (iout,30) i,name(i),k,name(k), & rv,sqrt(rik2),e 30 format (' VDW-Hal',4x,i5,'-',a3,1x,i5,'-', & a3,12x,2f10.4,f12.4) else write (iout,40) i,name(i),k,name(k), & rv,sqrt(rik2),e 40 format (' VDW-Hal',4x,i5,'-',a3,1x,i5,'-', & a3,3x,'(XTAL)',3x,2f10.4,f12.4) end if end if end if end if 50 continue end do if (repeat) then repeat = .false. start = kbx(ii) + 1 stop = nlight goto 10 end if end do return end c c c ############################################################## c ## ## c ## subroutine ehal3c -- buffered 14-7 analysis via list ## c ## ## c ############################################################## c c c "ehal3c" calculates the buffered 14-7 van der Waals energy c and also partitions the energy among the atoms using a c pairwise neighbor list c c subroutine ehal3c (esol) implicit none include 'sizes.i' include 'action.i' include 'analyz.i' include 'atmtyp.i' include 'atoms.i' include 'bound.i' include 'couple.i' include 'energi.i' include 'group.i' include 'inform.i' include 'inter.i' include 'iounit.i' include 'molcul.i' include 'mutant.i' include 'neigh.i' include 'shunt.i' include 'usage.i' include 'vdw.i' include 'vdwpot.i' integer i,j,k integer ii,iv,it integer kk,kv,kt integer im,km,ikmax,ikmin integer iv14(maxatm) real*8 e,eps,rdn real*8 fgrp,rv,rv7 real*8 xi,yi,zi real*8 xr,yr,zr real*8 rho,tau,taper real*8 rik,rik2,rik3 real*8 rik4,rik5,rik7 real*8 esol,erep real*8 xred(maxatm) real*8 yred(maxatm) real*8 zred(maxatm) real*8 vscale(maxatm) logical proceed,usei logical header,huge c c c zero out the van der Waals energy and partitioning terms c nev = 0 ev = 0.0d0 do i = 1, n aev(i) = 0.0d0 end do header = .true. c c zero out the array marking 1-4 vdw interactions c do i = 1, n iv14(i) = 0 end do c c set the coefficients for the switching function c call switch ('VDW') c c apply any reduction factor to the atomic coordinates c do k = 1, nvdw i = ivdw(k) iv = ired(i) rdn = kred(i) xred(i) = rdn*(x(i)-x(iv)) + x(iv) yred(i) = rdn*(y(i)-y(iv)) + y(iv) zred(i) = rdn*(z(i)-z(iv)) + z(iv) end do c c find the van der Waals energy via neighbor list search c do ii = 1, nvdw i = ivdw(ii) iv = ired(i) it = class(i) im = molcule(i) xi = xred(i) yi = yred(i) zi = zred(i) usei = (use(i) .or. use(iv)) do j = 1, nvlst(ii) vscale(ivdw(vlst(j,ii))) = 1.0d0 end do do j = 1, n12(i) vscale(i12(j,i)) = v2scale end do do j = 1, n13(i) vscale(i13(j,i)) = v3scale end do do j = 1, n14(i) vscale(i14(j,i)) = v4scale iv14(i14(j,i)) = i end do do j = 1, n15(i) vscale(i15(j,i)) = v5scale end do c c decide whether to compute the current interaction c do kk = 1, nvlst(ii) k = ivdw(vlst(kk,ii)) kv = ired(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 = class(k) km = molcule(k) xr = xi - xred(k) yr = yi - yred(k) zr = zi - zred(k) if (use_image) call image (xr,yr,zr,0) rik2 = xr*xr + yr*yr + zr*zr c c check for an interaction distance less than the cutoff c if (rik2 .le. off2) 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) rv7 = rv**7 rik7 = rik**7 rho = rik7 + ghal*rv7 tau = (dhal+1.0d0) / (rik + dhal*rv) e = eps * rv7 * tau**7 & * ((ghal+1.0d0)*rv7/rho-2.0d0) c c use Weeks-Chandler-Andersen for solute-solvent interaction c ikmax = max(im,km) ikmin = min(im,km) if (ikmax.eq.nmol .and. ikmin.ne.nmol) then if (rik .lt. rv) then erep = e + eps else erep = 0.0d0 end if e = e + lambda*(erep-e) esol = esol + e 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 e = e * taper end if c c scale the interaction based on its group membership c if (use_group) e = e * fgrp c c increment the overall van der Waals energy components c if (e .ne. 0.0d0) then nev = nev + 1 ev = ev + e aev(i) = aev(i) + 0.5d0*e aev(k) = aev(k) + 0.5d0*e end if c c increment the total intermolecular energy c if (molcule(i) .ne. molcule(k)) then einter = einter + e end if c c print a message if the energy of this interaction is large c huge = (e .gt. 10.0d0) if ((debug.and.e.ne.0.0d0) & .or. (verbose.and.huge)) then if (header) then header = .false. write (iout,10) 10 format (/,' Individual van der Waals', & ' Interactions :', & //,' Type',13x,'Atom Names', & 18x,'Minimum',4x,'Actual', & 6x,'Energy',/) end if write (iout,20) i,name(i),k,name(k), & rv,sqrt(rik2),e 20 format (' VDW-Hal',4x,i5,'-',a3,1x,i5,'-', & a3,12x,2f10.4,f12.4) end if end if end if end do end do return end