PairPotential:abstractsupertype for pair potentials

SitePotential:abstractsupertype for generic site potentials. Concrete subtypes should overload the cutoff_radius, eval_site and eval_grad_site methods.

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)

write docs...

If V <: SitePotential then it should implement the method

dv = eval_grad_site(V, Rs, Zs, z0)

where Rs::AbstractVector{<: SVector{3}} and Zs::AbstractVector of atom ids (e.g., atomic numbers), while z0 is a single atom id.

The output dv should be an AbstractVector containing SVector{3,T} blocks.

If V <: SitePotential then it should implement the method

val = eval_site(V, Rs, Zs, z0)

where Rs::AbstractVector{<: SVector{3}} and Zs::AbstractVector of atom ids (e.g. atomic numbers), while z0 is a single atom id.

The output val should be a single number, namely the site energy.

get_neighbours(nlist::PairList, at, i::Integer) -> Js, Rs, Zs, z0

If V <: SitePotential then it can implement the method

Pblock = precon(V, Rs, Zs, z0)

where Rs::AbstractVector{<: SVector{3}} and Zs::AbstractVector of atom ids (e.g., atomic numbers), while z0 is a single atom id. The output Pblock should be an AbstractMatrix containing SMatrix{3,3,T} blocks.

Unlike eval_site and eval_grad_site, this method is optional. It can be used to speedup geometry optimization, sampling and related tasks.

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.