Skip to content

Commit

Permalink
changes after Calvados PR
Browse files Browse the repository at this point in the history
  • Loading branch information
@jgreener64
jgreener64 committed Nov 7, 2024
1 parent 8de9e89 commit d7cb3aa
Show file tree
Hide file tree
Showing 3 changed files with 156 additions and 168 deletions.
61 changes: 28 additions & 33 deletions src/interactions/coulomb.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -342,10 +342,6 @@ end
end
end





@doc raw"""
Yukawa(; cutoff, use_neighbors, weight_special, coulomb_const, kappa)
Expand All @@ -355,30 +351,28 @@ The potential energy is defined as
```math
V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}} \exp(-\kappa r_{ij})
```
and the force on each atom by
```math
F(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}^2} \exp(-\kappa r_{ij})\left(\kappa r_{ij} + 1\right) \vec{r}_{ij}
```
"""
@kwdef struct Yukawa{C, W, T, D}
cutoff::C= NoCutoff()
use_neighbors::Bool= false
weight_special::W =1
@kwdef struct Yukawa{C, W, T, K}
cutoff::C = NoCutoff()
use_neighbors::Bool = false
weight_special::W = 1
coulomb_const::T = coulomb_const
kappa::D = 1.0*u"nm^-1"
kappa::K = 1.0*u"nm^-1"
end

use_neighbors(inter::Yukawa) = inter.use_neighbors

function Base.zero(yukawa::Yukawa{C, W, T ,D}) where {C, W, T, D}
return Yukawa{C, W, T, D}(
function Base.zero(yukawa::Yukawa{C, W, T, K}) where {C, W, T, K}
return Yukawa(
yukawa.cutoff,
yukawa.use_neighbors,
yukawa.weight_special,
yukawa.coulomb_const,
yukawa.kappa
zero(W),
zero(T),
zero(K),
)
end

Expand All @@ -387,18 +381,18 @@ function Base.:+(c1::Yukawa, c2::Yukawa)
c1.cutoff,
c1.use_neighbors,
c1.weight_special + c2.weight_special,
c1.coulomb_const,
c1.kappa
c1.coulomb_const + c2.coulomb_const,
c1.kappa + c2.kappa,
)
end

@inline function force(inter::Yukawa{C},
dr,
atom_i,
atom_j,
force_units=u"kJ * mol^-1 * nm^-1",
special=false,
args...) where C
@inline function force(inter::Yukawa,
dr,
atom_i,
atom_j,
force_units=u"kJ * mol^-1 * nm^-1",
special=false,
args...)
r2 = sum(abs2, dr)
cutoff = inter.cutoff
coulomb_const = inter.coulomb_const
Expand All @@ -416,15 +410,16 @@ end

function force_divr(::Yukawa, r2, invr2, (coulomb_const, qi, qj, kappa))
r = sqrt(r2)
return (coulomb_const * qi * qj) / r^3 *exp(-kappa*r) * (kappa*r+1)
return (coulomb_const * qi * qj) * exp(-kappa * r) * (kappa * r + 1) / r^3
end

@inline function potential_energy(inter::Yukawa{C},
dr,
atom_i,
atom_j,
energy_units=u"kJ * mol^-1",
special::Bool=false, args...) where C
@inline function potential_energy(inter::Yukawa,
dr,
atom_i,
atom_j,
energy_units=u"kJ * mol^-1",
special::Bool=false,
args...)
r2 = sum(abs2, dr)
cutoff = inter.cutoff
coulomb_const = inter.coulomb_const
Expand All @@ -440,5 +435,5 @@ end
end

function potential(::Yukawa, r2, invr2, (coulomb_const, qi, qj, kappa))
return (coulomb_const * qi * qj) * √invr2 * exp(-kappa*sqrt(r2))
return (coulomb_const * qi * qj) * √invr2 * exp(-kappa * sqrt(r2))
end
93 changes: 45 additions & 48 deletions src/interactions/lennard_jones.jl
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
export
LennardJones,
LennardJonesSoftCore,
AshbaughHatch,
AshbaughHatchAtom
AshbaughHatch

function lj_zero_shortcut(atom_i, atom_j)
return iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) ||
Expand Down Expand Up @@ -293,23 +292,23 @@ end


@doc raw"""
AshbaughHatch(; cutoff,use_neighbors,shortcut, ϵ_mixing,σ_mixing, λ_mixing,weight_special)
AshbaughHatch(; cutoff, use_neighbors, shortcut, ϵ_mixing, σ_mixing,
λ_mixing, weight_special)
The AshbaughHatch ($V_{\text{AH}}$) is a modified Lennard-Jones ($V_{\text{LJ}}$) 6-12 interaction between two atoms.
The Ashbaugh-Hatch potential ($V_{\text{AH}}$) is a modified Lennard-Jones ($V_{\text{LJ}}$)
6-12 interaction between two atoms.
The potential energy is defined by
The potential energy is defined as
```math
V_{\text{LJ}}(r_{ij}) = 4\varepsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right] \\
```
```math
V_{\text{AH}}(r_{ij}) =
\begin{cases}
V_{\text{LJ}}(r_{ij}) +\varepsilon_{ij}(1-λ_{ij}) &, r_{ij}\leq 2^{1/6}σ \\
λ_{ij}V_{\text{LJ}}(r_{ij}) &, 2^{1/6}σ \leq r_{ij}
\end{cases}
```
and the force on each atom by
```math
\vec{F}_{\text{AH}} =
Expand All @@ -328,30 +327,20 @@ where
If ``\lambda`` is one this gives the standard [`LennardJones`](@ref) potential.
"""
@kwdef struct AshbaughHatch{C, H, E, F,L,W}
@kwdef struct AshbaughHatch{C, H, S, E, L, W}
cutoff::C = NoCutoff()
use_neighbors::Bool = false
shortcut::H = lj_zero_shortcut
σ_mixing::S = lorentz_σ_mixing
ϵ_mixing::E = lorentz_ϵ_mixing
σ_mixing::F = lorentz_σ_mixing
λ_mixing::L = lorentz_λ_mixing
weight_special::W = 1
end


