Merge branch 'OpenBioSim:devel' into feature_DMR

Author: akalpokas

.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/SireCAS/lambdaschedule.cpp

@@ -65,7 +65,7 @@ static RegisterMetaType<LambdaSchedule> r_schedule;
 
 QDataStream &operator<<(QDataStream &ds, const LambdaSchedule &schedule)
 {
-    writeHeader(ds, r_schedule, 3);
+    writeHeader(ds, r_schedule, 4);
 
     SharedDataStream sds(ds);
 
@@ -75,6 +75,7 @@ QDataStream &operator<<(QDataStream &ds, const LambdaSchedule &schedule)
         << schedule.default_equations
         << schedule.stage_equations
         << schedule.mol_schedules
+        << schedule.coupled_levers
         << static_cast<const Property &>(schedule);
 
     return ds;
@@ -91,21 +92,30 @@ QDataStream &operator>>(QDataStream &ds, LambdaSchedule &schedule)
 {
     VersionID v = readHeader(ds, r_schedule);
 
-    if (v == 1 or v == 2 or v == 3)
+    if (v == 1 or v == 2 or v == 3 or v == 4)
     {
         SharedDataStream sds(ds);
 
         sds >> schedule.constant_values;
 
-        if (v == 3)
+        if (v >= 3)
             sds >> schedule.force_names;
 
         sds >> schedule.lever_names >> schedule.stage_names >>
             schedule.default_equations >> schedule.stage_equations;
 
-        if (v == 2 or v == 3)
+        if (v >= 2)
             sds >> schedule.mol_schedules;
 
+        if (v >= 4)
+            sds >> schedule.coupled_levers;
+        else
+        {
+            // Populate the default coupling for streams written before v4
+            schedule.coupled_levers[_get_lever_name("cmap", "cmap_grid")] =
+                _get_lever_name("torsion", "torsion_k");
+        }
+
         if (v < 3)
         {
             // need to make sure that the lever names are namespaced
@@ -146,13 +156,17 @@ QDataStream &operator>>(QDataStream &ds, LambdaSchedule &schedule)
         }
     }
     else
-        throw version_error(v, "1, 2, 3", r_schedule, CODELOC);
+        throw version_error(v, "1, 2, 3, 4", r_schedule, CODELOC);
 
     return ds;
 }
 
 LambdaSchedule::LambdaSchedule() : ConcreteProperty<LambdaSchedule, Property>()
 {
+    // By default, cmap_grid falls back to torsion_k so that customising
+    // torsion scaling automatically keeps the CMAP correction in sync.
+    coupled_levers[_get_lever_name("cmap", "cmap_grid")] =
+        _get_lever_name("torsion", "torsion_k");
 }
 
 LambdaSchedule::LambdaSchedule(const LambdaSchedule &other)
@@ -162,7 +176,8 @@ LambdaSchedule::LambdaSchedule(const LambdaSchedule &other)
       force_names(other.force_names),
       lever_names(other.lever_names), stage_names(other.stage_names),
       default_equations(other.default_equations),
-      stage_equations(other.stage_equations)
+      stage_equations(other.stage_equations),
+      coupled_levers(other.coupled_levers)
 {
 }
 
@@ -1006,6 +1021,32 @@ void LambdaSchedule::removeEquation(const QString &stage,
     this->stage_equations[idx].remove(_get_lever_name(force, lever));
 }
 
+/** Couple the lever 'force::lever' to 'fallback_force::fallback_lever'.
+ *  If no custom equation has been set for 'force::lever' at any stage,
+ *  the equation for 'fallback_force::fallback_lever' will be used instead
+ *  of the stage default. This is a single level of indirection — the
+ *  fallback lever is not itself followed further.
+ *
+ *  By default, 'cmap::cmap_grid' is coupled to 'torsion::torsion_k' so
+ *  that custom torsion schedules automatically keep the CMAP correction
+ *  in sync.
+ */
+void LambdaSchedule::coupleLever(const QString &force, const QString &lever,
+                                  const QString &fallback_force,
+                                  const QString &fallback_lever)
+{
+    coupled_levers[_get_lever_name(force, lever)] =
+        _get_lever_name(fallback_force, fallback_lever);
+}
+
+/** Remove any coupling for the lever 'force::lever', reverting it to
+ *  use the stage default equation when no custom equation is set.
+ */
+void LambdaSchedule::removeCoupledLever(const QString &force, const QString &lever)
+{
+    coupled_levers.remove(_get_lever_name(force, lever));
+}
+
 /** Return whether or not the specified 'lever' in the specified 'force'
  *  at the specified 'stage' has a custom equation set for it
  */
@@ -1077,6 +1118,32 @@ SireCAS::Expression LambdaSchedule::_getEquation(int stage,
         return it.value();
     }
 
+    // Check coupled levers: if this lever is coupled to another, try to
+    // find a custom equation for the coupled lever before falling back to
+    // the stage default. This is a single level of indirection (no recursion)
+    // to prevent loops.
+    auto coupled_it = this->coupled_levers.find(lever_name);
+
+    if (coupled_it != this->coupled_levers.end())
+    {
+        const auto &coupled_name = coupled_it.value();
+        const int sep = coupled_name.indexOf("::");
+        const QString coupled_force = sep >= 0 ? coupled_name.left(sep) : "*";
+        const QString coupled_lever_part = sep >= 0 ? coupled_name.mid(sep + 2) : coupled_name;
+
+        it = equations.find(coupled_name);
+        if (it != equations.end())
+            return it.value();
+
+        it = equations.find(_get_lever_name(coupled_force, "*"));
+        if (it != equations.end())
+            return it.value();
+
+        it = equations.find(_get_lever_name("*", coupled_lever_part));
+        if (it != equations.end())
+            return it.value();
+    }
+
     // we don't have any match, so return the default equation for this stage
     return this->default_equations[stage];
 }

corelib/src/libs/SireCAS/lambdaschedule.h

@@ -164,6 +164,12 @@ namespace SireCAS
                                       const QString &force = "*",
                                       const QString &lever = "*") const;
 
+        void coupleLever(const QString &force, const QString &lever,
+                         const QString &fallback_force,
+                         const QString &fallback_lever);
+
+        void removeCoupledLever(const QString &force, const QString &lever);
+
         SireCAS::Expression getEquation(const QString &stage = "*",
                                         const QString &force = "*",
                                         const QString &lever = "*") const;
@@ -248,6 +254,12 @@ namespace SireCAS
             particular stage */
         QVector<QHash<QString, SireCAS::Expression>> stage_equations;
 
+        /** Coupled lever fallbacks: if a lever (force::lever) has no custom
+            equation set, use the equation for the paired lever instead of
+            falling straight back to the stage default. Stored as
+            _get_lever_name(force, lever) → _get_lever_name(fb_force, fb_lever). */
+        QHash<QString, QString> coupled_levers;
+
         /** The symbol used to represent the lambda value */
         static SireCAS::Symbol lambda_symbol;
 

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