Constant potential simulations on a mesh Open Access

Molecular dynamics (MD) simulations are a popular method to investigate electrochemical processes of ionic liquids,^1–4 supercapacitors,^5–7 water in salt electrolytes,^8,9 and batteries^10,11 at the fundamental atomic level. For such cases, the constant potential method (CPM)^1,12 has become a valuable tool in simulations involving an electrolyte between metallic electrodes. In contrast to the rather simple approach of using a homogeneous and constant charge on all electrode atoms, the CPM adjusts the charges of the electrode atoms as a function of a constant electric potential applied between the electrodes. This method has been used to calculate the differential capacitance of various capacitor designs containing aqueous electrolytes or mixtures of ionic liquids and solvents as electrolytes.^3,5,13,14 Moreover, the mechanisms of electric double layer formation in supercapacitors have been extensively studied in the past, leading to different ion density profiles and realistic formation timescales compared to a fixed charge model for both flat and nanoporous electrodes.^7,15–17 Simulations without an electrolyte furthermore illustrated an influence of the electrode surface structure even for empty capacitors.³ Other applications include charge transfers at an electrochemical interface¹⁸ and the hydration of metal surfaces.¹⁹ Recently, even a Thomas–Fermi model has been incorporated to account for imperfect shielding of metals.²⁰ Though, this approach necessitates electrodes with rather large thicknesses comparable to the characteristic screening length of the specific electrode material.

A major shortcoming is that CPMs are only feasible for relatively small systems due to their higher computational cost—particularly when used in a different context in ab initio approaches.^21,22 Coarse-grained models are a useful approach to circumvent this and are commonly employed.² However, if the molecules have not been parameterized, an all-atomic molecular representation is required. Apart from this, alternative approaches using image charges^10,11 or potentiostats with the electric potential as a thermodynamic degree of freedom²³ have been developed. Furthermore, the finite field method allows the use of a CPM with periodic boundary conditions in three dimensions for slab-shaped geometries, requiring no artificial vacuum.^24,25 Until now, these approaches only allow planar electrodes or a homogeneous distribution of charge.

In this work, a more efficient algorithm based on a particle–particle particle–mesh (P³M) method is presented to compute the Coulomb interactions within a CPM on considerably smaller timescales as an alternative to the commonly used Ewald summation methods. The resulting expressions and their differences to previous implementations^15,26 of particle–mesh algorithms in combination with the CPM are discussed in detail. Upon application to moderate system sizes, as introduced in more detail later in this work, a noticeable speed-up is observed by using P³M in combination with a CPM.

Although there are already several CPM implementations available, e.g., in the MetalWalls code²⁷ or as an add-on for the molecular dynamics code Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)²⁸ by Wang et al.,¹⁷ we present here a new implementation of the CPM for LAMMPS that is conceptually more in line with the fundamental programming idea of LAMMPS by taking full advantage of its modularity. In order to describe the Coulomb interactions in capacitors periodic in two dimensions, two forms of the Ewald summation are commonly used: a summation in three dimensions with a dipole correction (EW3Dc)^29,30 or a summation explicitly over two dimensions (EW2D).^31,32 This implementation includes both EW2D and EW3Dc as well as an extension with the P³M method. Furthermore, if the EW3Dc approach is used within the CPM, close attention should be paid to the charge neutrality, as will be demonstrated for a simple model capacitor containing an aqueous electrolyte between two metal electrodes.

Molecular dynamics simulations with charged particles require the calculation of charge interactions, which are often idealized by simple Coulomb interactions between point charges. Consider a system of N point charges q_i at positions r_i which satisfies the charge neutrality condition $∑iNqi=0$ in a periodically repeated rectangular simulation cell with lengths in the x-, y-, and z-direction of L_x, L_y, and L_z, respectively. The total Coulomb energy is then given by the following expression:

Due to the slowly decaying 1/r term, the interactions between all particles and their periodic images must be taken into account for an accurate result. The first sum runs over the periodic images of the system, with the prime indicating that we skip the i = j interaction for n = 0. Note that we have omitted the 1/(4πε₀) factor for simplicity. Calculating this interaction poses a major challenge due to its long-range character, and to tackle this issue, the Ewald summation is frequently employed. In this formalism, the Coulomb interaction is divided into two sums, commonly denoted as a short- and a long-range part, that are calculated in real-space and (Fourier transformed) k-space, respectively. However, splitting the Coulomb sum into two parts requires a special handling of the self-interaction term originating from the long-range treatment of compensating Gaussian charge distributions.

