# COMB3FieldCorrection¶

class COMB3FieldCorrection(particleType1, particleType2, Pchi, PJ, r_i=None, r_cut=None)

Constructor of the potential.

Parameters: particleType1 (ParticleType or ParticleIdentifier) – Identifier of the first particle type. particleType2 (ParticleType or ParticleIdentifier) – Identifier of the second particle type. Pchi (PhysicalQuantity of length**2 * charge) – Potential parameter. PJ (PhysicalQuantity of length**4) – Potential parameter. r_i (PhysicalQuantity of type length) – Inner cutoff radius (distance where the smoothing starts) r_cut (PhysicalQuantity of type length) – The cutoff radius of this potential.
getAllParameterNames()

Return the names of all used parameters as a list.

getAllParameters()

Return all parameters of this potential and their current values as a <parameterName / parameterValue> dictionary.

static getDefaults()

Get the default parameters of this potential and return them in form of a dictionary of <parameter name, default value> key-value pairs.

getParameter(parameterName)

Get the current value of the parameter parameterName.

setCutoff(r_cut)

Set the cutoff radius for this potential.

Parameters: r_cut (PhysicalQuantity of type length) – The cutoff radius of this potential.
setPJ(PJ)

Set the parameter PJ.

Parameters: PJ (PhysicalQuantity of length**4) – Potential parameter.
setParameter(parameterName, value)

Set the parameter parameterName to the given value.

Parameters: parameterName (str) – The name of the parameter that will be modified. value – The new value that will be assigned to the parameter parameterName.
setPchi(Pchi)

Set the parameter Pchi.

Parameters: Pchi (PhysicalQuantity of length**2 * charge) – Potential parameter.

## Usage Examples¶

Define a COMB3 potential for a TiN rocksalt crystal by adding particle types and interaction functions to the TremoloXPotentialSet.

# Set up a Titanium-Nitride cell

vector_a = [4.235, 0.0, 0.0]*Angstrom
vector_b = [0.0, 4.235, 0.0]*Angstrom
vector_c = [0.0, 0.0, 4.235]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
elements = [Titanium, Nitrogen, Titanium, Nitrogen, Titanium, Nitrogen,
Titanium, Nitrogen]

# Define coordinates
fractional_coordinates = [[ 0. ,  0. ,  0. ],
[ 0.5,  0.5,  0.5],
[ 0.5,  0.5,  0. ],
[ 0. ,  0. ,  0.5],
[ 0.5,  0. ,  0.5],
[ 0. ,  0.5,  0. ],
[ 0. ,  0.5,  0.5],
[ 0.5,  0. ,  0. ]]

# Set up configuration
bulk_configuration = BulkConfiguration(
bravais_lattice=lattice,
elements=elements,
fractional_coordinates=fractional_coordinates
)

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------

potentialSet = TremoloXPotentialSet(name = 'COMB_NTi_2014')
symbol='N',
mass=14.0067 * atomic_mass_unit,
atomicNumber=7
))
symbol='Ti',
mass=47.867 * atomic_mass_unit,
atomicNumber=22
))

