diff --git a/src/interactions/coulomb.jl b/src/interactions/coulomb.jl
index 6c040864..c6b1466a 100644
--- a/src/interactions/coulomb.jl
+++ b/src/interactions/coulomb.jl
@@ -342,10 +342,6 @@ end
     end
 end
 
-
-
-
-
 @doc raw"""
     Yukawa(; cutoff, use_neighbors, weight_special, coulomb_const, kappa)
 
@@ -355,30 +351,28 @@ The potential energy is defined as
 ```math
 V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}} \exp(-\kappa r_{ij})
 ```
-
 and the force on each atom by 
-
 ```math
 F(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}^2} \exp(-\kappa r_{ij})\left(\kappa r_{ij} + 1\right) \vec{r}_{ij}
 ```
 """
-@kwdef struct Yukawa{C, W, T, D} 
-    cutoff::C= NoCutoff()
-    use_neighbors::Bool= false
-    weight_special::W =1
+@kwdef struct Yukawa{C, W, T, K} 
+    cutoff::C = NoCutoff()
+    use_neighbors::Bool = false
+    weight_special::W = 1
     coulomb_const::T = coulomb_const
-    kappa::D = 1.0*u"nm^-1"
+    kappa::K = 1.0*u"nm^-1"
 end
 
 use_neighbors(inter::Yukawa) = inter.use_neighbors
 
-function Base.zero(yukawa::Yukawa{C, W, T ,D}) where {C, W, T, D}
-    return Yukawa{C, W, T, D}(
+function Base.zero(yukawa::Yukawa{C, W, T, K}) where {C, W, T, K}
+    return Yukawa(
         yukawa.cutoff,
         yukawa.use_neighbors,
-        yukawa.weight_special,
-        yukawa.coulomb_const,
-        yukawa.kappa
+        zero(W),
+        zero(T),
+        zero(K),
     )
 end
 
@@ -387,18 +381,18 @@ function Base.:+(c1::Yukawa, c2::Yukawa)
         c1.cutoff,
         c1.use_neighbors,
         c1.weight_special + c2.weight_special,
-        c1.coulomb_const,
-        c1.kappa
+        c1.coulomb_const + c2.coulomb_const,
+        c1.kappa + c2.kappa,
     )
 end
 
-@inline function force(inter::Yukawa{C},
-                                    dr,
-                                    atom_i,
-                                    atom_j,
-                                    force_units=u"kJ * mol^-1 * nm^-1",
-                                    special=false,
-                                    args...) where C
+@inline function force(inter::Yukawa,
+                       dr,
+                       atom_i,
+                       atom_j,
+                       force_units=u"kJ * mol^-1 * nm^-1",
+                       special=false,
+                       args...)
     r2 = sum(abs2, dr)
     cutoff = inter.cutoff
     coulomb_const = inter.coulomb_const
@@ -416,15 +410,16 @@ end
 
 function force_divr(::Yukawa, r2, invr2, (coulomb_const, qi, qj, kappa))
     r = sqrt(r2)
-    return (coulomb_const * qi * qj) / r^3 *exp(-kappa*r) * (kappa*r+1)
+    return (coulomb_const * qi * qj) * exp(-kappa * r) * (kappa * r + 1) / r^3
 end
 
-@inline function potential_energy(inter::Yukawa{C},
-                                            dr,
-                                            atom_i,
-                                            atom_j,
-                                            energy_units=u"kJ * mol^-1",
-                                            special::Bool=false, args...) where C
+@inline function potential_energy(inter::Yukawa,
+                                  dr,
+                                  atom_i,
+                                  atom_j,
+                                  energy_units=u"kJ * mol^-1",
+                                  special::Bool=false,
+                                  args...)
     r2 = sum(abs2, dr)
     cutoff = inter.cutoff
     coulomb_const = inter.coulomb_const
@@ -440,5 +435,5 @@ end
 end
 
 function potential(::Yukawa, r2, invr2, (coulomb_const, qi, qj, kappa))
-    return (coulomb_const * qi * qj) * √invr2 * exp(-kappa*sqrt(r2))
+    return (coulomb_const * qi * qj) * √invr2 * exp(-kappa * sqrt(r2))
 end
diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl
index 0aa8ef35..d3a6252e 100644
--- a/src/interactions/lennard_jones.jl
+++ b/src/interactions/lennard_jones.jl
@@ -1,8 +1,7 @@
 export
     LennardJones,
     LennardJonesSoftCore,
-    AshbaughHatch,
-    AshbaughHatchAtom
+    AshbaughHatch
 
 function lj_zero_shortcut(atom_i, atom_j)
     return iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) ||