The short-range term using purely point charges (denoted by qq in the subscript) is defined as

and is cut off by a complementary error function with the convergence parameter α. The long-range interaction

sums over the reciprocal space vector k^T = 2π(u_x/L_x, u_y/L_y, u_z/L_z). u_x, u_y, and u_z are integers and V denotes the volume of the simulation box. The k → 0 part of the long-range interaction is treated separately as follows:

Furthermore, there exists an infinite boundary term (which arises from the regularization of the lattice sum) that is related to the square of the total dipole moment of the simulation cell (see, e.g., Refs. 32 and 33). In contrast to the usage of point charges for the electrolyte, we consider Gaussian shaped charge densities for the electrode atoms with a reciprocal width η_i centered at electrode atom positions R_i,

Interactions with Gaussian shaped charge densities modify the short-range term³⁴ 

with $ηij=ηiηj/ηi2+ηj2$ for two interacting Gaussian charges (denoted by the subscript “QQ”) and η_ij = η_i if a Gaussian charge interacts with a point charge (later denoted by the subscript “Qq”). It is noteworthy that the $ηiQi2/2π$ term does not compensate for the self-interaction between real-space Gaussian charges in k-space (as one might anticipate from the long-range term), but quite to the opposite it is the actual self-interaction of Gaussian charges in real-space. In fact, the interaction in reciprocal space is the same as for a system of point charges in Eqs. (3) and (4) and independent of η_i.³⁴ 

A constant potential applied between two infinite plates requires an estimation of the electrostatic potential at each electrode atom and a treatment of the long-range electrostatic interactions of all atoms in only two periodic dimensions. Here, we discuss two approaches to tackle this challenge. First, a commonly used dipole correction is added to the EW3D-method yielding the well-known EW3Dc method.³⁰ In a second approach, the energy is derived for a system that is explicitly non-periodic in one direction leading to the analytically correct EW2D method.³² As the order of summation is important in Eq. (3), Smith²⁹ and more recently Hu³² showed that a slab-like summation geometry, e.g., the xy plane first then extending over z, results in a net dipole for the infinite boundary term

with $Mz=∑iNqizi$ denoting the z-component of the total dipole moment. As shown by Yeh and Berkowitz,³⁰ adding $Uk→0EW3Dc$ from Eq. (7) as the infinite boundary term in a regular EW3D summation using otherwise so-called “tin-foil” conditions (i.e., $Uk→0EW3D=0$⁠) approximates an infinite slab reasonably well if sufficient empty space is introduced in the z-direction.

In the EW2D method, on the contrary, the Fourier summation is explicitly evaluated only for the periodic dimensions (a more detailed derivation is found in Ref. 32). This two-dimensional summation results in

for the infinite boundary term of an infinitely extended slab. The short-range term in the EW2D formalism is the same as in Eq. (6) but with n_z = 0 in the sum over |n| because the system is not periodic in this direction.

In contrast to a constant charge approach, an interface between the electrode and electrolyte leads to non-uniform and dynamic charge distributions in the electrode.^13,17 Initially developed by Siepmann and Sprik,¹² the constant potential method (CPM) was later modified by Reed et al.¹ to adapt the electrode charges Q_i at every time step in an extended Hamiltonian approach and was further improved by Gingrich and Wilson³⁴ who presented a rigorous derivation for the long-range interaction of Gaussian charge distributions. Assuming perfect conductivity, the charge Q_i of each electrode atom is obtained by minimizing the energy U({Q_i}) − ∑_iψ_iQ_i for an electric potential ψ_i applied on each electrode atom i. We now deviate a bit from the established notation of matrices and vectors and distinguish between charges and positions of electrode and electrolyte atoms by using upper and lower case letters, respectively. For the CPM approach used here, the electric potential is cast into a set of linear equations

with b containing all terms involving q_n, i.e., the electrolyte charges.^3,17 Given an external potential and a set of particle positions, the equation can then be solved for Q_m, e.g., by inverting the A-matrix,

or alternatively with a conjugate gradient (CG) method. Conveniently, the A-matrix is independent of the charges and electrolyte positions. Therefore, in principle, the inversion has to be carried out only once, provided the electrode positions are frozen in place. The derivative of the electrostatic energy [cf. long- and short-range terms in Eqs. (3) and (6)] with respect to Q_i is the electric potential,

