NPT ensemble
============

The following integration scheme is used to sample the NPT ensemble. Note that both sections thermostat and barostat should be present in the *runtime.inpt* file. It is not possible to use the NPT ensemble with 2D PBCs or for constant potential simulations.

Metalwalls uses an isotropic barostat and two Nosé-Hoover [[Nosé1984][1] & [Hoover1985][2]] thermostat chains (of length $`n_1`$ and $`n_2`$ provided by the user, default values are 5 and 5) as described by Martyna *et al.* [[Martyna1994][3]]. The first one is associated with the particles and the second one with the barostat. The additional variables are noted $`\zeta_{pj}`$, $`\zeta_{bj}`$ and $`\epsilon`$ and their momentum $`v_{\zeta p j}`$, $`v_{\zeta b j}`$ (($`j\in[1,n_1]`$) or $`[1,n_2]`$, respectively) and $`v_{\epsilon}`$.The fictitious mass of the first thermostat is given by 
```math
Q_{p1}=N_{f}kT\tau_t^2  
```
where $`N_f`$ is the number of degrees of freedom (see [here](dof)) and $`\tau_t`$ 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_{pj}=kT\tau_t^2
```
The barostat has a mass
```math
W=(N_f+3)kT\tau_b
```
where $`\tau_b`$ is the barostat relaxation time provided by the user (default value is 20670.5 a.u. $`\approx`$ 500 fs). The thermostats in the barostat chain have a mass
```math
Q_{bj}=kT\tau_t^2
```
Note that they have the same relaxation time as the ones involved in the thermostat for the atoms. The time integration involves the following steps:


(NPT1) 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}}
```
(NPT2) Update of the positions and of the velocities due to the barostat/thermostats at $`t=0`$
```math
\mathbf{r}''_{i}(\Delta t) \leftarrow \mathbf{r}'_{i}(\Delta t)- [\mathbf{v}_{i}(0)v_{\zeta 1}(0)+(2+\frac{3}{N_f})(\mathbf{v}_{i}(0)v_{\epsilon}(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)+(2+\frac{3}{N_f})(\mathbf{v}_{i}(0)v_{\epsilon}(0))]\frac{\Delta t}{2}
```
(NPT3) Update of the barostat/thermostats variables
```math
\zeta_{pj}(\Delta t) \leftarrow \zeta_{pj}(0)+v_{\zeta p j}(0)\Delta t +G_{\zeta p j}(0)\frac{\Delta t^2}{2}\\
\zeta_{bj}(\Delta t) \leftarrow \zeta_{bj}(0)+v_{\zeta b j}(0)\Delta t +G_{\zeta b j}(0)\frac{\Delta t^2}{2}\\
\epsilon(\Delta t) \leftarrow \epsilon(0)+v_{\epsilon}(0)\Delta t +[\frac{F_{\epsilon}(0)}{W}-v_{\epsilon}(0)v_{\zeta b 1}(0)]\frac{\Delta t^2}{2}
```
where $`G_{\zeta p j}`$, $`G_{\zeta b j}`$ and $`F_{\epsilon}`$ are the ''forces'' acting on the barostat/thermostats variables (defined below)

(NPT4) Update the cell parameters ($`\alpha=x,y,z`$) and rescale the position and velocities accordingly
```math
L_\alpha(\Delta t) \leftarrow e^{\epsilon(\Delta t)-\epsilon(0)}  L_\alpha(0)\\
\mathbf{r}_{i}(\Delta t) \leftarrow e^{\epsilon(\Delta t)-\epsilon(0)}  \mathbf{r}''_{i}(\Delta t)\\
\mathbf{v}'''_{i}(\Delta t) \leftarrow e^{\epsilon(\Delta t)-\epsilon(0)} \mathbf{v}''_{i}(\Delta t)
```
Note that at this point the number of $`k`$-vectors (see the Ewald section) and their repartition between the various cores have to be redefined.

(NPT5) Enforce the constraints on the new positions and velocities with the first part of RATTLE

(NPT6) Compute the new electrolyte atom dipoles (polarizable force field) and radii (for the DAIM force field)

(NPT7) Compute the new forces at $`\Delta t`$

(NPT8) 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}}\\
```

(NPT9) Enforce the constraints on the new velocities with the second part of RATTLE

(NPT10) Make an initial guess for the barostat/thermostats variables velocities at time $`\Delta t`$
```math
v_{\zeta p j}^{(0)}(\Delta t) \leftarrow v_{\zeta p j}(-\Delta t)+2G_{\zeta p j}(0)\Delta t\\
v_{\zeta b j}^{(0)}(\Delta t) \leftarrow v_{\zeta b j}(-\Delta t)+2G_{\zeta b j}(0)\Delta t\\
v_{\epsilon}^{(0)}(\Delta t) \leftarrow v_{\epsilon}(-\Delta t)+2\frac{F_{\epsilon}(0)}{W}\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`$)

(NPT11a) 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)+(2+\frac{3}{N_f})v_{\epsilon}^{(k-1)}(\Delta t)]}
```
(NPT11b) Rescaling the velocities
```math
\mathbf{v}_{i}^{(k)}(\Delta t) \leftarrow \alpha^{(k)}\mathbf{v}_{i}^{(k-1)}(\Delta t)
```
(NPT11c) Update the internal pressure $`P^{(k)}`$

