From 05e38a9d45d402e071eb6e0b4779b897519bd954 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Wed, 10 Jul 2024 11:14:00 +0100 Subject: [PATCH 01/15] Initial update for angle restraints --- corelib/src/libs/SireMM/anglerestraint.cpp | 1108 ++++++++++++++------ corelib/src/libs/SireMM/anglerestraint.h | 155 +-- 2 files changed, 863 insertions(+), 400 deletions(-) diff --git a/corelib/src/libs/SireMM/anglerestraint.cpp b/corelib/src/libs/SireMM/anglerestraint.cpp index a765049c0..06e6d536e 100644 --- a/corelib/src/libs/SireMM/anglerestraint.cpp +++ b/corelib/src/libs/SireMM/anglerestraint.cpp @@ -27,16 +27,16 @@ #include "anglerestraint.h" -#include "SireFF/forcetable.h" +// #include "SireFF/forcetable.h" -#include "SireCAS/conditional.h" -#include "SireCAS/power.h" -#include "SireCAS/symbols.h" -#include "SireCAS/values.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/angle.h" #include "SireUnits/units.h" #include "SireStream/datastream.h" @@ -44,12 +44,14 @@ #include "SireCAS/errors.h" +#include + using namespace SireMM; -using namespace SireMol; -using namespace SireFF; +// using namespace SireMol; +// using namespace SireFF; using namespace SireID; using namespace SireBase; -using namespace SireCAS; +// using namespace SireCAS; using namespace SireMaths; using namespace SireStream; using namespace SireUnits; @@ -62,14 +64,14 @@ using namespace SireUnits::Dimension; 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); + sds << angrest.atms << angrest._theta0 << angrest._ktheta; return ds; } @@ -83,11 +85,7 @@ QDataStream &operator>>(QDataStream &ds, AngleRestraint &angrest) { 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]); + sds >> angrest.atms >> angrest._theta0 >> angrest._ktheta; } else throw version_error(v, "1", r_angrest, CODELOC); @@ -95,461 +93,903 @@ QDataStream &operator>>(QDataStream &ds, AngleRestraint &angrest) 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() +/** Null constructor */ +AngleRestraint::AngleRestraint() + : ConcreteProperty(), + _ktheta(0), _theta0(0) { } -void AngleRestraint::calculateTheta() -{ - if (this->restraintFunction().isFunction(theta())) +/** 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) +{ + // Need to think here about validating the angle and force constant values + // if (atoms.count() != 3) + // { + // throw SireError::invalid_arg(QObject::tr( + // "Wrong number of inputs for an Angle restraint. You need to " + // "provide 3 atoms (%1).") + // .arg(atoms.count()), + // // .arg(theta0.count()) + // // .arg(ktheta.count()), + // CODELOC); + // } + + // make sure that we have 3 distinct atoms + QSet distinct; + distinct.reserve(3); + + for (const auto &atom : atoms) { - 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); + if (atom >= 0) + distinct.insert(atom); } -} - -/** 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 (distinct.count() != 3) + // throw SireError::invalid_arg(QObject::tr( + // "There is something wrong with the atoms provided. " + // "They should all be unique and all greater than or equal to 0."), + // CODELOC); - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); + atms = atoms.toVector(); +} - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]); +/* Copy constructor*/ +AngleRestraint::AngleRestraint(const AngleRestraint &other) + : ConcreteProperty(other), + atms(other.atms), _theta0(other._theta0), _ktheta(other._ktheta) - 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) +/* Destructor */ +AngleRestraint::~AngleRestraint() { - p[0] = point0; - p[1] = point1; - p[2] = point2; +} - force_expression = this->restraintFunction().differentiate(theta()); +AngleRestraint &AngleRestraint::operator=(const AngleRestraint &other) +{ + if (this != &other) + { + Property::operator=(other); + atms = other.atms; + _theta0 = other._theta0; + _ktheta = other._ktheta; + } - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); + return *this; +} - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]); +bool AngleRestraint::operator==(const AngleRestraint &other) const +{ + return atms == other.atms and + _theta0 == other._theta0 and + _ktheta == other._ktheta; +} - this->calculateTheta(); +bool AngleRestraint::operator!=(const AngleRestraint &other) const +{ + return not operator==(other); } -/** 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) +AngleRestraints AngleRestraint::operator+(const AngleRestraint &other) const { - p[0] = point0; - p[1] = point1; - p[2] = point2; + return AngleRestraints(*this) + other; +} - 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); - } +AngleRestraints AngleRestraint::operator+(const AngleRestraints &other) const +{ + return AngleRestraints(*this) + other; +} - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]); +const char *AngleRestraint::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); +} - this->calculateTheta(); +const char *AngleRestraint::what() const +{ + return AngleRestraint::typeName(); } -/** Copy constructor */ -AngleRestraint::AngleRestraint(const AngleRestraint &other) - : ConcreteProperty(other), force_expression(other.force_expression), - intra_molecule_points(other.intra_molecule_points) +AngleRestraint *AngleRestraint::clone() const { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } + return new AngleRestraint(*this); } -/** Destructor */ -AngleRestraint::~AngleRestraint() +bool AngleRestraint::isNull() const { + return atms.isEmpty(); } -/** Copy assignment operator */ -AngleRestraint &AngleRestraint::operator=(const AngleRestraint &other) +QString AngleRestraint::toString() const { - if (this != &other) + if (this->isNull()) + return QObject::tr("AngleRestraint::null"); + else { - ExpressionRestraint3D::operator=(other); + QStringList a; - for (int i = 0; i < this->nPoints(); ++i) + for (const auto &atom : atms) { - p[i] = other.p[i]; + a.append(QString::number(atom)); } - - force_expression = other.force_expression; - intra_molecule_points = other.intra_molecule_points; + return QString("AngleRestraint( [%1], theta0=%2, ktheta=%3 )") + .arg(a.join(", ")) + .arg(_theta0.toString()) + .arg(_ktheta.toString()); } - - return *this; } -/** Comparison operator */ -bool AngleRestraint::operator==(const AngleRestraint &other) const +/** Return the force constant for the restraint */ +SireUnits::Dimension::HarmonicAngleConstant AngleRestraint::ktheta() 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); + return this->_ktheta; } -/** Comparison operator */ -bool AngleRestraint::operator!=(const AngleRestraint &other) const +/** Return the equilibrium angle for the restraint */ +SireUnits::Dimension::Angle AngleRestraint::theta0() const { - return not AngleRestraint::operator==(other); + return this->_theta0; } -const char *AngleRestraint::typeName() +/** Return the atoms involved in the restraint */ +QVector AngleRestraint::atoms() const { - return QMetaType::typeName(qMetaTypeId()); + return this->atms; } -/** This restraint involves three points */ -int AngleRestraint::nPoints() const +/////// +/////// Implementation of AngleRestraints +/////// + +/** Serialise to a binary datastream */ + +static const RegisterMetaType r_angrests; + +QDataStream &operator<<(QDataStream &ds, const AngleRestraints &angrests) { - return 3; + writeHeader(ds, r_angrests, 1); + + SharedDataStream sds(ds); + + sds << angrests.r + << static_cast(angrests); + + return ds; } -/** Return the ith point */ -const Point &AngleRestraint::point(int i) const +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, AngleRestraints &angrests) { - i = Index(i).map(this->nPoints()); + VersionID v = readHeader(ds, r_angrests); + + if (v == 1) + { + SharedDataStream sds(ds); - return p[i].read(); + sds >> angrests.r >> + static_cast(angrests); + } + else + throw version_error(v, "1", r_angrests, CODELOC); + + return ds; } -/** Return the first point */ -const Point &AngleRestraint::point0() const +/** Null constructor */ +AngleRestraints::AngleRestraints() + : ConcreteProperty() { - return p[0].read(); } -/** Return the second point */ -const Point &AngleRestraint::point1() const +AngleRestraints::AngleRestraints(const QString &name) + : ConcreteProperty(name) { - return p[1].read(); } -/** Return the third point */ -const Point &AngleRestraint::point2() const +AngleRestraints::AngleRestraints(const AngleRestraint &restraint) + : ConcreteProperty() { - return p[2].read(); + if (not restraint.isNull()) + r.append(restraint); } -/** Return the built-in symbols of this restraint */ -Symbols AngleRestraint::builtinSymbols() const +AngleRestraints::AngleRestraints(const QList &restraints) + : ConcreteProperty() { - if (this->restraintFunction().isFunction(theta())) - return theta(); - else - return Symbols(); + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } } -/** Return the values of the built-in symbols of this restraint */ -Values AngleRestraint::builtinValues() const +AngleRestraints::AngleRestraints(const QString &name, + const AngleRestraint &restraint) + : ConcreteProperty(name) { - if (this->restraintFunction().isFunction(theta())) - return theta() == this->values()[theta()]; - else - return Values(); + if (not restraint.isNull()) + r.append(restraint); } -/** Return the differential of this restraint with respect to - the symbol 'symbol' +AngleRestraints::AngleRestraints(const QString &name, + const QList &restraints) + : ConcreteProperty(name) +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } +} - \throw SireCAS::unavailable_differential -*/ -RestraintPtr AngleRestraint::differentiate(const Symbol &symbol) const +AngleRestraints::AngleRestraints(const AngleRestraints &other) + : ConcreteProperty(other), r(other.r) { - 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 +/* Desctructor */ +AngleRestraints::~AngleRestraints() +{ +} - \throw SireVol::incompatible_space -*/ -void AngleRestraint::setSpace(const Space &new_space) +AngleRestraints &AngleRestraints::operator=(const AngleRestraints &other) { - if (not this->space().equals(new_space)) + if (this != &other) { - 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; - } + Restraints::operator=(other); + r = other.r; } + + return *this; } -/** Return the function used to calculate the restraint force */ -const Expression &AngleRestraint::differentialRestraintFunction() const +bool AngleRestraints::operator==(const AngleRestraints &other) const { - return force_expression; + return r == other.r and Restraints::operator==(other); } -/** 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 AngleRestraints::operator!=(const AngleRestraints &other) 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); + return not operator==(other); } -/** 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 +const char *AngleRestraints::typeName() { - 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; + return QMetaType::typeName(qMetaTypeId()); +} - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by an angle restraint."), - CODELOC); +const char *AngleRestraints::what() const +{ + return AngleRestraints::typeName(); } -/** Update the points of this restraint using new molecule data from 'moldata' +AngleRestraints *AngleRestraints::clone() const +{ + return new AngleRestraints(*this); +} - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void AngleRestraint::update(const MoleculeData &moldata) +QString AngleRestraints::toString() const { - if (this->contains(moldata.number())) - { - AngleRestraint old_state(*this); + if (this->isEmpty()) + return QObject::tr("AngleRestraints::null"); - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(moldata); - } + QStringList parts; - this->calculateTheta(); - } - catch (...) + const auto n = this->count(); + + if (n <= 10) + { + for (int i = 0; i < n; i++) { - AngleRestraint::operator=(old_state); - throw; + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); } } -} - -/** 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)) + else { - AngleRestraint old_state(*this); - - try + for (int i = 0; i < 5; i++) { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(molecules); - } - - this->calculateTheta(); + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); } - catch (...) + + parts.append("..."); + + for (int i = n - 5; i < n; i++) { - AngleRestraint::operator=(old_state); - throw; + 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 the molecules used in this restraint */ -Molecules AngleRestraint::molecules() const +/** Return whether or not this is empty */ +bool AngleRestraints::isEmpty() const { - Molecules mols; + return this->r.isEmpty(); +} - for (int i = 0; i < this->nPoints(); ++i) - { - mols += p[i].read().molecules(); - } +/** Return whether or not this is empty */ +bool AngleRestraints::isNull() const +{ + return this->isEmpty(); +} - return mols; +/** Return the number of restraints */ +int AngleRestraints::nRestraints() const +{ + return this->r.count(); } -/** Return whether or not this restraint affects the molecule - with number 'molnum' */ -bool AngleRestraint::contains(MolNum molnum) const +/** Return the number of restraints */ +int AngleRestraints::count() const { - return p[0].read().contains(molnum) or p[1].read().contains(molnum) or p[2].read().contains(molnum); + return this->nRestraints(); } -/** Return whether or not this restraint affects the molecule - with ID 'molid' */ -bool AngleRestraint::contains(const MolID &molid) const +/** Return the number of restraints */ +int AngleRestraints::size() const { - return p[0].read().contains(molid) or p[1].read().contains(molid) or p[2].read().contains(molid); + return this->nRestraints(); } -/** 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 the ith restraint */ +const AngleRestraint &AngleRestraints::at(int i) const { - return p[0].read().usesMoleculesIn(forcetable) or p[1].read().usesMoleculesIn(forcetable) or - p[2].read().usesMoleculesIn(forcetable); + i = SireID::Index(i).map(this->r.count()); + + return this->r.at(i); } -/** Return whether or not this restraint involves any of the molecules - in 'molecules' */ -bool AngleRestraint::usesMoleculesIn(const Molecules &molecules) const +/** Return the ith restraint */ +const AngleRestraint &AngleRestraints::operator[](int i) const { - return p[0].read().usesMoleculesIn(molecules) or p[1].read().usesMoleculesIn(molecules) or - p[2].read().usesMoleculesIn(molecules); + return this->at(i); } -static Expression harmonicFunction(double force_constant) +/** Return all of the restraints */ +QList AngleRestraints::restraints() const { - if (SireMaths::isZero(force_constant)) - return 0; - else - return force_constant * pow(AngleRestraint::theta(), 2); + return this->r; } -static Expression diffHarmonicFunction(double force_constant) +/** Add a restraints onto the list */ +void AngleRestraints::add(const AngleRestraint &restraint) { - if (SireMaths::isZero(force_constant)) - return 0; - else - return (2 * force_constant) * AngleRestraint::theta(); + if (not restraint.isNull()) + r.append(restraint); } -/** 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) +/** Add a restraint onto the list */ +void AngleRestraints::add(const AngleRestraints &restraints) { - return AngleRestraint(point0, point1, point2, ::harmonicFunction(force_constant), - ::diffHarmonicFunction(force_constant)); + this->r += restraints.r; } -static Expression halfHarmonicFunction(double force_constant, double angle) +/** Add a restraint onto the list */ +AngleRestraints &AngleRestraints::operator+=(const AngleRestraint &restraint) { - if (SireMaths::isZero(force_constant)) - return 0; + this->add(restraint); + return *this; +} - else if (angle <= 0) - // this is just a harmonic function - return ::harmonicFunction(force_constant); +/** Add a restraint onto the list */ +AngleRestraints AngleRestraints::operator+(const AngleRestraint &restraint) const +{ + AngleRestraints ret(*this); + ret += restraint; + return *this; +} - else - { - const Symbol &theta = AngleRestraint::theta(); - return Conditional(GreaterThan(theta, angle), force_constant * pow(theta - angle, 2), 0); - } +/** Add restraints onto the list */ +AngleRestraints &AngleRestraints::operator+=(const AngleRestraints &restraints) +{ + this->add(restraints); + return *this; } -static Expression diffHalfHarmonicFunction(double force_constant, double angle) +/** Add restraints onto the list */ +AngleRestraints AngleRestraints::operator+(const AngleRestraints &restraints) const { - if (SireMaths::isZero(force_constant)) - return 0; + AngleRestraints ret(*this); + ret += restraints; + return *this; +} - else if (angle <= 0) - // this is just a harmonic function - return ::diffHarmonicFunction(force_constant); +// QDataStream &operator<<(QDataStream &ds, const AngleRestraint &angrest) +// { +// writeHeader(ds, r_angrest, 1); - else - { - const Symbol &theta = AngleRestraint::theta(); - return Conditional(GreaterThan(theta, angle), (2 * force_constant) * (theta - angle), 0); - } -} +// SharedDataStream sds(ds); -/** 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); +// sds << angrest.p[0] << angrest.p[1] << angrest.p[2] << angrest.force_expression +// << static_cast(angrest); - return AngleRestraint(point0, point1, point2, ::halfHarmonicFunction(force_constant, acute_angle), - ::diffHalfHarmonicFunction(force_constant, acute_angle)); -} +// 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); +// } +// } + +// 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 index 6901c78ca..c8a4f228e 100644 --- a/corelib/src/libs/SireMM/anglerestraint.h +++ b/corelib/src/libs/SireMM/anglerestraint.h @@ -28,56 +28,49 @@ #ifndef SIREMM_ANGLERESTRAINT_H #define SIREMM_ANGLERESTRAINT_H -#include "SireFF/point.h" +// #include "SireFF/point.h" -#include "restraint.h" +#include "restraints.h" -#include "SireCAS/expression.h" -#include "SireCAS/symbol.h" +// #include "SireCAS/expression.h" +// #include "SireCAS/symbol.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 { - 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 + /** 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 AngleRestraint &); - friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, AngleRestraint &); + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::AngleRestraint &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::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 QList &atoms, + const SireUnits::Dimension::Angle &theta0, + const SireUnits::Dimension::HarmonicAngleConstant &ktheta); AngleRestraint(const AngleRestraint &other); @@ -88,72 +81,102 @@ namespace SireMM 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; - const Point &point(int i) const; + AngleRestraint *clone() const; - const Point &point0() const; - const Point &point1() const; - const Point &point2() const; + QString toString() const; - int nPoints() const; + bool isNull() const; - static const Symbol &theta(); + QVector atoms() const; - Symbols builtinSymbols() const; - Values builtinValues() const; + SireUnits::Dimension::Angle theta0() const; + SireUnits::Dimension::HarmonicAngleConstant ktheta() const; - RestraintPtr differentiate(const Symbol &symbol) const; + private: + /** Atoms involved in the angle restraint */ + QVector atms; - void setSpace(const Space &space); + /** Equilibrium angle */ + SireUnits::Dimension::Angle _theta0; - const Expression &differentialRestraintFunction() const; + /** Harmonic angle constant */ + SireUnits::Dimension::HarmonicAngleConstant _ktheta; + }; - void force(MolForceTable &forcetable, double scale_force = 1) const; - void force(ForceTable &forcetable, double scale_force = 1) const; + /** 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 &); - void update(const MoleculeData &moldata); + public: + AngleRestraints(); + AngleRestraints(const QString &name); - void update(const Molecules &molecules); + AngleRestraints(const AngleRestraint &restraint); + AngleRestraints(const QList &restraints); - Molecules molecules() const; + AngleRestraints(const QString &name, + const AngleRestraint &restraint); - bool contains(MolNum molnum) const; - bool contains(const MolID &molid) const; + AngleRestraints(const QString &name, + const QList &restraints); - bool usesMoleculesIn(const ForceTable &forcetable) const; - bool usesMoleculesIn(const Molecules &molecules) const; + AngleRestraints(const AngleRestraints &other); - static AngleRestraint harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const HarmonicAngleForceConstant &force_constant); + ~AngleRestraints(); - static AngleRestraint halfHarmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const SireUnits::Dimension::Angle &angle, - const HarmonicAngleForceConstant &force_constant); + AngleRestraints &operator=(const AngleRestraints &other); - protected: - AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Expression &nrg_restraint, const Expression &force_restraint); + bool operator==(const AngleRestraints &other) const; + bool operator!=(const AngleRestraints &other) const; - private: - void calculateTheta(); + static const char *typeName(); + const char *what() const; - /** The three points between which the restraint is calculated */ - SireFF::PointPtr p[3]; + AngleRestraints *clone() const; - /** The expression used to calculate the force */ - Expression force_expression; + QString toString() const; - /** Whether or not all three points are within the same molecule */ - bool intra_molecule_points; - }; + 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); -} // namespace SireMM + 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 From 12af94a60299cda9aa7b90aff0eed9cf547dafa8 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Wed, 10 Jul 2024 18:14:45 +0100 Subject: [PATCH 02/15] Intial addition of alchemical dihedral restraints and angle restraints update --- corelib/src/libs/SireMM/anglerestraint.cpp | 2 +- corelib/src/libs/SireMM/dihedralrestraint.cpp | 639 ++++++++---------- corelib/src/libs/SireMM/dihedralrestraint.h | 152 +++-- wrapper/MM/AngleRestraints.pypp.cpp | 254 +++++++ wrapper/MM/AngleRestraints.pypp.hpp | 10 + wrapper/MM/DihedralRestraints.pypp.cpp | 254 +++++++ wrapper/MM/DihedralRestraints.pypp.hpp | 10 + 7 files changed, 910 insertions(+), 411 deletions(-) create mode 100644 wrapper/MM/AngleRestraints.pypp.cpp create mode 100644 wrapper/MM/AngleRestraints.pypp.hpp create mode 100644 wrapper/MM/DihedralRestraints.pypp.cpp create mode 100644 wrapper/MM/DihedralRestraints.pypp.hpp diff --git a/corelib/src/libs/SireMM/anglerestraint.cpp b/corelib/src/libs/SireMM/anglerestraint.cpp index 06e6d536e..00e93dde5 100644 --- a/corelib/src/libs/SireMM/anglerestraint.cpp +++ b/corelib/src/libs/SireMM/anglerestraint.cpp @@ -770,7 +770,7 @@ AngleRestraints AngleRestraints::operator+(const AngleRestraints &restraints) co // p[i].edit().setSpace(new_space); // } -// Restraint3D::setSpace(new_space); +// Restraint::setSpace(new_space); // this->calculateTheta(); // } diff --git a/corelib/src/libs/SireMM/dihedralrestraint.cpp b/corelib/src/libs/SireMM/dihedralrestraint.cpp index bea13ed94..68ed915f5 100644 --- a/corelib/src/libs/SireMM/dihedralrestraint.cpp +++ b/corelib/src/libs/SireMM/dihedralrestraint.cpp @@ -27,16 +27,16 @@ #include "dihedralrestraint.h" -#include "SireFF/forcetable.h" +// #include "SireFF/forcetable.h" -#include "SireCAS/conditional.h" -#include "SireCAS/power.h" -#include "SireCAS/symbols.h" -#include "SireCAS/values.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/angle.h" #include "SireUnits/units.h" #include "SireStream/datastream.h" @@ -44,12 +44,14 @@ #include "SireCAS/errors.h" +#include + using namespace SireMM; -using namespace SireMol; -using namespace SireFF; +// using namespace SireMol; +// using namespace SireFF; using namespace SireID; using namespace SireBase; -using namespace SireCAS; +// using namespace SireCAS; using namespace SireMaths; using namespace SireStream; using namespace SireUnits; @@ -62,14 +64,14 @@ using namespace SireUnits::Dimension; 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); + sds << dihrest.atms << dihrest._phi0 << dihrest._kphi; return ds; } @@ -83,12 +85,7 @@ QDataStream &operator>>(QDataStream &ds, DihedralRestraint &dihrest) { 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]); + sds >> dihrest.atms >> dihrest._phi0 >> dihrest._kphi; } else throw version_error(v, "1", r_dihrest, CODELOC); @@ -96,163 +93,101 @@ QDataStream &operator>>(QDataStream &ds, DihedralRestraint &dihrest) 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() +/** Null constructor */ +DihedralRestraint::DihedralRestraint() + : ConcreteProperty(), + _phi0(0), _kphi(0) { - 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 + 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) +{ + // Need to think here about validating the angle and force constant values + // if (atoms.count() != 3) + // { + // throw SireError::invalid_arg(QObject::tr( + // "Wrong number of inputs for an Angle restraint. You need to " + // "provide 3 atoms (%1).") + // .arg(atoms.count()), + // // .arg(phi0.count()) + // // .arg(kphi.count()), + // CODELOC); + // } + + // make sure that we have 3 distinct atoms + QSet distinct; + distinct.reserve(4); + + for (const auto &atom : atoms) { - 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); + if (atom >= 0) + distinct.insert(atom); } - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]) and - Point::areIntraMoleculePoints(p[0], p[3]); + // if (distinct.count() != 3) + // throw SireError::invalid_arg(QObject::tr( + // "There is something wrong with the atoms provided. " + // "They should all be unique and all greater than or equal to 0."), + // CODELOC); - this->calculatePhi(); + atms = atoms.toVector(); } -/** Copy constructor */ +/* Copy constructor*/ DihedralRestraint::DihedralRestraint(const DihedralRestraint &other) - : ConcreteProperty(other), force_expression(other.force_expression), - intra_molecule_points(other.intra_molecule_points) + : ConcreteProperty(other), + atms(other.atms), _phi0(other._phi0), _kphi(other._kphi) + { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } } -/** Destructor */ +/* 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; + Property::operator=(other); + atms = other.atms; + _phi0 = other._phi0; + _kphi = other._kphi; } 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); + return atms == other.atms and + _phi0 == other._phi0 and + _kphi == other._kphi; } -/** Comparison operator */ bool DihedralRestraint::operator!=(const DihedralRestraint &other) const { - return not DihedralRestraint::operator==(other); + 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() @@ -260,317 +195,315 @@ const char *DihedralRestraint::typeName() return QMetaType::typeName(qMetaTypeId()); } -/** This restraint involves four points */ -int DihedralRestraint::nPoints() const +const char *DihedralRestraint::what() const { - return 4; + return DihedralRestraint::typeName(); } -/** Return the ith point */ -const Point &DihedralRestraint::point(int i) const +DihedralRestraint *DihedralRestraint::clone() const { - i = Index(i).map(this->nPoints()); + return new DihedralRestraint(*this); +} - return p[i].read(); +bool DihedralRestraint::isNull() const +{ + return atms.isEmpty(); } -/** Return the first point */ -const Point &DihedralRestraint::point0() const +QString DihedralRestraint::toString() const { - return p[0].read(); + 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 second point */ -const Point &DihedralRestraint::point1() const +/** Return the force constant for the restraint */ +SireUnits::Dimension::HarmonicAngleConstant DihedralRestraint::kphi() const { - return p[1].read(); + return this->_kphi; } -/** Return the third point */ -const Point &DihedralRestraint::point2() const +/** Return the equilibrium angle for the restraint */ +SireUnits::Dimension::Angle DihedralRestraint::phi0() const { - return p[2].read(); + return this->_phi0; } -/** Return the fourth point */ -const Point &DihedralRestraint::point3() const +/** Return the atoms involved in the restraint */ +QVector DihedralRestraint::atoms() const { - return p[3].read(); + return this->atms; } -/** Return the built-in symbols for this restraint */ -Symbols DihedralRestraint::builtinSymbols() const +/////// +/////// Implementation of DihedralRestraints +/////// + +/** Serialise to a binary datastream */ + +static const RegisterMetaType r_dihrests; + +QDataStream &operator<<(QDataStream &ds, const DihedralRestraints &dihrests) { - if (this->restraintFunction().isFunction(phi())) - return phi(); - else - return Symbols(); + writeHeader(ds, r_dihrests, 1); + + SharedDataStream sds(ds); + + sds << dihrests.r + << static_cast(dihrests); + + return ds; } -/** Return the values of the built-in symbols of this restraint */ -Values DihedralRestraint::builtinValues() const +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, DihedralRestraints &dihrests) { - if (this->restraintFunction().isFunction(phi())) - return phi() == this->values()[phi()]; + VersionID v = readHeader(ds, r_dihrests); + + if (v == 1) + { + SharedDataStream sds(ds); + + sds >> dihrests.r >> + static_cast(dihrests); + } else - return Values(); + throw version_error(v, "1", r_dihrests, CODELOC); + + return ds; } -/** Return the differential of this restraint with respect to - the symbol 'symbol' +/** Null constructor */ +DihedralRestraints::DihedralRestraints() + : ConcreteProperty() +{ +} - \throw SireCAS::unavailable_differential -*/ -RestraintPtr DihedralRestraint::differentiate(const Symbol &symbol) const +DihedralRestraints::DihedralRestraints(const QString &name) + : ConcreteProperty(name) { - 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 +DihedralRestraints::DihedralRestraints(const DihedralRestraint &restraint) + : ConcreteProperty() +{ + if (not restraint.isNull()) + r.append(restraint); +} - \throw SireVol::incompatible_space -*/ -void DihedralRestraint::setSpace(const Space &new_space) +DihedralRestraints::DihedralRestraints(const QList &restraints) + : ConcreteProperty() { - if (not this->space().equals(new_space)) + for (const auto &restraint : restraints) { - DihedralRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().setSpace(new_space); - } + if (not restraint.isNull()) + r.append(restraint); + } +} - Restraint3D::setSpace(new_space); +DihedralRestraints::DihedralRestraints(const QString &name, + const DihedralRestraint &restraint) + : ConcreteProperty(name) +{ + if (not restraint.isNull()) + r.append(restraint); +} - this->calculatePhi(); - } - catch (...) - { - DihedralRestraint::operator=(old_state); - throw; - } +DihedralRestraints::DihedralRestraints(const QString &name, + const QList &restraints) + : ConcreteProperty(name) +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); } } -/** Return the function used to calculate the restraint force */ -const Expression &DihedralRestraint::differentialRestraintFunction() const +DihedralRestraints::DihedralRestraints(const DihedralRestraints &other) + : ConcreteProperty(other), r(other.r) { - 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 +/* Desctructor */ +DihedralRestraints::~DihedralRestraints() { - 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; +DihedralRestraints &DihedralRestraints::operator=(const DihedralRestraints &other) +{ + if (this != &other) + { + Restraints::operator=(other); + r = other.r; + } - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by a dihedral restraint."), - CODELOC); + return *this; } -/** 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 DihedralRestraints::operator==(const DihedralRestraints &other) 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); + return r == other.r and Restraints::operator==(other); +} - if (not(in_p0 or in_p1 or in_p2 or in_p3)) - // this molecule is not affected by the restraint - return; +bool DihedralRestraints::operator!=(const DihedralRestraints &other) const +{ + return not operator==(other); +} - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by a dihedral restraint."), - CODELOC); +const char *DihedralRestraints::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); } -/** Update the points of this restraint using new molecule data from 'moldata' +const char *DihedralRestraints::what() const +{ + return DihedralRestraints::typeName(); +} - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void DihedralRestraint::update(const MoleculeData &moldata) +DihedralRestraints *DihedralRestraints::clone() const { - if (this->contains(moldata.number())) - { - DihedralRestraint old_state(*this); + return new DihedralRestraints(*this); +} - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(moldata); - } +QString DihedralRestraints::toString() const +{ + if (this->isEmpty()) + return QObject::tr("DihedralRestraints::null"); - this->calculatePhi(); - } - catch (...) - { - DihedralRestraint::operator=(old_state); - throw; - } - } -} + QStringList parts; -/** Update the points of this restraint using new molecule data from 'molecules' + const auto n = this->count(); - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void DihedralRestraint::update(const Molecules &molecules) -{ - if (this->usesMoleculesIn(molecules)) + if (n <= 10) { - DihedralRestraint old_state(*this); - - try + for (int i = 0; i < n; i++) { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(molecules); - } - - this->calculatePhi(); + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); } - catch (...) + } + else + { + for (int i = 0; i < 5; i++) { - DihedralRestraint::operator=(old_state); - throw; + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); } - } -} -/** Return the molecules used in this restraint */ -Molecules DihedralRestraint::molecules() const -{ - Molecules mols; + parts.append("..."); - for (int i = 0; i < this->nPoints(); ++i) - { - mols += p[i].read().molecules(); + for (int i = n - 5; i < n; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } } - return mols; + 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 restraint affects the molecule - with number 'molnum' */ -bool DihedralRestraint::contains(MolNum molnum) const +/** Return whether or not this is empty */ +bool DihedralRestraints::isEmpty() 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 this->r.isEmpty(); } -/** Return whether or not this restraint affects the molecule - with ID 'molid' */ -bool DihedralRestraint::contains(const MolID &molid) const +/** Return whether or not this is empty */ +bool DihedralRestraints::isNull() 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 this->isEmpty(); } -/** 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 the number of restraints */ +int DihedralRestraints::nRestraints() 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 this->r.count(); } -/** Return whether or not this restraint involves any of the molecules - in 'molecules' */ -bool DihedralRestraint::usesMoleculesIn(const Molecules &molecules) const +/** Return the number of restraints */ +int DihedralRestraints::count() 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); + return this->nRestraints(); } -static Expression harmonicFunction(double force_constant) +/** Return the number of restraints */ +int DihedralRestraints::size() const { - if (SireMaths::isZero(force_constant)) - return 0; - else - return force_constant * pow(DihedralRestraint::phi(), 2); + return this->nRestraints(); } -static Expression diffHarmonicFunction(double force_constant) +/** Return the ith restraint */ +const DihedralRestraint &DihedralRestraints::at(int i) const { - if (SireMaths::isZero(force_constant)) - return 0; - else - return (2 * force_constant) * DihedralRestraint::phi(); + i = SireID::Index(i).map(this->r.count()); + + return this->r.at(i); } -/** 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 the ith restraint */ +const DihedralRestraint &DihedralRestraints::operator[](int i) const { - return DihedralRestraint(point0, point1, point2, point3, ::harmonicFunction(force_constant), - ::diffHarmonicFunction(force_constant)); + return this->at(i); } -static Expression halfHarmonicFunction(double force_constant, double angle) +/** Return all of the restraints */ +QList DihedralRestraints::restraints() const { - 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); - } + return this->r; } -static Expression diffHalfHarmonicFunction(double force_constant, double angle) +/** Add a restraints onto the list */ +void DihedralRestraints::add(const DihedralRestraint &restraint) { - if (SireMaths::isZero(force_constant)) - return 0; + if (not restraint.isNull()) + r.append(restraint); +} - else if (angle <= 0) - // this is just a harmonic function - return ::diffHarmonicFunction(force_constant); +/** Add a restraint onto the list */ +void DihedralRestraints::add(const DihedralRestraints &restraints) +{ + this->r += restraints.r; +} - else - { - const Symbol &phi = DihedralRestraint::phi(); - return Conditional(GreaterThan(phi, angle), (2 * force_constant) * (phi - angle), 0); - } +/** Add a restraint onto the list */ +DihedralRestraints &DihedralRestraints::operator+=(const DihedralRestraint &restraint) +{ + this->add(restraint); + return *this; } -/** 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) +/** Add a restraint onto the list */ +DihedralRestraints DihedralRestraints::operator+(const DihedralRestraint &restraint) const { - double ang = angle.to(radians); + DihedralRestraints ret(*this); + ret += restraint; + return *this; +} - return DihedralRestraint(point0, point1, point2, point3, ::halfHarmonicFunction(force_constant, ang), - ::diffHalfHarmonicFunction(force_constant, ang)); +/** 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/dihedralrestraint.h b/corelib/src/libs/SireMM/dihedralrestraint.h index 606ee3b78..b340a3411 100644 --- a/corelib/src/libs/SireMM/dihedralrestraint.h +++ b/corelib/src/libs/SireMM/dihedralrestraint.h @@ -28,41 +28,49 @@ #ifndef SIREMM_DIHEDRALRESTRAINT_H #define SIREMM_DIHEDRALRESTRAINT_H -#include "anglerestraint.h" -#include "restraint.h" +// #include "SireFF/point.h" + +#include "restraints.h" + +// #include "SireCAS/expression.h" +// #include "SireCAS/symbol.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 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 + /** This class represents a single angle restraint between any three + * atoms in a system + * @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 &); + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::DihedralRestraint &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::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 QList &atoms, + const SireUnits::Dimension::Angle &phi0, + const SireUnits::Dimension::HarmonicAngleConstant &kphi); DihedralRestraint(const DihedralRestraint &other); @@ -73,72 +81,102 @@ namespace SireMM 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; - const Point &point(int i) const; + DihedralRestraint *clone() const; - const Point &point0() const; - const Point &point1() const; - const Point &point2() const; - const Point &point3() const; + QString toString() const; - int nPoints() const; + bool isNull() const; - static const Symbol &phi(); + QVector atoms() const; - Symbols builtinSymbols() const; - Values builtinValues() const; + SireUnits::Dimension::Angle phi0() const; + SireUnits::Dimension::HarmonicAngleConstant kphi() const; - RestraintPtr differentiate(const Symbol &symbol) const; + private: + /** Atoms involved in the angle restraint */ + QVector atms; - void setSpace(const Space &space); + /** Equilibrium angle */ + SireUnits::Dimension::Angle _phi0; - const Expression &differentialRestraintFunction() const; + /** Harmonic angle constant */ + SireUnits::Dimension::HarmonicAngleConstant _kphi; + }; - void force(MolForceTable &forcetable, double scale_force = 1) const; - void force(ForceTable &forcetable, double scale_force = 1) const; + /** 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 &); - void update(const MoleculeData &moldata); - void update(const Molecules &molecules); + public: + DihedralRestraints(); + DihedralRestraints(const QString &name); - Molecules molecules() const; + DihedralRestraints(const DihedralRestraint &restraint); + DihedralRestraints(const QList &restraints); - bool contains(MolNum molnum) const; - bool contains(const MolID &molid) const; + DihedralRestraints(const QString &name, + const DihedralRestraint &restraint); - bool usesMoleculesIn(const ForceTable &forcetable) const; - bool usesMoleculesIn(const Molecules &molecules) const; + DihedralRestraints(const QString &name, + const QList &restraints); - static DihedralRestraint harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const HarmonicAngleForceConstant &force_constant); + DihedralRestraints(const DihedralRestraints &other); - static DihedralRestraint halfHarmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const SireUnits::Dimension::Angle &angle, - const HarmonicAngleForceConstant &force_constant); + ~DihedralRestraints(); - protected: - DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const PointRef &point3, - const Expression &nrg_restraint, const Expression &force_restraint); + DihedralRestraints &operator=(const DihedralRestraints &other); - private: - void calculatePhi(); + bool operator==(const DihedralRestraints &other) const; + bool operator!=(const DihedralRestraints &other) const; - /** The four points between which the restraint is calculated */ - SireFF::PointPtr p[4]; + static const char *typeName(); + const char *what() const; - /** The expression used to calculate the force */ - Expression force_expression; + DihedralRestraints *clone() const; - /** Whether or not all four points are within the same molecule */ - bool intra_molecule_points; - }; + 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); -} // namespace SireMM + 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 +#endif \ No newline at end of file diff --git a/wrapper/MM/AngleRestraints.pypp.cpp b/wrapper/MM/AngleRestraints.pypp.cpp new file mode 100644 index 000000000..fdf557c21 --- /dev/null +++ b/wrapper/MM/AngleRestraints.pypp.cpp @@ -0,0 +1,254 @@ +// 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/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 + +#include "anglerestraint.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/DihedralRestraints.pypp.cpp b/wrapper/MM/DihedralRestraints.pypp.cpp new file mode 100644 index 000000000..21f56e7e0 --- /dev/null +++ b/wrapper/MM/DihedralRestraints.pypp.cpp @@ -0,0 +1,254 @@ +// 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/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 + +#include "dihedralrestraint.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 From 8809ac5a4e9163fad41fbb65cef4030d11673930 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Tue, 16 Jul 2024 15:45:14 +0100 Subject: [PATCH 03/15] Update wrappers for angle and dihedral restraints classes --- wrapper/MM/AngleRestraint.pypp.cpp | 296 ++++--------------------- wrapper/MM/DihedralRestraint.pypp.cpp | 308 ++++---------------------- wrapper/MM/_MM.main.cpp | 38 +++- 3 files changed, 125 insertions(+), 517 deletions(-) diff --git a/wrapper/MM/AngleRestraint.pypp.cpp b/wrapper/MM/AngleRestraint.pypp.cpp index 052916bb4..33954da56 100644 --- a/wrapper/MM/AngleRestraint.pypp.cpp +++ b/wrapper/MM/AngleRestraint.pypp.cpp @@ -3,7 +3,6 @@ // (C) Christopher Woods, GPL >= 3 License #include "boost/python.hpp" -#include "Helpers/clone_const_reference.hpp" #include "AngleRestraint.pypp.hpp" namespace bp = boost::python; @@ -32,6 +31,8 @@ namespace bp = boost::python; #include "anglerestraint.h" +#include + #include "anglerestraint.h" SireMM::AngleRestraint __copy__(const SireMM::AngleRestraint &other){ return SireMM::AngleRestraint(other); } @@ -47,162 +48,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::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 ); - - AngleRestraint_exposer.def( - "halfHarmonic" - , halfHarmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), 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::AngleRestraint::atoms - } - { //::SireMM::AngleRestraint::harmonic - - 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 ::QVector< long long > ( ::SireMM::AngleRestraint::*atoms_function_type)( ) const; + atoms_function_type atoms_function_value( &::SireMM::AngleRestraint::atoms ); AngleRestraint_exposer.def( - "harmonic" - , harmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("force_constant") ) + "atoms" + , atoms_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" ); + , "Return the atoms involved in the restraint" ); } - { //::SireMM::AngleRestraint::molecules + { //::SireMM::AngleRestraint::isNull - typedef ::SireMol::Molecules ( ::SireMM::AngleRestraint::*molecules_function_type)( ) const; - molecules_function_type molecules_function_value( &::SireMM::AngleRestraint::molecules ); + typedef bool ( ::SireMM::AngleRestraint::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::AngleRestraint::isNull ); AngleRestraint_exposer.def( - "molecules" - , molecules_function_value + "isNull" + , isNull_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 +106,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 + { //::SireMM::AngleRestraint::theta0 - typedef ::SireFF::Point const & ( ::SireMM::AngleRestraint::*point2_function_type)( ) const; - point2_function_type point2_function_value( &::SireMM::AngleRestraint::point2 ); + typedef ::SireUnits::Dimension::Angle ( ::SireMM::AngleRestraint::*theta0_function_type)( ) const; + theta0_function_type theta0_function_value( &::SireMM::AngleRestraint::theta0 ); AngleRestraint_exposer.def( - "point2" - , point2_function_value - , bp::return_value_policy() - , "Return the third point" ); - - } - { //::SireMM::AngleRestraint::setSpace - - typedef void ( ::SireMM::AngleRestraint::*setSpace_function_type)( ::SireVol::Space const & ) ; - setSpace_function_type setSpace_function_value( &::SireMM::AngleRestraint::setSpace ); - - 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 +142,26 @@ void register_AngleRestraint_class(){ , "" ); } - { //::SireMM::AngleRestraint::update - - typedef void ( ::SireMM::AngleRestraint::*update_function_type)( ::SireMol::MoleculeData const & ) ; - update_function_type update_function_value( &::SireMM::AngleRestraint::update ); - - AngleRestraint_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::AngleRestraint::update + { //::SireMM::AngleRestraint::what - typedef void ( ::SireMM::AngleRestraint::*update_function_type)( ::SireMol::Molecules 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("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::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/DihedralRestraint.pypp.cpp b/wrapper/MM/DihedralRestraint.pypp.cpp index 65ed2feb1..accddf41c 100644 --- a/wrapper/MM/DihedralRestraint.pypp.cpp +++ b/wrapper/MM/DihedralRestraint.pypp.cpp @@ -3,7 +3,6 @@ // (C) Christopher Woods, GPL >= 3 License #include "boost/python.hpp" -#include "Helpers/clone_const_reference.hpp" #include "DihedralRestraint.pypp.hpp" namespace bp = boost::python; @@ -32,6 +31,8 @@ namespace bp = boost::python; #include "dihedralrestraint.h" +#include + #include "dihedralrestraint.h" SireMM::DihedralRestraint __copy__(const SireMM::DihedralRestraint &other){ return SireMM::DihedralRestraint(other); } @@ -47,162 +48,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 - - typedef ::SireCAS::Symbols ( ::SireMM::DihedralRestraint::*builtinSymbols_function_type)( ) const; - builtinSymbols_function_type builtinSymbols_function_value( &::SireMM::DihedralRestraint::builtinSymbols ); - - DihedralRestraint_exposer.def( - "builtinSymbols" - , builtinSymbols_function_value - , bp::release_gil_policy() - , "Return the built-in symbols for this restraint" ); - - } - { //::SireMM::DihedralRestraint::builtinValues - - typedef ::SireCAS::Values ( ::SireMM::DihedralRestraint::*builtinValues_function_type)( ) const; - builtinValues_function_type builtinValues_function_value( &::SireMM::DihedralRestraint::builtinValues ); - - DihedralRestraint_exposer.def( - "builtinValues" - , builtinValues_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 + { //::SireMM::DihedralRestraint::atoms - typedef ::SireCAS::Expression const & ( ::SireMM::DihedralRestraint::*differentialRestraintFunction_function_type)( ) const; - differentialRestraintFunction_function_type differentialRestraintFunction_function_value( &::SireMM::DihedralRestraint::differentialRestraintFunction ); + typedef ::QVector< long long > ( ::SireMM::DihedralRestraint::*atoms_function_type)( ) const; + atoms_function_type atoms_function_value( &::SireMM::DihedralRestraint::atoms ); 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") ) + "atoms" + , atoms_function_value , bp::release_gil_policy() - , "Return the differential of this restraint with respect to\nthe symbol symbol\nThrow: SireCAS::unavailable_differential\n" ); + , "Return the atoms involved in the restraint" ); } - { //::SireMM::DihedralRestraint::force + { //::SireMM::DihedralRestraint::isNull - typedef void ( ::SireMM::DihedralRestraint::*force_function_type)( ::SireFF::MolForceTable &,double ) const; - force_function_type force_function_value( &::SireMM::DihedralRestraint::force ); + typedef bool ( ::SireMM::DihedralRestraint::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::DihedralRestraint::isNull ); 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") ) + "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::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 +106,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 + { //::SireMM::DihedralRestraint::phi0 - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point0_function_type)( ) const; - point0_function_type point0_function_value( &::SireMM::DihedralRestraint::point0 ); + typedef ::SireUnits::Dimension::Angle ( ::SireMM::DihedralRestraint::*phi0_function_type)( ) const; + phi0_function_type phi0_function_value( &::SireMM::DihedralRestraint::phi0 ); 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 - - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point2_function_type)( ) const; - point2_function_type point2_function_value( &::SireMM::DihedralRestraint::point2 ); - - 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 +142,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 - - typedef void ( ::SireMM::DihedralRestraint::*update_function_type)( ::SireMol::Molecules const & ) ; - update_function_type update_function_value( &::SireMM::DihedralRestraint::update ); - - DihedralRestraint_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::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 + { //::SireMM::DihedralRestraint::what - typedef bool ( ::SireMM::DihedralRestraint::*usesMoleculesIn_function_type)( ::SireMol::Molecules const & ) const; - usesMoleculesIn_function_type usesMoleculesIn_function_value( &::SireMM::DihedralRestraint::usesMoleculesIn ); + typedef char const * ( ::SireMM::DihedralRestraint::*what_function_type)( ) const; + what_function_type what_function_value( &::SireMM::DihedralRestraint::what ); DihedralRestraint_exposer.def( - "usesMoleculesIn" - , usesMoleculesIn_function_value - , ( bp::arg("molecules") ) + "what" + , what_function_value , 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/_MM.main.cpp b/wrapper/MM/_MM.main.cpp index 5a98d9755..b8559a1bd 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" @@ -546,11 +550,11 @@ BOOST_PYTHON_MODULE(_MM){ register_AngleParameterName_class(); - register_Restraint_class(); + register_AngleRestraint_class(); - register_Restraint3D_class(); + register_Restraints_class(); - register_AngleRestraint_class(); + register_AngleRestraints_class(); register_InternalSymbolsBase_class(); @@ -586,8 +590,6 @@ BOOST_PYTHON_MODULE(_MM){ register_BondRestraint_class(); - register_Restraints_class(); - register_BondRestraints_class(); register_BondSymbols_class(); @@ -646,6 +648,14 @@ BOOST_PYTHON_MODULE(_MM){ register_CLJParameterNames3D_class(); + register_CLJPotentialInterface_InterCLJPotential__class(); + + register_CLJPotentialInterface_InterSoftCLJPotential__class(); + + register_CLJPotentialInterface_IntraCLJPotential__class(); + + register_CLJPotentialInterface_IntraSoftCLJPotential__class(); + register_CLJProbe_class(); register_CLJRFFunction_class(); @@ -672,6 +682,10 @@ BOOST_PYTHON_MODULE(_MM){ register_CoulombNBPairs_class(); + register_CoulombPotentialInterface_InterCoulombPotential__class(); + + register_CoulombPotentialInterface_IntraCoulombPotential__class(); + register_CoulombProbe_class(); register_Dihedral_class(); @@ -682,8 +696,14 @@ BOOST_PYTHON_MODULE(_MM){ register_DihedralRestraint_class(); + register_DihedralRestraints_class(); + register_DihedralSymbols_class(); + register_Restraint_class(); + + register_Restraint3D_class(); + register_DistanceRestraint_class(); register_DoubleDistanceRestraint_class(); @@ -776,6 +796,10 @@ BOOST_PYTHON_MODULE(_MM){ register_LJPerturbation_class(); + register_LJPotentialInterface_InterLJPotential__class(); + + register_LJPotentialInterface_IntraLJPotential__class(); + register_LJProbe_class(); register_MMDetail_class(); @@ -820,6 +844,10 @@ BOOST_PYTHON_MODULE(_MM){ register_SoftCLJComponent_class(); + register_SoftCLJPotentialInterface_InterSoftCLJPotential__class(); + + register_SoftCLJPotentialInterface_IntraSoftCLJPotential__class(); + register_StretchBendComponent_class(); register_StretchBendSymbols_class(); From 386912d8b56e4613414c49542f9375a909c7b0e7 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Tue, 16 Jul 2024 16:05:36 +0100 Subject: [PATCH 04/15] Added intial python code for implementing alchemical dihedral restraints --- src/sire/mm/__init__.py | 5 ++ src/sire/restraints/__init__.py | 2 +- src/sire/restraints/_restraints.py | 100 +++++++++++++++++++++++++++++ 3 files changed, 106 insertions(+), 1 deletion(-) diff --git a/src/sire/mm/__init__.py b/src/sire/mm/__init__.py index 2255fd1ac..047d4a585 100644 --- a/src/sire/mm/__init__.py +++ b/src/sire/mm/__init__.py @@ -11,6 +11,8 @@ "Improper", "PositionRestraint", "PositionRestraints", + "DihedralRestraint", + "DihedralRestraints", "SelectorAngle", "SelectorBond", "SelectorDihedral", @@ -40,6 +42,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/restraints/__init__.py b/src/sire/restraints/__init__.py index 5f5af1465..f8283671f 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"] -from ._restraints import bond, boresch, distance, positional +from ._restraints import 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 d4c64774e..0d505bd1e 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -1,6 +1,7 @@ __all__ = [ "boresch", "bond", + "dihedral", "distance", "positional", ] @@ -374,6 +375,105 @@ def _check_stability_boresch_restraint(restraint_components, temperature=u("298 " values further from 0 or pi radians." ) +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 'ktheta' and equilibrium + bond length theta0. + + These restraints will be per atom-atom distance. If a list of k and/or r0 + values are passed, then different values could be used for + different atom-atom distances (assuming the same number as the number of + atom-atom distances). Otherwise, all atom-atom distances will use the + same parameters. + + If r0 is None, then the current atom-atom distance for + each atom-atom pair will be used as the equilibium value. + + If k is None, then a default value of 150 kcal mol-1 A-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 list of str or SireUnits::Dimension::GeneralUnit, optional + The force constants for the torsion restraints. + If None, this will default to 100 kcal mol-1 rad-2. + If a list, then this should be a. + Default is None. + + theta0 : list of 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 + ------- + 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]. + + Examples + -------- + Create a set of Dihedral restraints for the ligand in the system + 'system', specifying all of the force constants and equilibrium + values: + """ + from .. import u + from ..base import create_map + # from ..mm import DihedralRestraint + 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"but only {len(atoms)} atoms were provided." + ) + + from .. import measure + + if kphi is None: + kphi = u("100 kcal mol-1 rad-2") + + # TODO: Add support for multiple dihedral restraints + if phi0 is None: + # calculate all of the current angles + from .. import measure + # only support 1 dihedral restraint at the moment + phi0 = measure(atoms[0], atoms[1], atoms[2], atoms[3]) + + elif type(phi0) is list: + phi0 = [u(x) for x in phi0] + else: + phi0 = u(phi0) + + mols = mols.atoms() + + if name is None: + restraints = DihedralRestraints() + else: + restraints = DihedralRestraints(name=name) + + print(f"mols.find(atoms): {mols.find(atoms)}") + print(f"phi0: {phi0}") + print(f"kphi: {kphi}") + restraints.add(DihedralRestraint(mols.find(atoms), phi0, kphi)) + return restraints + def distance(mols, atoms0, atoms1, r0=None, k=None, name=None, map=None): """ From 69a7ee282c3a26e57f5160f8df0a17ebbc9a1a3a Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Tue, 16 Jul 2024 16:11:58 +0100 Subject: [PATCH 05/15] Added prototype code for adding alchemical dihedral restraints to OpenMM system --- .../SireOpenMM/sire_to_openmm_system.cpp | 95 ++++++++++++++++--- 1 file changed, 83 insertions(+), 12 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 03e2ad242..67720812b 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -8,37 +8,38 @@ #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/atomljs.h" #include "SireMM/bondrestraints.h" -#include "SireMM/positionalrestraints.h" #include "SireMM/boreschrestraints.h" +#include "SireMM/dihedralrestraint.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" @@ -393,6 +394,71 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, } } +/** Add all of the dihedral 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_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) + + // TODO: Add support for multiple dihedral restraints + + const auto energy_expression = QString( + "rho*k*delta*delta;" + "delta=(phi-phi0)") + .toStdString(); + + auto *restraintff = new OpenMM::CustomTorsionForce(energy_expression); + + restraintff->addPerTorsionParameter("rho"); + restraintff->addPerTorsionParameter("k"); + restraintff->addPerTorsionParameter("phi0"); + + 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 < 3; ++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; // kphi + parameters[2] = restraint.phi0().value(); // phi0 (already in radians) + + // 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 @@ -1625,7 +1691,12 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, CODELOC); // we now need to choose what to do based on the type of restraint... - if (prop.read().isA()) + if (prop.read().isA()) + { + _add_dihedral_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); From c29c4eb7ebac487ff2a9a39f6edcc206f1e2b20b Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Tue, 16 Jul 2024 17:04:30 +0100 Subject: [PATCH 06/15] Fix conversion of dihedral restraints to OpenMM system --- .../SireOpenMM/sire_to_openmm_system.cpp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 67720812b..10089d82d 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -418,15 +418,24 @@ void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, // TODO: Add support for multiple dihedral restraints const auto energy_expression = QString( - "rho*k*delta*delta;" - "delta=(phi-phi0)") + "rho*k*min(dtheta, 2*pi-dtheta)^2;" + "dtheta = abs(theta-theta0);" + "pi = 3.1415926535;") .toStdString(); auto *restraintff = new OpenMM::CustomTorsionForce(energy_expression); + // OLD CODE + // restraintff->addPerTorsionParameter("rho"); + // restraintff->addPerTorsionParameter("k"); + // restraintff->addPerTorsionParameter("phi0"); + + // it seems that OpenMM wants to call the torsion angle theta rather than phi + // we need to rename our parameters accordingly + // NEW CODE restraintff->addPerTorsionParameter("rho"); restraintff->addPerTorsionParameter("k"); - restraintff->addPerTorsionParameter("phi0"); + restraintff->addPerTorsionParameter("theta0"); restraintff->setUsesPeriodicBoundaryConditions(true); @@ -450,7 +459,7 @@ void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, std::vector parameters; parameters.resize(3); - parameters[0] = 1.0; // rho + parameters[0] = 1.0; // rho parameters[1] = restraint.kphi().value() * internal_to_ktheta; // kphi parameters[2] = restraint.phi0().value(); // phi0 (already in radians) From 4a0726d5c2a36f7e0954252cab75432bf8f6f195 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Wed, 17 Jul 2024 11:47:06 +0100 Subject: [PATCH 07/15] Added full end-to-end alchemical dihedral restraints implementation - bugs might still be lurking somewhere --- .../Convert/SireOpenMM/LambdaLever.pypp.cpp | 2 +- wrapper/Convert/SireOpenMM/lambdalever.cpp | 55 ++++++++++++++++++- .../SireOpenMM/sire_to_openmm_system.cpp | 2 +- 3 files changed, 56 insertions(+), 3 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp b/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp index 85786b4a2..152f5fc2f 100644 --- a/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp @@ -199,7 +199,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( diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index a70a0c18f..e7d4e62fd 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -28,8 +28,8 @@ #include "lambdalever.h" -#include "SireBase/propertymap.h" #include "SireBase/arrayproperty.hpp" +#include "SireBase/propertymap.h" #include "SireCAS/values.h" @@ -1825,6 +1825,53 @@ void _update_restraint_in_context(OpenMM::CustomCompoundBondForce *ff, double rh 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, @@ -1886,6 +1933,12 @@ void LambdaLever::updateRestraintInContext(OpenMM::Force &ff, double rho, dynamic_cast(&ff), rho, context); } + else if (ff_type == "CustomTorsionForce") + { + _update_restraint_in_context( + dynamic_cast(&ff), + rho, context); + } else if (ff_type == "CustomCompoundBondForce") { _update_restraint_in_context( diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 10089d82d..d660d5dc7 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -451,7 +451,7 @@ void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, std::vector particles; particles.resize(4); - for (int i = 0; i < 3; ++i) + for (int i = 0; i < 4; ++i) { particles[i] = restraint.atoms()[i]; } From 9eeb82435b366ef833bfcfb5fd4363430131ba11 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Fri, 19 Jul 2024 16:55:18 +0100 Subject: [PATCH 08/15] Added full end-to-end alchemical angle restraints implementation, and also renamed the core alchemical angle/dihedral restraint c++ files --- corelib/src/libs/SireMM/CMakeLists.txt | 8 +- corelib/src/libs/SireMM/anglerestraint.cpp | 995 ------------------ corelib/src/libs/SireMM/anglerestraint.h | 182 ---- ...alrestraint.cpp => dihedralrestraints.cpp} | 15 +- ...hedralrestraint.h => dihedralrestraints.h} | 9 +- src/sire/mm/__init__.py | 5 + src/sire/restraints/__init__.py | 4 +- src/sire/restraints/_restraints.py | 130 ++- wrapper/Convert/SireOpenMM/lambdalever.cpp | 53 + .../SireOpenMM/sire_to_openmm_system.cpp | 84 +- wrapper/MM/AngleRestraint.pypp.cpp | 16 +- wrapper/MM/AngleRestraints.pypp.cpp | 16 +- wrapper/MM/CMakeAutogenFile.txt | 2 + wrapper/MM/DihedralRestraint.pypp.cpp | 16 +- wrapper/MM/DihedralRestraints.pypp.cpp | 16 +- wrapper/MM/SireMM_registrars.cpp | 6 +- wrapper/MM/_MM.main.cpp | 28 +- wrapper/MM/active_headers.h | 4 +- 18 files changed, 262 insertions(+), 1327 deletions(-) delete mode 100644 corelib/src/libs/SireMM/anglerestraint.cpp delete mode 100644 corelib/src/libs/SireMM/anglerestraint.h rename corelib/src/libs/SireMM/{dihedralrestraint.cpp => dihedralrestraints.cpp} (97%) rename corelib/src/libs/SireMM/{dihedralrestraint.h => dihedralrestraints.h} (97%) 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 00e93dde5..000000000 --- a/corelib/src/libs/SireMM/anglerestraint.cpp +++ /dev/null @@ -1,995 +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" - -#include - -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.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) -{ - // Need to think here about validating the angle and force constant values - // if (atoms.count() != 3) - // { - // throw SireError::invalid_arg(QObject::tr( - // "Wrong number of inputs for an Angle restraint. You need to " - // "provide 3 atoms (%1).") - // .arg(atoms.count()), - // // .arg(theta0.count()) - // // .arg(ktheta.count()), - // CODELOC); - // } - - // make sure that we have 3 distinct atoms - QSet distinct; - distinct.reserve(3); - - for (const auto &atom : atoms) - { - if (atom >= 0) - distinct.insert(atom); - } - - // if (distinct.count() != 3) - // throw SireError::invalid_arg(QObject::tr( - // "There is something wrong with the atoms provided. " - // "They should all be unique and all greater than or equal to 0."), - // CODELOC); - - 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; -} - -// 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); -// } -// } - -// 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); -// } - -// Restraint::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 c8a4f228e..000000000 --- a/corelib/src/libs/SireMM/anglerestraint.h +++ /dev/null @@ -1,182 +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 "restraints.h" - -// #include "SireCAS/expression.h" -// #include "SireCAS/symbol.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/dihedralrestraints.cpp similarity index 97% rename from corelib/src/libs/SireMM/dihedralrestraint.cpp rename to corelib/src/libs/SireMM/dihedralrestraints.cpp index 68ed915f5..61a086e15 100644 --- a/corelib/src/libs/SireMM/dihedralrestraint.cpp +++ b/corelib/src/libs/SireMM/dihedralrestraints.cpp @@ -25,20 +25,10 @@ * \*********************************************/ -#include "dihedralrestraint.h" - -// #include "SireFF/forcetable.h" - -// #include "SireCAS/conditional.h" -// #include "SireCAS/power.h" -// #include "SireCAS/symbols.h" -// #include "SireCAS/values.h" +#include "dihedralrestraints.h" #include "SireID/index.h" - -// #include "SireUnits/angle.h" #include "SireUnits/units.h" - #include "SireStream/datastream.h" #include "SireStream/shareddatastream.h" @@ -47,11 +37,8 @@ #include 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; diff --git a/corelib/src/libs/SireMM/dihedralrestraint.h b/corelib/src/libs/SireMM/dihedralrestraints.h similarity index 97% rename from corelib/src/libs/SireMM/dihedralrestraint.h rename to corelib/src/libs/SireMM/dihedralrestraints.h index b340a3411..d3f6b8a35 100644 --- a/corelib/src/libs/SireMM/dihedralrestraint.h +++ b/corelib/src/libs/SireMM/dihedralrestraints.h @@ -25,16 +25,11 @@ * \*********************************************/ -#ifndef SIREMM_DIHEDRALRESTRAINT_H -#define SIREMM_DIHEDRALRESTRAINT_H - -// #include "SireFF/point.h" +#ifndef SIREMM_DIHEDRALRESTRAINTS_H +#define SIREMM_DIHEDRALRESTRAINTS_H #include "restraints.h" -// #include "SireCAS/expression.h" -// #include "SireCAS/symbol.h" - #include "SireUnits/dimensions.h" #include "SireUnits/generalunit.h" diff --git a/src/sire/mm/__init__.py b/src/sire/mm/__init__.py index 047d4a585..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", @@ -31,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 diff --git a/src/sire/restraints/__init__.py b/src/sire/restraints/__init__.py index f8283671f..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, dihedral, 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 0d505bd1e..f1cb8377f 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -19,6 +19,95 @@ 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"but only {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. You can setup one restraint at a time however." + ) + + # TODO: Add support for multiple angle restraints + if theta0 is None: + # calculate all of the current angles + 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. You can setup one restraint at a time however." + ) + 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, @@ -379,19 +468,13 @@ 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 'ktheta' and equilibrium - bond length theta0. + passed values of the force constant 'kphi' and equilibrium + torsion angle phi0. - These restraints will be per atom-atom distance. If a list of k and/or r0 - values are passed, then different values could be used for - different atom-atom distances (assuming the same number as the number of - atom-atom distances). Otherwise, all atom-atom distances will use the - same parameters. - - If r0 is None, then the current atom-atom distance for - each atom-atom pair will be used as the equilibium value. + If phi0 is None, then the current torsional angle for + the provided atoms will be used as the equilibium value. - If k is None, then a default value of 150 kcal mol-1 A-2 will be used + If kphi is None, then a default value of 100 kcal mol-1 rad-2 will be used Parameters ---------- @@ -401,14 +484,13 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): atoms : SireMol::Selector The atoms to restrain. - kphi : str or SireUnits::Dimension::GeneralUnit or list of str or SireUnits::Dimension::GeneralUnit, optional + 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. - If a list, then this should be a. Default is None. - theta0 : list of str or SireUnits::Dimension::GeneralUnit, optional - The equilibrium angles for the angle restraints. If None, these + 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. @@ -418,12 +500,6 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): A container of Dihedral restraints, where the first restraint is the DihedralRestraint created. The Dihedral restraint created can be extracted with DihedralRestraints[0]. - - Examples - -------- - Create a set of Dihedral restraints for the ligand in the system - 'system', specifying all of the force constants and equilibrium - values: """ from .. import u from ..base import create_map @@ -449,6 +525,11 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): 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. You can setup one restraint at a time however." + ) + # TODO: Add support for multiple dihedral restraints if phi0 is None: # calculate all of the current angles @@ -457,7 +538,9 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): phi0 = measure(atoms[0], atoms[1], atoms[2], atoms[3]) elif type(phi0) is list: - phi0 = [u(x) for x in phi0] + raise NotImplementedError( + "Setup of multiple dihedral restraints simultaneously is not currently supported. You can setup one restraint at a time however." + ) else: phi0 = u(phi0) @@ -468,9 +551,6 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): else: restraints = DihedralRestraints(name=name) - print(f"mols.find(atoms): {mols.find(atoms)}") - print(f"phi0: {phi0}") - print(f"kphi: {kphi}") restraints.add(DihedralRestraint(mols.find(atoms), phi0, kphi)) return restraints diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index e7d4e62fd..dde2e2eb0 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1825,6 +1825,53 @@ 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, @@ -1933,6 +1980,12 @@ void LambdaLever::updateRestraintInContext(OpenMM::Force &ff, double rho, dynamic_cast(&ff), rho, context); } + else if (ff_type == "CustomAngleForce") + { + _update_restraint_in_context( + dynamic_cast(&ff), + rho, context); + } else if (ff_type == "CustomTorsionForce") { _update_restraint_in_context( diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index d660d5dc7..551c1f0ba 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -21,10 +21,11 @@ #include "SireMol/moleditor.h" #include "SireMM/amberparams.h" +#include "SireMM/anglerestraints.h" #include "SireMM/atomljs.h" #include "SireMM/bondrestraints.h" #include "SireMM/boreschrestraints.h" -#include "SireMM/dihedralrestraint.h" +#include "SireMM/dihedralrestraints.h" #include "SireMM/positionalrestraints.h" #include "SireMM/selectorbond.h" @@ -394,10 +395,68 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, } } -/** Add all of the dihedral restraints from 'restraints' to the passed +/** 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->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) @@ -415,24 +474,16 @@ void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, // e_torsion = k_phi(dphi)^2 where // dphi = abs(phi - phi0) - // TODO: Add support for multiple dihedral restraints - const auto energy_expression = QString( - "rho*k*min(dtheta, 2*pi-dtheta)^2;" + "rho*k*min(dtheta, two_pi-dtheta)^2;" "dtheta = abs(theta-theta0);" - "pi = 3.1415926535;") + "two_pi=6.283185307179586;") .toStdString(); auto *restraintff = new OpenMM::CustomTorsionForce(energy_expression); - // OLD CODE - // restraintff->addPerTorsionParameter("rho"); - // restraintff->addPerTorsionParameter("k"); - // restraintff->addPerTorsionParameter("phi0"); - // it seems that OpenMM wants to call the torsion angle theta rather than phi // we need to rename our parameters accordingly - // NEW CODE restraintff->addPerTorsionParameter("rho"); restraintff->addPerTorsionParameter("k"); restraintff->addPerTorsionParameter("theta0"); @@ -460,8 +511,8 @@ void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, parameters.resize(3); parameters[0] = 1.0; // rho - parameters[1] = restraint.kphi().value() * internal_to_ktheta; // kphi - parameters[2] = restraint.phi0().value(); // phi0 (already in radians) + 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); @@ -1705,6 +1756,11 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, _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(), diff --git a/wrapper/MM/AngleRestraint.pypp.cpp b/wrapper/MM/AngleRestraint.pypp.cpp index 33954da56..b90320621 100644 --- a/wrapper/MM/AngleRestraint.pypp.cpp +++ b/wrapper/MM/AngleRestraint.pypp.cpp @@ -7,33 +7,21 @@ 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); } diff --git a/wrapper/MM/AngleRestraints.pypp.cpp b/wrapper/MM/AngleRestraints.pypp.cpp index fdf557c21..07edf38ca 100644 --- a/wrapper/MM/AngleRestraints.pypp.cpp +++ b/wrapper/MM/AngleRestraints.pypp.cpp @@ -8,33 +8,21 @@ 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::AngleRestraints __copy__(const SireMM::AngleRestraints &other){ return SireMM::AngleRestraints(other); } 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 accddf41c..a5d04277e 100644 --- a/wrapper/MM/DihedralRestraint.pypp.cpp +++ b/wrapper/MM/DihedralRestraint.pypp.cpp @@ -7,33 +7,21 @@ 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); } diff --git a/wrapper/MM/DihedralRestraints.pypp.cpp b/wrapper/MM/DihedralRestraints.pypp.cpp index 21f56e7e0..44d49c86a 100644 --- a/wrapper/MM/DihedralRestraints.pypp.cpp +++ b/wrapper/MM/DihedralRestraints.pypp.cpp @@ -8,33 +8,21 @@ 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::DihedralRestraints __copy__(const SireMM::DihedralRestraints &other){ return SireMM::DihedralRestraints(other); } 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 b8559a1bd..addf5ef01 100644 --- a/wrapper/MM/_MM.main.cpp +++ b/wrapper/MM/_MM.main.cpp @@ -550,6 +550,10 @@ BOOST_PYTHON_MODULE(_MM){ register_AngleParameterName_class(); + register_Restraint_class(); + + register_Restraint3D_class(); + register_AngleRestraint_class(); register_Restraints_class(); @@ -648,14 +652,6 @@ BOOST_PYTHON_MODULE(_MM){ register_CLJParameterNames3D_class(); - register_CLJPotentialInterface_InterCLJPotential__class(); - - register_CLJPotentialInterface_InterSoftCLJPotential__class(); - - register_CLJPotentialInterface_IntraCLJPotential__class(); - - register_CLJPotentialInterface_IntraSoftCLJPotential__class(); - register_CLJProbe_class(); register_CLJRFFunction_class(); @@ -682,10 +678,6 @@ BOOST_PYTHON_MODULE(_MM){ register_CoulombNBPairs_class(); - register_CoulombPotentialInterface_InterCoulombPotential__class(); - - register_CoulombPotentialInterface_IntraCoulombPotential__class(); - register_CoulombProbe_class(); register_Dihedral_class(); @@ -700,10 +692,6 @@ BOOST_PYTHON_MODULE(_MM){ register_DihedralSymbols_class(); - register_Restraint_class(); - - register_Restraint3D_class(); - register_DistanceRestraint_class(); register_DoubleDistanceRestraint_class(); @@ -796,10 +784,6 @@ BOOST_PYTHON_MODULE(_MM){ register_LJPerturbation_class(); - register_LJPotentialInterface_InterLJPotential__class(); - - register_LJPotentialInterface_IntraLJPotential__class(); - register_LJProbe_class(); register_MMDetail_class(); @@ -844,10 +828,6 @@ BOOST_PYTHON_MODULE(_MM){ register_SoftCLJComponent_class(); - register_SoftCLJPotentialInterface_InterSoftCLJPotential__class(); - - register_SoftCLJPotentialInterface_IntraSoftCLJPotential__class(); - register_StretchBendComponent_class(); register_StretchBendSymbols_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" From a500bc9eb5a04447c7e631685a6c26a1e70a2d98 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Fri, 19 Jul 2024 18:13:32 +0100 Subject: [PATCH 09/15] Add the missing corelib files for angle restraints --- corelib/src/libs/SireMM/anglerestraints.cpp | 496 ++++++++++++++++++++ corelib/src/libs/SireMM/anglerestraints.h | 177 +++++++ 2 files changed, 673 insertions(+) create mode 100644 corelib/src/libs/SireMM/anglerestraints.cpp create mode 100644 corelib/src/libs/SireMM/anglerestraints.h diff --git a/corelib/src/libs/SireMM/anglerestraints.cpp b/corelib/src/libs/SireMM/anglerestraints.cpp new file mode 100644 index 000000000..42cec09df --- /dev/null +++ b/corelib/src/libs/SireMM/anglerestraints.cpp @@ -0,0 +1,496 @@ +/********************************************\ + * + * 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 "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) +{ + // Need to think here about validating the angle and force constant values + // if (atoms.count() != 3) + // { + // throw SireError::invalid_arg(QObject::tr( + // "Wrong number of inputs for an Angle restraint. You need to " + // "provide 3 atoms (%1).") + // .arg(atoms.count()), + // // .arg(theta0.count()) + // // .arg(ktheta.count()), + // CODELOC); + // } + + // make sure that we have 3 distinct atoms + QSet distinct; + distinct.reserve(3); + + for (const auto &atom : atoms) + { + if (atom >= 0) + distinct.insert(atom); + } + + // if (distinct.count() != 3) + // throw SireError::invalid_arg(QObject::tr( + // "There is something wrong with the atoms provided. " + // "They should all be unique and all greater than or equal to 0."), + // CODELOC); + + 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..393d30849 --- /dev/null +++ b/corelib/src/libs/SireMM/anglerestraints.h @@ -0,0 +1,177 @@ +/********************************************\ + * + * 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_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 From 691223a0f2495f709a1e0296b8a486009d4704d7 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Fri, 10 Jan 2025 15:56:31 +0000 Subject: [PATCH 10/15] Updated restraints tutorial to include angle and dihedral restraints, also slightly tidied up python API for these restraints --- doc/source/tutorial/part06/03_restraints.rst | 68 ++++++++++++++++++++ src/sire/restraints/_restraints.py | 19 ++---- 2 files changed, 75 insertions(+), 12 deletions(-) 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/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 786163a73..83ea94717 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -1,4 +1,5 @@ __all__ = [ + "angle", "boresch", "bond", "dihedral", @@ -71,7 +72,7 @@ def angle(mols, atoms, theta0=None, ktheta=None, name=None, map=None): if len(atoms) != 3: raise ValueError( "You need to provide 3 atoms to create an angle restraint" - f"but only {len(atoms)} atoms were provided." + f"whereas {len(atoms)} atoms were provided." ) from .. import measure @@ -81,18 +82,16 @@ def angle(mols, atoms, theta0=None, ktheta=None, name=None, map=None): elif type(ktheta) is list: raise NotImplementedError( - "Setup of multiple angle restraints simultaneously is not currently supported. You can setup one restraint at a time however." + "Setup of multiple angle restraints simultaneously is not currently supported. Please set up each restraint individually and then combine them into multiple restraints." ) - # TODO: Add support for multiple angle restraints if theta0 is None: - # calculate all of the current angles 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. You can setup one restraint at a time however." + "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) @@ -505,7 +504,6 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): """ from .. import u from ..base import create_map - # from ..mm import DihedralRestraint from ..mm import DihedralRestraint, DihedralRestraints map = create_map(map) @@ -519,7 +517,7 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): if len(atoms) != 4: raise ValueError( "You need to provide 4 atoms to create a dihedral restraint" - f"but only {len(atoms)} atoms were provided." + f"whereas {len(atoms)} atoms were provided." ) from .. import measure @@ -529,19 +527,16 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): elif type(kphi) is list: raise NotImplementedError( - "Setup of multiple dihedral restraints simultaneously is not currently supported. You can setup one restraint at a time however." + "Setup of multiple dihedral restraints simultaneously is not currently supported. Please set up each restraint individually and then combine them into multiple restraints." ) - # TODO: Add support for multiple dihedral restraints if phi0 is None: - # calculate all of the current angles from .. import measure - # only support 1 dihedral restraint at the moment 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. You can setup one restraint at a time however." + "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) From ab2188dbb6e5ff51b3ba5fe10956996e126c6fbb Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Fri, 10 Jan 2025 16:29:09 +0000 Subject: [PATCH 11/15] Add basic unit tests for the new restraints --- .../test_angle_dihedral_restraints.py | 88 +++++++++++++++++++ 1 file changed, 88 insertions(+) create mode 100644 tests/restraints/test_angle_dihedral_restraints.py diff --git a/tests/restraints/test_angle_dihedral_restraints.py b/tests/restraints/test_angle_dihedral_restraints.py new file mode 100644 index 000000000..9bacd8069 --- /dev/null +++ b/tests/restraints/test_angle_dihedral_restraints.py @@ -0,0 +1,88 @@ +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 \ No newline at end of file From 35725d8d89e8de8ba2aa0f41d9ebcc5a6fee1f9a Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Fri, 10 Jan 2025 16:37:11 +0000 Subject: [PATCH 12/15] Updated changelog and contributor pages for the PR --- doc/source/changelog.rst | 2 ++ doc/source/contributors.rst | 1 + 2 files changed, 3 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 84df9d070..39b255b85 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -19,6 +19,8 @@ organisation on `GitHub `__. * Fixed update of triclinic box vectors in ``SOMD`` following ``OpenMM`` bug fix. +* Added support for angle and dihedral restraints which can be used in alchemical and standard simulations. + `2024.3.1 `__ - December 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 `__ From 65fbd034d2addf414690b1c1d7b86883d7fc9f21 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Tue, 14 Jan 2025 16:31:27 +0000 Subject: [PATCH 13/15] Removed whitespace that was added and tidied up python restraints files [ci skip] --- src/sire/restraints/_restraints.py | 19 ++++++++----- .../test_angle_dihedral_restraints.py | 27 +++++++++++++++---- wrapper/Convert/SireOpenMM/lambdalever.cpp | 20 +++++++------- 3 files changed, 45 insertions(+), 21 deletions(-) diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 83ea94717..6e4d5d050 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -49,7 +49,7 @@ def angle(mols, atoms, theta0=None, ktheta=None, name=None, map=None): 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 @@ -82,16 +82,19 @@ def angle(mols, atoms, theta0=None, ktheta=None, name=None, map=None): 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." + "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." + "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) @@ -465,6 +468,7 @@ def _check_stability_boresch_restraint(restraint_components, temperature=u("298 " values further from 0 or pi radians." ) + 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' @@ -494,7 +498,7 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): 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 @@ -527,16 +531,19 @@ def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): 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." + "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." + "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) diff --git a/tests/restraints/test_angle_dihedral_restraints.py b/tests/restraints/test_angle_dihedral_restraints.py index 9bacd8069..f6096ec21 100644 --- a/tests/restraints/test_angle_dihedral_restraints.py +++ b/tests/restraints/test_angle_dihedral_restraints.py @@ -1,6 +1,7 @@ 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() @@ -9,23 +10,30 @@ def test_default_angle_restraints_setup(ala_mols): 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") + 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") + 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() @@ -41,29 +49,37 @@ def test_default_dihedral_restraints_setup(ala_mols): 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") + 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") + 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() @@ -74,6 +90,7 @@ def test_multiple_angle_restraints(ala_mols): 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() @@ -85,4 +102,4 @@ def test_multiple_dihedral_restraints(ala_mols): restraint2 = sr.restraints.dihedral(mols=mols, atoms=dih2.atoms()) restraints.add(restraint1) restraints.add(restraint2) - assert restraints.num_restraints() == 3 \ No newline at end of file + assert restraints.num_restraints() == 3 diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 32ffcf599..ad6ccb580 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1473,8 +1473,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, if (nbidx < 0) throw SireError::program_bug(QObject::tr( - "Unset NB14 index for a ghost atom?"), - CODELOC); + "Unset NB14 index for a ghost atom?"), + CODELOC); coul_14_scale = morphed_ghost14_charge_scale[j]; lj_14_scale = morphed_ghost14_lj_scale[j]; @@ -1602,11 +1602,11 @@ double LambdaLever::setLambda(OpenMM::Context &context, double length, k; bondff->getBondParameters(index, particle1, particle2, - length, k); + length, k); bondff->setBondParameters(index, particle1, particle2, - morphed_bond_length[j], - morphed_bond_k[j]); + morphed_bond_length[j], + morphed_bond_k[j]); } } @@ -1651,13 +1651,13 @@ double LambdaLever::setLambda(OpenMM::Context &context, double size, k; angff->getAngleParameters(index, - particle1, particle2, particle3, - size, k); + particle1, particle2, particle3, + size, k); angff->setAngleParameters(index, - particle1, particle2, particle3, - morphed_angle_size[j], - morphed_angle_k[j]); + particle1, particle2, particle3, + morphed_angle_size[j], + morphed_angle_k[j]); } } From 945674cf361b4f933a56c3b006715b3350feb8a3 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Tue, 14 Jan 2025 17:38:31 +0000 Subject: [PATCH 14/15] Added force names to custom OpenMM angle and torsion forces --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 4 ++-- wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index ad6ccb580..a05af04dc 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -2012,13 +2012,13 @@ void LambdaLever::updateRestraintInContext(OpenMM::Force &ff, double rho, dynamic_cast(&ff), rho, context); } - else if (ff_type == "CustomAngleForce") + else if (ff_type == "AngleRestraintForce") { _update_restraint_in_context( dynamic_cast(&ff), rho, context); } - else if (ff_type == "CustomTorsionForce") + else if (ff_type == "TorsionRestraintForce") { _update_restraint_in_context( dynamic_cast(&ff), diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 3bee9af57..4c5b928d9 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -427,6 +427,7 @@ void _add_angle_restraints(const SireMM::AngleRestraints &restraints, auto *restraintff = new OpenMM::CustomAngleForce(energy_expression); + restraintff->setName("AngleRestraintForce"); restraintff->addPerAngleParameter("rho"); restraintff->addPerAngleParameter("k"); restraintff->addPerAngleParameter("theta0"); @@ -489,6 +490,7 @@ void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, // 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"); From 0e56c8eb18400da6797fc5169793c30f43744fe4 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Tue, 21 Jan 2025 15:50:59 +0000 Subject: [PATCH 15/15] Removed redundant comments for restraints code [ci skip] --- corelib/src/libs/SireMM/anglerestraints.cpp | 23 ++---------- corelib/src/libs/SireMM/anglerestraints.h | 2 +- .../src/libs/SireMM/dihedralrestraints.cpp | 35 +++++-------------- corelib/src/libs/SireMM/dihedralrestraints.h | 4 +-- 4 files changed, 14 insertions(+), 50 deletions(-) diff --git a/corelib/src/libs/SireMM/anglerestraints.cpp b/corelib/src/libs/SireMM/anglerestraints.cpp index 42cec09df..5eaf82fc1 100644 --- a/corelib/src/libs/SireMM/anglerestraints.cpp +++ b/corelib/src/libs/SireMM/anglerestraints.cpp @@ -2,7 +2,7 @@ * * Sire - Molecular Simulation Framework * - * Copyright (C) 2009 Christopher Woods + * 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 @@ -98,19 +98,8 @@ AngleRestraint::AngleRestraint(const QList &atoms, : ConcreteProperty(), _theta0(theta0), _ktheta(ktheta) { - // Need to think here about validating the angle and force constant values - // if (atoms.count() != 3) - // { - // throw SireError::invalid_arg(QObject::tr( - // "Wrong number of inputs for an Angle restraint. You need to " - // "provide 3 atoms (%1).") - // .arg(atoms.count()), - // // .arg(theta0.count()) - // // .arg(ktheta.count()), - // CODELOC); - // } - - // make sure that we have 3 distinct atoms + + // Make sure that we have 3 distinct atoms QSet distinct; distinct.reserve(3); @@ -120,12 +109,6 @@ AngleRestraint::AngleRestraint(const QList &atoms, distinct.insert(atom); } - // if (distinct.count() != 3) - // throw SireError::invalid_arg(QObject::tr( - // "There is something wrong with the atoms provided. " - // "They should all be unique and all greater than or equal to 0."), - // CODELOC); - atms = atoms.toVector(); } diff --git a/corelib/src/libs/SireMM/anglerestraints.h b/corelib/src/libs/SireMM/anglerestraints.h index 393d30849..b34de4822 100644 --- a/corelib/src/libs/SireMM/anglerestraints.h +++ b/corelib/src/libs/SireMM/anglerestraints.h @@ -2,7 +2,7 @@ * * Sire - Molecular Simulation Framework * - * Copyright (C) 2009 Christopher Woods + * 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 diff --git a/corelib/src/libs/SireMM/dihedralrestraints.cpp b/corelib/src/libs/SireMM/dihedralrestraints.cpp index 61a086e15..70054b34c 100644 --- a/corelib/src/libs/SireMM/dihedralrestraints.cpp +++ b/corelib/src/libs/SireMM/dihedralrestraints.cpp @@ -2,7 +2,7 @@ * * Sire - Molecular Simulation Framework * - * Copyright (C) 2009 Christopher Woods + * 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 @@ -28,9 +28,9 @@ #include "dihedralrestraints.h" #include "SireID/index.h" -#include "SireUnits/units.h" #include "SireStream/datastream.h" #include "SireStream/shareddatastream.h" +#include "SireUnits/units.h" #include "SireCAS/errors.h" @@ -80,7 +80,6 @@ QDataStream &operator>>(QDataStream &ds, DihedralRestraint &dihrest) return ds; } - /** Null constructor */ DihedralRestraint::DihedralRestraint() : ConcreteProperty(), @@ -88,29 +87,17 @@ DihedralRestraint::DihedralRestraint() { } - /** 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) + const SireUnits::Dimension::Angle &phi0, + const SireUnits::Dimension::HarmonicAngleConstant &kphi) : ConcreteProperty(), _phi0(phi0), _kphi(kphi) { - // Need to think here about validating the angle and force constant values - // if (atoms.count() != 3) - // { - // throw SireError::invalid_arg(QObject::tr( - // "Wrong number of inputs for an Angle restraint. You need to " - // "provide 3 atoms (%1).") - // .arg(atoms.count()), - // // .arg(phi0.count()) - // // .arg(kphi.count()), - // CODELOC); - // } - - // make sure that we have 3 distinct atoms + + // Make sure that we have 4 distinct atoms QSet distinct; distinct.reserve(4); @@ -120,12 +107,6 @@ DihedralRestraint::DihedralRestraint(const QList &atoms, distinct.insert(atom); } - // if (distinct.count() != 3) - // throw SireError::invalid_arg(QObject::tr( - // "There is something wrong with the atoms provided. " - // "They should all be unique and all greater than or equal to 0."), - // CODELOC); - atms = atoms.toVector(); } @@ -301,7 +282,7 @@ DihedralRestraints::DihedralRestraints(const QList &restraint } DihedralRestraints::DihedralRestraints(const QString &name, - const DihedralRestraint &restraint) + const DihedralRestraint &restraint) : ConcreteProperty(name) { if (not restraint.isNull()) @@ -309,7 +290,7 @@ DihedralRestraints::DihedralRestraints(const QString &name, } DihedralRestraints::DihedralRestraints(const QString &name, - const QList &restraints) + const QList &restraints) : ConcreteProperty(name) { for (const auto &restraint : restraints) diff --git a/corelib/src/libs/SireMM/dihedralrestraints.h b/corelib/src/libs/SireMM/dihedralrestraints.h index d3f6b8a35..f40b03e0c 100644 --- a/corelib/src/libs/SireMM/dihedralrestraints.h +++ b/corelib/src/libs/SireMM/dihedralrestraints.h @@ -2,7 +2,7 @@ * * Sire - Molecular Simulation Framework * - * Copyright (C) 2009 Christopher Woods + * 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 @@ -50,7 +50,7 @@ SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::DihedralRestraints namespace SireMM { - /** This class represents a single angle restraint between any three + /** This class represents a single torsion restraint between any four * atoms in a system * @author Christopher Woods */