Here, the sum of all terms containing electrolyte charges q_n equals the term b_i in Eq. (9) while the rest can be written as ∑_mA_imQ_m. Unfortunately, each element of A and b in Eq. (10) requires a summation over all k-vectors in the reciprocal space. Hence, calculating the energy from A and b would result in a large computational effort compared to the long-range energy in Eq. (3), which features only one sum over the reciprocal space. However, this step is necessary to satisfy the linear form of Eq. (9). One can think of the sums over the reciprocal space and electrolyte atoms in the long-range term as discrete Fourier transforms, which are the basis for the improvements in Sec. II D.

The charges obtained from Eq. (10) do not necessarily fulfill the charge neutrality constraint ∑_iQ_i = 0. In order to enforce this constraint, one can replace A⁻¹ with the symmetric matrix

with e^T = (1, ⋯, 1).³ In principle, one could explore the possibility of implementing the charge neutrality condition via the S-matrix formalism using a CG method. Such an approach, however, is far from trivial, and we do not pursue it here and instead solve Eq. (9) via matrix inversion and replace A⁻¹ with S unless noted otherwise.

In the Ewald formalism, the computational effort for the long-range terms of the A-matrix and b-vector in Eq. (10) quickly increases with the system size because for each element a summation over all k-vectors is required [cf. sum over k ≠ 0 in Eq. (11)]. The P³M method is a trade-off to the Ewald method in terms of speed, but at the expense of accuracy, which is facilitated by a Fourier transform into the k-space and mapping charges onto a regular mesh.^35,36 However, the error introduced by projecting charges onto the mesh is controllable by choosing an appropriate mesh and stencil size.

First, we will only consider electrode charges and extend this formalism later for electrode–electrolyte interactions by deriving an expression for the b-vector mentioned in Eq. (9). A weight function W(r_μ − R_i) is used to map the charges with position R_i onto the N_mesh points r_μ of the mesh, resulting in the charge density

A discrete Fourier transform of the charge density onto a mesh in reciprocal space gives

Thanks to the equal spacing of the mesh, the fast Fourier transform (FFT) algorithm can be used.³⁷ The long-range energy of electrode–electrode interactions

is then calculated in the k-space, similar to the energy in the Ewald formalism in Eq. (3). $Ĝ(k)$ is the lattice Green’s function, for which several optimizations have been developed in order to closely match the results of the Ewald method.^38,39 We use a lattice Green’s function that is optimized for the ik-differentiation scheme, for which a derivative is obtained via multiplication with ik in reciprocal space yielding the electric field. The electric field is then used to calculate the forces acting on the charged particles.³⁸ This method and the corresponding Green’s function are readily available in LAMMPS and give the most accurate forces.⁴⁰ In the following, the focus is on the specific implementation of a CPM via P³M; a more general description of P³M can be found in Ref. 40. By inserting Eqs. (13) and (14) into Eq. (15), we arrive at an expression featuring a sum over charge pairs as follows, which directly facilitates an efficient evaluation of the long-range part of the A-matrix required for the CPM:

We obtain an expression with a sum over the k-space for each pair of mesh points, similar to the sum for each charge pair in the Ewald formalism. However, in contrast to the A-matrix in the Ewald formalism, it is straightforward to calculate the A-matrix on the mesh,

exploiting the FFT algorithm, thanks to the equal spacing of distances r_μν on the grid. To project the A-matrix from the mesh back to the particle positions, a multiplication with the weight matrix W_μi = W(r_μ − R_i) is necessary,

For the projection between mesh and particle positions, a small subset of nearby mesh points is sufficient for each particle because the weighting function vanishes for large distances. Thus, the sums over μ and ν do not run over all mesh points but only those in the vicinity of particles i and j, respectively. In LAMMPS, this subset of the mesh, for which the weighting function is non-zero, is implemented as a cube-shaped stencil of mesh points with an adjustable edge length. To further reduce the computational cost, one can calculate the matrix in two steps by introducing an intermediate matrix

