This potential class is part of the third generation charge-optimized-many-body (COMB) potential [LDPS12]. The total potential energy consists of an electrostatic part \(E^{ES}\) and a short-ranged term \(E^{SR}\). Additional terms, such as van-der-Waals inteactions or angle correction terms, can be added later on.

As in the first and second generation COMB potentials, variable charges \(q_i\) are used. They are determined by minimizing \(E^{total}\) in each time step with respect to the charges. Additionally, induced dipoles \(\bf{\Delta}_i\) may be calculated for each particle, if polarization effects are allowed. The electrostatic part \(E^{ES}\) consists of four different terms: The self-energy term \(E^{self}\), the interactions between the variable charges \(E^{qq}\), the interactions between variable and fixed charges \(E^{qZ}\), and the polarization term \(E^{pol}\).

\[E^{ES} = E^{self} + E^{qq} + E^{qZ} + E^{pol}\]

These four terms are defined as follows:

\[E^{self} = \sum_i \left [ J_i^{(1)}(q_i -q_i^0) + J_i^{(2)} (q_i - q_i^0)^2 + J_i^{(3)} (q_i - q_i^0)^3 + J_i^{(4)} (q_i - q_i^0)^4 \right ]\]
\[\begin{split}E^{qq} = \sum_i \sum_{j>i} \frac{1}{4\pi \epsilon_0} q_i J_{ij}^{qq} q_j\end{split}\]
\[\begin{split}E^{qZ} = \sum_i \sum_{j>i} \frac{1}{4\pi \epsilon_0} q_i J_{ij}^{qZ} Z_j\end{split}\]
\[E^{pol} = \sum_i \left [\frac{1}{4\pi \epsilon_0} \frac{\Delta_i^T \Delta_i}{2\pi} + \frac{1}{4\pi \epsilon_0} \Delta_i^T \left[\sum_{j\neq i} q_j\frac{\partial J_{ij}^{qq}}{\partial r_{ij}}\frac{R_{ij}}{r_{ij}} \right ] + \frac{1}{2} \sum_{j\neq i} \Delta_i^T T_{ij} \Delta_j \right ]\]

All long-range interactions that occurr in \(E^{ES}\) are computed using the Wolf-summation [WKPE99]. The dipoles are calculated by minimizing \(E^{pol}\) with respect to the \(\Delta_i\), which corresponds to solving a set of linear equations.

The short-range interactions are given by

\[E^{SR} = E^{Rep} + E^{Attr}\]

which are defined by the following equations:

\[\begin{split}E^{Rep} = \sum_i \sum_{j>i} f_{ij}^S A_{ij}^{*} \exp(-\lambda_{ij} r_{ij})\end{split}\]
\[\begin{split}E^{Attr} = \sum_i \sum_{j>i}f_{ij}^S(b_{ij} + b_{ji})B_{ij}^{*}\sum_{k=1}^3\left[ B_{ij}^{(k)} \exp(-\mu_{ij}^{(k)} r_{ij}) \right]\end{split}\]
\[f_{ij}^S = f_C(r_{ij}, r_{ij}^{inner}, r_{ij}^{cut})\]
\[A_{ij}^{*} = A_{ij} \exp(0.5*\lambda_i D_i + 0.5 \lambda_j D_j)\]
\[B_{ij}^{*} = \exp(0.5\mu_i D_i + 0.5 \mu_j D_j) \sqrt{\left( a_i^B - | b_i^B(q_i-Q_i^0) |^{n_i^B}\right)\left( a_j^B - | b_j^B(q_j-Q_j^0) |^{n_j^B}\right)}\]
\[D_i = D_i^U + |b_i^D(Q_i^U - q_i)|^{n_i^D}\]
\[b_i^D = \frac{D_i^L - D_i^U)^{1/n_i^D}}{Q_i^U - Q_i^L}\]
\[n_i^D = \frac{\ln(-D_i^U/(D_i^U - D_i^L))}{\ln(Q_i^U/(Q_i^U - Q_i^L))}\]
\[b_{ij} = \left[1 + \left(\beta_i\sum_{k\neq i, k\neq j}\zeta_{ijk} \right)^n_i + P_{ij} \right]^{-1/(2n_i)}\]
\[\zeta_{ijk} = f_{ik}^S N_{ik} \exp(\beta_{ij}^{m_i}(r_{ij} - r_{ik})^{m_i})\sum_{l=0}^6 b_{ij}^{(l)} \cos(\theta_{ijk})^l\]
\[P_{ij} = c_{ij}^{(0)}\Omega_i + c_{ij}^{(1)} \exp(c_{ij}^{(2)}\Omega_i) + c_{ij}^{(3)}\]
\[\Omega_i = \sum_{j\neq i} f^S_{ij} N_{ij}\]
\[\Delta Q_i = 0.5 (Q^U_i - Q^L_i)\]
\[Q^O_i = 0.5(Q^U_i + Q^L_i)\]
\[a^B_i = \frac{1}{1 - |Q^O_i / \Delta Q_i|^{n^B_i}}\]
\[\begin{split}b^B_i &= \frac{|a^B_i|^{1/n^B_i}}{\Delta Q_i}\end{split}\]

A COMB3 potential is initiated by setting up a comb3particle object for each particle type. This will automatically activate all eletrostatic interactions that act on the given particle type. Note, that the constructor parameter p corresponds to the polarizability \(P_i\) in \(E^{pol}\).

The short-range interactions are activated by specifying the corresponding pair parameters in comb3pairpotential. Here, the parameters b0, b1, and b2 correspond to \(B_{ij}^0\), \(B_{ij}^1\), and \(B_{ij}^2\), respectively, while the constructor arguments di can be used to set \(b_{ij}^i\) (i=0,...,6). Note, that these pair interactions are non-symmetric, meaning that both particle-type combinations have to be specified, e.g. Titanium-Nitrogen and Nitrogen-Titanium in the above example.

One can also add a comb3fieldcorrection term \(E^{field}\) of the form

\[E^{field} = \frac{1}{4\pi\epsilon_0} \sum_i \sum_{i\neq j} f_C(r_{ij}, r_{ij}^{inner}, r_{ij}^{cut}) \left[\frac{P^\chi_{ij}q_j}{r_{ij}^3} + \frac{P^J q_j^2}{r_{ij}^5} \right] \, .\]

Similar to comb3pairpotential these interactions are non-symmetric and need to be specified for both particle-pair-combinations.

[LDPS12]Tao Liang, Bryce Devine, Simon R. Phillpot, and Susan B. Sinnott. Variable charge reactive potential for hydrocarbons to simulate organic-copper interactions. The Journal of Physical Chemistry A, 116(30):7976–7991, 2012. URL:, doi:10.1021/jp212083t.
[WKPE99]D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht. Exact method for the simulation of coulombic systems by spherically truncated pairwise r−1 summation. The Journal of Chemical Physics, 110(17):8254–8282, 1999. URL:, doi: