feat: Add on sphere (#2011)

disorderedmaterials · Nov 19, 2024 · 84bd77c · 84bd77c
1 parent 887387c
commit 84bd77c
Show file tree

Hide file tree

Showing 20 changed files with 543 additions and 80 deletions.
diff --git a/src/classes/molecule.cpp b/src/classes/molecule.cpp
@@ -6,6 +6,7 @@
 #include "classes/box.h"
 #include "classes/species.h"
 #include "classes/speciesBond.h"
+#include "data/atomicMasses.h"
 #include "templates/algorithms.h"
 
 /*
@@ -143,7 +144,7 @@ Vec3<double> Molecule::unFold(const Box *box)
 }
 
 // Set centre of geometry of molecule
-void Molecule::setCentreOfGeometry(const Box *box, const Vec3<double> newCentre)
+void Molecule::setCentreOfGeometry(const Box *box, const Vec3<double> &newCentre)
 {
     // Calculate Molecule centre of geometry
     Vec3<double> newR;
@@ -169,6 +170,30 @@ Vec3<double> Molecule::centreOfGeometry(const Box *box) const
     return cog / nAtoms();
 }
 
+// Calculate and return centre of geometry over supplied atom indices
+Vec3<double> Molecule::centreOfGeometry(const Box *box, const std::vector<int> &indices) const
+{
+    const auto ref = atoms_[indices.front()]->r();
+    return std::accumulate(std::next(indices.begin()), indices.end(), ref,
+                           [&](const auto &acc, const auto idx) { return acc + box->minimumImage(atoms_[idx]->r(), ref); }) /
+           indices.size();
+}
+
+// Calculate and return centre of mass over supplied atom indices
+Vec3<double> Molecule::centreOfMass(const Box *box, const std::vector<int> &indices) const
+{
+    auto mass = AtomicMass::mass(atoms_[indices.front()]->speciesAtom()->Z());
+    const auto ref = atoms_[indices.front()]->r();
+    auto sums = std::accumulate(std::next(indices.begin()), indices.end(), std::pair<Vec3<double>, double>(ref * mass, mass),
+                                [&](const auto &acc, const auto idx)
+                                {
+                                    auto mass = AtomicMass::mass(atoms_[idx]->speciesAtom()->Z());
+                                    return std::pair<Vec3<double>, double>(
+                                        acc.first + box->minimumImage(atoms_[idx]->r(), ref) * mass, acc.second + mass);
+                                });
+    return sums.first / sums.second;
+}
+
 // Transform molecule with supplied matrix, using centre of geometry as the origin
 void Molecule::transform(const Box *box, const Matrix3 &transformationMatrix)
 {
@@ -207,7 +232,7 @@ void Molecule::transform(const Box *box, const Matrix3 &transformationMatrix, co
 }
 
 // Translate whole molecule by the delta specified
-void Molecule::translate(const Vec3<double> delta)
+void Molecule::translate(const Vec3<double> &delta)
 {
     for (auto n = 0; n < nAtoms(); ++n)
         atom(n)->translateCoordinates(delta);

diff --git a/src/classes/molecule.h b/src/classes/molecule.h
@@ -73,9 +73,13 @@ class Molecule : public std::enable_shared_from_this<Molecule>
     // Un-fold molecule so it is not cut by box boundaries, returning the centre of geometry
     Vec3<double> unFold(const Box *box);
     // Set centre of geometry
-    void setCentreOfGeometry(const Box *box, const Vec3<double> newCentre);
+    void setCentreOfGeometry(const Box *box, const Vec3<double> &newCentre);
     // Calculate and return centre of geometry
     Vec3<double> centreOfGeometry(const Box *box) const;
+    // Calculate and return centre of geometry over supplied atom indices
+    Vec3<double> centreOfGeometry(const Box *box, const std::vector<int> &indices) const;
+    // Calculate and return centre of mass over supplied atom indices
+    Vec3<double> centreOfMass(const Box *box, const std::vector<int> &indices) const;
     // Transform molecule with supplied matrix, using centre of geometry as the origin
     void transform(const Box *box, const Matrix3 &transformationMatrix);
     // Transform molecule with supplied matrix about specified origin
@@ -84,7 +88,7 @@ class Molecule : public std::enable_shared_from_this<Molecule>
     void transform(const Box *box, const Matrix3 &transformationMatrix, const Vec3<double> &origin,
                    const std::vector<int> &targetAtoms);
     // Translate whole molecule by the delta specified
-    void translate(const Vec3<double> delta);
+    void translate(const Vec3<double> &delta);
     // Translate specified atoms by the delta specified
     void translate(const Vec3<double> &delta, const std::vector<int> &targetAtoms);
 };
diff --git a/src/classes/site.cpp b/src/classes/site.cpp
@@ -2,8 +2,11 @@
 // Copyright (c) 2024 Team Dissolve and contributors
 
 #include "classes/site.h"
-
 #include "base/messenger.h"
+#include "classes/box.h"
+#include "classes/molecule.h"
+#include "classes/speciesSite.h"
+#include "classes/speciesSiteInstance.h"
 #include <utility>
 
 Site::Site(const SpeciesSite *parent, std::optional<int> uniqueSiteIndex, std::shared_ptr<const Molecule> molecule,
@@ -19,6 +22,34 @@ Site::Site(const SpeciesSite *parent, std::optional<int> uniqueSiteIndex, std::s
 {
 }
 
+Site::Site(const SpeciesSite *parent, std::optional<int> uniqueSiteIndex, std::shared_ptr<const Molecule> molecule,
+           const SpeciesSiteInstance &instance, const Box *box)
+    : parent_(parent), molecule_(std::move(molecule))
+{
+    origin_ = parent_->originMassWeighted() ? molecule_->centreOfMass(box, instance.originIndices())
+                                            : molecule_->centreOfGeometry(box, instance.originIndices());
+
+    if (parent_->hasAxes())
+    {
+        // Get vector from site origin to x-axis reference point and normalise it
+        auto x = box->minimumVector(origin_, molecule_->centreOfGeometry(box, instance.xAxisIndices()));
+        x.normalise();
+
+        axes_.setColumn(0, x);
+
+        // Get vector from site origin to y-axis reference point, normalise it, and orthogonalise
+        auto y = box->minimumVector(origin_, molecule_->centreOfGeometry(box, instance.yAxisIndices()));
+        y.orthogonalise(x);
+        y.normalise();
+
+        axes_.setColumn(1, y);
+
+        axes_.setColumn(2, x * y);
+
+        hasAxes_ = true;
+    }
+}
+
 // Return enum options for SiteAxis
 EnumOptions<Site::SiteAxis> Site::siteAxis()
 {

diff --git a/src/classes/site.h b/src/classes/site.h
@@ -9,8 +9,10 @@
 #include <optional>
 
 // Forward Declarations
+class Box;
 class Molecule;
 class SpeciesSite;
+class SpeciesSiteInstance;
 
 // Site Definition
 class Site
@@ -20,6 +22,8 @@ class Site
          std::shared_ptr<const Molecule> molecule = {}, const Vec3<double> &origin = {});
     Site(const SpeciesSite *parent = nullptr, std::optional<int> uniqueSiteIndex = {},
          std::shared_ptr<const Molecule> molecule = {}, const Matrix3 &axes = {}, const Vec3<double> &origin = {});
+    Site(const SpeciesSite *parent, std::optional<int> uniqueSiteIndex, std::shared_ptr<const Molecule> molecule,
+         const SpeciesSiteInstance &instance, const Box *box);
     ~Site() = default;
     Site &operator=(const Site &source) = default;
     Site(const Site &source) = default;
@@ -66,6 +70,4 @@ class Site
     bool hasAxes() const;
     // Return local axes
     const Matrix3 &axes() const;
-    // Rotate about axis
-    void rotate(double angle, Site::SiteAxis axis);
 };
diff --git a/src/classes/siteStack.cpp b/src/classes/siteStack.cpp
@@ -6,7 +6,6 @@
 #include "classes/configuration.h"
 #include "classes/species.h"
 #include "classes/speciesSite.h"
-#include "data/atomicMasses.h"
 #include <algorithm>
 #include <numeric>
 
@@ -33,31 +32,6 @@ const SpeciesSite *SiteStack::speciesSite() const { return speciesSite_; }
  * Generation
  */
 
-// Calculate geometric centre of atoms in the given molecule
-Vec3<double> SiteStack::centreOfGeometry(const Molecule &mol, const Box *box, const std::vector<int> &indices)
-{
-    const auto ref = mol.atom(indices.front())->r();
-    return std::accumulate(std::next(indices.begin()), indices.end(), ref,
-                           [&ref, &mol, box](const auto &acc, const auto idx)
-                           { return acc + box->minimumImage(mol.atom(idx)->r(), ref); }) /
-           indices.size();
-}
-
-// Calculate (mass-weighted) coordinate centre of atoms in the given molecule
-Vec3<double> SiteStack::centreOfMass(const Molecule &mol, const Box *box, const std::vector<int> &indices)
-{
-    auto mass = AtomicMass::mass(mol.atom(indices.front())->speciesAtom()->Z());
-    const auto ref = mol.atom(indices.front())->r();
-    auto sums = std::accumulate(std::next(indices.begin()), indices.end(), std::pair<Vec3<double>, double>(ref * mass, mass),
-                                [&ref, &mol, box](const auto &acc, const auto idx)
-                                {
-                                    auto mass = AtomicMass::mass(mol.atom(idx)->speciesAtom()->Z());
-                                    return std::pair<Vec3<double>, double>(
-                                        acc.first + box->minimumImage(mol.atom(idx)->r(), ref) * mass, acc.second + mass);
-                                });
-    return sums.first / sums.second;
-}
-
 // Create stack for specified Configuration and site
 bool SiteStack::create(Configuration *cfg, const SpeciesSite *site)
 {
@@ -95,26 +69,7 @@ bool SiteStack::create(Configuration *cfg, const SpeciesSite *site)
         auto index = 0;
         for (const auto &instance : instances)
         {
-            origin = speciesSite_->originMassWeighted() ? centreOfMass(*molecule, box, instance.originIndices())
-                                                        : centreOfGeometry(*molecule, box, instance.originIndices());
-
-            if (sitesHaveOrientation_)
-            {
-                // Get vector from site origin to x-axis reference point and normalise it
-                x = box->minimumVector(origin, centreOfGeometry(*molecule, box, instance.xAxisIndices()));
-                x.normalise();
-
-                // Get vector from site origin to y-axis reference point, normalise it, and orthogonalise
-                y = box->minimumVector(origin, centreOfGeometry(*molecule, box, instance.yAxisIndices()));
-                y.orthogonalise(x);
-                y.normalise();
-
-                sites_.emplace_back(speciesSite_, index, molecule, Matrix3(x, y, x * y), origin);
-            }
-            else
-                sites_.emplace_back(speciesSite_, index, molecule, origin);
-
-            ++index;
+            sites_.emplace_back(speciesSite_, index++, molecule, instance, box);
         }
     }
     return true;

diff --git a/src/classes/siteStack.h b/src/classes/siteStack.h
@@ -38,12 +38,6 @@ class SiteStack
     /*
      * Creation
      */
-    private:
-    // Calculate geometric centre of atoms in the given molecule
-    Vec3<double> centreOfGeometry(const Molecule &mol, const Box *box, const std::vector<int> &indices);
-    // Calculate (mass-weighted) coordinate centre of atoms in the given molecule
-    Vec3<double> centreOfMass(const Molecule &mol, const Box *box, const std::vector<int> &indices);
-
     public:
     // Create stack for specified Configuration and site
     bool create(Configuration *cfg, const SpeciesSite *site);

diff --git a/src/generator/CMakeLists.txt b/src/generator/CMakeLists.txt
@@ -4,6 +4,8 @@ add_library(
   add.h
   addBase.cpp
   addBase.h
+  addOnSphere.cpp
+  addOnSphere.h
   addPair.cpp
   addPair.h
   box.cpp