Embedded Atom Method Potential

src/PotentialTypes/EmbeddedAtomPotentials/eam_potentials.jl

@@ -0,0 +1,67 @@
+# Contributions by: Stefan Bringuier
+################################################################################
+# Types of Potentials
+################################################################################
+
+include("suttonchen.jl")
+
+export SuttonChen
+export get_parameters, get_hyperparameters, serialize_parameters, serialize_hyperparameters
+
+################################################################################
+# InteratomicPotentials API implmentations for EAM potentials
+################################################################################
+
+function energy_and_force(A::AbstractSystem, p::EmbeddedAtomPotential)
+    nnlist = neighborlist(A, p.rcutoff)
+
+    e = 0.0
+    f = fill(SVector{3}(zeros(3)), length(A))
+    ρ = zeros(length(A))
+    dρ = fill(SVector{3}(zeros(3)),length(A))
+    for ii in 1:length(A)
+        for (jj, r, R) in zip(nnlist.j[ii], nnlist.r[ii], nnlist.R[ii])
+            species = unique([atomic_symbol(A, ii), atomic_symbol(A, jj)])
+            if (intersect(species, p.species) == species)
+
+                ρ[ii] += rho(R,p)
+                ρ[jj] += rho(R,p)
+                dρᵢⱼ = drho_dR(R,r,p)
+                dρ[ii] = dρ[ii] + dρᵢⱼ
+                dρ[jj] = dρ[jj] - dρᵢⱼ
+
+                e += potential_energy_repulsive(R, p)
+                fo = force_repulsive(R, r, p)
+                f[ii] = f[ii] + fo
+                f[jj] = f[jj] - fo
+            end
+        end
+        e -= potential_energy_embedding(ρ[ii],p)
+        f[ii] = f[ii] - force_embedding(ρ[ii],dρ[ii],p)
+    end
+    (; e, f)
+end
+
+
+function force(R::AbstractFloat,r::SVector{3,<:AbstractFloat}, p::EmbeddedAtomPotential)
+    ρ,dρ = rho(R,p), drho_dR(R,r,p)
+    fo = force_repulsive(R,r,p)
+    fo -= force_embedding(ρ,dρ,p)
+    return fo
+end
+
+force(r::SVector{3,<:AbstractFloat}, p::EmbeddedAtomPotential) = force(norm(r),r,p)
+
+function virial_stress(A::AbstractSystem, p::EmbeddedAtomPotential)
+    nnlist = neighborlist(A, p.rcutoff)
+
+    v = SVector{6}(zeros(6))
+    for ii in 1:length(A)
+        for (r, R) in zip(nnlist.r[ii], nnlist.R[ii])
+            fo = force(R, r, p)
+            vi = r * fo'
+            v = v + [vi[1, 1], vi[2, 2], vi[3, 3], vi[3, 2], vi[3, 1], vi[2, 1]]
+        end
+    end
+    v
+end

src/PotentialTypes/EmbeddedAtomPotentials/suttonchen.jl

@@ -0,0 +1,68 @@
+# Contributions by: Stefan Bringuier
+# Reference: A. P. Sutton and J. Chen, Phil. Mag. Lett. 61, 139 (1990)
+
+##############################   Sutton-Chen EAM  ###################################
+struct SuttonChen <:EmbeddedAtomPotential
+    ϵ
+    a
+    c
+    n
+    m
+    rcutoff
+    species :: AbstractVector
+end
+
+get_parameters(sc::SuttonChen) = Parameter{ (:ϵ,:a,:c,:n,:m) }((sc.ϵ,sc.a,sc.c,sc.n,sc.m))
+set_parameters(p::Parameter{(:ϵ,:a,:c,:n,:m)}, sc::SuttonChen) = SuttonChen(p.ϵ,p.a,p.c,p.n,p.m,sc.rcutoff,sc.species)
+
+deserialize_parameters(p::Parameter{(:ϵ, :a, :c, :n, :m)}, sc::SuttonChen) = [p.ϵ, p.a, p.c, p.n, p.m]
+serialize_parameters(p::Vector, sc::SuttonChen) = Parameter{(:ϵ, :a, :c, :n, :m)}( (p[1], p[2], p[3], p[4], p[5]) )
+
+get_hyperparameters(sc::SuttonChen) = Parameter{(:rcutoff,)}( (sc.rcutoff,) )
+set_hyperparameters(p::Parameter{(:rcutoff,)}, sc::SuttonChen) = SuttonChen(sc.ϵ, sc.a, sc.c, sc.n, sc.m, p.rcutoff, sc.species)
+
+deserialize_hyperparameters(p::Parameter{(:rcutoff,)}, sc::SuttonChen) = [p.rcutoff]
+serialize_hyperparameters(p::Vector, sc::SuttonChen) = Parameter{(:rcutoff, )}( (p[1],) )
+
+# ##############################   Cutoff  ###################################
+function fc(R::AbstractFloat, sc::SuttonChen)
+    return R-sc.rcutoff < 0.0 ? exp(1/(R-sc.rcutoff)) : 0.0
+end
+
+function dfc_dR(R::AbstractFloat, sc::SuttonChen)
+    return R-sc.rcutoff < 0.0 ? (-1/(R-sc.rcutoff))*exp(1/(R-sc.rcutoff)) : 0.0
+end
+
+
+# ##############################   Density  ###################################
+function rho(R::AbstractFloat, sc::SuttonChen)
+    ρ = (sc.a / R)^(sc.m) * fc(R,sc)
+    return ρ
+end
+
+function drho_dR(R::AbstractFloat,r::SVector{3,<:AbstractFloat},sc::SuttonChen)
+    fo = (-sc.m/R) * (sc.a/R)^(sc.m) * fc(R,sc) + (sc.a/R)^(sc.m) * dfc_dR(R,sc)
+    return SVector(fo .* r ./ R)
+end
+
+# ##############################   Energy  ###################################
+function potential_energy_repulsive(R::AbstractFloat, sc::SuttonChen)
+    return sc.ϵ * (sc.a/R)^(sc.n)  * fc(R,sc)
+end
+
+function potential_energy_embedding(ρ::AbstractFloat,sc::SuttonChen)
+    return sc.ϵ * sc.c * sqrt(ρ)
+end
+# ##############################   Force   ###################################
+
+function force_repulsive(R::AbstractFloat, r::SVector{3,<:AbstractFloat}, sc::SuttonChen)
+    fo = sc.ϵ * ((-sc.n/R) * (sc.a/R)^(sc.n) * fc(R,sc) + (sc.a/R)^(sc.n) * dfc_dR(R,sc))
+    return SVector(fo .* r ./ R)
+end
+
+function force_embedding(ρ::AbstractFloat, dρ::SVector{3,<:AbstractFloat}, sc::SuttonChen)
+    return sc.ϵ * (sc.c/(2*sqrt(ρ)) * dρ)
+end
+
+# ##############################   Gradients  ###################################
+# TO DO