@@ -293,15 +292,16 @@ end
 
 
 @doc raw"""
-    AshbaughHatch(; cutoff,use_neighbors,shortcut, ϵ_mixing,σ_mixing, λ_mixing,weight_special)
+    AshbaughHatch(; cutoff, use_neighbors, shortcut, ϵ_mixing, σ_mixing,
+                  λ_mixing, weight_special)
 
-The AshbaughHatch ($V_{\text{AH}}$) is a modified Lennard-Jones ($V_{\text{LJ}}$) 6-12 interaction between two atoms.
+The Ashbaugh-Hatch potential ($V_{\text{AH}}$) is a modified Lennard-Jones ($V_{\text{LJ}}$)
+6-12 interaction between two atoms.
 
-The potential energy is defined by
+The potential energy is defined as
 ```math
 V_{\text{LJ}}(r_{ij}) = 4\varepsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right] \\ 
 ```
-
 ```math
 V_{\text{AH}}(r_{ij}) =
     \begin{cases}
@@ -309,7 +309,6 @@ V_{\text{AH}}(r_{ij}) =
        λ_{ij}V_{\text{LJ}}(r_{ij})  &,  2^{1/6}σ \leq r_{ij}
     \end{cases}
 ```
-
 and the force on each atom by
 ```math
 \vec{F}_{\text{AH}} =
@@ -328,30 +327,20 @@ where
 
 If ``\lambda`` is one this gives the standard [`LennardJones`](@ref) potential.
 """
-@kwdef struct AshbaughHatch{C, H, E, F,L,W} 
+@kwdef struct AshbaughHatch{C, H, S, E, L, W} 
     cutoff::C = NoCutoff()
     use_neighbors::Bool = false
     shortcut::H = lj_zero_shortcut
+    σ_mixing::S = lorentz_σ_mixing
     ϵ_mixing::E = lorentz_ϵ_mixing
-    σ_mixing::F = lorentz_σ_mixing
     λ_mixing::L = lorentz_λ_mixing
     weight_special::W = 1
 end
 
-
-@kwdef struct AshbaughHatchAtom{C, M, S, E, L}
-    index::Int = 1
-    charge::C = 0.0
-    mass::M =1.0u"g/mol"
-    σ::S =0.0u"nm"
-    ϵ::E =0.0u"kJ * mol^-1"
-    λ::L=1.0
-end
-
 use_neighbors(inter::AshbaughHatch) = inter.use_neighbors
 
