COMBSMCoulomb

class COMBSMCoulomb(particleType1, particleType2, xii, xij)

Constructor of the potential.

Note: COMBParticle must be added for BOTH particle types!

Parameters:
  • particleType1 (ParticleType or ParticleIdentifier) – Identifier of the first particle type.
  • particleType2 (ParticleType or ParticleIdentifier) – Identifier of the second particle type.
  • xii (PhysicalQuantity of type length**-1) – Potential parameter.
  • xij (PhysicalQuantity of type length**-1) – Potential parameter.
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.
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.
setXii(xii)

Set the parameter xii.

Parameters:xii (PhysicalQuantity of type length**-1) – Potential parameter.
setXij(xij)

Set the parameter xij.

Parameters:xij (PhysicalQuantity of type length**-1) – Potential parameter.

Usage Examples

Define a COMB potential for a SiO2 quartz crystal by adding particle types and interaction functions to the TremoloXPotentialSet.

# -------------------------------------------------------------
# Set up a SiO2 Quartz crystal
# -------------------------------------------------------------

# Set up lattice
lattice = Hexagonal(4.916*Angstrom, 5.4054*Angstrom)

# Define elements
elements = [Silicon, Silicon, Silicon, Oxygen, Oxygen, Oxygen, Oxygen, Oxygen,
            Oxygen]

# Define coordinates
fractional_coordinates = [[ 0.4697,  0.0000,  0.0000    ],
                          [ 0.0000,  0.4697,  0.66666667],
                          [ 0.5303,  0.5303,  0.33333333],
                          [ 0.4135,  0.2669,  0.1191    ],
                          [ 0.2669,  0.4135,  0.547567  ],
                          [ 0.7331,  0.1466,  0.785767  ],
                          [ 0.5865,  0.8534,  0.214233  ],
                          [ 0.8534,  0.5865,  0.452433  ],
                          [ 0.1466,  0.7331,  0.8809    ]]

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

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

potentialSet = TremoloXPotentialSet(name='COMB_SiO')

# Add particle types for Si and O.
potentialSet.addParticleType(ParticleType(symbol='Si',
                                          mass=28.0855*atomic_mass_unit,
                                          atomicNumber=14))
potentialSet.addParticleType(ParticleType(symbol='O',
                                          mass=15.9994*atomic_mass_unit,
                                          atomicNumber=8))

# Add COMB potential for Si.
potential = COMBSingleTypePotential(particleType='Si',
                                    A = 1830.81*eV,
                                    B = 471.17*eV,
                                    R = 2.6*Angstrom,
                                    S = 3.0*Angstrom,
                                    l = 2.4799*1/Angstrom,
                                    mu = 1.7322*1/Angstrom,
                                    beta = 1.0999e-06,
                                    m = 3,
                                    n = 0.78734,
                                    c = 100390.0,
                                    d = 16.218,
                                    h = -0.59826,
                                    QL = -4.0,
                                    QU = 4.0,
                                    DL = 1.7982*Angstrom,
                                    DU = -1.7094*Angstrom,
                                    nB = 10,
                                    J1 = -0.3*eV/elementary_charge,
                                    J2 = -2.2*eV/elementary_charge**2,
                                    J3 = -1.25*eV/elementary_charge**3,
                                    J4 = 2.53*eV/elementary_charge**4,
                                    q0 = 0.6112*elementary_charge,
                                    uniformCharge = False)
potentialSet.addPotential(potential)

# Add COMB potential for O.
potential = COMBSingleTypePotential(particleType = 'O',
                                    A = 3326.7*eV,
                                    B = 260.89*eV,
                                    R = 2.6*Angstrom,
                                    S = 3.0*Angstrom,
                                    l = 5.36*1/Angstrom,
                                    mu = 2.68*1/Angstrom,
                                    beta = 2.0,
                                    m = 1,
                                    n = 1.0,
                                    c = 6.6,
                                    d = 1,
                                    h = -0.229,
                                    QL = -1.8349,
                                    QU = 5.5046,
                                    DL = 0.00148*Angstrom,
                                    DU = -0.00112*Angstrom,
                                    nB = 10,
                                    J1 = 12.006*eV/elementary_charge,
                                    J2 = 10.5205*eV/elementary_charge**2,
                                    J3 = 0.0*eV/elementary_charge**3,
                                    J4 = 0.0*eV/elementary_charge**4,
                                    q0 = 0.0*elementary_charge,
                                    uniformCharge = False)
potentialSet.addPotential(potential)

# Combination rules between Si and O.
potential = COMBMixitPotential(particleType1 = 'O',
                               particleType2 = 'Si',
                               K = 20.0,
                               r0 = 1.65*Angstrom)
