9. Electrostatics¶
The Coulomb (or electrostatic) interaction is defined as
follows. For a pair of particles at distance
where
is a prefactor which can be set by the user. The commonly used Bjerrum length
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 attached to the system object to become active. Only one electrostatics method can be active at any time.
Note that using the electrostatic interaction also requires assigning charges to
the particles via the particle property
q
.
All solvers need a prefactor and a set of other required parameters. This example shows the general usage of the electrostatic method P3M (often stylized as P3M). An instance of the solver is created and attached to the system, at which point it will be automatically activated. This activation will internally call a tuning function to achieve the requested accuracy:
import espressomd
import espressomd.electrostatics
system = espressomd.System(box_l=[10, 10, 10])
system.time_step = 0.01
system.part.add(pos=[[0, 0, 0], [1, 1, 1]], q=[-1, 1])
solver = espressomd.electrostatics.P3M(prefactor=2., accuracy=1e-3)
system.electrostatics.solver = solver
where the prefactor is defined as
The solver can be detached with either:
system.electrostatics.solver = None
or:
system.electrostatics.clear()
9.1. Coulomb P3M¶
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 (True, True, True)
. When using a non-metallic dielectric
constant (epsilon != 0.0
), the box must be cubic.
Make sure that you know the relevance of the P3M parameters before using P3M!
If you are not sure, read the following references:
[Cerdà et al., 2008, Deserno, 2000, Deserno and Holm, 1998, Deserno and Holm, 1998, Deserno et al., 2000, Ewald, 1921, Hockney and Eastwood, 1988, Kolafa and Perram, 1992].
9.1.1. Tuning Coulomb P3M¶
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.
The tuning method is called when the handle of the Coulomb P3M is added to
the actor list. Some parameters can be fixed (r_cut
, cao
, mesh
)
to speed up the tuning if the parameters are already known.
Please 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 [Hockney and Eastwood, 1988] and its real space error [Kolafa and Perram, 1992] 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.
During tuning, the algorithm reports 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.
9.1.2. Coulomb P3M on GPU¶
espressomd.electrostatics.P3MGPU
The GPU implementation of P3M calculates the far field contribution to the forces on the GPU. The near-field contribution to the forces, as well as the near- and far-field contributions to the energies are calculated on the CPU. 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 extension Dielectric interfaces with the ICC* algorithm.
The algorithm doesn’t have kernels to compute energies and pressures and therefore uses the respective CPU kernels with the parameters tuned for the GPU force kernel.
9.2. Debye-Hückel potential¶
The Debye-Hückel electrostatic potential is defined by
where
For
9.3. Reaction Field method¶
espressomd.electrostatics.ReactionField
The Reaction Field electrostatic potential is defined by
where
with
The term in
9.4. Dielectric interfaces with the ICC algorithm¶
espressomd.electrostatic_extensions.ICC
The ICC
import espressomd.electrostatics
import espressomd.electrostatic_extensions
p3m = espressomd.electrostatics.P3M(...)
icc = espressomd.electrostatic_extensions.ICC(...)
system.electrostatics.solver = p3m
system.electrostatics.extension = icc
The ICC particles are setup as normal ESPResSo particles. Note that they should be fixed in space and need an initial non-zero charge. The following 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 size
box_l = 9.
system.box_l = [box_l, box_l, 12.]
nicc = 3 # linear density
nicc_per_electrode = nicc**2 # surface density
nicc_tot = 2 * nicc_per_electrode
iccArea = box_l**2 / nicc_per_electrode
l = box_l / nicc
# Lists to collect required parameters
iccNormals = []
iccAreas = []
iccSigmas = []
iccEpsilons = []
# Add the fixed ICC particles:
# Left electrode (normal [0, 0, 1])
for xi in range(nicc):
for yi in range(nicc):
system.part.add(pos=[l * xi, l * yi, 0.], q=-0.0001,
type=icc_type, fix=[True, True, True])
iccNormals.extend([0, 0, 1] * nicc_per_electrode)
# Right electrode (normal [0, 0, -1])
for xi in range(nicc):
for yi in range(nicc):
system.part.add(pos=[l * xi, l * yi, box_l], q=0.0001,
type=icc_type, fix=[True, True, True])
iccNormals.extend([0, 0, -1] * nicc_per_electrode)
# Common area, sigma and metallic epsilon
iccAreas.extend([iccArea] * nicc_tot)
iccSigmas.extend([0] * nicc_tot)
iccEpsilons.extend([100000] * nicc_tot)
icc = espressomd.electrostatic_extensions.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.electrostatics.extension = 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 [Arnold et al., 2013, Kesselheim et al., 2011, Tyagi et al., 2010] before using it.
9.5. Electrostatic Layer Correction (ELC)¶
ELC is an extension of the P3M electrostatics solver for explicit 2D periodic
systems. It can account for different dielectric jumps on both sides of the
non-periodic direction. In more detail, it is a special procedure that
converts a 3D electrostatic method to a 2D method in computational order N.
The periodicity has to be set to (True, True, True)
. ELC cancels the electrostatic
contribution of the periodic replica in
Usage notes:
The non-periodic direction is always the
-direction.The method relies on a slab of the simulation box perpendicular to the
-direction not to contain particles. The size in -direction of this slab is controlled by thegap_size
parameter. The user has to ensure that no particles enter this region by means of constraints or by fixing the particles’ z-coordinate. When particles enter the slab of the specified size, an error will be thrown.
ELC is an ESPResSo actor and is used with:
import espressomd.electrostatics
p3m = espressomd.electrostatics.P3M(prefactor=1, accuracy=1e-4)
elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3)
system.electrostatics.solver = elc
Although it is technically feasible to detach elc
from the system
and then to attach the p3m
object, it is not recommended because the P3M
parameters are mutated by ELC, e.g. the epsilon
is made metallic.
It is safer to instantiate a new P3M object instead of recycling one that
has been adapted by ELC.
ELC can also be used to simulate 2D periodic systems with image charges,
specified by dielectric contrasts on the non-periodic boundaries
([Tyagi et al., 2008]). This is achieved by setting the dielectric jump from the
simulation region (middle) to bottom (at
elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3,
delta_mid_top=0.9, delta_mid_bot=0.1)
The fully metallic case const_pot
option.
Toggle const_pot
on to maintain a constant electric potential difference
pot_diff
between the xy-planes at
elc = espressomd.electrostatics.ELC(actor=p3m, gap_size=box_l * 0.2, maxPWerror=1e-3,
const_pot=True, delta_mid_bot=100.0)
This is done by countering the total dipole moment of the system with the
electric field pot_diff
.
9.6. MMM1D¶
espressomd.electrostatics.MMM1D
Please cite [Arnold and Holm, 2005] when using MMM1D. See MMM1D theory for the details.
MMM1D is used with:
import espressomd.electrostatics
mmm1d = espressomd.electrostatics.MMM1D(prefactor=C, far_switch_radius=fr,
maxPWerror=err, tune=False, bessel_cutoff=bc)
mmm1d = espressomd.electrostatics.MMM1D(prefactor=C, maxPWerror=err)
where the prefactor timings
argument of the
MMM1D
class,
which controls the number of test force calculations.
9.7. ScaFaCoS electrostatics¶
espressomd.electrostatics.Scafacos
ESPResSo can use the 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 with
espressomd.electrostatics.Scafacos.get_available_methods()
.
To use ScaFaCoS, create an instance of Scafacos
and attach it to the system. Three parameters have to be specified:
prefactor
(as defined in (1)), method_name
,
method_params
. The method-specific parameters are described in the
ScaFaCoS manual. In addition, methods supporting tuning have a parameter
tolerance_field
which sets the desired root mean square accuracy for
the electric field.
To use a specific electrostatics solver from ScaFaCoS for your system,
e.g. ewald
, set its cutoff to
import espressomd.electrostatics
scafacos = espressomd.electrostatics.Scafacos(
prefactor=1, method_name="ewald",
method_params={"ewald_r_cut": 1.5, "tolerance_field": 1e-3})
system.electrostatics.solver = scafacos
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.