diff --git a/docs/src/documentation.md b/docs/src/documentation.md index f125a92b..8353186e 100644 --- a/docs/src/documentation.md +++ b/docs/src/documentation.md @@ -543,12 +543,14 @@ In Molly there are three types of interactions: The available pairwise interactions are: - [`LennardJones`](@ref) - [`LennardJonesSoftCore`](@ref) +- [`AshbaughHatch`](@ref) - [`SoftSphere`](@ref) - [`Mie`](@ref) - [`Buckingham`](@ref) - [`Coulomb`](@ref) - [`CoulombSoftCore`](@ref) - [`CoulombReactionField`](@ref) +- [`Yukawa`](@ref) - [`Gravity`](@ref) The available specific interactions are: diff --git a/src/interactions/coulomb.jl b/src/interactions/coulomb.jl index 46815d98..6c040864 100644 --- a/src/interactions/coulomb.jl +++ b/src/interactions/coulomb.jl @@ -1,7 +1,8 @@ export Coulomb, CoulombSoftCore, - CoulombReactionField + CoulombReactionField, + Yukawa const coulomb_const = 138.93545764u"kJ * mol^-1 * nm" # 1 / 4πϵ0 @@ -340,3 +341,104 @@ end return pe end end + + + + + +@doc raw""" + Yukawa(; cutoff, use_neighbors, weight_special, coulomb_const, kappa) + +The Yukawa electrostatic interaction between two atoms. + +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 + coulomb_const::T = coulomb_const + kappa::D = 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}( + yukawa.cutoff, + yukawa.use_neighbors, + yukawa.weight_special, + yukawa.coulomb_const, + yukawa.kappa + ) +end + +function Base.:+(c1::Yukawa, c2::Yukawa) + return Yukawa( + c1.cutoff, + c1.use_neighbors, + c1.weight_special + c2.weight_special, + c1.coulomb_const, + c1.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 + r2 = sum(abs2, dr) + cutoff = inter.cutoff + coulomb_const = inter.coulomb_const + qi, qj = atom_i.charge, atom_j.charge + kappa = inter.kappa + params = (coulomb_const, qi, qj, kappa) + + f = force_divr_with_cutoff(inter, r2, params, cutoff, force_units) + if special + return f * dr * inter.weight_special + else + return f * dr + end +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) +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 + r2 = sum(abs2, dr) + cutoff = inter.cutoff + coulomb_const = inter.coulomb_const + qi, qj = atom_i.charge, atom_j.charge + params = (coulomb_const, qi, qj, inter.kappa) + + pe = potential_with_cutoff(inter, r2, params, cutoff, energy_units) + if special + return pe * inter.weight_special + else + return pe + end +end + +function potential(::Yukawa, r2, invr2, (coulomb_const, qi, qj, kappa)) + return (coulomb_const * qi * qj) * √invr2 * exp(-kappa*sqrt(r2)) +end diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 7cc84375..0aa8ef35 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -1,6 +1,8 @@ export LennardJones, - LennardJonesSoftCore + LennardJonesSoftCore, + AshbaughHatch, + AshbaughHatchAtom function lj_zero_shortcut(atom_i, atom_j) return iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) || @@ -13,6 +15,14 @@ function lorentz_σ_mixing(atom_i, atom_j) return (atom_i.σ + atom_j.σ) / 2 end +function lorentz_ϵ_mixing(atom_i, atom_j) + return (atom_i.ϵ + atom_j.ϵ) / 2 +end + +function lorentz_λ_mixing(atom_i, atom_j) + return (atom_i.λ + atom_j.λ) / 2 +end + function geometric_σ_mixing(atom_i, atom_j) return sqrt(atom_i.σ * atom_j.σ) end @@ -280,3 +290,159 @@ function potential(::LennardJonesSoftCore, r2, invr2, (σ2, ϵ, σ6_fac)) six_term = σ2^3 * inv(r2^3 + σ2^3 * σ6_fac) return 4ϵ * (six_term ^ 2 - six_term) end + + +@doc raw""" + 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 potential energy is defined by +```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}} = + \begin{cases} + F_{\text{LJ}}(r_{ij}) &, r_{ij} \leq 2^{1/6}σ \\ + λ_{ij}F_{\text{LJ}}(r_{ij}) &, 2^{1/6}σ \leq r_{ij} + \end{cases} +``` +where +```math +\begin{aligned} +\vec{F}_{\text{LJ}}\ +&= \frac{24\varepsilon_{ij}}{r_{ij}^2} \left[2\left(\frac{\sigma_{ij}^{6}}{r_{ij}^{6}}\right)^2 -\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right] \vec{r_{ij}} +\end{aligned} +``` + +If ``\lambda`` is one this gives the standard [`LennardJones`](@ref) potential. +""" +@kwdef struct AshbaughHatch{C, H, E, F,L,W} + cutoff::C = NoCutoff() + use_neighbors::Bool = false + shortcut::H = lj_zero_shortcut + ϵ_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}( + lj.cutoff, + lj,use_neighbors, + lj.shortcut, + lj.σ_mixing, + lj.ϵ_mixing, + lj.λ_mixing, + zero(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}( + l1.cutoff, + l1.use_neighbors, + l1.shortcut, + l1.σ_mixing, + l1.ϵ_mixing, + l1.λ_mixing, + l1.weight_special + l2.weight_special, + ) +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} + 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) + + cutoff = inter.cutoff + r2 = sum(abs2, dr) + σ2 = σ^2 + params = (σ2, ϵ, λ) + + f = force_divr_with_cutoff(inter, r2, params, cutoff,force_units) + if special + return f * dr * inter.weight_special + else + return f * dr + end +end + +@inline function force_divr(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ)) + if r2 < 2^(1/3)*σ2 + return force_divr(LennardJones(), r2, invr2, (σ2, ϵ)) + else + 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} + 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) + + cutoff = inter.cutoff + r2 = sum(abs2, dr) + σ2 = σ^2 + params = (σ2, ϵ, λ) + + pe = potential_with_cutoff(inter, r2, params, cutoff,energy_units) + if special + return pe * inter.weight_special + else + return pe + end +end + +@inline function potential(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ)) + if r2 < 2^(1/3)*σ2 + return potential(LennardJones(), r2, invr2, (σ2, ϵ)) + ϵ*(1-λ) + else + return λ* potential(LennardJones(), r2, invr2, (σ2, ϵ)) + end +end \ No newline at end of file diff --git a/test/interactions.jl b/test/interactions.jl index f6706604..2e857518 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -33,6 +33,29 @@ ) end + function do_shortcut(atom_i, atom_j) return true end + + for inter in (LennardJones(;shortcut=do_shortcut), + Mie(m=6, n=12;shortcut=do_shortcut), + LennardJonesSoftCore(α=1, λ=0, p=2;shortcut=do_shortcut), + SoftSphere(;shortcut=do_shortcut), + Buckingham(;shortcut=do_shortcut), + AshbaughHatch(;shortcut=do_shortcut), + ) + + @test isapprox( + force(inter, dr12, a1, a1), + SVector(0.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-9u"kJ * mol^-1 * nm^-1", + ) + + @test isapprox( + potential_energy(inter, dr12, a1, a1), + 0.0u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + end + inter = LennardJonesSoftCore(α=1, λ=0.5, p=2) @test isapprox( force(inter, dr12, a1, a1), @@ -83,7 +106,7 @@ B::TB C::TC end - + buck_c1 = SVector(10.0, 10.0, 10.0)u"Å" buck_c2 = SVector(13.0, 10.0, 10.0)u"Å" buck_c3 = SVector(14.0, 10.0, 10.0)u"Å" @@ -159,6 +182,37 @@ atol=1e-5u"kJ * mol^-1", ) + inter = Yukawa(; weight_special=0.5) + @test isapprox( + force(inter, dr12, a1, a1), + SVector(1486.7077156786308, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-5u"kJ * mol^-1 * nm^-1", + ) + @test isapprox( + force(inter, dr13, a1, a1), + SVector(814.8981977722481, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-5u"kJ * mol^-1 * nm^-1", + ) + @test isapprox( + force(inter, dr13, a1, a1,u"kJ * mol^-1 * nm^-1",true), + SVector(0.5*814.8981977722481, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-5u"kJ * mol^-1 * nm^-1", + ) + @test isapprox( + potential_energy(inter, dr12, a1, a1), + 343.08639592583785u"kJ * mol^-1"; + atol=1e-5u"kJ * mol^-1", + ) + @test isapprox( + potential_energy(inter, dr13, a1, a1), + 232.8280565063u"kJ * mol^-1"; + atol=1e-5u"kJ * mol^-1", + ) + @test isapprox( + potential_energy(inter, dr13, a1, a1,u"kJ * mol^-1 * nm^-1",true), + 0.5*232.8280565063u"kJ * mol^-1"; + atol=1e-5u"kJ * mol^-1", + ) c1_grav = SVector(1.0, 1.0, 1.0)u"m" c2_grav = SVector(6.0, 1.0, 1.0)u"m" a1_grav, a2_grav = Atom(mass=1e6u"kg"), Atom(mass=1e5u"kg") @@ -385,4 +439,67 @@ -108.166724117u"kJ * mol^-1"; atol=1e-7u"kJ * mol^-1", ) + AH_a1 = AshbaughHatchAtom(charge=1.0, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1", λ=1.0) + inter = AshbaughHatch(weight_special=0.5) + + @test isapprox( + force(inter, dr12, AH_a1, AH_a1, boundary), + SVector(16.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-9u"kJ * mol^-1 * nm^-1", + ) + @test isapprox( + force(inter, dr13, AH_a1, AH_a1, boundary), + SVector(-1.375509739, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-9u"kJ * mol^-1 * nm^-1", + ) + @test isapprox( + potential_energy(inter, dr12, AH_a1, AH_a1, boundary), + 0.0u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + @test isapprox( + potential_energy(inter, dr13, AH_a1, AH_a1, boundary), + -0.1170417309u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + + + AH_a1 = AshbaughHatchAtom(charge=1.0, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1", λ=0.5) + @test isapprox( + potential_energy(inter, dr13, AH_a1, AH_a1), + -0.058520865u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + + @test isapprox( + force(inter, dr13, AH_a1, AH_a1), + SVector(-0.68775486946, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-9u"kJ * mol^-1 * nm^-1", + ) + + c4 = SVector(1.28, 1.0, 1.0)u"nm" + dr14 = vector(c1, c4, boundary) + @test isapprox( + potential_energy(inter, dr14, AH_a1, AH_a1), + 0.7205987916u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + + @test isapprox( + force(inter, dr14, AH_a1, AH_a1), + SVector(52.5306754422, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-9u"kJ * mol^-1 * nm^-1", + ) + + @test isapprox( + potential_energy(inter, dr14, AH_a1, AH_a1,u"kJ * mol^-1 * nm^-1", true), + 0.5*0.7205987916u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + + @test isapprox( + force(inter, dr14, AH_a1, AH_a1,u"kJ * mol^-1 * nm^-1",true), + SVector(0.5*52.5306754422, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-9u"kJ * mol^-1 * nm^-1", + ) end