MeamElementPotential

class MeamElementPotential(particleType, latticeType, nearestNeighbors, alpha, beta, referenceDistance, referenceEnergy, scalingFactor, weightingFactors, rho, gamma, attrac=0.0, repuls=0.0, nn2=False, zbl=True)

Constructor of the potential.

Parameters:
  • particleType (ParticleType or ParticleIdentifier) – Identifier of the particle type.
  • latticeType (str) – Type of the reference structure.
  • nearestNeighbors (int) – Number of nearest neighbors in the reference structure.
  • alpha (float) – Exponential decay factor for pair potential.
  • beta (sequence of float) – Exponential decay factors for the background electron densities. Must be of length four.
  • referenceDistance ([type]) – [description].
  • referenceEnergy (PhysicalQuantity of type energy) – The cohesive energy.
  • scalingFactor (float) – Scaling Factor in embedding function.
  • weightingFactors (sequence of float) – Weighting factors in the reference configuration. Must be of length three.
  • rho ([type]) – [description].
  • gamma (int) – Choice for gamma function.
  • attrac (float) – The attractive factor in the pair potential.
  • repuls (float) – The repulsive factor in the pair potential.
  • nn2 (bool) – Triggers second nearest neighbor formulation for symmetric pair.
  • zbl (bool) – Toggles zbl blending for small distances.
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.

setAlpha(alpha)

Set the parameter alpha.

Parameters:alpha (float) – Exponential decay factor for pair potential.
setAttrac(attrac)

Set the parameter attrac.

Parameters:attrac (float) – The attractive factor in the pair potential.
setBeta(beta)

Set the parameter beta.

Parameters:beta (sequence of float) – Exponential decay factors for the background electron densities. Must be of length four.
setGamma(gamma)

Set the parameter gamma.

Parameters:gamma (int) – Choice for gamma function.
setLatticeType(latticeType)

Set the reference lattice type.

Parameters:latticeType (str) – Type of the reference structure.
setNN2(nn2)

Set the parameter nn2.

Parameters:nn2 (bool) – Triggers second nearest neighbor formulation for symmetric pair.
setNearestNeighbors(nearestNeighbors)

Set the parameter nearestNeighbors.

:param nearestNeighbors:Number of nearest neighbors in the reference structure. :type nearestNeighbors:int

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.
setReferenceDistance(distance)

Set the parameter distance.

Parameters:distance (PhysicalQuantity of type length) – The equilibrium nearest neighbor distance.
setReferenceEnergy(energy)

Set the parameter energy.

Parameters:energy (PhysicalQuantity of type energy) – The cohesive energy.
setRepuls(repuls)

Set the parameter repuls.

Parameters:repuls (float) – The repulsive factor in the pair potential.
setRho(rho)

Set the parameter rho.

Parameters:rho (float) – Relative electron scaling factor.
setScalingFactor(scalingFactor)

Set the parameter scalingFactor.

Parameters:scalingFactor (float) – Scaling Factor in embedding function.
setWeightingFactors(factors)

Set the parameter factors.

Parameters:factors (sequence of float) – Weighting factors in the reference configuration. Must be of length three.
setZBL(zbl)

Set the ZBL parameter.

Parameters:zbl (bool) – Toggles zbl blending for small distances.

Usage Examples

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

# Set up lattice
vector_a = [5.4306, 0.0, 0.0]*Angstrom
vector_b = [0.0, 5.4306, 0.0]*Angstrom
vector_c = [0.0, 0.0, 5.4306]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

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

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

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

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
potentialSet = TremoloXPotentialSet(name='MEAM_Si_2007')
potentialSet.addParticleType(ParticleType(
    symbol='Si',
    mass=28.0855*atomic_mass_unit,
    charge=None,
    sigma=None,
    sigma14=None,
    epsilon=None,
    epsilon14=None,
    atomicNumber=14,
    tags=[],
))


