testsuite: ELC-IC: added Madelung-energy testcase

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):