Skip to content

Backport fix from PR #372 #373

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lohedges merged 1 commit into from
Oct 22, 2025
Merged

Backport fix from PR #372 #373

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
93 changes: 50 additions & 43 deletions corelib/src/libs/SireMM/amberparams.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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));
}
}
}
}
Expand Down
6 changes: 6 additions & 0 deletions corelib/src/libs/SireMol/connectivity.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1066,6 +1066,9 @@ QList<AtomIdx> ConnectivityBase::findPath(AtomIdx atom0, AtomIdx atom1) const
{
QList<QList<AtomIdx>> paths = findPaths(atom0, atom1);

// Sort the paths so that the result is reproducible.
std::stable_sort(paths.begin(), paths.end());

QList<AtomIdx> shortest;

foreach (const QList<AtomIdx> &path, paths)
Expand All @@ -1087,6 +1090,9 @@ QList<AtomIdx> ConnectivityBase::findPath(AtomIdx atom0, AtomIdx atom1, int max_
{
QList<QList<AtomIdx>> paths = findPaths(atom0, atom1, max_length);

// Sort the paths so that the result is reproducible.
std::stable_sort(paths.begin(), paths.end());

QList<AtomIdx> shortest;

foreach (const QList<AtomIdx> &path, paths)
Expand Down
6 changes: 6 additions & 0 deletions doc/source/changelog.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,12 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.

* 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 <https://github.com/openbiosim/sire/compare/2025.1.0...2025.2.0>`__ - October 2025
--------------------------------------------------------------------------------------------

Expand Down