From 84554b0442aa72be72fa68c0ee1a4a9371e8bfae Mon Sep 17 00:00:00 2001 From: Andy Simmonett Date: Fri, 12 May 2017 10:14:18 -0600 Subject: [PATCH] Fix reciprocal polarization energy --- source/epolar.f | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/source/epolar.f b/source/epolar.f index 2c2056cb8..4b9a48a0f 100644 --- a/source/epolar.f +++ b/source/epolar.f @@ -1759,7 +1759,7 @@ subroutine eprecip integer m1,m2,m3 integer ntot,nff integer nf1,nf2,nf3 - real*8 e,r1,r2,r3 + real*8 e,r1,r2,r3,f real*8 h1,h2,h3 real*8 volterm,denom real*8 hsq,expterm @@ -1769,6 +1769,7 @@ subroutine eprecip real*8, allocatable :: fuind(:,:) real*8, allocatable :: fuinp(:,:) c + f = electric / dielec c c return if the Ewald coefficient is zero c @@ -1890,7 +1891,7 @@ subroutine eprecip call fphi_mpole (fphi) do i = 1, npole do j = 1, 20 - fphi(j,i) = electric * fphi(j,i) + fphi(j,i) = f * fphi(j,i) end do end do end if @@ -1926,7 +1927,7 @@ subroutine eprecip if (.not. use_bounds) then expterm = 0.5d0 * pi / xbox struc2 = qgrid(1,1,1,1)**2 + qgrid(2,1,1,1)**2 - e = 0.5d0 * electric * expterm * struc2 + e = 0.5d0 * f * expterm * struc2 ep = ep + e end if c @@ -1939,7 +1940,7 @@ subroutine eprecip e = e + fuind(k,i)*fphi(k+1,i) end do end do - e = 0.5d0 * e + e = 0.5d0 * f * e ep = ep + e c c perform deallocation of some local arrays