_option = MeamGlobalOption(
    delr=0.1*Angstrom,
    erose=2,
    wf_mixing=2,
    r_cut=4.5*Angstrom,
    augment_1st=0.0,
    embedding_negative=False,
    density_scaling=False,
)
potentialSet.addOption(_option)
_potential = MeamElementPotential(
    particleType=ParticleIdentifier('Si', []),
    latticeType='dia',
    nearestNeighbors=4,
    alpha=4.89890486934,
    beta=[3.55, 2.5, 0.0, 7.5, ],
    referenceDistance=2.35*Angstrom,
    referenceEnergy=4.63*eV,
    scalingFactor=0.58,
    weightingFactors=[1.8, 5.25, -2.61, ],
    rho=1.0,
    gamma=3,
    attrac=0.0,
    repuls=0.0,
    nn2=True,
    zbl=False,
)
potentialSet.addPotential(_potential)
_potential = MeamScreeningPotential(
    particleType1=ParticleIdentifier('Si', []),
    particleType2=ParticleIdentifier('Si', []),
    particleType3=ParticleIdentifier('Si', []),
    Cmin=1.41,
    Cmax=2.8,
)
potentialSet.addPotential(_potential)

calculator = TremoloXCalculator(parameters=potentialSet)
calculator.setVerletListsDelta(0.25*Angstrom)

bulk_configuration.setCalculator(calculator)
bulk_configuration.update()

Notes

The modified embedded-atom-method (MEAM) potential extends the conventional embedded-atom-method by including directional bonding [Bas92]. This allows the application to both metallic and covalent systems, as well as mixed systems.

The energy function for the Modified Embedded Atom Method as presented e.g. in [GWS03] is given as the sum of some pair term and some embedding term as follows:

\[E= \sum_i \left\{ F_i(\bar{\rho}_i) + \frac{1}{2} \sum_{j\ne i} \Phi_{ij}(r_{ij})\right\}\]

The embedding function is of the form

\[F_i(\bar{\rho}_i) = A_i E_i^0\bar{\rho}_i \log \bar{\rho}_i\]

for \(\bar{\rho}_i > 0\) and for \(\bar{\rho}_i \le 0\)

\[F_i(\bar{\rho}_i) = 0\]

if embedding_negative is 0 or

\[F_i(\bar{\rho}_i) = -A_i E_i^0 \bar{\rho}_i\]

if it is 1, where \(A_i\) and \(E_i^0\) are element dependent parameters. The background electron density \(\bar{\rho}_i\) is given by:

\[\bar{\rho}_i = \frac{\bar{\rho}_i^{(0)}}{\rho_i^0} G_i (\Gamma_i)\]

where

\[\Gamma_i = \sum_{k=1}^3 t_i^{(k)}\left(\frac{\bar{\rho}_i^{(k)}}{\bar{\rho}_i^{(0)}}\right)^2\]

and \(G_i\) is one of

(1)\[\begin{split}G_i(\Gamma) = \begin{cases} \sqrt{1+\Gamma} & gamma=0\\ \exp(\Gamma / 2) & gamma=1\\ \text{sign}(1+\Gamma) \sqrt{|1+\Gamma|} & gamma=2\\ \frac{2}{1+\exp(-\Gamma)} &gamma=3\\ \sqrt{1+\Gamma} & gamma=4 \end{cases}\end{split}\]

With the element dependent constant \(\rho_i^0\), the number of nearest neighbors in the reference structure \(Z_{i0}\) and \(\Gamma_i^\text{ref}\) \(\Gamma\) evaluated in the reference structure, the composition dependent density scaling can be chosen to be either

(2)\[\rho_i^0 = \rho_{i0}Z_{i0}G_i(\Gamma_i^\text{ref})\]

or

(3)\[\rho_i^0 = \rho_{i0} Z_{i0}\]

Note that for gamma=0 and gamma=2 \(G_i(\Gamma_i^\text{ref})\) is always set to 1. The partial electron densities are defined as

