Molecular dynamics simulations in a constant potential ensemble are an increasingly important tool to investigate charging mechanisms in next-generation energy storage devices. We present a highly efficient approach to compute electrostatic interactions in simulations employing a constant potential method (CPM) by introducing a particle–particle particle–mesh solver specifically designed for treating long-range interactions in a CPM. Moreover, we present evidence that a dipole correction term—commonly used in simulations with a slab-like geometry—must be used with caution if it is also to be used within a CPM. It is demonstrated that artifacts arising from the usage of the dipole correction term can be circumvented by enforcing a charge neutrality condition in the evaluation of the electrode charges at a given external potential.

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 batteries10,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.3 Other applications include charge transfers at an electrochemical interface18 and the hydration of metal surfaces.19 Recently, even a Thomas–Fermi model has been incorporated to account for imperfect shielding of metals.20 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.2 However, if the molecules have not been parameterized, an all-atomic molecular representation is required. Apart from this, alternative approaches using image charges10,11 or potentiostats with the electric potential as a thermodynamic degree of freedom23 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 (P3M) 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 implementations15,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 P3M in combination with a CPM.

Although there are already several CPM implementations available, e.g., in the MetalWalls code27 or as an add-on for the molecular dynamics code Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)28 by Wang et al.,17 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 P3M 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 qi at positions ri 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 Lx, Ly, and Lz, respectively. The total Coulomb energy is then given by the following expression:

U=12|n|i,j=1Nqiqj|rij+n|.
(1)

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πε0) 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

Usr,qqEW3D=12|n|i,j=1Nqiqjerfcα|rij+n||rij+n|απi=1Nqi2
(2)

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

UlrEW3D=12Vk04πk2ek24α2i=1Nqieikri2
(3)

sums over the reciprocal space vector kT = 2π(ux/Lx, uy/Ly, uz/Lz). ux, uy, and uz 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:

Uk0EW3D=12Vlimk04πk2ek24α2i=1Nqieikri2.
(4)

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 Ri,

χi(r)=Qiηi2π3/2expηi2rRi2.
(5)

Interactions with Gaussian shaped charge densities modify the short-range term34 

Usr,QQEW3D=i=1Nηi2παπQi2+12|n|i,j=1NQiQjerfcα|Rij+n||Rij+n|erfcηij|Rij+n||Rij+n|,
(6)

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.34 

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.30 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.32 As the order of summation is important in Eq. (3), Smith29 and more recently Hu32 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

Uk0EW3Dc=2πVMz2,
(7)

with Mz=iNqizi denoting the z-component of the total dipole moment. As shown by Yeh and Berkowitz,30 adding Uk0EW3Dc from Eq. (7) as the infinite boundary term in a regular EW3D summation using otherwise so-called “tin-foil” conditions (i.e., Uk0EW3D=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

Uk0EW2D=πLxLyi,jqiQjexpzij2α2α+πzijerfαzij,
(8)

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 nz = 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,12 the constant potential method (CPM) was later modified by Reed et al.1 to adapt the electrode charges Qi at every time step in an extended Hamiltonian approach and was further improved by Gingrich and Wilson34 who presented a rigorous derivation for the long-range interaction of Gaussian charge distributions. Assuming perfect conductivity, the charge Qi of each electrode atom is obtained by minimizing the energy U({Qi}) − iψiQi 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

ψi=mAimQm+bi,
(9)

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

Qm=iA1miψibi,
(10)

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 Qi is the electric potential,

UQi=mQmerfcαRmierfcηi2RmiRmi+nqnerfcαrnierfcηirnirni+1Vk04πk2expk24α2×mQmeikRmi+eikRinqneikrn+2ηi2απQi+Uk0Qi=ψi.
(11)

Here, the sum of all terms containing electrolyte charges qn equals the term bi in Eq. (9) while the rest can be written as mAimQm. 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 iQi = 0. In order to enforce this constraint, one can replace A−1 with the symmetric matrix

S=A1A1eeTA1eTA1e,
(12)

with eT = (1, ⋯, 1).3 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−1 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 P3M 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μRi) is used to map the charges with position Ri onto the Nmesh points rμ of the mesh, resulting in the charge density

ρQ(rμ)=1NmeshiQiW(rμRi).
(13)

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

