Force fields
==========

This section describes the various force fields implemented in MW. 
Each species type can interact through different potentials: interactomic potentials and intramolecular potentials. Each atom has a set of coordinates $`\mathbf{r}_i=(x_i, y_i, z_i)`$ and can have an electric charge, an electric dipole (characterized by an atomic polarizability when a polarizable force field is used) and an effective radius (characterized by an atomic deformability in the dipole aspherical ion model).

Throughout this page, we write the distance between two particles $`i`$ and $`j`$

```math
r_{ij} = \sqrt{(x_j-x_i)^2 + (y_j-y_i)^2 + (z_j-z_i)^2}
```

Atomic species are distinguished between electrode atoms and electrolyte atoms depending on their charge distribution. Electrode species are characterized by a Gaussian charge distribution of amplitude $`Q_i`$ and width $`\eta_i^{-1}`$

```math
\rho_i (\mathbf{r}) = Q_i \eta_i^3 \pi ^{-3/2} \exp(-\eta_i^2 |\mathbf{r} - \mathbf{r}_i|^2)
```

The value of the charge $`Q_i`$ can be kept constant or can fluctuate under a constant potential constraint. If so, the charges on the electrode atoms are computed at each timestep (for more information see [here](constant-potential)).

The electrodes are always treated as rigid bodies and can move only by considering them as pistons that translate along the $`z`$ direction, due to forces from the rest of the system and a pressure force $`F= P \times S`$ (see [here](NPzT-ensemble) for more information).

On the other hand, electrolyte species have a point charge distribution 

```math
\rho_i (\mathbf{r}) = q_i \delta (\mathbf{r} - \mathbf{r}_i)
```

They are usually mobile and can be polarizable (in polarizable force fields) and deformable (in the aspherical ion model).

In the following, we give potential energy expressions and the associated parameters that should be defined in the *runtime.inpt* file (see [here](system-configuration)). Given the interaction potential $`V(r_{ij})`$, the force felt by the particle $`i`$ due to particle $`j`$ is 

```math
\mathbf{F_{ij}} = - \mathbf{\nabla} V(r_{ij}) = - \frac{\partial V(r_{ij})}{\partial r_{ij}} \mathbf{\nabla} r_{ij} = \frac{\partial V(r_{ij})}{\partial r_{ij}} \frac{1}{r_{ij}}
\begin{pmatrix}
x_j - x_i \\
y_j - y_i \\
z_j - z_i 
\end{pmatrix}
```
and of course we have
```math
\mathbf{F_{ij}} = - \mathbf{F_{ji}}
```



