Fix 199

corelib/src/libs/SireMol/structureeditor.cpp

@@ -3210,6 +3210,134 @@ const MoleculeInfoData &StructureEditor::commitInfo()
             }
         }
 
+        // make sure that the AtomIdx order matches the CGAtomIdx order
+        // This is needed to remove confusion in code that assumes
+        // AtomIdx order is always correct. This is a big change in code
+        // behaviour, but it is necessary to make the code more robust,
+        // and also recognises the reality that CutGroups are now "under
+        // the hood" in the code and not really visible or known about by
+        // end users.
+        int atom_count = 0;
+        bool in_atomidx_order = true;
+
+        for (int i = 0; i < this->nCutGroupsInMolecule(); ++i)
+        {
+            const auto cgdata = this->getCGData(CGIdx(i));
+
+            for (const auto &atomidx : cgdata.get<1>())
+            {
+                if (atomidx != AtomIdx(atom_count))
+                {
+                    in_atomidx_order = false;
+                    break;
+                }
+
+                ++atom_count;
+            }
+        }
+
+        if (atom_count != this->nAtomsInMolecule())
+        {
+            in_atomidx_order = false;
+        }
+
+        if (not in_atomidx_order)
+        {
+            SireBase::Console::warning(QObject::tr(
+                "The atoms in the CutGroups are not in the same order as the atoms in the molecule. "
+                "Rebuilding CutGroups so that the CGAtomIdx order matches the AtomIdx order."));
+
+            // are the atoms contiguous in residues?
+            bool contiguous_residues = true;
+            QSet<ResIdx> seen_residues;
+            seen_residues.reserve(this->nResiduesInMolecule());
+            ResIdx prev_residx = ResIdx::null();
+
+            for (int i = 0; i < this->nAtomsInMolecule(); ++i)
+            {
+                const auto info = this->getAtomData(AtomIdx(i));
+
+                if (prev_residx.isNull())
+                {
+                    prev_residx = info.get<4>();
+                    seen_residues.insert(prev_residx);
+                }
+                else
+                {
+                    if (info.get<4>() != prev_residx)
+                    {
+                        prev_residx = info.get<4>();
+
+                        if (seen_residues.contains(prev_residx))
+                        {
+                            contiguous_residues = false;
+                            break;
+                        }
+                        else
+                        {
+                            seen_residues.insert(prev_residx);
+                        }
+                    }
+                }
+            }
+
+            // delete all existing CutGroups
+            this->removeAllCutGroups();
+
+            if (contiguous_residues)
+            {
+                prev_residx = ResIdx::null();
+                int cgcount = 0;
+
+                for (int i = 0; i < this->nAtomsInMolecule(); ++i)
+                {
+                    const auto info = this->getAtomData(AtomIdx(i));
+
+                    if (prev_residx.isNull() or info.get<4>() != prev_residx)
+                    {
+                        prev_residx = info.get<4>();
+                        this->addCutGroup().rename(CGName(QString::number(cgcount)));
+                        cgcount += 1;
+                    }
+
+                    this->reparentAtom(this->getUID(AtomIdx(i)), CGIdx(cgcount - 1));
+                }
+
+                if (cgcount != this->nResiduesInMolecule())
+                {
+                    SireBase::Console::warning(QObject::tr(
+                                                   "The number of CutGroups created (%1) does not match "
+                                                   "the number of residues in the molecule (%2)! Rebuilding "
+                                                   "into a single CutGroup.")
+                                                   .arg(cgcount)
+                                                   .arg(this->nResiduesInMolecule()));
+
+                    this->removeAllCutGroups();
+
+                    this->addCutGroup().rename(CGName("0"));
+
+                    for (int i = 0; i < this->nAtomsInMolecule(); ++i)
+                    {
+                        this->reparentAtom(this->getUID(AtomIdx(i)), CGIdx(0));
+                    }
+                }
+            }
+            else
+            {
+                // create a new CutGroup and add all atoms to it
+                SireBase::Console::warning(QObject::tr(
+                    "The atoms in the molecule are not contiguous in residues. "
+                    "Rebuilding into a single CutGroup."));
+
+                this->addCutGroup().rename(CGName("0"));
+
+                for (int i = 0; i < this->nAtomsInMolecule(); ++i)
+                {
+                    this->reparentAtom(this->getUID(AtomIdx(i)), CGIdx(0));
+                }
+            }
+        }
+
         d->cached_molinfo = new MoleculeInfoData(*this);
     }
 

