3.1. Energy Function

The ACT supports a range of different potentials for classical force field simulations. These potentials are described briefly in what follows.

3.1.1. Non-bonded interactions

3.1.1.1. Coulomb Interaction using Point charges

The Coulomb interaction can be computed with point charges:

(3.1)\[V_{coul}(r_{ij}) ~=~ \frac{q_i q_j} {4\pi\varepsilon_0\varepsilon_r r_{ij}}.\]

3.1.1.2. Coulomb Interaction using Gaussian charges

Alternatively, the Coulomb interaction can be described using Gaussian shielded charges

(3.2)\[V_{coul}(r_{ij}) ~=~ \frac{q_i q_j {\rm erf}\left(\zeta_{ij} r_{ij} \right)}{4\pi\varepsilon_0\varepsilon_r r_{ij}}, \quad \zeta_{ij} = \frac{\zeta_i \zeta_j}{\sqrt{\zeta_i^2 + \zeta_j^2}},\]

where \(\varepsilon_0\) is the permittivity of vacuum, \(q_{i,j}\) are the charges and \(\zeta_{i,j}\) are the charge distribution widths (screening factors). Note that Eqn. (3.2) includes a relative dielectric constant \(\varepsilon_r\) that can be used to parameterize a non-polarizable force field using charge-scaling [53].

3.1.1.3. Coulomb Interaction using Slater charges

In addition, Slater distributed charges can be used as described in ref. [30], from which the text below is an adapted version.

The spherical Slater orbital wavefunction is given by

(3.3)\[\psi_{n}({\bf r}) ~=~ \sqrt{\frac{(2\zeta)^{2n+1}}{4\pi(2n)!}} {\bf r}^{n-1}{\rm e}^{-\zeta {\bf r}}\]

where \(n\) is the highest quantum number of the element and \(\zeta\) is the exponent giving the width of the charge density. The distribution of charge can be described by the charge density, which is the square of the wave function ( (3.3) ). The Coulomb integral can then be written as [54, 55]:

(3.4)\[J_{ij}({\bf r}) \sim \int \int |\psi_n({\bf r_i})|^2\frac{q_i q_j}{|{\bf r_i}- {\bf r_j}|}|\psi_m({\bf r_j})|^2 d{\bf r_i} d{\bf r_j}\]

where \(\psi_n\) and \(\psi_m\) are the Slater wave functions of atoms \(i\) and \(j\) with quantum numbers \(n\) and \(m\) and with partial charges \(q_i\) and \(q_j\), respectively. Both (3.2) and (3.4) have a finite limit as \(r \rightarrow 0\). As a result, these functions are well behaved at small \(r\). A number of methods [56, 57] have been proposed to evaluate (3.4). Hentschke gives an analytical solution [55]:

(3.5)\[\begin{split}J_{ij}({\bf r}) = \begin{split} & \frac{1}{4 \pi \epsilon_0}\frac{q_iq_j}{\|{\bf r_i} - {\bf r_j}\|} \frac{4\zeta_i^{2n+1}\zeta_j^{2m+1}}{(2n)!(2m)!} \frac{\partial^{2n-2}\partial^{2m-2}}{\partial\zeta_i^{2n-2}\partial\zeta_j^{2m-2}} \left(\frac{1}{\zeta_i^3 \zeta_j^3}\right) \\ & \times\biggl[ 1 -\frac{(3\zeta_i^2 - \zeta_j^2)\zeta_j^4}{(\zeta_i - \zeta_j)^3(\zeta_i + \zeta_j)^3} e^{-2\zeta_ir_{ij}} -\frac{(\zeta_i^2 - 3\zeta_j^2)\zeta_i^4}{(\zeta_i - \zeta_j)^3(\zeta_i + \zeta_j)^3} e^{-2\zeta_jr_{ij}} \\ & - \frac{\zeta_i\zeta_j^4}{(\zeta_i-\zeta_j)^2(\zeta_i+\zeta_j)^2} r_{ij} e^{-2\zeta_ir_{ij}} - \frac{\zeta_i^4\zeta_j}{(\zeta_i-\zeta_j)^2(\zeta_i+\zeta_j)^2} r_{ij} e^{-2\zeta_jr_{ij}} \biggr]. \end{split}\end{split}\]

