is a prefactor which can be set by the user.
The commonly used Bdrm length \(l_B = e_o^2 / (4 \pi \epsilon_0 \epsilon_r k_B T)\) is the length at which the Coulomb energy between two unit charges is equal to the thermal energy \(k_B T\).
Based on the this length, the prefactor is given by \(C=l_B k_B T\).

Computing electrostatic interactions is computationally very expensive.
ESPResSo features some state-of-the-art algorithms to deal with these
interactions as efficiently as possible, but almost all of them require
some knowledge to use them properly. Uneducated use can result in
completely unphysical simulations.

Coulomb interactions have to be added to the list of active actors of the system object to become
active. This is done by calling the add-method of espressomd.system.System.actors.
Only one electrostatics method can be active at any time.

This example shows the general usage of an electrostatic method <SOLVER>.
All of them need the Bjerrum length and a set of other required parameters.
First, an instance of the solver is created and only after adding it to the actors
list, it is activated. Internally the method calls a tuning routine on
activation to achieve the given accuracy:

For this feature to work, you need to have the fftw3 library
installed on your system. In ESPResSo, you can check if it is compiled in by
checking for the feature FFTW with espressomd.features.
P3M requires full periodicity (1 1 1). Make sure that you know the relevance of the
P3M parameters before using P3M! If you are not sure, read the following
references
[Ewa21][HE88][KP92][DH98a][DH98b][DHL00][Des00][CBLCHolm08].

The tuning method is called when the handle of the Coulomb P3M is added to the
actor list. At this point, the system should already contain the charged
particles. Set parameters are fixed and not changed by the tuning algorithm.
This can be useful to speed up the tuning during testing or if the parameters
are already known.

To prevent the automatic tuning, set the tune parameter to False.
To manually tune or retune P3M, call espresso.electrostatics.P3M.Tune().
Note, however, that this is a method the P3M object inherited from
espressomd.electrostatics.ElectrostaticInteraction.
All parameters passed to the method are fixed in the tuning routine. If not
specified in the Tune() method, the parameters prefactor and
accuracy are reused.

It is not easy to calculate the various parameters of the P3M method
such that the method provides the desired accuracy at maximum speed. To
simplify this, it provides a function to automatically tune the algorithm.
Note that for this function to work properly, your system should already
contain an initial configuration of charges and the correct initial box
size. Also note that the provided tuning algorithms works very well on
homogeneous charge distributions, but might not achieve the requested
precision for highly inhomogeneous or symmetric systems. For example,
because of the nature of the P3M algorithm, systems are problematic
where most charges are placed in one plane, one small region, or on a
regular grid.

The function employs the analytical expression of the error estimate for
the P3M method [HE88] and its real space error [KP92] to
obtain sets of parameters that yield the desired accuracy, then it measures how
long it takes to compute the Coulomb interaction using these parameter sets and
chooses the set with the shortest run time.

After execution the tuning routines report the tested parameter sets,
the corresponding k-space and real-space errors and the timings needed
for force calculations. In the output, the timings are given in units of
milliseconds, length scales are in units of inverse box lengths.

The GPU implementation of P3M calculates the far field portion on the GPU.
It uses the same parameters and interface functionality as the CPU version of
the solver. It should be noted that this does not always provide significant
increase in performance. Furthermore it computes the far field interactions
with only single precision which limits the maximum precision. The algorithm
does not work in combination with the electrostatic extensions Dielectric interfaces with the ICC\star algorithm
and Electrostatic Layer Correction (ELC).

where \(C\) is defined as in Eqn. (1).
The Debye-Hückel potential is an approximate method for calculating
electrostatic interactions, but technically it is treated as other
short-ranged non-bonding potentials. For \(r>r_{\mathrm cut}\) it is
set to zero which introduces a step in energy. Therefore, it introduces
fluctuations in energy.

The ICC\(\star\) algorithm allows to take into account arbitrarily shaped
dielectric interfaces and dynamic charge induction. For instance, it can be
used to simulate a curved metallic boundary. This is done by iterating the
charge on a set of spatially fixed ICC particles until they correctly
represent the influence of the dielectric discontinuity. All ICC particles
need a certain area, normal vector and dielectric constant to specify the
surface. ICC relies on a Coulomb solver that is already initialized. So far, it
is implemented and well tested with the Coulomb solver P3M. ICC is an ESPResSo
actor and can be activated via:

icc=ICC(<SeethefollowinglistofICCparameters>)system.actors.add(icc)

Parameters are:

first_id:

ID of the first ICC Particle.

n_icc:

Total number of ICC Particles.

convergence:

Abort criteria of the iteration. It corresponds to the maximum relative
change of any of the interface particle’s charge.

relaxation:

SOR relaxation parameter.

ext_field:

Homogeneous electric field added to the calculation of dielectric boundary forces.

max_iterations:

Maximal number of iterations.

eps_out:

Relative permittivity of the outer region (where the particles are).