corelib/src/libs/SireSystem/merge.cpp

@@ -390,12 +390,26 @@ namespace SireSystem
 
                 auto res = mol.residue(residx);
 
+                // get the AtomIdx of the last atom in this residue
+                AtomIdx last_atomidx(0);
+
+                if (res.nAtoms() > 0)
+                {
+                    last_atomidx = res.atom(res.nAtoms() - 1).index();
+                }
+
                 // add the atom - it has the name "Xxx" as it doesn't exist
                 // in the reference state
                 auto atom = res.add(AtomName("Xxx"));
                 largest_atomnum = AtomNum(largest_atomnum.value() + 1);
                 atom.renumber(largest_atomnum);
 
+                // ensure that its index follows on from the index of the
+                // last atom in the residue - this is so that we keep
+                // the AtomIdx and CGAtomIdx orders in sync, and don't
+                // force a complex reordering of the atoms when we commit
+                atom.reindex(last_atomidx + 1);
+
                 // reparent this atom to the CutGroup for this residue
                 atom.reparent(cgidx);
 

corelib/src/libs/SireSystem/system.cpp

@@ -4162,6 +4162,65 @@ void System::makeWhole()
     this->makeWhole(PropertyMap());
 }
 
+void System::makeWhole(const Vector &center, const PropertyMap &map)
+{
+    if (this->needsAccepting())
+    {
+        this->accept();
+    }
+
+    if (not this->containsProperty(map["space"]))
+        return;
+
+    if (not this->property(map["space"]).isA<Space>())
+        return;
+
+    const auto &space = this->property(map["space"]).asA<Space>();
+
+    if (not space.isPeriodic())
+        return;
+
+    PropertyMap m = map;
+    m.set("space", space);
+
+    // get a list of all molecules in the system
+    const SelectorMol mols(*this);
+
+    SelectorMol changed_mols;
+
+    for (const auto &mol : mols)
+    {
+        auto new_mol = mol.move().makeWhole(center, m).commit();
+
+        if (new_mol.data().version() != mol.data().version())
+        {
+            changed_mols.append(new_mol);
+        }
+    }
+
+    if (not changed_mols.isEmpty())
+    {
+        Delta delta(*this, true);
+
+        // this ensures that only a single copy of System is used - prevents
+        // unnecessary copying
+        this->operator=(System());
+        delta.update(changed_mols.toMolecules());
+        this->operator=(delta.apply());
+
+        if (this->needsAccepting())
+        {
+            delta = Delta();
+            this->accept();
+        }
+    }
+}
+
+void System::makeWhole(const Vector &center)
+{
+    this->makeWhole(center, PropertyMap());
+}
+
 const char *System::typeName()
 {
     return QMetaType::typeName(qMetaTypeId<System>());

corelib/src/libs/SireSystem/system.h

@@ -443,6 +443,10 @@ namespace SireSystem
         void makeWhole();
         void makeWhole(const SireBase::PropertyMap &map);
 
+        void makeWhole(const SireMaths::Vector &center);
+        void makeWhole(const SireMaths::Vector &center,
+                       const SireBase::PropertyMap &map);
+
         static const System &null();
 
     protected:

doc/source/changelog.rst

@@ -55,6 +55,26 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 * Added a map option (fix_perturbable_zero_sigmas) to prevent perturbation of
   the Lennard-Jones sigma parameter for ghost atoms during alchemical free energy simulations.
 
+* [CHANGE IN BEHAVIOUR] - added code that ensures that, when editing molecules,
+  the CGAtomIdx order will always follow the AtomIdx order of atoms. This is
+  because a lot of code had implicitly assumed this, and so it was a cause
+  of bugs when this wasn't the case. Now, when you edit a molecule, on committing,
+  the orders will be checked. If they don't agree, then the CutGroups will be
+  reordered, with atoms reordered as necessary to make the CGAtomIdx order match
+  the AtomIdx order. If this isn't possible (e.g. because atoms in CutGroups
+  are not contiguous), then the molecule will be converted to a single-cutgroup
+  molecule, with the atoms placed in AtomIdx order. As part of this change,
+  the merge code will now also ensure that added atoms are added with the
+  correct AtomIdx, rather than added as the last atoms in the molecule. This
+  is also more natural. This fixes issue #202.
+
+* Added the "center" keyword argument to the ``make_whole`` functions of
+  :class:`~sire.mol.Cursors`, :class:`~sire.mol.CursorsM` and
+  :class:`~sire.system.System` (as well as to the legacy System class).
+  Also allowed the constructor of :class:`~sire.maths.Vector` to recognise
+  ``origin`` and ``zero`` as arguments, meaning you can write
+  ``cursor.make_whole(center="origin")``. This fixes issue #199.
+
 * Please add an item to this changelog when you create your PR
 
 `2024.1.0 <https://github.com/openbiosim/sire/compare/2023.5.2...2024.1.0>`__ - April 2024

