Merge pull request #22 from JuliaMolSim/acutils

Generic codes to AtomsCalculatorsUtils and Update to AC@0.2

.github/workflows/CI.yml

@@ -19,8 +19,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.8'
-          - '1'
+          - '1.9'
+          - '1.10'
           - 'nightly'
         os:
           - ubuntu-latest

.gitignore

@@ -2,3 +2,4 @@
 /docs/Manifest.toml
 /docs/build/
 tune.json
+.vscode

Project.toml

@@ -1,33 +1,26 @@
 name = "EmpiricalPotentials"
 uuid = "38527215-9240-4c91-a638-d4250620c9e2"
-authors = ["Christoph Ortner <christohortner@gmail.com> and contributors"]
-version = "0.1.3"
+version = "0.2.0-dev"
 
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
 AtomsCalculators = "a3e0e189-c65a-42c1-833c-339540406eb1"
-ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
-DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
-Folds = "41a02a25-b8f0-4f67-bc48-60067656b558"
+AtomsCalculatorsUtilities = "9855a07e-8816-4d1b-ac92-859c17475477"
+Bumper = "8ce10254-0962-460f-a3d8-1f77fea1446e"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NeighbourLists = "2fcf5ba9-9ed4-57cf-b73f-ff513e316b9c"
-ObjectPools = "658cac36-ff0f-48ad-967c-110375d98c9d"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 AtomsBase = "0.3"
-AtomsCalculators = "0.1.2, 0.2"
-ChunkSplitters = "2"
-DiffResults = "1"
-Folds = "0.2.8"
+AtomsCalculators = "0.2"
+AtomsCalculatorsUtilities = "0.1.0"
+Bumper = "0.6, 0.7"
 ForwardDiff = "0.10"
-JSON = "0.21"
 LinearAlgebra = "1"
 NeighbourLists = "0.5.6"
-ObjectPools = "0.3.1"
 StaticArrays = "1.7.0"
 Unitful = "1"
 julia = "1.8"

README.md