\[\begin{split}\begin{eqnarray*} \bar{\rho}_i^{(0)} &=& \sum_{j\ne i} \rho_j^{a(0)}(r_{ij})S_{ij} \\ (\bar{\rho}_i^{(1)})^2 &=& \sum_{\alpha = 1}^3 \left[\sum_{j \ne i} \rho_j^{a(1)} \frac{r_{ij\alpha}}{r_{ij}} S_{ij} \right]^2 \\ (\bar{\rho}_i^{(2)})^2 &=& \sum_{\alpha =1}^3 \sum_{\beta =1}^3 \left[\sum_{j\ne i} \rho_j^{a(2)} \frac{r_{ij\alpha}r_{ij\beta}}{r_{ij}^2} S_{ij} \right]^2 - \frac{1}{3} \left[ \sum_{j\ne i} \rho_j^{a(2)}S_{ij}\right]^2\\ (\bar{\rho}_i^{(3)})^2 &=& \sum_{\alpha =1}^3 \sum_{\beta =1}^3 \sum_{\gamma =1}^3 \left[\sum_{j\ne i} \rho_j^{a(3)} \frac{r_{ij\alpha}r_{ij\beta}r_{ij\gamma}}{r_{ij}^3} S_{ij}\right]^2 - \frac{3}{5}\sum_{\alpha =1}^3\left[\sum_{j \ne i} \rho_j^{a(3)} \frac{r_{ij\alpha}}{r_{ij}} S_{ij}\right]^2 \end{eqnarray*}\end{split}\]

with the atomic electron densities given by

\[\rho_i^{a(k)}(r_{ij}) = \rho_{i0} \exp \left[-\beta_i^{(k)} \left(\frac{r_{ij}}{r_i^0} - 1 \right)\right]\]

Finally, the average weighting factors are given by

(4)\[\begin{split}t_i^{(k)} = \begin{cases} \frac{1}{\bar{\rho}_i^{(0)}} \sum_{j\ne i} t_{0,j}^{(k)} \rho_j^{a(0)} S_{ij} & \text{wf-mixing=0} \\ \frac{ \sum_{j\ne i} t_{0,j}^{(k)} \rho_j^{a(0)} S_{ij}} { \sum_{j\ne i} \left(t_{0,j}^{(k)}\right)^2 \rho_j^{a(0)} S_{ij}}& \text{wf-mixing=1} \\ t_{0,i}^{(k)} & \text{wf-mixing=2} \end{cases}\end{split}\]

with element dependent parameters \(t_{0,i}^{(k)}\).

The pair term \(\Phi_{ij}\) is defined as