@kwdef struct AshbaughHatchAtom{C, M, S, E, L}
index::Int = 1
charge::C = 0.0
mass::M =1.0u"g/mol"
σ::S =0.0u"nm"
ϵ::E =0.0u"kJ * mol^-1"
λ::L=1.0
end

use_neighbors(inter::AshbaughHatch) = inter.use_neighbors

function Base.zero(lj::AshbaughHatch{C, H, E, F,L,W}) where {C, H, E, F,L,W}
return AshbaughHatch{C, H, E, F,L,W}(
function Base.zero(lj::AshbaughHatch{C, H, S, E, L, W}) where {C, H, S, E, L, W}
return AshbaughHatch(
lj.cutoff,
lj,use_neighbors,
lj.shortcut,
Expand All @@ -362,9 +351,8 @@ function Base.zero(lj::AshbaughHatch{C, H, E, F,L,W}) where {C, H, E, F,L,W}
)
end

function Base.:+(l1::AshbaughHatch{C, H, E, F,L,W},
l2::AshbaughHatch{C, H, E, F,L,W}) where {C, H, E, F,L,W}
return AshbaughHatch{C, H, E, F,L,W}(
function Base.:+(l1::AshbaughHatch, l2::AshbaughHatch)
return AshbaughHatch(
l1.cutoff,
l1.use_neighbors,
l1.shortcut,
Expand All @@ -375,18 +363,27 @@ function Base.:+(l1::AshbaughHatch{C, H, E, F,L,W},
)
end

@inline function force(inter::AshbaughHatch{S, C},
dr,
atom_i,
atom_j,
force_units=u"kJ * mol^-1 * nm^-1",
special::Bool=false, args...) where {S, C}
@kwdef struct AshbaughHatchAtom{T, M, C, S, E, L}
index::Int = 1
atom_type::T = 1
mass::M = 1.0u"g/mol"
charge::C = 0.0
σ::S = 0.0u"nm"
ϵ::E = 0.0u"kJ * mol^-1"
λ::L = 1.0
end

@inline function force(inter::AshbaughHatch,
dr,
atom_i,
atom_j,
force_units=u"kJ * mol^-1 * nm^-1",
special::Bool=false,
args...)
if inter.shortcut(atom_i, atom_j)
return ustrip.(zero(dr)) * force_units
end

# Lorentz-Berthelot mixing rules use the arithmetic average for σ
# Otherwise use the geometric average
ϵ = inter.ϵ_mixing(atom_i, atom_j)
σ = inter.σ_mixing(atom_i, atom_j)
λ = inter.λ_mixing(atom_i, atom_j)
Expand All @@ -396,7 +393,7 @@ end
σ2 = σ^2
params = (σ2, ϵ, λ)

f = force_divr_with_cutoff(inter, r2, params, cutoff,force_units)
f = force_divr_with_cutoff(inter, r2, params, cutoff, force_units)
if special
return f * dr * inter.weight_special
else
Expand All @@ -405,24 +402,24 @@ end
end

@inline function force_divr(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ))
if r2 < 2^(1/3)*σ2
if r2 < (2^(1/3) * σ2)
return force_divr(LennardJones(), r2, invr2, (σ2, ϵ))
else
return λ*force_divr(LennardJones(), r2, invr2, (σ2, ϵ))
return λ * force_divr(LennardJones(), r2, invr2, (σ2, ϵ))
end
end

@inline function potential_energy(inter::AshbaughHatch{S, C},
dr,
atom_i,
atom_j,
energy_units=u"kJ * mol^-1",
special::Bool=false,
args...) where {S, C}
@inline function potential_energy(inter::AshbaughHatch,
dr,
atom_i,
atom_j,
energy_units=u"kJ * mol^-1",
special::Bool=false,
args...)
if inter.shortcut(atom_i, atom_j)
return ustrip(zero(dr[1])) * energy_units
end
ϵ = inter.ϵ_mixing(atom_i, atom_j) ### probably not needed for Calvados2 use case
ϵ = inter.ϵ_mixing(atom_i, atom_j)
σ = inter.σ_mixing(atom_i, atom_j)
λ = inter.λ_mixing(atom_i, atom_j)

Expand All @@ -431,7 +428,7 @@ end
σ2 = σ^2
params = (σ2, ϵ, λ)

pe = potential_with_cutoff(inter, r2, params, cutoff,energy_units)
pe = potential_with_cutoff(inter, r2, params, cutoff, energy_units)
if special
return pe * inter.weight_special
else
Expand All @@ -440,9 +437,9 @@ end
end

@inline function potential(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ))
if r2 < 2^(1/3)*σ2
return potential(LennardJones(), r2, invr2, (σ2, ϵ)) + ϵ*(1-λ)
if r2 < (2^(1/3) * σ2)
return potential(LennardJones(), r2, invr2, (σ2, ϵ)) + ϵ * (1 - λ)
else
return λ* potential(LennardJones(), r2, invr2, (σ2, ϵ))
return λ * potential(LennardJones(), r2, invr2, (σ2, ϵ))
end
end
end
Loading

0 comments on commit d7cb3aa

Please sign in to comment.