ELC energy wrong for non-neutral systems

[KonradBreitsprecher]

Dragging a particle through the box in a non-neutral system (testscript below), the ELC energy jumps close to the boundaries (inside the elc_params.space_layer):

The current testcase (elc_vs_mmm2d_nonneutral.py) fails to cover this case because the particles don't enter the space layer (not dragged close enough to the walls).

I made the following observations:

• Forces are ok (agree with MMM2D)
• Neutral systems are ok
• The error has be related to one of the many parts, where additional contributions from the space layer appear (p.r.p[2] < elc_params.space_layer).
• I could not narrow down the error by finding the single contribution that causes a discontinuity. In fact, several parts of the energy calc. cause jumps at the space layer, and probably fail to cancel each other for non-neutral system.
• Another hot candidate is the part of the energy calculation in coulomb.cpp, where ELC messes with the system charges and P3M sums.

I don't see through the ELC energy calculation (both non-neutral and neutral) and
I'm afraid the best way out of this is to not support ELC for non-neutral systems.

Test script:

from __future__ import print_function
import unittest as ut
import espressomd
import numpy as np
import espressomd.electrostatics
from espressomd import electrostatic_extensions


class ELC_vs_MMM2D_neutral(ut.TestCase):
    # Handle to espresso system

    system = espressomd.System(box_l=[1.0, 1.0, 1.0])
    acc = 1e-6
    elc_gap = 8.0
    box_l = 10.0
    bl2 = box_l * 0.5
    system.time_step = 0.01
    system.cell_system.skin = 0.1

    def test_elc_vs_mmm2d(self):
        elc_param_sets = {
            "inert": {"gap_size": self.elc_gap, "maxPWerror": self.acc, "check_neutrality": False},
            "dielectric": {
                "gap_size": self.elc_gap,
                "maxPWerror": self.acc,
                "delta_mid_bot": 0.1,
                "delta_mid_top": 0.9,
                "check_neutrality": False},
            "const_pot_0": {
                "gap_size": self.elc_gap,
                "maxPWerror": self.acc,
                "const_pot": 1,
                "pot_diff": 0.0,
                "check_neutrality": False},
            "const_pot_1": {
                "gap_size": self.elc_gap,
                "maxPWerror": self.acc,
                "const_pot": 1,
                "pot_diff": 1.0,
                "check_neutrality": False},
            "const_pot_m1": {
                "gap_size": self.elc_gap, 
                "maxPWerror": self.acc, 
                "const_pot": 1, 
                "pot_diff": -1.0, 
                "check_neutrality": False},
        }

        mmm2d_param_sets = {
            "inert": {"prefactor": 1.0, "maxPWerror": self.acc, "check_neutrality": False},
            "dielectric": {
                "prefactor": 1.0,
                "maxPWerror": self.acc,
                "dielectric_contrast_on": 1,
                "delta_mid_bot": 0.1,
                "delta_mid_top": 0.9,
                "check_neutrality": False},
            "const_pot_0": {
                "prefactor": 1.0,
                "maxPWerror": self.acc,
                "const_pot": 1,
                "pot_diff": 0.0,
                "check_neutrality": False},
            "const_pot_1": {
                "prefactor": 1.0,
                "maxPWerror": self.acc,
                "const_pot": 1,
                "pot_diff": 1.0,
                "check_neutrality": False},
            "const_pot_m1": {
                "prefactor": 1.0, 
                "maxPWerror": self.acc, 
                "const_pot": 1, 
                "pot_diff": -1.0, 
                "check_neutrality": False},
        }

        self.system.box_l = [self.box_l, self.box_l, self.box_l]
        buf_node_grid = self.system.cell_system.node_grid
        self.system.cell_system.set_layered(
            n_layers=10, use_verlet_lists=False)

        self.system.periodicity = [1, 1, 0]

        q = 3.0
        non_neutral_fac = 3.0

        self.system.part.add(id=0, pos=(5.0, 5.0, 5.0), q=-non_neutral_fac*q)
        self.system.part.add(id=1, pos=(2.0, 2.0, 5.0), q=q / 3.0)
        self.system.part.add(id=2, pos=(2.0, 5.0, 2.0), q=q / 3.0)
        self.system.part.add(id=3, pos=(5.0, 2.0, 7.0), q=q / 3.0)

        #case="inert"
        #case="dielectric"
        case="const_pot_0"
        #case="const_pot_1"

        #MMM2D
        mmm2d = espressomd.electrostatics.MMM2D(**mmm2d_param_sets[case])
        self.system.actors.add(mmm2d)
        mmm2d_res = {}

        mmm2d_res[case] = self.scan()
        
        self.system.actors.remove(mmm2d)

        #ELC
        self.system.box_l = [self.box_l, self.box_l, self.box_l + self.elc_gap]
        self.system.cell_system.set_domain_decomposition(
            use_verlet_lists=True)
        self.system.cell_system.node_grid = buf_node_grid
        self.system.periodicity = [1, 1, 1]
        #p3m = espressomd.electrostatics.P3M(prefactor=1.0, accuracy=self.acc, mesh=[16, 16, 24], cao=6, check_neutrality=False)
        p3m = espressomd.electrostatics.P3M(prefactor=1.0, accuracy=self.acc, check_neutrality=False)
        self.system.actors.add(p3m)

        elc = electrostatic_extensions.ELC(**elc_param_sets[case])
        self.system.actors.add(elc)
        elc_res = {}

        elc_res[case] = self.scan()
        
        np.savetxt('data.dat',np.hstack((elc_res[case],mmm2d_res[case])))


    def scan(self):
        n = 100
        d = 0.05
        res = []
        for i in range(n + 1):
            z = self.box_l - d - 1.0 * i / n * (self.box_l - 2 * d)
            self.system.part[0].pos = [self.bl2, self.bl2, z]
            self.system.integrator.run(0)
            energy = self.system.analysis.energy()
            m = [z]
            m.extend(self.system.part[0].f)
            m.append(energy['coulomb'])
            res.append(m)

        return res

if __name__ == "__main__":
    ut.main()

[jonaslandsgesell]

In the following comes some speculation: To me the jump looks like being caused by a wrong scaling factor within the outer layers. The trend of the curve looks similar, but scaled.

How does this jump scale when you change the prefactor?

PS: thanks for reporting this problem! This energy problem (as long as it persists) prevents using ELC in a nonneutral system together with the cpH or reaction ensemble (because they rely on energy changes)

[RudolfWeeber]

Offline discussion with @reinaual:
No solution was found inspite of considerable effort.
We're going to block the affected combintaion of features.

[RudolfWeeber]

@reinaual, can you please open a Pr which blocks the activation of ELC with the non-working parameter combinations.
Please also add entries to the ES 4.1 release notes in the Wiki, explaining what you fixed and what wasn't fixed.

[RudolfWeeber]

acitvating the affected featues was blocked. De-milestoning this issue.