From ca51e4fcf9c714862c72edd3454afa1819265c5f Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Fri, 21 Jun 2024 16:43:45 +0100 Subject: [PATCH] fix: Correct tabulated energies and forces when applying internal scaling factors (#1908) --- src/classes/pairPotential.cpp | 143 +++++++++++-------- src/classes/pairPotential.h | 60 ++++---- src/classes/potentialMap.cpp | 94 ++++++------ src/classes/potentialMap.h | 28 ++-- src/classes/speciesAtom.cpp | 4 +- src/classes/speciesTorsion.cpp | 4 +- src/classes/speciesTorsion.h | 2 +- src/gui/models/pairPotentialModel.cpp | 2 +- src/io/export/pairPotential.cpp | 2 +- src/kernels/energy.cpp | 4 +- src/kernels/energy.h | 2 +- src/kernels/force.cpp | 8 +- src/kernels/force.h | 4 +- src/main/pairPotentials.cpp | 8 +- src/modules/energy/process.cpp | 6 +- src/modules/forces/process.cpp | 6 +- tests/pairPotentials/CMakeLists.txt | 1 + tests/pairPotentials/analytic.cpp | 16 +-- tests/pairPotentials/scaleFactors.cpp | 196 ++++++++++++++++++++++++++ 19 files changed, 411 insertions(+), 179 deletions(-) create mode 100644 tests/pairPotentials/scaleFactors.cpp diff --git a/src/classes/pairPotential.cpp b/src/classes/pairPotential.cpp index eab4ed0c67..d74b3c8dcd 100644 --- a/src/classes/pairPotential.cpp +++ b/src/classes/pairPotential.cpp @@ -13,17 +13,22 @@ PairPotential::CoulombTruncationScheme PairPotential::coulombTruncationScheme_ = PairPotential::ShiftedCoulombTruncation; PairPotential::ShortRangeTruncationScheme PairPotential::shortRangeTruncationScheme_ = PairPotential::ShiftedShortRangeTruncation; -bool PairPotential::includeAtomTypeCharges_ = false; +bool PairPotential::includeCoulombPotential_ = false; PairPotential::PairPotential(std::string_view nameI, std::string_view nameJ) - : nameI_(nameI), nameJ_(nameJ), totalPotentialInterpolation_(totalPotential_), derivativeInterpolation_(derivative_), - interactionPotential_{Functions1D::Form::None, ""}, potentialFunction_{Functions1D::Form::None, {}} + : nameI_(nameI), nameJ_(nameJ), totalShortRangePotentialInterpolation_(referenceShortRangePotential_), + coulombPotentialInterpolation_(coulombPotential_), totalPotentialInterpolation_(totalPotential_), + totalShortRangeDerivativeInterpolation_(totalShortRangeDerivative_), coulombDerivativeInterpolation_(coulombDerivative_), + totalDerivativeInterpolation_(totalDerivative_), interactionPotential_{Functions1D::Form::None, ""}, + potentialFunction_{Functions1D::Form::None, {}} { } PairPotential::PairPotential(std::string_view nameI, std::string_view nameJ, const InteractionPotential &potential) - : nameI_(nameI), nameJ_(nameJ), interactionPotential_(potential), totalPotentialInterpolation_(totalPotential_), - derivativeInterpolation_(derivative_) + : nameI_(nameI), nameJ_(nameJ), interactionPotential_(potential), + totalShortRangePotentialInterpolation_(referenceShortRangePotential_), coulombPotentialInterpolation_(coulombPotential_), + totalPotentialInterpolation_(totalPotential_), totalShortRangeDerivativeInterpolation_(totalShortRangeDerivative_), + coulombDerivativeInterpolation_(coulombDerivative_), totalDerivativeInterpolation_(totalDerivative_) { potentialFunction_.setFormAndParameters(interactionPotential_.form(), interactionPotential_.parameters()); } @@ -57,11 +62,11 @@ void PairPotential::setShortRangeTruncationScheme(PairPotential::ShortRangeTrunc // Return short-ranged truncation scheme PairPotential::ShortRangeTruncationScheme PairPotential::shortRangeTruncationScheme() { return shortRangeTruncationScheme_; } -// Set whether atom type charges should be included in the generated potential -void PairPotential::setIncludeAtomTypeCharges(bool b) { includeAtomTypeCharges_ = b; } +// Set whether Coulomb contributions should be included in the generated potential +void PairPotential::setIncludeCoulombPotential(bool b) { includeCoulombPotential_ = b; } -// Return whether atom type charges should be included in the generated potential -bool PairPotential::includeAtomTypeCharges() { return includeAtomTypeCharges_; } +// Return whether Coulomb contributions should be included in the generated potential +bool PairPotential::includeCoulombPotential() { return includeCoulombPotential_; } // Set Coulomb truncation scheme void PairPotential::setCoulombTruncationScheme(PairPotential::CoulombTruncationScheme scheme) @@ -81,12 +86,12 @@ void PairPotential::setData1DNames() { totalPotential_.setTag(fmt::format("{}-{}", nameI_, nameJ_)); - additionalPotential_.setTag(fmt::format("{}-{} (Add)", nameI_, nameJ_)); + additionalShortRangePotential_.setTag(fmt::format("{}-{} (Add SR)", nameI_, nameJ_)); + referenceShortRangePotential_.setTag(fmt::format("{}-{} (Ref SR)", nameI_, nameJ_)); - shortRangePotential_.setTag(fmt::format("{}-{} (SR)", nameI_, nameJ_)); - coulombPotential_.setTag(fmt::format("{}-{} (Coul)", nameI_, nameJ_)); + coulombPotential_.setTag(fmt::format("{}-{} (Elec)", nameI_, nameJ_)); - derivative_.setTag(fmt::format("{}-{} (dU/dr)", nameI_, nameJ_)); + totalDerivative_.setTag(fmt::format("{}-{} (dU/dr)", nameI_, nameJ_)); } // Set names reflecting target atom types for potential @@ -125,8 +130,8 @@ void PairPotential::setInteractionPotentialForm(Functions1D::Form form) // Return interaction potential const InteractionPotential &PairPotential::interactionPotential() const { return interactionPotential_; } -// Return atom typeCharge product (if including Coulomb terms) -double PairPotential::atomTypeChargeProduct() const { return atomTypeChargeProduct_; } +// Return local charge product (if including Coulomb terms) +double PairPotential::localChargeProduct() const { return localChargeProduct_; } /* * Tabulated PairPotential @@ -166,22 +171,32 @@ double PairPotential::analyticShortRangeForce(double r, PairPotential::ShortRang void PairPotential::updateTotals() { // Update total energy - for (auto &&[total, sr, coul, add] : zip(totalPotential_.values(), shortRangePotential_.values(), - coulombPotential_.values(), additionalPotential_.values())) - total = sr + coul + add; + for (auto &&[total, totalSR, refSR, addSR, coul] : + zip(totalPotential_.values(), totalShortRangePotential_.values(), referenceShortRangePotential_.values(), + additionalShortRangePotential_.values(), coulombPotential_.values())) + { + totalSR = refSR + addSR; + total = totalSR + coul; + } - // Recalculate interpolation + // Recalculate interpolations + totalShortRangePotentialInterpolation_.interpolate(Interpolator::ThreePointInterpolation); + coulombPotentialInterpolation_.interpolate(Interpolator::ThreePointInterpolation); totalPotentialInterpolation_.interpolate(Interpolator::ThreePointInterpolation); - // Calculate derivative of total potential - derivative_ = Derivative::derivative(totalPotential_); + // Calculate derivatives + totalShortRangeDerivative_ = Derivative::derivative(totalShortRangePotential_); + coulombDerivative_ = Derivative::derivative(coulombPotential_); + totalDerivative_ = Derivative::derivative(totalPotential_); - // Update interpolation - derivativeInterpolation_.interpolate(Interpolator::ThreePointInterpolation); + // Update interpolations for derivatives + totalShortRangeDerivativeInterpolation_.interpolate(Interpolator::ThreePointInterpolation); + coulombDerivativeInterpolation_.interpolate(Interpolator::ThreePointInterpolation); + totalDerivativeInterpolation_.interpolate(Interpolator::ThreePointInterpolation); } // Generate energy and force tables -void PairPotential::tabulate(double maxR, double delta, double qi, double qj) +void PairPotential::tabulate(double maxR, double delta, double chargeProduct) { // Determine delta_ = delta; @@ -191,43 +206,45 @@ void PairPotential::tabulate(double maxR, double delta, double qi, double qj) // Precalculate some quantities shortRangeEnergyAtCutoff_ = analyticShortRangeEnergy(range_, PairPotential::NoShortRangeTruncation); shortRangeForceAtCutoff_ = analyticShortRangeForce(range_, PairPotential::NoShortRangeTruncation); - atomTypeChargeProduct_ = qi * qj; + localChargeProduct_ = includeCoulombPotential_ ? chargeProduct : 0.0; // Set up containers const auto nPoints = int(range_ / delta); - shortRangePotential_.initialise(nPoints); + referenceShortRangePotential_.initialise(nPoints); for (auto n = 0; n < nPoints; ++n) - shortRangePotential_.xAxis()[n] = n * delta_; - coulombPotential_ = shortRangePotential_; - additionalPotential_ = shortRangePotential_; - totalPotential_ = shortRangePotential_; - derivative_ = shortRangePotential_; - - // Tabulate short-range and coulomb energies - for (auto &&[r, sr, coul] : zip(shortRangePotential_.xAxis(), shortRangePotential_.values(), coulombPotential_.values())) + referenceShortRangePotential_.xAxis()[n] = n * delta_; + coulombPotential_ = referenceShortRangePotential_; + additionalShortRangePotential_ = referenceShortRangePotential_; + totalShortRangePotential_ = referenceShortRangePotential_; + totalPotential_ = referenceShortRangePotential_; + totalDerivative_ = referenceShortRangePotential_; + + // Tabulate reference short-range and coulomb energies + for (auto &&[r, refSR, coul] : + zip(referenceShortRangePotential_.xAxis(), referenceShortRangePotential_.values(), coulombPotential_.values())) { - sr = analyticShortRangeEnergy(r); - coul = includeAtomTypeCharges_ ? analyticCoulombEnergy(atomTypeChargeProduct_, r) : 0.0; + refSR = analyticShortRangeEnergy(r); + coul = analyticCoulombEnergy(localChargeProduct_, r); } // Since the first point at r = 0.0 risks being a nan, set it to ten times the second point instead - shortRangePotential_.value(0) = 10.0 * shortRangePotential_.value(1); + referenceShortRangePotential_.value(0) = 10.0 * referenceShortRangePotential_.value(1); coulombPotential_.value(0) = 10.0 * coulombPotential_.value(1); // Ensure additional potential is set to zero and update full potential - std::fill(additionalPotential_.values().begin(), additionalPotential_.values().end(), 0); + std::fill(additionalShortRangePotential_.values().begin(), additionalShortRangePotential_.values().end(), 0); // Update totals updateTotals(); } -// Add supplied function to the short-range potential -void PairPotential::addShortRangePotential(const Function1DWrapper &potential, bool overwriteExisting) +// Add supplied function to the reference short-range potential +void PairPotential::addToReferenceShortRangePotential(const Function1DWrapper &potential, bool overwriteExisting) { if (overwriteExisting) - std::fill(shortRangePotential_.values().begin(), shortRangePotential_.values().end(), 0.0); + std::fill(referenceShortRangePotential_.values().begin(), referenceShortRangePotential_.values().end(), 0.0); - for (auto &&[r, sr] : zip(shortRangePotential_.xAxis(), shortRangePotential_.values())) + for (auto &&[r, sr] : zip(referenceShortRangePotential_.xAxis(), referenceShortRangePotential_.values())) sr += potential.y(r); // Update totals @@ -247,25 +264,32 @@ double PairPotential::energy(double r) return totalPotentialInterpolation_.y(r, r * rDelta_); } +double PairPotential::energy(double r, double elecScale, double srScale) +{ + assert(r >= 0); + + return totalShortRangePotentialInterpolation_.y(r, r * rDelta_) * srScale + + coulombPotentialInterpolation_.y(r, r * rDelta_) * elecScale; +} // Return analytic potential at specified r, including Coulomb term from local atomtype charges -double PairPotential::analyticEnergy(double r) const +double PairPotential::analyticEnergy(double r, double elecScale, double srScale) const { if (r > range_) return 0.0; // Short-range potential and Coulomb contribution - return analyticShortRangeEnergy(r) + analyticCoulombEnergy(atomTypeChargeProduct_, r); + return srScale * analyticShortRangeEnergy(r) + elecScale * analyticCoulombEnergy(localChargeProduct_, r); } // Return analytic potential at specified r, including Coulomb term from supplied charge product -double PairPotential::analyticEnergy(double qiqj, double r, double elecScale, double vdwScale, +double PairPotential::analyticEnergy(double qiqj, double r, double elecScale, double srScale, PairPotential::CoulombTruncationScheme truncation) const { if (r > range_) return 0.0; - return analyticShortRangeEnergy(r) * vdwScale + analyticCoulombEnergy(qiqj, r, truncation) * elecScale; + return analyticShortRangeEnergy(r) * srScale + analyticCoulombEnergy(qiqj, r, truncation) * elecScale; } // Return analytic coulomb potential energy of specified charges @@ -285,27 +309,34 @@ double PairPotential::force(double r) { assert(r >= 0); - return -derivativeInterpolation_.y(r, r * rDelta_); + return -totalDerivativeInterpolation_.y(r, r * rDelta_); +} +double PairPotential::force(double r, double elecScale, double srScale) +{ + assert(r >= 0); + + return -(totalShortRangeDerivativeInterpolation_.y(r, r * rDelta_) * srScale + + coulombDerivativeInterpolation_.y(r, r * rDelta_) * elecScale); } // Return analytic force at specified r -double PairPotential::analyticForce(double r) const +double PairPotential::analyticForce(double r, double elecScale, double srScale) const { if (r > range_) return 0.0; // Short-range potential and Coulomb contribution - return analyticShortRangeForce(r) + analyticCoulombForce(atomTypeChargeProduct_, r); + return srScale * analyticShortRangeForce(r) + elecScale * analyticCoulombForce(localChargeProduct_, r); } // Return analytic force at specified r, including Coulomb term from supplied charge product -double PairPotential::analyticForce(double qiqj, double r, double elecScale, double vdwScale, +double PairPotential::analyticForce(double qiqj, double r, double elecScale, double srScale, PairPotential::CoulombTruncationScheme truncation) const { if (r > range_) return 0.0; - return analyticShortRangeForce(r) * vdwScale + analyticCoulombForce(qiqj, r) * elecScale; + return analyticShortRangeForce(r) * srScale + analyticCoulombForce(qiqj, r) * elecScale; } // Return analytic coulomb force of specified charges @@ -336,21 +367,21 @@ double PairPotential::analyticCoulombForce(double qiqj, double r, PairPotential: const Data1D &PairPotential::totalPotential() const { return totalPotential_; } // Return full tabulated derivative -const Data1D &PairPotential::derivative() const { return derivative_; } +const Data1D &PairPotential::derivative() const { return totalDerivative_; } // Return short-range potential -const Data1D &PairPotential::shortRangePotential() const { return shortRangePotential_; } +const Data1D &PairPotential::shortRangePotential() const { return referenceShortRangePotential_; } // Return Coulomb potential const Data1D &PairPotential::coulombPotential() const { return coulombPotential_; } // Return additional potential -const Data1D &PairPotential::additionalPotential() const { return additionalPotential_; } +const Data1D &PairPotential::additionalPotential() const { return additionalShortRangePotential_; } // Zero additional potential void PairPotential::resetAdditionalPotential() { - std::fill(additionalPotential_.values().begin(), additionalPotential_.values().end(), 0.0); + std::fill(additionalShortRangePotential_.values().begin(), additionalShortRangePotential_.values().end(), 0.0); updateTotals(); } @@ -358,7 +389,7 @@ void PairPotential::resetAdditionalPotential() // Set additional potential void PairPotential::setAdditionalPotential(Data1D &newUAdditional) { - additionalPotential_ = newUAdditional; + additionalShortRangePotential_ = newUAdditional; updateTotals(); } diff --git a/src/classes/pairPotential.h b/src/classes/pairPotential.h index fec9644939..bf0123052f 100644 --- a/src/classes/pairPotential.h +++ b/src/classes/pairPotential.h @@ -47,8 +47,8 @@ class PairPotential : Serialisable<> double shortRangeEnergyAtCutoff_{0.0}; // Short-range force at cutoff distance (used by truncation scheme) double shortRangeForceAtCutoff_{0.0}; - // Whether atom type charges should be included in the generated potential - static bool includeAtomTypeCharges_; + // Whether Coulomb contributions should be included in the generated potential + static bool includeCoulombPotential_; // Truncation scheme to apply to Coulomb part of potential static CoulombTruncationScheme coulombTruncationScheme_; @@ -57,10 +57,10 @@ class PairPotential : Serialisable<> static void setShortRangeTruncationScheme(ShortRangeTruncationScheme scheme); // Return short-ranged truncation scheme static ShortRangeTruncationScheme shortRangeTruncationScheme(); - // Set whether atom type charges should be included in the generated potential - static void setIncludeAtomTypeCharges(bool b); - // Return whether atom type charges should be included in the generated potential - static bool includeAtomTypeCharges(); + // Set whether Coulomb contributions should be included in the generated potential + static void setIncludeCoulombPotential(bool b); + // Return whether Coulomb contributions should be included in the generated potential + static bool includeCoulombPotential(); // Set Coulomb truncation scheme static void setCoulombTruncationScheme(CoulombTruncationScheme scheme); // Return Coulomb truncation scheme @@ -76,7 +76,7 @@ class PairPotential : Serialisable<> InteractionPotential interactionPotential_; Function1DWrapper potentialFunction_; // Charge product (if including Coulomb terms) - double atomTypeChargeProduct_{0.0}; + double localChargeProduct_{0.0}; private: // Set Data1D names from source AtomTypes @@ -96,8 +96,8 @@ class PairPotential : Serialisable<> void setInteractionPotentialForm(Functions1D::Form form); // Return interaction potential const InteractionPotential &interactionPotential() const; - // Return atom typeCharge product (if including Coulomb terms) - double atomTypeChargeProduct() const; + // Return local charge product (if including Coulomb terms) + double localChargeProduct() const; /* * Tabulated Potential @@ -108,17 +108,19 @@ class PairPotential : Serialisable<> // Distance between points in tabulated potentials double delta_{-1.0}, rDelta_{0.0}; // Tabulated original potential, calculated from AtomType parameters - Data1D shortRangePotential_, coulombPotential_; - // Additional potential, generated by some means - Data1D additionalPotential_; - // Total potential (original plus additional), + Data1D referenceShortRangePotential_, coulombPotential_; + // Additional short-range potential, generated by some means + Data1D additionalShortRangePotential_; + // Total short-range potential: reference plus additional + Data1D totalShortRangePotential_; + // Total potential: reference and additional short-range plus Coulomb Data1D totalPotential_; - // Interpolation of full potential - Interpolator totalPotentialInterpolation_; - // Tabulated derivative of full potential - Data1D derivative_; - // Interpolation of derivative of full potential - Interpolator derivativeInterpolation_; + // Interpolators for potentials + Interpolator totalShortRangePotentialInterpolation_, coulombPotentialInterpolation_, totalPotentialInterpolation_; + // Tabulated derivative + Data1D totalDerivative_, totalShortRangeDerivative_, coulombDerivative_; + // Interpolators for derivatives + Interpolator totalShortRangeDerivativeInterpolation_, coulombDerivativeInterpolation_, totalDerivativeInterpolation_; private: // Return analytic short range potential energy @@ -132,19 +134,20 @@ class PairPotential : Serialisable<> public: // Generate energy and force tables - void tabulate(double maxR, double delta, double qi = 0.0, double qj = 0.0); - // Add supplied function to the short-range potential - void addShortRangePotential(const Function1DWrapper &potential, bool overwriteExisting = false); + void tabulate(double maxR, double delta, double chargeProduct = 0.0); + // Add supplied function to the reference short-range potential + void addToReferenceShortRangePotential(const Function1DWrapper &potential, bool overwriteExisting = false); // Return range of potential double range() const; // Return spacing between points double delta() const; // Return potential at specified r double energy(double r); - // Return analytic potential at specified r, including Coulomb term from local atomtype charges - double analyticEnergy(double r) const; + double energy(double r, double elecScale, double srScale); + // Return analytic potential at specified r, including Coulomb term from local charge product + double analyticEnergy(double r, double elecScale, double srScale) const; // Return analytic potential at specified r, including Coulomb term from supplied charge product - double analyticEnergy(double qiqj, double r, double elecScale = 1.0, double vdwScale = 1.0, + double analyticEnergy(double qiqj, double r, double elecScale, double srScale, PairPotential::CoulombTruncationScheme truncation = PairPotential::coulombTruncationScheme()) const; // Return analytic coulomb potential energy of specified charge product double @@ -152,10 +155,11 @@ class PairPotential : Serialisable<> PairPotential::CoulombTruncationScheme truncation = PairPotential::coulombTruncationScheme()) const; // Return derivative of potential at specified r double force(double r); - // Return analytic force at specified r, including Coulomb term from local atomtype charges - double analyticForce(double r) const; + double force(double r, double elecScale, double srScale); + // Return analytic force at specified r, including Coulomb term from local charge product + double analyticForce(double r, double elecScale, double srScale) const; // Return analytic force at specified r, including Coulomb term from supplied charge product - double analyticForce(double qiqj, double r, double elecScale = 1.0, double vdwScale = 1.0, + double analyticForce(double qiqj, double r, double elecScale, double srScale, PairPotential::CoulombTruncationScheme truncation = PairPotential::coulombTruncationScheme()) const; // Return analytic coulomb force of specified charge product double diff --git a/src/classes/potentialMap.cpp b/src/classes/potentialMap.cpp index b4f80f98c8..fb72073b76 100644 --- a/src/classes/potentialMap.cpp +++ b/src/classes/potentialMap.cpp @@ -76,13 +76,13 @@ double PotentialMap::energy(const Atom &i, const Atom &j, double r) const // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the // interpolated potential auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}]; - return pp->energy(r) + (PairPotential::includeAtomTypeCharges() + return pp->energy(r) + (PairPotential::includeCoulombPotential() ? 0 : pp->analyticCoulombEnergy(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r)); } -// Return energy between Atoms at distance specified, scaling electrostatic and van der Waals components -double PotentialMap::energy(const Atom &i, const Atom &j, double r, double elecScale, double vdwScale) const +// Return energy between Atoms at distance specified, scaling electrostatic and short-range components +double PotentialMap::energy(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const { assert(r >= 0.0); assert(i.speciesAtom() && j.speciesAtom()); @@ -90,9 +90,9 @@ double PotentialMap::energy(const Atom &i, const Atom &j, double r, double elecS // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the // interpolated potential auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}]; - return PairPotential::includeAtomTypeCharges() - ? pp->energy(r) * elecScale - : pp->energy(r) * vdwScale + + return PairPotential::includeCoulombPotential() + ? pp->energy(r, elecScale, srScale) + : pp->energy(r) * srScale + pp->analyticCoulombEnergy(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r) * elecScale; } @@ -106,11 +106,11 @@ double PotentialMap::energy(const SpeciesAtom *i, const SpeciesAtom *j, double r // interpolated potential auto *pp = potentialMatrix_[{i->atomType()->index(), j->atomType()->index()}]; return pp->energy(r) + - (PairPotential::includeAtomTypeCharges() ? 0 : pp->analyticCoulombEnergy(i->charge() * j->charge(), r)); + (PairPotential::includeCoulombPotential() ? 0 : pp->analyticCoulombEnergy(i->charge() * j->charge(), r)); } -// Return energy between SpeciesAtoms at distance specified, scaling electrostatic and van der Waals components -double PotentialMap::energy(const SpeciesAtom *i, const SpeciesAtom *j, double r, double elecScale, double vdwScale) const +// Return energy between SpeciesAtoms at distance specified, scaling electrostatic and short-range components +double PotentialMap::energy(const SpeciesAtom *i, const SpeciesAtom *j, double r, double elecScale, double srScale) const { assert(r >= 0.0); assert(i && j); @@ -118,37 +118,37 @@ double PotentialMap::energy(const SpeciesAtom *i, const SpeciesAtom *j, double r // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the // interpolated potential auto *pp = potentialMatrix_[{i->atomType()->index(), j->atomType()->index()}]; - return PairPotential::includeAtomTypeCharges() - ? pp->energy(r) * elecScale - : pp->energy(r) * vdwScale + pp->analyticCoulombEnergy(i->charge() * j->charge(), r) * elecScale; + return PairPotential::includeCoulombPotential() + ? pp->energy(r, elecScale, srScale) + : pp->energy(r) * srScale + pp->analyticCoulombEnergy(i->charge() * j->charge(), r) * elecScale; } // Return analytic energy between Atom types at distance specified -double PotentialMap::analyticEnergy(const Atom *i, const Atom *j, double r) const +double PotentialMap::analyticEnergy(const Atom &i, const Atom &j, double r) const { assert(r >= 0.0); assert(i && j); // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being local to the atom // types - auto *pp = potentialMatrix_[{i->masterTypeIndex(), j->masterTypeIndex()}]; - return PairPotential::includeAtomTypeCharges() - ? pp->analyticEnergy(r) - : pp->analyticEnergy(i->speciesAtom()->charge() * j->speciesAtom()->charge(), r); + auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}]; + return PairPotential::includeCoulombPotential() + ? pp->analyticEnergy(r, 1.0, 1.0) + : pp->analyticEnergy(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r, 1.0, 1.0); } -// Return analytic energy between Atom types at distance specified, scaling electrostatic and van der Waals components -double PotentialMap::analyticEnergy(const Atom *i, const Atom *j, double r, double elecScale, double vdwScale) const +// Return analytic energy between Atom types at distance specified, scaling electrostatic and short-range components +double PotentialMap::analyticEnergy(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const { assert(r >= 0.0); assert(i && j); // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being local to the atom // types - auto *pp = potentialMatrix_[{i->masterTypeIndex(), j->masterTypeIndex()}]; - return PairPotential::includeAtomTypeCharges() - ? pp->analyticEnergy(r) * elecScale - : pp->analyticEnergy(i->speciesAtom()->charge() * j->speciesAtom()->charge(), r, elecScale, vdwScale); + auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}]; + return PairPotential::includeCoulombPotential() + ? pp->analyticEnergy(r, elecScale, srScale) + : pp->analyticEnergy(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r, elecScale, srScale); } // Return force between Atoms at distance specified @@ -157,20 +157,20 @@ double PotentialMap::force(const Atom &i, const Atom &j, double r) const // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the // interpolated potential auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}]; - return PairPotential::includeAtomTypeCharges() + return PairPotential::includeCoulombPotential() ? pp->force(r) : pp->force(r) + pp->analyticCoulombForce(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r); } -// Return force between Atoms at distance specified, scaling electrostatic and van der Waals components -double PotentialMap::force(const Atom &i, const Atom &j, double r, double elecScale, double vdwScale) const +// Return force between Atoms at distance specified, scaling electrostatic and short-range components +double PotentialMap::force(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const { // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the // interpolated potential auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}]; - return PairPotential::includeAtomTypeCharges() - ? pp->force(r) * elecScale - : pp->force(r) * vdwScale + + return PairPotential::includeCoulombPotential() + ? pp->force(r, elecScale, srScale) + : pp->force(r) * srScale + pp->analyticCoulombForce(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r) * elecScale; } @@ -183,12 +183,12 @@ double PotentialMap::force(const SpeciesAtom *i, const SpeciesAtom *j, double r) // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the // interpolated potential auto *pp = potentialMatrix_[{i->atomType()->index(), j->atomType()->index()}]; - return PairPotential::includeAtomTypeCharges() ? pp->force(r) - : pp->force(r) + pp->analyticCoulombForce(i->charge() * j->charge(), r); + return PairPotential::includeCoulombPotential() ? pp->force(r) + : pp->force(r) + pp->analyticCoulombForce(i->charge() * j->charge(), r); } -// Return force between SpeciesAtoms at distance specified, scaling electrostatic and van der Waals components -double PotentialMap::force(const SpeciesAtom *i, const SpeciesAtom *j, double r, double elecScale, double vdwScale) const +// Return force between SpeciesAtoms at distance specified, scaling electrostatic and short-range components +double PotentialMap::force(const SpeciesAtom *i, const SpeciesAtom *j, double r, double elecScale, double srScale) const { assert(r >= 0.0); assert(i && j); @@ -196,13 +196,13 @@ double PotentialMap::force(const SpeciesAtom *i, const SpeciesAtom *j, double r, // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the // interpolated potential auto *pp = potentialMatrix_[{i->atomType()->index(), j->atomType()->index()}]; - return PairPotential::includeAtomTypeCharges() - ? pp->force(r) * elecScale - : pp->force(r) * vdwScale + pp->analyticCoulombForce(i->charge() * j->charge(), r) * elecScale; + return PairPotential::includeCoulombPotential() + ? pp->force(r, elecScale, srScale) + : pp->force(r) * srScale + pp->analyticCoulombForce(i->charge() * j->charge(), r) * elecScale; } // Return analytic force between Atom types at distance specified -double PotentialMap::analyticForce(const Atom *i, const Atom *j, double r) const +double PotentialMap::analyticForce(const Atom &i, const Atom &j, double r) const { assert(r >= 0.0); assert(i && j); @@ -210,14 +210,14 @@ double PotentialMap::analyticForce(const Atom *i, const Atom *j, double r) const // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the // interpolated potential - auto *pp = potentialMatrix_[{i->masterTypeIndex(), j->masterTypeIndex()}]; - return PairPotential::includeAtomTypeCharges() - ? pp->analyticForce(r) - : pp->analyticForce(i->speciesAtom()->charge() * j->speciesAtom()->charge(), r); + auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}]; + return PairPotential::includeCoulombPotential() + ? pp->analyticForce(r, 1.0, 1.0) + : pp->analyticForce(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r, 1.0, 1.0); } -// Return analytic force between Atom types at distance specified, scaling electrostatic and van der Waals components -double PotentialMap::analyticForce(const Atom *i, const Atom *j, double r, double elecScale, double vdwScale) const +// Return analytic force between Atom types at distance specified, scaling electrostatic and short-range components +double PotentialMap::analyticForce(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const { assert(r >= 0.0); assert(i && j); @@ -225,8 +225,8 @@ double PotentialMap::analyticForce(const Atom *i, const Atom *j, double r, doubl // Check to see whether Coulomb terms should be calculated from atomic charges, rather than them being included in the // interpolated potential - auto *pp = potentialMatrix_[{i->masterTypeIndex(), j->masterTypeIndex()}]; - return PairPotential::includeAtomTypeCharges() - ? pp->analyticForce(r) * elecScale - : pp->analyticForce(i->speciesAtom()->charge() * j->speciesAtom()->charge(), r, elecScale, vdwScale); + auto *pp = potentialMatrix_[{i.masterTypeIndex(), j.masterTypeIndex()}]; + return PairPotential::includeCoulombPotential() + ? pp->analyticForce(r, elecScale, srScale) + : pp->analyticForce(i.speciesAtom()->charge() * j.speciesAtom()->charge(), r, elecScale, srScale); } diff --git a/src/classes/potentialMap.h b/src/classes/potentialMap.h index 0287e3f65e..0723916b18 100644 --- a/src/classes/potentialMap.h +++ b/src/classes/potentialMap.h @@ -44,26 +44,26 @@ class PotentialMap public: // Return energy between Atoms at distance specified double energy(const Atom &i, const Atom &j, double r) const; - // Return energy between Atoms at distance specified, scaling electrostatic and van der Waals components - double energy(const Atom &i, const Atom &j, double r, double elecScale, double vdwScale) const; + // Return energy between Atoms at distance specified, scaling electrostatic and short-range components + double energy(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const; // Return energy between SpeciesAtoms at distance specified double energy(const SpeciesAtom *i, const SpeciesAtom *j, double r) const; - // Return energy between SpeciesAtoms at distance specified, scaling electrostatic and van der Waals components - double energy(const SpeciesAtom *i, const SpeciesAtom *j, double r, double elecScale, double vdwScale) const; + // Return energy between SpeciesAtoms at distance specified, scaling electrostatic and short-range components + double energy(const SpeciesAtom *i, const SpeciesAtom *j, double r, double elecScale, double srScale) const; // Return analytic energy between Atom types at distance specified - double analyticEnergy(const Atom *i, const Atom *j, double r) const; - // Return analytic energy between Atom types at distance specified, scaling electrostatic and van der Waals components - double analyticEnergy(const Atom *i, const Atom *j, double r, double elecScale, double vdwScale) const; + double analyticEnergy(const Atom &i, const Atom &j, double r) const; + // Return analytic energy between Atom types at distance specified, scaling electrostatic and short-range components + double analyticEnergy(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const; // Return force between Atoms at distance specified double force(const Atom &i, const Atom &j, double r) const; - // Return force between Atoms at distance specified, scaling electrostatic and van der Waals components - double force(const Atom &i, const Atom &j, double r, double elecScale, double vdwScale) const; + // Return force between Atoms at distance specified, scaling electrostatic and short-range components + double force(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const; // Return force between SpeciesAtoms at distance specified double force(const SpeciesAtom *i, const SpeciesAtom *j, double r) const; - // Return force between SpeciesAtoms at distance specified, scaling electrostatic and van der Waals components - double force(const SpeciesAtom *i, const SpeciesAtom *j, double r, double elecScale, double vdwScale) const; + // Return force between SpeciesAtoms at distance specified, scaling electrostatic and short-range components + double force(const SpeciesAtom *i, const SpeciesAtom *j, double r, double elecScale, double srScale) const; // Return analytic force between Atom types at distance specified - double analyticForce(const Atom *i, const Atom *j, double r) const; - // Return analytic force between Atom types at distance specified, scaling electrostatic and van der Waals components - double analyticForce(const Atom *i, const Atom *j, double r, double elecScale, double vdwScale) const; + double analyticForce(const Atom &i, const Atom &j, double r) const; + // Return analytic force between Atom types at distance specified, scaling electrostatic and short-range components + double analyticForce(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const; }; diff --git a/src/classes/speciesAtom.cpp b/src/classes/speciesAtom.cpp index 42cf7e1b08..ca5b71d4ea 100644 --- a/src/classes/speciesAtom.cpp +++ b/src/classes/speciesAtom.cpp @@ -220,12 +220,12 @@ void SpeciesAtom::setScaledInteractions() scaledInteractions_.clear(); std::function addInteractionFunction = - [&](SpeciesAtom *j, SpeciesAtom::ScaledInteraction scaledType, double elecScale, double vdwScale) + [&](SpeciesAtom *j, SpeciesAtom::ScaledInteraction scaledType, double elecScale, double srScale) { auto it = std::find_if(scaledInteractions_.begin(), scaledInteractions_.end(), [j](const auto &p) { return p.first == j; }); if (it == scaledInteractions_.end()) - scaledInteractions_.emplace_back(j, ScaledInteractionDefinition{scaledType, elecScale, vdwScale}); + scaledInteractions_.emplace_back(j, ScaledInteractionDefinition{scaledType, elecScale, srScale}); }; /* diff --git a/src/classes/speciesTorsion.cpp b/src/classes/speciesTorsion.cpp index 0a04864c08..b4b852a3a3 100644 --- a/src/classes/speciesTorsion.cpp +++ b/src/classes/speciesTorsion.cpp @@ -192,12 +192,12 @@ bool SpeciesTorsion::isSelected() const */ // Set 1-4 scaling factors -bool SpeciesTorsion::set14ScalingFactors(double elecScale, double vdwScale) +bool SpeciesTorsion::set14ScalingFactors(double elecScale, double srScale) { if (masterTerm_) return Messenger::error("Refused to set 1-4 scaling factors since master parameters are referenced.\n"); electrostatic14Scaling_ = elecScale; - vdw14Scaling_ = vdwScale; + vdw14Scaling_ = srScale; return true; } diff --git a/src/classes/speciesTorsion.h b/src/classes/speciesTorsion.h index edac397053..ce46c535db 100644 --- a/src/classes/speciesTorsion.h +++ b/src/classes/speciesTorsion.h @@ -83,7 +83,7 @@ class SpeciesTorsion : public SpeciesIntra public: // Set 1-4 scaling factors - bool set14ScalingFactors(double elecScale, double vdwScale); + bool set14ScalingFactors(double elecScale, double srScale); // Set electrostatic 1-4 scaling factor for the interaction bool setElectrostatic14Scaling(double scaling); // Return electrostatic 1-4 scaling factor for the interaction diff --git a/src/gui/models/pairPotentialModel.cpp b/src/gui/models/pairPotentialModel.cpp index d72e52c428..495d6d23c4 100644 --- a/src/gui/models/pairPotentialModel.cpp +++ b/src/gui/models/pairPotentialModel.cpp @@ -56,7 +56,7 @@ QVariant PairPotentialModel::data(const QModelIndex &index, int role) const case (Columns::NameJColumn): return QString::fromStdString(std::string(pp->nameJ())); case (Columns::ChargeProductColumn): - return PairPotential::includeAtomTypeCharges() ? QString::number(pp->atomTypeChargeProduct()) : QString("--"); + return PairPotential::includeCoulombPotential() ? QString::number(pp->localChargeProduct()) : QString("--"); case (Columns::ShortRangeFormColumn): return QString::fromStdString(Functions1D::forms().keyword(pp->interactionPotential().form())); case (Columns::ShortRangeParametersColumn): diff --git a/src/io/export/pairPotential.cpp b/src/io/export/pairPotential.cpp index 5b5fb48648..b889859acc 100644 --- a/src/io/export/pairPotential.cpp +++ b/src/io/export/pairPotential.cpp @@ -39,7 +39,7 @@ bool PairPotentialExportFileFormat::exportBlock(LineParser &parser, PairPotentia for (auto &&[r, tot, deriv, sr, coul, add] : zip(shortRange.xAxis(), total.values(), derivative.values(), shortRange.values(), coulomb.values(), additional.values())) if (!parser.writeLineF("{:10.6e} {:12.6e} {:12.6e} {:12.6e} {:12.6e} {:12.6e} {:12.6e}\n", r, tot, deriv, - sr + coul, add, pp->analyticEnergy(r), pp->analyticForce(r))) + sr + coul, add, pp->analyticEnergy(r, 1.0, 1.0), pp->analyticForce(r, 1.0, 1.0))) return false; return true; diff --git a/src/kernels/energy.cpp b/src/kernels/energy.cpp index aecd053854..d37cb96b45 100644 --- a/src/kernels/energy.cpp +++ b/src/kernels/energy.cpp @@ -67,9 +67,9 @@ EnergyKernel::EnergyKernel(const Configuration *cfg, const ProcessPool &procPool double EnergyKernel::pairPotentialEnergy(const Atom &i, const Atom &j, double r) const { return potentialMap_.energy(i, j, r); } // Return PairPotential energy between atoms, scaling electrostatic and van der Waals components -double EnergyKernel::pairPotentialEnergy(const Atom &i, const Atom &j, double r, double elecScale, double vdwScale) const +double EnergyKernel::pairPotentialEnergy(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const { - return potentialMap_.energy(i, j, r, elecScale, vdwScale); + return potentialMap_.energy(i, j, r, elecScale, srScale); } /* diff --git a/src/kernels/energy.h b/src/kernels/energy.h index c8f0f5ac63..5153c4e19b 100644 --- a/src/kernels/energy.h +++ b/src/kernels/energy.h @@ -85,7 +85,7 @@ class EnergyKernel : public GeometryKernel // Return PairPotential energy between atoms virtual double pairPotentialEnergy(const Atom &i, const Atom &j, double r) const; // Return PairPotential energy between atoms, scaling electrostatic and van der Waals components - virtual double pairPotentialEnergy(const Atom &i, const Atom &j, double r, double elecScale, double vdwScale) const; + virtual double pairPotentialEnergy(const Atom &i, const Atom &j, double r, double elecScale, double srScale) const; /* * PairPotential Terms diff --git a/src/kernels/force.cpp b/src/kernels/force.cpp index 32c5b59777..b3dd9c975a 100644 --- a/src/kernels/force.cpp +++ b/src/kernels/force.cpp @@ -43,7 +43,7 @@ void ForceKernel::forcesWithoutMim(const Atom &i, int indexI, const Atom &j, int // Calculate inter-particle forces between Atoms provided, scaling electrostatic and van der Waals components void ForceKernel::forcesWithoutMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f, double elecScale, - double vdwScale) const + double srScale) const { auto vij = j.r() - i.r(); auto distanceSq = vij.magnitudeSq(); @@ -51,7 +51,7 @@ void ForceKernel::forcesWithoutMim(const Atom &i, int indexI, const Atom &j, int return; auto r = sqrt(distanceSq); vij /= r; - vij *= potentialMap_.force(i, j, r, elecScale, vdwScale); + vij *= potentialMap_.force(i, j, r, elecScale, srScale); f[indexI] -= vij; f[indexJ] += vij; } @@ -72,7 +72,7 @@ void ForceKernel::forcesWithMim(const Atom &i, int indexI, const Atom &j, int in // Calculate inter-particle forces between Atoms provided, scaling electrostatic and van der Waals components void ForceKernel::forcesWithMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f, double elecScale, - double vdwScale) const + double srScale) const { auto vij = box_->minimumVector(i.r(), j.r()); auto distanceSq = vij.magnitudeSq(); @@ -80,7 +80,7 @@ void ForceKernel::forcesWithMim(const Atom &i, int indexI, const Atom &j, int in return; auto r = sqrt(distanceSq); vij /= r; - vij *= potentialMap_.force(i, j, r, elecScale, vdwScale); + vij *= potentialMap_.force(i, j, r, elecScale, srScale); f[indexI] -= vij; f[indexJ] += vij; } diff --git a/src/kernels/force.h b/src/kernels/force.h index a45a3773ce..ce977a07b3 100644 --- a/src/kernels/force.h +++ b/src/kernels/force.h @@ -50,12 +50,12 @@ class ForceKernel : public GeometryKernel void forcesWithoutMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f) const; // Calculate inter-particle forces between Atoms provided, scaling electrostatic and van der Waals components void forcesWithoutMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f, double elecScale, - double vdwScale) const; + double srScale) const; // Calculate inter-particle forces between Atoms provided void forcesWithMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f) const; // Calculate inter-particle forces between Atoms provided, scaling electrostatic and van der Waals components void forcesWithMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f, double elecScale, - double vdwScale) const; + double srScale) const; // Calculate forces between two cells void cellToCellPairPotentialForces(const Cell *cell, const Cell *otherCell, bool applyMim, ForceVector &f) const; diff --git a/src/main/pairPotentials.cpp b/src/main/pairPotentials.cpp index 3b006b579b..915988fcd4 100644 --- a/src/main/pairPotentials.cpp +++ b/src/main/pairPotentials.cpp @@ -105,7 +105,7 @@ bool Dissolve::updatePairPotentials(std::optional useCombinationRulesHint) auto useCombinationRules = useCombinationRulesHint.value_or(useCombinationRules_); // Set the charge handling for all pair potentials - PairPotential::setIncludeAtomTypeCharges(atomTypeChargeSource_); + PairPotential::setIncludeCoulombPotential(atomTypeChargeSource_); // First step - remove any pair potentials which reference non-existent atom types pairPotentials_.erase(std::remove_if(pairPotentials_.begin(), pairPotentials_.end(), @@ -156,7 +156,7 @@ bool Dissolve::updatePairPotentials(std::optional useCombinationRulesHint) // Re-tabulate the potentials to account for changes in charge inclusion/exclusion, range etc. as well as parameters for (auto &&[at1, at2, pot] : pairPotentials_) { - pot->tabulate(pairPotentialRange_, pairPotentialDelta_, at1->charge(), at2->charge()); + pot->tabulate(pairPotentialRange_, pairPotentialDelta_, at1->charge() * at2->charge()); } // Third step - apply any overrides @@ -196,10 +196,10 @@ bool Dissolve::updatePairPotentials(std::optional useCombinationRulesHint) case (PairPotentialOverride::PairPotentialOverrideType::Off): break; case (PairPotentialOverride::PairPotentialOverrideType::Add): - pp->addShortRangePotential(overridePotential); + pp->addToReferenceShortRangePotential(overridePotential); break; case (PairPotentialOverride::PairPotentialOverrideType::Replace): - pp->addShortRangePotential(overridePotential, true); + pp->addToReferenceShortRangePotential(overridePotential, true); break; } diff --git a/src/modules/energy/process.cpp b/src/modules/energy/process.cpp index 5f39361dff..d46a979634 100644 --- a/src/modules/energy/process.cpp +++ b/src/modules/energy/process.cpp @@ -182,9 +182,9 @@ Module::ExecutionResult EnergyModule::process(ModuleContext &moduleContext) // Get intramolecular scaling of atom pair auto &&[scalingType, elec14, vdw14] = i->scaling(j); if (scalingType == SpeciesAtom::ScaledInteraction::NotScaled) - correctSelfEnergy += potentialMap.analyticEnergy(i, j, r); + correctSelfEnergy += potentialMap.analyticEnergy(*i, *j, r); else if (scalingType == SpeciesAtom::ScaledInteraction::Scaled) - correctSelfEnergy += potentialMap.analyticEnergy(i, j, r, elec14, vdw14); + correctSelfEnergy += potentialMap.analyticEnergy(*i, *j, r, elec14, vdw14); } } @@ -207,7 +207,7 @@ Module::ExecutionResult EnergyModule::process(ModuleContext &moduleContext) if (r > cutoff) continue; - correctInterEnergy += potentialMap.analyticEnergy(i, j, r); + correctInterEnergy += potentialMap.analyticEnergy(*i, *j, r); } } } diff --git a/src/modules/forces/process.cpp b/src/modules/forces/process.cpp index 4a9962721a..db8c9f193f 100644 --- a/src/modules/forces/process.cpp +++ b/src/modules/forces/process.cpp @@ -104,9 +104,9 @@ Module::ExecutionResult ForcesModule::process(ModuleContext &moduleContext) vecij /= r; if (scalingType == SpeciesAtom::ScaledInteraction::NotScaled) - vecij *= potentialMap.analyticForce(molN->atom(ii), molN->atom(jj), r); + vecij *= potentialMap.analyticForce(*molN->atom(ii), *molN->atom(jj), r); else if (scalingType == SpeciesAtom::ScaledInteraction::Scaled) - vecij *= potentialMap.analyticForce(molN->atom(ii), molN->atom(jj), r, elec14, vdw14); + vecij *= potentialMap.analyticForce(*molN->atom(ii), *molN->atom(jj), r, elec14, vdw14); fInter[offsetN + ii] -= vecij; fInter[offsetN + jj] += vecij; @@ -136,7 +136,7 @@ Module::ExecutionResult ForcesModule::process(ModuleContext &moduleContext) auto r = sqrt(magjisq); vecij /= r; - vecij *= potentialMap.analyticForce(i, j, r); + vecij *= potentialMap.analyticForce(*i, *j, r); fInter[offsetN + ii] -= vecij; fInter[offsetM + jj] += vecij; diff --git a/tests/pairPotentials/CMakeLists.txt b/tests/pairPotentials/CMakeLists.txt index 4e37ef6ebc..9ee381fdf1 100644 --- a/tests/pairPotentials/CMakeLists.txt +++ b/tests/pairPotentials/CMakeLists.txt @@ -1,2 +1,3 @@ dissolve_add_test(SRC analytic.cpp) dissolve_add_test(SRC overrides.cpp) +dissolve_add_test(SRC scaleFactors.cpp) diff --git a/tests/pairPotentials/analytic.cpp b/tests/pairPotentials/analytic.cpp index 15b3e621ad..65aaa17154 100644 --- a/tests/pairPotentials/analytic.cpp +++ b/tests/pairPotentials/analytic.cpp @@ -7,10 +7,10 @@ namespace UnitTest { -class PairPotentialsTest : public ::testing::Test +class PairPotentialsAnalyticTest : public ::testing::Test { public: - PairPotentialsTest() + PairPotentialsAnalyticTest() { PairPotential::setShortRangeTruncationScheme(PairPotential::NoShortRangeTruncation); @@ -61,39 +61,39 @@ class PairPotentialsTest : public ::testing::Test void testEnergy(Functions1D::Form form, std::string_view parameters, double rStart = testRDelta_) { test( - form, parameters, rStart, [=](double r) { return pairPotential_->analyticEnergy(r); }, + form, parameters, rStart, [=](double r) { return pairPotential_->analyticEnergy(r, 1.0, 1.0); }, [=](double r) { return pairPotential_->energy(r); }); } // Test analytic vs tabulated force for specified form and parameters void testForce(Functions1D::Form form, std::string_view parameters, double rStart = testRDelta_) { test( - form, parameters, rStart, [=](double r) { return pairPotential_->analyticForce(r); }, + form, parameters, rStart, [=](double r) { return pairPotential_->analyticForce(r, 1.0, 1.0); }, [=](double r) { return pairPotential_->force(r); }); } }; -TEST_F(PairPotentialsTest, NoneForm) +TEST_F(PairPotentialsAnalyticTest, NoneForm) { testEnergy(Functions1D::Form::None, ""); testForce(Functions1D::Form::None, ""); } -TEST_F(PairPotentialsTest, LennardJonesForm) +TEST_F(PairPotentialsAnalyticTest, LennardJonesForm) { // Lennard-Jones is super steep at r -> 0 so start a little after that testEnergy(Functions1D::Form::LennardJones126, "epsilon=0.35 sigma=2.166", testRDelta_ * 5); testForce(Functions1D::Form::LennardJones126, "epsilon=0.35 sigma=2.166", testRDelta_ * 5); } -TEST_F(PairPotentialsTest, Buckingham) +TEST_F(PairPotentialsAnalyticTest, Buckingham) { // Values put in for TeO2 from https://pubs.rsc.org/en/content/articlelanding/2014/cp/c4cp01273a testEnergy(Functions1D::Form::Buckingham, "A=153919.5503 B=2.8912848 C=96.48514925", testRDelta_ * 5); testForce(Functions1D::Form::Buckingham, "A=8005439.257 B=6.211565 C=3025.962812", testRDelta_ * 5); } -TEST_F(PairPotentialsTest, Gaussian) +TEST_F(PairPotentialsAnalyticTest, Gaussian) { testEnergy(Functions1D::Form::GaussianPotential, "A=1.0 fwhm=2.4 x0=0", testRDelta_ * 5); testForce(Functions1D::Form::GaussianPotential, "A=1.0 fwhm=2.4 x0=0", testRDelta_ * 5); diff --git a/tests/pairPotentials/scaleFactors.cpp b/tests/pairPotentials/scaleFactors.cpp new file mode 100644 index 0000000000..55890e386c --- /dev/null +++ b/tests/pairPotentials/scaleFactors.cpp @@ -0,0 +1,196 @@ +// SPDX-License-Identifier: GPL-3.0-or-later +// Copyright (c) 2024 Team Dissolve and contributors + +#include "classes/coreData.h" +#include "classes/localMolecule.h" +#include "classes/potentialMap.h" +#include "classes/species.h" +#include + +namespace UnitTest +{ +class PairPotentialsScaleFactorsTest : public ::testing::Test +{ + public: + PairPotentialsScaleFactorsTest() + { + // Set up atom types + atC1_ = coreData_.addAtomType(Elements::C); + atC1_->setName("C1"); + atC1_->setCharge(-0.1); + atC1_->interactionPotential().setFormAndParameters(ShortRangeFunctions::Form::LennardJones, ljParameters); + atC2_ = coreData_.addAtomType(Elements::C); + atC2_->setCharge(0.1); + atC2_->setName("C2"); + atC2_->interactionPotential().setFormAndParameters(ShortRangeFunctions::Form::LennardJones, ljParameters); + + // Create Species + fourSixthsBenzene_.addAtom(Elements::C, {-1.390000, 0.000000, 0.000000}, -0.2, atC1_); + fourSixthsBenzene_.addAtom(Elements::C, {-0.695000, 1.203775, 0.000000}, 0.2, atC2_); + fourSixthsBenzene_.addAtom(Elements::C, {0.695000, 1.203775, 0.000000}, 0.2, atC2_); + fourSixthsBenzene_.addAtom(Elements::C, {1.390000, 0.000000, 0.000000}, -0.2, atC1_); + fourSixthsBenzene_.addMissingBonds(); + torsion_ = fourSixthsBenzene_.addTorsion(0, 1, 2, 3); + fourSixthsBenzene_.setUpScaledInteractions(); + + // Create a molecule based on the species + molecule_.setSpecies(&fourSixthsBenzene_); + for (auto &&[molAtom, spAtom] : zip(molecule_.localAtoms(), fourSixthsBenzene_.atoms())) + { + molAtom.setCoordinates(spAtom.r()); + molAtom.setMasterTypeIndex(spAtom.atomType()->index()); + } + + // Set up the function wrapper + potentialWrapper_.setFormAndParameters(interactionPotential_.form(), interactionPotential_.parameters()); + } + // Set up potentials + void setUpPotentials(bool useAtomTypeCharges) + { + const auto ppRange = 15.0, ppDelta = 0.001; + + potentialMap_.clear(); + pairPotentials_.clear(); + PairPotential::setIncludeCoulombPotential(useAtomTypeCharges); + PairPotential::setShortRangeTruncationScheme(PairPotential::NoShortRangeTruncation); + PairPotential::setCoulombTruncationScheme(PairPotential::NoCoulombTruncation); + + // Set up pair potentials + auto *pp11 = std::get<2>(pairPotentials_.emplace_back(atC1_, atC1_, + std::make_unique(atC1_->name(), atC1_->name()))) + .get(); + pp11->setInteractionPotential(interactionPotential_); + pp11->tabulate(ppRange, ppDelta, atC1_->charge() * atC1_->charge()); + auto *pp12 = std::get<2>(pairPotentials_.emplace_back(atC1_, atC2_, + std::make_unique(atC1_->name(), atC2_->name()))) + .get(); + pp12->setInteractionPotential(interactionPotential_); + pp12->tabulate(ppRange, ppDelta, atC1_->charge() * atC2_->charge()); + auto *pp22 = std::get<2>(pairPotentials_.emplace_back(atC2_, atC2_, + std::make_unique(atC2_->name(), atC2_->name()))) + .get(); + pp22->setInteractionPotential(interactionPotential_); + pp22->tabulate(ppRange, ppDelta, atC2_->charge() * atC2_->charge()); + + // Create PotentialMap + potentialMap_.initialise(coreData_.atomTypes(), pairPotentials_, ppRange); + } + // Return reference energy at distance r given specified charge product and scalings + double referenceEnergy(double r, double chargeProduct, double elecScale = 1.0, double srScale = 1.0) + { + return srScale * potentialWrapper_.y(r) + elecScale * COULCONVERT * chargeProduct / r; + } + // Return reference force at distance r given specified charge product and scalings + double referenceForce(double r, double chargeProduct, double elecScale = 1.0, double srScale = 1.0) + { + return -srScale * potentialWrapper_.dYdX(r) + elecScale * COULCONVERT * chargeProduct / (r * r); + } + // Perform scaling tests on production routines + template void testScalings(const Particle &i, const Particle &j, double r, double refChargeProduct) + { + // No scaling + EXPECT_NEAR(potentialMap_.energy(i, j, r), referenceEnergy(r, refChargeProduct), testTolerance_); + EXPECT_NEAR(potentialMap_.force(i, j, r), referenceForce(r, refChargeProduct), testTolerance_); + + auto elecScale = 0.5, srScale = 0.5; + + // Equal scaling for short-range and electrostatics + EXPECT_NEAR(potentialMap_.energy(i, j, r, elecScale, srScale), referenceEnergy(r, refChargeProduct, elecScale, srScale), + testTolerance_); + EXPECT_NEAR(potentialMap_.force(i, j, r, elecScale, srScale), referenceForce(r, refChargeProduct, elecScale, srScale), + testTolerance_); + + // Unlike scalings + elecScale = 0.25; + srScale = 0.75; + EXPECT_NEAR(potentialMap_.energy(i, j, r, elecScale, srScale), referenceEnergy(r, refChargeProduct, elecScale, srScale), + testTolerance_); + EXPECT_NEAR(potentialMap_.force(i, j, r, elecScale, srScale), referenceForce(r, refChargeProduct, elecScale, srScale), + testTolerance_); + } + // Perform scaling tests on analytic routines + void testAnalyticScalings(const Atom &i, const Atom &j, double r, double refChargeProduct) + { + // No scaling + EXPECT_NEAR(potentialMap_.analyticEnergy(i, j, r), referenceEnergy(r, refChargeProduct), testTolerance_); + EXPECT_NEAR(potentialMap_.analyticForce(i, j, r), referenceForce(r, refChargeProduct), testTolerance_); + auto elecScale = 0.5, srScale = 0.5; + + // Equal scaling for short-range and electrostatics + EXPECT_NEAR(potentialMap_.analyticEnergy(i, j, r, elecScale, srScale), + referenceEnergy(r, refChargeProduct, elecScale, srScale), testTolerance_); + EXPECT_NEAR(potentialMap_.analyticForce(i, j, r, elecScale, srScale), + referenceForce(r, refChargeProduct, elecScale, srScale), testTolerance_); + + // Unlike scalings + elecScale = 0.25; + srScale = 0.75; + EXPECT_NEAR(potentialMap_.analyticEnergy(i, j, r, elecScale, srScale), + referenceEnergy(r, refChargeProduct, elecScale, srScale), testTolerance_); + EXPECT_NEAR(potentialMap_.analyticForce(i, j, r, elecScale, srScale), + referenceForce(r, refChargeProduct, elecScale, srScale), testTolerance_); + } + + protected: + // Lennard-Jones parameters for test + static constexpr auto ljParameters = "epsilon=1.0 sigma=3.0"; + // Double value test tolerance + static constexpr auto testTolerance_ = 1.0e-8; + + CoreData coreData_; + std::shared_ptr atC1_, atC2_; + std::vector pairPotentials_; + PotentialMap potentialMap_; + Species fourSixthsBenzene_; + SpeciesTorsion torsion_; + LocalMolecule molecule_; + + // Test potential function and wrapper, equivalent to the one defined in the potential map + InteractionPotential interactionPotential_{Functions1D::Form::LennardJones126, ljParameters}; + Function1DWrapper potentialWrapper_; +}; + +TEST_F(PairPotentialsScaleFactorsTest, SpeciesEnergyWithAtomTypeCharges) +{ + setUpPotentials(true); + + auto &i = fourSixthsBenzene_.atom(0); + auto &j = fourSixthsBenzene_.atom(3); + + testScalings(&i, &j, (j.r() - i.r()).magnitude(), atC1_->charge() * atC1_->charge()); +} + +TEST_F(PairPotentialsScaleFactorsTest, SpeciesEnergyWithSpeciesCharges) +{ + setUpPotentials(false); + + auto &i = fourSixthsBenzene_.atom(0); + auto &j = fourSixthsBenzene_.atom(3); + + testScalings(&i, &j, (j.r() - i.r()).magnitude(), i.charge() * j.charge()); +} + +TEST_F(PairPotentialsScaleFactorsTest, MoleculeEnergyWithAtomTypeCharges) +{ + setUpPotentials(true); + + auto &i = molecule_.localAtoms()[0]; + auto &j = molecule_.localAtoms()[3]; + + testScalings(i, j, (j.r() - i.r()).magnitude(), atC1_->charge() * atC1_->charge()); + testAnalyticScalings(i, j, (j.r() - i.r()).magnitude(), atC1_->charge() * atC1_->charge()); +} + +TEST_F(PairPotentialsScaleFactorsTest, MoleculeEnergyWithSpeciesCharges) +{ + setUpPotentials(false); + + auto &i = molecule_.localAtoms()[0]; + auto &j = molecule_.localAtoms()[3]; + + testScalings(i, j, (j.r() - i.r()).magnitude(), fourSixthsBenzene_.atom(0).charge() * fourSixthsBenzene_.atom(3).charge()); + testAnalyticScalings(i, j, (j.r() - i.r()).magnitude(), + fourSixthsBenzene_.atom(0).charge() * fourSixthsBenzene_.atom(3).charge()); +} + +} // namespace UnitTest