refactor: Forcefield tidy up

src/classes/species.cpp

@@ -172,45 +172,6 @@ void Species::print() const
     }
 }
 
-/*
- * Forcefield
- */
-
-// Apply terms from specified Forcefield
-bool Species::applyForcefieldTerms(std::shared_ptr<Forcefield> ff)
-{
-    if (!ff)
-        return Messenger::error("No forcefield supplied!\n");
-
-    // Assign atom types to the species
-    if (!ff->assignAtomTypes(this, Forcefield::TypeAll, false).empty())
-        return false;
-
-    // Assign intramolecular terms
-    if (!ff->assignIntramolecular(this))
-        return false;
-
-    return true;
-}
-
-// Clear forcefield terms
-void Species::clearForcefieldTerms(bool nullifyAtomTypes)
-{
-    if (nullifyAtomTypes)
-        clearAtomTypes();
-
-    for (auto &b : bonds_)
-        b.setInteractionFormAndParameters(BondFunctions::Form::None, std::vector<double>());
-
-    for (auto &a : angles_)
-        a.setInteractionFormAndParameters(AngleFunctions::Form::None, std::vector<double>());
-
-    for (auto &t : torsions_)
-        t.setInteractionFormAndParameters(TorsionFunctions::Form::None, std::vector<double>());
-
-    impropers_.clear();
-}
-
 /*
  * Serialisation
  */

src/classes/species.h

@@ -107,40 +107,37 @@ class Species : public Serialisable<>
     // Return vector of SpeciesBond
     std::vector<SpeciesBond> &bonds();
     const std::vector<SpeciesBond> &bonds() const;
-    // Return whether SpeciesBond between SpeciesAtoms exists
-    bool hasBond(const SpeciesAtom *i, const SpeciesAtom *j) const;
-    bool hasBond(int i, int j) const;
-    // Return the SpeciesBond between the specified SpeciesAtom indices
-    OptionalReferenceWrapper<SpeciesBond> getBond(int indexI, int indexJ);
+    // Return the SpeciesBond between the specified SpeciesAtom indices, if it exists
+    OptionalReferenceWrapper<SpeciesBond> getBond(const SpeciesAtom *i, const SpeciesAtom *j);
     // Remove bonds crossing periodic boundaries
     void removePeriodicBonds();
     // Return vector of SpeciesAngle
     std::vector<SpeciesAngle> &angles();
     const std::vector<SpeciesAngle> &angles() const;
-    // Return whether SpeciesAngle between SpeciesAtoms exists
-    bool hasAngle(SpeciesAtom *i, SpeciesAtom *j, SpeciesAtom *k) const;
-    // Return the SpeciesAngle between the specified SpeciesAtom indices
-    OptionalReferenceWrapper<SpeciesAngle> getAngle(int indexI, int indexJ, int indexK);
+    // Return the SpeciesAngle between the specified SpeciesAtom indices, if it exists
+    OptionalReferenceWrapper<SpeciesAngle> getAngle(const SpeciesAtom *i, const SpeciesAtom *j, const SpeciesAtom *k);
     // Return vector of SpeciesTorsion
     std::vector<SpeciesTorsion> &torsions();
     const std::vector<SpeciesTorsion> &torsions() const;
-    // Return whether SpeciesTorsion between SpeciesAtoms exists
-    bool hasTorsion(SpeciesAtom *i, SpeciesAtom *j, SpeciesAtom *k, SpeciesAtom *l) const;
-    // Return the SpeciesTorsion between the specified SpeciesAtom indices
-    OptionalReferenceWrapper<SpeciesTorsion> getTorsion(int indexI, int indexJ, int indexK, int indexL);
+    // Return the SpeciesTorsion between the specified SpeciesAtom indices, if it exists
+    OptionalReferenceWrapper<SpeciesTorsion> getTorsion(const SpeciesAtom *i, const SpeciesAtom *j, const SpeciesAtom *k,
+                                                        const SpeciesAtom *l);
     // Return vector of SpeciesImproper
     std::vector<SpeciesImproper> &impropers();
     const std::vector<SpeciesImproper> &impropers() const;
