Many-body potentials

matscipy implements support for generic manybody potentials. The class matscipy.calculators.manybody.Manybody implements the functional form

\[\begin{split} U = \frac{1}{2} \sum_{\substack{ij\\ i\neq j}} U_2(r^2_{ij}) + U_\text{m}(r^2_{ij}, \xi_{ij})\end{split}\]

with

\[\begin{split} \xi_{ij} = \sum_{\substack{k\\ k\neq i,j}} \Xi(r^2_{ij}, r^2_{ik}, r^2_{jk})\end{split}\]

as described, e.g. by Müser et al. and Grießer et al.. On top of energies and forces, the calculator can compute second derivatives (with respect to positions and strain degrees of freedom). Explicit functional forms of \(U_2\), \(U_\text{m}\) and \(\Xi\) are implemented for a number of potential.

Kumagai potential

The following example computes the elastic constants of a small representation of amorphous silicon using the potential by Kumagai et al.. We first load the amorphous structure.

from ase.io import read

def interactive_view(system):
    from ase.visualize import view
    # Interactive view of the lattice
    v = view(system, viewer='ngl')

    # Resize widget
    v.view._remote_call("setSize", target="Widget", args=["300px", "300px"])
    v.view.center()
    return v

a = read('../../tests/aSi.cfg')
a.symbols[:] = 'Si'

interactive_view(a)

/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"

Now we setup the calculator with the Kumagai potential and its parametrization.

from matscipy.calculators.manybody import Manybody
from matscipy.calculators.manybody.explicit_forms import Kumagai
from matscipy.calculators.manybody.explicit_forms.kumagai import Kumagai_Comp_Mat_Sci_39_Si
from IPython.display import display, Markdown

a.calc = Manybody(**Kumagai(Kumagai_Comp_Mat_Sci_39_Si))

display(Markdown(f'Cohesive energy = {a.get_potential_energy() / len(a):.2f} eV/atom'))

Cohesive energy = -4.42 eV/atom

We can also compute the Born elastic constants, i.e. the affine elastic constants:

import numpy as np
from ase.units import GPa
from matscipy.elasticity import elastic_moduli, full_3x3x3x3_to_Voigt_6x6

# Born elastic constants (without nonaffine displacements)
C = a.calc.get_property('born_constants', a)

# Get isotropic elastic moduli
E, nu, G, B, K = elastic_moduli(full_3x3x3x3_to_Voigt_6x6(C))

display(Markdown(f"Young's modulus = {E.mean() / GPa:.1f} GPa"))
display(Markdown(f"Poisson number = {(nu + np.eye(3)).sum()/6:.3f}"))

Young’s modulus = 170.9 GPa

Poisson number = 0.212

Other many-body potentials

To runs the same code with Stillinger-Weber, replace the calculator by

from matscipy.calculators.manybody.explicit_forms import StillingerWeber
from matscipy.calculators.manybody.explicit_forms.stillinger_weber import Stillinger_Weber_PRB_31_5262_Si

a.calc = Manybody(**StillingerWeber(Stillinger_Weber_PRB_31_5262_Si))

display(Markdown(f'Cohesive energy = {a.get_potential_energy() / len(a):.2f} eV/atom'))

Cohesive energy = -4.11 eV/atom

Other examples of potentials that are implemented are Tersoff, Brenner (without the lookup tables), Erhart-Albe and others. Second derivatives are supported for all of these.