diff --git a/README.rst b/README.rst index eea059fd2..6c296fbd0 100644 --- a/README.rst +++ b/README.rst @@ -45,7 +45,7 @@ To create a new environment: mamba create -n openbiosim "python<3.11" mamba activate openbiosim - mamba install -c openbiosim sire + mamba install -c conda-forge -c openbiosim sire To install the latest development version you can use: @@ -53,7 +53,7 @@ To install the latest development version you can use: mamba create -n openbiosim-dev "python<3.11" mamba activate openbiosim-dev - mamba install -c openbiosim/label/dev sire + mamba install -c conda-forge -c openbiosim/label/dev sire However, as you are here, it is likely you want to download the latest, greatest version of the code, which you will need to compile. To compile diff --git a/corelib/src/libs/SireCAS/lambdaschedule.cpp b/corelib/src/libs/SireCAS/lambdaschedule.cpp index 78b9bed63..168a0a158 100644 --- a/corelib/src/libs/SireCAS/lambdaschedule.cpp +++ b/corelib/src/libs/SireCAS/lambdaschedule.cpp @@ -432,14 +432,23 @@ void LambdaSchedule::clear() this->constant_values = Values(); } +/** Append a morph stage onto this schedule. The morph stage is a + * standard stage that scales each forcefield parameter by + * (1-:lambda:).initial + :lambda:.final + */ +void LambdaSchedule::addMorphStage(const QString &name) +{ + this->addStage(name, (this->lam() * this->final()) + + ((1 - this->lam()) * this->initial())); +} + /** Append a morph stage onto this schedule. The morph stage is a * standard stage that scales each forcefield parameter by * (1-:lambda:).initial + :lambda:.final */ void LambdaSchedule::addMorphStage() { - this->addStage("morph", (this->lam() * this->final()) + - ((1 - this->lam()) * this->initial())); + this->addMorphStage("morph"); } /** Sandwich the current set of stages with a charge-descaling and @@ -450,7 +459,9 @@ void LambdaSchedule::addMorphStage() * by :gamma:. A final charge-rescaling stage is then appended that * scales the charge parameter from :gamma:.final to final. */ -void LambdaSchedule::addChargeScaleStages(double scale) +void LambdaSchedule::addChargeScaleStages(const QString &decharge_name, + const QString &recharge_name, + double scale) { auto scl = this->setConstant("γ", scale); @@ -462,11 +473,24 @@ void LambdaSchedule::addChargeScaleStages(double scale) } // now prepend the decharging stage, and append the recharging stage - this->prependStage("decharge", this->initial()); - this->appendStage("recharge", this->final()); + this->prependStage(decharge_name, this->initial()); + this->appendStage(recharge_name, this->final()); + + this->setEquation(decharge_name, "charge", (1.0 - ((1.0 - scl) * this->lam())) * this->initial()); + this->setEquation(recharge_name, "charge", (1.0 - ((1.0 - scl) * (1.0 - this->lam()))) * this->final()); +} - this->setEquation("decharge", "charge", (1.0 - ((1.0 - scl) * this->lam())) * this->initial()); - this->setEquation("recharge", "charge", (1.0 - ((1.0 - scl) * (1.0 - this->lam()))) * this->final()); +/** Sandwich the current set of stages with a charge-descaling and + * a charge-scaling stage. This prepends a charge-descaling stage + * that scales the charge parameter down from `initial` to + * :gamma:.initial (where :gamma:=`scale`). The charge parameter in all of + * the exising stages in this schedule are then multiplied + * by :gamma:. A final charge-rescaling stage is then appended that + * scales the charge parameter from :gamma:.final to final. + */ +void LambdaSchedule::addChargeScaleStages(double scale) +{ + this->addChargeScaleStages("decharge", "recharge", scale); } /** Prepend a stage called 'name' which uses the passed 'equation' diff --git a/corelib/src/libs/SireCAS/lambdaschedule.h b/corelib/src/libs/SireCAS/lambdaschedule.h index fc432e0d2..2026899f9 100644 --- a/corelib/src/libs/SireCAS/lambdaschedule.h +++ b/corelib/src/libs/SireCAS/lambdaschedule.h @@ -117,7 +117,12 @@ namespace SireCAS const SireCAS::Expression &equation); void addMorphStage(); + void addMorphStage(const QString &name); + void addChargeScaleStages(double scale = 0.2); + void addChargeScaleStages(const QString &decharge_name, + const QString &recharge_name, + double scale = 0.2); void setEquation(const QString &stage, const QString &lever, diff --git a/corelib/src/libs/SireIO/moleculeparser.cpp b/corelib/src/libs/SireIO/moleculeparser.cpp index 09c585e05..63f4b92bd 100644 --- a/corelib/src/libs/SireIO/moleculeparser.cpp +++ b/corelib/src/libs/SireIO/moleculeparser.cpp @@ -2323,17 +2323,24 @@ QStringList MoleculeParser::write(const System &system, const QString &filename, QStringList filenames; QStringList fileformats; - const auto format_property = map["fileformat"]; - - if (format_property.hasValue()) + if (map.specified("fileformat")) { - try + const auto format_property = map["fileformat"]; + + if (format_property.hasSource()) { - fileformats = format_property.value().asA().toString().split(","); + fileformats = format_property.source().split(","); } - catch (...) + else { - fileformats.append(format_property.value().asA().formatName()); + try + { + fileformats = format_property.value().asA().toString().split(","); + } + catch (...) + { + fileformats.append(format_property.value().asA().formatName()); + } } const auto fileinfo = QFileInfo(filename); @@ -2353,7 +2360,7 @@ QStringList MoleculeParser::write(const System &system, const QString &filename, // we need to find the format from the system try { - fileformats = system.property(format_property).asA().toString().split(","); + fileformats = system.property(map["fileformat"]).asA().toString().split(","); } catch (...) { @@ -2404,17 +2411,24 @@ QStringList MoleculeParser::write(const System &system, const QStringList &files QStringList filenames; QStringList fileformats; - const auto format_property = map["fileformat"]; - - if (format_property.hasValue()) + if (map.specified("fileformat")) { - try + const auto format_property = map["fileformat"]; + + if (format_property.hasSource()) { - fileformats = format_property.value().asA().toString().split(","); + fileformats = format_property.source().split(","); } - catch (...) + else { - fileformats.append(format_property.value().asA().formatName()); + try + { + fileformats = format_property.value().asA().toString().split(","); + } + catch (...) + { + fileformats.append(format_property.value().asA().formatName()); + } } if (files.count() != fileformats.count()) @@ -2449,7 +2463,7 @@ QStringList MoleculeParser::write(const System &system, const QStringList &files // we may need to find the format from the system try { - fileformats = system.property(format_property).asA().toString().split(","); + fileformats = system.property(map["fileformat"]).asA().toString().split(","); } catch (...) { diff --git a/corelib/src/libs/SireMaths/energytrajectory.cpp b/corelib/src/libs/SireMaths/energytrajectory.cpp index fa53fbdf8..ee68830a0 100644 --- a/corelib/src/libs/SireMaths/energytrajectory.cpp +++ b/corelib/src/libs/SireMaths/energytrajectory.cpp @@ -52,6 +52,7 @@ QDataStream &operator<<(QDataStream &ds, const EnergyTrajectory &etraj) SharedDataStream sds(ds); sds << etraj.time_values << etraj.energy_values + << etraj.label_values << etraj.props << static_cast(etraj); return ds; @@ -65,7 +66,7 @@ QDataStream &operator>>(QDataStream &ds, EnergyTrajectory &etraj) { SharedDataStream sds(ds); - sds >> etraj.time_values >> etraj.energy_values >> static_cast(etraj); + sds >> etraj.time_values >> etraj.energy_values >> etraj.label_values >> etraj.props >> static_cast(etraj); } else throw version_error(v, "1", r_etraj, CODELOC); @@ -80,7 +81,8 @@ EnergyTrajectory::EnergyTrajectory() EnergyTrajectory::EnergyTrajectory(const EnergyTrajectory &other) : ConcreteProperty(other), - time_values(other.time_values), energy_values(other.energy_values) + time_values(other.time_values), energy_values(other.energy_values), + label_values(other.label_values), props(other.props) { } @@ -94,6 +96,8 @@ EnergyTrajectory &EnergyTrajectory::operator=(const EnergyTrajectory &other) { time_values = other.time_values; energy_values = other.energy_values; + label_values = other.label_values; + props = other.props; Property::operator=(other); } @@ -103,7 +107,9 @@ EnergyTrajectory &EnergyTrajectory::operator=(const EnergyTrajectory &other) bool EnergyTrajectory::operator==(const EnergyTrajectory &other) const { return time_values == other.time_values and - energy_values == other.energy_values; + energy_values == other.energy_values and + label_values == other.label_values and + props == other.props; } bool EnergyTrajectory::operator!=(const EnergyTrajectory &other) const @@ -137,18 +143,32 @@ QString EnergyTrajectory::toString() const auto keys = this->energy_values.keys(); keys.sort(); - keys.insert(0, "time"); - parts.append(keys.join("\t")); + auto label_keys = this->label_values.keys(); + label_keys.sort(); + + auto headers = label_keys + keys; + + headers.insert(0, "time"); + + parts.append(headers.join("\t")); if (n <= 10) { for (int i = 0; i < n; ++i) { const auto vals = this->get(i, GeneralUnit(picosecond), GeneralUnit(kcal_per_mol)); + const auto label_vals = this->getLabels(i); QStringList v; + v.append(QString::number(vals["time"])); + + for (const auto &key : label_keys) + { + v.append(label_vals[key]); + } + for (const auto &key : keys) { v.append(QString::number(vals[key])); @@ -162,9 +182,17 @@ QString EnergyTrajectory::toString() const for (int i = 0; i < 5; ++i) { const auto vals = this->get(i, GeneralUnit(picosecond), GeneralUnit(kcal_per_mol)); + const auto label_vals = this->getLabels(i); QStringList v; + v.append(QString::number(vals["time"])); + + for (const auto &key : label_keys) + { + v.append(label_vals[key]); + } + for (const auto &key : keys) { v.append(QString::number(vals[key])); @@ -178,9 +206,17 @@ QString EnergyTrajectory::toString() const for (int i = n - 5; i < n; ++i) { const auto vals = this->get(i, GeneralUnit(picosecond), GeneralUnit(kcal_per_mol)); + const auto label_vals = this->getLabels(i); QStringList v; + v.append(QString::number(vals["time"])); + + for (const auto &key : label_keys) + { + v.append(label_vals[key]); + } + for (const auto &key : keys) { v.append(QString::number(vals[key])); @@ -190,16 +226,17 @@ QString EnergyTrajectory::toString() const } } - return QObject::tr("SelectorMol( size=%1\n%2\n)").arg(n).arg(parts.join("\n")); + return QObject::tr("EnergyTrajectory( size=%1\n%2\n)").arg(n).arg(parts.join("\n")); } /** Set the energies at time 'time' to the components contained - * in 'energies' + * in 'energies', and the labels to those contained in 'labels' */ void EnergyTrajectory::set(const GeneralUnit &time, - const QHash &energies) + const QHash &energies, + const QHash &labels) { - if (energies.isEmpty()) + if (energies.isEmpty() and labels.isEmpty()) return; auto t = Time(time); @@ -218,6 +255,21 @@ void EnergyTrajectory::set(const GeneralUnit &time, auto nrg = MolarEnergy(it.value()); } + for (auto it = labels.constBegin(); it != labels.constEnd(); ++it) + { + // make sure all of the values are valid + if (it.key() == "time") + { + throw SireError::invalid_key(QObject::tr( + "You cannot call a label component 'time'. This name " + "is reserved for the time component."), + CODELOC); + } + } + + // we won't check for duplication in the energies and labels as these + // can be resolved on printout + // make sure we have populated all of the columns - this // populates with NaNs if the column doesn't exist QVector null_column; @@ -235,6 +287,21 @@ void EnergyTrajectory::set(const GeneralUnit &time, } } + QVector null_label_column; + + for (const auto &key : labels.keys()) + { + if (not this->label_values.contains(key)) + { + if (null_label_column.isEmpty() and not time_values.isEmpty()) + { + null_label_column = QVector(time_values.count()); + } + + this->label_values.insert(key, null_label_column); + } + } + // find the index of the time, and add it as needed. // Times are sorted, so start from the end and work towards // the beginning (as we will normally be adding energies onto @@ -266,6 +333,12 @@ void EnergyTrajectory::set(const GeneralUnit &time, { it.value().append(NAN); } + + for (auto it = this->label_values.begin(); + it != this->label_values.end(); ++it) + { + it.value().append(QString()); + } } else if (must_create) { @@ -276,6 +349,12 @@ void EnergyTrajectory::set(const GeneralUnit &time, { it.value().insert(idx, NAN); } + + for (auto it = this->label_values.begin(); + it != this->label_values.end(); ++it) + { + it.value().insert(idx, QString()); + } } for (auto it = energies.constBegin(); @@ -283,6 +362,21 @@ void EnergyTrajectory::set(const GeneralUnit &time, { this->energy_values[it.key()][idx] = MolarEnergy(it.value()).value(); } + + for (auto it = labels.constBegin(); + it != labels.constEnd(); ++it) + { + this->label_values[it.key()][idx] = it.value(); + } +} + +/** Set the energies at time 'time' to the components contained + * in 'energies' + */ +void EnergyTrajectory::set(const GeneralUnit &time, + const QHash &energies) +{ + this->set(time, energies, QHash()); } /** Return whether or not this is empty */ @@ -362,12 +456,73 @@ QHash EnergyTrajectory::get(int i, return vals; } +/** Return the labels at the ith row */ +QHash EnergyTrajectory::getLabels(int i) const +{ + i = SireID::Index(i).map(this->count()); + + QHash vals; + vals.reserve(this->label_values.count()); + + for (auto it = this->label_values.constBegin(); + it != this->label_values.constEnd(); + ++it) + { + vals.insert(it.key(), it.value()[i]); + } + + return vals; +} + +/** Return the labels at the ith row as numbers */ +QHash EnergyTrajectory::getLabelsAsNumbers(int i) const +{ + i = SireID::Index(i).map(this->count()); + + QHash vals; + vals.reserve(this->label_values.count()); + + for (auto it = this->label_values.constBegin(); + it != this->label_values.constEnd(); + ++it) + { + if (it.value()[i].isEmpty()) + { + vals.insert(it.key(), NAN); + continue; + } + + double val; + bool ok = false; + + val = it.value()[i].toDouble(&ok); + + if (not ok) + throw SireError::invalid_cast(QObject::tr( + "Cannot convert label '%1' at index %2 to a number. Value is %3.") + .arg(it.key()) + .arg(i) + .arg(it.value()[i]), + CODELOC); + + vals.insert(it.key(), val); + } + + return vals; +} + /** Return all of the energy keys */ QStringList EnergyTrajectory::keys() const { return this->energy_values.keys(); } +/** Return all of the label keys */ +QStringList EnergyTrajectory::labelKeys() const +{ + return this->label_values.keys(); +} + /** Return all of the time values (the time column). This is * sorted from earliest to latest time, and is in the default internal * unit @@ -377,7 +532,70 @@ QVector EnergyTrajectory::times() const return this->time_values; } -/** Return all of the energy values for the passed key (the energy-key column). +/** Return all of the label values for the passed key (the label-key column). + * This is in the same order as the times. + */ +QVector EnergyTrajectory::labels(const QString &key) const +{ + auto it = this->label_values.constFind(key); + + if (it == this->label_values.constEnd()) + throw SireError::invalid_key(QObject::tr( + "There is no label component with key '%1'. Valid " + "keys are [ %2 ].") + .arg(key) + .arg(this->label_values.keys().join(", ")), + CODELOC); + + return it.value(); +} + +/** Return all of the label values for the passed key (the label-key column) + * as numbers. This is in the same order as the times. + */ +QVector EnergyTrajectory::labelsAsNumbers(const QString &key) const +{ + auto it = this->label_values.constFind(key); + + if (it == this->label_values.constEnd()) + throw SireError::invalid_key(QObject::tr( + "There is no label component with key '%1'. Valid " + "keys are [ %2 ].") + .arg(key) + .arg(this->label_values.keys().join(", ")), + CODELOC); + + const auto &l = it.value(); + + QVector vals; + vals.reserve(l.count()); + + for (const auto &label : l) + { + if (label.isEmpty()) + { + vals.append(NAN); + continue; + } + + bool ok = false; + double val = label.toDouble(&ok); + + if (not ok) + throw SireError::invalid_cast(QObject::tr( + "Cannot convert label '%1' to a number. Value is %2.") + .arg(key) + .arg(label), + CODELOC); + + vals.append(val); + } + + return vals; +} + +/** Return all of the + * energy values for the passed key (the energy-key column). * This is in the same order as the times, and is in the default internal * unit */ @@ -434,3 +652,49 @@ QVector EnergyTrajectory::energies(const QString &key, return e; } + +/** Set a property on the trajectory. This is used to store additional + * information about the trajectory, such as the simulation temperature + */ +void EnergyTrajectory::setProperty(const QString &key, const SireBase::Property &value) +{ + props.setProperty(key, value); +} + +/** Get a property on the trajectory. This could be additional + * information about the trajectory, such as the simulation temperature + */ +const SireBase::Property &EnergyTrajectory::property(const SireBase::PropertyName &key) const +{ + return props.property(key); +} + +/** Return whether or not the trajectory has a property with the passed key */ +bool EnergyTrajectory::hasProperty(const SireBase::PropertyName &key) +{ + return props.hasProperty(key); +} + +/** Return all of the properties on the trajectory */ +const SireBase::Properties &EnergyTrajectory::properties() const +{ + return props; +} + +/** Return all of the property keys on the trajectory */ +QStringList EnergyTrajectory::propertyKeys() const +{ + return props.propertyKeys(); +} + +/** Remove a property from the trajectory */ +void EnergyTrajectory::removeProperty(const QString &key) +{ + props.removeProperty(key); +} + +/** Clear all of the properties on the trajectory */ +void EnergyTrajectory::clearProperties() +{ + props = Properties(); +} diff --git a/corelib/src/libs/SireMaths/energytrajectory.h b/corelib/src/libs/SireMaths/energytrajectory.h index 2b7152206..66509a970 100644 --- a/corelib/src/libs/SireMaths/energytrajectory.h +++ b/corelib/src/libs/SireMaths/energytrajectory.h @@ -33,6 +33,7 @@ #include "SireUnits/dimensions.h" #include "SireBase/property.h" +#include "SireBase/properties.h" SIRE_BEGIN_HEADER @@ -87,24 +88,54 @@ namespace SireMaths const SireUnits::Dimension::GeneralUnit &time_unit, const SireUnits::Dimension::GeneralUnit &energy_unit) const; + QHash getLabels(int i) const; + QHash getLabelsAsNumbers(int i) const; + void set(const SireUnits::Dimension::GeneralUnit &time, const QHash &energies); + void set(const SireUnits::Dimension::GeneralUnit &time, + const QHash &energies, + const QHash &labels); + QStringList keys() const; + QStringList labelKeys() const; QVector times() const; QVector energies(const QString &key) const; + QVector labels(const QString &key) const; + QVector labelsAsNumbers(const QString &key) const; QVector times(const SireUnits::Dimension::GeneralUnit &time_unit) const; QVector energies(const QString &key, const SireUnits::Dimension::GeneralUnit &energy_unit) const; + void setProperty(const QString &key, const SireBase::Property &value); + + const SireBase::Property &property(const SireBase::PropertyName &key) const; + + bool hasProperty(const SireBase::PropertyName &key); + + const SireBase::Properties &properties() const; + + QStringList propertyKeys() const; + + void removeProperty(const QString &key); + + void clearProperties(); + private: /** All of the time values */ QVector time_values; /** All of the energy values */ QHash> energy_values; + + /** All of the label values */ + QHash> label_values; + + /** Additional properties for the table (e.g. temperature) */ + SireBase::Properties props; }; } diff --git a/corelib/src/libs/SireMol/findmcs.cpp b/corelib/src/libs/SireMol/findmcs.cpp index f62a0dc3c..0a3c841de 100644 --- a/corelib/src/libs/SireMol/findmcs.cpp +++ b/corelib/src/libs/SireMol/findmcs.cpp @@ -30,15 +30,58 @@ #include "sireglobal.h" +#ifdef WIN32 + +#include "evaluator.h" + +#include "SireError/errors.h" + +using namespace SireMol; +using namespace SireBase; +using namespace SireUnits::Dimension; + +QHash Evaluator::findMCS(const MoleculeView &other, const AtomMatcher &matcher, const Time &timeout, + bool match_light_atoms, const PropertyMap &map0, const PropertyMap &map1, + int min_heavy_protons, bool verbose) const +{ + throw SireError::unsupported(QObject::tr( + "MCS searches are not supported on Windows because of incompatibilities " + "between boost::describe and the compiler used to build sire on Windows."), + CODELOC); + return QHash(); +} + +QVector> Evaluator::findMCSmatches(const MoleculeView &other, const AtomMatcher &matcher, + const Time &timeout, bool match_light_atoms, + const PropertyMap &map0, const PropertyMap &map1, + int min_heavy_protons, bool verbose) const +{ + throw SireError::unsupported(QObject::tr( + "MCS searches are not supported on Windows because of incompatibilities " + "between boost::describe and the compiler used to build sire on Windows."), + CODELOC); + return QVector>(); +} + +#else + #ifdef __clang__ #pragma clang diagnostic push #pragma clang diagnostic ignored "-Wdeprecated-declarations" #pragma clang diagnostic ignored "-Wdeprecated-builtins" #endif -#include -#include -#include + #include + +// this header causes a compile error on Windows because of missing +// boost_base_descriptor_fn, 'boost_public_member_descriptor_fn +// boost_protected_member_descriptor_fn and +// boost_private_member_descriptor_fn. This is because +// it includes boost/describe.hpp which seems not to build... +#include + +#include + #ifdef __clang__ #pragma clang diagnostic pop #endif @@ -1001,3 +1044,5 @@ QVector> Evaluator::findMCSmatches(const MoleculeView &o return pvt_findMCSmatches(*this, other, matcher, timeout, false, map0, map1, false, true, min_heavy_protons, verbose); } + +#endif // ifdef WIN32 diff --git a/corelib/src/libs/SireMove/ensemble.cpp b/corelib/src/libs/SireMove/ensemble.cpp index 8555a1fe4..ea3c3b960 100644 --- a/corelib/src/libs/SireMove/ensemble.cpp +++ b/corelib/src/libs/SireMove/ensemble.cpp @@ -358,7 +358,7 @@ QString Ensemble::toString() const if (this->isConstantTemperature()) { - parts.append(QObject::tr("temperature = %1 C").arg(this->temperature().to(Celsius()))); + parts.append(QObject::tr("temperature = %1").arg(Celsius(this->temperature()).toString())); } if (parts.isEmpty()) diff --git a/corelib/src/libs/SireSearch/ast.h b/corelib/src/libs/SireSearch/ast.h index 0fcd29f0b..0b4a97e61 100644 --- a/corelib/src/libs/SireSearch/ast.h +++ b/corelib/src/libs/SireSearch/ast.h @@ -50,17 +50,8 @@ SIRE_BEGIN_HEADER #define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS #define BOOST_MPL_LIMIT_LIST_SIZE 40 -#include -#include - -#include -#include -#include -#include -#include - #include -#include +#include #include diff --git a/corelib/src/libs/SireUnits/ast.h b/corelib/src/libs/SireUnits/ast.h index fffec8003..506dcf39d 100644 --- a/corelib/src/libs/SireUnits/ast.h +++ b/corelib/src/libs/SireUnits/ast.h @@ -35,17 +35,7 @@ SIRE_BEGIN_HEADER -#include -#include - -#include -#include -#include -#include -#include - #include -#include #include diff --git a/corelib/src/libs/SireUnits/dimensions.cpp b/corelib/src/libs/SireUnits/dimensions.cpp index de224090b..0721c93ca 100644 --- a/corelib/src/libs/SireUnits/dimensions.cpp +++ b/corelib/src/libs/SireUnits/dimensions.cpp @@ -1,6 +1,7 @@ #include "dimensions.h" #include "temperature.h" +#include "generalunit.h" namespace SireUnits { diff --git a/doc/source/blogposts.rst b/doc/source/blogposts.rst new file mode 100644 index 000000000..ab856d8c8 --- /dev/null +++ b/doc/source/blogposts.rst @@ -0,0 +1,45 @@ +======= +Diaries +======= + +Development of :mod:`sire` is an ongoing process involving many people. +This section of the documentation is a collection of diary-style +blog posts that describe the development of :mod:`sire` and +new features as they are added. + +What's in a file? Using sire search +----------------------------------- + +* `What's in a file? Using sire search `__ + +In this first post, written around the 2023.1 release, +we describe the :doc:`sire search ` +functionality that enables :mod:`sire` to act like a molecular information +search engine. + +Working with RDKit +------------------ + +* `Working with RDKit `__ + +In this second post, written around the 2023.2 release, +we show how :doc:`integration with RDKit ` +enables :mod:`sire` to work with SMILES strings and RDKit molecules. + +Working with OpenMM +------------------- + +* `Working with OpenMM `__ + +In this third post, also written around the 2023.2 release, +we show how :doc:`integration with OpenMM ` +enables :mod:`sire` to perform GPU-accelerated molecular dynamics simulations. + +GPU-accelerated free energies +----------------------------- + +* `GPU-accelerated free energies `__ + +In this fourth post, written around the 2023.4 release, +we show how the OpenMM integration has been extended to +enable the :doc:`calculation of free energies `. diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 80e55db12..126be7082 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -12,8 +12,8 @@ Development was migrated into the `OpenBioSim `__ organisation on `GitHub `__. -`2023.4.0 `__ - September 2023 ----------------------------------------------------------------------------------------------- +`2023.4.0 `__ - October 2023 +-------------------------------------------------------------------------------------------- * Added ``closest`` and ``furthest`` keywords to enable searching for the n closest or furthest views. This is very general, and is described in the @@ -36,12 +36,6 @@ organisation on `GitHub `__. Added a new ``sire.morph`` module that includes functions that should make it easier to set up, view and control morphs (perturbations). -* Added an ``EnergyTrajectory`` class that lets us record the energy - trajectory along a dynamics simulation. This includes recording energies - at different λ-windows to that being simulated, thereby providing - the raw data for free energy calculations. By default the - ``EnergyTrajectory`` is returned to the user as a pandas DataFrame. - * Forced all new-style modules to import when `sr.use_new_api()` is called. This will make it easier to use sire with multiprocessing. @@ -89,17 +83,39 @@ organisation on `GitHub `__. Restraints can be named, meaning that you can scale different restraints at different stages and by different values across the λ-coordinate. -* Added support for one or more "permanent" distance restraints. These are - distance restraints that are always applied, and are never scaled by λ. +* Added an :class:`~sire.maths.EnergyTrajectory` class that lets us record the + energy trajectory along a dynamics simulation. This includes recording + energies at different λ-windows to that being simulated, thereby providing + the raw data for free energy calculations. By default the + ``EnergyTrajectory`` is returned to the user as a pandas DataFrame. + +* Added the ability to export an :class:`~sire.maths.EnergyTrajectory` as + an alchemlyb-compatible data frame. Added :func:`sire.morph.to_alchemlyb` + to convert lots of ``EnergyTrajectory`` objects (or files containing + s3 streams) into a single alchemlyb-compatible data frame that is + ready for analysis. You can now calculate relative hydration and binding + free energies and analyse the results using alchemlyb. This is documented + in the :doc:`tutorial `. + +* Added a :func:`sire.morph.repartition_hydrogen_masses` to make it easier to + repartition hydrogen masses during alchemical free energy simulations. + Set the default mass factor to 1.5 to support a 4 fs timestep with the + default ``LangevinMiddleIntegrator``. + +* Added support for an Andersen thermostat in the OpenMM dynamics layer. + +* Added support for scaling intramolecular non-bonded scale factors to the + ``LambdaLever``, so that we have rudimentary support for perturbations + that involve bond breaking and forming. + +* Added support to somd for one or more "permanent" distance restraints. These + are distance restraints that are always applied, and are never scaled by λ. This allows the release of all other distance restraints to a single harmonic or flat-bottomed restraint. When the ligand is fully decoupled, the free energy of release of the single remaining restraint can be computed without simulation. See for more details. -* Please add the changelog entry for your PR here. We will add the link to your PR - during the code review :-) - `2023.3.2 `__ - September 2023 ---------------------------------------------------------------------------------------------- @@ -143,9 +159,9 @@ organisation on `GitHub `__. is now roughly O(N), using a hash to lookup at the molecule level, so that we don't have to test individual atoms. -* Fixed :module:`sire.wrapper.tools.standardstatecorrection`. This stopped working after - the commit https://github.com/OpenBioSim/sire/commit/e2e370940894315838fb8f65e141baaf07050ce0, - because not all required changes were included. +* Fixed ``StandardStateCorrection``. This stopped working after + the commit https://github.com/OpenBioSim/sire/commit/e2e370940894315838fb8f65e141baaf07050ce0, + because not all required changes were included. * Fix for crash when not passing a map to the SelectorImproper constructor diff --git a/doc/source/cheatsheet/openmm.rst b/doc/source/cheatsheet/openmm.rst index 95bed8f23..d5345681c 100644 --- a/doc/source/cheatsheet/openmm.rst +++ b/doc/source/cheatsheet/openmm.rst @@ -76,6 +76,11 @@ Available keys and allowable values are listed below. +---------------------------+----------------------------------------------------------+ | Key | Valid values | +===========================+==========================================================+ +| com_reset_frequency | The frequency at which the ``CMMotionRemover`` acts to | +| | remove center of mass relative motion. If this is not | +| | set (the default) then center of mass motion is not | +| | removed. | ++---------------------------+----------------------------------------------------------+ | constraint | Type of constraint to use for bonds and/or angles. | | | Valid strings are ``none``, ``h-bonds``, ``bonds``, | | | ``h-bonds-h-angles`` and ``bonds-h-angles``. | @@ -109,7 +114,7 @@ Available keys and allowable values are listed below. | integrator | The MD integrator to use, e.g. | | | ``verlet``, ``leapfrog``, ``langevin``, | | | ``langevin_middle``, ``nose_hoover``, | -| | ``brownian`` | +| | ``brownian``, ``andersen`` | +---------------------------+----------------------------------------------------------+ | friction | Friction value for the integrator, in inverse time, e.g. | | | ``5.0 / sr.units.picosecond`` | diff --git a/doc/source/features.rst b/doc/source/features.rst index 65f2a91a2..5fa08b231 100644 --- a/doc/source/features.rst +++ b/doc/source/features.rst @@ -5,6 +5,8 @@ Features :mod:`sire` is a powerful Python/C++ module for loading and manipulating molecular (predominantly biomolecular) systems. +* :doc:`Calculate alchemical free energies ` using + GPU-accelerated molecular dynamics simulations. * :doc:`Load and save ` molecules from a number of :doc:`molecular file formats `. diff --git a/doc/source/index.rst b/doc/source/index.rst index 0af92d3c5..ce6e90c6a 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -30,9 +30,10 @@ Tutorial Blog Posts ========== -* `What's in a file? Using sire search `__ -* `Working with RDKit `__ -* `Working with OpenMM `__ +.. toctree:: + :maxdepth: 2 + + blogposts Detailed Guides =============== diff --git a/doc/source/install.rst b/doc/source/install.rst index 6f1b42fd2..4c5e2f2bc 100644 --- a/doc/source/install.rst +++ b/doc/source/install.rst @@ -108,28 +108,30 @@ And then... Install sire into a new environment We recommend that :mod:`sire` is installed into a new (clean) environment. This minimises the risk of failures caused by incompatible dependencies. -Sire is currently packaged for Python 3.8, 3.9 and Python 3.10. We will start -by creating a Python 3.10 environment that we will call ``openbiosim``. +Sire is currently packaged for Python 3.9, 3.10 and Python 3.11. We will start +by creating a Python 3.11 environment that we will call ``openbiosim``. .. code-block:: bash - $ mamba create -n openbiosim "python<3.11" + $ mamba create -n openbiosim "python<3.12" .. note:: - We use ``python<3.11`` as this will install the most recent 3.10 + We use ``python<3.12`` as this will install the most recent 3.11 release of python. We can now install :mod:`sire` into that environment by typing .. code-block:: bash - $ mamba install -n openbiosim -c openbiosim sire + $ mamba install -n openbiosim -c conda-forge -c openbiosim sire .. note:: The option ``-n openbiosim`` tells ``mamba`` to install :mod:`sire` - into the ``openbiosim`` environment. The option ``-c openbiosim`` + into the ``openbiosim`` environment. The option ``-c conda-forge`` + tells ``mamba`` to use the ``conda-forge`` channel for all + dependencies. The option ``-c openbiosim`` tells ``mamba`` to install :mod:`sire` from the ``openbiosim`` conda channel. @@ -137,7 +139,7 @@ If you want the latest development release, then install by typing .. code-block:: bash - $ mamba install -n openbiosim -c "openbiosim/label/dev" sire + $ mamba install -n openbiosim -c conda-forge -c "openbiosim/label/dev" sire You may (optionally) want to install additional tools such as ``ipython`` and ``jupyterlab``. To do this, type @@ -167,6 +169,14 @@ to learn how to use :mod:`sire` or the 3. Also easy installation - Run in a container ============================================== +.. warning:: + + Because of low demand, pre-built containers are created only for + major releases, and may be out of date compared to the newest release. + Please `get in touch `__ + if you want to use a container and would like us to build the latest + release. + Another route to install :mod:`sire` is to download and run our pre-built containers. These can be run via `docker `__ (on Linux, MacOS and Windows) @@ -233,7 +243,7 @@ branch if you are a developer). You compile :mod:`sire` into an existing anaconda / miniconda environment. Please create and activate an environment, e.g. by following `the instructions <_Install_Mambaforge>` to install a fresh ``mambaforge`` and -then creating and activating Python 3.10 environment called +then creating and activating Python 3.11 environment called ``openbiosim``. Next, download the source code. You could download the latest development @@ -434,8 +444,8 @@ edit the pins in this file too, if you want to do some fine-tuning. You may need to edit the recipe to fix version inconsistencies. This is especially the case for ``rdkit`` - you need to to make - sure that if you specify a version for ``rdkit`` in your - ``environment.yml`` that you also use the same version + sure that if you specify a version for ``rdkit`` in your + ``environment.yml`` that you also use the same version for the ``rdkit-devel`` package. E. Building the package @@ -453,8 +463,8 @@ sire conda package, e.g. .. note:: The above command assumes that you don't need any other channels included - to install all of the packages included in your ``environment.yml``. - The ``actions/update_recipe.py`` script will print out the correct + to install all of the packages included in your ``environment.yml``. + The ``actions/update_recipe.py`` script will print out the correct ``conda mambabuild`` command at the end, which includes any extra channels that are needed. @@ -487,10 +497,10 @@ to conda, etc.). .. note:: A full set of tests will be run on the package after it has been built. - Some of these tests may fail if you have edited the recipe to remove - some of the dependencies. If this happens, you can decide to ignore + Some of these tests may fail if you have edited the recipe to remove + some of the dependencies. If this happens, you can decide to ignore the tests, e.g. by removing them from the conda recipe (``meta.yml``) - or by just copying the file that is produced and has been placed into + or by just copying the file that is produced and has been placed into the ``conda-bld/broken`` directory. You can then install it, either via the channel you've uploaded to, or by diff --git a/doc/source/quickstart/index.rst b/doc/source/quickstart/index.rst index a204ff54f..81e78d2d2 100644 --- a/doc/source/quickstart/index.rst +++ b/doc/source/quickstart/index.rst @@ -235,6 +235,54 @@ And you can even run molecular dynamics using the integration with OpenMM. 18 18 19.0 60.048626 1364.033346 9976.885470 33.771207 9.353054e-08 -59483.125256 144.011481 -64.081234 -838.157968 -48806.614328 19 19 20.0 61.609644 1363.227511 9935.113811 32.689491 2.968531e-07 -59452.910379 144.467384 -64.420695 -835.410744 -48815.633975 +This integration extends to running GPU-accelerated alchemical molecular +dynamics free energy calculations. For example, load up this merged +molecule that uses a λ-coordinate to morph between ethane and methanol. + +>>> mols = sr.load(sr.expand(sr.tutorial_url, "merged_molecule.s3")) + +We'll now select the reference state (ethane)... + +>>> for mol in mols.molecules("molecule property is_perturbable"): +... mols.update(mol.perturbation().link_to_reference().commit()) + +To calculate a free energy, we would need to run multiple simulations +across λ, calculating the difference in energy between neighbouring +λ-windows. How to do this is +:doc:`described in full here <../tutorial/index_part06>`. For this +quick start guide, we'll just run at λ=0.5, calculating the difference +in energy between λ=0.5 and λ=0.0 and λ=1.0. + +First, lets minimise at λ=0.5. + +>>> mols = mols.minimisation(lambda_value=0.5).run().commit() + +And not lets run a short dynamics simulation at λ=0.5, calculating the +energy differences to λ=0.0 and λ=1.0. + +>>> d = mols.dynamics(lambda_value=0.5, timestep="1fs", temperature="25oC") +>>> d.run("5ps", energy_frequency="0.5ps", lambda_windows=[0.0, 1.0]) + +We can now extract the energy differences in an alchemlyb-compatible +pandas DataFrame. + +>>> print(d.commit().energy_trajectory().to_alchemlyb()) + 0.0 0.5 1.0 +time fep-lambda +0.5 0.5 -46062.581449 -46066.727158 -46064.969441 +1.0 0.5 -43335.286994 -43351.114108 -43348.340617 +1.5 0.5 -41676.617061 -41685.901393 -41684.741190 +2.0 0.5 -40392.020981 -40449.196087 -40448.185263 +2.5 0.5 -39720.265484 -39743.531651 -39740.788035 +3.0 0.5 -39067.899327 -39094.033563 -39092.723982 +3.5 0.5 -38659.468399 -38677.805074 -38676.375989 +4.0 0.5 -38354.706210 -38366.022091 -38363.577232 +4.5 0.5 -38160.932310 -38175.773525 -38172.820779 +5.0 0.5 -37844.757596 -37850.187961 -37848.280865 + +See :ref:`the tutorial ` +for a complete script to run and analyse relative free energy calculations. + This is just the beginning of what :mod:`sire` can do! To learn more, please take a look at :doc:`the detailed guides <../cheatsheet/index>` or :doc:`the tutorial <../tutorial/index>`. diff --git a/doc/source/tutorial/index_part06.rst b/doc/source/tutorial/index_part06.rst index 6dc6e0dc1..31d8c0778 100644 --- a/doc/source/tutorial/index_part06.rst +++ b/doc/source/tutorial/index_part06.rst @@ -26,3 +26,5 @@ restraints as you perturb between different molecules. part06/02_alchemical_dynamics part06/03_restraints part06/04_alchemical_restraints + part06/05_free_energy_perturbation + part06/06_faster_fep diff --git a/doc/source/tutorial/part06/01_merge.rst b/doc/source/tutorial/part06/01_merge.rst index 21fcecb3d..9d265901a 100644 --- a/doc/source/tutorial/part06/01_merge.rst +++ b/doc/source/tutorial/part06/01_merge.rst @@ -1,45 +1,225 @@ -================= -Merging Molecules -================= +================ +Merged Molecules +================ -This page is a WORK IN PROGRESS. +The concept of a "merged molecule" is central to the way that free energy +calculations are implemented in sire. A merged molecule is one that +represents both a "reference" state and a "perturbed" state. These are +the two states that the free energy simulation will morph between, and +for which the free energy difference will be calculated. -This page will talk about how we can create merge molecules. +For example, here we have pre-prepared a merged molecule that represents +the perturbation from ethane to methanol. -This uses BioSimSpace, so examples will start with - ->>> import BioSimSpace as BSS >>> import sire as sr +>>> mols = sr.load(sr.expand(sr.tutorial_url, "merged_molecule.s3")) +>>> print(mols) +System( name=BioSimSpace_System num_molecules=4054 num_residues=4054 num_atoms=12167 ) -This page will then detail additional functionality that can be used -to control or edit the merge molecule. - -Examples include shrinking ghost atoms. +.. note:: ->>> sr.morph.shrink_ghost_atoms(mols, length="0.6 A") + The ``.s3`` file format is an internal binary format used by sire to + store any object. You can create ``.s3`` files using the + :func:`sire.stream.save` function, and load them using the + :func:`sire.stream.load` function. These files have the ``.s3`` + suffix when created by sire, and the ``.bss`` suffix when created + by BioSimSpace. The ``.s3`` format is designed to + be portable and backwards compatible, but standard file formats, + e.g. ``.pdb``, ``.mol2``, ``.sdf``, etc. should be preferred for + long-term storage of data. -Or doing extra work like adding restraints, or adding extra controls -or forcefield parameters that represent best practice for complex -morphs. +.. note:: -Also show how to visualise the morph, e.g. + This merged molecule was created using `BioSimSpace `__, + and then saved using :func:`sire.stream.save`. + You should use BioSimSpace if you want to create merged molecules + yourself. ->>> merged_mol.perturbation().view() +This system contains a single merged molecule in a box of water. Merged +molecules are idenfitied by the molecule property ``is_perturbable``, which +will be ``True``. We can extract the merged molecule from this system using -Also show the interface for checking if a molecule is perturbable +>>> mol = mols["molecule property is_perturbable"] +>>> print(mol) +Molecule( Merged_Molecule:6 num_atoms=8 num_residues=1 ) + +A merged molecule is exactly the same as all other molecules in sire. The +difference is that it contains two set of the molecular properties; +one that represents the reference state, and one that represents the +perturbed state. These are identified by the ``0`` and ``1`` suffixes. + +For example, the reference state coordinates are in the ``coordinates0`` +property; + +>>> print(mol.property("coordinates0")) +AtomCoords( size=8 +0: (25.71278789630378, 24.93752746353058, 25.253932968775896) +1: (24.28721210369622, 25.062578848992636, 24.746067031224108) +2: (25.911542134040474, 23.88985958968847, 25.56394874111798) +3: (26.425045597240814, 25.22062162179178, 24.45094936681738) +4: (25.86160072175123, 25.60943815695771, 26.125882372721225) +5: (24.13839927824877, 24.3906681555655, 23.87411762727878) +6: (24.088796970412663, 26.110140410311534, 24.43511653656108) +7: (23.574954402759186, 24.779484690731437, 25.549050633182624) +) + +while the perturbed state coordinates are in the ``coordinates1`` property; + +>>> print(mol.property("coordinates1")) +AtomCoords( size=8 +0: (25.65553270521631, 24.945670198487242, 25.22503902796385) +1: (24.34753270521631, 25.064670198487246, 24.744039027963847) +2: (25.911542134040474, 23.88985958968847, 25.56394874111798) +3: (26.247532705216308, 25.207670198487243, 24.474039027963848) +4: (25.86160072175123, 25.60943815695771, 26.125882372721225) +5: (24.194532705216307, 24.386670198487245, 23.877039027963846) +6: (24.14453270521631, 26.114670198487243, 24.441039027963846) +7: (23.63753270521631, 24.781670198487245, 25.548039027963846) +) + +Similarly the reference state atomic charges are in the ``charge0`` property; + +>>> print(mol.property("charge0")) +SireMol::AtomCharges( size=8 +0: -0.09435 |e| +1: -0.09435 |e| +2: 0.03145 |e| +3: 0.03145 |e| +4: 0.03145 |e| +5: 0.03145 |e| +6: 0.03145 |e| +7: 0.03145 |e| +) + +while the perturbed state atomic charges are in the ``charge1`` property; + +>>> print(mol.property("charge1")) +SireMol::AtomCharges( size=8 +0: -0.5988 |e| +1: 0.1167 |e| +2: 0 |e| +3: 0.396 |e| +4: 0 |e| +5: 0.0287 |e| +6: 0.0287 |e| +7: 0.0287 |e| +) + +The atomic elements are in the ``element0`` and ``element1`` properties; + +>>> print(mol.property("element0")) +SireMol::AtomElements( size=8 +0: Carbon (C, 6) +1: Carbon (C, 6) +2: Hydrogen (H, 1) +3: Hydrogen (H, 1) +4: Hydrogen (H, 1) +5: Hydrogen (H, 1) +6: Hydrogen (H, 1) +7: Hydrogen (H, 1) +) +>>> print(mol.property("element1")) +SireMol::AtomElements( size=8 +0: Oxygen (O, 8) +1: Carbon (C, 6) +2: dummy (Xx, 0) +3: Hydrogen (H, 1) +4: dummy (Xx, 0) +5: Hydrogen (H, 1) +6: Hydrogen (H, 1) +7: Hydrogen (H, 1) +) + +Here we can see that the atoms at indexes 2 and 4 go from being hydrogens +in the reference state (with charges of 0.03145 |e|) to being ghost +(or dummy) atoms in the perturbed state, with charges of zero. + +Viewing merged molecules +------------------------ + +The standard :func:`~sire.mol.SelectorMol.view` function uses the standard +``coordinates``, ``element`` and other properties to view molecules. These +properties don't exist in our merged molecule, as instead we have +``coordinates0``, ``coordinates1``, ``element0``, ``element1``, etc. + +To view the molecule, we need to choose which of the reference or perturbed +states we want to view. We do this by linking the standard properties to +either the reference or perturbed versions, e.g. linking ``coordinates`` +to ``coordinates0`` if we want to view the reference state. + +We could do this manually, but it would be a bit tedious. To save typing, +sire provides a :class:`sire.morph.Perturbation` class that makes it easier +to work with merged molecules. You can access this via the +:meth:`~sire.mol.Molecule.perturbation` method. + +>>> pert = mol.perturbation() +>>> print(pert) +Perturbation( Molecule( Merged_Molecule:6 num_atoms=8 num_residues=1 ) ) + +.. note:: + + Calling the :meth:`~sire.mol.Molecule.perturbation` method on a molecule + that is not a merged molecule will raise an exception. + +The :class:`~sire.morph.Perturbation` class provides the +:meth:`~sire.morph.Perturbation.link_to_reference` and +:meth:`~sire.morph.Perturbation.link_to_perturbed` methods. These can +be used to link all of the standard properties to either the reference +or perturbed values. + +>>> mol = pert.link_to_reference().commit() +>>> mol.view() + +.. image:: images/06_01_01.jpg + :alt: View of the reference state (ethane) + +has viewed the reference state (ethane), while + +>>> mol = pert.link_to_perturbed().commit() +>>> mol["not element Xx"].view() + +.. image:: images/06_01_02.jpg + :alt: View of the perturbed state (methanol) + +has viewed the perturbed state (methanol). + +.. note:: + + The perturbed state includes two ghost (dummy) atoms, which should + normally be invisible. However, the ``view`` function will show all + atoms, including ghosts. To hide the ghost atoms, we have chosen + to view all non-ghost atoms, i.e. all atoms that are not element + ``Xx``. + +Viewing merged molecules in their environment +--------------------------------------------- + +So far we have just viewed the merged molecule in isolation. However, the +molecule exists in a system, in this case, a box of water. We can view the +merged molecule in its environment by updating the system with the result +of linking the molecule to either the reference or perturbed states, +e.g. + +>>> mols = mols.update(pert.link_to_reference().commit()) ->>> merged_mol.is_perturbable() +has updated the system with a copy of the merged molecule where all of +its standard properties are linked to the reference state. While -And to search for perturbable molecules in a collection. +>>> mols = mols.update(pert.link_to_perturbed().commit()) ->>> perturbable_mols = mols.molecules("property is_perturbable") +updates the system with a copy of the merged molecule where all of its +standard properties are linked to the perturbed state. -or +In general, a system could contain many merged molecules. To link all of them +to the reference state you could use ->>> mol = mols["molecule property is_perturbable"] +>>> for mol in mols.molecules("molecule property is_perturbable"): +... mols.update(mol.perturbation().link_to_reference().commit()) -Also show how to use links to link to the 0 or 1 properties +or to link all of them to the perturbed state you could use ->>> mol = mol.edit().add_link("coordinates0", "coordinates").commit() +>>> for mol in mols.molecules("molecule property is_perturbable"): +... mols.update(mol.perturbation().link_to_perturbed().commit()) -(would be good to add this to the Cursor API) +Now you could view and manipulate them as normal, e.g. using +``mols.view()`` etc. diff --git a/doc/source/tutorial/part06/02_alchemical_dynamics.rst b/doc/source/tutorial/part06/02_alchemical_dynamics.rst index b03c63b8d..2368c1c71 100644 --- a/doc/source/tutorial/part06/02_alchemical_dynamics.rst +++ b/doc/source/tutorial/part06/02_alchemical_dynamics.rst @@ -6,9 +6,7 @@ You can create an alchemical molecular dynamics simulation in exactly the same way as you would a normal molecular dynamics simulation. There are two options: -* Use the high-level interface based on the -:func:`~sire.system.System.minimisation` and -:func:`~sire.system.System.dynamics` functions. +* Use the high-level interface based on the :func:`~sire.system.System.minimisation` and :func:`~sire.system.System.dynamics` functions. * Use the low-level interface that works directly with native OpenMM objects. High level interface @@ -17,36 +15,35 @@ High level interface The simplest route is to use the high-level interface. Calling :func:`~sire.system.System.minimisation` or :func:`~sire.system.System.dynamics` on any collection of molecules (or -view) that contains any perturbable molecule will automatically create +view) that contains one or more merged molecules will automatically create an alchemical simulation. There are extra options that you can pass to the simulation that will control how it is run: * ``lambda_value`` - this sets the global λ-value for the simulation. - λ is a parameter that controls the morph from the reference molecule - (at λ=0) to the perturbed molecule (at λ=1). + λ is a parameter that controls the morph from the reference state + (at λ=0) to the perturbed state (at λ=1). * ``swap_end_states`` - if set to ``True``, this will swap the end states - of the perturbation. The morph will run from the perturbed molecule - (at λ=0) to the reference molecule (at λ=1). Note that the coordinates + of the perturbation. The morph will run from the perturbed state + (at λ=0) to the reference state (at λ=1). Note that the coordinates of the perturbed molecule will be used in this case to start the - simulation. + simulation. This can be useful to calculate the reverse of the free + energy potential of mean force (PMF), to check that the reverse + free energy equals to forwards free energy. -* ``schedule`` - set the λ-schedule which specifies how λ morphs between - the reference and perturbed molecules. - -* ``shift_delta`` - set the ``shift_delta`` parameter which is used to - control the electrostatic and van der Waals softening potential that - smooths the creation or deletion of ghost atoms. This is a floating - point number that defaults to ``1.0``, which should be good for - most perturbations. - -* ``coulomb_power`` - set the ``coulomb_power`` parameter which is used - to control the electrostatic softening potential that smooths the - creation and deletion of ghost atoms. This is an integer that defaults - to ``0``, which should be good for most perturbations. +* ``ignore_perturbations`` - if set to ``True``, this will ignore any + perturbations, and will run the dynamics without a λ-coordinate, using + the properties from the reference state (or from the perturbed state + if ``swap_end_states`` is ``True``). This is useful if you want to + run a standard dynamics simulation of the reference or perturbed states, + without the machinery of alchemical dynamics. For example, we could minimise our alchemical system at λ=0.5 using +>>> import sire as sr +>>> mols = sr.load(sr.expand(sr.tutorial_url, "merged_molecule.s3")) +>>> for mol in mols.molecules("molecule property is_perturbable"): +... mols.update(mol.perturbation().link_to_reference().commit()) >>> mols = mols.minimisation(lambda_value=0.5).run().commit() We can then run some dynamics at this λ-value using @@ -54,15 +51,11 @@ We can then run some dynamics at this λ-value using >>> d = mols.dynamics(lambda_value=0.5, timestep="4fs", temperature="25oC") >>> d.run("10ps", save_frequency="1ps") >>> mols = d.commit() +>>> print(d) +Dynamics(completed=10 ps, energy=-33741.4 kcal mol-1, speed=90.7 ns day-1) The result of dynamics is a trajectory run at this λ-value. You can view the -trajectory as you would any other, e.g. - ->>> mols.view() - -You could view specifically the perturbable molecules using - ->>> mols["molecule property is_perturbable"].view() +trajectory as you would any other, using ``mols.view()``. Energy Trajectories ------------------- @@ -73,18 +66,53 @@ sampled by the molecules during the trajectory. You can access this energy trajectory via the :func:`~sire.system.System.energy_trajectory` function, e.g. ->>> df = mols.energy_trajectory() +>>> t = mols.energy_trajectory() +>>> print(t) +EnergyTrajectory( size=10 +time lambda 0.5 kinetic potential +1 0.5 -45509.8 4402.32 -45509.8 +2 0.5 -43644.2 6105.82 -43644.2 +3 0.5 -42466.9 6635.76 -42466.9 +4 0.5 -41678.9 7006.96 -41678.9 +5 0.5 -41306.2 7019.09 -41306.2 +6 0.5 -41201.9 7098.73 -41201.9 +7 0.5 -41046.8 7229.69 -41046.8 +8 0.5 -40949.4 7172.58 -40949.4 +9 0.5 -40928.5 7234.13 -40928.5 +10 0.5 -40855.4 7226.52 -40855.4 +) + +.. note:: + + Dynamics uses a random number generator for the initial velocities + and temperature control. The exact energies you get from dynamics will + be different to what is shown here. + +The result is a :class:`sire.maths.EnergyTrajectory` object. If you want, +you can also ask for the energy trajetory to be returned as a +`pandas DataFrame `__, +by using the ``to_pandas`` argument, e.g. + +>>> df = mols.energy_trajectory(to_pandas=True) >>> print(df) -PANDAS DATAFRAME + lambda 0.5 kinetic potential +time +1.0 0.5 -45509.784041 4402.317703 -45509.784041 +2.0 0.5 -43644.165016 6105.819473 -43644.165016 +3.0 0.5 -42466.942263 6635.760728 -42466.942263 +4.0 0.5 -41678.910475 7006.955030 -41678.910475 +5.0 0.5 -41306.210905 7019.094675 -41306.210905 +6.0 0.5 -41201.944653 7098.727370 -41201.944653 +7.0 0.5 -41046.829930 7229.685671 -41046.829930 +8.0 0.5 -40949.435093 7172.578640 -40949.435093 +9.0 0.5 -40928.462340 7234.130710 -40928.462340 +10.0 0.5 -40855.386336 7226.516672 -40855.386336 -By default, this trajectory is returned as a -`pandas DataFrame `__. -You can get the underlying :class:`sire.maths.EnergyTrajectory` function -by passing ``to_pandas=False``, e.g. +.. note:: ->>> t = mols.energy_trajectory(to_pandas=False) ->>> print(t) -EnergyTrajectory() + You could also have obtained a DataFrame by calling the + :meth:`~sire.maths.EnergyTrajectory.to_pandas` function of + the :class:`~sire.maths.EnergyTrajectory` object. You calculate free energies by evaluating the potential energy for different values of λ during dynamics. You can control which values of λ are used @@ -93,12 +121,41 @@ values of λ during dynamics. You can control which values of λ are used >>> d = mols.dynamics(lambda_value=0.5, timestep="4fs", temperature="25oC") >>> d.run("10ps", save_frequency="1ps", lambda_windows=[0.4, 0.6]) >>> mols = d.commit() ->>> print(mols.get_energy_trajectory()) +>>> print(mols.energy_trajectory().to_pandas()) + lambda 0.4 0.5 0.6 kinetic potential +time +1.0 0.5 NaN -45509.784041 NaN 4402.317703 -45509.784041 +2.0 0.5 NaN -43644.165016 NaN 6105.819473 -43644.165016 +3.0 0.5 NaN -42466.942263 NaN 6635.760728 -42466.942263 +4.0 0.5 NaN -41678.910475 NaN 7006.955030 -41678.910475 +5.0 0.5 NaN -41306.210905 NaN 7019.094675 -41306.210905 +6.0 0.5 NaN -41201.944653 NaN 7098.727370 -41201.944653 +7.0 0.5 NaN -41046.829930 NaN 7229.685671 -41046.829930 +8.0 0.5 NaN -40949.435093 NaN 7172.578640 -40949.435093 +9.0 0.5 NaN -40928.462340 NaN 7234.130710 -40928.462340 +10.0 0.5 NaN -40855.386336 NaN 7226.516672 -40855.386336 +11.0 0.5 -43210.778752 -43211.624385 -43211.937516 6219.582732 -43211.624385 +12.0 0.5 -42405.359296 -42405.219032 -42404.934648 6807.796036 -42405.219032 +13.0 0.5 -41724.163073 -41724.470944 -41724.664572 6942.723457 -41724.470944 +14.0 0.5 -41282.749354 -41282.638965 -41282.324706 7174.870863 -41282.638965 +15.0 0.5 -41091.395386 -41090.926489 -41090.283597 7121.938171 -41090.926489 +16.0 0.5 -41020.530186 -41020.300294 -41019.896407 7243.171748 -41020.300294 +17.0 0.5 -41027.939363 -41027.739347 -41027.365337 7112.568694 -41027.739347 +18.0 0.5 -40947.962069 -40948.210188 -40948.254438 7086.990231 -40948.210188 +19.0 0.5 -41063.162834 -41064.038343 -41064.470977 7173.505327 -41064.038343 +20.0 0.5 -40997.466132 -40997.774003 -40997.907880 7211.324644 -40997.774003 This has run a new trajectory, evaluating the potential energy at the simulation λ-value (0.5) as well as at λ-windows 0.4 and 0.6. You can pass in as many or as few λ-windows as you want. +.. note:: + + Notice how the potential energy is evaluated at λ=0.4, λ=0.5 and λ=0.6 + only from 11ps onwards. The first 10ps was the first block of dynamics + where we only evaluted the energy at the simulated λ-value. The second + block of 10ps also evaluated the energy at λ=0.4 and λ=0.6. + Controlling the trajectory frequency ------------------------------------ @@ -116,40 +173,269 @@ but saving energies every 20 femtoseconds. >>> d = mols.dynamics(lambda_value=0.5, timestep="4fs", temperature="25oC") >>> d.run("10ps", frame_frequency="1ps", energy_frequency="20fs", -... lambda_windows=[0.4, 0.6], save_velocities=False) +... lambda_windows=[0.4, 0.6]) >>> mols = d.commit() ->>> print(mols.get_energy_trajectory()) +>>> print(mols.energy_trajectory().to_pandas()) + lambda 0.4 0.5 0.6 kinetic potential +time +1.00 0.5 NaN -45509.784041 NaN 4402.317703 -45509.784041 +2.00 0.5 NaN -43644.165016 NaN 6105.819473 -43644.165016 +3.00 0.5 NaN -42466.942263 NaN 6635.760728 -42466.942263 +4.00 0.5 NaN -41678.910475 NaN 7006.955030 -41678.910475 +5.00 0.5 NaN -41306.210905 NaN 7019.094675 -41306.210905 +... ... ... ... ... ... ... +29.92 0.5 -40892.930998 -40892.850485 -40892.506350 7306.237779 -40892.850485 +29.94 0.5 -40890.720195 -40891.326823 -40891.729581 7314.180669 -40891.326823 +29.96 0.5 -40869.209679 -40868.920036 -40868.426522 7305.782760 -40868.920036 +29.98 0.5 -40828.250071 -40828.378688 -40828.243683 7270.175960 -40828.378688 +30.00 0.5 -40764.913551 -40764.504405 -40763.921264 7140.896560 -40764.504405 + +Controlling perturbations with a λ-schedule +------------------------------------------- + +So far the perturbation from the reference to the perturbed state has been +linear. λ has acted on each of the perturbable properties of the molecule +by scaling the ``initial`` value from the reference state to the ``final`` +value in the perturbed state using the equation + +.. math:: + + (1 - \lambda) \times \mathrm{initial} + \lambda \times \mathrm{final} + +This shows that at λ=0, the perturbable properties are set to the +``initial`` value, and at λ=1, the perturbable properties are set to the +``final`` value. At intermediate values of λ, the perturbable properties +are linearly interpolated between the ``initial`` and ``final`` values, +e.g. at λ=0.5, the perturbable properties are set to half-way between the +``initial`` and ``final`` values. + +The perturbation of the parameters is controlled in the code using +a :class:`sire.cas.LambdaSchedule`. + +You can get the λ-schedule used by the dynamics simulation using the +:func:`~sire.mol.Dynamics.lambda_schedule` function, e.g. + +>>> s = d.lambda_schedule() +>>> print(s) +LambdaSchedule( + morph: λ * final + initial * (-λ + 1) +) + +You can plot how this schedule would morph the perturbable properties +using the :func:`~sire.cas.LambdaSchedule.plot` function, e.g. + +>>> df = get_lever_values(initial=2.0, final=3.0) +>>> df.plot() + +.. image:: images/06_02_01.jpg + :alt: View of the default λ-schedule + +This shows how the different levers available in this schedule would morph +a hyperthetical parameter that has an ``initial`` value of ``2.0`` and a +``final`` value of ``3.0``. + +In this case the levers are all identical, so would change the parameter +in the same way. You can choose your own equation for the λ-schedule. +For example, maybe we want to scale the charge by the square of λ. + +>>> s.set_equation("morph", "charge", +... s.lam()**2 * s.final() + s.initial() * (1 - s.lam()**2)) +>>> print(s) +LambdaSchedule( + morph: λ * final + initial * (-λ + 1) + charge: λ^2 * final + initial * (-λ^2 + 1) +) + +This shows that in the default ``morph`` stage of this schedule, we will +scale the ``charge`` parameters by λ^2, while all other parameters (the +default) will scale using λ. We can plot the result of this using; + +>>> s.get_lever_values(initial=2.0, final=3.0).plot() + +.. image:: images/06_02_02.jpg + :alt: View of the λ-schedule where charge is scaled by λ^2 + +The above affected the default ``morph`` stage of the schedule. You can +append or prepend additional stages to the schedule using the +:meth:`~sire.cas.LambdaSchedule.append_stage` or +:meth:`~sire.cas.LambdaSchedule.prepend_stage` functions, e.g. + +>>> s.append_stage("scale", s.final()) +>>> print(s) + +would append a second stage, called ``scale``, which by default would +use the ``final`` value of the parameter. We could then add a lever to +this stage that scales down the charge to 0, + +>>> s.set_equation("scale", "charge", (1-s.lam()) * s.final()) +>>> print(s) +LambdaSchedule( + morph: λ * final + initial * (-λ + 1) + charge: λ^2 * final + initial * (-λ^2 + 1) + scale: final + charge: final * (-λ + 1) +) + +Again, it is worth plotting the impact of this schedule on a hyperthetical +parameter. + +>>> s.get_lever_values(initial=2.0, final=3.0).plot() + +.. image:: images/06_02_03.jpg + :alt: View of the λ-schedule where charge is scaled in a second stage to zero. + +.. note:: + + The value of the stage-local λ in each stage moves from 0 to 1. + This is different to the global λ, which moves from 0 to 1 evenly + over all stages. For example, in this 2-stage schedule, values of + global λ between 0 and 0.5 will be in the ``morph`` stage, while + values of global λ between 0.5 and 1 will be in the ``scale`` stage. + Within the ``morph`` stage, the local λ will move from 0 to 1 + (corresponding to global λ values of 0 to 0.5), while + within the ``scale`` stage, the local λ will move from 0 to 1 + (corresponding to gloabl λ values of 0.5 to 1). + +Through the combination of adding stages and specifyig different equations +for levers, you can have a lot of control over how the perturbable properties +are morphed from the reference to the perturbed states. + +To make things easier, there are some simple functions that let you add +some common stages. + +>>> s = sr.cas.LambdaSchedule() +>>> print(s) +LambdaSchedule::null +>>> s.add_morph_stage("morph") +>>> print(s) +LambdaSchedule( + morph: final * λ + initial * (-λ + 1) +) +>>> s.add_charge_scale_stages("decharge", "recharge", scale=0.2) +>>> print(s) +LambdaSchedule( + decharge: initial + charge: initial * (-λ * (-γ + 1) + 1) + morph: final * λ + initial * (-λ + 1) + charge: γ * (final * λ + initial * (-λ + 1)) + recharge: final + charge: final * (-(-γ + 1) * (-λ + 1) + 1) + γ == 0.2 +) + +has created a null schedule, and then added a ``morph`` stage using the +default perturbation equation. This is then sandwiched by two stages; +a ``decharge`` stage that scales the charge lever from the ``initial`` +value to γ times that value, and a ``recharge`` stage that scales +the charge lever from γ times the ``final`` value to the full +``final`` value. It also scales the charge lever in the ``morph`` stage +by γ, which is set to 0.2 for all stages. + +We can see how this would affect a hyperthetical parameter that goes +from an ``initial`` value of 2.0 to a ``final`` value of 3.0 via + +>>> s.get_lever_values(initial=2.0, final=3.0).plot() + +.. image:: images/06_02_04.jpg + :alt: View of the λ-schedule that sandwiches a standard morph stage + between two stages that scale the charge lever. .. note:: - The ``save_velocities`` option turns on or off the saving of atomic - velocities to the frame trajectory. Typically you don't need to - record velocities, so it is safe to switch them off. This can - reduce memory consumption and also slightly speed up your simulation. + Schedules constructed outside of the dynamics simulation do not have + the full set of levers (e.g. torsion_k, dih_scale etc) as + levers are only added as they are needed (hence why only + ``default`` and ``charge`` are shown here). The additional levers + are added when the schedule is added to the simulation. + +Once you have created your schedule you can add it via the +:meth:`~sire.mol.Dynamics.set_schedule` function of the +:class:`~sire.mol.Dynamics` object, e.g. + +>>> d.set_schedule(s) + +Alternatively, you can set the schedule when you call the +:meth:`~sire.mol.SelectorMol.dynamics` function, e.g. + +>>> d = mols.dynamics(lambda_value=0.5, timestep="4fs", temperature="25oC", +... schedule=s) +>>> print(d.get_schedule()) +LambdaSchedule( + decharge: initial + charge: (-(-γ + 1) * λ + 1) * initial + morph: final * λ + (-λ + 1) * initial + charge: γ * (final * λ + (-λ + 1) * initial) + recharge: final + charge: (-(-γ + 1) * (-λ + 1) + 1) * final + γ == 0.2 +) + +Ghost Atoms and Softening potentials +------------------------------------ + +Internally the alchemical dynamics simulation works by calculating morphed +forcefield parameters whenever λ is changed, and then calling the +OpenMM `updateParametersInContext function `__ +function to update those parameters in all of the OpenMM Force objects that +are used to calculate atomic forces. This is a very efficient way of +performing a perturbation, as it allows vanilla (standard) OpenMM forces +to be used for the dynamics. However, there are challenges with how +we handle atoms which are created or deleted during the perturbation. + +These atoms, which we call "ghost atoms", provide either a space into which +a new atom is grown, or a space from which an atom is deleted. To ensure +these ghost atoms don't cause dynamics instabilities, we use a softening +potential to model their interactions with all other atoms. These +softening potentials (also called "soft-core" potentials) soften the +charge and Lennard-Jones interactions between the ghost atoms and all +other atoms using an α (alpha) parameter. This is a perturbable parameter +of the atoms, which is equal to 1 when the atom is in a ghost state, +and 0 when it is not. For example, an atom which exists in the reference +state but becomes a ghost in the perturbed state would have an α value +that would go from 0 to 1. Alternatively, an atom which does not exist +in the reference state, and that appears in the perturbed state, would +have an α value that would go from 1 to 0. -Setting up a λ-schedule ------------------------ +.. note:: -A λ-schedule (represented using the :class:`sire.cas.LambdaSchedule` class) -specifies how the λ-parameter morphs from the reference to the perturbed -molecules. The λ-schedule achieves this... + You can have as many or few ghost atoms as you want in your merged + molecule. If all atoms become ghosts, then this is the same as + completely decoupling the molecule, as you would do in an + absolute binding free energy calculation. Equally, you could + run a "dual topology" calculation by two sets of atoms in your + merged molecule - the first set start as the reference state atoms, + and all become ghosts, while the second set start all as ghosts + and become the perturbed state atoms. In "single topology" calculations + you would only use ghost atoms for those which don't exist in either + of the end states. -WRITE MORE ABOUT THE λ-schedule +There are two parameters that control the softening potential: + +* ``shift_delta`` - set the ``shift_delta`` parameter which is used to + control the electrostatic and van der Waals softening potential that + smooths the creation or deletion of ghost atoms. This is a floating + point number that defaults to ``1.0``, which should be good for + most perturbations. + +* ``coulomb_power`` - set the ``coulomb_power`` parameter which is used + to control the electrostatic softening potential that smooths the + creation and deletion of ghost atoms. This is an integer that defaults + to ``0``, which should be good for most perturbations. Low level interface ------------------- The high-level interface is just a set of convienient wrappers around the OpenMM objects which are used to run the simulation. If you convert -any set of views (or view) that contains perturbable molecules, then an +any set of views (or view) that contains merged molecules, then an alchemical OpenMM context will be returned. >>> context = sr.convert.to(mols, "openmm") >>> print(context) -OUTPUT +openmm::Context( num_atoms=12167 integrator=VerletIntegrator timestep=1.0 fs platform=HIP ) The context is held in a low-level class, -:class:`~sire.Convert.SireOpenMM.SOMMContext`, inherits from the +:class:`~sire.Convert.SireOpenMM.SOMMContext`, which inherits from the standard `OpenMM Context `__ class. @@ -174,3 +460,40 @@ are; * :func:`~sire.Convert.SireOpenMM.SOMMContext.get_energy` - return the current potential energy of the context. This will be in :mod:`sire` units if ``to_sire_units`` is ``True`` (the default). + +Note that you can also set the ``lambda_value`` and ``lambda_schedule`` +when you create the context using the ``map``, e.g. + +>>> context = sr.convert.to(mols, "openmm", +... map={"lambda_value": 0.5, "schedule": s}) +>>> print(context) +openmm::Context( num_atoms=12167 integrator=VerletIntegrator timestep=1.0 fs platform=HIP ) +>>> print(context.get_lambda()) +0.5 +>>> print(context.get_lambda_schedule()) +LambdaSchedule( + decharge: initial + charge: (-(-γ + 1) * λ + 1) * initial + morph: final * λ + (-λ + 1) * initial + charge: γ * (final * λ + (-λ + 1) * initial) + recharge: final + charge: (-(-γ + 1) * (-λ + 1) + 1) * final + γ == 0.2 +) + +You can then run dynamics as you would do normally using the standard +OpenMM python API, e.g. + +>>> integrator = context.getIntegrator() +>>> integrator.step(100) + +You can then call ``get_potential_energy()`` and ``set_lambda()`` to +get the energy during dynamics for different values of λ, e.g. + +>>> context.set_lambda(0.0) +>>> print(context.get_lambda(), context.get_potential_energy()) +0.0 -38727.2 kcal mol-1 +>>> context.set_lambda(0.5) +>>> print(context.get_lambda(), context.get_potential_energy()) +0.5 -38743.8 kcal mol-1 + diff --git a/doc/source/tutorial/part06/03_restraints.rst b/doc/source/tutorial/part06/03_restraints.rst index 084e83401..e715aca9e 100644 --- a/doc/source/tutorial/part06/03_restraints.rst +++ b/doc/source/tutorial/part06/03_restraints.rst @@ -311,18 +311,7 @@ of the atoms. Here, >>> d.run("10ps") >>> mols = d.commit() -we have run dynamics where all of the carbon atoms are fixed. - -.. note:: - - Be careful using constraints with fixed atoms. If a constraint involves - a fixed atom and a mobile atom, then there is a high risk that the - constraint won't be able to be satisfied during dynamics, and - ``Particle coordinate is NaN`` error or similar will occur. - It it safest to use a small timestep (e.g. 1 fs) when constraining - only parts of molecules. - -while +we have run dynamics where all of the carbon atoms are fixed, while >>> d = mols.dynamics(timestep="1fs", temperature="25oC", ... fixed=[0, 2, 4, 6, 8]) @@ -345,3 +334,24 @@ molecules within the sphere of mobile atoms, e.g. ... fixed=f"not molecules within {radius} of {center}") >>> d.run("10ps") >>> mols = d.commit() + +This implements the restraints in OpenMM by setting the masses of the fixed +atoms to zero. This signals the OpenMM integrator to skip over these +atoms and not move them. + +.. note:: + + Be careful using constraints with fixed atoms. If a constraint involves + a fixed atom and a mobile atom, then there is a high risk that the + constraint won't be able to be satisfied during dynamics, and + ``Particle coordinate is NaN`` error or similar will occur. + It it safest to use a small timestep (e.g. 1 fs) when constraining + only parts of molecules. + +.. note:: + + While the fixed atoms are not moved by the integrator, they are still + included in the energy calculation. This means that the computational + cost of the simulation will still be high, and also that energies + will not be conserved. Make sure you use an integrator that provides + a heat bath (e.g. NVT or NPT integrators). diff --git a/doc/source/tutorial/part06/05_free_energy_perturbation.rst b/doc/source/tutorial/part06/05_free_energy_perturbation.rst new file mode 100644 index 000000000..92514ea53 --- /dev/null +++ b/doc/source/tutorial/part06/05_free_energy_perturbation.rst @@ -0,0 +1,316 @@ +============================================= +Running a Free Energy Perturbation Simulation +============================================= + +Now that you have your (possibly restrained) alchemical system, you can +run a free energy perturbation simulation. There are many ways to do this, +and the exact protocol will depend on what type of free energy you are +calculating, what level of accuracy you require and what compute resource +you have available. :mod:`sire` provides the building blocks to run a +simulation in any way you like, but in this tutorial we will use the +a very simple protocol. Take a look at +`BioSimSpace `__ or the upcoming +`somd2 software `__ if you are +interested in higher-level interfaces to this functionality that +automatically run more complex protocols. + +At its simplest, a free energy perturbation simulation consists of a series +of alchemical molecular dynamics simulations run at different values +of the alchemical parameter, λ. Let's start with the ethane to methanol +merged molecule from the previous section. + +>>> import sire as sr +>>> mols = sr.load(sr.expand(sr.tutorial_url, "merged_molecule.s3")) +>>> print(mols) +System( name=BioSimSpace_System num_molecules=4054 num_residues=4054 num_atoms=12167 ) + +And lets link the properties to the reference state, so that we start the +simulation using the reference state coordinates. + +>>> for mol in mols.molecules("molecule property is_perturbable"): +... mols.update(mol.perturbation().link_to_reference().commit()) + +Next we will run through 21 evenly-space λ values from 0 to 1. We've picked +21 because it is a reasonable number that should over-sample the λ-coordinate. + +For each λ-value, we will minimise the system, generate random velocities +from a Boltzmann distribution for the simulation temperature, +equilibrate the molecules for 2 ps, then +run dynamics for 25 ps. We will calculate the energy of each simulation at +each λ-value, plus at the neighbouring λ-values. This will allow us to +calculate energy differences which we will later use to calculate the +free energy difference. + +We will calculate these energies every 0.1 ps, and will write them at the +end of each block of dynamics via the :class:`~sire.maths.EnergyTrajectory` +output which we will stream as a :mod:`sire` object. + +This will let us subsequently calculate the free energy across λ using +`alchemlyb `__. + +>>> for l in range(0, 105, 5): +... # turn l into the lambda value by dividing by 100 +... lambda_value = l / 100.0 +... print(f"Simulating lambda={lambda_value:.2f}") +... # minimise the system at this lambda value +... min_mols = mols.minimisation(lambda_value=lambda_value).run().commit() +... # create a dynamics object for the system +... d = min_mols.dynamics(timestep="1fs", temperature="25oC", +... lambda_value=lambda_value) +... # generate random velocities +... d.randomise_velocities() +... # equilibrate, not saving anything +... d.run("2ps", save_frequency=0) +... print("Equilibration complete") +... print(d) +... # get the values of lambda for neighbouring windows +... lambda_windows = [lambda_value] +... if lambda_value > 0: +... lambda_windows.insert(0, (l-5)/100.0) +... if lambda_value < 1: +... lambda_windows.append((l+5)/100.0) +... # run the dynamics, saving the energy every 0.1 ps +... d.run("25ps", energy_frequency="0.1ps", frame_frequency=0, +... lambda_windows=lambda_windows) +... print("Dynamics complete") +... print(d) +... # stream the EnergyTrajectory to a sire save stream object +... sr.stream.save(d.commit().energy_trajectory(), +... f"energy_{lambda_value:.2f}.s3") + +This is a very simple protocol for a free energy simulation, with simulation +times chosen so that it will run relatively quickly. For example, with a +modern GPU this should run at about 60 nanoseconds of sampling per day, +so take only 45 seconds or so for each of the 21 λ-windows. + +.. note:: + + This is not enough sampling to calculate a properly converged free energy + for most molecular systems. You would need to experiment with different + equilibration and simulation times, or use a more sophisticated algorithm + implemented via `BioSimSpace `__ or + `somd2 `__. + +.. note:: + + We set the forwards and backwards λ-values to ``(l-5)/100.0`` and + ``(l+5)/100.0`` so that they exactly match the way we set the reference + λ-value. This is because alchemlyb matches energies using a + floating point representation of λ, so small numerical differences + between, e.g. ``(l+5)/100.0`` and ``lambda_value + 0.05`` will lead + to errors (as in this case, ``(l+5)/100.0`` is exactly every 0.05, while + ``lambda_value + 0.05`` has numerical error, meaning some values are + at, e.g. ``0.15000000000000002``). + +At the end of the simulation, you should have 21 energy files, one for each +λ-window. These are called ``energy_0.00.s3``, ``energy_0.05.s3``, ..., +``energy_1.00.s3``. + +Processing Free Energy Data using alchemlyb +-------------------------------------------- + +`alchemlyb `__ is a simple +alchemistry python package for easily analysing the results of free energy +simulations. They have excellent documentation on their website that you +can use if you want to go into the details of how to calculate free +energies. + +Here, we will show a simple BAR analysis of the data that we have just +generated. We can do this because we have calculated data which +alchemlyb can convert into reduced potentials for each λ-window. + +First, we need to import alchemlyb + +>>> import alchemlyb + +.. note:: + + If you see an error then you may need to install (or reinstall) + alchemlyb. You can do this using conda or mamba, e.g. + ``mamba install -c conda-forge alchemlyb``. + +Next, we will load all of the :class:`~sire.maths.EnergyTrajectory` objects +for each λ-window, and will convert them into pandas DataFrames arranged +into an alchemlyb-compatible format. We could do this manually by first +loading all of the s3 files containing the :class:`~sire.maths.EnergyTrajectory` +objects... + +>>> import sire as sr +>>> from glob import glob +>>> dfs = [] +>>> energy_files = glob("energy*.s3") +>>> energy_files.sort() +>>> for energy_file in energy_files: +... dfs.append(sr.stream.load(energy_file).to_alchemlyb()) + +.. note:: + + If you wanted, you could have put the dataframes generated above + directly into the ``dfs`` list here, and not saved them to disk + via the ``.s3`` files. However, this would risk you having to re-run + all of the simulation if you wanted to change the analysis below. + +.. note:: + + Be careful to load the DataFrames in λ-order. The ``glob`` function + can return the files in a random order, hence why we need to sort + this list. This sort only works because we have used a naming convention + for the files that puts them in λ-order. They must be in the right + order or else alchemlyb will calculate the free energy incorrectly + (it uses the column-order rather than the λ-order when calculating + free energies). + +...then joining them together all of these DataFrames into a single +DataFrame... + +>>> import pandas as pd +>>> df = pd.concat(dfs) +>>> print(df) + 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 ... 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 +time fep-lambda ... +2.1 0.0 -39300.750665 -39301.636585 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN +2.2 0.0 -39231.737759 -39232.504175 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN +2.3 0.0 -38982.663906 -38982.683430 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN +2.4 0.0 -38950.308505 -38950.357904 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN +2.5 0.0 -38732.275522 -38733.221193 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN +... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... +26.6 1.0 NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN -37208.485229 -37208.880770 +26.7 1.0 NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN -37332.051194 -37332.476611 +26.8 1.0 NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN -37276.392733 -37276.609020 +26.9 1.0 NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN -37234.238097 -37234.603762 +27.0 1.0 NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN -37268.684798 -37269.319345 + +.. note:: + + Do not worry about the large number of ``NaN`` values. These just show that + we have only calculated free energy differences along the diagonal of this + DataFrame, i.e. only between the simulated and neighbouring λ-windows. + +...or we can use the in-build :func:`sire.morph.to_alchemlyb` function to +do all of the above for us. + +>>> df = sr.morph.to_alchemlyb("energy*.s3") + +This function is not only quicker, but it will automatically sort the +data by λ-value, meaning that you don't need to worry about the naming +convention of your files. + +Now we can tell alchemlyb to calculate the free energy using the BAR method. + +>>> from alchemlyb.estimators import BAR +>>> b = BAR() +>>> b.fit(df) +>>> print(b.delta_f_.loc[0.00, 1.00]) +-3.009332126465953 + +You can get a convergence plot, showing how the free energy changes as +a function of the simulation length using the ``convergence_plot`` function. + +>>> from alchemlyb.convergence import forward_backward_convergence +>>> from alchemlyb.visualisation import plot_convergence +>>> f = forward_backward_convergence(dfs, "bar") +>>> plot_convergence(f) + +.. image:: images/06_05_01.jpg + :alt: Convergence of the free energy estimate as a function of the fraction of simulation length + +All of this shows that the relative free energy for the perturbation of +ethane to methanol in water is about -3.0 kcal mol-1. + +To get the relative hydration free energy, we would need to complete the +cycle by calculating the relative free energy for the perturbation in the +gas phase. We could do this using this code (which is almost identical to +above, except we only simulate the perturbable molecule, and save +the :class:`~sire.maths.EnergyTrajectory` objects to ``energy_gas_{lambda}.s3`` +instead of ``energy_{lambda}.s3``). + +>>> import sire as sr +>>> mols = sr.load(sr.expand(sr.tutorial_url, "merged_molecule.s3")) +>>> mol = mols.molecule("molecule property is_perturbable") +>>> mol.update(mol.perturbation().link_to_reference().commit()) +>>> for l in range(0, 105, 5): +... # turn l into the lambda value by dividing by 100 +... lambda_value = l / 100.0 +... print(f"Simulating lambda={lambda_value:.2f}") +... # minimise the system at this lambda value +... min_mol = mol.minimisation(lambda_value=lambda_value, +... vacuum=True).run().commit() +... # create a dynamics object for the system +... d = min_mol.dynamics(timestep="1fs", temperature="25oC", +... lambda_value=lambda_value, +... vacuum=True) +... # generate random velocities +... d.randomise_velocities() +... # equilibrate, not saving anything +... d.run("2ps", save_frequency=0) +... print("Equilibration complete") +... print(d) +... # get the values of lambda for neighbouring windows +... lambda_windows = [lambda_value] +... if lambda_value > 0: +... lambda_windows.insert(0, (l-5)/100.0) +... if lambda_value < 1: +... lambda_windows.append((l+5)/100.0) +... # run the dynamics, saving the energy every 0.1 ps +... d.run("25ps", energy_frequency="0.1ps", frame_frequency=0, +... lambda_windows=lambda_windows) +... print("Dynamics complete") +... print(d) +... # stream the EnergyTrajectory to a sire save stream object +... sr.stream.save(d.commit().energy_trajectory(), +... f"energy_gas_{lambda_value:.2f}.s3") + +.. note:: + + The option ``vacuum=True`` tells the minimisation and dynamics to + remove any simulation space that may be attached to the molecule(s), + and instead set the space to a :class:`~sire.space.Cartesian` space. + This has the affect of simulating the molecules in vacuum. + +This should run more quickly than the simulation in water, e.g. about +15 seconds per window (at about 150 nanoseconds per day of sampling). + +We can then analyse the results using the same analysis code, except we +switch to analysing the ``energy_gas_{lambda}.s3`` files instead. + +>>> df = sr.morph.to_alchemlyb("energy_gas*.s3") +>>> print(df) + 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 ... 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 +time fep-lambda ... +2.1 0.0 6.966376 7.115377 NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN +2.2 0.0 6.614100 6.434525 NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN +2.3 0.0 6.782235 6.831925 NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN +2.4 0.0 4.196341 4.227257 NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN +2.5 0.0 4.141862 4.205784 NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN +... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... +26.6 1.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN 8.246601 8.707174 +26.7 1.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN 7.169596 7.100342 +26.8 1.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN 6.422485 6.466972 +26.9 1.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN 5.892950 6.129184 +27.0 1.0 NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN 7.487893 7.618827 +>>> from alchemlyb.estimators import BAR +>>> b = BAR() +>>> b.fit(df) +>>> print(b.delta_f_.loc[0.00, 1.00]) +3.186998316374265 + +This shows that the relative free energy for the perturbation of ethane +to methanol in the gas phase is about 3.2 kcal mol-1. Subtracting this +from the free energy in water gives a relative hydration free energy of +about -6.2 kcal mol-1, which is in reasonable agreement with +`published results from other codes `__ +which are in the range of -5.99 kcal mol-1 to -6.26 kcal mol-1. + +.. note:: + + The quoted published results are the difference in computed + absolute hydration free energies calculated using difference codes + and protocols for ethane and methanol, as reported in table 2 + of the linked paper. + +.. note:: + + There will be some variation between different codes and different + protocols, as the convergence of the free energy estimate is sensitive + to the length of the dynamics simulation at each λ-value. In this case, + we used very short simulations. diff --git a/doc/source/tutorial/part06/06_faster_fep.rst b/doc/source/tutorial/part06/06_faster_fep.rst new file mode 100644 index 000000000..fa4b5495a --- /dev/null +++ b/doc/source/tutorial/part06/06_faster_fep.rst @@ -0,0 +1,477 @@ +====================================== +Running Faster Free Energy Simulations +====================================== + +The previous section showed how to calculate the relative hydration free +energy of ethane and methanol using alchemical dynamics. Short dynamics +simulation were run for each λ-value, with the energy differences +between neighbouring λ-values used to calculate the free energy differences. + +The simulations used a timestep of 1 fs, and calculated energy differences +every 0.1 ps (100 steps). This calculated a free energy that should be +accurate, but the simulation took a long time to run. In this section, we +will show how to calculate the same result, but with a much faster +simulation. + +Timesteps and Constraints +------------------------- + +The easiest way to speed up a simulation is to increase the dynamic +timestep. This is the amount of time between each step of the simulation, +i.e. the amount of time between each calculation of the forces. The longer +the timestep, the faster the simulation will run, but at the cost of +more unstable dynamics. At an extreme, the simulation will become +totally unstable and an +`OpenMM Particle Exception `__ +will be raised. + +As a rule of thumb, the timestep should be less than the fastest vibrational +motion in the simulation. Since the fastest vibrations will likely involve +bonds with hydrogen atoms, we can make the simulation more stable by +contraining the bonds that involve hydrogen atoms. + +We can choose the constraint to use via the ``constraint`` keyword, e.g. +after loading the molecules and minimising, + +>>> import sire as sr +>>> mols = sr.load(sr.expand(sr.tutorial_url, "merged_molecule.s3")) +>>> for mol in mols.molecules("molecule property is_perturbable"): +... mols.update(mol.perturbation().link_to_reference().commit()) +>>> mols = mols.minimisation().run().commit() + +...we can turn on constraints of the bonds involving hydrogen atoms by +setting ``constraint`` to ``h-bonds``. + +>>> d = mols.dynamics(timestep="2fs", temperature="25oC", constraint="h-bonds") +>>> d.run("5ps") +>>> print(d) +Dynamics(completed=5 ps, energy=-32140.6 kcal mol-1, speed=51.2 ns day-1) + +We would hope that this is faster than running a simulation with no constraints +and a 1 fs timestep... + +>>> d = mols.dynamics(timestep="1fs", temperature="25oC", constraint="none") +>>> d.run("5ps") +>>> print(d) +Dynamics(completed=5 ps, energy=-27197.8 kcal mol-1, speed=67.3 ns day-1) + +...but this is not the case. We can see here that the cost of the constraints +has out-weighed the benefit of having half the simulation steps. + +The simulation can go a lot faster with a 4 fs timestep... + +>>> d = mols.dynamics(timestep="4fs", temperature="25oC", constraint="h-bonds") +>>> d.run("4ps") +>>> print(d) +Dynamics(completed=4 ps, energy=-32807 kcal mol-1, speed=100.9 ns day-1) + +However, turning on ``h-bonds`` constraints would not be enough to keep the +simulation stable for larger timesteps. For example, if we use a timestep +of 5 fs... + +>>> d = mols.dynamics(timestep="5fs", temperature="25oC", constraint="h-bonds") +>>> d.run("5ps") +OpenMMException: Particle coordinate is NaN. For more information, see +https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#nan + +For such timesteps, we would need to constrain all bonds and angles +involving hydrogen, using the ``h-bonds-h-angles`` constraint. + +>>> d = mols.dynamics(timestep="5fs", temperature="25oC", +... constraint="h-bonds-h-angles") +>>> d.run("5ps") +>>> print(d) +Dynamics(completed=5 ps, energy=-34450.6 kcal mol-1, speed=116.2 ns day-1) + +You can go even further by constraining all bonds, and all angles involving +hydrogen using the ``bonds-h-angles`` constraint... + +>>> d = mols.dynamics(timestep="8fs", temperature="25oC", +... constraint="bonds-h-angles") +>>> d.run("5ps") +>>> print(d) +Dynamics(completed=5 ps, energy=-34130.9 kcal mol-1, speed=152.6 ns day-1) + +.. note:: + + There is also the ``bonds`` constraint which constrains + all chemical bonds, but this doesn't really help unless + angles involving hydrogen aren't also constrained. + +.. note:: + + Different molecular systems will behave differently, and may produce + more, or less stable dynamics than that shown above. As a rule of thumb, + start using a 4 fs timestep with ``h-bonds-h-angles`` constraints, and + then either increase the timestep and dial up the constraints if you can, + or reduce the timestep and reduce the timestep if the simulation becomes + unstable. + +To make things easy, :mod:`sire` automatically chooses a suitable constraint +based on the timestep. You can see the constraint chosen using the +:meth:`~sire.mol.Dynamics.constraint` function. + +>>> d = mols.dynamics(timestep="8fs", temperature="25oC") +>>> print(d.constraint()) +bonds-h-angles + +You can disable all constraints by setting ``constraint`` equal to ``none``, +e.g. + +>>> d = mols.dynamics(timestep="4fs", temperature="25oC", constraint="none") +>>> print(d.constraint()) +none +>>> d.run("5ps") +RuntimeError: The kinetic energy has exceeded 1000 kcal mol-1 per atom (it is 2.2202087996265908e+16 kcal mol-1 atom-1, and +2.701328046505673e+20 kcal mol-1 total). This suggests that the simulation has become unstable. Try reducing the timestep and/or minimising +the system and run again. + +but do expect to see a ``ParticleException`` or other ``RuntimeError`` +exceptions raised at some point! + +.. note:: + + Setting ``constraint`` to the Python ``None`` type indicates that you + don't have a preference, and so :mod:`sire` will choose a suitable + constraint for you. Use the string ``"none"`` to explicitly disable + constraints. + +Constraints and Perturbable Molecules +------------------------------------- + +While constraints are useful for speeding up simulations, they can cause +problems when used with perturbable (merged) molecules. This is because the +constraints hold the bonds and/or angles at fixed values based on the +starting coordinates of the simulation. Changes in λ, which may change the +equilibrium bond and angle parameters, will not be reflected in the +free energy. This is because the constraints will stop dynamics from sampling +these perturbing bonds and angles. + +For example, changing the value of λ to 1.0, and sampling with bond and angle +constraints on would force methanol to adopt ethane's internal geometry. + +>>> print(mols[0].bond("element C", "element C").length()) +1.53625 Å +>>> d = mols.dynamics(timestep="1fs", temperature="25oC", lambda_value=1.0) +>>> d.run("5ps") +>>> print(d.commit()[0].bond("element C", "element C").length()) +1.48737 Å +>>> d = mols.dynamics(timestep="2fs", constraint="bonds-h-angles", +... temperature="25oC", lambda_value=1.0) +>>> d.run("5ps") +>>> print(d.commit()[0].bond("element C", "element C").length()) +1.53625 Å + +As seen here, the C-C bond length for ethane is 1.54 Å, while the C-O +bond length for methanol is 1.49 Å. Using ``bonds-h-angles`` constrains +this bond, meaning that the simulation at λ=1 uses ethane's bond length +(1.54 Å) rather than methanol's (1.49 Å). + +.. note:: + + The C-C bond morphs into the C-O bond during the perturbation from ethane + to methanol. Earlier, we mapped the default parameters to those of + ethane, meaning that the elements property of ethane is used by + default. This is why we searched for ``bond("element C", "element C")`` + rather than ``bond("element C", "element O")``. We would use + ``bond("element C", "element O")`` if we had used + :meth:`~sire.mol.Perturbation.link_to_perturbed` to set the + default properties. + +One solution is to choose a different constraint for perturbable molecules +than for the rest of the system. You can do this using the +``perturbable_constraint`` keyword, e.g. + +>>> d = mols.dynamics(timestep="2fs", +... constraint="bonds-h-angles", +... perturbable_constraint="none", +... temperature="25oC", +... lambda_value=1.0) +>>> d.run("5ps") +>>> print(d.commit()[0].bond("element C", "element C").length()) +1.42569 Å + +has run dynamics using no constraints on the perturbable molecules, +and ``bonds-h-angles`` constraints on all other molecules. This has +allowed sampling of the C-O bond in methanol, so that it was able to +vibrate around its equilibrium bond length (1.42 Å). + +The ``perturbable_constraint`` argument accepts the same values as +``constraint``, i.e. ``none``, ``h-bonds``, ``h-bond-h-angles`` etc. + +.. note:: + + By default, ``perturbable_constraint`` will have the same value + as ``constraint``. + +.. note:: + + You must set ``perturbable_constraint`` to the string ``none`` if you want + to disable constraints. Setting ``perturbable_constraint`` to the Python + ``None`` indicates you are providing no preference for any constraint, + and so the code automatically assigns ``perturbable_constraint`` as + equal to ``constraint``. + +Unfortunately, not constraining the bonds and/or angles of the perturbable +molecules will impact the stability of dynamics, and thus the size of +timestep that will be achievable. For example, + +>>> d = mols.dynamics(timestep="4fs", +... constraint="h-bonds-h-angles", +... perturbable_constraint="none", +... temperature="25oC", +... lambda_value=1.0) +>>> d.run("5ps") +OpenMMException: Particle coordinate is NaN. For more information, see +https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#nan + +.. note:: + + You can still use constraints on perturbable molecules. Just be careful + to minimise and then equilibrate the molecule(s) at the desired value + of λ without using constraints, so that the perturbable bonds and angles + have the right size for that value of λ. You can then run longer simulations + with constraints applied, as they will use the bond / angle sizes + measured from these starting coordinates as the constrained values. + +Hydrogen Mass Repartitioning +---------------------------- + +Bonds involving hydrogen atoms vibrate quickly because vibrational frequency +is related to atomic mass - the lighter the atom, the faster it will +vibrate. We can reduce the frequency of these vibrations by increasing the +mass of the hydrogens. Fortunately, free energy is derived from the +potential energy of the molecules, which is independent of their mass. +So, we are free to magically move mass from heavy atoms such as carbon to +their bonded hydrogen atoms without affecting the free energy. + +This method, called "hydrogen mass repartitioning", is implemented in +the :func:`sire.morph.repartition_hydrogen_masses` function. It takes a +molecule as argument, and returns that same molecule with its hydrogen +masses repartitioned. + +>>> mol = mols.molecule("molecule property is_perturbable") +>>> repartitioned_mol = sr.morph.repartition_hydrogen_masses(mol) +>>> for atom0, atom1 in zip(mol.atoms(), repartitioned_mol.atoms()): +... print(atom0, atom0.property("mass"), atom1.property("mass")) +Atom( C1:1 [ 25.71, 24.94, 25.25] ) 12.01 g mol-1 10.498 g mol-1 +Atom( C2:2 [ 24.29, 25.06, 24.75] ) 12.01 g mol-1 10.498 g mol-1 +Atom( H3:3 [ 25.91, 23.89, 25.56] ) 1.008 g mol-1 1.512 g mol-1 +Atom( H4:4 [ 26.43, 25.22, 24.45] ) 1.008 g mol-1 1.512 g mol-1 +Atom( H5:5 [ 25.86, 25.61, 26.13] ) 1.008 g mol-1 1.512 g mol-1 +Atom( H6:6 [ 24.14, 24.39, 23.87] ) 1.008 g mol-1 1.512 g mol-1 +Atom( H7:7 [ 24.09, 26.11, 24.44] ) 1.008 g mol-1 1.512 g mol-1 +Atom( H8:8 [ 23.57, 24.78, 25.55] ) 1.008 g mol-1 1.512 g mol-1 + +This has repartioned the hydrogen masses using a ``mass_factor`` of 1.5. +This factor is known to be good for using the (default) +``LangevinMiddleIntegrator`` with a timestep of 4 fs. You may need to use +a different factor for a different integrator or timestep, e.g. + +>>> rep4_mol = sr.morph.repartition_hydrogen_masses(mol, mass_factor=4) +>>> for atom0, atom1 in zip(mol.atoms(), rep4_mol.atoms()): +... print(atom0, atom0.property("mass"), atom1.property("mass")) +Atom( C1:1 [ 25.71, 24.94, 25.25] ) 12.01 g mol-1 2.938 g mol-1 +Atom( C2:2 [ 24.29, 25.06, 24.75] ) 12.01 g mol-1 2.938 g mol-1 +Atom( H3:3 [ 25.91, 23.89, 25.56] ) 1.008 g mol-1 4.032 g mol-1 +Atom( H4:4 [ 26.43, 25.22, 24.45] ) 1.008 g mol-1 4.032 g mol-1 +Atom( H5:5 [ 25.86, 25.61, 26.13] ) 1.008 g mol-1 4.032 g mol-1 +Atom( H6:6 [ 24.14, 24.39, 23.87] ) 1.008 g mol-1 4.032 g mol-1 +Atom( H7:7 [ 24.09, 26.11, 24.44] ) 1.008 g mol-1 4.032 g mol-1 +Atom( H8:8 [ 23.57, 24.78, 25.55] ) 1.008 g mol-1 4.032 g mol-1 + +would multiply the hydrogen masses by a factor of 4. + +.. note:: + + By default, this function will repartition the masses of the hydrogens + in both the standard mass property (``mass``) and also for the end-state + mass properties (``mass0`` and ``mass1`` if they exist). You can disable + repartitioning of the end-state masses by setting ``include_end_states`` + to ``False``. You can choose to repartition specific mass properties by + passing in a property map, e.g. ``map={"mass": "mass0"}``. + +The repartitioned molecule has the same mass as the original molecule, but +the hydrogens have been made heavier. The mass of the carbon atoms has been +reduced to compensate. + +>>> print(mol.mass(), repartitioned_mol.mass()) +30.068 g mol-1 30.068 g mol-1 + +It is normal to only repartition the hydrogen masses of perturbable molecules. +This lets us disable the constraints on the perturbable molecules, but +still use a large timestep (e.g. 3 fs or 4 fs). This is because +the hydrogens on the perturbable molecules would become heavier, +and so the vibrations of those atoms should have a lower frequency. + +>>> mols.update(repartitioned_mol) +>>> d = mols.dynamics(timestep="4fs", temperature="25oC", +... constraint="h-bonds", +... perturbable_constraint="none") +>>> d.run("5ps") +>>> print(d) +Dynamics(completed=5 ps, energy=-31950.1 kcal mol-1, speed=100.8 ns day-1) + +This has given us the best of both worlds - a fast simulation with a larger +timestep, plus no constraints on the perturbable molecules. + +Using this protocol, we can now recalculate the relative free energy of +ethane and methanol. + +>>> for l in range(0, 105, 5): +... # turn l into the lambda value by dividing by 100 +... lambda_value = l / 100.0 +... print(f"Simulating lambda={lambda_value:.2f}") +... # minimise the system at this lambda value +... min_mols = mols.minimisation(lambda_value=lambda_value).run().commit() +... # create a dynamics object for the system +... d = min_mols.dynamics(timestep="4fs", temperature="25oC", +... lambda_value=lambda_value, +... constraint="h-bonds", +... perturbable_constraint="none") +... # generate random velocities +... d.randomise_velocities() +... # equilibrate, not saving anything +... d.run("2ps", save_frequency=0) +... print("Equilibration complete") +... print(d) +... # get the values of lambda for neighbouring windows +... lambda_windows = [lambda_value] +... if lambda_value > 0: +... lambda_windows.insert(0, (l-5)/100.0) +... if lambda_value < 1: +... lambda_windows.append((l+5)/100.0) +... # run the dynamics, saving the energy every 0.1 ps +... d.run("25ps", energy_frequency="0.1ps", frame_frequency=0, +... lambda_windows=lambda_windows) +... print("Dynamics complete") +... print(d) +... # stream the EnergyTrajectory to a sire save stream object +... sr.stream.save(d.commit().energy_trajectory(), +... f"energy_fast_{lambda_value:.2f}.s3") + +The simulation runs 33% faster than before, taking about 30 seconds per +λ-window. + +We can now calculate the free energy using alchemlyb as before; + +>>> df = sr.morph.to_alchemlyb("energy_fast_*.s3") +>>> from alchemlyb.estimators import BAR +>>> b = BAR() +>>> b.fit(df) +>>> print(b.delta_f_.loc[0.00, 1.00]) +-2.9972763590251836 + +Within error, this is the same free energy as before, but calculated in +a little less time. + +.. note:: + + The random error on these calculations can be seen in ``b.d_delta_f``, + which shows the error per λ-window is about 0.02 kcal mol-1. Summed over + the whole λ-coordinate suggest an error of about 0.4 kcal mol-1. You + could reduce this error by running the simulation for longer + (e.g. 250 ps per λ-window) or by running multiple repeats and taking + and average. + +.. _ExampleFEPScript: + +Complete Example Script +----------------------- + +Putting everything together, here is a simple script that does all of the +work of calculating the relative hydration free energy of ethane +and methanol. The key parameters (e.g. timestep, constraint, run time, +λ-values etc) are pulled out as variables at the top. The script then +runs a dynamics simulation for each λ-window for the water leg, then +a dynamics simulation using the same parameters and λ-windows for the +vacuum leg. The free energies are collected and then calculated using BAR +from alchemlyb. This is a good starting point for you to adapt for your +own simulations. Or take a look at +`BioSimSpace `__ or the upcoming +`somd2 software `__ if you are +interested in higher-level interfaces to this functionality that +automatically run more complex protocols. + +.. literalinclude:: scripts/run_md.py + +The relative hydration free energy calculated using the above script is: + +:: + + Water phase: -3.2587718667371606 kcal mol-1 + Vacuum phase: 2.980451221216327 kcal mol-1 + Hydration free energy: -6.239223087953487 kcal mol-1 + +This is in good agreement with +`published results from other codes `__ +which are typically -5.99 kcal mol-1 to -6.26 kcal mol-1. + +.. note:: + + There will be some variation between different codes and different + protocols, as the convergence of the free energy estimate is sensitive + to the length of the dynamics simulation at each λ-value. In this case, + we used very short simulations. Edit the script to increase the amount + of simulation time per λ-value to get a more accurate result. + +This script can be edited to run a longer, more statistically accurate +by setting ``run_time`` to a larger value, e.g. 250 ps. On my machine, +this results in a simulation that takes 5 minutes per λ-window for the +water phase at 100 ns day-1, and 45 seconds per λ-window for the vacuum +phase, at 950 ns day-1 (total run time about 120 minutes). The resulting +free energy is in remarkable agreement with that from the shorter +simulation. + +:: + + Water phase: -3.3027026336540715 kcal mol-1 + Vacuum phase: 2.9761399609273567 kcal mol-1 + Hydration free energy: -6.278842594581429 kcal mol-1 + +We can check the convergence and get an error estimate as before. First +the water phase; + +>>> import sire as sr +>>> from glob import glob +>>> from alchemlyb.convergence import forward_backward_convergence +>>> from alchemlyb.visualisation import plot_convergence +>>> dfs = [] +>>> s3files = glob("energy_water*.s3") +>>> s3files.sort() +>>> for s3 in s3files: +... dfs.append(sr.stream.load(s3).to_alchemlyb()) +>>> f = forward_backward_convergence(dfs, "bar") +>>> plot_convergence(f) + +.. note:: + + It is important that the dataframes are passed to alchemlyb in + ascending order of λ-value, or else it will calculate the wrong + free energy. This is why we had to sort the result of globbing + the files, and why a file naming scheme was chosen that names + the files in ascending λ order. + +.. image:: images/06_06_01.jpg + :alt: Convergence of the free energy estimate for the water phase. + +And now the vacuum phase; + +>>> dfs = [] +>>> s3files = glob("energy_vacuum*.s3") +>>> s3files.sort() +>>> for s3 in s3files: +... dfs.append(sr.stream.load(s3).to_alchemlyb()) +>>> f = forward_backward_convergence(dfs, "bar") +>>> plot_convergence(f) + +.. image:: images/06_06_02.jpg + :alt: Convergence of the free energy estimate for the vacuum phase. + +From the convergence plots and the output printed by alchemlyb, we can +see that the error on the water leg is 0.02 kcal mol-1, and the error +on the vacuum leg is 0.01 kcal mol-1. Summed in quadrature, this gives +an error of 0.022 kcal mol-1 on the hydration free energy. This gives +a final result of -6.279 ± 0.022 kcal mol-1, which is in good agreement with +`published results from other codes `__ +which are in the range of -5.99 kcal mol-1 to -6.26 kcal mol-1. diff --git a/doc/source/tutorial/part06/images/06_01_01.jpg b/doc/source/tutorial/part06/images/06_01_01.jpg new file mode 100644 index 000000000..3b89213fe Binary files /dev/null and b/doc/source/tutorial/part06/images/06_01_01.jpg differ diff --git a/doc/source/tutorial/part06/images/06_01_02.jpg b/doc/source/tutorial/part06/images/06_01_02.jpg new file mode 100644 index 000000000..604a84aaf Binary files /dev/null and b/doc/source/tutorial/part06/images/06_01_02.jpg differ diff --git a/doc/source/tutorial/part06/images/06_02_01.jpg b/doc/source/tutorial/part06/images/06_02_01.jpg new file mode 100644 index 000000000..8674ec7aa Binary files /dev/null and b/doc/source/tutorial/part06/images/06_02_01.jpg differ diff --git a/doc/source/tutorial/part06/images/06_02_02.jpg b/doc/source/tutorial/part06/images/06_02_02.jpg new file mode 100644 index 000000000..5618f084b Binary files /dev/null and b/doc/source/tutorial/part06/images/06_02_02.jpg differ diff --git a/doc/source/tutorial/part06/images/06_02_03.jpg b/doc/source/tutorial/part06/images/06_02_03.jpg new file mode 100644 index 000000000..44694b4c6 Binary files /dev/null and b/doc/source/tutorial/part06/images/06_02_03.jpg differ diff --git a/doc/source/tutorial/part06/images/06_02_04.jpg b/doc/source/tutorial/part06/images/06_02_04.jpg new file mode 100644 index 000000000..56af4e349 Binary files /dev/null and b/doc/source/tutorial/part06/images/06_02_04.jpg differ diff --git a/doc/source/tutorial/part06/images/06_05_01.jpg b/doc/source/tutorial/part06/images/06_05_01.jpg new file mode 100644 index 000000000..7ec45085b Binary files /dev/null and b/doc/source/tutorial/part06/images/06_05_01.jpg differ diff --git a/doc/source/tutorial/part06/images/06_06_01.jpg b/doc/source/tutorial/part06/images/06_06_01.jpg new file mode 100644 index 000000000..39df8cc89 Binary files /dev/null and b/doc/source/tutorial/part06/images/06_06_01.jpg differ diff --git a/doc/source/tutorial/part06/images/06_06_02.jpg b/doc/source/tutorial/part06/images/06_06_02.jpg new file mode 100644 index 000000000..c1046088d Binary files /dev/null and b/doc/source/tutorial/part06/images/06_06_02.jpg differ diff --git a/doc/source/tutorial/part06/scripts/run_md.py b/doc/source/tutorial/part06/scripts/run_md.py new file mode 100644 index 000000000..f30b8ad0a --- /dev/null +++ b/doc/source/tutorial/part06/scripts/run_md.py @@ -0,0 +1,136 @@ +import sire as sr + +mols = sr.load(sr.expand(sr.tutorial_url, "merged_molecule.s3")) + +for mol in mols.molecules("molecule property is_perturbable"): + mol = mol.perturbation().link_to_reference().commit() + mol = sr.morph.repartition_hydrogen_masses(mol, mass_factor=1.5) + mols.update(mol) + +mols = mols.minimisation().run().commit() + +timestep = "4fs" +energy_frequency = "0.1ps" + +constraint = "h-bonds" +perturbable_constraint = None + +equil_time = "2ps" +run_time = "25ps" + +lambda_values = [x / 100.0 for x in range(0, 101, 5)] + +print("\nWater leg") + +for i, lambda_value in enumerate(lambda_values): + print(f"Simulating lambda={lambda_value:.2f}") + # minimise the system at this lambda value + min_mols = mols.minimisation(lambda_value=lambda_value).run().commit() + + # create a dynamics object for the system + d = min_mols.dynamics( + timestep=timestep, + temperature="25oC", + lambda_value=lambda_value, + constraint=constraint, + perturbable_constraint=perturbable_constraint, + ) + + # generate random velocities + d.randomise_velocities() + + # equilibrate, not saving anything + d.run(equil_time, save_frequency=0) + print("Equilibration complete") + print(d) + + # get the values of lambda for neighbouring windows + lambda_windows = lambda_values[ + max(i - 1, 0) : min(len(lambda_values), i + 2) + ] + + # run the dynamics, saving the energy every 0.1 ps + d.run( + run_time, + energy_frequency=energy_frequency, + frame_frequency=0, + lambda_windows=lambda_windows, + ) + print("Dynamics complete") + print(d) + + # stream the EnergyTrajectory to a sire save stream object + sr.stream.save( + d.commit().energy_trajectory(), f"energy_water_{lambda_value:.2f}.s3" + ) + +print("\nVacuum leg") + +for i, lambda_value in enumerate(lambda_values): + print(f"Simulating lambda={lambda_value:.2f}") + # minimise the system at this lambda value + min_mols = ( + mols[0] + .minimisation(lambda_value=lambda_value, vacuum=True) + .run() + .commit() + ) + + # create a dynamics object for the system + d = min_mols.dynamics( + timestep=timestep, + temperature="25oC", + lambda_value=lambda_value, + constraint=constraint, + perturbable_constraint=perturbable_constraint, + vacuum=True, + ) + + # generate random velocities + d.randomise_velocities() + + # equilibrate, not saving anything + d.run(equil_time, save_frequency=0) + print("Equilibration complete") + print(d) + + # get the values of lambda for neighbouring windows + lambda_windows = lambda_values[ + max(i - 1, 0) : min(len(lambda_values), i + 2) + ] + + # run the dynamics, saving the energy every 0.1 ps + d.run( + run_time, + energy_frequency=energy_frequency, + frame_frequency=0, + lambda_windows=lambda_windows, + ) + print("Dynamics complete") + print(d) + + # stream the EnergyTrajectory to a sire save stream object + sr.stream.save( + d.commit().energy_trajectory(), + f"energy_vacuum_{lambda_value:.2f}.s3", + ) + +from alchemlyb.estimators import BAR + +df = sr.morph.to_alchemlyb("energy_water_*.s3") + +b = BAR() +b.fit(df) + +dG_solv = b.delta_f_.loc[0.00, 1.00] + +df = sr.morph.to_alchemlyb("energy_vacuum_*.s3") + +b = BAR() +b.fit(df) + +dG_vac = b.delta_f_.loc[0.00, 1.00] + +print(f"Water phase: {dG_solv} kcal mol-1") +print(f"Vacuum phase: {dG_vac} kcal mol-1") +print(f"Hydration free energy: {dG_solv - dG_vac} kcal mol-1") diff --git a/src/sire/_pythonize.py b/src/sire/_pythonize.py index bbfbe2612..5fe58c5dc 100644 --- a/src/sire/_pythonize.py +++ b/src/sire/_pythonize.py @@ -155,8 +155,10 @@ def _pythonize_modules(modules, delete_old: bool = True): _is_using_old_api = None _is_using_new_api = None +_is_in_loading_process = False -def _load_new_api_modules(delete_old: bool = True): + +def _load_new_api_modules(delete_old: bool = True, is_base: bool = False): """ Internal function to load the new API modules, pythonizing the function names as we go. If `delete_old` is True, then @@ -167,8 +169,16 @@ def _load_new_api_modules(delete_old: bool = True): _is_using_new_api = True + global _is_in_loading_process + + if _is_in_loading_process: + return + + _is_in_loading_process = True + # call Pythonize on all of the new modules from .legacy import ( + Base, Move, IO, System, @@ -177,7 +187,6 @@ def _load_new_api_modules(delete_old: bool = True): FF, Mol, Analysis, - Base, CAS, Cluster, Convert, # does not need pythonizing, but importing will make it visible @@ -192,8 +201,8 @@ def _load_new_api_modules(delete_old: bool = True): _pythonize_modules( [ - Analysis._Analysis, Base._Base, + Analysis._Analysis, CAS._CAS, Cluster._Cluster, Error._Error, @@ -223,6 +232,18 @@ def _load_new_api_modules(delete_old: bool = True): if have_lazy_import: # Now make sure that all new modules have been loaded + # (we need to import base first) + from . import base + + if lazy_import.LazyModule in type(base).mro(): + # this module is lazily loaded - use 'dir' to load it + dir(base) + + if is_base: + # return, as we will only import base here + _is_in_loading_process = False + return + from . import ( move, io, @@ -233,7 +254,6 @@ def _load_new_api_modules(delete_old: bool = True): ff, mol, analysis, - base, cas, cluster, error, @@ -257,7 +277,6 @@ def _load_new_api_modules(delete_old: bool = True): ff, mol, analysis, - base, cas, cluster, error, @@ -274,6 +293,8 @@ def _load_new_api_modules(delete_old: bool = True): # this module is lazily loaded - use 'dir' to load it dir(M) + _is_in_loading_process = False + def use_mixed_api(support_old_module_names: bool = False): """Load Sire using both the new (python-style) and old APIs. This @@ -313,10 +334,19 @@ def use_mixed_api(support_old_module_names: bool = False): _load_new_api_modules(delete_old=False) -def use_new_api(): - """Load Sire using the new (python-style) API. This will be called +def use_new_api(is_base: bool = False): + """ + Load Sire using the new (python-style) API. This will be called automatically when you load any of the new Python modules, so you shouldn't need to call this yourself. + + Parameters + ---------- + + is_base: bool (defaults to False) + Whether or not this is being called by the sire.base module. + This triggers a special case where we only load sire.base, + and thus avoid circular imports """ global _is_using_new_api, _is_using_old_api @@ -340,7 +370,7 @@ def use_new_api(): raise ImportError(msg) # Now bring in the new API - _load_new_api_modules(delete_old=True) + _load_new_api_modules(delete_old=True, is_base=is_base) def use_old_api(): diff --git a/src/sire/base/__init__.py b/src/sire/base/__init__.py index e9b2a17eb..8d67cddac 100644 --- a/src/sire/base/__init__.py +++ b/src/sire/base/__init__.py @@ -8,7 +8,7 @@ from ._progressbar import * -_use_new_api() +_use_new_api(is_base=True) wrap = _Base.wrap @@ -107,11 +107,13 @@ def __str__(obj): def __propertymap_set(obj, key, value): try: obj.__orig__set(key, value) + return except Exception: pass try: obj.__orig__set(key, _Base.PropertyName(value)) + return except Exception: pass diff --git a/src/sire/maths/__init__.py b/src/sire/maths/__init__.py index cbec5949c..48efd90aa 100644 --- a/src/sire/maths/__init__.py +++ b/src/sire/maths/__init__.py @@ -170,23 +170,133 @@ def create_quaternion(angle=None, axis=None, matrix=None, quaternion=None): if not hasattr(EnergyTrajectory, "to_pandas"): - def _to_pandas(obj): + def _to_pandas(obj, temperature=None, to_alchemlyb: bool = False): """ Return the energy trajectory as a pandas DataFrame + + Parameters + ---------- + + temperature: temperature + The temperature of the simulation. If this is + not set then the temperature from this table's + `ensemble` or `temperature` property will be + used. Note that you only need a temperature + if you are converting to alchemlyb format. + + to_alchemlyb: bool + This will format the DataFrame in a way that is + compatible with alchemlyb. This will allow the + DataFrame to be used as part of an alchemlyb + free energy calculation. """ import pandas as pd from ..units import picosecond, kcal_per_mol data = {} - data["time"] = obj.times(picosecond.get_default()) + if to_alchemlyb: + time_unit = picosecond + time_unit_string = "ps" + + energy_unit = kcal_per_mol + energy_unit_string = "kcal/mol" + + if temperature is None: + # look for the temperature in the ensemble property + if obj.has_property("ensemble"): + temperature = obj.property("ensemble").temperature() + + # ok, try the temperature property + if temperature is None and obj.has_property("temperature"): + temperature = obj.property("temperature") + + if temperature is None: + raise ValueError( + "You must specify the temperature of the simulation " + "when converting to alchemlyb format, or ensure that " + "the trajectory has an ensemble or temperature " + "property." + ) + else: + time_unit = picosecond.get_default() + time_unit_string = time_unit.unit_string() + + energy_unit = kcal_per_mol.get_default() + energy_unit_string = energy_unit.unit_string() + + data["time"] = obj.times(time_unit) + + keys = obj.label_keys() + keys.sort() + + for key in keys: + if to_alchemlyb and key == "lambda": + data["fep-lambda"] = obj.labels_as_numbers(key) + else: + # use float keys if possible + try: + column_header = float(key) + except Exception: + column_header = key + + try: + data[column_header] = obj.labels_as_numbers(key) + except Exception: + data[column_header] = obj.labels(key) keys = obj.keys() keys.sort() + if to_alchemlyb: + keys.remove("kinetic") + keys.remove("potential") + for key in keys: - data[str(key)] = obj.energies(key, kcal_per_mol.get_default()) + # use float keys if possible + try: + column_header = float(key) + except Exception: + column_header = key + + data[column_header] = obj.energies(key, energy_unit) - return pd.DataFrame(data).set_index("time") + if to_alchemlyb: + df = pd.DataFrame(data).set_index(["time", "fep-lambda"]) + else: + df = pd.DataFrame(data).set_index("time") + + if temperature is not None: + from .. import u + from ..units import kelvin + + df.attrs["temperature"] = u(temperature).to(kelvin) + + df.attrs["energy_unit"] = energy_unit_string + df.attrs["time_unit"] = time_unit_string + + return df + + def _to_alchemlyb(obj, temperature=None): + """ + Return the energy trajectory as an alchemlyb-formatted pandas DataFrame + + Parameters + ---------- + + temperature: temperature + The temperature of the simulation. If this is + not set then the temperature from this table's + `ensemble` or `temperature` property will be + used. + + Returns + ------- + + pandas.DataFrame + A pandas DataFrame that is compatible with alchemlyb. + """ + return obj.to_pandas(temperature=temperature, to_alchemlyb=True) EnergyTrajectory.to_pandas = _to_pandas + EnergyTrajectory.to_alchemlyb = _to_alchemlyb diff --git a/src/sire/mol/__init__.py b/src/sire/mol/__init__.py index 29de1a105..5c6129cea 100644 --- a/src/sire/mol/__init__.py +++ b/src/sire/mol/__init__.py @@ -1405,19 +1405,25 @@ def _dynamics( save_frequency=None, energy_frequency=None, frame_frequency=None, + save_velocities=None, constraint=None, + perturbable_constraint=None, + include_constrained_energies: bool = True, schedule=None, lambda_value=None, swap_end_states=None, + ignore_perturbations=None, temperature=None, pressure=None, + vacuum=None, shift_delta=None, coulomb_power=None, restraints=None, fixed=None, + platform=None, device=None, + precision=None, map=None, - **kwargs, ): """ Return a Dynamics object that can be used to perform @@ -1452,6 +1458,12 @@ def _dynamics( to save frames during the trajectory. This can be overridden by setting frame_frequency during an individual run. + save_velocities: bool + Whether or not to save velocities when saving trajectory frames + during the simulation. This defaults to False, as velocity + trajectories aren't often needed, and they double the amount + of space that is required for a trajectory. + constraint: str The type of constraint to use for bonds and/or angles, e.g. `h-bonds`, `bonds` etc. @@ -1459,6 +1471,19 @@ def _dynamics( for the full list of options. This will be automatically guessed from the timestep if it isn't set. + perturbable_constraint: str + The type of constraint to use for perturbable bonds and/or angles, + e.g. `h-bonds`, `bonds` etc. + See https://sire.openbiosim.org/cheatsheet/openmm.html#choosing-options + for the full list of options. This equal the value of `constraint` + if it isn't set. + + include_constrained_energies: bool + Whether or not to include the energies of the constrained bonds + and angles. If this is False, then the internal bond or angle + energy of the constrained degrees of freedom are not included + in the total energy, and their forces are not evaluated. + schedule: sire.cas.LambdaSchedule The schedule used to control how perturbable forcefield parameters should be morphed as a function of lambda. If this is not set @@ -1478,6 +1503,14 @@ def _dynamics( use the coordinates of the perturbed molecule as the starting point. + ignore_perturbations: bool + Whether or not to ignore perturbations. If this is True, then + the perturbation will be ignored, and the simulation will + be run using the properties of the reference molecule + (or the perturbed molecule if swap_end_states is True). This + is useful if you just want to run standard molecular dynamics + of the reference or perturbed states. + temperature: temperature The temperature at which to run the simulation. A microcanonical (NVE) simulation will be run if you don't @@ -1488,6 +1521,12 @@ def _dynamics( microcanonical (NVE) or canonical (NVT) simulation will be run if the pressure is not set. + vacuum: bool (optional) + Whether or not to run the simulation in vacuum. If this is + set to `True`, then the simulation space automatically be + replaced by a `sire.vol.Cartesian` space, and the + simulation run in vacuum. + shift_delta: length The shift_delta parameter that controls the electrostatic and van der Waals softening potential that smooths the @@ -1509,11 +1548,19 @@ def _dynamics( that should be fixed in place during the simulation. These atoms will not be moved by dynamics. + platform: str + The name of the OpenMM platform on which to run the dynamics, + e.g. "CUDA", "OpenCL", "Metal" etc. + device: str or int The ID of the GPU (or accelerator) used to accelerate the simulation. This would be CUDA_DEVICE_ID or similar if CUDA was used. This can be any valid OpenMM device string + precision: str + The desired precision for the simulation (e.g. `single`, + `mixed` or `double`) + map: dict A dictionary of additional options. Note that any options set in this dictionary that are also specified via one of @@ -1524,7 +1571,12 @@ def _dynamics( from ..system import System from .. import u - map = create_map(map, kwargs) + map = create_map(map) + + if vacuum is True: + from ..vol import Cartesian + + map.set("space", Cartesian()) if not map.specified("space"): map = create_map(map, {"space": "space"}) @@ -1538,7 +1590,12 @@ def _dynamics( # space is not a shared property, so may be lost when we # convert to molecules. Make sure this doens't happen by # adding the space directly to the property map - map.set("space", view.property(map["space"])) + try: + map.set("space", view.property(map["space"])) + except Exception: + from ..vol import Cartesian + + map.set("space", Cartesian()) # Set default values if these have not been set if cutoff is None and not map.specified("cutoff"): @@ -1566,15 +1623,18 @@ def _dynamics( if save_frequency is None and not map.specified("save_frequency"): from ..units import picosecond - save_frequency = 25 * picosecond - else: - save_frequency = u(save_frequency) + map.set("save_frequency", 25 * picosecond) + elif save_frequency is not None: + map.set("save_frequency", u(save_frequency)) if energy_frequency is not None: - map.set("energy_frequency", energy_frequency) + map.set("energy_frequency", u(energy_frequency)) if frame_frequency is not None: - map.set("frame_frequency", frame_frequency) + map.set("frame_frequency", u(frame_frequency)) + + if save_velocities is not None: + map.set("save_velocities", save_velocities) if constraint is None and not map.specified("constraint"): from ..units import femtosecond @@ -1583,18 +1643,25 @@ def _dynamics( # it must be in the map timestep = map["timestep"].value() - if timestep > 2 * femtosecond: + if timestep > 4 * femtosecond: # need constraint on everything constraint = "bonds-h-angles" + elif timestep > 2 * femtosecond: + # need constraint on everything + constraint = "h-bonds-h-angles" + elif timestep > 1 * femtosecond: # need it just on H bonds and angles - constraint = "h-bonds-h-angles" + constraint = "h-bonds" else: # can get away with no constraints constraint = "none" + if perturbable_constraint is not None: + perturbable_constraint = str(perturbable_constraint).lower() + if temperature is not None: temperature = u(temperature) map.set("temperature", temperature) @@ -1606,18 +1673,28 @@ def _dynamics( if device is not None: map.set("device", str(device)) + if precision is not None: + map.set("precision", str(precision).lower()) + + if include_constrained_energies is not None: + map.set("include_constrained_energies", include_constrained_energies) + + if platform is not None: + map.set("platform", str(platform)) + return Dynamics( view, cutoff=cutoff, cutoff_type=cutoff_type, timestep=timestep, - save_frequency=save_frequency, constraint=str(constraint), + perturbable_constraint=perturbable_constraint, schedule=schedule, lambda_value=lambda_value, shift_delta=shift_delta, coulomb_power=coulomb_power, swap_end_states=swap_end_states, + ignore_perturbations=ignore_perturbations, restraints=restraints, fixed=fixed, map=map, @@ -1629,16 +1706,21 @@ def _minimisation( cutoff=None, cutoff_type=None, constraint=None, + perturbable_constraint=None, + include_constrained_energies: bool = True, schedule=None, lambda_value=None, swap_end_states=None, + ignore_perturbations=None, + vacuum=None, shift_delta=None, coulomb_power=None, + platform=None, device=None, + precision=None, restraints=None, fixed=None, map=None, - **kwargs, ): """ Return a Minimisation object that can be used to perform @@ -1656,8 +1738,20 @@ def _minimisation( The type of constraint to use for bonds and/or angles, e.g. `h-bonds`, `bonds` etc. See https://sire.openbiosim.org/cheatsheet/openmm.html#choosing-options - for the full list of options. This will be automatically - guessed from the timestep if it isn't set. + for the full list of options. This is `none` if it is not set. + + perturbable_constraint: str + The type of constraint to use for perturbable bonds and/or angles, + e.g. `h-bonds`, `bonds` etc. + See https://sire.openbiosim.org/cheatsheet/openmm.html#choosing-options + for the full list of options. This equal the value of `constraint` + if it isn't set. + + include_constrained_energies: bool + Whether or not to include the energies of the perturbable bonds + and angles. If this is False, then the internal bond or angle + energy of the perturbable degrees of freedom are not included + in the total energy, and their forces are not evaluated. schedule: sire.cas.LambdaSchedule The schedule used to control how perturbable forcefield parameters @@ -1678,6 +1772,20 @@ def _minimisation( use the coordinates of the perturbed molecule as the starting point. + ignore_perturbations: bool + Whether or not to ignore perturbations. If this is True, then + the perturbation will be ignored, and the simulation will + be run using the properties of the reference molecule + (or the perturbed molecule if swap_end_states is True). This + is useful if you just want to run standard molecular dynamics + of the reference or perturbed states. + + vacuum: bool (optional) + Whether or not to run the simulation in vacuum. If this is + set to `True`, then the simulation space automatically be + replaced by a `sire.vol.Cartesian` space, and the + simulation run in vacuum. + shift_delta: length The shift_delta parameter that controls the electrostatic and van der Waals softening potential that smooths the @@ -1699,11 +1807,19 @@ def _minimisation( that should be fixed in place during the simulation. These atoms will not be moved by minimisation. + platform: str + The name of the OpenMM platform on which to run the dynamics, + e.g. "CUDA", "OpenCL", "Metal" etc. + device: str or int The ID of the GPU (or accelerator) used to accelerate minimisation. This would be CUDA_DEVICE_ID or similar if CUDA was used. This can be any valid OpenMM device string + precision: str + The desired precision for the simulation (e.g. `single`, + `mixed` or `double`) + map: dict A dictionary of additional options. Note that any options set in this dictionary that are also specified via one of @@ -1714,7 +1830,12 @@ def _minimisation( from ..system import System from .. import u - map = create_map(map, kwargs) + map = create_map(map) + + if vacuum is True: + from ..vol import Cartesian + + map.set("space", Cartesian()) if not map.specified("space"): map = create_map(map, {"space": "space"}) @@ -1728,7 +1849,12 @@ def _minimisation( # space is not a shared property, so may be lost when we # convert to molecules. Make sure this doens't happen by # adding the space directly to the property map - map.set("space", view.property(map["space"])) + try: + map.set("space", view.property(map["space"])) + except Exception: + from ..vol import Cartesian + + map.set("space", Cartesian()) # Set default values if these have not been set if cutoff is None and not map.specified("cutoff"): @@ -1749,8 +1875,20 @@ def _minimisation( if device is not None: map.set("device", str(device)) + if precision is not None: + map.set("precision", str(precision).lower()) + if constraint is not None: - map.set("constraint", str(constraint)) + map.set("constraint", str(constraint).lower()) + + if perturbable_constraint is not None: + map.set("perturbable_constraint", str(perturbable_constraint).lower()) + + if include_constrained_energies is not None: + map.set("include_constrained_energies", include_constrained_energies) + + if platform is not None: + map.set("platform", str(platform)) return Minimisation( view, @@ -1759,6 +1897,7 @@ def _minimisation( schedule=schedule, lambda_value=lambda_value, swap_end_states=swap_end_states, + ignore_perturbations=ignore_perturbations, shift_delta=shift_delta, coulomb_power=coulomb_power, restraints=restraints, diff --git a/src/sire/mol/_dynamics.py b/src/sire/mol/_dynamics.py index 3bd1c8e75..17cad4472 100644 --- a/src/sire/mol/_dynamics.py +++ b/src/sire/mol/_dynamics.py @@ -29,7 +29,7 @@ def __init__(self, mols=None, map=None, **kwargs): from . import selection_to_atoms # turn the fixed atoms into a list of atoms - map.set( + self._map.set( "fixed", mols.atoms().find(selection_to_atoms(mols, fixed_atoms)), ) @@ -60,6 +60,9 @@ def __init__(self, mols=None, map=None, **kwargs): to_pandas=False, map=self._map ) + # make sure that the ensemble is recorded in the trajectory + self._energy_trajectory.set_property("ensemble", self.ensemble()) + self._num_atoms = len(self._sire_mols.atoms()) from ..units import nanosecond @@ -209,10 +212,10 @@ def _exit_dynamics_block( * kcal_per_mol ) - if lambda_windows is not None: - sim_lambda_value = self._omm_mols.get_lambda() - nrgs[str(sim_lambda_value)] = nrgs["potential"] + sim_lambda_value = self._omm_mols.get_lambda() + nrgs[str(sim_lambda_value)] = nrgs["potential"] + if lambda_windows is not None: for lambda_value in lambda_windows: if lambda_value != sim_lambda_value: self._omm_mols.set_lambda(lambda_value) @@ -225,7 +228,9 @@ def _exit_dynamics_block( self._omm_mols.set_lambda(sim_lambda_value) - self._energy_trajectory.set(self._current_time, nrgs) + self._energy_trajectory.set( + self._current_time, nrgs, {"lambda": str(sim_lambda_value)} + ) self._is_running = False @@ -273,13 +278,43 @@ def ensemble(self): return Ensemble(map=self._map) + def randomise_velocities(self, temperature=None, random_seed: int = None): + """ + Set the velocities to random values, drawn from the Boltzmann + distribution for the current temperature. + + Parameters + ---------- + + - temperature (temperature): The temperature to use. If None, then + the current temperature will be used + - random_seed (int): The random seed to use. If None, then + a random seed will be generated + """ + if self.is_null(): + return + + if temperature is not None: + self.set_temperature(temperature, rescale_velocities=False) + + from ..units import kelvin + + if random_seed is None: + self._omm_mols.setVelocitiesToTemperature( + self.ensemble().temperature().to(kelvin) + ) + else: + self._omm_mols.setVelocitiesToTemperature( + self.ensemble().temperature().to(kelvin), random_seed + ) + def set_ensemble(self, ensemble, rescale_velocities: bool = True): """ Set the ensemble for the dynamics. Note that this will only let you change the temperature and/or pressure of the ensemble. You can't change its fundemental nature. - If rescalse_velocities is True, then the velocities will + If rescale_velocities is True, then the velocities will be rescaled to the new temperature. """ if self.is_null(): @@ -303,6 +338,41 @@ def set_ensemble(self, ensemble, rescale_velocities: bool = True): self._map["pressure"] = ensemble.pressure() self._omm_mols.set_pressure(ensemble.pressure()) + def set_temperature(self, temperature, rescale_velocities: bool = True): + """ + Set the temperature for the dynamics. Note that this will only + let you change the temperature of the ensemble. + You can't change its fundemental nature. + + If rescale_velocities is True, then the velocities will be + rescaled to the new temperature. + """ + if self.is_null(): + return + + ensemble = self.ensemble() + + ensemble.set_temperature(temperature) + + self.set_ensemble( + ensemble=ensemble, rescale_velocities=rescale_velocities + ) + + def set_pressure(self, pressure): + """ + Set the pressure for the dynamics. Note that this will only + let you change the pressure of the ensemble. + You can't change its fundemental nature. + """ + if self.is_null(): + return + + ensemble = self.ensemble() + + ensemble.set_pressure(pressure) + + self.set_ensemble(ensemble=ensemble) + def constraint(self): if self.is_null(): return None @@ -312,6 +382,15 @@ def constraint(self): else: return "none" + def perturbable_constraint(self): + if self.is_null(): + return None + else: + if self._map.specified("perturbable_constraint"): + return self._map["perturbable_constraint"].source() + else: + return self.constraint() + def get_schedule(self): if self.is_null(): return None @@ -634,10 +713,14 @@ def get_steps_till_save(completed: int, total: int): completed = 0 if completed >= total: - return (0, True, True) + return ( + 0, + frame_frequency_steps > 0, + energy_frequency_steps > 0, + ) elif frame_frequency_steps <= 0 and energy_frequency_steps <= 0: - return (total, True, True) + return (total, False, False) n_to_end = total - completed @@ -825,6 +908,8 @@ def commit(self): self._energy_trajectory, map=self._map ) + self._sire_mols.set_ensemble(self.ensemble()) + def _add_extra(extras, key, value): if value is not None: @@ -846,11 +931,12 @@ def __init__( cutoff=None, cutoff_type=None, timestep=None, - save_frequency=None, constraint=None, + perturbable_constraint=None, schedule=None, lambda_value=None, swap_end_states=None, + ignore_perturbations=None, shift_delta=None, coulomb_power=None, restraints=None, @@ -867,11 +953,12 @@ def __init__( if timestep is not None: _add_extra(extras, "timestep", u(timestep)) - _add_extra(extras, "save_frequency", save_frequency) _add_extra(extras, "constraint", constraint) + _add_extra(extras, "perturbable_constraint", perturbable_constraint) _add_extra(extras, "schedule", schedule) _add_extra(extras, "lambda", lambda_value) _add_extra(extras, "swap_end_states", swap_end_states) + _add_extra(extras, "ignore_perturbations", ignore_perturbations) if shift_delta is not None: _add_extra(extras, "shift_delta", u(shift_delta)) @@ -908,7 +995,10 @@ def __str__(self): def __repr__(self): return self.__str__() - def minimise(self, max_iterations: int = 10000): + def minimise( + self, + max_iterations: int = 10000, + ): """ Perform minimisation on the molecules, running a maximum of max_iterations iterations. @@ -918,7 +1008,9 @@ def minimise(self, max_iterations: int = 10000): - max_iterations (int): The maximum number of iterations to run """ if not self._d.is_null(): - self._d.run_minimisation(max_iterations=max_iterations) + self._d.run_minimisation( + max_iterations=max_iterations, + ) return self @@ -929,7 +1021,7 @@ def run( frame_frequency=None, energy_frequency=None, lambda_windows=None, - save_velocities: bool = True, + save_velocities: bool = None, auto_fix_minimise: bool = True, ): """ @@ -989,7 +1081,7 @@ def run( save_velocities: bool Whether or not to save the velocities when running dynamics. - By default this is True. Set this to False if you aren't + By default this is False. Set this to True if you are interested in saving the velocities. auto_fix_minimise: bool @@ -1000,6 +1092,14 @@ def run( in these cases, and then runs the requested dynamics. """ if not self._d.is_null(): + if save_velocities is None: + if self._d._map.specified("save_velocities"): + save_velocities = ( + self._d._map["save_velocities"].value().as_bool() + ) + else: + save_velocities = False + self._d.run( time=time, save_frequency=save_frequency, @@ -1060,6 +1160,44 @@ def set_ensemble(self, ensemble): """ self._d.set_ensemble(ensemble) + def set_temperature(self, temperature, rescale_velocities: bool = True): + """ + Set the temperature for the dynamics. Note that this will only + let you change the temperature of the ensemble. + You can't change its fundemental nature. + + If rescale_velocities is True, then the velocities will be + rescaled to the new temperature. + """ + self._d.set_temperature( + temperature=temperature, rescale_velocities=rescale_velocities + ) + + def set_pressure(self, pressure): + """ + Set the pressure for the dynamics. Note that this will only + let you change the pressure of the ensemble. + You can't change its fundemental nature. + """ + self._d.set_pressure(pressure=pressure) + + def randomise_velocities(self, temperature=None, random_seed: int = None): + """ + Set the velocities to random values, drawn from the Boltzmann + distribution for the current temperature. + + Parameters + ---------- + + - temperature (temperature): The temperature to use. If None, then + the current temperature will be used + - random_seed (int): The random seed to use. If None, then + a random seed will be generated + """ + self._d.randomise_velocities( + temperature=temperature, random_seed=random_seed + ) + def constraint(self): """ Return the constraint used for the dynamics (e.g. constraining @@ -1067,6 +1205,13 @@ def constraint(self): """ return self._d.constraint() + def perturbable_constraint(self): + """ + Return the perturbable constraint used for the dynamics (e.g. + constraining bonds involving hydrogens etc.) + """ + return self._d.perturbable_constraint() + def info(self): """ Return the information that describes the forcefield that will @@ -1204,7 +1349,9 @@ def current_kinetic_energy(self): """ return self._d.current_kinetic_energy() - def energy_trajectory(self, to_pandas: bool = True): + def energy_trajectory( + self, to_pandas: bool = False, to_alchemlyb: bool = False + ): """ Return the energy trajectory. This is the trajectory of energy values that have been captured during dynamics. @@ -1215,8 +1362,10 @@ def energy_trajectory(self, to_pandas: bool = True): """ t = self._d.energy_trajectory() - if to_pandas: - return t.to_pandas() + if to_pandas or to_alchemlyb: + return t.to_pandas( + to_alchemlyb=to_alchemlyb, + ) else: return t diff --git a/src/sire/mol/_minimisation.py b/src/sire/mol/_minimisation.py index 2d7ae2eca..4d5cebaa7 100644 --- a/src/sire/mol/_minimisation.py +++ b/src/sire/mol/_minimisation.py @@ -18,6 +18,7 @@ def __init__( schedule=None, lambda_value=None, swap_end_states=None, + ignore_perturbations=None, shift_delta=None, coulomb_power=None, restraints=None, @@ -33,6 +34,7 @@ def __init__( _add_extra(extras, "schedule", schedule) _add_extra(extras, "lambda", lambda_value) _add_extra(extras, "swap_end_states", swap_end_states) + _add_extra(extras, "ignore_perturbations", ignore_perturbations) _add_extra(extras, "shift_delta", shift_delta) _add_extra(extras, "coulomb_power", coulomb_power) _add_extra(extras, "restraints", restraints) diff --git a/src/sire/morph/CMakeLists.txt b/src/sire/morph/CMakeLists.txt index ea5a9f23e..0cc85b497 100644 --- a/src/sire/morph/CMakeLists.txt +++ b/src/sire/morph/CMakeLists.txt @@ -7,7 +7,9 @@ # Add your script to this list set ( SCRIPTS __init__.py + _alchemy.py _ghost_atoms.py + _hmr.py _perturbation.py _repex.py ) diff --git a/src/sire/morph/__init__.py b/src/sire/morph/__init__.py index f39200ac5..13596f702 100644 --- a/src/sire/morph/__init__.py +++ b/src/sire/morph/__init__.py @@ -1,7 +1,17 @@ -__all__ = ["shrink_ghost_atoms", "replica_exchange", "Perturbation"] +__all__ = [ + "shrink_ghost_atoms", + "replica_exchange", + "repartition_hydrogen_masses", + "to_alchemlyb", + "Perturbation", +] from ._perturbation import * from ._ghost_atoms import * from ._repex import * + +from ._hmr import * + +from ._alchemy import * diff --git a/src/sire/morph/_alchemy.py b/src/sire/morph/_alchemy.py new file mode 100644 index 000000000..f895212ac --- /dev/null +++ b/src/sire/morph/_alchemy.py @@ -0,0 +1,81 @@ +__all__ = ["to_alchemlyb"] + + +def to_alchemlyb(energy_trajectories, temperature=None): + """ + Convert the passed list of energy trajectories into + a single alchemlyb-compatible DataFrame, ready for + free energy analysis via alchemlyb. + + Parameters + ---------- + + energy_trajectories : list of sire.maths.EnergyTrajectory objects, + list of filenames of s3 files, + globbed expression for list of filenames etc. + A list of EnergyTrajectory objects, each containing the + energy trajectory for a single simulation at a single + lambda value. + + temperature : temperature, optional + The temperature of the simulation. If not provided, + the temperature will be taken from the values in + each EnergyTrajectory + + Returns + ------- + + pandas.DataFrame + A single DataFrame containing the energy trajectories + from all simulations, ready for free energy analysis + via alchemlyb. + """ + if type(energy_trajectories) is str: + # this could be a globbed file path + import glob + + energy_trajectories = glob.glob(energy_trajectories) + elif type(energy_trajectories) is not list: + energy_trajectories = [energy_trajectories] + + # convert all of the energy trajectories to pandas DataFrames + dataframes = {} + + for energy_trajectory in energy_trajectories: + if type(energy_trajectory) is str: + # load the energy trajectory from the s3 file + from ..stream import load + + energy_trajectory = load(energy_trajectory) + + df = energy_trajectory.to_alchemlyb(temperature=temperature) + + # get the lambda value of the first row + try: + lambda_value = df.index[0][1] + except IndexError: + # this is a corrupt energy trajectory... + from ..utils import Console + + Console.warning("Corrupt energy trajectory detected. Skipping...") + continue + + if lambda_value in dataframes: + dataframes[lambda_value].append(df) + else: + dataframes[lambda_value] = [df] + + # get the list of lambda values and sort... + lambda_values = list(dataframes.keys()) + lambda_values.sort() + + # now create a list of dataframes in lambda order + dfs = [] + + for lambda_value in lambda_values: + dfs.extend(dataframes[lambda_value]) + + # now concatenate the dataframes + import pandas as pd + + return pd.concat(dfs) diff --git a/src/sire/morph/_hmr.py b/src/sire/morph/_hmr.py new file mode 100644 index 000000000..d51b2a1da --- /dev/null +++ b/src/sire/morph/_hmr.py @@ -0,0 +1,98 @@ +__all__ = ["repartition_hydrogen_masses"] + + +def repartition_hydrogen_masses( + mols, + mass_factor=1.5, + ignore_water: bool = False, + ignore_non_water: bool = False, + include_end_states: bool = True, + map=None, +): + """ + Increase the mass of hydrogen atoms to hmass times * amu, and subtract the + mass increase from the heavy atom the hydrogen is bonded to. + + (based heavily on the repartitionMasses function in + Sire.Tools.OpenMMMD) + + Parameters + ---------- + + mol : sire.mol.Molecule or list of molecules, or System + The molecule(s) whose hydrogen masses should be repartitioned + + mass_factor : float + The factor to multiply the mass of hydrogen atoms by. Using + the default of 1.5 is known to be a good value to use to + achieve a 4 fs timestep with the (default) LangevinMiddle integrator + + ignore_water : bool + Whether or not to ignore water molecules (default False) + + ignore_non_water : bool + Whether or not to ignore non-water molecules (default False) + + include_end_states : bool + Whether or not to repartition the masses of the end states + of perturbable molecules (default True) (i.e. this will + automatically repartition `mass0` and `mass1` in addition + to `mass`) + + map : dict + The property map used to identify molecular properties + + Returns + ------- + + sire.mol.Molecule, list of molecules or System + The repartitioned molecule(s) + """ + + try: + from ..legacy.IO import ( + repartitionHydrogenMass as _repartition_hydrogen_mass, + ) + except Exception: + from ..legacy.IO import ( + repartition_hydrogen_mass as _repartition_hydrogen_mass, + ) + + if ignore_water and ignore_non_water: + # ignoring everything ;-) + return mols + + from ..base import create_map + + map = create_map(map) + + # convert the flag into the integer defined in the BioSimSpace function + if ignore_water: + water = 0 + elif ignore_non_water: + water = 2 + else: + water = 1 + + mols = mols.clone() + + for mol in mols.molecules(): + mol = _repartition_hydrogen_mass(mol, mass_factor, water, map) + + if include_end_states: + mass0 = f"{map['mass'].source()}0" + mass1 = f"{map['mass'].source()}1" + + if mol.has_property(mass0): + map0 = map.clone() + map0.set("mass", mass0) + mol = _repartition_hydrogen_mass(mol, mass_factor, water, map0) + + if mol.has_property(mass1): + map1 = map.clone() + map1.set("mass", mass1) + mol = _repartition_hydrogen_mass(mol, mass_factor, water, map1) + + mols.update(mol) + + return mols diff --git a/src/sire/system/_system.py b/src/sire/system/_system.py index 56836417c..c5b8ace8f 100644 --- a/src/sire/system/_system.py +++ b/src/sire/system/_system.py @@ -374,8 +374,21 @@ def minimisation(self, *args, **kwargs): The type of constraint to use for bonds and/or angles, e.g. `h-bonds`, `bonds` etc. See https://sire.openbiosim.org/cheatsheet/openmm.html#choosing-options - for the full list of options. This will be automatically - guessed from the timestep if it isn't set. + for the full list of options. This will be `none` if it hasn't + been set. + + perturbable_constraint: str + The type of constraint to use for perturbable bonds and/or angles, + e.g. `h-bonds`, `bonds` etc. + See https://sire.openbiosim.org/cheatsheet/openmm.html#choosing-options + for the full list of options. This equal the value of `constraint` + if it isn't set. + + include_constrained_energies: bool + Whether or not to include the energies of the perturbable bonds + and angles. If this is False, then the internal bond or angle + energy of the perturbable degrees of freedom are not included + in the total energy, and their forces are not evaluated. schedule: sire.cas.LambdaSchedule The schedule used to control how perturbable forcefield parameters @@ -396,6 +409,14 @@ def minimisation(self, *args, **kwargs): use the coordinates of the perturbed molecule as the starting point. + ignore_perturbations: bool + Whether or not to ignore perturbations. If this is True, then + the perturbation will be ignored, and the simulation will + be run using the properties of the reference molecule + (or the perturbed molecule if swap_end_states is True). This + is useful if you just want to run standard molecular dynamics + of the reference or perturbed states. + shift_delta: length The shift_delta parameter that controls the electrostatic and van der Waals softening potential that smooths the @@ -407,6 +428,12 @@ def minimisation(self, *args, **kwargs): softening potential that smooths the creation and deletion of ghost atoms during a potential. This defaults to 0. + vacuum: bool + Whether or not to run the simulation in vacuum. If this is + set to `True`, then the simulation space automatically be + replaced by a `sire.vol.Cartesian` space, and the + simulation run in vacuum. + restraints: sire.mm.Restraints or list[sire.mm.Restraints] A single set of restraints, or a list of sets of restraints that will be applied to the atoms during @@ -417,6 +444,10 @@ def minimisation(self, *args, **kwargs): that should be fixed in place during the simulation. These atoms will not be moved by minimisation. + platform: str + The name of the OpenMM platform on which to run the dynamics, + e.g. "CUDA", "OpenCL", "Metal" etc. + device: str or int The ID of the GPU (or accelerator) used to accelerate minimisation. This would be CUDA_DEVICE_ID or similar @@ -466,6 +497,12 @@ def dynamics(self, *args, **kwargs): to save frames during the trajectory. This can be overridden by setting frame_frequency during an individual run. + save_velocities: bool + Whether or not to save velocities when saving trajectory frames + during the simulation. This defaults to False, as velocity + trajectories aren't often needed, and they double the amount + of space that is required for a trajectory. + constraint: str The type of constraint to use for bonds and/or angles, e.g. `h-bonds`, `bonds` etc. @@ -473,6 +510,19 @@ def dynamics(self, *args, **kwargs): for the full list of options. This will be automatically guessed from the timestep if it isn't set. + perturbable_constraint: str + The type of constraint to use for perturbable bonds and/or angles, + e.g. `h-bonds`, `bonds` etc. + See https://sire.openbiosim.org/cheatsheet/openmm.html#choosing-options + for the full list of options. This equal the value of `constraint` + if it isn't set. + + include_constrained_energies: bool + Whether or not to include the energies of the perturbable bonds + and angles. If this is False, then the internal bond or angle + energy of the perturbable degrees of freedom are not included + in the total energy, and their forces are not evaluated. + schedule: sire.cas.LambdaSchedule The schedule used to control how perturbable forcefield parameters should be morphed as a function of lambda. If this is not set @@ -492,6 +542,14 @@ def dynamics(self, *args, **kwargs): use the coordinates of the perturbed molecule as the starting point. + ignore_perturbations: bool + Whether or not to ignore perturbations. If this is True, then + the perturbation will be ignored, and the simulation will + be run using the properties of the reference molecule + (or the perturbed molecule if swap_end_states is True). This + is useful if you just want to run standard molecular dynamics + of the reference or perturbed states. + temperature: temperature The temperature at which to run the simulation. A microcanonical (NVE) simulation will be run if you don't @@ -502,6 +560,12 @@ def dynamics(self, *args, **kwargs): microcanonical (NVE) or canonical (NVT) simulation will be run if the pressure is not set. + vacuum: bool + Whether or not to run the simulation in vacuum. If this is + set to `True`, then the simulation space automatically be + replaced by a `sire.vol.Cartesian` space, and the + simulation run in vacuum. + shift_delta: length The shift_delta parameter that controls the electrostatic and van der Waals softening potential that smooths the @@ -523,11 +587,19 @@ def dynamics(self, *args, **kwargs): that should be fixed in place during the simulation. These atoms will not be moved by dynamics. + platform: str + The name of the OpenMM platform on which to run the dynamics, + e.g. "CUDA", "OpenCL", "Metal" etc. + device: str or int The ID of the GPU (or accelerator) used to accelerate the simulation. This would be CUDA_DEVICE_ID or similar if CUDA was used. This can be any valid OpenMM device string + precision: str + The desired precision for the simulation (e.g. `single`, + `mixed` or `double`) + map: dict A dictionary of additional options. Note that any options set in this dictionary that are also specified via one of @@ -538,6 +610,33 @@ def dynamics(self, *args, **kwargs): return _dynamics(self, *args, **kwargs) + def ensemble(self, map=None): + """ + Return the last ensemble that was used to simulate this system. + This returns a microcanonical ensemble if None was used + """ + from ..base import create_map + + map = create_map(map) + + try: + return self._system.property(map["ensemble"]) + except Exception: + from ..move import Ensemble + + return Ensemble.microcanonical() + + def set_ensemble(self, ensemble, map=None): + """ + Set the ensemble that was last used to simulate this system. + This is copied into the map["ensemble"] property + """ + from ..base import create_map + + map = create_map(map) + + self._system.set_property(map["ensemble"].source(), ensemble) + def energy(self, *args, **kwargs): """Calculate and return the energy of this System (or of the matching index/search subset of this System) @@ -669,11 +768,23 @@ def set_time(self, time, map=None): self._molecules = None - def energy_trajectory(self, to_pandas=True, map=None): + def energy_trajectory( + self, to_pandas: bool = False, to_alchemlyb: bool = False, map=None + ): """ Return the energy trajectory for this System. This is the history of energies evaluate during any dynamics runs. It could include - energies calculated at different values of lambda + energies calculated at different values of lambda. + + Parameters + ---------- + + to_pandas: bool + Whether or not to return the energy trajectory as a pandas DataFrame. + + to_alchemlyb: bool + Whether or not to return the energy trajectory as a pandas DataFrame + that is formatted to usable in alchemlyb """ from ..base import create_map @@ -706,8 +817,20 @@ def energy_trajectory(self, to_pandas=True, map=None): traj = self._system.property(traj_propname) - if to_pandas: - return traj.to_pandas() + if to_pandas or to_alchemlyb: + try: + return traj.to_pandas(to_alchemlyb=to_alchemlyb) + except Exception: + ensemble = self.ensemble() + + if ensemble.is_constant_temperature(): + temperature = ensemble.temperature() + else: + temperature = None + + return traj.to_pandas( + to_alchemlyb=to_alchemlyb, temperature=temperature + ) else: return traj diff --git a/tests/conftest.py b/tests/conftest.py index 62dda5fdd..88f267b8a 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -140,3 +140,13 @@ def merged_ethane_methanol(): @pytest.fixture(scope="session") def merged_zan_ose(): return sr.load_test_files("merged_ligand.s3") + + +@pytest.fixture(scope="session") +def ethane_12dichloroethane(): + return sr.load_test_files("ethane_12dichloroethane.bss") + + +@pytest.fixture(scope="session") +def pentane_cyclopentane(): + return sr.load_test_files("pentane_cyclopentane.bss") diff --git a/tests/convert/test_openmm.py b/tests/convert/test_openmm.py index 7e6001e78..2ea8bf90b 100644 --- a/tests/convert/test_openmm.py +++ b/tests/convert/test_openmm.py @@ -11,7 +11,11 @@ def test_openmm_single_energy(kigaki_mols): mol = mols[0] - map = {"space": sr.vol.Cartesian(), "platform": "Reference"} + map = { + "space": sr.vol.Cartesian(), + "platform": "Reference", + "constraint": "bonds-h-angles", + } omm = sr.convert.to(mol, "openmm", map=map) @@ -35,7 +39,9 @@ def test_openmm_single_energy(kigaki_mols): energy = energy.value_in_unit(energy.unit) # these won't be exactly the same - this is 5227 +/- 0.1 kJ mol-1 - assert mol.energy(map=map).to(sr.units.kJ_per_mol) == pytest.approx(energy, abs=0.1) + assert mol.energy(map=map).to(sr.units.kJ_per_mol) == pytest.approx( + energy, abs=0.1 + ) @pytest.mark.skipif( @@ -46,7 +52,11 @@ def test_openmm_multi_energy_small_cart(kigaki_mols): # first, try just 50 molecules in a cartesian space mols = kigaki_mols[0:50] - map = {"space": sr.vol.Cartesian(), "platform": "Reference"} + map = { + "space": sr.vol.Cartesian(), + "platform": "Reference", + "constraint": "bonds-h-angles", + } omm = sr.convert.to(mols, "openmm", map=map) @@ -78,6 +88,7 @@ def test_openmm_multi_energy_all_cart(kigaki_mols): "cutoff_type": "REACTION_FIELD", "dielectric": 1.0, "platform": "Reference", + "constraint": "bonds-h-angles", } omm = sr.convert.to(mols, "openmm", map=map) @@ -109,6 +120,7 @@ def test_openmm_multi_energy_all_cart_cutoff(kigaki_mols): "cutoff_type": "REACTION_FIELD", "dielectric": 78.0, "platform": "Reference", + "constraint": "bonds-h-angles", } omm = sr.convert.to(mols, "openmm", map=map) @@ -139,6 +151,7 @@ def test_openmm_multi_energy_all_periodic_cutoff(kigaki_mols): "cutoff_type": "REACTION_FIELD", "dielectric": 78.0, "platform": "Reference", + "constraint": "bonds-h-angles", } omm = sr.convert.to(mols, "openmm", map=map) @@ -169,7 +182,8 @@ def test_openmm_dynamics(ala_mols): "cutoff_type": "REACTION_FIELD", "dielectric": 78.0, "temperature": 25 * sr.units.celsius, - "platform": "Reference" + "platform": "Reference", + "constraint": "bonds-h-angles", # "pressure": 1 * sr.units.atm, # currently disagree with energies for NPT... } @@ -228,6 +242,7 @@ def test_openmm_options(ala_mols): "pressure": 1 * sr.units.atm, "friction": 5 / sr.units.picosecond, "platform": "Reference", + "constraint": "bonds-h-angles", } omm = sr.convert.to(mol, "openmm", map=m) @@ -240,3 +255,33 @@ def test_openmm_options(ala_mols): except ValueError: # maybe OpenCL or CUDA are not supported pass + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_openmm_ignore_constrained(ala_mols): + mols = ala_mols + + mol = mols[0] + + d = mol.dynamics( + constraint="bonds-h-angles", + include_constrained_energies=True, + platform="Reference", + ) + + nrg1 = d.current_potential_energy() + + d = mol.dynamics( + constraint="bonds-h-angles", + include_constrained_energies=False, + platform="Reference", + ) + + nrg2 = d.current_potential_energy() + + # these two energies should be different, because + # we should be ignoring the constrained bonds and angles + assert abs(nrg2.value() - nrg1.value()) > 1.0 diff --git a/tests/convert/test_openmm_lambda.py b/tests/convert/test_openmm_lambda.py index bc8e7bb31..1131f9807 100644 --- a/tests/convert/test_openmm_lambda.py +++ b/tests/convert/test_openmm_lambda.py @@ -11,14 +11,17 @@ def _run_test(mols, is_slow=False, use_taylor=False): mols = c.commit() + has_water = len(mols) > 1 + # in this case the perturbable molecule is the first one merge = mols[0] # and everything else is water - if is_slow: - water = mols["water"] - else: - water = mols["closest 5 waters to molidx 0"] + if has_water: + if is_slow: + water = mols["water"] + else: + water = mols["closest 5 waters to molidx 0"] # this function will extract the lambda0 or lambda1 end state def get_end_state(mol, state, remove_state): @@ -34,18 +37,29 @@ def get_end_state(mol, state, remove_state): return c.commit() # this is the perturbable system - mols = merge + water + if has_water: + mols = merge + water - # these are the non-perturbable systems at lambda 0 and lambda 1 - mols0 = get_end_state(merge, "0", "1") + water - mols1 = get_end_state(merge, "1", "0") + water + # these are the non-perturbable systems at lambda 0 and lambda 1 + mols0 = get_end_state(merge, "0", "1") + water + mols1 = get_end_state(merge, "1", "0") + water + else: + mols = merge + + # these are the non-perturbable systems at lambda 0 and lambda 1 + mols0 = get_end_state(merge, "0", "1") + mols1 = get_end_state(merge, "1", "0") # a very basic lambda schedule l = sr.cas.LambdaSchedule() l.add_stage("morph", (1 - l.lam()) * l.initial() + l.lam() * l.final()) # need to use the reference platform on GH Actions - map = {"platform": "Reference", "schedule": l} + map = { + "platform": "Reference", + "schedule": l, + "constraint": "bonds-h-angles", + } if use_taylor: map["use_taylor_softening"] = True @@ -137,3 +151,25 @@ def test_openmm_scale_lambda_ligand(merged_zan_ose): ) def test_big_openmm_scale_lambda_ligand(merged_zan_ose): _run_test(merged_zan_ose.clone(), True) + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_openmm_scale_lambda_dichloroethane(ethane_12dichloroethane): + _run_test(ethane_12dichloroethane.clone(), False) + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_openmm_scale_lambda_cyclopentane(pentane_cyclopentane): + mols = pentane_cyclopentane.clone() + + for mol in mols.molecules("property is_perturbable"): + mol = mol.edit().add_link("connectivity", "connectivity0").commit() + mols.update(mol) + + _run_test(mols, False) diff --git a/tests/mol/test_mapping.py b/tests/mol/test_mapping.py index 705765f1e..f38302c31 100644 --- a/tests/mol/test_mapping.py +++ b/tests/mol/test_mapping.py @@ -1,6 +1,10 @@ import sire as sr +import pytest +import sys + +@pytest.mark.skipif(sys.platform == "win32", reason="Does not run on windows") def test_mapping(ala_mols, ala_traj): mols = ala_mols.clone() mols2 = ala_traj.clone() diff --git a/tests/mol/test_rmsd.py b/tests/mol/test_rmsd.py index be35c4dc7..a359fa304 100644 --- a/tests/mol/test_rmsd.py +++ b/tests/mol/test_rmsd.py @@ -1,7 +1,10 @@ import sire as sr + import pytest +import sys +@pytest.mark.skipif(sys.platform == "win32", reason="Does not run on windows") def test_rmsd(ala_trr, ala_mols): mols = ala_trr.clone() mols2 = ala_mols.clone() diff --git a/tests/morph/test_hmr.py b/tests/morph/test_hmr.py new file mode 100644 index 000000000..3a30f34a0 --- /dev/null +++ b/tests/morph/test_hmr.py @@ -0,0 +1,17 @@ +import pytest +import sire as sr + + +def test_repartition_hydrogen_masses(ala_mols): + mol = ala_mols.molecules()[0] + + newmol = sr.morph.repartition_hydrogen_masses(mol) + + for atom0, atom1 in zip(mol.atoms(), newmol.atoms()): + if atom0.element().num_protons() == 1: + assert atom1.mass() > atom0.mass() + assert atom1.mass() >= sr.u("1.5 g mol-1") + else: + assert atom1.mass() <= atom0.mass() + + assert mol.mass().value() == pytest.approx(newmol.mass().value()) diff --git a/wrapper/CAS/CMakeAutogenFile.txt b/wrapper/CAS/CMakeAutogenFile.txt index 2b2e30254..2a7b25a3c 100644 --- a/wrapper/CAS/CMakeAutogenFile.txt +++ b/wrapper/CAS/CMakeAutogenFile.txt @@ -1,72 +1,72 @@ # WARNING - AUTOGENERATED FILE - CONTENTS WILL BE OVERWRITTEN! set ( PYPP_SOURCES - Abs.pypp.cpp - AlwaysFalse.pypp.cpp + I.pypp.cpp AlwaysTrue.pypp.cpp - ArcCos.pypp.cpp - ArcCosh.pypp.cpp - ArcCot.pypp.cpp - ArcCoth.pypp.cpp - ArcCsc.pypp.cpp - ArcCsch.pypp.cpp - ArcSec.pypp.cpp - ArcSech.pypp.cpp - ArcSin.pypp.cpp + Exp.pypp.cpp + Power.pypp.cpp + _CAS_free_functions.pypp.cpp + ExBase.pypp.cpp + NotEqualTo.pypp.cpp + IntegerPower.pypp.cpp + Csch.pypp.cpp + Tanh.pypp.cpp ArcSinh.pypp.cpp - ArcTan.pypp.cpp - ArcTanh.pypp.cpp - ComplexPower.pypp.cpp - ComplexValues.pypp.cpp - Condition.pypp.cpp + Tan.pypp.cpp + ArcCoth.pypp.cpp + Values.pypp.cpp Conditional.pypp.cpp + Factor.pypp.cpp Constant.pypp.cpp - Cos.pypp.cpp - Cosh.pypp.cpp - Cot.pypp.cpp - Coth.pypp.cpp - Csc.pypp.cpp - Csch.pypp.cpp - DoubleFunc.pypp.cpp - EqualTo.pypp.cpp - ExBase.pypp.cpp - Exp.pypp.cpp - Expression.pypp.cpp ExpressionBase.pypp.cpp - ExpressionProperty.pypp.cpp - Factor.pypp.cpp - GreaterOrEqualThan.pypp.cpp - GreaterThan.pypp.cpp - I.pypp.cpp - Identities.pypp.cpp - IntegerPower.pypp.cpp - IntegrationConstant.pypp.cpp - LambdaSchedule.pypp.cpp LessOrEqualThan.pypp.cpp + ArcCosh.pypp.cpp + RationalPower.pypp.cpp + ComplexPower.pypp.cpp + Sech.pypp.cpp + ArcTan.pypp.cpp LessThan.pypp.cpp + SymbolExpression.pypp.cpp + GreaterOrEqualThan.pypp.cpp + LambdaSchedule.pypp.cpp + SymbolComplex.pypp.cpp + GreaterThan.pypp.cpp + Cosh.pypp.cpp Ln.pypp.cpp - Max.pypp.cpp - Min.pypp.cpp - NotEqualTo.pypp.cpp - Power.pypp.cpp - PowerConstant.pypp.cpp - PowerFunction.pypp.cpp + Cot.pypp.cpp + ArcSech.pypp.cpp Product.pypp.cpp - RationalPower.pypp.cpp RealPower.pypp.cpp Sec.pypp.cpp - Sech.pypp.cpp + DoubleFunc.pypp.cpp + Csc.pypp.cpp + SymbolValue.pypp.cpp + ComplexValues.pypp.cpp + Sum.pypp.cpp + IntegrationConstant.pypp.cpp + EqualTo.pypp.cpp + Min.pypp.cpp + Abs.pypp.cpp + ArcCos.pypp.cpp Sin.pypp.cpp + Symbol.pypp.cpp + Identities.pypp.cpp + ArcSec.pypp.cpp + PowerConstant.pypp.cpp + AlwaysFalse.pypp.cpp + ArcCot.pypp.cpp + Max.pypp.cpp + ExpressionProperty.pypp.cpp + ArcSin.pypp.cpp + Coth.pypp.cpp + Cos.pypp.cpp SingleFunc.pypp.cpp + Condition.pypp.cpp + ArcCsch.pypp.cpp + Expression.pypp.cpp + ArcCsc.pypp.cpp + PowerFunction.pypp.cpp Sinh.pypp.cpp - Sum.pypp.cpp - Symbol.pypp.cpp - SymbolComplex.pypp.cpp - SymbolExpression.pypp.cpp - SymbolValue.pypp.cpp - Tan.pypp.cpp - Tanh.pypp.cpp - Values.pypp.cpp - _CAS_free_functions.pypp.cpp + ArcTanh.pypp.cpp SireCAS_containers.cpp SireCAS_registrars.cpp ) diff --git a/wrapper/CAS/LambdaSchedule.pypp.cpp b/wrapper/CAS/LambdaSchedule.pypp.cpp index bff2a2557..6fd39eb88 100644 --- a/wrapper/CAS/LambdaSchedule.pypp.cpp +++ b/wrapper/CAS/LambdaSchedule.pypp.cpp @@ -45,6 +45,18 @@ void register_LambdaSchedule_class(){ , ( bp::arg("scale")=0.20000000000000001 ) , "Sandwich the current set of stages with a charge-descaling and\n a charge-scaling stage. This prepends a charge-descaling stage\n that scales the charge parameter down from `initial` to\n :gamma:.initial (where :gamma:=`scale`). The charge parameter in all of\n the exising stages in this schedule are then multiplied\n by :gamma:. A final charge-rescaling stage is then appended that\n scales the charge parameter from :gamma:.final to final.\n" ); + } + { //::SireCAS::LambdaSchedule::addChargeScaleStages + + typedef void ( ::SireCAS::LambdaSchedule::*addChargeScaleStages_function_type)( ::QString const &,::QString const &,double ) ; + addChargeScaleStages_function_type addChargeScaleStages_function_value( &::SireCAS::LambdaSchedule::addChargeScaleStages ); + + LambdaSchedule_exposer.def( + "addChargeScaleStages" + , addChargeScaleStages_function_value + , ( bp::arg("decharge_name"), bp::arg("recharge_name"), bp::arg("scale")=0.20000000000000001 ) + , "" ); + } { //::SireCAS::LambdaSchedule::addLever @@ -83,6 +95,19 @@ void register_LambdaSchedule_class(){ , bp::release_gil_policy() , "Append a morph stage onto this schedule. The morph stage is a\n standard stage that scales each forcefield parameter by\n (1-:lambda:).initial + :lambda:.final\n" ); + } + { //::SireCAS::LambdaSchedule::addMorphStage + + typedef void ( ::SireCAS::LambdaSchedule::*addMorphStage_function_type)( ::QString const & ) ; + addMorphStage_function_type addMorphStage_function_value( &::SireCAS::LambdaSchedule::addMorphStage ); + + LambdaSchedule_exposer.def( + "addMorphStage" + , addMorphStage_function_value + , ( bp::arg("name") ) + , bp::release_gil_policy() + , "" ); + } { //::SireCAS::LambdaSchedule::addStage @@ -382,7 +407,7 @@ void register_LambdaSchedule_class(){ , morph_function_value , ( bp::arg("lever"), bp::arg("initial"), bp::arg("final"), bp::arg("lambda_value") ) , bp::release_gil_policy() - , "Return the parameters for the specified lever called `lever_name`\n that have been morphed from the passed list of initial values\n (in `initial`) to the passed list of final values (in `final`)\n for the specified global value of :lambda: (in `lambda_value`).\n\n The morphed parameters will be returned in the matching\n order to `initial` and `final`.\n\n This morphs floating point parameters. There is an overload\n of this function that morphs integer parameters, in which\n case the result would be rounded to the nearest integer.\n" ); + , "Return the parameters for the specified lever called `lever_name`\n that have been morphed from the passed list of initial values\n (in `initial`) to the passed list of final values (in `final`)\n for the specified global value of :lambda: (in `lambda_value`).\n\n The morphed parameters will be returned in the matching\n order to `initial` and `final`.\n\n This morphs a single floating point parameters.\n" ); } { //::SireCAS::LambdaSchedule::morph diff --git a/wrapper/Convert/SireOpenMM/_sommcontext.py b/wrapper/Convert/SireOpenMM/_sommcontext.py index 39c3bd074..5e27ed254 100644 --- a/wrapper/Convert/SireOpenMM/_sommcontext.py +++ b/wrapper/Convert/SireOpenMM/_sommcontext.py @@ -43,6 +43,8 @@ def __init__( # turn this system into a context if map.specified("lambda"): lambda_value = map["lambda"].value().as_double() + elif map.specified("lambda_value"): + lambda_value = map["lambda_value"].value().as_double() else: lambda_value = 0.0 diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 55529e26d..358889699 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -107,13 +107,6 @@ bool LambdaLever::hasLever(const QString &lever_name) void LambdaLever::addLever(const QString &lever_name) { - if (this->hasLever(lever_name)) - throw SireError::invalid_key(QObject::tr( - "You cannot add a new lever called '%1' as a lever with " - "this name already exists.") - .arg(lever_name), - CODELOC); - this->lambda_schedule.addLever(lever_name); } @@ -269,71 +262,75 @@ double LambdaLever::setLambda(OpenMM::Context &context, // calculate the new parameters for this lambda value const auto morphed_charges = this->lambda_schedule.morph( "charge", - perturbable_mol.getCharges(), - perturbable_mol.perturbed->getCharges(), + perturbable_mol.getCharges0(), + perturbable_mol.getCharges1(), lambda_value); const auto morphed_sigmas = this->lambda_schedule.morph( "sigma", - perturbable_mol.getSigmas(), - perturbable_mol.perturbed->getSigmas(), + perturbable_mol.getSigmas0(), + perturbable_mol.getSigmas1(), lambda_value); const auto morphed_epsilons = this->lambda_schedule.morph( "epsilon", - perturbable_mol.getEpsilons(), - perturbable_mol.perturbed->getEpsilons(), + perturbable_mol.getEpsilons0(), + perturbable_mol.getEpsilons1(), lambda_value); const auto morphed_alphas = this->lambda_schedule.morph( "alpha", - perturbable_mol.getAlphas(), - perturbable_mol.perturbed->getAlphas(), + perturbable_mol.getAlphas0(), + perturbable_mol.getAlphas1(), lambda_value); const auto morphed_bond_k = this->lambda_schedule.morph( "bond_k", - perturbable_mol.getBondKs(), - perturbable_mol.perturbed->getBondKs(), + perturbable_mol.getBondKs0(), + perturbable_mol.getBondKs1(), lambda_value); const auto morphed_bond_length = this->lambda_schedule.morph( "bond_length", - perturbable_mol.getBondLengths(), - perturbable_mol.perturbed->getBondLengths(), + perturbable_mol.getBondLengths0(), + perturbable_mol.getBondLengths1(), lambda_value); const auto morphed_angle_k = this->lambda_schedule.morph( "angle_k", - perturbable_mol.getAngleKs(), - perturbable_mol.perturbed->getAngleKs(), + perturbable_mol.getAngleKs0(), + perturbable_mol.getAngleKs1(), lambda_value); const auto morphed_angle_size = this->lambda_schedule.morph( "angle_size", - perturbable_mol.getAngleSizes(), - perturbable_mol.perturbed->getAngleSizes(), + perturbable_mol.getAngleSizes0(), + perturbable_mol.getAngleSizes1(), lambda_value); const auto morphed_torsion_phase = this->lambda_schedule.morph( "torsion_phase", - perturbable_mol.getTorsionPhases(), - perturbable_mol.perturbed->getTorsionPhases(), - lambda_value); - - const auto morphed_torsion_periodicity = this->lambda_schedule.morph( - "torsion_periodicity", - perturbable_mol.getTorsionPeriodicities(), - perturbable_mol.perturbed->getTorsionPeriodicities(), + perturbable_mol.getTorsionPhases0(), + perturbable_mol.getTorsionPhases1(), lambda_value); const auto morphed_torsion_k = this->lambda_schedule.morph( "torsion_k", - perturbable_mol.getTorsionKs(), - perturbable_mol.perturbed->getTorsionKs(), + perturbable_mol.getTorsionKs0(), + perturbable_mol.getTorsionKs1(), lambda_value); - // will eventually morph the NB14 / exception parameters :-) + const auto morphed_charge_scale = this->lambda_schedule.morph( + "charge_scale", + perturbable_mol.getChargeScales0(), + perturbable_mol.getChargeScales1(), + lambda_value); + + const auto morphed_lj_scale = this->lambda_schedule.morph( + "lj_scale", + perturbable_mol.getLJScales0(), + perturbable_mol.getLJScales1(), + lambda_value); // now update the forcefields int start_index = start_idxs.value("clj", -1); @@ -346,8 +343,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, { for (int j = 0; j < nparams; ++j) { - const bool is_from_ghost = perturbable_mol.from_ghost_idxs.contains(j); - const bool is_to_ghost = perturbable_mol.to_ghost_idxs.contains(j); + const bool is_from_ghost = perturbable_mol.getFromGhostIdxs().contains(j); + const bool is_to_ghost = perturbable_mol.getToGhostIdxs().contains(j); // reduced charge custom_params[0] = morphed_charges[j]; @@ -386,24 +383,26 @@ double LambdaLever::setLambda(OpenMM::Context &context, } } - const auto idxs = perturbable_mol.exception_idxs.value("clj"); + const auto idxs = perturbable_mol.getExceptionIndicies("clj"); if (not idxs.isEmpty()) { - for (int j = 0; j < perturbable_mol.exception_params.count(); ++j) + const auto exception_atoms = perturbable_mol.getExceptionAtoms(); + + for (int j = 0; j < exception_atoms.count(); ++j) { - const auto ¶m = perturbable_mol.exception_params[j]; + const auto &atoms = exception_atoms[j]; - const auto atom0 = std::get<0>(param); - const auto atom1 = std::get<1>(param); + const auto atom0 = std::get<0>(atoms); + const auto atom1 = std::get<1>(atoms); - const auto coul_14_scale = std::get<2>(param); - const auto lj_14_scale = std::get<3>(param); + const auto coul_14_scale = morphed_charge_scale[j]; + const auto lj_14_scale = morphed_lj_scale[j]; const bool atom0_is_ghost = perturbable_mol.isGhostAtom(atom0); const bool atom1_is_ghost = perturbable_mol.isGhostAtom(atom1); - const auto p = get_exception(std::get<0>(param), std::get<1>(param), + const auto p = get_exception(atom0, atom1, start_index, coul_14_scale, lj_14_scale, morphed_charges, morphed_sigmas, morphed_epsilons, morphed_alphas); @@ -521,7 +520,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, dihff->setTorsionParameters(index, particle1, particle2, particle3, particle4, - morphed_torsion_periodicity[j], + periodicity, morphed_torsion_phase[j], morphed_torsion_k[j]); } @@ -762,7 +761,7 @@ int LambdaLever::addPerturbableMolecule(const OpenMMMolecule &molecule, const QHash &starts) { // should add in some sanity checks for these inputs - this->perturbable_mols.append(molecule); + this->perturbable_mols.append(PerturbableOpenMMMolecule(molecule)); this->start_indicies.append(starts); this->perturbable_maps.insert(molecule.number, molecule.perturtable_map); return this->perturbable_mols.count() - 1; @@ -777,8 +776,7 @@ void LambdaLever::setExceptionIndicies(int mol_idx, { mol_idx = SireID::Index(mol_idx).map(this->perturbable_mols.count()); - this->perturbable_mols[mol_idx].exception_idxs.insert( - name, exception_idxs); + this->perturbable_mols[mol_idx].setExceptionIndicies(name, exception_idxs); } /** Return all of the property maps used to find the perturbable properties @@ -796,5 +794,12 @@ LambdaSchedule LambdaLever::getSchedule() const void LambdaLever::setSchedule(const LambdaSchedule &sched) { + QStringList levers = lambda_schedule.getLevers(); + lambda_schedule = sched; + + for (const auto &lever : levers) + { + lambda_schedule.addLever(lever); + } } diff --git a/wrapper/Convert/SireOpenMM/lambdalever.h b/wrapper/Convert/SireOpenMM/lambdalever.h index ecb4a4de1..ee158ce25 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.h +++ b/wrapper/Convert/SireOpenMM/lambdalever.h @@ -109,7 +109,7 @@ namespace SireOpenMM SireCAS::LambdaSchedule lambda_schedule; /** The list of perturbable molecules */ - QVector perturbable_mols; + QVector perturbable_mols; /** The start indicies of the parameters in each named forcefield for each perturbable moleucle */ diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 9f3d02cb8..23cd19bc9 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -38,6 +38,10 @@ using namespace SireMM; using namespace SireMol; using namespace SireBase; +//////// +//////// Implementation of OpenMMMolecule +//////// + OpenMMMolecule::OpenMMMolecule() { } @@ -95,11 +99,11 @@ OpenMMMolecule::OpenMMMolecule(const Molecule &mol, } else if (c == "h-bonds-h-angles") { - constraint_type = CONSTRAIN_HBONDS | CONSTRAIN_ANGLES; + constraint_type = CONSTRAIN_HBONDS | CONSTRAIN_HANGLES; } else if (c == "bonds-h-angles") { - constraint_type = CONSTRAIN_BONDS | CONSTRAIN_ANGLES; + constraint_type = CONSTRAIN_BONDS | CONSTRAIN_HANGLES; } else { @@ -116,13 +120,56 @@ OpenMMMolecule::OpenMMMolecule(const Molecule &mol, constraint_type = CONSTRAIN_NONE; } + if (map.specified("perturbable_constraint")) + { + const auto c = map["perturbable_constraint"].source().toLower().simplified(); + + if (c == "none") + { + perturbable_constraint_type = CONSTRAIN_NONE; + } + else if (c == "h-bonds") + { + perturbable_constraint_type = CONSTRAIN_HBONDS; + } + else if (c == "bonds") + { + perturbable_constraint_type = CONSTRAIN_BONDS; + } + else if (c == "h-bonds-h-angles") + { + perturbable_constraint_type = CONSTRAIN_HBONDS | CONSTRAIN_HANGLES; + } + else if (c == "bonds-h-angles") + { + perturbable_constraint_type = CONSTRAIN_BONDS | CONSTRAIN_HANGLES; + } + else + { + throw SireError::invalid_key(QObject::tr( + "Unrecognised perturbable constraint type '%1'. Valid values are " + "'none', 'h-bonds', 'bonds', 'h-bonds-h-angles' or " + "'bonds-h-angles',") + .arg(c), + CODELOC); + } + } + else + { + perturbable_constraint_type = constraint_type; + } + if (ffinfo.isAmberStyle()) { if (is_perturbable) { // update the map to find the lambda=0 properties + // Note that we don't use the coordinates0 or coordinates1 + // properties, because we need to build the molecule from + // its current coordinates (which should represent the + // current lambda state) QStringList props = {"LJ", "ambertype", "angle", "atomtype", - "bond", "charge", "coordinates", + "bond", "charge", "dihedral", "element", "forcefield", "gb_radii", "gb_screening", "improper", "intrascale", "mass", "name", @@ -289,18 +336,14 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, } // extract the coordinates and convert to OpenMM units - const auto coords_prop = map["coordinates"]; + const auto &c = moldata.property(map["coordinates"]).asA(); + this->coords = QVector(nats, OpenMM::Vec3(0, 0, 0)); + auto coords_data = coords.data(); - if (moldata.hasProperty(coords_prop)) + for (int i = 0; i < nats; ++i) { - const auto &c = moldata.property(coords_prop).asA(); - auto coords_data = coords.data(); - - for (int i = 0; i < nats; ++i) - { - coords_data[i] = to_vec3(c.at(idx_to_cgatomidx_data[i])); - } + coords_data[i] = to_vec3(c.at(idx_to_cgatomidx_data[i])); } // extract the velocities and convert to OpenMM units @@ -398,7 +441,6 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, cljs_data[i] = std::make_tuple(chg, sig, eps); } - this->bond_pairs.clear(); this->bond_params.clear(); this->constraints.clear(); @@ -408,13 +450,19 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, QSet constrained_pairs; + bool include_constrained_energies = true; + + if (map.specified("include_constrained_energies")) + { + include_constrained_energies = map["include_constrained_energies"].value().asABoolean(); + } + for (auto it = params.bonds().constBegin(); it != params.bonds().constEnd(); ++it) { const auto bondid = it.key().map(molinfo); const auto &bondparam = it.value().first; - const auto includes_h = it.value().second; int atom0 = bondid.get<0>().value(); int atom1 = bondid.get<1>().value(); @@ -425,20 +473,33 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, const double k = bondparam.k() * bond_k_to_openmm; const double r0 = bondparam.r0() * bond_r0_to_openmm; - const bool has_light_atom = includes_h or (masses_data[atom0] < 2.5 or masses_data[atom1] < 2.5); + const bool has_light_atom = (masses_data[atom0] < 2.5 or masses_data[atom1] < 2.5); const bool has_massless_atom = masses_data[atom0] < 0.5 or masses_data[atom1] < 0.5; - this->bond_pairs.append(std::make_pair(atom0, atom1)); + auto this_constraint_type = constraint_type; - if ((not has_massless_atom) and ((constraint_type & CONSTRAIN_BONDS) or (has_light_atom and (constraint_type & CONSTRAIN_HBONDS)))) + if (is_perturbable) { - this->constraints.append(std::make_tuple(atom0, atom1, r0)); - constrained_pairs.insert(to_pair(atom0, atom1)); + this_constraint_type = perturbable_constraint_type; } - else + + bool bond_is_not_constrained = true; + + if ((not has_massless_atom) and ((this_constraint_type & CONSTRAIN_BONDS) or (has_light_atom and (this_constraint_type & CONSTRAIN_HBONDS)))) { - this->bond_params.append(std::make_tuple(atom0, atom1, r0, k)); + // add the constraint - this constrains the bond to whatever length it has now + const auto delta = coords[atom1] - coords[atom0]; + const auto constraint_length = std::sqrt((delta[0] * delta[0]) + + (delta[1] * delta[1]) + + (delta[2] * delta[2])); + + this->constraints.append(std::make_tuple(atom0, atom1, constraint_length)); + constrained_pairs.insert(to_pair(atom0, atom1)); + bond_is_not_constrained = false; } + + if (include_constrained_energies or bond_is_not_constrained) + this->bond_params.append(std::make_tuple(atom0, atom1, r0, k)); } // now the angles @@ -467,24 +528,36 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, const auto key = to_pair(atom0, atom2); + bool angle_is_not_constrained = true; + if (not constrained_pairs.contains(key)) { + auto this_constraint_type = constraint_type; + + if (is_perturbable) + { + this_constraint_type = perturbable_constraint_type; + } + // only include the angle X-y-Z if X-Z are not already constrained - if ((constraint_type & CONSTRAIN_ANGLES) and is_h_x_h) + if ((this_constraint_type & CONSTRAIN_HANGLES) and is_h_x_h) { const auto delta = coords[atom2] - coords[atom0]; - const auto length = std::sqrt((delta[0] * delta[0]) + - (delta[1] * delta[1]) + - (delta[2] * delta[2])); - constraints.append(std::make_tuple(atom0, atom2, length)); + const auto constraint_length = std::sqrt((delta[0] * delta[0]) + + (delta[1] * delta[1]) + + (delta[2] * delta[2])); + constraints.append(std::make_tuple(atom0, atom2, + constraint_length)); constrained_pairs.insert(key); - } - else - { - ang_params.append(std::make_tuple(atom0, atom1, atom2, - theta0, k)); + angle_is_not_constrained = false; } } + else + angle_is_not_constrained = false; + + if (include_constrained_energies or angle_is_not_constrained) + ang_params.append(std::make_tuple(atom0, atom1, atom2, + theta0, k)); } // now the dihedrals @@ -498,7 +571,6 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, { const auto dihid = it.key().map(molinfo); const auto &dihparam = it.value().first; - const auto includes_h = it.value().second; int atom0 = dihid.get<0>().value(); int atom1 = dihid.get<1>().value(); @@ -545,7 +617,6 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, { const auto impid = it.key().map(molinfo); const auto &impparam = it.value().first; - const auto includes_h = it.value().second; const int atom0 = impid.get<0>().value(); const int atom1 = impid.get<1>().value(); @@ -579,22 +650,7 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, } } - this->buildPerturbableExceptions(mol, constrained_pairs, map); - - /*if (is_perturbable) - { - // build all of the exceptions for perturbable molecules. - // This is because the exceptions could change during the - // perturbation, so we need more control over them - this->buildPerturbableExceptions(mol, constrained_pairs, map); - } - else - { - // just find any exceptions that would be missing based on just - // a cursary look at the bonding - this->buildStandardExceptions(mol, params, - constrained_pairs, map); - }*/ + this->buildExceptions(mol, constrained_pairs, map); } bool is_ghost(const std::tuple &clj) @@ -616,10 +672,6 @@ bool is_ghost_lj(const std::tuple &clj) * state. Ensure that there is a one-to-one mapping, with them all * in the same order. Any that are missing are added as nulls in * the correct end state. - * - * Note that bonds and angles involved in constraints have been - * removed from the list of potentials. This is because we don't - * evaluate the energy of constrained degrees of freedom */ void OpenMMMolecule::alignInternals(const PropertyMap &map) { @@ -872,82 +924,90 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) perturbed->dih_params = dih_params_1; - // we may need to align the the standard_14_pairs - // and the custom_14_pairs - this is the work that will be needed - // if the bonding changes across the perturbation! -} + // now align all of the exceptions - this should allow the bonding + // to change during the perturbation + QVector> exception_params_1; + exception_params_1.reserve(exception_params.count()); -void OpenMMMolecule::buildStandardExceptions(const Molecule &mol, - const AmberParams ¶ms, - QSet &constrained_pairs, - const PropertyMap &map) -{ - exception_params.clear(); + found_index_0 = QVector(exception_params.count(), false); + found_index_1 = QVector(perturbed->exception_params.count(), false); - // add all of the 1-4 terms - for (auto it = params.nb14s().constBegin(); - it != params.nb14s().constEnd(); - ++it) + for (int i = 0; i < exception_params.count(); ++i) { - const double cscl = it.value().cscl(); - const double ljscl = it.value().ljscl(); + const auto &ex0 = exception_params.at(i); - const auto nbid = it.key().map(molinfo); + int atom0 = std::get<0>(ex0); + int atom1 = std::get<1>(ex0); - const int atom0 = nbid.get<0>().value(); - const int atom1 = nbid.get<1>().value(); - - exception_params.append(std::make_tuple( - atom0, atom1, cscl, ljscl)); - } - - // add in all of the excluded atoms - const auto excl_pairs = ExcludedPairs(mol, map); + bool found = false; - if (excl_pairs.count() > 0) - { - for (int i = 0; i < excl_pairs.count(); ++i) + for (int j = 0; j < perturbed->exception_params.count(); ++j) { - auto pair = excl_pairs[i]; + const auto &ex1 = perturbed->exception_params.at(j); - exception_params.append(std::make_tuple( - std::get<0>(pair).value(), std::get<1>(pair).value(), - 0.0, 0.0)); + if (std::get<0>(ex1) == atom0 and std::get<1>(ex1) == atom1) + { + // we have found the matching exception! + exception_params_1.append(ex1); + found_index_0[i] = true; + found_index_1[j] = true; + found = true; + break; + } } - // and finally (finally!) find any atoms that are not bonded to - // anything else and make sure that they are constrained. These - // atoms will be excluded atoms (by definition) so just look - // through those - const auto &connectivity = mol.property(map["connectivity"]).asA(); + if (not found) + { + // add a null exception with the same scale factors + exception_params_1.append(std::tuple(atom0, atom1, 1.0, 1.0)); + found_index_0[i] = true; + } + } - for (int i = 0; i < excl_pairs.count(); ++i) + for (int j = 0; j < perturbed->exception_params.count(); ++j) + { + if (not found_index_1[j]) { - auto pair = excl_pairs[i]; + // need to add an exception missing in the reference state + const auto &ex1 = perturbed->exception_params.at(j); - const int atom0 = std::get<0>(pair).value(); - const int atom1 = std::get<1>(pair).value(); + int atom0 = std::get<0>(ex1); + int atom1 = std::get<1>(ex1); - if (not constrained_pairs.contains(to_pair(atom0, atom1))) - { - if (connectivity.nConnections(std::get<0>(pair)) == 0 or - connectivity.nConnections(std::get<1>(pair)) == 0) - { - const auto delta = coords[atom1] - coords[atom0]; - const auto length = std::sqrt((delta[0] * delta[0]) + - (delta[1] * delta[1]) + - (delta[2] * delta[2])); - constraints.append(std::make_tuple(atom0, atom1, length)); - constrained_pairs.insert(to_pair(atom0, atom1)); - } - } + // add a null exception + exception_params.append(std::tuple(atom0, atom1, 1.0, 1.0)); + exception_params_1.append(ex1); + found_index_1[j] = true; } } + + // all of found_index_0 and found_index_1 should be true... + if (found_index_0.indexOf(false) != -1 or found_index_1.indexOf(false) != -1) + { + throw SireError::program_bug(QObject::tr( + "Failed to align the exceptions!"), + CODELOC); + } + + perturbed->exception_params = exception_params_1; + + if (exception_params.count() != perturbed->exception_params.count()) + { + throw SireError::program_bug(QObject::tr( + "Different number of exceptions between the reference " + "(%1) and perturbed (%2) states.") + .arg(exception_params.count()) + .arg(perturbed->exception_params.count()), + CODELOC); + } } -void OpenMMMolecule::buildPerturbableExceptions(const Molecule &mol, - QSet &constrained_pairs, - const PropertyMap &map) +/** Internal function that builds all of the exceptions for all of the + atoms in the molecule +*/ +void OpenMMMolecule::buildExceptions(const Molecule &mol, + QSet &constrained_pairs, + const PropertyMap &map) { // we will build the complete exception list, and will not rely // on this list being built by an OpenMM forcefield. This is because @@ -1183,11 +1243,15 @@ void OpenMMMolecule::copyInCoordsAndVelocities(OpenMM::Vec3 *c, OpenMM::Vec3 *v) } } +/** Return the alpha parameters of all atoms in atom order for + * this molecule + */ QVector OpenMMMolecule::getAlphas() const { return this->alphas; } +/** Return all of the atom charges in atom order for this molecule */ QVector OpenMMMolecule::getCharges() const { const int natoms = this->cljs.count(); @@ -1205,6 +1269,7 @@ QVector OpenMMMolecule::getCharges() const return charges; } +/** Return all of the LJ sigma parameters in atom order for this molecule */ QVector OpenMMMolecule::getSigmas() const { const int natoms = this->cljs.count(); @@ -1222,6 +1287,7 @@ QVector OpenMMMolecule::getSigmas() const return sigmas; } +/** Return all of the LJ epsilon parameters in atom order for this molecule */ QVector OpenMMMolecule::getEpsilons() const { const int natoms = this->cljs.count(); @@ -1239,6 +1305,7 @@ QVector OpenMMMolecule::getEpsilons() const return epsilons; } +/** Return all of the bond k parameters in bond order for this molecule */ QVector OpenMMMolecule::getBondKs() const { const int nbonds = this->bond_params.count(); @@ -1256,6 +1323,7 @@ QVector OpenMMMolecule::getBondKs() const return bond_ks; } +/** Return all of the bond length parameters in bond order for this molecule */ QVector OpenMMMolecule::getBondLengths() const { const int nbonds = this->bond_params.count(); @@ -1273,6 +1341,7 @@ QVector OpenMMMolecule::getBondLengths() const return bond_lengths; } +/** Return all of the angle K parameters in angle order for this molecule */ QVector OpenMMMolecule::getAngleKs() const { const int nangs = this->ang_params.count(); @@ -1290,6 +1359,7 @@ QVector OpenMMMolecule::getAngleKs() const return ang_ks; } +/** Return all of the angle size parameters in angle order for this molecule */ QVector OpenMMMolecule::getAngleSizes() const { const int nangs = this->ang_params.count(); @@ -1307,6 +1377,9 @@ QVector OpenMMMolecule::getAngleSizes() const return ang_sizes; } +/** Return all of the dihedral torsion periodicities in dihedral order + * for this molecule + */ QVector OpenMMMolecule::getTorsionPeriodicities() const { const int ndihs = this->dih_params.count(); @@ -1324,6 +1397,9 @@ QVector OpenMMMolecule::getTorsionPeriodicities() const return dih_periodicities; } +/** Return all of the torsion phase parameters for all of the dihedrals + * in dihedral order for this molecule + */ QVector OpenMMMolecule::getTorsionPhases() const { const int ndihs = this->dih_params.count(); @@ -1341,6 +1417,9 @@ QVector OpenMMMolecule::getTorsionPhases() const return dih_phases; } +/** Return all of the torsion K values for all of the dihedrals in + * dihedral order for this molecule + */ QVector OpenMMMolecule::getTorsionKs() const { const int ndihs = this->dih_params.count(); @@ -1357,3 +1436,379 @@ QVector OpenMMMolecule::getTorsionKs() const return dih_ks; } + +/** Return the atom indexes of the atoms in the exceptions, in + * exception order for this molecule + */ +QVector> OpenMMMolecule::getExceptionAtoms() const +{ + const int nexceptions = this->exception_params.count(); + + QVector> exception_atoms(nexceptions); + + auto exception_atoms_data = exception_atoms.data(); + const auto exception_params_data = this->exception_params.constData(); + + for (int i = 0; i < nexceptions; ++i) + { + exception_atoms_data[i] = std::make_pair(std::get<0>(exception_params_data[i]), + std::get<1>(exception_params_data[i])); + } + + return exception_atoms; +} + +/** Return all of the coulomb intramolecular scale factors (nbscl) + * in exception order for this molecule + */ +QVector OpenMMMolecule::getChargeScales() const +{ + const int nexceptions = this->exception_params.count(); + + QVector charge_scales(nexceptions); + + auto charge_scales_data = charge_scales.data(); + const auto exception_params_data = this->exception_params.constData(); + + for (int i = 0; i < nexceptions; ++i) + { + charge_scales_data[i] = std::get<2>(exception_params_data[i]); + } + + return charge_scales; +} + +/** Return all of the LJ intramolecular scale factors (nbscl) + * in exception order for this molecule + */ +QVector OpenMMMolecule::getLJScales() const +{ + const int nexceptions = this->exception_params.count(); + + QVector lj_scales(nexceptions); + + auto lj_scales_data = lj_scales.data(); + const auto exception_params_data = this->exception_params.constData(); + + for (int i = 0; i < nexceptions; ++i) + { + lj_scales_data[i] = std::get<3>(exception_params_data[i]); + } + + return lj_scales; +} + +//////// +//////// Implementation of PerturbableOpenMMMolecule +//////// + +/** Null constructor */ +PerturbableOpenMMMolecule::PerturbableOpenMMMolecule() +{ +} + +/** Construct from the passed OpenMMMolecule */ +PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const OpenMMMolecule &mol) +{ + if (mol.perturbed.get() == 0) + throw SireError::incompatible_error(QObject::tr( + "You cannot construct a PerturbableOpenMMMolecule from an " + "OpenMMMolecule that has not been perturbed!"), + CODELOC); + + alpha0 = mol.getAlphas(); + alpha1 = mol.perturbed->getAlphas(); + + chg0 = mol.getCharges(); + chg1 = mol.perturbed->getCharges(); + + sig0 = mol.getSigmas(); + sig1 = mol.perturbed->getSigmas(); + + eps0 = mol.getEpsilons(); + eps1 = mol.perturbed->getEpsilons(); + + bond_k0 = mol.getBondKs(); + bond_k1 = mol.perturbed->getBondKs(); + + bond_r0 = mol.getBondLengths(); + bond_r1 = mol.perturbed->getBondLengths(); + + ang_k0 = mol.getAngleKs(); + ang_k1 = mol.perturbed->getAngleKs(); + + ang_t0 = mol.getAngleSizes(); + ang_t1 = mol.perturbed->getAngleSizes(); + + tors_k0 = mol.getTorsionKs(); + tors_k1 = mol.perturbed->getTorsionKs(); + + tors_periodicity0 = mol.getTorsionPeriodicities(); + tors_periodicity1 = mol.perturbed->getTorsionPeriodicities(); + + tors_phase0 = mol.getTorsionPhases(); + tors_phase1 = mol.perturbed->getTorsionPhases(); + + charge_scl0 = mol.getChargeScales(); + charge_scl1 = mol.perturbed->getChargeScales(); + + exception_atoms = mol.getExceptionAtoms(); + + lj_scl0 = mol.getLJScales(); + lj_scl1 = mol.perturbed->getLJScales(); + + to_ghost_idxs = mol.to_ghost_idxs; + from_ghost_idxs = mol.from_ghost_idxs; +} + +/** Copy constructor */ +PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const PerturbableOpenMMMolecule &other) + : alpha0(other.alpha0), alpha1(other.alpha1), + chg0(other.chg0), chg1(other.chg1), + sig0(other.sig0), sig1(other.sig1), + eps0(other.eps0), eps1(other.eps1), + bond_k0(other.bond_k0), bond_k1(other.bond_k1), + bond_r0(other.bond_r0), bond_r1(other.bond_r1), + ang_k0(other.ang_k0), ang_k1(other.ang_k1), + ang_t0(other.ang_t0), ang_t1(other.ang_t1), + tors_k0(other.tors_k0), tors_k1(other.tors_k1), + tors_periodicity0(other.tors_periodicity0), + tors_periodicity1(other.tors_periodicity1), + tors_phase0(other.tors_phase0), tors_phase1(other.tors_phase1), + charge_scl0(other.charge_scl0), charge_scl1(other.charge_scl1), + lj_scl0(other.lj_scl0), lj_scl1(other.lj_scl1), + to_ghost_idxs(other.to_ghost_idxs), from_ghost_idxs(other.from_ghost_idxs), + exception_atoms(other.exception_atoms), exception_idxs(other.exception_idxs) +{ +} + +/** Destructor */ +PerturbableOpenMMMolecule::~PerturbableOpenMMMolecule() +{ +} + +/** Comparison operator */ +bool PerturbableOpenMMMolecule::operator==(const PerturbableOpenMMMolecule &other) const +{ + return alpha0 == alpha1 and chg0 == chg1 and sig0 == sig1 and eps0 == eps1 and + bond_k0 == bond_k1 and bond_r0 == bond_r1 and ang_k0 == ang_k1 and + ang_t0 == ang_t1 and tors_k0 == tors_k1 and tors_periodicity0 == tors_periodicity1 and + tors_phase0 == tors_phase1 and charge_scl0 == charge_scl1 and lj_scl0 == lj_scl1 and + to_ghost_idxs == other.to_ghost_idxs and from_ghost_idxs == other.from_ghost_idxs and + exception_atoms == other.exception_atoms and exception_idxs == other.exception_idxs; +} + +/** Comparison operator */ +bool PerturbableOpenMMMolecule::operator!=(const PerturbableOpenMMMolecule &other) const +{ + return not this->operator==(other); +} + +/** Return the alpha parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getAlphas0() const +{ + return this->alpha0; +} + +/** Return the alpha parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getAlphas1() const +{ + return this->alpha1; +} + +/** Return the atom charges of the reference state */ +QVector PerturbableOpenMMMolecule::getCharges0() const +{ + return this->chg0; +} + +/** Return the atom charges of the perturbed state */ +QVector PerturbableOpenMMMolecule::getCharges1() const +{ + return this->chg1; +} + +/** Return the LJ sigma parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getSigmas0() const +{ + return this->sig0; +} + +/** Return the LJ sigma parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getSigmas1() const +{ + return this->sig1; +} + +/** Return the LJ epsilon parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getEpsilons0() const +{ + return this->eps0; +} + +/** Return the LJ epsilon parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getEpsilons1() const +{ + return this->eps1; +} + +/** Return the bond k parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getBondKs0() const +{ + return this->bond_k0; +} + +/** Return the bond k parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getBondKs1() const +{ + return this->bond_k1; +} + +/** Return the bond length parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getBondLengths0() const +{ + return this->bond_r0; +} + +/** Return the bond length parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getBondLengths1() const +{ + return this->bond_r1; +} + +/** Return the angle k parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getAngleKs0() const +{ + return this->ang_k0; +} + +/** Return the angle k parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getAngleKs1() const +{ + return this->ang_k1; +} + +/** Return the angle size parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getAngleSizes0() const +{ + return this->ang_t0; +} + +/** Return the angle size parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getAngleSizes1() const +{ + return this->ang_t1; +} + +/** Return the torsion k parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getTorsionKs0() const +{ + return this->tors_k0; +} + +/** Return the torsion k parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getTorsionKs1() const +{ + return this->tors_k1; +} + +/** Return the torsion periodicity parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getTorsionPeriodicities0() const +{ + return this->tors_periodicity0; +} + +/** Return the torsion periodicity parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getTorsionPeriodicities1() const +{ + return this->tors_periodicity1; +} + +/** Return the torsion phase parameters of the reference state */ +QVector PerturbableOpenMMMolecule::getTorsionPhases0() const +{ + return this->tors_phase0; +} + +/** Return the torsion phase parameters of the perturbed state */ +QVector PerturbableOpenMMMolecule::getTorsionPhases1() const +{ + return this->tors_phase1; +} + +/** Return the coulomb intramolecular scale factors of the reference state */ +QVector PerturbableOpenMMMolecule::getChargeScales0() const +{ + return this->charge_scl0; +} + +/** Return the coulomb intramolecular scale factors of the perturbed state */ +QVector PerturbableOpenMMMolecule::getChargeScales1() const +{ + return this->charge_scl1; +} + +/** Return the LJ intramolecular scale factors of the reference state */ +QVector PerturbableOpenMMMolecule::getLJScales0() const +{ + return this->lj_scl0; +} + +/** Return the LJ intramolecular scale factors of the perturbed state */ +QVector PerturbableOpenMMMolecule::getLJScales1() const +{ + return this->lj_scl1; +} + +/** Return the indexes of the atoms that are to be ghosted in the + * perturbed state + */ +QSet PerturbableOpenMMMolecule::getToGhostIdxs() const +{ + return this->to_ghost_idxs; +} + +/** Return the indexes of the atoms that were ghosts in the + * reference state + */ +QSet PerturbableOpenMMMolecule::getFromGhostIdxs() const +{ + return this->from_ghost_idxs; +} + +/** Return true if the atom is a ghost atom in the + * referenece or perturbed states */ +bool PerturbableOpenMMMolecule::isGhostAtom(int atom) const +{ + return from_ghost_idxs.contains(atom) or to_ghost_idxs.contains(atom); +} + +/** Return the indices of the atoms in the exceptions */ +QVector> PerturbableOpenMMMolecule::getExceptionAtoms() const +{ + return this->exception_atoms; +} + +/** Return the global indexes of the exceptions in the non-bonded and + * ghost-14 forces + */ +QVector> PerturbableOpenMMMolecule::getExceptionIndicies(const QString &name) const +{ + return this->exception_idxs.value(name); +} + +/** Set the global indexes of the exceptions in the non-bonded and + * ghost-14 forces + */ +void PerturbableOpenMMMolecule::setExceptionIndicies(const QString &name, + const QVector> &exception_idxs) +{ + if (exception_idxs.count() != this->exception_atoms.count()) + throw SireError::incompatible_error(QObject::tr( + "The number of exception indicies (%1) does not match the number of exceptions (%2)") + .arg(exception_idxs.count()) + .arg(this->exception_atoms.count()), + CODELOC); + + this->exception_idxs.insert(name, exception_idxs); +} diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.h b/wrapper/Convert/SireOpenMM/openmmmolecule.h index 022876c65..8158e5c64 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.h +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.h @@ -18,7 +18,9 @@ namespace SireOpenMM { /** Internal class used to hold all of the extracted information - * of an OpenMM Molecule + * of an OpenMM Molecule. You should not use this outside + * of the sire_to_openmm_system function. It holds lots of scratch + * data that may change in the future. */ class OpenMMMolecule { @@ -28,7 +30,7 @@ namespace SireOpenMM CONSTRAIN_NONE = 0x0000, CONSTRAIN_BONDS = 0x0001, CONSTRAIN_HBONDS = 0x0010, - CONSTRAIN_ANGLES = 0x0100 + CONSTRAIN_HANGLES = 0x1000 }; OpenMMMolecule(); @@ -64,6 +66,11 @@ namespace SireOpenMM QVector getTorsionPhases() const; QVector getTorsionKs() const; + QVector> getExceptionAtoms() const; + + QVector getChargeScales() const; + QVector getLJScales() const; + bool isGhostAtom(int atom) const; std::tuple @@ -109,10 +116,7 @@ namespace SireOpenMM /** Charge and LJ parameters (sigma / epsilon) */ QVector> cljs; - /** Indexes of all bond pairs */ - QVector> bond_pairs; - - /** Set of 1-4 or excluded pairs + /** Set of 1-4 or excluded pairs (with coulomb and LJ scaling factors) */ QVector> exception_params; @@ -158,6 +162,10 @@ namespace SireOpenMM /** What type of constraint to use */ qint32 constraint_type; + /** What type of constraint to use when the bond/angle involves + * perturbable atoms */ + qint32 perturbable_constraint_type; + private: void constructFromAmber(const SireMol::Molecule &mol, const SireMM::AmberParams ¶ms, @@ -165,18 +173,112 @@ namespace SireOpenMM const SireBase::PropertyMap &map, bool is_perturbable); - void buildStandardExceptions(const SireMol::Molecule &mol, - const SireMM::AmberParams ¶ms, - QSet &constrained_pairs, - const SireBase::PropertyMap &map); - - void buildPerturbableExceptions(const SireMol::Molecule &mol, - QSet &constrained_pairs, - const SireBase::PropertyMap &map); + void buildExceptions(const SireMol::Molecule &mol, + QSet &constrained_pairs, + const SireBase::PropertyMap &map); void alignInternals(const SireBase::PropertyMap &map); }; + /** This class holds all of the information of an OpenMM molecule + * that can be perturbed using a LambdaSchedule. The data is held + * in easy-to-access arrays, with guarantees that the arrays are + * compatible and the data is aligned. + */ + class PerturbableOpenMMMolecule + { + public: + PerturbableOpenMMMolecule(); + + PerturbableOpenMMMolecule(const OpenMMMolecule &mol); + PerturbableOpenMMMolecule(const PerturbableOpenMMMolecule &other); + + ~PerturbableOpenMMMolecule(); + + bool operator==(const PerturbableOpenMMMolecule &other) const; + bool operator!=(const PerturbableOpenMMMolecule &other) const; + + QVector getAlphas0() const; + QVector getAlphas1() const; + + QVector getCharges0() const; + QVector getCharges1() const; + + QVector getSigmas0() const; + QVector getSigmas1() const; + QVector getEpsilons0() const; + QVector getEpsilons1() const; + + QVector getBondKs0() const; + QVector getBondKs1() const; + QVector getBondLengths0() const; + QVector getBondLengths1() const; + + QVector getAngleKs0() const; + QVector getAngleKs1() const; + QVector getAngleSizes0() const; + QVector getAngleSizes1() const; + + QVector getTorsionKs0() const; + QVector getTorsionKs1() const; + QVector getTorsionPeriodicities0() const; + QVector getTorsionPeriodicities1() const; + QVector getTorsionPhases0() const; + QVector getTorsionPhases1() const; + + QVector getChargeScales0() const; + QVector getChargeScales1() const; + QVector getLJScales0() const; + QVector getLJScales1() const; + + QSet getToGhostIdxs() const; + QSet getFromGhostIdxs() const; + + bool isGhostAtom(int atom) const; + + QVector> getExceptionAtoms() const; + + QVector> getExceptionIndicies(const QString &name) const; + + void setExceptionIndicies(const QString &name, + const QVector> &exception_idxs); + + private: + /** The array of parameters for the two end states, aligned + * so that they can be morphed via the LambdaLever + */ + QVector alpha0, alpha1; + QVector chg0, chg1; + QVector sig0, sig1; + QVector eps0, eps1; + QVector bond_k0, bond_k1; + QVector bond_r0, bond_r1; + QVector ang_k0, ang_k1; + QVector ang_t0, ang_t1; + QVector tors_k0, tors_k1; + QVector tors_periodicity0, tors_periodicity1; + QVector tors_phase0, tors_phase1; + QVector charge_scl0, charge_scl1; + QVector lj_scl0, lj_scl1; + + /** The indexes of atoms that become ghosts in the + * perturbed state + */ + QSet to_ghost_idxs; + + /** The indexes of atoms that are ghosts in the reference + * state and are real in the perturbed state + */ + QSet from_ghost_idxs; + + /** The indicies of the atoms in the exceptions, in exception order */ + QVector> exception_atoms; + + /** The indicies of the added exceptions - only populated + * if this is a peturbable molecule */ + QHash>> exception_idxs; + }; + } SIRE_END_HEADER diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 9c775561e..724bd861c 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -752,6 +752,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, lambda_lever.addLever("sigma"); lambda_lever.addLever("epsilon"); + // and the exceptions + lambda_lever.addLever("charge_scale"); + lambda_lever.addLever("lj_scale"); + // Do the same for the bond, angle and torsion forces lambda_lever.setForceIndex("bond", system.addForce(bondff)); lambda_lever.addLever("bond_length"); @@ -763,7 +767,6 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, lambda_lever.setForceIndex("torsion", system.addForce(dihff)); lambda_lever.addLever("torsion_phase"); - lambda_lever.addLever("torsion_periodicity"); lambda_lever.addLever("torsion_k"); /// @@ -1256,65 +1259,23 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, std::get<4>(dih), std::get<5>(dih), std::get<6>(dih)); } - // now add the constraints - if (any_perturbable and mol.isPerturbable()) + for (const auto &constraint : mol.constraints) { - // we may want to select out constraints that involve - // perturbing atoms... For now we will include them all, - // but will use the current coordinates to set the - // constrained values, rather than using the - // forcefield parameters. This will work as long - // as the user has minimised the system at the desired - // lambda value before running dynamics... - for (const auto &constraint : mol.constraints) - { - const auto atom0 = std::get<0>(constraint); - const auto atom1 = std::get<1>(constraint); - - const auto mass0 = system.getParticleMass(atom0 + start_index); - const auto mass1 = system.getParticleMass(atom1 + start_index); + const auto atom0 = std::get<0>(constraint); + const auto atom1 = std::get<1>(constraint); - if (mass0 != 0 and mass1 != 0) - { - // we can add the constraint - const auto coords0 = mol.coords[atom0]; - const auto coords1 = mol.coords[atom1]; - - const auto delta = coords1 - coords0; + const auto mass0 = system.getParticleMass(atom0 + start_index); + const auto mass1 = system.getParticleMass(atom1 + start_index); - const auto length = std::sqrt((delta[0] * delta[0]) + - (delta[1] * delta[1]) + - (delta[2] * delta[2])); - - system.addConstraint( - atom0 + start_index, - atom1 + start_index, - length); - } - // else we will need to think about how to constrain bonds - // involving fixed atoms. Could we fix the other atom too? - } - } - else - { - for (const auto &constraint : mol.constraints) + if (mass0 != 0 and mass1 != 0) { - const auto atom0 = std::get<0>(constraint); - const auto atom1 = std::get<1>(constraint); - const auto mass0 = system.getParticleMass(atom0 + start_index); - const auto mass1 = system.getParticleMass(atom1 + start_index); - - if (mass0 != 0 and mass1 != 0) - { - - system.addConstraint(atom0 + start_index, - atom1 + start_index, - std::get<2>(constraint)); - } - // else we will need to think about how to constrain bonds - // involving fixed atoms. Could we fix the other atom too? + system.addConstraint(atom0 + start_index, + atom1 + start_index, + std::get<2>(constraint)); } + // else we will need to think about how to constrain bonds + // involving fixed atoms. Could we fix the other atom too? } start_index += mol.masses.count(); @@ -1437,10 +1398,16 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, std::get<4>(p), true); } + // these are the indexes of the exception in the + // non-bonded forcefields and also the ghost-14 forcefield exception_idxs[j] = std::make_pair(idx, nbidx); - // remove this interaction from the ghost forcefields - if (ghost_ghostff != 0) + // remove this interaction from the ghost forcefields, only + // if it involves ghost atoms. If it doens't involve ghost atoms + // then we cannot remove it, because OpenMM doesn't support + // having a non-zero exception in NonbondedForce while having + // exclusions between the same atoms in CustomNonbondedForce + if (ghost_ghostff != 0 and (atom0_is_ghost or atom1_is_ghost)) { ghost_ghostff->addExclusion(std::get<0>(p), std::get<1>(p)); ghost_nonghostff->addExclusion(std::get<0>(p), std::get<1>(p)); diff --git a/wrapper/Convert/__init__.py b/wrapper/Convert/__init__.py index 677d3a298..f85c0f4e6 100644 --- a/wrapper/Convert/__init__.py +++ b/wrapper/Convert/__init__.py @@ -149,6 +149,9 @@ def sire_to_openmm(mols, map): friction = friction.to(1.0 / picosecond) / openmm.unit.picosecond + use_andersen = False + temperature = None + if integrator is None: if ensemble.is_nve(): integrator = openmm.VerletIntegrator(timestep) @@ -206,6 +209,12 @@ def sire_to_openmm(mols, map): temperature, friction, timestep ) + elif integrator == "andersen": + # use a verlet integrator and switch on the + # andersen thermostat with the specified frequency + integrator = openmm.VerletIntegrator(timestep) + use_andersen = True + else: raise ValueError(f"Unrecognised integrator {integrator}") @@ -220,6 +229,11 @@ def sire_to_openmm(mols, map): # system must be an openmm.System() or else we will crash! openmm_metadata = _sire_to_openmm_system(system, mols, map) + # If we want temperature controlled by an Andersen thermostat + # then add this here + if use_andersen: + system.addForce(openmm.AndersenThermostat(temperature, friction)) + # If we want NPT and this is periodic then we have to # add the barostat to the system if ensemble.is_npt(): diff --git a/wrapper/IO/_IO_load.cpp b/wrapper/IO/_IO_load.cpp index 916f2cfea..963f1b5b0 100644 --- a/wrapper/IO/_IO_load.cpp +++ b/wrapper/IO/_IO_load.cpp @@ -46,7 +46,7 @@ System load_molecules(const QStringList &files, { const auto &f = files.at(0); - if (f.endsWith(".s3")) + if (f.endsWith(".s3") or f.endsWith(".bss")) { // try to load this as a Sire s3 file try @@ -55,8 +55,22 @@ System load_molecules(const QStringList &files, } catch (...) { - // it couldn't be loaded - try one of the - // parsers below to try to read it + } + + // try to load this as a Sire s3 file + try + { + auto m = SireStream::loadType(f); + + auto s = System(); + s.setName(m.name()); + auto g = MoleculeGroup("all"); + g.add(m); + s.add(g); + return s; + } + catch (...) + { } } } diff --git a/wrapper/Maths/CMakeAutogenFile.txt b/wrapper/Maths/CMakeAutogenFile.txt index a525364c8..048bb92b8 100644 --- a/wrapper/Maths/CMakeAutogenFile.txt +++ b/wrapper/Maths/CMakeAutogenFile.txt @@ -1,54 +1,54 @@ # WARNING - AUTOGENERATED FILE - CONTENTS WILL BE OVERWRITTEN! set ( PYPP_SOURCES - Accumulator.pypp.cpp + AxisSet.pypp.cpp + ExpAverage.pypp.cpp + NVector.pypp.cpp + MultiUInt.pypp.cpp + Sphere.pypp.cpp + RecordValues.pypp.cpp Array2D_Matrix_.pypp.cpp - Array2D_NMatrix_.pypp.cpp + NullAccumulator.pypp.cpp + SphereArrayProperty.pypp.cpp + N4Matrix.pypp.cpp + AverageAndStddev.pypp.cpp + _Maths_free_functions.pypp.cpp + MultiFixed.pypp.cpp Array2D_SireMaths_AccumulatorPtr_.pypp.cpp - Array2D_Vector_.pypp.cpp + TrigArray2D_Matrix_.pypp.cpp ArrayProperty_Vector_.pypp.cpp + VectorArrayProperty.pypp.cpp Average.pypp.cpp - AverageAndStddev.pypp.cpp - AxisSet.pypp.cpp - BennettsFreeEnergyAverage.pypp.cpp - Complex.pypp.cpp + FreeEnergyAverage.pypp.cpp + SphereProperty.pypp.cpp + Median.pypp.cpp + VectorProperty.pypp.cpp DistVector.pypp.cpp + Vector.pypp.cpp EnergyTrajectory.pypp.cpp - ExpAverage.pypp.cpp - FreeEnergyAverage.pypp.cpp - Histogram.pypp.cpp - HistogramBin.pypp.cpp - HistogramValue.pypp.cpp Line.pypp.cpp + HistogramBin.pypp.cpp + Triangle.pypp.cpp + Quaternion.pypp.cpp Matrix.pypp.cpp - Median.pypp.cpp + Accumulator.pypp.cpp + Rational.pypp.cpp + Transform.pypp.cpp + HistogramValue.pypp.cpp + Array2D_Vector_.pypp.cpp + BennettsFreeEnergyAverage.pypp.cpp + TrigArray2D_Vector_.pypp.cpp MultiDouble.pypp.cpp - MultiFixed.pypp.cpp + Array2D_NMatrix_.pypp.cpp MultiFloat.pypp.cpp + Histogram.pypp.cpp MultiInt.pypp.cpp - MultiUInt.pypp.cpp - MultiVector.pypp.cpp - N4Matrix.pypp.cpp - NMatrix.pypp.cpp - NVector.pypp.cpp - NullAccumulator.pypp.cpp - Plane.pypp.cpp - Quaternion.pypp.cpp - RanGenerator.pypp.cpp - Rational.pypp.cpp - RecordValues.pypp.cpp - Sphere.pypp.cpp - SphereArrayProperty.pypp.cpp - SphereProperty.pypp.cpp Torsion.pypp.cpp - Transform.pypp.cpp - Triangle.pypp.cpp - TrigArray2D_Matrix_.pypp.cpp - TrigArray2D_Vector_.pypp.cpp + Complex.pypp.cpp TrigMatrix.pypp.cpp - Vector.pypp.cpp - VectorArrayProperty.pypp.cpp - VectorProperty.pypp.cpp - _Maths_free_functions.pypp.cpp + NMatrix.pypp.cpp + MultiVector.pypp.cpp + RanGenerator.pypp.cpp + Plane.pypp.cpp SireMaths_containers.cpp SireMaths_properties.cpp SireMaths_registrars.cpp diff --git a/wrapper/Maths/EnergyTrajectory.pypp.cpp b/wrapper/Maths/EnergyTrajectory.pypp.cpp index c0cc4410d..0bad36b1a 100644 --- a/wrapper/Maths/EnergyTrajectory.pypp.cpp +++ b/wrapper/Maths/EnergyTrajectory.pypp.cpp @@ -3,6 +3,7 @@ // (C) Christopher Woods, GPL >= 3 License #include "boost/python.hpp" +#include "Helpers/clone_const_reference.hpp" #include "EnergyTrajectory.pypp.hpp" namespace bp = boost::python; @@ -38,6 +39,18 @@ void register_EnergyTrajectory_class(){ EnergyTrajectory_exposer_t EnergyTrajectory_exposer = EnergyTrajectory_exposer_t( "EnergyTrajectory", "This class holds the trajectory of energies, organised by\ntimestep the energy was recorded and the types of energy\n(e.g. kinetic, potential, values at different lambda windows)\n", bp::init< >("") ); bp::scope EnergyTrajectory_scope( EnergyTrajectory_exposer ); EnergyTrajectory_exposer.def( bp::init< SireMaths::EnergyTrajectory const & >(( bp::arg("other") ), "") ); + { //::SireMaths::EnergyTrajectory::clearProperties + + typedef void ( ::SireMaths::EnergyTrajectory::*clearProperties_function_type)( ) ; + clearProperties_function_type clearProperties_function_value( &::SireMaths::EnergyTrajectory::clearProperties ); + + EnergyTrajectory_exposer.def( + "clearProperties" + , clearProperties_function_value + , bp::release_gil_policy() + , "" ); + + } { //::SireMaths::EnergyTrajectory::count typedef int ( ::SireMaths::EnergyTrajectory::*count_function_type)( ) const; @@ -86,7 +99,7 @@ void register_EnergyTrajectory_class(){ , get_function_value , ( bp::arg("i") ) , bp::release_gil_policy() - , "" ); + , "Return the time and energy components at the ith row.\n Values are returned in internal units\n" ); } { //::SireMaths::EnergyTrajectory::get @@ -99,6 +112,45 @@ void register_EnergyTrajectory_class(){ , get_function_value , ( bp::arg("i"), bp::arg("time_unit"), bp::arg("energy_unit") ) , bp::release_gil_policy() + , "Return the time and energy components at the ith row.\n Values are returned in the specified units\n" ); + + } + { //::SireMaths::EnergyTrajectory::getLabels + + typedef ::QHash< QString, QString > ( ::SireMaths::EnergyTrajectory::*getLabels_function_type)( int ) const; + getLabels_function_type getLabels_function_value( &::SireMaths::EnergyTrajectory::getLabels ); + + EnergyTrajectory_exposer.def( + "getLabels" + , getLabels_function_value + , ( bp::arg("i") ) + , bp::release_gil_policy() + , "" ); + + } + { //::SireMaths::EnergyTrajectory::getLabelsAsNumbers + + typedef ::QHash< QString, double > ( ::SireMaths::EnergyTrajectory::*getLabelsAsNumbers_function_type)( int ) const; + getLabelsAsNumbers_function_type getLabelsAsNumbers_function_value( &::SireMaths::EnergyTrajectory::getLabelsAsNumbers ); + + EnergyTrajectory_exposer.def( + "getLabelsAsNumbers" + , getLabelsAsNumbers_function_value + , ( bp::arg("i") ) + , bp::release_gil_policy() + , "" ); + + } + { //::SireMaths::EnergyTrajectory::hasProperty + + typedef bool ( ::SireMaths::EnergyTrajectory::*hasProperty_function_type)( ::SireBase::PropertyName const & ) ; + hasProperty_function_type hasProperty_function_value( &::SireMaths::EnergyTrajectory::hasProperty ); + + EnergyTrajectory_exposer.def( + "hasProperty" + , hasProperty_function_value + , ( bp::arg("key") ) + , bp::release_gil_policy() , "" ); } @@ -137,6 +189,44 @@ void register_EnergyTrajectory_class(){ , bp::release_gil_policy() , "Return all of the energy keys" ); + } + { //::SireMaths::EnergyTrajectory::labelKeys + + typedef ::QStringList ( ::SireMaths::EnergyTrajectory::*labelKeys_function_type)( ) const; + labelKeys_function_type labelKeys_function_value( &::SireMaths::EnergyTrajectory::labelKeys ); + + EnergyTrajectory_exposer.def( + "labelKeys" + , labelKeys_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMaths::EnergyTrajectory::labels + + typedef ::QVector< QString > ( ::SireMaths::EnergyTrajectory::*labels_function_type)( ::QString const & ) const; + labels_function_type labels_function_value( &::SireMaths::EnergyTrajectory::labels ); + + EnergyTrajectory_exposer.def( + "labels" + , labels_function_value + , ( bp::arg("key") ) + , bp::release_gil_policy() + , "" ); + + } + { //::SireMaths::EnergyTrajectory::labelsAsNumbers + + typedef ::QVector< double > ( ::SireMaths::EnergyTrajectory::*labelsAsNumbers_function_type)( ::QString const & ) const; + labelsAsNumbers_function_type labelsAsNumbers_function_value( &::SireMaths::EnergyTrajectory::labelsAsNumbers ); + + EnergyTrajectory_exposer.def( + "labelsAsNumbers" + , labelsAsNumbers_function_value + , ( bp::arg("key") ) + , bp::release_gil_policy() + , "" ); + } EnergyTrajectory_exposer.def( bp::self != bp::self ); { //::SireMaths::EnergyTrajectory::operator= @@ -164,6 +254,56 @@ void register_EnergyTrajectory_class(){ , ( bp::arg("i") ) , "" ); + } + { //::SireMaths::EnergyTrajectory::properties + + typedef ::SireBase::Properties const & ( ::SireMaths::EnergyTrajectory::*properties_function_type)( ) const; + properties_function_type properties_function_value( &::SireMaths::EnergyTrajectory::properties ); + + EnergyTrajectory_exposer.def( + "properties" + , properties_function_value + , bp::return_value_policy< bp::copy_const_reference >() + , "" ); + + } + { //::SireMaths::EnergyTrajectory::property + + typedef ::SireBase::Property const & ( ::SireMaths::EnergyTrajectory::*property_function_type)( ::SireBase::PropertyName const & ) const; + property_function_type property_function_value( &::SireMaths::EnergyTrajectory::property ); + + EnergyTrajectory_exposer.def( + "property" + , property_function_value + , ( bp::arg("key") ) + , bp::return_value_policy() + , "" ); + + } + { //::SireMaths::EnergyTrajectory::propertyKeys + + typedef ::QStringList ( ::SireMaths::EnergyTrajectory::*propertyKeys_function_type)( ) const; + propertyKeys_function_type propertyKeys_function_value( &::SireMaths::EnergyTrajectory::propertyKeys ); + + EnergyTrajectory_exposer.def( + "propertyKeys" + , propertyKeys_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMaths::EnergyTrajectory::removeProperty + + typedef void ( ::SireMaths::EnergyTrajectory::*removeProperty_function_type)( ::QString const & ) ; + removeProperty_function_type removeProperty_function_value( &::SireMaths::EnergyTrajectory::removeProperty ); + + EnergyTrajectory_exposer.def( + "removeProperty" + , removeProperty_function_value + , ( bp::arg("key") ) + , bp::release_gil_policy() + , "" ); + } { //::SireMaths::EnergyTrajectory::set @@ -177,6 +317,32 @@ void register_EnergyTrajectory_class(){ , bp::release_gil_policy() , "Set the energies at time time to the components contained\n in energies\n" ); + } + { //::SireMaths::EnergyTrajectory::set + + typedef void ( ::SireMaths::EnergyTrajectory::*set_function_type)( ::SireUnits::Dimension::GeneralUnit const &,::QHash< QString, SireUnits::Dimension::GeneralUnit > const &,::QHash< QString, QString > const & ) ; + set_function_type set_function_value( &::SireMaths::EnergyTrajectory::set ); + + EnergyTrajectory_exposer.def( + "set" + , set_function_value + , ( bp::arg("time"), bp::arg("energies"), bp::arg("labels") ) + , bp::release_gil_policy() + , "" ); + + } + { //::SireMaths::EnergyTrajectory::setProperty + + typedef void ( ::SireMaths::EnergyTrajectory::*setProperty_function_type)( ::QString const &,::SireBase::Property const & ) ; + setProperty_function_type setProperty_function_value( &::SireMaths::EnergyTrajectory::setProperty ); + + EnergyTrajectory_exposer.def( + "setProperty" + , setProperty_function_value + , ( bp::arg("key"), bp::arg("value") ) + , bp::release_gil_policy() + , "" ); + } { //::SireMaths::EnergyTrajectory::size