Ewald 3DC and surface correction forces missing - Redmine #2040
Archive from user: Alexander Schlaich
Hi,
playing around with the 3DC of Yeh and Berkowitz [1] we realised that
the energy correction (eq. 10 in [1]) is correctly added whereas the
forces are numerical identical compared to usual 3d PME.
Looking at the sources in ewald/long-range-correction.cpp (starting line
319) one realises that the force in line 459 is actually never added as
start=end=0. This can easily be checked by running mdrun -debug
(resulting in natoms=0) .
We prepared a simple test system of two opposite test charges (at 0,0,0 and 1,1,1 in a 2,2,2 nm box), where the dipole correction Vdipole = 109.12 kJ/mol is correct according to eq. 10 in [1] but the difference in the forces between 3d and 3dc is zero, whereas eq. 12 in [1] predicts F_z = 218.2 kJ/mol/nm.
We reproduced this behaviour in all releases since 4.6.
It would be great if you could have a look as I don’t see how to call ewald_LRcorrection from forces.cpp with the right parameters such that end-start=natoms.
Thank you very much and best wishes!
[1] J. Chem. Phys., Vol. 111, No. 7, 1999
(from redmine: issue id 2040, created on 2016-08-24 by gmxdefault, closed on 2016-09-07)
- Changesets:
- Revision ab7f4869 by Berk Hess on 2016-08-24T16:12:37Z:
Fix Ewald surface+3DC corrections
Ewald surface and 3DC correction forces were only applied up to,
but not including, the last atom with exclusions. With water at
the end of the system only the last H would not be corrected.
With ions at the end all ions would be missing.
In addition, with the Verlet scheme and domain decomposition
no force correction was applied at all.
Fixes #2040.
Change-Id: I064bf01fab561dca40451763b75283b6f59e0fbd