NVT ensemble
============

Two integration schemes can be used to sample the NVT ensemble. 

* [Born Oppenheimer-type MD for the electrode atom charges or the electrolyte ion dipoles](#BO)
* [MaZe algorithm for the electrode atom charges or the electrolyte ion dipoles](#MaZe)

The first one is also used when no electrodes or dipoles are present. The latter cannot be used for the calculation of the ion radii in the DAIM force field.<br>
Note that the section barostat should not be present in the *runtime.inpt* file.

Metalwalls uses the Nosé-Hoover [[Nosé1984][1] & [Hoover1985][2]] thermostat chain (of length $`n`$ provided by the user, default value is 5) as described by Martyna *et al.* [[Martyna1992][3]]. The additional variables are noted $`\zeta_j`$ and their momentum $`v_{\zeta j}`$ ($`j\in[1,n]`$). The fictitious mass of the first thermostat is given by 
```math
Q_1=N_{f}kT\tau^2  
```
where $`N_f`$ is the number of degrees of freedom (see [here](dof)) and $`\tau`$ the thermostat relaxation time provided by the user (default value is 4134.1 a.u. $`\approx`$ 100 fs). The other thermostats in the chain have a mass
```math
Q_j=kT\tau^2
```

<a name="BO">

Born-Oppenheimer 
==========

</a>

(NVT-BO1) Update of the positions and of the velocities due to $`\mathbf{F}(0)`$
```math
\mathbf{r}'_{i}(\Delta t) \leftarrow \mathbf{r}_{i}(0)+\mathbf{v}_{i}(0)\Delta t + \mathbf{F}_{i}(0)\frac{\Delta t^2}{2 m_{i}}\\
\mathbf{v}'_{i}(\Delta t) \leftarrow \mathbf{v}_{i}(0)+\mathbf{F}_{i}(0)\frac{\Delta t}{2 m_{i}}
```
(NVT-BO2) Update of the positions and of the velocities due to the thermostat at $`t=0`$
```math
\mathbf{r}_{i}(\Delta t) \leftarrow \mathbf{r}'_{i}(\Delta t)- \mathbf{v}_{i}(0)v_{\zeta 1}(0)\frac{\Delta t^2}{2}\\
\mathbf{v}''_{i}(\Delta t) \leftarrow \mathbf{v}'_{i}(\Delta t)- \mathbf{v}_{i}(0)v_{\zeta 1}(0)\frac{\Delta t}{2}
```
(NVT-BO3) Update of the thermostat variables
```math
\zeta_j(\Delta t) \leftarrow \zeta_j(0)+v_{\zeta j}(0)\Delta t +G_{\zeta j}(0)\frac{\Delta t^2}{2}
```
where $`G_{\zeta j}`$ are the ''forces'' acting on the thermostat variables (defined below)

(NVT-BO4) Enforce the constraints on the new positions and velocities with the first part of RATTLE

(NVT-BO5) Compute the new electrolyte atom dipoles (polarizable force field), electrode atom charges (in constant potential simulations) and electrolyte atom radii (for the DAIM force field)

(NVT-BO6) Compute the new forces at $`\Delta t`$

(NVT-BO7) Make an initial guess for the velocities at time $`\Delta t`$ by adding the contribution due to $`\mathbf{F}(\Delta t)`$
```math
\mathbf{v}_{i}^{(0)}(\Delta t) \leftarrow \mathbf{v}''_{i}(\Delta t)+\mathbf{F}_{i}(\Delta t)\frac{\Delta t}{2 m_{i}}
```
(NVT-BO8) Enforce the constraints on the new velocities with the second part of RATTLE

(NVT-BO9) Make an initial guess for the thermostat variable velocities at time $`\Delta t`$
```math
v_{\zeta j}^{(0)}(\Delta t) \leftarrow v_{\zeta j}(-\Delta t)+2G_{\zeta j}(0)\Delta t
```
There is then an iterative procedure to determine the velocities (of both the atoms and of the thermostat variables) in a self-consistent way. The step of the procedure is noted $`k`$ ($`k > 1`$)

(NVT-BO10a) Determination of the scaling factor for the velocities
```math
\alpha^{(k)} \leftarrow \frac{1}{1+\frac{\Delta t}{2} v_{\zeta 1}^{(k-1)}(\Delta t)}
```
(NVT-BO10b) Rescaling the velocities
```math
\mathbf{v}_{i}^{(k)}(\Delta t) \leftarrow \alpha^{(k)}\mathbf{v}_{i}^{(k-1)}(\Delta t)
```
(NVT-BO10c) Update the first thermostat in the chain (force and velocity)
```math
G_{\zeta 1}^{(k)}(\Delta t) \leftarrow \frac{1}{Q_1}\left[\sum m_{i}\mathbf{v}^{(k)}_{i}(\Delta t)\mathbf{v}^{(k)}_{i}(\Delta t)-N_{f}kT \right]-v_{\zeta 1}^{(k-1)}(\Delta t)v_{\zeta 2}^{(k-1)}(\Delta t)\\
v_{\zeta 1}^{(k)}(\Delta t) \leftarrow v_{\zeta 1}(0)+\left[G_{\zeta 1}(0) +G_{\zeta 1}^{(k)}(\Delta t)\right]\frac{\Delta t}{2}
```
(NVT-BO10d) Update the other thermostats (forces and velocities)
```math
G_{\zeta j}^{(k)}(\Delta t) \leftarrow \frac{1}{Q_j}\left[Q_{j-1}v_{\zeta j-1}^{(k)}(\Delta t)v_{\zeta j-1}^{(k)}(\Delta t)-kT \right]-v_{\zeta j}^{(k-1)}(\Delta t)v_{\zeta j+1}^{(k-1)}(\Delta t)\\
v_{\zeta j}^{(k)}(\Delta t) \leftarrow v_{\zeta j}(0)+\left[G_{\zeta j}(0) +G_{\zeta j}^{(k)}(\Delta t)\right]\frac{\Delta t}{2}
```
(NVT-BO10e) The convergence is obtained when the kinetic energy difference is below $`N_f kT \times tolerance`$ where the $`tolerance`$ can be provided in the input (default value: 1.0e-15). The maximal number of iterations can also be provided in the input (default value: 50). If the velocities are not converged at this stage the program exits with an error.

<a name="MaZe">

MaZe 
==========

</a>

(NVT-MaZe1) Update of the positions and of the velocities due to $`\mathbf{F}(0)`$
```math
\mathbf{r}'_{i}(\Delta t) \leftarrow \mathbf{r}_{i}(0)+\mathbf{v}_{i}(0)\Delta t + \mathbf{F}_{i}(0)\frac{\Delta t^2}{2 m_{i}}\\
\mathbf{v}'_{i}(\Delta t) \leftarrow \mathbf{v}_{i}(0)+\mathbf{F}_{i}(0)\frac{\Delta t}{2 m_{i}}
```
(NVT-MaZe2) Update of the positions and of the velocities due to the thermostat at $`t=0`$
```math
\mathbf{r}_{i}(\Delta t) \leftarrow \mathbf{r}'_{i}(\Delta t)- \mathbf{v}_{i}(0)v_{\zeta 1}(0)\frac{\Delta t^2}{2}\\
\mathbf{v}''_{i}(\Delta t) \leftarrow \mathbf{v}'_{i}(\Delta t)- \mathbf{v}_{i}(0)v_{\zeta 1}(0)\frac{\Delta t}{2}
```
(NVT-MaZe3) Update of the thermostat variables
```math
\zeta_j(\Delta t) \leftarrow \zeta_j(0)+v_{\zeta j}(0)\Delta t +G_{\zeta j}(0)\frac{\Delta t^2}{2}
```
where $`G_{\zeta j}`$ are the ''forces'' acting on the thermostat variables (defined below)

(NVT-MaZe4) Enforce the constraints on the new positions and velocities with the first part of RATTLE

(NVT-MaZe5) The Verlet algorithm is used to provide an expression of the additional variables $`s`$ (electrode charges for constant potential simulations or electrolyte dipoles for polarizable force fields) at $`\Delta t`$ as a function of the Lagrange multipliers $`\gamma`$:
```math
s_i(\Delta t) \leftarrow 2s_i(0) - s_i(-\Delta t) + \Delta t^2\sum_k\gamma_k\frac{\partial^2 U}{\partial s_i\partial s_k}
```
Initialization details for the variables $`s`$ are provided [here](maze-initialization).

(NVT-MaZe6) The Lagrange multipliers $`\gamma`$ are found solving the equation of the adiabatic constraints $`\sigma(\mathbf{r}(\Delta t), s(\Delta t)) = 0`$ where $`\sigma = \frac{\partial U}{\partial s}`$.

(NVT-MaZe7) Compute the new forces at $`\Delta t`$

(NVT-MaZe8) Make an initial guess for the velocities at time $`\Delta t`$ by adding the contribution due to $`\mathbf{F}(\Delta t)`$
```math
\mathbf{v}_{i}^{(0)}(\Delta t) \leftarrow \mathbf{v}''_{i}(\Delta t)+\mathbf{F}_{i}(\Delta t)\frac{\Delta t}{2 m_{i}}
```
(NVT-MaZe9) Enforce the constraints on the new velocities with the second part of RATTLE

(NVT-MaZe10) Make an initial guess for the thermostat variable velocities at time $`\Delta t`$
```math
v_{\zeta j}^{(0)}(\Delta t) \leftarrow v_{\zeta j}(-\Delta t)+2G_{\zeta j}(0)\Delta t
```
There is then an iterative procedure to determine the velocities (of both the atoms and of the thermostat variables) in a self-consistent way. The step of the procedure is noted $`k`$ ($`k > 1`$)

(NVT-MaZe11a) Determination of the scaling factor for the velocities
```math
\alpha^{(k)} \leftarrow \frac{1}{1+\frac{\Delta t}{2} v_{\zeta 1}^{(k-1)}(\Delta t)}
```
(NVT-MaZe11b) Rescaling the velocities
```math
\mathbf{v}_{i}^{(k)}(\Delta t) \leftarrow \alpha^{(k)}\mathbf{v}_{i}^{(k-1)}(\Delta t)
```
(NVT-MaZe11c) Update the first thermostat in the chain (force and velocity)
```math
G_{\zeta 1}^{(k)}(\Delta t) \leftarrow \frac{1}{Q_1}\left[\sum m_{i}\mathbf{v}^{(k)}_{i}(\Delta t)\mathbf{v}^{(k)}_{i}(\Delta t)-N_{f}kT \right]-v_{\zeta 1}^{(k-1)}(\Delta t)v_{\zeta 2}^{(k-1)}(\Delta t)\\
v_{\zeta 1}^{(k)}(\Delta t) \leftarrow v_{\zeta 1}(0)+\left[G_{\zeta 1}(0) +G_{\zeta 1}^{(k)}(\Delta t)\right]\frac{\Delta t}{2}
```
(NVT-MaZe11d) Update the other thermostats (forces and velocities)
```math
G_{\zeta j}^{(k)}(\Delta t) \leftarrow \frac{1}{Q_j}\left[Q_{j-1}v_{\zeta j-1}^{(k)}(\Delta t)v_{\zeta j-1}^{(k)}(\Delta t)-kT \right]-v_{\zeta j}^{(k-1)}(\Delta t)v_{\zeta j+1}^{(k-1)}(\Delta t)\\
v_{\zeta j}^{(k)}(\Delta t) \leftarrow v_{\zeta j}(0)+\left[G_{\zeta j}(0) +G_{\zeta j}^{(k)}(\Delta t)\right]\frac{\Delta t}{2}
```
(NVT-MaZe11e) The convergence is obtained when the kinetic energy difference is below $`N_f kT \times tolerance`$ where the $`tolerance`$ can be provided in the input (default value: 1.0e-15). The maximal number of iterations can also be provided in the input (default value: 50). If the velocities are not converged at this stage the program exits with an error.

References
----------
[1]:https://doi.org/10.1063/1.447334
[Nosé1984]
S. Nosé.
A unified formulation of the constant temperature molecular dynamics methods 
J. Chem. Phys. 81, 511 (1984)

[2]:https://doi.org/10.1103/PhysRevA.31.1695
[Hoover1985]
W. G. Hoover
Canonical dynamics: Equilibrium phase-space distributions
Phys. Rev. A 31, 1695 (1985)

[3]:https://doi.org/10.1063/1.463940
[Martyna1992]
G. J. Martyna, M. L. Klein and M. Tuckerman.
Nosé-Hoover chains: The canonical ensemble via continuous dynamics
J. Chem. Phys. 97, 2635 (1992)