normals:

List of size n_icc with normal vectors pointing into the outer region.

areas

List of size n_icc with areas of the discretized surface.

sigmas

List of size n_icc with an additional surface charge density in
absence of any charge induction

epsilons

List of size n_icc with the dielectric constant associated to the area.

The ICC particles are setup as normal ESPResSo particles. Note that they should be
fixed in space and need an initial nonzero charge. The following usage example
sets up parallel metallic plates and activates ICC:

# Set the ICC line density and calculate the number of# ICC particles according to the box sizel=3.2nicc=int(box_l/l)nicc_per_electrode=nicc*niccnicc_tot=2*nicc_per_electrodeiccArea=box_l*box_l/nicc_per_electrodel=box_l/nicc# Lists to collect required parametersiccNormals=[]iccAreas=[]iccSigmas=[]iccEpsilons=[]# Add the fixed ICC particles:# Left electrode (normal [0,0,1])forxiinxrange(nicc):foryiinxrange(nicc):system.part.add(pos=[l*xi,l*yi,0],q=-0.0001,fix=[1,1,1],type=icc_type)iccNormals.extend([0,0,1]*nicc_per_electrode)# Right electrode (normal [0,0,-1])forxiinxrange(nicc):foryiinxrange(nicc):system.part.add(pos=[l*xi,l*yi,box_l],q=0.0001,fix=[1,1,1],type=icc_type)iccNormals.extend([0,0,-1]*nicc_per_electrode)# Common area, sigma and metallic epsiloniccAreas.extend([iccArea]*nicc_tot)iccSigmas.extend([0]*nicc_tot)iccEpsilons.extend([100000]*nicc_tot)icc=ICC(first_id=0,n_icc=nicc_tot,convergence=1e-4,relaxation=0.75,ext_field=[0,0,0],max_iterations=100,eps_out=1.0,normals=iccNormals,areas=iccAreas,sigmas=iccSigmas,epsilons=iccEpsilons)system.actors.add(icc)

With each iteration, ICC has to solve electrostatics which can severely slow
down the integration. The performance can be improved by using multiple cores,
a minimal set of ICC particles and convergence and relaxation parameters that
result in a minimal number of iterations. Also please make sure to read the
corresponding articles, mainly [ALK+ed][TSuzenS+10][KSH11] before
using it.

MMM2D is an electrostatics solver for explicit 2D periodic systems.
It can account for different dielectric jumps on both sides of the
non-periodic direction. MMM2D Coulomb method needs periodicity 1 1 0 and the
layered cell system. The performance of the method depends on the number of
slices of the cell system, which has to be tuned manually. It is
automatically ensured that the maximal pairwise error is smaller than
the given bound. Note that the user has to take care that the particles don’t
leave the box in the non-periodic z-direction e.g. with constraints. By default,
no dielectric contrast is set and it is used as:

where the prefactor \(C\) is defined in Eqn. (1).
For a detailed list of parameters see espressomd.electrostatics.MMM2D.
The last two, mutually exclusive parameters dielectric and
dielectric_constants_on allow to specify dielectric contrasts at the
upper and lower boundaries of the simulation box. The first form
specifies the respective dielectric constants in the media, which
however is only used to calculate the contrasts. That is, specifying
\(\epsilon_t=\epsilon_m=\epsilon_b=\text{const}\) is always
identical to \(\epsilon_t=\epsilon_m=\epsilon_b=1\):

The second form specifies only the dielectric contrasts at the boundaries,
that is \(\Delta_t=\frac{\epsilon_m-\epsilon_t}{\epsilon_m+\epsilon_t}\)
and \(\Delta_b=\frac{\epsilon_m-\epsilon_b}{\epsilon_m+\epsilon_b}\).
Using this form allows to choose \(\Delta_{t/b}=-1\), corresponding
to metallic boundary conditions:

This is done by countering the total dipole moment of the system with the
electric field \(E_{induced}\) and superposing a homogeneous electric field
\(E_{applied} = \frac{U}{L}\) to retain \(U\). This mimics the
induction of surface charges \(\pm\sigma = E_{induced} \cdot \epsilon_0\)
for planar electrodes at \(z=0\) and \(z=L\) in a capacitor connected
to a battery with voltage pot_diff. Using 0 is equivalent to
\(\Delta_{t/b}=-1\).

Finally, the far cutoff setting should only be used for testing reasons,
otherwise you are more safe with the automatic tuning. If you even don’t know
what it is, do not even think of touching the far cutoff. For details on the
MMM family of algorithms, refer to appendix The MMM family of algorithms.
Please cite [AH02b] when using MMM2D.

A complete (but unphysical) sample script for a plate capacitor simulated with MMM2D
can be found in /samples/visualization_mmm2d.py.

