Merge pull request #189 from ywitzky/YW_Calvados2

Added the force fields necessary for Calvados2

docs/src/documentation.md

@@ -543,12 +543,14 @@ In Molly there are three types of interactions:
 The available pairwise interactions are:
 - [`LennardJones`](@ref)
 - [`LennardJonesSoftCore`](@ref)
+- [`AshbaughHatch`](@ref)
 - [`SoftSphere`](@ref)
 - [`Mie`](@ref)
 - [`Buckingham`](@ref)
 - [`Coulomb`](@ref)
 - [`CoulombSoftCore`](@ref)
 - [`CoulombReactionField`](@ref)
+- [`Yukawa`](@ref)
 - [`Gravity`](@ref)
 
 The available specific interactions are:

src/interactions/coulomb.jl

@@ -1,7 +1,8 @@
 export
     Coulomb,
     CoulombSoftCore,
-    CoulombReactionField
+    CoulombReactionField,
+    Yukawa
 
 const coulomb_const = 138.93545764u"kJ * mol^-1 * nm" # 1 / 4πϵ0
 
@@ -340,3 +341,104 @@ end
         return pe
     end
 end
+
+
+
+
+
+@doc raw"""
+    Yukawa(; cutoff, use_neighbors, weight_special, coulomb_const, kappa)
+
+The Yukawa electrostatic interaction between two atoms.
+
+The potential energy is defined as
+```math
+V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}} \exp(-\kappa r_{ij})
+```
+
+and the force on each atom by 
+
+```math
+F(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}^2} \exp(-\kappa r_{ij})\left(\kappa r_{ij} + 1\right) \vec{r}_{ij}
+```
+"""
+@kwdef struct Yukawa{C, W, T, D} 
+    cutoff::C= NoCutoff()
+    use_neighbors::Bool= false
+    weight_special::W =1
+    coulomb_const::T = coulomb_const
+    kappa::D = 1.0*u"nm^-1"
+end
+
+use_neighbors(inter::Yukawa) = inter.use_neighbors
+
+function Base.zero(yukawa::Yukawa{C, W, T ,D}) where {C, W, T, D}
+    return Yukawa{C, W, T, D}(
+        yukawa.cutoff,
+        yukawa.use_neighbors,
+        yukawa.weight_special,
+        yukawa.coulomb_const,
+        yukawa.kappa
+    )
+end
+
+function Base.:+(c1::Yukawa, c2::Yukawa)
+    return Yukawa(
+        c1.cutoff,
+        c1.use_neighbors,
+        c1.weight_special + c2.weight_special,
+        c1.coulomb_const,
+        c1.kappa
+    )
+end
+
+@inline function force(inter::Yukawa{C},
+                                    dr,
+                                    atom_i,
+                                    atom_j,
+                                    force_units=u"kJ * mol^-1 * nm^-1",
+                                    special=false,
+                                    args...) where C
+    r2 = sum(abs2, dr)
+    cutoff = inter.cutoff
+    coulomb_const = inter.coulomb_const
+    qi, qj = atom_i.charge, atom_j.charge
+    kappa = inter.kappa
+    params = (coulomb_const, qi, qj, kappa)
+
+    f = force_divr_with_cutoff(inter, r2, params, cutoff, force_units)
+    if special
+        return f * dr * inter.weight_special
+    else
+        return f * dr
+    end
+end
+
+function force_divr(::Yukawa, r2, invr2, (coulomb_const, qi, qj, kappa))
+    r = sqrt(r2)
+    return (coulomb_const * qi * qj) / r^3 *exp(-kappa*r) * (kappa*r+1)
+end
+
+@inline function potential_energy(inter::Yukawa{C},
+                                            dr,
+                                            atom_i,
+                                            atom_j,
+                                            energy_units=u"kJ * mol^-1",
+                                            special::Bool=false, args...) where C
+    r2 = sum(abs2, dr)
+    cutoff = inter.cutoff
+    coulomb_const = inter.coulomb_const
+    qi, qj = atom_i.charge, atom_j.charge
+    params = (coulomb_const, qi, qj, inter.kappa)
+
+    pe = potential_with_cutoff(inter, r2, params, cutoff, energy_units)
+    if special
+        return pe * inter.weight_special
+    else
+        return pe
+    end
+end
+
+function potential(::Yukawa, r2, invr2, (coulomb_const, qi, qj, kappa))
+    return (coulomb_const * qi * qj) * √invr2 * exp(-kappa*sqrt(r2))
+end