potentialSet.addPotential(potential)

# Three-body-terms.
potential = COMBTriplePotential(particleType1 = 'Si',
                                particleType2 = 'O',
                                particleType3 = 'O',
                                cosTheta = -0.333313248,
                                K0 = 3.16666666667*eV,
                                K1 = 0.0*eV,
                                K2 = 6.33333333333*eV,
                                K3 = 0.0*eV,
                                K4 = 0.0*eV,
                                K5 = 0.0*eV,
                                K6 = 0.0*eV)
potentialSet.addPotential(potential)
potential = COMBTriplePotential(particleType1 = 'O',
                                particleType2 = 'Si',
                                particleType3 = 'Si',
                                cosTheta = -0.806238149,
                                K0 = 1.5*eV,
                                K1 = 0.0*eV,
                                K2 = 3.0*eV,
                                K3 = 0.0*eV,
                                K4 = 0.0*eV,
                                K5 = 0.0*eV,
                                K6 = 0.0*eV)
potentialSet.addPotential(potential)

# Coulomb interactions.
potential = COMBCoulomb(particleType1 = 'Si',
                        particleType2 = 'Si',
                        R = 1.1*Angstrom,
                        S = 5.45*Angstrom,
                        eta = 0.380689)
potentialSet.addPotential(potential)
potential = COMBCoulomb(particleType1 = 'Si',
                        particleType2 = 'O',
                        R = 1.1*Angstrom,
                        S = 5.45*Angstrom,
                        eta = 0.2817839)
potentialSet.addPotential(potential)
potential = COMBCoulomb(particleType1 = 'O',
                        particleType2 = 'O',
                        R = 1.1*Angstrom,
                        S = 5.45*Angstrom,
                        eta = 0.20857489)
potentialSet.addPotential(potential)
calculator = TremoloXCalculator(parameters=potentialSet)

bulk_configuration.setCalculator(calculator)

Notes

The charge-optimized multibody (COMB) potential is an extended version of the Tersoff potential (cf. TersoffSingleTypePotential), which additionally allows to take interatomic charge transfer into account [Yas96], [YSP07].

The total energy is composed of the following contributions.

\[E = \sum_i \phi_i + \frac{1}{2} \sum_{j \neq i} \left [ U_{ij}^\mathrm{rep} + U_{ij}^\mathrm{sht} + U_{ij}^\mathrm{ion} + U_{ij}^\mathrm{corr} \right ] + \frac{1}{2} \sum_{j \neq i} \sum_{k\neq i\atop k\neq j} U^\mathrm{triple}_{ijk}\]

Using the smoothed cutoff function \(f_C\) from TersoffSingleTypePotential we can write the potential components as follows:

\[\phi_i = J1_i(q_i - q_i^0) + \frac{1}{2} J2_i(q_i - q_i^0)^2 + J3_i (q_i - q_i^0)^3 + J4_i (q_i - q_i^0)^4 \, ,\]
\[U^\mathrm{rep}_{ij} = f^S_{ij} A_{ij} \exp(-\lambda_{ij} r_{ij}) \left(1 + K_{ij}(1 - r_{ij} / r_{ij}^0)^2\right) \, ,\]
\[U^\mathrm{sht}_{ij} = -f^S_{ij} b_{ij} B_{ij} \exp(-\mu_{ij} r_{ij}) \, ,\]
\[U^\mathrm{ion}_{ij} = \frac{f^L_{ij} \eta_i\eta_j q_i q_j}{4\pi\varepsilon_0 r_{ij}} \, ,\]
\[f^S_{ij} = f_C(r_{ij}, R_{ij}, S_{ij}) \, ,\]
\[A^S_i = A_i \exp(\lambda_i D_i) \, ,\]
\[B^S_i = B_i \exp(\mu_i D_i) \left(a^B_i - |b^B_i(q_i - Q^O_i)|^{n^B_i} \right) \, ,\]
\[D_i = D^U_i + |b^D_i(Q^U_i - q_i)|^{n^D_i} \, ,\]
\[b^D_i = \frac{(D^L_i - D^U_i)^{1 / n^D_i}}{Q^U_i - Q^L_i} \, ,\]
\[n^D_i = \frac{\ln(-D^U_i / (D^U_i - D^L_i))}{\ln(Q^U_i / (Q^U_i - Q^L_i)} \, ,\]
\[b_{ij} = \left( 1 + \left(\beta_i \sum_{k\neq i\atop k\neq j} \zeta_{ijk}\right)^{n_i} \right)^{-1/(2n_i)} \, ,\]
\[\zeta_{ijk} = f^S_{ik} \exp(\mu_{ij}^{m_i} (r_{ij} - r_{ik})^{m_i}) \left(1 + \frac{c_i^2}{d_i^2} - \frac{c_i^2}{d_i^2 + (h_i - \cos \theta_{i,j,k})^2} \right) \, ,\]
\[\Delta Q_i = \frac{1}{2} (Q^U_i - Q^L_i) \, ,\]
\[Q^O_i = \frac{1}{2} (Q^U_i + Q^L_i) \, ,\]
\[a^B_i = \frac{1}{1 - |Q^O_i / \Delta Q_i|^{n^B_i}} \, ,\]
\[b^B_i = \frac{|a^B_i|^{1/n^B_i}}{\Delta Q_i} \, .\]

These parameters are specified in COMBParticle, or, for interactions between different particle types, in COMBPairPotential.

The charges \(q_i\) are not static. Instead, they are recalculated in each simulation step by minimizing \(E\) with respect to the charges.

There are several additions to the basic COMB potential. One can for example add additional bond bending terms (cf. COMBTriplePotential):

\[U^\mathrm{triple}_{ijk} = f^S_{ij} f^S_{ik} \sum_{d = 0}^6 K_{ijk}^d L^d\left(\cos \theta_{ijk} - \cos \theta_{ijk}^0\right) ,\]

where \(L^d\) denotes the Legendre polynomial of degree \(d\).

In [DSC+11], a self-energy correction term was introduced (cf. COMBSelfEnergyCorrection):

\[U^\mathrm{corr}_{ij} = \frac{q_i^2}{4\pi\varepsilon_0} \left( \frac{a_i q_j}{r_{ij}^3} + \frac{b_i q_j^2}{r_{ij}^5}\right) .\]

Alternatively, the modification proposed in [SDK+10],

\[U^\mathrm{corr}_{ij} = \frac{1}{4\pi\varepsilon_0} \left( \frac{a_i q_j}{r_{ij}^5} + \frac{b_i q_j^2}{r_{ij}^5}\right) ,\]

can be used by specifying COMBSelfEnergyCorrectionShan.

To add electrostatic interactions, you either use damped point-charge interactions,

\[U^\mathrm{ion}_{ij} = \frac{f^S_{ij} \eta_{ij} q_i q_j}{4\pi\varepsilon_0 r_{ij}} \, ,\]

with parameters specified in COMBPointWiseCoulomb, or resort to the Streitz-Mintmire model (cf. COMBSMCoulomb),

\[U^\mathrm{ion}_{ij} = \frac{1}{4\pi\varepsilon_0 r_{ij}} J_{ij} q_i q_j \, ,\]

where \(J_{ij}\) is calculated by integration over the charge densities \(\rho_i\) [DSC+11],

\[\rho_i(x) = \frac{\xi^3}{\pi} \exp(-2\xi_i |x - r_i|) .\]

Instead of specifying all interaction parameters between different particle types manually via COMBPairPotential, these are often determined by employing the following combination rules:

\[A_{ij} = \sqrt{A_i A_j} \, ,\]
\[B_{ij} = \sqrt{B_i B_j} \, ,\]
\[\lambda_{ij} = 0.5(\lambda_i + \lambda_j) \, ,\]
\[\mu_{ij} = 0.5(\mu_i + \mu_j) \, .\]

These combination rules are enabled by using the COMBMixitPotential for the desired particle type combination.

[DSC+11](1, 2) B. Devine, T.-R. Shan, Y.-T. Cheng, A. J. H. McGaughey, M. Lee, S. R. Phillpot, and S. B. Sinnott. Atomistic simulations of copper oxidation and cu/cu2o interfaces using charge-optimized many-body potentials. Phys. Rev. B, 84:125308, Sep 2011. doi:10.1103/PhysRevB.84.125308.
[SDK+10]T.-R. Shan, B. D. Devine, T. W. Kemper, S. B. Sinnott, and S. R. Phillpot. Charge-optimized many-body potential for the hafnium/hafnium oxide system. Phys. Rev. B, 81:125328, Mar 2010. doi:10.1103/PhysRevB.81.125328.
[Yas96]Akio Yasukawa. Using An Extended Tersoff Interatomic Potential to Analyze The Static-Fatigue Strength of SiO2 under Atmospheric Influence. JSME international journal. Ser. A, Mechanics and material engineering, 39(3):313–320, jul 1996. URL: http://ci.nii.ac.jp/naid/110002964467/en/.
[YSP07]J. Yu, S. B. Sinnott, and S. R. Phillpot. Charge optimized many-body potential for the si∕sio2 system. Phys. Rev. B, 75:085311, Feb 2007. doi:10.1103/PhysRevB.75.085311.