This results in a linear rather than a quadratic scaling of the required calculations with the stencil size. However, in terms of memory efficiency, this scales with the mesh size because X_μj has to be calculated for all μ. This means that the intermediate matrix can become prohibitively large as it has a size of N_mesh × N_electrode. In the following, the intermediate matrix is used because of its lower computational effort for the investigated systems. The A-matrix resulting from Eq. (18) consequently replaces the traditional long-range part of the CPM and is used along with the short-range and the other correction terms required in Eq. (10). Only the long-range term is hence approximated with P³M, and all other terms, including the k → 0 correction term in Eq. (7), are unchanged and similar to a regular three-dimensional P³M with the dipole correction. Consequently, we denote the usage of P³M in the CPM here as P³M3Dc. The dipole correction from EW3Dc is much more convenient, readily implemented in LAMMPS, and (under most circumstances) computationally less expensive than the correction from the EW2D method.³⁰ 

A similar approach is used for the electrostatic potential of the electrode–electrolyte interaction (i.e., the b-vector). The long-range part of the electrolyte–electrode interactions in P³M is given by

The Fourier transform of the electrolyte charge density $ρ̂q(k)$ is defined akin to the electrode density in Eqs. (13) and (14). Again, Eqs. (13) and (14) are inserted for the charge densities of the electrodes,

Following the spirit of P³M, the electrode particles interact with a quantity on the mesh, namely, the Fourier transform of $Ĝ(k)ρ̂q*(k)$⁠. In fact, a very similar pattern is used to obtain forces from the interaction of particles with the electric field in the ik-differentiation scheme. Hence, calculating $blrP3M$ is fast, as it merely requires an FFT followed by a sum over the stencil of the respective electrode atom for each element of the vector.

It is important to realize that an efficient algorithm for the b-vector is essential, as it needs to be evaluated for each charge update that is performed after a fixed interval, typically the molecular dynamics time step. While particle–mesh algorithms have been used for this approach before, they required a sum over electrolyte charges for each element of b_lr.^15,26 In contrast to this, Eq. (21) necessitates just a sum over the stencil, leading to a better scaling with respect to the number of electrolyte particles. Adapting the expression in Eq. (21) for other methods that use a CPM could be an interesting application in the future. In particular, for the finite field method, which is proposed as another alternative to implement a CPM,^24,25 the adaptation should be straightforward since the above expressions are also valid for a system with periodic boundary conditions in three dimensions.

In contrast to the EW2D method in MetalWalls,²⁷ Wang et al.¹⁷ have implemented a CPM using the EW3Dc approach. Here, we present a new implementation of CPM in LAMMPS featuring both Ewald methods as well as P³M3Dc. To the best of our knowledge, this is the first implementation of the CPM using the P³M3Dc method in the form presented above. Our implementation is based on the KSPACE-package available in LAMMPS, which provides both Ewald and P³M methods for regular MD simulations. As a model system, we choose a gold capacitor with ten atomic layers in each electrode, adapted to the LAMMPS format from an example in the MetalWalls repository.⁴¹ 

We start by discussing the use of the dipole correction in the CPM and investigate the effect of the associated slab factor as well as the impact of the matrix symmetrization [cf. Eq. (12)] on the charge density of a vacuum capacitor. Our model system consisted of an empty gold capacitor with 6480 atoms per electrode (cf. inset in Fig. 1), to which a voltage of Δψ = 2 V was applied between the electrodes and which were separated by a vacuum of 100 Å. The bottom of the box is set to z = 0, i.e., the system is not centrosymmetric along the z-axis.

Additional vacuum is required in the z-direction between periodically repeated slabs in order to calculate appropriately the long-range terms for slab-like systems using the EW3Dc formalism and is added via the usual slab factor V/V₀. An increased slab factor gives rise to increased computational efforts due to a larger number of required k-points. Nevertheless, using sufficient vacuum is necessary to combat artifacts that stem from carrying out the summation in the k-space in three dimensions for a system only periodic in two. Looking more closely at Fig. 1, we observe a rather strong dependence of the resulting surface charge densities σ on the electrodes on the slab factor if the S-matrix is not used. This is indeed expected but quite strong and surprising for an asymmetric simulation setup like this. However, by enforcing charge neutrality using the S-matrix as defined in Eq. (12), charge densities converge almost immediately to the reference values using EW2D (crosses in Fig. 1) even for very small slab factors and agree well with only a difference of about 0.2% between both methods. Since the electrodes are very homogeneous in the x- and y-direction, a small distance between periodic images is sufficient and the dipole in the z-direction is appropriately corrected using Eq. (7).

In an attempt to explain the large deviations from calculations without enforced charge neutrality, one can look at the derivative of the correction term

