Polarization effects
====================

MetalWalls allows 3D bulk simulations of polarizable electrolytes. This is interesting for simulating systems where terms of the multipole expansion higher than the first, play an important role in the interactions among particles. In MW this is done in the framework of the Polarizable Ion Model (PIM) up to the dipolar order, as described, for example, in [[Aguado2003][1]]. This model considers point dipoles $`\bm{\mu}`$ associated to every polarizable particle.
Computation of induced dipoles was first introduced using a Car-Parrinello scheme [[Sprik1988][2]] and was later adapted to a Born-Oppenheimer-like approach. The method was also discussed in the framework of Mass-Zero (MaZe) constrained molecular dynamics in [[Coretti2018][3]].

The induced dipoles are calculated at each timestep of the simulation minimizing the energy with respect to the additional variables representing the point dipoles. The energy to minimize includes the coulomb interactions between different species due to the presence of the dipoles $`V_{pol}`$ and the amount of energy necessary to form the dipole. In particular, the form of the energy to be minimized is given by

```math
V(\bm{\mu})  = V_{pol}(\bm{\mu}) + \sum_{i=1}^{N}\frac{|\bm{\mu}|^2}{2\alpha_i}
``` 

where $`N`$ is the number of polarizable electrolyte atoms, $`\alpha`$ is the polarizability of the atom and $`\bm{\mu}=(\bm{\mu}_1,...,\bm{\mu}_N)`$. More detail on the model and the explicit form of the polarization energy used can be found in [[Aguado2003][1]] and [here](force-field).

The dipoles $`\bm{\mu}`$ that minimize the polarization energy are found solving the system of $`N`$ equations

```math
\bm{\nabla}_{\bm{\mu}_i}V(\bm{\mu}) = 0 \quad \Rightarrow \quad \frac{\bm{\mu}_i}{\alpha_i} = -\bm{\nabla}_{\bm{\mu}_i}V_{pol}(\bm{\mu})
```

These sets of equations rule the dynamics of the system for the two methods available in the code to perform simulations of polarizable force fields.

Born-Oppenheimer approach
-------------------------

A first method to solve the dynamics of the electrolyte atoms with a polarizable force field is to solve directly the equation

```math
\frac{\bm{\mu}_i}{\alpha_i} = -\bm{\nabla}_{\bm{\mu}_i}V_{pol}(\bm{\mu})
```

for each electrolyte configuration. It is possible to recognize in the RHS of the equation above the definition of the electric field generated by the surrounding particles on particle $`i`$, so that a more straightforward way to write the equation is
```math
\bm{\mu}_i = \alpha_i \bm{E}_i
```
Being the polarization energy quadratic in the induced dipoles (the electric field is therefore linear), this set of equations can be solved efficiently by means of iterative procedures such as the conjugate gradient method. For the same reason this equation could be solved also by direct matrix inversion, but this method becomes inefficient rapidly with the increase of the number of polarizable particles, as the matrix inversion has to be performed at each timestep.

Mass-Zero constrained MD (MaZe)
-------------------------------

Another approach, based on the method introduced in [[Coretti2018][3]], is to consider the motion of the electrolyte particles subject to the set of holonomic constraints given by 
```math
\bm{\nabla}_{\bm{\mu}_i}V(\bm{\mu}) = 0
```
This approach is based on an extended Lagrangian, where the induced dipoles are considered additional dynamical variables for the system. Contrary to already existing extended-Lagrangian approaches *à la* Car-Parrinello, the mass-zero limit is enforced rigorously in this method and the motion of the additional degrees of freedom is properly separated by 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 rigorously the adiabatic separation. In [[Bonella2020][4]] it has been proven that the dynamics also samples the Born-Oppenheimer distribution exactly. 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.1605941
[Aguado2003]
A. Aguado and P. A. Madden
Ewald summation of electrostatic multipole interactions up to the quadrupolar level
J. Chem. Phys. 119, 7471 (2003)

[2]:https://doi.org/10.1063/1.455722
[Sprik1988]
M. Sprik and M. L. Klein
A polarizable model for water using distributed charge sites 
J. Chem. Phys. 89, 7556 (1988)

[3]: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).

[4]: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)