(NPT11d) Update the first thermostat in the chain for particles (force and velocity)
```math
G_{\zeta p 1}^{(k)}(\Delta t) \leftarrow \frac{1}{Q_{p1}}\left[\sum m_{i}\mathbf{v}^{(k)}_{i}(\Delta t)\mathbf{v}^{(k)}_{i}(\Delta t)-N_{f}kT \right]-v_{\zeta p 1}^{(k-1)}(\Delta t)v_{\zeta p 2}^{(k-1)}(\Delta t)\\
v_{\zeta p 1}^{(k)}(\Delta t) \leftarrow v_{\zeta p 1}(0)+\left[G_{\zeta p 1}(0) +G_{\zeta p 1}^{(k)}(\Delta t)\right]\frac{\Delta t}{2}
```
(NPT11e) Update the other thermostats in the chain for particles (forces and velocities)
```math
G_{\zeta p j}^{(k)}(\Delta t) \leftarrow \frac{1}{Q_{pj}}\left[Q_{p(j-1)}v_{\zeta p (j-1)}^{(k)}(\Delta t)v_{\zeta p (j-1)}^{(k)}(\Delta t)-kT \right]-v_{\zeta p j}^{(k-1)}(\Delta t)v_{\zeta p (j+1)}^{(k-1)}(\Delta t)\\
v_{\zeta p j}^{(k)}(\Delta t) \leftarrow v_{\zeta p j}(0)+\left[G_{\zeta p j}(0) +G_{\zeta p j}^{(k)}(\Delta t)\right]\frac{\Delta t}{2}
```
(NPT11f) Update the barostat force and velocity
```math
F_{\epsilon}^{(k)}(\Delta t) \leftarrow 3V(P^{(k)}-P_{external})+\frac{3}{N_f}\sum m_i \mathbf{v}^{(k)}_{i}(\Delta t)\mathbf{v}^{(k)}_{i}(\Delta t)-W v_{\epsilon}^{(k-1)}(\Delta t)v_{\zeta b 1}^{(k-1)}(\Delta t)\\
v_{\epsilon}^{(k)}(\Delta t) \leftarrow v_{\epsilon}(0)+[F_{\epsilon}(0)+F_{\epsilon}^{(k)}(\Delta t)]\frac{\Delta t}{2W}
```
(NPT11g) Update the first thermostat in the chain for barostat (force and velocity)
```math
G_{\zeta b 1}^{(k)}(\Delta t) \leftarrow \frac{1}{Q_{b1}}\left[W v^{(k)}_{\epsilon}(\Delta t)v^{(k)}_{\epsilon}(\Delta t)-kT \right]-v_{\zeta b 1}^{(k-1)}(\Delta t)v_{\zeta b 2}^{(k-1)}(\Delta t)\\
v_{\zeta b 1}^{(k)}(\Delta t) \leftarrow v_{\zeta b 1}(0)+\left[G_{\zeta b 1}(0) +G_{\zeta b 1}^{(k)}(\Delta t)\right]\frac{\Delta t}{2}
```
(NPT11h) Update the other thermostats in the chain for barostat (forces and velocities)
```math
G_{\zeta b j}^{(k)}(\Delta t) \leftarrow \frac{1}{Q_{bj}}\left[Q_{b(j-1)}v_{\zeta b (j-1)}^{(k)}(\Delta t)v_{\zeta b (j-1)}^{(k)}(\Delta t)-kT \right]-v_{\zeta b j}^{(k-1)}(\Delta t)v_{\zeta b (j+1)}^{(k-1)}(\Delta t)\\
v_{\zeta b j}^{(k)}(\Delta t) \leftarrow v_{\zeta b j}(0)+\left[G_{\zeta b j}(0) +G_{\zeta b j}^{(k)}(\Delta t)\right]\frac{\Delta t}{2}
```
(NPT11i) 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 (in the thermostat section, default value: 1.0e-15). The maximal number of iterations can also be provided in the input (in the thermostat section, 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.467468
[Martyna1994]
G. J. Martyna, D. J. Tobias and M. L. Klein.
Constant pressure molecular dynamics algorithms
J. Chem. Phys. 101, 4177 (1994)