src/PotentialTypes/types.jl

@@ -6,6 +6,7 @@
 
 abstract type ArbitraryPotential end
 abstract type EmpiricalPotential <: ArbitraryPotential end
+abstract type EmbeddedAtomPotential <: ArbitraryPotential end
 abstract type MixedPotential <: ArbitraryPotential end
 
 ################################################################################
@@ -20,9 +21,15 @@ export ArbitraryPotential
 include("EmpiricalPotentials/empirical_potentials.jl")
 export EmpiricalPotential
 
+
+################################################################################
+# Manybody Potentials
+################################################################################
+include("EmbeddedAtomPotentials/eam_potentials.jl")
+export EmbeddedAtomPotential
+
 ################################################################################
 # Mixed Potentials
 ################################################################################
 include("MixedPotentials/mixed_potentials.jl")
 export MixedPotential
-

test/runtests.jl

@@ -2,4 +2,5 @@ using Test
 
 @time @testset "Potentials.jl" begin
     include("lj_test.jl")
+    include("suttonchen.jl")
 end

test/suttonchen_test.jl

@@ -0,0 +1,69 @@
+# Contributions by: Stefan Bringuier
+
+using AtomsBase
+using InteratomicPotentials
+using StaticArrays
+using Unitful
+using UnitfulAtomic
+
+# Reference: A. P. Sutton and J. Chen, Phil. Mag. Lett. 61, 139 (1990)
+ϵ = 1.2832e-2
+a,c = 3.61,39.432
+nn,m = 9,6
+rcutoff = 10.0
+rc = 5.5
+sc = SuttonChen(ϵ,a,c,nn,m,rc, [:Cu,:Cu])
+
+pair_dist,energy = [], []
+lcell = 15.000
+for i in 0.0:0.05:5.0
+    atom1 = Atom(element, (@SVector [0.0, 0.0, 0.0]) * u"Å")
+    atom2 = Atom(element, (@SVector [1.0+i, 0.0, 0.0]) * u"Å")
+    atoms = [atom1, atom2]
+    box = [[lcell, 0.0, 0.0], [0.0, lcell, 0.0], [0.0, 0.0, lcell]]
+    bcs = [Periodic(), Periodic(), Periodic()]
+    system = FlexibleSystem(atoms, box * u"Å", bcs)
+    n = length(system)
+    nnlist = neighborlist(system, rcutoff)
+    push!(pair_dist,nnlist.R[1][1])
+    push!(energy,potential_energy(system,sc))
+end
+
+min_energy, iloc = findmin(energy/length(system))
+min_dist = pair_dist[iloc]
+@test -1.1 < min_energy < -1.3
+@test 2.05e0 < min_dist < 2.2e0
+
+e_binding = -3.54 # eV/atom
+eqpair_dist = 2.578 # Å
+lat_param = 3.615 # Å
+
+function fcc_cu(;latp=3.615)
+    atom1 = Atom(element, (@SVector [0.0, 0.0, 0.0]) * latp*u"Å")
+    atom2 = Atom(element, (@SVector [0.5, 0.5, 0.0]) * latp*u"Å")
+    atom3 = Atom(element, (@SVector [0.5, 0.0, 0.5]) * latp*u"Å")
+    atom4 = Atom(element, (@SVector [0.0, 0.5, 0.5]) * latp*u"Å")
+    atoms = [atom1,atom2,atom3,atom4]
+    box = [[latp, 0.0, 0.0], [0.0, latp, 0.0], [0.0, 0.0, latp]]
+    bcs = [Periodic(), Periodic(), Periodic()]
+    system = FlexibleSystem(atoms, box * u"Å", bcs)
+    return system
+end
+
+pair_dist,fcc_energy  = [], []
+for latp in 2.25:0.005:8.00
+    system = fcc_cu(latp=latp)
+    n = length(system)
+    nnlist = neighborlist(system, rcutoff)
+    push!(pair_dist,nnlist.R[1][1])
+    push!(fcc_energy,potential_energy(system,sc))
+end
+
+min_fcc_energy, iloc = findmin(fcc_energy/length(system))
+min_dist = pair_dist[iloc]
+@test -3.545 < min_fcc_energy < -3.550
+@test 2.54e0 < min_dist < 2.58e0
+
+
+force_equil_fcc  = force(fcc_cu(),sc)
+@test sum(force_equil_fcc) ≈ SVector(0.0,0.0,0.0)