From 1009e4da4d338ee41f82aed39fdb09ebac087608 Mon Sep 17 00:00:00 2001 From: Immanuel Felix Sulzer <100942429+Immi000@users.noreply.github.com> Date: Mon, 3 Nov 2025 16:24:07 +0100 Subject: [PATCH 1/8] update LJBeutler to more general implementation --- src/interactions/lennard_jones.jl | 308 +++++++++++++++++++++--------- 1 file changed, 214 insertions(+), 94 deletions(-) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 87f7630a..9ca368e2 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -140,142 +140,262 @@ function pairwise_pe(::LennardJones, r, (σ2, ϵ)) end @doc raw""" - LennardJonesSoftCoreBeutler(; cutoff, α, λ, use_neighbors, shortcut, σ_mixing, - ϵ_mixing, weight_special) - -The Lennard-Jones 6-12 interaction between two atoms with a soft core, used for -the appearing and disappearing of atoms. - -See [Beutler et al. 1994](https://doi.org/10.1016/0009-2614(94)00397-1). -The potential energy is defined as -```math -V(r_{ij}) = \frac{C^{(12)}}{r_{LJ}^{12}} - \frac{C^{(6)}}{r_{LJ}^{6}} -``` -and the force on each atom by -```math -\vec{F}_i = ((\frac{12C^{(12)}}{r_{LJ}^{13}} - \frac{6C^{(6)}}{r_{LJ}^7})(\frac{r_{ij}}{r_{LJ}})^5) \frac{\vec{r_{ij}}}{r_{ij}} -``` -where -```math -r_{LJ} = (\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}}+r^6)^{1/6} -``` -and -```math -C^{(12)} = 4\epsilon\sigma^{12} -C^{(6)} = 4\epsilon\sigma^{6} -``` - -If ``\lambda`` is 1.0, this gives the standard [`LennardJones`](@ref) potential and means -the atom is fully turned on. -If ``\lambda`` is zero the interaction is turned off. -``\alpha`` determines the strength of softening the function. + LennardJonesSoftCoreBeutler(; λ, α, p, cutoff, use_neighbors, inter_state_a, inter_state_b, shortcut, weight_special) + + The Lennard-Jones soft core interaction as described by + [Beutler et al. 1994](https://doi.org/10.1016/0009-2614(94)00397-1). This interaction is a linear + interpolation between two Lennard-Jones interactions A and B controlled by a parameter λ. The + parameters `inter_state_a` and `inter_state_b` can be a named tuple of `σ_mixing` and `ϵ_mixing` + for the interaction in the respective state or `nothing` for no interaction in that state. + + The potential energy is defined as + ````math + V(r)&=(1-\lambda)V^A(r_A)+λV^B(r_B) \\ + r_A&=\left(\alpha\sigma_A^6\lambda^p+r^6\right)^\frac{1}{6} \\ + r_B&=\left(\alpha\sigma_B^6(1-\lambda)^p+r^6\right)^\frac{1}{6} + ```` + with $V^A$ and $V^B$ being the regular Lennard-Jones potentials in states $A$ and $B$ respectively. + + The force on each atom is + ```math + F(r)=(1-\lambda)F^A(r_A)\left(\frac{r}{r_A}\right)^5 + \lambda F^B(r_B)\left(\frac{r}{r_B}\right)^5 + ``` + with $F^A$ and $F^B$ being the regular Lennard-Jones forces in states $A$ and $B$ respectively. """ -@kwdef struct LennardJonesSoftCoreBeutler{C, A, L, H, S, E, W, R} <: PairwiseInteraction - cutoff::C = NoCutoff() - α::A = 1 - λ::L = 0 - use_neighbors::Bool = false - shortcut::H = lj_zero_shortcut - σ_mixing::S = lorentz_σ_mixing - ϵ_mixing::E = geometric_ϵ_mixing - weight_special::W = 1 - σ6_fac::R = (α * (1-λ)) +struct LennardJonesSoftCoreBeutler{L, A, P, C, F, W, SA, SB} <: + ThermodynamicIntegrationInteraction + λ::L + α::A + p::P + cutoff::C + use_neighbors::Bool + inter_state_a::SA + inter_state_b::SB + shortcut::F + weight_special::W +end + +function LennardJonesSoftCoreBeutler(; + λ, + α, + p, + cutoff = NoCutoff(), + use_neighbors = false, + inter_state_a = nothing, + inter_state_b = (σ_mixing = lorentz_σ_mixing, ϵ_mixing = geometric_ϵ_mixing), + shortcut = lj_zero_shortcut, + weight_special = 1) + return LennardJonesSoftCoreBeutler{typeof(λ), typeof(α), typeof(p), typeof(cutoff), + typeof(shortcut), typeof(weight_special), typeof(inter_state_a), typeof(inter_state_b)}( + λ, α, p, cutoff, use_neighbors, inter_state_a, + inter_state_b, shortcut, weight_special) end use_neighbors(inter::LennardJonesSoftCoreBeutler) = inter.use_neighbors -function Base.zero(lj::LennardJonesSoftCoreBeutler{C, A, L, H, S, E, W, R}) where {C, A, L, H, S, E, W, R} - return LennardJonesSoftCoreBeutler( - lj.cutoff, - zero(A), - zero(L), - lj.use_neighbors, - lj.shortcut, - lj.σ_mixing, - lj.ϵ_mixing, - zero(W), - zero(R), +function Base.zero(inter::LennardJonesSoftCoreBeutler{L, A, P, C, F, W, SA, SB}) where {L, A, P, C, F, W, SA, SB} + return LennardJonesSoftCoreBeutler{L, A, P, C, F, W, SA, SB}( + zero(inter.λ), + zero(inter.α), + zero(inter.p), + inter.cutoff, + inter.use_neighbors, + inter.inter_state_a, + inter.inter_state_b, + inter.shortcut, + zero(inter.weight_special) ) end function Base.:+(l1::LennardJonesSoftCoreBeutler, l2::LennardJonesSoftCoreBeutler) return LennardJonesSoftCoreBeutler( - l1.cutoff, - l1.α + l2.α, l1.λ + l2.λ, + l1.α + l2.α, + l1.p + l2.p, + l1.cutoff, l1.use_neighbors, + l1.inter_state_a, + l1.inter_state_b, l1.shortcut, - l1.σ_mixing, - l1.ϵ_mixing, - l1.weight_special + l2.weight_special, - l1.σ6_fac + l2.σ6_fac, + l1.weight_special + l2.weight_special ) end -@inline function force(inter::LennardJonesSoftCoreBeutler, - dr, - atom_i, - atom_j, - force_units=u"kJ * mol^-1 * nm^-1", - special=false, - args...) - if inter.shortcut(atom_i, atom_j) +const lennard_jones = LennardJones() + +@inline function Molly.force( + sc::LennardJonesSoftCoreBeutler, + dr, + atom_i, + atom_j, + force_units = u"kJ * mol^-1 * nm^-1", + special = false, + args...) + zero_force = ustrip(zero(dr[1])) * force_units + if sc.shortcut(atom_i, atom_j) return ustrip.(zero(dr)) * force_units end - σ6 = inter.σ_mixing(atom_i, atom_j)^6 - ϵ = inter.ϵ_mixing(atom_i, atom_j) - cutoff = inter.cutoff r = norm(dr) - C6 = 4 * ϵ * σ6 - params = (C6 * σ6, C6, inter.σ6_fac) + params = (sc.λ, sc.α, sc.p, sc.inter_state_a, sc.inter_state_b, atom_i, atom_j, zero_force) - f = force_cutoff(cutoff, inter, r, params) + f = force_cutoff(sc.cutoff, sc, r, params) fdr = (f / r) * dr if special - return fdr * inter.weight_special + return fdr * sc.weight_special else return fdr end end -function pairwise_force(::LennardJonesSoftCoreBeutler, r, (C12, C6, σ6_fac)) - R = sqrt(cbrt((σ6_fac*(C12/C6))+r^6)) - R6 = R^6 - return (((12*C12)/(R6*R6*R)) - ((6*C6)/(R6*R)))*((r/R)^5) +@inline function get_ra(r6, σ2, α, λ, p) + return cbrt(sqrt(α * σ2^3 * λ^p + r6)) end -@inline function potential_energy(inter::LennardJonesSoftCoreBeutler, - dr, - atom_i, - atom_j, - energy_units=u"kJ * mol^-1", - special=false, - args...) - if inter.shortcut(atom_i, atom_j) - return ustrip(zero(dr[1])) * energy_units +@inline function get_rb(r6, σ2, α, λ, p) + return cbrt(sqrt(α * σ2^3 * (1 - λ)^p + r6)) +end + +@inline function pairwise_force(::LennardJonesSoftCoreBeutler, r, + (λ, α, p, inter_state_a, inter_state_b, atom_i, atom_j, zero_force)) + r6 = r^6 + + force_term_a, force_term_b = zero_force, zero_force + + # force for system in state A + if !isnothing(inter_state_a) + σ_a = inter_state_a.σ_mixing(atom_i, atom_j) + ϵ_a = inter_state_a.ϵ_mixing(atom_i, atom_j) + + σ_a2 = σ_a^2 + r_a = get_ra(r6, σ_a2, α, λ, p) + + force_term_a = (1 - λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) * r6 / r_a^5 + end + + # force for system in state B + if !isnothing(inter_state_b) + σ_b = inter_state_b.σ_mixing(atom_i, atom_j) + ϵ_b = inter_state_b.ϵ_mixing(atom_i, atom_j) + + σ_b2 = σ_b^2 + r_b = get_rb(r6, σ_b2, α, λ, p) + + force_term_b = λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) * r6 / r_b^5 + end + + return force_term_a + force_term_b +end + +@inline function Molly.potential_energy( + sc::LennardJonesSoftCoreBeutler, + dr, + atom_i, + atom_j, + energy_units = u"kJ*mol^-1", + special = false, + args...) + zero_energy = ustrip(zero(dr[1])) * energy_units + if sc.shortcut(atom_i, atom_j) + return zero_energy end - σ6 = inter.σ_mixing(atom_i, atom_j)^6 - ϵ = inter.ϵ_mixing(atom_i, atom_j) - cutoff = inter.cutoff r = norm(dr) - C6 = 4 * ϵ * σ6 - params = (C6 * σ6, C6, inter.σ6_fac) + params = ( + sc.λ, sc.α, sc.p, sc.inter_state_a, sc.inter_state_b, atom_i, atom_j, zero_energy) - pe = pe_cutoff(cutoff, inter, r, params) + pe = pe_cutoff(sc.cutoff, sc, r, params) if special - return pe * inter.weight_special + return pe * sc.weight_special else return pe end end -function pairwise_pe(::LennardJonesSoftCoreBeutler, r, (C12, C6, σ6_fac)) - R6 = (σ6_fac*(C12/C6))+r^6 - return ((C12/(R6*R6)) - (C6/(R6))) +@inline function pairwise_pe( + ::LennardJonesSoftCoreBeutler, r, ( + λ, α, p, inter_state_a, inter_state_b, atom_i, atom_j, zero_energy)) + r6 = r^6 + + energy_term_a, energy_term_b = zero_energy, zero_energy + + # energy for system in state A + if !isnothing(inter_state_a) + σ_a = inter_state_a.σ_mixing(atom_i, atom_j) + ϵ_a = inter_state_a.ϵ_mixing(atom_i, atom_j) + + σ_a2 = σ_a^2 + r_a = cbrt(sqrt(α * σ_a2^3 * λ^p + r6)) + + energy_term_a = (1 - λ) * pairwise_pe(lennard_jones, r_a, (σ_a2, ϵ_a)) + end + + # energy for system in state B + if !isnothing(inter_state_b) + σ_b = inter_state_b.σ_mixing(atom_i, atom_j) + ϵ_b = inter_state_b.ϵ_mixing(atom_i, atom_j) + + σ_b2 = σ_b^2 + r_b = cbrt(sqrt(α * σ_b2^3 * (1 - λ)^p + r6)) + + energy_term_b = λ * pairwise_pe(lennard_jones, r_b, (σ_b2, ϵ_b)) + end + + return energy_term_a + energy_term_b end +@inline function ∂H_∂λ( + sc::LennardJonesSoftCoreBeutler, + dr, + atom_i, + atom_j, + energy_units = u"kJ*mol^-1", + special = false) + zero_energy = ustrip(zero(dr[1])) * energy_units + if sc.shortcut(atom_i, atom_j) + return zero_energy + end + + r = norm(dr) + r6 = r^6 + pα_6 = sc.p * sc.α / 6 + + term_state_a, term_state_b = zero_energy, zero_energy + + # ∂V/∂λ for system in state A + if !isnothing(sc.inter_state_a) + σ_a = sc.inter_state_a.σ_mixing(atom_i, atom_j) + ϵ_a = sc.inter_state_a.ϵ_mixing(atom_i, atom_j) + + σ_a2 = σ_a^2 + r_a = cbrt(sqrt(sc.α * σ_a2^3 * sc.λ^p + r6)) + + term_state_a = pairwise_pe(lennard_jones, r_a, (σ_a2, ϵ_a)) + + (1 - sc.λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) / + r_a^5 * σ_a2^3 * sc.λ^(p - 1) + end + + # ∂V/∂λ for system in state B + if !isnothing(sc.inter_state_b) + σ_b = sc.inter_state_b.σ_mixing(atom_i, atom_j) + ϵ_b = sc.inter_state_b.ϵ_mixing(atom_i, atom_j) + + σ_b2 = σ_b^2 + r_b = cbrt(sqrt(sc.α * σ_b2^3 * (1 - sc.λ)^sc.p + r6)) + + term_state_b = pairwise_pe(lennard_jones, r_b, (σ_b2, ϵ_b)) + + pα_6 * sc.λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) / + r_b^5 * σ_b2^3 * (1 - sc.λ)^(sc.p - 1) + end + + if special + return (term_state_b - term_state_a) * sc.weight_special + else + return term_state_b - term_state_a + end +end + + @doc raw""" LennardJonesSoftCoreGapsys(; cutoff, α, λ, use_neighbors, shortcut, σ_mixing, ϵ_mixing, weight_special) From 13540661ade5890035661309485916988ba6c85a Mon Sep 17 00:00:00 2001 From: Immanuel Felix Sulzer <100942429+Immi000@users.noreply.github.com> Date: Mon, 3 Nov 2025 16:25:22 +0100 Subject: [PATCH 2/8] pairwise interaction --- src/interactions/lennard_jones.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 9ca368e2..9b1cf5cf 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -162,8 +162,7 @@ end ``` with $F^A$ and $F^B$ being the regular Lennard-Jones forces in states $A$ and $B$ respectively. """ -struct LennardJonesSoftCoreBeutler{L, A, P, C, F, W, SA, SB} <: - ThermodynamicIntegrationInteraction +struct LennardJonesSoftCoreBeutler{L, A, P, C, F, W, SA, SB} <: PairwiseInteraction λ::L α::A p::P From 5118d70074f0918f7a48838d8460b61aca2ff7b9 Mon Sep 17 00:00:00 2001 From: Immanuel Sulzer Date: Mon, 3 Nov 2025 17:17:49 +0100 Subject: [PATCH 3/8] default value for p and remove wrong NoUnits --- docs/src/examples.md | 2 +- src/interactions/lennard_jones.jl | 9 ++++----- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/docs/src/examples.md b/docs/src/examples.md index aff06a14..b479d2d6 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -988,7 +988,7 @@ energies(inter) = map(dists) do dist c1 = SVector(1.0, 1.0, 1.0)u"nm" c2 = SVector(dist + 1.0, 1.0, 1.0)u"nm" vec = vector(c1, c2, boundary) - potential_energy(inter, vec, a1, a2, NoUnits) + potential_energy(inter, vec, a1, a2) end function plot_interactions(ax, title, xlabel, ylabel, data, ylims_range) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 9b1cf5cf..725f8417 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -177,7 +177,7 @@ end function LennardJonesSoftCoreBeutler(; λ, α, - p, + p = 1, cutoff = NoCutoff(), use_neighbors = false, inter_state_a = nothing, @@ -291,7 +291,7 @@ end dr, atom_i, atom_j, - energy_units = u"kJ*mol^-1", + energy_units = u"kJ * mol^-1", special = false, args...) zero_energy = ustrip(zero(dr[1])) * energy_units @@ -300,8 +300,7 @@ end end r = norm(dr) - params = ( - sc.λ, sc.α, sc.p, sc.inter_state_a, sc.inter_state_b, atom_i, atom_j, zero_energy) + params = (sc.λ, sc.α, sc.p, sc.inter_state_a, sc.inter_state_b, atom_i, atom_j, zero_energy) pe = pe_cutoff(sc.cutoff, sc, r, params) if special @@ -348,7 +347,7 @@ end dr, atom_i, atom_j, - energy_units = u"kJ*mol^-1", + energy_units = u"kJ * mol^-1", special = false) zero_energy = ustrip(zero(dr[1])) * energy_units if sc.shortcut(atom_i, atom_j) From dc54df4fc9e3bcc577e19d6b9f30e0303b5db1bb Mon Sep 17 00:00:00 2001 From: Immanuel Date: Mon, 10 Nov 2025 14:19:25 +0100 Subject: [PATCH 4/8] prettification and final bugfixes --- src/interactions/lennard_jones.jl | 52 +++++++++++-------------------- test/interactions.jl | 8 ++--- 2 files changed, 23 insertions(+), 37 deletions(-) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 725f8417..62632a5e 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -162,32 +162,16 @@ end ``` with $F^A$ and $F^B$ being the regular Lennard-Jones forces in states $A$ and $B$ respectively. """ -struct LennardJonesSoftCoreBeutler{L, A, P, C, F, W, SA, SB} <: PairwiseInteraction - λ::L - α::A - p::P - cutoff::C - use_neighbors::Bool - inter_state_a::SA - inter_state_b::SB - shortcut::F - weight_special::W -end - -function LennardJonesSoftCoreBeutler(; - λ, - α, - p = 1, - cutoff = NoCutoff(), - use_neighbors = false, - inter_state_a = nothing, - inter_state_b = (σ_mixing = lorentz_σ_mixing, ϵ_mixing = geometric_ϵ_mixing), - shortcut = lj_zero_shortcut, - weight_special = 1) - return LennardJonesSoftCoreBeutler{typeof(λ), typeof(α), typeof(p), typeof(cutoff), - typeof(shortcut), typeof(weight_special), typeof(inter_state_a), typeof(inter_state_b)}( - λ, α, p, cutoff, use_neighbors, inter_state_a, - inter_state_b, shortcut, weight_special) +@kwdef struct LennardJonesSoftCoreBeutler{L, A, P, C, F, W, SA, SB} <: PairwiseInteraction + λ::L = 1 + α::A = 1 + p::P = 1 + cutoff::C = NoCutoff() + use_neighbors::Bool = false + inter_state_a::SA = nothing + inter_state_b::SB = (σ_mixing = lorentz_σ_mixing, ϵ_mixing = geometric_ϵ_mixing) + shortcut::F = lj_zero_shortcut + weight_special::W = 1 end use_neighbors(inter::LennardJonesSoftCoreBeutler) = inter.use_neighbors @@ -230,11 +214,13 @@ const lennard_jones = LennardJones() force_units = u"kJ * mol^-1 * nm^-1", special = false, args...) - zero_force = ustrip(zero(dr[1])) * force_units + if sc.shortcut(atom_i, atom_j) return ustrip.(zero(dr)) * force_units end + zero_force = ustrip(zero(dr[1])) * force_units + r = norm(dr) params = (sc.λ, sc.α, sc.p, sc.inter_state_a, sc.inter_state_b, atom_i, atom_j, zero_force) @@ -269,7 +255,7 @@ end σ_a2 = σ_a^2 r_a = get_ra(r6, σ_a2, α, λ, p) - force_term_a = (1 - λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) * r6 / r_a^5 + force_term_a = (1 - λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) * (r / r_a)^5 end # force for system in state B @@ -280,7 +266,7 @@ end σ_b2 = σ_b^2 r_b = get_rb(r6, σ_b2, α, λ, p) - force_term_b = λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) * r6 / r_b^5 + force_term_b = λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) * (r / r_b)^5 end return force_term_a + force_term_b @@ -323,7 +309,7 @@ end ϵ_a = inter_state_a.ϵ_mixing(atom_i, atom_j) σ_a2 = σ_a^2 - r_a = cbrt(sqrt(α * σ_a2^3 * λ^p + r6)) + r_a = get_ra(r6, σ_a2, α, λ, p) energy_term_a = (1 - λ) * pairwise_pe(lennard_jones, r_a, (σ_a2, ϵ_a)) end @@ -334,7 +320,7 @@ end ϵ_b = inter_state_b.ϵ_mixing(atom_i, atom_j) σ_b2 = σ_b^2 - r_b = cbrt(sqrt(α * σ_b2^3 * (1 - λ)^p + r6)) + r_b = get_rb(r6, σ_b2, α, λ, p) energy_term_b = λ * pairwise_pe(lennard_jones, r_b, (σ_b2, ϵ_b)) end @@ -366,7 +352,7 @@ end ϵ_a = sc.inter_state_a.ϵ_mixing(atom_i, atom_j) σ_a2 = σ_a^2 - r_a = cbrt(sqrt(sc.α * σ_a2^3 * sc.λ^p + r6)) + r_a = get_ra(r6, σ_a2, α, λ, p) term_state_a = pairwise_pe(lennard_jones, r_a, (σ_a2, ϵ_a)) + (1 - sc.λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) / @@ -379,7 +365,7 @@ end ϵ_b = sc.inter_state_b.ϵ_mixing(atom_i, atom_j) σ_b2 = σ_b^2 - r_b = cbrt(sqrt(sc.α * σ_b2^3 * (1 - sc.λ)^sc.p + r6)) + r_b = get_rb(r6, σ_b2, α, λ, p) term_state_b = pairwise_pe(lennard_jones, r_b, (σ_b2, ϵ_b)) + pα_6 * sc.λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) / diff --git a/test/interactions.jl b/test/interactions.jl index 67a05696..79669ef9 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -47,22 +47,22 @@ inter = LennardJonesSoftCoreBeutler(α=0.3, λ=0.5) @test isapprox( force(inter, dr14, a1, a1), - SVector(35.093676538737824, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + SVector(17.546838269368916, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; atol=1e-9u"kJ * mol^-1 * nm^-1", ) @test isapprox( force(inter, dr13, a1, a1), - SVector(-1.3236594727846438, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + SVector(-0.6618297363923222, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; atol=1e-9u"kJ * mol^-1 * nm^-1", ) @test isapprox( potential_energy(inter, dr14, a1, a1), - 29.629058917654785u"kJ * mol^-1"; + 14.814529458827394u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) @test isapprox( potential_energy(inter, dr13, a1, a1), - -0.11464014233084913u"kJ * mol^-1"; + -0.05732007116542457u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) From 8e22d435e6858aeee607be394ad6bf62cc1e1c30 Mon Sep 17 00:00:00 2001 From: Immanuel Date: Mon, 10 Nov 2025 14:53:29 +0100 Subject: [PATCH 5/8] formatting and clarification --- src/interactions/lennard_jones.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 62632a5e..c64b7900 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -144,9 +144,10 @@ end The Lennard-Jones soft core interaction as described by [Beutler et al. 1994](https://doi.org/10.1016/0009-2614(94)00397-1). This interaction is a linear - interpolation between two Lennard-Jones interactions A and B controlled by a parameter λ. The - parameters `inter_state_a` and `inter_state_b` can be a named tuple of `σ_mixing` and `ϵ_mixing` - for the interaction in the respective state or `nothing` for no interaction in that state. + interpolation between two soft-core Lennard-Jones interactions A and B controlled by a parameter λ. + The system is in state A for λ = 0 and in state B for λ = 1. The parameters `inter_state_a` and + `inter_state_b` can be a named tuple of `σ_mixing` and `ϵ_mixing` for the interaction in the + respective state or `nothing` for no interaction in that state. The potential energy is defined as ````math @@ -243,6 +244,7 @@ end @inline function pairwise_force(::LennardJonesSoftCoreBeutler, r, (λ, α, p, inter_state_a, inter_state_b, atom_i, atom_j, zero_force)) + r6 = r^6 force_term_a, force_term_b = zero_force, zero_force @@ -280,6 +282,7 @@ end energy_units = u"kJ * mol^-1", special = false, args...) + zero_energy = ustrip(zero(dr[1])) * energy_units if sc.shortcut(atom_i, atom_j) return zero_energy @@ -296,9 +299,8 @@ end end end -@inline function pairwise_pe( - ::LennardJonesSoftCoreBeutler, r, ( - λ, α, p, inter_state_a, inter_state_b, atom_i, atom_j, zero_energy)) +@inline function pairwise_pe(::LennardJonesSoftCoreBeutler, r, + (λ, α, p, inter_state_a, inter_state_b, atom_i, atom_j, zero_energy)) r6 = r^6 energy_term_a, energy_term_b = zero_energy, zero_energy @@ -335,6 +337,7 @@ end atom_j, energy_units = u"kJ * mol^-1", special = false) + zero_energy = ustrip(zero(dr[1])) * energy_units if sc.shortcut(atom_i, atom_j) return zero_energy From c4f18ac5aa185b585d2526f62411b6bf0a47b2fb Mon Sep 17 00:00:00 2001 From: Immanuel Date: Mon, 10 Nov 2025 17:11:27 +0100 Subject: [PATCH 6/8] missing factor --- src/interactions/lennard_jones.jl | 55 ++++++++++++++++++++++++++++++- 1 file changed, 54 insertions(+), 1 deletion(-) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index c64b7900..b2a9f7d4 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -358,7 +358,7 @@ end r_a = get_ra(r6, σ_a2, α, λ, p) term_state_a = pairwise_pe(lennard_jones, r_a, (σ_a2, ϵ_a)) + - (1 - sc.λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) / + pα_6 * (1 - sc.λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) / r_a^5 * σ_a2^3 * sc.λ^(p - 1) end @@ -541,6 +541,59 @@ function pairwise_pe(inter::LennardJonesSoftCoreGapsys, r, (C12, C6)) end end +@inline function ∂H_∂λ( + sc::LennardJonesSoftCoreGapsys, + dr, + atom_i, + atom_j, + energy_units = u"kJ * mol^-1", + special = false) + + zero_energy = ustrip(zero(dr[1])) * energy_units + if sc.shortcut(atom_i, atom_j) + return zero_energy + end + + r = norm(dr) + r6 = r^6 + pα_6 = sc.p * sc.α / 6 + + term_state_a, term_state_b = zero_energy, zero_energy + + # ∂V/∂λ for system in state A + if !isnothing(sc.inter_state_a) + σ_a = sc.inter_state_a.σ_mixing(atom_i, atom_j) + ϵ_a = sc.inter_state_a.ϵ_mixing(atom_i, atom_j) + + σ_a2 = σ_a^2 + r_a = get_ra(r6, σ_a2, α, λ, p) + + term_state_a = pairwise_pe(lennard_jones, r_a, (σ_a2, ϵ_a)) + + (1 - sc.λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) / + r_a^5 * σ_a2^3 * sc.λ^(p - 1) + end + + # ∂V/∂λ for system in state B + if !isnothing(sc.inter_state_b) + σ_b = sc.inter_state_b.σ_mixing(atom_i, atom_j) + ϵ_b = sc.inter_state_b.ϵ_mixing(atom_i, atom_j) + + σ_b2 = σ_b^2 + r_b = get_rb(r6, σ_b2, α, λ, p) + + term_state_b = pairwise_pe(lennard_jones, r_b, (σ_b2, ϵ_b)) + + pα_6 * sc.λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) / + r_b^5 * σ_b2^3 * (1 - sc.λ)^(sc.p - 1) + end + + if special + return (term_state_b - term_state_a) * sc.weight_special + else + return term_state_b - term_state_a + end + +end + @doc raw""" AshbaughHatch(; cutoff, use_neighbors, shortcut, ϵ_mixing, σ_mixing, λ_mixing, weight_special) From 5968cee7191054eef2010433c29d136c7b44a374 Mon Sep 17 00:00:00 2001 From: Immanuel Date: Mon, 10 Nov 2025 18:04:16 +0100 Subject: [PATCH 7/8] =?UTF-8?q?pairwise=5F=E2=88=82H=5F=E2=88=82=CE=BB=20f?= =?UTF-8?q?unction=20for=20future=20cutoff=20support=20and=20=E2=88=82H=5F?= =?UTF-8?q?=E2=88=82=CE=BB-tests?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/interactions/lennard_jones.jl | 44 +++++++++++++++++++------------ test/interactions.jl | 35 ++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 17 deletions(-) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index b2a9f7d4..4a5ad3b7 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -344,45 +344,55 @@ end end r = norm(dr) + params = (sc.λ, sc.α, sc.p, sc.inter_state_a, sc.inter_state_b, atom_i, atom_j, zero_energy) + + dh_dl = pairwise_∂H_∂λ(sc, r, params) + + if special + return dh_dl * sc.weight_special + else + return dh_dl + end +end + +@inline function pairwise_∂H_∂λ(::LennardJonesSoftCoreBeutler, r, + (λ, α, p, inter_state_a, inter_state_b, atom_i, atom_j, zero_energy)) + r6 = r^6 - pα_6 = sc.p * sc.α / 6 + pα_6 = p * α / 6 term_state_a, term_state_b = zero_energy, zero_energy # ∂V/∂λ for system in state A - if !isnothing(sc.inter_state_a) - σ_a = sc.inter_state_a.σ_mixing(atom_i, atom_j) - ϵ_a = sc.inter_state_a.ϵ_mixing(atom_i, atom_j) + if !isnothing(inter_state_a) + σ_a = inter_state_a.σ_mixing(atom_i, atom_j) + ϵ_a = inter_state_a.ϵ_mixing(atom_i, atom_j) σ_a2 = σ_a^2 r_a = get_ra(r6, σ_a2, α, λ, p) term_state_a = pairwise_pe(lennard_jones, r_a, (σ_a2, ϵ_a)) + - pα_6 * (1 - sc.λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) / - r_a^5 * σ_a2^3 * sc.λ^(p - 1) + pα_6 * (1 - λ) * pairwise_force(lennard_jones, r_a, (σ_a2, ϵ_a)) / + r_a^5 * σ_a2^3 * λ^(p - 1) end # ∂V/∂λ for system in state B - if !isnothing(sc.inter_state_b) - σ_b = sc.inter_state_b.σ_mixing(atom_i, atom_j) - ϵ_b = sc.inter_state_b.ϵ_mixing(atom_i, atom_j) + if !isnothing(inter_state_b) + σ_b = inter_state_b.σ_mixing(atom_i, atom_j) + ϵ_b = inter_state_b.ϵ_mixing(atom_i, atom_j) σ_b2 = σ_b^2 r_b = get_rb(r6, σ_b2, α, λ, p) term_state_b = pairwise_pe(lennard_jones, r_b, (σ_b2, ϵ_b)) + - pα_6 * sc.λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) / - r_b^5 * σ_b2^3 * (1 - sc.λ)^(sc.p - 1) + pα_6 * λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) / + r_b^5 * σ_b2^3 * (1 - λ)^(p - 1) end - if special - return (term_state_b - term_state_a) * sc.weight_special - else - return term_state_b - term_state_a - end + return term_state_a + term_state_b + end - @doc raw""" LennardJonesSoftCoreGapsys(; cutoff, α, λ, use_neighbors, shortcut, σ_mixing, ϵ_mixing, weight_special) diff --git a/test/interactions.jl b/test/interactions.jl index 79669ef9..f995b6b8 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -66,6 +66,41 @@ atol=1e-9u"kJ * mol^-1", ) + @test isapprox( + Molly.∂H_∂λ(inter, dr12, a1, a2), + -0.12598966767217717u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + @test isapprox( + Molly.∂H_∂λ(inter, dr13, a1, a2), + -0.0317055834215144u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + @test isapprox( + Molly.∂H_∂λ(inter, dr14, a1, a2), + 62.957442851334115u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + + inter = LennardJonesSoftCoreBeutler(α=0.5, λ=0.8, + inter_state_a=(σ_mixing=Molly.lorentz_σ_mixing, ϵ_mixing=Molly.geometric_ϵ_mixing), + inter_state_b=nothing) + @test isapprox( + Molly.∂H_∂λ(inter, dr12, a1, a2), + -0.1197457491693602u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + @test isapprox( + Molly.∂H_∂λ(inter, dr13, a1, a2), + -0.031184500618766355u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + @test isapprox( + Molly.∂H_∂λ(inter, dr14, a1, a2), + 3.4324768629314484u"kJ * mol^-1"; + atol=1e-9u"kJ * mol^-1", + ) + inter = LennardJonesSoftCoreGapsys(α=0.85, λ=0.5) @test isapprox( force(inter, dr14, a1, a1), From 1e3e86ccadf4878837787e0ef0f4cf1a287b40b6 Mon Sep 17 00:00:00 2001 From: Immanuel Sulzer Date: Wed, 26 Nov 2025 12:41:14 +0100 Subject: [PATCH 8/8] explicit non-interaction instead of nothing --- src/interactions/lennard_jones.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 4a5ad3b7..5e406e41 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -169,12 +169,14 @@ end p::P = 1 cutoff::C = NoCutoff() use_neighbors::Bool = false - inter_state_a::SA = nothing + inter_state_a::SA = NoInterState() inter_state_b::SB = (σ_mixing = lorentz_σ_mixing, ϵ_mixing = geometric_ϵ_mixing) shortcut::F = lj_zero_shortcut weight_special::W = 1 end +struct NoInterState end + use_neighbors(inter::LennardJonesSoftCoreBeutler) = inter.use_neighbors function Base.zero(inter::LennardJonesSoftCoreBeutler{L, A, P, C, F, W, SA, SB}) where {L, A, P, C, F, W, SA, SB} @@ -250,7 +252,7 @@ end force_term_a, force_term_b = zero_force, zero_force # force for system in state A - if !isnothing(inter_state_a) + if !isa(inter_state_a, NoInterState) σ_a = inter_state_a.σ_mixing(atom_i, atom_j) ϵ_a = inter_state_a.ϵ_mixing(atom_i, atom_j) @@ -261,8 +263,8 @@ end end # force for system in state B - if !isnothing(inter_state_b) - σ_b = inter_state_b.σ_mixing(atom_i, atom_j) + if !isa(inter_state_b, NoInterState) + σ_b = inter_state_b.σ_mixing(atom_i, atom_j) ϵ_b = inter_state_b.ϵ_mixing(atom_i, atom_j) σ_b2 = σ_b^2 @@ -306,7 +308,7 @@ end energy_term_a, energy_term_b = zero_energy, zero_energy # energy for system in state A - if !isnothing(inter_state_a) + if !isa(inter_state_a, NoInterState) σ_a = inter_state_a.σ_mixing(atom_i, atom_j) ϵ_a = inter_state_a.ϵ_mixing(atom_i, atom_j) @@ -317,7 +319,7 @@ end end # energy for system in state B - if !isnothing(inter_state_b) + if !isa(inter_state_b, NoInterState) σ_b = inter_state_b.σ_mixing(atom_i, atom_j) ϵ_b = inter_state_b.ϵ_mixing(atom_i, atom_j) @@ -364,7 +366,7 @@ end term_state_a, term_state_b = zero_energy, zero_energy # ∂V/∂λ for system in state A - if !isnothing(inter_state_a) + if !isa(inter_state_a, NoInterState) σ_a = inter_state_a.σ_mixing(atom_i, atom_j) ϵ_a = inter_state_a.ϵ_mixing(atom_i, atom_j) @@ -377,7 +379,7 @@ end end # ∂V/∂λ for system in state B - if !isnothing(inter_state_b) + if !isa(inter_state_b, NoInterState) σ_b = inter_state_b.σ_mixing(atom_i, atom_j) ϵ_b = inter_state_b.ϵ_mixing(atom_i, atom_j)