ELC can be used to simulate charged system with 2D periodicity. In more
detail, is a special procedure that converts a 3D electrostatic method to a 2D
method in computational order N. Currently, it only supports P3M. This means,
that you will first have to set up the P3M algorithm before using ELC. The
algorithm is definitely faster than MMM2D for larger numbers of particles
(\(>400\) at reasonable accuracy requirements). The periodicity has to be
set to 111 still, ELC cancels the electrostatic contribution of the
periodic replica in z-direction. Make sure that you read the papers on ELC
([AdeJoannisH02][TAH08]) before using it. ELC is an ESPResSo actor and is used
with:

The gap size gives the height of the empty region between the system box
and the neighboring artificial images. ESPResSo does not
make sure that the gap is actually empty, this is the users
responsibility. The method will compute fine if the condition is not
fulfilled, however, the error bound will not be reached. Therefore you
should really make sure that the gap region is empty (e.g. with wall
constraints).

maxPWerror:

The maximal pairwise error sets the least upper bound (LUB) error of
the force between any two charges without prefactors (see the papers).
The algorithm tries to find parameters to meet this LUB requirements or
will throw an error if there are none.

delta_mid_top/delta_mid_bot:

ELC can also be used to simulate 2D periodic systems with image charges,
specified by dielectric contrasts on the non-periodic boundaries
([TAH08]). Similar to MMM2D, these can be set with the
keywords delta_mid_bot and delta_mid_top, setting the dielectric
jump from the simulation region (middle) to bottom (at z<0) and
from middle to top (z>box_l[2]-gap_size). The fully metallic case
delta_mid_top=delta_mid_bot=-1 would lead to divergence of the
forces/energies in ELC and is therefore only possible with the
const_pot option.

const_pot:

As described, setting this to 1 leads to fully metallic boundaries and
behaves just like the mmm2d parameter of the same name: It maintains a
constant potential pot_diff by countering the total dipole moment of
the system and adding a homogeneous electric field according to
pot_diff.

pot_diff:

Used in conjunction with const_pot set to 1, this sets the potential difference
between the boundaries in the z-direction between z=0 and
z=box_l[2]-gap_size.

far_cut:

The setting of the far cutoff is only intended for testing and allows to
directly set the cutoff. In this case, the maximal pairwise error is
ignored.

neutralize:

By default, ELC just as P3M adds a homogeneous neutralizing background
to the system in case of a net charge. However, unlike in three dimensions,
this background adds a parabolic potential across the
slab [BAC09]. Therefore, under normal circumstance, you will
probably want to disable the neutralization for non-neutral systems.
This corresponds then to a formal regularization of the forces and
energies [BAC09]. Also, if you add neutralizing walls
explicitly as constraints, you have to disable the neutralization.
When using a dielectric contrast or full metallic walls
(delta_mid_top!=0 or delta_mid_bot!=0 or
const_pot=1), neutralize is overwritten and switched off internally.
Note that the special case of non-neutral systems with a non-metallic dielectric jump (eg.
delta_mid_top or delta_mid_bot in ]-1,1[) is not covered by the
algorithm and will throw an error.

where the prefactor \(C\) is defined in Eqn. (1).
MMM1D Coulomb method for systems with periodicity 0 0 1. Needs the
nsquared cell system (see section Cellsystems). The first form sets parameters
manually. The switch radius determines at which xy-distance the force
calculation switches from the near to the far formula. The Bessel cutoff
does not need to be specified as it is automatically determined from the
particle distances and maximal pairwise error. The second tuning form
just takes the maximal pairwise error and tries out a lot of switching
radii to find out the fastest one. If this takes too long, you can
change the value of the setmd variable timings, which controls the number of
test force calculations.

MMM1D is also available in a GPU implementation. Unlike its CPU
counterpart, it does not need the nsquared cell system. The first form
sets parameters manually. The switch radius determines at which
xy-distance the force calculation switches from the near to the far
formula. If the Bessel cutoff is not explicitly given, it is determined
from the maximal pairwise error, otherwise this error only counts for
the near formula. The second tuning form just takes the maximal pairwise
error and tries out a lot of switching radii to find out the fastest
one.

Espresso can use the electrostatics methods from the SCAFACOS Scalable
fast Coulomb solvers library. The specific methods available depend on the compile-time options of the library, and can be queried using espressomd.scafacos.available_methods

To use SCAFACOS, create an instance of espressomd.electrostatics.Scafacos and add it to the list of active actors. Three parameters have to be specified:

The method-specific parameters are described in the SCAFACOS manual.
Additionally, methods supporting tuning have the parameter tolerance_field which sets the desired root mean square accuracy for the electric field

To use the, e.g., ewald solver from SCAFACOS as electrostatics solver for your system, set its
cutoff to \(1.5\) and tune the other parameters for an accuracy of
\(10^{-3}\), use:

For details of the various methods and their parameters please refer to
the SCAFACOS manual. To use this feature, SCAFACOS has to be built as a shared library. SCAFACOS can be used only once, either for Coulomb or for dipolar interactions.