feat: Phantom atoms (#1953)

Co-authored-by: Adam Washington <[email protected]>

src/classes/atom.cpp

@@ -54,6 +54,9 @@ void Atom::setSpeciesAtom(const SpeciesAtom *spAtom) { speciesAtom_ = spAtom; }
 // Return SpeciesAtom that this Atom represents
 const SpeciesAtom *Atom::speciesAtom() const { return speciesAtom_; }
 
+// Return whether the atom's presence is that specified
+bool Atom::isPresence(SpeciesAtom::Presence presence) const { return speciesAtom_->isPresence(presence); }
+
 // Set Molecule in which this Atom exists
 void Atom::setMolecule(std::shared_ptr<Molecule> mol) { molecule_ = std::move(mol); }
 

src/classes/atom.h

@@ -67,6 +67,8 @@ class Atom
     void setSpeciesAtom(const SpeciesAtom *spAtom);
     // Return SpeciesAtom that this Atom represents
     const SpeciesAtom *speciesAtom() const;
+    // Return whether the atom's presence is that specified
+    bool isPresence(SpeciesAtom::Presence presence) const;
     // Set Molecule in which this Atom exists
     void setMolecule(std::shared_ptr<Molecule> mol);
     // Return Molecule in which this Atom exists

src/classes/atomType.h

@@ -51,6 +51,11 @@ class AtomType : public Serialisable<>
     int index_{-1};
 
     public:
+    // Enumeration for special type indices
+    enum SpecialTypeIndex
+    {
+        Ignore = -1
+    };
     // Return short-range interaction potential
     InteractionPotential<ShortRangeFunctions> &interactionPotential();
     const InteractionPotential<ShortRangeFunctions> &interactionPotential() const;
@@ -65,6 +70,10 @@ class AtomType : public Serialisable<>
     // Return whether our parameters are the same as those provided
     bool sameParametersAs(const AtomType *other, bool checkCharge = false);
 
+    /*
+     * I/O
+     */
+    public:
     // Express as a serialisable value
     SerialisedValue serialise() const override;
     // Read values from a serialisable value

src/classes/atomTypeMix.cpp

@@ -78,8 +78,7 @@ void AtomTypeMix::addIsotope(std::shared_ptr<AtomType> atomType, Sears91::Isotop
 // Finalise list, calculating fractional populations etc.
 void AtomTypeMix::finalise()
 {
-    // Finalise AtomTypeData
-    double total = totalPopulation();
+    auto total = totalPopulation();
     for (auto &atd : types_)
         atd.finalise(total);
 }
@@ -120,19 +119,15 @@ void AtomTypeMix::finalise(const std::vector<std::shared_ptr<AtomType>> &exchang
 // Make all AtomTypeData in the list reference only their natural isotope
 void AtomTypeMix::naturalise()
 {
-    // Loop over AtomTypes in the source list
     for (auto &atd : types_)
         atd.naturalise();
 }
 
 // Check for presence of AtomType
 bool AtomTypeMix::contains(const std::shared_ptr<AtomType> &atomType) const
 {
-    for (auto &atd : types_)
-        if (atd.atomType() == atomType)
-            return true;
-
-    return false;
+    return std::find_if(types_.begin(), types_.end(), [atomType](const auto &atd) { return atd.atomType() == atomType; }) !=
+           types_.end();
 }
 
 // Check for presence of AtomType/Isotope pair
@@ -156,35 +151,29 @@ std::vector<AtomTypeData>::const_iterator AtomTypeMix::begin() const { return ty
 std::vector<AtomTypeData>::const_iterator AtomTypeMix::end() const { return types_.end(); }
 
 // Return index of AtomType
-int AtomTypeMix::indexOf(const std::shared_ptr<AtomType> &atomtype) const
+std::optional<int> AtomTypeMix::indexOf(const std::shared_ptr<AtomType> &atomType) const
 {
-    auto count = 0;
-    for (auto &atd : types_)
-    {
-        if (atd.atomType() == atomtype)
-            return count;
-        ++count;
-    }
-
-    return -1;
+    auto it = std::find_if(types_.begin(), types_.end(), [atomType](const auto &atd) { return atd.atomType() == atomType; });
+    if (it == types_.end())
+        return {};
+    else
+        return it - types_.begin();
 }
 
 // Return index of names AtomType
-int AtomTypeMix::indexOf(std::string_view name) const
+std::optional<int> AtomTypeMix::indexOf(std::string_view name) const
 {
-    auto count = 0;
-    for (auto &atd : types_)
-    {
-        if (DissolveSys::sameString(atd.atomType()->name(), name))
-            return count;
-        ++count;
-    }
-
-    return -1;
+    auto it = std::find_if(types_.begin(), types_.end(),
+                           [name](const auto &atd) { return DissolveSys::sameString(atd.atomType()->name(), name); });
+    if (it == types_.end())
+        return {};
+    else
+        return it - types_.begin();
 }
 
 // Return indices of AtomType pair
-std::pair<int, int> AtomTypeMix::indexOf(const std::shared_ptr<AtomType> &at1, const std::shared_ptr<AtomType> &at2) const
+std::optional<std::pair<int, int>> AtomTypeMix::indexOf(const std::shared_ptr<AtomType> &at1,
+                                                        const std::shared_ptr<AtomType> &at2) const
 {
     auto count = 0, index = -1;
     for (auto &atd : types_)
@@ -194,19 +183,19 @@ std::pair<int, int> AtomTypeMix::indexOf(const std::shared_ptr<AtomType> &at1, c
             if (index == -1)
                 index = count;
             else
-                return {count, index};
+                return {{count, index}};
         }
         if (atd.atomType() == at2)
         {
             if (index == -1)
                 index = count;
             else
-                return {index, count};
+                return {{index, count}};
         }
         ++count;
     }
 
-    return {-1, -1};
+    return {};
 }
 
 // Return total population of all types

src/classes/atomTypeMix.h

@@ -61,11 +61,12 @@ class AtomTypeMix
     // Return ending iterator
     std::vector<AtomTypeData>::const_iterator end() const;
     // Return index of AtomType
-    int indexOf(const std::shared_ptr<AtomType> &atomtype) const;
+    std::optional<int> indexOf(const std::shared_ptr<AtomType> &atomType) const;
     // Return indices of AtomType pair
-    std::pair<int, int> indexOf(const std::shared_ptr<AtomType> &at1, const std::shared_ptr<AtomType> &at2) const;
+    std::optional<std::pair<int, int>> indexOf(const std::shared_ptr<AtomType> &at1,
+                                               const std::shared_ptr<AtomType> &at2) const;
     // Return index of named AtomType
-    int indexOf(std::string_view name) const;
+    std::optional<int> indexOf(std::string_view name) const;
     // Return total population of all types
     double totalPopulation() const;
     // Return nth referenced AtomType

src/classes/configuration.cpp

@@ -76,8 +76,8 @@ bool Configuration::generate(const GeneratorContext &generatorContext)
     // Make sure all objects know about each other
     updateObjectRelationships();
 
-    // Finalise used AtomType list
-    atomTypes_.finalise();
+    // Finalise atom type populations
+    atomTypePopulations_.finalise();
 
     // Link targeted potentials to atoms
     linkTargetedPotentials();

src/classes/configuration.h

@@ -81,8 +81,8 @@ class Configuration : public Serialisable<const CoreData &>
     private:
     // Species populations present in the Configuration
     std::vector<std::pair<const Species *, int>> speciesPopulations_;
-    // AtomType mix, containing unique (non-isotopic) atom types over all Species used in this configuration
-    AtomTypeMix atomTypes_;
+    // AtomType populations in the configuration
+    AtomTypeMix atomTypePopulations_;
     // Contents version, incremented whenever Configuration content or Atom positions change
     VersionCounter contentsVersion_;
     // Molecule vector
@@ -95,10 +95,8 @@ class Configuration : public Serialisable<const CoreData &>
     void empty();
     // Return specified used type
     std::shared_ptr<AtomType> atomTypes(int index);
-    // Return AtomTypeMix for this Configuration
-    const AtomTypeMix &atomTypes() const;
-    // Return number of atom types used in this Configuration
-    int nAtomTypes() const;
+    // Return atom type populations for this Configuration
+    const AtomTypeMix &atomTypePopulations() const;
     // Adjust population of specified Species in the Configuration
     void adjustSpeciesPopulation(const Species *sp, int delta);
     // Return Species populations within the Configuration
@@ -137,8 +135,8 @@ class Configuration : public Serialisable<const CoreData &>
     std::shared_ptr<Molecule> molecule(int n);
     // Add new Atom to Configuration
     Atom &addAtom(const SpeciesAtom *sourceAtom, const std::shared_ptr<Molecule> &molecule, Vec3<double> r = Vec3<double>());
-    // Return number of Atoms in Configuration
-    int nAtoms() const;
+    // Return the number of atoms in the configuration (or only those with the specified presence)
+    int nAtoms(SpeciesAtom::Presence withPresence = SpeciesAtom::Presence::Any) const;
     // Return Atom array
     std::vector<Atom> &atoms();
     const std::vector<Atom> &atoms() const;

src/classes/configuration_contents.cpp

@@ -14,7 +14,7 @@ void Configuration::empty()
 {
     molecules_.clear();
     atoms_.clear();
-    atomTypes_.clear();
+    atomTypePopulations_.clear();
     appliedSizeFactor_ = std::nullopt;
     speciesPopulations_.clear();
     globalPotentials_.clear();
@@ -24,11 +24,8 @@ void Configuration::empty()
     ++contentsVersion_;
 }
 
-// Return AtomTypeMix for this Configuration
-const AtomTypeMix &Configuration::atomTypes() const { return atomTypes_; }
-
-// Return number of atom types used in this Configuration
-int Configuration::nAtomTypes() const { return atomTypes_.nItems(); }
+// Return atom type populations for this Configuration
+const AtomTypeMix &Configuration::atomTypePopulations() const { return atomTypePopulations_; }
 
 // Adjust population of specified Species in the Configuration
 void Configuration::adjustSpeciesPopulation(const Species *sp, int delta)
@@ -91,7 +88,7 @@ std::optional<double> Configuration::atomicDensity() const
     if (nAtoms() == 0)
         return {};
 
-    return nAtoms() / box_->volume();
+    return nAtoms(SpeciesAtom::Presence::Physical) / box_->volume();
 }
 
 // Return the chemical density (g/cm3) of the Configuration
@@ -233,19 +230,30 @@ Atom &Configuration::addAtom(const SpeciesAtom *sourceAtom, const std::shared_pt
     // Set the position
     newAtom.setCoordinates(r);
 
-    // Update our typeIndex (non-isotopic) and set local and master type indices
-    if (sourceAtom->atomType() != nullptr)
+    // Update atom type population and set local type index
+    if (sourceAtom->isPresence(SpeciesAtom::Presence::Physical))
     {
-        auto &&[atd, atdIndex] = atomTypes_.add(sourceAtom->atomType(), 1);
+        auto &&[atd, atdIndex] = atomTypePopulations_.add(sourceAtom->atomType(), 1);
         newAtom.setLocalTypeIndex(atdIndex);
-        newAtom.setMasterTypeIndex(sourceAtom->atomType()->index());
     }
+    else
+        newAtom.setLocalTypeIndex(AtomType::Ignore);
+
+    // Set master index for pair potential lookup
+    newAtom.setMasterTypeIndex(sourceAtom->atomType()->index());
 
     return newAtom;
 }
 
-// Return number of Atoms in Configuration
-int Configuration::nAtoms() const { return atoms_.size(); }
+// Return the number of atoms in the configuration (or only those with the specified presence)
+int Configuration::nAtoms(SpeciesAtom::Presence withPresence) const
+{
+    if (withPresence == SpeciesAtom::Presence::Any)
+    {
+        return atoms_.size();
+    }
+    return std::count_if(atoms_.begin(), atoms_.end(), [withPresence](const auto &i) { return i.isPresence(withPresence); });
+}
 
 // Return Atom array
 std::vector<Atom> &Configuration::atoms() { return atoms_; }

src/classes/configuration_io.cpp

@@ -157,7 +157,7 @@ bool Configuration::deserialise(LineParser &parser, const CoreData &coreData, do
     }
 
     // Finalise used AtomType list
-    atomTypes_.finalise();
+    atomTypePopulations_.finalise();
 
     // Scale box and cells according to the applied size factor
     auto appliedSF = appliedSizeFactor_.value_or(defaultSizeFactor_);