ρ̂Q(k)=NmeshμρQ(rμ)expikrμ.
(14)

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

Ulr,QQP3M=12Vk0Ĝ(k)ρ̂Q(k)2,
(15)

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.38 This method and the corresponding Green’s function are readily available in LAMMPS and give the most accurate forces.40 In the following, the focus is on the specific implementation of a CPM via P3M; a more general description of P3M 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:

Ulr,QQP3M=12Vi,jQiQjμ,νW(rμRi)W(rνRj)×k0expik(rμrν)Ĝ(k)=12VQTWTAmeshWQ.
(16)

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,

Ameshμν=k0expikrμνĜ(k),
(17)

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μRi) is necessary,

AlrP3Mij=1Vμ,νWiμAmeshμνWνj.
(18)

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

Xμj=νAmeshμνWνj.
(19)

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 Nmesh × Nelectrode. 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 P3M, and all other terms, including the k → 0 correction term in Eq. (7), are unchanged and similar to a regular three-dimensional P3M with the dipole correction. Consequently, we denote the usage of P3M in the CPM here as P3M3Dc. 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.30 

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 P3M is given by

Ulr,QqP3M=1Vk0Ĝ(k)ρ̂Q(k)ρ̂q*(k).
(20)

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,

Ulr,QqP3M=1ViQiμW(rμRi)k0expikrμĜ(k)ρ̂q*(k)=QTblrP3M.
(21)

Following the spirit of P3M, 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 blr.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,27 Wang et al.17 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 P3M3Dc. To the best of our knowledge, this is the first implementation of the CPM using the P3M3Dc method in the form presented above. Our implementation is based on the KSPACE-package available in LAMMPS, which provides both Ewald and P3M 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.41 

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.

FIG. 1.

Surface charges for Δψ = 2 V using EW3Dc with or without enforcing charge neutrality (i.e., with S-matrix) plotted against the slab factor V/V0. The inset shows the simulation box for V/V0 = 2.

FIG. 1.

Surface charges for Δψ = 2 V using EW3Dc with or without enforcing charge neutrality (i.e., with S-matrix) plotted against the slab factor V/V0. The inset shows the simulation box for V/V0 = 2.

Close modal

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/V0. 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

Uk0EW3DcQj=4πVzjMz
(22)

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 Uk0EW3Dc/Qj of 4π/VΔzMz, 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 Smith29 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.3 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 P3M 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 P3M3Dc 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.

FIG. 2.

Convergence of the long-range part of A and b for different k-space settings. The employed model capacitor for Alr is shown in the inset (for blr cf. Fig. 3). k-points for EW3Dc are gradually increased in order to achieve convergence at the cost of increasing computation time. For Alr, with the P3M3Dc method (blue lines), both the number of k-points (specified next to lines) and the edge length of the stencil is varied between 2 and 5 (indicated by blue markers). For blr, the stencil edge length is fixed at 5 and meshes are increased gradually. Computation times for the b-vector are obtained from molecular dynamics and are given as an average per simulation step.

FIG. 2.

Convergence of the long-range part of A and b for different k-space settings. The employed model capacitor for Alr is shown in the inset (for blr cf. Fig. 3). k-points for EW3Dc are gradually increased in order to achieve convergence at the cost of increasing computation time. For Alr, with the P3M3Dc method (blue lines), both the number of k-points (specified next to lines) and the edge length of the stencil is varied between 2 and 5 (indicated by blue markers). For blr, the stencil edge length is fixed at 5 and meshes are increased gradually. Computation times for the b-vector are obtained from molecular dynamics and are given as an average per simulation step.

Close modal

Using our new P3M 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 P3M 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−4 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 P3M 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 P3M implementation. Here, the distance to a reference bref from an EW3Dc calculation with a large k-space is used to assess the accuracy of our P3M 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 P3M3Dc 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 P3M3Dc 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 P3M3Dc 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 P3M 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.

FIG. 3.

Comparison of differential capacitances obtained using different implementations and methods. The model system containing an aqueous NaCl solution between two gold electrodes is shown in the inset. The electrode atoms are colored according to the charge and electrolyte atoms according to the atom type. The vacuum contribution to the capacitance is calculated from the A-matrix (blue bars), the electrolyte contribution from electrode charge variances.

FIG. 3.

