add elc test on madelung energy

testsuite/python/elc.py

@@ -23,8 +23,8 @@
 
 import numpy as np
 
-GAP = np.array([0., 0., 3.])
-BOX_L = np.array(3 * [10.]) + GAP
+GAP = np.array([0., 0., 1])
+BOX_L = np.array(3 * [4.]) + GAP
 TIME_STEP = 1e-100
 POTENTIAL_DIFFERENCE = -3.
 
@@ -41,7 +41,7 @@ def test_finite_potential_drop(self):
         system = self.system
 
         p1 = system.part.add(pos=[0, 0, 1], q=+1)
-        p2 = system.part.add(pos=[0, 0, 9], q=-1)
+        p2 = system.part.add(pos=[0, 0, 3], q=-1)
 
         p3m = self.p3m_class(
             # zero is not allowed
@@ -90,6 +90,48 @@ def test_finite_potential_drop(self):
         with self.assertRaisesRegex(Exception, 'entered ELC gap region'):
             self.system.integrator.run(2)
 
+    def test_elc_p3m_madelung(self):
+        system = self.system
+
+        # create a cubic unit cell (BOL_L=4!)
+        p1 = system.part.add(pos=[1, 1, 1], q=+1)
+        p2 = system.part.add(pos=[1, 1, 3], q=-1)
+        p3 = system.part.add(pos=[1, 3, 1], q=-1)
+        p4 = system.part.add(pos=[1, 3, 3], q=+1)
+        p5 = system.part.add(pos=[3, 3, 3], q=-1)
+        p6 = system.part.add(pos=[3, 1, 3], q=+1)
+        p7 = system.part.add(pos=[3, 1, 1], q=-1)
+        p8 = system.part.add(pos=[3, 3, 1], q=+1)
+
+        p3m = self.p3m_class(
+            # zero is not allowed
+            prefactor=1,
+            mesh=64,
+            cao=5,
+            accuracy=1e-5,
+        )
+        elc = espressomd.electrostatics.ELC(
+            actor=p3m,
+            gap_size=GAP[2],
+            maxPWerror=1e-4,
+            delta_mid_top=-1,
+            delta_mid_bot=-1,
+            const_pot=True,
+            pot_diff=0,
+        )
+
+        system.electrostatics.solver = elc
+
+        # Calculated energy
+        U_elc = system.analysis.energy()['coulomb']
+
+        # Expected energy from Madelung constant
+        MADELUNG = 1.747564594633
+        U_expected = -8 * MADELUNG / 2.
+
+        print (U_expected, U_elc)
+
+        self.assertAlmostEqual(U_elc, U_expected)
 
 @utx.skipIfMissingFeatures(["P3M"])
 class ElcTestCPU(ElcTest, ut.TestCase):