option = COMBCoulombOption(alpha = 0.2*1/Angstrom, r_cut = 11.0*Angstrom)
potential = COMB3Particle(
particleType = 'Ti',
l = 2.136011*1/Angstrom,
mu = 1.178831*1/Angstrom,
m = 1,
n = 0.566048,
QL = -4.0*elementary_charge,
QU = 4.0*elementary_charge,
DL = 0.005*Angstrom,
DU = -0.5*Angstrom,
nB = 10,
J1 = 3.095768*eV/elementary_charge,
J2 = 4.23028*eV/elementary_charge**2,
J3 = -1.039759*eV/elementary_charge**3,
J4 = 0.357428*eV/elementary_charge**4,
q0 = 0.0*elementary_charge,
xi = 0.724352*1/Angstrom,
Z = 3.022932*elementary_charge,
P = 0.335*Angstrom**3,
uniformCharge = False
)
potential = COMB3Particle(
particleType = 'N',
l = 5.218037*1/Angstrom,
mu = 3.738549*1/Angstrom,
m = 1,
n = 1.0,
QL = -2.0*elementary_charge,
QU = 6.0*elementary_charge,
DL = 0.007664*Angstrom,
DU = -1.213951*Angstrom,
nB = 10,
J1 = 6.59963*eV/elementary_charge,
J2 = 5.955097*eV/elementary_charge**2,
J3 = 0.760433*eV/elementary_charge**3,
J4 = 0.009388*eV/elementary_charge**4,
q0 = 0.0*elementary_charge,
xi = 1.371794*1/Angstrom,
Z = -1.53917*elementary_charge,
P = 0.5167863369*Angstrom**3,
uniformCharge = False
)
potential = COMB3PairPotential(
particleType1 = 'Ti',
particleType2 = 'Ti',
A = 516.5873248*eV,
B0 = 117.0421345*eV,
B1 = 0.0*eV,
B2 = 0.0*eV,
l = 2.136011*1/Angstrom,
mu0 = 1.178831*1/Angstrom,
mu1 = 0.0*1/Angstrom,
mu2 = 0.0*1/Angstrom,
beta = 0.545629*1/Angstrom,
b0 = 0.077183,
b1 = 0.146606,
b2 = 0.226674,
b3 = 0.093054,
b4 = -0.130433,
b5 = -0.065549,
b6 = 0.5504,
c0 = -0.012318,
c1 = 0.175079,
c2 = 0.066085,
c3 = -0.076109,
N = 1.0,
r_i = 3.9*Angstrom,
r_cut = 4.1*Angstrom
)
potential = COMB3PairPotential(
particleType1 = 'N',
particleType2 = 'N',
A = 7654.972656*eV,
B0 = 2102.295654*eV,
B1 = 0.0*eV,
B2 = 0.0*eV,
l = 5.218037*1/Angstrom,
mu0 = 3.738549*1/Angstrom,
mu1 = 0.0*1/Angstrom,
mu2 = 0.0*1/Angstrom,
beta = 3.738549*1/Angstrom,
b0 = 1.691325,
b1 = 0.283832,
b2 = 0.698596,
b3 = 1.724015,
b4 = 0.812665,
b5 = 0.0,
b6 = 0.0,
c0 = 0.0,
c1 = 0.0,
c2 = 0.0,
c3 = 0.0,
N = 1.0,
r_i = 2.0*Angstrom,
r_cut = 2.3*Angstrom
)
potential = COMB3PairPotential(
particleType1 = 'Ti',
particleType2 = 'N',
A = 1661.765935*eV,
B0 = 104.805547*eV,
B1 = 135.992262*eV,
B2 = 30.112007*eV,
l = 3.17815*1/Angstrom,
mu0 = 2.453223*1/Angstrom,
mu1 = 1.766886*1/Angstrom,
mu2 = 2.439896*1/Angstrom,
beta = 2.297286*1/Angstrom,
b0 = 0.225186,
b1 = 0.269113,
b2 = 0.261613,
b3 = 0.033596,
b4 = -0.113159,
b5 = 2.2e-05,
b6 = 3e-06,
c0 = -0.032039,
c1 = -0.082354,
c2 = 0.024479,
c3 = -0.073187,
N = 1.0,
r_i = 2.6*Angstrom,
r_cut = 3.0*Angstrom
)
potential = COMB3PairPotential(
particleType1 = 'N',
particleType2 = 'Ti',
A = 1661.765935*eV,
B0 = 104.805547*eV,
B1 = 135.992262*eV,
B2 = 30.112007*eV,
l = 3.17815*1/Angstrom,
mu0 = 2.453223*1/Angstrom,
mu1 = 1.766886*1/Angstrom,
mu2 = 2.439896*1/Angstrom,
beta = 5.066036*1/Angstrom,
b0 = 0.024877,
b1 = -0.103268,
b2 = 0.263335,
b3 = 0.412664,
b4 = 0.017287,
b5 = 2.3e-05,
b6 = 0.0,
c0 = -0.056517,
c1 = -0.036331,
c2 = 0.007137,
c3 = 0.06462,
N = 1.0,
r_i = 2.6*Angstrom,
r_cut = 3.0*Angstrom
)
potential = COMB3FieldCorrection(
particleType1 = 'Ti',
particleType2 = 'Ti',
Pchi = 0.0253919125995*Angstrom**2*elementary_charge,
PJ = 0.0531158828648*Angstrom**4,
r_i = 10.0*Angstrom,
r_cut = 11.0*Angstrom
)
potential = COMB3FieldCorrection(
particleType1 = 'Ti',
particleType2 = 'N',
Pchi = 0.0253919125995*Angstrom**2*elementary_charge,
PJ = 0.0531158828648*Angstrom**4,
r_i = 10.0*Angstrom,
r_cut = 11.0*Angstrom
)
potential = COMB3FieldCorrection(
particleType1 = 'N',
particleType2 = 'Ti',
Pchi = 0.0728019460398*Angstrom**2*elementary_charge,
PJ = 0.271659870233*Angstrom**4,
r_i = 10.0*Angstrom,
r_cut = 11.0*Angstrom
)
potential = COMB3FieldCorrection(
particleType1 = 'N',
particleType2 = 'N',
Pchi = 0.0728019460398*Angstrom**2*elementary_charge,
PJ = 0.271659870233*Angstrom**4,
r_i = 10.0*Angstrom,
r_cut = 11.0*Angstrom
)
potential = AngleCorrection6Potential(
particleType1 = 'Ti',
particleType2 = 'Ti',
particleType3 = 'Ti',
k0 = 0.0*eV,
k1 = 0.043561*eV,
k2 = 0.0*eV,
k3 = 0.049195*eV,
k4 = 0.0*eV,
k5 = 0.0*eV,
k6 = 0.0*eV,
costheta0 = 0.0,
r_i = 3.9*Angstrom,
r_cut = 4.1*Angstrom
)
calculator = TremoloXCalculator(parameters=potentialSet)
calculator.setInternalOrdering("default")

bulk_configuration.setCalculator(calculator)


### Notes¶

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: http://dx.doi.org/10.1021/jp212083t, 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: http://scitation.aip.org/content/aip/journal/jcp/110/17/10.1063/1.478738, doi:http://dx.doi.org/10.1063/1.478738.