employed to obtain the electric potential of a charge neutral system in the EW3Dc formalism. According to this, a system shifted by Δz along the z-axis has an offset in the derivative $∂Uk→0EW3Dc/∂Qj$ of 4π/VΔzM_z, while the energy term $Uk=0EW3Dc$⁠, which is sufficient in regular MD simulations that do not make use of the CPM, is invariant to this shift. The offset decreases with increasing gap introduced between the periodic images of the non-periodic dimension because the volume increases while Δz remains constant. As demonstrated in Fig. 1, this shift leads to a net charge of the system that decreases for large vacuum gaps between the periodically repeated images of the slabs. It is noteworthy that Smith²⁹ already pointed out that a suitable choice of coordinates is necessary for using Eq. (7), but they do not discuss the implications in too much detail, particularly not with respect to a CPM. Fortunately, by implementing the charge neutrality condition as shown by Scalfi et al.³ in our algorithm, only the difference between the two potentials is considered in the simulations. Hence, the obtained charge densities sum up to zero and resemble the values obtained with the EW2D approach. As an alternative, one could try to find an expression for the correction in Eq. (22) that is invariant with respect to shifts along the z-axis.

First, we compare the convergence between the P³M and a regular Ewald method for a vacuum capacitor, which is again composed of 6480 gold atoms in each electrode with an electrode separation of 100 Å (cf. the inset in Fig. 2). The A-matrix is generally independent of the electrolyte, and in the absence of an electrolyte, the charge on the electrodes is solely determined by the A-matrix. In this case of a static calculation, obtaining and inverting the A-matrix make up the main part of the calculations. Figure 2 illustrates the charge densities obtained using P³M3Dc with varying convergence parameters and compares them to the results obtained using a traditional EW3Dc that is used as reference here. Resulting calculation times for the long-range term of A are given on the x-axis in the logarithmic scale. It should be mentioned that the time for the matrix inversion is not included because it was mostly unaffected and about 21 min in all our calculations. For both methods, the number of k-points was set proportional to the box length in each direction, taking also the added vacuum in the non-periodic z-direction into account. Charge densities converge quite quickly by employing the traditional EW3Dc approach as the size of the reciprocal space is increased. However, the calculations required for this become prohibitively expensive and lead to quite long overall calculation times.

Using our new P³M implementation for the CPM in combination with the intermediate matrix of Eq. (19) results in orders of magnitude shorter wall times as shown in Fig. 2. In the P³M implementation of LAMMPS, each processor is assigned a subsection of the mesh, which allows an efficient distribution of the intermediate matrix across processors. Thus, the memory required for the intermediate matrix was manageable for reasonable mesh sizes. In the frame of this work, we investigated the behavior with an increasing stencil width ranging from 2 to 5 for different mesh sizes. For a stencil width of 5 and meshes with sizes (24, 24, 101) and (32, 32, 135), relative differences below 10⁻⁴ are obtained for the charge densities compared to the converged value from the EW3Dc method. This results in a substantial reduction in computation time by more than an order of magnitude.

Next, we also apply the P³M method to calculate the b-vector. For this, a capacitor with an aqueous NaCl electrolyte confined by two gold electrodes is used (cf. the inset of Fig. 3). The area of the electrodes is one-fourth that of the capacitor in the previous part, as this smaller system size is sufficient to illustrate the advantages of the P³M implementation. Here, the distance to a reference b_ref from an EW3Dc calculation with a large k-space is used to assess the accuracy of our P³M implementation. Resulting long-range computation times of each setting with respect to the accuracy (in relation to the reference) are indicated by points in Fig. 2 (the respective number of k-points used for each point in Fig. 2 are found in Table S1 in the supplementary material). Figure 2 shows that P³M3Dc converges quickly to the reference as the k-space is increased with a fixed edge length of 5 for the stencil. The mesh size for P³M3Dc ranges from (6, 6, 52) to (32, 32, 280) and the size of the reciprocal space for EW3Dc from (3, 3, 26) to (12, 12, 105). Even though larger reciprocal spaces are needed for P³M3Dc than for EW3Dc to obtain results with a comparable accuracy, significantly lower computation times are required. The computation times for the matrix multiplication [cf. Eq. (10)] as well as the short-range terms and the dipole correction are independent of the long-range solver and are about 10 ms for all calculations (cf. Table S2 in the supplementary material). Of these contributions, the matrix multiplication is the most expensive part and while this specific implementation could be further optimized, its scaling with the square of the electrode size will always lead to a large contribution for large systems. While the long-range term requires over 80% of the time in the CPM routine for the more expensive Ewald settings, this value drops to 30% with P³M at a similar accuracy. Thus, the achieved speed-up of the long-range term has a significant impact on the total run time. Thanks to the improved scaling, one can expect the difference to be even more notable for larger systems, and naturally, this aspect will be investigated in our future works.

