Merge pull request #434 from OpenBioSim/fix_ring_break_merge

Fix mergeIntrascale using per-state connectivity

Author: lohedges

corelib/src/libs/SireIO/biosimspace.cpp

@@ -1760,76 +1760,58 @@ namespace SireIO
 
     SireBase::PropertyList mergeIntrascale(const CLJNBPairs &nb0,
                                            const CLJNBPairs &nb1,
-                                           const MoleculeInfoData &merged_info,
+                                           const Connectivity &conn0,
+                                           const Connectivity &conn1,
+                                           const CLJScaleFactor &sf14,
                                            const QHash<AtomIdx, AtomIdx> &mol0_merged_mapping,
                                            const QHash<AtomIdx, AtomIdx> &mol1_merged_mapping)
     {
-        // Helper lambda: copy scaling factors from 'nb' to 'nb_merged' according
-        // to the provided mapping. Takes nb_merged by reference to avoid copies.
-        // When copy_all is true, ALL pairs between mapped atoms are written
-        // (including the default (1,1)), so that the correct end-state values
-        // always overwrite any values set by a prior ghost-topology pass.
-        // When copy_all is false, only non-(1,1) pairs are written (sufficient
-        // for the first/ghost pass, where unset pairs default to (1,1) anyway).
-        auto copyIntrascale = [&](const CLJNBPairs &nb, CLJNBPairs &nb_merged,
-                                  const QHash<AtomIdx, AtomIdx> &mapping,
-                                  bool copy_all = false)
+        // Build base CLJNBPairs from the per-state merged connectivity. This
+        // correctly captures bonded distances (1-2, 1-3, 1-4) in the merged
+        // atom space, including paths that only exist in one end state (e.g.
+        // the ring-closure bond for ring-breaking perturbations).
+        CLJNBPairs intra0(conn0, sf14);
+        CLJNBPairs intra1(conn1, sf14);
+
+        // Override with per-pair values from the individual molecule intrascales.
+        // For standard AMBER molecules these are identical to the connectivity-
+        // derived values and the override is a no-op. For force fields with
+        // non-default per-pair scale factors (e.g. GLYCAM funct=2 with (1,1)
+        // instead of global sf14 for certain 1-4 pairs) the override replaces
+        // the global-sf14 value set by the base construction with the correct
+        // per-pair value.
+        auto overrideIntrascale = [&](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);
+                    const AtomIdx merged_aj = mapping.value(aj, AtomIdx(-1));
 
-                    // Get the scaling factor for this pair of atoms.
-                    const CLJScaleFactor sf = nb.get(ai, aj);
+                    if (merged_aj == AtomIdx(-1))
+                        continue;
 
-                    // Copy if this is a non-default scaling factor, or if
-                    // copy_all is set (second pass: must overwrite ghost-topology
-                    // values even when the correct end-state value is (1,1)).
-                    if (copy_all or 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);
-                    }
+                    const CLJScaleFactor nb_sf = nb.get(ai, aj);
+                    const CLJScaleFactor base_sf = nb_merged.get(merged_ai, merged_aj);
+
+                    if (nb_sf.coulomb() != base_sf.coulomb() or nb_sf.lj() != base_sf.lj())
+                        nb_merged.set(merged_ai, merged_aj, nb_sf);
                 }
             }
         };
 
-        // Create the intrascale objects for the merged end-states.
-        CLJNBPairs intra0(merged_info);
-        CLJNBPairs intra1(merged_info);
-
-        // Copy scaling factors from the original intrascale objects to the
-        // merged intrascale objects. For each end state, the ghost molecule's
-        // topology is written first (non-default pairs only), then the correct
-        // end-state topology overwrites with copy_all=true so that (1,1) pairs
-        // are also written. This handles ring-breaking perturbations where a
-        // pair that is excluded/1-4 in one state becomes fully interacting (1,1)
-        // in the other state and must not be left at the ghost state's value.
-        copyIntrascale(nb1, intra0, mol1_merged_mapping);
-        copyIntrascale(nb0, intra0, mol0_merged_mapping, true);
-        copyIntrascale(nb0, intra1, mol0_merged_mapping);
-        copyIntrascale(nb1, intra1, mol1_merged_mapping, true);
+        overrideIntrascale(nb0, intra0, mol0_merged_mapping);
+        overrideIntrascale(nb1, intra1, mol1_merged_mapping);
 
         // Assemble the intrascale objects into a property list to return.
         SireBase::PropertyList ret;

corelib/src/libs/SireIO/biosimspace.h

@@ -41,6 +41,7 @@
 #include "SireMM/cljnbpairs.h"
 
 #include "SireMol/atomidxmapping.h"
+#include "SireMol/connectivity.h"
 #include "SireMol/moleculeinfodata.h"
 #include "SireMol/select.h"
 
@@ -411,7 +412,9 @@ namespace SireIO
     SIREIO_EXPORT SireBase::PropertyList mergeIntrascale(
         const SireMM::CLJNBPairs &nb0,
         const SireMM::CLJNBPairs &nb1,
-        const SireMol::MoleculeInfoData &merged_info,
+        const SireMol::Connectivity &conn0,
+        const SireMol::Connectivity &conn1,
+        const SireMM::CLJScaleFactor &sf14,
         const QHash<SireMol::AtomIdx, SireMol::AtomIdx> &mol0_merged_mapping,
         const QHash<SireMol::AtomIdx, SireMol::AtomIdx> &mol1_merged_mapping);