(5)\[\begin{split}\begin{eqnarray} \Phi_{ij}(r_{ij}) &=& \bar{\Phi}_{ij}(r_{ij}) S_{ij}\\ \bar{\Phi}_{ij}(r_{ij}) &=& \frac{1}{Z_{ij0}} [2E_i^u(r_{ij}) - F_i(\hat{\rho}_i (r_{ij})) - F_j(\hat{\rho}_j(r_{ij})] \end{eqnarray}\end{split}\]

where \(Z_{ij0}\) is the number of nearest neighbors in the reference structure, \(F_i(\hat{\rho}_{ij})\) the embedding function evaluated in the reference structure and

(6)\[\begin{split}E_i^u(r_{ij}) = \begin{cases} -E_{ij}^c \left(1+a_{ij}^{*} + a_{ij}^{(3)} {a_{ij}^{*}}^3 \frac{r^0_{ij}}{r_{ij}}\right) \exp(-a_{ij}^*) & \text{erose=0} \\ -E_{ij}^c \left(1+a_{ij}^* + \left(-\text{attrac}+\frac{\text{repuls}}{r_{ij}}\right){a_{ij}^*}^3\right)\exp(-a_{ij}^*) & \text{erose=1}\\ -E_{ij}^c\left(1+a_{ij}^* + a_{ij}^{(3)}{a_{ij}^*}^3\right)\exp(-a_{ij}^*) & \text{erose=2} \end{cases}\end{split}\]

with

\[a_{ij}^*(r_{ij}) = \alpha_{ij} \left(\frac{r_{ij}}{r^0_{ij}} - 1\right)\]

and

\[\begin{split}a_{ij}^{(3)} = \begin{cases} \text{repuls} & a_{ij}^* < 0\\ \text{attrac} & a_{ij}^* \ge 0 \end{cases}\end{split}\]

Some reference structures use a more complex form of (5). For small \(r_{ij}\), this pair term is blended with a ZBL potential by default.newline Last, the screening function \(S_{ij}\) is defined by

\[\begin{split}\begin{eqnarray*} S_{ij}&=&\bar{S}_{ij} f_c\left(\frac{r_c - r_{ij}}{\Delta r}\right)\\ \bar{S}_{ij} &=& \prod_{k\ne i,j} S_{ikj}\\ S_{ikj} &=& f_c\left(\frac{C_{ikj}-C_{\text{min},ikj}}{C_{\text{max},ikj} -C_{\text{min},ikj}}\right)\\ C_{ikj}&=& 1+2\frac{r_{ij}^2 r_{ik}^2 + r_{ij}^2r_{jk}^2-r_{ij}^4}{r_{ij}^4-(r_{ik}^2-r_{jk}^2)^2}\\ f_c(x)&=&\begin{cases} 1 & x \ge 1 \\ [1-(1-x)^4]^2 & 0<x<1\\ 0 & x \le 0 \end{cases} \end{eqnarray*}\end{split}\]

where \(r_c\) is the global cutoff radius, \(\Delta r\) is the length of a smoothing region for \(r\) near \(r_c\) and \(C_\text{min}\) and \(C_\text{max}\) are element-triple dependent parameters.

Note that the cutoff radius for this term needs to be larger than \(r_c\), depending on \(C_{\text{max}}\), because even triples with some particle distance \(> r_c\) may have a screening effect.

The parameter lattice defines the reference structure for the given particle type. The parameters alpha, beta_k, re, Ec, scaling_factor, weighting_factor_k, rho, gamma, attrac and repuls correspond to the variables \(\alpha\), \(\beta^{(k)}\), \(r^0\), \(E^c\), \(A\), \(t_0^{(k)}\), \(\rho_{0}\), attrac and repuls from the previous section respectively. nn2 enables the second nearest neighbor formulation as described in cite{lee2000second}, zbl enables ZBL blending for small distances if set to 1.

re and Ec correspond to \(r^0\) and \(E^c\). Note that all element_data entries must be given before the first element_pair entry.newline If augment_1st is set to 1, we set

\[t_{0,i}^{(1)} = t_{0,i}^{(1)} + \frac{3}{5} t_{0,i}^{(3)}\]

for all \(i\). This is usually only needed for older parameter sets, modern ones already include this correction. The default density_scaling uses (2), setting it to 1 uses (3) instead.

The latticeType parameter must be one of the following strings:

  • dia = diamond
  • fcc = face centered cubic
  • bcc = body centered cubic
  • dim = dimer
  • b1 = rock salt
  • hcp = hexagonal close-packed
  • c11 = MoSi2 structure
  • l12 = Cu3Au structure
  • b2 = CsCl structure

[Bas92]M. I. Baskes. Modified embedded-atom potentials for cubic materials and impurities. Phys. Rev. B, 46:2727–2742, Aug 1992. URL: http://link.aps.org/doi/10.1103/PhysRevB.46.2727, doi:10.1103/PhysRevB.46.2727.
[GWS03]PM Gullet, G Wagner, and A Slepoy. Numerical tools for atomistic simulations. SANDIA Report, 8782:2003, 2003.