Skip to content

Commit b281b88

Browse files
lohedgeslohedges
committed
Add fast pre-compute of genuine 1-4 pairs.
1 parent ac18e40 commit b281b88

1 file changed

Lines changed: 168 additions & 194 deletions

File tree

‎corelib/src/libs/SireIO/grotop.cpp‎

Lines changed: 168 additions & 194 deletions
Original file line numberDiff line numberDiff line change
@@ -4093,6 +4093,34 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype,
40934093
// Store the molinfo object;
40944094
const auto molinfo = mol.info();
40954095

4096+
// Precompute the set of genuine 1-4 bonded atom pairs once. This lets
4097+
// us distinguish GLYCAM-style 1-4 pairs (CLJScaleFactor(1,1)) from
4098+
// 1-5+ pairs that return the same default value but must not appear in
4099+
// [pairs]. Empty if connectivity is unavailable (fallback: write all).
4100+
QSet<QPair<AtomIdx, AtomIdx>> pairs_14;
4101+
4102+
try {
4103+
const auto conn = mol.property("connectivity").asA<Connectivity>();
4104+
const int natoms = molinfo.nAtoms();
4105+
for (int a = 0; a < natoms; ++a) {
4106+
const AtomIdx aidx(a);
4107+
for (const auto &b : conn.connectionsTo(aidx)) {
4108+
for (const auto &c : conn.connectionsTo(b)) {
4109+
if (c == aidx)
4110+
continue;
4111+
for (const auto &d : conn.connectionsTo(c)) {
4112+
if (d == b or d == aidx)
4113+
continue;
4114+
// aidx-b-c-d is a dihedral: aidx and d are 1-4 bonded.
4115+
pairs_14.insert(QPair<AtomIdx, AtomIdx>(aidx, d));
4116+
pairs_14.insert(QPair<AtomIdx, AtomIdx>(d, aidx));
4117+
}
4118+
}
4119+
}
4120+
}
4121+
} catch (...) {
4122+
}
4123+
40964124
if (is_perturbable) {
40974125
CLJNBPairs scl0;
40984126
CLJNBPairs scl1;
@@ -4130,132 +4158,95 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype,
41304158
} catch (...) {
41314159
}
41324160

4133-
// A set of recorded 1-4 pairs.
4134-
QSet<QPair<AtomIdx, AtomIdx>> recorded_pairs;
4135-
41364161
bool fix_null_perturbable_14s = false;
41374162

41384163
if (mol.hasProperty("fix_null_perturbable_14s"))
41394164
fix_null_perturbable_14s = mol.property("fix_null_perturbable_14s")
41404165
.asA<BooleanProperty>()
41414166
.value();
41424167