Comparison of differential capacitances obtained using different implementations and methods. The model system containing an aqueous NaCl solution between two gold electrodes is shown in the inset. The electrode atoms are colored according to the charge and electrolyte atoms according to the atom type. The vacuum contribution to the capacitance is calculated from the A-matrix (blue bars), the electrolyte contribution from electrode charge variances.

Close modal

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.42 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.3 have derived a rigorous expression for the differential capacitance Cdiff 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.13 

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.3Figure 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. MetalWalls27 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 (P3M) 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 P3M is used for the long-range part already for moderate system sizes with about 104 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/V0) 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.3 when dipole corrections in the various Ewald summation techniques are employed in constant potential simulations. In the future, the implementation of P3M 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.

1.
S. K.
Reed
,
O. J.
Lanning
, and
P. A.
Madden
,
J. Chem. Phys.
126
,
084704
(
2007
).
2.
C.
Merlet
,
M.
Salanne
, and
B.
Rotenberg
,
J. Phys. Chem. C
116
,
7687
(
2012
).
3.
L.
Scalfi
,
D. T.
Limmer
,
A.
Coretti
,
S.
Bonella
,
P. A.
Madden
,
M.
Salanne
, and
B.
Rotenberg
,
Phys. Chem. Chem. Phys.
22
,
10480
(
2020
); arXiv:1911.09351.
4.
D. J.
Eyckens
,
L.
Servinis
,
C.
Scheffler
,
E.
Wölfel
,
B.
Demir
,
T. R.
Walsh
, and
L. C.
Henderson
,
J. Mater. Chem. A
6
,
4504
(
2018
).
5.
C.
Merlet
,
B.
Rotenberg
,
P. A.
Madden
,
P.-L.
Taberna
,
P.
Simon
,
Y.
Gogotsi
, and
M.
Salanne
,
Nat. Mater.
11
,
306
(
2012
).
6.
S. M.
Mahurin
,
E.
Mamontov
,
M. W.
Thompson
,
P.
Zhang
,
C. H.
Turner
,
P. T.
Cummings
, and
S.
Dai
,
Appl. Phys. Lett.
109
,
143111
(
2016
).
7.
T.
Mo
,
S.
Bi
,
Y.
Zhang
,
V.
Presser
,
X.
Wang
,
Y.
Gogotsi
, and
G.
Feng
,
ACS Nano
14
,
2395
(
2020
).
8.
G.
Feng
,
X.
Jiang
,
R.
Qiao
, and
A. A.
Kornyshev
,
ACS Nano
8
,
11685
(
2014
).
9.
J.
Vatamanu
and
O.
Borodin
,
J. Phys. Chem. Lett.
8
,
4362
(
2017
).
10.
M. K.
Petersen
,
R.
Kumar
,
H. S.
White
, and
G. A.
Voth
,
J. Phys. Chem. C
116
,
4903
(
2012
).
11.
K. A.
Dwelle
and
A. P.
Willard
,
J. Phys. Chem. C
123
,
24095
(
2019
).
12.
J. I.
Siepmann
and
M.
Sprik
,
J. Chem. Phys.
102
,
511
(
1995
).
13.
D. T.
Limmer
,
C.
Merlet
,
M.
Salanne
,
D.
Chandler
,
P. A.
Madden
,
R.
Van Roij
, and
B.
Rotenberg
,
Phys. Rev. Lett.
111
,
106102
(
2013
); arXiv:1306.6903.
14.
B.
Uralcan
,
I. A.
Aksay
,
P. G.
Debenedetti
, and
D. T.
Limmer
,
J. Phys. Chem. Lett.
7
,
2333
(
2016
); arXiv:1604.03995.
15.
J.
Vatamanu
,
O.
Borodin
, and
G. D.
Smith
,
J. Phys. Chem. B
115
,
3073
(
2011
).
16.
C.
Merlet
,
C.
Péan
,
B.
Rotenberg
,
P. A.
Madden
,
P.
Simon
, and
M.
Salanne
,
J. Phys. Chem. Lett.
4
,
264
(
2013
).
17.
Z.
Wang
,
Y.
Yang
,
D. L.
Olmsted
,
M.
Asta
, and
B. B.
Laird
,
J. Chem. Phys.
141
,
184102
(
2014
); arXiv:1408.0839.
18.
S. K.
Reed
,
P. A.
Madden
, and
A.
Papadopoulos
,
J. Chem. Phys.
128
,
124701
(
2008
).
19.
D. T.
Limmer
,
A. P.
Willard
,
P.
Madden
, and
D.
Chandler
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
4200
(
2013
).
20.
L.
Scalfi
,
T.
Dufils
,
K. G.
Reeves
,
B.
Rotenberg
, and
M.
Salanne
,
J. Chem. Phys.
153
,
174704
(
2020
); arXiv:1910.13341.
21.
S.
Surendralal
,
M.
Todorova
,
M. W.
Finnis
, and
J.
Neugebauer
,
Phys. Rev. Lett.
120
,
246801
(
2018
).
22.
J.
Rossmeisl
,
J. K.
Nørskov
,
C. D.
Taylor
,
M. J.
Janik
, and
M.
Neurock
,
J. Phys. Chem. B
110
,
21833
(
2006
).
23.
F.
Deißenbeck
,
C.
Freysoldt
,
M.
Todorova
,
J.
Neugebauer
, and
S.
Wippermann
,
Phys. Rev. Lett.
126
,
136803
(
2021
); arXiv:2003.08156.
24.
T.
Dufils
,
G.
Jeanmairet
,
B.
Rotenberg
,
M.
Sprik
, and
M.
Salanne
,
Phys. Rev. Lett.
123
,
195501
(
2019
); arXiv:1907.00622.
25.
T.
Dufils
,
M.
Sprik
, and
M.
Salanne
,
J. Phys. Chem. Lett.
12
,
4357
(
2021
); arXiv:2104.03177.
26.
J.
Vatamanu
,
O.
Borodin
, and
G. D.
Smith
,
Phys. Chem. Chem. Phys.
12
,
170
(
2010
).
27.
A.
Marin-Laflèche
,
M.
Haefele
,
L.
Scalfi
,
A.
Coretti
,
T.
Dufils
,
G.
Jeanmairet
,
S.
Reed
,
A.
Serva
,
R.
Berthin
,
C.
Bacon
,
S.
Bonella
,
B.
Rotenberg
,
P.
Madden
, and
M.
Salanne
,
J. Open Source Softw.
5
,
2373
(
2020
).
28.
29.
E. R.
Smith
,
Proc. R. Soc. London, Ser. A
375
,
475
(
1981
).
30.
I.-C.
Yeh
and
M. L.
Berkowitz
,
J. Chem. Phys.
111
,
3155
(
1999
).
31.
D. E.
Parry
,
Surf. Sci.
49
,
433
(
1975
).
32.
Z.
Hu
,
J. Chem. Theory Comput.
10
,
5254
(
2014
).
33.
V.
Ballenegger
,
J. Chem. Phys.
140
,
161102
(
2014
).
34.
T. R.
Gingrich
and
M.
Wilson
,
Chem. Phys. Lett.
500
,
178
(
2010
).
35.
R. W.
Hockney
,
S. P.
Goel
, and
J. W.
Eastwood
,
Chem. Phys. Lett.
21
,
589
(
1973
).
36.
J. W.
Eastwood
,
R. W.
Hockney
, and
D. N.
Lawrence
,
Comput. Phys. Commun.
19
,
215
(
1980
).
37.
J. W.
Cooley
and
J. W.
Tukey
,
Math. Comput.
19
,
297
(
1965
).
38.
J. W.
Eastwood
,
J. Comput. Phys.
18
,
1
(
1975
).
39.
V.
Ballenegger
,
J. J.
Cerdà
, and
C.
Holm
,
J. Chem. Theory Comput.
8
,
936
(
2012
).
40.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
, 2nd ed. (
Oxford University Press
,
2017
), Vol. 1, Chap. 6, pp.
216
257
.
41.
See Gitlab.com/ampere2/metalwalls/-/tree/20.05/example/thomas-fermi-gold-capacitor for the MetalWalls input files of the original model capacitor system adopted in this work.
42.
C.
Merlet
,
D. T.
Limmer
,
M.
Salanne
,
R.
Van Roij
,
P. A.
Madden
,
D.
Chandler
, and
B.
Rotenberg
,
J. Phys. Chem. C
118
,
18291
(
2014
); arXiv:1404.0343.

Supplementary Material