@@ -3,3 +3,15 @@
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaMolSim.github.io/EmpiricalPotentials.jl/stable/)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaMolSim.github.io/EmpiricalPotentials.jl/dev/)
 [![Build Status](https://github.com/JuliaMolSim/EmpiricalPotentials.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaMolSim/EmpiricalPotentials.jl/actions/workflows/CI.yml?query=branch%3Amain)
+
+Implementation of `AtomsBase` and `AtomsCalculators` compatible 
+empirical interatomic potentials. At the moment, the following 
+potentials are provided: 
+- LennardJones (multi-species) 
+- Morse (multi-species)
+- ZBL
+- StillingerWeber (Si)
+
+EAM is planned, but there is no ETA. At the moment potentials have fixed parameters. Extension for parameterized potentials (low-level AtomsCalculators interface) are planned, but also no ETA. 
+
+Issues to request adding other potentials or functionality or PRs implementing them are welcome.

src/EmpiricalPotentials.jl

@@ -1,10 +1,30 @@
 module EmpiricalPotentials
 
-# Write your package code here.
+using Unitful, ForwardDiff, StaticArrays, Bumper 
+
+using ForwardDiff: Dual
+
+import AtomsCalculators: energy_unit, length_unit, force_unit 
+
+import AtomsCalculatorsUtilities
+import AtomsCalculatorsUtilities.SitePotentials
+import AtomsCalculatorsUtilities.SitePotentials: eval_site, eval_grad_site, 
+                                 hessian_site, block_hessian_site,  
+                                 cutoff_radius, SitePotential, 
+                                 ad_block_hessian_site, ad_hessian_site, 
+                                 cutoff_radius
+
+import AtomsCalculatorsUtilities.PairPotentials: PairPotential, 
+                                                 eval_pair 
+
+include("utils.jl")                                                 
+
+include("lennardjones.jl")
+include("morse.jl")
+include("zbl.jl")
+
+include("stillingerweber.jl")
 
-include("sitepotentials.jl")
-include("pairpotentials.jl")
-include("stillingerweber/stillingerweber.jl")
 
 
 end

src/lennardjones.jl

@@ -0,0 +1,93 @@
+
+
+export LennardJones
+
+
+# ----------- Lennard-Jones potential
+# to make this parametric, all we have to do is move the emins and rmins 
+# into a `ps` parameter NamedTuple. I suggest to do this in the next iteration 
+# and first stabilize the non-parametric potentials. 
+
+""" 
+    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.
+```julia 
+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. 
+"""
+mutable struct LennardJones{NZ, TZ, T, UL, UE} <: PairPotential 
+    zlist::NTuple{NZ, TZ}
+    emins::SMatrix{NZ, NZ, T} 
+    rmins::SMatrix{NZ, NZ, T}
+    rcut::T
+end
+
+_fltype(::LennardJones{NZ, TZ, T}) where {NZ, TZ, T} = T 
+
+length_unit(::LennardJones{NZ, TZ, T, UL, UE}) where {NZ, TZ, T, UL, UE} = UL 
+
+energy_unit(::LennardJones{NZ, TZ, T, UL, UE}) where {NZ, TZ, T, UL, UE} = UE 
+
+cutoff_radius(V::LennardJones) = V.rcut * length_unit(V)
+
+
+function eval_pair(V::LennardJones, r, z1, z2)
+    iz1 = _z2i(V.zlist, z1)
+    iz2 = _z2i(V.zlist, z2)
+    if iz1 == 0 || iz1 == 0 || r > V.rcut
+        return zero(_fltype(V))
+    end
+
+    _lj(s) = s^12 - 2 * s^6
+    _dlj(s) = 12 * s^11 - 12 * s^5 
+    _lj_tilt(s, scut) = _lj(s) - _lj(scut) - (s - scut) * _dlj(scut)
+
+    rmin = V.rmins[iz1, iz2]
+    emin = V.emins[iz1, iz2]
+    s = rmin / r 
+    scut = rmin / V.rcut 
+
+    return emin * _lj_tilt(s, scut)
+end
+
+
+function LennardJones(emins::Dict, rmins::Dict, rcut::Unitful.Length) 
+
+    zlist = tuple( unique(reduce(vcat, collect.(keys(emins))))... )
+    NZ = length(zlist)
+    TZ = typeof(first(zlist))
+    UL = unit(rcut) 
+    UE = unit(first(values(emins)))
+    T = typeof(ustrip(rcut))
+    @assert all(typeof(v) == TZ for v in zlist)
+    @assert all(unit(v) == UL for v in values(rmins))
+    @assert all(unit(v) == UE for v in values(emins)) 
+    @assert all(typeof(ustrip(v)) == T for v in values(rmins))
+    @assert all(typeof(ustrip(v)) == T for v in values(emins))
+
+    _emin(z1, z2) = haskey(emins, (z1, z2)) ? emins[(z1, z2)] : emins[(z2, z1)]
+    _rmin(z1, z2) = haskey(rmins, (z1, z2)) ? rmins[(z1, z2)] : rmins[(z2, z1)]
+
+    emins = SMatrix{NZ, NZ}([ ustrip(_emin(z1, z2)) 
+                              for z1 in zlist, z2 in zlist ])
+    rmins = SMatrix{NZ, NZ}([ ustrip(_rmin(z1, z2)) 
+                              for z1 in zlist, z2 in zlist ])
+
+
+    return LennardJones{NZ, TZ, T, UL, UE}(zlist, emins, rmins, ustrip(rcut))
+end
+
+

src/morse.jl

@@ -0,0 +1,87 @@
+
+
+export Morse
+
+
+
+# ----------- Lennard-Jones potential
+# to make this parametric, all we have to do is move the emins and rmins 
+# into a `ps` parameter NamedTuple. I suggest to do this in the next iteration 
+# and first stabilize the non-parametric potentials. 
+
+""" 
+    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)
+```julia 
+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)]`.
+"""
+mutable struct Morse{NZ, TZ, T, UL, UE} <: PairPotential 
+    zlist::NTuple{NZ, TZ}
+    params::SMatrix{NZ, NZ, Tuple{T, T, T}}
+    rcut::T
+end
+
+_fltype(::Morse{NZ, TZ, T}) where {NZ, TZ, T} = T 
+
+length_unit(::Morse{NZ, TZ, T, UL, UE}) where {NZ, TZ, T, UL, UE} = UL 
+
+energy_unit(::Morse{NZ, TZ, T, UL, UE}) where {NZ, TZ, T, UL, UE} = UE 
+
+cutoff_radius(V::Morse) = V.rcut * length_unit(V)
+
+
+function eval_pair(V::Morse, r, z1, z2)
+    iz1 = _z2i(V.zlist, z1)
+    iz2 = _z2i(V.zlist, z2)
+    if iz1 == 0 || iz1 == 0 || r > V.rcut
+        return zero(_fltype(V))
+    end
+
+    _m(s) = exp(-2 * s) - 2 * exp(-s)
+    _dm(s) = - 2 * exp(-2 * s) + 2 * exp(-s)
+    _m_tilt(s, scut) = _m(s) - _m(scut) - (s - scut) * _dm(scut)
+
+    emin, r0, α = V.params[iz1, iz2]
+    s = α * (r / r0 - 1)
+    scut = α * (V.rcut / r0 - 1)
+
+    return emin * _m_tilt(s, scut)
+end
+
+
+function Morse(params::Dict, rcut::Unitful.Length)
+
+    zlist = tuple( unique(reduce(vcat, collect.(keys(params))))... )
+    NZ = length(zlist)
+    TZ = typeof(first(zlist))
+    UL = unit(rcut) 
+    UE = unit(first(values(params))[1])
+    T = typeof(ustrip(rcut))
+    @assert all(typeof(v) == TZ for v in zlist)
+    @assert all(unit(v[1]) == UE for v in values(params)) 
+    @assert all(unit(v[2]) == UL for v in values(params))
+    @assert all(typeof(ustrip(v[1])) == T for v in values(params))
+    @assert all(typeof(ustrip(v[2])) == T for v in values(params))
+    @assert all(typeof(v[3]) == T for v in values(params))
+
+    _params(z1, z2) = haskey(params, (z1, z2)) ? params[(z1, z2)] : params[(z2, z1)]
+
+    P = SMatrix{NZ, NZ, Tuple{T, T, T}}(
+                    [ustrip.(_params(z1, z2)) for z1 in zlist, z2 in zlist])
+
+    return Morse{NZ, TZ, T, UL, UE}(zlist, P, ustrip(rcut))
+end
+
+