Merge pull request #407 from OpenBioSim/fix_glycam

Fix parsing of GLYCAM topology files

Author: lohedges

.clang-format

@@ -0,0 +1,11 @@
+BasedOnStyle: LLVM
+UseTab: Never
+IndentWidth: 4
+TabWidth: 4
+BreakBeforeBraces: Allman
+AllowShortIfStatementsOnASingleLine: false
+IndentCaseLabels: false
+ColumnLimit: 0
+AccessModifierOffset: -4
+NamespaceIndentation: All
+FixNamespaceComments: false

actions/generate_recipe.py

@@ -375,6 +375,9 @@ def generate_recipe(data, features, git_remote, git_branch, git_version, git_num
 
     # Script test (pytest)
     lines.append("  - script:")
+    lines.append("      - if: win")
+    lines.append("        then: set PYTHONUTF8=1")
+    lines.append("        else: export PYTHONUTF8=1")
     lines.append("      - pytest -vvv --color=yes --runveryslow ./tests")
     lines.append("    files:")
     lines.append("      source:")

corelib/src/libs/SireIO/biosimspace.cpp

@@ -32,11 +32,14 @@
 
 #include "SireError/errors.h"
 
+#include "SireMM/cljnbpairs.h"
+
 #include "SireMol/atomelements.h"
 #include "SireMol/atommasses.h"
 #include "SireMol/connectivity.h"
 #include "SireMol/core.h"
 #include "SireMol/mgname.h"
+#include "SireMol/moleculeinfodata.h"
 #include "SireMol/moleditor.h"
 #include "SireMol/molidx.h"
 
@@ -49,6 +52,7 @@
 
 using namespace SireBase;
 using namespace SireMaths;
+using namespace SireMM;
 using namespace SireMol;
 using namespace SireUnits;
 using namespace SireVol;
@@ -1627,7 +1631,7 @@ namespace SireIO
 
     Molecule createSodiumIon(const Vector &coords, const QString model, const PropertyMap &map)
     {
-		// Strip all whitespace from the model name and convert to upper case.
+        // Strip all whitespace from the model name and convert to upper case.
         auto _model = model.simplified().replace(" ", "").toUpper();
 
         // Create a hash between the allowed model names and their templace files.
@@ -1658,7 +1662,7 @@ namespace SireIO
 
     Molecule createChlorineIon(const Vector &coords, const QString model, const PropertyMap &map)
     {
-		// Strip all whitespace from the model name and convert to upper case.
+        // Strip all whitespace from the model name and convert to upper case.
         auto _model = model.simplified().replace(" ", "").toUpper();
 
         // Create a hash between the allowed model names and their templace files.
@@ -1730,9 +1734,9 @@ namespace SireIO
 
                 // Set the new coordinates.
                 molecule = molecule.edit()
-                                   .atom(AtomIdx(j))
-                                   .setProperty(coord_prop, coord)
-                                   .molecule();
+                               .atom(AtomIdx(j))
+                               .setProperty(coord_prop, coord)
+                               .molecule();
             }
 
             // Update the molecule in the system.
@@ -1754,4 +1758,75 @@ namespace SireIO
         return Vector(nx, ny, nz);
     }
 
+    SireBase::PropertyList mergeIntrascale(const CLJNBPairs &nb0,
+                                           const CLJNBPairs &nb1,
+                                           const MoleculeInfoData &merged_info,
+                                           const QHash<AtomIdx, AtomIdx> &mol0_merged_mapping,
+                                           const QHash<AtomIdx, AtomIdx> &mol1_merged_mapping)
+    {
+        // Helper lambda: copy the non-default scaling factors from 'nb' to
+        // 'nb_merged' according to the provided mapping. Takes nb_merged by
+        // reference to avoid copies.
+        auto copyIntrascale = [&](const CLJNBPairs &nb, CLJNBPairs &nb_merged,
+                                  const QHash<AtomIdx, AtomIdx> &mapping)
+        {
+            const int n = nb.nAtoms();
+
+            for (int i = 0; i < n; ++i)
+            {
+                const AtomIdx ai(i);
+
+                // Get the index of this atom in the merged system.
+                const AtomIdx merged_ai = mapping.value(ai, AtomIdx(-1));
+
+                // If this atom hasn't been mapped to the merged system, then we
+                // can skip it, as any scaling factors involving this atom will
+                // just use the default.
+                if (merged_ai == AtomIdx(-1))
+                    continue;
+
+                for (int j = i; j < n; ++j)
+                {
+                    const AtomIdx aj(j);
+
+                    // Get the scaling factor for this pair of atoms.
+                    const CLJScaleFactor sf = nb.get(ai, aj);
+
+                    // This is a non-default scaling factor, so we need to copy
+                    // it across to the merged intrascale object according to
+                    // the mapping.
+                    if (sf.coulomb() != 1.0 or sf.lj() != 1.0)
+                    {
+                        // Get the index of the second atom in the merged system.
+                        const AtomIdx merged_aj = mapping.value(aj, AtomIdx(-1));
+
+                        // Only set the scaling factor if both atoms have been
+                        // mapped to the merged system. If one of the atoms
+                        // hasn't been mapped, then we can just use the default.
+                        if (merged_aj != AtomIdx(-1))
+                            nb_merged.set(merged_ai, merged_aj, sf);
+                    }
+                }
+            }
+        };
+
+        // Create the intrascale objects for the merged end-states.
+        CLJNBPairs intra0(merged_info);
+        CLJNBPairs intra1(merged_info);
+
+        // Copy the non-default scaling factors from the original intrascale
+        // objects to the merged intrascale objects according to the provided
+        // mappings.
+        copyIntrascale(nb1, intra0, mol1_merged_mapping);
+        copyIntrascale(nb0, intra0, mol0_merged_mapping);
+        copyIntrascale(nb0, intra1, mol0_merged_mapping);
+        copyIntrascale(nb1, intra1, mol1_merged_mapping);
+
+        // Assemble the intrascale objects into a property list to return.
+        SireBase::PropertyList ret;
+        ret.append(intra0);
+        ret.append(intra1);
+        return ret;
+    }
+
 } // namespace SireIO

