Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Feature grotop14s #173

Merged
merged 3 commits into from
Mar 8, 2024
Merged

Feature grotop14s #173

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
65bd51d
Added a flag to prevent LJ 14 terms being scaled to zero if this (non…
chryswoods Mar 6, 2024
adbcd03
Moved `fix_null_perturbable_14s` to be a molecule property, rather th…
chryswoods Mar 7, 2024
4a4cbe6
Merge branch 'devel' into feature_grotop14s
chryswoods Mar 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 73 additions & 5 deletions corelib/src/libs/SireIO/grotop.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3157,7 +3157,8 @@ static QStringList writeAtomTypes(QMap<QPair<int, QString>, GroMolType> &moltyps
}

/** 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)
static QStringList writeMolType(const QString &name, const GroMolType &moltype, const Molecule &mol,
bool uses_parallel)
{
QStringList lines;

Expand Down Expand Up @@ -4036,9 +4037,33 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype,
return;
}

AtomLJs ljs0;
AtomLJs ljs1;

try
{
ljs0 = mol.property("LJ0").asA<AtomLJs>();
}
catch (...)
{
}

try
{
ljs1 = mol.property("LJ1").asA<AtomLJs>();
}
catch (...)
{
}

// A set of recorded 1-4 pairs.
QSet<QPair<AtomIdx, AtomIdx>> 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<BooleanProperty>().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)
Expand Down Expand Up @@ -4076,6 +4101,47 @@ static QStringList writeMolType(const QString &name, const GroMolType &moltype,
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));
}
}
Expand Down Expand Up @@ -4243,7 +4309,8 @@ static QStringList writeMolTypes(const QMap<QPair<int, QString>, GroMolType> &mo
for (int i = r.begin(); i < r.end(); ++i)
{
QStringList typlines =
::writeMolType(keys[i].second, moltyps[keys[i]], examples[keys[i]], uses_parallel);
::writeMolType(keys[i].second, moltyps[keys[i]], examples[keys[i]],
uses_parallel);

QMutexLocker lkr(&mutex);
typs.insert(keys[i].second, typlines);
Expand All @@ -4254,7 +4321,8 @@ static QStringList writeMolTypes(const QMap<QPair<int, QString>, GroMolType> &mo
for (auto it = moltyps.constBegin(); it != moltyps.constEnd(); ++it)
{
typs.insert(it.key().second,
::writeMolType(it.key().second, it.value(), examples[it.key()], uses_parallel));
::writeMolType(it.key().second, it.value(), examples[it.key()],
uses_parallel));
}
}

Expand Down Expand Up @@ -4538,8 +4606,8 @@ GroTop::GroTop(const SireSystem::System &system, const PropertyMap &map)
// the molecules
lines += ::writeAtomTypes(idx_name_to_mtyp, idx_name_to_example, ffield, map);

// now convert these into text lines that can be written as the file
lines += ::writeMolTypes(idx_name_to_mtyp, idx_name_to_example, usesParallel(), isSorted);
lines += ::writeMolTypes(idx_name_to_mtyp, idx_name_to_example, usesParallel(),
isSorted);

// now write the system part
lines += ::writeSystem(system.name(), mol_to_moltype);
Expand Down
Loading