diff --git a/.clang-format b/.clang-format new file mode 100644 index 000000000..45eee0732 --- /dev/null +++ b/.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 diff --git a/actions/generate_recipe.py b/actions/generate_recipe.py index 8010a0e0c..09be4c500 100644 --- a/actions/generate_recipe.py +++ b/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:") diff --git a/corelib/src/libs/SireIO/biosimspace.cpp b/corelib/src/libs/SireIO/biosimspace.cpp index 89dd01ca5..f0a2d6222 100644 --- a/corelib/src/libs/SireIO/biosimspace.cpp +++ b/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 &mol0_merged_mapping, + const QHash &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 &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 diff --git a/corelib/src/libs/SireIO/biosimspace.h b/corelib/src/libs/SireIO/biosimspace.h index 081735cff..6e7f64612 100644 --- a/corelib/src/libs/SireIO/biosimspace.h +++ b/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 &mol0_merged_mapping, + const QHash &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 diff --git a/corelib/src/libs/SireIO/gro87.cpp b/corelib/src/libs/SireIO/gro87.cpp index adb44d753..55c9d6484 100644 --- a/corelib/src/libs/SireIO/gro87.cpp +++ b/corelib/src/libs/SireIO/gro87.cpp @@ -1889,43 +1889,59 @@ System Gro87::startSystem(const PropertyMap &map) const int ncg = 0; - QSet 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 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 diff --git a/corelib/src/libs/SireIO/grotop.cpp b/corelib/src/libs/SireIO/grotop.cpp index ef84a1e70..fd44d6ed9 100644 --- a/corelib/src/libs/SireIO/grotop.cpp +++ b/corelib/src/libs/SireIO/grotop.cpp @@ -46,11 +46,11 @@ #include "SireMM/atomljs.h" #include "SireMM/cljnbpairs.h" +#include "SireMM/cmapfunctions.h" #include "SireMM/fouratomfunctions.h" #include "SireMM/internalff.h" #include "SireMM/threeatomfunctions.h" #include "SireMM/twoatomfunctions.h" -#include "SireMM/cmapfunctions.h" #include "SireBase/booleanproperty.h" #include "SireBase/parallel.h" @@ -135,7 +135,7 @@ QDataStream &operator>>(QDataStream &ds, GroAtom &atom) atom.mss = mass * g_per_mol; } else - throw version_error(v, "1,2", r_groatom, CODELOC); + throw version_error(v, "1,2,3", r_groatom, CODELOC); return ds; } @@ -351,7 +351,7 @@ static const RegisterMetaType r_gromoltyp(NO_ROOT); QDataStream &operator<<(QDataStream &ds, const GroMolType &moltyp) { - writeHeader(ds, r_gromoltyp, 3); + writeHeader(ds, r_gromoltyp, 4); SharedDataStream sds(ds); @@ -367,8 +367,22 @@ QDataStream &operator>>(QDataStream &ds, GroMolType &moltyp) { VersionID v = readHeader(ds, r_gromoltyp); - if (v == 1 or v == 2 or v == 3) + if (v == 4) + { + // v4: all fields in correct order, fixing the v3 write/read misalignment + // where ffield was written but not read back. + SharedDataStream sds(ds); + + sds >> moltyp.nme >> moltyp.warns >> moltyp.atms0 >> moltyp.atms1 >> moltyp.bnds0 >> moltyp.bnds1 >> + moltyp.angs0 >> moltyp.angs1 >> moltyp.dihs0 >> moltyp.dihs1 >> + moltyp.cmaps0 >> moltyp.cmaps1 >> moltyp.first_atoms0 >> moltyp.first_atoms1 >> + moltyp.ffield0 >> moltyp.ffield1 >> moltyp.nexcl0 >> moltyp.nexcl1 >> moltyp.is_perturbable; + } + else if (v == 1 or v == 2 or v == 3) { + // v1/v2/v3: preserved as-is for backward compatibility. + // Note: v3 streams have a write/read misalignment (ffield was written + // but not read back); existing v3 caches will be corrupt. v4 fixes this. SharedDataStream sds(ds); sds >> moltyp.nme >> moltyp.warns >> moltyp.atms0 >> moltyp.atms1 >> moltyp.bnds0 >> moltyp.bnds1 >> @@ -397,7 +411,7 @@ QDataStream &operator>>(QDataStream &ds, GroMolType &moltyp) sds >> moltyp.nexcl0 >> moltyp.nexcl1 >> moltyp.is_perturbable; } else - throw version_error(v, "1,2,3", r_gromoltyp, CODELOC); + throw version_error(v, "1,2,3,4", r_gromoltyp, CODELOC); return ds; } @@ -1617,7 +1631,8 @@ GroMolType::GroMolType(const GroMolType &other) : nme(other.nme), warns(other.warns), atms0(other.atms0), atms1(other.atms1), first_atoms0(other.first_atoms0), first_atoms1(other.first_atoms1), bnds0(other.bnds0), bnds1(other.bnds1), angs0(other.angs0), angs1(other.angs1), dihs0(other.dihs0), dihs1(other.dihs1), cmaps0(other.cmaps0), cmaps1(other.cmaps1), - ffield0(other.ffield0), ffield1(other.ffield1), nexcl0(other.nexcl0), + ffield0(other.ffield0), ffield1(other.ffield1), + explicit_pairs(other.explicit_pairs), nexcl0(other.nexcl0), nexcl1(other.nexcl1), is_perturbable(other.is_perturbable) { } @@ -1646,10 +1661,11 @@ GroMolType &GroMolType::operator=(const GroMolType &other) dihs1 = other.dihs1; cmaps0 = other.cmaps0; cmaps1 = other.cmaps1; + explicit_pairs = other.explicit_pairs; ffield0 = other.ffield0; ffield1 = other.ffield1; nexcl0 = other.nexcl0; - nexcl0 = other.nexcl1; + nexcl1 = other.nexcl1; is_perturbable = other.is_perturbable; } @@ -1663,6 +1679,7 @@ bool GroMolType::operator==(const GroMolType &other) const first_atoms0 == other.first_atoms0 and first_atoms1 == other.first_atoms1 and bnds0 == other.bnds0 and bnds1 == other.bnds1 and angs0 == other.angs0 and angs1 == other.angs1 and dihs0 == other.dihs0 and dihs1 == other.dihs1 and cmaps0 == other.cmaps0 and cmaps1 == other.cmaps1 and + explicit_pairs == other.explicit_pairs and nexcl0 == other.nexcl0 and nexcl1 == other.nexcl1 and is_perturbable == other.is_perturbable; } @@ -2483,6 +2500,19 @@ QHash GroMolType::cmaps(bool is_lambda1) const return cmaps0; } +/** Add an explicit 1-4 pair (from a [pairs] funct=2 line) with given + * coulomb and LJ scale factors */ +void GroMolType::addExplicitPair(const BondID &pair, double cscl, double ljscl) +{ + explicit_pairs.insert(pair, qMakePair(cscl, ljscl)); +} + +/** Return the explicit 1-4 pair scale factors (from [pairs] funct=2) */ +QHash> GroMolType::explicitPairs() const +{ + return explicit_pairs; +} + /** Sanitise all of the CMAP terms - this sets the string equal to "1", * as the information contained previously has already been read */ @@ -2613,7 +2643,8 @@ QStringList GroMolType::settlesLines(bool is_lambda1) const // lambda function to check whether a four point water model // is OPC water, which is determined by the virtual site charge // value being < -1.1 - auto is_opc = [this, is_lambda1]() -> bool { + auto is_opc = [this, is_lambda1]() -> bool + { if (is_lambda1) { for (const auto &atm : atms1) @@ -3699,7 +3730,7 @@ static QStringList writeCMAPTypes(const QHash &cmap_para /** Internal function used to convert a Gromacs Moltyp to a set of lines */ static QStringList writeMolType(const QString &name, const GroMolType &moltype, const Molecule &mol, - bool uses_parallel) + bool uses_parallel, int combining_rules = 2) { QStringList lines; @@ -4605,6 +4636,41 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype, // Store the molinfo object; const auto molinfo = mol.info(); + // Precompute the set of genuine 1-4 bonded atom pairs once. This lets + // us distinguish GLYCAM-style 1-4 pairs (CLJScaleFactor(1,1)) from + // 1-5+ pairs that return the same default value but must not appear in + // [pairs]. Empty if connectivity is unavailable (fallback: write all). + QSet> pairs_14; + + try + { + const auto conn = mol.property("connectivity").asA(); + const int natoms = molinfo.nAtoms(); + for (int a = 0; a < natoms; ++a) + { + const AtomIdx aidx(a); + for (const auto &b : conn.connectionsTo(aidx)) + { + for (const auto &c : conn.connectionsTo(b)) + { + if (c == aidx) + continue; + for (const auto &d : conn.connectionsTo(c)) + { + if (d == b or d == aidx) + continue; + // aidx-b-c-d is a dihedral: aidx and d are 1-4 bonded. + pairs_14.insert(QPair(aidx, d)); + pairs_14.insert(QPair(d, aidx)); + } + } + } + } + } + catch (...) + { + } + if (is_perturbable) { CLJNBPairs scl0; @@ -4630,10 +4696,16 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype, AtomLJs ljs0; AtomLJs ljs1; + AtomCharges charges0; + AtomCharges charges1; + bool has_ljs = false; + bool has_charges = false; try { ljs0 = mol.property("LJ0").asA(); + ljs1 = mol.property("LJ1").asA(); + has_ljs = true; } catch (...) { @@ -4641,102 +4713,117 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype, try { - ljs1 = mol.property("LJ1").asA(); + charges0 = mol.property("charge0").asA(); + charges1 = mol.property("charge1").asA(); + has_charges = true; } catch (...) { } - // A set of recorded 1-4 pairs. - QSet> recorded_pairs; - bool fix_null_perturbable_14s = false; if (mol.hasProperty("fix_null_perturbable_14s")) fix_null_perturbable_14s = mol.property("fix_null_perturbable_14s").asA().value(); - // Must record every pair that has a non-default scaling factor. - // Loop over intrascale matrix by cut-groups to avoid N^2 loop. - for (int i = 0; i < scl0.nGroups(); ++i) + // When connectivity is available, iterate over genuine 1-4 bonded pairs + // — O(N_dihedrals). Otherwise use nonDefaultElements() on each CG pair + // — O(N_bonded). Both avoid the O(N^2) atom-pair loop. + auto write_pair14 = [&](AtomIdx idx0, AtomIdx idx1) { - for (int j = 0; j < scl0.nGroups(); ++j) + const auto s0 = scl0.get(idx0, idx1); + const auto s1 = scl1.get(idx0, idx1); + if (s0.coulomb() == 1 and s0.lj() == 1 and + s1.coulomb() == 1 and s1.lj() == 1) { - const auto s0 = scl0.get(CGIdx(i), CGIdx(j)); - const auto s1 = scl1.get(CGIdx(i), CGIdx(j)); - - if (not s0.isEmpty() and not s1.isEmpty()) + // Both end states have full 1-4 interaction (GLYCAM). Write as + // funct=2 with explicit LJ from state 0. + if (has_ljs and has_charges) { - const auto idxs0 = molinfo.getAtomsIn(CGIdx(i)); - const auto idxs1 = molinfo.getAtomsIn(CGIdx(j)); - - for (const auto &idx0 : idxs0) + const auto cgidx0 = molinfo.cgAtomIdx(idx0); + const auto cgidx1 = molinfo.cgAtomIdx(idx1); + const auto &lj0 = ljs0.at(cgidx0); + const auto &lj1 = ljs0.at(cgidx1); + LJParameter lj_ij; + if (combining_rules == 2) + lj_ij = lj0.combineArithmetic(lj1); + else + lj_ij = lj0.combineGeometric(lj1); + const double qi = charges0.at(cgidx0).to(mod_electron); + const double qj = charges0.at(cgidx1).to(mod_electron); + scllines.append( + QString("%1 %2 2 1.0 %3 %4 %5 %6") + .arg(idx0 + 1, 6) + .arg(idx1 + 1, 6) + .arg(qi, 11, 'f', 6) + .arg(qj, 11, 'f', 6) + .arg(lj_ij.sigma().to(nanometer), 18, 'e', 11) + .arg(lj_ij.epsilon().to(kJ_per_mol), 18, 'e', 11)); + } + else + { + scllines.append(QString("%1 %2 1") + .arg(idx0 + 1, 6) + .arg(idx1 + 1, 6)); + } + } + else if (not(s0.coulomb() == 0 and s0.lj() == 0 and + s1.coulomb() == 0 and s1.lj() == 0)) + { + // Standard partial 1-4 scaling, or a mixed perturbation. + if (fix_null_perturbable_14s) + { + const auto &lj0_0 = ljs0.get(idx0); + const auto &lj0_1 = ljs1.get(idx0); + const auto &lj1_0 = ljs0.get(idx1); + const auto &lj1_1 = ljs1.get(idx1); + if (lj0_0.epsilon().value() == 0 or lj0_1.epsilon().value() == 0 or + lj1_0.epsilon().value() == 0 or lj1_1.epsilon().value() == 0) { - for (const auto &idx1 : idxs1) - { - QPair pair = QPair(idx0, idx1); + LJParameter lj0, lj1; + lj0 = (lj0_0.epsilon().value() == 0) ? lj0_1 : lj0_0; + lj1 = (lj1_0.epsilon().value() == 0) ? lj1_1 : lj1_0; + auto lj = (combining_rules == 2) ? lj0.combineArithmetic(lj1) + : lj0.combineGeometric(lj1); + double scl = (s0.lj() != 0) ? s0.lj() : s1.lj(); + scllines.append( + QString("%1 %2 1 %3 %4 %3 %4") + .arg(idx0 + 1, 6) + .arg(idx1 + 1, 6) + .arg(lj.sigma().to(nanometer), 11, 'f', 5) + .arg(scl * lj.epsilon().to(kJ_per_mol), 11, 'f', 5)); + return; + } + } + scllines.append(QString("%1 %2 1") + .arg(idx0 + 1, 6) + .arg(idx1 + 1, 6)); + } + }; - // Make sure this is a new atom pair. - if (not recorded_pairs.contains(pair)) - { - // Insert the pair and its inverse. - recorded_pairs.insert(pair); - pair = QPair(idx1, idx0); - recorded_pairs.insert(pair); - - const auto s0 = scl0.get(idx0, idx1); - const auto s1 = scl1.get(idx0, idx1); - - if (not((s0.coulomb() == 0 and s0.lj() == 0 and s1.coulomb() == 0 and - s1.lj() == 0) or - (s0.coulomb() == 1 and s0.lj() == 1 and s1.coulomb() == 1 and - s1.lj() == 1))) - { - // This is a non-default pair. - if (fix_null_perturbable_14s) - { - // get the initial and perturbed charge and LJ parameters - const auto &lj0_0 = ljs0.get(idx0); - const auto &lj0_1 = ljs1.get(idx0); - const auto &lj1_0 = ljs0.get(idx1); - const auto &lj1_1 = ljs1.get(idx1); - - if (lj0_0.epsilon().value() == 0 or lj0_1.epsilon().value() == 0 or - lj1_0.epsilon().value() == 0 or lj1_1.epsilon().value() == 0) - { - // we need to avoid having a null 1-4 LJ parameter, so use the non-dummy state - LJParameter lj0, lj1; - - if (lj0_0.epsilon().value() == 0) - lj0 = lj0_1; - else - lj0 = lj0_0; - - if (lj1_0.epsilon().value() == 0) - lj1 = lj1_1; - else - lj1 = lj1_0; - - auto lj = lj0.combineArithmetic(lj1); - - double scl = s0.lj(); - - if (scl == 0) - scl = s1.lj(); - - scllines.append(QString("%1 %2 1 %3 %4 %3 %4") - .arg(idx0 + 1, 6) - .arg(idx1 + 1, 6) - .arg(lj.sigma().to(nanometer), 11, 'f', 5) - .arg(scl * lj.epsilon().to(kJ_per_mol), 11, 'f', 5)); - - continue; - } - } - - scllines.append(QString("%1 %2 1").arg(idx0 + 1, 6).arg(idx1 + 1, 6)); - } - } - } + if (not pairs_14.isEmpty()) + { + for (const auto &pair14 : pairs_14) + { + if (pair14.first >= pair14.second) + continue; + write_pair14(pair14.first, pair14.second); + } + } + else + { + // No connectivity: iterate over non-default CG atom pair entries. + // GLYCAM-style (1,1) pairs cannot be identified here and are omitted. + for (int i = 0; i < scl0.nGroups(); ++i) + { + for (int j = 0; j < scl0.nGroups(); ++j) + { + for (const auto &[row, col, s0] : + scl0.get(CGIdx(i), CGIdx(j)).nonDefaultElements()) + { + if (row >= col) + continue; + write_pair14(AtomIdx(row), AtomIdx(col)); } } } @@ -4755,45 +4842,106 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype, return; } - // A set of recorded 1-4 pairs. - QSet> recorded_pairs; + // Get LJ and charge properties for writing funct=2 explicit pairs. + AtomLJs ljs; + AtomCharges charges; + bool has_ljs = false; + bool has_charges = false; - // Must record every pair that has a non-default scaling factor. - // Loop over intrascale matrix by cut-groups to avoid N^2 loop. - for (int i = 0; i < scl.nGroups(); ++i) + try { - for (int j = 0; j < scl.nGroups(); ++j) - { - const auto s = scl.get(CGIdx(i), CGIdx(j)); + ljs = mol.property("LJ").asA(); + has_ljs = true; + } + catch (...) + { + } - if (not s.isEmpty()) + try + { + charges = mol.property("charge").asA(); + has_charges = true; + } + catch (...) + { + } + + // When connectivity is available, iterate over genuine 1-4 bonded pairs + // — O(N_dihedrals). Otherwise use nonDefaultElements() on each CG pair + // — O(N_bonded). Both avoid the O(N^2) atom-pair loop. + auto write_pair14 = [&](AtomIdx idx0, AtomIdx idx1) + { + const auto s = scl.get(idx0, idx1); + if (s.coulomb() == 0 and s.lj() == 0) + { + // excluded — skip + } + else if (s.coulomb() == 1 and s.lj() == 1) + { + // Full 1-4 interaction (GLYCAM). Write as funct=2 with explicit LJ + // parameters because funct=1 would apply fudgeLJ scaling. + if (has_ljs and has_charges) + { + const auto cgidx0 = molinfo.cgAtomIdx(idx0); + const auto cgidx1 = molinfo.cgAtomIdx(idx1); + const auto &lj0 = ljs.at(cgidx0); + const auto &lj1 = ljs.at(cgidx1); + LJParameter lj_ij; + if (combining_rules == 2) + lj_ij = lj0.combineArithmetic(lj1); + else + lj_ij = lj0.combineGeometric(lj1); + const double qi = charges.at(cgidx0).to(mod_electron); + const double qj = charges.at(cgidx1).to(mod_electron); + scllines.append( + QString("%1 %2 2 1.0 %3 %4 %5 %6") + .arg(idx0 + 1, 6) + .arg(idx1 + 1, 6) + .arg(qi, 11, 'f', 6) + .arg(qj, 11, 'f', 6) + .arg(lj_ij.sigma().to(nanometer), 18, 'e', 11) + .arg(lj_ij.epsilon().to(kJ_per_mol), 18, 'e', 11)); + } + else { - const auto idxs0 = molinfo.getAtomsIn(CGIdx(i)); - const auto idxs1 = molinfo.getAtomsIn(CGIdx(j)); + // Fall back to funct=1; energy will be wrong if fudgeLJ != 1.0. + scllines.append(QString("%1 %2 1") + .arg(idx0 + 1, 6) + .arg(idx1 + 1, 6)); + } + } + else + { + // Standard partial 1-4 scaling (e.g. AMBER). Write as funct=1. + scllines.append(QString("%1 %2 1") + .arg(idx0 + 1, 6) + .arg(idx1 + 1, 6)); + } + }; - for (const auto &idx0 : idxs0) + if (not pairs_14.isEmpty()) + { + for (const auto &pair14 : pairs_14) + { + if (pair14.first >= pair14.second) + continue; + write_pair14(pair14.first, pair14.second); + } + } + else + { + // No connectivity: iterate over non-default CG atom pair entries. + // GLYCAM-style (1,1) pairs cannot be identified here and are omitted. + for (int i = 0; i < scl.nGroups(); ++i) + { + for (int j = 0; j < scl.nGroups(); ++j) + { + for (const auto &[row, col, s] : + scl.get(CGIdx(i), CGIdx(j)).nonDefaultElements()) { - for (const auto &idx1 : idxs1) - { - QPair pair = QPair(idx0, idx1); - - // Make sure this is a new atom pair. - if (not recorded_pairs.contains(pair)) - { - // Insert the pair and its inverse. - recorded_pairs.insert(pair); - pair = QPair(idx1, idx0); - recorded_pairs.insert(pair); - - const auto s = scl.get(idx0, idx1); - - if (not((s.coulomb() == 0 and s.lj() == 0) or (s.coulomb() == 1 and s.lj() == 1))) - { - // This is a non-default pair. - scllines.append(QString("%1 %2 1").arg(idx0 + 1, 6).arg(idx1 + 1, 6)); - } - } - } + if (row >= col) + continue; + write_pair14(AtomIdx(row), AtomIdx(col)); } } } @@ -7836,6 +7984,72 @@ QStringList GroTop::processDirectives(const QMap &taglocs, const Q }; // function that extracts all of the information from the 'cmap' lines + // function that extracts explicit 1-4 pair scale factors from the 'pairs' lines. + // funct=1 pairs are standard (use global fudge_qq/fudge_lj) and are handled + // automatically by gen-pairs, so we only need to store funct=2 explicit pairs. + // funct=2 format: ai aj 2 fudgeQQ qi qj sigma epsilon + // The LJ scale is 1.0 for funct=2 because sigma/epsilon are the full combined values. + auto addPairsTo = [&](GroMolType &moltype, int linenum) + { + QStringList lines = getDirectiveLines(linenum); + + for (const auto &line : lines) + { + const auto words = line.split(" "); + + if (words.count() < 3) + { + moltype.addWarning(QObject::tr("Cannot extract pair information " + "from the line '%1' as it should contain at least three words.") + .arg(line)); + continue; + } + + bool ok0, ok1, ok2; + + int atm0 = words[0].toInt(&ok0); + int atm1 = words[1].toInt(&ok1); + int funct = words[2].toInt(&ok2); + + if (not(ok0 and ok1 and ok2)) + { + moltype.addWarning(QObject::tr("Cannot extract pair information " + "from the line '%1' as the first three words need to be integers.") + .arg(line)); + continue; + } + + if (funct == 1) + { + // Standard pair: uses global fudge_qq/fudge_lj. + // The gen-pairs mechanism already handles these, so no explicit storage needed. + continue; + } + else if (funct == 2) + { + // Explicit pair: ai aj 2 fudgeQQ qi qj sigma epsilon + // The fudgeQQ is the coulomb scale factor; LJ params are used directly (lj_scl = 1.0). + double cscl = fudge_qq; // default to global fudge_qq if not specified + if (words.count() > 3) + { + bool ok; + double val = words[3].toDouble(&ok); + if (ok) + cscl = val; + } + + moltype.addExplicitPair(BondID(AtomNum(atm0), AtomNum(atm1)), cscl, 1.0); + } + else + { + moltype.addWarning(QObject::tr("Unsupported pair function type %1 in line '%2'. " + "Only funct=1 and funct=2 are supported.") + .arg(funct) + .arg(line)); + } + } + }; + auto addCMAPsTo = [&](GroMolType &moltype, int linenum) { QStringList lines = getDirectiveLines(linenum); @@ -7917,9 +8131,13 @@ QStringList GroTop::processDirectives(const QMap &taglocs, const Q addCMAPsTo(moltype, linenum); } + for (auto linenum : moltag.values("pairs")) + { + addPairsTo(moltype, linenum); + } + // now print out warnings for any lines that are missed... - const QStringList missed_tags = {"pairs", - "pairs_nb", + const QStringList missed_tags = {"pairs_nb", "exclusions", "contraints", "settles", @@ -8327,6 +8545,33 @@ Molecule GroTop::createMolecule(const GroMolType &moltype, QStringList &errors, ChainStructureEditor chain; CGStructureEditor cgroup; + // 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 within a molecule cause duplicate_residue errors. + // next_unique is always > every number assigned so far, so conflict + // resolution is O(1) rather than a linear scan. + QSet 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 (topology) residue number of the current residue + // separately, because unique_resnum may assign a different number. The + // "is this a new residue?" check must use the original topology number. + int current_orig_resnum = -1; + ResName current_resname; + auto different_chain = [&](const ChainName &name) { if (name.isNull() and chain.isEmpty()) @@ -8348,14 +8593,16 @@ Molecule GroTop::createMolecule(const GroMolType &moltype, QStringList &errors, if (not atom.chainName().isNull()) { chain = mol.add(ChainName(atom.chainName())); - res = chain.add(ResNum(atom.residueNumber())); + res = chain.add(unique_resnum(atom.residueNumber())); } else { - res = mol.add(ResNum(atom.residueNumber())); + res = mol.add(unique_resnum(atom.residueNumber())); } res = res.rename(atom.residueName()); + current_orig_resnum = atom.residueNumber(); + current_resname = atom.residueName(); } else if (different_chain(atom.chainName())) { @@ -8367,25 +8614,30 @@ Molecule GroTop::createMolecule(const GroMolType &moltype, QStringList &errors, { // residue is not in a chain chain = ChainStructureEditor(); - res = mol.add(ResNum(atom.residueNumber())); + res = mol.add(unique_resnum(atom.residueNumber())); } else { // residue is in a chain chain = mol.add(ChainName(atom.chainName())); - res = chain.add(ResNum(atom.residueNumber())); + res = chain.add(unique_resnum(atom.residueNumber())); } res = res.rename(atom.residueName()); + current_orig_resnum = atom.residueNumber(); + current_resname = atom.residueName(); } - else if (atom.residueNumber() != res.number() or atom.residueName() != res.name()) + else if (atom.residueNumber() != current_orig_resnum or + atom.residueName() != current_resname) { // this atom is in a different residue cgroup = mol.add(CGName(QString::number(cgidx))); cgidx += 1; - res = mol.add(ResNum(atom.residueNumber())); + res = mol.add(unique_resnum(atom.residueNumber())); res = res.rename(atom.residueName()); + current_orig_resnum = atom.residueNumber(); + current_resname = atom.residueName(); } // add the atom to the residue @@ -8618,6 +8870,27 @@ GroTop::PropsAndErrors GroTop::getBondProperties(const MoleculeInfo &molinfo, co else { CLJNBPairs nbpairs(conn, CLJScaleFactor(fudge_qq, fudge_lj)); + + // Override with any explicitly specified [pairs] funct=2 entries. + // These carry their own fudgeQQ (coulomb scale) and use lj_scl=1.0 + // (sigma/epsilon in funct=2 are the full combined values, not scaled by fudgeLJ). + const auto explicit_pairs = moltype.explicitPairs(); + for (auto it = explicit_pairs.constBegin(); it != explicit_pairs.constEnd(); ++it) + { + const auto &pair = it.key(); + const auto &scl = it.value(); + try + { + AtomIdx idx0 = molinfo.atomIdx(pair.atom0()); + AtomIdx idx1 = molinfo.atomIdx(pair.atom1()); + nbpairs.set(idx0, idx1, CLJScaleFactor(scl.first, scl.second)); + } + catch (...) + { + // atom not found — skip silently (already warned during parsing) + } + } + props.setProperty("intrascale", nbpairs); } } diff --git a/corelib/src/libs/SireIO/grotop.h b/corelib/src/libs/SireIO/grotop.h index 2cbf2c41a..679de9182 100644 --- a/corelib/src/libs/SireIO/grotop.h +++ b/corelib/src/libs/SireIO/grotop.h @@ -43,6 +43,7 @@ #include "SireMM/cmapparameter.h" #include +#include SIRE_BEGIN_HEADER @@ -247,6 +248,9 @@ namespace SireIO QMultiHash dihedrals(bool is_lambda1 = false) const; QHash cmaps(bool is_lambda1 = false) const; + void addExplicitPair(const SireMol::BondID &pair, double cscl, double ljscl); + QHash> explicitPairs() const; + bool isWater(bool is_lambda1 = false) const; QStringList settlesLines(bool is_lambda1 = false) const; @@ -296,6 +300,12 @@ namespace SireIO SireMM::MMDetail ffield0; SireMM::MMDetail ffield1; + /** Explicit 1-4 pair scale factors from [pairs] funct=2 lines. + * Key: the atom pair. Value: (coulomb_scl, lj_scl) where + * coulomb_scl comes from the fudgeQQ column and lj_scl = 1.0 + * (since funct=2 LJ parameters are used directly, not further scaled). */ + QHash> explicit_pairs; + /** The number of excluded atoms */ qint64 nexcl0; qint64 nexcl1; diff --git a/corelib/src/libs/SireMM/amberparams.cpp b/corelib/src/libs/SireMM/amberparams.cpp index b102966f4..5312666f6 100644 --- a/corelib/src/libs/SireMM/amberparams.cpp +++ b/corelib/src/libs/SireMM/amberparams.cpp @@ -103,9 +103,7 @@ QDataStream &operator>>(QDataStream &ds, AmberBond &bond) } /** Construct with the passed bond constant and equilibrium bond length */ -AmberBond::AmberBond(double k, double r0) : _k(k), _r0(r0) -{ -} +AmberBond::AmberBond(double k, double r0) : _k(k), _r0(r0) {} /** Construct from the passed expression */ AmberBond::AmberBond(const Expression &f, const Symbol &R) : _k(0), _r0(0) @@ -136,15 +134,17 @@ AmberBond::AmberBond(const Expression &f, const Symbol &R) : _k(0), _r0(0) { if (not factor.power().isConstant()) { - errors.append(QObject::tr("Power of R must be constant, not %1").arg(factor.power().toString())); + errors.append(QObject::tr("Power of R must be constant, not %1") + .arg(factor.power().toString())); continue; } if (not factor.factor().isConstant()) { - errors.append(QObject::tr("The value of K in K (R - R0)^2 must be constant. " - "Here it is %1") - .arg(factor.factor().toString())); + errors.append( + QObject::tr("The value of K in K (R - R0)^2 must be constant. " + "Here it is %1") + .arg(factor.factor().toString())); continue; } @@ -175,14 +175,17 @@ AmberBond::AmberBond(const Expression &f, const Symbol &R) : _k(0), _r0(0) } else { - errors.append(QObject::tr("Power of R^2 must equal 2.0, 1.0 or 0.0, not %1").arg(power)); + errors.append( + QObject::tr("Power of R^2 must equal 2.0, 1.0 or 0.0, not %1") + .arg(power)); continue; } } else { errors.append( - QObject::tr("Cannot have a factor that does not include R. %1").arg(factor.symbol().toString())); + QObject::tr("Cannot have a factor that does not include R. %1") + .arg(factor.symbol().toString())); } } @@ -192,17 +195,20 @@ AmberBond::AmberBond(const Expression &f, const Symbol &R) : _k(0), _r0(0) // kr0 should be equal to -2 k r0 if (std::abs(_k * _r0 + 0.5 * kr0) > 0.001) { - errors.append(QObject::tr("How can the power of R be %1. It should be 2 x %2 x %3 = %4.") - .arg(kr0) - .arg(_k) - .arg(_r0) - .arg(2 * _k * _r0)); + errors.append( + QObject::tr( + "How can the power of R be %1. It should be 2 x %2 x %3 = %4.") + .arg(kr0) + .arg(_k) + .arg(_r0) + .arg(2 * _k * _r0)); } if (not errors.isEmpty()) { throw SireError::incompatible_error( - QObject::tr("Cannot extract an AmberBond with function K ( %1 - R0 )^2 from the " + QObject::tr("Cannot extract an AmberBond with function K ( %1 - R0 )^2 " + "from the " "expression %2, because\n%3") .arg(R.toString()) .arg(f.toString()) @@ -211,13 +217,9 @@ AmberBond::AmberBond(const Expression &f, const Symbol &R) : _k(0), _r0(0) } } -AmberBond::AmberBond(const AmberBond &other) : _k(other._k), _r0(other._r0) -{ -} +AmberBond::AmberBond(const AmberBond &other) : _k(other._k), _r0(other._r0) {} -AmberBond::~AmberBond() -{ -} +AmberBond::~AmberBond() {} double AmberBond::operator[](int i) const { @@ -271,10 +273,7 @@ const char *AmberBond::typeName() return QMetaType::typeName(qMetaTypeId()); } -const char *AmberBond::what() const -{ - return AmberBond::typeName(); -} +const char *AmberBond::what() const { return AmberBond::typeName(); } /** Return the energy evaluated from this bond for the passed bond length */ double AmberBond::energy(double r) const @@ -321,11 +320,10 @@ QDataStream &operator>>(QDataStream &ds, AmberAngle &angle) return ds; } -AmberAngle::AmberAngle(double k, double theta0) : _k(k), _theta0(theta0) -{ -} +AmberAngle::AmberAngle(double k, double theta0) : _k(k), _theta0(theta0) {} -AmberAngle::AmberAngle(const Expression &f, const Symbol &theta) : _k(0), _theta0(0) +AmberAngle::AmberAngle(const Expression &f, const Symbol &theta) + : _k(0), _theta0(0) { if (f.isZero()) { @@ -353,15 +351,17 @@ AmberAngle::AmberAngle(const Expression &f, const Symbol &theta) : _k(0), _theta { if (not factor.power().isConstant()) { - errors.append(QObject::tr("Power of theta must be constant, not %1").arg(factor.power().toString())); + errors.append(QObject::tr("Power of theta must be constant, not %1") + .arg(factor.power().toString())); continue; } if (not factor.factor().isConstant()) { - errors.append(QObject::tr("The value of K in K (theta - theta0)^2 must be " - "constant. Here it is %1") - .arg(factor.factor().toString())); + errors.append( + QObject::tr("The value of K in K (theta - theta0)^2 must be " + "constant. Here it is %1") + .arg(factor.factor().toString())); continue; } @@ -392,14 +392,17 @@ AmberAngle::AmberAngle(const Expression &f, const Symbol &theta) : _k(0), _theta } else { - errors.append(QObject::tr("Power of theta^2 must equal 2.0, 1.0 or 0.0, not %1").arg(power)); + errors.append( + QObject::tr("Power of theta^2 must equal 2.0, 1.0 or 0.0, not %1") + .arg(power)); continue; } } else { errors.append( - QObject::tr("Cannot have a factor that does not include theta. %1").arg(factor.symbol().toString())); + QObject::tr("Cannot have a factor that does not include theta. %1") + .arg(factor.symbol().toString())); } } @@ -409,17 +412,20 @@ AmberAngle::AmberAngle(const Expression &f, const Symbol &theta) : _k(0), _theta // ktheta0 should be equal to -k theta0 if (std::abs(_k * _theta0 + 0.5 * ktheta0) > 0.001) { - errors.append(QObject::tr("How can the power of theta be %1. It should be 2 x %2 x %3 = %4.") - .arg(ktheta0) - .arg(_k) - .arg(_theta0) - .arg(2 * _k * _theta0)); + errors.append( + QObject::tr( + "How can the power of theta be %1. It should be 2 x %2 x %3 = %4.") + .arg(ktheta0) + .arg(_k) + .arg(_theta0) + .arg(2 * _k * _theta0)); } if (not errors.isEmpty()) { throw SireError::incompatible_error( - QObject::tr("Cannot extract an AmberAngle with function K ( %1 - theta0 )^2 from the " + QObject::tr("Cannot extract an AmberAngle with function K ( %1 - " + "theta0 )^2 from the " "expression %2, because\n%3") .arg(theta.toString()) .arg(f.toString()) @@ -428,13 +434,10 @@ AmberAngle::AmberAngle(const Expression &f, const Symbol &theta) : _k(0), _theta } } -AmberAngle::AmberAngle(const AmberAngle &other) : _k(other._k), _theta0(other._theta0) -{ -} +AmberAngle::AmberAngle(const AmberAngle &other) + : _k(other._k), _theta0(other._theta0) {} -AmberAngle::~AmberAngle() -{ -} +AmberAngle::~AmberAngle() {} double AmberAngle::operator[](int i) const { @@ -486,10 +489,7 @@ const char *AmberAngle::typeName() return QMetaType::typeName(qMetaTypeId()); } -const char *AmberAngle::what() const -{ - return AmberAngle::typeName(); -} +const char *AmberAngle::what() const { return AmberAngle::typeName(); } double AmberAngle::energy(double theta) const { @@ -533,18 +533,13 @@ QDataStream &operator>>(QDataStream &ds, AmberDihPart &dih) return ds; } -AmberDihPart::AmberDihPart(double k, double periodicity, double phase) : _k(k), _periodicity(periodicity), _phase(phase) -{ -} +AmberDihPart::AmberDihPart(double k, double periodicity, double phase) + : _k(k), _periodicity(periodicity), _phase(phase) {} AmberDihPart::AmberDihPart(const AmberDihPart &other) - : _k(other._k), _periodicity(other._periodicity), _phase(other._phase) -{ -} + : _k(other._k), _periodicity(other._periodicity), _phase(other._phase) {} -AmberDihPart::AmberDihPart::~AmberDihPart() -{ -} +AmberDihPart::AmberDihPart::~AmberDihPart() {} double AmberDihPart::operator[](int i) const { @@ -568,7 +563,8 @@ AmberDihPart &AmberDihPart::operator=(const AmberDihPart &other) bool AmberDihPart::operator==(const AmberDihPart &other) const { - return _k == other._k and _periodicity == other._periodicity and _phase == other._phase; + return _k == other._k and _periodicity == other._periodicity and + _phase == other._phase; } bool AmberDihPart::operator!=(const AmberDihPart &other) const @@ -599,10 +595,7 @@ const char *AmberDihPart::typeName() return QMetaType::typeName(qMetaTypeId()); } -const char *AmberDihPart::what() const -{ - return AmberDihPart::typeName(); -} +const char *AmberDihPart::what() const { return AmberDihPart::typeName(); } double AmberDihPart::energy(double phi) const { @@ -611,7 +604,10 @@ double AmberDihPart::energy(double phi) const QString AmberDihPart::toString() const { - return QObject::tr("AmberDihPart( k = %1, periodicity = %2, phase = %3 )").arg(_k).arg(_periodicity).arg(_phase); + return QObject::tr("AmberDihPart( k = %1, periodicity = %2, phase = %3 )") + .arg(_k) + .arg(_periodicity) + .arg(_phase); } /////////// @@ -641,16 +637,15 @@ QDataStream &operator>>(QDataStream &ds, AmberDihedral &dih) return ds; } -AmberDihedral::AmberDihedral() -{ -} +AmberDihedral::AmberDihedral() {} AmberDihedral::AmberDihedral(AmberDihPart part) { _parts = QVector(1, part); } -AmberDihedral::AmberDihedral(const Expression &f, const Symbol &phi, bool test_ryckaert_bellemans) +AmberDihedral::AmberDihedral(const Expression &f, const Symbol &phi, + bool test_ryckaert_bellemans) { if (f.isZero()) { @@ -686,8 +681,9 @@ AmberDihedral::AmberDihedral(const Expression &f, const Symbol &phi, bool test_r auto F1 = nrg_factor * (3.0 * F3 - 2.0 * params[1]); // Convert the expression to a Fourier series. - f_copy = Expression(0.5 * (F1 * (1 + Cos(phi)) + F2 * (1 - Cos(2 * phi)) + F3 * (1 + Cos(3 * phi)) + - F4 * (1 - Cos(4 * phi)))); + f_copy = Expression( + 0.5 * (F1 * (1 + Cos(phi)) + F2 * (1 - Cos(2 * phi)) + + F3 * (1 + Cos(3 * phi)) + F4 * (1 - Cos(4 * phi)))); } } catch (...) @@ -800,16 +796,18 @@ AmberDihedral::AmberDihedral(const Expression &f, const Symbol &phi, bool test_r auto cos_term = neg_terms[i].second.base().asA(); // Store the negated factor and shifted cosine term. - neg_terms[i] = QPair(factors[i], Cos(cos_term.argument() + SireMaths::pi)); + neg_terms[i] = QPair( + factors[i], Cos(cos_term.argument() + SireMaths::pi)); } } } else { throw SireError::incompatible_error( - QObject::tr("Cannot extract an Amber-format dihedral expression from '%1' as " - "the expression must be a series of terms of type " - "'k{ 1 + cos[ per %2 - phase ] }'. Errors include\n%3") + QObject::tr( + "Cannot extract an Amber-format dihedral expression from '%1' as " + "the expression must be a series of terms of type " + "'k{ 1 + cos[ per %2 - phase ] }'. Errors include\n%3") .arg(f.toString()) .arg(phi.toString()) .arg(errors.join("\n")), @@ -847,7 +845,8 @@ AmberDihedral::AmberDihedral(const Expression &f, const Symbol &phi, bool test_r { if (not factor.power().isConstant()) { - errors.append(QObject::tr("Power of phi must be constant, not %1").arg(factor.power().toString())); + errors.append(QObject::tr("Power of phi must be constant, not %1") + .arg(factor.power().toString())); ok = false; continue; } @@ -875,7 +874,8 @@ AmberDihedral::AmberDihedral(const Expression &f, const Symbol &phi, bool test_r } else { - errors.append(QObject::tr("Power of phi must equal 1.0 or 0.0, not %1").arg(power)); + errors.append(QObject::tr("Power of phi must equal 1.0 or 0.0, not %1") + .arg(power)); ok = false; continue; } @@ -890,9 +890,10 @@ AmberDihedral::AmberDihedral(const Expression &f, const Symbol &phi, bool test_r if (not errors.isEmpty()) { throw SireError::incompatible_error( - QObject::tr("Cannot extract an Amber-format dihedral expression from '%1' as " - "the expression must be a series of terms of type " - "'k{ 1 + cos[ per %2 - phase ] }'. Errors include\n%3") + QObject::tr( + "Cannot extract an Amber-format dihedral expression from '%1' as " + "the expression must be a series of terms of type " + "'k{ 1 + cos[ per %2 - phase ] }'. Errors include\n%3") .arg(f.toString()) .arg(phi.toString()) .arg(errors.join("\n")), @@ -907,18 +908,16 @@ AmberDihedral::AmberDihedral(const Expression &f, const Symbol &phi, bool test_r for (const auto &term : terms) { // Remember that the expression uses the negative of the phase ;-) - _parts.append(AmberDihPart(std::get<0>(term), std::get<1>(term), -std::get<2>(term))); + _parts.append(AmberDihPart(std::get<0>(term), std::get<1>(term), + -std::get<2>(term))); } } } -AmberDihedral::AmberDihedral(const AmberDihedral &other) : _parts(other._parts) -{ -} +AmberDihedral::AmberDihedral(const AmberDihedral &other) + : _parts(other._parts) {} -AmberDihedral::~AmberDihedral() -{ -} +AmberDihedral::~AmberDihedral() {} AmberDihedral &AmberDihedral::operator+=(const AmberDihPart &part) { @@ -954,10 +953,7 @@ const char *AmberDihedral::typeName() return QMetaType::typeName(qMetaTypeId()); } -const char *AmberDihedral::what() const -{ - return AmberDihedral::typeName(); -} +const char *AmberDihedral::what() const { return AmberDihedral::typeName(); } AmberDihPart AmberDihedral::operator[](int i) const { @@ -1043,17 +1039,12 @@ QDataStream &operator>>(QDataStream &ds, AmberNB14 &nb) return ds; } -AmberNB14::AmberNB14(double cscl, double ljscl) : _cscl(cscl), _ljscl(ljscl) -{ -} +AmberNB14::AmberNB14(double cscl, double ljscl) : _cscl(cscl), _ljscl(ljscl) {} -AmberNB14::AmberNB14(const AmberNB14 &other) : _cscl(other._cscl), _ljscl(other._ljscl) -{ -} +AmberNB14::AmberNB14(const AmberNB14 &other) + : _cscl(other._cscl), _ljscl(other._ljscl) {} -AmberNB14::~AmberNB14() -{ -} +AmberNB14::~AmberNB14() {} double AmberNB14::operator[](int i) const { @@ -1124,14 +1115,13 @@ const char *AmberNB14::typeName() return QMetaType::typeName(qMetaTypeId()); } -const char *AmberNB14::what() const -{ - return AmberNB14::typeName(); -} +const char *AmberNB14::what() const { return AmberNB14::typeName(); } QString AmberNB14::toString() const { - return QObject::tr("AmberNB14( cscl = %1, ljscl = %2 )").arg(_cscl).arg(_ljscl); + return QObject::tr("AmberNB14( cscl = %1, ljscl = %2 )") + .arg(_cscl) + .arg(_ljscl); } /** Return the value converted to a CLJScaleFactor */ @@ -1168,24 +1158,19 @@ QDataStream &operator>>(QDataStream &ds, AmberNBDihPart ¶m) } /** Constructor */ -AmberNBDihPart::AmberNBDihPart() -{ -} +AmberNBDihPart::AmberNBDihPart() {} /** Construct from the passed parameters */ -AmberNBDihPart::AmberNBDihPart(const AmberDihPart &dihedral, const AmberNB14 &nb14) : dih(dihedral), nbscl(nb14) -{ -} +AmberNBDihPart::AmberNBDihPart(const AmberDihPart &dihedral, + const AmberNB14 &nb14) + : dih(dihedral), nbscl(nb14) {} /** Copy constructor */ -AmberNBDihPart::AmberNBDihPart(const AmberNBDihPart &other) : dih(other.dih), nbscl(other.nbscl) -{ -} +AmberNBDihPart::AmberNBDihPart(const AmberNBDihPart &other) + : dih(other.dih), nbscl(other.nbscl) {} /** Destructor */ -AmberNBDihPart::~AmberNBDihPart() -{ -} +AmberNBDihPart::~AmberNBDihPart() {} /** Copy assignment operator */ AmberNBDihPart &AmberNBDihPart::operator=(const AmberNBDihPart &other) @@ -1247,14 +1232,13 @@ const char *AmberNBDihPart::typeName() return QMetaType::typeName(qMetaTypeId()); } -const char *AmberNBDihPart::what() const -{ - return AmberNBDihPart::typeName(); -} +const char *AmberNBDihPart::what() const { return AmberNBDihPart::typeName(); } QString AmberNBDihPart::toString() const { - return QObject::tr("AmberNBDihPart( param == %1, scl == %2 )").arg(dih.toString()).arg(nbscl.toString()); + return QObject::tr("AmberNBDihPart( param == %1, scl == %2 )") + .arg(dih.toString()) + .arg(nbscl.toString()); } /////////// @@ -1270,12 +1254,15 @@ QDataStream &operator<<(QDataStream &ds, const AmberParams &amberparam) SharedDataStream sds(ds); - sds << amberparam.molinfo << amberparam.amber_charges << amberparam.amber_ljs << amberparam.amber_masses - << amberparam.amber_elements << amberparam.amber_types << amberparam.born_radii << amberparam.amber_screens - << amberparam.amber_treechains << amberparam.exc_atoms << amberparam.amber_bonds << amberparam.amber_angles - << amberparam.amber_dihedrals << amberparam.amber_impropers << amberparam.amber_nb14s - << amberparam.cmap_funcs << amberparam.radius_set - << amberparam.propmap << static_cast(amberparam); + sds << amberparam.molinfo << amberparam.amber_charges << amberparam.amber_ljs + << amberparam.amber_masses << amberparam.amber_elements + << amberparam.amber_types << amberparam.born_radii + << amberparam.amber_screens << amberparam.amber_treechains + << amberparam.exc_atoms << amberparam.amber_bonds + << amberparam.amber_angles << amberparam.amber_dihedrals + << amberparam.amber_impropers << amberparam.amber_nb14s + << amberparam.cmap_funcs << amberparam.radius_set << amberparam.propmap + << static_cast(amberparam); return ds; } @@ -1289,12 +1276,16 @@ QDataStream &operator>>(QDataStream &ds, AmberParams &amberparam) { SharedDataStream sds(ds); - sds >> amberparam.molinfo >> amberparam.amber_charges >> amberparam.amber_ljs >> amberparam.amber_masses >> - amberparam.amber_elements >> amberparam.amber_types >> amberparam.born_radii >> amberparam.amber_screens >> - amberparam.amber_treechains >> amberparam.exc_atoms >> amberparam.amber_bonds >> amberparam.amber_angles >> - amberparam.amber_dihedrals >> amberparam.amber_impropers >> amberparam.amber_nb14s >> - amberparam.cmap_funcs >> - amberparam.radius_set >> amberparam.propmap >> static_cast(amberparam); + sds >> amberparam.molinfo >> amberparam.amber_charges >> + amberparam.amber_ljs >> amberparam.amber_masses >> + amberparam.amber_elements >> amberparam.amber_types >> + amberparam.born_radii >> amberparam.amber_screens >> + amberparam.amber_treechains >> amberparam.exc_atoms >> + amberparam.amber_bonds >> amberparam.amber_angles >> + amberparam.amber_dihedrals >> amberparam.amber_impropers >> + amberparam.amber_nb14s >> amberparam.cmap_funcs >> + amberparam.radius_set >> amberparam.propmap >> + static_cast(amberparam); } else if (v == 2) { @@ -1302,11 +1293,15 @@ QDataStream &operator>>(QDataStream &ds, AmberParams &amberparam) amberparam.cmap_funcs.clear(); - sds >> amberparam.molinfo >> amberparam.amber_charges >> amberparam.amber_ljs >> amberparam.amber_masses >> - amberparam.amber_elements >> amberparam.amber_types >> amberparam.born_radii >> amberparam.amber_screens >> - amberparam.amber_treechains >> amberparam.exc_atoms >> amberparam.amber_bonds >> amberparam.amber_angles >> - amberparam.amber_dihedrals >> amberparam.amber_impropers >> amberparam.amber_nb14s >> - amberparam.radius_set >> amberparam.propmap >> static_cast(amberparam); + sds >> amberparam.molinfo >> amberparam.amber_charges >> + amberparam.amber_ljs >> amberparam.amber_masses >> + amberparam.amber_elements >> amberparam.amber_types >> + amberparam.born_radii >> amberparam.amber_screens >> + amberparam.amber_treechains >> amberparam.exc_atoms >> + amberparam.amber_bonds >> amberparam.amber_angles >> + amberparam.amber_dihedrals >> amberparam.amber_impropers >> + amberparam.amber_nb14s >> amberparam.radius_set >> amberparam.propmap >> + static_cast(amberparam); } else throw version_error(v, "2,3", r_amberparam, CODELOC); @@ -1315,9 +1310,8 @@ QDataStream &operator>>(QDataStream &ds, AmberParams &amberparam) } /** Null Constructor */ -AmberParams::AmberParams() : ConcreteProperty() -{ -} +AmberParams::AmberParams() + : ConcreteProperty() {} /** Constructor for the passed molecule*/ AmberParams::AmberParams(const MoleculeView &mol, const PropertyMap &map) @@ -1353,27 +1347,26 @@ AmberParams::AmberParams(const MoleculeView &mol, const PropertyMap &map) } /** Constructor for the passed molecule*/ -AmberParams::AmberParams(const MoleculeInfo &info) : ConcreteProperty(), molinfo(info) -{ -} +AmberParams::AmberParams(const MoleculeInfo &info) + : ConcreteProperty(), molinfo(info) {} /** Constructor for the passed molecule*/ AmberParams::AmberParams(const MoleculeInfoData &info) - : ConcreteProperty(), molinfo(info) -{ -} + : ConcreteProperty(), molinfo(info) {} /** Copy constructor */ AmberParams::AmberParams(const AmberParams &other) - : ConcreteProperty(), molinfo(other.molinfo), amber_charges(other.amber_charges), - amber_ljs(other.amber_ljs), amber_masses(other.amber_masses), amber_elements(other.amber_elements), - amber_types(other.amber_types), born_radii(other.born_radii), amber_screens(other.amber_screens), - amber_treechains(other.amber_treechains), exc_atoms(other.exc_atoms), amber_bonds(other.amber_bonds), - amber_angles(other.amber_angles), amber_dihedrals(other.amber_dihedrals), amber_impropers(other.amber_impropers), - amber_nb14s(other.amber_nb14s), cmap_funcs(other.cmap_funcs), - radius_set(other.radius_set), propmap(other.propmap) -{ -} + : ConcreteProperty(), molinfo(other.molinfo), + amber_charges(other.amber_charges), amber_ljs(other.amber_ljs), + amber_masses(other.amber_masses), amber_elements(other.amber_elements), + amber_types(other.amber_types), born_radii(other.born_radii), + amber_screens(other.amber_screens), + amber_treechains(other.amber_treechains), exc_atoms(other.exc_atoms), + amber_bonds(other.amber_bonds), amber_angles(other.amber_angles), + amber_dihedrals(other.amber_dihedrals), + amber_impropers(other.amber_impropers), amber_nb14s(other.amber_nb14s), + cmap_funcs(other.cmap_funcs), radius_set(other.radius_set), + propmap(other.propmap), is_perturbable(other.is_perturbable) {} /** Copy assignment operator */ AmberParams &AmberParams::operator=(const AmberParams &other) @@ -1399,27 +1392,31 @@ AmberParams &AmberParams::operator=(const AmberParams &other) cmap_funcs = other.cmap_funcs; radius_set = other.radius_set; propmap = other.propmap; + is_perturbable = other.is_perturbable; } return *this; } /** Destructor */ -AmberParams::~AmberParams() -{ -} +AmberParams::~AmberParams() {} /** Comparison operator */ bool AmberParams::operator==(const AmberParams &other) const { - return (molinfo == other.molinfo and amber_charges == other.amber_charges and amber_ljs == other.amber_ljs and - amber_masses == other.amber_masses and amber_elements == other.amber_elements and - amber_types == other.amber_types and born_radii == other.born_radii and - amber_screens == other.amber_screens and amber_treechains == other.amber_treechains and - exc_atoms == other.exc_atoms and amber_bonds == other.amber_bonds and amber_angles == other.amber_angles and - amber_dihedrals == other.amber_dihedrals and amber_impropers == other.amber_impropers and - amber_nb14s == other.amber_nb14s and cmap_funcs == other.cmap_funcs and - radius_set == other.radius_set and propmap == other.propmap); + return ( + molinfo == other.molinfo and amber_charges == other.amber_charges and + amber_ljs == other.amber_ljs and amber_masses == other.amber_masses and + amber_elements == other.amber_elements and + amber_types == other.amber_types and born_radii == other.born_radii and + amber_screens == other.amber_screens and + amber_treechains == other.amber_treechains and + exc_atoms == other.exc_atoms and amber_bonds == other.amber_bonds and + amber_angles == other.amber_angles and + amber_dihedrals == other.amber_dihedrals and + amber_impropers == other.amber_impropers and + amber_nb14s == other.amber_nb14s and cmap_funcs == other.cmap_funcs and + radius_set == other.radius_set and propmap == other.propmap); } /** Comparison operator */ @@ -1430,34 +1427,22 @@ bool AmberParams::operator!=(const AmberParams &other) const /** Return the layout of the molecule whose flexibility is contained in this object */ -MoleculeInfo AmberParams::info() const -{ - return molinfo; -} +MoleculeInfo AmberParams::info() const { return molinfo; } /** Set the property map that should be used to find and update properties of the molecule */ -void AmberParams::setPropertyMap(const PropertyMap &map) -{ - propmap = map; -} +void AmberParams::setPropertyMap(const PropertyMap &map) { propmap = map; } /** Return the property map that is used to find and update properties of the molecule */ -const PropertyMap &AmberParams::propertyMap() const -{ - return propmap; -} +const PropertyMap &AmberParams::propertyMap() const { return propmap; } /** Validate this set of parameters. This checks that all of the requirements for an Amber set of parameters are met, e.g. that all Atom indicies are contiguous and in-order, and that all atoms contiguously fill all residues etc. This returns any errors as strings. An empty set of strings indicates that there are no errors */ -QStringList AmberParams::validate() const -{ - return QStringList(); -} +QStringList AmberParams::validate() const { return QStringList(); } /** Validate this set of parameters. In addition to checking that the requirements are met, this also does any work needed to fix problems, @@ -1479,113 +1464,181 @@ QStringList AmberParams::validateAndFix() // All 1-4 scaling factors should match up with actual dihedrals - validate // that this is the case and fix any problems if we can - tbb::parallel_for(tbb::blocked_range(0, exc_atoms.nGroups()), [&](const tbb::blocked_range &r) - { - for (int icg = r.begin(); icg < r.end(); ++icg) + tbb::parallel_for( + tbb::blocked_range(0, exc_atoms.nGroups()), + [&](const tbb::blocked_range &r) { - const int nats0 = molinfo.nAtoms(CGIdx(icg)); - - for (int jcg = 0; jcg < exc_atoms.nGroups(); ++jcg) + for (int icg = r.begin(); icg < r.end(); ++icg) { - auto group_pairs = exc_atoms.get(CGIdx(icg), CGIdx(jcg)); + const int nats0 = molinfo.nAtoms(CGIdx(icg)); - if (group_pairs.isEmpty() and group_pairs.defaultValue() == CLJScaleFactor(1, 1)) + for (int jcg = 0; jcg < exc_atoms.nGroups(); ++jcg) { - // non of the pairs of atoms between these two groups are - // bonded - continue; - } - - const int nats1 = molinfo.nAtoms(CGIdx(jcg)); + auto group_pairs = exc_atoms.get(CGIdx(icg), CGIdx(jcg)); - // compare all pairs of atoms - for (int i = 0; i < nats0; ++i) - { - for (int j = 0; j < nats1; ++j) + if (group_pairs.isEmpty() and + group_pairs.defaultValue() == CLJScaleFactor(1, 1)) { - const auto s = group_pairs.get(i, j); + // non of the pairs of atoms between these two groups are + // bonded + continue; + } - if ((not(s.coulomb() == 0 or s.coulomb() == 1)) or (not(s.lj() == 0 or s.lj() == 1))) + const int nats1 = molinfo.nAtoms(CGIdx(jcg)); + + // compare all pairs of atoms + for (int i = 0; i < nats0; ++i) + { + for (int j = 0; j < nats1; ++j) { - const auto atm0 = molinfo.atomIdx(CGAtomIdx(CGIdx(icg), Index(i))); - const auto atm3 = molinfo.atomIdx(CGAtomIdx(CGIdx(jcg), Index(j))); + const auto s = group_pairs.get(i, j); - if (not has_connectivity) + // Process any non-zero 1-4 pair that isn't purely excluded + // (0,0). This includes both partial-scaling pairs (e.g. + // 0.833, 0.5 for standard AMBER) and full-interaction pairs + // (1.0, 1.0 for GLYCAM SCNB=1.0/SCEE=1.0). + if (not(s.coulomb() == 0.0 and s.lj() == 0.0)) { - // have to use the connectivity that is implied by the bonds - QMutexLocker lkr(&mutex); + const auto atm0 = + molinfo.atomIdx(CGAtomIdx(CGIdx(icg), Index(i))); + const auto atm3 = + molinfo.atomIdx(CGAtomIdx(CGIdx(jcg), Index(j))); + if (not has_connectivity) { - conn = this->connectivity(); - has_connectivity = true; + // have to use the connectivity that is implied by the + // bonds + QMutexLocker lkr(&mutex); + if (not has_connectivity) + { + conn = this->connectivity(); + has_connectivity = true; + } } - } - // find the shortest bonded paths between these two atoms - const auto paths = conn.findPaths(atm0, atm3, 4); - - for (const auto &path : paths) - { - if (path.count() != 4) + // find the shortest bonded paths between these two atoms + const auto paths = conn.findPaths(atm0, atm3, 4); + + // If the shortest bonded path between these two atoms is + // fewer than 4 atoms (i.e. they are 1-2 or 1-3), + // connectivity always enforces their exclusion from the + // non-bonded calculation regardless of what the intrascale + // says. For perturbable molecules this is expected (a + // ring-closure bond in one end-state can turn a 1-4 pair + // into a 1-3 pair in the other). For non-perturbable + // molecules it likely indicates a topology issue, so warn. { - 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; + bool has_short_path = false; + for (const auto &path : paths) + { + if (path.count() < 4) + { + has_short_path = true; + break; + } + } + if (has_short_path) + { + if (not is_perturbable) + { + QMutexLocker lkr(&mutex); + qWarning().noquote() + << QObject::tr( + "WARNING: Have a 1-4 scaling factor " + "(%1/%2) between atoms %3:%4 and %5:%6 " + "that are fewer than 4 bonds apart. " + "Skipping — connectivity will enforce " + "their exclusion. This may indicate a " + "topology issue.") + .arg(s.coulomb()) + .arg(s.lj()) + .arg(molinfo.name(atm0).value()) + .arg(atm0.value()) + .arg(molinfo.name(atm3).value()) + .arg(atm3.value()); + } + continue; + } } - // convert the atom IDs into a canonical form - auto dih = this->convert(DihedralID(path[0], path[1], path[2], path[3])); - - // skip if we already have this dihedral - if (new_dihedrals.contains(dih)) - continue; - - // qDebug() << "ADDING NULL DIHEDRAL FOR" << dih.toString(); - - // does this bond involve hydrogen? - //- this relies on "AtomElements" being full - bool contains_hydrogen = false; - - if (not amber_elements.isEmpty()) + for (const auto &path : paths) { - 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); + 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])); + + // skip if we already have this dihedral + if (new_dihedrals.contains(dih)) + continue; + + // qDebug() << "ADDING NULL DIHEDRAL FOR" << + // dih.toString(); + + // does this bond involve hydrogen? + //- this relies on "AtomElements" being full + bool contains_hydrogen = false; + + 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); + } + + // 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)); + + // now add in the 1-4 pair + BondID nb14pair = + this->convert(BondID(dih.atom0(), dih.atom3())); + + // 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)); } - - // 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)); - - // now add in the 1-4 pair - BondID nb14pair = this->convert(BondID(dih.atom0(), dih.atom3())); - - // 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)); } } } } } - } }); + }); amber_dihedrals = new_dihedrals; amber_nb14s = new_nb14s; @@ -1601,8 +1654,9 @@ QString AmberParams::toString() const if (molinfo.nAtoms() == 0) return QObject::tr("AmberParams::null"); - return QObject::tr("AmberParams( nAtoms()=%6 nBonds=%1, nAngles=%2, nDihedrals=%3 " - "nImpropers=%4 n14s=%5 )") + return QObject::tr( + "AmberParams( nAtoms()=%6 nBonds=%1, nAngles=%2, nDihedrals=%3 " + "nImpropers=%4 n14s=%5 )") .arg(amber_bonds.count()) .arg(amber_angles.count()) .arg(amber_dihedrals.count()) @@ -1674,7 +1728,8 @@ ImproperID AmberParams::convert(const ImproperID &improper) const /** Return whether or not this flexibility is compatible with the molecule whose info is in 'molinfo' */ -bool AmberParams::isCompatibleWith(const SireMol::MoleculeInfoData &molinfo) const +bool AmberParams::isCompatibleWith( + const SireMol::MoleculeInfoData &molinfo) const { return info().UID() == molinfo.UID(); } @@ -1685,52 +1740,28 @@ const char *AmberParams::typeName() } /** Return the charges on the atoms */ -AtomCharges AmberParams::charges() const -{ - return amber_charges; -} +AtomCharges AmberParams::charges() const { return amber_charges; } /** Return the atom masses */ -AtomMasses AmberParams::masses() const -{ - return amber_masses; -} +AtomMasses AmberParams::masses() const { return amber_masses; } /** Return the atom elements */ -AtomElements AmberParams::elements() const -{ - return amber_elements; -} +AtomElements AmberParams::elements() const { return amber_elements; } /** Return the atom LJ parameters */ -AtomLJs AmberParams::ljs() const -{ - return amber_ljs; -} +AtomLJs AmberParams::ljs() const { return amber_ljs; } /** Return all of the amber atom types */ -AtomStringProperty AmberParams::amberTypes() const -{ - return amber_types; -} +AtomStringProperty AmberParams::amberTypes() const { return amber_types; } /** Return all of the Born radii of the atoms */ -AtomRadii AmberParams::gbRadii() const -{ - return born_radii; -} +AtomRadii AmberParams::gbRadii() const { return born_radii; } /** Return all of the Born screening parameters for the atoms */ -AtomFloatProperty AmberParams::gbScreening() const -{ - return amber_screens; -} +AtomFloatProperty AmberParams::gbScreening() const { return amber_screens; } /** Return all of the Amber treechain classification for all of the atoms */ -AtomStringProperty AmberParams::treeChains() const -{ - return amber_treechains; -} +AtomStringProperty AmberParams::treeChains() const { return amber_treechains; } void AmberParams::createContainers() { @@ -1749,9 +1780,13 @@ void AmberParams::createContainers() } /** Set the atom parameters for the specified atom to the provided values */ -void AmberParams::add(const AtomID &atom, SireUnits::Dimension::Charge charge, SireUnits::Dimension::MolarMass mass, - const SireMol::Element &element, const SireMM::LJParameter &ljparam, const QString &amber_type, - SireUnits::Dimension::Length born_radius, double screening_parameter, const QString &treechain) +void AmberParams::add(const AtomID &atom, SireUnits::Dimension::Charge charge, + SireUnits::Dimension::MolarMass mass, + const SireMol::Element &element, + const SireMM::LJParameter &ljparam, + const QString &amber_type, + SireUnits::Dimension::Length born_radius, + double screening_parameter, const QString &treechain) { createContainers(); @@ -1770,7 +1805,8 @@ void AmberParams::add(const AtomID &atom, SireUnits::Dimension::Charge charge, S /** Set the LJ exceptions for the specified atom - this replaces any * existing exceptions */ -void AmberParams::set(const AtomID &atom, const QList &exceptions) +void AmberParams::set(const AtomID &atom, + const QList &exceptions) { createContainers(); amber_ljs.set(molinfo.atomIdx(atom).value(), exceptions); @@ -1786,12 +1822,12 @@ void AmberParams::set(const AtomID &atom0, const AtomID &atom1, other.createContainers(); amber_ljs.set(molinfo.atomIdx(atom0).value(), - other.molinfo.atomIdx(atom1).value(), - other.amber_ljs, ljparam); + other.molinfo.atomIdx(atom1).value(), other.amber_ljs, ljparam); } /** Set the LJ exception between atom0 and atom1 to 'ljparam' */ -void AmberParams::set(const AtomID &atom0, const AtomID &atom1, const LJ1264Parameter &ljparam) +void AmberParams::set(const AtomID &atom0, const AtomID &atom1, + const LJ1264Parameter &ljparam) { this->set(atom0, atom1, *this, ljparam); } @@ -1813,16 +1849,10 @@ Connectivity AmberParams::connectivity() const /** Set the radius set used by LEAP to assign the Born radii of the atoms. This is just a string that is used to label the radius set in the PRM file */ -void AmberParams::setRadiusSet(const QString &rset) -{ - radius_set = rset; -} +void AmberParams::setRadiusSet(const QString &rset) { radius_set = rset; } /** Return the radius set used by LEAP to assign the Born radii */ -QString AmberParams::radiusSet() const -{ - return radius_set; -} +QString AmberParams::radiusSet() const { return radius_set; } /** Set the excluded atoms of the molecule. This should be a CLJNBPairs with the value equal to 0 for atom0-atom1 pairs @@ -1873,10 +1903,12 @@ CLJNBPairs AmberParams::cljScaleFactors() const return nbpairs; } -void AmberParams::add(const BondID &bond, double k, double r0, bool includes_h) +void AmberParams::add(const BondID &bond, double k, double r0, + bool includes_h) { BondID b = convert(bond); - amber_bonds.insert(this->convert(bond), qMakePair(AmberBond(k, r0), includes_h)); + amber_bonds.insert(this->convert(bond), + qMakePair(AmberBond(k, r0), includes_h)); } void AmberParams::remove(const BondID &bond) @@ -1908,9 +1940,11 @@ TwoAtomFunctions AmberParams::bondFunctions() const return bondFunctions(Symbol("r")); } -void AmberParams::add(const AngleID &angle, double k, double theta0, bool includes_h) +void AmberParams::add(const AngleID &angle, double k, double theta0, + bool includes_h) { - amber_angles.insert(this->convert(angle), qMakePair(AmberAngle(k, theta0), includes_h)); + amber_angles.insert(this->convert(angle), + qMakePair(AmberAngle(k, theta0), includes_h)); } void AmberParams::remove(const AngleID &angle) @@ -1923,12 +1957,14 @@ AmberAngle AmberParams::getParameter(const AngleID &angle) const return amber_angles.value(this->convert(angle)).first; } -/** Return all of the angle parameters converted to a set of ThreeAtomFunctions */ +/** Return all of the angle parameters converted to a set of ThreeAtomFunctions + */ ThreeAtomFunctions AmberParams::angleFunctions(const Symbol &THETA) const { ThreeAtomFunctions funcs(molinfo); - for (auto it = amber_angles.constBegin(); it != amber_angles.constEnd(); ++it) + for (auto it = amber_angles.constBegin(); it != amber_angles.constEnd(); + ++it) { funcs.set(it.key(), it.value().first.toExpression(THETA)); } @@ -1936,13 +1972,15 @@ ThreeAtomFunctions AmberParams::angleFunctions(const Symbol &THETA) const return funcs; } -/** Return all of the angle parameters converted to a set of ThreeAtomFunctions */ +/** Return all of the angle parameters converted to a set of ThreeAtomFunctions + */ ThreeAtomFunctions AmberParams::angleFunctions() const { return angleFunctions(Symbol("theta")); } -void AmberParams::add(const DihedralID &dihedral, double k, double periodicity, double phase, bool includes_h) +void AmberParams::add(const DihedralID &dihedral, double k, double periodicity, + double phase, bool includes_h) { // convert the dihedral into AtomIdx indicies DihedralID d = this->convert(dihedral); @@ -1954,7 +1992,9 @@ void AmberParams::add(const DihedralID &dihedral, double k, double periodicity, } else { - amber_dihedrals.insert(d, qMakePair(AmberDihedral(AmberDihPart(k, periodicity, phase)), includes_h)); + amber_dihedrals.insert( + d, qMakePair(AmberDihedral(AmberDihPart(k, periodicity, phase)), + includes_h)); } } @@ -1968,12 +2008,14 @@ AmberDihedral AmberParams::getParameter(const DihedralID &dihedral) const return amber_dihedrals.value(this->convert(dihedral)).first; } -/** Return all of the dihedral parameters converted to a set of FourAtomFunctions */ +/** Return all of the dihedral parameters converted to a set of + * FourAtomFunctions */ FourAtomFunctions AmberParams::dihedralFunctions(const Symbol &PHI) const { FourAtomFunctions funcs(molinfo); - for (auto it = amber_dihedrals.constBegin(); it != amber_dihedrals.constEnd(); ++it) + for (auto it = amber_dihedrals.constBegin(); it != amber_dihedrals.constEnd(); + ++it) { funcs.set(it.key(), it.value().first.toExpression(PHI)); } @@ -1981,13 +2023,15 @@ FourAtomFunctions AmberParams::dihedralFunctions(const Symbol &PHI) const return funcs; } -/** Return all of the dihedral parameters converted to a set of FourAtomFunctions */ +/** Return all of the dihedral parameters converted to a set of + * FourAtomFunctions */ FourAtomFunctions AmberParams::dihedralFunctions() const { return dihedralFunctions(Symbol("phi")); } -void AmberParams::add(const ImproperID &improper, double k, double periodicity, double phase, bool includes_h) +void AmberParams::add(const ImproperID &improper, double k, double periodicity, + double phase, bool includes_h) { ImproperID imp = this->convert(improper); @@ -1997,7 +2041,9 @@ void AmberParams::add(const ImproperID &improper, double k, double periodicity, } else { - amber_impropers.insert(imp, qMakePair(AmberDihedral(AmberDihPart(k, periodicity, phase)), includes_h)); + amber_impropers.insert( + imp, qMakePair(AmberDihedral(AmberDihPart(k, periodicity, phase)), + includes_h)); } } @@ -2011,12 +2057,14 @@ AmberDihedral AmberParams::getParameter(const ImproperID &improper) const return amber_impropers.value(this->convert(improper)).first; } -/** Return all of the improper parameters converted to a set of FourAtomFunctions */ +/** Return all of the improper parameters converted to a set of + * FourAtomFunctions */ FourAtomFunctions AmberParams::improperFunctions(const Symbol &PHI) const { FourAtomFunctions funcs(molinfo); - for (auto it = amber_impropers.constBegin(); it != amber_impropers.constEnd(); ++it) + for (auto it = amber_impropers.constBegin(); it != amber_impropers.constEnd(); + ++it) { funcs.set(it.key(), it.value().first.toExpression(PHI)); } @@ -2024,7 +2072,8 @@ FourAtomFunctions AmberParams::improperFunctions(const Symbol &PHI) const return funcs; } -/** Return all of the improper parameters converted to a set of FourAtomFunctions */ +/** Return all of the improper parameters converted to a set of + * FourAtomFunctions */ FourAtomFunctions AmberParams::improperFunctions() const { return improperFunctions(Symbol("phi")); @@ -2033,17 +2082,14 @@ FourAtomFunctions AmberParams::improperFunctions() const /** Return all of the CMAP functions for the molecule. This will be empty * if there are no CMAP functions for this molecule */ -CMAPFunctions AmberParams::cmapFunctions() const -{ - return this->cmap_funcs; -} +CMAPFunctions AmberParams::cmapFunctions() const { return this->cmap_funcs; } /** Add the passed CMAP parameter for this set of 5 atoms. This will replace * any existing CMAP parameter for this set of atoms */ -void AmberParams::add(const AtomID &atom0, const AtomID &atom1, const AtomID &atom2, - const AtomID &atom3, const AtomID &atom4, - const CMAPParameter &cmap) +void AmberParams::add(const AtomID &atom0, const AtomID &atom1, + const AtomID &atom2, const AtomID &atom3, + const AtomID &atom4, const CMAPParameter &cmap) { if (cmap_funcs.isEmpty()) { @@ -2054,8 +2100,9 @@ void AmberParams::add(const AtomID &atom0, const AtomID &atom1, const AtomID &at } /** Remove the CMAP function from the passed set of 5 atoms */ -void AmberParams::removeCMAP(const AtomID &atom0, const AtomID &atom1, const AtomID &atom2, - const AtomID &atom3, const AtomID &atom4) +void AmberParams::removeCMAP(const AtomID &atom0, const AtomID &atom1, + const AtomID &atom2, const AtomID &atom3, + const AtomID &atom4) { if (cmap_funcs.isEmpty()) { @@ -2073,8 +2120,9 @@ void AmberParams::removeCMAP(const AtomID &atom0, const AtomID &atom1, const Ato /** Return the CMAP parameter for the passed 5 atoms. This returns a null * parameter if there is no matching CMAP parameter for this set of atoms */ -CMAPParameter AmberParams::getCMAP(const AtomID &atom0, const AtomID &atom1, const AtomID &atom2, - const AtomID &atom3, const AtomID &atom4) const +CMAPParameter AmberParams::getCMAP(const AtomID &atom0, const AtomID &atom1, + const AtomID &atom2, const AtomID &atom3, + const AtomID &atom4) const { if (cmap_funcs.isEmpty()) { @@ -2105,7 +2153,9 @@ AmberParams &AmberParams::operator+=(const AmberParams &other) if (not this->isCompatibleWith(other.info()) or propmap != other.propmap) { throw SireError::incompatible_error( - QObject::tr("Cannot combine Amber parameters, as the two sets are incompatible!"), CODELOC); + QObject::tr("Cannot combine Amber parameters, as the two sets are " + "incompatible!"), + CODELOC); } if (not other.amber_charges.isEmpty()) @@ -2168,7 +2218,8 @@ AmberParams &AmberParams::operator+=(const AmberParams &other) } else if (not other.amber_bonds.isEmpty()) { - for (auto it = other.amber_bonds.constBegin(); it != other.amber_bonds.constEnd(); ++it) + for (auto it = other.amber_bonds.constBegin(); + it != other.amber_bonds.constEnd(); ++it) { amber_bonds.insert(it.key(), it.value()); } @@ -2180,7 +2231,8 @@ AmberParams &AmberParams::operator+=(const AmberParams &other) } else if (not other.amber_angles.isEmpty()) { - for (auto it = other.amber_angles.constBegin(); it != other.amber_angles.constEnd(); ++it) + for (auto it = other.amber_angles.constBegin(); + it != other.amber_angles.constEnd(); ++it) { amber_angles.insert(it.key(), it.value()); } @@ -2192,7 +2244,8 @@ AmberParams &AmberParams::operator+=(const AmberParams &other) } else if (not other.amber_dihedrals.isEmpty()) { - for (auto it = other.amber_dihedrals.constBegin(); it != other.amber_dihedrals.constEnd(); ++it) + for (auto it = other.amber_dihedrals.constBegin(); + it != other.amber_dihedrals.constEnd(); ++it) { amber_dihedrals.insert(it.key(), it.value()); } @@ -2204,7 +2257,8 @@ AmberParams &AmberParams::operator+=(const AmberParams &other) } else if (not other.amber_impropers.isEmpty()) { - for (auto it = other.amber_impropers.constBegin(); it != other.amber_impropers.constEnd(); ++it) + for (auto it = other.amber_impropers.constBegin(); + it != other.amber_impropers.constEnd(); ++it) { amber_impropers.insert(it.key(), it.value()); } @@ -2228,7 +2282,8 @@ AmberParams &AmberParams::operator+=(const AmberParams &other) } else if (not other.amber_nb14s.isEmpty()) { - for (auto it = other.amber_nb14s.constBegin(); it != other.amber_nb14s.constEnd(); ++it) + for (auto it = other.amber_nb14s.constBegin(); + it != other.amber_nb14s.constEnd(); ++it) { amber_nb14s.insert(it.key(), it.value()); } @@ -2261,10 +2316,11 @@ void AmberParams::updateFrom(const MoleculeView &molview) this->_pvt_updateFrom(molview.data()); } -/** Internal function used to grab the property, catching errors and signalling if - the correct property has been found */ +/** Internal function used to grab the property, catching errors and signalling + if the correct property has been found */ template -T getProperty(const PropertyName &prop, const MoleculeData &moldata, bool *found) +T getProperty(const PropertyName &prop, const MoleculeData &moldata, + bool *found) { if (moldata.hasProperty(prop)) { @@ -2281,8 +2337,10 @@ T getProperty(const PropertyName &prop, const MoleculeData &moldata, bool *found return T(); } -/** Internal function used to guess the masses of atoms based on their element */ -void guessMasses(AtomMasses &masses, const AtomElements &elements, bool *has_masses) +/** Internal function used to guess the masses of atoms based on their element + */ +void guessMasses(AtomMasses &masses, const AtomElements &elements, + bool *has_masses) { for (int i = 0; i < elements.nCutGroups(); ++i) { @@ -2300,7 +2358,8 @@ void guessMasses(AtomMasses &masses, const AtomElements &elements, bool *has_mas } /** Internal function used to guess the element of atoms based on their name */ -AtomElements guessElements(const MoleculeInfoData &molinfo, bool *has_elements) +AtomElements guessElements(const MoleculeInfoData &molinfo, + bool *has_elements) { AtomElements elements(molinfo); @@ -2321,8 +2380,10 @@ AtomElements guessElements(const MoleculeInfoData &molinfo, bool *has_elements) } /** Construct the hash of bonds */ -void AmberParams::getAmberBondsFrom(const TwoAtomFunctions &funcs, const MoleculeData &moldata, - const QVector &is_hydrogen, const PropertyMap &map) +void AmberParams::getAmberBondsFrom(const TwoAtomFunctions &funcs, + const MoleculeData &moldata, + const QVector &is_hydrogen, + const PropertyMap &map) { // get the set of all bond functions const auto potentials = funcs.potentials(); @@ -2341,21 +2402,29 @@ void AmberParams::getAmberBondsFrom(const TwoAtomFunctions &funcs, const Molecul throw SireError::program_bug(QObject::tr("is_hydrogen is wrong!"), CODELOC); // convert each of these into an AmberBond - tbb::parallel_for(tbb::blocked_range(0, potentials.count()), [&](const tbb::blocked_range &r) - { - for (int i = r.begin(); i < r.end(); ++i) + tbb::parallel_for( + tbb::blocked_range(0, potentials.count()), + [&](const tbb::blocked_range &r) { - const auto &potential = potentials_data[i]; + for (int i = r.begin(); i < r.end(); ++i) + { + const auto &potential = potentials_data[i]; - // convert the atom IDs into a canonical form - BondID bond = this->convert(BondID(potential.atom0(), potential.atom1())); + // convert the atom IDs into a canonical form + BondID bond = + this->convert(BondID(potential.atom0(), potential.atom1())); - // does this bond involve hydrogen? - this relies on "AtomElements" being full - bool contains_hydrogen = is_hydrogen_data[molinfo.atomIdx(potential.atom0()).value()] or - is_hydrogen_data[molinfo.atomIdx(potential.atom1()).value()]; + // does this bond involve hydrogen? - this relies on "AtomElements" + // being full + bool contains_hydrogen = + is_hydrogen_data[molinfo.atomIdx(potential.atom0()).value()] or + is_hydrogen_data[molinfo.atomIdx(potential.atom1()).value()]; - bonds_data[i] = std::make_tuple(bond, AmberBond(potential.function(), Symbol("r")), contains_hydrogen); - } }); + bonds_data[i] = std::make_tuple( + bond, AmberBond(potential.function(), Symbol("r")), + contains_hydrogen); + } + }); // finally add all of these into the amber_bonds hash amber_bonds.clear(); @@ -2378,7 +2447,8 @@ void AmberParams::getAmberBondsFrom(const TwoAtomFunctions &funcs, const Molecul if (amberbond.k() != 0) { - amber_bonds.insert(std::get<0>(bonds_data[i]), qMakePair(amberbond, std::get<2>(bonds_data[i]))); + amber_bonds.insert(std::get<0>(bonds_data[i]), + qMakePair(amberbond, std::get<2>(bonds_data[i]))); } else if (keep_null_bonds) { @@ -2386,14 +2456,17 @@ void AmberParams::getAmberBondsFrom(const TwoAtomFunctions &funcs, const Molecul // to the current bond length const auto &bondid = std::get<0>(bonds_data[i]); double r0 = Bond(moldata, bondid).length(map).to(angstrom); - amber_bonds.insert(bondid, qMakePair(AmberBond(0, r0), std::get<2>(bonds_data[i]))); + amber_bonds.insert( + bondid, qMakePair(AmberBond(0, r0), std::get<2>(bonds_data[i]))); } } } /** Construct the hash of angles */ -void AmberParams::getAmberAnglesFrom(const ThreeAtomFunctions &funcs, const MoleculeData &moldata, - const QVector &is_hydrogen, const PropertyMap &map) +void AmberParams::getAmberAnglesFrom(const ThreeAtomFunctions &funcs, + const MoleculeData &moldata, + const QVector &is_hydrogen, + const PropertyMap &map) { // get the set of all angle functions const auto potentials = funcs.potentials(); @@ -2409,23 +2482,30 @@ void AmberParams::getAmberAnglesFrom(const ThreeAtomFunctions &funcs, const Mole throw SireError::program_bug(QObject::tr("is_hydrogen is wrong!"), CODELOC); // convert each of these into an AmberAngle - tbb::parallel_for(tbb::blocked_range(0, potentials.count()), [&](const tbb::blocked_range &r) - { - for (int i = r.begin(); i < r.end(); ++i) + tbb::parallel_for( + tbb::blocked_range(0, potentials.count()), + [&](const tbb::blocked_range &r) { - const auto potential = potentials.constData()[i]; - - // convert the atom IDs into a canonical form - AngleID angle = this->convert(AngleID(potential.atom0(), potential.atom1(), potential.atom2())); - - // does this angle involve hydrogen? - this relies on "AtomElements" being full - bool contains_hydrogen = is_hydrogen_data[molinfo.atomIdx(potential.atom0()).value()] or - is_hydrogen_data[molinfo.atomIdx(potential.atom1()).value()] or - is_hydrogen_data[molinfo.atomIdx(potential.atom2()).value()]; - - angles_data[i] = - std::make_tuple(angle, AmberAngle(potential.function(), Symbol("theta")), contains_hydrogen); - } }); + for (int i = r.begin(); i < r.end(); ++i) + { + const auto potential = potentials.constData()[i]; + + // convert the atom IDs into a canonical form + AngleID angle = this->convert( + AngleID(potential.atom0(), potential.atom1(), potential.atom2())); + + // does this angle involve hydrogen? - this relies on "AtomElements" + // being full + bool contains_hydrogen = + is_hydrogen_data[molinfo.atomIdx(potential.atom0()).value()] or + is_hydrogen_data[molinfo.atomIdx(potential.atom1()).value()] or + is_hydrogen_data[molinfo.atomIdx(potential.atom2()).value()]; + + angles_data[i] = std::make_tuple( + angle, AmberAngle(potential.function(), Symbol("theta")), + contains_hydrogen); + } + }); // finally add all of these into the amber_angles hash amber_angles.clear(); @@ -2448,7 +2528,8 @@ void AmberParams::getAmberAnglesFrom(const ThreeAtomFunctions &funcs, const Mole if (amberangle.k() != 0) { - amber_angles.insert(std::get<0>(angles_data[i]), qMakePair(amberangle, std::get<2>(angles_data[i]))); + amber_angles.insert(std::get<0>(angles_data[i]), + qMakePair(amberangle, std::get<2>(angles_data[i]))); } else if (keep_null_angles) { @@ -2456,20 +2537,24 @@ void AmberParams::getAmberAnglesFrom(const ThreeAtomFunctions &funcs, const Mole // to the current angle size const auto &angid = std::get<0>(angles_data[i]); double theta0 = Angle(moldata, angid).size(map).to(radians); - amber_angles.insert(angid, qMakePair(AmberAngle(0, theta0), std::get<2>(angles_data[i]))); + amber_angles.insert( + angid, qMakePair(AmberAngle(0, theta0), std::get<2>(angles_data[i]))); } } } /** Construct the hash of dihedrals */ -void AmberParams::getAmberDihedralsFrom(const FourAtomFunctions &funcs, const MoleculeData &moldata, - const QVector &is_hydrogen, const PropertyMap &map) +void AmberParams::getAmberDihedralsFrom(const FourAtomFunctions &funcs, + const MoleculeData &moldata, + const QVector &is_hydrogen, + const PropertyMap &map) { // get the set of all dihedral functions const auto potentials = funcs.potentials(); // create temporary space to hold the converted dihedrals - QVector> dihedrals(potentials.count()); + QVector> dihedrals( + potentials.count()); auto dihedrals_data = dihedrals.data(); const auto &molinfo = moldata.info(); @@ -2479,25 +2564,32 @@ void AmberParams::getAmberDihedralsFrom(const FourAtomFunctions &funcs, const Mo throw SireError::program_bug(QObject::tr("is_hydrogen is wrong!"), CODELOC); // convert each of these into an AmberDihedral - tbb::parallel_for(tbb::blocked_range(0, potentials.count()), [&](const tbb::blocked_range &r) - { - for (int i = r.begin(); i < r.end(); ++i) + tbb::parallel_for( + tbb::blocked_range(0, potentials.count()), + [&](const tbb::blocked_range &r) { - const auto potential = potentials.constData()[i]; - - // convert the atom IDs into a canonical form - DihedralID dihedral = - this->convert(DihedralID(potential.atom0(), potential.atom1(), potential.atom2(), potential.atom3())); - - // does this bond involve hydrogen? - this relies on "AtomElements" being full - bool contains_hydrogen = is_hydrogen_data[molinfo.atomIdx(potential.atom0()).value()] or - is_hydrogen_data[molinfo.atomIdx(potential.atom1()).value()] or - is_hydrogen_data[molinfo.atomIdx(potential.atom2()).value()] or - is_hydrogen_data[molinfo.atomIdx(potential.atom3()).value()]; - - dihedrals_data[i] = - std::make_tuple(dihedral, AmberDihedral(potential.function(), Symbol("phi")), contains_hydrogen); - } }); + for (int i = r.begin(); i < r.end(); ++i) + { + const auto potential = potentials.constData()[i]; + + // convert the atom IDs into a canonical form + DihedralID dihedral = + this->convert(DihedralID(potential.atom0(), potential.atom1(), + potential.atom2(), potential.atom3())); + + // does this bond involve hydrogen? - this relies on "AtomElements" + // being full + bool contains_hydrogen = + is_hydrogen_data[molinfo.atomIdx(potential.atom0()).value()] or + is_hydrogen_data[molinfo.atomIdx(potential.atom1()).value()] or + is_hydrogen_data[molinfo.atomIdx(potential.atom2()).value()] or + is_hydrogen_data[molinfo.atomIdx(potential.atom3()).value()]; + + dihedrals_data[i] = std::make_tuple( + dihedral, AmberDihedral(potential.function(), Symbol("phi")), + contains_hydrogen); + } + }); // finally add all of these into the amber_dihedrals hash amber_dihedrals.clear(); @@ -2506,19 +2598,23 @@ void AmberParams::getAmberDihedralsFrom(const FourAtomFunctions &funcs, const Mo for (int i = 0; i < dihedrals.count(); ++i) { amber_dihedrals.insert(std::get<0>(dihedrals_data[i]), - qMakePair(std::get<1>(dihedrals_data[i]), std::get<2>(dihedrals_data[i]))); + qMakePair(std::get<1>(dihedrals_data[i]), + std::get<2>(dihedrals_data[i]))); } } /** Construct the hash of impropers */ -void AmberParams::getAmberImpropersFrom(const FourAtomFunctions &funcs, const MoleculeData &moldata, - const QVector &is_hydrogen, const PropertyMap &map) +void AmberParams::getAmberImpropersFrom(const FourAtomFunctions &funcs, + const MoleculeData &moldata, + const QVector &is_hydrogen, + const PropertyMap &map) { // get the set of all improper functions const auto potentials = funcs.potentials(); // create temporary space to hold the converted dihedrals - QVector> impropers(potentials.count()); + QVector> impropers( + potentials.count()); auto impropers_data = impropers.data(); const auto &molinfo = moldata.info(); @@ -2528,25 +2624,32 @@ void AmberParams::getAmberImpropersFrom(const FourAtomFunctions &funcs, const Mo throw SireError::program_bug(QObject::tr("is_hydrogen is wrong!"), CODELOC); // convert each of these into an AmberDihedral - tbb::parallel_for(tbb::blocked_range(0, potentials.count()), [&](const tbb::blocked_range &r) - { - for (int i = r.begin(); i < r.end(); ++i) + tbb::parallel_for( + tbb::blocked_range(0, potentials.count()), + [&](const tbb::blocked_range &r) { - const auto potential = potentials.constData()[i]; - - // convert the atom IDs into a canonical form - ImproperID improper = - this->convert(ImproperID(potential.atom0(), potential.atom1(), potential.atom2(), potential.atom3())); - - // does this bond involve hydrogen? - this relies on "AtomElements" being full - bool contains_hydrogen = is_hydrogen_data[molinfo.atomIdx(potential.atom0()).value()] or - is_hydrogen_data[molinfo.atomIdx(potential.atom1()).value()] or - is_hydrogen_data[molinfo.atomIdx(potential.atom2()).value()] or - is_hydrogen_data[molinfo.atomIdx(potential.atom3()).value()]; - - impropers_data[i] = - std::make_tuple(improper, AmberDihedral(potential.function(), Symbol("phi")), contains_hydrogen); - } }); + for (int i = r.begin(); i < r.end(); ++i) + { + const auto potential = potentials.constData()[i]; + + // convert the atom IDs into a canonical form + ImproperID improper = + this->convert(ImproperID(potential.atom0(), potential.atom1(), + potential.atom2(), potential.atom3())); + + // does this bond involve hydrogen? - this relies on "AtomElements" + // being full + bool contains_hydrogen = + is_hydrogen_data[molinfo.atomIdx(potential.atom0()).value()] or + is_hydrogen_data[molinfo.atomIdx(potential.atom1()).value()] or + is_hydrogen_data[molinfo.atomIdx(potential.atom2()).value()] or + is_hydrogen_data[molinfo.atomIdx(potential.atom3()).value()]; + + impropers_data[i] = std::make_tuple( + improper, AmberDihedral(potential.function(), Symbol("phi")), + contains_hydrogen); + } + }); // finally add all of these into the amber_dihedrals hash amber_impropers.clear(); @@ -2555,12 +2658,14 @@ void AmberParams::getAmberImpropersFrom(const FourAtomFunctions &funcs, const Mo for (int i = 0; i < impropers.count(); ++i) { amber_impropers.insert(std::get<0>(impropers_data[i]), - qMakePair(std::get<1>(impropers_data[i]), std::get<2>(impropers_data[i]))); + qMakePair(std::get<1>(impropers_data[i]), + std::get<2>(impropers_data[i]))); } } /** Construct the excluded atom set and 14 NB parameters */ -void AmberParams::getAmberNBsFrom(const CLJNBPairs &nbpairs, const FourAtomFunctions &dihedrals) +void AmberParams::getAmberNBsFrom(const CLJNBPairs &nbpairs, + const FourAtomFunctions &dihedrals) { // first, copy in the CLJNBPairs from the molecule exc_atoms = nbpairs; @@ -2578,7 +2683,8 @@ void AmberParams::getAmberNBsFrom(const CLJNBPairs &nbpairs, const FourAtomFunct const auto potential = potentials.constData()[i]; // convert the atom IDs into a canonical form - BondID nb14pair = this->convert(BondID(potential.atom0(), potential.atom3())); + BondID nb14pair = + this->convert(BondID(potential.atom0(), potential.atom3())); const auto molinfo = info(); @@ -2587,24 +2693,15 @@ void AmberParams::getAmberNBsFrom(const CLJNBPairs &nbpairs, const FourAtomFunct // extract the nb14 term from exc_atoms auto nbscl = nbpairs.get(nb14pair.atom0(), nb14pair.atom1()); - if (nbscl.coulomb() != 1.0 or nbscl.lj() != 1.0) + if (nbscl.coulomb() != 0.0 or nbscl.lj() != 0.0) { - if (nbscl.coulomb() != 0.0 or nbscl.lj() != 0.0) - { - // add them to the list of 14 scale factors - new_nb14s.insert(nb14pair, AmberNB14(nbscl.coulomb(), nbscl.lj())); + // add them to the list of 14 scale factors. + // This handles both standard AMBER (e.g. 0.833, 0.5) and + // GLYCAM-style (1.0, 1.0) where SCEE=1.0 and SCNB=1.0. + new_nb14s.insert(nb14pair, AmberNB14(nbscl.coulomb(), nbscl.lj())); - // and remove them from the excluded atoms list - exc_atoms.set(nb14pair.atom0(), nb14pair.atom1(), CLJScaleFactor(0)); - } - else - { - const auto tscl = nbpairs.get(nb14pair.atom0(), nb14pair.atom1()); - } - } - else - { - const auto tscl = nbpairs.get(nb14pair.atom0(), nb14pair.atom1()); + // and remove them from the excluded atoms list + exc_atoms.set(nb14pair.atom0(), nb14pair.atom1(), CLJScaleFactor(0)); } } } @@ -2618,24 +2715,33 @@ void AmberParams::_pvt_createFrom(const MoleculeData &moldata) // pull out all of the molecular properties const PropertyMap &map = propmap; + // check if this is a perturbable molecule + is_perturbable = moldata.hasProperty("is_perturbable") and + moldata.property("is_perturbable").asABoolean(); + // first, all of the atom-based properties bool has_charges, has_ljs, has_masses, has_elements, has_ambertypes; - amber_charges = getProperty(map["charge"], moldata, &has_charges); + amber_charges = + getProperty(map["charge"], moldata, &has_charges); amber_ljs = getProperty(map["LJ"], moldata, &has_ljs); amber_masses = getProperty(map["mass"], moldata, &has_masses); - amber_elements = getProperty(map["element"], moldata, &has_elements); - amber_types = getProperty(map["ambertype"], moldata, &has_ambertypes); + amber_elements = + getProperty(map["element"], moldata, &has_elements); + amber_types = getProperty(map["ambertype"], moldata, + &has_ambertypes); if (not has_ambertypes) { // first look for the bondtypes property, as OPLS uses this - amber_types = getProperty(map["bondtype"], moldata, &has_ambertypes); + amber_types = getProperty(map["bondtype"], moldata, + &has_ambertypes); if (not has_ambertypes) { // look for the atomtypes property - amber_types = getProperty(map["atomtype"], moldata, &has_ambertypes); + amber_types = getProperty(map["atomtype"], moldata, + &has_ambertypes); } } @@ -2655,14 +2761,17 @@ void AmberParams::_pvt_createFrom(const MoleculeData &moldata) } } - if (not(has_charges and has_ljs and has_masses and has_elements and has_ambertypes)) + if (not(has_charges and has_ljs and has_masses and has_elements and + has_ambertypes)) { // it is not possible to create the parameter object if we don't have // these atom-based parameters throw SireBase::missing_property( - QObject::tr("Cannot create an AmberParams object for molecule %1 as it is missing " - "necessary atom based properties: has_charges = %2, has_ljs = %3, " - "has_masses = %4, has_elements = %5, has_ambertypes = %6.") + QObject::tr( + "Cannot create an AmberParams object for molecule %1 as it is " + "missing " + "necessary atom based properties: has_charges = %2, has_ljs = %3, " + "has_masses = %4, has_elements = %5, has_ambertypes = %6.") .arg(Molecule(moldata).toString()) .arg(has_charges) .arg(has_ljs) @@ -2676,15 +2785,19 @@ void AmberParams::_pvt_createFrom(const MoleculeData &moldata) bool has_radii, has_screening, has_treechains; born_radii = getProperty(map["gb_radii"], moldata, &has_radii); - amber_screens = getProperty(map["gb_screening"], moldata, &has_screening); - amber_treechains = getProperty(map["treechain"], moldata, &has_treechains); + amber_screens = getProperty(map["gb_screening"], moldata, + &has_screening); + amber_treechains = getProperty(map["treechain"], moldata, + &has_treechains); if (has_radii) { // see if there is a label for the source of the GB parameters bool has_source; - radius_set = getProperty(map["gb_radius_set"], moldata, &has_source).value(); + radius_set = + getProperty(map["gb_radius_set"], moldata, &has_source) + .value(); if (not has_source) { @@ -2697,15 +2810,23 @@ void AmberParams::_pvt_createFrom(const MoleculeData &moldata) } // now lets get the bonded parameters (if they exist...) - bool has_bonds, has_ubs, has_angles, has_dihedrals, has_impropers, has_nbpairs, has_cmaps; - - const auto bonds = getProperty(map["bond"], moldata, &has_bonds); - const auto ub_bonds = getProperty(map["urey_bradley"], moldata, &has_ubs); - const auto angles = getProperty(map["angle"], moldata, &has_angles); - const auto dihedrals = getProperty(map["dihedral"], moldata, &has_dihedrals); - const auto impropers = getProperty(map["improper"], moldata, &has_impropers); - const auto nbpairs = getProperty(map["intrascale"], moldata, &has_nbpairs); - const auto cmaps = getProperty(map["cmap"], moldata, &has_cmaps); + bool has_bonds, has_ubs, has_angles, has_dihedrals, has_impropers, + has_nbpairs, has_cmaps; + + const auto bonds = + getProperty(map["bond"], moldata, &has_bonds); + const auto ub_bonds = + getProperty(map["urey_bradley"], moldata, &has_ubs); + const auto angles = + getProperty(map["angle"], moldata, &has_angles); + const auto dihedrals = + getProperty(map["dihedral"], moldata, &has_dihedrals); + const auto impropers = + getProperty(map["improper"], moldata, &has_impropers); + const auto nbpairs = + getProperty(map["intrascale"], moldata, &has_nbpairs); + const auto cmaps = + getProperty(map["cmap"], moldata, &has_cmaps); // get all of the atoms that contain hydrogen QVector is_hydrogen; @@ -2723,7 +2844,8 @@ void AmberParams::_pvt_createFrom(const MoleculeData &moldata) auto elements = amber_elements.toVector(); if (elements.count() != natoms) - throw SireError::program_bug(QObject::tr("Wrong elements count!"), CODELOC); + throw SireError::program_bug(QObject::tr("Wrong elements count!"), + CODELOC); const auto *elements_data = elements.constData(); @@ -2747,32 +2869,37 @@ void AmberParams::_pvt_createFrom(const MoleculeData &moldata) if (has_bonds) { - nb_functions.append([&]() - { getAmberBondsFrom(bonds, moldata, is_hydrogen, map); }); + nb_functions.append( + [&]() + { getAmberBondsFrom(bonds, moldata, is_hydrogen, map); }); } if (has_ubs) { - nb_functions.append([&]() - { getAmberBondsFrom(ub_bonds, moldata, is_hydrogen, map); }); + nb_functions.append( + [&]() + { getAmberBondsFrom(ub_bonds, moldata, is_hydrogen, map); }); } if (has_angles) { - nb_functions.append([&]() - { getAmberAnglesFrom(angles, moldata, is_hydrogen, map); }); + nb_functions.append( + [&]() + { getAmberAnglesFrom(angles, moldata, is_hydrogen, map); }); } if (has_dihedrals) { - nb_functions.append([&]() - { getAmberDihedralsFrom(dihedrals, moldata, is_hydrogen, map); }); + nb_functions.append( + [&]() + { getAmberDihedralsFrom(dihedrals, moldata, is_hydrogen, map); }); } if (has_impropers) { - nb_functions.append([&]() - { getAmberImpropersFrom(impropers, moldata, is_hydrogen, map); }); + nb_functions.append( + [&]() + { getAmberImpropersFrom(impropers, moldata, is_hydrogen, map); }); } if (has_nbpairs) @@ -2788,12 +2915,14 @@ void AmberParams::_pvt_createFrom(const MoleculeData &moldata) if (not errors.isEmpty()) { - throw SireError::io_error(QObject::tr("Problem creating the AmberParams object for molecule %1 : %2. " - "Errors include:\n%3") - .arg(moldata.name().value()) - .arg(moldata.number().value()) - .arg(errors.join("\n")), - CODELOC); + throw SireError::io_error( + QObject::tr( + "Problem creating the AmberParams object for molecule %1 : %2. " + "Errors include:\n%3") + .arg(moldata.name().value()) + .arg(moldata.number().value()) + .arg(errors.join("\n")), + CODELOC); } } @@ -2814,7 +2943,9 @@ void AmberParams::_pvt_updateFrom(const MoleculeData &moldata) this->_pvt_createFrom(moldata); } -PropertyPtr AmberParams::_pvt_makeCompatibleWith(const MoleculeInfoData &newinfo, const AtomMatcher &atommatcher) const +PropertyPtr +AmberParams::_pvt_makeCompatibleWith(const MoleculeInfoData &newinfo, + const AtomMatcher &atommatcher) const { try { @@ -2826,7 +2957,8 @@ PropertyPtr AmberParams::_pvt_makeCompatibleWith(const MoleculeInfoData &newinfo return ret; } - QHash matched_atoms = atommatcher.match(this->info(), molinfo); + QHash matched_atoms = + atommatcher.match(this->info(), molinfo); return this->_pvt_makeCompatibleWith(molinfo, matched_atoms); } @@ -2837,8 +2969,9 @@ PropertyPtr AmberParams::_pvt_makeCompatibleWith(const MoleculeInfoData &newinfo } } -PropertyPtr AmberParams::_pvt_makeCompatibleWith(const MoleculeInfoData &newinfo, - const QHash &map) const +PropertyPtr +AmberParams::_pvt_makeCompatibleWith(const MoleculeInfoData &newinfo, + const QHash &map) const { if (map.isEmpty()) { @@ -2848,7 +2981,8 @@ PropertyPtr AmberParams::_pvt_makeCompatibleWith(const MoleculeInfoData &newinfo return ret; } - throw SireError::incomplete_code("Cannot make compatible if atom order has changed!", CODELOC); + throw SireError::incomplete_code( + "Cannot make compatible if atom order has changed!", CODELOC); } /** Merge this property with another property */ @@ -2859,14 +2993,17 @@ PropertyList AmberParams::merge(const MolViewProperty &other, { if (not other.isA()) { - throw SireError::incompatible_error(QObject::tr("Cannot merge %1 with %2 as they are different types.") - .arg(this->what()) - .arg(other.what()), - CODELOC); + throw SireError::incompatible_error( + QObject::tr("Cannot merge %1 with %2 as they are different types.") + .arg(this->what()) + .arg(other.what()), + CODELOC); } - SireBase::Console::warning(QObject::tr("Merging %1 properties is not yet implemented. Returning two copies of the original property.") - .arg(this->what())); + SireBase::Console::warning( + QObject::tr("Merging %1 properties is not yet implemented. Returning two " + "copies of the original property.") + .arg(this->what())); SireBase::PropertyList ret; diff --git a/corelib/src/libs/SireMM/amberparams.h b/corelib/src/libs/SireMM/amberparams.h index 0e6ffc6db..e20d4b46f 100644 --- a/corelib/src/libs/SireMM/amberparams.h +++ b/corelib/src/libs/SireMM/amberparams.h @@ -605,6 +605,11 @@ namespace SireMM /** The property map used (if any) to identify the properties that hold the amber parameters */ SireBase::PropertyMap propmap; + + /** Whether this molecule is perturbable (has an is_perturbable property). + Used in validateAndFix to handle ring-closure 1-3 pairs that may + carry incorrect (1,1) CLJScaleFactors from old-style merges. */ + bool is_perturbable = false; }; #ifndef SIRE_SKIP_INLINE_FUNCTIONS diff --git a/corelib/src/libs/SireMM/atompairs.hpp b/corelib/src/libs/SireMM/atompairs.hpp index 4cf92c2ab..82d2ee880 100644 --- a/corelib/src/libs/SireMM/atompairs.hpp +++ b/corelib/src/libs/SireMM/atompairs.hpp @@ -125,6 +125,8 @@ namespace SireMM const T &defaultValue() const; + QList> nonDefaultElements() const; + private: /** The matrix of objects associated with each pair of atoms */ SireBase::SparseMatrix data; @@ -352,6 +354,13 @@ namespace SireMM return data.defaultValue(); } + /** Return all non-default (i, j, value) elements of this atom pairs matrix */ + template + SIRE_INLINE_TEMPLATE QList> CGAtomPairs::nonDefaultElements() const + { + return data.nonDefaultElements(); + } + /** Make sure that this container can hold dim_x by dim_y atom pairs */ template SIRE_OUTOFLINE_TEMPLATE void CGAtomPairs::reserve(int dim_x, int dim_y) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index e3842d002..bfe5a1bff 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -23,6 +23,8 @@ organisation on `GitHub `__. * Add convenience function to ``sire.mol.dynamics`` to get current energy trajectory records. +* Fixed parsing of AMBER and GROMACS GLYCAM force field topologies. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/tests/io/test_grotop.py b/tests/io/test_grotop.py index d1399d475..4dd5bc91e 100644 --- a/tests/io/test_grotop.py +++ b/tests/io/test_grotop.py @@ -1,3 +1,5 @@ +import sys + import sire as sr import pytest @@ -375,3 +377,108 @@ def test_grotop_water(tmpdir, water_model): is_settles = True if line.startswith("[ vsite3 ]"): is_vsite = True + + +@pytest.mark.slow +def test_glycam(tmpdir): + """Test that a topology using the GLYCAM force field (SCEE=1.0, SCNB=1.0) + is read and written correctly. + + GLYCAM uses full 1-4 interactions (no scaling), unlike standard AMBER + which uses SCEE=1.2, SCNB=2.0. In GROMACS this is represented by funct=2 + pairs with explicit LJ parameters rather than funct=1 which would apply + the global fudgeLJ/fudgeQQ scaling. + """ + + # Load the GLYCAM topology and coordinates. + mols = sr.load_test_files("glycam.top", "glycam.gro") + + # The system contains a protein (PROA) and a glycan (CARB). + # Find them by molecule name. + glycan = mols["molname CARB"] + protein = mols["molname PROA"] + + # The glycan should use GLYCAM-style full 1-4 interactions. + # Its LJ property should be readable and non-trivial. + glycan_lj = glycan.property("LJ") + assert glycan_lj is not None + + protein_lj = protein.property("LJ") + assert protein_lj is not None + + # Write the whole system to GROMACS topology + coordinate files. + d = tmpdir.mkdir("test_glycam") + f = sr.save(mols, d.join("glycam_out"), format=["GroTop", "Gro87"]) + + # Parse the written file to verify the [pairs] sections are correct. + # The glycan (CARB) molecule type must use funct=2 (explicit LJ pairs, + # fudgeQQ=1.0) because its 1-4 interactions are unscaled. + # The protein (PROA) molecule type uses funct=1 (standard AMBER scaling). + current_moltype = None + expect_moltype_name = False + in_pairs = False + carb_funct2_count = 0 + carb_funct1_count = 0 + proa_funct1_count = 0 + + with open(f[0], "r") as fh: + for line in fh: + stripped = line.strip() + if not stripped or stripped.startswith(";"): + continue + if stripped == "[ moleculetype ]": + expect_moltype_name = True + in_pairs = False + continue + if expect_moltype_name: + current_moltype = stripped.split()[0] + expect_moltype_name = False + continue + if stripped.startswith("["): + in_pairs = stripped == "[ pairs ]" + continue + if in_pairs: + words = stripped.split() + if len(words) >= 3: + funct = int(words[2]) + if current_moltype == "CARB": + if funct == 2: + carb_funct2_count += 1 + else: + carb_funct1_count += 1 + elif current_moltype == "PROA": + if funct == 1: + proa_funct1_count += 1 + + # The CARB glycan uses GLYCAM (SCEE=1.0, SCNB=1.0), so the vast majority + # of its 1-4 pairs must be funct=2. The original topology has exactly one + # funct=1 pair (atoms 8-14), which is preserved as funct=1. + assert carb_funct2_count > 0, "No funct=2 pairs found for CARB glycan" + assert ( + carb_funct2_count > carb_funct1_count + ), "CARB glycan has more funct=1 pairs than funct=2; GLYCAM scaling not applied" + + # Protein pairs must use funct=1 (standard AMBER scaling). + assert proa_funct1_count > 0, "No funct=1 pairs found for PROA protein" + + # Reload and verify the energy is self-consistent after the GROMACS roundtrip. + # Before the fix, funct=2 pairs were written as funct=1, which would apply + # fudgeLJ=0.5 to the glycan 1-4 LJ interactions and give a different energy. + mols2 = sr.load(f, show_warnings=False) + glycan2 = mols2["molname CARB"] + + glycan_lj2 = glycan2.property("LJ") + + # Check a representative atom (index 1, type Oh) whose sigma and epsilon + # survive the roundtrip within floating-point formatting precision. + assert glycan_lj2[1].sigma().value() == pytest.approx( + glycan_lj[1].sigma().value(), rel=1e-3 + ) + assert glycan_lj2[1].epsilon().value() == pytest.approx( + glycan_lj[1].epsilon().value(), rel=1e-3 + ) + + # The CLJ energy calculation returns NaN for large solvated systems on Windows + # (a separate Windows-specific bug). Skip the energy check on that platform. + if sys.platform != "win32": + assert mols2.energy().value() == pytest.approx(mols.energy().value(), rel=1e-3) diff --git a/tests/io/test_prmtop.py b/tests/io/test_prmtop.py index 3c790ca3d..07e27f16d 100644 --- a/tests/io/test_prmtop.py +++ b/tests/io/test_prmtop.py @@ -1,3 +1,5 @@ +import sys + import sire as sr import pytest @@ -271,3 +273,80 @@ def test_reorder_prmtop(reordered_protein): assert mols[0]["element Zn"].num_atoms() == 2 assert mols[1]["element Zn"].num_atoms() == 2 + + +@pytest.mark.slow +def test_glycam(tmpdir): + """Test that a topology using the GLYCAM force field (SCEE=1.0, SCNB=1.0) + round-trips correctly through AMBER prm7 format. + + GLYCAM uses full 1-4 interactions (no scaling), so SCEE=1.0 and SCNB=1.0 + for glycan dihedrals. The protein dihedrals use standard AMBER scaling + (SCEE=1.2, SCNB=2.0). Before the fix, CLJScaleFactor(1.0, 1.0) pairs were + silently dropped when building AmberParams, so glycan dihedrals were written + with SCEE=0 and SCNB=0, giving zero 1-4 interactions. + """ + + # Load the GLYCAM topology and coordinates. + mols = sr.load_test_files("glycam.top", "glycam.gro") + + # Write to AMBER prm7 + rst7 format. + d = tmpdir.mkdir("test_glycam_amber") + f = sr.save(mols, d.join("glycam_out"), format=["PRM7", "RST7"]) + + # Parse SCEE_SCALE_FACTOR and SCNB_SCALE_FACTOR from the written prm7. + # Format: 5 values per line, 16 chars each (AmberFormat FLOAT 5 16 8). + scee_values = [] + scnb_values = [] + reading = None + + with open(f[0], "r") as fh: + for line in fh: + if line.startswith("%FLAG SCEE_SCALE_FACTOR"): + reading = "scee" + continue + elif line.startswith("%FLAG SCNB_SCALE_FACTOR"): + reading = "scnb" + continue + elif line.startswith("%FLAG") or line.startswith("%FORMAT"): + if line.startswith("%FLAG"): + reading = None + continue + if reading == "scee": + scee_values.extend(float(x) for x in line.split()) + elif reading == "scnb": + scnb_values.extend(float(x) for x in line.split()) + + assert len(scee_values) > 0, "No SCEE_SCALE_FACTOR entries found" + assert len(scnb_values) > 0, "No SCNB_SCALE_FACTOR entries found" + + # The system has both glycan (SCEE=1.0, SCNB=1.0) and protein + # (SCEE=1.2, SCNB=2.0) dihedrals. Both values must be present. + # Before the fix, all glycan entries would be 0.0. + assert any( + v == pytest.approx(1.0) for v in scee_values + ), "No SCEE=1.0 entries found; GLYCAM glycan dihedrals were not written correctly" + assert any( + v == pytest.approx(1.2, rel=1e-3) for v in scee_values + ), "No SCEE=1.2 entries found; standard AMBER protein dihedrals are missing" + assert any( + v == pytest.approx(1.0) for v in scnb_values + ), "No SCNB=1.0 entries found; GLYCAM glycan dihedrals were not written correctly" + assert any( + v == pytest.approx(2.0, rel=1e-3) for v in scnb_values + ), "No SCNB=2.0 entries found; standard AMBER protein dihedrals are missing" + + # SCEE/SCNB=0.0 is valid and expected — it marks dihedrals that share + # terminal atoms with another dihedral and should not contribute a 1-4 term. + # The pre-fix bug was that ALL glycan dihedral entries were 0.0 because + # CLJScaleFactor(1.0, 1.0) pairs were silently dropped. Having both 1.0 and + # 1.2 present (checked above) confirms the fix is working correctly. + + # Reload and verify the energy is self-consistent after the AMBER roundtrip. + # Before the fix, glycan 1-4 pairs had SCEE=0/SCNB=0, giving zero 1-4 + # interactions and a different energy. + mols2 = sr.load(f) + # The CLJ energy calculation returns NaN for large solvated systems on Windows + # (a separate Windows-specific bug). Skip the energy check on that platform. + if sys.platform != "win32": + assert mols2.energy().value() == pytest.approx(mols.energy().value(), rel=1e-3) diff --git a/wrapper/IO/_IO_free_functions.pypp.cpp b/wrapper/IO/_IO_free_functions.pypp.cpp index 4a5b57a62..e307165f0 100644 --- a/wrapper/IO/_IO_free_functions.pypp.cpp +++ b/wrapper/IO/_IO_free_functions.pypp.cpp @@ -7,6 +7,13 @@ namespace bp = boost::python; +#include "SireBase/propertylist.h" + +#include "SireMM/cljnbpairs.h" + +#include "SireMol/atomidxmapping.h" +#include "SireMol/moleculeinfodata.h" + #include "SireBase/getinstalldir.h" #include "SireBase/parallel.h" @@ -871,16 +878,45 @@ void register_free_functions(){ } { //::SireIO::updateCoordinatesAndVelocities - + typedef ::boost::tuples::tuple< SireSystem::System, QHash< SireMol::MolIdx, SireMol::MolIdx >, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > ( *updateCoordinatesAndVelocities_function_type )( ::SireSystem::System const &,::SireSystem::System const &,::SireSystem::System const &,::QHash< SireMol::MolIdx, SireMol::MolIdx > const &,bool const,::SireBase::PropertyMap const &,::SireBase::PropertyMap const & ); updateCoordinatesAndVelocities_function_type updateCoordinatesAndVelocities_function_value( &::SireIO::updateCoordinatesAndVelocities ); - - bp::def( + + bp::def( "updateCoordinatesAndVelocities" , updateCoordinatesAndVelocities_function_value , ( bp::arg("original_system"), bp::arg("renumbered_system"), bp::arg("updated_system"), bp::arg("molecule_mapping"), bp::arg("is_lambda1")=(bool const)(false), bp::arg("map0")=SireBase::PropertyMap(), bp::arg("map1")=SireBase::PropertyMap() ) , "Update the coordinates and velocities of original_system with those from\nupdated_system.\n\nPar:am system_original\nThe original system.\n\nPar:am system_renumbered\nThe original system, atoms and residues have been renumbered to be\nunique and in ascending order.\n\nPar:am system_updated\nThe updated system, where molecules may not be in the same order.\n\nPar:am map0\nA dictionary of user-defined molecular property names for system0.\n\nPar:am map1\nA dictionary of user-defined molecular property names for system1.\n\nRetval: system, mapping\nThe system with updated coordinates and velocities and a mapping\nbetween the molecule indices in both systems.\n" ); - + + } + + { //::SireIO::mergeIntrascale + + typedef ::SireBase::PropertyList ( *mergeIntrascale_function_type )( + ::SireMM::CLJNBPairs const &, + ::SireMM::CLJNBPairs const &, + ::SireMol::MoleculeInfoData const &, + ::QHash< SireMol::AtomIdx, SireMol::AtomIdx > const &, + ::QHash< SireMol::AtomIdx, SireMol::AtomIdx > const & ); + mergeIntrascale_function_type mergeIntrascale_function_value( &::SireIO::mergeIntrascale ); + + bp::def( + "mergeIntrascale" + , mergeIntrascale_function_value + , ( bp::arg("nb0"), bp::arg("nb1"), bp::arg("merged_info"), + bp::arg("mol0_merged_mapping"), bp::arg("mol1_merged_mapping") ) + , "Merge the CLJNBPairs (intrascale) of two molecules into the end-state\n" + "intrascales of a perturbable merged molecule.\n" + "\n" + "Uses a two-pass approach to correctly preserve actual per-pair scale factors\n" + "(including GLYCAM funct=2 overrides of 1.0/1.0) without relying on global\n" + "forcefield scale factors. For intrascale0, nb1 is copied first (so ghost\n" + "atoms from mol1 get correct exclusions), then nb0 overwrites (so mapped and\n" + "mol0-unique atoms get correct state-0 values). intrascale1 is built with the\n" + "same logic in reverse.\n" + "\n" + "Returns a PropertyList [intrascale0, intrascale1].\n" ); + } }