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 87f7630a..5e406e41 100644
--- a/src/interactions/lennard_jones.jl
+++ b/src/interactions/lennard_jones.jl
@@ -140,140 +140,259 @@ 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 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
+    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()
+@kwdef struct LennardJonesSoftCoreBeutler{L, A, P, C, F, W, SA, SB} <: PairwiseInteraction
+    λ::L = 1
     α::A = 1
-    λ::L = 0
+    p::P = 1
+    cutoff::C = NoCutoff()
     use_neighbors::Bool = false
-    shortcut::H = lj_zero_shortcut
-    σ_mixing::S = lorentz_σ_mixing
-    ϵ_mixing::E = geometric_ϵ_mixing
+    inter_state_a::SA = NoInterState()
+    inter_state_b::SB = (σ_mixing = lorentz_σ_mixing, ϵ_mixing = geometric_ϵ_mixing)
+    shortcut::F = lj_zero_shortcut
     weight_special::W = 1
-    σ6_fac::R = (α * (1-λ))
 end
 
+struct NoInterState 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...)
+    
+    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
+    zero_force = ustrip(zero(dr[1])) * force_units
+
     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 !isa(inter_state_a, NoInterState)
+        σ_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)) * (r / r_a)^5
+    end
+
+    # force for system in 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)
+
+        σ_b2 = σ_b^2
+        r_b = get_rb(r6, σ_b2, α, λ, p)
+
+        force_term_b = λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) * (r / 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 !isa(inter_state_a, NoInterState)
+        σ_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)
+
+        energy_term_a = (1 - λ) * pairwise_pe(lennard_jones, r_a, (σ_a2, ϵ_a))
+    end
+
+    # energy for system in 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)
+
+        σ_b2 = σ_b^2
+        r_b = get_rb(r6, σ_b2, α, λ, p)
+
+        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)
+    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 = p * α / 6
+
+    term_state_a, term_state_b = zero_energy, zero_energy
+
+    # ∂V/∂λ for system in 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)
+
+        σ_a2 = σ_a^2
+        r_a = get_ra(r6, σ_a2, α, λ, p)
+
+        term_state_a = pairwise_pe(lennard_jones, r_a, (σ_a2, ϵ_a)) +
+                       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 !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
+        r_b = get_rb(r6, σ_b2, α, λ, p)
+
+        term_state_b = pairwise_pe(lennard_jones, r_b, (σ_b2, ϵ_b)) +
+                       pα_6 * λ * pairwise_force(lennard_jones, r_b, (σ_b2, ϵ_b)) /
+                       r_b^5 * σ_b2^3 * (1 - λ)^(p - 1)
+    end
+
+    return term_state_a + term_state_b
+    
 end
 
 @doc raw"""
@@ -434,6 +553,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)
diff --git a/test/interactions.jl b/test/interactions.jl
index 67a05696..f995b6b8 100644
--- a/test/interactions.jl
+++ b/test/interactions.jl
@@ -47,22 +47,57 @@
     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",
+    )
+
+    @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",
     )