-function Base.zero(lj::AshbaughHatch{C, H, E, F,L,W}) where {C, H, E, F,L,W}
-    return AshbaughHatch{C, H, E, F,L,W}(
+function Base.zero(lj::AshbaughHatch{C, H, S, E, L, W}) where {C, H, S, E, L, W}
+    return AshbaughHatch(
         lj.cutoff,
         lj,use_neighbors,
         lj.shortcut,
@@ -362,9 +351,8 @@ function Base.zero(lj::AshbaughHatch{C, H, E, F,L,W}) where {C, H, E, F,L,W}
     )
 end
 
-function Base.:+(l1::AshbaughHatch{C, H, E, F,L,W},
-                 l2::AshbaughHatch{C, H, E, F,L,W}) where {C, H, E, F,L,W}
-    return AshbaughHatch{C, H, E, F,L,W}(
+function Base.:+(l1::AshbaughHatch, l2::AshbaughHatch)
+    return AshbaughHatch(
         l1.cutoff,
         l1.use_neighbors,
         l1.shortcut,
@@ -375,18 +363,27 @@ function Base.:+(l1::AshbaughHatch{C, H, E, F,L,W},
     )
 end
 
-@inline function force(inter::AshbaughHatch{S, C},
-                                    dr,
-                                    atom_i,
-                                    atom_j,
-                                    force_units=u"kJ * mol^-1 * nm^-1",
-                                    special::Bool=false, args...) where {S, C}
+@kwdef struct AshbaughHatchAtom{T, M, C, S, E, L}
+    index::Int = 1
+    atom_type::T = 1
+    mass::M = 1.0u"g/mol"
+    charge::C = 0.0
+    σ::S = 0.0u"nm"
+    ϵ::E = 0.0u"kJ * mol^-1"
+    λ::L = 1.0
+end
+
+@inline function force(inter::AshbaughHatch,
+                       dr,
+                       atom_i,
+                       atom_j,
+                       force_units=u"kJ * mol^-1 * nm^-1",
+                       special::Bool=false,
+                       args...)
     if inter.shortcut(atom_i, atom_j)
         return ustrip.(zero(dr)) * force_units
     end
 
-    # Lorentz-Berthelot mixing rules use the arithmetic average for σ
-    # Otherwise use the geometric average
     ϵ = inter.ϵ_mixing(atom_i, atom_j)
     σ = inter.σ_mixing(atom_i, atom_j)
     λ = inter.λ_mixing(atom_i, atom_j)
@@ -396,7 +393,7 @@ end
     σ2 = σ^2
     params = (σ2, ϵ, λ)
 
-    f = force_divr_with_cutoff(inter, r2, params, cutoff,force_units)
+    f = force_divr_with_cutoff(inter, r2, params, cutoff, force_units)
     if special
         return f * dr * inter.weight_special
     else
@@ -405,24 +402,24 @@ end
 end
 
 @inline function force_divr(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ))
-    if r2 < 2^(1/3)*σ2
+    if r2 < (2^(1/3) * σ2)
         return  force_divr(LennardJones(), r2, invr2, (σ2, ϵ))
     else
-        return λ*force_divr(LennardJones(), r2, invr2, (σ2, ϵ)) 
+        return λ * force_divr(LennardJones(), r2, invr2, (σ2, ϵ)) 
     end
 end
 
-@inline function potential_energy(inter::AshbaughHatch{S, C},
-                                            dr,
-                                            atom_i,
-                                            atom_j,
-                                            energy_units=u"kJ * mol^-1",
-                                            special::Bool=false,
-                                            args...) where {S, C}
+@inline function potential_energy(inter::AshbaughHatch,
+                                  dr,
+                                  atom_i,
+                                  atom_j,
+                                  energy_units=u"kJ * mol^-1",
+                                  special::Bool=false,
+                                  args...)
     if inter.shortcut(atom_i, atom_j)
         return ustrip(zero(dr[1])) * energy_units
     end
-    ϵ = inter.ϵ_mixing(atom_i, atom_j)     ### probably not needed for Calvados2 use case
+    ϵ = inter.ϵ_mixing(atom_i, atom_j)
     σ = inter.σ_mixing(atom_i, atom_j)
     λ = inter.λ_mixing(atom_i, atom_j)
 
@@ -431,7 +428,7 @@ end
     σ2 = σ^2
     params = (σ2, ϵ, λ)
 
-    pe = potential_with_cutoff(inter, r2, params, cutoff,energy_units)
+    pe = potential_with_cutoff(inter, r2, params, cutoff, energy_units)
     if special
         return pe * inter.weight_special
     else
@@ -440,9 +437,9 @@ end
 end
 
 @inline function potential(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ))
-    if r2 < 2^(1/3)*σ2
-        return   potential(LennardJones(), r2, invr2, (σ2, ϵ)) + ϵ*(1-λ) 
+    if r2 < (2^(1/3) * σ2)
+        return potential(LennardJones(), r2, invr2, (σ2, ϵ)) + ϵ * (1 - λ) 
     else
-        return λ* potential(LennardJones(), r2, invr2, (σ2, ϵ))
+        return λ * potential(LennardJones(), r2, invr2, (σ2, ϵ))
     end
-end
\ No newline at end of file
+end
diff --git a/test/interactions.jl b/test/interactions.jl
index 2e857518..ef4c898b 100644
--- a/test/interactions.jl
+++ b/test/interactions.jl
@@ -33,29 +33,6 @@
         )
     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),
@@ -78,6 +55,64 @@
         atol=1e-9u"kJ * mol^-1",
     )
 