Structural variability in the electrolyte at the atomic level—rarely if ever taken into account in standard models of double layer theory, such as the Gouy–Chapman model—has a strong influence on the differential capacitance.⁴² Thus, the differential capacitance allows comparisons of the dynamics with different CPM implementations; exploring the influence of the explicit electric double layer structure on the differential capacitance with atomistic simulations is moreover instrumental to a better understanding of electrochemical phenomena at the solid–liquid interface. Recently, Scalfi et al.³ have derived a rigorous expression for the differential capacitance C_diff as a function of the electrode charge fluctuations using an NVT ensemble in combination with a CPM. It was found that the capacitance comprises two parts: a contribution from an empty capacitor, which is solely a function of the A-matrix, and an electrolyte contribution. The latter can be calculated from the electrode charge fluctuations evolving during the simulation at a constant potential.¹³ 

For a rigorous comparison of different implementations of the CPM and the various potentially useful long-range approaches, we used the same gold capacitor with an aqueous NaCl electrolyte as previously used by Scalfi et al.³ Figure 3 summarizes both contributions to the differential capacitance and the standard error resulting in the electrolytic contributions stemming from the various long-range methods implemented in a CPM. MetalWalls²⁷ results are given along with those obtained with our CPM implementation in LAMMPS, denoted as USER-CONP, using the various Ewald methods mentioned above, which are also indicated. Charge fluctuations for calculating the electrolytic contribution have been obtained after a 2 ns equilibration in an NVT simulation at 298 K comprising a simulation time of 3 ns using a time step of 1 fs and at a potential difference of 2 V. While the electrolytic capacitance contributions do show some negligible differences due to stochastic effects and limited simulations length, vacuum capacitances only depend on the entries in the A-matrix and agree almost perfectly. Evaluating the vacuum capacitances from the matrices of the various approaches as given in Fig. 3 resulted in relative differences between the vacuum capacitances that are smaller than 0.1%.

We demonstrated that the developed particle–particle particle–mesh (P³M) method facilitates the incorporation of long-range electrostatic interactions based on a constant potential approach at considerably smaller timescales in comparison to the previously employed techniques. In particular, and in contrast to describing the electrode–electrolyte interactions by a homogeneous charge on the electrodes in so-called constant charge approaches, a more physically realistic approach, such as a constant potential based on a dipole corrected three-dimensional Ewald (EW3Dc) method for the long-range electrostatic computations, can achieve significantly more efficient computations when P³M is used for the long-range part already for moderate system sizes with about 10⁴ particles. This was achieved by efficiently calculating not only the A-matrix—originating from the linear equations used to determine the electrode charges—on a regular mesh, but also compute the b-vector on the same mesh, which includes all electrostatic interactions between electrodes and electrolytes.

In addition, we demonstrated that when using a dipole correction in a constant potential simulation, which is usually required for the EW3Dc method, a charge neutrality condition in the evaluation of the electrostatic potential must be enforced. Otherwise, an artificial potential shift will be introduced if the simulation is not perfectly symmetric even if an additional vacuum (as given by the slab factor V/V₀) commonly found in the literature is incorporated. Apparently, even doubling the introduced vacuum volume by twice the distance between the electrodes (which greatly increases the computational cost) cannot completely eliminate this artifact. Hence, we strongly encourage the use of the S-matrix in Eq. (12) as introduced in Scalfi et al.³ when dipole corrections in the various Ewald summation techniques are employed in constant potential simulations. In the future, the implementation of P³M in constant potential simulations may allow not only larger systems but also moving electrode atoms that require frequent matrix updates or (not yet implemented) adapted conjugate gradient methods.

See the supplementary material for a discussion on the scaling for the long-range term of the b-vector as well as detailed computation times for the b-vector calculations discussed in Sec. III B.

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) (Projektnummer 192346071—SFB 986 and Projektnummer 390794421—GRK 2462).

The source code and the here investigated examples are made available at https://github.com/USER-CONP/lammps/tree/v1.1 and are intended to be incorporated in the main branch of the LAMMPS distribution.