From bf4890824cc7f3fb1a63506b1eb012cbac531aa4 Mon Sep 17 00:00:00 2001
From: Alexander Reinauer <areinauer@icp.uni-stuttgart.de>
Date: Fri, 8 Dec 2023 15:51:45 +0100
Subject: [PATCH] testsuite: ELC-IC: added Madelung-energy testcase

Co-authored-by: Alexander Schlaich <schlaich@icp.uni-stuttgart.de>
---
 testsuite/python/elc.py | 56 +++++++++++++++++++++++++++++++++++++----
 1 file changed, 51 insertions(+), 5 deletions(-)

diff --git a/testsuite/python/elc.py b/testsuite/python/elc.py
index ec35086b96..f13d940e2f 100644
--- a/testsuite/python/elc.py
+++ b/testsuite/python/elc.py
@@ -1,5 +1,5 @@
 #
-# Copyright (C) 2010-2022 The ESPResSo project
+# Copyright (C) 2010-2023 The ESPResSo project
 #
 # This file is part of ESPResSo.
 #
@@ -22,15 +22,13 @@
 import espressomd.electrostatics
 
 import numpy as np
+import itertools
 
-GAP = np.array([0., 0., 3.])
-BOX_L = np.array(3 * [10.]) + GAP
 TIME_STEP = 1e-100
-POTENTIAL_DIFFERENCE = -3.
 
 
 class ElcTest:
-    system = espressomd.System(box_l=BOX_L, time_step=TIME_STEP)
+    system = espressomd.System(box_l=[1.] * 3, time_step=TIME_STEP)
     system.cell_system.skin = 0.0
 
     def tearDown(self):
@@ -40,6 +38,12 @@ def tearDown(self):
     def test_finite_potential_drop(self):
         system = self.system
 
+        GAP = np.array([0., 0., 3.])
+        BOX_L = np.array(3 * [10.]) + GAP
+        POTENTIAL_DIFFERENCE = -3.
+
+        system.box_l = BOX_L
+
         p1 = system.part.add(pos=[0, 0, 1], q=+1)
         p2 = system.part.add(pos=[0, 0, 9], q=-1)
 
@@ -90,6 +94,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
+
+        n_pairs = 6
+
+        BOX_L = 2 * n_pairs
+        ELC_GAP = 0.4 * BOX_L
+        system.box_l = [BOX_L, BOX_L, BOX_L + ELC_GAP]
+
+        for j, k, l in itertools.product(range(2 * n_pairs), repeat=3):
+            system.part.add(pos=[j + 0.5, k + 0.5, l + 0.5],
+                            q=(-1)**(j + k + l))
+
+        p3m = self.p3m_class(
+            prefactor=1,
+            mesh=[60, 60, 84],
+            cao=7,
+            r_cut=3.314,
+            alpha=1.18,
+            accuracy=3e-7,
+            tune=False,
+        )
+        elc = espressomd.electrostatics.ELC(
+            actor=p3m,
+            gap_size=ELC_GAP,
+            maxPWerror=3e-7,
+            delta_mid_top=-1,
+            delta_mid_bot=-1,
+            const_pot=True,
+            pot_diff=0,
+        )
+
+        system.electrostatics.solver = elc
+
+        MADELUNG = -1.74756459463318219
+        U_expected = MADELUNG
+
+        U_elc = 2. * \
+            system.analysis.energy()['coulomb'] / len(system.part.all())
+
+        np.testing.assert_allclose(U_elc, U_expected, atol=0., rtol=1e-6)
+
 
 @utx.skipIfMissingFeatures(["P3M"])
 class ElcTestCPU(ElcTest, ut.TestCase):