corelib/src/libs/SireIO/biosimspace.h

@@ -36,6 +36,12 @@
 
 #include "SireMaths/vector.h"
 
+#include "SireBase/propertylist.h"
+
+#include "SireMM/cljnbpairs.h"
+
+#include "SireMol/atomidxmapping.h"
+#include "SireMol/moleculeinfodata.h"
 #include "SireMol/select.h"
 
 SIRE_BEGIN_HEADER
@@ -372,6 +378,43 @@ namespace SireIO
         const bool is_lambda1 = false, const PropertyMap &map = PropertyMap());
 
     Vector cross(const Vector &v0, const Vector &v1);
+
+    //! Merge the CLJNBPairs (intrascale) of two molecules into a perturbable
+    /*! merged molecule's end-state intrascales.
+
+        Expands nb0 from molecule0's atom index space into the merged
+        molecule's space (preserving actual per-pair scale factors, including
+        force-field-specific overrides such as GLYCAM funct=2 pairs), then
+        calls CLJNBPairs::merge to produce intrascale0 and intrascale1.
+
+        \param nb0
+            The CLJNBPairs for molecule0 in its original atom index space.
+
+        \param nb1
+            The CLJNBPairs for molecule1 in its original atom index space.
+
+        \param merged_info
+            The MoleculeInfoData for the merged molecule.
+
+        \param mol0_merged_mapping
+            A hash mapping each AtomIdx in molecule0's original space to
+            the corresponding AtomIdx in the merged molecule's space.
+
+        \param atom_mapping
+            The AtomIdxMapping describing how merged-molecule atom indices
+            correspond to molecule1 atom indices (used by CLJNBPairs::merge).
+
+        \retval [intrascale0, intrascale1]
+            A PropertyList containing the CLJNBPairs for the lambda=0 and
+            lambda=1 end states of the merged molecule.
+     */
+    SIREIO_EXPORT SireBase::PropertyList mergeIntrascale(
+        const SireMM::CLJNBPairs &nb0,
+        const SireMM::CLJNBPairs &nb1,
+        const SireMol::MoleculeInfoData &merged_info,
+        const QHash<SireMol::AtomIdx, SireMol::AtomIdx> &mol0_merged_mapping,
+        const QHash<SireMol::AtomIdx, SireMol::AtomIdx> &mol1_merged_mapping);
+
 } // namespace SireIO
 
 SIRE_EXPOSE_FUNCTION(SireIO::isAmberWater)