-    // Return whether SpeciesImproper between SpeciesAtoms exists
-    bool hasImproper(SpeciesAtom *i, SpeciesAtom *j, SpeciesAtom *k, SpeciesAtom *l) const;
-    // Return the SpeciesImproper between the specified SpeciesAtom indices
-    OptionalReferenceWrapper<SpeciesImproper> getImproper(int indexI, int indexJ, int indexK, int indexL);
+    // Add a new improper term between the specified atoms
+    SpeciesImproper &addImproper(SpeciesAtom *i, SpeciesAtom *j, SpeciesAtom *k, SpeciesAtom *l);
+    // Return the SpeciesImproper between the specified SpeciesAtom indices, if it exists
+    OptionalReferenceWrapper<SpeciesImproper> getImproper(const SpeciesAtom *i, const SpeciesAtom *j, const SpeciesAtom *k,
+                                                          const SpeciesAtom *l);
     // Return whether the attached atoms lists have been created
     bool attachedAtomListsGenerated() const;
     // Determine angles and torsions from bond connectivity
     void determineAnglesAndTorsions();
     // Finalise internal relationships related to geometry once it is defined
     void finaliseGeometry();
+    // Clear forcefield data from intramolecular terms
+    void clearIntramolecularForcefieldTerms();
 
     /*
      * Intramolecular Common Terms
@@ -218,15 +215,6 @@ class Species : public Serialisable<>
     // Create Box definition with specified lengths and angles
     void createBox(const Vector3 lengths, const Vector3 angles, bool nonPeriodic = false);
 
-    /*
-     * Forcefield
-     */
-    public:
-    // Apply terms from specified Forcefield
-    bool applyForcefieldTerms(std::shared_ptr<Forcefield> ff);
-    // Clear forcefield terms
-    void clearForcefieldTerms(bool nullifyAtomTypes = true);
-
     /*
      * Isotopologues
      */

src/classes/speciesIntra.h

@@ -34,12 +34,16 @@ template <class Intra, class Functions> class SpeciesIntra : public Serialisable
     SpeciesIntra &operator=(SpeciesIntra &&source) = delete;
 
     /*
-     * SpeciesAtom Information
+     * Parent Species
      */
     protected:
     // Parent species
     Species *parent_{nullptr};
 
+    public:
+    // Return parent species
+    Species *parent() const { return parent_; }
+
     /*
      * SpeciesAtom Information
      */

src/classes/species_intra.cpp

@@ -15,19 +15,9 @@ std::vector<SpeciesBond> &Species::bonds() { return bonds_; }
 
 const std::vector<SpeciesBond> &Species::bonds() const { return bonds_; }
 