Eqn. (3.5) was implemented in a Mathematica ^TM program from which C++ code was generated for the analytical computation of \(J_{ij}\) and its analytical derivatives with respect to \(r\), which are necessary for computing forces. Due to the nature of Eqn. (3.5), there are many terms with large powers, particularly for \(n>3\). Thus, the equations have to be implemented using the arbitrary precision arithmetic library Class Library for Numbers ( CLN ). to avoid numerical instabilities. It should be noted that the arbitrary precision library significantly increases the computational cost to analytically solve Eqn. (3.5)

3.1.1.4. Long range Coulomb interactions

It is good practice [58] to use the particle-mesh Ewald [59, 60] method to treat long-range electrostatics interactions in the condensed phase. In short, the method splits the Coulomb potential, either Eqn. (3.1), Eqn. (3.2) or something similar into a short and a long-range part as follows

(3.6)\[V_{coul}(r_{ij}) ~=~ V_{coul} \,{\rm erfc}(\alpha r_{ij}) + V_{coul}\, {\rm erf}(\alpha r_{ij})\]

using the fact that the complementary error function {rm erfc} equal 1 - {rm erf}. The first term here is the short-range part, that is computed in real space, and the second term is the long-range part that is computed in Fourier (reciprocal) space by moving charges on a regular grid [60]. The constant \(\alpha\) determines how quickly the short-range potential decays to zero, and it can be computed from a user-specified relative error tolerance related to moving charges on a grid. Given the cut-off distance \(rc\) and an error tolerance \(\epsilon\), typically \(10^{-4}\), we have

(3.7)\[\alpha ~=~ \frac{\sqrt{-{\rm ln}\,2\epsilon}}{rc}\]

where ln is the natural logarithm. This shows that \(\alpha\) has the dimension of 1 over distance in the units of rc. Fig. 3.1 shows how the division over short and long-range interactions works in practice, and how the long range contribution should be incorporated into simulations using a modified Coulomb function.

Fig. 3.1 Total Coulomb energy V as a function of distance r for two like unit charges, Gaussian-shielded Coulomb and short and long-range parts of the energy function. Ewald-tolerance \(\epsilon\) 0.0001, rc 0.9 nm, \(\alpha\) 3.24269/nm and Gaussian screening \(\zeta\) 5/nm.

3.1.1.5. Polarization

Polarization can be treated explicitly as well in the ACT using the core-shell model [32, 33]

(3.8)\[V_{pol}(r_{cs}) ~=~ \frac{1}{4\pi\varepsilon_0} \frac{q_s^2}{2\alpha_c}r_{cs}^2,\]

where \(q_s\) is the shell charge, \(\alpha_c\) is the polarizability and \(r_{cs}\) the distance between core and shell particles.

3.1.1.6. Induction correction

A term to add extra attraction between atoms was proposed in ref. [61]

(3.9)\[V_{ic}(r_{ij}) ~=~ -A_{ij} {\rm e}^{-b_{ij} r_{ij}}\]

with \(A_{ij}\) and \(b_{ij}\) constants to be trained. The geometric combination rule is used for \(A_{ij}\) and the arithmetic combination rule for \(b_{ij}\). It is not certain whether the exponential functional form is optimal to describe this interaction and there is no theory behind this.

3.1.1.7. Pauli Repulsion and London Dispersion

The repulsion and dispersion interactions acting between atoms \(i\) and \(j\) can be described by a number of potentials. In a recent study on noble gases [25] we found that the 14-7 potential due to Halgren [62] was one of the most accurate ones:

(3.10)\[V_{14-7}(r_{ij}) ~=~ \epsilon_{ij} \left(\frac{1+\delta_{ij}}{\frac{r_{ij}}{\sigma_{ij}} + \delta_{ij}} \right)^7\left( \frac{1+\gamma_{ij}}{(\frac{r_{ij}}{\sigma_{ij}})^7 +\gamma_{ij}} -2 \right)\]