src/sire/maths/_vector.py

@@ -237,6 +237,13 @@ class containing 3 double precision values. These values
     """
 
     def __init__(self, *args, **kwargs):
+        if len(args) == 1:
+            # check for "zero" or "origin"
+            arg0 = str(args[0]).strip().lower()
+
+            if arg0 == "zero" or arg0 == "origin":
+                args = [0.0, 0.0, 0.0]
+
         from ..units import angstrom
         from .. import u
 

src/sire/mol/_cursor.py

@@ -2624,7 +2624,7 @@ def delete_frame(self, *args, **kwargs):
 
         return self
 
-    def make_whole(self, *args, map=None):
+    def make_whole(self, center=None, map=None):
         """
         Make all of the atoms operated on by this cursor whole
         (they won't be broken across a periodic box boundary)
@@ -2634,7 +2634,7 @@ def make_whole(self, *args, map=None):
         which they should be wrapped.
         """
         for cursor in self._cursors:
-            cursor.make_whole(*args, map=map)
+            cursor.make_whole(center=center, map=map)
 
         return self
 
@@ -3779,7 +3779,7 @@ def delete_frame(self, *args, **kwargs):
 
         return self
 
-    def make_whole(self, *args, map=None):
+    def make_whole(self, center=None, map=None):
         """
         Make all of the atoms operated on by this cursor whole
         (they won't be broken across a periodic box boundary)
@@ -3789,7 +3789,7 @@ def make_whole(self, *args, map=None):
         which they should be wrapped.
         """
         for cursor in self._cursors:
-            cursor.make_whole(*args, map=map)
+            cursor.make_whole(center=center, map=map)
 
         return self
 

src/sire/system/_system.py

@@ -145,18 +145,30 @@ def numbers(self):
         """Return the numbers of all of the molecules in this System"""
         return self.molecules().numbers()
 
-    def make_whole(self, map=None):
+    def make_whole(self, center=None, map=None):
         """
         Make all of the molecules in this system whole. This
         maps each molecule into the current space, such that no
         molecule is broken across a periodic box boundary
         """
-        if map is None:
-            self._system.make_whole()
+        if center is None:
+            if map is None:
+                self._system.make_whole()
+            else:
+                from ..base import create_map
+
+                self._system.make_whole(map=create_map(map))
         else:
-            from ..base import create_map
+            from ..maths import Vector
+
+            center = Vector(center)
+
+            if map is None:
+                self._system.make_whole(center=center)
+            else:
+                from ..base import create_map
 
-            self._system.make_whole(map=create_map(map))
+                self._system.make_whole(center=center, map=create_map(map))
 
         self._molecules = None
 

tests/mol/test_make_whole.py

@@ -59,9 +59,7 @@ def test_auto_make_whole_on_load_frame(wrapped_mols):
 
 
 def test_auto_make_whole_on_load():
-    mols = sr.load_test_files(
-        "wrapped.rst7", "wrapped.prm7", map={"make_whole": True}
-    )
+    mols = sr.load_test_files("wrapped.rst7", "wrapped.prm7", map={"make_whole": True})
 
     _assert_correct_com(mols[0].evaluate().center_of_mass())
 
@@ -84,3 +82,19 @@ def test_auto_make_whole_on_load_no_breakage(kigaki_mols):
         kigaki_mols[0].evaluate().center_of_mass()
         == mols[0].evaluate().center_of_mass()
     )
+
+
+def test_make_whole_center_args(ala_mols):
+    mols = ala_mols
+
+    c = mols[0].cursor()
+    c.make_whole(center="origin")
+
+    c = mols.cursor()
+    c.make_whole(center=0)
+
+    c = mols[0].atoms().cursor()
+    c.make_whole(center=(1, 2, 3))
+
+    mols = mols.clone()
+    mols.make_whole(center=("1A", "2A", "3A"))