* [Electrostatic interactions](#elec) 
   - [Coulomb potential](#coulomb)
   - [Polarization effects](#polarization)
* [van der Waals interactions](#vdW)
   - [Lennard-Jones](#lj)
   - [Fumi-Tosi](#ft)
   - [Extended Fumi-Tosi](#xft)
   - [Dipole Aspherical Ion Model](#daim)
   - [Steele wall](#steele)
* [Intramolecular interactions](#intra)
   - [Harmonic bond](intramolecular-potential)
   - [Harmonic angle](intramolecular-potential)
   - [Dihedral angle](intramolecular-potential)
   - [Improper angle](intramolecular-potential)


<a name="elec">

# Electrostatic interactions

</a>

<a name="coulomb">

## Coulomb potential

</a>

The Coulomb potential for a set of point charges and Gaussian charges $`\rho=\sum_i \rho_i`$ is given by
```math
V[\rho]=\frac{1}{2}\int d^3 \mathbf{r} d^3 \mathbf{r}'\frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{\mid \mathbf{r}-\mathbf{r}'\mid}
```
This expression leads to three different pair potentials:

- For point charge - point charge pairs
```math
V(r_{ij})=\frac{q_iq_j}{r_{ij}}
```
- For point charge - Gaussian charge pairs
```math
V(r_{ij})=\frac{Q_iq_j}{r_{ij}}{\rm erf}[\eta r_{ij}]
```
- For Gaussian charge - Gaussian charge pairs
```math
V(r_{ij})=\frac{Q_iQ_j}{r_{ij}}{\rm erf}[\sqrt{2}\eta r_{ij}]
```
And a self-interaction term for the Gaussian charges:
```math
V_{\rm self}=\frac{\eta}{\sqrt{2\pi}}Q_i^2
```

This potential is long-ranged and it cannot be truncated at a cut-off distance. It is calculated using the Ewald summation method, for which the expressions depend on the periodicity (2D or 3D) of the system. More informations on the implementation are given [here](ewald).

<a name="polarization">

## Polarization effects

</a>

These terms can be computed for systems of point charges only. The polarizable atoms carry an induced dipole $`\mathbf{\mu}_i`$ which interacts with all the point charges via the potential [[Wilson1993][1]]
```math
V(r_{ij}) =f^{ji}_4(r_{ij})\frac{q_j \mathbf{\mu}_i\cdot\mathbf{r}_{ij}}{r_{ij}^3}
```
and with all the other induced dipoles via the potential
```math
V(r_{ij}) =\mathbf{\mu}_i\mathbf{T}_{ij}\mathbf{\mu}_j
```
where $`\mathbf{T}_{ij}`$ is the dipole-dipole interaction tensor:
```math
\mathbf{T}_{ij}=\frac{1}{r_{ij}^3}\mathbf{I}-3\frac{1}{r_{ij}^5}\begin{bmatrix} x^2 & xy & xz \\ yx & y^2 & yz \\ zx & zy & z^2\end{bmatrix}
```
in which $`\mathbf{I}`$ is the identity matrix and $`x`$, $`y`$ and $`z`$ are the Cartesian components of the vector $`\mathbf{r}_{ij}`$.

The charge-dipole interactions are damped at short-range with Tang-Toennies functions [[Tang1984][2]] defined as
```math
f^{ij}_n(r_{ij}) = 1 - c^{ij}_n  \exp(-b^{ij}_n r_{ij}) \sum_{k=0}^{n} \frac{(b^{ij}_n r_{ij})^k}{k!}
```
The parameters needed for this damping function are $`n`$, $`b^{ij}_n`$ and $`c^{ij}_n`$. Note that generally $`b^{ij}_n = b^{ji}_n`$ but $`c^{ij}_n`$ differs from $`c^{ji}_n`$. 

In addition to the pair potentials, polarization effects include an additional term which represents the work needed to polarize the atom
```math
V=\frac{1}{2\alpha_i}\mu_i^2
```
where $`\alpha_i`$ is the polarizability of the atom, which is a parameter of the force field. The induced dipoles are additional degrees of freedom which are determined at each time step of the simulation following a procedure described [here](polarization-effects).

<a name="vdW">

# van der Waals interactions

</a>

Several potentials are implemented to describe the van der Waals interactions

<a name="lj">

## Lennard-Jones 

</a>

The Lennard-Jones potential is a simple pair potential that accounts for repulsion and dispersion [[Lennard-Jones1924][3]].

```math
V(r_{ij}) = 4 \epsilon_{ij} \left( \frac{\sigma_{ij}^{12}}{r_{ij}^{12}} - \frac{\sigma_{ij}^6}{r_{ij}^{6}} \right) = \frac{B_{ij}}{r_{ij}^{12}} - \frac{A_{ij}}{r_{ij}^{6}}
```

The $`\epsilon_{ij}`$ parameters define the depth of the energy minimum and the $`\sigma_{ij}`$ parameters the radial position of this minimum for each pair. There are mixing rules to derive cross parameters $`ij`$ from direct parameters $`ii`$ and $`jj`$, such as the Lorentz-Berthelot rule:

```math
\epsilon_{ij} = \sqrt{\epsilon_{ii}\epsilon_{jj}} \\
\sigma_{ij} = (\sigma_{ii} + \sigma_{jj})/2
```

<a name="ft">

## Fumi-Tosi

</a>

This potential is a simplified version of the Born-Huggins-Mayer potential (the asymptotic expansion for the dispersion term includes dipole-dipole and dipole-quadrupole terms) [[Tosi1964][4]]

```math
V_{ij}(r_{ij}) = B_{ij} \exp(-\eta _{ij} r_{ij}) - f_{dd}^{ij}(r_{ij}) \frac{C_{ij}}{r_{ij}^{6}} - f_{dq}^{ij}(r_{ij}) \frac{D_{ij}}{r_{ij}^{8}}
```
where Tang-Toennies functions are used to damp the dispersion terms at short range:
```math
f_{dd}^{ij}(r_{ij}) = 1 - \exp(-d_{dd}^{ij} r_{ij}) \sum_{k=0}^6 \frac{(d_{dd}^{ij} r_{ij})^k}{k!}
```
and 
```math
f_{dq}^{ij}(r_{ij}) = 1 - \exp(-d_{dq}^{ij} r_{ij}) \sum_{k=0}^8 \frac{(d_{dq}^{ij} r_{ij})^k}{k!}
```

The parameters of this potential are $`\eta_{ij}`$, $`B_{ij}`$, $`C_{ij}`$, $`D_{ij}`$, $`d_{dd}^{ij}`$, $`d_{dq}^{ij}`$ for each pair. There are mixing rules to derive cross parameters $`ij`$ from direct parameters $`ii`$ and $`jj`$. The one implemented in MW states that
```math
\eta_{ij} = (\eta_{ii}+\eta_{jj}) \frac{\eta_{ii}\eta_{jj}}{\eta_{ii}^2+\eta_{jj}^2} \\
B_{ij} = \sqrt{B_{ii}B_{jj}} 
```
and the rules for $`d_{dd}^{ij}`$ and $`d_{dq}^{ij}`$ are the same as for $`\eta_{ij}`$, while the ones for $`C_{ij}`$ and $`D_{ij}`$ are the same as for $`B_{ij}`$.

<a name="xft">

## Extended Fumi-Tosi (XFT)

</a>

This potential is an extension of the Fumi-Tosi potential, which was introduced to have a more repulsive wall at short-range

```math
V_{ij}(r_{ij}) =\frac{B_{ij} \exp(-\eta_{ij}r_{ij})}{r_{ij}^{n}}  + B'_{ij} \exp(-\eta'_{ij} r^2_{ij}) - f_{dd}^{ij}(r_{ij}) \frac{C_{ij}}{r_{ij}^{6}} - f_{dq}^{ij}(r_{ij}) \frac{D_{ij}}{r_{ij}^{8}}
```
where the same damping functions as for the Tosi-Fumi potential are used.

The parameters of this potential are $`\eta_{ij}`$, $`B_{ij}`$, $`n`$, $`\eta'_{ij}`$, $`B'_{ij}`$, $`C_{ij}`$, $`D_{ij}`$, $`d_{dd}^{ij}`$, $`d_{dq}^{ij}`$ for each pair. Mixing rules are similar to the ones of the Fumi-Tosi potential (but $`n`$ has to be the same for all the pairs).

<a name="daim">

## Dipole Aspherical Ion Model

</a>

This potential is an extension of the Fumi-Tosi potential to account for the impact of the deformation of the electronic clouds in condensed phases on the repulsion energy [[Madden1996][5]]

```math
V(r_{ij}) = B^1_{ij}\exp(-\eta^1_{ij} r_{ij}) + B^2_{ij}\exp(-\eta^2_{ij} \rho_{ij}) + B^3_{ij}\exp(-\eta^3_{ij} \rho_{ij}) - f_{dd}^{ij}(r_{ij})\frac{C_{ij}}{r_{ij}^6} - f_{dq}^{ij}(r_{ij})\frac{D_{ij}}{r_{ij}^8}
```
where

```math
\rho_{ij}=r_{ij}-\delta\sigma_{i}-\delta\sigma_{j} -\frac{1}{r_{ij}}\bm{r_{ij}}\cdot(\bm{\nu}_{i}-\bm{\nu}_{j})
```
in which $`\sigma_{i}`$ represents the deviation of the ionic radius of the particle $`i`$ and $`\bm{\nu}_i`$ expresses its distorsion in a dipolar shape. The set of quantities $`\{\sigma_{i},\nu_{ix},\nu_{iy},\nu_{iz}\}_{i\in N}`$ are additional degrees of freedom which are computed at each timestep by minimization of the total energy through a conjugate-gradient approach described [here](ion-deformation-effects). The same damping functions as for the Tosi-Fumi potential are used.

The parameters needed for this potential are  $`\eta^1_{ij}`$, $`\eta^2_{ij}`$, $`\eta^3_{ij}`$, $`B^1_{ij}`$, $`B^2_{ij}`$, $`B^3_{ij}`$, $`C_{ij}`$, $`D_{ij}`$, $`d_{dd}^{ij}`$, $`d_{dq}^{ij}`$. No mixing rules are implemented for this potential.

When this potential is used there is an additional self-energy term corresponding to the cost of the deformation:

```math
V_{self}=\sum_i D_i[\exp(\beta_i \delta\sigma_i)+\exp(-\beta_i\delta\sigma_i)]+\exp(\zeta_i^2\mid\mathbf{\nu}_{i}\mid^2)-1
```

in which $`D_i`$, $`\beta_i`$ and $`\zeta_i`$ are atomic parameters that control the deformability of the ion.

<a name="steele">

## Steele wall

</a>

The Steele potential is a one-dimensional potential used to model non-structured walls, or a 10-4-3 Lennard-Jones-type model. It is defined as

```math
V(|z_{i}-z_0|) = 2 \pi \rho_w \epsilon _{fw} \sigma _{fw} \Delta [2/5 \left(\frac{\sigma _{fw}}{z}\right) ^{10} - \left(\frac{\sigma _{fw}}{z}\right)^4 - \frac{\sigma_{fw}^4}{3 \Delta (z + 0.61 \Delta)^3}]
```
where $`z`$ is the distance to the wall. The parameters are $`\rho_w`$ (the solid wall density), $`\epsilon_{fw}`$ (the fluid-wall interaction energy), $`\sigma_{fw}`$ (the fluid-wall interaction length) and $`\Delta`$ (that corresponds to the distance between cristallographic planes in the solid).

<a name="intra">

# Intramolecular interactions

</a>

When molecular species are defined, additional intramolecular terms are included as described [here](intramolecular-potential).

References
----------
[1]:https://doi.org/10.1088/0953-8984/5/17/004
[Wilson1993]
M. Wilson and P. A. Madden.
Polarization effects in ionic systems from first principles.
J. Phys.: Condens. Matter 5, 2687 (1993)

[2]:https://doi.org/10.1063/1.447150
[Tang1984]
K. T. Tang and J. P. Toennies.
An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients.
J. Chem. Phys. 80, 3726 (1984)

[3]:https://doi.org/10.1098/rspa.1924.0082
[Lennard-Jones1924]
J. E. Lennard-Jones.
On the Determination of Molecular Fields.
Proc. R. Soc. Lond. A 106, 463 (1924)

[4]:https://doi.org/10.1016/0022-3697(64)90160-X
[Tosi1964]
M. P. Tosi and F. G. Fumi
Ionic sizes and born repulsive parameters in the NaCl-type alkali halides - II: The generalized Huggins-Mayer form.
J. Phys. Chem. Solids 25, 45 (1964)

[5]:https://doi.org/10.1039/CS9962500339
[Madden1996]
P. A. Madden and M. Wilson.
'Covalent' effects in 'ionic' systems
Chem. Soc. Rev. 25, 339 (1996)