4143-
// Must record every pair that has a non-default scaling factor.
4144-
// Loop over intrascale matrix by cut-groups to avoid N^2 loop.
4145-
for (int i = 0; i < scl0.nGroups(); ++i) {
4146-
for (int j = 0; j < scl0.nGroups(); ++j) {
4147-
const auto s0 = scl0.get(CGIdx(i), CGIdx(j));
4148-
const auto s1 = scl1.get(CGIdx(i), CGIdx(j));
4149-
4150-
if (not s0.isEmpty() and not s1.isEmpty()) {
4151-
const auto idxs0 = molinfo.getAtomsIn(CGIdx(i));
4152-
const auto idxs1 = molinfo.getAtomsIn(CGIdx(j));
4153-
4154-
for (const auto &idx0 : idxs0) {
4155-
for (const auto &idx1 : idxs1) {
4156-
QPair<AtomIdx, AtomIdx> pair =
4157-
QPair<AtomIdx, AtomIdx>(idx0, idx1);
4158-
4159-
// Make sure this is a new atom pair.
4160-
if (not recorded_pairs.contains(pair)) {
4161-
// Insert the pair and its inverse.
4162-
recorded_pairs.insert(pair);
4163-
pair = QPair<AtomIdx, AtomIdx>(idx1, idx0);
4164-
recorded_pairs.insert(pair);
4165-
4166-
const auto s0 = scl0.get(idx0, idx1);
4167-
const auto s1 = scl1.get(idx0, idx1);
4168-
4169-
if (s0.coulomb() == 1 and s0.lj() == 1 and
4170-
s1.coulomb() == 1 and s1.lj() == 1) {
4171-
// Both endstates have full 1-4 interaction (e.g. GLYCAM
4172-
// SCNB=1.0/SCEE=1.0). Write as funct=2 with explicit LJ
4173-
// parameters from state 0 (identical in both states).
4174-
if (has_ljs and has_charges) {
4175-
const auto cgidx0 = molinfo.cgAtomIdx(idx0);
4176-
const auto cgidx1 = molinfo.cgAtomIdx(idx1);
4177-
4178-
const auto &lj0 = ljs0.at(cgidx0);
4179-
const auto &lj1 = ljs0.at(cgidx1);
4180-
4181-
LJParameter lj_ij;
4182-
if (combining_rules == 2)
4183-
lj_ij = lj0.combineArithmetic(lj1);
4184-
else
4185-
lj_ij = lj0.combineGeometric(lj1);
4186-
4187-
const double qi = charges0.at(cgidx0).to(mod_electron);
4188-
const double qj = charges0.at(cgidx1).to(mod_electron);
4189-
4190-
scllines.append(
4191-
QString("%1 %2 2 1.0 %3 %4 %5 %6")
4192-
.arg(idx0 + 1, 6)
4193-
.arg(idx1 + 1, 6)
4194-
.arg(qi, 11, 'f', 6)
4195-
.arg(qj, 11, 'f', 6)
4196-
.arg(lj_ij.sigma().to(nanometer), 18, 'e', 11)
4197-
.arg(lj_ij.epsilon().to(kJ_per_mol), 18, 'e',
4198-
11));
4199-
} else {
4200-
scllines.append(QString("%1 %2 1")
4201-
.arg(idx0 + 1, 6)
4202-
.arg(idx1 + 1, 6));
4203-
}
4204-
} else if (not(s0.coulomb() == 0 and s0.lj() == 0 and
4205-
s1.coulomb() == 0 and s1.lj() == 0)) {
4206-
// This is a non-default, non-full pair (e.g. standard AMBER
4207-
// partial scaling, or a mixed perturbation).
4208-
if (fix_null_perturbable_14s) {
4209-
// get the initial and perturbed charge and LJ parameters
4210-
const auto &lj0_0 = ljs0.get(idx0);
4211-
const auto &lj0_1 = ljs1.get(idx0);
4212-
const auto &lj1_0 = ljs0.get(idx1);
4213-
const auto &lj1_1 = ljs1.get(idx1);
4214-
4215-
if (lj0_0.epsilon().value() == 0 or
4216-
lj0_1.epsilon().value() == 0 or
4217-
lj1_0.epsilon().value() == 0 or
4218-
lj1_1.epsilon().value() == 0) {
4219-
// we need to avoid having a null 1-4 LJ parameter, so
4220-
// use the non-dummy state
4221-
LJParameter lj0, lj1;
4222-
4223-
if (lj0_0.epsilon().value() == 0)
4224-
lj0 = lj0_1;
4225-
else
4226-
lj0 = lj0_0;
4227-
4228-
if (lj1_0.epsilon().value() == 0)
4229-
lj1 = lj1_1;
4230-
else
4231-
lj1 = lj1_0;
4232-
4233-
auto lj = (combining_rules == 2) ? lj0.combineArithmetic(lj1)
4234-
: lj0.combineGeometric(lj1);
4235-
4236-
double scl = s0.lj();
4237-
4238-
if (scl == 0)
4239-
scl = s1.lj();
4240-
4241-
scllines.append(
4242-
QString("%1 %2 1 %3 %4 %3 %4")
4168+
// When connectivity is available, iterate over genuine 1-4 bonded pairs
4169+
// — O(N_dihedrals). Otherwise use nonDefaultElements() on each CG pair
4170+
// — O(N_bonded). Both avoid the O(N^2) atom-pair loop.
4171+
const auto write_pair14 = [&](AtomIdx idx0, AtomIdx idx1) {
4172+
const auto s0 = scl0.get(idx0, idx1);
4173+
const auto s1 = scl1.get(idx0, idx1);
4174+
if (s0.coulomb() == 1 and s0.lj() == 1 and
4175+
s1.coulomb() == 1 and s1.lj() == 1) {
4176+
// Both end states have full 1-4 interaction (GLYCAM). Write as
4177+
// funct=2 with explicit LJ from state 0.
4178+
if (has_ljs and has_charges) {
4179+
const auto cgidx0 = molinfo.cgAtomIdx(idx0);
4180+
const auto cgidx1 = molinfo.cgAtomIdx(idx1);
4181+
const auto &lj0 = ljs0.at(cgidx0);
4182+
const auto &lj1 = ljs0.at(cgidx1);
4183+
LJParameter lj_ij;
4184+
if (combining_rules == 2)
4185+
lj_ij = lj0.combineArithmetic(lj1);
4186+
else
4187+
lj_ij = lj0.combineGeometric(lj1);
4188+
const double qi = charges0.at(cgidx0).to(mod_electron);
4189+
const double qj = charges0.at(cgidx1).to(mod_electron);
4190+
scllines.append(
4191+
QString("%1 %2 2 1.0 %3 %4 %5 %6")
4192+
.arg(idx0 + 1, 6)
4193+
.arg(idx1 + 1, 6)
4194+
.arg(qi, 11, 'f', 6)
4195+
.arg(qj, 11, 'f', 6)
4196+
.arg(lj_ij.sigma().to(nanometer), 18, 'e', 11)
4197+
.arg(lj_ij.epsilon().to(kJ_per_mol), 18, 'e', 11));
4198+
} else {
4199+
scllines.append(QString("%1 %2 1")
42434200
.arg(idx0 + 1, 6)
4244-
.arg(idx1 + 1, 6)
4245-
.arg(lj.sigma().to(nanometer), 11, 'f', 5)
4246-
.arg(scl * lj.epsilon().to(kJ_per_mol), 11, 'f',
4247-
5));
4248-
4249-
continue;
4250-
}
4251-
}
4252-
4253-
scllines.append(QString("%1 %2 1")
4254-
.arg(idx0 + 1, 6)
4255-
.arg(idx1 + 1, 6));
4256-
}
4257-
}
4258-
}
4201+
.arg(idx1 + 1, 6));
4202+
}
4203+
} else if (not(s0.coulomb() == 0 and s0.lj() == 0 and
4204+
s1.coulomb() == 0 and s1.lj() == 0)) {
4205+
// Standard partial 1-4 scaling, or a mixed perturbation.
4206+
if (fix_null_perturbable_14s) {
4207+
const auto &lj0_0 = ljs0.get(idx0);
4208+
const auto &lj0_1 = ljs1.get(idx0);
4209+
const auto &lj1_0 = ljs0.get(idx1);
4210+
const auto &lj1_1 = ljs1.get(idx1);
4211+
if (lj0_0.epsilon().value() == 0 or lj0_1.epsilon().value() == 0 or
4212+
lj1_0.epsilon().value() == 0 or lj1_1.epsilon().value() == 0) {
4213+
LJParameter lj0, lj1;
4214+
lj0 = (lj0_0.epsilon().value() == 0) ? lj0_1 : lj0_0;
4215+
lj1 = (lj1_0.epsilon().value() == 0) ? lj1_1 : lj1_0;
4216+
auto lj = (combining_rules == 2) ? lj0.combineArithmetic(lj1)
4217+
: lj0.combineGeometric(lj1);
4218+
double scl = (s0.lj() != 0) ? s0.lj() : s1.lj();
4219+
scllines.append(
4220+
QString("%1 %2 1 %3 %4 %3 %4")
4221+
.arg(idx0 + 1, 6)
4222+
.arg(idx1 + 1, 6)
4223+
.arg(lj.sigma().to(nanometer), 11, 'f', 5)
4224+
.arg(scl * lj.epsilon().to(kJ_per_mol), 11, 'f', 5));
4225+
return;
4226+
}
4227+
}
4228+
scllines.append(QString("%1 %2 1")
4229+
.arg(idx0 + 1, 6)
4230+
.arg(idx1 + 1, 6));
4231+
}
4232+
};
4233+
4234+
if (not pairs_14.isEmpty()) {
4235+
for (const auto &pair14 : pairs_14) {
4236+
if (pair14.first >= pair14.second)
4237+
continue;
4238+
write_pair14(pair14.first, pair14.second);
4239+
}
4240+
} else {
4241+
// No connectivity: iterate over non-default CG atom pair entries.
4242+
// GLYCAM-style (1,1) pairs cannot be identified here and are omitted.
4243+
for (int i = 0; i < scl0.nGroups(); ++i) {
4244+
for (int j = 0; j < scl0.nGroups(); ++j) {
4245+
for (const auto &[row, col, s0] :
4246+
scl0.get(CGIdx(i), CGIdx(j)).nonDefaultElements()) {
4247+
if (row >= col)
4248+
continue;
4249+
write_pair14(AtomIdx(row), AtomIdx(col));
42594250
}
42604251
}
42614252
}
@@ -4275,7 +4266,6 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype,
42754266
bool has_ljs = false;
42764267
bool has_charges = false;
42774268

4278-
42794269
try {
42804270
ljs = mol.property("LJ").asA<AtomLJs>();
42814271
has_ljs = true;
@@ -4288,82 +4278,66 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype,
42884278
} catch (...) {
42894279
}
42904280

4291-
// A set of recorded 1-4 pairs.
4292-
QSet<QPair<AtomIdx, AtomIdx>> recorded_pairs;
4293-
4294-
// Must record every pair that has a non-default scaling factor.
4295-
// Loop over intrascale matrix by cut-groups to avoid N^2 loop.
4296-
for (int i = 0; i < scl.nGroups(); ++i) {
4297-
for (int j = 0; j < scl.nGroups(); ++j) {
4298-
const auto s = scl.get(CGIdx(i), CGIdx(j));
4299-
4300-
if (not s.isEmpty()) {
4301-
const auto idxs0 = molinfo.getAtomsIn(CGIdx(i));
4302-
const auto idxs1 = molinfo.getAtomsIn(CGIdx(j));
4303-
4304-
for (const auto &idx0 : idxs0) {
4305-
for (const auto &idx1 : idxs1) {
4306-
QPair<AtomIdx, AtomIdx> pair =
4307-
QPair<AtomIdx, AtomIdx>(idx0, idx1);
4308-
4309-
// Make sure this is a new atom pair.
4310-
if (not recorded_pairs.contains(pair)) {
4311-
// Insert the pair and its inverse.
4312-
recorded_pairs.insert(pair);
4313-
pair = QPair<AtomIdx, AtomIdx>(idx1, idx0);
4314-
recorded_pairs.insert(pair);
4315-
4316-
const auto s = scl.get(idx0, idx1);
4317-
4318-
if (s.coulomb() == 0 and s.lj() == 0) {
4319-
// Fully excluded: don't write.
4320-
} else if (s.coulomb() == 1 and s.lj() == 1) {
4321-
// Full 1-4 interaction (e.g. GLYCAM with SCNB=1.0,
4322-
// SCEE=1.0). Must write as funct=2 with explicit LJ
4323-
// parameters because funct=1 with gen-pairs would apply
4324-
// fudgeLJ and reduce the interaction, and not listing the
4325-
// pair would give zero interaction.
4326-
if (has_ljs and has_charges) {
4327-
const auto cgidx0 = molinfo.cgAtomIdx(idx0);
4328-
const auto cgidx1 = molinfo.cgAtomIdx(idx1);
4329-
4330-
const auto &lj0 = ljs.at(cgidx0);
4331-
const auto &lj1 = ljs.at(cgidx1);
4332-
4333-
LJParameter lj_ij;
4334-
if (combining_rules == 2)
4335-
lj_ij = lj0.combineArithmetic(lj1);
4336-
else
4337-
lj_ij = lj0.combineGeometric(lj1);
4338-
4339-
const double qi = charges.at(cgidx0).to(mod_electron);
4340-
const double qj = charges.at(cgidx1).to(mod_electron);
4341-
4342-
scllines.append(
4343-
QString("%1 %2 2 1.0 %3 %4 %5 %6")
4281+
// When connectivity is available, iterate over genuine 1-4 bonded pairs
4282+
// — O(N_dihedrals). Otherwise use nonDefaultElements() on each CG pair
4283+
// — O(N_bonded). Both avoid the O(N^2) atom-pair loop.
4284+
const auto write_pair14 = [&](AtomIdx idx0, AtomIdx idx1) {
4285+
const auto s = scl.get(idx0, idx1);
4286+
if (s.coulomb() == 0 and s.lj() == 0) {
4287+
// excluded — skip
4288+
} else if (s.coulomb() == 1 and s.lj() == 1) {
4289+
// Full 1-4 interaction (GLYCAM). Write as funct=2 with explicit LJ
4290+
// parameters because funct=1 would apply fudgeLJ scaling.
4291+
if (has_ljs and has_charges) {
4292+
const auto cgidx0 = molinfo.cgAtomIdx(idx0);
4293+
const auto cgidx1 = molinfo.cgAtomIdx(idx1);
4294+
const auto &lj0 = ljs.at(cgidx0);
4295+
const auto &lj1 = ljs.at(cgidx1);
4296+
LJParameter lj_ij;
4297+
if (combining_rules == 2)
4298+
lj_ij = lj0.combineArithmetic(lj1);
4299+
else
4300+
lj_ij = lj0.combineGeometric(lj1);
4301+
const double qi = charges.at(cgidx0).to(mod_electron);
4302+
const double qj = charges.at(cgidx1).to(mod_electron);
4303+
scllines.append(
4304+
QString("%1 %2 2 1.0 %3 %4 %5 %6")
4305+
.arg(idx0 + 1, 6)
4306+
.arg(idx1 + 1, 6)
4307+
.arg(qi, 11, 'f', 6)
4308+
.arg(qj, 11, 'f', 6)
4309+
.arg(lj_ij.sigma().to(nanometer), 18, 'e', 11)
4310+
.arg(lj_ij.epsilon().to(kJ_per_mol), 18, 'e', 11));
4311+
} else {
4312+
// Fall back to funct=1; energy will be wrong if fudgeLJ != 1.0.
4313+
scllines.append(QString("%1 %2 1")
4314+
.arg(idx0 + 1, 6)
4315+
.arg(idx1 + 1, 6));
4316+
}
4317+
} else {
4318+
// Standard partial 1-4 scaling (e.g. AMBER). Write as funct=1.
4319+
scllines.append(QString("%1 %2 1")
43444320
.arg(idx0 + 1, 6)
4345-
.arg(idx1 + 1, 6)
4346-
.arg(qi, 11, 'f', 6)
4347-
.arg(qj, 11, 'f', 6)
4348-
.arg(lj_ij.sigma().to(nanometer), 18, 'e', 11)
4349-
.arg(lj_ij.epsilon().to(kJ_per_mol), 18, 'e',
4350-
11));
4351-
} else {
4352-
// Fall back to funct=1; the energy will be wrong if
4353-
// fudgeLJ != 1.0, but we have no LJ parameters to use.
4354-
scllines.append(QString("%1 %2 1")
4355-
.arg(idx0 + 1, 6)
4356-
.arg(idx1 + 1, 6));
4357-
}
4358-
} else {
4359-
// Standard partial 1-4 (e.g. fudgeQQ/fudgeLJ): write as
4360-
// funct=1.
4361-
scllines.append(QString("%1 %2 1")
4362-
.arg(idx0 + 1, 6)
4363-
.arg(idx1 + 1, 6));
4364-
}
4365-
}
4366-
}
4321+
.arg(idx1 + 1, 6));
4322+
}
4323+
};
4324+
4325+
if (not pairs_14.isEmpty()) {
4326+
for (const auto &pair14 : pairs_14) {
4327+
if (pair14.first >= pair14.second)
4328+
continue;
4329+
write_pair14(pair14.first, pair14.second);
4330+
}
4331+
} else {
4332+
// No connectivity: iterate over non-default CG atom pair entries.
4333+
// GLYCAM-style (1,1) pairs cannot be identified here and are omitted.
4334+
for (int i = 0; i < scl.nGroups(); ++i) {
4335+
for (int j = 0; j < scl.nGroups(); ++j) {
4336+
for (const auto &[row, col, s] :
4337+
scl.get(CGIdx(i), CGIdx(j)).nonDefaultElements()) {
4338+
if (row >= col)
4339+
continue;
4340+
write_pair14(AtomIdx(row), AtomIdx(col));
43674341
}
43684342
}
43694343
}

0 commit comments

Comments
 (0)