src/interactions/lennard_jones.jl

@@ -1,6 +1,8 @@
 export
     LennardJones,
-    LennardJonesSoftCore
+    LennardJonesSoftCore,
+    AshbaughHatch,
+    AshbaughHatchAtom
 
 function lj_zero_shortcut(atom_i, atom_j)
     return iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) ||
@@ -13,6 +15,14 @@ function lorentz_σ_mixing(atom_i, atom_j)
     return (atom_i.σ + atom_j.σ) / 2
 end
 
+function lorentz_ϵ_mixing(atom_i, atom_j)
+    return (atom_i.ϵ + atom_j.ϵ) / 2
+end
+
+function lorentz_λ_mixing(atom_i, atom_j)
+    return (atom_i.λ + atom_j.λ) / 2
+end
+
 function geometric_σ_mixing(atom_i, atom_j)
     return sqrt(atom_i.σ * atom_j.σ)
 end
@@ -280,3 +290,159 @@ function potential(::LennardJonesSoftCore, r2, invr2, (σ2, ϵ, σ6_fac))
     six_term = σ2^3 * inv(r2^3 + σ2^3 * σ6_fac)
     return 4ϵ * (six_term ^ 2 - six_term)
 end
+
+
+@doc raw"""
+    AshbaughHatch(; cutoff,use_neighbors,shortcut, ϵ_mixing,σ_mixing, λ_mixing,weight_special)
+
+The AshbaughHatch ($V_{\text{AH}}$) is a modified Lennard-Jones ($V_{\text{LJ}}$) 6-12 interaction between two atoms.
+
+The potential energy is defined by
+```math
+V_{\text{LJ}}(r_{ij}) = 4\varepsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right] \\ 
+```
+
+```math
+V_{\text{AH}}(r_{ij}) =
+    \begin{cases}
+      V_{\text{LJ}}(r_{ij}) +\varepsilon_{ij}(1-λ_{ij}) &,  r_{ij}\leq  2^{1/6}σ  \\
+       λ_{ij}V_{\text{LJ}}(r_{ij})  &,  2^{1/6}σ \leq r_{ij}
+    \end{cases}
+```
+
+and the force on each atom by
+```math
+\vec{F}_{\text{AH}} =
+    \begin{cases}
+      F_{\text{LJ}}(r_{ij})  &,  r_{ij} \leq  2^{1/6}σ  \\
+       λ_{ij}F_{\text{LJ}}(r_{ij})  &,  2^{1/6}σ \leq r_{ij}
+    \end{cases}
+```
+where
+```math
+\begin{aligned}
+\vec{F}_{\text{LJ}}\
+&= \frac{24\varepsilon_{ij}}{r_{ij}^2} \left[2\left(\frac{\sigma_{ij}^{6}}{r_{ij}^{6}}\right)^2 -\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]  \vec{r_{ij}}
+\end{aligned}
+```
+
+If ``\lambda`` is one this gives the standard [`LennardJones`](@ref) potential.
+"""
+@kwdef struct AshbaughHatch{C, H, E, F,L,W} 
+    cutoff::C = NoCutoff()
+    use_neighbors::Bool = false
+    shortcut::H = lj_zero_shortcut
+    ϵ_mixing::E = lorentz_ϵ_mixing
+    σ_mixing::F = lorentz_σ_mixing
+    λ_mixing::L = lorentz_λ_mixing
+    weight_special::W = 1
+end
+
+
+@kwdef struct AshbaughHatchAtom{C, M, S, E, L}
+    index::Int = 1
+    charge::C = 0.0
+    mass::M =1.0u"g/mol"
+    σ::S =0.0u"nm"
+    ϵ::E =0.0u"kJ * mol^-1"
+    λ::L=1.0
+end
+
+use_neighbors(inter::AshbaughHatch) = inter.use_neighbors
+
+function Base.zero(lj::AshbaughHatch{C, H, E, F,L,W}) where {C, H, E, F,L,W}
+    return AshbaughHatch{C, H, E, F,L,W}(
+        lj.cutoff,
+        lj,use_neighbors,
+        lj.shortcut,
+        lj.σ_mixing,
+        lj.ϵ_mixing,
+        lj.λ_mixing,
+        zero(W),
+    )
+end
+
+function Base.:+(l1::AshbaughHatch{C, H, E, F,L,W},
+                 l2::AshbaughHatch{C, H, E, F,L,W}) where {C, H, E, F,L,W}
+    return AshbaughHatch{C, H, E, F,L,W}(
+        l1.cutoff,
+        l1.use_neighbors,
+        l1.shortcut,
+        l1.σ_mixing,
+        l1.ϵ_mixing,
+        l1.λ_mixing,
+        l1.weight_special + l2.weight_special,
+    )
+end
+
+@inline function force(inter::AshbaughHatch{S, C},
+                                    dr,
+                                    atom_i,
+                                    atom_j,
+                                    force_units=u"kJ * mol^-1 * nm^-1",
+                                    special::Bool=false, args...) where {S, C}
+    if inter.shortcut(atom_i, atom_j)
+        return ustrip.(zero(dr)) * force_units
+    end
+
+    # Lorentz-Berthelot mixing rules use the arithmetic average for σ
+    # Otherwise use the geometric average
+    ϵ = inter.ϵ_mixing(atom_i, atom_j)
+    σ = inter.σ_mixing(atom_i, atom_j)
+    λ = inter.λ_mixing(atom_i, atom_j)
+
+    cutoff = inter.cutoff
+    r2 = sum(abs2, dr)
+    σ2 = σ^2
+    params = (σ2, ϵ, λ)
+
+    f = force_divr_with_cutoff(inter, r2, params, cutoff,force_units)
+    if special
+        return f * dr * inter.weight_special
+    else
+        return f * dr
+    end
+end
+
+@inline function force_divr(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ))
+    if r2 < 2^(1/3)*σ2
+        return  force_divr(LennardJones(), r2, invr2, (σ2, ϵ))
+    else
+        return λ*force_divr(LennardJones(), r2, invr2, (σ2, ϵ)) 
+    end
+end
+
+@inline function potential_energy(inter::AshbaughHatch{S, C},
+                                            dr,
+                                            atom_i,
+                                            atom_j,
+                                            energy_units=u"kJ * mol^-1",
+                                            special::Bool=false,
+                                            args...) where {S, C}
+    if inter.shortcut(atom_i, atom_j)
+        return ustrip(zero(dr[1])) * energy_units
+    end
+    ϵ = inter.ϵ_mixing(atom_i, atom_j)     ### probably not needed for Calvados2 use case
+    σ = inter.σ_mixing(atom_i, atom_j)
+    λ = inter.λ_mixing(atom_i, atom_j)
+
+    cutoff = inter.cutoff
+    r2 = sum(abs2, dr)
+    σ2 = σ^2
+    params = (σ2, ϵ, λ)
+
+    pe = potential_with_cutoff(inter, r2, params, cutoff,energy_units)
+    if special
+        return pe * inter.weight_special
+    else
+        return pe
+    end
+end
+
+@inline function potential(::AshbaughHatch, r2, invr2, (σ2, ϵ, λ))
+    if r2 < 2^(1/3)*σ2
+        return   potential(LennardJones(), r2, invr2, (σ2, ϵ)) + ϵ*(1-λ) 
+    else
+        return λ* potential(LennardJones(), r2, invr2, (σ2, ϵ))
+    end
+end