From 3b085287bb232dbfe677f4c938f2927c9546c07d Mon Sep 17 00:00:00 2001 From: Stefan Bringuier Date: Tue, 1 Mar 2022 12:41:28 -0800 Subject: [PATCH 1/2] Initial implementation of EAM energy and force routines. Implemented the Sutton-Chen formulation. --- .../EmbeddedAtomPotentials/eam_potentials.jl | 67 ++++++++++++++++++ .../EmbeddedAtomPotentials/suttonchen.jl | 68 +++++++++++++++++++ src/PotentialTypes/types.jl | 9 ++- 3 files changed, 143 insertions(+), 1 deletion(-) create mode 100644 src/PotentialTypes/EmbeddedAtomPotentials/eam_potentials.jl create mode 100644 src/PotentialTypes/EmbeddedAtomPotentials/suttonchen.jl diff --git a/src/PotentialTypes/EmbeddedAtomPotentials/eam_potentials.jl b/src/PotentialTypes/EmbeddedAtomPotentials/eam_potentials.jl new file mode 100644 index 0000000..eaa17f2 --- /dev/null +++ b/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 diff --git a/src/PotentialTypes/EmbeddedAtomPotentials/suttonchen.jl b/src/PotentialTypes/EmbeddedAtomPotentials/suttonchen.jl new file mode 100644 index 0000000..9ef2d3a --- /dev/null +++ b/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 diff --git a/src/PotentialTypes/types.jl b/src/PotentialTypes/types.jl index 9d3448b..154af02 100644 --- a/src/PotentialTypes/types.jl +++ b/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 - From 944171247b67c3143285224f679d25cb591a75c8 Mon Sep 17 00:00:00 2001 From: Stefan Bringuier Date: Tue, 1 Mar 2022 12:42:20 -0800 Subject: [PATCH 2/2] Test for EAM Sutton-Chen. Not passing, i.e., not matching expected values. --- test/runtests.jl | 1 + test/suttonchen_test.jl | 69 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 70 insertions(+) create mode 100644 test/suttonchen_test.jl diff --git a/test/runtests.jl b/test/runtests.jl index b8bf3fd..b03c03a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,4 +2,5 @@ using Test @time @testset "Potentials.jl" begin include("lj_test.jl") + include("suttonchen.jl") end diff --git a/test/suttonchen_test.jl b/test/suttonchen_test.jl new file mode 100644 index 0000000..43e6e4e --- /dev/null +++ b/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)