+    AH_a1 = Molly.AshbaughHatchAtom(charge=1.0, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1", λ=1.0)
+    inter = AshbaughHatch(weight_special=0.5)
+    @test isapprox(
+        force(inter, dr12, AH_a1, AH_a1, boundary),
+        SVector(16.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
+        atol=1e-9u"kJ * mol^-1 * nm^-1",
+    )
+    @test isapprox(
+        force(inter, dr13, AH_a1, AH_a1, boundary),
+        SVector(-1.375509739, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
+        atol=1e-9u"kJ * mol^-1 * nm^-1",
+    )
+    @test isapprox(
+        potential_energy(inter, dr12, AH_a1, AH_a1, boundary),
+        0.0u"kJ * mol^-1";
+        atol=1e-9u"kJ * mol^-1",
+    )
+    @test isapprox(
+        potential_energy(inter, dr13, AH_a1, AH_a1, boundary),
+        -0.1170417309u"kJ * mol^-1";
+        atol=1e-9u"kJ * mol^-1",
+    )
+
+    AH_a1 = Molly.AshbaughHatchAtom(charge=1.0, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1", λ=0.5)
+    @test isapprox(
+        potential_energy(inter, dr13, AH_a1, AH_a1),
+        -0.058520865u"kJ * mol^-1";
+        atol=1e-9u"kJ * mol^-1",
+    )
+    @test isapprox(
+        force(inter, dr13, AH_a1, AH_a1),
+        SVector(-0.68775486946, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
+        atol=1e-9u"kJ * mol^-1 * nm^-1",
+    )
+
+    c4 = SVector(1.28, 1.0, 1.0)u"nm"
+    dr14 = vector(c1, c4, boundary)
+    @test isapprox(
+        potential_energy(inter, dr14, AH_a1, AH_a1),
+        0.7205987916u"kJ * mol^-1";
+        atol=1e-9u"kJ * mol^-1",
+    )
+    @test isapprox(
+        force(inter, dr14, AH_a1, AH_a1),
+        SVector(52.5306754422, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
+        atol=1e-9u"kJ * mol^-1 * nm^-1",
+    )
+    @test isapprox(
+        potential_energy(inter, dr14,  AH_a1, AH_a1, u"kJ * mol^-1 * nm^-1", true),
+        0.5 * 0.7205987916u"kJ * mol^-1";
+        atol=1e-9u"kJ * mol^-1",
+    )
+    @test isapprox(
+        force(inter, dr14,  AH_a1, AH_a1, u"kJ * mol^-1 * nm^-1", true),
+        SVector(0.5 * 52.5306754422, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
+        atol=1e-9u"kJ * mol^-1 * nm^-1",
+    )
+
     inter = SoftSphere()
     @test isapprox(
         force(inter, dr12, a1, a1),
@@ -106,7 +141,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"Å"
@@ -194,8 +229,8 @@
         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";
+        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(
@@ -209,8 +244,8 @@
         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";
+        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"
@@ -439,67 +474,28 @@
         -108.166724117u"kJ * mol^-1";
         atol=1e-7u"kJ * mol^-1",
     )
-    AH_a1 = AshbaughHatchAtom(charge=1.0, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1", λ=1.0)
-    inter = AshbaughHatch(weight_special=0.5)
-
-    @test isapprox(
-        force(inter, dr12, AH_a1, AH_a1, boundary),
-        SVector(16.0, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
-        atol=1e-9u"kJ * mol^-1 * nm^-1",
-    )
-    @test isapprox(
-        force(inter, dr13, AH_a1, AH_a1, boundary),
-        SVector(-1.375509739, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
-        atol=1e-9u"kJ * mol^-1 * nm^-1",
-    )
-    @test isapprox(
-        potential_energy(inter, dr12, AH_a1, AH_a1, boundary),
-        0.0u"kJ * mol^-1";
-        atol=1e-9u"kJ * mol^-1",
-    )
-    @test isapprox(
-        potential_energy(inter, dr13, AH_a1, AH_a1, boundary),
-        -0.1170417309u"kJ * mol^-1";
-        atol=1e-9u"kJ * mol^-1",
-    )
-
-
-    AH_a1 = AshbaughHatchAtom(charge=1.0, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1", λ=0.5)
-    @test isapprox(
-        potential_energy(inter, dr13,  AH_a1, AH_a1),
-        -0.058520865u"kJ * mol^-1";
-        atol=1e-9u"kJ * mol^-1",
-    )
-
-    @test isapprox(
-        force(inter, dr13,  AH_a1, AH_a1),
-        SVector(-0.68775486946, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
-        atol=1e-9u"kJ * mol^-1 * nm^-1",
-    )
 
-    c4 = SVector(1.28, 1.0, 1.0)u"nm"
-    dr14 = vector(c1, c4, boundary)
-    @test isapprox(
-        potential_energy(inter, dr14,  AH_a1, AH_a1),
-        0.7205987916u"kJ * mol^-1";
-        atol=1e-9u"kJ * mol^-1",
-    )
+    do_shortcut(atom_i, atom_j) = true
 
-    @test isapprox(
-        force(inter, dr14,  AH_a1, AH_a1),
-        SVector(52.5306754422, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
-        atol=1e-9u"kJ * mol^-1 * nm^-1",
-    )
+    for inter in (
+            LennardJones(; shortcut=do_shortcut),
+            LennardJonesSoftCore(α=1, λ=0, p=2; shortcut=do_shortcut),
+            AshbaughHatch(; shortcut=do_shortcut),
+            SoftSphere(; shortcut=do_shortcut),
+            Mie(m=6, n=12; shortcut=do_shortcut),
+            Buckingham(; shortcut=do_shortcut),
+        )
 
-    @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, 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(
-        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",
-    )
+        @test isapprox(
+            potential_energy(inter, dr12, a1, a1),
+            0.0u"kJ * mol^-1";
+            atol=1e-9u"kJ * mol^-1",
+        )
+    end
 end