Constant potential simulations
==============================

The originality of MetalWalls lies in the possibility of simulating materials 
at constant electric potential instead of constant charge. This is interesting for 
simulating conductors where the atom charge is not localized on the atom but 
can spread in the material and adapt to the external medium. 
In MW, this is done by allowing the charge on the material atoms to fluctuate 
under the condition of constant potential. The method was first introduced by Siepmann and Sprik [[Siepmann1995][1]] in the framework of the Car-Parrinello algorithm. It was later adapted to the case of 2D periodic boundary conditions and extended to a Born-Oppenheimer-like approach in the first version of MW [[Reed2007][2]]. The method was also discussed in the framework of Mass-Zero (MaZe) constrained molecular dynamics in [[Coretti2020][3]]. An updated description of the statistical mechanics of the constant potential ensemble can be found in the reference [[Scalfi2020][4]].

In MW the electrode charges are modeled as Gaussian charge distributions centered on the electrode charge site, with integrated charges $`Q_i`$ (and $`Q=(Q_1,...Q_N)`$). The electric potential generated by the electrode atoms is $`\frac{\partial V_c(Q)}{\partial Q_i}`$, where $`V_c(Q)`$ is the Coulomb energy of the system (see [here](ewald) for details). It follows that the constant electric potential condition within the electrode materials is given by
```math
\frac{\partial V_c(Q)}{\partial Q_i} = \Psi_i 
```

where $`\Psi_i`$ is the applied potential difference on the electrode sites. This is equivalent to
```math
\frac{\partial V(Q)}{\partial Q_i} = 0
```

where the total potential energy $`V(Q)`$ is defined as
```math
V(Q) = V_c(Q) - \sum_{i=1}^{N}\Psi_iQ_i
```

and includes the work necessary to charge the electrodes. The charges $`Q`$ that satisfy the constant potential condition are found solving the system of $`N`$ equations given above.

It is also frequent, in these kinds of simulations, to enforce the so-called electroneutrality constraint, defined by the equation
```math
\sum_{i=1}^{N}Q_i = 0
```

As the electrolyte ions undergo thermal motion, thereby changing the external potential, these equations must be solved to determine the electrode charges consistent with the constant potential condition. Two methods, described below, have been applied to achieve this.

## Born-Oppenheimer approach

A first method to simulate the dynamics of the electrolyte atoms under the constant potential condition is to solve directly the equation
```math
\frac{\partial V_c(Q)}{\partial Q_i} = \Psi_i
```

for each electrolyte configuration. This method relies on the quadratic form of the Coulomb energy with respect to the electrode charges to solve efficiently the equation for the $`Q`$. In particular, $`V_c`$ can be written in the form
```math
V_c(Q) = \frac{1}{2}Q^{T}AQ - Q^{T}b + c
```

which yields a set of equations for the charges given by $`AQ - b = \Psi`$, where the matrix $`A`$ depends only on the geometry of the electrodes while the vector $`b`$ depends also on the positions of the electrolyte atoms. The set of equations can then be solved by direct numerical matrix inversion yielding the charges 
```math
Q = A^{-1}(b+\Psi) 
```

This is very efficient when the electrodes are fixed as the matrix $`A`$ is constant and the inversion has to be performed only once at the beginning of the simulation. The MW implementation also allows constant pressure $`P_z`$ simulations to be run by treating the electrodes as [rigid pistons](NPzT-ensemble). In this case, a more efficient way to solve the linear system is by using the conjugate gradient iterative method.

The fulfillment of the electroneutrality constraint in this framework and its repercussions on the statistical mechanics of the constant-potential ensemble are discussed at length in [[Scalfi2020][4]].

## Mass-Zero constrained MD (MaZe)

Another approach, based on the method introduced in [[Coretti2018][5]], is to consider the motion of the electrolyte particles subject to the set of holonomic constraints given by 
```math
\frac{\partial V(Q)}{\partial Q_i} = 0
```

This approach is based on an extended Lagrangian, where the charges $`Q`$ are considered additional dynamical variables for the system. Unlike the already existing extended-Lagrangian approaches *à la* Car and Parrinello, the mass-zero limit is enforced rigorously in this method and the additional degrees of freedom are properly separated from the physical ones in an adiabatic sense. This yields dynamical equations of motion for the additional variables which can be solved numerically by symplectic and time-reversible algorithms, while keeping the adiabatic separation. In [[Bonella2020][6]] it has been proven that the dynamics also samples the Born-Oppenheimer distribution exactly. The application of the method to constant potential simulations is described in detail in [[Coretti2020][3]] together with the possibility of enforcing the electroneutrality constraint. The equations of motion for the additional variables are given for the [NVE](nve-ensemble) and [NVT](nvt-ensemble) ensembles.

References
----------
[1]:https://doi.org/10.1063/1.469429
[Siepmann1995]
J. I. Siepmann and M. Sprik.
Influence of Surface-Topology and Electrostatic Potential on Water Electrode Systems.
J. Chem. Phys. 102, 511 (1995)

[2]:https://doi.org/10.1063/1.2464084
[Reed2007]
S. K. Reed, O. J. Lanning and P. A. Madden.
Electrochemical Interface Between an Ionic Liquid and a Model Metallic Electrode.
J. Chem. Phys. 126, 084704 (2007)

[3]:https://doi.org/10.1063/5.0007192
[Coretti2020]
A. Coretti, L. Scalfi, C. Bacon, B. Rotenberg, R. Vuilleumier, G. Ciccotti, M. Salanne, and S. Bonella. 
Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems.
J. Chem. Phys. 152 (19), 194701 (2020).

[4]:https://doi.org/10.1039/C9CP06285H
[Scalfi2020]
L. Scalfi, D. T. Limmer, A. Coretti, S. Bonella, P. A. Madden, M. Salanne and B. Rotenberg.
Charge fluctuations from molecular simulations in the constant-potential ensemble.
Phys. Chem. Chem. Phys. 22, 10480 (2020)

[5]:https://doi.org/10.1063/1.5055704
[Coretti2018]
A. Coretti, S. Bonella, and G. Ciccotti. 
Communication: Constrained molecular dynamics for polarizable models.
J. Chem. Phys. 149 (19), 191102 (2018).

[6]:https://doi.org/10.1039/D0CP00163E
[Bonella2020]
S. Bonella, A. Coretti, R. Vuilleumier, and G. Ciccotti. 
Adiabatic motion and statistical mechanics *via* mass-zero constrained dynamics
Phys. Chem. Chem. Phys. 22, 10775 (2020)
 
