From 8ef2d6afd916ab150b56cd54e0845d1169a5773a Mon Sep 17 00:00:00 2001 From: Thomas Kemmer Date: Thu, 28 Nov 2024 17:02:13 +0100 Subject: [PATCH] FF: don't cache atom positions Fixes: #134 Signed-off-by: Thomas Kemmer --- src/forcefields/common/bend_component.jl | 17 +++++------ src/forcefields/common/nonbonded_component.jl | 22 ++++++--------- src/forcefields/common/stretch_component.jl | 14 ++++------ src/forcefields/common/torsion_component.jl | 28 ++++++++----------- 4 files changed, 34 insertions(+), 47 deletions(-) diff --git a/src/forcefields/common/bend_component.jl b/src/forcefields/common/bend_component.jl index 08497a34..fbb3bb74 100644 --- a/src/forcefields/common/bend_component.jl +++ b/src/forcefields/common/bend_component.jl @@ -6,11 +6,8 @@ export θ₀::T k::T a1::Atom{T} - a1r::Vector3{T} a2::Atom{T} - a2r::Vector3{T} a3::Atom{T} - a3r::Vector3{T} end @auto_hash_equals mutable struct QuadraticBendComponent{T<:Real} <: AbstractForceFieldComponent{T} @@ -100,9 +97,9 @@ function setup!(qbc::QuadraticBendComponent{T}) where {T<:Real} QuadraticAngleBend( T(θ₀_factor*only(qab.theta0)), T(k_factor*only(qab.k)), - a1, a1.r, - a2, a2.r, - a3, a3.r + a1, + a2, + a3, )) end end @@ -117,8 +114,8 @@ function update!(qbc::QuadraticBendComponent{T}) where {T<:Real} end @inline function compute_energy(qab::QuadraticAngleBend{T})::T where {T<:Real} - v1 = qab.a1r .- qab.a2r - v2 = qab.a3r .- qab.a2r + v1 = qab.a1.r .- qab.a2.r + v2 = qab.a3.r .- qab.a2.r sq_length = squared_norm(v1) * squared_norm(v2) @@ -149,8 +146,8 @@ end function compute_forces!(qab::QuadraticAngleBend{T}) where {T<:Real} # calculate the vectors between the atoms and normalize if possible - v1 = qab.a1r .- qab.a2r - v2 = qab.a3r .- qab.a2r + v1 = qab.a1.r .- qab.a2.r + v2 = qab.a3.r .- qab.a2.r v1_length = norm(v1) v2_length = norm(v2) diff --git a/src/forcefields/common/nonbonded_component.jl b/src/forcefields/common/nonbonded_component.jl index 1058e9d7..bd2fed8b 100644 --- a/src/forcefields/common/nonbonded_component.jl +++ b/src/forcefields/common/nonbonded_component.jl @@ -86,9 +86,7 @@ end distance::T scaling_factor::T a1::Atom{T} - a1r::Vector3{T} a2::Atom{T} - a2r::Vector3{T} switching_function::CubicSwitchingFunction{T} end @@ -98,9 +96,7 @@ end scaling_factor::T distance_dependent_dielectric::Bool a1::Atom{T} - a1r::Vector3{T} a2::Atom{T} - a2r::Vector3{T} switching_function::CubicSwitchingFunction{T} end @@ -162,8 +158,8 @@ end params.B_ij, T(distance), scaling_factor, - atom_1, atom_1.r, - atom_2, atom_2.r, + atom_1, + atom_2, switching_function ) ) @@ -416,8 +412,8 @@ function update!(nbc::NonBondedComponent{T}) where {T<:Real} lj_candidate[3], vicinal_pair ? scaling_es_1_4 : T(1.0), distance_dependent_dielectric, - atom_1, atom_1.r, - atom_2, atom_2.r, + atom_1, + atom_2, es_switching_function ) push!(electrostatic_interactions, es) @@ -443,8 +439,8 @@ function update!(nbc::NonBondedComponent{T}) where {T<:Real} only(h_params.B), T(lj_candidate[3]), T(1.0), - atom_1, atom_1.r, - atom_2, atom_2.r, + atom_1, + atom_2, vdw_switching_function ) ) @@ -516,7 +512,7 @@ function compute_energy!(nbc::NonBondedComponent{T})::T where {T<:Real} end function compute_forces!(lji::LennardJonesInteraction{T, 12, 6}) where {T<:Real} - direction = lji.a1r .- lji.a2r + direction = lji.a1.r .- lji.a2.r sq_distance = squared_norm(direction) @@ -552,7 +548,7 @@ function compute_forces!(lji::LennardJonesInteraction{T, 12, 6}) where {T<:Real} end function compute_forces!(hb::LennardJonesInteraction{T, 12, 10}) where {T<:Real} - direction = hb.a1r .- hb.a2r + direction = hb.a1.r .- hb.a2.r sq_distance = squared_norm(direction) @@ -592,7 +588,7 @@ function compute_forces!(hb::LennardJonesInteraction{T, 12, 10}) where {T<:Real} end function compute_forces!(esi::ElectrostaticInteraction{T}) where {T<:Real} - direction = esi.a1r .- esi.a2r + direction = esi.a1.r .- esi.a2.r sq_distance = squared_norm(direction) diff --git a/src/forcefields/common/stretch_component.jl b/src/forcefields/common/stretch_component.jl index 60c2c63e..85226063 100644 --- a/src/forcefields/common/stretch_component.jl +++ b/src/forcefields/common/stretch_component.jl @@ -6,9 +6,7 @@ export r0::T k::T a1::Atom{T} - a1r::Vector3{T} a2::Atom{T} - a2r::Vector3{T} end @auto_hash_equals mutable struct QuadraticStretchComponent{T<:Real} <: AbstractForceFieldComponent{T} @@ -60,7 +58,7 @@ function setup!(qsc::QuadraticStretchComponent{T}) where {T} if has_flag(bond, :TYPE__HYDROGEN) # skip hydrogen bonds - stretches[i] = QuadraticBondStretch(one(T), zero(T), a1, a1.r, a2, a2.r) + stretches[i] = QuadraticBondStretch(one(T), zero(T), a1, a2) end qbs = coalesce( @@ -82,13 +80,13 @@ function setup!(qsc::QuadraticStretchComponent{T}) where {T} end # we don't want to get any force or energy component from this stretch - QuadraticBondStretch(one(T), zero(T), a1, a1.r, a2, a2.r) + QuadraticBondStretch(one(T), zero(T), a1, a2) else QuadraticBondStretch( T(r0_factor*only(qbs.r0)), T(k_factor*only(qbs.k)), - a1, a1.r, - a2, a2.r + a1, + a2 ) end end @@ -101,7 +99,7 @@ function update!(qsc::QuadraticStretchComponent{T}) where {T<:Real} end @inline function compute_energy(qbs::QuadraticBondStretch{T})::T where {T<:Real} - d = distance(qbs.a1r, qbs.a2r) + d = distance(qbs.a1.r, qbs.a2.r) qbs.k * (d - qbs.r0)^2 end @@ -117,7 +115,7 @@ function compute_energy!(qsc::QuadraticStretchComponent{T})::T where {T<:Real} end function compute_forces!(qbs::QuadraticBondStretch{T}) where {T<:Real} - direction = qbs.a1r .- qbs.a2r + direction = qbs.a1.r .- qbs.a2.r distance = norm(direction) if distance == zero(T) diff --git a/src/forcefields/common/torsion_component.jl b/src/forcefields/common/torsion_component.jl index 4081f78f..7da81392 100644 --- a/src/forcefields/common/torsion_component.jl +++ b/src/forcefields/common/torsion_component.jl @@ -8,13 +8,9 @@ export f::Vector{Int} div::Vector{Int} a1::Atom{T} - a1r::Vector3{T} a2::Atom{T} - a2r::Vector3{T} a3::Atom{T} - a3r::Vector3{T} a4::Atom{T} - a4r::Vector3{T} end @auto_hash_equals mutable struct TorsionComponent{T<:Real} <: AbstractForceFieldComponent{T} @@ -113,10 +109,10 @@ function _try_assign_torsion!( ϕ₀_factor .* getproperty.(pts, :phi0), getproperty.(pts, :f), getproperty.(pts, :div), - a1, a1.r, - a2, a2.r, - a3, a3.r, - a4, a4.r + a1, + a2, + a3, + a4 ) ) end @@ -262,10 +258,10 @@ end @inline function compute_energy(pt::CosineTorsion{T})::T where {T<:Real} energy = zero(T) - a23 = pt.a3r .- pt.a2r + a23 = pt.a3.r .- pt.a2.r - cross2321 = normalize(cross(a23, pt.a1r .- pt.a2r)) - cross2334 = normalize(cross(a23, pt.a4r .- pt.a3r)) + cross2321 = normalize(cross(a23, pt.a1.r .- pt.a2.r)) + cross2334 = normalize(cross(a23, pt.a4.r .- pt.a3.r)) if !isnan(cross2321[1]) && !isnan(cross2334[1]) cos_ϕ = clamp(dot(cross2321, cross2334), T(-1.0), T(1.0)) @@ -296,9 +292,9 @@ function compute_energy!(tc::TorsionComponent{T})::T where {T<:Real} end function compute_forces!(ct::CosineTorsion{T}) where {T<:Real} - a21 = ct.a1r .- ct.a2r - a23 = ct.a3r .- ct.a2r - a34 = ct.a4r .- ct.a3r + a21 = ct.a1.r .- ct.a2.r + a23 = ct.a3.r .- ct.a2.r + a34 = ct.a4.r .- ct.a3.r cross2321 = cross(a23, a21) cross2334 = cross(a23, a34) @@ -325,8 +321,8 @@ function compute_forces!(ct::CosineTorsion{T}) where {T<:Real} ∂E∂ϕ *= -1 end - a13 = ct.a3r .- ct.a1r - a24 = ct.a4r .- ct.a2r + a13 = ct.a3.r .- ct.a1.r + a24 = ct.a4.r .- ct.a2.r dEdt = (∂E∂ϕ / (length_cross2321^2 * norm(a23)) * cross(cross2321, a23)) dEdu = -(∂E∂ϕ / (length_cross2334^2 * norm(a23)) * cross(cross2334, a23))