LennardJones

Basic implementation of a (multi-species) Lennard-Jones potential with finite cutoff radius that is imposed by "shifting and tilting" the potential at the cutoff. It can be constructed as follows.

emins = Dict( (z1, z1) => -1.0u"eV", 
              (z1, z2) => -0.5u"eV", 
              (z2, z2) => -0.25u"eV" )
rmins = Dict( (z1, z1) => 2.7u"Å", 
              (z1, z2) => 3.2u"Å", 
              (z2, z2) => 3.0u"Å" )
rcut = 6.0u"Å"              
lj = LennardJones(emins, rmins, rcut)

It is assumed that the potential is symmetric, i.e. emins[(z1, z2)] == emins[(z2, z1)] and so forth.

Morse

Basic implementation of a (multi-species) Morse potential with finite cutoff radius that is imposed by "shifting and tilting" the potential at the cutoff. It can be constructed as follows. The parameters are (energy scale, equilibrium bond length, stiffness parameter)

params = Dict( (z1, z1) => ( -1.0u"eV",  2.7u"Å", 4.1 ),     
               (z1, z2) => ( -0.5u"eV",  3.2u"Å", 3.5 ),
               (z2, z2) => ( -0.25u"eV", 3.0u"Å", 4.3 ) )
rcut = 6.0u"Å"              
V = Morse(params, rcut) 

It is assumed that the potential is symmetric, i.e. params[(z1, z2)] == params[(z2, z1)].

Stillinger-Weber potential with parameters for Si.

Functional form and default parameters match the original SW potential from [Stillinger/Weber, PRB 1985].

The StillingerWeber type can also by "abused" to generate arbitrary bond-angle potentials of the form Σᵢⱼ V₂(rᵢⱼ) + Σᵢⱼₖ V₃(rᵢⱼ) V₃(rᵢₖ) (cos Θᵢⱼₖ + 1/3)²

Constructor admits the following key-word parameters: ϵ=2.1675, σ = 2.0951, A=7.049556277, B=0.6022245584, p = 4, a = 1.8, λ=21.0, γ=1.20, atom_number=14

which enter the potential as follows:

V2(r) = 0.5 * ϵ * A * (B * (r/σ)^(-p) - 1.0) * exp(1.0 / (r/σ - a))
V3(r) = sqrt(ϵ * λ) * exp(γ / (r/σ - a))

The brittle keyword can be used to switch to the parameters λ to 42.0, which is appropriate to simulate brittle fracture. (need reference for this)

ZBL

Basic implementation of a ZBL potential (TODO insert reference). The original ZBL has not cutoff but the interface we use enforces a cutoff. This is the only parameter. It can be constructed as follows.

rcut = 6.0u"Å"              
zbl = ZBL(rcut) 

The current version of this potential assumes eV and Å as energy and force units. A PR to generalize this is welcome.

_z2i(zlist, z)

return an index in zlist or 0 if z was not found. If nothing should be returned instead of 0, then use findfirst instead.

sw_bondangle(S1, S2) -> (dot(S1, S2) + 1.0/3.0)^2

• not this assumes that S1, S2 are normalised
• see sw_bondangle_d for the derivative

b := sw_bondangle(S1, S2) then

sw_bondangle_d(S1, S2, r1, r2) -> b, db1, db2

where dbi is the derivative of b w.r.t. Ri where Si= Ri/ri.