diff --git a/src/interactions/coulomb.jl b/src/interactions/coulomb.jl index 9cc424e1..e2baa9f4 100644 --- a/src/interactions/coulomb.jl +++ b/src/interactions/coulomb.jl @@ -118,15 +118,15 @@ 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{1}{4\pi\epsilon_0} \frac{q_iq_j}{r_Q^{1/6}} +V(r_{ij}) = \lambda\frac{1}{4\pi\epsilon_0} \frac{q_iq_j}{r_Q^{1/6}} ``` and the force on each atom by ```math -\vec{F}_i = \frac{1}{4\pi\epsilon_0} \frac{q_iq_jr_{ij}^5}{r_Q^{7/6}}\frac{\vec{r_{ij}}}{r_{ij}} +\vec{F}_i = \lambda\frac{1}{4\pi\epsilon_0} \frac{q_iq_jr_{ij}^5}{r_Q^{7/6}}\frac{\vec{r_{ij}}}{r_{ij}} ``` where ```math -r_{Q} = (\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}})+r_{ij}^6) +r_{Q} = \left(\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}}\right)+r_{ij}^6 ``` and ```math @@ -194,7 +194,7 @@ end σ6 = inter.σ_mixing(atom_i, atom_j)^6 ϵ = inter.ϵ_mixing(atom_i, atom_j) C6 = 4 * ϵ * σ6 - params = (ke, qi, qj, C6 * σ6, C6, inter.σ6_fac) + params = (ke, qi, qj, C6 * σ6, C6, inter.σ6_fac, inter.λ) f = force_cutoff(cutoff, inter, r, params) fdr = (f / r) * dr @@ -205,10 +205,10 @@ end end end -function pairwise_force(::CoulombSoftCoreBeutler, r, (ke, qi, qj, C12, C6, σ6_fac)) +function pairwise_force(::CoulombSoftCoreBeutler, r, (ke, qi, qj, C12, C6, σ6_fac, λ)) r3 = r^3 R = ((σ6_fac*(C12/C6))+(r3*r3))*sqrt(cbrt(((σ6_fac*(C12/C6))+(r3*r3)))) - return ke * ((qi*qj)/R) * (r3*r*r) + return λ * ke * ((qi*qj)/R) * (r3*r*r) end @inline function potential_energy(inter::CoulombSoftCoreBeutler, @@ -225,7 +225,7 @@ end σ6 = inter.σ_mixing(atom_i, atom_j)^6 ϵ = inter.ϵ_mixing(atom_i, atom_j) C6 = 4 * ϵ * σ6 - params = (ke, qi, qj, C6 *σ6, C6, inter.σ6_fac) + params = (ke, qi, qj, C6 *σ6, C6, inter.σ6_fac, inter.λ) pe = pe_cutoff(cutoff, inter, r, params) if special @@ -235,9 +235,9 @@ end end end -function pairwise_pe(::CoulombSoftCoreBeutler, r, (ke, qi, qj, C12, C6, σ6_fac)) +function pairwise_pe(::CoulombSoftCoreBeutler, r, (ke, qi, qj, C12, C6, σ6_fac, λ)) R = sqrt(cbrt((σ6_fac*(C12/C6))+r^6)) - return ke * ((qi * qj)/R) + return λ * ke * ((qi * qj)/R) end @doc raw""" @@ -321,7 +321,7 @@ end cutoff = inter.cutoff ke = inter.coulomb_const qi, qj = atom_i.charge, atom_j.charge - params = (ke, qi, qj, inter.σQ, inter.σ6_fac) + params = (ke, qi, qj, inter.σQ, inter.σ6_fac, inter.λ) f = force_cutoff(cutoff, inter, r, params) fdr = (f / r) * dr @@ -332,13 +332,13 @@ end end end -function pairwise_force(::CoulombSoftCoreGapsys, r, (ke, qi, qj, σQ, σ6_fac)) +function pairwise_force(::CoulombSoftCoreGapsys, r, (ke, qi, qj, σQ, σ6_fac, λ)) qij = qi * qj R = σ6_fac*(oneunit(r)+(σQ*abs(qij))) if r >= R - return ke * (qij/(r^2)) + return λ * ke * (qij/(r^2)) elseif r < R - return ke * (-(((2*qij)/(R^3)) * r) + ((3*qij)/(R^2))) + return λ * ke * (-(((2*qij)/(R^3)) * r) + ((3*qij)/(R^2))) end end @@ -353,7 +353,7 @@ end cutoff = inter.cutoff ke = inter.coulomb_const qi, qj = atom_i.charge, atom_j.charge - params = (ke, qi, qj, inter.σQ, inter.σ6_fac) + params = (ke, qi, qj, inter.σQ, inter.σ6_fac, inter.λ) pe = pe_cutoff(cutoff, inter, r, params) if special @@ -363,13 +363,13 @@ end end end -function pairwise_pe(::CoulombSoftCoreGapsys, r, (ke, qi, qj, σQ, σ6_fac)) +function pairwise_pe(::CoulombSoftCoreGapsys, r, (ke, qi, qj, σQ, σ6_fac, λ)) qij = qi * qj R = σ6_fac*(oneunit(r)+(σQ*abs(qij))) if r >= R - return ke * (qij/r) + return λ * ke * (qij/r) elseif r < R - return ke * (((qij/(R^3))*(r^2))-(((3*qij)/(R^2))*r)+((3*qij)/R)) + return λ * ke * (((qij/(R^3))*(r^2))-(((3*qij)/(R^2))*r)+((3*qij)/R)) end end diff --git a/src/interactions/lennard_jones.jl b/src/interactions/lennard_jones.jl index 87f7630a..1e598c23 100644 --- a/src/interactions/lennard_jones.jl +++ b/src/interactions/lennard_jones.jl @@ -149,15 +149,15 @@ 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}} +V(r_{ij}) = \lambda\left(\frac{C^{(12)}}{r_{LJ}^{12}} - \frac{C^{(6)}}{r_{LJ}^{6}}\right ``` 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}} +\vec{F}_i = \lambda\left(\left(\frac{12C^{(12)}}{r_{LJ}^{13}} - \frac{6C^{(6)}}{r_{LJ}^7}\right)\left(\frac{r_{ij}}{r_{LJ}}\right)^5\right) \frac{\vec{r_{ij}}}{r_{ij}} ``` where ```math -r_{LJ} = (\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}}+r^6)^{1/6} +r_{LJ} = \left(\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}}+r^6\right)^{1/6} ``` and ```math @@ -228,7 +228,7 @@ end cutoff = inter.cutoff r = norm(dr) C6 = 4 * ϵ * σ6 - params = (C6 * σ6, C6, inter.σ6_fac) + params = (C6 * σ6, C6, inter.σ6_fac, inter.λ) f = force_cutoff(cutoff, inter, r, params) fdr = (f / r) * dr @@ -239,10 +239,10 @@ end end end -function pairwise_force(::LennardJonesSoftCoreBeutler, r, (C12, C6, σ6_fac)) +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) + return λ*(((12*C12)/(R6*R6*R)) - ((6*C6)/(R6*R)))*((r/R)^5) end @inline function potential_energy(inter::LennardJonesSoftCoreBeutler, @@ -261,7 +261,7 @@ end cutoff = inter.cutoff r = norm(dr) C6 = 4 * ϵ * σ6 - params = (C6 * σ6, C6, inter.σ6_fac) + params = (C6 * σ6, C6, inter.σ6_fac, inter.λ) pe = pe_cutoff(cutoff, inter, r, params) if special @@ -271,9 +271,9 @@ end end end -function pairwise_pe(::LennardJonesSoftCoreBeutler, r, (C12, C6, σ6_fac)) +function pairwise_pe(::LennardJonesSoftCoreBeutler, r, (C12, C6, σ6_fac, λ)) R6 = (σ6_fac*(C12/C6))+r^6 - return ((C12/(R6*R6)) - (C6/(R6))) + return λ*((C12/(R6*R6)) - (C6/(R6))) end @doc raw""" @@ -368,7 +368,7 @@ end cutoff = inter.cutoff r = norm(dr) C6 = 4 * ϵ * σ6 - params = (C6 * σ6, C6) + params = (C6 * σ6, C6, inter.λ) f = force_cutoff(cutoff, inter, r, params) fdr = (f / r) * dr @@ -379,16 +379,16 @@ end end end -function pairwise_force(inter::LennardJonesSoftCoreGapsys, r, (C12, C6)) +function pairwise_force(inter::LennardJonesSoftCoreGapsys, r, (C12, C6, λ)) R = inter.α*sqrt(cbrt((26*(C12/C6)*(1-inter.λ)/7))) r6 = r^6 invR = inv(R) invR2 = invR^2 invR6 = invR^6 if r >= R - return (((12*C12)/(r6*r6*r))-((6*C6)/(r6*r))) + return λ * (((12*C12)/(r6*r6*r))-((6*C6)/(r6*r))) elseif r < R - return (((-156*C12*(invR6*invR6*invR2)) + (42*C6*(invR2*invR6)))*r + + return λ * (((-156*C12*(invR6*invR6*invR2)) + (42*C6*(invR2*invR6)))*r + (168*C12*(invR6*invR6*invR)) - (48*C6*(invR6*invR))) end end @@ -409,7 +409,7 @@ end cutoff = inter.cutoff r = norm(dr) C6 = 4 * ϵ * σ6 - params = (C6 * σ6, C6) + params = (C6 * σ6, C6, inter.λ) pe = pe_cutoff(cutoff, inter, r, params) if special @@ -419,16 +419,16 @@ end end end -function pairwise_pe(inter::LennardJonesSoftCoreGapsys, r, (C12, C6)) +function pairwise_pe(inter::LennardJonesSoftCoreGapsys, r, (C12, C6, λ)) R = inter.α*sqrt(cbrt((26*(C12/C6)*(1-inter.λ)/7))) r6 = r^6 invR = inv(R) invR2 = invR^2 invR6 = invR^6 if r >= R - return (C12/(r6*r6))-(C6/(r6)) + return λ * (C12/(r6*r6))-(C6/(r6)) elseif r < R - return ((78*C12*(invR6*invR6*invR2)) - (21*C6*(invR2*invR6)))*(r^2) - + return λ * ((78*C12*(invR6*invR6*invR2)) - (21*C6*(invR2*invR6)))*(r^2) - ((168*C12*(invR6*invR6*invR)) - (48*C6*(invR6*invR)))*r + (91*C12*(invR6*invR6)) - (28*C6*(invR6)) end diff --git a/test/interactions.jl b/test/interactions.jl index 67a05696..ec8702ef 100644 --- a/test/interactions.jl +++ b/test/interactions.jl @@ -47,44 +47,44 @@ 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", ) inter = LennardJonesSoftCoreGapsys(α=0.85, λ=0.5) @test isapprox( force(inter, dr14, a1, a1), - SVector(516.8457758610879, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + SVector(258.42288793054365, 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.37550973892212, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + SVector(-0.68775486946106, 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), - 51.814929518108826u"kJ * mol^-1"; + 45.35853723577515u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) @test isapprox( potential_energy(inter, dr13, a1, a1), - -0.1170417308807374u"kJ * mol^-1"; + -0.12971227169036878u"kJ * mol^-1"; atol=1e-9u"kJ * mol^-1", ) @@ -231,44 +231,44 @@ inter = CoulombSoftCoreBeutler(α=0.3, λ=0.5) @test isapprox( force(inter, dr13, a1, a1), - SVector(842.06163558501, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + SVector(421.030817792505, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; atol=1e-5u"kJ * mol^-1 * nm^-1", ) @test isapprox( force(inter, dr14, a1, a1), - SVector(57.48819886473348, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + SVector(28.74409943236674, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; atol=1e-5u"kJ * mol^-1 * nm^-1", ) @test isapprox( potential_energy(inter, dr13, a1, a1), - 345.81678703197457u"kJ * mol^-1"; + 172.90839351598729u"kJ * mol^-1"; atol=1e-5u"kJ * mol^-1", ) @test isapprox( potential_energy(inter, dr14, a1, a1), - 634.3822744723304u"kJ * mol^-1"; + 317.1911372361652u"kJ * mol^-1"; atol=1e-5u"kJ * mol^-1", ) inter = CoulombSoftCoreGapsys(α=0.3, λ=0.5, σQ=1u"nm") @test isapprox( force(inter, dr13, a1, a1), - SVector(731.0108657043696, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + SVector(365.5054328521848, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; atol=1e-5u"kJ * mol^-1 * nm^-1", ) @test isapprox( force(inter, dr14, a1, a1), - SVector(1276.8008892849266, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; + SVector(638.4004446424633, 0.0, 0.0)u"kJ * mol^-1 * nm^-1"; atol=1e-5u"kJ * mol^-1 * nm^-1", ) @test isapprox( potential_energy(inter, dr13, a1, a1), - 341.80053925707637u"kJ * mol^-1"; + 170.90026962853818u"kJ * mol^-1"; atol=1e-5u"kJ * mol^-1", ) @test isapprox( potential_energy(inter, dr14, a1, a1), - 642.9723025054708u"kJ * mol^-1"; + 321.4861512527354u"kJ * mol^-1"; atol=1e-5u"kJ * mol^-1", )