diff --git a/corelib/src/libs/SireMM/amberparams.cpp b/corelib/src/libs/SireMM/amberparams.cpp index 541d2c2e7..b102966f4 100644 --- a/corelib/src/libs/SireMM/amberparams.cpp +++ b/corelib/src/libs/SireMM/amberparams.cpp @@ -1521,59 +1521,66 @@ QStringList AmberParams::validateAndFix() } } - // find the shortest bonded path between these two atoms - const auto path = conn.findPath(atm0, atm3, 4); + // find the shortest bonded paths between these two atoms + const auto paths = conn.findPaths(atm0, atm3, 4); - if (path.count() != 4) + for (const auto &path : paths) { - QMutexLocker lkr(&mutex); - errors.append( - QObject::tr("Have a 1-4 scaling factor (%1/%2) " - "between atoms %3:%4 and %5:%6 despite there being no physical " - "dihedral between these two atoms. All 1-4 scaling factors MUST " - "be associated with " - "physical dihedrals. The shortest path is %7") - .arg(s.coulomb()) - .arg(s.lj()) - .arg(molinfo.name(atm0).value()) - .arg(atm0.value()) - .arg(molinfo.name(atm3).value()) - .arg(atm3.value()) - .arg(Sire::toString(path))); - continue; - } + if (path.count() != 4) + { + QMutexLocker lkr(&mutex); + errors.append( + QObject::tr("Have a 1-4 scaling factor (%1/%2) " + "between atoms %3:%4 and %5:%6 despite there being no physical " + "dihedral between these two atoms. All 1-4 scaling factors MUST " + "be associated with " + "physical dihedrals. The shortest path is %7") + .arg(s.coulomb()) + .arg(s.lj()) + .arg(molinfo.name(atm0).value()) + .arg(atm0.value()) + .arg(molinfo.name(atm3).value()) + .arg(atm3.value()) + .arg(Sire::toString(path))); + continue; + } - // convert the atom IDs into a canonical form - auto dih = this->convert(DihedralID(path[0], path[1], path[2], path[3])); + // convert the atom IDs into a canonical form + auto dih = this->convert(DihedralID(path[0], path[1], path[2], path[3])); - // qDebug() << "ADDING NULL DIHEDRAL FOR" << dih.toString(); + // skip if we already have this dihedral + if (new_dihedrals.contains(dih)) + continue; - // does this bond involve hydrogen? - //- this relies on "AtomElements" being full - bool contains_hydrogen = false; + // qDebug() << "ADDING NULL DIHEDRAL FOR" << dih.toString(); - if (not amber_elements.isEmpty()) - { - contains_hydrogen = - (amber_elements.at(molinfo.cgAtomIdx(dih.atom0())).nProtons() < 2) or - (amber_elements.at(molinfo.cgAtomIdx(dih.atom1())).nProtons() < 2) or - (amber_elements.at(molinfo.cgAtomIdx(dih.atom2())).nProtons() < 2) or - (amber_elements.at(molinfo.cgAtomIdx(dih.atom3())).nProtons() < 2); - } + // does this bond involve hydrogen? + //- this relies on "AtomElements" being full + bool contains_hydrogen = false; - // create a null dihedral parameter and add this to the set - QMutexLocker lkr(&mutex); - new_dihedrals.insert( - dih, qMakePair(AmberDihedral(Expression(0), Symbol("phi")), contains_hydrogen)); + if (not amber_elements.isEmpty()) + { + contains_hydrogen = + (amber_elements.at(molinfo.cgAtomIdx(dih.atom0())).nProtons() < 2) or + (amber_elements.at(molinfo.cgAtomIdx(dih.atom1())).nProtons() < 2) or + (amber_elements.at(molinfo.cgAtomIdx(dih.atom2())).nProtons() < 2) or + (amber_elements.at(molinfo.cgAtomIdx(dih.atom3())).nProtons() < 2); + } - // now add in the 1-4 pair - BondID nb14pair = this->convert(BondID(dih.atom0(), dih.atom3())); + // create a null dihedral parameter and add this to the set + QMutexLocker lkr(&mutex); + new_dihedrals.insert( + dih, qMakePair(AmberDihedral(Expression(0), Symbol("phi")), contains_hydrogen)); - // add them to the list of 14 scale factors - new_nb14s.insert(nb14pair, AmberNB14(s.coulomb(), s.lj())); + // now add in the 1-4 pair + BondID nb14pair = this->convert(BondID(dih.atom0(), dih.atom3())); - // and remove them from the excluded atoms list - new_exc.set(nb14pair.atom0(), nb14pair.atom1(), CLJScaleFactor(0)); + // add them to the list of 14 scale factors + new_nb14s.insert(nb14pair, AmberNB14(s.coulomb(), s.lj())); + + // and remove them from the excluded atoms list + new_exc.set(nb14pair.atom0(), nb14pair.atom1(), CLJScaleFactor(0)); + } } } } diff --git a/corelib/src/libs/SireMol/connectivity.cpp b/corelib/src/libs/SireMol/connectivity.cpp index e884e9caf..0f927c374 100644 --- a/corelib/src/libs/SireMol/connectivity.cpp +++ b/corelib/src/libs/SireMol/connectivity.cpp @@ -1066,6 +1066,9 @@ QList ConnectivityBase::findPath(AtomIdx atom0, AtomIdx atom1) const { QList> paths = findPaths(atom0, atom1); + // Sort the paths so that the result is reproducible. + std::stable_sort(paths.begin(), paths.end()); + QList shortest; foreach (const QList &path, paths) @@ -1087,6 +1090,9 @@ QList ConnectivityBase::findPath(AtomIdx atom0, AtomIdx atom1, int max_ { QList> paths = findPaths(atom0, atom1, max_length); + // Sort the paths so that the result is reproducible. + std::stable_sort(paths.begin(), paths.end()); + QList shortest; foreach (const QList &path, paths) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 046c0f454..8a2759e55 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -21,6 +21,12 @@ organisation on `GitHub `__. * Set appropriate default OpenMM box vectors when a system has no periodic space. +* Make the result of ``Sire::Mol::Connectivity::findPath`` reproducible by sorting the + output paths. + +* Add missing dihedrals for all 1-4 atom paths in ``Sire::MM::AmberParams::validateAndFix``, + not just the first one found. + `2025.2.0 `__ - October 2025 --------------------------------------------------------------------------------------------