From 36ed01e90d6d6c0c1be0cff31ed0b06175ea4413 Mon Sep 17 00:00:00 2001 From: ywitzky <134135819+ywitzky@users.noreply.github.com> Date: Thu, 26 Sep 2024 11:15:52 +0200 Subject: [PATCH 1/8] added docs and first ashbaugh hatch implementation --- src/interactions/lennard_jones.jl | 191 +++++++++++++++++++++++++++++- 1 file changed, 190 insertions(+), 1 deletion(-) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 3dc1e83a..279b001a 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -1,6 +1,7 @@ export LennardJones, - LennardJonesSoftCore + LennardJonesSoftCore, + AshbaughHatch @doc raw""" LennardJones(; cutoff, use_neighbors, lorentz_mixing, weight_special, weight_solute_solvent, @@ -322,3 +323,191 @@ 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, α, λ, p, use_neighbors, lorentz_mixing, weight_special, + weight_solute_solvent, force_units, energy_units, skip_shortcut) + +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}) -λ_{ij}V_{\text{LJ}}(r_{c})+\varepsilon_{ij}(1-λ_{ij}) &, r_{ij}\leq 2^{1/6}σ \\ + λ_{ij}\left(V_{\text{LJ}}(r_{ij}) -V_{\text{LJ}}(r_{c} \right) &, 2^{1/6}σ \leq r_{ij}< r_{c}\\ + 0 &, r_{ij}\gt r_{c} + \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}\left(F_{\text{LJ}}(r_{ij}) &, 2^{1/6}σ \leq r_{ij}< r_{c}\\ + 0 &, r_{ij}\gt r_{c} + \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] +\end{aligned} +``` + +If ``\lambda`` is one this gives the standard [`LennardJones`](@ref) potential. +""" +struct AshbaughHatch{S, C, W, WS, F, E} <: PairwiseInteraction + cutoff::C + use_neighbors::Bool + lorentz_mixing::Bool + weight_special::W + weight_solute_solvent::WS + force_units::F + energy_units::E +end + +function AshbaughHatch(; + cutoff=ShiftedPotentialCutoff(), + use_neighbors=false, + lorentz_mixing=true, + weight_special=1, + weight_solute_solvent=1, + force_units=u"kJ * mol^-1 * nm^-1", + energy_units=u"kJ * mol^-1", + skip_shortcut=false) + return AshbaughHatch{skip_shortcut, typeof(cutoff), typeof(weight_special), + typeof(weight_solute_solvent), typeof(force_units), typeof(energy_units)}( + cutoff, use_neighbors, lorentz_mixing, weight_special, weight_solute_solvent, + force_units, energy_units) +end + +use_neighbors(inter::AshbaughHatch) = inter.use_neighbors + +is_solute(at::Atom) = at.solute +is_solute(at) = false + +function Base.zero(lj::AshbaughHatch{S, C, W, WS, F, E}) where {S, C, W, WS, F, E} + return LennardJones{S, C, W, WS, F, E}( + lj.cutoff, + false, + false, + zero(W), + zero(WS), + lj.force_units, + lj.energy_units, + ) +end + +function Base.:+(l1::AshbaughHatch{S, C, W, WS, F, E}, + l2::AshbaughHatch{S, C, W, WS, F, E}) where {S, C, W, WS, F, E} + return AshbaughHatch{S, C, W, WS, F, E}( + l1.cutoff, + l1.use_neighbors, + l1.lorentz_mixing, + l1.weight_special + l2.weight_special, + l1.weight_solute_solvent + l2.weight_solute_solvent, + l1.force_units, + l1.energy_units, + ) +end + +@inline function force(inter::AshbaughHatch{S, C}, + dr, + coord_i, + coord_j, + atom_i, + atom_j, + boundary, + special::Bool=false) where {S, C} + if !S && (iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) || + iszero_value(atom_i.σ) || iszero_value(atom_j.σ) || + iszero_value(atom_i.λ) || iszero_value(atom_j.λ)) + return ustrip.(zero(coord_i)) * inter.force_units + end + + # Lorentz-Berthelot mixing rules use the arithmetic average for σ + # Otherwise use the geometric average + σ = inter.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ) + λ = inter.lorentz_mixing ? (atom_i.λ + atom_j.λ) / 2 : sqrt(atom_i.λ * atom_j.λ) + + if (is_solute(atom_i) && !is_solute(atom_j)) || (is_solute(atom_j) && !is_solute(atom_i)) + ϵ = inter.weight_solute_solvent * sqrt(atom_i.ϵ * atom_j.ϵ) + else + ϵ = sqrt(atom_i.ϵ * atom_j.ϵ) + end + + cutoff = inter.cutoff + r2 = sum(abs2, dr) + σ2 = σ^2 + params = (σ2, ϵ, λ) + + f = force_divr_with_cutoff(inter, r2, params, cutoff, coord_i, inter.force_units) + if special + return f * dr * inter.weight_special + else + return f * dr + end +end + +function force_divr(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ)) + #six_term = (σ2 * invr2) ^ 3 + if r2 < 2^(1/3)*σ2 + return force_divr(::LennardJones, r2, invr2, (σ2, ϵ))# (24ϵ * invr2) * (2 * six_term ^ 2 - six_term) + else + return λ*force_divr(::LennardJones, r2, invr2, (σ2, ϵ)) #*(24ϵ * invr2) * (2 * six_term ^ 2 - six_term) + end +end + +@inline function potential_energy(inter::AshbaughHatch{S, C}, + dr, + coord_i, + coord_j, + atom_i, + atom_j, + boundary, + special::Bool=false) where {S, C} + if !S && (iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) || + iszero_value(atom_i.σ) || iszero_value(atom_j.σ) || + iszero_value(atom_i.λ) || iszero_value(atom_j.λ)) + return ustrip(zero(coord_i[1])) * inter.energy_units + end + + σ = inter.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ) + λ = inter.lorentz_mixing ? (atom_i.λ + atom_j.λ) / 2 : sqrt(atom_i.λ * atom_j.λ) + + ### probably not needed for Calvados2 use case + if (is_solute(atom_i) && !is_solute(atom_j)) || (is_solute(atom_j) && !is_solute(atom_i)) + ϵ = inter.weight_solute_solvent * sqrt(atom_i.ϵ * atom_j.ϵ) + else + ϵ = sqrt(atom_i.ϵ * atom_j.ϵ) + end + + cutoff = inter.cutoff + r2 = sum(abs2, dr) + σ2 = σ^2 + params = (σ2, ϵ, λ) + + pe = potential_with_cutoff(inter, r2, params, cutoff, coord_i, inter.energy_units) + if special + return pe * inter.weight_special + else + return pe + end +end + +function potential(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ)) + #six_term = (σ2 * invr2) ^ 3 + if r2 < 2^(1/3)*σ2 + return potential(::LennardJones, r2, invr2, (σ2, ϵ)) #4ϵ * (six_term ^ 2 - six_term) + ϵ*(1-λ) + else + return λ* potential(::LennardJones, r2, invr2, (σ2, ϵ)) #4ϵ * (six_term ^ 2 - six_term) + end +end \ No newline at end of file From a7525dd0455baaa645d855ca31d5ccea81ee3c8a Mon Sep 17 00:00:00 2001 From: ywitzky <134135819+ywitzky@users.noreply.github.com> Date: Thu, 26 Sep 2024 14:50:49 +0200 Subject: [PATCH 2/8] added fixed and test for ashbaugh hatch potential --- docs/src/documentation.md | 1 + src/interactions/lennard_jones.jl | 52 ++++++++++++++++++++---------- test/interactions.jl | 53 ++++++++++++++++++++++++++++++- 3 files changed, 88 insertions(+), 18 deletions(-) diff --git a/docs/src/documentation.md b/docs/src/documentation.md index c8731c6f..d25ec664 100644 --- a/docs/src/documentation.md +++ b/docs/src/documentation.md @@ -540,6 +540,7 @@ 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) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 279b001a..29816f63 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -1,7 +1,8 @@ export LennardJones, LennardJonesSoftCore, - AshbaughHatch + AshbaughHatch, + AshbaughHatchAtom @doc raw""" LennardJones(; cutoff, use_neighbors, lorentz_mixing, weight_special, weight_solute_solvent, @@ -340,7 +341,7 @@ V_{\text{LJ}}(r_{ij}) = 4\varepsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}} V_{\text{AH}}(r_{ij}) = \begin{cases} V_{\text{LJ}}(r_{ij}) -λ_{ij}V_{\text{LJ}}(r_{c})+\varepsilon_{ij}(1-λ_{ij}) &, r_{ij}\leq 2^{1/6}σ \\ - λ_{ij}\left(V_{\text{LJ}}(r_{ij}) -V_{\text{LJ}}(r_{c} \right) &, 2^{1/6}σ \leq r_{ij}< r_{c}\\ + λ_{ij}\left(V_{\text{LJ}}(r_{ij}) -V_{\text{LJ}}(r_{c}) \right) &, 2^{1/6}σ \leq r_{ij}< r_{c}\\ 0 &, r_{ij}\gt r_{c} \end{cases} ``` @@ -349,8 +350,8 @@ 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}\left(F_{\text{LJ}}(r_{ij}) &, 2^{1/6}σ \leq r_{ij}< r_{c}\\ + F_{\text{LJ}}(r_{ij}) &, r_{ij} \leq 2^{1/6}σ \\ + λ_{ij}F_{\text{LJ}}(r_{ij}) &, 2^{1/6}σ \leq r_{ij}< r_{c}\\ 0 &, r_{ij}\gt r_{c} \end{cases} ``` @@ -358,7 +359,7 @@ 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] +&= \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} ``` @@ -374,8 +375,30 @@ struct AshbaughHatch{S, C, W, WS, F, E} <: PairwiseInteraction energy_units::E end +struct AshbaughHatchAtom{C, M, S, E, L} + index::Int + charge::C + mass::M + σ::S + ϵ::E + λ::L + solute::Bool +end + +function AshbaughHatchAtom(; + index=1, + charge=0.0, + mass=1.0u"g/mol", + σ=0.0u"nm", + ϵ=0.0u"kJ * mol^-1", + λ=1.0, + solute=false) +return AshbaughHatchAtom(index, charge, mass, σ, ϵ,λ, solute) +end + + function AshbaughHatch(; - cutoff=ShiftedPotentialCutoff(), + cutoff=NoCutoff(), use_neighbors=false, lorentz_mixing=true, weight_special=1, @@ -391,9 +414,6 @@ end use_neighbors(inter::AshbaughHatch) = inter.use_neighbors -is_solute(at::Atom) = at.solute -is_solute(at) = false - function Base.zero(lj::AshbaughHatch{S, C, W, WS, F, E}) where {S, C, W, WS, F, E} return LennardJones{S, C, W, WS, F, E}( lj.cutoff, @@ -457,12 +477,11 @@ end end end -function force_divr(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ)) - #six_term = (σ2 * invr2) ^ 3 +@inline function force_divr(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ)) if r2 < 2^(1/3)*σ2 - return force_divr(::LennardJones, r2, invr2, (σ2, ϵ))# (24ϵ * invr2) * (2 * six_term ^ 2 - six_term) + return force_divr(LennardJones(), r2, invr2, (σ2, ϵ)) else - return λ*force_divr(::LennardJones, r2, invr2, (σ2, ϵ)) #*(24ϵ * invr2) * (2 * six_term ^ 2 - six_term) + return λ*force_divr(LennardJones(), r2, invr2, (σ2, ϵ)) end end @@ -503,11 +522,10 @@ end end end -function potential(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ)) - #six_term = (σ2 * invr2) ^ 3 +@inline function potential(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ)) if r2 < 2^(1/3)*σ2 - return potential(::LennardJones, r2, invr2, (σ2, ϵ)) #4ϵ * (six_term ^ 2 - six_term) + ϵ*(1-λ) + return potential(LennardJones(), r2, invr2, (σ2, ϵ)) + ϵ*(1-λ) else - return λ* potential(::LennardJones, r2, invr2, (σ2, ϵ)) #4ϵ * (six_term ^ 2 - six_term) + 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 e8f24897..7b9e5909 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -83,7 +83,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"Å" @@ -386,4 +386,55 @@ -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() + + @test isapprox( + force(inter, dr12, c1, c2, 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, c1, c3, 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, c1, c2, AH_a1, AH_a1, boundary), + 0.0u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + @test isapprox( + potential_energy(inter, dr13, c1, c3, 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, c1, c3, AH_a1, AH_a1, boundary), + -0.058520865u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + + @test isapprox( + force(inter, dr13, c1, c3, AH_a1, AH_a1, boundary), + SVector(-0.68775486946, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-9u"kJ * mol^-1 * nm^-1", + ) + + c2 = SVector(1.28, 1.0, 1.0)u"nm" + dr12 = vector(c1, c2, boundary) + @test isapprox( + potential_energy(inter, dr12, c1, c2, AH_a1, AH_a1, boundary), + 0.7205987916u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + + @test isapprox( + force(inter, dr12, c1, c2, AH_a1, AH_a1, boundary), + SVector(52.5306754422, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + atol=1e-9u"kJ * mol^-1 * nm^-1", + ) end From c48dfddf465fd47d94b100d16d06ab284e8f6dd2 Mon Sep 17 00:00:00 2001 From: ywitzky <134135819+ywitzky@users.noreply.github.com> Date: Tue, 1 Oct 2024 16:47:32 +0200 Subject: [PATCH 3/8] Added Yukawa potential --- docs/src/documentation.md | 1 + src/interactions/coulomb.jl | 127 +++++++++++++++++++++++++++++++++++- test/interactions.jl | 23 +++++++ 3 files changed, 150 insertions(+), 1 deletion(-) diff --git a/docs/src/documentation.md b/docs/src/documentation.md index c8731c6f..3bb4ee6c 100644 --- a/docs/src/documentation.md +++ b/docs/src/documentation.md @@ -546,6 +546,7 @@ The available pairwise interactions are: - [`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 595d160e..8d74ea4a 100644 --- a/src/interactions/coulomb.jl +++ b/src/interactions/coulomb.jl @@ -1,7 +1,8 @@ export Coulomb, CoulombSoftCore, - CoulombReactionField + CoulombReactionField, + Yukawa @doc raw""" Coulomb(; cutoff, use_neighbors, weight_special, coulomb_const, force_units, energy_units) @@ -377,3 +378,127 @@ end return pe end end + + + + + +@doc raw""" + Yukawa(; cutoff, use_neighbors, weight_special, coulomb_const, force_units, energy_units) + +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} +``` +""" +struct Yukawa{C, W, T, F, E, D} <: PairwiseInteraction + cutoff::C + use_neighbors::Bool + weight_special::W + coulomb_const::T + force_units::F + energy_units::E + kappa::D +end + +#const coulombconst = 138.93545764u"kJ * mol^-1 * nm" # 1 / 4πϵ0 + +function Yukawa(; + cutoff=NoCutoff(), + use_neighbors=false, + weight_special=1, + coulomb_const=coulombconst, + force_units=u"kJ * mol^-1 * nm^-1", + energy_units=u"kJ * mol^-1", + kappa=1.0*u"nm^-1") + return Yukawa{typeof(cutoff), typeof(weight_special), typeof(coulomb_const), typeof(force_units), typeof(energy_units),typeof(kappa)}( + cutoff, use_neighbors, weight_special, coulomb_const, force_units, energy_units, kappa) +end + +use_neighbors(inter::Yukawa) = inter.use_neighbors + +function Base.zero(yukawa::Yukawa{C, W, T, F, E,D}) where {C, W, T, F, E,D} + return Yukawa{C, W, T, F, E,D}( + yukawa.cutoff, + false, + zero(W), + zero(T), + yukawa.force_units, + yukawa.energy_units, + 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 + c2.coulomb_const, + c1.force_units, + c1.energy_units, + c1.kappa + ) +end + +@inline function force(inter::Yukawa{C}, + dr, + coord_i, + coord_j, + atom_i, + atom_j, + boundary, + special::Bool=false) 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, coord_i, inter.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.0) +end + +@inline function potential_energy(inter::Yukawa{C}, + dr, + coord_i, + coord_j, + atom_i, + atom_j, + boundary, + special::Bool=false) 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, coord_i, inter.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/test/interactions.jl b/test/interactions.jl index e8f24897..f207985e 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -159,6 +159,29 @@ atol=1e-5u"kJ * mol^-1", ) + + inter = Yukawa(; kappa=1.0u"nm^-1") + @test isapprox( + force(inter, dr12, c1, c2, a1, a1, boundary), + 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, c1, c3, a1, a1, boundary), + SVector(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, c1, c2, a1, a1, boundary), + 343.08639592583785u"kJ * mol^-1"; + atol=1e-5u"kJ * mol^-1", + ) + @test isapprox( + potential_energy(inter, dr13, c1, c3, a1, a1, boundary), + 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") From 3f885867c3477d3590622be43006742fdbeb22a0 Mon Sep 17 00:00:00 2001 From: ywitzky <134135819+ywitzky@users.noreply.github.com> Date: Tue, 1 Oct 2024 16:47:32 +0200 Subject: [PATCH 4/8] Added Yukawa potential Signed-off-by: ywitzky <134135819+ywitzky@users.noreply.github.com> --- docs/src/documentation.md | 1 + src/interactions/coulomb.jl | 127 +++++++++++++++++++++++++++++++++++- test/interactions.jl | 23 +++++++ 3 files changed, 150 insertions(+), 1 deletion(-) diff --git a/docs/src/documentation.md b/docs/src/documentation.md index d25ec664..a3e01e8f 100644 --- a/docs/src/documentation.md +++ b/docs/src/documentation.md @@ -547,6 +547,7 @@ The available pairwise interactions are: - [`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 595d160e..8d74ea4a 100644 --- a/src/interactions/coulomb.jl +++ b/src/interactions/coulomb.jl @@ -1,7 +1,8 @@ export Coulomb, CoulombSoftCore, - CoulombReactionField + CoulombReactionField, + Yukawa @doc raw""" Coulomb(; cutoff, use_neighbors, weight_special, coulomb_const, force_units, energy_units) @@ -377,3 +378,127 @@ end return pe end end + + + + + +@doc raw""" + Yukawa(; cutoff, use_neighbors, weight_special, coulomb_const, force_units, energy_units) + +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} +``` +""" +struct Yukawa{C, W, T, F, E, D} <: PairwiseInteraction + cutoff::C + use_neighbors::Bool + weight_special::W + coulomb_const::T + force_units::F + energy_units::E + kappa::D +end + +#const coulombconst = 138.93545764u"kJ * mol^-1 * nm" # 1 / 4πϵ0 + +function Yukawa(; + cutoff=NoCutoff(), + use_neighbors=false, + weight_special=1, + coulomb_const=coulombconst, + force_units=u"kJ * mol^-1 * nm^-1", + energy_units=u"kJ * mol^-1", + kappa=1.0*u"nm^-1") + return Yukawa{typeof(cutoff), typeof(weight_special), typeof(coulomb_const), typeof(force_units), typeof(energy_units),typeof(kappa)}( + cutoff, use_neighbors, weight_special, coulomb_const, force_units, energy_units, kappa) +end + +use_neighbors(inter::Yukawa) = inter.use_neighbors + +function Base.zero(yukawa::Yukawa{C, W, T, F, E,D}) where {C, W, T, F, E,D} + return Yukawa{C, W, T, F, E,D}( + yukawa.cutoff, + false, + zero(W), + zero(T), + yukawa.force_units, + yukawa.energy_units, + 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 + c2.coulomb_const, + c1.force_units, + c1.energy_units, + c1.kappa + ) +end + +@inline function force(inter::Yukawa{C}, + dr, + coord_i, + coord_j, + atom_i, + atom_j, + boundary, + special::Bool=false) 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, coord_i, inter.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.0) +end + +@inline function potential_energy(inter::Yukawa{C}, + dr, + coord_i, + coord_j, + atom_i, + atom_j, + boundary, + special::Bool=false) 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, coord_i, inter.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/test/interactions.jl b/test/interactions.jl index 7b9e5909..03ab0d7a 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -159,6 +159,29 @@ atol=1e-5u"kJ * mol^-1", ) + + inter = Yukawa(; kappa=1.0u"nm^-1") + @test isapprox( + force(inter, dr12, c1, c2, a1, a1, boundary), + 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, c1, c3, a1, a1, boundary), + SVector(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, c1, c2, a1, a1, boundary), + 343.08639592583785u"kJ * mol^-1"; + atol=1e-5u"kJ * mol^-1", + ) + @test isapprox( + potential_energy(inter, dr13, c1, c3, a1, a1, boundary), + 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") From c51c48ffa5a62b6850dd8c23ad35218b4f7785e5 Mon Sep 17 00:00:00 2001 From: ywitzky <134135819+ywitzky@users.noreply.github.com> Date: Thu, 31 Oct 2024 20:14:38 +0100 Subject: [PATCH 5/8] updated implementation to match recent changes --- src/interactions/coulomb.jl | 54 ++++++------- src/interactions/lennard_jones.jl | 130 ++++++++++++++---------------- test/interactions.jl | 28 ++++--- 3 files changed, 100 insertions(+), 112 deletions(-) diff --git a/src/interactions/coulomb.jl b/src/interactions/coulomb.jl index 1fe696f4..55703e68 100644 --- a/src/interactions/coulomb.jl +++ b/src/interactions/coulomb.jl @@ -362,18 +362,19 @@ and the force on each atom by 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} ``` """ -struct Yukawa{C, W, T, F, E, D} <: PairwiseInteraction - cutoff::C - use_neighbors::Bool - weight_special::W - coulomb_const::T - force_units::F - energy_units::E - kappa::D -end #const coulombconst = 138.93545764u"kJ * mol^-1 * nm" # 1 / 4πϵ0 +@kwdef struct Yukawa{C, W, T, D} # <: PairwiseInteraction + cutoff::C= NoCutoff() + use_neighbors::Bool= false + weight_special::W =1 + coulomb_const::T = coulomb_const + kappa::D = kappa=1.0*u"nm^-1" +end + + +#= function Yukawa(; cutoff=NoCutoff(), use_neighbors=false, @@ -384,18 +385,16 @@ function Yukawa(; kappa=1.0*u"nm^-1") return Yukawa{typeof(cutoff), typeof(weight_special), typeof(coulomb_const), typeof(force_units), typeof(energy_units),typeof(kappa)}( cutoff, use_neighbors, weight_special, coulomb_const, force_units, energy_units, kappa) -end +end=# use_neighbors(inter::Yukawa) = inter.use_neighbors -function Base.zero(yukawa::Yukawa{C, W, T, F, E,D}) where {C, W, T, F, E,D} - return Yukawa{C, W, T, F, E,D}( +function Base.zero(yukawa::Yukawa{C, W, T ,D}) where {C, W, T, D} + return Yukawa{C, W, T, D}( yukawa.cutoff, - false, - zero(W), - zero(T), - yukawa.force_units, - yukawa.energy_units, + yukawa.use_neighbors, + yukawa.weight_special, + yukawa.coulomb_const, yukawa.kappa ) end @@ -405,21 +404,18 @@ function Base.:+(c1::Yukawa, c2::Yukawa) c1.cutoff, c1.use_neighbors, c1.weight_special + c2.weight_special, - c1.coulomb_const + c2.coulomb_const, - c1.force_units, - c1.energy_units, + c1.coulomb_const,# + c2.coulomb_const, c1.kappa ) end @inline function force(inter::Yukawa{C}, dr, - coord_i, - coord_j, atom_i, atom_j, - boundary, - special::Bool=false) where C + 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 @@ -427,7 +423,7 @@ end kappa = inter.kappa params = (coulomb_const, qi, qj, kappa) - f = force_divr_with_cutoff(inter, r2, params, cutoff, coord_i, inter.force_units) + f = force_divr_with_cutoff(inter, r2, params, cutoff, force_units) if special return f * dr * inter.weight_special else @@ -442,19 +438,17 @@ end @inline function potential_energy(inter::Yukawa{C}, dr, - coord_i, - coord_j, atom_i, atom_j, - boundary, - special::Bool=false) where C + 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, coord_i, inter.energy_units) + pe = potential_with_cutoff(inter, r2, params, cutoff, energy_units) if special return pe * inter.weight_special else diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index c0a38120..2082b84c 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -15,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 @@ -323,26 +331,29 @@ where If ``\lambda`` is one this gives the standard [`LennardJones`](@ref) potential. """ -struct AshbaughHatch{S, C, W, WS, F, E} <: PairwiseInteraction - cutoff::C - use_neighbors::Bool - lorentz_mixing::Bool - weight_special::W - weight_solute_solvent::WS - force_units::F - energy_units::E +@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 -struct AshbaughHatchAtom{C, M, S, E, L} - index::Int - charge::C - mass::M - σ::S - ϵ::E - λ::L - solute::Bool + + +@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 + solute::Bool = false end +#= function AshbaughHatchAtom(; index=1, charge=0.0, @@ -352,9 +363,9 @@ function AshbaughHatchAtom(; λ=1.0, solute=false) return AshbaughHatchAtom(index, charge, mass, σ, ϵ,λ, solute) -end - +end =# +#= function AshbaughHatch(; cutoff=NoCutoff(), use_neighbors=false, @@ -368,66 +379,57 @@ function AshbaughHatch(; typeof(weight_solute_solvent), typeof(force_units), typeof(energy_units)}( cutoff, use_neighbors, lorentz_mixing, weight_special, weight_solute_solvent, force_units, energy_units) -end +end=# use_neighbors(inter::AshbaughHatch) = inter.use_neighbors -function Base.zero(lj::AshbaughHatch{S, C, W, WS, F, E}) where {S, C, W, WS, F, E} - return LennardJones{S, C, W, WS, F, E}( +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, - false, - false, + lj,use_neighbors, + lj.shortcut, + lj.σ_mixing, + lj.ϵ_mixing, + lj.λ_mixing, zero(W), - zero(WS), - lj.force_units, - lj.energy_units, ) end -function Base.:+(l1::AshbaughHatch{S, C, W, WS, F, E}, - l2::AshbaughHatch{S, C, W, WS, F, E}) where {S, C, W, WS, F, E} - return AshbaughHatch{S, C, W, WS, F, E}( +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.lorentz_mixing, + l1.shortcut, + l1.σ_mixing, + l1.ϵ_mixing, + l1.λ_mixing, l1.weight_special + l2.weight_special, - l1.weight_solute_solvent + l2.weight_solute_solvent, - l1.force_units, - l1.energy_units, ) end @inline function force(inter::AshbaughHatch{S, C}, dr, - coord_i, - coord_j, atom_i, atom_j, - boundary, - special::Bool=false) where {S, C} - if !S && (iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) || - iszero_value(atom_i.σ) || iszero_value(atom_j.σ) || - iszero_value(atom_i.λ) || iszero_value(atom_j.λ)) - return ustrip.(zero(coord_i)) * inter.force_units + 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.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ) - λ = inter.lorentz_mixing ? (atom_i.λ + atom_j.λ) / 2 : sqrt(atom_i.λ * atom_j.λ) - - if (is_solute(atom_i) && !is_solute(atom_j)) || (is_solute(atom_j) && !is_solute(atom_i)) - ϵ = inter.weight_solute_solvent * sqrt(atom_i.ϵ * atom_j.ϵ) - else - ϵ = sqrt(atom_i.ϵ * atom_j.ϵ) - end + ϵ = 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, coord_i, inter.force_units) + f = force_divr_with_cutoff(inter, r2, params, cutoff,force_units) if special return f * dr * inter.weight_special else @@ -445,34 +447,24 @@ end @inline function potential_energy(inter::AshbaughHatch{S, C}, dr, - coord_i, - coord_j, atom_i, atom_j, - boundary, - special::Bool=false) where {S, C} - if !S && (iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) || - iszero_value(atom_i.σ) || iszero_value(atom_j.σ) || - iszero_value(atom_i.λ) || iszero_value(atom_j.λ)) - return ustrip(zero(coord_i[1])) * inter.energy_units - end - - σ = inter.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ) - λ = inter.lorentz_mixing ? (atom_i.λ + atom_j.λ) / 2 : sqrt(atom_i.λ * atom_j.λ) - - ### probably not needed for Calvados2 use case - if (is_solute(atom_i) && !is_solute(atom_j)) || (is_solute(atom_j) && !is_solute(atom_i)) - ϵ = inter.weight_solute_solvent * sqrt(atom_i.ϵ * atom_j.ϵ) - else - ϵ = sqrt(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, coord_i, inter.energy_units) + pe = potential_with_cutoff(inter, r2, params, cutoff,energy_units) if special return pe * inter.weight_special else diff --git a/test/interactions.jl b/test/interactions.jl index a3139874..66649d13 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -76,6 +76,7 @@ 0.0253410816u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) + struct BuckinghamAtom{M, TA, TB, TC} mass::M @@ -160,24 +161,24 @@ ) - inter = Yukawa(; kappa=1.0u"nm^-1") + inter = Yukawa() @test isapprox( - force(inter, dr12, c1, c2, a1, a1, boundary), + 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, c1, c3, a1, a1, boundary), + 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( - potential_energy(inter, dr12, c1, c2, a1, a1, boundary), + potential_energy(inter, dr12, a1, a1), 343.08639592583785u"kJ * mol^-1"; atol=1e-5u"kJ * mol^-1", ) @test isapprox( - potential_energy(inter, dr13, c1, c3, a1, a1, boundary), + potential_energy(inter, dr13, a1, a1), 232.8280565063u"kJ * mol^-1"; atol=1e-5u"kJ * mol^-1", ) @@ -413,35 +414,36 @@ inter = AshbaughHatch() @test isapprox( - force(inter, dr12, c1, c2, AH_a1, AH_a1, boundary), + 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, c1, c3, AH_a1, AH_a1, boundary), + 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, c1, c2, AH_a1, AH_a1, boundary), + 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, c1, c3, AH_a1, AH_a1, boundary), + 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, c1, c3, AH_a1, AH_a1, boundary), + potential_energy(inter, dr13, AH_a1, AH_a1), -0.058520865u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) @test isapprox( - force(inter, dr13, c1, c3, AH_a1, AH_a1, boundary), + 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", ) @@ -449,13 +451,13 @@ c2 = SVector(1.28, 1.0, 1.0)u"nm" dr12 = vector(c1, c2, boundary) @test isapprox( - potential_energy(inter, dr12, c1, c2, AH_a1, AH_a1, boundary), + potential_energy(inter, dr12, AH_a1, AH_a1), 0.7205987916u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) @test isapprox( - force(inter, dr12, c1, c2, AH_a1, AH_a1, boundary), + force(inter, dr12, 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", ) From 6952e99607c46c62bd2c14c675339b28a454dc08 Mon Sep 17 00:00:00 2001 From: ywitzky <134135819+ywitzky@users.noreply.github.com> Date: Thu, 31 Oct 2024 20:14:38 +0100 Subject: [PATCH 6/8] updated implementation to match recent changes --- src/interactions/coulomb.jl | 61 ++++--------- src/interactions/lennard_jones.jl | 147 +++++++++++------------------- test/interactions.jl | 28 +++--- 3 files changed, 89 insertions(+), 147 deletions(-) diff --git a/src/interactions/coulomb.jl b/src/interactions/coulomb.jl index 1fe696f4..e8c5925e 100644 --- a/src/interactions/coulomb.jl +++ b/src/interactions/coulomb.jl @@ -362,40 +362,22 @@ and the force on each atom by 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} ``` """ -struct Yukawa{C, W, T, F, E, D} <: PairwiseInteraction - cutoff::C - use_neighbors::Bool - weight_special::W - coulomb_const::T - force_units::F - energy_units::E - kappa::D -end - -#const coulombconst = 138.93545764u"kJ * mol^-1 * nm" # 1 / 4πϵ0 - -function Yukawa(; - cutoff=NoCutoff(), - use_neighbors=false, - weight_special=1, - coulomb_const=coulombconst, - force_units=u"kJ * mol^-1 * nm^-1", - energy_units=u"kJ * mol^-1", - kappa=1.0*u"nm^-1") - return Yukawa{typeof(cutoff), typeof(weight_special), typeof(coulomb_const), typeof(force_units), typeof(energy_units),typeof(kappa)}( - cutoff, use_neighbors, weight_special, coulomb_const, force_units, energy_units, kappa) +@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 = kappa=1.0*u"nm^-1" end use_neighbors(inter::Yukawa) = inter.use_neighbors -function Base.zero(yukawa::Yukawa{C, W, T, F, E,D}) where {C, W, T, F, E,D} - return Yukawa{C, W, T, F, E,D}( +function Base.zero(yukawa::Yukawa{C, W, T ,D}) where {C, W, T, D} + return Yukawa{C, W, T, D}( yukawa.cutoff, - false, - zero(W), - zero(T), - yukawa.force_units, - yukawa.energy_units, + yukawa.use_neighbors, + yukawa.weight_special, + yukawa.coulomb_const, yukawa.kappa ) end @@ -405,21 +387,18 @@ function Base.:+(c1::Yukawa, c2::Yukawa) c1.cutoff, c1.use_neighbors, c1.weight_special + c2.weight_special, - c1.coulomb_const + c2.coulomb_const, - c1.force_units, - c1.energy_units, + c1.coulomb_const, c1.kappa ) end @inline function force(inter::Yukawa{C}, dr, - coord_i, - coord_j, atom_i, atom_j, - boundary, - special::Bool=false) where C + 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 @@ -427,7 +406,7 @@ end kappa = inter.kappa params = (coulomb_const, qi, qj, kappa) - f = force_divr_with_cutoff(inter, r2, params, cutoff, coord_i, inter.force_units) + f = force_divr_with_cutoff(inter, r2, params, cutoff, force_units) if special return f * dr * inter.weight_special else @@ -442,19 +421,17 @@ end @inline function potential_energy(inter::Yukawa{C}, dr, - coord_i, - coord_j, atom_i, atom_j, - boundary, - special::Bool=false) where C + 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, coord_i, inter.energy_units) + pe = potential_with_cutoff(inter, r2, params, cutoff, energy_units) if special return pe * inter.weight_special else diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index c0a38120..11b66438 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -15,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 @@ -323,111 +331,76 @@ where If ``\lambda`` is one this gives the standard [`LennardJones`](@ref) potential. """ -struct AshbaughHatch{S, C, W, WS, F, E} <: PairwiseInteraction - cutoff::C - use_neighbors::Bool - lorentz_mixing::Bool - weight_special::W - weight_solute_solvent::WS - force_units::F - energy_units::E -end - -struct AshbaughHatchAtom{C, M, S, E, L} - index::Int - charge::C - mass::M - σ::S - ϵ::E - λ::L - solute::Bool -end - -function AshbaughHatchAtom(; - index=1, - charge=0.0, - mass=1.0u"g/mol", - σ=0.0u"nm", - ϵ=0.0u"kJ * mol^-1", - λ=1.0, - solute=false) -return AshbaughHatchAtom(index, charge, mass, σ, ϵ,λ, solute) +@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 -function AshbaughHatch(; - cutoff=NoCutoff(), - use_neighbors=false, - lorentz_mixing=true, - weight_special=1, - weight_solute_solvent=1, - force_units=u"kJ * mol^-1 * nm^-1", - energy_units=u"kJ * mol^-1", - skip_shortcut=false) - return AshbaughHatch{skip_shortcut, typeof(cutoff), typeof(weight_special), - typeof(weight_solute_solvent), typeof(force_units), typeof(energy_units)}( - cutoff, use_neighbors, lorentz_mixing, weight_special, weight_solute_solvent, - force_units, energy_units) +@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 + solute::Bool = false end use_neighbors(inter::AshbaughHatch) = inter.use_neighbors -function Base.zero(lj::AshbaughHatch{S, C, W, WS, F, E}) where {S, C, W, WS, F, E} - return LennardJones{S, C, W, WS, F, E}( +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, - false, - false, + lj,use_neighbors, + lj.shortcut, + lj.σ_mixing, + lj.ϵ_mixing, + lj.λ_mixing, zero(W), - zero(WS), - lj.force_units, - lj.energy_units, ) end -function Base.:+(l1::AshbaughHatch{S, C, W, WS, F, E}, - l2::AshbaughHatch{S, C, W, WS, F, E}) where {S, C, W, WS, F, E} - return AshbaughHatch{S, C, W, WS, F, E}( +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.lorentz_mixing, + l1.shortcut, + l1.σ_mixing, + l1.ϵ_mixing, + l1.λ_mixing, l1.weight_special + l2.weight_special, - l1.weight_solute_solvent + l2.weight_solute_solvent, - l1.force_units, - l1.energy_units, ) end @inline function force(inter::AshbaughHatch{S, C}, dr, - coord_i, - coord_j, atom_i, atom_j, - boundary, - special::Bool=false) where {S, C} - if !S && (iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) || - iszero_value(atom_i.σ) || iszero_value(atom_j.σ) || - iszero_value(atom_i.λ) || iszero_value(atom_j.λ)) - return ustrip.(zero(coord_i)) * inter.force_units + 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.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ) - λ = inter.lorentz_mixing ? (atom_i.λ + atom_j.λ) / 2 : sqrt(atom_i.λ * atom_j.λ) - - if (is_solute(atom_i) && !is_solute(atom_j)) || (is_solute(atom_j) && !is_solute(atom_i)) - ϵ = inter.weight_solute_solvent * sqrt(atom_i.ϵ * atom_j.ϵ) - else - ϵ = sqrt(atom_i.ϵ * atom_j.ϵ) - end + ϵ = 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, coord_i, inter.force_units) + f = force_divr_with_cutoff(inter, r2, params, cutoff,force_units) if special return f * dr * inter.weight_special else @@ -445,34 +418,24 @@ end @inline function potential_energy(inter::AshbaughHatch{S, C}, dr, - coord_i, - coord_j, atom_i, atom_j, - boundary, - special::Bool=false) where {S, C} - if !S && (iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) || - iszero_value(atom_i.σ) || iszero_value(atom_j.σ) || - iszero_value(atom_i.λ) || iszero_value(atom_j.λ)) - return ustrip(zero(coord_i[1])) * inter.energy_units - end - - σ = inter.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ) - λ = inter.lorentz_mixing ? (atom_i.λ + atom_j.λ) / 2 : sqrt(atom_i.λ * atom_j.λ) - - ### probably not needed for Calvados2 use case - if (is_solute(atom_i) && !is_solute(atom_j)) || (is_solute(atom_j) && !is_solute(atom_i)) - ϵ = inter.weight_solute_solvent * sqrt(atom_i.ϵ * atom_j.ϵ) - else - ϵ = sqrt(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, coord_i, inter.energy_units) + pe = potential_with_cutoff(inter, r2, params, cutoff,energy_units) if special return pe * inter.weight_special else diff --git a/test/interactions.jl b/test/interactions.jl index a3139874..66649d13 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -76,6 +76,7 @@ 0.0253410816u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) + struct BuckinghamAtom{M, TA, TB, TC} mass::M @@ -160,24 +161,24 @@ ) - inter = Yukawa(; kappa=1.0u"nm^-1") + inter = Yukawa() @test isapprox( - force(inter, dr12, c1, c2, a1, a1, boundary), + 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, c1, c3, a1, a1, boundary), + 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( - potential_energy(inter, dr12, c1, c2, a1, a1, boundary), + potential_energy(inter, dr12, a1, a1), 343.08639592583785u"kJ * mol^-1"; atol=1e-5u"kJ * mol^-1", ) @test isapprox( - potential_energy(inter, dr13, c1, c3, a1, a1, boundary), + potential_energy(inter, dr13, a1, a1), 232.8280565063u"kJ * mol^-1"; atol=1e-5u"kJ * mol^-1", ) @@ -413,35 +414,36 @@ inter = AshbaughHatch() @test isapprox( - force(inter, dr12, c1, c2, AH_a1, AH_a1, boundary), + 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, c1, c3, AH_a1, AH_a1, boundary), + 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, c1, c2, AH_a1, AH_a1, boundary), + 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, c1, c3, AH_a1, AH_a1, boundary), + 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, c1, c3, AH_a1, AH_a1, boundary), + potential_energy(inter, dr13, AH_a1, AH_a1), -0.058520865u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) @test isapprox( - force(inter, dr13, c1, c3, AH_a1, AH_a1, boundary), + 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", ) @@ -449,13 +451,13 @@ c2 = SVector(1.28, 1.0, 1.0)u"nm" dr12 = vector(c1, c2, boundary) @test isapprox( - potential_energy(inter, dr12, c1, c2, AH_a1, AH_a1, boundary), + potential_energy(inter, dr12, AH_a1, AH_a1), 0.7205987916u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) @test isapprox( - force(inter, dr12, c1, c2, AH_a1, AH_a1, boundary), + force(inter, dr12, 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", ) From 3ec0feabd7e7452b558a9d63d383e8922ca5cad4 Mon Sep 17 00:00:00 2001 From: ywitzky <134135819+ywitzky@users.noreply.github.com> Date: Wed, 6 Nov 2024 17:33:49 +0100 Subject: [PATCH 7/8] added review requests --- src/interactions/coulomb.jl | 6 +++--- src/interactions/lennard_jones.jl | 12 ++++-------- test/interactions.jl | 9 ++++----- 3 files changed, 11 insertions(+), 16 deletions(-) diff --git a/src/interactions/coulomb.jl b/src/interactions/coulomb.jl index e8c5925e..6c040864 100644 --- a/src/interactions/coulomb.jl +++ b/src/interactions/coulomb.jl @@ -347,7 +347,7 @@ end @doc raw""" - Yukawa(; cutoff, use_neighbors, weight_special, coulomb_const, force_units, energy_units) + Yukawa(; cutoff, use_neighbors, weight_special, coulomb_const, kappa) The Yukawa electrostatic interaction between two atoms. @@ -367,7 +367,7 @@ F(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}^2} \exp(-\kappa r_{ij})\le use_neighbors::Bool= false weight_special::W =1 coulomb_const::T = coulomb_const - kappa::D = kappa=1.0*u"nm^-1" + kappa::D = 1.0*u"nm^-1" end use_neighbors(inter::Yukawa) = inter.use_neighbors @@ -416,7 +416,7 @@ 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.0) + return (coulomb_const * qi * qj) / r^3 *exp(-kappa*r) * (kappa*r+1) end @inline function potential_energy(inter::Yukawa{C}, diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 11b66438..0aa8ef35 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -293,8 +293,7 @@ end @doc raw""" - AshbaughHatch(; cutoff, α, λ, p, use_neighbors, lorentz_mixing, weight_special, - weight_solute_solvent, force_units, energy_units, skip_shortcut) + 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. @@ -306,9 +305,8 @@ V_{\text{LJ}}(r_{ij}) = 4\varepsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}} ```math V_{\text{AH}}(r_{ij}) = \begin{cases} - V_{\text{LJ}}(r_{ij}) -λ_{ij}V_{\text{LJ}}(r_{c})+\varepsilon_{ij}(1-λ_{ij}) &, r_{ij}\leq 2^{1/6}σ \\ - λ_{ij}\left(V_{\text{LJ}}(r_{ij}) -V_{\text{LJ}}(r_{c}) \right) &, 2^{1/6}σ \leq r_{ij}< r_{c}\\ - 0 &, r_{ij}\gt r_{c} + 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} ``` @@ -317,8 +315,7 @@ and the force on each atom by \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}< r_{c}\\ - 0 &, r_{ij}\gt r_{c} + λ_{ij}F_{\text{LJ}}(r_{ij}) &, 2^{1/6}σ \leq r_{ij} \end{cases} ``` where @@ -349,7 +346,6 @@ end σ::S =0.0u"nm" ϵ::E =0.0u"kJ * mol^-1" λ::L=1.0 - solute::Bool = false end use_neighbors(inter::AshbaughHatch) = inter.use_neighbors diff --git a/test/interactions.jl b/test/interactions.jl index 66649d13..4eb9cb3d 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -76,7 +76,6 @@ 0.0253410816u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) - struct BuckinghamAtom{M, TA, TB, TC} mass::M @@ -448,16 +447,16 @@ atol=1e-9u"kJ * mol^-1 * nm^-1", ) - c2 = SVector(1.28, 1.0, 1.0)u"nm" - dr12 = vector(c1, c2, boundary) + c4 = SVector(1.28, 1.0, 1.0)u"nm" + dr14 = vector(c1, c2, boundary) @test isapprox( - potential_energy(inter, dr12, AH_a1, AH_a1), + potential_energy(inter, dr14, AH_a1, AH_a1), 0.7205987916u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) @test isapprox( - force(inter, dr12, AH_a1, AH_a1), + 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", ) From 3c4f6ac087018b8cba3f8260cb6f017f72a83898 Mon Sep 17 00:00:00 2001 From: ywitzky <134135819+ywitzky@users.noreply.github.com> Date: Wed, 6 Nov 2024 19:43:26 +0100 Subject: [PATCH 8/8] increase test coverage for forces and potentials --- test/interactions.jl | 56 ++++++++++++++++++++++++++++++++++++++------ 1 file changed, 49 insertions(+), 7 deletions(-) diff --git a/test/interactions.jl b/test/interactions.jl index 4eb9cb3d..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), @@ -159,8 +182,7 @@ atol=1e-5u"kJ * mol^-1", ) - - inter = Yukawa() + 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"; @@ -171,6 +193,11 @@ 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"; @@ -181,7 +208,11 @@ 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") @@ -408,9 +439,8 @@ -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() + inter = AshbaughHatch(weight_special=0.5) @test isapprox( force(inter, dr12, AH_a1, AH_a1, boundary), @@ -432,7 +462,7 @@ -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( @@ -448,7 +478,7 @@ ) c4 = SVector(1.28, 1.0, 1.0)u"nm" - dr14 = vector(c1, c2, boundary) + dr14 = vector(c1, c4, boundary) @test isapprox( potential_energy(inter, dr14, AH_a1, AH_a1), 0.7205987916u"kJ * mol^-1"; @@ -460,4 +490,16 @@ 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