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')
potentialSet.addParticleType(ParticleType(
    symbol='N',
    mass=14.0067 * atomic_mass_unit,
    atomicNumber=7
))
potentialSet.addParticleType(ParticleType(
    symbol='Ti',
    mass=47.867 * atomic_mass_unit,
    atomicNumber=22
))

option = COMBCoulombOption(alpha = 0.2*1/Angstrom, r_cut = 11.0*Angstrom)
potentialSet.addOption(option)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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
)
potentialSet.addPotential(potential)
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.