Merge #3310

3310: core: Fixed ELC potential difference prefactor + test r=RudolfWeeber a=fweik

Description of changes:
 - When using ELC/MMM2d with potential boundary conditions, one
   of the correction factors was over-counted resulting in wrong energies.
 
The other term in the correction in L 348 of elc.cpp also looks suspicious,
but I could not find a straight forward way to test that. All of this code should
only be used with the utmost caution.

Co-authored-by: Florian Weik <[email protected]>

src/core/electrostatics_magnetostatics/elc.cpp

@@ -347,7 +347,7 @@ static double dipole_energy(const ParticleRange &particles) {
       // zero potential difference contribution
       eng += pref * height_inverse / uz * Utils::sqr(gblcblk[6]);
       // external potential shift contribution
-      eng -= elc_params.pot_diff * height_inverse * gblcblk[6];
+      eng -= 2 * elc_params.pot_diff * height_inverse * gblcblk[6];
     }
 
     /* counter the P3M homogeneous background contribution to the

src/core/electrostatics_magnetostatics/mmm2d.cpp

@@ -641,7 +641,7 @@ static double z_energy(const ParticleRange &particles) {
       // zero potential difference contribution
       eng += gbl_dm_z * gbl_dm_z * coulomb.prefactor * 2 * M_PI * ux * uy * uz;
       // external potential shift contribution
-      eng -= mmm2d_params.pot_diff * uz * gbl_dm_z;
+      eng -= 2. * mmm2d_params.pot_diff * uz * gbl_dm_z;
     }
   }
 

testsuite/python/CMakeLists.txt

@@ -206,6 +206,7 @@ python_test(FILE time_series.py MAX_NUM_PROC 1)
 python_test(FILE linear_momentum.py)
 python_test(FILE linear_momentum_lb.py MAX_NUM_PROC 2 LABELS gpu)
 python_test(FILE mmm1d.py MAX_NUM_PROC 2 LABELS gpu)
+python_test(FILE elc.py MAX_NUM_PROC 2)
 
 if(PY_H5PY)
   python_test(FILE h5md.py MAX_NUM_PROC 2)

testsuite/python/elc.py

@@ -0,0 +1,76 @@
+# Copyright (C) 2010-2019 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+import unittest as ut
+import unittest_decorators as utx
+import espressomd
+from espressomd import electrostatics, electrostatic_extensions
+
+import numpy as np
+
+GAP = np.array([0, 0, 3.])
+BOX_L = np.array(3 * [10]) + GAP
+TIME_STEP = 1e-100
+POTENTIAL_DIFFERENCE = -3.
+
+
+@utx.skipIfMissingFeatures(["P3M"])
+class ElcTest(ut.TestCase):
+    system = espressomd.System(box_l=BOX_L, time_step=TIME_STEP)
+    system.cell_system.skin = 0.0
+
+    def test_finite_potential_drop(self):
+        s = self.system
+
+        p1 = s.part.add(pos=GAP + [0, 0, 1], q=+1)
+        p2 = s.part.add(pos=GAP + [0, 0, 9], q=-1)
+
+        s.actors.add(
+            electrostatics.P3M(
+                # zero is not allowed
+                prefactor=1e-100,
+                mesh=32,
+                cao=5,
+                accuracy=1e-3,
+            ))
+
+        s.actors.add(
+            electrostatic_extensions.ELC(
+                gap_size=GAP[2],
+                maxPWerror=1e-3,
+                delta_mid_top=-1,
+                delta_mid_bot=-1,
+                const_pot=1,
+                pot_diff=POTENTIAL_DIFFERENCE,
+            ))
+
+        # Calculated energy
+        U_elc = s.analysis.energy()['coulomb']
+
+        # Expected E-Field is voltage drop over the box
+        E_expected = POTENTIAL_DIFFERENCE / (BOX_L[2] - GAP[2])
+        # Expected potential is -E_expected * z, so
+        U_expected = -E_expected * (p1.pos[2] * p1.q + p2.pos[2] * p2.q)
+
+        self.assertAlmostEqual(U_elc, U_expected)
+
+        s.integrator.run(0)
+        self.assertAlmostEqual(E_expected, p1.f[2] / p1.q)
+        self.assertAlmostEqual(E_expected, p2.f[2] / p2.q)
+
+
+if __name__ == "__main__":
+    ut.main()