@@ -387,6 +430,7 @@ SIRE_EXPOSE_FUNCTION(SireIO::updateCoordinatesAndVelocities)
 SIRE_EXPOSE_FUNCTION(SireIO::createSodiumIon)
 SIRE_EXPOSE_FUNCTION(SireIO::createChlorineIon)
 SIRE_EXPOSE_FUNCTION(SireIO::setCoordinates)
+SIRE_EXPOSE_FUNCTION(SireIO::mergeIntrascale)
 
 SIRE_END_HEADER
 

corelib/src/libs/SireIO/gro87.cpp

@@ -1889,43 +1889,59 @@ System Gro87::startSystem(const PropertyMap &map) const
 
         int ncg = 0;
 
-        QSet<ResNum> completed_residues;
+        // Track used residue numbers to handle topologies where numbering
+        // restarts (e.g. glycan residues numbered from 1 after a protein
+        // chain also starting from 1). Duplicate ResNums cause errors when
+        // looking up residues by number.
+        // next_unique is always > every number assigned so far, so conflict
+        // resolution is O(1) rather than a linear scan.
+        QSet<int> used_resnums;
+        int next_unique = 0;
+        auto unique_resnum = [&](int resnum) -> ResNum
+        {
+            if (!used_resnums.contains(resnum))
+            {
+                used_resnums.insert(resnum);
+                if (resnum >= next_unique)
+                    next_unique = resnum + 1;
+                return ResNum(resnum);
+            }
+            // conflict: assign next number beyond all previously seen
+            used_resnums.insert(next_unique);
+            return ResNum(next_unique++);
+        };
+
+        // Track the original residue number/name to detect residue boundaries.
+        // We compare against the original topology values (not remapped numbers).
+        // Store editors for the current residue and CutGroup so reparenting uses
+        // O(1) index lookups rather than O(n) name scans.
+        int current_orig_resnum = -1;
+        QString current_resname;
+        ResStructureEditor current_res;
+        CGStructureEditor current_cg;
 
         for (int i = 0; i < atmnams.count(); ++i)
         {
             auto atom = moleditor.add(AtomNum(atmnums[i]));
             atom = atom.rename(AtomName(atmnams[i]));
 
-            const ResNum resnum(resnums[i]);
-
-            if (completed_residues.contains(resnum))
-            {
-                auto res = moleditor.residue(resnum);
-
-                if (res.name().value() != resnams[i])
-                {
-                    // different residue
-                    res = moleditor.add(resnum);
-                    res = res.rename(ResName(resnams[i]));
-                    ncg += 1;
-                    moleditor.add(CGName(QString::number(ncg)));
-                }
+            const int orig_resnum = resnums[i];
+            const QString &resnam = resnams[i];
 
-                atom = atom.reparent(res.index());
-                atom = atom.reparent(CGName(QString::number(ncg)));
-            }
-            else
+            if (orig_resnum != current_orig_resnum || resnam != current_resname)
             {
-                auto res = moleditor.add(resnum);
-                res = res.rename(ResName(resnams[i]));
-                atom = atom.reparent(res.index());
+                // new residue
+                current_res = moleditor.add(unique_resnum(orig_resnum));
+                current_res = current_res.rename(ResName(resnam));
+                current_orig_resnum = orig_resnum;
+                current_resname = resnam;
 
                 ncg += 1;
-                auto cg = moleditor.add(CGName(QString::number(ncg)));
-                atom = atom.reparent(cg.index());
-
-                completed_residues.insert(resnum);
+                current_cg = moleditor.add(CGName(QString::number(ncg)));
             }
+
+            atom = atom.reparent(current_res.index());
+            atom = atom.reparent(current_cg.index());
         }
 
         // we have created the molecule - now add in the coordinates/velocities as needed