-// Return whether SpeciesBond between specified SpeciesAtoms exists
-bool Species::hasBond(const SpeciesAtom *i, const SpeciesAtom *j) const
+// Return the SpeciesBond between the specified SpeciesAtoms, if it exists
+OptionalReferenceWrapper<SpeciesBond> Species::getBond(const SpeciesAtom *i, const SpeciesAtom *j)
 {
-    return std::any_of(bonds_.cbegin(), bonds_.cend(), [i, j](const auto &bond) { return bond.matches(i, j); });
-}
-bool Species::hasBond(int i, int j) const { return hasBond(&atom(i), &atom(j)); }
-
-// Return the SpeciesBond between the specified SpeciesAtom indices
-OptionalReferenceWrapper<SpeciesBond> Species::getBond(int indexI, int indexJ)
-{
-    auto *i = &atoms_[indexI];
-    auto *j = &atoms_[indexJ];
-
     auto it = std::find_if(bonds_.begin(), bonds_.end(), [i, j](auto &bond) { return bond.matches(i, j); });
     if (it == bonds_.end())
         return {};
@@ -57,19 +47,9 @@ std::vector<SpeciesAngle> &Species::angles() { return angles_; }
 
 const std::vector<SpeciesAngle> &Species::angles() const { return angles_; }
 
-// Return whether SpeciesAngle between SpeciesAtoms exists
-bool Species::hasAngle(SpeciesAtom *i, SpeciesAtom *j, SpeciesAtom *k) const
-{
-    return std::any_of(angles_.cbegin(), angles_.cend(), [i, j, k](const auto &angle) { return angle.matches(i, j, k); });
-}
-
-// Return the SpeciesAngle between the specified SpeciesAtom indic
-OptionalReferenceWrapper<SpeciesAngle> Species::getAngle(int indexI, int indexJ, int indexK)
+// Return the SpeciesAngle between the specified SpeciesAtoms, if it exists
+OptionalReferenceWrapper<SpeciesAngle> Species::getAngle(const SpeciesAtom *i, const SpeciesAtom *j, const SpeciesAtom *k)
 {
-    auto *i = &atoms_[indexI];
-    auto *j = &atoms_[indexJ];
-    auto *k = &atoms_[indexK];
-
     auto it = std::find_if(angles_.begin(), angles_.end(), [i, j, k](auto &angle) { return angle.matches(i, j, k); });
     if (it == angles_.end())
         return {};
@@ -82,20 +62,10 @@ std::vector<SpeciesTorsion> &Species::torsions() { return torsions_; }
 
 const std::vector<SpeciesTorsion> &Species::torsions() const { return torsions_; }
 
-// Return whether SpeciesTorsion between SpeciesAtoms exists
-bool Species::hasTorsion(SpeciesAtom *i, SpeciesAtom *j, SpeciesAtom *k, SpeciesAtom *l) const
-{
-    return std::any_of(torsions_.cbegin(), torsions_.cend(), [&](const auto &t) { return t.matches(i, j, k, l); });
-}
-
-// Return the SpeciesTorsion between the specified SpeciesAtom indices
-OptionalReferenceWrapper<SpeciesTorsion> Species::getTorsion(int indexI, int indexJ, int indexK, int indexL)
+// Return the SpeciesTorsion between the specified SpeciesAtoms, if it exists
+OptionalReferenceWrapper<SpeciesTorsion> Species::getTorsion(const SpeciesAtom *i, const SpeciesAtom *j, const SpeciesAtom *k,
+                                                             const SpeciesAtom *l)
 {
-    auto *i = &atoms_[indexI];
-    auto *j = &atoms_[indexJ];
-    auto *k = &atoms_[indexK];
-    auto *l = &atoms_[indexL];
-
     auto it =
         std::find_if(torsions_.begin(), torsions_.end(), [i, j, k, l](auto &torsion) { return torsion.matches(i, j, k, l); });
     if (it == torsions_.end())
@@ -109,21 +79,20 @@ std::vector<SpeciesImproper> &Species::impropers() { return impropers_; }
 
 const std::vector<SpeciesImproper> &Species::impropers() const { return impropers_; }
 
-// Return whether SpeciesImproper between SpeciesAtoms exists
-bool Species::hasImproper(SpeciesAtom *i, SpeciesAtom *j, SpeciesAtom *k, SpeciesAtom *l) const
+// Add a new improper term between the specified atoms
+SpeciesImproper &Species::addImproper(SpeciesAtom *i, SpeciesAtom *j, SpeciesAtom *k, SpeciesAtom *l)
 {
-    return std::any_of(impropers_.cbegin(), impropers_.cend(),
-                       [i, j, k, l](auto &improper) { return improper.matches(i, j, k, l); });
+    auto optImproper = getImproper(i, j, k, l);
+    if (optImproper)
+        return *optImproper;
+
+    return impropers_.emplace_back(this, i, j, k, l);
 }
 
-// Return the SpeciesImproper between the specified SpeciesAtom indices
-OptionalReferenceWrapper<SpeciesImproper> Species::getImproper(int indexI, int indexJ, int indexK, int indexL)
+// Return the SpeciesImproper between the specified SpeciesAtoms, if it exists
+OptionalReferenceWrapper<SpeciesImproper> Species::getImproper(const SpeciesAtom *i, const SpeciesAtom *j, const SpeciesAtom *k,
+                                                               const SpeciesAtom *l)
 {
-    auto *i = &atoms_[indexI];
-    auto *j = &atoms_[indexJ];
-    auto *k = &atoms_[indexK];
-    auto *l = &atoms_[indexL];
-
     auto it = std::find_if(impropers_.begin(), impropers_.end(),
                            [i, j, k, l](auto &improper) { return improper.matches(i, j, k, l); });
     if (it == impropers_.end())
@@ -332,6 +301,21 @@ void Species::finaliseGeometry()
     attachedAtomListsGenerated_ = true;
 }
 
+// Clear intramolecular forcefield terms
+void Species::clearIntramolecularForcefieldTerms()
+{
+    for (auto &b : bonds_)
+        b.setInteractionFormAndParameters(BondFunctions::Form::None, std::vector<double>());
+
+    for (auto &a : angles_)
+        a.setInteractionFormAndParameters(AngleFunctions::Form::None, std::vector<double>());
+
+    for (auto &t : torsions_)
+        t.setInteractionFormAndParameters(TorsionFunctions::Form::None, std::vector<double>());
+
+    impropers_.clear();
+}
+
 // Return periodic box
 const Box *Species::box() const { return box_.get(); }
 

src/classes/species_io.cpp

@@ -95,7 +95,7 @@ bool Species::read(LineParser &parser, CoreData &coreData)
                 // Create a new angle definition between the specified atoms
                 if (angleVectorFixed && angleIndex < angles_.size())
                 {
-                    if (hasAngle(&atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1], &atoms_[parser.argi(3) - 1]))
+                    if (getAngle(&atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1], &atoms_[parser.argi(3) - 1]))
                         return false;
                     angles_[angleIndex] = SpeciesAngle(this, &atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1],
                                                        &atoms_[parser.argi(3) - 1]);
@@ -194,7 +194,7 @@ bool Species::read(LineParser &parser, CoreData &coreData)
                 // Create a new bond definition between the specified atoms
                 if (bondVectorFixed && bondIndex < bonds_.size())
                 {
-                    if (hasBond(&atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1]))
+                    if (getBond(&atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1]))
                         return false;
                     bonds_[bondIndex] = SpeciesBond(this, &atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1]);
                     b = bonds_[bondIndex++];