where \(\epsilon_{ij}\) is the well depth at minimum, \(\gamma_{ij}\) and \(\delta_{ij}\) are dimensionless numbers that were originally shared for all elements. However, in our previous work [25], we treated \(\gamma\) and \(\delta\) as free atom-specific parameters subject to optimization and combination rules.

Furthermore, the generalized 4-parameter Buckingham potential with an adjustable long-range attraction is implemented as [63]

(3.11)\[V_{GBH}(r_{ij}) = \epsilon_{ij}\frac{\delta_{ij} + 2\gamma_{ij} + 6}{2\gamma_{ij}} \frac{1}{1+ (\frac{r}{\sigma_{ij}})^6 } \left[ \frac{6+\delta_{ij}}{\delta_{ij} +2\gamma_{ij}+6} e^{\gamma_{ij} (1-\frac{r}{\sigma_{ij}})} - 1 \right] -\frac{\epsilon_{ij}}{1+(\frac{r_{ij}}{\sigma})^{\delta}}\]

where \(\gamma\) and \(\delta\) again are dimensionless constants.

Another potential that is supported is the buffered (Wang) Buckingham potential [64]:

(3.12)\[V_{WBH}(r_{ij}) = \left(\frac{2\epsilon_{ij}}{1-\frac{3}{\gamma_{ij}+3}}\right) \left( \frac{\sigma_{ij}^6}{\sigma_{ij}^6+r_{ij}^6} \right) \left[\frac{3}{\gamma_{ij}+3}e^{\displaystyle{\gamma_{ij}\left(1-\frac{r_{ij}}{\sigma_{ij}}\right)}} - 1 \right]\]

where \(\epsilon_{ij}\) is the well depth at minimum, \(\sigma_{ij}\) is the Van der Waals radius and \(\gamma_{ij}\) determines the steepness of the potential and \(r_{ij}\) is the distance between particles \(i\) and \(j\). The buffered Buckingham potential mentioned above have been used to develop an accurate phase-transferable model for alkali-halides [19, 20, 21, 22, 23] and the original Buckingham potential [65] is implemented as

(3.13)\[V_{BH}(r_{ij}) ~=~ A_{ij}{\rm e}^{\displaystyle -b_{ij}r_{ij}} - \frac{C_{ij}}{r_{ij}^6}\]

where \(A_{ij}\), \(b_{ij}\) and \(C_{ij}\) are parameters to be trained. An alternative way of formulating the Buckingham potential [66] is

(3.14)\[E_{MBH} = \frac{\epsilon}{1-\frac{6}{\gamma}} \left[ \frac{6}{\gamma} e^{\gamma (1-\frac{r}{\sigma})} - \left(\frac{\sigma}{r}\right)^6 \right]\]

where \(\epsilon\) is the well-depth \(\sigma\) reflects the position of the minimum and \(\gamma\) is again a dimensionless constant. Although mathematically identical to Eqn. (3.13) we have shown that the potentials behave differently when combination rules are applied [25].

The well-known Tang-Toennies potential [67, 68] containing five parameters is implemented according to:

(3.15)\[V_{TT}(x) = Ae^{-b x} - \sum_{n=3}^5\left[1-e^{-b x}\sum_{k=0}^{2n}\frac{(b x)^k}{k!}\right]\frac{C_{2n}}{x^{2n}}\]

where repulsion and dispersion terms shared the parameter \(b\). Here, all five parameters can be trained by the ACT. This potential was used in our recent work on noble gases published where we investigated combination rules [25].

In other works, the Tang-Toennies potential is extended with an additional parameter [69] giving the exponential functions a different decay length \(b\):

(3.16)\[V_{TT2}(x) = Ae^{-b_{exp} x} - \sum_{n=3}^5\left[1-e^{-b_{disp} x}\sum_{k=0}^{2n}\frac{(b_{disp} x)^k}{k!}\right]\frac{C_{2n}}{x^{2n}}\]

and the ACT supports both these functions. Van Vleet et al. described a more accurate formula for the Pauli repulsion [70] that was applied to derive a force field later [71]. In this model the exchange repulsion is given by the Slater-ISA formula

