diff --git a/corelib/src/libs/SireCAS/lambdaschedule.cpp b/corelib/src/libs/SireCAS/lambdaschedule.cpp index 11eae8e35..e1cc50913 100644 --- a/corelib/src/libs/SireCAS/lambdaschedule.cpp +++ b/corelib/src/libs/SireCAS/lambdaschedule.cpp @@ -1032,8 +1032,12 @@ bool LambdaSchedule::hasForceSpecificEquation(const QString &stage, if (force == "*") return false; - else - return this->stage_equations[idx].contains(_get_lever_name(force, lever)); + + // If there is a default equation for this force, always return true + if (this->stage_equations[idx].contains(_get_lever_name(force, "*"))) + return true; + + return this->stage_equations[idx].contains(_get_lever_name(force, lever)); } SireCAS::Expression LambdaSchedule::_getEquation(int stage, diff --git a/corelib/src/libs/SireIO/sdf.cpp b/corelib/src/libs/SireIO/sdf.cpp index c5043dee0..44b9cd193 100644 --- a/corelib/src/libs/SireIO/sdf.cpp +++ b/corelib/src/libs/SireIO/sdf.cpp @@ -959,6 +959,13 @@ SDFMolecule parseMolecule(const Molecule &molecule, QStringList &errors, const P sdfmol.addProperty("CHG", i + 1, QString::number(charge)); charge = 0; } + else if (charge > 0) + { + if (charge < 4) + { + charge = 4 - charge; + } + } else if (charge < 0) { switch (charge) diff --git a/corelib/src/libs/SireMM/CMakeLists.txt b/corelib/src/libs/SireMM/CMakeLists.txt index f228580aa..f0e650fba 100644 --- a/corelib/src/libs/SireMM/CMakeLists.txt +++ b/corelib/src/libs/SireMM/CMakeLists.txt @@ -16,7 +16,7 @@ include_directories(${TBB_INCLUDE_DIR}) # Define the headers in SireMM set ( SIREMM_HEADERS amberparams.h - anglerestraint.h + anglerestraints.h angle.h atomfunctions.h atomljs.h @@ -43,7 +43,7 @@ set ( SIREMM_HEADERS cljworkspace.h coulombpotential.h dihedral.h - dihedralrestraint.h + dihedralrestraints.h distancerestraint.h errors.h excludedpairs.h @@ -110,7 +110,7 @@ set ( SIREMM_SOURCES amberparams.cpp angle.cpp - anglerestraint.cpp + anglerestraints.cpp atomfunctions.cpp atomljs.cpp bond.cpp @@ -135,7 +135,7 @@ set ( SIREMM_SOURCES cljworkspace.cpp coulombpotential.cpp dihedral.cpp - dihedralrestraint.cpp + dihedralrestraints.cpp distancerestraint.cpp errors.cpp excludedpairs.cpp diff --git a/corelib/src/libs/SireMM/anglerestraint.cpp b/corelib/src/libs/SireMM/anglerestraint.cpp deleted file mode 100644 index a765049c0..000000000 --- a/corelib/src/libs/SireMM/anglerestraint.cpp +++ /dev/null @@ -1,555 +0,0 @@ -/********************************************\ - * - * Sire - Molecular Simulation Framework - * - * Copyright (C) 2009 Christopher Woods - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * - * For full details of the license please see the COPYING file - * that should have come with this distribution. - * - * You can contact the authors at https://sire.openbiosim.org - * -\*********************************************/ - -#include "anglerestraint.h" - -#include "SireFF/forcetable.h" - -#include "SireCAS/conditional.h" -#include "SireCAS/power.h" -#include "SireCAS/symbols.h" -#include "SireCAS/values.h" - -#include "SireID/index.h" - -#include "SireUnits/angle.h" -#include "SireUnits/units.h" - -#include "SireStream/datastream.h" -#include "SireStream/shareddatastream.h" - -#include "SireCAS/errors.h" - -using namespace SireMM; -using namespace SireMol; -using namespace SireFF; -using namespace SireID; -using namespace SireBase; -using namespace SireCAS; -using namespace SireMaths; -using namespace SireStream; -using namespace SireUnits; -using namespace SireUnits::Dimension; - -//////////// -//////////// Implementation of AngleRestraint -//////////// - -static const RegisterMetaType r_angrest; - -/** Serialise to a binary datastream */ -QDataStream &operator<<(QDataStream &ds, const AngleRestraint &angrest) -{ - writeHeader(ds, r_angrest, 1); - - SharedDataStream sds(ds); - - sds << angrest.p[0] << angrest.p[1] << angrest.p[2] << angrest.force_expression - << static_cast(angrest); - - return ds; -} - -/** Extract from a binary datastream */ -QDataStream &operator>>(QDataStream &ds, AngleRestraint &angrest) -{ - VersionID v = readHeader(ds, r_angrest); - - if (v == 1) - { - SharedDataStream sds(ds); - - sds >> angrest.p[0] >> angrest.p[1] >> angrest.p[2] >> angrest.force_expression >> - static_cast(angrest); - - angrest.intra_molecule_points = Point::areIntraMoleculePoints(angrest.p[0], angrest.p[1]) and - Point::areIntraMoleculePoints(angrest.p[0], angrest.p[2]); - } - else - throw version_error(v, "1", r_angrest, CODELOC); - - return ds; -} - -Q_GLOBAL_STATIC_WITH_ARGS(Symbol, getThetaSymbol, ("theta")); - -/** Return the symbol that represents the angle between the points (theta) */ -const Symbol &AngleRestraint::theta() -{ - return *(getThetaSymbol()); -} - -/** Constructor */ -AngleRestraint::AngleRestraint() : ConcreteProperty() -{ -} - -void AngleRestraint::calculateTheta() -{ - if (this->restraintFunction().isFunction(theta())) - { - SireUnits::Dimension::Angle angle; - - if (intra_molecule_points) - // we don't use the space when calculating intra-molecular angles - angle = Vector::angle(p[0].read().point(), p[1].read().point(), p[2].read().point()); - else - angle = this->space().calcAngle(p[0].read().point(), p[1].read().point(), p[2].read().point()); - - ExpressionRestraint3D::_pvt_setValue(theta(), angle); - } -} - -/** Construct a restraint that acts on the angle within the - three points 'point0', 'point1' and 'point2' (theta == a(012)), - restraining the angle within these points using the expression - 'restraint' */ -AngleRestraint::AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Expression &restraint) - : ConcreteProperty(restraint) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - - force_expression = this->restraintFunction().differentiate(theta()); - - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]); - - this->calculateTheta(); -} - -/** Construct a restraint that acts on the angle within the - three points 'point0', 'point1' and 'point2' (theta == a(012)), - restraining the angle within these points using the expression - 'restraint' */ -AngleRestraint::AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Expression &restraint, const Values &values) - : ConcreteProperty(restraint, values) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - - force_expression = this->restraintFunction().differentiate(theta()); - - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]); - - this->calculateTheta(); -} - -/** Internal constructor used to construct a restraint using the specified - points, energy expression and force expression */ -AngleRestraint::AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Expression &nrg_restraint, const Expression &force_restraint) - : ConcreteProperty(nrg_restraint), force_expression(force_restraint) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - - if (force_expression.isConstant()) - { - force_expression = force_expression.evaluate(Values()); - } - else - { - if (not this->restraintFunction().symbols().contains(force_expression.symbols())) - throw SireError::incompatible_error( - QObject::tr("You cannot use a force function which uses more symbols " - "(%1) than the energy function (%2).") - .arg(Sire::toString(force_expression.symbols()), Sire::toString(restraintFunction().symbols())), - CODELOC); - } - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]); - - this->calculateTheta(); -} - -/** Copy constructor */ -AngleRestraint::AngleRestraint(const AngleRestraint &other) - : ConcreteProperty(other), force_expression(other.force_expression), - intra_molecule_points(other.intra_molecule_points) -{ - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } -} - -/** Destructor */ -AngleRestraint::~AngleRestraint() -{ -} - -/** Copy assignment operator */ -AngleRestraint &AngleRestraint::operator=(const AngleRestraint &other) -{ - if (this != &other) - { - ExpressionRestraint3D::operator=(other); - - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } - - force_expression = other.force_expression; - intra_molecule_points = other.intra_molecule_points; - } - - return *this; -} - -/** Comparison operator */ -bool AngleRestraint::operator==(const AngleRestraint &other) const -{ - return this == &other or (ExpressionRestraint3D::operator==(other) and p[0] == other.p[0] and p[1] == other.p[1] and - p[2] == other.p[2] and force_expression == other.force_expression); -} - -/** Comparison operator */ -bool AngleRestraint::operator!=(const AngleRestraint &other) const -{ - return not AngleRestraint::operator==(other); -} - -const char *AngleRestraint::typeName() -{ - return QMetaType::typeName(qMetaTypeId()); -} - -/** This restraint involves three points */ -int AngleRestraint::nPoints() const -{ - return 3; -} - -/** Return the ith point */ -const Point &AngleRestraint::point(int i) const -{ - i = Index(i).map(this->nPoints()); - - return p[i].read(); -} - -/** Return the first point */ -const Point &AngleRestraint::point0() const -{ - return p[0].read(); -} - -/** Return the second point */ -const Point &AngleRestraint::point1() const -{ - return p[1].read(); -} - -/** Return the third point */ -const Point &AngleRestraint::point2() const -{ - return p[2].read(); -} - -/** Return the built-in symbols of this restraint */ -Symbols AngleRestraint::builtinSymbols() const -{ - if (this->restraintFunction().isFunction(theta())) - return theta(); - else - return Symbols(); -} - -/** Return the values of the built-in symbols of this restraint */ -Values AngleRestraint::builtinValues() const -{ - if (this->restraintFunction().isFunction(theta())) - return theta() == this->values()[theta()]; - else - return Values(); -} - -/** Return the differential of this restraint with respect to - the symbol 'symbol' - - \throw SireCAS::unavailable_differential -*/ -RestraintPtr AngleRestraint::differentiate(const Symbol &symbol) const -{ - if (this->restraintFunction().isFunction(symbol)) - return AngleRestraint(p[0], p[1], p[2], restraintFunction().differentiate(symbol), this->values()); - else - return NullRestraint(); -} - -/** Set the space used to evaluate the energy of this restraint - - \throw SireVol::incompatible_space -*/ -void AngleRestraint::setSpace(const Space &new_space) -{ - if (not this->space().equals(new_space)) - { - AngleRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().setSpace(new_space); - } - - Restraint3D::setSpace(new_space); - - this->calculateTheta(); - } - catch (...) - { - AngleRestraint::operator=(old_state); - throw; - } - } -} - -/** Return the function used to calculate the restraint force */ -const Expression &AngleRestraint::differentialRestraintFunction() const -{ - return force_expression; -} - -/** Calculate the force acting on the molecule in the forcetable 'forcetable' - caused by this restraint, and add it on to the forcetable scaled by - 'scale_force' */ -void AngleRestraint::force(MolForceTable &forcetable, double scale_force) const -{ - bool in_p0 = p[0].read().contains(forcetable.molNum()); - bool in_p1 = p[1].read().contains(forcetable.molNum()); - bool in_p2 = p[2].read().contains(forcetable.molNum()); - - if (not(in_p0 or in_p1 or in_p2)) - // this molecule is not affected by the restraint - return; - - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by an angle restraint."), - CODELOC); -} - -/** Calculate the force acting on the molecules in the forcetable 'forcetable' - caused by this restraint, and add it on to the forcetable scaled by - 'scale_force' */ -void AngleRestraint::force(ForceTable &forcetable, double scale_force) const -{ - bool in_p0 = p[0].read().usesMoleculesIn(forcetable); - bool in_p1 = p[1].read().usesMoleculesIn(forcetable); - bool in_p2 = p[2].read().usesMoleculesIn(forcetable); - - if (not(in_p0 or in_p1 or in_p2)) - // this molecule is not affected by the restraint - return; - - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by an angle restraint."), - CODELOC); -} - -/** Update the points of this restraint using new molecule data from 'moldata' - - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void AngleRestraint::update(const MoleculeData &moldata) -{ - if (this->contains(moldata.number())) - { - AngleRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(moldata); - } - - this->calculateTheta(); - } - catch (...) - { - AngleRestraint::operator=(old_state); - throw; - } - } -} - -/** Update the points of this restraint using new molecule data from 'molecules' - - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void AngleRestraint::update(const Molecules &molecules) -{ - if (this->usesMoleculesIn(molecules)) - { - AngleRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(molecules); - } - - this->calculateTheta(); - } - catch (...) - { - AngleRestraint::operator=(old_state); - throw; - } - } -} - -/** Return the molecules used in this restraint */ -Molecules AngleRestraint::molecules() const -{ - Molecules mols; - - for (int i = 0; i < this->nPoints(); ++i) - { - mols += p[i].read().molecules(); - } - - return mols; -} - -/** Return whether or not this restraint affects the molecule - with number 'molnum' */ -bool AngleRestraint::contains(MolNum molnum) const -{ - return p[0].read().contains(molnum) or p[1].read().contains(molnum) or p[2].read().contains(molnum); -} - -/** Return whether or not this restraint affects the molecule - with ID 'molid' */ -bool AngleRestraint::contains(const MolID &molid) const -{ - return p[0].read().contains(molid) or p[1].read().contains(molid) or p[2].read().contains(molid); -} - -/** Return whether or not this restraint involves any of the molecules - that are in the forcetable 'forcetable' */ -bool AngleRestraint::usesMoleculesIn(const ForceTable &forcetable) const -{ - return p[0].read().usesMoleculesIn(forcetable) or p[1].read().usesMoleculesIn(forcetable) or - p[2].read().usesMoleculesIn(forcetable); -} - -/** Return whether or not this restraint involves any of the molecules - in 'molecules' */ -bool AngleRestraint::usesMoleculesIn(const Molecules &molecules) const -{ - return p[0].read().usesMoleculesIn(molecules) or p[1].read().usesMoleculesIn(molecules) or - p[2].read().usesMoleculesIn(molecules); -} - -static Expression harmonicFunction(double force_constant) -{ - if (SireMaths::isZero(force_constant)) - return 0; - else - return force_constant * pow(AngleRestraint::theta(), 2); -} - -static Expression diffHarmonicFunction(double force_constant) -{ - if (SireMaths::isZero(force_constant)) - return 0; - else - return (2 * force_constant) * AngleRestraint::theta(); -} - -/** Return a distance restraint that applies a harmonic potential between - the points 'point0' and 'point1' using a force constant 'force_constant' */ -AngleRestraint AngleRestraint::harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const HarmonicAngleForceConstant &force_constant) -{ - return AngleRestraint(point0, point1, point2, ::harmonicFunction(force_constant), - ::diffHarmonicFunction(force_constant)); -} - -static Expression halfHarmonicFunction(double force_constant, double angle) -{ - if (SireMaths::isZero(force_constant)) - return 0; - - else if (angle <= 0) - // this is just a harmonic function - return ::harmonicFunction(force_constant); - - else - { - const Symbol &theta = AngleRestraint::theta(); - return Conditional(GreaterThan(theta, angle), force_constant * pow(theta - angle, 2), 0); - } -} - -static Expression diffHalfHarmonicFunction(double force_constant, double angle) -{ - if (SireMaths::isZero(force_constant)) - return 0; - - else if (angle <= 0) - // this is just a harmonic function - return ::diffHarmonicFunction(force_constant); - - else - { - const Symbol &theta = AngleRestraint::theta(); - return Conditional(GreaterThan(theta, angle), (2 * force_constant) * (theta - angle), 0); - } -} - -/** Return a distance restraint that applied a half-harmonic potential - between the points 'point0' and 'point1' above a distance 'distance' - using a force constant 'force_constant' */ -AngleRestraint AngleRestraint::halfHarmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Angle &angle, const HarmonicAngleForceConstant &force_constant) -{ - double acute_angle = acute(angle).to(radians); - - return AngleRestraint(point0, point1, point2, ::halfHarmonicFunction(force_constant, acute_angle), - ::diffHalfHarmonicFunction(force_constant, acute_angle)); -} diff --git a/corelib/src/libs/SireMM/anglerestraint.h b/corelib/src/libs/SireMM/anglerestraint.h deleted file mode 100644 index 6901c78ca..000000000 --- a/corelib/src/libs/SireMM/anglerestraint.h +++ /dev/null @@ -1,159 +0,0 @@ -/********************************************\ - * - * Sire - Molecular Simulation Framework - * - * Copyright (C) 2009 Christopher Woods - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * - * For full details of the license please see the COPYING file - * that should have come with this distribution. - * - * You can contact the authors at https://sire.openbiosim.org - * -\*********************************************/ - -#ifndef SIREMM_ANGLERESTRAINT_H -#define SIREMM_ANGLERESTRAINT_H - -#include "SireFF/point.h" - -#include "restraint.h" - -#include "SireCAS/expression.h" -#include "SireCAS/symbol.h" - -#include "SireUnits/dimensions.h" - -SIRE_BEGIN_HEADER - -namespace SireMM -{ - class AngleRestraint; -} - -SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::AngleRestraint &); -SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::AngleRestraint &); - -namespace SireMM -{ - - using SireFF::Point; - using SireFF::PointRef; - - using SireCAS::Expression; - using SireCAS::Symbol; - using SireCAS::Symbols; - - // typedef the unit of a harmonic force constant ( MolarEnergy / Angle^2 ) - typedef SireUnits::Dimension::PhysUnit<1, 2, -2, 0, 0, -1, -2> HarmonicAngleForceConstant; - - /** This is a restraint that operates on the angle between - three SireMM::Point objects (e.g. three atoms in a molecule) - - @author Christopher Woods - */ - class SIREMM_EXPORT AngleRestraint : public SireBase::ConcreteProperty - { - - friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const AngleRestraint &); - friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, AngleRestraint &); - - public: - AngleRestraint(); - - AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const Expression &restraint); - - AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const Expression &restraint, - const Values &values); - - AngleRestraint(const AngleRestraint &other); - - ~AngleRestraint(); - - AngleRestraint &operator=(const AngleRestraint &other); - - bool operator==(const AngleRestraint &other) const; - bool operator!=(const AngleRestraint &other) const; - - static const char *typeName(); - - const Point &point(int i) const; - - const Point &point0() const; - const Point &point1() const; - const Point &point2() const; - - int nPoints() const; - - static const Symbol &theta(); - - Symbols builtinSymbols() const; - Values builtinValues() const; - - RestraintPtr differentiate(const Symbol &symbol) const; - - void setSpace(const Space &space); - - const Expression &differentialRestraintFunction() const; - - void force(MolForceTable &forcetable, double scale_force = 1) const; - void force(ForceTable &forcetable, double scale_force = 1) const; - - void update(const MoleculeData &moldata); - - void update(const Molecules &molecules); - - Molecules molecules() const; - - bool contains(MolNum molnum) const; - bool contains(const MolID &molid) const; - - bool usesMoleculesIn(const ForceTable &forcetable) const; - bool usesMoleculesIn(const Molecules &molecules) const; - - static AngleRestraint harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const HarmonicAngleForceConstant &force_constant); - - static AngleRestraint halfHarmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const SireUnits::Dimension::Angle &angle, - const HarmonicAngleForceConstant &force_constant); - - protected: - AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Expression &nrg_restraint, const Expression &force_restraint); - - private: - void calculateTheta(); - - /** The three points between which the restraint is calculated */ - SireFF::PointPtr p[3]; - - /** The expression used to calculate the force */ - Expression force_expression; - - /** Whether or not all three points are within the same molecule */ - bool intra_molecule_points; - }; - -} // namespace SireMM - -Q_DECLARE_METATYPE(SireMM::AngleRestraint) - -SIRE_EXPOSE_CLASS(SireMM::AngleRestraint) - -SIRE_END_HEADER - -#endif diff --git a/corelib/src/libs/SireMM/anglerestraints.cpp b/corelib/src/libs/SireMM/anglerestraints.cpp new file mode 100644 index 000000000..5eaf82fc1 --- /dev/null +++ b/corelib/src/libs/SireMM/anglerestraints.cpp @@ -0,0 +1,479 @@ +/********************************************\ + * + * Sire - Molecular Simulation Framework + * + * Copyright (C) 2025 Christopher Woods + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * For full details of the license please see the COPYING file + * that should have come with this distribution. + * + * You can contact the authors at https://sire.openbiosim.org + * +\*********************************************/ + +#include "anglerestraints.h" + +#include "SireID/index.h" + +#include "SireUnits/units.h" + +#include "SireStream/datastream.h" +#include "SireStream/shareddatastream.h" + +#include "SireCAS/errors.h" + +#include + +using namespace SireMM; +using namespace SireID; +using namespace SireBase; +using namespace SireMaths; +using namespace SireStream; +using namespace SireUnits; +using namespace SireUnits::Dimension; + +//////////// +//////////// Implementation of AngleRestraint +//////////// + +static const RegisterMetaType r_angrest; + +/** Serialise to a binary datastream */ + +QDataStream &operator<<(QDataStream &ds, const AngleRestraint &angrest) +{ + writeHeader(ds, r_angrest, 1); + + SharedDataStream sds(ds); + + sds << angrest.atms << angrest._theta0 << angrest._ktheta; + + return ds; +} + +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, AngleRestraint &angrest) +{ + VersionID v = readHeader(ds, r_angrest); + + if (v == 1) + { + SharedDataStream sds(ds); + + sds >> angrest.atms >> angrest._theta0 >> angrest._ktheta; + } + else + throw version_error(v, "1", r_angrest, CODELOC); + + return ds; +} + +/** Null constructor */ +AngleRestraint::AngleRestraint() + : ConcreteProperty(), + _ktheta(0), _theta0(0) +{ +} + +/** Construct a restraint that acts on the angle within the + three atoms 'atom0', 'atom1' and 'atom2' (theta == a(012)), + restraining the angle within these atoms */ +AngleRestraint::AngleRestraint(const QList &atoms, + const SireUnits::Dimension::Angle &theta0, + const SireUnits::Dimension::HarmonicAngleConstant &ktheta) + : ConcreteProperty(), + _theta0(theta0), _ktheta(ktheta) +{ + + // Make sure that we have 3 distinct atoms + QSet distinct; + distinct.reserve(3); + + for (const auto &atom : atoms) + { + if (atom >= 0) + distinct.insert(atom); + } + + atms = atoms.toVector(); +} + +/* Copy constructor*/ +AngleRestraint::AngleRestraint(const AngleRestraint &other) + : ConcreteProperty(other), + atms(other.atms), _theta0(other._theta0), _ktheta(other._ktheta) + +{ +} + +/* Destructor */ +AngleRestraint::~AngleRestraint() +{ +} + +AngleRestraint &AngleRestraint::operator=(const AngleRestraint &other) +{ + if (this != &other) + { + Property::operator=(other); + atms = other.atms; + _theta0 = other._theta0; + _ktheta = other._ktheta; + } + + return *this; +} + +bool AngleRestraint::operator==(const AngleRestraint &other) const +{ + return atms == other.atms and + _theta0 == other._theta0 and + _ktheta == other._ktheta; +} + +bool AngleRestraint::operator!=(const AngleRestraint &other) const +{ + return not operator==(other); +} + +AngleRestraints AngleRestraint::operator+(const AngleRestraint &other) const +{ + return AngleRestraints(*this) + other; +} + +AngleRestraints AngleRestraint::operator+(const AngleRestraints &other) const +{ + return AngleRestraints(*this) + other; +} + +const char *AngleRestraint::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); +} + +const char *AngleRestraint::what() const +{ + return AngleRestraint::typeName(); +} + +AngleRestraint *AngleRestraint::clone() const +{ + return new AngleRestraint(*this); +} + +bool AngleRestraint::isNull() const +{ + return atms.isEmpty(); +} + +QString AngleRestraint::toString() const +{ + if (this->isNull()) + return QObject::tr("AngleRestraint::null"); + else + { + QStringList a; + + for (const auto &atom : atms) + { + a.append(QString::number(atom)); + } + return QString("AngleRestraint( [%1], theta0=%2, ktheta=%3 )") + .arg(a.join(", ")) + .arg(_theta0.toString()) + .arg(_ktheta.toString()); + } +} + +/** Return the force constant for the restraint */ +SireUnits::Dimension::HarmonicAngleConstant AngleRestraint::ktheta() const +{ + return this->_ktheta; +} + +/** Return the equilibrium angle for the restraint */ +SireUnits::Dimension::Angle AngleRestraint::theta0() const +{ + return this->_theta0; +} + +/** Return the atoms involved in the restraint */ +QVector AngleRestraint::atoms() const +{ + return this->atms; +} + +/////// +/////// Implementation of AngleRestraints +/////// + +/** Serialise to a binary datastream */ + +static const RegisterMetaType r_angrests; + +QDataStream &operator<<(QDataStream &ds, const AngleRestraints &angrests) +{ + writeHeader(ds, r_angrests, 1); + + SharedDataStream sds(ds); + + sds << angrests.r + << static_cast(angrests); + + return ds; +} + +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, AngleRestraints &angrests) +{ + VersionID v = readHeader(ds, r_angrests); + + if (v == 1) + { + SharedDataStream sds(ds); + + sds >> angrests.r >> + static_cast(angrests); + } + else + throw version_error(v, "1", r_angrests, CODELOC); + + return ds; +} + +/** Null constructor */ +AngleRestraints::AngleRestraints() + : ConcreteProperty() +{ +} + +AngleRestraints::AngleRestraints(const QString &name) + : ConcreteProperty(name) +{ +} + +AngleRestraints::AngleRestraints(const AngleRestraint &restraint) + : ConcreteProperty() +{ + if (not restraint.isNull()) + r.append(restraint); +} + +AngleRestraints::AngleRestraints(const QList &restraints) + : ConcreteProperty() +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } +} + +AngleRestraints::AngleRestraints(const QString &name, + const AngleRestraint &restraint) + : ConcreteProperty(name) +{ + if (not restraint.isNull()) + r.append(restraint); +} + +AngleRestraints::AngleRestraints(const QString &name, + const QList &restraints) + : ConcreteProperty(name) +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } +} + +AngleRestraints::AngleRestraints(const AngleRestraints &other) + : ConcreteProperty(other), r(other.r) +{ +} + +/* Desctructor */ +AngleRestraints::~AngleRestraints() +{ +} + +AngleRestraints &AngleRestraints::operator=(const AngleRestraints &other) +{ + if (this != &other) + { + Restraints::operator=(other); + r = other.r; + } + + return *this; +} + +bool AngleRestraints::operator==(const AngleRestraints &other) const +{ + return r == other.r and Restraints::operator==(other); +} + +bool AngleRestraints::operator!=(const AngleRestraints &other) const +{ + return not operator==(other); +} + +const char *AngleRestraints::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); +} + +const char *AngleRestraints::what() const +{ + return AngleRestraints::typeName(); +} + +AngleRestraints *AngleRestraints::clone() const +{ + return new AngleRestraints(*this); +} + +QString AngleRestraints::toString() const +{ + if (this->isEmpty()) + return QObject::tr("AngleRestraints::null"); + + QStringList parts; + + const auto n = this->count(); + + if (n <= 10) + { + for (int i = 0; i < n; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + } + else + { + for (int i = 0; i < 5; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + + parts.append("..."); + + for (int i = n - 5; i < n; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + } + + return QObject::tr("AngleRestraints( name=%1, size=%2\n%3\n )") + .arg(this->name()) + .arg(n) + .arg(parts.join("\n")); +} + +/** Return whether or not this is empty */ +bool AngleRestraints::isEmpty() const +{ + return this->r.isEmpty(); +} + +/** Return whether or not this is empty */ +bool AngleRestraints::isNull() const +{ + return this->isEmpty(); +} + +/** Return the number of restraints */ +int AngleRestraints::nRestraints() const +{ + return this->r.count(); +} + +/** Return the number of restraints */ +int AngleRestraints::count() const +{ + return this->nRestraints(); +} + +/** Return the number of restraints */ +int AngleRestraints::size() const +{ + return this->nRestraints(); +} + +/** Return the ith restraint */ +const AngleRestraint &AngleRestraints::at(int i) const +{ + i = SireID::Index(i).map(this->r.count()); + + return this->r.at(i); +} + +/** Return the ith restraint */ +const AngleRestraint &AngleRestraints::operator[](int i) const +{ + return this->at(i); +} + +/** Return all of the restraints */ +QList AngleRestraints::restraints() const +{ + return this->r; +} + +/** Add a restraints onto the list */ +void AngleRestraints::add(const AngleRestraint &restraint) +{ + if (not restraint.isNull()) + r.append(restraint); +} + +/** Add a restraint onto the list */ +void AngleRestraints::add(const AngleRestraints &restraints) +{ + this->r += restraints.r; +} + +/** Add a restraint onto the list */ +AngleRestraints &AngleRestraints::operator+=(const AngleRestraint &restraint) +{ + this->add(restraint); + return *this; +} + +/** Add a restraint onto the list */ +AngleRestraints AngleRestraints::operator+(const AngleRestraint &restraint) const +{ + AngleRestraints ret(*this); + ret += restraint; + return *this; +} + +/** Add restraints onto the list */ +AngleRestraints &AngleRestraints::operator+=(const AngleRestraints &restraints) +{ + this->add(restraints); + return *this; +} + +/** Add restraints onto the list */ +AngleRestraints AngleRestraints::operator+(const AngleRestraints &restraints) const +{ + AngleRestraints ret(*this); + ret += restraints; + return *this; +} diff --git a/corelib/src/libs/SireMM/anglerestraints.h b/corelib/src/libs/SireMM/anglerestraints.h new file mode 100644 index 000000000..b34de4822 --- /dev/null +++ b/corelib/src/libs/SireMM/anglerestraints.h @@ -0,0 +1,177 @@ +/********************************************\ + * + * Sire - Molecular Simulation Framework + * + * Copyright (C) 2025 Christopher Woods + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * For full details of the license please see the COPYING file + * that should have come with this distribution. + * + * You can contact the authors at https://sire.openbiosim.org + * +\*********************************************/ + +#ifndef SIREMM_ANGLERESTRAINTS_H +#define SIREMM_ANGLERESTRAINTS_H + + +#include "restraints.h" +#include "SireUnits/dimensions.h" +#include "SireUnits/generalunit.h" + +SIRE_BEGIN_HEADER + +namespace SireMM +{ + class AngleRestraint; + class AngleRestraints; +} + +SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::AngleRestraint &); +SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::AngleRestraint &); + +SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::AngleRestraints &); +SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::AngleRestraints &); + +namespace SireMM +{ + + /** This class represents a single angle restraint between any three + * atoms in a system + * @author Christopher Woods + */ + class SIREMM_EXPORT AngleRestraint + : public SireBase::ConcreteProperty + { + + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::AngleRestraint &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::AngleRestraint &); + + public: + AngleRestraint(); + AngleRestraint(const QList &atoms, + const SireUnits::Dimension::Angle &theta0, + const SireUnits::Dimension::HarmonicAngleConstant &ktheta); + + AngleRestraint(const AngleRestraint &other); + + ~AngleRestraint(); + + AngleRestraint &operator=(const AngleRestraint &other); + + bool operator==(const AngleRestraint &other) const; + bool operator!=(const AngleRestraint &other) const; + + AngleRestraints operator+(const AngleRestraint &other) const; + AngleRestraints operator+(const AngleRestraints &other) const; + + static const char *typeName(); + const char *what() const; + + AngleRestraint *clone() const; + + QString toString() const; + + bool isNull() const; + + QVector atoms() const; + + SireUnits::Dimension::Angle theta0() const; + SireUnits::Dimension::HarmonicAngleConstant ktheta() const; + + private: + /** Atoms involved in the angle restraint */ + QVector atms; + + /** Equilibrium angle */ + SireUnits::Dimension::Angle _theta0; + + /** Harmonic angle constant */ + SireUnits::Dimension::HarmonicAngleConstant _ktheta; + }; + + /** This class represents a collection of angle restraints */ + class SIREMM_EXPORT AngleRestraints + : public SireBase::ConcreteProperty + { + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::AngleRestraints &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::AngleRestraints &); + + public: + AngleRestraints(); + AngleRestraints(const QString &name); + + AngleRestraints(const AngleRestraint &restraint); + AngleRestraints(const QList &restraints); + + AngleRestraints(const QString &name, + const AngleRestraint &restraint); + + AngleRestraints(const QString &name, + const QList &restraints); + + AngleRestraints(const AngleRestraints &other); + + ~AngleRestraints(); + + AngleRestraints &operator=(const AngleRestraints &other); + + bool operator==(const AngleRestraints &other) const; + bool operator!=(const AngleRestraints &other) const; + + static const char *typeName(); + const char *what() const; + + AngleRestraints *clone() const; + + QString toString() const; + + bool isEmpty() const; + bool isNull() const; + + int count() const; + int size() const; + int nRestraints() const; + + const AngleRestraint &at(int i) const; + const AngleRestraint &operator[](int i) const; + + QList restraints() const; + + void add(const AngleRestraint &restraint); + void add(const AngleRestraints &restraints); + + AngleRestraints &operator+=(const AngleRestraint &restraint); + AngleRestraints &operator+=(const AngleRestraints &restraints); + + AngleRestraints operator+(const AngleRestraint &restraint) const; + AngleRestraints operator+(const AngleRestraints &restraints) const; + + private: + /** List of restraints */ + QList r; + }; +} + +Q_DECLARE_METATYPE(SireMM::AngleRestraint) +Q_DECLARE_METATYPE(SireMM::AngleRestraints) + +SIRE_EXPOSE_CLASS(SireMM::AngleRestraint) +SIRE_EXPOSE_CLASS(SireMM::AngleRestraints) +SIRE_END_HEADER + +#endif diff --git a/corelib/src/libs/SireMM/dihedralrestraint.cpp b/corelib/src/libs/SireMM/dihedralrestraint.cpp deleted file mode 100644 index bea13ed94..000000000 --- a/corelib/src/libs/SireMM/dihedralrestraint.cpp +++ /dev/null @@ -1,576 +0,0 @@ -/********************************************\ - * - * Sire - Molecular Simulation Framework - * - * Copyright (C) 2009 Christopher Woods - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * - * For full details of the license please see the COPYING file - * that should have come with this distribution. - * - * You can contact the authors at https://sire.openbiosim.org - * -\*********************************************/ - -#include "dihedralrestraint.h" - -#include "SireFF/forcetable.h" - -#include "SireCAS/conditional.h" -#include "SireCAS/power.h" -#include "SireCAS/symbols.h" -#include "SireCAS/values.h" - -#include "SireID/index.h" - -#include "SireUnits/angle.h" -#include "SireUnits/units.h" - -#include "SireStream/datastream.h" -#include "SireStream/shareddatastream.h" - -#include "SireCAS/errors.h" - -using namespace SireMM; -using namespace SireMol; -using namespace SireFF; -using namespace SireID; -using namespace SireBase; -using namespace SireCAS; -using namespace SireMaths; -using namespace SireStream; -using namespace SireUnits; -using namespace SireUnits::Dimension; - -//////////// -//////////// Implementation of DihedralRestraint -//////////// - -static const RegisterMetaType r_dihrest; - -/** Serialise to a binary datastream */ -QDataStream &operator<<(QDataStream &ds, const DihedralRestraint &dihrest) -{ - writeHeader(ds, r_dihrest, 1); - - SharedDataStream sds(ds); - - sds << dihrest.p[0] << dihrest.p[1] << dihrest.p[2] << dihrest.p[3] << dihrest.force_expression - << static_cast(dihrest); - - return ds; -} - -/** Extract from a binary datastream */ -QDataStream &operator>>(QDataStream &ds, DihedralRestraint &dihrest) -{ - VersionID v = readHeader(ds, r_dihrest); - - if (v == 1) - { - SharedDataStream sds(ds); - - sds >> dihrest.p[0] >> dihrest.p[1] >> dihrest.p[2] >> dihrest.p[3] >> dihrest.force_expression >> - static_cast(dihrest); - - dihrest.intra_molecule_points = Point::areIntraMoleculePoints(dihrest.p[0], dihrest.p[1]) and - Point::areIntraMoleculePoints(dihrest.p[0], dihrest.p[2]) and - Point::areIntraMoleculePoints(dihrest.p[0], dihrest.p[3]); - } - else - throw version_error(v, "1", r_dihrest, CODELOC); - - return ds; -} - -Q_GLOBAL_STATIC_WITH_ARGS(Symbol, getPhiSymbol, ("phi")); - -/** Return the symbol that represents the dihedral angle between the points (phi) */ -const Symbol &DihedralRestraint::phi() -{ - return *(getPhiSymbol()); -} - -/** Constructor */ -DihedralRestraint::DihedralRestraint() : ConcreteProperty() -{ -} - -void DihedralRestraint::calculatePhi() -{ - if (this->restraintFunction().isFunction(phi())) - { - SireUnits::Dimension::Angle angle; - - if (intra_molecule_points) - // we don't use the space when calculating intra-molecular angles - angle = - Vector::dihedral(p[0].read().point(), p[1].read().point(), p[2].read().point(), p[3].read().point()); - else - angle = this->space().calcDihedral(p[0].read().point(), p[1].read().point(), p[2].read().point(), - p[3].read().point()); - - ExpressionRestraint3D::_pvt_setValue(phi(), angle); - } -} - -/** Construct a restraint that acts on the angle within the - three points 'point0', 'point1' and 'point2' (theta == a(012)), - restraining the angle within these points using the expression - 'restraint' */ -DihedralRestraint::DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const Expression &restraint) - : ConcreteProperty(restraint) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - p[3] = point3; - - force_expression = this->restraintFunction().differentiate(phi()); - - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]) and - Point::areIntraMoleculePoints(p[0], p[3]); - - this->calculatePhi(); -} - -/** Construct a restraint that acts on the angle within the - three points 'point0', 'point1' and 'point2' (theta == a(012)), - restraining the angle within these points using the expression - 'restraint' */ -DihedralRestraint::DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const Expression &restraint, const Values &values) - : ConcreteProperty(restraint, values) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - p[3] = point3; - - force_expression = this->restraintFunction().differentiate(phi()); - - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]) and - Point::areIntraMoleculePoints(p[0], p[3]); - - this->calculatePhi(); -} - -/** Internal constructor used to construct a restraint using the specified - points, energy expression and force expression */ -DihedralRestraint::DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const Expression &nrg_restraint, - const Expression &force_restraint) - : ConcreteProperty(nrg_restraint), force_expression(force_restraint) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - p[3] = point3; - - if (force_expression.isConstant()) - { - force_expression = force_expression.evaluate(Values()); - } - else - { - if (not this->restraintFunction().symbols().contains(force_expression.symbols())) - throw SireError::incompatible_error( - QObject::tr("You cannot use a force function which uses more symbols " - "(%1) than the energy function (%2).") - .arg(Sire::toString(force_expression.symbols()), Sire::toString(restraintFunction().symbols())), - CODELOC); - } - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]) and - Point::areIntraMoleculePoints(p[0], p[3]); - - this->calculatePhi(); -} - -/** Copy constructor */ -DihedralRestraint::DihedralRestraint(const DihedralRestraint &other) - : ConcreteProperty(other), force_expression(other.force_expression), - intra_molecule_points(other.intra_molecule_points) -{ - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } -} - -/** Destructor */ -DihedralRestraint::~DihedralRestraint() -{ -} - -/** Copy assignment operator */ -DihedralRestraint &DihedralRestraint::operator=(const DihedralRestraint &other) -{ - if (this != &other) - { - ExpressionRestraint3D::operator=(other); - - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } - - force_expression = other.force_expression; - intra_molecule_points = other.intra_molecule_points; - } - - return *this; -} - -/** Comparison operator */ -bool DihedralRestraint::operator==(const DihedralRestraint &other) const -{ - return this == &other or (ExpressionRestraint3D::operator==(other) and p[0] == other.p[0] and p[1] == other.p[1] and - p[2] == other.p[2] and p[3] == other.p[3] and force_expression == other.force_expression); -} - -/** Comparison operator */ -bool DihedralRestraint::operator!=(const DihedralRestraint &other) const -{ - return not DihedralRestraint::operator==(other); -} - -const char *DihedralRestraint::typeName() -{ - return QMetaType::typeName(qMetaTypeId()); -} - -/** This restraint involves four points */ -int DihedralRestraint::nPoints() const -{ - return 4; -} - -/** Return the ith point */ -const Point &DihedralRestraint::point(int i) const -{ - i = Index(i).map(this->nPoints()); - - return p[i].read(); -} - -/** Return the first point */ -const Point &DihedralRestraint::point0() const -{ - return p[0].read(); -} - -/** Return the second point */ -const Point &DihedralRestraint::point1() const -{ - return p[1].read(); -} - -/** Return the third point */ -const Point &DihedralRestraint::point2() const -{ - return p[2].read(); -} - -/** Return the fourth point */ -const Point &DihedralRestraint::point3() const -{ - return p[3].read(); -} - -/** Return the built-in symbols for this restraint */ -Symbols DihedralRestraint::builtinSymbols() const -{ - if (this->restraintFunction().isFunction(phi())) - return phi(); - else - return Symbols(); -} - -/** Return the values of the built-in symbols of this restraint */ -Values DihedralRestraint::builtinValues() const -{ - if (this->restraintFunction().isFunction(phi())) - return phi() == this->values()[phi()]; - else - return Values(); -} - -/** Return the differential of this restraint with respect to - the symbol 'symbol' - - \throw SireCAS::unavailable_differential -*/ -RestraintPtr DihedralRestraint::differentiate(const Symbol &symbol) const -{ - if (this->restraintFunction().isFunction(symbol)) - return DihedralRestraint(p[0], p[1], p[2], p[3], restraintFunction().differentiate(symbol), this->values()); - else - return NullRestraint(); -} - -/** Set the space used to evaluate the energy of this restraint - - \throw SireVol::incompatible_space -*/ -void DihedralRestraint::setSpace(const Space &new_space) -{ - if (not this->space().equals(new_space)) - { - DihedralRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().setSpace(new_space); - } - - Restraint3D::setSpace(new_space); - - this->calculatePhi(); - } - catch (...) - { - DihedralRestraint::operator=(old_state); - throw; - } - } -} - -/** Return the function used to calculate the restraint force */ -const Expression &DihedralRestraint::differentialRestraintFunction() const -{ - return force_expression; -} - -/** Calculate the force acting on the molecule in the forcetable 'forcetable' - caused by this restraint, and add it on to the forcetable scaled by - 'scale_force' */ -void DihedralRestraint::force(MolForceTable &forcetable, double scale_force) const -{ - bool in_p0 = p[0].read().contains(forcetable.molNum()); - bool in_p1 = p[1].read().contains(forcetable.molNum()); - bool in_p2 = p[2].read().contains(forcetable.molNum()); - bool in_p3 = p[3].read().contains(forcetable.molNum()); - - if (not(in_p0 or in_p1 or in_p2 or in_p3)) - // this molecule is not affected by the restraint - return; - - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by a dihedral restraint."), - CODELOC); -} - -/** Calculate the force acting on the molecules in the forcetable 'forcetable' - caused by this restraint, and add it on to the forcetable scaled by - 'scale_force' */ -void DihedralRestraint::force(ForceTable &forcetable, double scale_force) const -{ - bool in_p0 = p[0].read().usesMoleculesIn(forcetable); - bool in_p1 = p[1].read().usesMoleculesIn(forcetable); - bool in_p2 = p[2].read().usesMoleculesIn(forcetable); - bool in_p3 = p[3].read().usesMoleculesIn(forcetable); - - if (not(in_p0 or in_p1 or in_p2 or in_p3)) - // this molecule is not affected by the restraint - return; - - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by a dihedral restraint."), - CODELOC); -} - -/** Update the points of this restraint using new molecule data from 'moldata' - - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void DihedralRestraint::update(const MoleculeData &moldata) -{ - if (this->contains(moldata.number())) - { - DihedralRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(moldata); - } - - this->calculatePhi(); - } - catch (...) - { - DihedralRestraint::operator=(old_state); - throw; - } - } -} - -/** Update the points of this restraint using new molecule data from 'molecules' - - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void DihedralRestraint::update(const Molecules &molecules) -{ - if (this->usesMoleculesIn(molecules)) - { - DihedralRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(molecules); - } - - this->calculatePhi(); - } - catch (...) - { - DihedralRestraint::operator=(old_state); - throw; - } - } -} - -/** Return the molecules used in this restraint */ -Molecules DihedralRestraint::molecules() const -{ - Molecules mols; - - for (int i = 0; i < this->nPoints(); ++i) - { - mols += p[i].read().molecules(); - } - - return mols; -} - -/** Return whether or not this restraint affects the molecule - with number 'molnum' */ -bool DihedralRestraint::contains(MolNum molnum) const -{ - return p[0].read().contains(molnum) or p[1].read().contains(molnum) or p[2].read().contains(molnum) or - p[3].read().contains(molnum); -} - -/** Return whether or not this restraint affects the molecule - with ID 'molid' */ -bool DihedralRestraint::contains(const MolID &molid) const -{ - return p[0].read().contains(molid) or p[1].read().contains(molid) or p[2].read().contains(molid) or - p[3].read().contains(molid); -} - -/** Return whether or not this restraint involves any of the molecules - that are in the forcetable 'forcetable' */ -bool DihedralRestraint::usesMoleculesIn(const ForceTable &forcetable) const -{ - return p[0].read().usesMoleculesIn(forcetable) or p[1].read().usesMoleculesIn(forcetable) or - p[2].read().usesMoleculesIn(forcetable) or p[3].read().usesMoleculesIn(forcetable); -} - -/** Return whether or not this restraint involves any of the molecules - in 'molecules' */ -bool DihedralRestraint::usesMoleculesIn(const Molecules &molecules) const -{ - return p[0].read().usesMoleculesIn(molecules) or p[1].read().usesMoleculesIn(molecules) or - p[2].read().usesMoleculesIn(molecules) or p[3].read().usesMoleculesIn(molecules); -} - -static Expression harmonicFunction(double force_constant) -{ - if (SireMaths::isZero(force_constant)) - return 0; - else - return force_constant * pow(DihedralRestraint::phi(), 2); -} - -static Expression diffHarmonicFunction(double force_constant) -{ - if (SireMaths::isZero(force_constant)) - return 0; - else - return (2 * force_constant) * DihedralRestraint::phi(); -} - -/** Return a distance restraint that applies a harmonic potential between - the points 'point0' and 'point1' using a force constant 'force_constant' */ -DihedralRestraint DihedralRestraint::harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const HarmonicAngleForceConstant &force_constant) -{ - return DihedralRestraint(point0, point1, point2, point3, ::harmonicFunction(force_constant), - ::diffHarmonicFunction(force_constant)); -} - -static Expression halfHarmonicFunction(double force_constant, double angle) -{ - if (SireMaths::isZero(force_constant)) - return 0; - - else if (angle <= 0) - // this is just a harmonic function - return ::harmonicFunction(force_constant); - - else - { - const Symbol &phi = DihedralRestraint::phi(); - return Conditional(GreaterThan(phi, angle), force_constant * pow(phi - angle, 2), 0); - } -} - -static Expression diffHalfHarmonicFunction(double force_constant, double angle) -{ - if (SireMaths::isZero(force_constant)) - return 0; - - else if (angle <= 0) - // this is just a harmonic function - return ::diffHarmonicFunction(force_constant); - - else - { - const Symbol &phi = DihedralRestraint::phi(); - return Conditional(GreaterThan(phi, angle), (2 * force_constant) * (phi - angle), 0); - } -} - -/** Return a distance restraint that applied a half-harmonic potential - between the points 'point0' and 'point1' above a distance 'distance' - using a force constant 'force_constant' */ -DihedralRestraint DihedralRestraint::halfHarmonic(const PointRef &point0, const PointRef &point1, - const PointRef &point2, const PointRef &point3, const Angle &angle, - const HarmonicAngleForceConstant &force_constant) -{ - double ang = angle.to(radians); - - return DihedralRestraint(point0, point1, point2, point3, ::halfHarmonicFunction(force_constant, ang), - ::diffHalfHarmonicFunction(force_constant, ang)); -} diff --git a/corelib/src/libs/SireMM/dihedralrestraint.h b/corelib/src/libs/SireMM/dihedralrestraint.h deleted file mode 100644 index 606ee3b78..000000000 --- a/corelib/src/libs/SireMM/dihedralrestraint.h +++ /dev/null @@ -1,144 +0,0 @@ -/********************************************\ - * - * Sire - Molecular Simulation Framework - * - * Copyright (C) 2009 Christopher Woods - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * - * For full details of the license please see the COPYING file - * that should have come with this distribution. - * - * You can contact the authors at https://sire.openbiosim.org - * -\*********************************************/ - -#ifndef SIREMM_DIHEDRALRESTRAINT_H -#define SIREMM_DIHEDRALRESTRAINT_H - -#include "anglerestraint.h" -#include "restraint.h" - -SIRE_BEGIN_HEADER - -namespace SireMM -{ - class DihedralRestraint; -} - -SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::DihedralRestraint &); -SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::DihedralRestraint &); - -namespace SireMM -{ - - /** This is a restraint that operates on the dihedral angle between - four SireMM::Point objects (e.g. four atoms in a molecule) - - @author Christopher Woods - */ - class SIREMM_EXPORT DihedralRestraint : public SireBase::ConcreteProperty - { - - friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const DihedralRestraint &); - friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, DihedralRestraint &); - - public: - DihedralRestraint(); - - DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const PointRef &point3, - const Expression &restraint); - - DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const PointRef &point3, - const Expression &restraint, const Values &values); - - DihedralRestraint(const DihedralRestraint &other); - - ~DihedralRestraint(); - - DihedralRestraint &operator=(const DihedralRestraint &other); - - bool operator==(const DihedralRestraint &other) const; - bool operator!=(const DihedralRestraint &other) const; - - static const char *typeName(); - - const Point &point(int i) const; - - const Point &point0() const; - const Point &point1() const; - const Point &point2() const; - const Point &point3() const; - - int nPoints() const; - - static const Symbol &phi(); - - Symbols builtinSymbols() const; - Values builtinValues() const; - - RestraintPtr differentiate(const Symbol &symbol) const; - - void setSpace(const Space &space); - - const Expression &differentialRestraintFunction() const; - - void force(MolForceTable &forcetable, double scale_force = 1) const; - void force(ForceTable &forcetable, double scale_force = 1) const; - - void update(const MoleculeData &moldata); - void update(const Molecules &molecules); - - Molecules molecules() const; - - bool contains(MolNum molnum) const; - bool contains(const MolID &molid) const; - - bool usesMoleculesIn(const ForceTable &forcetable) const; - bool usesMoleculesIn(const Molecules &molecules) const; - - static DihedralRestraint harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const HarmonicAngleForceConstant &force_constant); - - static DihedralRestraint halfHarmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const SireUnits::Dimension::Angle &angle, - const HarmonicAngleForceConstant &force_constant); - - protected: - DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const PointRef &point3, - const Expression &nrg_restraint, const Expression &force_restraint); - - private: - void calculatePhi(); - - /** The four points between which the restraint is calculated */ - SireFF::PointPtr p[4]; - - /** The expression used to calculate the force */ - Expression force_expression; - - /** Whether or not all four points are within the same molecule */ - bool intra_molecule_points; - }; - -} // namespace SireMM - -Q_DECLARE_METATYPE(SireMM::DihedralRestraint) - -SIRE_EXPOSE_CLASS(SireMM::DihedralRestraint) - -SIRE_END_HEADER - -#endif diff --git a/corelib/src/libs/SireMM/dihedralrestraints.cpp b/corelib/src/libs/SireMM/dihedralrestraints.cpp new file mode 100644 index 000000000..70054b34c --- /dev/null +++ b/corelib/src/libs/SireMM/dihedralrestraints.cpp @@ -0,0 +1,477 @@ +/********************************************\ + * + * Sire - Molecular Simulation Framework + * + * Copyright (C) 2025 Christopher Woods + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * For full details of the license please see the COPYING file + * that should have come with this distribution. + * + * You can contact the authors at https://sire.openbiosim.org + * +\*********************************************/ + +#include "dihedralrestraints.h" + +#include "SireID/index.h" +#include "SireStream/datastream.h" +#include "SireStream/shareddatastream.h" +#include "SireUnits/units.h" + +#include "SireCAS/errors.h" + +#include + +using namespace SireMM; +using namespace SireID; +using namespace SireBase; +using namespace SireMaths; +using namespace SireStream; +using namespace SireUnits; +using namespace SireUnits::Dimension; + +//////////// +//////////// Implementation of DihedralRestraint +//////////// + +static const RegisterMetaType r_dihrest; + +/** Serialise to a binary datastream */ + +QDataStream &operator<<(QDataStream &ds, const DihedralRestraint &dihrest) +{ + writeHeader(ds, r_dihrest, 1); + + SharedDataStream sds(ds); + + sds << dihrest.atms << dihrest._phi0 << dihrest._kphi; + + return ds; +} + +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, DihedralRestraint &dihrest) +{ + VersionID v = readHeader(ds, r_dihrest); + + if (v == 1) + { + SharedDataStream sds(ds); + + sds >> dihrest.atms >> dihrest._phi0 >> dihrest._kphi; + } + else + throw version_error(v, "1", r_dihrest, CODELOC); + + return ds; +} + +/** Null constructor */ +DihedralRestraint::DihedralRestraint() + : ConcreteProperty(), + _phi0(0), _kphi(0) +{ +} + +/** Construct a restraint that acts on the angle within the + four atoms 'atom0', 'atom1', 'atom2' 'atom3' (phi == a(0123)), + restraining the angle within these atoms */ +DihedralRestraint::DihedralRestraint(const QList &atoms, + const SireUnits::Dimension::Angle &phi0, + const SireUnits::Dimension::HarmonicAngleConstant &kphi) + : ConcreteProperty(), + _phi0(phi0), _kphi(kphi) +{ + + // Make sure that we have 4 distinct atoms + QSet distinct; + distinct.reserve(4); + + for (const auto &atom : atoms) + { + if (atom >= 0) + distinct.insert(atom); + } + + atms = atoms.toVector(); +} + +/* Copy constructor*/ +DihedralRestraint::DihedralRestraint(const DihedralRestraint &other) + : ConcreteProperty(other), + atms(other.atms), _phi0(other._phi0), _kphi(other._kphi) + +{ +} + +/* Destructor */ +DihedralRestraint::~DihedralRestraint() +{ +} + +DihedralRestraint &DihedralRestraint::operator=(const DihedralRestraint &other) +{ + if (this != &other) + { + Property::operator=(other); + atms = other.atms; + _phi0 = other._phi0; + _kphi = other._kphi; + } + + return *this; +} + +bool DihedralRestraint::operator==(const DihedralRestraint &other) const +{ + return atms == other.atms and + _phi0 == other._phi0 and + _kphi == other._kphi; +} + +bool DihedralRestraint::operator!=(const DihedralRestraint &other) const +{ + return not operator==(other); +} + +DihedralRestraints DihedralRestraint::operator+(const DihedralRestraint &other) const +{ + return DihedralRestraints(*this) + other; +} + +DihedralRestraints DihedralRestraint::operator+(const DihedralRestraints &other) const +{ + return DihedralRestraints(*this) + other; +} + +const char *DihedralRestraint::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); +} + +const char *DihedralRestraint::what() const +{ + return DihedralRestraint::typeName(); +} + +DihedralRestraint *DihedralRestraint::clone() const +{ + return new DihedralRestraint(*this); +} + +bool DihedralRestraint::isNull() const +{ + return atms.isEmpty(); +} + +QString DihedralRestraint::toString() const +{ + if (this->isNull()) + return QObject::tr("DihedralRestraint::null"); + else + { + QStringList a; + + for (const auto &atom : atms) + { + a.append(QString::number(atom)); + } + return QString("DihedralRestraint( [%1], phi0=%2, kphi=%3 )") + .arg(a.join(", ")) + .arg(_phi0.toString()) + .arg(_kphi.toString()); + } +} + +/** Return the force constant for the restraint */ +SireUnits::Dimension::HarmonicAngleConstant DihedralRestraint::kphi() const +{ + return this->_kphi; +} + +/** Return the equilibrium angle for the restraint */ +SireUnits::Dimension::Angle DihedralRestraint::phi0() const +{ + return this->_phi0; +} + +/** Return the atoms involved in the restraint */ +QVector DihedralRestraint::atoms() const +{ + return this->atms; +} + +/////// +/////// Implementation of DihedralRestraints +/////// + +/** Serialise to a binary datastream */ + +static const RegisterMetaType r_dihrests; + +QDataStream &operator<<(QDataStream &ds, const DihedralRestraints &dihrests) +{ + writeHeader(ds, r_dihrests, 1); + + SharedDataStream sds(ds); + + sds << dihrests.r + << static_cast(dihrests); + + return ds; +} + +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, DihedralRestraints &dihrests) +{ + VersionID v = readHeader(ds, r_dihrests); + + if (v == 1) + { + SharedDataStream sds(ds); + + sds >> dihrests.r >> + static_cast(dihrests); + } + else + throw version_error(v, "1", r_dihrests, CODELOC); + + return ds; +} + +/** Null constructor */ +DihedralRestraints::DihedralRestraints() + : ConcreteProperty() +{ +} + +DihedralRestraints::DihedralRestraints(const QString &name) + : ConcreteProperty(name) +{ +} + +DihedralRestraints::DihedralRestraints(const DihedralRestraint &restraint) + : ConcreteProperty() +{ + if (not restraint.isNull()) + r.append(restraint); +} + +DihedralRestraints::DihedralRestraints(const QList &restraints) + : ConcreteProperty() +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } +} + +DihedralRestraints::DihedralRestraints(const QString &name, + const DihedralRestraint &restraint) + : ConcreteProperty(name) +{ + if (not restraint.isNull()) + r.append(restraint); +} + +DihedralRestraints::DihedralRestraints(const QString &name, + const QList &restraints) + : ConcreteProperty(name) +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } +} + +DihedralRestraints::DihedralRestraints(const DihedralRestraints &other) + : ConcreteProperty(other), r(other.r) +{ +} + +/* Desctructor */ +DihedralRestraints::~DihedralRestraints() +{ +} + +DihedralRestraints &DihedralRestraints::operator=(const DihedralRestraints &other) +{ + if (this != &other) + { + Restraints::operator=(other); + r = other.r; + } + + return *this; +} + +bool DihedralRestraints::operator==(const DihedralRestraints &other) const +{ + return r == other.r and Restraints::operator==(other); +} + +bool DihedralRestraints::operator!=(const DihedralRestraints &other) const +{ + return not operator==(other); +} + +const char *DihedralRestraints::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); +} + +const char *DihedralRestraints::what() const +{ + return DihedralRestraints::typeName(); +} + +DihedralRestraints *DihedralRestraints::clone() const +{ + return new DihedralRestraints(*this); +} + +QString DihedralRestraints::toString() const +{ + if (this->isEmpty()) + return QObject::tr("DihedralRestraints::null"); + + QStringList parts; + + const auto n = this->count(); + + if (n <= 10) + { + for (int i = 0; i < n; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + } + else + { + for (int i = 0; i < 5; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + + parts.append("..."); + + for (int i = n - 5; i < n; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + } + + return QObject::tr("DihedralRestraints( name=%1, size=%2\n%3\n )") + .arg(this->name()) + .arg(n) + .arg(parts.join("\n")); +} + +/** Return whether or not this is empty */ +bool DihedralRestraints::isEmpty() const +{ + return this->r.isEmpty(); +} + +/** Return whether or not this is empty */ +bool DihedralRestraints::isNull() const +{ + return this->isEmpty(); +} + +/** Return the number of restraints */ +int DihedralRestraints::nRestraints() const +{ + return this->r.count(); +} + +/** Return the number of restraints */ +int DihedralRestraints::count() const +{ + return this->nRestraints(); +} + +/** Return the number of restraints */ +int DihedralRestraints::size() const +{ + return this->nRestraints(); +} + +/** Return the ith restraint */ +const DihedralRestraint &DihedralRestraints::at(int i) const +{ + i = SireID::Index(i).map(this->r.count()); + + return this->r.at(i); +} + +/** Return the ith restraint */ +const DihedralRestraint &DihedralRestraints::operator[](int i) const +{ + return this->at(i); +} + +/** Return all of the restraints */ +QList DihedralRestraints::restraints() const +{ + return this->r; +} + +/** Add a restraints onto the list */ +void DihedralRestraints::add(const DihedralRestraint &restraint) +{ + if (not restraint.isNull()) + r.append(restraint); +} + +/** Add a restraint onto the list */ +void DihedralRestraints::add(const DihedralRestraints &restraints) +{ + this->r += restraints.r; +} + +/** Add a restraint onto the list */ +DihedralRestraints &DihedralRestraints::operator+=(const DihedralRestraint &restraint) +{ + this->add(restraint); + return *this; +} + +/** Add a restraint onto the list */ +DihedralRestraints DihedralRestraints::operator+(const DihedralRestraint &restraint) const +{ + DihedralRestraints ret(*this); + ret += restraint; + return *this; +} + +/** Add restraints onto the list */ +DihedralRestraints &DihedralRestraints::operator+=(const DihedralRestraints &restraints) +{ + this->add(restraints); + return *this; +} + +/** Add restraints onto the list */ +DihedralRestraints DihedralRestraints::operator+(const DihedralRestraints &restraints) const +{ + DihedralRestraints ret(*this); + ret += restraints; + return *this; +} \ No newline at end of file diff --git a/corelib/src/libs/SireMM/dihedralrestraints.h b/corelib/src/libs/SireMM/dihedralrestraints.h new file mode 100644 index 000000000..f40b03e0c --- /dev/null +++ b/corelib/src/libs/SireMM/dihedralrestraints.h @@ -0,0 +1,177 @@ +/********************************************\ + * + * Sire - Molecular Simulation Framework + * + * Copyright (C) 2025 Christopher Woods + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * For full details of the license please see the COPYING file + * that should have come with this distribution. + * + * You can contact the authors at https://sire.openbiosim.org + * +\*********************************************/ + +#ifndef SIREMM_DIHEDRALRESTRAINTS_H +#define SIREMM_DIHEDRALRESTRAINTS_H + +#include "restraints.h" + +#include "SireUnits/dimensions.h" +#include "SireUnits/generalunit.h" + +SIRE_BEGIN_HEADER + +namespace SireMM +{ + class DihedralRestraint; + class DihedralRestraints; +} + +SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::DihedralRestraint &); +SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::DihedralRestraint &); + +SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::DihedralRestraints &); +SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::DihedralRestraints &); + +namespace SireMM +{ + + /** This class represents a single torsion restraint between any four + * atoms in a system + * @author Christopher Woods + */ + class SIREMM_EXPORT DihedralRestraint + : public SireBase::ConcreteProperty + { + + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::DihedralRestraint &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::DihedralRestraint &); + + public: + DihedralRestraint(); + DihedralRestraint(const QList &atoms, + const SireUnits::Dimension::Angle &phi0, + const SireUnits::Dimension::HarmonicAngleConstant &kphi); + + DihedralRestraint(const DihedralRestraint &other); + + ~DihedralRestraint(); + + DihedralRestraint &operator=(const DihedralRestraint &other); + + bool operator==(const DihedralRestraint &other) const; + bool operator!=(const DihedralRestraint &other) const; + + DihedralRestraints operator+(const DihedralRestraint &other) const; + DihedralRestraints operator+(const DihedralRestraints &other) const; + + static const char *typeName(); + const char *what() const; + + DihedralRestraint *clone() const; + + QString toString() const; + + bool isNull() const; + + QVector atoms() const; + + SireUnits::Dimension::Angle phi0() const; + SireUnits::Dimension::HarmonicAngleConstant kphi() const; + + private: + /** Atoms involved in the angle restraint */ + QVector atms; + + /** Equilibrium angle */ + SireUnits::Dimension::Angle _phi0; + + /** Harmonic angle constant */ + SireUnits::Dimension::HarmonicAngleConstant _kphi; + }; + + /** This class represents a collection of angle restraints */ + class SIREMM_EXPORT DihedralRestraints + : public SireBase::ConcreteProperty + { + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::DihedralRestraints &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::DihedralRestraints &); + + public: + DihedralRestraints(); + DihedralRestraints(const QString &name); + + DihedralRestraints(const DihedralRestraint &restraint); + DihedralRestraints(const QList &restraints); + + DihedralRestraints(const QString &name, + const DihedralRestraint &restraint); + + DihedralRestraints(const QString &name, + const QList &restraints); + + DihedralRestraints(const DihedralRestraints &other); + + ~DihedralRestraints(); + + DihedralRestraints &operator=(const DihedralRestraints &other); + + bool operator==(const DihedralRestraints &other) const; + bool operator!=(const DihedralRestraints &other) const; + + static const char *typeName(); + const char *what() const; + + DihedralRestraints *clone() const; + + QString toString() const; + + bool isEmpty() const; + bool isNull() const; + + int count() const; + int size() const; + int nRestraints() const; + + const DihedralRestraint &at(int i) const; + const DihedralRestraint &operator[](int i) const; + + QList restraints() const; + + void add(const DihedralRestraint &restraint); + void add(const DihedralRestraints &restraints); + + DihedralRestraints &operator+=(const DihedralRestraint &restraint); + DihedralRestraints &operator+=(const DihedralRestraints &restraints); + + DihedralRestraints operator+(const DihedralRestraint &restraint) const; + DihedralRestraints operator+(const DihedralRestraints &restraints) const; + + private: + /** List of restraints */ + QList r; + }; +} + +Q_DECLARE_METATYPE(SireMM::DihedralRestraint) +Q_DECLARE_METATYPE(SireMM::DihedralRestraints) + +SIRE_EXPOSE_CLASS(SireMM::DihedralRestraint) +SIRE_EXPOSE_CLASS(SireMM::DihedralRestraints) +SIRE_END_HEADER + +#endif \ No newline at end of file diff --git a/corelib/src/libs/SireMol/stereochemistry.cpp b/corelib/src/libs/SireMol/stereochemistry.cpp index e2570b8ef..cb9228e83 100644 --- a/corelib/src/libs/SireMol/stereochemistry.cpp +++ b/corelib/src/libs/SireMol/stereochemistry.cpp @@ -75,6 +75,8 @@ Stereochemistry::Stereochemistry(const QString &str) : ConcretePropertystereo_type = 1; + else if (s == "cis or trans") + this->stereo_type = 3; else if (s == "down") this->stereo_type = 6; else if (s == "not stereo") @@ -90,14 +92,14 @@ Stereochemistry Stereochemistry::fromSDF(int value) { Stereochemistry ret; - if (value == 0 or value == 1 or value == -1 or value == 6) + if (value == 0 or value == 1 or value == -1 or value == 3 or value == 6) { ret.stereo_type = value; } else { throw SireError::invalid_arg(QObject::tr("Invalid stereo type '%1'. Should be an integer in " - "[-1, 0, 1, 6]") + "[-1, 0, 1, 3, 6]") .arg(value), CODELOC); } @@ -194,6 +196,8 @@ QString Stereochemistry::toString() const return "not stereo"; case 1: return "up"; + case 3: + return "cis or trans"; case 6: return "down"; case -1: @@ -220,6 +224,8 @@ int Stereochemistry::toSDF() const { case 1: return 1; + case 3: + return 3; case 6: return 6; default: @@ -256,6 +262,12 @@ Stereochemistry Stereochemistry::up() return Stereochemistry::fromSDF(1); } +/** Return a cis or trans Stereochemistry */ +Stereochemistry Stereochemistry::cisOrTrans() +{ + return Stereochemistry::fromSDF(3); +} + /** Return a down Stereochemistry */ Stereochemistry Stereochemistry::down() { @@ -286,6 +298,12 @@ bool Stereochemistry::isUp() const return this->stereo_type == 1; } +/** Return whether or not this is a cis or trans bond */ +bool Stereochemistry::isCisOrTrans() const +{ + return this->stereo_type == 3; +} + /** Return whether or not this is a down bond */ bool Stereochemistry::isDown() const { diff --git a/corelib/src/libs/SireMol/stereochemistry.h b/corelib/src/libs/SireMol/stereochemistry.h index 0e45a81c9..60d368482 100644 --- a/corelib/src/libs/SireMol/stereochemistry.h +++ b/corelib/src/libs/SireMol/stereochemistry.h @@ -73,6 +73,7 @@ namespace SireMol static const char *typeName(); static Stereochemistry up(); + static Stereochemistry cisOrTrans(); static Stereochemistry down(); static Stereochemistry notStereo(); static Stereochemistry undefined(); @@ -85,6 +86,7 @@ namespace SireMol bool isDefined() const; bool isUp() const; + bool isCisOrTrans() const; bool isDown() const; bool isNotStereo() const; diff --git a/corelib/src/libs/SireMove/openmmfrenergydt.cpp b/corelib/src/libs/SireMove/openmmfrenergydt.cpp index 8d2cee119..9a8695d1b 100644 --- a/corelib/src/libs/SireMove/openmmfrenergydt.cpp +++ b/corelib/src/libs/SireMove/openmmfrenergydt.cpp @@ -129,7 +129,7 @@ QDataStream &operator<<(QDataStream &ds, const OpenMMFrEnergyDT &velver) sds << velver.frequent_save_velocities << velver.molgroup << velver.solutegroup << velver.CutoffType << velver.cutoff_distance << velver.field_dielectric << velver.Andersen_flag << velver.Andersen_frequency - << velver.MCBarostat_flag << velver.MCBarostat_frequency << velver.ConstraintType << velver.Pressure + << velver.MCBarostat_flag << velver.MCBarostat_membrane_flag << velver.MCBarostat_frequency << velver.ConstraintType << velver.Pressure << velver.Temperature << velver.platform_type << velver.Restraint_flag << velver.CMMremoval_frequency << velver.energy_frequency << velver.device_index << velver.Alchemical_value << velver.coulomb_power << velver.shift_delta << velver.delta_alchemical << velver.buffer_coords << velver.gradients @@ -151,7 +151,7 @@ QDataStream &operator>>(QDataStream &ds, OpenMMFrEnergyDT &velver) sds >> velver.frequent_save_velocities >> velver.molgroup >> velver.solutegroup >> velver.CutoffType >> velver.cutoff_distance >> velver.field_dielectric >> velver.Andersen_flag >> velver.Andersen_frequency >> - velver.MCBarostat_flag >> velver.MCBarostat_frequency >> velver.ConstraintType >> velver.Pressure >> + velver.MCBarostat_flag >> velver.MCBarostat_membrane_flag >> velver.MCBarostat_frequency >> velver.ConstraintType >> velver.Pressure >> velver.Temperature >> velver.platform_type >> velver.Restraint_flag >> velver.CMMremoval_frequency >> velver.energy_frequency >> velver.device_index >> velver.Alchemical_value >> velver.coulomb_power >> velver.shift_delta >> velver.delta_alchemical >> velver.buffer_coords >> velver.gradients >> @@ -175,7 +175,7 @@ OpenMMFrEnergyDT::OpenMMFrEnergyDT(bool frequent_save) : ConcreteProperty(), frequent_save_velocities(frequent_save), molgroup(MoleculeGroup()), solutegroup(MoleculeGroup()), openmm_system(0), isInitialised(false), CutoffType("nocutoff"), cutoff_distance(1.0 * nanometer), field_dielectric(78.3), Andersen_flag(false), - Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_frequency(25), ConstraintType("none"), + Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_membrane_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), Temperature(300.0 * kelvin), platform_type("Reference"), Restraint_flag(false), CMMremoval_frequency(0), energy_frequency(100), device_index("0"), Alchemical_value(0.5), coulomb_power(0), shift_delta(2.0), delta_alchemical(0.001), buffer_coords(false), gradients() @@ -188,7 +188,7 @@ OpenMMFrEnergyDT::OpenMMFrEnergyDT(const MoleculeGroup &molecule_group, const Mo : ConcreteProperty(), frequent_save_velocities(frequent_save), molgroup(molecule_group), solutegroup(solute_group), openmm_system(0), isInitialised(false), CutoffType("nocutoff"), cutoff_distance(1.0 * nanometer), field_dielectric(78.3), Andersen_flag(false), - Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_frequency(25), ConstraintType("none"), + Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_membrane_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), Temperature(300.0 * kelvin), platform_type("Reference"), Restraint_flag(false), CMMremoval_frequency(0), energy_frequency(100), device_index("0"), Alchemical_value(0.5), coulomb_power(0), shift_delta(2.0), delta_alchemical(0.001), buffer_coords(false), gradients() @@ -201,7 +201,7 @@ OpenMMFrEnergyDT::OpenMMFrEnergyDT(const OpenMMFrEnergyDT &other) molgroup(other.molgroup), solutegroup(other.solutegroup), openmm_system(other.openmm_system), isInitialised(other.isInitialised), CutoffType(other.CutoffType), cutoff_distance(other.cutoff_distance), field_dielectric(other.field_dielectric), Andersen_flag(other.Andersen_flag), - Andersen_frequency(other.Andersen_frequency), MCBarostat_flag(other.MCBarostat_flag), + Andersen_frequency(other.Andersen_frequency), MCBarostat_flag(other.MCBarostat_flag), MCBarostat_membrane_flag(other.MCBarostat_membrane_flag), MCBarostat_frequency(other.MCBarostat_frequency), ConstraintType(other.ConstraintType), Pressure(other.Pressure), Temperature(other.Temperature), platform_type(other.platform_type), Restraint_flag(other.Restraint_flag), CMMremoval_frequency(other.CMMremoval_frequency), energy_frequency(other.energy_frequency), @@ -232,6 +232,7 @@ OpenMMFrEnergyDT &OpenMMFrEnergyDT::operator=(const OpenMMFrEnergyDT &other) Andersen_flag = other.Andersen_flag; Andersen_frequency = other.Andersen_frequency; MCBarostat_flag = other.MCBarostat_flag; + MCBarostat_membrane_flag = other.MCBarostat_membrane_flag; MCBarostat_frequency = other.MCBarostat_frequency; ConstraintType = other.ConstraintType; Pressure = other.Pressure; @@ -521,9 +522,21 @@ void OpenMMFrEnergyDT::initialise() const double converted_Temperature = convertTo(Temperature.value(), kelvin); const double converted_Pressure = convertTo(Pressure.value(), bar); - OpenMM::MonteCarloBarostat *barostat = - new OpenMM::MonteCarloBarostat(converted_Pressure, converted_Temperature, MCBarostat_frequency); - system_openmm->addForce(barostat); + if (MCBarostat_membrane_flag) + { + // Simple options for now: zero surface tension, XY isotropic, Z free + const double surface_Tension = 0; + OpenMM::MonteCarloMembraneBarostat::XYMode xymode = OpenMM::MonteCarloMembraneBarostat::XYIsotropic; + OpenMM::MonteCarloMembraneBarostat::ZMode zmode = OpenMM::MonteCarloMembraneBarostat::ZFree; + OpenMM::MonteCarloMembraneBarostat * barostat = new OpenMM::MonteCarloMembraneBarostat(converted_Pressure, surface_Tension, converted_Temperature, xymode, zmode, MCBarostat_frequency); + system_openmm->addForce(barostat); + } + else + { + OpenMM::MonteCarloBarostat *barostat = + new OpenMM::MonteCarloBarostat(converted_Pressure, converted_Temperature, MCBarostat_frequency); + system_openmm->addForce(barostat); + } if (Debug) { @@ -531,6 +544,10 @@ void OpenMMFrEnergyDT::initialise() qDebug() << "Temperature = " << converted_Temperature << " K\n"; qDebug() << "Pressure = " << converted_Pressure << " bar\n"; qDebug() << "Frequency every " << MCBarostat_frequency << " steps\n"; + if (MCBarostat_membrane_flag) + { + qDebug() << "Membrane barostat, surface tension 0, XY isotropic, Z free\n"; + } } } @@ -1795,6 +1812,18 @@ bool OpenMMFrEnergyDT::getMCBarostat(void) return MCBarostat_flag; } +/** Set Monte Carlo membrane Barostat on/off */ +void OpenMMFrEnergyDT::setMCBarostat_membrane(bool MCBarostat_membrane) +{ + MCBarostat_membrane_flag = MCBarostat_membrane; +} + +bool OpenMMFrEnergyDT::getMCBarostat_membrane(void) +{ + + return MCBarostat_membrane_flag; +} + /** Get the Monte Carlo Barostat frequency in time speps */ int OpenMMFrEnergyDT::getMCBarostat_frequency(void) { diff --git a/corelib/src/libs/SireMove/openmmfrenergydt.h b/corelib/src/libs/SireMove/openmmfrenergydt.h index 88332e308..dded8a7da 100644 --- a/corelib/src/libs/SireMove/openmmfrenergydt.h +++ b/corelib/src/libs/SireMove/openmmfrenergydt.h @@ -114,6 +114,9 @@ namespace SireMove void setMCBarostat_frequency(int); int getMCBarostat_frequency(void); + bool getMCBarostat_membrane(void); + void setMCBarostat_membrane(bool); + QString getConstraintType(void); void setConstraintType(QString); @@ -180,6 +183,7 @@ namespace SireMove double Andersen_frequency; bool MCBarostat_flag; + bool MCBarostat_membrane_flag; int MCBarostat_frequency; QString ConstraintType; diff --git a/corelib/src/libs/SireMove/openmmfrenergyst.cpp b/corelib/src/libs/SireMove/openmmfrenergyst.cpp index ded8dcf8f..a27495e36 100644 --- a/corelib/src/libs/SireMove/openmmfrenergyst.cpp +++ b/corelib/src/libs/SireMove/openmmfrenergyst.cpp @@ -130,7 +130,7 @@ QDataStream &operator<<(QDataStream &ds, const OpenMMFrEnergyST &velver) sds << velver.frequent_save_velocities << velver.molgroup << velver.solute << velver.solutehard << velver.solutetodummy << velver.solutefromdummy << velver.combiningRules << velver.CutoffType << velver.cutoff_distance << velver.field_dielectric << velver.Andersen_flag << velver.Andersen_frequency - << velver.MCBarostat_flag << velver.MCBarostat_frequency << velver.ConstraintType << velver.Pressure + << velver.MCBarostat_flag << velver.MCBarostat_membrane_flag << velver.MCBarostat_frequency << velver.ConstraintType << velver.Pressure << velver.Temperature << velver.platform_type << velver.Restraint_flag << velver.CMMremoval_frequency << velver.buffer_frequency << velver.energy_frequency << velver.device_index << velver.precision << velver.Alchemical_value << velver.coulomb_power << velver.shift_delta << velver.delta_alchemical @@ -152,7 +152,7 @@ QDataStream &operator>>(QDataStream &ds, OpenMMFrEnergyST &velver) sds >> velver.frequent_save_velocities >> velver.molgroup >> velver.solute >> velver.solutehard >> velver.solutetodummy >> velver.solutefromdummy >> velver.combiningRules >> velver.CutoffType >> velver.cutoff_distance >> velver.field_dielectric >> velver.Andersen_flag >> velver.Andersen_frequency >> - velver.MCBarostat_flag >> velver.MCBarostat_frequency >> velver.ConstraintType >> velver.Pressure >> + velver.MCBarostat_flag >> velver.MCBarostat_membrane_flag >> velver.MCBarostat_frequency >> velver.ConstraintType >> velver.Pressure >> velver.Temperature >> velver.platform_type >> velver.Restraint_flag >> velver.CMMremoval_frequency >> velver.buffer_frequency >> velver.energy_frequency >> velver.device_index >> velver.precision >> velver.Alchemical_value >> velver.coulomb_power >> velver.shift_delta >> velver.delta_alchemical >> @@ -180,7 +180,7 @@ OpenMMFrEnergyST::OpenMMFrEnergyST(bool frequent_save) solutefromdummy(MoleculeGroup()), openmm_system(0), openmm_context(0), isSystemInitialised(false), isContextInitialised(false), combiningRules("arithmetic"), CutoffType("nocutoff"), cutoff_distance(1.0 * nanometer), field_dielectric(78.3), Andersen_flag(false), Andersen_frequency(90.0), - MCBarostat_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), + MCBarostat_flag(false), MCBarostat_membrane_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), Temperature(300.0 * kelvin), platform_type("Reference"), Restraint_flag(false), CMMremoval_frequency(0), buffer_frequency(0), energy_frequency(100), device_index("0"), precision("single"), Alchemical_value(0.5), coulomb_power(0), shift_delta(2.0), delta_alchemical(0.001), alchemical_array(), finite_diff_gradients(), @@ -199,7 +199,7 @@ OpenMMFrEnergyST::OpenMMFrEnergyST(const MoleculeGroup &molecule_group, const Mo solutefromdummy(solute_fromdummy), openmm_system(0), openmm_context(0), isSystemInitialised(false), isContextInitialised(false), combiningRules("arithmetic"), CutoffType("nocutoff"), cutoff_distance(1.0 * nanometer), field_dielectric(78.3), Andersen_flag(false), Andersen_frequency(90.0), - MCBarostat_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), + MCBarostat_flag(false), MCBarostat_membrane_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), Temperature(300.0 * kelvin), platform_type("Reference"), Restraint_flag(false), CMMremoval_frequency(0), buffer_frequency(0), energy_frequency(100), device_index("0"), precision("single"), Alchemical_value(0.5), coulomb_power(0), shift_delta(2.0), delta_alchemical(0.001), alchemical_array(), finite_diff_gradients(), @@ -217,7 +217,7 @@ OpenMMFrEnergyST::OpenMMFrEnergyST(const OpenMMFrEnergyST &other) isSystemInitialised(other.isSystemInitialised), isContextInitialised(other.isContextInitialised), combiningRules(other.combiningRules), CutoffType(other.CutoffType), cutoff_distance(other.cutoff_distance), field_dielectric(other.field_dielectric), Andersen_flag(other.Andersen_flag), - Andersen_frequency(other.Andersen_frequency), MCBarostat_flag(other.MCBarostat_flag), + Andersen_frequency(other.Andersen_frequency), MCBarostat_flag(other.MCBarostat_flag), MCBarostat_membrane_flag(other.MCBarostat_membrane_flag), MCBarostat_frequency(other.MCBarostat_frequency), ConstraintType(other.ConstraintType), Pressure(other.Pressure), Temperature(other.Temperature), platform_type(other.platform_type), Restraint_flag(other.Restraint_flag), CMMremoval_frequency(other.CMMremoval_frequency), buffer_frequency(other.buffer_frequency), @@ -259,6 +259,7 @@ OpenMMFrEnergyST &OpenMMFrEnergyST::operator=(const OpenMMFrEnergyST &other) Andersen_flag = other.Andersen_flag; Andersen_frequency = other.Andersen_frequency; MCBarostat_flag = other.MCBarostat_flag; + MCBarostat_membrane_flag = other.MCBarostat_membrane_flag; MCBarostat_frequency = other.MCBarostat_frequency; ConstraintType = other.ConstraintType; Pressure = other.Pressure; @@ -302,7 +303,7 @@ bool OpenMMFrEnergyST::operator==(const OpenMMFrEnergyST &other) const combiningRules == other.combiningRules and CutoffType == other.CutoffType and cutoff_distance == other.cutoff_distance and field_dielectric == other.field_dielectric and Andersen_flag == other.Andersen_flag and Andersen_frequency == other.Andersen_frequency and - MCBarostat_flag == other.MCBarostat_flag and MCBarostat_frequency == other.MCBarostat_frequency and + MCBarostat_flag == other.MCBarostat_flag and MCBarostat_membrane_flag == other.MCBarostat_membrane_flag and MCBarostat_frequency == other.MCBarostat_frequency and ConstraintType == other.ConstraintType and Pressure == other.Pressure and Temperature == other.Temperature and platform_type == other.platform_type and Restraint_flag == other.Restraint_flag and CMMremoval_frequency == other.CMMremoval_frequency and @@ -1241,13 +1242,29 @@ void OpenMMFrEnergyST::initialise() const double converted_Temperature = convertTo(Temperature.value(), kelvin); const double converted_Pressure = convertTo(Pressure.value(), bar); - OpenMM::MonteCarloBarostat *barostat = - new OpenMM::MonteCarloBarostat(converted_Pressure, converted_Temperature, MCBarostat_frequency); + if (MCBarostat_membrane_flag == true) + { + // Simple options for now: zero surface tension, XY isotropic, Z free + const double surface_Tension = 0; + OpenMM::MonteCarloMembraneBarostat::XYMode xymode = OpenMM::MonteCarloMembraneBarostat::XYIsotropic; + OpenMM::MonteCarloMembraneBarostat::ZMode zmode = OpenMM::MonteCarloMembraneBarostat::ZFree; + OpenMM::MonteCarloMembraneBarostat * barostat = new OpenMM::MonteCarloMembraneBarostat(converted_Pressure, surface_Tension, converted_Temperature, xymode, zmode, MCBarostat_frequency); - // Set The random seed - barostat->setRandomNumberSeed(random_seed); + //Set The random seed + barostat->setRandomNumberSeed(random_seed); + + system_openmm->addForce(barostat); + } + else + { + OpenMM::MonteCarloBarostat *barostat = + new OpenMM::MonteCarloBarostat(converted_Pressure, converted_Temperature, MCBarostat_frequency); - system_openmm->addForce(barostat); + // Set The random seed + barostat->setRandomNumberSeed(random_seed); + + system_openmm->addForce(barostat); + } if (Debug) { @@ -1255,6 +1272,10 @@ void OpenMMFrEnergyST::initialise() qDebug() << "Temperature = " << converted_Temperature << " K\n"; qDebug() << "Pressure = " << converted_Pressure << " bar\n"; qDebug() << "Frequency every " << MCBarostat_frequency << " steps\n"; + if (MCBarostat_membrane_flag) + { + qDebug() << "Membrane barostat, surface tension 0, XY isotropic, Z free\n"; + } } } /*******************************************************BONDED @@ -4351,6 +4372,16 @@ bool OpenMMFrEnergyST::getMCBarostat(void) return MCBarostat_flag; } +void OpenMMFrEnergyST::setMCBarostatMembrane(bool MCBarostat_membrane) +{ + MCBarostat_membrane_flag = MCBarostat_membrane; +} + +bool OpenMMFrEnergyST::getMCBarostatMembrane(void) +{ + return MCBarostat_membrane_flag; +} + /** Get the Monte Carlo Barostat frequency in time speps */ int OpenMMFrEnergyST::getMCBarostatFrequency(void) { diff --git a/corelib/src/libs/SireMove/openmmfrenergyst.h b/corelib/src/libs/SireMove/openmmfrenergyst.h index fded8e087..15fd80027 100644 --- a/corelib/src/libs/SireMove/openmmfrenergyst.h +++ b/corelib/src/libs/SireMove/openmmfrenergyst.h @@ -122,6 +122,9 @@ namespace SireMove bool getMCBarostat(void); void setMCBarostat(bool); + bool getMCBarostatMembrane(void); + void setMCBarostatMembrane(bool); + void setMCBarostatFrequency(int); int getMCBarostatFrequency(void); @@ -247,6 +250,7 @@ namespace SireMove double Andersen_frequency; bool MCBarostat_flag; + bool MCBarostat_membrane_flag; int MCBarostat_frequency; QString ConstraintType; diff --git a/corelib/src/libs/SireMove/openmmmdintegrator.cpp b/corelib/src/libs/SireMove/openmmmdintegrator.cpp index ee6200919..094952fa9 100644 --- a/corelib/src/libs/SireMove/openmmmdintegrator.cpp +++ b/corelib/src/libs/SireMove/openmmmdintegrator.cpp @@ -124,7 +124,7 @@ QDataStream &operator<<(QDataStream &ds, const OpenMMMDIntegrator &velver) sds << velver.frequent_save_velocities << velver.molgroup << velver.Integrator_type << velver.friction << velver.CutoffType << velver.cutoff_distance << velver.field_dielectric << velver.tolerance_ewald_pme - << velver.Andersen_flag << velver.Andersen_frequency << velver.MCBarostat_flag << velver.MCBarostat_frequency + << velver.Andersen_flag << velver.Andersen_frequency << velver.MCBarostat_flag << velver.MCBarostat_membrane_flag << velver.MCBarostat_frequency << velver.ConstraintType << velver.Pressure << velver.Temperature << velver.platform_type << velver.Restraint_flag << velver.CMMremoval_frequency << velver.buffer_frequency << velver.device_index << velver.LJ_dispersion << velver.precision << velver.integration_tol << velver.timeskip @@ -147,7 +147,7 @@ QDataStream &operator>>(QDataStream &ds, OpenMMMDIntegrator &velver) sds >> velver.frequent_save_velocities >> velver.molgroup >> velver.Integrator_type >> velver.friction >> velver.CutoffType >> velver.cutoff_distance >> velver.field_dielectric >> velver.tolerance_ewald_pme >> - velver.Andersen_flag >> velver.Andersen_frequency >> velver.MCBarostat_flag >> + velver.Andersen_flag >> velver.Andersen_frequency >> velver.MCBarostat_flag >> velver.MCBarostat_membrane_flag >> velver.MCBarostat_frequency >> velver.ConstraintType >> velver.Pressure >> velver.Temperature >> velver.platform_type >> velver.Restraint_flag >> velver.CMMremoval_frequency >> velver.buffer_frequency >> velver.device_index >> velver.LJ_dispersion >> velver.precision >> velver.integration_tol >> @@ -195,7 +195,7 @@ OpenMMMDIntegrator::OpenMMMDIntegrator(bool frequent_save) molgroup(MoleculeGroup()), openmm_system(0), openmm_context(0), isSystemInitialised(false), isContextInitialised(false), Integrator_type("leapfrogverlet"), friction(1.0 / picosecond), CutoffType("nocutoff"), cutoff_distance(1.0 * nanometer), field_dielectric(78.3), tolerance_ewald_pme(0.0001), - Andersen_flag(false), Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_frequency(25), + Andersen_flag(false), Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_membrane_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), Temperature(300.0 * kelvin), platform_type("Reference"), Restraint_flag(false), CMMremoval_frequency(0), buffer_frequency(0), device_index("0"), LJ_dispersion(true), precision("single"), reinetialise_context(false), integration_tol(0.001), timeskip(0.0 * picosecond), @@ -211,7 +211,7 @@ OpenMMMDIntegrator::OpenMMMDIntegrator(const MoleculeGroup &molecule_group, bool molgroup(molecule_group), openmm_system(0), openmm_context(0), isSystemInitialised(false), isContextInitialised(false), Integrator_type("leapfrogverlet"), friction(1.0 / picosecond), CutoffType("nocutoff"), cutoff_distance(1.0 * nanometer), field_dielectric(78.3), tolerance_ewald_pme(0.0001), - Andersen_flag(false), Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_frequency(25), + Andersen_flag(false), Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_membrane_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), Temperature(300.0 * kelvin), platform_type("Reference"), Restraint_flag(false), CMMremoval_frequency(0), buffer_frequency(0), device_index("0"), LJ_dispersion(true), precision("single"), reinetialise_context(false), integration_tol(0.001), timeskip(0.0 * picosecond), @@ -229,7 +229,7 @@ OpenMMMDIntegrator::OpenMMMDIntegrator(const OpenMMMDIntegrator &other) Integrator_type(other.Integrator_type), friction(other.friction), CutoffType(other.CutoffType), cutoff_distance(other.cutoff_distance), field_dielectric(other.field_dielectric), tolerance_ewald_pme(other.tolerance_ewald_pme), Andersen_flag(other.Andersen_flag), - Andersen_frequency(other.Andersen_frequency), MCBarostat_flag(other.MCBarostat_flag), + Andersen_frequency(other.Andersen_frequency), MCBarostat_flag(other.MCBarostat_flag), MCBarostat_membrane_flag(other.MCBarostat_membrane_flag), MCBarostat_frequency(other.MCBarostat_frequency), ConstraintType(other.ConstraintType), Pressure(other.Pressure), Temperature(other.Temperature), platform_type(other.platform_type), Restraint_flag(other.Restraint_flag), CMMremoval_frequency(other.CMMremoval_frequency), buffer_frequency(other.buffer_frequency), @@ -264,6 +264,7 @@ OpenMMMDIntegrator &OpenMMMDIntegrator::operator=(const OpenMMMDIntegrator &othe Andersen_flag = other.Andersen_flag; Andersen_frequency = other.Andersen_frequency; MCBarostat_flag = other.MCBarostat_flag; + MCBarostat_membrane_flag = other.MCBarostat_membrane_flag; MCBarostat_frequency = other.MCBarostat_frequency; ConstraintType = other.ConstraintType; Pressure = other.Pressure; @@ -291,7 +292,7 @@ bool OpenMMMDIntegrator::operator==(const OpenMMMDIntegrator &other) const isSystemInitialised == other.isSystemInitialised and isContextInitialised == other.isContextInitialised and CutoffType == other.CutoffType and cutoff_distance == other.cutoff_distance and field_dielectric == other.field_dielectric and Andersen_flag == other.Andersen_flag and - Andersen_frequency == other.Andersen_frequency and MCBarostat_flag == other.MCBarostat_flag and + Andersen_frequency == other.Andersen_frequency and MCBarostat_flag == other.MCBarostat_flag and MCBarostat_membrane_flag == other.MCBarostat_membrane_flag and MCBarostat_frequency == other.MCBarostat_frequency and ConstraintType == other.ConstraintType and Pressure == other.Pressure and Temperature == other.Temperature and platform_type == other.platform_type and Restraint_flag == other.Restraint_flag and CMMremoval_frequency == other.CMMremoval_frequency and @@ -487,9 +488,21 @@ void OpenMMMDIntegrator::initialise() const double converted_Temperature = convertTo(Temperature.value(), kelvin); const double converted_Pressure = convertTo(Pressure.value(), bar); - OpenMM::MonteCarloBarostat *barostat = - new OpenMM::MonteCarloBarostat(converted_Pressure, converted_Temperature, MCBarostat_frequency); - system_openmm->addForce(barostat); + if (MCBarostat_membrane_flag == true) + { + // Simple options for now: zero surface tension, XY isotropic, Z free + const double surface_Tension = 0; + OpenMM::MonteCarloMembraneBarostat::XYMode xymode = OpenMM::MonteCarloMembraneBarostat::XYIsotropic; + OpenMM::MonteCarloMembraneBarostat::ZMode zmode = OpenMM::MonteCarloMembraneBarostat::ZFree; + OpenMM::MonteCarloMembraneBarostat * barostat = new OpenMM::MonteCarloMembraneBarostat(converted_Pressure, surface_Tension, converted_Temperature, xymode, zmode, MCBarostat_frequency); + system_openmm->addForce(barostat); + } + else // normal barostat + { + OpenMM::MonteCarloBarostat *barostat = + new OpenMM::MonteCarloBarostat(converted_Pressure, converted_Temperature, MCBarostat_frequency); + system_openmm->addForce(barostat); + } if (Debug) { @@ -498,6 +511,10 @@ void OpenMMMDIntegrator::initialise() qDebug() << "Pressure = " << converted_Pressure << " bar\n"; qDebug() << "Frequency every " << MCBarostat_frequency << " steps\n"; qDebug() << "Lennard Jones Dispersion term is set to " << LJ_dispersion << "\n"; + if (MCBarostat_membrane_flag) + { + qDebug() << "Membrane barostat, surface tension 0, XY isotropic, Z free\n"; + } } } @@ -1928,6 +1945,16 @@ bool OpenMMMDIntegrator::getMCBarostat(void) return MCBarostat_flag; } +void OpenMMMDIntegrator::setMCBarostatMembrane(bool MCBarostatMembrane) +{ + MCBarostat_membrane_flag = MCBarostatMembrane; +} + +bool OpenMMMDIntegrator::getMCBarostatMembrane(void) +{ + return MCBarostat_membrane_flag; +} + /** Get the Monte Carlo Barostat frequency in time speps */ int OpenMMMDIntegrator::getMCBarostatFrequency(void) { diff --git a/corelib/src/libs/SireMove/openmmmdintegrator.h b/corelib/src/libs/SireMove/openmmmdintegrator.h index 0e026ca9f..48953cb7b 100644 --- a/corelib/src/libs/SireMove/openmmmdintegrator.h +++ b/corelib/src/libs/SireMove/openmmmdintegrator.h @@ -127,6 +127,9 @@ namespace SireMove bool getMCBarostat(void); void setMCBarostat(bool); + bool getMCBarostatMembrane(void); + void setMCBarostatMembrane(bool); + void setMCBarostatFrequency(int); int getMCBarostatFrequency(void); @@ -203,6 +206,7 @@ namespace SireMove double Andersen_frequency; bool MCBarostat_flag; + bool MCBarostat_membrane_flag; int MCBarostat_frequency; QString ConstraintType; diff --git a/corelib/src/libs/SireMove/openmmpmefep.cpp b/corelib/src/libs/SireMove/openmmpmefep.cpp index f85086460..9b00f0b0d 100644 --- a/corelib/src/libs/SireMove/openmmpmefep.cpp +++ b/corelib/src/libs/SireMove/openmmpmefep.cpp @@ -144,7 +144,7 @@ QDataStream &operator<<(QDataStream &ds, const OpenMMPMEFEP &velver) sds << velver.frequent_save_velocities << velver.molgroup << velver.solute << velver.solutehard << velver.solutetodummy << velver.solutefromdummy << velver.combiningRules << velver.CutoffType << velver.cutoff_distance << velver.field_dielectric << velver.Andersen_flag << velver.Andersen_frequency - << velver.MCBarostat_flag << velver.MCBarostat_frequency << velver.ConstraintType << velver.Pressure + << velver.MCBarostat_flag << velver.MCBarostat_membrane_flag << velver.MCBarostat_frequency << velver.ConstraintType << velver.Pressure << velver.Temperature << velver.platform_type << velver.Restraint_flag << velver.CMMremoval_frequency << velver.buffer_frequency << velver.energy_frequency << velver.device_index << velver.precision << velver.current_lambda << velver.coulomb_power << velver.shift_delta << velver.delta_alchemical @@ -166,7 +166,7 @@ QDataStream &operator>>(QDataStream &ds, OpenMMPMEFEP &velver) sds >> velver.frequent_save_velocities >> velver.molgroup >> velver.solute >> velver.solutehard >> velver.solutetodummy >> velver.solutefromdummy >> velver.combiningRules >> velver.CutoffType >> velver.cutoff_distance >> velver.field_dielectric >> velver.Andersen_flag >> velver.Andersen_frequency >> - velver.MCBarostat_flag >> velver.MCBarostat_frequency >> velver.ConstraintType >> velver.Pressure >> + velver.MCBarostat_flag >> velver.MCBarostat_membrane_flag >> velver.MCBarostat_frequency >> velver.ConstraintType >> velver.Pressure >> velver.Temperature >> velver.platform_type >> velver.Restraint_flag >> velver.CMMremoval_frequency >> velver.buffer_frequency >> velver.energy_frequency >> velver.device_index >> velver.precision >> velver.current_lambda >> velver.coulomb_power >> velver.shift_delta >> velver.delta_alchemical >> @@ -194,7 +194,7 @@ OpenMMPMEFEP::OpenMMPMEFEP(bool frequent_save) solutefromdummy(MoleculeGroup()), openmm_system(0), openmm_context(0), isSystemInitialised(false), isContextInitialised(false), combiningRules("arithmetic"), CutoffType("PME"), cutoff_distance(1.0 * nanometer), field_dielectric(78.3), Andersen_flag(false), Andersen_frequency(90.0), MCBarostat_flag(false), - MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), Temperature(300.0 * kelvin), + MCBarostat_membrane_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), Temperature(300.0 * kelvin), platform_type("Reference"), Restraint_flag(false), CMMremoval_frequency(0), buffer_frequency(0), energy_frequency(100), device_index("0"), precision("single"), current_lambda(0.5), coulomb_power(0), shift_delta(2.0), delta_alchemical(0.001), alchemical_array(), finite_diff_gradients(), pot_energies(), @@ -212,7 +212,7 @@ OpenMMPMEFEP::OpenMMPMEFEP(const MoleculeGroup &molecule_group, const MoleculeGr solute(solute_group), solutehard(solute_hard), solutetodummy(solute_todummy), solutefromdummy(solute_fromdummy), openmm_system(0), openmm_context(0), isSystemInitialised(false), isContextInitialised(false), combiningRules("arithmetic"), CutoffType("PME"), cutoff_distance(1.0 * nanometer), field_dielectric(78.3), - Andersen_flag(false), Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_frequency(25), + Andersen_flag(false), Andersen_frequency(90.0), MCBarostat_flag(false), MCBarostat_membrane_flag(false), MCBarostat_frequency(25), ConstraintType("none"), Pressure(1.0 * bar), Temperature(300.0 * kelvin), platform_type("Reference"), Restraint_flag(false), CMMremoval_frequency(0), buffer_frequency(0), energy_frequency(100), device_index("0"), precision("single"), current_lambda(0.5), coulomb_power(0), shift_delta(2.0), delta_alchemical(0.001), @@ -230,7 +230,7 @@ OpenMMPMEFEP::OpenMMPMEFEP(const OpenMMPMEFEP &other) isSystemInitialised(other.isSystemInitialised), isContextInitialised(other.isContextInitialised), combiningRules(other.combiningRules), CutoffType(other.CutoffType), cutoff_distance(other.cutoff_distance), field_dielectric(other.field_dielectric), Andersen_flag(other.Andersen_flag), - Andersen_frequency(other.Andersen_frequency), MCBarostat_flag(other.MCBarostat_flag), + Andersen_frequency(other.Andersen_frequency), MCBarostat_flag(other.MCBarostat_flag), MCBarostat_membrane_flag(other.MCBarostat_membrane_flag), MCBarostat_frequency(other.MCBarostat_frequency), ConstraintType(other.ConstraintType), Pressure(other.Pressure), Temperature(other.Temperature), platform_type(other.platform_type), Restraint_flag(other.Restraint_flag), CMMremoval_frequency(other.CMMremoval_frequency), buffer_frequency(other.buffer_frequency), @@ -272,6 +272,7 @@ OpenMMPMEFEP &OpenMMPMEFEP::operator=(const OpenMMPMEFEP &other) Andersen_flag = other.Andersen_flag; Andersen_frequency = other.Andersen_frequency; MCBarostat_flag = other.MCBarostat_flag; + MCBarostat_membrane_flag = other.MCBarostat_membrane_flag; MCBarostat_frequency = other.MCBarostat_frequency; ConstraintType = other.ConstraintType; Pressure = other.Pressure; @@ -315,7 +316,7 @@ bool OpenMMPMEFEP::operator==(const OpenMMPMEFEP &other) const combiningRules == other.combiningRules and CutoffType == other.CutoffType and cutoff_distance == other.cutoff_distance and field_dielectric == other.field_dielectric and Andersen_flag == other.Andersen_flag and Andersen_frequency == other.Andersen_frequency and - MCBarostat_flag == other.MCBarostat_flag and MCBarostat_frequency == other.MCBarostat_frequency and + MCBarostat_flag == other.MCBarostat_flag and MCBarostat_membrane_flag == other.MCBarostat_membrane_flag and MCBarostat_frequency == other.MCBarostat_frequency and ConstraintType == other.ConstraintType and Pressure == other.Pressure and Temperature == other.Temperature and platform_type == other.platform_type and Restraint_flag == other.Restraint_flag and CMMremoval_frequency == other.CMMremoval_frequency and @@ -386,17 +387,39 @@ void OpenMMPMEFEP::addMCBarostat(OpenMM::System &system) const double converted_Temperature = convertTo(Temperature.value(), kelvin); const double converted_Pressure = convertTo(Pressure.value(), bar); - auto barostat = new OpenMM::MonteCarloBarostat(converted_Pressure, converted_Temperature, MCBarostat_frequency); + if (MCBarostat_membrane_flag == true) + { + // Simple options for now: zero surface tension, XY isotropic, Z free + const double surface_Tension = 0; + OpenMM::MonteCarloMembraneBarostat::XYMode xymode = OpenMM::MonteCarloMembraneBarostat::XYIsotropic; + OpenMM::MonteCarloMembraneBarostat::ZMode zmode = OpenMM::MonteCarloMembraneBarostat::ZFree; + auto barostat = new OpenMM::MonteCarloMembraneBarostat(converted_Pressure, surface_Tension, converted_Temperature, xymode, zmode, MCBarostat_frequency); + + //Set The random seed + barostat->setRandomNumberSeed(random_seed); + + system.addForce(barostat); + } + else + { + auto barostat = new OpenMM::MonteCarloBarostat(converted_Pressure, converted_Temperature, MCBarostat_frequency); - barostat->setRandomNumberSeed(random_seed); - system.addForce(barostat); + // Set The random seed + barostat->setRandomNumberSeed(random_seed); + + system.addForce(barostat); + } if (Debug) { - qDebug() << "Monte Carlo Barostat set"; - qDebug() << "Temperature =" << converted_Temperature << " K"; - qDebug() << "Pressure =" << converted_Pressure << " bar"; - qDebug() << "Frequency every" << MCBarostat_frequency << " steps"; + qDebug() << "Monte Carlo Barostat set\n"; + qDebug() << "Temperature = " << converted_Temperature << " K\n"; + qDebug() << "Pressure = " << converted_Pressure << " bar\n"; + qDebug() << "Frequency every " << MCBarostat_frequency << " steps\n"; + if (MCBarostat_membrane_flag) + { + qDebug() << "Membrane barostat, surface tension 0, XY isotropic, Z free\n"; + } } } @@ -3665,8 +3688,8 @@ boost::tuples::tuple OpenMMPMEFEP::calculateGradient(dou { double double_increment = incr_plus - incr_minus; double gradient = 0; - double potential_energy_lambda_plus_delta; - double potential_energy_lambda_minus_delta; + double potential_energy_lambda_plus_delta = 0; + double potential_energy_lambda_minus_delta = 0; double forward_m; double backward_m; if (incr_plus < 1.0) @@ -3887,6 +3910,16 @@ bool OpenMMPMEFEP::getMCBarostat(void) return MCBarostat_flag; } +void OpenMMPMEFEP::setMCBarostatMembrane(bool MCBarostat_membrane) +{ + MCBarostat_membrane_flag = MCBarostat_membrane; +} + +bool OpenMMPMEFEP::getMCBarostatMembrane(void) +{ + return MCBarostat_membrane_flag; +} + /** Get the Monte Carlo Barostat frequency in time speps */ int OpenMMPMEFEP::getMCBarostatFrequency(void) { diff --git a/corelib/src/libs/SireMove/openmmpmefep.h b/corelib/src/libs/SireMove/openmmpmefep.h index a6c57ae3c..c73b7267e 100644 --- a/corelib/src/libs/SireMove/openmmpmefep.h +++ b/corelib/src/libs/SireMove/openmmpmefep.h @@ -127,6 +127,9 @@ namespace SireMove bool getMCBarostat(void); void setMCBarostat(bool); + bool getMCBarostatMembrane(void); + void setMCBarostatMembrane(bool); + void setMCBarostatFrequency(int); int getMCBarostatFrequency(void); @@ -250,6 +253,7 @@ namespace SireMove double Andersen_frequency; bool MCBarostat_flag; + bool MCBarostat_membrane_flag; int MCBarostat_frequency; QString ConstraintType; diff --git a/corelib/src/libs/SireSystem/system.cpp b/corelib/src/libs/SireSystem/system.cpp index 874d0abb1..868a2b883 100644 --- a/corelib/src/libs/SireSystem/system.cpp +++ b/corelib/src/libs/SireSystem/system.cpp @@ -4090,6 +4090,7 @@ void System::deleteAllFrames(const SireBase::PropertyMap &map) { this->accept(); this->mustNowRecalculateFromScratch(); + this->system_trajectory = 0; MolGroupsBase::deleteAllFrames(map); } diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 2cda5c669..1ee8d41cc 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -12,15 +12,38 @@ Development was migrated into the `OpenBioSim `__ organisation on `GitHub `__. -`2024.4.0 `__ - December 2024 --------------------------------------------------------------------------------------------- - -* Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR. +`2024.4.0 `__ - Feb 2025 +---------------------------------------------------------------------------------------- * Fixed update of triclinic box vectors in ``SOMD`` following ``OpenMM`` bug fix. +* Don't automatically save energies and frames when ``dynamics.run()`` returns. + +* Improved handling of ``OpenMM`` NaN errors during dynamics. + +* Restore thread state before raising exceptions in the Sire OpenMM minimiser. + +* Add support for Replica Exchange with Solute Tempering (REST2) simulations. + +* Null SystemTrajectory pointer when all frames are deleted. + +* Handle cis/trans double bond stereochemistry values in SDF bond blocks. + +* Fix handling of positive formal charge when writing SDF files. + +* Fix lambda schedule discussion and plots in the tutorial. + +* Added support for angle and dihedral restraints which can be used in alchemical and standard simulations. + +* Allow user to compute energy trajectory over a subset of the lambda windows for each lambda. + +* Fix the ``hasForceSpecificEquation`` function in the ``LambdaLever`` class so that it returns true if + there is a default equation for the force. + +* Added support for the ``OpenMM`` ``MonteCarloMembraneBarostat`` to ``SOMD``. + `2024.3.1 `__ - December 2024 ---------------------------------------------------------------------------------------------- +-------------------------------------------------------------------------------------------- * Fixed instantiaton of ``QByteArray`` in ``Sire::Mol::Frame::toByteArray`` and count bytes with ``QByteArray::size``. @@ -32,17 +55,17 @@ organisation on `GitHub `__. * Make minimisation function settings consistent across API. -* Reload ``TorchQMForce`` module if OpenMM platform changes. +* Reload TorchQMForce module if OpenMM platform changes. -* Fixed calculation of non-equilibrium work when starting from QM state. +* Fix calculation of non-equilibrium work when starting from QM state. -* Fixed stereochemistry in RDKit molecule conversion. +* Fix stereochemistry in RDKit molecule conversion. * Fixed :func:`sire.restraints.get_standard_state_correction` to be consistent with the definition of the "force constanst" as ``F = 2 ** k ** x`` (rather than ``F = k ** x``). Updated docstrings and restraints documentation to make it clear how the force constants are defined. -* Fixed thread safety issue in Sire OpenMM minimiser. +* Fix thread safety issue in Sire OpenMM minimiser. `2024.3.0 `__ - October 2024 -------------------------------------------------------------------------------------------- diff --git a/doc/source/contributors.rst b/doc/source/contributors.rst index 495575326..4622d0ba3 100644 --- a/doc/source/contributors.rst +++ b/doc/source/contributors.rst @@ -22,3 +22,4 @@ can be recognised. * `@msuruzon `__ * `@kexul `__ * `@mb2055 `__ +* `@akalpokas `__ diff --git a/doc/source/tutorial/part06/02_alchemical_dynamics.rst b/doc/source/tutorial/part06/02_alchemical_dynamics.rst index 52c4861ab..15765f992 100644 --- a/doc/source/tutorial/part06/02_alchemical_dynamics.rst +++ b/doc/source/tutorial/part06/02_alchemical_dynamics.rst @@ -23,6 +23,20 @@ the simulation that will control how it is run: λ is a parameter that controls the morph from the reference state (at λ=0) to the perturbed state (at λ=1). +* ``rest2_scale`` - this sets the scale factor for Replica Exchange + with Solute Tempering (REST2) simulations. This is a floating point + number that defaults to ``1.0``. The scale factor defines the temperature + of the REST2 region relative to the rest of the system, which can be used + to soften dihedral potentials and non-bonded interactions in the REST2 + region to aid conformational sampling. + +* ``rest2_selection`` - this is a selection string that specifies additional + atoms to include in the REST2 region, e.g. relevant atoms within the protein + binding site. (By default, the REST2 region comprises all atoms in perturbable + molecules.) If the selection includes atoms from a perturbable molecule, then + only those atoms within the perturbable molecule will be used, i.e. this + allows a user to specify a sub-set of atoms within a perturbable molecule. + * ``swap_end_states`` - if set to ``True``, this will swap the end states of the perturbation. The morph will run from the perturbed state (at λ=0) to the reference state (at λ=1). Note that the coordinates @@ -214,7 +228,7 @@ 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() +>>> s = d.get_schedule() >>> print(s) LambdaSchedule( morph: λ * final + initial * (-λ + 1) @@ -224,13 +238,13 @@ 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() +>>> df.plot(subplots=True) .. 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 +a hypothetical 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 @@ -249,7 +263,7 @@ 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() +>>> s.get_lever_values(initial=2.0, final=3.0).plot(subplots=True) .. image:: images/06_02_02.jpg :alt: View of the λ-schedule where charge is scaled by λ^2 @@ -276,10 +290,10 @@ LambdaSchedule( charge: final * (-λ + 1) ) -Again, it is worth plotting the impact of this schedule on a hyperthetical +Again, it is worth plotting the impact of this schedule on a hypothetical parameter. ->>> s.get_lever_values(initial=2.0, final=3.0).plot() +>>> s.get_lever_values(initial=2.0, final=3.0).plot(subplots=True) .. image:: images/06_02_03.jpg :alt: View of the λ-schedule where charge is scaled in a second stage to zero. @@ -331,23 +345,15 @@ 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 +We can see how this would affect a hypothetical 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() +>>> s.get_lever_values(initial=2.0, final=3.0).plot(subplots=True) .. 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:: - - 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. @@ -452,7 +458,8 @@ are; new value of λ for the context. Note that this should only really be used to change λ to evaluate energies at different λ-windows. It is better to re-create the context if you want to simulate - at a different λ-value. + at a different λ-value. This function can also be used to set the + ``rest2_scale`` parameter for the context. * :func:`~sire.Convert.SireOpenMM.SOMMContext.get_lambda_schedule` - return the λ-schedule used to control the morph. * :func:`~sire.Convert.SireOpenMM.SOMMContext.set_lambda_schedule` - set the diff --git a/doc/source/tutorial/part06/03_restraints.rst b/doc/source/tutorial/part06/03_restraints.rst index 341dca762..e510098f1 100644 --- a/doc/source/tutorial/part06/03_restraints.rst +++ b/doc/source/tutorial/part06/03_restraints.rst @@ -236,6 +236,74 @@ BondRestraints( size=630 the first atom in the first molecule and the oxygen atoms in all of the water molecules. +Angle or Dihedral Restraints +--------------------------- + +The :func:`sire.restraints.angle` or :func:`sire.restraints.dihedral` functions +are used to create angle or distance restraints. + +Just like the other restraint functions, these functions take +the set of molecules or system you want to simulate, +plus a search string, lists of atom indexes, or molecule views +holding the atoms that you want to restrain., e.g. + +>>> restraints = sr.restraints.angle(mols=mols, atoms=[0, 1, 2]) +>>> print(restraints) +AngleRestraints( name=restraint, size=1 +0: AngleRestraint( [0, 1, 2], theta0=112.8°, ktheta=0.0304617 kcal mol-1 °-2 ) + ) + +or + +>>> restraints = sr.restraints.dihedral(mols=mols, atoms=[0, 1, 2, 3]) +>>> print(restraints) +DihedralRestraints( name=restraint, size=1 +0: DihedralRestraint( [0, 1, 2, 3], phi0=244.528°, kphi=0.0304617 kcal mol-1 °-2 ) + ) + +creates a single harmonic angle or dihedral restraint that acts between +the specified atoms. By default, the equilibrium angle (theta0 or phi0) +is the current angle between the atoms (112.8° or 244.528°), +and the force constant (ktheta or kphi) is 100 kcal mol-1 rad-2. + +You can set these via the ``ktheta`` or ``kphi`` and ``theta0`` or ``phi0`` arguments depending +on the restraint used, e.g. + +>>> restraints = sr.restraints.angle(mols=mols, atoms=[0, 1, 2], +... theta0="1.5 rad", ktheta="50 kcal mol-1 rad-2") +>>> print(restraints) +AngleRestraints( name=restraint, size=1 +0: AngleRestraint( [0, 1, 2], theta0=85.9437°, ktheta=0.0152309 kcal mol-1 °-2 ) + ) + +or + +>>> restraints = sr.restraints.dihedral(mols=mols, atoms=[0, 1, 2, 3], +... phi0="2 rad", kphi="10 kcal mol-1 rad-2") +>>> print(restraints) +DihedralRestraints( name=restraint, size=1 +0: DihedralRestraint( [0, 1, 2, 3], phi0=114.592°, kphi=0.00304617 kcal mol-1 °-2 ) + ) + +You can specify the atoms using a search string, passing the atoms themselves, +or passing the atoms from a molecular container. + +>>> ang = mols.angles()[0] +>>> restraints = sr.restraints.angle(mols=mols, atoms=ang.atoms()) +>>> print(restraints) +AngleRestraints( name=restraint, size=1 +0: AngleRestraint( [0, 1, 2], theta0=112.8°, ktheta=0.0304617 kcal mol-1 °-2 ) + ) + +or + +>>> dih = mols.dihedrals()[0] +>>> restraints = sr.restraints.dihedral(mols=mols, atoms=dih.atoms()) +>>> print(restraints) +DihedralRestraints( name=restraint, size=1 +0: DihedralRestraint( [0, 1, 4, 5], phi0=243.281°, kphi=0.0304617 kcal mol-1 °-2 ) + ) + Boresch Restraints --------------------------- diff --git a/doc/source/tutorial/part06/images/06_02_01.jpg b/doc/source/tutorial/part06/images/06_02_01.jpg index 8674ec7aa..461cb16eb 100644 Binary files a/doc/source/tutorial/part06/images/06_02_01.jpg 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 index 5618f084b..1198cde9e 100644 Binary files a/doc/source/tutorial/part06/images/06_02_02.jpg 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 index 44694b4c6..fe39f9ab0 100644 Binary files a/doc/source/tutorial/part06/images/06_02_03.jpg 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 index 56af4e349..e7b0348c3 100644 Binary files a/doc/source/tutorial/part06/images/06_02_04.jpg and b/doc/source/tutorial/part06/images/06_02_04.jpg differ diff --git a/src/sire/mm/__init__.py b/src/sire/mm/__init__.py index 2255fd1ac..9df64e5d1 100644 --- a/src/sire/mm/__init__.py +++ b/src/sire/mm/__init__.py @@ -4,6 +4,8 @@ "AmberDihPart", "AmberDihedral", "Angle", + "AngleRestraint", + "AngleRestraints", "Bond", "BondRestraint", "BondRestraints", @@ -11,6 +13,8 @@ "Improper", "PositionRestraint", "PositionRestraints", + "DihedralRestraint", + "DihedralRestraints", "SelectorAngle", "SelectorBond", "SelectorDihedral", @@ -29,6 +33,9 @@ _use_new_api() +AngleRestraint = _MM.AngleRestraint +AngleRestraints = _MM.AngleRestraints + # It would be better if these were called "DistanceRestraints", # but there is already a legacy Sire.MM class with this name BondRestraint = _MM.BondRestraint @@ -40,6 +47,9 @@ PositionalRestraint = _MM.PositionalRestraint PositionalRestraints = _MM.PositionalRestraints +DihedralRestraint = _MM.DihedralRestraint +DihedralRestraints = _MM.DihedralRestraints + AmberBond = _MM.AmberBond AmberAngle = _MM.AmberAngle AmberDihPart = _MM.AmberDihPart diff --git a/src/sire/mol/__init__.py b/src/sire/mol/__init__.py index ad8004fa9..f3a7ded7d 100644 --- a/src/sire/mol/__init__.py +++ b/src/sire/mol/__init__.py @@ -1567,6 +1567,8 @@ def _dynamics( ignore_perturbations=None, temperature=None, pressure=None, + rest2_scale=None, + rest2_selection=None, vacuum=None, shift_delta=None, shift_coulomb=None, @@ -1686,6 +1688,30 @@ def _dynamics( microcanonical (NVE) or canonical (NVT) simulation will be run if the pressure is not set. + rest2_scale: float + The scaling factor to apply when running a REST2 simulation + (Replica Exchange with Solute Tempering). This defaults to 1.0. + This specifies the ratio of the temperature of the REST2 region + to the temperature of the rest of the system. + + rest2_selection: str + A selection string for atoms to include in the REST2 region + in addition to any perturbable molecules. For example, + "molidx 0 and residx 0,1,2" would select atoms from the first + three residues of the first molecule. If None, then all atoms + within perturbable molecules will be included in the REST2 + region. When atoms within a perturbable molecule are also + included in the selection, then only those atoms will be + considered as part of the REST2 region. This allows REST2 to + be applied to protein mutations. + + A selection string for atoms to include in the REST2 region in + addition to the perturbable molecule. For example, + "molidx 0 and residx 0,1,2" would select atoms from the first + three residues of the first molecule. Selected atoms should be + within the same molecule. If None, then only the perturbable + molecule will be used for the REST2 region. + vacuum: bool (optional) Whether or not to run the simulation in vacuum. If this is set to `True`, then the simulation space automatically be @@ -1852,6 +1878,12 @@ def _dynamics( pressure = u(pressure) map.set("pressure", pressure) + if rest2_scale is not None: + map.set("rest2_scale", rest2_scale) + + if rest2_selection is not None: + map.set("rest2_selection", rest2_selection) + if device is not None: map.set("device", str(device)) diff --git a/src/sire/mol/_dynamics.py b/src/sire/mol/_dynamics.py index fe6eb7d2d..5403075e9 100644 --- a/src/sire/mol/_dynamics.py +++ b/src/sire/mol/_dynamics.py @@ -179,6 +179,84 @@ def __init__(self, mols=None, map=None, **kwargs): self._is_running = False self._schedule_changed = False + # Check for a REST2 scaling factor. + if map.specified("rest2_scale"): + try: + rest2_scale = map["rest2_scale"].value().as_double() + except: + raise ValueError("'rest2_scale' must be of type 'float'") + if rest2_scale < 1.0: + raise ValueError("'rest2_scale' must be >= 1.0") + else: + rest2_scale = 1.0 + + # Check for an additional REST2 selection. + if map.specified("rest2_selection"): + try: + rest2_selection = str(map["rest2_selection"]) + except: + raise ValueError("'rest2_selection' must be of type 'str'") + + try: + from . import selection_to_atoms + + # Try to find the REST2 selection. + atoms = selection_to_atoms(mols, rest2_selection) + except: + raise ValueError( + "Invalid 'rest2_selection' string. Please check the selection syntax." + ) + + # Store all the perturbable molecules associated with the selection + # and remove perturbable atoms from the selection. Remove alchemical ions + # from the selection. + pert_mols = {} + non_pert_atoms = atoms.to_list() + for atom in atoms: + mol = atom.molecule() + if mol.has_property("is_alchemical_ion"): + non_pert_atoms.remove(atom) + elif mol.has_property("is_perturbable"): + non_pert_atoms.remove(atom) + if mol.number() not in pert_mols: + pert_mols[mol.number()] = [atom] + else: + pert_mols[mol.number()].append(atom) + + # Now create a boolean is_rest2 mask for the atoms in the perturbable molecules. + # Only do this if there are perturbable atoms in the selection. + if len(non_pert_atoms) != len(atoms): + for num in pert_mols: + mol = self._sire_mols[num] + is_rest2 = [False] * mol.num_atoms() + for atom in pert_mols[num]: + is_rest2[atom.index().value()] = True + + # Set the is_rest2 property for each perturbable molecule. + mol = ( + mol.edit() + .set_property("is_rest2", is_rest2) + .molecule() + .commit() + ) + + # Update the system. + self._sire_mols.update(mol) + + # Search for alchemical ions and exclude them via a REST2 mask. + try: + for mol in self._sire_mols.molecules("property is_alchemical_ion"): + is_rest2 = [False] * mol.num_atoms() + mol = ( + mol.edit() + .set_property("is_rest2", is_rest2) + .molecule() + .commit() + ) + self._sire_mols.update(mol) + except: + pass + from ..convert import to self._omm_mols = to(self._sire_mols, "openmm", map=self._map) @@ -194,6 +272,10 @@ def __init__(self, mols=None, map=None, **kwargs): else: self._enforce_periodic_box = False + # Prepare the OpenMM REST2 data structures. + if map.specified("rest2_selection"): + if len(non_pert_atoms) > 0: + self._omm_mols._prepare_rest2(self._sire_mols, non_pert_atoms) else: self._sire_mols = None self._energy_trajectory = None @@ -209,6 +291,8 @@ def _update_from(self, state, state_has_cv, nsteps_completed): if self.is_null(): return + self._current_step = nsteps_completed + if not state_has_cv[0]: # there is no information to update return @@ -247,8 +331,6 @@ def _update_from(self, state, state_has_cv, nsteps_completed): map=self._map, ) - self._current_step = nsteps_completed - self._sire_mols.update(mols_to_update.molecules()) if self._ffinfo.space().is_periodic(): @@ -273,8 +355,11 @@ def _exit_dynamics_block( save_frame: bool = False, save_energy: bool = False, lambda_windows=[], + rest2_scale_factors=[], save_velocities: bool = False, delta_lambda: float = None, + num_energy_neighbours: int = None, + null_energy: float = None, ): if not self._is_running: raise SystemError("Cannot stop dynamics that is not running!") @@ -304,6 +389,14 @@ def _exit_dynamics_block( self._elapsed_time = current_time self._current_time += delta + # store the number of lambda windows + if lambda_windows is not None: + num_lambda_windows = len(lambda_windows) + + # compute energies for all windows + if num_energy_neighbours is None: + num_energy_neighbours = num_lambda_windows + if save_energy: # should save energy here nrgs = {} @@ -323,8 +416,9 @@ def _exit_dynamics_block( ) sim_lambda_value = self._omm_mols.get_lambda() + sim_rest2_scale = self._omm_mols.get_rest2_scale() - # Store the potential energy and accumulated non-equilibrium work. + # store the potential energy and accumulated non-equilibrium work if self._is_interpolate: nrg = nrgs["potential"] @@ -339,19 +433,38 @@ def _exit_dynamics_block( nrgs[str(sim_lambda_value)] = nrgs["potential"] if lambda_windows is not None: - for lambda_value in lambda_windows: + # get the index of the simulation lambda value in the + # lambda windows list + try: + lambda_index = lambda_windows.index(sim_lambda_value) + except: + num_energy_neighbours = num_lambda_windows + lambda_index = i + + for i, (lambda_value, rest2_scale) in enumerate( + zip(lambda_windows, rest2_scale_factors) + ): if lambda_value != sim_lambda_value: - self._omm_mols.set_lambda( - lambda_value, update_constraints=False - ) - nrgs[str(lambda_value)] = ( - self._omm_mols.get_potential_energy( - to_sire_units=False - ).value_in_unit(openmm.unit.kilocalorie_per_mole) - * kcal_per_mol - ) + if abs(lambda_index - i) <= num_energy_neighbours: + self._omm_mols.set_lambda( + lambda_value, + rest2_scale=rest2_scale, + update_constraints=False, + ) + nrgs[str(lambda_value)] = ( + self._omm_mols.get_potential_energy( + to_sire_units=False + ).value_in_unit(openmm.unit.kilocalorie_per_mole) + * kcal_per_mol + ) + else: + nrgs[str(lambda_value)] = null_energy * kcal_per_mol - self._omm_mols.set_lambda(sim_lambda_value, update_constraints=False) + self._omm_mols.set_lambda( + sim_lambda_value, + rest2_scale=sim_rest2_scale, + update_constraints=False, + ) if self._is_interpolate: self._energy_trajectory.set( @@ -545,7 +658,10 @@ def set_schedule(self, schedule): if not self.is_null(): self._omm_mols.set_lambda_schedule(schedule) self._schedule_changed = True - self.set_lambda(self._omm_mols.get_lambda()) + self.set_lambda( + self._omm_mols.get_lambda(), + rest2_scale=self._omm_mols.get_rest2_scale(), + ) def get_lambda(self): if self.is_null(): @@ -553,7 +669,12 @@ def get_lambda(self): else: return self._omm_mols.get_lambda() - def set_lambda(self, lambda_value: float, update_constraints: bool = True): + def set_lambda( + self, + lambda_value: float, + rest2_scale: float = 1.0, + update_constraints: bool = True, + ): if not self.is_null(): s = self.get_schedule() @@ -562,14 +683,18 @@ def set_lambda(self, lambda_value: float, update_constraints: bool = True): lambda_value = s.clamp(lambda_value) - if (not self._schedule_changed) and ( - lambda_value == self._omm_mols.get_lambda() + if ( + (not self._schedule_changed) + and (lambda_value == self._omm_mols.get_lambda()) + and (rest2_scale == self._omm_mols.get_rest2_scale()) ): # nothing to do return self._omm_mols.set_lambda( - lambda_value=lambda_value, update_constraints=update_constraints + lambda_value=lambda_value, + rest2_scale=rest2_scale, + update_constraints=update_constraints, ) self._schedule_changed = False self._clear_state() @@ -780,13 +905,11 @@ def _rebuild_and_minimise(self): from ..utils import Console Console.warning( - "Something went wrong when running dynamics. Since no steps " - "were completed, it is likely that the system needs minimising. " - "The system will be minimised, and then dynamics will be " - "attempted again. If an error still occurs, then it is likely " - "that the step size is too large, the molecules are " - "over-constrained, or there is something more fundemental " - "going wrong..." + "Something went wrong when running dynamics. The system will be " + "minimised, and then dynamics will be attempted again. If an " + "error still occurs, then it is likely that the step size is too " + "large, the molecules are over-constrained, or there is something " + "more fundemental going wrong..." ) # rebuild the molecules @@ -803,8 +926,13 @@ def run( frame_frequency=None, energy_frequency=None, lambda_windows=None, + rest2_scale_factors=None, save_velocities: bool = None, + save_frame_on_exit: bool = False, + save_energy_on_exit: bool = False, auto_fix_minimise: bool = True, + num_energy_neighbours: int = None, + null_energy: str = None, ): if self.is_null(): return @@ -815,8 +943,13 @@ def run( "frame_frequency": frame_frequency, "energy_frequency": energy_frequency, "lambda_windows": lambda_windows, + "rest2_scale_factors": rest2_scale_factors, "save_velocities": save_velocities, + "save_frame_on_exit": save_frame_on_exit, + "save_energy_on_exit": save_energy_on_exit, "auto_fix_minimise": auto_fix_minimise, + "num_energy_neighbours": num_energy_neighbours, + "null_energy": null_energy, } from concurrent.futures import ThreadPoolExecutor @@ -836,6 +969,17 @@ def run( if energy_frequency is not None: energy_frequency = u(energy_frequency) + if null_energy is not None: + null_energy = u(null_energy) + else: + null_energy = u("10000 kcal/mol") + + if num_energy_neighbours is not None: + try: + num_energy_neighbours = int(num_energy_neighbours) + except: + num_energy_neighbours = len(lambda_windows) + try: steps_to_run = int(time.to(picosecond) / self.timestep().to(picosecond)) except Exception: @@ -857,8 +1001,13 @@ def run( save_frequency = 25 else: save_frequency = save_frequency.to(picosecond) + no_save = False + else: + no_save = True + save_frequency = steps_to_run + 1 if energy_frequency != 0: + no_save_energy = False if energy_frequency is None: if self._map.specified("energy_frequency"): energy_frequency = ( @@ -866,10 +1015,15 @@ def run( ) else: energy_frequency = save_frequency + no_save_energy = no_save else: energy_frequency = energy_frequency.to(picosecond) + else: + energy_frequency = steps_to_run + 1 + no_save_energy = True if frame_frequency != 0: + no_save_frame = False if frame_frequency is None: if self._map.specified("frame_frequency"): frame_frequency = ( @@ -877,8 +1031,12 @@ def run( ) else: frame_frequency = save_frequency + no_save_frame = no_save else: frame_frequency = frame_frequency.to(picosecond) + else: + frame_frequency = steps_to_run + 1 + no_save_frame = True completed = 0 @@ -899,6 +1057,9 @@ def run( # Fixed lambda value. else: delta_lambda = None + + # Set the REST2 scaling factors. + rest2_scale_factors = [1.0, 1.0] else: delta_lambda = None if lambda_windows is not None: @@ -908,6 +1069,15 @@ def run( if self._map.specified("lambda_windows"): lambda_windows = self._map["lambda_windows"].value() + if rest2_scale_factors is not None: + if len(rest2_scale_factors) != len(lambda_windows): + raise ValueError( + "len(rest2_scale_factors) must be equal to len(lambda_windows)" + ) + else: + if lambda_windows is not None: + rest2_scale_factors = [1.0] * len(lambda_windows) + def runfunc(num_steps): try: integrator = self._omm_mols.getIntegrator() @@ -927,60 +1097,50 @@ def process_block(state, state_has_cv, nsteps_completed): } ) - def get_steps_till_save(completed: int, total: int): - """Internal function to calculate the number of steps - to run before the next save. This returns a tuple - of number of steps, and then if a frame should be - saved and if the energy should be saved - """ - if completed < 0: - completed = 0 - - if completed >= total: - return ( - 0, - frame_frequency_steps > 0, - energy_frequency_steps > 0, - ) - - elif frame_frequency_steps <= 0 and energy_frequency_steps <= 0: - return (total, False, False) - - n_to_end = total - completed - - if frame_frequency_steps > 0: - n_to_frame = min( - frame_frequency_steps - (completed % frame_frequency_steps), - n_to_end, - ) - else: - n_to_frame = total - completed - - if energy_frequency_steps > 0: - n_to_energy = min( - energy_frequency_steps - (completed % energy_frequency_steps), - n_to_end, - ) - else: - n_to_energy = total - completed - - if n_to_frame == n_to_energy: - return (n_to_frame, True, True) - elif n_to_frame < n_to_energy: - return (n_to_frame, True, False) - else: - return (n_to_energy, False, True) - - block_size = 50 - state = None state_has_cv = (False, False) saved_last_frame = False + # whether the energy or frame were saved after the current block + have_saved_frame = False + have_saved_energy = False + class NeedsMinimiseError(Exception): pass nsteps_before_run = self._current_step + # if this is the first call, then set the save frequencies + if nsteps_before_run == 0: + self._next_save_frame = frame_frequency_steps + self._next_save_energy = energy_frequency_steps + self._prev_frame_frequency_steps = frame_frequency_steps + self._prev_energy_frequency_steps = energy_frequency_steps + self._prev_no_frame = no_save_frame + self._prev_no_energy = no_save_energy + # handle adjustments to the save frequencies + else: + if frame_frequency_steps != self._prev_frame_frequency_steps: + if self._prev_no_frame: + self._next_save_frame = nsteps_before_run + frame_frequency_steps + else: + self._next_save_frame = ( + self._next_save_frame + + frame_frequency_steps + - self._prev_frame_frequency_steps + ) + if energy_frequency_steps != self._prev_energy_frequency_steps: + if self._prev_no_energy: + self._next_save_energy = nsteps_before_run + energy_frequency_steps + else: + self._next_save_energy = ( + self._next_save_energy + + energy_frequency_steps + - self._prev_energy_frequency_steps + ) + self._prev_no_frame = no_save_frame + self._prev_frame_frequency_steps = frame_frequency_steps + self._prev_no_energy = no_save_energy + self._prev_energy_frequency_steps = energy_frequency_steps from ..base import ProgressBar from ..units import second @@ -995,13 +1155,51 @@ class NeedsMinimiseError(Exception): with ThreadPoolExecutor() as pool: while completed < steps_to_run: - ( - nrun_till_save, - save_frame, - save_energy, - ) = get_steps_till_save(completed, steps_to_run) + block_size = 50 - assert nrun_till_save > 0 + steps_till_frame = self._next_save_frame - ( + completed + nsteps_before_run + ) + if steps_till_frame <= 0 or ( + steps_till_frame <= block_size + and steps_till_frame <= steps_to_run - completed + ): + save_frame = True + self._next_save_frame += frame_frequency_steps + if frame_frequency_steps < block_size: + block_size = frame_frequency_steps + else: + save_frame = False + + steps_till_energy = self._next_save_energy - ( + completed + nsteps_before_run + ) + if steps_till_energy <= 0 or ( + steps_till_energy <= block_size + and steps_till_energy <= steps_to_run - completed + ): + save_energy = True + self._next_save_energy += energy_frequency_steps + if energy_frequency_steps < block_size: + block_size = energy_frequency_steps + else: + save_energy = False + + # save the last frame if we're about to exit and the user + # has requested it + if ( + save_frame_on_exit + and completed + block_size >= steps_to_run + ): + save_frame = True + + # save the last energy if we're about to exit and the user + # has requested it + if ( + save_energy_on_exit + and completed + block_size >= steps_to_run + ): + save_energy = True self._enter_dynamics_block() @@ -1016,47 +1214,62 @@ class NeedsMinimiseError(Exception): else: process_promise = None - while nrun_till_save > 0: - nrun = block_size + nrun = block_size - if 2 * nrun > nrun_till_save: - nrun = nrun_till_save + # this block will exceed the run time so reduce the size + if nrun > steps_to_run - completed: + nrun = steps_to_run - completed - # run the current block in the background - run_promise = pool.submit(runfunc, nrun) + # run the current block in the background + run_promise = pool.submit(runfunc, nrun) - result = None + result = None - while not run_promise.done(): - try: - result = run_promise.result(timeout=1.0) - except Exception: - pass - - # catch rare edge case where the promise timed - # out, but then completed before the .done() - # test in the next loop iteration - if result is None: - result = run_promise.result() - - if result == 0: - completed += nrun - nrun_till_save -= nrun - progress.set_progress(completed) - run_promise = None - else: - # make sure we finish processing the last block - if process_promise is not None: - try: - process_promise.result() - except Exception: - pass - - if completed == 0 and auto_fix_minimise: + while not run_promise.done(): + try: + result = run_promise.result(timeout=1.0) + except Exception as e: + if ( + "NaN" in str(e) + and not have_saved_frame + and not have_saved_energy + and auto_fix_minimise + ): raise NeedsMinimiseError() - # something went wrong - re-raise the exception - raise result + # catch rare edge case where the promise timed + # out, but then completed before the .done() + # test in the next loop iteration + if result is None: + result = run_promise.result() + + if result == 0: + completed += nrun + progress.set_progress(completed) + run_promise = None + else: + # make sure we finish processing the last block + if process_promise is not None: + try: + process_promise.result() + except Exception as e: + if ( + "NaN" in str(e) + and not have_saved_frame + and not have_saved_energy + and auto_fix_minimise + ): + raise NeedsMinimiseError() + + if ( + not have_saved_frame + and not have_saved_energy + and auto_fix_minimise + ): + raise NeedsMinimiseError() + + # something went wrong - re-raise the exception + raise result # make sure we've finished processing the last block if process_promise is not None: @@ -1069,12 +1282,18 @@ class NeedsMinimiseError(Exception): save_frame=save_frame, save_energy=save_energy, lambda_windows=lambda_windows, + rest2_scale_factors=rest2_scale_factors, save_velocities=save_velocities, delta_lambda=delta_lambda, + num_energy_neighbours=num_energy_neighbours, + null_energy=null_energy.value(), ) saved_last_frame = False + have_saved_frame = save_frame + have_saved_energy = save_energy + kinetic_energy = state.getKineticEnergy().value_in_unit( openmm.unit.kilocalorie_per_mole ) @@ -1086,7 +1305,11 @@ class NeedsMinimiseError(Exception): state = None saved_last_frame = True - if completed == 0 and auto_fix_minimise: + if ( + not have_saved_frame + and not have_saved_energy + and auto_fix_minimise + ): raise NeedsMinimiseError() raise RuntimeError( @@ -1331,8 +1554,13 @@ def run( frame_frequency=None, energy_frequency=None, lambda_windows=None, + rest2_scale_factors=None, save_velocities: bool = None, + save_frame_on_exit: bool = False, + save_energy_on_exit: bool = False, auto_fix_minimise: bool = True, + num_energy_neighbours: int = None, + null_energy: str = None, ): """ Perform dynamics on the molecules. @@ -1389,17 +1617,51 @@ def run( we always save the potential energy of the simulated lambda value, even if it is not in the list of lambda windows. + rest2_scale_factors: list[float] + The scaling factors for the REST2 region for each lambda + window. + save_velocities: bool Whether or not to save the velocities when running dynamics. By default this is False. Set this to True if you are interested in saving the velocities. + save_frame_on_exit: bool + Whether to save a trajectory frame on exit, regardless of + whether the frame frequency has been reached. + + save_energy_on_exit: bool + Whether to save the energy on exit, regardless of whether + the energy frequency has been reached. + auto_fix_minimise: bool Whether or not to automatically run minimisation if the trajectory exits with an error in the first few steps. Such failures often indicate that the system needs minimsing. This automatically runs the minimisation in these cases, and then runs the requested dynamics. + + num_energy_neighbours: int + The number of neighbouring windows to use when computing + the energy trajectory for the simulation lambda value. + This can be used to compute energies over a subset of the + values in 'lambda_windows', hence reducing the cost of + computing the energy trajectory. Note that the simulation + lambda value must be contained in 'lambda_windows', so it + is recommended that the values are rounded. A value of + 'null_energy' will be added to the energy trajectory for the + lambda windows that are omitted. Note that a similar result + can be achieved by simply removing any lambda values from + 'lambda_windows' that you don't want, but this will result + in an energy trajectory that only contains results for + the specified lambda values. If None, then all lambda windows + will be used. + + null_energy: str + The energy value to use for lambda windows that are not + being computed as part of the energy trajectory, i.e. when + 'num_energy_neighbours' is less than len(lambda_windows). + By default, a value of '10000 kcal mol-1' is used. """ if not self._d.is_null(): if save_velocities is None: @@ -1414,8 +1676,13 @@ def run( frame_frequency=frame_frequency, energy_frequency=energy_frequency, lambda_windows=lambda_windows, + rest2_scale_factors=rest2_scale_factors, save_velocities=save_velocities, + save_frame_on_exit=save_frame_on_exit, + save_energy_on_exit=save_energy_on_exit, auto_fix_minimise=auto_fix_minimise, + num_energy_neighbours=num_energy_neighbours, + null_energy=null_energy, ) return self @@ -1444,13 +1711,22 @@ def get_lambda(self): """ return self._d.get_lambda() - def set_lambda(self, lambda_value: float, update_constraints: bool = True): + def set_lambda( + self, + lambda_value: float, + rest2_scale: float = 1.0, + update_constraints: bool = True, + ): """ Set the current value of lambda for this system. This will update the forcefield parameters in the context according to the data in the LambdaSchedule. This does nothing if this isn't a perturbable system. + The `rest2_scale` parameter specifies the temperature of the + REST2 region relative to the rest of the system at the specified + lambda value. + If `update_constraints` is True, then this will also update the constraint length of any constrained perturbable bonds. These will be set to the r0 value for that bond at this @@ -1458,9 +1734,25 @@ def set_lambda(self, lambda_value: float, update_constraints: bool = True): the constraint will not be changed. """ self._d.set_lambda( - lambda_value=lambda_value, update_constraints=update_constraints + lambda_value=lambda_value, + rest2_scale=rest2_scale, + update_constraints=update_constraints, ) + def get_rest2_scale(self): + """ + Return the current REST2 scaling factor. + """ + if self.is_null(): + return None + return self._d.get_rest2_scale() + + def set_rest2_scale(self, rest2_scale: float): + """ + Set the current REST2 scaling factor. + """ + self._d.set_rest2_scale(rest2_scale=rest2_scale) + def ensemble(self): """ Return the ensemble in which the simulation is being performed @@ -1646,35 +1938,56 @@ def current_energy(self): """ return self._d.current_energy() - def current_potential_energy(self, lambda_values=None): + def current_potential_energy(self, lambda_values=None, rest2_scale_factors=None): """ Return the current potential energy. If `lambda_values` is passed (which should be a list of lambda values) then this will return the energies (as a list) at the requested lambda values + + If `rest2_scale_factors` is passed, then these will be + used to scale the temperature of the REST2 region at each + lambda value. """ if lambda_values is None: return self._d.current_potential_energy() else: if not isinstance(lambda_values, list): lambda_values = [lambda_values] + if rest2_scale_factors is None: + rest2_scale_factors = [1.0] * len(lambda_values) + else: + if not isinstance(rest2_scale_factors, list): + rest2_scale_factors = [rest2_scale_factors] + else: + if len(rest2_scale_factors) != len(lambda_values): + raise ValueError( + "len(rest2_scale_factors) must be equal to len(lambda_values)" + ) # save the current value of lambda so we # can restore it old_lambda = self.get_lambda() + old_rest2_scale = self.get_rest2_scale() nrgs = [] try: - for lambda_value in lambda_values: - self.set_lambda(lambda_value, update_constraints=False) + for lambda_value, rest2_scale in zip( + lambda_values, rest2_scale_factors + ): + self.set_lambda( + lambda_value, rest2_scale=rest2_scale, update_constraints=False + ) nrgs.append(self._d.current_potential_energy()) except Exception: - self.set_lambda(old_lambda, update_constraints=False) + self.set_lambda(old_lambda, old_rest2_scale, update_constraints=False) raise - self.set_lambda(old_lambda, update_constraints=False) + self.set_lambda( + old_lambda, rest2_scale=old_rest2_scale, update_constraints=False + ) return nrgs diff --git a/src/sire/mol/_minimisation.py b/src/sire/mol/_minimisation.py index 32e3dc8f9..b6c89b9fe 100644 --- a/src/sire/mol/_minimisation.py +++ b/src/sire/mol/_minimisation.py @@ -17,6 +17,7 @@ def __init__( cutoff_type=None, schedule=None, lambda_value=None, + rest2_scale=None, swap_end_states=None, ignore_perturbations=None, shift_delta=None, @@ -37,6 +38,7 @@ def __init__( _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, "rest2_scale", rest2_scale) if shift_delta is not None: _add_extra(extras, "shift_delta", u(shift_delta)) diff --git a/src/sire/morph/_perturbation.py b/src/sire/morph/_perturbation.py index 6f6698a5e..4c7bab201 100644 --- a/src/sire/morph/_perturbation.py +++ b/src/sire/morph/_perturbation.py @@ -840,7 +840,7 @@ def to_openmm( from ..convert.openmm import PerturbableOpenMMMolecule try: - return PerturbableOpenMMMolecule(self._mol.molecule(), map=map) + return PerturbableOpenMMMolecule(self._mol.molecule(), 0, map=map) except KeyError: # probably need to choose an end-state - default to reference return PerturbableOpenMMMolecule( diff --git a/src/sire/qm/_emle.py b/src/sire/qm/_emle.py index 591ac599d..a47ebaefa 100644 --- a/src/sire/qm/_emle.py +++ b/src/sire/qm/_emle.py @@ -38,7 +38,10 @@ def get_forces(self): ) # Get the OpenMM EMLE force. - emle_force = _deepcopy(d._d._omm_mols.getSystem().getForce(0)) + for force in d._d._omm_mols.getSystem().getForces(): + if "QMForce" in force.getName(): + break + emle_force = _deepcopy(force) # Create a null CustomBondForce to add the EMLE interpolation # parameter. @@ -84,7 +87,10 @@ def get_forces(self): ) # Get the OpenMM EMLE force. - emle_force = _deepcopy(d._d._omm_mols.getSystem().getForce(0)) + for force in d._d._omm_mols.getSystem().getForces(): + if "QMForce" in force.getName(): + break + emle_force = _deepcopy(force) # Create a null CustomBondForce to add the EMLE interpolation # parameter. @@ -225,23 +231,10 @@ def emle( # Create an engine from an EMLE calculator. if isinstance(calculator, _EMLECalculator): - # Determine the callback name. Use an optimised version of the callback - # if the user has specified "torchani" as the backend and is using - # "electrostatic" embedding. - if calculator._backend == "torchani" and calculator._method == "electrostatic": - try: - from emle.models import ANI2xEMLE as _ANI2xEMLE - - callback = "_sire_callback_optimised" - except: - callback = "_sire_callback" - else: - callback = "_sire_callback" - # Create the EMLE engine. engine = EMLEEngine( calculator, - callback, + "_sire_callback", cutoff, neighbour_list_frequency, False, diff --git a/src/sire/restraints/__init__.py b/src/sire/restraints/__init__.py index 5f5af1465..166cb985c 100644 --- a/src/sire/restraints/__init__.py +++ b/src/sire/restraints/__init__.py @@ -1,4 +1,4 @@ -__all__ = ["positional", "bond", "distance", "boresch", "get_standard_state_correction"] +__all__ = ["angle", "positional", "bond", "dihedral", "distance", "boresch", "get_standard_state_correction"] -from ._restraints import bond, boresch, distance, positional +from ._restraints import angle, bond, boresch, dihedral, distance, positional from ._standard_state_correction import get_standard_state_correction diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 077cb97c0..6e4d5d050 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -1,6 +1,8 @@ __all__ = [ + "angle", "boresch", "bond", + "dihedral", "distance", "positional", ] @@ -18,6 +20,96 @@ def _to_atoms(mols, atoms): return selection_to_atoms(mols, atoms) +def angle(mols, atoms, theta0=None, ktheta=None, name=None, map=None): + """ + Create a set of anglel restraints from all of the atoms in 'atoms' + where all atoms are contained in the container 'mols', using the + passed values of the force constant 'ktheta' and equilibrium + angle value theta0. + + If theta0 is None, then the current angle for + provided atoms will be used as the equilibium value. + + If ktheta is None, then a default value of 100 kcal mol-1 rad-2 will be used + + Parameters + ---------- + mols : sire.system._system.System + The system containing the atoms. + + atoms : SireMol::Selector + The atoms to restrain. + + ktheta : str or SireUnits::Dimension::GeneralUnit or, optional + The force constants for the angle restraints. + If None, this will default to 100 kcal mol-1 rad-2. + Default is None. + + theta0 : str or SireUnits::Dimension::GeneralUnit, optional + The equilibrium angles for the angle restraints. If None, these + will be measured from the current coordinates of the atoms. + Default is None. + + Returns + ------- + AnglelRestraints : SireMM::AngleRestraints + A container of angle restraints, where the first restraint is + the AngleRestraint created. The angle restraint created can be + extracted with AngleRestraints[0]. + """ + from .. import u + from ..base import create_map + from ..mm import AngleRestraint, AngleRestraints + + map = create_map(map) + map_dict = map.to_dict() + ktheta = ktheta if ktheta is not None else map_dict.get("ktheta", None) + theta0 = theta0 if theta0 is not None else map_dict.get("theta0", None) + name = name if name is not None else map_dict.get("name", None) + + atoms = _to_atoms(mols, atoms) + + if len(atoms) != 3: + raise ValueError( + "You need to provide 3 atoms to create an angle restraint" + f"whereas {len(atoms)} atoms were provided." + ) + + from .. import measure + + if ktheta is None: + ktheta = u("100 kcal mol-1 rad-2") + + elif type(ktheta) is list: + raise NotImplementedError( + "Setup of multiple angle restraints simultaneously is not currently supported. " + "Please set up each restraint individually and then combine them into multiple restraints." + ) + + if theta0 is None: + from .. import measure + + theta0 = measure(atoms[0], atoms[1], atoms[2]) + + elif type(theta0) is list: + raise NotImplementedError( + "Setup of multiple angle restraints simultaneously is not currently supported. " + "Please set up each restraint individually and then combine them into multiple restraints." + ) + else: + theta0 = u(theta0) + + mols = mols.atoms() + + if name is None: + restraints = AngleRestraints() + else: + restraints = AngleRestraints(name=name) + + restraints.add(AngleRestraint(mols.find(atoms), theta0, ktheta)) + return restraints + + def boresch( mols, receptor, @@ -377,6 +469,96 @@ def _check_stability_boresch_restraint(restraint_components, temperature=u("298 ) +def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): + """ + Create a set of dihedral restraints from all of the atoms in 'atoms' + where all atoms are contained in the container 'mols', using the + passed values of the force constant 'kphi' and equilibrium + torsion angle phi0. + + If phi0 is None, then the current torsional angle for + the provided atoms will be used as the equilibium value. + + If kphi is None, then a default value of 100 kcal mol-1 rad-2 will be used + + Parameters + ---------- + mols : sire.system._system.System + The system containing the atoms. + + atoms : SireMol::Selector + The atoms to restrain. + + kphi : str or SireUnits::Dimension::GeneralUnit or, optional + The force constants for the torsion restraints. + If None, this will default to 100 kcal mol-1 rad-2. + Default is None. + + phi0 : str or SireUnits::Dimension::GeneralUnit, optional + The equilibrium torsional angle for restraints. If None, these + will be measured from the current coordinates of the atoms. + Default is None. + + Returns + ------- + DihedralRestraints : SireMM::DihedralRestraints + A container of Dihedral restraints, where the first restraint is + the DihedralRestraint created. The Dihedral restraint created can be + extracted with DihedralRestraints[0]. + """ + from .. import u + from ..base import create_map + from ..mm import DihedralRestraint, DihedralRestraints + + map = create_map(map) + map_dict = map.to_dict() + kphi = kphi if kphi is not None else map_dict.get("kphi", None) + phi0 = phi0 if phi0 is not None else map_dict.get("phi0", None) + name = name if name is not None else map_dict.get("name", None) + + atoms = _to_atoms(mols, atoms) + + if len(atoms) != 4: + raise ValueError( + "You need to provide 4 atoms to create a dihedral restraint" + f"whereas {len(atoms)} atoms were provided." + ) + + from .. import measure + + if kphi is None: + kphi = u("100 kcal mol-1 rad-2") + + elif type(kphi) is list: + raise NotImplementedError( + "Setup of multiple dihedral restraints simultaneously is not currently supported. " + "Please set up each restraint individually and then combine them into multiple restraints." + ) + + if phi0 is None: + from .. import measure + + phi0 = measure(atoms[0], atoms[1], atoms[2], atoms[3]) + + elif type(phi0) is list: + raise NotImplementedError( + "Setup of multiple dihedral restraints simultaneously is not currently supported. " + "Please set up each restraint individually and then combine them into multiple restraints." + ) + else: + phi0 = u(phi0) + + mols = mols.atoms() + + if name is None: + restraints = DihedralRestraints() + else: + restraints = DihedralRestraints(name=name) + + restraints.add(DihedralRestraint(mols.find(atoms), phi0, kphi)) + return restraints + + def distance(mols, atoms0, atoms1, r0=None, k=None, name=None, map=None): """ Create a set of distance restraints from all of the atoms in 'atoms0' diff --git a/src/sire/system/_system.py b/src/sire/system/_system.py index e2bd72a7f..0846f535b 100644 --- a/src/sire/system/_system.py +++ b/src/sire/system/_system.py @@ -592,6 +592,23 @@ def dynamics(self, *args, **kwargs): microcanonical (NVE) or canonical (NVT) simulation will be run if the pressure is not set. + rest2_scale: float + The scaling factor to apply when running a REST2 simulation + (Replica Exchange with Solute Tempering). This defaults to 1.0. + This specifies the ratio of the temperature of the REST2 region + to the temperature of the rest of the system. + + rest2_selection: str + A selection string for atoms to include in the REST2 region + in addition to any perturbable molecules. For example, + "molidx 0 and residx 0,1,2" would select atoms from the first + three residues of the first molecule. If None, then all atoms + within perturbable molecules will be included in the REST2 + region. When atoms within a perturbable molecule are also + included in the selection, then only those atoms will be + considered as part of the REST2 region. This allows REST2 to + be applied to protein mutations. + vacuum: bool Whether or not to run the simulation in vacuum. If this is set to `True`, then the simulation space automatically be diff --git a/tests/cas/test_lambdaschedule.py b/tests/cas/test_lambdaschedule.py index 92dddb70b..5348f4611 100644 --- a/tests/cas/test_lambdaschedule.py +++ b/tests/cas/test_lambdaschedule.py @@ -145,3 +145,11 @@ def test_lambdaschedule(): ) _assert_same_equation(l.lam(), l.get_equation(stage="scale_up"), morph5_equation) + +@pytest.mark.parametrize("force, lever, contained", [ + ("ghost-14", "kappa", True), + ("ghost-14", "epsilon", True) +]) +def test_has_force_specific_equation(force, lever, contained): + l = sr.cas.LambdaSchedule.standard_decouple() + assert l.has_force_specific_equation("decouple", force, lever) == contained \ No newline at end of file diff --git a/tests/convert/test_openmm_constraints.py b/tests/convert/test_openmm_constraints.py index c7427391a..ef98bb276 100644 --- a/tests/convert/test_openmm_constraints.py +++ b/tests/convert/test_openmm_constraints.py @@ -467,7 +467,6 @@ def test_asymmetric_constraints(merged_ethane_methanol): # different constraints. from math import isclose - from openmm import XmlSerializer from tempfile import NamedTemporaryFile # Extract the molecule. @@ -497,9 +496,9 @@ def test_asymmetric_constraints(merged_ethane_methanol): xml0 = NamedTemporaryFile() xml1 = NamedTemporaryFile() with open(xml0.name, "w") as f: - f.write(XmlSerializer.serialize(d_forwards._d._omm_mols.getSystem())) + f.write(d_forwards.to_xml()) with open(xml1.name, "w") as f: - f.write(XmlSerializer.serialize(d_backwards._d._omm_mols.getSystem())) + f.write(d_backwards.to_xml()) # Load the serialised systems and sort. with open(xml0.name, "r") as f: diff --git a/tests/convert/test_openmm_rest2.py b/tests/convert/test_openmm_rest2.py new file mode 100644 index 000000000..7ac498a52 --- /dev/null +++ b/tests/convert/test_openmm_rest2.py @@ -0,0 +1,268 @@ +from math import isclose + +import pytest + +import sire as sr + + +@pytest.fixture(scope="module") +def toluene_methane(): + """ + Load the toluene methane perturbable system. + """ + return sr.load_test_files("toluene_methane.s3") + + +@pytest.mark.parametrize( + ["mols", "rest2_selection", "excluded_atoms"], + [ + ("toluene_methane", None, []), + ("toluene_methane", "not atomidx 0,1,2,3", [0, 1, 2, 3]), + ("ala_mols", "molidx 0", []), + ], +) +def test_rest2(mols, rest2_selection, excluded_atoms, request): + """ + Test that REST2 modifications are correctly applied to the system. + """ + + # Load the test perturbation. + mol = request.getfixturevalue(mols)[0] + + # Link to the reference state. + try: + mol = sr.morph.link_to_reference(mol) + is_perturbable = True + except: + is_perturbable = False + pass + + # Work out the number of dihedrals in the system. + num_dihedrals = len(mol.dihedrals()) + + # Create a dynamics object. + d = mol.dynamics(platform="Reference", rest2_selection=rest2_selection) + + # Extract the OpenMM system. + omm_system = d.context().getSystem() + + # Find the PeriodicTorsionForce. + for force in omm_system.getForces(): + if force.getName() == "PeriodicTorsionForce": + break + + # Store the initial parameters. + torsion_params_initial = [] + excluded_torsion_indices = [] + for x in range(num_dihedrals): + i, j, k, l, periodicity, phase, force_constant = force.getTorsionParameters(x) + torsion_params_initial.append(force_constant) + # This torsion is not in the REST2 region and should be excluded. + if ( + i in excluded_atoms + or j in excluded_atoms + or k in excluded_atoms + or l in excluded_atoms + ): + excluded_torsion_indices.append(x) + + # Find the NonbondedForce. + for force in omm_system.getForces(): + if force.getName() == "NonbondedForce": + break + + # Store the initial parameters. + nonbonded_params_initial = [] + excluded_nonbonded_indices = [] + for i in range(force.getNumParticles()): + charge, sigma, epsilon = force.getParticleParameters(i) + nonbonded_params_initial.append((charge, epsilon)) + if i in excluded_atoms: + excluded_nonbonded_indices.append(i) + exception_params_initial = [] + excluded_exceptions = [] + for i in range(force.getNumExceptions()): + x, y, chargeProd, sigma, epsilon = force.getExceptionParameters(i) + exception_params_initial.append((chargeProd, epsilon)) + if x in excluded_atoms or y in excluded_atoms: + excluded_exceptions.append(i) + + # Handle custom forces for pertubable molecules. + if is_perturbable: + # Find the ghost/ghost nonbonded interaction. + for force in omm_system.getForces(): + if force.getName() == "GhostGhostNonbondedForce": + break + + # Store the initial parameters. + ghost_ghost_params_initial = [] + excluded_ghost_ghost_indices = [] + for i in range(force.getNumParticles()): + charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( + force.getParticleParameters(i) + ) + ghost_ghost_params_initial.append((charge, two_sqrt_epsilon)) + if i in excluded_atoms: + excluded_ghost_ghost_indices.append(i) + + # Find the ghost/non-ghost nonbonded interaction. + for force in omm_system.getForces(): + if force.getName() == "GhostNonGhostNonbondedForce": + break + + # Store the initial parameters. + ghost_non_ghost_params_initial = [] + excluded_ghost_non_ghost_indices = [] + for i in range(force.getNumParticles()): + charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( + force.getParticleParameters(i) + ) + ghost_non_ghost_params_initial.append((charge, two_sqrt_epsilon)) + if i in excluded_atoms: + excluded_ghost_non_ghost_indices.append(i) + + # Update the REST2 scaling factor. + d.set_lambda(0.0, rest2_scale=2.0) + + # Extract the OpenMM system. + omm_system = d.context().getSystem() + + # Find the PeriodicTorsionForce. + for force in omm_system.getForces(): + if force.getName() == "PeriodicTorsionForce": + break + + # Store the modified parameters. + torsion_params_modified = [] + for i in range(force.getNumTorsions()): + torsion_params_modified.append(force.getTorsionParameters(i)[-1]) + + # Find the NonbondedForce. + for force in omm_system.getForces(): + if force.getName() == "NonbondedForce": + break + + # Store the modified parameters. + nonbonded_params_modified = [] + for i in range(force.getNumParticles()): + charge, sigma, epsilon = force.getParticleParameters(i) + nonbonded_params_modified.append((charge, epsilon)) + exception_params_modified = [] + for i in range(force.getNumExceptions()): + exception_params_modified.append(force.getExceptionParameters(i)[-3::2]) + + # Find the ghost/ghost nonbonded interaction. + for force in omm_system.getForces(): + if force.getName() == "GhostGhostNonbondedForce": + break + + # Handle custom forces for pertubable molecules. + if is_perturbable: + # Store the modified parameters. + ghost_ghost_params_modified = [] + for i in range(force.getNumParticles()): + charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( + force.getParticleParameters(i) + ) + ghost_ghost_params_modified.append((charge, two_sqrt_epsilon)) + + # Find the ghost/non-ghost nonbonded interaction. + for force in omm_system.getForces(): + if force.getName() == "GhostNonGhostNonbondedForce": + break + + # Store the modified parameters. + ghost_non_ghost_params_modified = [] + for i in range(force.getNumParticles()): + charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( + force.getParticleParameters(i) + ) + ghost_non_ghost_params_modified.append((charge, two_sqrt_epsilon)) + + # Store the scaling factor. + scale = 0.5 + + # All dihedral force constants should be halved. + for i in range(num_dihedrals): + if i in excluded_torsion_indices: + # This torsion is not in the REST2 region so is unchanged. + assert isclose( + torsion_params_modified[i]._value, torsion_params_initial[i]._value + ) + else: + assert isclose( + torsion_params_modified[i]._value, + torsion_params_initial[i]._value * scale, + ) + + # All impropers should be unchanged. + for i in range(num_dihedrals, len(torsion_params_initial)): + assert isclose( + torsion_params_modified[i]._value, torsion_params_initial[i]._value + ) + + # Nonbonded charges should be scaled by the square root of the scaling + # factor and epsilon should be scaled by the scaling factor. + for i in range(len(nonbonded_params_initial)): + charge, epsilon = nonbonded_params_initial[i] + charge_modified, epsilon_modified = nonbonded_params_modified[i] + if i in excluded_nonbonded_indices: + # This atom is not in the REST2 region so is unchanged. + assert isclose(charge_modified._value, charge._value) + assert isclose(epsilon_modified._value, epsilon._value) + else: + assert isclose(charge_modified._value, charge._value * scale**0.5) + assert isclose(epsilon_modified._value, epsilon._value * scale) + + # For exceptions, both the charge and epsilon should be scaled by the + # scaling factor. + for i in range(len(exception_params_initial)): + charge, epsilon = exception_params_initial[i] + charge_modified, epsilon_modified = exception_params_modified[i] + if i in excluded_exceptions: + # This exception is not in the REST2 region so is unchanged. + assert isclose(charge_modified._value, charge._value) + assert isclose(epsilon_modified._value, epsilon._value) + else: + assert isclose(charge_modified._value, charge._value * scale) + if epsilon._value > 1e-6: + assert isclose(epsilon_modified._value, epsilon._value * scale) + + if is_perturbable: + # Ghost/ghost nonbonded charges should be scaled by the square root of the + # scaling factor and epsilon should be scaled by the scaling factor. + # (Note that epsilon is stored as sqrt(epsilon) so the scaling factor is + # also square rooted.) + for i in range(len(ghost_ghost_params_initial)): + charge, two_sqrt_epsilon = ghost_ghost_params_initial[i] + charge_modified, two_sqrt_epsilon_modified = ghost_ghost_params_modified[i] + if i in excluded_ghost_ghost_indices: + # This atom is not in the REST2 region so is unchanged. + assert isclose(charge_modified, charge) + assert isclose(two_sqrt_epsilon_modified, two_sqrt_epsilon) + else: + assert isclose(charge_modified, charge * scale**0.5) + if two_sqrt_epsilon > 1e-6: + assert isclose( + two_sqrt_epsilon_modified, two_sqrt_epsilon * scale**0.5 + ) + + # Ghost/non-ghost nonbonded charges should be scaled by the square root of + # the scaling factor and epsilon should be scaled by the scaling factor. + # (Note that epsilon is stored as sqrt(epsilon) so the scaling factor is + # also square rooted.) + for i in range(len(ghost_non_ghost_params_initial)): + charge, two_sqrt_epsilon = ghost_non_ghost_params_initial[i] + charge_modified, two_sqrt_epsilon_modified = ( + ghost_non_ghost_params_modified[i] + ) + if i in excluded_ghost_non_ghost_indices: + # This atom is not in the REST2 region so is unchanged. + assert isclose(charge_modified, charge) + assert isclose(two_sqrt_epsilon_modified, two_sqrt_epsilon) + else: + assert isclose(charge_modified, charge * scale**0.5) + if two_sqrt_epsilon > 1e-6: + assert isclose( + two_sqrt_epsilon_modified, two_sqrt_epsilon * scale**0.5 + ) diff --git a/tests/io/test_sdf.py b/tests/io/test_sdf.py new file mode 100644 index 000000000..5631caedc --- /dev/null +++ b/tests/io/test_sdf.py @@ -0,0 +1,8 @@ +import sire as sr + + +def test_cis_or_trans(): + """ + Make sure that we can read an SDF file containing a cis or trans double bond. + """ + mols = sr.load_test_files("cis_trans_double_bond.sdf") diff --git a/tests/mol/test_dynamics.py b/tests/mol/test_dynamics.py index 78842ee75..7b86fad4a 100644 --- a/tests/mol/test_dynamics.py +++ b/tests/mol/test_dynamics.py @@ -75,3 +75,73 @@ def test_cutoff_options(ala_mols): assert d.info().cutoff() == sr.u("7.5 A") assert d.info().cutoff_type() == "PME" + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_sample_frequency(ala_mols): + """ + Test that energies and frames are saved at the correct frequency. + """ + + from sire.base import ProgressBar + + ProgressBar.set_silent() + + mols = ala_mols + + d = mols.dynamics(platform="Reference", timestep="1 fs") + + # Create a list of lambda windows. + lambdas = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] + + # Run 10 cycles of dynamics, saving energies every 2 fs and frames every 5 fs. + for i in range(10): + d.run( + "1 fs", + energy_frequency="2 fs", + frame_frequency="5 fs", + lambda_windows=lambdas, + ) + + # Get the energy trajectory. + nrg_traj = d.energy_trajectory() + + # Make sure the trajectory has 5 frames. + assert len(nrg_traj) == 5 + + # Get the updated system. + new_mols = d.commit() + + # Check that the trajectory has 2 frames. + assert new_mols.num_frames() == 2 + + # Now check when we request that a trajectory frame is saved when run exits. + + # Recreate the dynamics object. + d = mols.dynamics(platform="Reference", timestep="1 fs") + + # Run 10 cycles of dynamics, saving energies frames on exit. + for i in range(10): + d.run( + "1 fs", + energy_frequency="2 fs", + frame_frequency="5 fs", + lambda_windows=lambdas, + save_frame_on_exit=True, + save_energy_on_exit=True, + ) + + # Get the energy trajectory. + nrg_traj = d.energy_trajectory() + + # Make sure the trajectory has 10 frames. + assert len(nrg_traj) == 10 + + # Get the updated system. + new_mols = d.commit() + + # Check that the trajectory has 10 frames. + assert new_mols.num_frames() == 10 diff --git a/tests/restraints/test_angle_dihedral_restraints.py b/tests/restraints/test_angle_dihedral_restraints.py new file mode 100644 index 000000000..f6096ec21 --- /dev/null +++ b/tests/restraints/test_angle_dihedral_restraints.py @@ -0,0 +1,105 @@ +import pytest +import sire as sr + + +def test_default_angle_restraints_setup(ala_mols): + """Tests that angle restraints can be set up correctly with default parameters.""" + mols = ala_mols.clone() + restraints = sr.restraints.angle(mols=mols, atoms=[0, 1, 2]) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2] + assert restraints[0].ktheta().value() == 100.0 + + +def test_angle_restraint_custom_ktheta(ala_mols): + """Tests that angle restraints can be set up correctly with custom ktheta.""" + mols = ala_mols.clone() + restraints = sr.restraints.angle( + mols=mols, atoms=[0, 1, 2], ktheta="10 kcal mol-1 rad-2" + ) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2] + assert restraints[0].ktheta().value() == 10.0 + + +def test_angle_restraint_custom_ktheta_and_theta0(ala_mols): + """Tests that angle restraints can be set up correctly with custom ktheta and theta0.""" + mols = ala_mols.clone() + restraints = sr.restraints.angle( + mols=mols, atoms=[0, 1, 2], ktheta="10 kcal mol-1 rad-2", theta0="2 rad" + ) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2] + assert restraints[0].ktheta().value() == 10.0 + assert restraints[0].theta0().value() == 2.0 + + +def test_angle_restraint_molecular_container(ala_mols): + """Tests that angle restraints can be set up correctly with a molecular container (i.e. passing the angle of the molecule).""" + mols = ala_mols.clone() + ang = mols.angles()[0] + restraints = sr.restraints.angle(mols=mols, atoms=ang.atoms()) + + +def test_default_dihedral_restraints_setup(ala_mols): + """Tests that dihedral restraints can be set up correctly with default parameters.""" + mols = ala_mols.clone() + restraints = sr.restraints.dihedral(mols=mols, atoms=[0, 1, 2, 3]) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2, 3] + assert restraints[0].kphi().value() == 100.0 + + +def test_dihedral_restraint_custom_kphi(ala_mols): + """Tests that dihedral restraints can be set up correctly with custom kphi.""" + mols = ala_mols.clone() + restraints = sr.restraints.dihedral( + mols=mols, atoms=[0, 1, 2, 3], kphi="10 kcal mol-1 rad-2" + ) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2, 3] + assert restraints[0].kphi().value() == 10.0 + + +def test_dihedral_restraint_custom_kphi_and_phi0(ala_mols): + """Tests that dihedral restraints can be set up correctly with custom kphi and phi0.""" + mols = ala_mols.clone() + restraints = sr.restraints.dihedral( + mols=mols, atoms=[0, 1, 2, 3], kphi="10 kcal mol-1 rad-2", phi0="2 rad" + ) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2, 3] + assert restraints[0].kphi().value() == 10.0 + assert restraints[0].phi0().value() == 2.0 + + +def test_dihedral_restraint_molecular_container(ala_mols): + """Tests that dihedral restraints can be set up correctly with a molecular container (i.e. passing the dihedral of the molecule).""" + mols = ala_mols.clone() + dih = mols.dihedrals()[0] + restraints = sr.restraints.dihedral(mols=mols, atoms=dih.atoms()) + + +def test_multiple_angle_restraints(ala_mols): + """Tests that multiple angle restraints can be set up correctly.""" + mols = ala_mols.clone() + ang0 = mols.angles()[0] + ang1 = mols.angles()[1] + restraints = sr.restraints.angle(mols=mols, atoms=ang0.atoms()) + restraint1 = sr.restraints.angle(mols=mols, atoms=ang1.atoms()) + restraints.add(restraint1) + assert restraints.num_restraints() == 2 + + +def test_multiple_dihedral_restraints(ala_mols): + """Tests that multiple dihedral restraints can be set up correctly.""" + mols = ala_mols.clone() + dih0 = mols.dihedrals()[0] + dih1 = mols.dihedrals()[1] + dih2 = mols.dihedrals()[2] + restraints = sr.restraints.dihedral(mols=mols, atoms=dih0.atoms()) + restraint1 = sr.restraints.dihedral(mols=mols, atoms=dih1.atoms()) + restraint2 = sr.restraints.dihedral(mols=mols, atoms=dih2.atoms()) + restraints.add(restraint1) + restraints.add(restraint2) + assert restraints.num_restraints() == 3 diff --git a/version.txt b/version.txt index c60e0c0d6..3f841c887 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -2024.3.1 +2024.4.0 diff --git a/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp b/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp index 54e828a53..891398949 100644 --- a/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp @@ -205,7 +205,7 @@ void register_LambdaLever_class(){ } { //::SireOpenMM::LambdaLever::setExceptionIndicies - typedef void ( ::SireOpenMM::LambdaLever::*setExceptionIndicies_function_type)( int,::QString const &,::QVector< boost::tuples::tuple< int, int, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > > const & ) ; + typedef void ( ::SireOpenMM::LambdaLever::*setExceptionIndicies_function_type)( int,::QString const &,::QVector< boost::tuples::tuple< int, int > > const & ) ; setExceptionIndicies_function_type setExceptionIndicies_function_value( &::SireOpenMM::LambdaLever::setExceptionIndicies ); LambdaLever_exposer.def( @@ -231,13 +231,13 @@ void register_LambdaLever_class(){ } { //::SireOpenMM::LambdaLever::setLambda - typedef double ( ::SireOpenMM::LambdaLever::*setLambda_function_type)( ::OpenMM::Context &,double,bool ) const; + typedef double ( ::SireOpenMM::LambdaLever::*setLambda_function_type)( ::OpenMM::Context &,double,double,bool ) const; setLambda_function_type setLambda_function_value( &::SireOpenMM::LambdaLever::setLambda ); LambdaLever_exposer.def( "setLambda" , setLambda_function_value - , ( bp::arg("system"), bp::arg("lambda_value"), bp::arg("update_constraints")=(bool)(true) ) + , ( bp::arg("system"), bp::arg("lambda_value"), bp::arg("rest2_scale")=(double)(1.0), bp::arg("update_constraints")=(bool)(true) ) , "Set the value of lambda in the passed context. Returns the\n actual value of lambda set.\n" ); } diff --git a/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp b/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp index 2c8453127..95cd479f1 100644 --- a/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp @@ -126,7 +126,7 @@ void register_PerturbableOpenMMMolecule_class(){ PerturbableOpenMMMolecule_exposer_t PerturbableOpenMMMolecule_exposer = PerturbableOpenMMMolecule_exposer_t( "PerturbableOpenMMMolecule", "This class holds all of the information of an OpenMM molecule\nthat can be perturbed using a LambdaSchedule. The data is held\nin easy-to-access arrays, with guarantees that the arrays are\ncompatible and the data is aligned.\n", bp::init< >("Null constructor") ); bp::scope PerturbableOpenMMMolecule_scope( PerturbableOpenMMMolecule_exposer ); PerturbableOpenMMMolecule_exposer.def( bp::init< SireOpenMM::OpenMMMolecule const &, bp::optional< SireBase::PropertyMap const & > >(( bp::arg("mol"), bp::arg("map")=SireBase::PropertyMap() ), "Construct from the passed OpenMMMolecule") ); - PerturbableOpenMMMolecule_exposer.def( bp::init< SireMol::Molecule const &, bp::optional< SireBase::PropertyMap const & > >(( bp::arg("mol"), bp::arg("map")=SireBase::PropertyMap() ), "Construct from a passed molecule and map") ); + PerturbableOpenMMMolecule_exposer.def( bp::init< SireMol::Molecule const &, int, bp::optional< SireBase::PropertyMap const & > >(( bp::arg("mol"), bp::arg("map")=SireBase::PropertyMap() ), "Construct from a passed molecule and map") ); PerturbableOpenMMMolecule_exposer.def( bp::init< SireOpenMM::PerturbableOpenMMMolecule const & >(( bp::arg("other") ), "Copy constructor") ); { //::SireOpenMM::PerturbableOpenMMMolecule::angles diff --git a/wrapper/Convert/SireOpenMM/_sommcontext.py b/wrapper/Convert/SireOpenMM/_sommcontext.py index cceab24ab..15fcba546 100644 --- a/wrapper/Convert/SireOpenMM/_sommcontext.py +++ b/wrapper/Convert/SireOpenMM/_sommcontext.py @@ -69,8 +69,30 @@ def __init__( # place the coordinates and velocities into the context set_openmm_coordinates_and_velocities(self, metadata) + # Check for a REST2 scaling factor. + if map.specified("rest2_scale"): + try: + rest2_scale = map["rest2_scale"].value().as_double() + except: + raise ValueError("'rest2_scale' must be of type 'float'") + if rest2_scale < 1.0: + raise ValueError("'rest2_scale' must be >= 1.0") + self._rest2_scale = rest2_scale + + else: + self._rest2_scale = 1.0 + + # Check for a REST2 selection. + if map.specified("rest2_selection"): + self._has_rest2_selection = True + else: + self._has_rest2_selection = False + self._lambda_value = self._lambda_lever.set_lambda( - self, lambda_value=lambda_value, update_constraints=True + self, + lambda_value=lambda_value, + rest2_scale=self._rest2_scale, + update_constraints=True, ) self._map = map @@ -80,6 +102,8 @@ def __init__( self._lambda_value = 0.0 self._map = None + self._is_non_pert_rest2 = False + def __str__(self): p = self.getPlatform() s = self.getSystem() @@ -224,18 +248,52 @@ def set_lambda_schedule(self, schedule): self._lambda_lever.set_schedule(schedule) - def set_lambda(self, lambda_value: float, update_constraints: bool = True): + def set_lambda( + self, + lambda_value: float, + rest2_scale: float = None, + update_constraints: bool = True, + ): """ Update the parameters in the context to set the lambda value - to 'lamval'. If update_constraints is True then the constraints - will be updated to match the new value of lambda + to 'lambda_value'. The 'rest2_scale' defines the temperature of + the REST2 region relative to the rest of the system. If + 'update_constraints' is True then the constraints will be updated + to match the new value of lambda. """ - if self._lambda_lever is None: - return - self._lambda_value = self._lambda_lever.set_lambda( - self, lambda_value=lambda_value, update_constraints=update_constraints - ) + # If not provided, use the REST2 scaling factor used to initalise + # the context. + if rest2_scale is None: + rest2_scale = self._rest2_scale + else: + if rest2_scale < 1.0: + raise ValueError("'rest2_scale' must be >= 1.0") + + if self._lambda_lever is not None: + self._lambda_value = self._lambda_lever.set_lambda( + self, + lambda_value=lambda_value, + rest2_scale=rest2_scale, + update_constraints=update_constraints, + ) + + # Update any additional parameters in the REST2 region. + if self._is_non_pert_rest2 and rest2_scale != self._rest2_scale: + self._update_rest2(lambda_value, rest2_scale) + self._rest2_scale = rest2_scale + + def get_rest2_scale(self): + """ + Return the temperature scale factor for the REST2 region. + """ + return self._rest2_scale + + def set_rest2_scale(self, rest2_scale): + """ + Set the temperature scale factor for the REST2 region. + """ + self._set_lambda(self._lambda_value, rest2_scale=rest2_scale) def set_temperature(self, temperature, rescale_velocities=True): """ @@ -314,3 +372,148 @@ def to_xml(self, f=None): handle.write(_XmlSerializer.serialize(self.getSystem())) else: f.write(_XmlSerializer.serialize(self.getSystem())) + + def _prepare_rest2(self, system, atoms): + """ + Internal method to prepare the REST2 data structures. + """ + + # Adapted from code in meld: https://github.com/maccallumlab/meld + + # Flag that REST2 modifications are being applied to non-perturbable atoms. + self._is_non_pert_rest2 = True + + import openmm + from ..Mol import AtomIdx, Connectivity + + # Work out the molecules to which the atoms belong. + mols = [] + for atom in atoms: + mol = atom.molecule() + # Perturbable molecules are handled separately. + if mol not in mols and not mol.has_property("is_perturbable"): + mols.append(mol) + + # Store the OpenMM system. + omm_system = self.getSystem() + + # Get the NonBonded force. + nonbonded_force = None + for force in omm_system.getForces(): + if isinstance(force, openmm.NonbondedForce): + nonbonded_force = force + break + if nonbonded_force is None: + raise ValueError("No NonbondedForce found in the OpenMM system.") + self._nonbonded_force = nonbonded_force + + # Get the PeriodicTorsionForce. + periodic_torsion_force = None + for force in omm_system.getForces(): + if isinstance(force, openmm.PeriodicTorsionForce): + periodic_torsion_force = force + break + if periodic_torsion_force is None: + raise ValueError("No PeriodicTorsionForce found in the OpenMM system.") + self._periodic_torsion_force = periodic_torsion_force + + # Initialise the parameter dictionaries. + self._nonbonded_params = {} + self._exception_params = {} + self._torsion_params = {} + + # Store the molecules in the system. + system_mols = system.molecules() + + # Process each of the molecules. + for mol in mols: + # Create the connectivity object for the molecule. + connectivity = Connectivity(mol.info()).edit() + + # Loop over the bonds in the molecule and connect the atoms. + for bond in mol.bonds(): + connectivity.connect(bond.atom0().index(), bond.atom1().index()) + connectivity = connectivity.commit() + + # Find the index of the molecule in the system. + mol_idx = system._system.mol_nums().index(mol.number()) + + # Work out the offset to apply to the atom indices to convert to system indices. + num_atoms = 0 + for i in range(mol_idx): + num_atoms += system_mols[i].num_atoms() + + # Create a list of atom indices. + atom_idxs = [atom.index().value() + num_atoms for atom in atoms] + + # Gather the nonbonded parameters for the atoms in the selection. + for idx in atom_idxs: + self._nonbonded_params[idx] = nonbonded_force.getParticleParameters(idx) + + # Store the exception parameters. + for param_index in range(nonbonded_force.getNumExceptions()): + params = nonbonded_force.getExceptionParameters(param_index) + if params[0] in atom_idxs and params[1] in atom_idxs: + self._exception_params[param_index] = params + + # Gather the torsion parameters for the atoms in the selection. + for param_index in range(periodic_torsion_force.getNumTorsions()): + params = periodic_torsion_force.getTorsionParameters(param_index) + i, j, k, l, _, _, _ = params + + # Don't modify non-REST2 torsions. + if ( + i not in atom_idxs + or j not in atom_idxs + or k not in atom_idxs + or l not in atom_idxs + ): + continue + + # Convert to AtomIdx objects. + idx_i = AtomIdx(i - num_atoms) + idx_l = AtomIdx(l - num_atoms) + + if connectivity.are_dihedraled(idx_i, idx_l): + self._torsion_params[param_index] = params + + def _update_rest2(self, lambda_value, rest2_scale): + """ + Internal method to update the REST2 parameters. + """ + + from math import sqrt + + # This is the temperature scale factor, so we need to invert to get the energy + # scale factor. + scale = 1.0 / rest2_scale + + # Store the REST2 charge scaling factor for non-bonded interactions. + sqrt_scale = sqrt(scale) + + # Update the non-bonded parameters. + for index, params in self._nonbonded_params.items(): + q, sigma, epsilon = params + self._nonbonded_force.setParticleParameters( + index, q * sqrt_scale, sigma, epsilon * scale + ) + + # Update the exception parameters. + for index, params in self._exception_params.items(): + i, j, q, sigma, epsilon = params + self._nonbonded_force.setExceptionParameters( + index, i, j, q * scale, sigma, epsilon * scale + ) + + # Update the parameters in the context. + self._nonbonded_force.updateParametersInContext(self) + + # Update the torsion parameters. + for index, params in self._torsion_params.items(): + i, j, k, l, periodicity, phase, fc = params + self._periodic_torsion_force.setTorsionParameters( + index, i, j, k, l, periodicity, phase, fc * scale + ) + + # Update the parameters in the context. + self._periodic_torsion_force.updateParametersInContext(self) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 67707a398..acc3d9de6 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -26,11 +26,11 @@ * \*********************************************/ -#include "pyqm.h" #include "lambdalever.h" +#include "pyqm.h" -#include "SireBase/propertymap.h" #include "SireBase/arrayproperty.hpp" +#include "SireBase/propertymap.h" #include "SireCAS/values.h" @@ -1126,6 +1126,7 @@ PropertyList LambdaLever::getLeverValues(const QVector &lambda_values, */ double LambdaLever::setLambda(OpenMM::Context &context, double lambda_value, + double rest2_scale, bool update_constraints) const { // go over each forcefield and update the parameters in the forcefield, @@ -1135,6 +1136,13 @@ double LambdaLever::setLambda(OpenMM::Context &context, lambda_value = this->lambda_schedule.clamp(lambda_value); + // This is the temperature scale factor, so we need to invert to get the energy + // scale factor. + rest2_scale = 1.0 / rest2_scale; + + // Store the REST charge scaling factor for non-bonded interactions. + const auto sqrt_rest2_scale = std::sqrt(rest2_scale); + // we need an editable reference to the system to get editable // pointers to the forces... OpenMM::System &system = const_cast(context.getSystem()); @@ -1188,6 +1196,10 @@ double LambdaLever::setLambda(OpenMM::Context &context, // now update the forcefields int start_index = start_idxs.value("clj", -1); + // record the index of the first atom in the perturbable molecule within + // the original Sire system + const auto start_atom_idx = perturbable_mol.getStartAtomIdx(); + if (start_index != -1 and cljff != 0) { const auto morphed_charges = cache.morph( @@ -1353,12 +1365,22 @@ double LambdaLever::setLambda(OpenMM::Context &context, const bool is_from_ghost = perturbable_mol.getFromGhostIdxs().contains(j); const bool is_to_ghost = perturbable_mol.getToGhostIdxs().contains(j); + double scale = 1.0; + double sqrt_scale = 1.0; + + // apply the REST2 scaling. + if (perturbable_mol.isRest2(j)) + { + scale = rest2_scale; + sqrt_scale = sqrt_rest2_scale; + } + // reduced charge - custom_params[0] = morphed_ghost_charges[j]; + custom_params[0] = sqrt_scale * morphed_ghost_charges[j]; // half_sigma custom_params[1] = 0.5 * morphed_ghost_sigmas[j]; // two_sqrt_epsilon - custom_params[2] = 2.0 * std::sqrt(morphed_ghost_epsilons[j]); + custom_params[2] = 2.0 * sqrt_scale * std::sqrt(morphed_ghost_epsilons[j]); // alpha custom_params[3] = morphed_ghost_alphas[j]; // kappa @@ -1379,11 +1401,11 @@ double LambdaLever::setLambda(OpenMM::Context &context, ghost_ghostff->setParticleParameters(start_index + j, custom_params); // reduced charge - custom_params[0] = morphed_nonghost_charges[j]; + custom_params[0] = sqrt_scale * morphed_nonghost_charges[j]; // half_sigma custom_params[1] = 0.5 * morphed_nonghost_sigmas[j]; // two_sqrt_epsilon - custom_params[2] = 2.0 * std::sqrt(morphed_nonghost_epsilons[j]); + custom_params[2] = 2.0 * sqrt_scale * std::sqrt(morphed_nonghost_epsilons[j]); // alpha custom_params[3] = morphed_nonghost_alphas[j]; // kappa @@ -1406,11 +1428,13 @@ double LambdaLever::setLambda(OpenMM::Context &context, if (is_from_ghost or is_to_ghost) { // don't set the LJ parameters in the cljff - cljff->setParticleParameters(start_index + j, morphed_charges[j], 0.0, 0.0); + cljff->setParticleParameters(start_index + j, sqrt_scale * morphed_charges[j], 0.0, 0.0); } else { - cljff->setParticleParameters(start_index + j, morphed_charges[j], morphed_sigmas[j], morphed_epsilons[j]); + cljff->setParticleParameters( + start_index + j, sqrt_scale* morphed_charges[j], + morphed_sigmas[j], scale * morphed_epsilons[j]); } } } @@ -1418,7 +1442,18 @@ double LambdaLever::setLambda(OpenMM::Context &context, { for (int j = 0; j < nparams; ++j) { - cljff->setParticleParameters(start_index + j, morphed_charges[j], morphed_sigmas[j], morphed_epsilons[j]); + double scale = 1.0; + double sqrt_scale = 1.0; + + // apply the REST2 scaling. + if (perturbable_mol.isRest2(j)) + { + scale = rest2_scale; + sqrt_scale = sqrt_rest2_scale; + } + + cljff->setParticleParameters(start_index + j, sqrt_scale * morphed_charges[j], + morphed_sigmas[j], scale * morphed_epsilons[j]); } } @@ -1446,6 +1481,14 @@ double LambdaLever::setLambda(OpenMM::Context &context, morphed_charges, morphed_sigmas, morphed_epsilons, morphed_alphas, morphed_kappas); + double scale = 1.0; + + // apply the REST2 scaling. + if (perturbable_mol.isRest2(atom0) and perturbable_mol.isRest2(atom1)) + { + scale = rest2_scale; + } + // don't set LJ terms for ghost atoms if (atom0_is_ghost or atom1_is_ghost) { @@ -1462,7 +1505,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, cljff->setExceptionParameters( boost::get<0>(idxs[j]), boost::get<0>(p), boost::get<1>(p), - boost::get<2>(p), 1e-9, 1e-9); + boost::get<2>(p) * scale, 1e-9, 1e-9); // exclude 14s for to/from ghost interactions if (not to_from_ghost and ghost_14ff != 0) @@ -1488,10 +1531,13 @@ double LambdaLever::setLambda(OpenMM::Context &context, morphed_ghost14_kappas); // parameters are q, sigma, four_epsilon and alpha - std::vector params14 = - {boost::get<2>(p), boost::get<3>(p), - 4.0 * boost::get<4>(p), boost::get<5>(p), - boost::get<6>(p)}; + std::vector params14 = { + boost::get<2>(p) * scale, + boost::get<3>(p), + 4.0 * boost::get<4>(p) * scale, + boost::get<5>(p), + boost::get<6>(p) + }; if (start_change_14 == -1) { @@ -1518,8 +1564,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, cljff->setExceptionParameters( boost::get<0>(idxs[j]), boost::get<0>(p), boost::get<1>(p), - boost::get<2>(p), boost::get<3>(p), - boost::get<4>(p)); + boost::get<2>(p) * scale, boost::get<3>(p), + boost::get<4>(p) * scale); } } } @@ -1679,6 +1725,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, const int nparams = morphed_torsion_k.count(); + const auto is_improper = perturbable_mol.getIsImproper(); + if (start_change_torsion == -1) { start_change_torsion = start_index; @@ -1707,12 +1755,30 @@ double LambdaLever::setLambda(OpenMM::Context &context, particle3, particle4, periodicity, phase, k); + // get the indices of the particles in the Sire molecule + const auto atom1 = particle1 - start_atom_idx; + const auto atom2 = particle2 - start_atom_idx; + const auto atom3 = particle3 - start_atom_idx; + const auto atom4 = particle4 - start_atom_idx; + + // Only apply the REST2 scaling factor if this is a proper dihedral + // and all atoms are in the REST2 region. + double scale = 1.0; + if (not is_improper[j] and + perturbable_mol.isRest2(atom1) and + perturbable_mol.isRest2(atom2) and + perturbable_mol.isRest2(atom3) and + perturbable_mol.isRest2(atom4)) + { + scale = rest2_scale; + } + dihff->setTorsionParameters(index, particle1, particle2, particle3, particle4, periodicity, morphed_torsion_phase[j], - morphed_torsion_k[j]); + morphed_torsion_k[j] * scale); } } } @@ -1857,6 +1923,100 @@ void _update_restraint_in_context(OpenMM::CustomCompoundBondForce *ff, double rh ff->updateParametersInContext(context); } +/** Update the parameters for a CustomAngleForce for scale factor 'rho' + * in the passed context */ +void _update_restraint_in_context(OpenMM::CustomAngleForce *ff, double rho, + OpenMM::Context &context) +{ + if (ff == 0) + throw SireError::invalid_cast(QObject::tr( + "Unable to cast the restraint force to an OpenMM::CustomAngleForce, " + "despite it reporting that is was an object of this type..."), + CODELOC); + + const int nangles = ff->getNumAngles(); + + if (nangles == 0) + // nothing to update + return; + + const int nparams = ff->getNumPerAngleParameters(); + + if (nparams == 0) + throw SireError::incompatible_error(QObject::tr( + "Unable to set 'rho' for this restraint as it has no custom parameters!"), + CODELOC); + + // we set the first parameter - we can see what the current value + // is from the first restraint. This is because rho should be the + // first parameter and have the same value for all restraints + std::vector custom_params; + custom_params.resize(nparams); + int atom0, atom1, atom2; + + ff->getAngleParameters(0, atom0, atom1, atom2, custom_params); + + if (custom_params[0] == rho) + // nothing to do - it is already equal to this value + return; + + for (int i = 0; i < nangles; ++i) + { + ff->getAngleParameters(i, atom0, atom1, atom2, custom_params); + custom_params[0] = rho; + ff->setAngleParameters(i, atom0, atom1, atom2, custom_params); + } + + ff->updateParametersInContext(context); +} + +/** Update the parameters for a CustomTorsionForce for scale factor 'rho' + * in the passed context */ +void _update_restraint_in_context(OpenMM::CustomTorsionForce *ff, double rho, + OpenMM::Context &context) +{ + if (ff == 0) + throw SireError::invalid_cast(QObject::tr( + "Unable to cast the restraint force to an OpenMM::CustomTorsionForce, " + "despite it reporting that is was an object of this type..."), + CODELOC); + + const int ntorsions = ff->getNumTorsions(); + + if (ntorsions == 0) + // nothing to update + return; + + const int nparams = ff->getNumPerTorsionParameters(); + + if (nparams == 0) + throw SireError::incompatible_error(QObject::tr( + "Unable to set 'rho' for this restraint as it has no custom parameters!"), + CODELOC); + + // we set the first parameter - we can see what the current value + // is from the first restraint. This is because rho should be the + // first parameter and have the same value for all restraints + std::vector custom_params; + custom_params.resize(nparams); + int atom0, atom1, atom2, atom3; + + ff->getTorsionParameters(0, atom0, atom1, atom2, atom3, custom_params); + + if (custom_params[0] == rho) + // nothing to do - it is already equal to this value + return; + + for (int i = 0; i < ntorsions; ++i) + { + ff->getTorsionParameters(i, atom0, atom1, atom2, atom3, custom_params); + custom_params[0] = rho; + ff->setTorsionParameters(i, atom0, atom1, atom2, atom3, custom_params); + } + + ff->updateParametersInContext(context); +} + /** Update the parameters for a CustomBondForce for scale factor 'rho' * in the passed context */ void _update_restraint_in_context(OpenMM::CustomBondForce *ff, double rho, @@ -1918,6 +2078,24 @@ void LambdaLever::updateRestraintInContext(OpenMM::Force &ff, double rho, dynamic_cast(&ff), rho, context); } + else if (ff_type == "AngleRestraintForce") + { + _update_restraint_in_context( + dynamic_cast(&ff), + rho, context); + } + else if (ff_type == "TorsionRestraintForce") + { + _update_restraint_in_context( + dynamic_cast(&ff), + rho, context); + } + else if (ff_type == "CustomCompoundBondForce") + { + _update_restraint_in_context( + dynamic_cast(&ff), + rho, context); + } else if (ff_type == "BoreschRestraintForce") { _update_restraint_in_context( @@ -2007,7 +2185,6 @@ int LambdaLever::addPerturbableMolecule(const OpenMMMolecule &molecule, return this->perturbable_mols.count() - 1; } - /** Set the exception indices for the perturbable molecule at * index 'mol_idx' */ diff --git a/wrapper/Convert/SireOpenMM/lambdalever.h b/wrapper/Convert/SireOpenMM/lambdalever.h index 513025800..c163ceca6 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.h +++ b/wrapper/Convert/SireOpenMM/lambdalever.h @@ -117,7 +117,7 @@ namespace SireOpenMM const PerturbableOpenMMMolecule &mol) const; double setLambda(OpenMM::Context &system, double lambda_value, - bool update_constraints = true) const; + double rest2_scale = 1.0, bool update_constraints = true) const; void setForceIndex(const QString &force, int index); diff --git a/wrapper/Convert/SireOpenMM/openmmminimise.cpp b/wrapper/Convert/SireOpenMM/openmmminimise.cpp index 500401b4b..d42a038ed 100644 --- a/wrapper/Convert/SireOpenMM/openmmminimise.cpp +++ b/wrapper/Convert/SireOpenMM/openmmminimise.cpp @@ -1090,6 +1090,8 @@ namespace SireOpenMM } } + PyEval_RestoreThread(_save); + if (is_success) { data.addLog("Minimisation successful!"); @@ -1105,8 +1107,6 @@ namespace SireOpenMM CODELOC); } - PyEval_RestoreThread(_save); - return data.getLog().join("\n"); } diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 40979610d..4552afb66 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -51,11 +51,15 @@ OpenMMMolecule::OpenMMMolecule() } OpenMMMolecule::OpenMMMolecule(const Molecule &mol, + int start_atom_idx, const PropertyMap &map) { molinfo = mol.info(); number = mol.number(); + // store the starting index of the atoms in the OpenMM system + this->start_atom_idx = start_atom_idx; + if (molinfo.nAtoms() == 0) { // nothing to extract @@ -480,6 +484,20 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, atoms = mol.atoms(); const int nats = atoms.count(); + // Create the REST2 region mask. + if (mol.hasProperty("is_rest2")) + { + const auto &is_rest2_prop = moldata.property("is_rest2").asA(); + for (int i = 0; i < nats; ++i) + { + is_rest2.append(is_rest2_prop.at(i) == 1); + } + } + else + { + is_rest2 = QVector(nats, true); + } + if (nats <= 0) { return; @@ -1057,6 +1075,7 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, const double dihedral_k_to_openmm = (SireUnits::kcal_per_mol).to(SireUnits::kJ_per_mol); dih_params.clear(); + is_improper.clear(); for (auto it = params.dihedrals().constBegin(); it != params.dihedrals().constEnd(); @@ -1087,6 +1106,7 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, dih_params.append(boost::make_tuple(atom0, atom1, atom2, atom3, periodicity, phase, v)); + is_improper.append(false); } else if (periodicity == 0 and v == 0) { @@ -1095,6 +1115,7 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, dih_params.append(boost::make_tuple(atom0, atom1, atom2, atom3, 1, phase, v)); + is_improper.append(false); } else { @@ -1127,6 +1148,7 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, dih_params.append(boost::make_tuple(atom0, atom1, atom2, atom3, periodicity, phase, v)); + is_improper.append(true); } else if (periodicity == 0 and v == 0) { @@ -1135,6 +1157,7 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, dih_params.append(boost::make_tuple(atom0, atom1, atom2, atom3, 1, phase, v)); + is_improper.append(true); } else { @@ -1369,7 +1392,9 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) perturbed->ang_params = ang_params_1; QVector> dih_params_1; + QVector is_improper_1; dih_params_1.reserve(dih_params.count()); + is_improper_1.reserve(dih_params.count()); found_index_0 = QVector(dih_params.count(), false); found_index_1 = QVector(perturbed->dih_params.count(), false); @@ -1377,6 +1402,7 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) for (int i = 0; i < dih_params.count(); ++i) { const auto &dih0 = dih_params.at(i); + const bool is_improper0 = is_improper.at(i); int atom0 = boost::get<0>(dih0); int atom1 = boost::get<1>(dih0); @@ -1390,6 +1416,7 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) if (not found_index_1[j]) { const auto &dih1 = perturbed->dih_params.at(j); + const bool is_improper1 = perturbed->is_improper.at(j); // we need to match all of the atoms AND the periodicity if (boost::get<0>(dih1) == atom0 and boost::get<1>(dih1) == atom1 and @@ -1398,6 +1425,7 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) { // we have found the matching torsion! dih_params_1.append(dih1); + is_improper_1.append(is_improper1); found_index_0[i] = true; found_index_1[j] = true; found = true; @@ -1410,6 +1438,7 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) { // add a null dihedral with the same periodicity and phase, but null k dih_params_1.append(boost::tuple(atom0, atom1, atom2, atom3, boost::get<4>(dih0), boost::get<5>(dih0), 0.0)); + is_improper_1.append(is_improper0); found_index_0[i] = true; } } @@ -1420,6 +1449,7 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) { // need to add a dihedral missing in the reference state const auto &dih1 = perturbed->dih_params.at(j); + const auto is_improper1 = perturbed->is_improper.at(j); int atom0 = boost::get<0>(dih1); int atom1 = boost::get<1>(dih1); @@ -1428,7 +1458,9 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) // add a null dihedral with the same periodicity and phase, but null k dih_params.append(boost::tuple(atom0, atom1, atom2, atom3, boost::get<4>(dih1), boost::get<5>(dih1), 0.0)); + is_improper.append(is_improper1); dih_params_1.append(dih1); + is_improper_1.append(is_improper1); found_index_1[j] = true; } } @@ -1442,6 +1474,8 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) } perturbed->dih_params = dih_params_1; + perturbed->is_improper = is_improper_1; + is_improper = is_improper_1; // now align all of the exceptions - this should allow the bonding // to change during the perturbation @@ -2089,10 +2123,11 @@ PerturbableOpenMMMolecule::PerturbableOpenMMMolecule() /** Construct from a passed molecule and map */ PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const Molecule &mol, + int start_atom_idx, const PropertyMap &map) : ConcreteProperty() { - this->operator=(PerturbableOpenMMMolecule(OpenMMMolecule(mol, map), map)); + this->operator=(PerturbableOpenMMMolecule(OpenMMMolecule(mol, start_atom_idx, map), map)); } /** Return whether or not this is null */ @@ -2116,6 +2151,9 @@ PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const OpenMMMolecule &mol, auto molecule = mol.atoms.molecule(); + // store the index of the first atom in the OpenMM system + this->start_atom_idx = mol.start_atom_idx; + for (int i = 0; i < mol.atoms.count(); ++i) { perturbed_atoms.append(mol.atoms(i)); @@ -2196,6 +2234,9 @@ PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const OpenMMMolecule &mol, perturbable_constraints = mol.perturbable_constraints; + is_improper = mol.is_improper; + is_rest2 = mol.is_rest2; + bool fix_perturbable_zero_sigmas = false; if (map.specified("fix_perturbable_zero_sigmas")) @@ -2261,7 +2302,10 @@ PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const PerturbableOpenMMMole to_ghost_idxs(other.to_ghost_idxs), from_ghost_idxs(other.from_ghost_idxs), exception_atoms(other.exception_atoms), exception_idxs(other.exception_idxs), perturbable_constraints(other.perturbable_constraints), - constraint_idxs(other.constraint_idxs) + constraint_idxs(other.constraint_idxs), + start_atom_idx(other.start_atom_idx), + is_improper(other.is_improper), + is_rest2(other.is_rest2) { } @@ -2705,3 +2749,21 @@ QList PerturbableOpenMMMolecule::torsions() const { return perturbed_dihs; } + +/** Return the improper flag vector for the perturbed dihedrals */ +QVector PerturbableOpenMMMolecule::getIsImproper() const +{ + return is_improper; +} + +/** Return the index of the first atom in the Sire system */ +int PerturbableOpenMMMolecule::getStartAtomIdx() const +{ + return start_atom_idx; +} + +/** Whether the atom is in the REST2 region. */ +bool PerturbableOpenMMMolecule::isRest2(int atom) const +{ + return is_rest2[atom]; +} diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.h b/wrapper/Convert/SireOpenMM/openmmmolecule.h index f21c44746..a903048cb 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.h +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.h @@ -44,6 +44,7 @@ namespace SireOpenMM OpenMMMolecule(); OpenMMMolecule(const SireMol::Molecule &mol, + int start_atom_idx, const SireBase::PropertyMap &map); ~OpenMMMolecule(); @@ -183,6 +184,15 @@ namespace SireOpenMM * perturbable atoms */ qint32 perturbable_constraint_type; + /** Whether each dihedral is an improper */ + QVector is_improper; + + /** Whether each atom is within the REST2 region */ + QVector is_rest2; + + /** The starting index of the first OpenMM atom in the original Sire system. */ + int start_atom_idx; + private: void constructFromAmber(const SireMol::Molecule &mol, const SireMM::AmberParams ¶ms, @@ -211,6 +221,7 @@ namespace SireOpenMM const SireBase::PropertyMap &map = SireBase::PropertyMap()); PerturbableOpenMMMolecule(const SireMol::Molecule &mol, + int start_atom_idx, const SireBase::PropertyMap &map = SireBase::PropertyMap()); PerturbableOpenMMMolecule(const PerturbableOpenMMMolecule &other); @@ -271,6 +282,9 @@ namespace SireOpenMM QSet getToGhostIdxs() const; QSet getFromGhostIdxs() const; + int getStartAtomIdx() const; + bool isRest2(int atom) const; + bool isGhostAtom(int atom) const; QVector> getExceptionAtoms() const; @@ -293,6 +307,8 @@ namespace SireOpenMM QList angles() const; QList torsions() const; + QVector getIsImproper() const; + private: /** The atoms that are perturbed, in the order they appear * in the arrays below @@ -358,8 +374,16 @@ namespace SireOpenMM * to the number of perturbable constraints in the molecule */ QVector constraint_idxs; - }; + /** Whether each dihedral is an improper */ + QVector is_improper; + + /** Whether each atom is within the REST2 region */ + QVector is_rest2; + + /** The starting index of the first OpenMM atom in the original Sire system. */ + int start_atom_idx; + }; } Q_DECLARE_METATYPE(SireOpenMM::PerturbableOpenMMMolecule) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index e3c46a03e..7016e7c7a 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -8,37 +8,39 @@ #include "SireSystem/forcefieldinfo.h" -#include "SireMol/core.h" -#include "SireMol/moleditor.h" -#include "SireMol/atomelements.h" #include "SireMol/atomcharges.h" #include "SireMol/atomcoords.h" +#include "SireMol/atomelements.h" #include "SireMol/atommasses.h" #include "SireMol/atomproperty.hpp" -#include "SireMol/connectivity.h" +#include "SireMol/atomvelocities.h" #include "SireMol/bondid.h" #include "SireMol/bondorder.h" -#include "SireMol/atomvelocities.h" +#include "SireMol/connectivity.h" +#include "SireMol/core.h" +#include "SireMol/moleditor.h" -#include "SireMM/atomljs.h" -#include "SireMM/selectorbond.h" #include "SireMM/amberparams.h" +#include "SireMM/anglerestraints.h" +#include "SireMM/atomljs.h" #include "SireMM/bondrestraints.h" -#include "SireMM/positionalrestraints.h" #include "SireMM/boreschrestraints.h" +#include "SireMM/dihedralrestraints.h" +#include "SireMM/positionalrestraints.h" +#include "SireMM/selectorbond.h" #include "SireVol/periodicbox.h" #include "SireVol/triclinicbox.h" #include "SireCAS/lambdaschedule.h" -#include "SireMaths/vector.h" #include "SireMaths/maths.h" +#include "SireMaths/vector.h" +#include "SireBase/generalunitproperty.h" +#include "SireBase/lengthproperty.h" #include "SireBase/parallel.h" #include "SireBase/propertylist.h" -#include "SireBase/lengthproperty.h" -#include "SireBase/generalunitproperty.h" #include "SireUnits/units.h" @@ -398,6 +400,132 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, } } +/** Add all of the angle restraints from 'restraints' to the passed + * system, which is acted on by the passed LambdaLever. The number + * of real (non-anchor) atoms in the OpenMM::System is 'natoms' + */ + +void _add_angle_restraints(const SireMM::AngleRestraints &restraints, + OpenMM::System &system, LambdaLever &lambda_lever, + int natoms) +{ + if (restraints.isEmpty()) + return; + + // energy expression of the angle restraint, which acts over three atoms + // + // theta = angle(P1, P2, P3) + // + // The energies are + // + // e_restraint = rho * (e_angle) + // e_angle = ktheta(theta - theta0)^2 + + const auto energy_expression = QString( + "rho*k*(theta-theta0)^2;") + .toStdString(); + + auto *restraintff = new OpenMM::CustomAngleForce(energy_expression); + + restraintff->setName("AngleRestraintForce"); + restraintff->addPerAngleParameter("rho"); + restraintff->addPerAngleParameter("k"); + restraintff->addPerAngleParameter("theta0"); + + restraintff->setUsesPeriodicBoundaryConditions(true); + + lambda_lever.addRestraintIndex(restraints.name(), + system.addForce(restraintff)); + + const double internal_to_ktheta = (1 * SireUnits::kcal_per_mol / (SireUnits::radian2)).to(SireUnits::kJ_per_mol / SireUnits::radian2); + + const auto atom_restraints = restraints.restraints(); + + for (const auto &restraint : atom_restraints) + { + std::vector particles; + particles.resize(3); + + for (int i = 0; i < 3; ++i) + { + particles[i] = restraint.atoms()[i]; + } + + std::vector parameters; + parameters.resize(3); + + parameters[0] = 1.0; // rho + parameters[1] = restraint.ktheta().value() * internal_to_ktheta; // k + parameters[2] = restraint.theta0().value(); // theta0 (already in radians) + + // restraintff->addTorsion(particles, parameters); + restraintff->addAngle(particles[0], particles[1], particles[2], parameters); + } +} + +void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, + OpenMM::System &system, LambdaLever &lambda_lever, + int natoms) +{ + if (restraints.isEmpty()) + return; + + // energy expression of the dihedral restraint, which acts over four atoms + // + // phi = dihedral(P1, P2, P3, P4) + // + // The energies are + // + // e_restraint = rho * (e_torsion) + // e_torsion = k_phi(dphi)^2 where + // dphi = abs(phi - phi0) + + const auto energy_expression = QString( + "rho*k*min(dtheta, two_pi-dtheta)^2;" + "dtheta = abs(theta-theta0);" + "two_pi=6.283185307179586;") + .toStdString(); + + auto *restraintff = new OpenMM::CustomTorsionForce(energy_expression); + + // it seems that OpenMM wants to call the torsion angle theta rather than phi + // we need to rename our parameters accordingly + restraintff->setName("TorsionRestraintForce"); + restraintff->addPerTorsionParameter("rho"); + restraintff->addPerTorsionParameter("k"); + restraintff->addPerTorsionParameter("theta0"); + + restraintff->setUsesPeriodicBoundaryConditions(true); + + lambda_lever.addRestraintIndex(restraints.name(), + system.addForce(restraintff)); + + const double internal_to_ktheta = (1 * SireUnits::kcal_per_mol / (SireUnits::radian2)).to(SireUnits::kJ_per_mol / SireUnits::radian2); + + const auto atom_restraints = restraints.restraints(); + + for (const auto &restraint : atom_restraints) + { + std::vector particles; + particles.resize(4); + + for (int i = 0; i < 4; ++i) + { + particles[i] = restraint.atoms()[i]; + } + + std::vector parameters; + parameters.resize(3); + + parameters[0] = 1.0; // rho + parameters[1] = restraint.kphi().value() * internal_to_ktheta; // k + parameters[2] = restraint.phi0().value(); // phi0 (already in radians) --> theta0 + + // restraintff->addTorsion(particles, parameters); + restraintff->addTorsion(particles[0], particles[1], particles[2], particles[3], parameters); + } +} + /** Set the coulomb and LJ cutoff in the passed NonbondedForce, * based on the information in the passed ForceFieldInfo. * This sets the cutoff type (e.g. PME) and the actual @@ -665,20 +793,28 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, QVector openmm_mols(nmols); auto openmm_mols_data = openmm_mols.data(); + // Create a vector containing the start index for the atoms in each molecule. + QVector start_atom_index(nmols); + start_atom_index[0] = 0; + for (int i = 1; i < nmols; ++i) + { + start_atom_index[i] = start_atom_index[i-1] + mols[i-1].nAtoms(); + } + if (SireBase::should_run_in_parallel(nmols, map)) { tbb::parallel_for(tbb::blocked_range(0, mols.count()), [&](const tbb::blocked_range &r) { for (int i=r.begin(); i()) + if (prop.read().isA()) + { + _add_dihedral_restraints(prop.read().asA(), + system, lambda_lever, start_index); + } + else if (prop.read().isA()) + { + _add_angle_restraints(prop.read().asA(), + system, lambda_lever, start_index); + } + else if (prop.read().isA()) { _add_positional_restraints(prop.read().asA(), system, lambda_lever, anchor_coords, start_index); diff --git a/wrapper/MM/AngleRestraint.pypp.cpp b/wrapper/MM/AngleRestraint.pypp.cpp index 052916bb4..b90320621 100644 --- a/wrapper/MM/AngleRestraint.pypp.cpp +++ b/wrapper/MM/AngleRestraint.pypp.cpp @@ -3,36 +3,25 @@ // (C) Christopher Woods, GPL >= 3 License #include "boost/python.hpp" -#include "Helpers/clone_const_reference.hpp" #include "AngleRestraint.pypp.hpp" namespace bp = boost::python; -#include "SireCAS/conditional.h" - #include "SireCAS/errors.h" -#include "SireCAS/power.h" - -#include "SireCAS/symbols.h" - -#include "SireCAS/values.h" - -#include "SireFF/forcetable.h" - #include "SireID/index.h" #include "SireStream/datastream.h" #include "SireStream/shareddatastream.h" -#include "SireUnits/angle.h" - #include "SireUnits/units.h" -#include "anglerestraint.h" +#include "anglerestraints.h" + +#include -#include "anglerestraint.h" +#include "anglerestraints.h" SireMM::AngleRestraint __copy__(const SireMM::AngleRestraint &other){ return SireMM::AngleRestraint(other); } @@ -47,162 +36,50 @@ SireMM::AngleRestraint __copy__(const SireMM::AngleRestraint &other){ return Sir void register_AngleRestraint_class(){ { //::SireMM::AngleRestraint - typedef bp::class_< SireMM::AngleRestraint, bp::bases< SireMM::Restraint3D, SireMM::Restraint, SireBase::Property > > AngleRestraint_exposer_t; - AngleRestraint_exposer_t AngleRestraint_exposer = AngleRestraint_exposer_t( "AngleRestraint", "This is a restraint that operates on the angle between\nthree SireMM::Point objects (e.g. three atoms in a molecule)\n\nAuthor: Christopher Woods\n", bp::init< >("Constructor") ); + typedef bp::class_< SireMM::AngleRestraint, bp::bases< SireBase::Property > > AngleRestraint_exposer_t; + AngleRestraint_exposer_t AngleRestraint_exposer = AngleRestraint_exposer_t( "AngleRestraint", "This class represents a single angle restraint between any three\natoms in a system\nAuthor: Christopher Woods\n", bp::init< >("Null constructor") ); bp::scope AngleRestraint_scope( AngleRestraint_exposer ); - AngleRestraint_exposer.def( bp::init< SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireCAS::Expression const & >(( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("restraint") ), "Construct a restraint that acts on the angle within the\nthree points point0, point1 and point2 (theta == a(012)),\nrestraining the angle within these points using the expression\nrestraint") ); - AngleRestraint_exposer.def( bp::init< SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireCAS::Expression const &, SireCAS::Values const & >(( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("restraint"), bp::arg("values") ), "Construct a restraint that acts on the angle within the\nthree points point0, point1 and point2 (theta == a(012)),\nrestraining the angle within these points using the expression\nrestraint") ); + AngleRestraint_exposer.def( bp::init< QList< long long > const &, SireUnits::Dimension::Angle const &, SireUnits::Dimension::HarmonicAngleConstant const & >(( bp::arg("atoms"), bp::arg("theta0"), bp::arg("ktheta") ), "Construct a restraint that acts on the angle within the\nthree atoms atom0, atom1 and atom2 (theta == a(012)),\nrestraining the angle within these atoms") ); AngleRestraint_exposer.def( bp::init< SireMM::AngleRestraint const & >(( bp::arg("other") ), "Copy constructor") ); - { //::SireMM::AngleRestraint::builtinSymbols - - typedef ::SireCAS::Symbols ( ::SireMM::AngleRestraint::*builtinSymbols_function_type)( ) const; - builtinSymbols_function_type builtinSymbols_function_value( &::SireMM::AngleRestraint::builtinSymbols ); - - AngleRestraint_exposer.def( - "builtinSymbols" - , builtinSymbols_function_value - , bp::release_gil_policy() - , "Return the built-in symbols of this restraint" ); - - } - { //::SireMM::AngleRestraint::builtinValues - - typedef ::SireCAS::Values ( ::SireMM::AngleRestraint::*builtinValues_function_type)( ) const; - builtinValues_function_type builtinValues_function_value( &::SireMM::AngleRestraint::builtinValues ); - - AngleRestraint_exposer.def( - "builtinValues" - , builtinValues_function_value - , bp::release_gil_policy() - , "Return the values of the built-in symbols of this restraint" ); - - } - { //::SireMM::AngleRestraint::contains - - typedef bool ( ::SireMM::AngleRestraint::*contains_function_type)( ::SireMol::MolNum ) const; - contains_function_type contains_function_value( &::SireMM::AngleRestraint::contains ); - - AngleRestraint_exposer.def( - "contains" - , contains_function_value - , ( bp::arg("molnum") ) - , bp::release_gil_policy() - , "Return whether or not this restraint affects the molecule\nwith number molnum" ); - - } - { //::SireMM::AngleRestraint::contains - - typedef bool ( ::SireMM::AngleRestraint::*contains_function_type)( ::SireMol::MolID const & ) const; - contains_function_type contains_function_value( &::SireMM::AngleRestraint::contains ); - - AngleRestraint_exposer.def( - "contains" - , contains_function_value - , ( bp::arg("molid") ) - , bp::release_gil_policy() - , "Return whether or not this restraint affects the molecule\nwith ID molid" ); - - } - { //::SireMM::AngleRestraint::differentialRestraintFunction - - typedef ::SireCAS::Expression const & ( ::SireMM::AngleRestraint::*differentialRestraintFunction_function_type)( ) const; - differentialRestraintFunction_function_type differentialRestraintFunction_function_value( &::SireMM::AngleRestraint::differentialRestraintFunction ); - - AngleRestraint_exposer.def( - "differentialRestraintFunction" - , differentialRestraintFunction_function_value - , bp::return_value_policy< bp::copy_const_reference >() - , "Return the function used to calculate the restraint force" ); - - } - { //::SireMM::AngleRestraint::differentiate - - typedef ::SireMM::RestraintPtr ( ::SireMM::AngleRestraint::*differentiate_function_type)( ::SireCAS::Symbol const & ) const; - differentiate_function_type differentiate_function_value( &::SireMM::AngleRestraint::differentiate ); - - AngleRestraint_exposer.def( - "differentiate" - , differentiate_function_value - , ( bp::arg("symbol") ) - , bp::release_gil_policy() - , "Return the differential of this restraint with respect to\nthe symbol symbol\nThrow: SireCAS::unavailable_differential\n" ); - - } - { //::SireMM::AngleRestraint::force - - typedef void ( ::SireMM::AngleRestraint::*force_function_type)( ::SireFF::MolForceTable &,double ) const; - force_function_type force_function_value( &::SireMM::AngleRestraint::force ); - - AngleRestraint_exposer.def( - "force" - , force_function_value - , ( bp::arg("forcetable"), bp::arg("scale_force")=1 ) - , "Calculate the force acting on the molecule in the forcetable forcetable\ncaused by this restraint, and add it on to the forcetable scaled by\nscale_force" ); - - } - { //::SireMM::AngleRestraint::force - - typedef void ( ::SireMM::AngleRestraint::*force_function_type)( ::SireFF::ForceTable &,double ) const; - force_function_type force_function_value( &::SireMM::AngleRestraint::force ); - - AngleRestraint_exposer.def( - "force" - , force_function_value - , ( bp::arg("forcetable"), bp::arg("scale_force")=1 ) - , "Calculate the force acting on the molecules in the forcetable forcetable\ncaused by this restraint, and add it on to the forcetable scaled by\nscale_force" ); + { //::SireMM::AngleRestraint::atoms - } - { //::SireMM::AngleRestraint::halfHarmonic - - typedef ::SireMM::AngleRestraint ( *halfHarmonic_function_type )( ::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireUnits::Dimension::Angle const &,::SireMM::HarmonicAngleForceConstant const & ); - halfHarmonic_function_type halfHarmonic_function_value( &::SireMM::AngleRestraint::halfHarmonic ); + typedef ::QVector< long long > ( ::SireMM::AngleRestraint::*atoms_function_type)( ) const; + atoms_function_type atoms_function_value( &::SireMM::AngleRestraint::atoms ); AngleRestraint_exposer.def( - "halfHarmonic" - , halfHarmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("angle"), bp::arg("force_constant") ) + "atoms" + , atoms_function_value , bp::release_gil_policy() - , "Return a distance restraint that applied a half-harmonic potential\nbetween the points point0 and point1 above a distance distance\nusing a force constant force_constant" ); + , "Return the atoms involved in the restraint" ); } - { //::SireMM::AngleRestraint::harmonic + { //::SireMM::AngleRestraint::isNull - typedef ::SireMM::AngleRestraint ( *harmonic_function_type )( ::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireMM::HarmonicAngleForceConstant const & ); - harmonic_function_type harmonic_function_value( &::SireMM::AngleRestraint::harmonic ); + typedef bool ( ::SireMM::AngleRestraint::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::AngleRestraint::isNull ); AngleRestraint_exposer.def( - "harmonic" - , harmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("force_constant") ) + "isNull" + , isNull_function_value , bp::release_gil_policy() - , "Return a distance restraint that applies a harmonic potential between\nthe points point0 and point1 using a force constant force_constant" ); - - } - { //::SireMM::AngleRestraint::molecules - - typedef ::SireMol::Molecules ( ::SireMM::AngleRestraint::*molecules_function_type)( ) const; - molecules_function_type molecules_function_value( &::SireMM::AngleRestraint::molecules ); - - AngleRestraint_exposer.def( - "molecules" - , molecules_function_value - , bp::release_gil_policy() - , "Return the molecules used in this restraint" ); + , "" ); } - { //::SireMM::AngleRestraint::nPoints + { //::SireMM::AngleRestraint::ktheta - typedef int ( ::SireMM::AngleRestraint::*nPoints_function_type)( ) const; - nPoints_function_type nPoints_function_value( &::SireMM::AngleRestraint::nPoints ); + typedef ::SireUnits::Dimension::HarmonicAngleConstant ( ::SireMM::AngleRestraint::*ktheta_function_type)( ) const; + ktheta_function_type ktheta_function_value( &::SireMM::AngleRestraint::ktheta ); AngleRestraint_exposer.def( - "nPoints" - , nPoints_function_value + "ktheta" + , ktheta_function_value , bp::release_gil_policy() - , "This restraint involves three points" ); + , "Return the force constant for the restraint" ); } AngleRestraint_exposer.def( bp::self != bp::self ); + AngleRestraint_exposer.def( bp::self + bp::self ); + AngleRestraint_exposer.def( bp::self + bp::other< SireMM::AngleRestraints >() ); { //::SireMM::AngleRestraint::operator= typedef ::SireMM::AngleRestraint & ( ::SireMM::AngleRestraint::*assign_function_type)( ::SireMM::AngleRestraint const & ) ; @@ -217,78 +94,28 @@ void register_AngleRestraint_class(){ } AngleRestraint_exposer.def( bp::self == bp::self ); - { //::SireMM::AngleRestraint::point - - typedef ::SireFF::Point const & ( ::SireMM::AngleRestraint::*point_function_type)( int ) const; - point_function_type point_function_value( &::SireMM::AngleRestraint::point ); - - AngleRestraint_exposer.def( - "point" - , point_function_value - , ( bp::arg("i") ) - , bp::return_value_policy() - , "Return the ith point" ); - - } - { //::SireMM::AngleRestraint::point0 - - typedef ::SireFF::Point const & ( ::SireMM::AngleRestraint::*point0_function_type)( ) const; - point0_function_type point0_function_value( &::SireMM::AngleRestraint::point0 ); - - AngleRestraint_exposer.def( - "point0" - , point0_function_value - , bp::return_value_policy() - , "Return the first point" ); - - } - { //::SireMM::AngleRestraint::point1 - - typedef ::SireFF::Point const & ( ::SireMM::AngleRestraint::*point1_function_type)( ) const; - point1_function_type point1_function_value( &::SireMM::AngleRestraint::point1 ); - - AngleRestraint_exposer.def( - "point1" - , point1_function_value - , bp::return_value_policy() - , "Return the second point" ); - - } - { //::SireMM::AngleRestraint::point2 - - typedef ::SireFF::Point const & ( ::SireMM::AngleRestraint::*point2_function_type)( ) const; - point2_function_type point2_function_value( &::SireMM::AngleRestraint::point2 ); - - AngleRestraint_exposer.def( - "point2" - , point2_function_value - , bp::return_value_policy() - , "Return the third point" ); - - } - { //::SireMM::AngleRestraint::setSpace + { //::SireMM::AngleRestraint::theta0 - typedef void ( ::SireMM::AngleRestraint::*setSpace_function_type)( ::SireVol::Space const & ) ; - setSpace_function_type setSpace_function_value( &::SireMM::AngleRestraint::setSpace ); + typedef ::SireUnits::Dimension::Angle ( ::SireMM::AngleRestraint::*theta0_function_type)( ) const; + theta0_function_type theta0_function_value( &::SireMM::AngleRestraint::theta0 ); AngleRestraint_exposer.def( - "setSpace" - , setSpace_function_value - , ( bp::arg("space") ) + "theta0" + , theta0_function_value , bp::release_gil_policy() - , "Set the space used to evaluate the energy of this restraint\nThrow: SireVol::incompatible_space\n" ); + , "Return the equilibrium angle for the restraint" ); } - { //::SireMM::AngleRestraint::theta + { //::SireMM::AngleRestraint::toString - typedef ::SireCAS::Symbol const & ( *theta_function_type )( ); - theta_function_type theta_function_value( &::SireMM::AngleRestraint::theta ); + typedef ::QString ( ::SireMM::AngleRestraint::*toString_function_type)( ) const; + toString_function_type toString_function_value( &::SireMM::AngleRestraint::toString ); AngleRestraint_exposer.def( - "theta" - , theta_function_value - , bp::return_value_policy() - , "Return the symbol that represents the angle between the points (theta)" ); + "toString" + , toString_function_value + , bp::release_gil_policy() + , "" ); } { //::SireMM::AngleRestraint::typeName @@ -303,69 +130,26 @@ void register_AngleRestraint_class(){ , "" ); } - { //::SireMM::AngleRestraint::update + { //::SireMM::AngleRestraint::what - typedef void ( ::SireMM::AngleRestraint::*update_function_type)( ::SireMol::MoleculeData const & ) ; - update_function_type update_function_value( &::SireMM::AngleRestraint::update ); + typedef char const * ( ::SireMM::AngleRestraint::*what_function_type)( ) const; + what_function_type what_function_value( &::SireMM::AngleRestraint::what ); AngleRestraint_exposer.def( - "update" - , update_function_value - , ( bp::arg("moldata") ) + "what" + , what_function_value , bp::release_gil_policy() - , "Update the points of this restraint using new molecule data from moldata\nThrow: SireBase::missing_property\nThrow: SireError::invalid_cast\nThrow: SireError::incompatible_error\n" ); - - } - { //::SireMM::AngleRestraint::update - - typedef void ( ::SireMM::AngleRestraint::*update_function_type)( ::SireMol::Molecules const & ) ; - update_function_type update_function_value( &::SireMM::AngleRestraint::update ); - - AngleRestraint_exposer.def( - "update" - , update_function_value - , ( bp::arg("molecules") ) - , bp::release_gil_policy() - , "Update the points of this restraint using new molecule data from molecules\nThrow: SireBase::missing_property\nThrow: SireError::invalid_cast\nThrow: SireError::incompatible_error\n" ); - - } - { //::SireMM::AngleRestraint::usesMoleculesIn - - typedef bool ( ::SireMM::AngleRestraint::*usesMoleculesIn_function_type)( ::SireFF::ForceTable const & ) const; - usesMoleculesIn_function_type usesMoleculesIn_function_value( &::SireMM::AngleRestraint::usesMoleculesIn ); - - AngleRestraint_exposer.def( - "usesMoleculesIn" - , usesMoleculesIn_function_value - , ( bp::arg("forcetable") ) - , bp::release_gil_policy() - , "Return whether or not this restraint involves any of the molecules\nthat are in the forcetable forcetable" ); - - } - { //::SireMM::AngleRestraint::usesMoleculesIn - - typedef bool ( ::SireMM::AngleRestraint::*usesMoleculesIn_function_type)( ::SireMol::Molecules const & ) const; - usesMoleculesIn_function_type usesMoleculesIn_function_value( &::SireMM::AngleRestraint::usesMoleculesIn ); - - AngleRestraint_exposer.def( - "usesMoleculesIn" - , usesMoleculesIn_function_value - , ( bp::arg("molecules") ) - , bp::release_gil_policy() - , "Return whether or not this restraint involves any of the molecules\nin molecules" ); + , "" ); } - AngleRestraint_exposer.staticmethod( "halfHarmonic" ); - AngleRestraint_exposer.staticmethod( "harmonic" ); - AngleRestraint_exposer.staticmethod( "theta" ); AngleRestraint_exposer.staticmethod( "typeName" ); AngleRestraint_exposer.def( "__copy__", &__copy__); AngleRestraint_exposer.def( "__deepcopy__", &__copy__); AngleRestraint_exposer.def( "clone", &__copy__); AngleRestraint_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::AngleRestraint >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); AngleRestraint_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::AngleRestraint >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); AngleRestraint_exposer.def_pickle(sire_pickle_suite< ::SireMM::AngleRestraint >()); AngleRestraint_exposer.def( "__str__", &__str__< ::SireMM::AngleRestraint > ); AngleRestraint_exposer.def( "__repr__", &__str__< ::SireMM::AngleRestraint > ); diff --git a/wrapper/MM/AngleRestraints.pypp.cpp b/wrapper/MM/AngleRestraints.pypp.cpp new file mode 100644 index 000000000..07edf38ca --- /dev/null +++ b/wrapper/MM/AngleRestraints.pypp.cpp @@ -0,0 +1,242 @@ +// This file has been generated by Py++. + +// (C) Christopher Woods, GPL >= 3 License + +#include "boost/python.hpp" +#include "Helpers/clone_const_reference.hpp" +#include "AngleRestraints.pypp.hpp" + +namespace bp = boost::python; + +#include "SireCAS/errors.h" + +#include "SireID/index.h" + +#include "SireStream/datastream.h" + +#include "SireStream/shareddatastream.h" + +#include "SireUnits/units.h" + +#include "anglerestraints.h" + +#include + +#include "anglerestraints.h" + +SireMM::AngleRestraints __copy__(const SireMM::AngleRestraints &other){ return SireMM::AngleRestraints(other); } + +#include "Helpers/copy.hpp" + +#include "Qt/qdatastream.hpp" + +#include "Helpers/str.hpp" + +#include "Helpers/release_gil_policy.hpp" + +#include "Helpers/len.hpp" + +void register_AngleRestraints_class(){ + + { //::SireMM::AngleRestraints + typedef bp::class_< SireMM::AngleRestraints, bp::bases< SireMM::Restraints, SireBase::Property > > AngleRestraints_exposer_t; + AngleRestraints_exposer_t AngleRestraints_exposer = AngleRestraints_exposer_t( "AngleRestraints", "This class represents a collection of angle restraints", bp::init< >("Null constructor") ); + bp::scope AngleRestraints_scope( AngleRestraints_exposer ); + AngleRestraints_exposer.def( bp::init< QString const & >(( bp::arg("name") ), "") ); + AngleRestraints_exposer.def( bp::init< SireMM::AngleRestraint const & >(( bp::arg("restraint") ), "") ); + AngleRestraints_exposer.def( bp::init< QList< SireMM::AngleRestraint > const & >(( bp::arg("restraints") ), "") ); + AngleRestraints_exposer.def( bp::init< QString const &, SireMM::AngleRestraint const & >(( bp::arg("name"), bp::arg("restraint") ), "") ); + AngleRestraints_exposer.def( bp::init< QString const &, QList< SireMM::AngleRestraint > const & >(( bp::arg("name"), bp::arg("restraints") ), "") ); + AngleRestraints_exposer.def( bp::init< SireMM::AngleRestraints const & >(( bp::arg("other") ), "") ); + { //::SireMM::AngleRestraints::add + + typedef void ( ::SireMM::AngleRestraints::*add_function_type)( ::SireMM::AngleRestraint const & ) ; + add_function_type add_function_value( &::SireMM::AngleRestraints::add ); + + AngleRestraints_exposer.def( + "add" + , add_function_value + , ( bp::arg("restraint") ) + , bp::release_gil_policy() + , "Add a restraint onto the list" ); + + } + { //::SireMM::AngleRestraints::add + + typedef void ( ::SireMM::AngleRestraints::*add_function_type)( ::SireMM::AngleRestraints const & ) ; + add_function_type add_function_value( &::SireMM::AngleRestraints::add ); + + AngleRestraints_exposer.def( + "add" + , add_function_value + , ( bp::arg("restraints") ) + , bp::release_gil_policy() + , "Add a restraint onto the list" ); + + } + { //::SireMM::AngleRestraints::at + + typedef ::SireMM::AngleRestraint const & ( ::SireMM::AngleRestraints::*at_function_type)( int ) const; + at_function_type at_function_value( &::SireMM::AngleRestraints::at ); + + AngleRestraints_exposer.def( + "at" + , at_function_value + , ( bp::arg("i") ) + , bp::return_value_policy() + , "Return the ith restraint" ); + + } + { //::SireMM::AngleRestraints::count + + typedef int ( ::SireMM::AngleRestraints::*count_function_type)( ) const; + count_function_type count_function_value( &::SireMM::AngleRestraints::count ); + + AngleRestraints_exposer.def( + "count" + , count_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + { //::SireMM::AngleRestraints::isEmpty + + typedef bool ( ::SireMM::AngleRestraints::*isEmpty_function_type)( ) const; + isEmpty_function_type isEmpty_function_value( &::SireMM::AngleRestraints::isEmpty ); + + AngleRestraints_exposer.def( + "isEmpty" + , isEmpty_function_value + , bp::release_gil_policy() + , "Return whether or not this is empty" ); + + } + { //::SireMM::AngleRestraints::isNull + + typedef bool ( ::SireMM::AngleRestraints::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::AngleRestraints::isNull ); + + AngleRestraints_exposer.def( + "isNull" + , isNull_function_value + , bp::release_gil_policy() + , "Return whether or not this is empty" ); + + } + { //::SireMM::AngleRestraints::nRestraints + + typedef int ( ::SireMM::AngleRestraints::*nRestraints_function_type)( ) const; + nRestraints_function_type nRestraints_function_value( &::SireMM::AngleRestraints::nRestraints ); + + AngleRestraints_exposer.def( + "nRestraints" + , nRestraints_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + AngleRestraints_exposer.def( bp::self != bp::self ); + AngleRestraints_exposer.def( bp::self + bp::other< SireMM::AngleRestraint >() ); + AngleRestraints_exposer.def( bp::self + bp::self ); + { //::SireMM::AngleRestraints::operator= + + typedef ::SireMM::AngleRestraints & ( ::SireMM::AngleRestraints::*assign_function_type)( ::SireMM::AngleRestraints const & ) ; + assign_function_type assign_function_value( &::SireMM::AngleRestraints::operator= ); + + AngleRestraints_exposer.def( + "assign" + , assign_function_value + , ( bp::arg("other") ) + , bp::return_self< >() + , "" ); + + } + AngleRestraints_exposer.def( bp::self == bp::self ); + { //::SireMM::AngleRestraints::operator[] + + typedef ::SireMM::AngleRestraint const & ( ::SireMM::AngleRestraints::*__getitem___function_type)( int ) const; + __getitem___function_type __getitem___function_value( &::SireMM::AngleRestraints::operator[] ); + + AngleRestraints_exposer.def( + "__getitem__" + , __getitem___function_value + , ( bp::arg("i") ) + , bp::return_value_policy() + , "" ); + + } + { //::SireMM::AngleRestraints::restraints + + typedef ::QList< SireMM::AngleRestraint > ( ::SireMM::AngleRestraints::*restraints_function_type)( ) const; + restraints_function_type restraints_function_value( &::SireMM::AngleRestraints::restraints ); + + AngleRestraints_exposer.def( + "restraints" + , restraints_function_value + , bp::release_gil_policy() + , "Return all of the restraints" ); + + } + { //::SireMM::AngleRestraints::size + + typedef int ( ::SireMM::AngleRestraints::*size_function_type)( ) const; + size_function_type size_function_value( &::SireMM::AngleRestraints::size ); + + AngleRestraints_exposer.def( + "size" + , size_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + { //::SireMM::AngleRestraints::toString + + typedef ::QString ( ::SireMM::AngleRestraints::*toString_function_type)( ) const; + toString_function_type toString_function_value( &::SireMM::AngleRestraints::toString ); + + AngleRestraints_exposer.def( + "toString" + , toString_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMM::AngleRestraints::typeName + + typedef char const * ( *typeName_function_type )( ); + typeName_function_type typeName_function_value( &::SireMM::AngleRestraints::typeName ); + + AngleRestraints_exposer.def( + "typeName" + , typeName_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMM::AngleRestraints::what + + typedef char const * ( ::SireMM::AngleRestraints::*what_function_type)( ) const; + what_function_type what_function_value( &::SireMM::AngleRestraints::what ); + + AngleRestraints_exposer.def( + "what" + , what_function_value + , bp::release_gil_policy() + , "" ); + + } + AngleRestraints_exposer.staticmethod( "typeName" ); + AngleRestraints_exposer.def( "__copy__", &__copy__); + AngleRestraints_exposer.def( "__deepcopy__", &__copy__); + AngleRestraints_exposer.def( "clone", &__copy__); + AngleRestraints_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::AngleRestraints >, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + AngleRestraints_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::AngleRestraints >, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + AngleRestraints_exposer.def_pickle(sire_pickle_suite< ::SireMM::AngleRestraints >()); + AngleRestraints_exposer.def( "__str__", &__str__< ::SireMM::AngleRestraints > ); + AngleRestraints_exposer.def( "__repr__", &__str__< ::SireMM::AngleRestraints > ); + AngleRestraints_exposer.def( "__len__", &__len_size< ::SireMM::AngleRestraints > ); + } + +} diff --git a/wrapper/MM/AngleRestraints.pypp.hpp b/wrapper/MM/AngleRestraints.pypp.hpp new file mode 100644 index 000000000..b1292a5cd --- /dev/null +++ b/wrapper/MM/AngleRestraints.pypp.hpp @@ -0,0 +1,10 @@ +// This file has been generated by Py++. + +// (C) Christopher Woods, GPL >= 3 License + +#ifndef AngleRestraints_hpp__pyplusplus_wrapper +#define AngleRestraints_hpp__pyplusplus_wrapper + +void register_AngleRestraints_class(); + +#endif//AngleRestraints_hpp__pyplusplus_wrapper diff --git a/wrapper/MM/CMakeAutogenFile.txt b/wrapper/MM/CMakeAutogenFile.txt index 3eeb5da11..3c3a4985a 100644 --- a/wrapper/MM/CMakeAutogenFile.txt +++ b/wrapper/MM/CMakeAutogenFile.txt @@ -68,6 +68,7 @@ set ( PYPP_SOURCES PositionalRestraints.pypp.cpp IntraGroupFF.pypp.cpp DihedralRestraint.pypp.cpp + DihedralRestraints.pypp.cpp IntraFF.pypp.cpp DihedralComponent.pypp.cpp ScaledChargeParameterNames3D.pypp.cpp @@ -181,6 +182,7 @@ set ( PYPP_SOURCES CoulombNBPairs.pypp.cpp AmberDihedral.pypp.cpp AngleRestraint.pypp.cpp + AngleRestraints.pypp.cpp TwoAtomFunctions.pypp.cpp IntraSoftCLJFFBase.pypp.cpp BondSymbols.pypp.cpp diff --git a/wrapper/MM/DihedralRestraint.pypp.cpp b/wrapper/MM/DihedralRestraint.pypp.cpp index 65ed2feb1..a5d04277e 100644 --- a/wrapper/MM/DihedralRestraint.pypp.cpp +++ b/wrapper/MM/DihedralRestraint.pypp.cpp @@ -3,36 +3,25 @@ // (C) Christopher Woods, GPL >= 3 License #include "boost/python.hpp" -#include "Helpers/clone_const_reference.hpp" #include "DihedralRestraint.pypp.hpp" namespace bp = boost::python; -#include "SireCAS/conditional.h" - #include "SireCAS/errors.h" -#include "SireCAS/power.h" - -#include "SireCAS/symbols.h" - -#include "SireCAS/values.h" - -#include "SireFF/forcetable.h" - #include "SireID/index.h" #include "SireStream/datastream.h" #include "SireStream/shareddatastream.h" -#include "SireUnits/angle.h" - #include "SireUnits/units.h" -#include "dihedralrestraint.h" +#include "dihedralrestraints.h" + +#include -#include "dihedralrestraint.h" +#include "dihedralrestraints.h" SireMM::DihedralRestraint __copy__(const SireMM::DihedralRestraint &other){ return SireMM::DihedralRestraint(other); } @@ -47,162 +36,50 @@ SireMM::DihedralRestraint __copy__(const SireMM::DihedralRestraint &other){ retu void register_DihedralRestraint_class(){ { //::SireMM::DihedralRestraint - typedef bp::class_< SireMM::DihedralRestraint, bp::bases< SireMM::Restraint3D, SireMM::Restraint, SireBase::Property > > DihedralRestraint_exposer_t; - DihedralRestraint_exposer_t DihedralRestraint_exposer = DihedralRestraint_exposer_t( "DihedralRestraint", "This is a restraint that operates on the dihedral angle between\nfour SireMM::Point objects (e.g. four atoms in a molecule)\n\nAuthor: Christopher Woods\n", bp::init< >("Constructor") ); + typedef bp::class_< SireMM::DihedralRestraint, bp::bases< SireBase::Property > > DihedralRestraint_exposer_t; + DihedralRestraint_exposer_t DihedralRestraint_exposer = DihedralRestraint_exposer_t( "DihedralRestraint", "This class represents a single angle restraint between any three\natoms in a system\nAuthor: Christopher Woods\n", bp::init< >("Null constructor") ); bp::scope DihedralRestraint_scope( DihedralRestraint_exposer ); - DihedralRestraint_exposer.def( bp::init< SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireCAS::Expression const & >(( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("point3"), bp::arg("restraint") ), "Construct a restraint that acts on the angle within the\nthree points point0, point1 and point2 (theta == a(012)),\nrestraining the angle within these points using the expression\nrestraint") ); - DihedralRestraint_exposer.def( bp::init< SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireCAS::Expression const &, SireCAS::Values const & >(( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("point3"), bp::arg("restraint"), bp::arg("values") ), "Construct a restraint that acts on the angle within the\nthree points point0, point1 and point2 (theta == a(012)),\nrestraining the angle within these points using the expression\nrestraint") ); + DihedralRestraint_exposer.def( bp::init< QList< long long > const &, SireUnits::Dimension::Angle const &, SireUnits::Dimension::HarmonicAngleConstant const & >(( bp::arg("atoms"), bp::arg("phi0"), bp::arg("kphi") ), "Construct a restraint that acts on the angle within the\nfour atoms atom0, atom1, atom2 atom3 (phi == a(0123)),\nrestraining the angle within these atoms") ); DihedralRestraint_exposer.def( bp::init< SireMM::DihedralRestraint const & >(( bp::arg("other") ), "Copy constructor") ); - { //::SireMM::DihedralRestraint::builtinSymbols + { //::SireMM::DihedralRestraint::atoms - typedef ::SireCAS::Symbols ( ::SireMM::DihedralRestraint::*builtinSymbols_function_type)( ) const; - builtinSymbols_function_type builtinSymbols_function_value( &::SireMM::DihedralRestraint::builtinSymbols ); + typedef ::QVector< long long > ( ::SireMM::DihedralRestraint::*atoms_function_type)( ) const; + atoms_function_type atoms_function_value( &::SireMM::DihedralRestraint::atoms ); DihedralRestraint_exposer.def( - "builtinSymbols" - , builtinSymbols_function_value + "atoms" + , atoms_function_value , bp::release_gil_policy() - , "Return the built-in symbols for this restraint" ); + , "Return the atoms involved in the restraint" ); } - { //::SireMM::DihedralRestraint::builtinValues + { //::SireMM::DihedralRestraint::isNull - typedef ::SireCAS::Values ( ::SireMM::DihedralRestraint::*builtinValues_function_type)( ) const; - builtinValues_function_type builtinValues_function_value( &::SireMM::DihedralRestraint::builtinValues ); + typedef bool ( ::SireMM::DihedralRestraint::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::DihedralRestraint::isNull ); DihedralRestraint_exposer.def( - "builtinValues" - , builtinValues_function_value + "isNull" + , isNull_function_value , bp::release_gil_policy() - , "Return the values of the built-in symbols of this restraint" ); - - } - { //::SireMM::DihedralRestraint::contains - - typedef bool ( ::SireMM::DihedralRestraint::*contains_function_type)( ::SireMol::MolNum ) const; - contains_function_type contains_function_value( &::SireMM::DihedralRestraint::contains ); - - DihedralRestraint_exposer.def( - "contains" - , contains_function_value - , ( bp::arg("molnum") ) - , bp::release_gil_policy() - , "Return whether or not this restraint affects the molecule\nwith number molnum" ); - - } - { //::SireMM::DihedralRestraint::contains - - typedef bool ( ::SireMM::DihedralRestraint::*contains_function_type)( ::SireMol::MolID const & ) const; - contains_function_type contains_function_value( &::SireMM::DihedralRestraint::contains ); - - DihedralRestraint_exposer.def( - "contains" - , contains_function_value - , ( bp::arg("molid") ) - , bp::release_gil_policy() - , "Return whether or not this restraint affects the molecule\nwith ID molid" ); - - } - { //::SireMM::DihedralRestraint::differentialRestraintFunction - - typedef ::SireCAS::Expression const & ( ::SireMM::DihedralRestraint::*differentialRestraintFunction_function_type)( ) const; - differentialRestraintFunction_function_type differentialRestraintFunction_function_value( &::SireMM::DihedralRestraint::differentialRestraintFunction ); - - DihedralRestraint_exposer.def( - "differentialRestraintFunction" - , differentialRestraintFunction_function_value - , bp::return_value_policy< bp::copy_const_reference >() - , "Return the function used to calculate the restraint force" ); - - } - { //::SireMM::DihedralRestraint::differentiate - - typedef ::SireMM::RestraintPtr ( ::SireMM::DihedralRestraint::*differentiate_function_type)( ::SireCAS::Symbol const & ) const; - differentiate_function_type differentiate_function_value( &::SireMM::DihedralRestraint::differentiate ); - - DihedralRestraint_exposer.def( - "differentiate" - , differentiate_function_value - , ( bp::arg("symbol") ) - , bp::release_gil_policy() - , "Return the differential of this restraint with respect to\nthe symbol symbol\nThrow: SireCAS::unavailable_differential\n" ); - - } - { //::SireMM::DihedralRestraint::force - - typedef void ( ::SireMM::DihedralRestraint::*force_function_type)( ::SireFF::MolForceTable &,double ) const; - force_function_type force_function_value( &::SireMM::DihedralRestraint::force ); - - DihedralRestraint_exposer.def( - "force" - , force_function_value - , ( bp::arg("forcetable"), bp::arg("scale_force")=1 ) - , "Calculate the force acting on the molecule in the forcetable forcetable\ncaused by this restraint, and add it on to the forcetable scaled by\nscale_force" ); - - } - { //::SireMM::DihedralRestraint::force - - typedef void ( ::SireMM::DihedralRestraint::*force_function_type)( ::SireFF::ForceTable &,double ) const; - force_function_type force_function_value( &::SireMM::DihedralRestraint::force ); - - DihedralRestraint_exposer.def( - "force" - , force_function_value - , ( bp::arg("forcetable"), bp::arg("scale_force")=1 ) - , "Calculate the force acting on the molecules in the forcetable forcetable\ncaused by this restraint, and add it on to the forcetable scaled by\nscale_force" ); - - } - { //::SireMM::DihedralRestraint::halfHarmonic - - typedef ::SireMM::DihedralRestraint ( *halfHarmonic_function_type )( ::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireUnits::Dimension::Angle const &,::SireMM::HarmonicAngleForceConstant const & ); - halfHarmonic_function_type halfHarmonic_function_value( &::SireMM::DihedralRestraint::halfHarmonic ); - - DihedralRestraint_exposer.def( - "halfHarmonic" - , halfHarmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("point3"), bp::arg("angle"), bp::arg("force_constant") ) - , bp::release_gil_policy() - , "Return a distance restraint that applied a half-harmonic potential\nbetween the points point0 and point1 above a distance distance\nusing a force constant force_constant" ); - - } - { //::SireMM::DihedralRestraint::harmonic - - typedef ::SireMM::DihedralRestraint ( *harmonic_function_type )( ::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireMM::HarmonicAngleForceConstant const & ); - harmonic_function_type harmonic_function_value( &::SireMM::DihedralRestraint::harmonic ); - - DihedralRestraint_exposer.def( - "harmonic" - , harmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("point3"), bp::arg("force_constant") ) - , bp::release_gil_policy() - , "Return a distance restraint that applies a harmonic potential between\nthe points point0 and point1 using a force constant force_constant" ); - - } - { //::SireMM::DihedralRestraint::molecules - - typedef ::SireMol::Molecules ( ::SireMM::DihedralRestraint::*molecules_function_type)( ) const; - molecules_function_type molecules_function_value( &::SireMM::DihedralRestraint::molecules ); - - DihedralRestraint_exposer.def( - "molecules" - , molecules_function_value - , bp::release_gil_policy() - , "Return the molecules used in this restraint" ); + , "" ); } - { //::SireMM::DihedralRestraint::nPoints + { //::SireMM::DihedralRestraint::kphi - typedef int ( ::SireMM::DihedralRestraint::*nPoints_function_type)( ) const; - nPoints_function_type nPoints_function_value( &::SireMM::DihedralRestraint::nPoints ); + typedef ::SireUnits::Dimension::HarmonicAngleConstant ( ::SireMM::DihedralRestraint::*kphi_function_type)( ) const; + kphi_function_type kphi_function_value( &::SireMM::DihedralRestraint::kphi ); DihedralRestraint_exposer.def( - "nPoints" - , nPoints_function_value + "kphi" + , kphi_function_value , bp::release_gil_policy() - , "This restraint involves four points" ); + , "Return the force constant for the restraint" ); } DihedralRestraint_exposer.def( bp::self != bp::self ); + DihedralRestraint_exposer.def( bp::self + bp::self ); + DihedralRestraint_exposer.def( bp::self + bp::other< SireMM::DihedralRestraints >() ); { //::SireMM::DihedralRestraint::operator= typedef ::SireMM::DihedralRestraint & ( ::SireMM::DihedralRestraint::*assign_function_type)( ::SireMM::DihedralRestraint const & ) ; @@ -217,90 +94,28 @@ void register_DihedralRestraint_class(){ } DihedralRestraint_exposer.def( bp::self == bp::self ); - { //::SireMM::DihedralRestraint::phi - - typedef ::SireCAS::Symbol const & ( *phi_function_type )( ); - phi_function_type phi_function_value( &::SireMM::DihedralRestraint::phi ); - - DihedralRestraint_exposer.def( - "phi" - , phi_function_value - , bp::return_value_policy() - , "Return the symbol that represents the dihedral angle between the points (phi)" ); - - } - { //::SireMM::DihedralRestraint::point - - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point_function_type)( int ) const; - point_function_type point_function_value( &::SireMM::DihedralRestraint::point ); - - DihedralRestraint_exposer.def( - "point" - , point_function_value - , ( bp::arg("i") ) - , bp::return_value_policy() - , "Return the ith point" ); - - } - { //::SireMM::DihedralRestraint::point0 - - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point0_function_type)( ) const; - point0_function_type point0_function_value( &::SireMM::DihedralRestraint::point0 ); - - DihedralRestraint_exposer.def( - "point0" - , point0_function_value - , bp::return_value_policy() - , "Return the first point" ); - - } - { //::SireMM::DihedralRestraint::point1 - - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point1_function_type)( ) const; - point1_function_type point1_function_value( &::SireMM::DihedralRestraint::point1 ); - - DihedralRestraint_exposer.def( - "point1" - , point1_function_value - , bp::return_value_policy() - , "Return the second point" ); - - } - { //::SireMM::DihedralRestraint::point2 + { //::SireMM::DihedralRestraint::phi0 - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point2_function_type)( ) const; - point2_function_type point2_function_value( &::SireMM::DihedralRestraint::point2 ); + typedef ::SireUnits::Dimension::Angle ( ::SireMM::DihedralRestraint::*phi0_function_type)( ) const; + phi0_function_type phi0_function_value( &::SireMM::DihedralRestraint::phi0 ); DihedralRestraint_exposer.def( - "point2" - , point2_function_value - , bp::return_value_policy() - , "Return the third point" ); - - } - { //::SireMM::DihedralRestraint::point3 - - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point3_function_type)( ) const; - point3_function_type point3_function_value( &::SireMM::DihedralRestraint::point3 ); - - DihedralRestraint_exposer.def( - "point3" - , point3_function_value - , bp::return_value_policy() - , "Return the fourth point" ); + "phi0" + , phi0_function_value + , bp::release_gil_policy() + , "Return the equilibrium angle for the restraint" ); } - { //::SireMM::DihedralRestraint::setSpace + { //::SireMM::DihedralRestraint::toString - typedef void ( ::SireMM::DihedralRestraint::*setSpace_function_type)( ::SireVol::Space const & ) ; - setSpace_function_type setSpace_function_value( &::SireMM::DihedralRestraint::setSpace ); + typedef ::QString ( ::SireMM::DihedralRestraint::*toString_function_type)( ) const; + toString_function_type toString_function_value( &::SireMM::DihedralRestraint::toString ); DihedralRestraint_exposer.def( - "setSpace" - , setSpace_function_value - , ( bp::arg("space") ) + "toString" + , toString_function_value , bp::release_gil_policy() - , "Set the space used to evaluate the energy of this restraint\nThrow: SireVol::incompatible_space\n" ); + , "" ); } { //::SireMM::DihedralRestraint::typeName @@ -315,69 +130,26 @@ void register_DihedralRestraint_class(){ , "" ); } - { //::SireMM::DihedralRestraint::update - - typedef void ( ::SireMM::DihedralRestraint::*update_function_type)( ::SireMol::MoleculeData const & ) ; - update_function_type update_function_value( &::SireMM::DihedralRestraint::update ); - - DihedralRestraint_exposer.def( - "update" - , update_function_value - , ( bp::arg("moldata") ) - , bp::release_gil_policy() - , "Update the points of this restraint using new molecule data from moldata\nThrow: SireBase::missing_property\nThrow: SireError::invalid_cast\nThrow: SireError::incompatible_error\n" ); - - } - { //::SireMM::DihedralRestraint::update + { //::SireMM::DihedralRestraint::what - typedef void ( ::SireMM::DihedralRestraint::*update_function_type)( ::SireMol::Molecules const & ) ; - update_function_type update_function_value( &::SireMM::DihedralRestraint::update ); + typedef char const * ( ::SireMM::DihedralRestraint::*what_function_type)( ) const; + what_function_type what_function_value( &::SireMM::DihedralRestraint::what ); DihedralRestraint_exposer.def( - "update" - , update_function_value - , ( bp::arg("molecules") ) + "what" + , what_function_value , bp::release_gil_policy() - , "Update the points of this restraint using new molecule data from molecules\nThrow: SireBase::missing_property\nThrow: SireError::invalid_cast\nThrow: SireError::incompatible_error\n" ); - - } - { //::SireMM::DihedralRestraint::usesMoleculesIn - - typedef bool ( ::SireMM::DihedralRestraint::*usesMoleculesIn_function_type)( ::SireFF::ForceTable const & ) const; - usesMoleculesIn_function_type usesMoleculesIn_function_value( &::SireMM::DihedralRestraint::usesMoleculesIn ); - - DihedralRestraint_exposer.def( - "usesMoleculesIn" - , usesMoleculesIn_function_value - , ( bp::arg("forcetable") ) - , bp::release_gil_policy() - , "Return whether or not this restraint involves any of the molecules\nthat are in the forcetable forcetable" ); - - } - { //::SireMM::DihedralRestraint::usesMoleculesIn - - typedef bool ( ::SireMM::DihedralRestraint::*usesMoleculesIn_function_type)( ::SireMol::Molecules const & ) const; - usesMoleculesIn_function_type usesMoleculesIn_function_value( &::SireMM::DihedralRestraint::usesMoleculesIn ); - - DihedralRestraint_exposer.def( - "usesMoleculesIn" - , usesMoleculesIn_function_value - , ( bp::arg("molecules") ) - , bp::release_gil_policy() - , "Return whether or not this restraint involves any of the molecules\nin molecules" ); + , "" ); } - DihedralRestraint_exposer.staticmethod( "halfHarmonic" ); - DihedralRestraint_exposer.staticmethod( "harmonic" ); - DihedralRestraint_exposer.staticmethod( "phi" ); DihedralRestraint_exposer.staticmethod( "typeName" ); DihedralRestraint_exposer.def( "__copy__", &__copy__); DihedralRestraint_exposer.def( "__deepcopy__", &__copy__); DihedralRestraint_exposer.def( "clone", &__copy__); DihedralRestraint_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::DihedralRestraint >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); DihedralRestraint_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::DihedralRestraint >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); DihedralRestraint_exposer.def_pickle(sire_pickle_suite< ::SireMM::DihedralRestraint >()); DihedralRestraint_exposer.def( "__str__", &__str__< ::SireMM::DihedralRestraint > ); DihedralRestraint_exposer.def( "__repr__", &__str__< ::SireMM::DihedralRestraint > ); diff --git a/wrapper/MM/DihedralRestraints.pypp.cpp b/wrapper/MM/DihedralRestraints.pypp.cpp new file mode 100644 index 000000000..44d49c86a --- /dev/null +++ b/wrapper/MM/DihedralRestraints.pypp.cpp @@ -0,0 +1,242 @@ +// This file has been generated by Py++. + +// (C) Christopher Woods, GPL >= 3 License + +#include "boost/python.hpp" +#include "Helpers/clone_const_reference.hpp" +#include "DihedralRestraints.pypp.hpp" + +namespace bp = boost::python; + +#include "SireCAS/errors.h" + +#include "SireID/index.h" + +#include "SireStream/datastream.h" + +#include "SireStream/shareddatastream.h" + +#include "SireUnits/units.h" + +#include "dihedralrestraints.h" + +#include + +#include "dihedralrestraints.h" + +SireMM::DihedralRestraints __copy__(const SireMM::DihedralRestraints &other){ return SireMM::DihedralRestraints(other); } + +#include "Helpers/copy.hpp" + +#include "Qt/qdatastream.hpp" + +#include "Helpers/str.hpp" + +#include "Helpers/release_gil_policy.hpp" + +#include "Helpers/len.hpp" + +void register_DihedralRestraints_class(){ + + { //::SireMM::DihedralRestraints + typedef bp::class_< SireMM::DihedralRestraints, bp::bases< SireMM::Restraints, SireBase::Property > > DihedralRestraints_exposer_t; + DihedralRestraints_exposer_t DihedralRestraints_exposer = DihedralRestraints_exposer_t( "DihedralRestraints", "This class represents a collection of angle restraints", bp::init< >("Null constructor") ); + bp::scope DihedralRestraints_scope( DihedralRestraints_exposer ); + DihedralRestraints_exposer.def( bp::init< QString const & >(( bp::arg("name") ), "") ); + DihedralRestraints_exposer.def( bp::init< SireMM::DihedralRestraint const & >(( bp::arg("restraint") ), "") ); + DihedralRestraints_exposer.def( bp::init< QList< SireMM::DihedralRestraint > const & >(( bp::arg("restraints") ), "") ); + DihedralRestraints_exposer.def( bp::init< QString const &, SireMM::DihedralRestraint const & >(( bp::arg("name"), bp::arg("restraint") ), "") ); + DihedralRestraints_exposer.def( bp::init< QString const &, QList< SireMM::DihedralRestraint > const & >(( bp::arg("name"), bp::arg("restraints") ), "") ); + DihedralRestraints_exposer.def( bp::init< SireMM::DihedralRestraints const & >(( bp::arg("other") ), "") ); + { //::SireMM::DihedralRestraints::add + + typedef void ( ::SireMM::DihedralRestraints::*add_function_type)( ::SireMM::DihedralRestraint const & ) ; + add_function_type add_function_value( &::SireMM::DihedralRestraints::add ); + + DihedralRestraints_exposer.def( + "add" + , add_function_value + , ( bp::arg("restraint") ) + , bp::release_gil_policy() + , "Add a restraints onto the list" ); + + } + { //::SireMM::DihedralRestraints::add + + typedef void ( ::SireMM::DihedralRestraints::*add_function_type)( ::SireMM::DihedralRestraints const & ) ; + add_function_type add_function_value( &::SireMM::DihedralRestraints::add ); + + DihedralRestraints_exposer.def( + "add" + , add_function_value + , ( bp::arg("restraints") ) + , bp::release_gil_policy() + , "Add a restraint onto the list" ); + + } + { //::SireMM::DihedralRestraints::at + + typedef ::SireMM::DihedralRestraint const & ( ::SireMM::DihedralRestraints::*at_function_type)( int ) const; + at_function_type at_function_value( &::SireMM::DihedralRestraints::at ); + + DihedralRestraints_exposer.def( + "at" + , at_function_value + , ( bp::arg("i") ) + , bp::return_value_policy() + , "Return the ith restraint" ); + + } + { //::SireMM::DihedralRestraints::count + + typedef int ( ::SireMM::DihedralRestraints::*count_function_type)( ) const; + count_function_type count_function_value( &::SireMM::DihedralRestraints::count ); + + DihedralRestraints_exposer.def( + "count" + , count_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + { //::SireMM::DihedralRestraints::isEmpty + + typedef bool ( ::SireMM::DihedralRestraints::*isEmpty_function_type)( ) const; + isEmpty_function_type isEmpty_function_value( &::SireMM::DihedralRestraints::isEmpty ); + + DihedralRestraints_exposer.def( + "isEmpty" + , isEmpty_function_value + , bp::release_gil_policy() + , "Return whether or not this is empty" ); + + } + { //::SireMM::DihedralRestraints::isNull + + typedef bool ( ::SireMM::DihedralRestraints::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::DihedralRestraints::isNull ); + + DihedralRestraints_exposer.def( + "isNull" + , isNull_function_value + , bp::release_gil_policy() + , "Return whether or not this is empty" ); + + } + { //::SireMM::DihedralRestraints::nRestraints + + typedef int ( ::SireMM::DihedralRestraints::*nRestraints_function_type)( ) const; + nRestraints_function_type nRestraints_function_value( &::SireMM::DihedralRestraints::nRestraints ); + + DihedralRestraints_exposer.def( + "nRestraints" + , nRestraints_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + DihedralRestraints_exposer.def( bp::self != bp::self ); + DihedralRestraints_exposer.def( bp::self + bp::other< SireMM::DihedralRestraint >() ); + DihedralRestraints_exposer.def( bp::self + bp::self ); + { //::SireMM::DihedralRestraints::operator= + + typedef ::SireMM::DihedralRestraints & ( ::SireMM::DihedralRestraints::*assign_function_type)( ::SireMM::DihedralRestraints const & ) ; + assign_function_type assign_function_value( &::SireMM::DihedralRestraints::operator= ); + + DihedralRestraints_exposer.def( + "assign" + , assign_function_value + , ( bp::arg("other") ) + , bp::return_self< >() + , "" ); + + } + DihedralRestraints_exposer.def( bp::self == bp::self ); + { //::SireMM::DihedralRestraints::operator[] + + typedef ::SireMM::DihedralRestraint const & ( ::SireMM::DihedralRestraints::*__getitem___function_type)( int ) const; + __getitem___function_type __getitem___function_value( &::SireMM::DihedralRestraints::operator[] ); + + DihedralRestraints_exposer.def( + "__getitem__" + , __getitem___function_value + , ( bp::arg("i") ) + , bp::return_value_policy() + , "" ); + + } + { //::SireMM::DihedralRestraints::restraints + + typedef ::QList< SireMM::DihedralRestraint > ( ::SireMM::DihedralRestraints::*restraints_function_type)( ) const; + restraints_function_type restraints_function_value( &::SireMM::DihedralRestraints::restraints ); + + DihedralRestraints_exposer.def( + "restraints" + , restraints_function_value + , bp::release_gil_policy() + , "Return all of the restraints" ); + + } + { //::SireMM::DihedralRestraints::size + + typedef int ( ::SireMM::DihedralRestraints::*size_function_type)( ) const; + size_function_type size_function_value( &::SireMM::DihedralRestraints::size ); + + DihedralRestraints_exposer.def( + "size" + , size_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + { //::SireMM::DihedralRestraints::toString + + typedef ::QString ( ::SireMM::DihedralRestraints::*toString_function_type)( ) const; + toString_function_type toString_function_value( &::SireMM::DihedralRestraints::toString ); + + DihedralRestraints_exposer.def( + "toString" + , toString_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMM::DihedralRestraints::typeName + + typedef char const * ( *typeName_function_type )( ); + typeName_function_type typeName_function_value( &::SireMM::DihedralRestraints::typeName ); + + DihedralRestraints_exposer.def( + "typeName" + , typeName_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMM::DihedralRestraints::what + + typedef char const * ( ::SireMM::DihedralRestraints::*what_function_type)( ) const; + what_function_type what_function_value( &::SireMM::DihedralRestraints::what ); + + DihedralRestraints_exposer.def( + "what" + , what_function_value + , bp::release_gil_policy() + , "" ); + + } + DihedralRestraints_exposer.staticmethod( "typeName" ); + DihedralRestraints_exposer.def( "__copy__", &__copy__); + DihedralRestraints_exposer.def( "__deepcopy__", &__copy__); + DihedralRestraints_exposer.def( "clone", &__copy__); + DihedralRestraints_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::DihedralRestraints >, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + DihedralRestraints_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::DihedralRestraints >, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + DihedralRestraints_exposer.def_pickle(sire_pickle_suite< ::SireMM::DihedralRestraints >()); + DihedralRestraints_exposer.def( "__str__", &__str__< ::SireMM::DihedralRestraints > ); + DihedralRestraints_exposer.def( "__repr__", &__str__< ::SireMM::DihedralRestraints > ); + DihedralRestraints_exposer.def( "__len__", &__len_size< ::SireMM::DihedralRestraints > ); + } + +} diff --git a/wrapper/MM/DihedralRestraints.pypp.hpp b/wrapper/MM/DihedralRestraints.pypp.hpp new file mode 100644 index 000000000..0aab6519a --- /dev/null +++ b/wrapper/MM/DihedralRestraints.pypp.hpp @@ -0,0 +1,10 @@ +// This file has been generated by Py++. + +// (C) Christopher Woods, GPL >= 3 License + +#ifndef DihedralRestraints_hpp__pyplusplus_wrapper +#define DihedralRestraints_hpp__pyplusplus_wrapper + +void register_DihedralRestraints_class(); + +#endif//DihedralRestraints_hpp__pyplusplus_wrapper diff --git a/wrapper/MM/SireMM_registrars.cpp b/wrapper/MM/SireMM_registrars.cpp index b7c70f668..9dabd1bd9 100644 --- a/wrapper/MM/SireMM_registrars.cpp +++ b/wrapper/MM/SireMM_registrars.cpp @@ -16,7 +16,7 @@ #include "selectordihedral.h" #include "clj14group.h" #include "improper.h" -#include "anglerestraint.h" +#include "anglerestraints.h" #include "bondrestraints.h" #include "twoatomfunctions.h" #include "boreschrestraints.h" @@ -54,7 +54,7 @@ #include "cljboxes.h" #include "internalgroupff.h" #include "internalcomponent.h" -#include "dihedralrestraint.h" +#include "dihedralrestraints.h" #include "selectorimproper.h" #include "intergroupff.h" #include "dihedral.h" @@ -112,6 +112,7 @@ void register_SireMM_objects() ObjectRegistry::registerConverterFor< SireMM::Improper >(); ObjectRegistry::registerConverterFor< SireMol::Mover >(); ObjectRegistry::registerConverterFor< SireMM::AngleRestraint >(); + ObjectRegistry::registerConverterFor< SireMM::AngleRestraints >(); ObjectRegistry::registerConverterFor< SireMM::BondRestraint >(); ObjectRegistry::registerConverterFor< SireMM::BondRestraints >(); ObjectRegistry::registerConverterFor< SireMM::TwoAtomFunctions >(); @@ -200,6 +201,7 @@ void register_SireMM_objects() ObjectRegistry::registerConverterFor< SireMM::Intra14Component >(); ObjectRegistry::registerConverterFor< SireMM::InternalComponent >(); ObjectRegistry::registerConverterFor< SireMM::DihedralRestraint >(); + ObjectRegistry::registerConverterFor< SireMM::DihedralRestraints >(); ObjectRegistry::registerConverterFor< SireMM::SelectorImproper >(); ObjectRegistry::registerConverterFor< SireMol::Mover >(); ObjectRegistry::registerConverterFor< SireMM::InterGroupFF >(); diff --git a/wrapper/MM/_MM.main.cpp b/wrapper/MM/_MM.main.cpp index 5a98d9755..addf5ef01 100644 --- a/wrapper/MM/_MM.main.cpp +++ b/wrapper/MM/_MM.main.cpp @@ -27,6 +27,8 @@ #include "AngleRestraint.pypp.hpp" +#include "AngleRestraints.pypp.hpp" + #include "AngleSymbols.pypp.hpp" #include "AtomFunction.pypp.hpp" @@ -161,6 +163,8 @@ #include "DihedralRestraint.pypp.hpp" +#include "DihedralRestraints.pypp.hpp" + #include "DihedralSymbols.pypp.hpp" #include "DistanceRestraint.pypp.hpp" @@ -552,6 +556,10 @@ BOOST_PYTHON_MODULE(_MM){ register_AngleRestraint_class(); + register_Restraints_class(); + + register_AngleRestraints_class(); + register_InternalSymbolsBase_class(); register_AngleSymbols_class(); @@ -586,8 +594,6 @@ BOOST_PYTHON_MODULE(_MM){ register_BondRestraint_class(); - register_Restraints_class(); - register_BondRestraints_class(); register_BondSymbols_class(); @@ -682,6 +688,8 @@ BOOST_PYTHON_MODULE(_MM){ register_DihedralRestraint_class(); + register_DihedralRestraints_class(); + register_DihedralSymbols_class(); register_DistanceRestraint_class(); diff --git a/wrapper/MM/active_headers.h b/wrapper/MM/active_headers.h index 3cb17a25a..c50bbe8c8 100644 --- a/wrapper/MM/active_headers.h +++ b/wrapper/MM/active_headers.h @@ -5,7 +5,7 @@ #include "amberparams.h" #include "angle.h" -#include "anglerestraint.h" +#include "anglerestraints.h" #include "atomfunctions.h" #include "atomljs.h" #include "bond.h" @@ -30,7 +30,7 @@ #include "cljworkspace.h" #include "coulombpotential.h" #include "dihedral.h" -#include "dihedralrestraint.h" +#include "dihedralrestraints.h" #include "distancerestraint.h" #include "excludedpairs.h" #include "fouratomfunctions.h" diff --git a/wrapper/Move/OpenMMFrEnergyDT.pypp.cpp b/wrapper/Move/OpenMMFrEnergyDT.pypp.cpp index 98270db10..2967c7709 100644 --- a/wrapper/Move/OpenMMFrEnergyDT.pypp.cpp +++ b/wrapper/Move/OpenMMFrEnergyDT.pypp.cpp @@ -309,6 +309,16 @@ void register_OpenMMFrEnergyDT_class(){ , "" ); } + { //::SireMove::OpenMMFrEnergyDT::getMCBarostat_membrane + + typedef bool ( ::SireMove::OpenMMFrEnergyDT::*getMCBarostat_membrane_function_type)( ) ; + getMCBarostat_membrane_function_type getMCBarostat_membrane_function_value( &::SireMove::OpenMMFrEnergyDT::getMCBarostat_membrane ); + + OpenMMFrEnergyDT_exposer.def( + "getMCBarostat_membrane" + , getMCBarostat_membrane_function_value + , "" ); + } { //::SireMove::OpenMMFrEnergyDT::getMCBarostat_frequency typedef int ( ::SireMove::OpenMMFrEnergyDT::*getMCBarostat_frequency_function_type)( ) ; @@ -614,6 +624,18 @@ void register_OpenMMFrEnergyDT_class(){ , bp::release_gil_policy() , "Set Monte Carlo Barostat onoff" ); + } + { //::SireMove::OpenMMFrEnergyDT::setMCBarostat_membrane + + typedef void ( ::SireMove::OpenMMFrEnergyDT::*setMCBarostat_membrane_function_type)( bool ) ; + setMCBarostat_membrane_function_type setMCBarostat_membrane_function_value( &::SireMove::OpenMMFrEnergyDT::setMCBarostat_membrane ); + + OpenMMFrEnergyDT_exposer.def( + "setMCBarostat_membrane" + , setMCBarostat_membrane_function_value + , ( bp::arg("arg0") ) + , "Set Monte Carlo Membrane Barostat onoff (only has an effect if Monte Carlo Barostat is on)" ); + } { //::SireMove::OpenMMFrEnergyDT::setMCBarostat_frequency diff --git a/wrapper/Move/OpenMMFrEnergyST.pypp.cpp b/wrapper/Move/OpenMMFrEnergyST.pypp.cpp index 7bd2c4502..a7bf2cf5d 100644 --- a/wrapper/Move/OpenMMFrEnergyST.pypp.cpp +++ b/wrapper/Move/OpenMMFrEnergyST.pypp.cpp @@ -411,6 +411,17 @@ void register_OpenMMFrEnergyST_class(){ , bp::release_gil_policy() , "" ); + } + { //::SireMove::OpenMMFrEnergyST::getMCBarostatMembrane + + typedef bool ( ::SireMove::OpenMMFrEnergyST::*getMCBarostatMembrane_function_type)( ) ; + getMCBarostatMembrane_function_type getMCBarostatMembrane_function_value( &::SireMove::OpenMMFrEnergyST::getMCBarostatMembrane ); + + OpenMMFrEnergyST_exposer.def( + "getMCBarostatMembrane" + , getMCBarostatMembrane_function_value + , "" ); + } { //::SireMove::OpenMMFrEnergyST::getMCBarostatFrequency @@ -869,6 +880,18 @@ void register_OpenMMFrEnergyST_class(){ , bp::release_gil_policy() , "Set Monte Carlo Barostat onoff" ); + } + { //::SireMove::OpenMMFrEnergyST::setMCBarostatMembrane + + typedef void ( ::SireMove::OpenMMFrEnergyST::*setMCBarostatMembrane_function_type)( bool ) ; + setMCBarostatMembrane_function_type setMCBarostatMembrane_function_value( &::SireMove::OpenMMFrEnergyST::setMCBarostatMembrane ); + + OpenMMFrEnergyST_exposer.def( + "setMCBarostatMembrane" + , setMCBarostatMembrane_function_value + , ( bp::arg("arg0") ) + , "Set Monte Carlo Membrane Barostat onoff (only has an effect if Monte Carlo Barostat is on)" ); + } { //::SireMove::OpenMMFrEnergyST::setMCBarostatFrequency @@ -880,7 +903,7 @@ void register_OpenMMFrEnergyST_class(){ , setMCBarostatFrequency_function_value , ( bp::arg("arg0") ) , bp::release_gil_policy() - , "Set the Monte Carlo Barostat frequency in time speps" ); + , "Set the Monte Carlo Barostat frequency in time steps" ); } { //::SireMove::OpenMMFrEnergyST::setPlatform diff --git a/wrapper/Move/OpenMMMDIntegrator.pypp.cpp b/wrapper/Move/OpenMMMDIntegrator.pypp.cpp index d8a3e6b89..e577211bf 100644 --- a/wrapper/Move/OpenMMMDIntegrator.pypp.cpp +++ b/wrapper/Move/OpenMMMDIntegrator.pypp.cpp @@ -322,6 +322,17 @@ void register_OpenMMMDIntegrator_class(){ , "" ); } + { //::SireMove::OpenMMMDIntegrator::getMCBarostatMembrane + + typedef bool ( ::SireMove::OpenMMMDIntegrator::*getMCBarostatMembrane_function_type)( ) ; + getMCBarostatMembrane_function_type getMCBarostatMembrane_function_value( &::SireMove::OpenMMMDIntegrator::getMCBarostatMembrane ); + + OpenMMMDIntegrator_exposer.def( + "getMCBarostatMembrane" + , getMCBarostatMembrane_function_value + , "" ); + } + { //::SireMove::OpenMMMDIntegrator::getMCBarostatFrequency typedef int ( ::SireMove::OpenMMMDIntegrator::*getMCBarostatFrequency_function_type)( ) ; @@ -677,6 +688,18 @@ void register_OpenMMMDIntegrator_class(){ , bp::release_gil_policy() , "Set Monte Carlo Barostat onoff" ); + } + { //::SireMove::OpenMMMDIntegrator::setMCBarostatMembrane + + typedef void ( ::SireMove::OpenMMMDIntegrator::*setMCBarostatMembrane_function_type)( bool ) ; + setMCBarostatMembrane_function_type setMCBarostatMembrane_function_value( &::SireMove::OpenMMMDIntegrator::setMCBarostatMembrane ); + + OpenMMMDIntegrator_exposer.def( + "setMCBarostatMembrane" + , setMCBarostatMembrane_function_value + , ( bp::arg("arg0") ) + , "Set Monte Carlo Membrane Barostat onoff (only has an effect if Monte Carlo Barostat is on)" ); + } { //::SireMove::OpenMMMDIntegrator::setMCBarostatFrequency diff --git a/wrapper/Move/OpenMMPMEFEP.pypp.cpp b/wrapper/Move/OpenMMPMEFEP.pypp.cpp index c70629b6f..7b786da06 100644 --- a/wrapper/Move/OpenMMPMEFEP.pypp.cpp +++ b/wrapper/Move/OpenMMPMEFEP.pypp.cpp @@ -414,6 +414,17 @@ void register_OpenMMPMEFEP_class(){ , "" ); } + { //::SireMove::OpenMMPMEFEP::getMCBarostatMembrane + + typedef bool ( ::SireMove::OpenMMPMEFEP::*getMCBarostatMembrane_function_type)( ) ; + getMCBarostatMembrane_function_type getMCBarostatMembrane_function_value( &::SireMove::OpenMMPMEFEP::getMCBarostatMembrane ); + + OpenMMPMEFEP_exposer.def( + "getMCBarostatMembrane" + , getMCBarostatMembrane_function_value + , bp::release_gil_policy() + , "" ); + } { //::SireMove::OpenMMPMEFEP::getMCBarostatFrequency typedef int ( ::SireMove::OpenMMPMEFEP::*getMCBarostatFrequency_function_type)( ) ; @@ -858,6 +869,19 @@ void register_OpenMMPMEFEP_class(){ , bp::release_gil_policy() , "Set Monte Carlo Barostat onoff" ); + } + { //::SireMove::OpenMMPMEFEP::setMCBarostatMembrane + + typedef void ( ::SireMove::OpenMMPMEFEP::*setMCBarostatMembrane_function_type)( bool ) ; + setMCBarostatMembrane_function_type setMCBarostatMembrane_function_value( &::SireMove::OpenMMPMEFEP::setMCBarostatMembrane ); + + OpenMMPMEFEP_exposer.def( + "setMCBarostatMembrane" + , setMCBarostatMembrane_function_value + , ( bp::arg("arg0") ) + , bp::release_gil_policy() + , "Set Monte Carlo Membrane Barostat onoff (only has an effect if Monte Carlo Barostat is on)" ); + } { //::SireMove::OpenMMPMEFEP::setMCBarostatFrequency diff --git a/wrapper/Tools/OpenMMMD.py b/wrapper/Tools/OpenMMMD.py index 45f7017ee..97fe9e401 100644 --- a/wrapper/Tools/OpenMMMD.py +++ b/wrapper/Tools/OpenMMMD.py @@ -268,6 +268,12 @@ """Whether or not to use a barostat (needed for NPT simulation).""", ) +barostat_membrane = Parameter( + "membrane barostat", + False, + """Whether the barostat is a membrane barostat (needed for membrane simulation).""", +) + andersen_frequency = Parameter( "andersen frequency", 10.0, """Collision frequency in units of (1/ps)""" ) @@ -824,6 +830,8 @@ def setupMoves(system, debug_seed, GPUS): if barostat.val: Integrator_OpenMM.setPressure(pressure.val) Integrator_OpenMM.setMCBarostat(barostat.val) + if barostat_membrane.val: + Integrator_OpenMM.setMCBarostatMembrane(barostat_membrane.val) Integrator_OpenMM.setMCBarostatFrequency(barostat_frequency.val) # print Integrator_OpenMM.getDeviceIndex() @@ -2145,6 +2153,8 @@ def setupMovesFreeEnergy(system, debug_seed, gpu_idx, lam_val): if barostat.val: Integrator_OpenMM.setPressure(pressure.val) Integrator_OpenMM.setMCBarostat(barostat.val) + if barostat_membrane.val: + Integrator_OpenMM.setMCBarostatMembrane(barostat_membrane.val) Integrator_OpenMM.setMCBarostatFrequency(barostat_frequency.val) # Choose a random seed for Sire if a debugging seed hasn't been set.