@@ -204,7 +204,7 @@ bool Species::read(LineParser &parser, CoreData &coreData)
                     auto i = parser.argi(1) - 1;
                     auto j = parser.argi(2) - 1;
                     // Check for existence of Bond already
-                    auto bondRef = getBond(i, j);
+                    auto bondRef = getBond(&atoms_[i], &atoms_[j]);
                     if (bondRef)
                     {
                         Messenger::warn("Refused to add a new SpeciesBond between atoms {} and {} in Species '{}' since it "
@@ -266,7 +266,7 @@ bool Species::read(LineParser &parser, CoreData &coreData)
                 break;
             case (Species::SpeciesKeyword::BondType):
                 // Find the specified bond
-                b = getBond(parser.argi(1) - 1, parser.argi(2) - 1);
+                b = getBond(&atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1]);
                 if (!b)
                 {
                     Messenger::error("Tried to set the bond type of bond between atoms {} and {}, but this bond "
@@ -316,7 +316,7 @@ bool Species::read(LineParser &parser, CoreData &coreData)
                 // Create a new improper definition
                 if (improperVectorFixed && improperIndex < impropers_.size())
                 {
-                    if (hasImproper(&atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1], &atoms_[parser.argi(3) - 1],
+                    if (getImproper(&atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1], &atoms_[parser.argi(3) - 1],
                                     &atoms_[parser.argi(4) - 1]))
                         return false;
                     impropers_[improperIndex] = SpeciesImproper(this, &atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1],
@@ -481,7 +481,7 @@ bool Species::read(LineParser &parser, CoreData &coreData)
                 // Create a new angle definition between the specified atoms
                 if (torsionVectorFixed && torsionIndex < torsions_.size())
                 {
-                    if (hasTorsion(&atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1], &atoms_[parser.argi(3) - 1],
+                    if (getTorsion(&atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1], &atoms_[parser.argi(3) - 1],
                                    &atoms_[parser.argi(4) - 1]))
                         return false;
                     torsions_[torsionIndex] = SpeciesTorsion(this, &atoms_[parser.argi(1) - 1], &atoms_[parser.argi(2) - 1],

src/data/ff/angleTerm.cpp

@@ -5,6 +5,12 @@
 #include "data/ff/atomType.h"
 #include "data/ff/ff.h"
 
+ForcefieldAngleTerm::ForcefieldAngleTerm(std::string_view typeI, std::string_view typeJ, std::string_view typeK,
+                                         AngleFunctions::Form form, const std::vector<double> &parameters)
+    : typeI_(typeI), typeJ_(typeJ), typeK_(typeK), form_(form), parameters_(parameters)
+{
+}
+
 ForcefieldAngleTerm::ForcefieldAngleTerm(std::string_view typeI, std::string_view typeJ, std::string_view typeK,
                                          AngleFunctions::Form form, std::string_view parameterString)
     : typeI_(typeI), typeJ_(typeJ), typeK_(typeK), form_(form)

src/data/ff/angleTerm.h

@@ -6,15 +6,17 @@
 #include "classes/speciesAngle.h"
 
 // Forward Declarations
-class Forcefield;
 class ForcefieldAtomType;
 
 // Forcefield Angle Term
 class ForcefieldAngleTerm
 {
     public:
-    ForcefieldAngleTerm(std::string_view typeI = "", std::string_view typeJ = "", std::string_view typeK = "",
-                        AngleFunctions::Form form = AngleFunctions::Form::None, std::string_view parameterString = "");
+    ForcefieldAngleTerm() = default;
+    ForcefieldAngleTerm(std::string_view typeI, std::string_view typeJ, std::string_view typeK, AngleFunctions::Form form,
+                        const std::vector<double> &parameters);
+    ForcefieldAngleTerm(std::string_view typeI, std::string_view typeJ, std::string_view typeK, AngleFunctions::Form form,
+                        std::string_view parameterString);
     ~ForcefieldAngleTerm() = default;
 
     /*

src/data/ff/bondTerm.cpp

@@ -5,6 +5,12 @@
 #include "data/ff/atomType.h"
 #include "data/ff/ff.h"
 
+ForcefieldBondTerm::ForcefieldBondTerm(std::string_view typeI, std::string_view typeJ, BondFunctions::Form form,
+                                       const std::vector<double> &parameters)
+    : typeI_(typeI), typeJ_(typeJ), form_(form), parameters_(parameters)
+{
+}
+
 ForcefieldBondTerm::ForcefieldBondTerm(std::string_view typeI, std::string_view typeJ, BondFunctions::Form form,
                                        std::string_view parameterString)
     : typeI_(typeI), typeJ_(typeJ), form_(form)

src/data/ff/bondTerm.h

@@ -6,15 +6,17 @@
 #include "classes/speciesBond.h"
 
 // Forward Declarations
-class Forcefield;
 class ForcefieldAtomType;
 
 // Forcefield Bond Term
 class ForcefieldBondTerm
 {
     public:
-    ForcefieldBondTerm(std::string_view typeI = "", std::string_view typeJ = "",
-                       BondFunctions::Form form = BondFunctions::Form::None, std::string_view parameterString = "");
+    ForcefieldBondTerm() = default;
+    ForcefieldBondTerm(std::string_view typeI, std::string_view typeJ, BondFunctions::Form form,
+                       const std::vector<double> &parameters);
+    ForcefieldBondTerm(std::string_view typeI, std::string_view typeJ, BondFunctions::Form form,
+                       std::string_view parameterString);
     ~ForcefieldBondTerm() = default;
 
     /*