(3.17)\[V_{exch} ~=~ A\left(\frac{1}{3}(br)^2 + br + 1\right){\rm exp}^{-br}\]

where it can be noted, that no additional parameters are needed for equation~:eq:slater_isa compared to the repulsion part in the Tang-Toennies potential (Eqn. (3.15)). In the ACT the Slater-ISA exchange can be combined with the damped dispersion part of the latter potential.

The Lennard-Jones 12-6 potential [72] is available as well:

(3.18)\[V_{LJ}(r_{ij}) ~=~ 4\varepsilon_{ij}\left(\frac{\sigma_{ij}^{12}}{r_{ij}^{12}}-\frac{\sigma_{ij}^{6}}{r_{ij}^6}\right)\]

and, in addition, the 12-6-4 potential

(3.19)\[V_{LJ}(r_{ij}) ~=~ 4\varepsilon_{ij}\left(\frac{\sigma_{ij}^{12}}{r_{ij}^{12}}-\frac{\sigma_{ij}^{6}}{r_{ij}^6}\right) - \frac{\gamma}{r^4}\]

which is meant to model ion-dipole interactions is supported too.

It is worth noting that in all these potentials (except Buckingham, Eqn. (3.13)) parameters describe both the exchange and dispersion interactions at the same time, which necessitates simultaneous training of these interaction. To account for anisotropy in exchange, a correction can be implemented using a virtual site particle on the atom that has a \(\sigma\)-hole interacting with other particles [26]:

(3.20)\[V_{EXCH,Corr}(r_{ij}) = -A_{ij}{\rm e}^{\displaystyle -b_{ij}r_{ij}}\]

where \(A_{ij}\) and \(b_{ij}\) are parameters to be optimized. This correction term (Eqn. (3.20)) does not need to be applied to all possible atom pairs, therefore a new combination rule, dubbed “Kronecker” is introduced

\[A_{ij} ~=~ (1-\delta_{ij})\frac{A_i + A_j}{2}\]

where \(\delta_{ij}\) is the Kronecker delta operating on particle types \(i\) and \(j\). This means that only interactions between \(\sigma\)-holes (represented by a virtual site) and atoms are non-zero. To avoid parameter explosion, the atomic \(A_i\) are set to zero and only the virtual site \(A_j\) is non-zero. This is reasonable since the virtual site represents a property of the \(\sigma\)-hole of the atom it is connected to. The arithmetic combination rule is used for \(b_{ij}\).

Multiple different combination rules can be used for parameters involved in van der Waals interactions, both within ACT and through the interface to OpenMM. For details, see the paper by Kříž et al. [25].

3.1.1.8. Treatment of long range dispersion

The different Van der Waals potentials described above can be used with Lennard-Jones PME [59, 73] if treated carefully. LJ-PME assumes a simple dispersion interaction given by

(3.21)\[V_{LJPME}(r_{ij}) ~=~ -\frac{C_{ij}}{r_{ij}^6}\]

for atoms \(i\) and \(j\), which differs from the potentials mentioned above. For LJ-PME, it is beneficial to use the geometric combination rule, hence

\[V_{LJPME}(r_{ij}) ~=~ -\frac{C_iC_j}{r_{ij}^6}\]

where \(C_i\) and \(C_j\), in contrast to Eqn. (3.21), have units of distance to the third power.

To minimize the difference between the potential used and the LJ dispersion over the whole volume outside the cut-off, and to specify correct inputs to OpenMM, starting from e.g. Eqn. (3.12) we have to solve the following equation:

(3.22)\[0 ~=~ \int_{rc}^{\infty} 4\pi r_{ij}^2\left[\left(-\frac{2\epsilon_{ij}}{1-\frac{3}{\gamma_{ij}+3}}\right) \left( \frac{\sigma_{ij}^6}{\sigma_{ij}^6+r_{ij}^6} \right) + \frac{C_iC_j}{r_{ij}^6} \right] {\rm d}r_{ij}\]

where we integrate from the cut-off \(rc\) to infinity. In this notation, we have to insert the combination rules for \(\epsilon\), \(\sigma\) and \(\gamma\). For the special case that \(i\) = \(j\), the constant \(C_{ii}\) can be derived using Mathematica ^TM as

\[C_{ii} ~=~ \epsilon_{i}\frac{3+\gamma_{i}}{\gamma_{i}}rc^3\sigma_{i}^3\left(\pi-2 {\rm arctan}\left[\left(\frac{rc}{\sigma_{i}}\right)^3\right]\right)\]

which we can then convert back to a \(\sigma'\) compatible with standard Lennard-Jones by

\[\sigma'_{i} = \left(\frac{C_{ii}}{4\epsilon_{i}}\right)^{1/6}.\]

Note however, that we need to get effective \(\sigma'\) for all atoms, which means the problem turns into a minimization problem where

(3.23)\[\varepsilon^2 ~=~ \sum_{i,j}\left[ \int_{rc}^{\infty} 4\pi r_{ij}^2\left[\left(-\frac{2\epsilon_{ij}}{1-\frac{3}{\gamma_{ij}+3}}\right) \left( \frac{\sigma_{ij}^6}{\sigma_{ij}^6+r_{ij}^6} \right) + \frac{C_iC_j}{r_{ij}^6} \right] {\rm d}r_{ij}\right]^2\]

has to be minimized with respect to \(C_i\) independently for each set of combination rules. In the Alexandria alkali-halide model we used the Wang-Buckingham potential (Eqn. (3.12)) in conjunction with the combination rules due to Kong [74] in which \(\sigma_{ij}\) depends on all the \(\sigma\), \(\epsilon\) and \(\gamma\). This makes it cumbersome to analytically derive integrals like Eqn. (3.23) for many combinations of potentials and combination rules.

Another problem is, that the parameters \(C_i\) that we need to compute the long-range dispersion interaction now depend on the cut-off, which means they should preferably not be computed beforehand. In summary, the \(C_i\) should be derived numerically given the potential and combination rules at the start of a simulation.

3.1.1.9. Combination Rules for Van der Waals Potentials

In a recent paper [25] we studied the effect of combination rules on potentials for noble gases and introduced two new combination rules (Eqn. (3.34) and Eqn. (3.35) ). Here is a brief summary, copied (with permission) from that paper.

Combination rules reduce the number of parameters required for the pair-wise potentials introduced above because the parameters describing the interaction between dissimilar X–Y atoms are reconstructed from parameters of homodimers X–X and Y–Y. In this way, it is necessary only to fit atomic parameters on data for homodimers. Different sets of mathematical expression were considered, as detailed below.

The two simplest expressions that have historically been used as combination rules are the geometric [75] and arithmetic [76] averages, respectively:

(3.24)\[X_{12} = \sqrt{x_{1} x_{2}}\]

%was used for all three parameters \(\epsilon\), \(\gamma\), \(\sigma\) in “geometric” rules. It is also used for \(\epsilon\) in “arithmetic” rules and for \(\sigma\) in “Kong-Mason” rules [74].

(3.25)\[X_{12} = \frac{x_{1}+ x_{2}}{2}\]

where \(x_1\) and \(x_2\) are the atomic parameters. Both these rules can in principle be applied to all parameters in the Van der Waals potentials described above and the rules are well-behaved mathematically.

Hogervorst introduced a set of combination rules for 12-6 Lennard-Jones and exp-6, modified Buckingham, potential (Eqn. (3.14)) [77]. He proposed using Eqn. (3.26) for \(\epsilon\) for both potential forms. For the \(\sigma\) of 12-6 potential (Eqn. (3.18)), the arithmetic mean (Eqn. (3.25)) was used. In the case of the modified Buckingham potential (eqn. (3.14)), expression~:eq:cr_sigma5, which depends on the combined parameters \(\epsilon_{1,2}\) and \(\gamma_{1,2}\) was advocated for the \(\sigma\), along with arithmetic mean for \(\gamma\).

(3.26)\[X_{12} = \frac{2 x_{1} x_{2}}{x_{1} + x_{2}}\]

(3.27)\[\sigma_{12}^6 = \sqrt{\frac{\epsilon_1\gamma_1\sigma_1^6}{\gamma_1 - 6} \frac{\epsilon_2\gamma_2\sigma_2^6}{\gamma_2 - 6}} \frac{(\gamma_{1,2}-6)}{\gamma_{1,2}\epsilon_{1,2}}\]

Where the \(\gamma\)^1,2 used Eqn. (3.25) and \(\epsilon\)^1,2 Eqn. (3.26). It should be noted that equation (3.26) is ill-behaved if both \(x_1\) and \(x_2\) are zero, while Eqn. (3.27) is ill-behaved if either \(\gamma_{1,2}\) or \(\epsilon_{1,2}\) is zero. Yang et al. [78] introduced an expression for the Morse potential [79], using Eqn. (3.26) for \(\epsilon\), while eqn. (3.28) is used for \(\sigma\) and \(\gamma\).

(3.28)\[X_{12} = \frac{x_{1} x_{2}(x_1 + x_2)}{x_{1}^2 + x_{2}^2}\]

The expression for \(\gamma\) proposed by Mason [80] for the exp-6 potential has the form %red{while used other rels for epsilon, sigma in alexandria}

(3.29)\[\gamma_{12} = \sqrt{\sigma_1 \sigma_2} \left(\frac{\gamma_1}{2\sigma_1}+ \frac{\gamma_2}{2\sigma_2}\right)\]

Waldman and Hagler [81] introduced expressions (3.30) for \(\epsilon\) and (3.31) for \(\sigma\)

to reproduce experimental well-depths and interaction distances.

(3.30)\[\begin{split}\epsilon_{12} = \sqrt{\epsilon_1 \epsilon_2}\frac{2\sigma_1^3\sigma_2^3}{\sigma_1^6+\sigma_2^6}\\\end{split}\]

(3.31)\[X_{12} = \left(\frac{X_1^6+X_2^6}{2}\right)^{1/6}\]

Although Eqn. (3.31) was devised for \(\sigma\) we have evaluated it for other parameters here as well, hence the notation with \(X\). Qi and coworkers advocated the use of buffered 14-7 Lennard-Jones potential (Eqn. (3.10)) due to Halgren [62], alongside combination expressions (3.30) for \(\epsilon\) and (3.32) for \(\sigma\): [82]

(3.32)\[X_{12} = \frac{x_1^3+x_2^3}{x_1^2+x_2^2}\]

A further relation, the harmonic mean rule, was proposed by Halgren [62]:

(3.33)\[X_{12} = \frac{4 x_1 x_2}{(x_1^{1/2}+x_2^{1/2})^2}\]

Finally, we introduce two new combination rules that we have not seen published previously. Since the \(\sigma\) in most potentials can be interpreted as a Van der Waals radius, we introduced a relation averaging third powers, corresponding to an atomic volume [25]:

(3.34)\[\sigma_{12} = \left(\frac{\sigma_1^3+\sigma_2^3}{2}\right)^{1/3}.\]

In addition, we applied the following rule for in particular \(\epsilon\) since it yield an \(X_{12}\) that is smaller than the geometric one (Eqn. (3.24)):

(3.35)\[X_{12} = \left(\frac{2}{x_1^{-2}+x_2^{-2}}\right)^{1/2}\]

The combination relations described above were permuted with each other into new combination rules. In this way, relations that depend on only one parameter type were used for any parameter. Relations depending on multiple parameters were used only for the specific parameter type combination they depend on (e.g. Eqn. (3.30) was only used for \(\epsilon\), using homodimer \(\epsilon\) and \(\sigma\)). In our previous work on alkali halides [19] we used combination rules according to eqn. (3.27) for \(\sigma\), eqn. (3.26) for \(\epsilon\) and eqn. (3.25) for \(\gamma\) with the Wang-Buckingham potential (Eqn. (3.12)).

A number of these combination rules can be written using the generalized mean equation [83]

(3.36)\[M_p(x_1,x_2,...,x_N) = \left(\frac{1}{N}\sum_{i=1}^N x_i^p\right)^{1/p}\]

for \(p \ne 0\). For \(p = 0\), Eqn. (3.36) turns into the geometric rule. Table 3.1 lists other well-known combination rules and their respective exponent \(p\). Hohm also describes combinations of the generalized mean and other similar expressions, but the most important observation he made was that the exponent \(p\) can be varied at will [83]. In the ACT it is possible to train this parameter along with the Van der Waals parameters.

Table 3.1 Correspondence between well-known combination rules and the generalized mean equation (3.36).

p	Name	Equation
-2	Inverse Square	(3.35)
-1	Hogervorst \(\epsilon\)	(3.26)
-1/2	Harmonic Mean	(3.33)
0	Geometric	(3.24)
1	Arithmetic	(3.25)
3	Volumetric	(3.34)
6	Waldman-Hagler \(\sigma\)	(3.31)

3.1.2. Bonded interactions

3.1.2.1. Harmonic potential

Bond vibrations can be described using a harmonic term based on the bond length \(r_{ij}\)

(3.37)\[V_b(r_{ij}) ~=~ \frac{k_{ij}^{b}}{2}\left(r_{ij}-r_{ij}^0\right)^2,\]

where \(k_{ij}^{b}\) is the force constant, and \(r_{ij}^0\) is the equilibrium bond length.

3.1.2.2. Morse potential

A Morse potential [79] can be used, with one addition:

(3.38)\[V_M(r_{ij}) ~=~ D_{ij}^e\left[{\rm e}^{\displaystyle{-2\beta_{ij}(r_{ij}-r_{ij}^0)}}-2{\rm e}^{\displaystyle{-\beta_{ij}(r_{ij}-r_{ij}^0)}}\right] + D_{ij}^0\]

The term \(D_{ij}^e\) roughly corresponds to a dissociation energy, however since there are Coulomb and/or Buckingham interactions between the atoms as well, the total “bond” potential is given by the sum of three terms and a correction term \(D_{ij}^0\) is needed to get the correct energy minimum.

3.1.2.3. Wei-Hua potential

In a very recent study we found the potential due to Hua [84, 85] to be the best compromise between accuracy of the vibrational frequencies and computational cost [27]. It is given by:

(3.39)\[U(r) ~=~ D_e\left(\left[ \frac{1-e^{-b(r-r_e)}}{1-ce^{-b(r-r_e)}} \right]^2 - 1 \right)\]

where \(D_e\) is the well-depth, \(r_e\) the equilibrium bond length, and \(b\) and \(c\) are constants with \(\|c\| < 1\).

3.1.2.4. Angle potential

Angle vibrations are described using a harmonic term based on the angle \(\theta_{ijk}\)

(3.40)\[V_a(\theta_{ijk}) ~=~ \frac{k_{ijk}^{\theta}}{2}\left(\theta_{ijk}-\theta_{ijk}^0\right)^2,\]

where \(k_{ijk}^{\theta}\) is the force constant, and \(\theta_{ijk}^0\) is the equilibrium angle.

3.1.2.5. Angle potential for linear compounds

The reference position, corresponding to a minimum energy structure, \(\mathbf{x}_j^0\) for a central atom \(j\) in a linear triplet of atoms \(i,j,k\) is given by

(3.41)\[\mathbf{x}_j^0 ~=~ {\rm a}\, \mathbf{x}_i + (1-{\rm a})\, \mathbf{x}_k\]

where \(a\) is a constant defined by the bond-lengths \(i-j\) and \(j-k\). In a group with bonds \(i-j\) and \(j-k\) with lengths \(b_{ij}\) and \(b_{jk}\) respectively, the constant is

(3.42)\[a ~=~ \frac{b_{jk}}{b_{ij}+b_{jk}}.\]

If the order of atoms is flipped \(a\) will change to \(1-a\). The potential \(V_{lin}\) is then given by

(3.43)\[V_{lin} ~=~ \frac{{\rm k}_{lin}}{2} \left(\mathbf{x}_j - \mathbf{x}_j^0\right)^2\]

with k _lin the force constant [24].

3.1.2.6. Out-of-plane vibrations

Finally, out-of-plane vibrations are treated by another harmonic potential

(3.44)\[V_i(\phi_{ijkl}) ~=~ \frac{k_{ijkl}^{\phi}}{2}\phi_{ijkl}^2,\]

where \(k_{ijkl}^{\phi}\) is the force constant and \(\phi_{ijkl}\) is defined by the angle between the two planes \(i,j,k\) and \(j,k,l\). This potential was historically termed improper dihedral.

3.1.2.7. Torsion potential

A torsion potential is implemented using a Fourier series:

(3.45)\[V_{d}(\phi_{ijkl}) ~=~ \sum_{n=0}^5 c_n {\rm cos}^n(\pi+\phi_{ijkl})\]

where \(c_n\) are constants and the torsion angle is defined as above. The constant \(\pi\) is added to be compatible with the Ryckaert-Bellemans potential [86] that is implemented in simulation codes like GROMACS [87] and OpenMM [88].

3.1.2.8. Proper dihedral

A simpler torsion potential is implemented for backward compatibility as:

(3.46)\[V_{d}(\phi_{ijkl}) ~=~ {\rm cos}(n \phi_{ijkl} + \phi_0)\]

where \(n\) is the multiplicity (number of minima in 360 degrees) \(\phi_0\) is an offset angle.

3.1.3. Special potentials

The ACT includes a flat-bottom position restraint potential according to

(3.47)\[\begin{split}V_{fbpr}(r) = \begin{cases} 0 & \hspace{2cm} r \leq r_0 \\ \frac{k}{2}(r-r_0)^2 & \hspace{2cm} r > r_0 \end{cases}\end{split}\]

where \(k\) is the force constant and \(r_0\) the radius of the sphere (centered at the origin) in which the potential is zero. The flat-bottom potential is activated by flags to the alexandria simulate command. It is useful mainly to keep molecules close to the origin and prevent them from flying into outer space.

3.1.4. Virtual Sites

A virtual site is an extra point, located at a defined position in a molecule. A variety of virtual site options is currently implemented within the ACT framework:

• a virtual site on top of an atom (VSITE1) [89]

• a virtual sites along the bond (VSITE2) for the description of anisotropic charge distribution and exchange [26] such as encountered in \(\sigma\) holes

• a virtual site on the bisector of a angle, like in the TIP4P water model [52] (VSITE3S) or in alcohol (VSITE3)

• off-plane virtual sites for modeling lone-pairs in \({\mathrm sp}^3\) hybridized compounds, such as water (VSITE3SOUT, symmetric) or asymmetric for compounds like alcohols [26, 90].

• Four-particle virtual sites (VSITE4) to model a lone-pair on an amine group.

3.1.5. Total energy

The total energy \(E\) of a compound then follows from

(3.48)\[E ~=~ V_{vdw}(r_{ij})+V_{coul}(r_{ij})+V_{pol}(r_{cs})+V_b(r_{ij})+V_a(\theta_{ijk})+V_i(\phi_{ijkl}) + V_{d}(\phi_{ijkl}).\]

Finally, it should be noted that the number of excluded neighbors is user-configurable. That means that atoms that are covalently bonded can interact both through the Buckingham (or Lennard-Jones) and Coulomb potentials, and through the bonded potentials. The main reason for this is that the short-range Coulomb interactions yield polarization anisotropy that is difficult to reproduce by a non-interacting model. To make sure that the forces on the atoms in a molecule are zero in the reference minimum-energy structure from quantum chemistry, both bond lengths \(r_{ij}^0\) and angles \(\theta_{ijk}^0\) can be treated as free parameters, that may differ substantially from the reference geometry. The number of exclusions can be selected separately for Coulomb and Van der Waals forces.

In total there are up to seven “atom” parameter types (\(\epsilon\), \(\sigma\), \(\gamma\), \(\delta\), \(\zeta\), \(q_s\), and \(\alpha\)) and 7-14 “bond” parameter types (\(k^b\), \(D^e\), \(r^0\), \(r^{max}\), \(k^{\theta}\), \(\theta^0\), \(k^{lin}\) and \(c_{n}\)) where \(n\) is the dihedral term index running from 0 to 6 to determine. All the atomic parameters are taken to be hybridization state dependent (corresponding to, for instance, sp ¹, sp ² and sp ³ carbon atoms). %It is straightforward to add support for other potentials, such as proper dihedrals.