diff --git a/corelib/src/libs/SireMM/CMakeLists.txt b/corelib/src/libs/SireMM/CMakeLists.txt index f228580aa..f0e650fba 100644 --- a/corelib/src/libs/SireMM/CMakeLists.txt +++ b/corelib/src/libs/SireMM/CMakeLists.txt @@ -16,7 +16,7 @@ include_directories(${TBB_INCLUDE_DIR}) # Define the headers in SireMM set ( SIREMM_HEADERS amberparams.h - anglerestraint.h + anglerestraints.h angle.h atomfunctions.h atomljs.h @@ -43,7 +43,7 @@ set ( SIREMM_HEADERS cljworkspace.h coulombpotential.h dihedral.h - dihedralrestraint.h + dihedralrestraints.h distancerestraint.h errors.h excludedpairs.h @@ -110,7 +110,7 @@ set ( SIREMM_SOURCES amberparams.cpp angle.cpp - anglerestraint.cpp + anglerestraints.cpp atomfunctions.cpp atomljs.cpp bond.cpp @@ -135,7 +135,7 @@ set ( SIREMM_SOURCES cljworkspace.cpp coulombpotential.cpp dihedral.cpp - dihedralrestraint.cpp + dihedralrestraints.cpp distancerestraint.cpp errors.cpp excludedpairs.cpp diff --git a/corelib/src/libs/SireMM/anglerestraint.cpp b/corelib/src/libs/SireMM/anglerestraint.cpp deleted file mode 100644 index a765049c0..000000000 --- a/corelib/src/libs/SireMM/anglerestraint.cpp +++ /dev/null @@ -1,555 +0,0 @@ -/********************************************\ - * - * Sire - Molecular Simulation Framework - * - * Copyright (C) 2009 Christopher Woods - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * - * For full details of the license please see the COPYING file - * that should have come with this distribution. - * - * You can contact the authors at https://sire.openbiosim.org - * -\*********************************************/ - -#include "anglerestraint.h" - -#include "SireFF/forcetable.h" - -#include "SireCAS/conditional.h" -#include "SireCAS/power.h" -#include "SireCAS/symbols.h" -#include "SireCAS/values.h" - -#include "SireID/index.h" - -#include "SireUnits/angle.h" -#include "SireUnits/units.h" - -#include "SireStream/datastream.h" -#include "SireStream/shareddatastream.h" - -#include "SireCAS/errors.h" - -using namespace SireMM; -using namespace SireMol; -using namespace SireFF; -using namespace SireID; -using namespace SireBase; -using namespace SireCAS; -using namespace SireMaths; -using namespace SireStream; -using namespace SireUnits; -using namespace SireUnits::Dimension; - -//////////// -//////////// Implementation of AngleRestraint -//////////// - -static const RegisterMetaType r_angrest; - -/** Serialise to a binary datastream */ -QDataStream &operator<<(QDataStream &ds, const AngleRestraint &angrest) -{ - writeHeader(ds, r_angrest, 1); - - SharedDataStream sds(ds); - - sds << angrest.p[0] << angrest.p[1] << angrest.p[2] << angrest.force_expression - << static_cast(angrest); - - return ds; -} - -/** Extract from a binary datastream */ -QDataStream &operator>>(QDataStream &ds, AngleRestraint &angrest) -{ - VersionID v = readHeader(ds, r_angrest); - - if (v == 1) - { - SharedDataStream sds(ds); - - sds >> angrest.p[0] >> angrest.p[1] >> angrest.p[2] >> angrest.force_expression >> - static_cast(angrest); - - angrest.intra_molecule_points = Point::areIntraMoleculePoints(angrest.p[0], angrest.p[1]) and - Point::areIntraMoleculePoints(angrest.p[0], angrest.p[2]); - } - else - throw version_error(v, "1", r_angrest, CODELOC); - - return ds; -} - -Q_GLOBAL_STATIC_WITH_ARGS(Symbol, getThetaSymbol, ("theta")); - -/** Return the symbol that represents the angle between the points (theta) */ -const Symbol &AngleRestraint::theta() -{ - return *(getThetaSymbol()); -} - -/** Constructor */ -AngleRestraint::AngleRestraint() : ConcreteProperty() -{ -} - -void AngleRestraint::calculateTheta() -{ - if (this->restraintFunction().isFunction(theta())) - { - SireUnits::Dimension::Angle angle; - - if (intra_molecule_points) - // we don't use the space when calculating intra-molecular angles - angle = Vector::angle(p[0].read().point(), p[1].read().point(), p[2].read().point()); - else - angle = this->space().calcAngle(p[0].read().point(), p[1].read().point(), p[2].read().point()); - - ExpressionRestraint3D::_pvt_setValue(theta(), angle); - } -} - -/** Construct a restraint that acts on the angle within the - three points 'point0', 'point1' and 'point2' (theta == a(012)), - restraining the angle within these points using the expression - 'restraint' */ -AngleRestraint::AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Expression &restraint) - : ConcreteProperty(restraint) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - - force_expression = this->restraintFunction().differentiate(theta()); - - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]); - - this->calculateTheta(); -} - -/** Construct a restraint that acts on the angle within the - three points 'point0', 'point1' and 'point2' (theta == a(012)), - restraining the angle within these points using the expression - 'restraint' */ -AngleRestraint::AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Expression &restraint, const Values &values) - : ConcreteProperty(restraint, values) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - - force_expression = this->restraintFunction().differentiate(theta()); - - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]); - - this->calculateTheta(); -} - -/** Internal constructor used to construct a restraint using the specified - points, energy expression and force expression */ -AngleRestraint::AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Expression &nrg_restraint, const Expression &force_restraint) - : ConcreteProperty(nrg_restraint), force_expression(force_restraint) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - - if (force_expression.isConstant()) - { - force_expression = force_expression.evaluate(Values()); - } - else - { - if (not this->restraintFunction().symbols().contains(force_expression.symbols())) - throw SireError::incompatible_error( - QObject::tr("You cannot use a force function which uses more symbols " - "(%1) than the energy function (%2).") - .arg(Sire::toString(force_expression.symbols()), Sire::toString(restraintFunction().symbols())), - CODELOC); - } - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]); - - this->calculateTheta(); -} - -/** Copy constructor */ -AngleRestraint::AngleRestraint(const AngleRestraint &other) - : ConcreteProperty(other), force_expression(other.force_expression), - intra_molecule_points(other.intra_molecule_points) -{ - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } -} - -/** Destructor */ -AngleRestraint::~AngleRestraint() -{ -} - -/** Copy assignment operator */ -AngleRestraint &AngleRestraint::operator=(const AngleRestraint &other) -{ - if (this != &other) - { - ExpressionRestraint3D::operator=(other); - - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } - - force_expression = other.force_expression; - intra_molecule_points = other.intra_molecule_points; - } - - return *this; -} - -/** Comparison operator */ -bool AngleRestraint::operator==(const AngleRestraint &other) const -{ - return this == &other or (ExpressionRestraint3D::operator==(other) and p[0] == other.p[0] and p[1] == other.p[1] and - p[2] == other.p[2] and force_expression == other.force_expression); -} - -/** Comparison operator */ -bool AngleRestraint::operator!=(const AngleRestraint &other) const -{ - return not AngleRestraint::operator==(other); -} - -const char *AngleRestraint::typeName() -{ - return QMetaType::typeName(qMetaTypeId()); -} - -/** This restraint involves three points */ -int AngleRestraint::nPoints() const -{ - return 3; -} - -/** Return the ith point */ -const Point &AngleRestraint::point(int i) const -{ - i = Index(i).map(this->nPoints()); - - return p[i].read(); -} - -/** Return the first point */ -const Point &AngleRestraint::point0() const -{ - return p[0].read(); -} - -/** Return the second point */ -const Point &AngleRestraint::point1() const -{ - return p[1].read(); -} - -/** Return the third point */ -const Point &AngleRestraint::point2() const -{ - return p[2].read(); -} - -/** Return the built-in symbols of this restraint */ -Symbols AngleRestraint::builtinSymbols() const -{ - if (this->restraintFunction().isFunction(theta())) - return theta(); - else - return Symbols(); -} - -/** Return the values of the built-in symbols of this restraint */ -Values AngleRestraint::builtinValues() const -{ - if (this->restraintFunction().isFunction(theta())) - return theta() == this->values()[theta()]; - else - return Values(); -} - -/** Return the differential of this restraint with respect to - the symbol 'symbol' - - \throw SireCAS::unavailable_differential -*/ -RestraintPtr AngleRestraint::differentiate(const Symbol &symbol) const -{ - if (this->restraintFunction().isFunction(symbol)) - return AngleRestraint(p[0], p[1], p[2], restraintFunction().differentiate(symbol), this->values()); - else - return NullRestraint(); -} - -/** Set the space used to evaluate the energy of this restraint - - \throw SireVol::incompatible_space -*/ -void AngleRestraint::setSpace(const Space &new_space) -{ - if (not this->space().equals(new_space)) - { - AngleRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().setSpace(new_space); - } - - Restraint3D::setSpace(new_space); - - this->calculateTheta(); - } - catch (...) - { - AngleRestraint::operator=(old_state); - throw; - } - } -} - -/** Return the function used to calculate the restraint force */ -const Expression &AngleRestraint::differentialRestraintFunction() const -{ - return force_expression; -} - -/** Calculate the force acting on the molecule in the forcetable 'forcetable' - caused by this restraint, and add it on to the forcetable scaled by - 'scale_force' */ -void AngleRestraint::force(MolForceTable &forcetable, double scale_force) const -{ - bool in_p0 = p[0].read().contains(forcetable.molNum()); - bool in_p1 = p[1].read().contains(forcetable.molNum()); - bool in_p2 = p[2].read().contains(forcetable.molNum()); - - if (not(in_p0 or in_p1 or in_p2)) - // this molecule is not affected by the restraint - return; - - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by an angle restraint."), - CODELOC); -} - -/** Calculate the force acting on the molecules in the forcetable 'forcetable' - caused by this restraint, and add it on to the forcetable scaled by - 'scale_force' */ -void AngleRestraint::force(ForceTable &forcetable, double scale_force) const -{ - bool in_p0 = p[0].read().usesMoleculesIn(forcetable); - bool in_p1 = p[1].read().usesMoleculesIn(forcetable); - bool in_p2 = p[2].read().usesMoleculesIn(forcetable); - - if (not(in_p0 or in_p1 or in_p2)) - // this molecule is not affected by the restraint - return; - - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by an angle restraint."), - CODELOC); -} - -/** Update the points of this restraint using new molecule data from 'moldata' - - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void AngleRestraint::update(const MoleculeData &moldata) -{ - if (this->contains(moldata.number())) - { - AngleRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(moldata); - } - - this->calculateTheta(); - } - catch (...) - { - AngleRestraint::operator=(old_state); - throw; - } - } -} - -/** Update the points of this restraint using new molecule data from 'molecules' - - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void AngleRestraint::update(const Molecules &molecules) -{ - if (this->usesMoleculesIn(molecules)) - { - AngleRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(molecules); - } - - this->calculateTheta(); - } - catch (...) - { - AngleRestraint::operator=(old_state); - throw; - } - } -} - -/** Return the molecules used in this restraint */ -Molecules AngleRestraint::molecules() const -{ - Molecules mols; - - for (int i = 0; i < this->nPoints(); ++i) - { - mols += p[i].read().molecules(); - } - - return mols; -} - -/** Return whether or not this restraint affects the molecule - with number 'molnum' */ -bool AngleRestraint::contains(MolNum molnum) const -{ - return p[0].read().contains(molnum) or p[1].read().contains(molnum) or p[2].read().contains(molnum); -} - -/** Return whether or not this restraint affects the molecule - with ID 'molid' */ -bool AngleRestraint::contains(const MolID &molid) const -{ - return p[0].read().contains(molid) or p[1].read().contains(molid) or p[2].read().contains(molid); -} - -/** Return whether or not this restraint involves any of the molecules - that are in the forcetable 'forcetable' */ -bool AngleRestraint::usesMoleculesIn(const ForceTable &forcetable) const -{ - return p[0].read().usesMoleculesIn(forcetable) or p[1].read().usesMoleculesIn(forcetable) or - p[2].read().usesMoleculesIn(forcetable); -} - -/** Return whether or not this restraint involves any of the molecules - in 'molecules' */ -bool AngleRestraint::usesMoleculesIn(const Molecules &molecules) const -{ - return p[0].read().usesMoleculesIn(molecules) or p[1].read().usesMoleculesIn(molecules) or - p[2].read().usesMoleculesIn(molecules); -} - -static Expression harmonicFunction(double force_constant) -{ - if (SireMaths::isZero(force_constant)) - return 0; - else - return force_constant * pow(AngleRestraint::theta(), 2); -} - -static Expression diffHarmonicFunction(double force_constant) -{ - if (SireMaths::isZero(force_constant)) - return 0; - else - return (2 * force_constant) * AngleRestraint::theta(); -} - -/** Return a distance restraint that applies a harmonic potential between - the points 'point0' and 'point1' using a force constant 'force_constant' */ -AngleRestraint AngleRestraint::harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const HarmonicAngleForceConstant &force_constant) -{ - return AngleRestraint(point0, point1, point2, ::harmonicFunction(force_constant), - ::diffHarmonicFunction(force_constant)); -} - -static Expression halfHarmonicFunction(double force_constant, double angle) -{ - if (SireMaths::isZero(force_constant)) - return 0; - - else if (angle <= 0) - // this is just a harmonic function - return ::harmonicFunction(force_constant); - - else - { - const Symbol &theta = AngleRestraint::theta(); - return Conditional(GreaterThan(theta, angle), force_constant * pow(theta - angle, 2), 0); - } -} - -static Expression diffHalfHarmonicFunction(double force_constant, double angle) -{ - if (SireMaths::isZero(force_constant)) - return 0; - - else if (angle <= 0) - // this is just a harmonic function - return ::diffHarmonicFunction(force_constant); - - else - { - const Symbol &theta = AngleRestraint::theta(); - return Conditional(GreaterThan(theta, angle), (2 * force_constant) * (theta - angle), 0); - } -} - -/** Return a distance restraint that applied a half-harmonic potential - between the points 'point0' and 'point1' above a distance 'distance' - using a force constant 'force_constant' */ -AngleRestraint AngleRestraint::halfHarmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Angle &angle, const HarmonicAngleForceConstant &force_constant) -{ - double acute_angle = acute(angle).to(radians); - - return AngleRestraint(point0, point1, point2, ::halfHarmonicFunction(force_constant, acute_angle), - ::diffHalfHarmonicFunction(force_constant, acute_angle)); -} diff --git a/corelib/src/libs/SireMM/anglerestraint.h b/corelib/src/libs/SireMM/anglerestraint.h deleted file mode 100644 index 6901c78ca..000000000 --- a/corelib/src/libs/SireMM/anglerestraint.h +++ /dev/null @@ -1,159 +0,0 @@ -/********************************************\ - * - * Sire - Molecular Simulation Framework - * - * Copyright (C) 2009 Christopher Woods - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * - * For full details of the license please see the COPYING file - * that should have come with this distribution. - * - * You can contact the authors at https://sire.openbiosim.org - * -\*********************************************/ - -#ifndef SIREMM_ANGLERESTRAINT_H -#define SIREMM_ANGLERESTRAINT_H - -#include "SireFF/point.h" - -#include "restraint.h" - -#include "SireCAS/expression.h" -#include "SireCAS/symbol.h" - -#include "SireUnits/dimensions.h" - -SIRE_BEGIN_HEADER - -namespace SireMM -{ - class AngleRestraint; -} - -SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::AngleRestraint &); -SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::AngleRestraint &); - -namespace SireMM -{ - - using SireFF::Point; - using SireFF::PointRef; - - using SireCAS::Expression; - using SireCAS::Symbol; - using SireCAS::Symbols; - - // typedef the unit of a harmonic force constant ( MolarEnergy / Angle^2 ) - typedef SireUnits::Dimension::PhysUnit<1, 2, -2, 0, 0, -1, -2> HarmonicAngleForceConstant; - - /** This is a restraint that operates on the angle between - three SireMM::Point objects (e.g. three atoms in a molecule) - - @author Christopher Woods - */ - class SIREMM_EXPORT AngleRestraint : public SireBase::ConcreteProperty - { - - friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const AngleRestraint &); - friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, AngleRestraint &); - - public: - AngleRestraint(); - - AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const Expression &restraint); - - AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const Expression &restraint, - const Values &values); - - AngleRestraint(const AngleRestraint &other); - - ~AngleRestraint(); - - AngleRestraint &operator=(const AngleRestraint &other); - - bool operator==(const AngleRestraint &other) const; - bool operator!=(const AngleRestraint &other) const; - - static const char *typeName(); - - const Point &point(int i) const; - - const Point &point0() const; - const Point &point1() const; - const Point &point2() const; - - int nPoints() const; - - static const Symbol &theta(); - - Symbols builtinSymbols() const; - Values builtinValues() const; - - RestraintPtr differentiate(const Symbol &symbol) const; - - void setSpace(const Space &space); - - const Expression &differentialRestraintFunction() const; - - void force(MolForceTable &forcetable, double scale_force = 1) const; - void force(ForceTable &forcetable, double scale_force = 1) const; - - void update(const MoleculeData &moldata); - - void update(const Molecules &molecules); - - Molecules molecules() const; - - bool contains(MolNum molnum) const; - bool contains(const MolID &molid) const; - - bool usesMoleculesIn(const ForceTable &forcetable) const; - bool usesMoleculesIn(const Molecules &molecules) const; - - static AngleRestraint harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const HarmonicAngleForceConstant &force_constant); - - static AngleRestraint halfHarmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const SireUnits::Dimension::Angle &angle, - const HarmonicAngleForceConstant &force_constant); - - protected: - AngleRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const Expression &nrg_restraint, const Expression &force_restraint); - - private: - void calculateTheta(); - - /** The three points between which the restraint is calculated */ - SireFF::PointPtr p[3]; - - /** The expression used to calculate the force */ - Expression force_expression; - - /** Whether or not all three points are within the same molecule */ - bool intra_molecule_points; - }; - -} // namespace SireMM - -Q_DECLARE_METATYPE(SireMM::AngleRestraint) - -SIRE_EXPOSE_CLASS(SireMM::AngleRestraint) - -SIRE_END_HEADER - -#endif diff --git a/corelib/src/libs/SireMM/anglerestraints.cpp b/corelib/src/libs/SireMM/anglerestraints.cpp new file mode 100644 index 000000000..5eaf82fc1 --- /dev/null +++ b/corelib/src/libs/SireMM/anglerestraints.cpp @@ -0,0 +1,479 @@ +/********************************************\ + * + * Sire - Molecular Simulation Framework + * + * Copyright (C) 2025 Christopher Woods + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * For full details of the license please see the COPYING file + * that should have come with this distribution. + * + * You can contact the authors at https://sire.openbiosim.org + * +\*********************************************/ + +#include "anglerestraints.h" + +#include "SireID/index.h" + +#include "SireUnits/units.h" + +#include "SireStream/datastream.h" +#include "SireStream/shareddatastream.h" + +#include "SireCAS/errors.h" + +#include + +using namespace SireMM; +using namespace SireID; +using namespace SireBase; +using namespace SireMaths; +using namespace SireStream; +using namespace SireUnits; +using namespace SireUnits::Dimension; + +//////////// +//////////// Implementation of AngleRestraint +//////////// + +static const RegisterMetaType r_angrest; + +/** Serialise to a binary datastream */ + +QDataStream &operator<<(QDataStream &ds, const AngleRestraint &angrest) +{ + writeHeader(ds, r_angrest, 1); + + SharedDataStream sds(ds); + + sds << angrest.atms << angrest._theta0 << angrest._ktheta; + + return ds; +} + +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, AngleRestraint &angrest) +{ + VersionID v = readHeader(ds, r_angrest); + + if (v == 1) + { + SharedDataStream sds(ds); + + sds >> angrest.atms >> angrest._theta0 >> angrest._ktheta; + } + else + throw version_error(v, "1", r_angrest, CODELOC); + + return ds; +} + +/** Null constructor */ +AngleRestraint::AngleRestraint() + : ConcreteProperty(), + _ktheta(0), _theta0(0) +{ +} + +/** Construct a restraint that acts on the angle within the + three atoms 'atom0', 'atom1' and 'atom2' (theta == a(012)), + restraining the angle within these atoms */ +AngleRestraint::AngleRestraint(const QList &atoms, + const SireUnits::Dimension::Angle &theta0, + const SireUnits::Dimension::HarmonicAngleConstant &ktheta) + : ConcreteProperty(), + _theta0(theta0), _ktheta(ktheta) +{ + + // Make sure that we have 3 distinct atoms + QSet distinct; + distinct.reserve(3); + + for (const auto &atom : atoms) + { + if (atom >= 0) + distinct.insert(atom); + } + + atms = atoms.toVector(); +} + +/* Copy constructor*/ +AngleRestraint::AngleRestraint(const AngleRestraint &other) + : ConcreteProperty(other), + atms(other.atms), _theta0(other._theta0), _ktheta(other._ktheta) + +{ +} + +/* Destructor */ +AngleRestraint::~AngleRestraint() +{ +} + +AngleRestraint &AngleRestraint::operator=(const AngleRestraint &other) +{ + if (this != &other) + { + Property::operator=(other); + atms = other.atms; + _theta0 = other._theta0; + _ktheta = other._ktheta; + } + + return *this; +} + +bool AngleRestraint::operator==(const AngleRestraint &other) const +{ + return atms == other.atms and + _theta0 == other._theta0 and + _ktheta == other._ktheta; +} + +bool AngleRestraint::operator!=(const AngleRestraint &other) const +{ + return not operator==(other); +} + +AngleRestraints AngleRestraint::operator+(const AngleRestraint &other) const +{ + return AngleRestraints(*this) + other; +} + +AngleRestraints AngleRestraint::operator+(const AngleRestraints &other) const +{ + return AngleRestraints(*this) + other; +} + +const char *AngleRestraint::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); +} + +const char *AngleRestraint::what() const +{ + return AngleRestraint::typeName(); +} + +AngleRestraint *AngleRestraint::clone() const +{ + return new AngleRestraint(*this); +} + +bool AngleRestraint::isNull() const +{ + return atms.isEmpty(); +} + +QString AngleRestraint::toString() const +{ + if (this->isNull()) + return QObject::tr("AngleRestraint::null"); + else + { + QStringList a; + + for (const auto &atom : atms) + { + a.append(QString::number(atom)); + } + return QString("AngleRestraint( [%1], theta0=%2, ktheta=%3 )") + .arg(a.join(", ")) + .arg(_theta0.toString()) + .arg(_ktheta.toString()); + } +} + +/** Return the force constant for the restraint */ +SireUnits::Dimension::HarmonicAngleConstant AngleRestraint::ktheta() const +{ + return this->_ktheta; +} + +/** Return the equilibrium angle for the restraint */ +SireUnits::Dimension::Angle AngleRestraint::theta0() const +{ + return this->_theta0; +} + +/** Return the atoms involved in the restraint */ +QVector AngleRestraint::atoms() const +{ + return this->atms; +} + +/////// +/////// Implementation of AngleRestraints +/////// + +/** Serialise to a binary datastream */ + +static const RegisterMetaType r_angrests; + +QDataStream &operator<<(QDataStream &ds, const AngleRestraints &angrests) +{ + writeHeader(ds, r_angrests, 1); + + SharedDataStream sds(ds); + + sds << angrests.r + << static_cast(angrests); + + return ds; +} + +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, AngleRestraints &angrests) +{ + VersionID v = readHeader(ds, r_angrests); + + if (v == 1) + { + SharedDataStream sds(ds); + + sds >> angrests.r >> + static_cast(angrests); + } + else + throw version_error(v, "1", r_angrests, CODELOC); + + return ds; +} + +/** Null constructor */ +AngleRestraints::AngleRestraints() + : ConcreteProperty() +{ +} + +AngleRestraints::AngleRestraints(const QString &name) + : ConcreteProperty(name) +{ +} + +AngleRestraints::AngleRestraints(const AngleRestraint &restraint) + : ConcreteProperty() +{ + if (not restraint.isNull()) + r.append(restraint); +} + +AngleRestraints::AngleRestraints(const QList &restraints) + : ConcreteProperty() +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } +} + +AngleRestraints::AngleRestraints(const QString &name, + const AngleRestraint &restraint) + : ConcreteProperty(name) +{ + if (not restraint.isNull()) + r.append(restraint); +} + +AngleRestraints::AngleRestraints(const QString &name, + const QList &restraints) + : ConcreteProperty(name) +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } +} + +AngleRestraints::AngleRestraints(const AngleRestraints &other) + : ConcreteProperty(other), r(other.r) +{ +} + +/* Desctructor */ +AngleRestraints::~AngleRestraints() +{ +} + +AngleRestraints &AngleRestraints::operator=(const AngleRestraints &other) +{ + if (this != &other) + { + Restraints::operator=(other); + r = other.r; + } + + return *this; +} + +bool AngleRestraints::operator==(const AngleRestraints &other) const +{ + return r == other.r and Restraints::operator==(other); +} + +bool AngleRestraints::operator!=(const AngleRestraints &other) const +{ + return not operator==(other); +} + +const char *AngleRestraints::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); +} + +const char *AngleRestraints::what() const +{ + return AngleRestraints::typeName(); +} + +AngleRestraints *AngleRestraints::clone() const +{ + return new AngleRestraints(*this); +} + +QString AngleRestraints::toString() const +{ + if (this->isEmpty()) + return QObject::tr("AngleRestraints::null"); + + QStringList parts; + + const auto n = this->count(); + + if (n <= 10) + { + for (int i = 0; i < n; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + } + else + { + for (int i = 0; i < 5; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + + parts.append("..."); + + for (int i = n - 5; i < n; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + } + + return QObject::tr("AngleRestraints( name=%1, size=%2\n%3\n )") + .arg(this->name()) + .arg(n) + .arg(parts.join("\n")); +} + +/** Return whether or not this is empty */ +bool AngleRestraints::isEmpty() const +{ + return this->r.isEmpty(); +} + +/** Return whether or not this is empty */ +bool AngleRestraints::isNull() const +{ + return this->isEmpty(); +} + +/** Return the number of restraints */ +int AngleRestraints::nRestraints() const +{ + return this->r.count(); +} + +/** Return the number of restraints */ +int AngleRestraints::count() const +{ + return this->nRestraints(); +} + +/** Return the number of restraints */ +int AngleRestraints::size() const +{ + return this->nRestraints(); +} + +/** Return the ith restraint */ +const AngleRestraint &AngleRestraints::at(int i) const +{ + i = SireID::Index(i).map(this->r.count()); + + return this->r.at(i); +} + +/** Return the ith restraint */ +const AngleRestraint &AngleRestraints::operator[](int i) const +{ + return this->at(i); +} + +/** Return all of the restraints */ +QList AngleRestraints::restraints() const +{ + return this->r; +} + +/** Add a restraints onto the list */ +void AngleRestraints::add(const AngleRestraint &restraint) +{ + if (not restraint.isNull()) + r.append(restraint); +} + +/** Add a restraint onto the list */ +void AngleRestraints::add(const AngleRestraints &restraints) +{ + this->r += restraints.r; +} + +/** Add a restraint onto the list */ +AngleRestraints &AngleRestraints::operator+=(const AngleRestraint &restraint) +{ + this->add(restraint); + return *this; +} + +/** Add a restraint onto the list */ +AngleRestraints AngleRestraints::operator+(const AngleRestraint &restraint) const +{ + AngleRestraints ret(*this); + ret += restraint; + return *this; +} + +/** Add restraints onto the list */ +AngleRestraints &AngleRestraints::operator+=(const AngleRestraints &restraints) +{ + this->add(restraints); + return *this; +} + +/** Add restraints onto the list */ +AngleRestraints AngleRestraints::operator+(const AngleRestraints &restraints) const +{ + AngleRestraints ret(*this); + ret += restraints; + return *this; +} diff --git a/corelib/src/libs/SireMM/anglerestraints.h b/corelib/src/libs/SireMM/anglerestraints.h new file mode 100644 index 000000000..b34de4822 --- /dev/null +++ b/corelib/src/libs/SireMM/anglerestraints.h @@ -0,0 +1,177 @@ +/********************************************\ + * + * Sire - Molecular Simulation Framework + * + * Copyright (C) 2025 Christopher Woods + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * For full details of the license please see the COPYING file + * that should have come with this distribution. + * + * You can contact the authors at https://sire.openbiosim.org + * +\*********************************************/ + +#ifndef SIREMM_ANGLERESTRAINTS_H +#define SIREMM_ANGLERESTRAINTS_H + + +#include "restraints.h" +#include "SireUnits/dimensions.h" +#include "SireUnits/generalunit.h" + +SIRE_BEGIN_HEADER + +namespace SireMM +{ + class AngleRestraint; + class AngleRestraints; +} + +SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::AngleRestraint &); +SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::AngleRestraint &); + +SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::AngleRestraints &); +SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::AngleRestraints &); + +namespace SireMM +{ + + /** This class represents a single angle restraint between any three + * atoms in a system + * @author Christopher Woods + */ + class SIREMM_EXPORT AngleRestraint + : public SireBase::ConcreteProperty + { + + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::AngleRestraint &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::AngleRestraint &); + + public: + AngleRestraint(); + AngleRestraint(const QList &atoms, + const SireUnits::Dimension::Angle &theta0, + const SireUnits::Dimension::HarmonicAngleConstant &ktheta); + + AngleRestraint(const AngleRestraint &other); + + ~AngleRestraint(); + + AngleRestraint &operator=(const AngleRestraint &other); + + bool operator==(const AngleRestraint &other) const; + bool operator!=(const AngleRestraint &other) const; + + AngleRestraints operator+(const AngleRestraint &other) const; + AngleRestraints operator+(const AngleRestraints &other) const; + + static const char *typeName(); + const char *what() const; + + AngleRestraint *clone() const; + + QString toString() const; + + bool isNull() const; + + QVector atoms() const; + + SireUnits::Dimension::Angle theta0() const; + SireUnits::Dimension::HarmonicAngleConstant ktheta() const; + + private: + /** Atoms involved in the angle restraint */ + QVector atms; + + /** Equilibrium angle */ + SireUnits::Dimension::Angle _theta0; + + /** Harmonic angle constant */ + SireUnits::Dimension::HarmonicAngleConstant _ktheta; + }; + + /** This class represents a collection of angle restraints */ + class SIREMM_EXPORT AngleRestraints + : public SireBase::ConcreteProperty + { + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::AngleRestraints &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::AngleRestraints &); + + public: + AngleRestraints(); + AngleRestraints(const QString &name); + + AngleRestraints(const AngleRestraint &restraint); + AngleRestraints(const QList &restraints); + + AngleRestraints(const QString &name, + const AngleRestraint &restraint); + + AngleRestraints(const QString &name, + const QList &restraints); + + AngleRestraints(const AngleRestraints &other); + + ~AngleRestraints(); + + AngleRestraints &operator=(const AngleRestraints &other); + + bool operator==(const AngleRestraints &other) const; + bool operator!=(const AngleRestraints &other) const; + + static const char *typeName(); + const char *what() const; + + AngleRestraints *clone() const; + + QString toString() const; + + bool isEmpty() const; + bool isNull() const; + + int count() const; + int size() const; + int nRestraints() const; + + const AngleRestraint &at(int i) const; + const AngleRestraint &operator[](int i) const; + + QList restraints() const; + + void add(const AngleRestraint &restraint); + void add(const AngleRestraints &restraints); + + AngleRestraints &operator+=(const AngleRestraint &restraint); + AngleRestraints &operator+=(const AngleRestraints &restraints); + + AngleRestraints operator+(const AngleRestraint &restraint) const; + AngleRestraints operator+(const AngleRestraints &restraints) const; + + private: + /** List of restraints */ + QList r; + }; +} + +Q_DECLARE_METATYPE(SireMM::AngleRestraint) +Q_DECLARE_METATYPE(SireMM::AngleRestraints) + +SIRE_EXPOSE_CLASS(SireMM::AngleRestraint) +SIRE_EXPOSE_CLASS(SireMM::AngleRestraints) +SIRE_END_HEADER + +#endif diff --git a/corelib/src/libs/SireMM/dihedralrestraint.cpp b/corelib/src/libs/SireMM/dihedralrestraint.cpp deleted file mode 100644 index bea13ed94..000000000 --- a/corelib/src/libs/SireMM/dihedralrestraint.cpp +++ /dev/null @@ -1,576 +0,0 @@ -/********************************************\ - * - * Sire - Molecular Simulation Framework - * - * Copyright (C) 2009 Christopher Woods - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * - * For full details of the license please see the COPYING file - * that should have come with this distribution. - * - * You can contact the authors at https://sire.openbiosim.org - * -\*********************************************/ - -#include "dihedralrestraint.h" - -#include "SireFF/forcetable.h" - -#include "SireCAS/conditional.h" -#include "SireCAS/power.h" -#include "SireCAS/symbols.h" -#include "SireCAS/values.h" - -#include "SireID/index.h" - -#include "SireUnits/angle.h" -#include "SireUnits/units.h" - -#include "SireStream/datastream.h" -#include "SireStream/shareddatastream.h" - -#include "SireCAS/errors.h" - -using namespace SireMM; -using namespace SireMol; -using namespace SireFF; -using namespace SireID; -using namespace SireBase; -using namespace SireCAS; -using namespace SireMaths; -using namespace SireStream; -using namespace SireUnits; -using namespace SireUnits::Dimension; - -//////////// -//////////// Implementation of DihedralRestraint -//////////// - -static const RegisterMetaType r_dihrest; - -/** Serialise to a binary datastream */ -QDataStream &operator<<(QDataStream &ds, const DihedralRestraint &dihrest) -{ - writeHeader(ds, r_dihrest, 1); - - SharedDataStream sds(ds); - - sds << dihrest.p[0] << dihrest.p[1] << dihrest.p[2] << dihrest.p[3] << dihrest.force_expression - << static_cast(dihrest); - - return ds; -} - -/** Extract from a binary datastream */ -QDataStream &operator>>(QDataStream &ds, DihedralRestraint &dihrest) -{ - VersionID v = readHeader(ds, r_dihrest); - - if (v == 1) - { - SharedDataStream sds(ds); - - sds >> dihrest.p[0] >> dihrest.p[1] >> dihrest.p[2] >> dihrest.p[3] >> dihrest.force_expression >> - static_cast(dihrest); - - dihrest.intra_molecule_points = Point::areIntraMoleculePoints(dihrest.p[0], dihrest.p[1]) and - Point::areIntraMoleculePoints(dihrest.p[0], dihrest.p[2]) and - Point::areIntraMoleculePoints(dihrest.p[0], dihrest.p[3]); - } - else - throw version_error(v, "1", r_dihrest, CODELOC); - - return ds; -} - -Q_GLOBAL_STATIC_WITH_ARGS(Symbol, getPhiSymbol, ("phi")); - -/** Return the symbol that represents the dihedral angle between the points (phi) */ -const Symbol &DihedralRestraint::phi() -{ - return *(getPhiSymbol()); -} - -/** Constructor */ -DihedralRestraint::DihedralRestraint() : ConcreteProperty() -{ -} - -void DihedralRestraint::calculatePhi() -{ - if (this->restraintFunction().isFunction(phi())) - { - SireUnits::Dimension::Angle angle; - - if (intra_molecule_points) - // we don't use the space when calculating intra-molecular angles - angle = - Vector::dihedral(p[0].read().point(), p[1].read().point(), p[2].read().point(), p[3].read().point()); - else - angle = this->space().calcDihedral(p[0].read().point(), p[1].read().point(), p[2].read().point(), - p[3].read().point()); - - ExpressionRestraint3D::_pvt_setValue(phi(), angle); - } -} - -/** Construct a restraint that acts on the angle within the - three points 'point0', 'point1' and 'point2' (theta == a(012)), - restraining the angle within these points using the expression - 'restraint' */ -DihedralRestraint::DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const Expression &restraint) - : ConcreteProperty(restraint) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - p[3] = point3; - - force_expression = this->restraintFunction().differentiate(phi()); - - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]) and - Point::areIntraMoleculePoints(p[0], p[3]); - - this->calculatePhi(); -} - -/** Construct a restraint that acts on the angle within the - three points 'point0', 'point1' and 'point2' (theta == a(012)), - restraining the angle within these points using the expression - 'restraint' */ -DihedralRestraint::DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const Expression &restraint, const Values &values) - : ConcreteProperty(restraint, values) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - p[3] = point3; - - force_expression = this->restraintFunction().differentiate(phi()); - - if (force_expression.isConstant()) - force_expression = force_expression.evaluate(Values()); - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]) and - Point::areIntraMoleculePoints(p[0], p[3]); - - this->calculatePhi(); -} - -/** Internal constructor used to construct a restraint using the specified - points, energy expression and force expression */ -DihedralRestraint::DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const Expression &nrg_restraint, - const Expression &force_restraint) - : ConcreteProperty(nrg_restraint), force_expression(force_restraint) -{ - p[0] = point0; - p[1] = point1; - p[2] = point2; - p[3] = point3; - - if (force_expression.isConstant()) - { - force_expression = force_expression.evaluate(Values()); - } - else - { - if (not this->restraintFunction().symbols().contains(force_expression.symbols())) - throw SireError::incompatible_error( - QObject::tr("You cannot use a force function which uses more symbols " - "(%1) than the energy function (%2).") - .arg(Sire::toString(force_expression.symbols()), Sire::toString(restraintFunction().symbols())), - CODELOC); - } - - intra_molecule_points = Point::areIntraMoleculePoints(p[0], p[1]) and Point::areIntraMoleculePoints(p[0], p[2]) and - Point::areIntraMoleculePoints(p[0], p[3]); - - this->calculatePhi(); -} - -/** Copy constructor */ -DihedralRestraint::DihedralRestraint(const DihedralRestraint &other) - : ConcreteProperty(other), force_expression(other.force_expression), - intra_molecule_points(other.intra_molecule_points) -{ - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } -} - -/** Destructor */ -DihedralRestraint::~DihedralRestraint() -{ -} - -/** Copy assignment operator */ -DihedralRestraint &DihedralRestraint::operator=(const DihedralRestraint &other) -{ - if (this != &other) - { - ExpressionRestraint3D::operator=(other); - - for (int i = 0; i < this->nPoints(); ++i) - { - p[i] = other.p[i]; - } - - force_expression = other.force_expression; - intra_molecule_points = other.intra_molecule_points; - } - - return *this; -} - -/** Comparison operator */ -bool DihedralRestraint::operator==(const DihedralRestraint &other) const -{ - return this == &other or (ExpressionRestraint3D::operator==(other) and p[0] == other.p[0] and p[1] == other.p[1] and - p[2] == other.p[2] and p[3] == other.p[3] and force_expression == other.force_expression); -} - -/** Comparison operator */ -bool DihedralRestraint::operator!=(const DihedralRestraint &other) const -{ - return not DihedralRestraint::operator==(other); -} - -const char *DihedralRestraint::typeName() -{ - return QMetaType::typeName(qMetaTypeId()); -} - -/** This restraint involves four points */ -int DihedralRestraint::nPoints() const -{ - return 4; -} - -/** Return the ith point */ -const Point &DihedralRestraint::point(int i) const -{ - i = Index(i).map(this->nPoints()); - - return p[i].read(); -} - -/** Return the first point */ -const Point &DihedralRestraint::point0() const -{ - return p[0].read(); -} - -/** Return the second point */ -const Point &DihedralRestraint::point1() const -{ - return p[1].read(); -} - -/** Return the third point */ -const Point &DihedralRestraint::point2() const -{ - return p[2].read(); -} - -/** Return the fourth point */ -const Point &DihedralRestraint::point3() const -{ - return p[3].read(); -} - -/** Return the built-in symbols for this restraint */ -Symbols DihedralRestraint::builtinSymbols() const -{ - if (this->restraintFunction().isFunction(phi())) - return phi(); - else - return Symbols(); -} - -/** Return the values of the built-in symbols of this restraint */ -Values DihedralRestraint::builtinValues() const -{ - if (this->restraintFunction().isFunction(phi())) - return phi() == this->values()[phi()]; - else - return Values(); -} - -/** Return the differential of this restraint with respect to - the symbol 'symbol' - - \throw SireCAS::unavailable_differential -*/ -RestraintPtr DihedralRestraint::differentiate(const Symbol &symbol) const -{ - if (this->restraintFunction().isFunction(symbol)) - return DihedralRestraint(p[0], p[1], p[2], p[3], restraintFunction().differentiate(symbol), this->values()); - else - return NullRestraint(); -} - -/** Set the space used to evaluate the energy of this restraint - - \throw SireVol::incompatible_space -*/ -void DihedralRestraint::setSpace(const Space &new_space) -{ - if (not this->space().equals(new_space)) - { - DihedralRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().setSpace(new_space); - } - - Restraint3D::setSpace(new_space); - - this->calculatePhi(); - } - catch (...) - { - DihedralRestraint::operator=(old_state); - throw; - } - } -} - -/** Return the function used to calculate the restraint force */ -const Expression &DihedralRestraint::differentialRestraintFunction() const -{ - return force_expression; -} - -/** Calculate the force acting on the molecule in the forcetable 'forcetable' - caused by this restraint, and add it on to the forcetable scaled by - 'scale_force' */ -void DihedralRestraint::force(MolForceTable &forcetable, double scale_force) const -{ - bool in_p0 = p[0].read().contains(forcetable.molNum()); - bool in_p1 = p[1].read().contains(forcetable.molNum()); - bool in_p2 = p[2].read().contains(forcetable.molNum()); - bool in_p3 = p[3].read().contains(forcetable.molNum()); - - if (not(in_p0 or in_p1 or in_p2 or in_p3)) - // this molecule is not affected by the restraint - return; - - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by a dihedral restraint."), - CODELOC); -} - -/** Calculate the force acting on the molecules in the forcetable 'forcetable' - caused by this restraint, and add it on to the forcetable scaled by - 'scale_force' */ -void DihedralRestraint::force(ForceTable &forcetable, double scale_force) const -{ - bool in_p0 = p[0].read().usesMoleculesIn(forcetable); - bool in_p1 = p[1].read().usesMoleculesIn(forcetable); - bool in_p2 = p[2].read().usesMoleculesIn(forcetable); - bool in_p3 = p[3].read().usesMoleculesIn(forcetable); - - if (not(in_p0 or in_p1 or in_p2 or in_p3)) - // this molecule is not affected by the restraint - return; - - throw SireError::incomplete_code(QObject::tr("Haven't yet written the code to calculate forces caused " - "by a dihedral restraint."), - CODELOC); -} - -/** Update the points of this restraint using new molecule data from 'moldata' - - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void DihedralRestraint::update(const MoleculeData &moldata) -{ - if (this->contains(moldata.number())) - { - DihedralRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(moldata); - } - - this->calculatePhi(); - } - catch (...) - { - DihedralRestraint::operator=(old_state); - throw; - } - } -} - -/** Update the points of this restraint using new molecule data from 'molecules' - - \throw SireBase::missing_property - \throw SireError::invalid_cast - \throw SireError::incompatible_error -*/ -void DihedralRestraint::update(const Molecules &molecules) -{ - if (this->usesMoleculesIn(molecules)) - { - DihedralRestraint old_state(*this); - - try - { - for (int i = 0; i < this->nPoints(); ++i) - { - p[i].edit().update(molecules); - } - - this->calculatePhi(); - } - catch (...) - { - DihedralRestraint::operator=(old_state); - throw; - } - } -} - -/** Return the molecules used in this restraint */ -Molecules DihedralRestraint::molecules() const -{ - Molecules mols; - - for (int i = 0; i < this->nPoints(); ++i) - { - mols += p[i].read().molecules(); - } - - return mols; -} - -/** Return whether or not this restraint affects the molecule - with number 'molnum' */ -bool DihedralRestraint::contains(MolNum molnum) const -{ - return p[0].read().contains(molnum) or p[1].read().contains(molnum) or p[2].read().contains(molnum) or - p[3].read().contains(molnum); -} - -/** Return whether or not this restraint affects the molecule - with ID 'molid' */ -bool DihedralRestraint::contains(const MolID &molid) const -{ - return p[0].read().contains(molid) or p[1].read().contains(molid) or p[2].read().contains(molid) or - p[3].read().contains(molid); -} - -/** Return whether or not this restraint involves any of the molecules - that are in the forcetable 'forcetable' */ -bool DihedralRestraint::usesMoleculesIn(const ForceTable &forcetable) const -{ - return p[0].read().usesMoleculesIn(forcetable) or p[1].read().usesMoleculesIn(forcetable) or - p[2].read().usesMoleculesIn(forcetable) or p[3].read().usesMoleculesIn(forcetable); -} - -/** Return whether or not this restraint involves any of the molecules - in 'molecules' */ -bool DihedralRestraint::usesMoleculesIn(const Molecules &molecules) const -{ - return p[0].read().usesMoleculesIn(molecules) or p[1].read().usesMoleculesIn(molecules) or - p[2].read().usesMoleculesIn(molecules) or p[3].read().usesMoleculesIn(molecules); -} - -static Expression harmonicFunction(double force_constant) -{ - if (SireMaths::isZero(force_constant)) - return 0; - else - return force_constant * pow(DihedralRestraint::phi(), 2); -} - -static Expression diffHarmonicFunction(double force_constant) -{ - if (SireMaths::isZero(force_constant)) - return 0; - else - return (2 * force_constant) * DihedralRestraint::phi(); -} - -/** Return a distance restraint that applies a harmonic potential between - the points 'point0' and 'point1' using a force constant 'force_constant' */ -DihedralRestraint DihedralRestraint::harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const HarmonicAngleForceConstant &force_constant) -{ - return DihedralRestraint(point0, point1, point2, point3, ::harmonicFunction(force_constant), - ::diffHarmonicFunction(force_constant)); -} - -static Expression halfHarmonicFunction(double force_constant, double angle) -{ - if (SireMaths::isZero(force_constant)) - return 0; - - else if (angle <= 0) - // this is just a harmonic function - return ::harmonicFunction(force_constant); - - else - { - const Symbol &phi = DihedralRestraint::phi(); - return Conditional(GreaterThan(phi, angle), force_constant * pow(phi - angle, 2), 0); - } -} - -static Expression diffHalfHarmonicFunction(double force_constant, double angle) -{ - if (SireMaths::isZero(force_constant)) - return 0; - - else if (angle <= 0) - // this is just a harmonic function - return ::diffHarmonicFunction(force_constant); - - else - { - const Symbol &phi = DihedralRestraint::phi(); - return Conditional(GreaterThan(phi, angle), (2 * force_constant) * (phi - angle), 0); - } -} - -/** Return a distance restraint that applied a half-harmonic potential - between the points 'point0' and 'point1' above a distance 'distance' - using a force constant 'force_constant' */ -DihedralRestraint DihedralRestraint::halfHarmonic(const PointRef &point0, const PointRef &point1, - const PointRef &point2, const PointRef &point3, const Angle &angle, - const HarmonicAngleForceConstant &force_constant) -{ - double ang = angle.to(radians); - - return DihedralRestraint(point0, point1, point2, point3, ::halfHarmonicFunction(force_constant, ang), - ::diffHalfHarmonicFunction(force_constant, ang)); -} diff --git a/corelib/src/libs/SireMM/dihedralrestraint.h b/corelib/src/libs/SireMM/dihedralrestraint.h deleted file mode 100644 index 606ee3b78..000000000 --- a/corelib/src/libs/SireMM/dihedralrestraint.h +++ /dev/null @@ -1,144 +0,0 @@ -/********************************************\ - * - * Sire - Molecular Simulation Framework - * - * Copyright (C) 2009 Christopher Woods - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - * - * For full details of the license please see the COPYING file - * that should have come with this distribution. - * - * You can contact the authors at https://sire.openbiosim.org - * -\*********************************************/ - -#ifndef SIREMM_DIHEDRALRESTRAINT_H -#define SIREMM_DIHEDRALRESTRAINT_H - -#include "anglerestraint.h" -#include "restraint.h" - -SIRE_BEGIN_HEADER - -namespace SireMM -{ - class DihedralRestraint; -} - -SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::DihedralRestraint &); -SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::DihedralRestraint &); - -namespace SireMM -{ - - /** This is a restraint that operates on the dihedral angle between - four SireMM::Point objects (e.g. four atoms in a molecule) - - @author Christopher Woods - */ - class SIREMM_EXPORT DihedralRestraint : public SireBase::ConcreteProperty - { - - friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const DihedralRestraint &); - friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, DihedralRestraint &); - - public: - DihedralRestraint(); - - DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const PointRef &point3, - const Expression &restraint); - - DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const PointRef &point3, - const Expression &restraint, const Values &values); - - DihedralRestraint(const DihedralRestraint &other); - - ~DihedralRestraint(); - - DihedralRestraint &operator=(const DihedralRestraint &other); - - bool operator==(const DihedralRestraint &other) const; - bool operator!=(const DihedralRestraint &other) const; - - static const char *typeName(); - - const Point &point(int i) const; - - const Point &point0() const; - const Point &point1() const; - const Point &point2() const; - const Point &point3() const; - - int nPoints() const; - - static const Symbol &phi(); - - Symbols builtinSymbols() const; - Values builtinValues() const; - - RestraintPtr differentiate(const Symbol &symbol) const; - - void setSpace(const Space &space); - - const Expression &differentialRestraintFunction() const; - - void force(MolForceTable &forcetable, double scale_force = 1) const; - void force(ForceTable &forcetable, double scale_force = 1) const; - - void update(const MoleculeData &moldata); - void update(const Molecules &molecules); - - Molecules molecules() const; - - bool contains(MolNum molnum) const; - bool contains(const MolID &molid) const; - - bool usesMoleculesIn(const ForceTable &forcetable) const; - bool usesMoleculesIn(const Molecules &molecules) const; - - static DihedralRestraint harmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const HarmonicAngleForceConstant &force_constant); - - static DihedralRestraint halfHarmonic(const PointRef &point0, const PointRef &point1, const PointRef &point2, - const PointRef &point3, const SireUnits::Dimension::Angle &angle, - const HarmonicAngleForceConstant &force_constant); - - protected: - DihedralRestraint(const PointRef &point0, const PointRef &point1, const PointRef &point2, const PointRef &point3, - const Expression &nrg_restraint, const Expression &force_restraint); - - private: - void calculatePhi(); - - /** The four points between which the restraint is calculated */ - SireFF::PointPtr p[4]; - - /** The expression used to calculate the force */ - Expression force_expression; - - /** Whether or not all four points are within the same molecule */ - bool intra_molecule_points; - }; - -} // namespace SireMM - -Q_DECLARE_METATYPE(SireMM::DihedralRestraint) - -SIRE_EXPOSE_CLASS(SireMM::DihedralRestraint) - -SIRE_END_HEADER - -#endif diff --git a/corelib/src/libs/SireMM/dihedralrestraints.cpp b/corelib/src/libs/SireMM/dihedralrestraints.cpp new file mode 100644 index 000000000..70054b34c --- /dev/null +++ b/corelib/src/libs/SireMM/dihedralrestraints.cpp @@ -0,0 +1,477 @@ +/********************************************\ + * + * Sire - Molecular Simulation Framework + * + * Copyright (C) 2025 Christopher Woods + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * For full details of the license please see the COPYING file + * that should have come with this distribution. + * + * You can contact the authors at https://sire.openbiosim.org + * +\*********************************************/ + +#include "dihedralrestraints.h" + +#include "SireID/index.h" +#include "SireStream/datastream.h" +#include "SireStream/shareddatastream.h" +#include "SireUnits/units.h" + +#include "SireCAS/errors.h" + +#include + +using namespace SireMM; +using namespace SireID; +using namespace SireBase; +using namespace SireMaths; +using namespace SireStream; +using namespace SireUnits; +using namespace SireUnits::Dimension; + +//////////// +//////////// Implementation of DihedralRestraint +//////////// + +static const RegisterMetaType r_dihrest; + +/** Serialise to a binary datastream */ + +QDataStream &operator<<(QDataStream &ds, const DihedralRestraint &dihrest) +{ + writeHeader(ds, r_dihrest, 1); + + SharedDataStream sds(ds); + + sds << dihrest.atms << dihrest._phi0 << dihrest._kphi; + + return ds; +} + +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, DihedralRestraint &dihrest) +{ + VersionID v = readHeader(ds, r_dihrest); + + if (v == 1) + { + SharedDataStream sds(ds); + + sds >> dihrest.atms >> dihrest._phi0 >> dihrest._kphi; + } + else + throw version_error(v, "1", r_dihrest, CODELOC); + + return ds; +} + +/** Null constructor */ +DihedralRestraint::DihedralRestraint() + : ConcreteProperty(), + _phi0(0), _kphi(0) +{ +} + +/** Construct a restraint that acts on the angle within the + four atoms 'atom0', 'atom1', 'atom2' 'atom3' (phi == a(0123)), + restraining the angle within these atoms */ +DihedralRestraint::DihedralRestraint(const QList &atoms, + const SireUnits::Dimension::Angle &phi0, + const SireUnits::Dimension::HarmonicAngleConstant &kphi) + : ConcreteProperty(), + _phi0(phi0), _kphi(kphi) +{ + + // Make sure that we have 4 distinct atoms + QSet distinct; + distinct.reserve(4); + + for (const auto &atom : atoms) + { + if (atom >= 0) + distinct.insert(atom); + } + + atms = atoms.toVector(); +} + +/* Copy constructor*/ +DihedralRestraint::DihedralRestraint(const DihedralRestraint &other) + : ConcreteProperty(other), + atms(other.atms), _phi0(other._phi0), _kphi(other._kphi) + +{ +} + +/* Destructor */ +DihedralRestraint::~DihedralRestraint() +{ +} + +DihedralRestraint &DihedralRestraint::operator=(const DihedralRestraint &other) +{ + if (this != &other) + { + Property::operator=(other); + atms = other.atms; + _phi0 = other._phi0; + _kphi = other._kphi; + } + + return *this; +} + +bool DihedralRestraint::operator==(const DihedralRestraint &other) const +{ + return atms == other.atms and + _phi0 == other._phi0 and + _kphi == other._kphi; +} + +bool DihedralRestraint::operator!=(const DihedralRestraint &other) const +{ + return not operator==(other); +} + +DihedralRestraints DihedralRestraint::operator+(const DihedralRestraint &other) const +{ + return DihedralRestraints(*this) + other; +} + +DihedralRestraints DihedralRestraint::operator+(const DihedralRestraints &other) const +{ + return DihedralRestraints(*this) + other; +} + +const char *DihedralRestraint::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); +} + +const char *DihedralRestraint::what() const +{ + return DihedralRestraint::typeName(); +} + +DihedralRestraint *DihedralRestraint::clone() const +{ + return new DihedralRestraint(*this); +} + +bool DihedralRestraint::isNull() const +{ + return atms.isEmpty(); +} + +QString DihedralRestraint::toString() const +{ + if (this->isNull()) + return QObject::tr("DihedralRestraint::null"); + else + { + QStringList a; + + for (const auto &atom : atms) + { + a.append(QString::number(atom)); + } + return QString("DihedralRestraint( [%1], phi0=%2, kphi=%3 )") + .arg(a.join(", ")) + .arg(_phi0.toString()) + .arg(_kphi.toString()); + } +} + +/** Return the force constant for the restraint */ +SireUnits::Dimension::HarmonicAngleConstant DihedralRestraint::kphi() const +{ + return this->_kphi; +} + +/** Return the equilibrium angle for the restraint */ +SireUnits::Dimension::Angle DihedralRestraint::phi0() const +{ + return this->_phi0; +} + +/** Return the atoms involved in the restraint */ +QVector DihedralRestraint::atoms() const +{ + return this->atms; +} + +/////// +/////// Implementation of DihedralRestraints +/////// + +/** Serialise to a binary datastream */ + +static const RegisterMetaType r_dihrests; + +QDataStream &operator<<(QDataStream &ds, const DihedralRestraints &dihrests) +{ + writeHeader(ds, r_dihrests, 1); + + SharedDataStream sds(ds); + + sds << dihrests.r + << static_cast(dihrests); + + return ds; +} + +/** Extract from a binary datastream */ +QDataStream &operator>>(QDataStream &ds, DihedralRestraints &dihrests) +{ + VersionID v = readHeader(ds, r_dihrests); + + if (v == 1) + { + SharedDataStream sds(ds); + + sds >> dihrests.r >> + static_cast(dihrests); + } + else + throw version_error(v, "1", r_dihrests, CODELOC); + + return ds; +} + +/** Null constructor */ +DihedralRestraints::DihedralRestraints() + : ConcreteProperty() +{ +} + +DihedralRestraints::DihedralRestraints(const QString &name) + : ConcreteProperty(name) +{ +} + +DihedralRestraints::DihedralRestraints(const DihedralRestraint &restraint) + : ConcreteProperty() +{ + if (not restraint.isNull()) + r.append(restraint); +} + +DihedralRestraints::DihedralRestraints(const QList &restraints) + : ConcreteProperty() +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } +} + +DihedralRestraints::DihedralRestraints(const QString &name, + const DihedralRestraint &restraint) + : ConcreteProperty(name) +{ + if (not restraint.isNull()) + r.append(restraint); +} + +DihedralRestraints::DihedralRestraints(const QString &name, + const QList &restraints) + : ConcreteProperty(name) +{ + for (const auto &restraint : restraints) + { + if (not restraint.isNull()) + r.append(restraint); + } +} + +DihedralRestraints::DihedralRestraints(const DihedralRestraints &other) + : ConcreteProperty(other), r(other.r) +{ +} + +/* Desctructor */ +DihedralRestraints::~DihedralRestraints() +{ +} + +DihedralRestraints &DihedralRestraints::operator=(const DihedralRestraints &other) +{ + if (this != &other) + { + Restraints::operator=(other); + r = other.r; + } + + return *this; +} + +bool DihedralRestraints::operator==(const DihedralRestraints &other) const +{ + return r == other.r and Restraints::operator==(other); +} + +bool DihedralRestraints::operator!=(const DihedralRestraints &other) const +{ + return not operator==(other); +} + +const char *DihedralRestraints::typeName() +{ + return QMetaType::typeName(qMetaTypeId()); +} + +const char *DihedralRestraints::what() const +{ + return DihedralRestraints::typeName(); +} + +DihedralRestraints *DihedralRestraints::clone() const +{ + return new DihedralRestraints(*this); +} + +QString DihedralRestraints::toString() const +{ + if (this->isEmpty()) + return QObject::tr("DihedralRestraints::null"); + + QStringList parts; + + const auto n = this->count(); + + if (n <= 10) + { + for (int i = 0; i < n; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + } + else + { + for (int i = 0; i < 5; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + + parts.append("..."); + + for (int i = n - 5; i < n; i++) + { + parts.append(QObject::tr("%1: %2").arg(i).arg(this->r.at(i).toString())); + } + } + + return QObject::tr("DihedralRestraints( name=%1, size=%2\n%3\n )") + .arg(this->name()) + .arg(n) + .arg(parts.join("\n")); +} + +/** Return whether or not this is empty */ +bool DihedralRestraints::isEmpty() const +{ + return this->r.isEmpty(); +} + +/** Return whether or not this is empty */ +bool DihedralRestraints::isNull() const +{ + return this->isEmpty(); +} + +/** Return the number of restraints */ +int DihedralRestraints::nRestraints() const +{ + return this->r.count(); +} + +/** Return the number of restraints */ +int DihedralRestraints::count() const +{ + return this->nRestraints(); +} + +/** Return the number of restraints */ +int DihedralRestraints::size() const +{ + return this->nRestraints(); +} + +/** Return the ith restraint */ +const DihedralRestraint &DihedralRestraints::at(int i) const +{ + i = SireID::Index(i).map(this->r.count()); + + return this->r.at(i); +} + +/** Return the ith restraint */ +const DihedralRestraint &DihedralRestraints::operator[](int i) const +{ + return this->at(i); +} + +/** Return all of the restraints */ +QList DihedralRestraints::restraints() const +{ + return this->r; +} + +/** Add a restraints onto the list */ +void DihedralRestraints::add(const DihedralRestraint &restraint) +{ + if (not restraint.isNull()) + r.append(restraint); +} + +/** Add a restraint onto the list */ +void DihedralRestraints::add(const DihedralRestraints &restraints) +{ + this->r += restraints.r; +} + +/** Add a restraint onto the list */ +DihedralRestraints &DihedralRestraints::operator+=(const DihedralRestraint &restraint) +{ + this->add(restraint); + return *this; +} + +/** Add a restraint onto the list */ +DihedralRestraints DihedralRestraints::operator+(const DihedralRestraint &restraint) const +{ + DihedralRestraints ret(*this); + ret += restraint; + return *this; +} + +/** Add restraints onto the list */ +DihedralRestraints &DihedralRestraints::operator+=(const DihedralRestraints &restraints) +{ + this->add(restraints); + return *this; +} + +/** Add restraints onto the list */ +DihedralRestraints DihedralRestraints::operator+(const DihedralRestraints &restraints) const +{ + DihedralRestraints ret(*this); + ret += restraints; + return *this; +} \ No newline at end of file diff --git a/corelib/src/libs/SireMM/dihedralrestraints.h b/corelib/src/libs/SireMM/dihedralrestraints.h new file mode 100644 index 000000000..f40b03e0c --- /dev/null +++ b/corelib/src/libs/SireMM/dihedralrestraints.h @@ -0,0 +1,177 @@ +/********************************************\ + * + * Sire - Molecular Simulation Framework + * + * Copyright (C) 2025 Christopher Woods + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * For full details of the license please see the COPYING file + * that should have come with this distribution. + * + * You can contact the authors at https://sire.openbiosim.org + * +\*********************************************/ + +#ifndef SIREMM_DIHEDRALRESTRAINTS_H +#define SIREMM_DIHEDRALRESTRAINTS_H + +#include "restraints.h" + +#include "SireUnits/dimensions.h" +#include "SireUnits/generalunit.h" + +SIRE_BEGIN_HEADER + +namespace SireMM +{ + class DihedralRestraint; + class DihedralRestraints; +} + +SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::DihedralRestraint &); +SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::DihedralRestraint &); + +SIREMM_EXPORT QDataStream &operator<<(QDataStream &, const SireMM::DihedralRestraints &); +SIREMM_EXPORT QDataStream &operator>>(QDataStream &, SireMM::DihedralRestraints &); + +namespace SireMM +{ + + /** This class represents a single torsion restraint between any four + * atoms in a system + * @author Christopher Woods + */ + class SIREMM_EXPORT DihedralRestraint + : public SireBase::ConcreteProperty + { + + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::DihedralRestraint &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::DihedralRestraint &); + + public: + DihedralRestraint(); + DihedralRestraint(const QList &atoms, + const SireUnits::Dimension::Angle &phi0, + const SireUnits::Dimension::HarmonicAngleConstant &kphi); + + DihedralRestraint(const DihedralRestraint &other); + + ~DihedralRestraint(); + + DihedralRestraint &operator=(const DihedralRestraint &other); + + bool operator==(const DihedralRestraint &other) const; + bool operator!=(const DihedralRestraint &other) const; + + DihedralRestraints operator+(const DihedralRestraint &other) const; + DihedralRestraints operator+(const DihedralRestraints &other) const; + + static const char *typeName(); + const char *what() const; + + DihedralRestraint *clone() const; + + QString toString() const; + + bool isNull() const; + + QVector atoms() const; + + SireUnits::Dimension::Angle phi0() const; + SireUnits::Dimension::HarmonicAngleConstant kphi() const; + + private: + /** Atoms involved in the angle restraint */ + QVector atms; + + /** Equilibrium angle */ + SireUnits::Dimension::Angle _phi0; + + /** Harmonic angle constant */ + SireUnits::Dimension::HarmonicAngleConstant _kphi; + }; + + /** This class represents a collection of angle restraints */ + class SIREMM_EXPORT DihedralRestraints + : public SireBase::ConcreteProperty + { + friend SIREMM_EXPORT QDataStream & ::operator<<(QDataStream &, const SireMM::DihedralRestraints &); + friend SIREMM_EXPORT QDataStream & ::operator>>(QDataStream &, SireMM::DihedralRestraints &); + + public: + DihedralRestraints(); + DihedralRestraints(const QString &name); + + DihedralRestraints(const DihedralRestraint &restraint); + DihedralRestraints(const QList &restraints); + + DihedralRestraints(const QString &name, + const DihedralRestraint &restraint); + + DihedralRestraints(const QString &name, + const QList &restraints); + + DihedralRestraints(const DihedralRestraints &other); + + ~DihedralRestraints(); + + DihedralRestraints &operator=(const DihedralRestraints &other); + + bool operator==(const DihedralRestraints &other) const; + bool operator!=(const DihedralRestraints &other) const; + + static const char *typeName(); + const char *what() const; + + DihedralRestraints *clone() const; + + QString toString() const; + + bool isEmpty() const; + bool isNull() const; + + int count() const; + int size() const; + int nRestraints() const; + + const DihedralRestraint &at(int i) const; + const DihedralRestraint &operator[](int i) const; + + QList restraints() const; + + void add(const DihedralRestraint &restraint); + void add(const DihedralRestraints &restraints); + + DihedralRestraints &operator+=(const DihedralRestraint &restraint); + DihedralRestraints &operator+=(const DihedralRestraints &restraints); + + DihedralRestraints operator+(const DihedralRestraint &restraint) const; + DihedralRestraints operator+(const DihedralRestraints &restraints) const; + + private: + /** List of restraints */ + QList r; + }; +} + +Q_DECLARE_METATYPE(SireMM::DihedralRestraint) +Q_DECLARE_METATYPE(SireMM::DihedralRestraints) + +SIRE_EXPOSE_CLASS(SireMM::DihedralRestraint) +SIRE_EXPOSE_CLASS(SireMM::DihedralRestraints) +SIRE_END_HEADER + +#endif \ No newline at end of file diff --git a/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 `__ 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/mm/__init__.py b/src/sire/mm/__init__.py index 2255fd1ac..9df64e5d1 100644 --- a/src/sire/mm/__init__.py +++ b/src/sire/mm/__init__.py @@ -4,6 +4,8 @@ "AmberDihPart", "AmberDihedral", "Angle", + "AngleRestraint", + "AngleRestraints", "Bond", "BondRestraint", "BondRestraints", @@ -11,6 +13,8 @@ "Improper", "PositionRestraint", "PositionRestraints", + "DihedralRestraint", + "DihedralRestraints", "SelectorAngle", "SelectorBond", "SelectorDihedral", @@ -29,6 +33,9 @@ _use_new_api() +AngleRestraint = _MM.AngleRestraint +AngleRestraints = _MM.AngleRestraints + # It would be better if these were called "DistanceRestraints", # but there is already a legacy Sire.MM class with this name BondRestraint = _MM.BondRestraint @@ -40,6 +47,9 @@ PositionalRestraint = _MM.PositionalRestraint PositionalRestraints = _MM.PositionalRestraints +DihedralRestraint = _MM.DihedralRestraint +DihedralRestraints = _MM.DihedralRestraints + AmberBond = _MM.AmberBond AmberAngle = _MM.AmberAngle AmberDihPart = _MM.AmberDihPart diff --git a/src/sire/restraints/__init__.py b/src/sire/restraints/__init__.py index 5f5af1465..166cb985c 100644 --- a/src/sire/restraints/__init__.py +++ b/src/sire/restraints/__init__.py @@ -1,4 +1,4 @@ -__all__ = ["positional", "bond", "distance", "boresch", "get_standard_state_correction"] +__all__ = ["angle", "positional", "bond", "dihedral", "distance", "boresch", "get_standard_state_correction"] -from ._restraints import bond, boresch, distance, positional +from ._restraints import angle, bond, boresch, dihedral, distance, positional from ._standard_state_correction import get_standard_state_correction diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 077cb97c0..6e4d5d050 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -1,6 +1,8 @@ __all__ = [ + "angle", "boresch", "bond", + "dihedral", "distance", "positional", ] @@ -18,6 +20,96 @@ def _to_atoms(mols, atoms): return selection_to_atoms(mols, atoms) +def angle(mols, atoms, theta0=None, ktheta=None, name=None, map=None): + """ + Create a set of anglel restraints from all of the atoms in 'atoms' + where all atoms are contained in the container 'mols', using the + passed values of the force constant 'ktheta' and equilibrium + angle value theta0. + + If theta0 is None, then the current angle for + provided atoms will be used as the equilibium value. + + If ktheta is None, then a default value of 100 kcal mol-1 rad-2 will be used + + Parameters + ---------- + mols : sire.system._system.System + The system containing the atoms. + + atoms : SireMol::Selector + The atoms to restrain. + + ktheta : str or SireUnits::Dimension::GeneralUnit or, optional + The force constants for the angle restraints. + If None, this will default to 100 kcal mol-1 rad-2. + Default is None. + + theta0 : str or SireUnits::Dimension::GeneralUnit, optional + The equilibrium angles for the angle restraints. If None, these + will be measured from the current coordinates of the atoms. + Default is None. + + Returns + ------- + AnglelRestraints : SireMM::AngleRestraints + A container of angle restraints, where the first restraint is + the AngleRestraint created. The angle restraint created can be + extracted with AngleRestraints[0]. + """ + from .. import u + from ..base import create_map + from ..mm import AngleRestraint, AngleRestraints + + map = create_map(map) + map_dict = map.to_dict() + ktheta = ktheta if ktheta is not None else map_dict.get("ktheta", None) + theta0 = theta0 if theta0 is not None else map_dict.get("theta0", None) + name = name if name is not None else map_dict.get("name", None) + + atoms = _to_atoms(mols, atoms) + + if len(atoms) != 3: + raise ValueError( + "You need to provide 3 atoms to create an angle restraint" + f"whereas {len(atoms)} atoms were provided." + ) + + from .. import measure + + if ktheta is None: + ktheta = u("100 kcal mol-1 rad-2") + + elif type(ktheta) is list: + raise NotImplementedError( + "Setup of multiple angle restraints simultaneously is not currently supported. " + "Please set up each restraint individually and then combine them into multiple restraints." + ) + + if theta0 is None: + from .. import measure + + theta0 = measure(atoms[0], atoms[1], atoms[2]) + + elif type(theta0) is list: + raise NotImplementedError( + "Setup of multiple angle restraints simultaneously is not currently supported. " + "Please set up each restraint individually and then combine them into multiple restraints." + ) + else: + theta0 = u(theta0) + + mols = mols.atoms() + + if name is None: + restraints = AngleRestraints() + else: + restraints = AngleRestraints(name=name) + + restraints.add(AngleRestraint(mols.find(atoms), theta0, ktheta)) + return restraints + + def boresch( mols, receptor, @@ -377,6 +469,96 @@ def _check_stability_boresch_restraint(restraint_components, temperature=u("298 ) +def dihedral(mols, atoms, phi0=None, kphi=None, name=None, map=None): + """ + Create a set of dihedral restraints from all of the atoms in 'atoms' + where all atoms are contained in the container 'mols', using the + passed values of the force constant 'kphi' and equilibrium + torsion angle phi0. + + If phi0 is None, then the current torsional angle for + the provided atoms will be used as the equilibium value. + + If kphi is None, then a default value of 100 kcal mol-1 rad-2 will be used + + Parameters + ---------- + mols : sire.system._system.System + The system containing the atoms. + + atoms : SireMol::Selector + The atoms to restrain. + + kphi : str or SireUnits::Dimension::GeneralUnit or, optional + The force constants for the torsion restraints. + If None, this will default to 100 kcal mol-1 rad-2. + Default is None. + + phi0 : str or SireUnits::Dimension::GeneralUnit, optional + The equilibrium torsional angle for restraints. If None, these + will be measured from the current coordinates of the atoms. + Default is None. + + Returns + ------- + DihedralRestraints : SireMM::DihedralRestraints + A container of Dihedral restraints, where the first restraint is + the DihedralRestraint created. The Dihedral restraint created can be + extracted with DihedralRestraints[0]. + """ + from .. import u + from ..base import create_map + from ..mm import DihedralRestraint, DihedralRestraints + + map = create_map(map) + map_dict = map.to_dict() + kphi = kphi if kphi is not None else map_dict.get("kphi", None) + phi0 = phi0 if phi0 is not None else map_dict.get("phi0", None) + name = name if name is not None else map_dict.get("name", None) + + atoms = _to_atoms(mols, atoms) + + if len(atoms) != 4: + raise ValueError( + "You need to provide 4 atoms to create a dihedral restraint" + f"whereas {len(atoms)} atoms were provided." + ) + + from .. import measure + + if kphi is None: + kphi = u("100 kcal mol-1 rad-2") + + elif type(kphi) is list: + raise NotImplementedError( + "Setup of multiple dihedral restraints simultaneously is not currently supported. " + "Please set up each restraint individually and then combine them into multiple restraints." + ) + + if phi0 is None: + from .. import measure + + phi0 = measure(atoms[0], atoms[1], atoms[2], atoms[3]) + + elif type(phi0) is list: + raise NotImplementedError( + "Setup of multiple dihedral restraints simultaneously is not currently supported. " + "Please set up each restraint individually and then combine them into multiple restraints." + ) + else: + phi0 = u(phi0) + + mols = mols.atoms() + + if name is None: + restraints = DihedralRestraints() + else: + restraints = DihedralRestraints(name=name) + + restraints.add(DihedralRestraint(mols.find(atoms), phi0, kphi)) + return restraints + + def distance(mols, atoms0, atoms1, r0=None, k=None, name=None, map=None): """ Create a set of distance restraints from all of the atoms in 'atoms0' diff --git a/tests/restraints/test_angle_dihedral_restraints.py b/tests/restraints/test_angle_dihedral_restraints.py new file mode 100644 index 000000000..f6096ec21 --- /dev/null +++ b/tests/restraints/test_angle_dihedral_restraints.py @@ -0,0 +1,105 @@ +import pytest +import sire as sr + + +def test_default_angle_restraints_setup(ala_mols): + """Tests that angle restraints can be set up correctly with default parameters.""" + mols = ala_mols.clone() + restraints = sr.restraints.angle(mols=mols, atoms=[0, 1, 2]) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2] + assert restraints[0].ktheta().value() == 100.0 + + +def test_angle_restraint_custom_ktheta(ala_mols): + """Tests that angle restraints can be set up correctly with custom ktheta.""" + mols = ala_mols.clone() + restraints = sr.restraints.angle( + mols=mols, atoms=[0, 1, 2], ktheta="10 kcal mol-1 rad-2" + ) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2] + assert restraints[0].ktheta().value() == 10.0 + + +def test_angle_restraint_custom_ktheta_and_theta0(ala_mols): + """Tests that angle restraints can be set up correctly with custom ktheta and theta0.""" + mols = ala_mols.clone() + restraints = sr.restraints.angle( + mols=mols, atoms=[0, 1, 2], ktheta="10 kcal mol-1 rad-2", theta0="2 rad" + ) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2] + assert restraints[0].ktheta().value() == 10.0 + assert restraints[0].theta0().value() == 2.0 + + +def test_angle_restraint_molecular_container(ala_mols): + """Tests that angle restraints can be set up correctly with a molecular container (i.e. passing the angle of the molecule).""" + mols = ala_mols.clone() + ang = mols.angles()[0] + restraints = sr.restraints.angle(mols=mols, atoms=ang.atoms()) + + +def test_default_dihedral_restraints_setup(ala_mols): + """Tests that dihedral restraints can be set up correctly with default parameters.""" + mols = ala_mols.clone() + restraints = sr.restraints.dihedral(mols=mols, atoms=[0, 1, 2, 3]) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2, 3] + assert restraints[0].kphi().value() == 100.0 + + +def test_dihedral_restraint_custom_kphi(ala_mols): + """Tests that dihedral restraints can be set up correctly with custom kphi.""" + mols = ala_mols.clone() + restraints = sr.restraints.dihedral( + mols=mols, atoms=[0, 1, 2, 3], kphi="10 kcal mol-1 rad-2" + ) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2, 3] + assert restraints[0].kphi().value() == 10.0 + + +def test_dihedral_restraint_custom_kphi_and_phi0(ala_mols): + """Tests that dihedral restraints can be set up correctly with custom kphi and phi0.""" + mols = ala_mols.clone() + restraints = sr.restraints.dihedral( + mols=mols, atoms=[0, 1, 2, 3], kphi="10 kcal mol-1 rad-2", phi0="2 rad" + ) + assert restraints.num_restraints() == 1 + assert restraints[0].atoms() == [0, 1, 2, 3] + assert restraints[0].kphi().value() == 10.0 + assert restraints[0].phi0().value() == 2.0 + + +def test_dihedral_restraint_molecular_container(ala_mols): + """Tests that dihedral restraints can be set up correctly with a molecular container (i.e. passing the dihedral of the molecule).""" + mols = ala_mols.clone() + dih = mols.dihedrals()[0] + restraints = sr.restraints.dihedral(mols=mols, atoms=dih.atoms()) + + +def test_multiple_angle_restraints(ala_mols): + """Tests that multiple angle restraints can be set up correctly.""" + mols = ala_mols.clone() + ang0 = mols.angles()[0] + ang1 = mols.angles()[1] + restraints = sr.restraints.angle(mols=mols, atoms=ang0.atoms()) + restraint1 = sr.restraints.angle(mols=mols, atoms=ang1.atoms()) + restraints.add(restraint1) + assert restraints.num_restraints() == 2 + + +def test_multiple_dihedral_restraints(ala_mols): + """Tests that multiple dihedral restraints can be set up correctly.""" + mols = ala_mols.clone() + dih0 = mols.dihedrals()[0] + dih1 = mols.dihedrals()[1] + dih2 = mols.dihedrals()[2] + restraints = sr.restraints.dihedral(mols=mols, atoms=dih0.atoms()) + restraint1 = sr.restraints.dihedral(mols=mols, atoms=dih1.atoms()) + restraint2 = sr.restraints.dihedral(mols=mols, atoms=dih2.atoms()) + restraints.add(restraint1) + restraints.add(restraint2) + assert restraints.num_restraints() == 3 diff --git a/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp b/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp index 54e828a53..517aa0ee4 100644 --- a/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp @@ -205,7 +205,7 @@ void register_LambdaLever_class(){ } { //::SireOpenMM::LambdaLever::setExceptionIndicies - typedef void ( ::SireOpenMM::LambdaLever::*setExceptionIndicies_function_type)( int,::QString const &,::QVector< boost::tuples::tuple< int, int, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > > const & ) ; + typedef void ( ::SireOpenMM::LambdaLever::*setExceptionIndicies_function_type)( int,::QString const &,::QVector< boost::tuples::tuple< int, int > > const & ) ; setExceptionIndicies_function_type setExceptionIndicies_function_value( &::SireOpenMM::LambdaLever::setExceptionIndicies ); LambdaLever_exposer.def( diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 67707a398..a05af04dc 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -26,11 +26,11 @@ * \*********************************************/ -#include "pyqm.h" #include "lambdalever.h" +#include "pyqm.h" -#include "SireBase/propertymap.h" #include "SireBase/arrayproperty.hpp" +#include "SireBase/propertymap.h" #include "SireCAS/values.h" @@ -1857,6 +1857,100 @@ void _update_restraint_in_context(OpenMM::CustomCompoundBondForce *ff, double rh ff->updateParametersInContext(context); } +/** Update the parameters for a CustomAngleForce for scale factor 'rho' + * in the passed context */ +void _update_restraint_in_context(OpenMM::CustomAngleForce *ff, double rho, + OpenMM::Context &context) +{ + if (ff == 0) + throw SireError::invalid_cast(QObject::tr( + "Unable to cast the restraint force to an OpenMM::CustomAngleForce, " + "despite it reporting that is was an object of this type..."), + CODELOC); + + const int nangles = ff->getNumAngles(); + + if (nangles == 0) + // nothing to update + return; + + const int nparams = ff->getNumPerAngleParameters(); + + if (nparams == 0) + throw SireError::incompatible_error(QObject::tr( + "Unable to set 'rho' for this restraint as it has no custom parameters!"), + CODELOC); + + // we set the first parameter - we can see what the current value + // is from the first restraint. This is because rho should be the + // first parameter and have the same value for all restraints + std::vector custom_params; + custom_params.resize(nparams); + int atom0, atom1, atom2; + + ff->getAngleParameters(0, atom0, atom1, atom2, custom_params); + + if (custom_params[0] == rho) + // nothing to do - it is already equal to this value + return; + + for (int i = 0; i < nangles; ++i) + { + ff->getAngleParameters(i, atom0, atom1, atom2, custom_params); + custom_params[0] = rho; + ff->setAngleParameters(i, atom0, atom1, atom2, custom_params); + } + + ff->updateParametersInContext(context); +} + +/** Update the parameters for a CustomTorsionForce for scale factor 'rho' + * in the passed context */ +void _update_restraint_in_context(OpenMM::CustomTorsionForce *ff, double rho, + OpenMM::Context &context) +{ + if (ff == 0) + throw SireError::invalid_cast(QObject::tr( + "Unable to cast the restraint force to an OpenMM::CustomTorsionForce, " + "despite it reporting that is was an object of this type..."), + CODELOC); + + const int ntorsions = ff->getNumTorsions(); + + if (ntorsions == 0) + // nothing to update + return; + + const int nparams = ff->getNumPerTorsionParameters(); + + if (nparams == 0) + throw SireError::incompatible_error(QObject::tr( + "Unable to set 'rho' for this restraint as it has no custom parameters!"), + CODELOC); + + // we set the first parameter - we can see what the current value + // is from the first restraint. This is because rho should be the + // first parameter and have the same value for all restraints + std::vector custom_params; + custom_params.resize(nparams); + int atom0, atom1, atom2, atom3; + + ff->getTorsionParameters(0, atom0, atom1, atom2, atom3, custom_params); + + if (custom_params[0] == rho) + // nothing to do - it is already equal to this value + return; + + for (int i = 0; i < ntorsions; ++i) + { + ff->getTorsionParameters(i, atom0, atom1, atom2, atom3, custom_params); + custom_params[0] = rho; + ff->setTorsionParameters(i, atom0, atom1, atom2, atom3, custom_params); + } + + ff->updateParametersInContext(context); +} + /** Update the parameters for a CustomBondForce for scale factor 'rho' * in the passed context */ void _update_restraint_in_context(OpenMM::CustomBondForce *ff, double rho, @@ -1918,6 +2012,24 @@ void LambdaLever::updateRestraintInContext(OpenMM::Force &ff, double rho, dynamic_cast(&ff), rho, context); } + else if (ff_type == "AngleRestraintForce") + { + _update_restraint_in_context( + dynamic_cast(&ff), + rho, context); + } + else if (ff_type == "TorsionRestraintForce") + { + _update_restraint_in_context( + dynamic_cast(&ff), + rho, context); + } + else if (ff_type == "CustomCompoundBondForce") + { + _update_restraint_in_context( + dynamic_cast(&ff), + rho, context); + } else if (ff_type == "BoreschRestraintForce") { _update_restraint_in_context( @@ -2007,7 +2119,6 @@ int LambdaLever::addPerturbableMolecule(const OpenMMMolecule &molecule, return this->perturbable_mols.count() - 1; } - /** Set the exception indices for the perturbable molecule at * index 'mol_idx' */ diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index e3c46a03e..4c5b928d9 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -8,37 +8,39 @@ #include "SireSystem/forcefieldinfo.h" -#include "SireMol/core.h" -#include "SireMol/moleditor.h" -#include "SireMol/atomelements.h" #include "SireMol/atomcharges.h" #include "SireMol/atomcoords.h" +#include "SireMol/atomelements.h" #include "SireMol/atommasses.h" #include "SireMol/atomproperty.hpp" -#include "SireMol/connectivity.h" +#include "SireMol/atomvelocities.h" #include "SireMol/bondid.h" #include "SireMol/bondorder.h" -#include "SireMol/atomvelocities.h" +#include "SireMol/connectivity.h" +#include "SireMol/core.h" +#include "SireMol/moleditor.h" -#include "SireMM/atomljs.h" -#include "SireMM/selectorbond.h" #include "SireMM/amberparams.h" +#include "SireMM/anglerestraints.h" +#include "SireMM/atomljs.h" #include "SireMM/bondrestraints.h" -#include "SireMM/positionalrestraints.h" #include "SireMM/boreschrestraints.h" +#include "SireMM/dihedralrestraints.h" +#include "SireMM/positionalrestraints.h" +#include "SireMM/selectorbond.h" #include "SireVol/periodicbox.h" #include "SireVol/triclinicbox.h" #include "SireCAS/lambdaschedule.h" -#include "SireMaths/vector.h" #include "SireMaths/maths.h" +#include "SireMaths/vector.h" +#include "SireBase/generalunitproperty.h" +#include "SireBase/lengthproperty.h" #include "SireBase/parallel.h" #include "SireBase/propertylist.h" -#include "SireBase/lengthproperty.h" -#include "SireBase/generalunitproperty.h" #include "SireUnits/units.h" @@ -398,6 +400,132 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, } } +/** Add all of the angle restraints from 'restraints' to the passed + * system, which is acted on by the passed LambdaLever. The number + * of real (non-anchor) atoms in the OpenMM::System is 'natoms' + */ + +void _add_angle_restraints(const SireMM::AngleRestraints &restraints, + OpenMM::System &system, LambdaLever &lambda_lever, + int natoms) +{ + if (restraints.isEmpty()) + return; + + // energy expression of the angle restraint, which acts over three atoms + // + // theta = angle(P1, P2, P3) + // + // The energies are + // + // e_restraint = rho * (e_angle) + // e_angle = ktheta(theta - theta0)^2 + + const auto energy_expression = QString( + "rho*k*(theta-theta0)^2;") + .toStdString(); + + auto *restraintff = new OpenMM::CustomAngleForce(energy_expression); + + restraintff->setName("AngleRestraintForce"); + restraintff->addPerAngleParameter("rho"); + restraintff->addPerAngleParameter("k"); + restraintff->addPerAngleParameter("theta0"); + + restraintff->setUsesPeriodicBoundaryConditions(true); + + lambda_lever.addRestraintIndex(restraints.name(), + system.addForce(restraintff)); + + const double internal_to_ktheta = (1 * SireUnits::kcal_per_mol / (SireUnits::radian2)).to(SireUnits::kJ_per_mol / SireUnits::radian2); + + const auto atom_restraints = restraints.restraints(); + + for (const auto &restraint : atom_restraints) + { + std::vector particles; + particles.resize(3); + + for (int i = 0; i < 3; ++i) + { + particles[i] = restraint.atoms()[i]; + } + + std::vector parameters; + parameters.resize(3); + + parameters[0] = 1.0; // rho + parameters[1] = restraint.ktheta().value() * internal_to_ktheta; // k + parameters[2] = restraint.theta0().value(); // theta0 (already in radians) + + // restraintff->addTorsion(particles, parameters); + restraintff->addAngle(particles[0], particles[1], particles[2], parameters); + } +} + +void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, + OpenMM::System &system, LambdaLever &lambda_lever, + int natoms) +{ + if (restraints.isEmpty()) + return; + + // energy expression of the dihedral restraint, which acts over four atoms + // + // phi = dihedral(P1, P2, P3, P4) + // + // The energies are + // + // e_restraint = rho * (e_torsion) + // e_torsion = k_phi(dphi)^2 where + // dphi = abs(phi - phi0) + + const auto energy_expression = QString( + "rho*k*min(dtheta, two_pi-dtheta)^2;" + "dtheta = abs(theta-theta0);" + "two_pi=6.283185307179586;") + .toStdString(); + + auto *restraintff = new OpenMM::CustomTorsionForce(energy_expression); + + // it seems that OpenMM wants to call the torsion angle theta rather than phi + // we need to rename our parameters accordingly + restraintff->setName("TorsionRestraintForce"); + restraintff->addPerTorsionParameter("rho"); + restraintff->addPerTorsionParameter("k"); + restraintff->addPerTorsionParameter("theta0"); + + restraintff->setUsesPeriodicBoundaryConditions(true); + + lambda_lever.addRestraintIndex(restraints.name(), + system.addForce(restraintff)); + + const double internal_to_ktheta = (1 * SireUnits::kcal_per_mol / (SireUnits::radian2)).to(SireUnits::kJ_per_mol / SireUnits::radian2); + + const auto atom_restraints = restraints.restraints(); + + for (const auto &restraint : atom_restraints) + { + std::vector particles; + particles.resize(4); + + for (int i = 0; i < 4; ++i) + { + particles[i] = restraint.atoms()[i]; + } + + std::vector parameters; + parameters.resize(3); + + parameters[0] = 1.0; // rho + parameters[1] = restraint.kphi().value() * internal_to_ktheta; // k + parameters[2] = restraint.phi0().value(); // phi0 (already in radians) --> theta0 + + // restraintff->addTorsion(particles, parameters); + restraintff->addTorsion(particles[0], particles[1], particles[2], particles[3], parameters); + } +} + /** Set the coulomb and LJ cutoff in the passed NonbondedForce, * based on the information in the passed ForceFieldInfo. * This sets the cutoff type (e.g. PME) and the actual @@ -1682,7 +1810,17 @@ 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_angle_restraints(prop.read().asA(), + system, lambda_lever, start_index); + } + else if (prop.read().isA()) { _add_positional_restraints(prop.read().asA(), system, lambda_lever, anchor_coords, start_index); diff --git a/wrapper/MM/AngleRestraint.pypp.cpp b/wrapper/MM/AngleRestraint.pypp.cpp index 052916bb4..b90320621 100644 --- a/wrapper/MM/AngleRestraint.pypp.cpp +++ b/wrapper/MM/AngleRestraint.pypp.cpp @@ -3,36 +3,25 @@ // (C) Christopher Woods, GPL >= 3 License #include "boost/python.hpp" -#include "Helpers/clone_const_reference.hpp" #include "AngleRestraint.pypp.hpp" namespace bp = boost::python; -#include "SireCAS/conditional.h" - #include "SireCAS/errors.h" -#include "SireCAS/power.h" - -#include "SireCAS/symbols.h" - -#include "SireCAS/values.h" - -#include "SireFF/forcetable.h" - #include "SireID/index.h" #include "SireStream/datastream.h" #include "SireStream/shareddatastream.h" -#include "SireUnits/angle.h" - #include "SireUnits/units.h" -#include "anglerestraint.h" +#include "anglerestraints.h" + +#include -#include "anglerestraint.h" +#include "anglerestraints.h" SireMM::AngleRestraint __copy__(const SireMM::AngleRestraint &other){ return SireMM::AngleRestraint(other); } @@ -47,162 +36,50 @@ SireMM::AngleRestraint __copy__(const SireMM::AngleRestraint &other){ return Sir void register_AngleRestraint_class(){ { //::SireMM::AngleRestraint - typedef bp::class_< SireMM::AngleRestraint, bp::bases< SireMM::Restraint3D, SireMM::Restraint, SireBase::Property > > AngleRestraint_exposer_t; - AngleRestraint_exposer_t AngleRestraint_exposer = AngleRestraint_exposer_t( "AngleRestraint", "This is a restraint that operates on the angle between\nthree SireMM::Point objects (e.g. three atoms in a molecule)\n\nAuthor: Christopher Woods\n", bp::init< >("Constructor") ); + typedef bp::class_< SireMM::AngleRestraint, bp::bases< SireBase::Property > > AngleRestraint_exposer_t; + AngleRestraint_exposer_t AngleRestraint_exposer = AngleRestraint_exposer_t( "AngleRestraint", "This class represents a single angle restraint between any three\natoms in a system\nAuthor: Christopher Woods\n", bp::init< >("Null constructor") ); bp::scope AngleRestraint_scope( AngleRestraint_exposer ); - AngleRestraint_exposer.def( bp::init< SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireCAS::Expression const & >(( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("restraint") ), "Construct a restraint that acts on the angle within the\nthree points point0, point1 and point2 (theta == a(012)),\nrestraining the angle within these points using the expression\nrestraint") ); - AngleRestraint_exposer.def( bp::init< SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireCAS::Expression const &, SireCAS::Values const & >(( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("restraint"), bp::arg("values") ), "Construct a restraint that acts on the angle within the\nthree points point0, point1 and point2 (theta == a(012)),\nrestraining the angle within these points using the expression\nrestraint") ); + AngleRestraint_exposer.def( bp::init< QList< long long > const &, SireUnits::Dimension::Angle const &, SireUnits::Dimension::HarmonicAngleConstant const & >(( bp::arg("atoms"), bp::arg("theta0"), bp::arg("ktheta") ), "Construct a restraint that acts on the angle within the\nthree atoms atom0, atom1 and atom2 (theta == a(012)),\nrestraining the angle within these atoms") ); AngleRestraint_exposer.def( bp::init< SireMM::AngleRestraint const & >(( bp::arg("other") ), "Copy constructor") ); - { //::SireMM::AngleRestraint::builtinSymbols - - typedef ::SireCAS::Symbols ( ::SireMM::AngleRestraint::*builtinSymbols_function_type)( ) const; - builtinSymbols_function_type builtinSymbols_function_value( &::SireMM::AngleRestraint::builtinSymbols ); - - AngleRestraint_exposer.def( - "builtinSymbols" - , builtinSymbols_function_value - , bp::release_gil_policy() - , "Return the built-in symbols of this restraint" ); - - } - { //::SireMM::AngleRestraint::builtinValues - - typedef ::SireCAS::Values ( ::SireMM::AngleRestraint::*builtinValues_function_type)( ) const; - builtinValues_function_type builtinValues_function_value( &::SireMM::AngleRestraint::builtinValues ); - - AngleRestraint_exposer.def( - "builtinValues" - , builtinValues_function_value - , bp::release_gil_policy() - , "Return the values of the built-in symbols of this restraint" ); - - } - { //::SireMM::AngleRestraint::contains - - typedef bool ( ::SireMM::AngleRestraint::*contains_function_type)( ::SireMol::MolNum ) const; - contains_function_type contains_function_value( &::SireMM::AngleRestraint::contains ); - - AngleRestraint_exposer.def( - "contains" - , contains_function_value - , ( bp::arg("molnum") ) - , bp::release_gil_policy() - , "Return whether or not this restraint affects the molecule\nwith number molnum" ); - - } - { //::SireMM::AngleRestraint::contains - - typedef bool ( ::SireMM::AngleRestraint::*contains_function_type)( ::SireMol::MolID const & ) const; - contains_function_type contains_function_value( &::SireMM::AngleRestraint::contains ); - - AngleRestraint_exposer.def( - "contains" - , contains_function_value - , ( bp::arg("molid") ) - , bp::release_gil_policy() - , "Return whether or not this restraint affects the molecule\nwith ID molid" ); - - } - { //::SireMM::AngleRestraint::differentialRestraintFunction - - typedef ::SireCAS::Expression const & ( ::SireMM::AngleRestraint::*differentialRestraintFunction_function_type)( ) const; - differentialRestraintFunction_function_type differentialRestraintFunction_function_value( &::SireMM::AngleRestraint::differentialRestraintFunction ); - - AngleRestraint_exposer.def( - "differentialRestraintFunction" - , differentialRestraintFunction_function_value - , bp::return_value_policy< bp::copy_const_reference >() - , "Return the function used to calculate the restraint force" ); - - } - { //::SireMM::AngleRestraint::differentiate - - typedef ::SireMM::RestraintPtr ( ::SireMM::AngleRestraint::*differentiate_function_type)( ::SireCAS::Symbol const & ) const; - differentiate_function_type differentiate_function_value( &::SireMM::AngleRestraint::differentiate ); - - AngleRestraint_exposer.def( - "differentiate" - , differentiate_function_value - , ( bp::arg("symbol") ) - , bp::release_gil_policy() - , "Return the differential of this restraint with respect to\nthe symbol symbol\nThrow: SireCAS::unavailable_differential\n" ); - - } - { //::SireMM::AngleRestraint::force - - typedef void ( ::SireMM::AngleRestraint::*force_function_type)( ::SireFF::MolForceTable &,double ) const; - force_function_type force_function_value( &::SireMM::AngleRestraint::force ); - - AngleRestraint_exposer.def( - "force" - , force_function_value - , ( bp::arg("forcetable"), bp::arg("scale_force")=1 ) - , "Calculate the force acting on the molecule in the forcetable forcetable\ncaused by this restraint, and add it on to the forcetable scaled by\nscale_force" ); - - } - { //::SireMM::AngleRestraint::force - - typedef void ( ::SireMM::AngleRestraint::*force_function_type)( ::SireFF::ForceTable &,double ) const; - force_function_type force_function_value( &::SireMM::AngleRestraint::force ); - - AngleRestraint_exposer.def( - "force" - , force_function_value - , ( bp::arg("forcetable"), bp::arg("scale_force")=1 ) - , "Calculate the force acting on the molecules in the forcetable forcetable\ncaused by this restraint, and add it on to the forcetable scaled by\nscale_force" ); + { //::SireMM::AngleRestraint::atoms - } - { //::SireMM::AngleRestraint::halfHarmonic - - typedef ::SireMM::AngleRestraint ( *halfHarmonic_function_type )( ::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireUnits::Dimension::Angle const &,::SireMM::HarmonicAngleForceConstant const & ); - halfHarmonic_function_type halfHarmonic_function_value( &::SireMM::AngleRestraint::halfHarmonic ); + typedef ::QVector< long long > ( ::SireMM::AngleRestraint::*atoms_function_type)( ) const; + atoms_function_type atoms_function_value( &::SireMM::AngleRestraint::atoms ); AngleRestraint_exposer.def( - "halfHarmonic" - , halfHarmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("angle"), bp::arg("force_constant") ) + "atoms" + , atoms_function_value , bp::release_gil_policy() - , "Return a distance restraint that applied a half-harmonic potential\nbetween the points point0 and point1 above a distance distance\nusing a force constant force_constant" ); + , "Return the atoms involved in the restraint" ); } - { //::SireMM::AngleRestraint::harmonic + { //::SireMM::AngleRestraint::isNull - typedef ::SireMM::AngleRestraint ( *harmonic_function_type )( ::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireMM::HarmonicAngleForceConstant const & ); - harmonic_function_type harmonic_function_value( &::SireMM::AngleRestraint::harmonic ); + typedef bool ( ::SireMM::AngleRestraint::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::AngleRestraint::isNull ); AngleRestraint_exposer.def( - "harmonic" - , harmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("force_constant") ) + "isNull" + , isNull_function_value , bp::release_gil_policy() - , "Return a distance restraint that applies a harmonic potential between\nthe points point0 and point1 using a force constant force_constant" ); - - } - { //::SireMM::AngleRestraint::molecules - - typedef ::SireMol::Molecules ( ::SireMM::AngleRestraint::*molecules_function_type)( ) const; - molecules_function_type molecules_function_value( &::SireMM::AngleRestraint::molecules ); - - AngleRestraint_exposer.def( - "molecules" - , molecules_function_value - , bp::release_gil_policy() - , "Return the molecules used in this restraint" ); + , "" ); } - { //::SireMM::AngleRestraint::nPoints + { //::SireMM::AngleRestraint::ktheta - typedef int ( ::SireMM::AngleRestraint::*nPoints_function_type)( ) const; - nPoints_function_type nPoints_function_value( &::SireMM::AngleRestraint::nPoints ); + typedef ::SireUnits::Dimension::HarmonicAngleConstant ( ::SireMM::AngleRestraint::*ktheta_function_type)( ) const; + ktheta_function_type ktheta_function_value( &::SireMM::AngleRestraint::ktheta ); AngleRestraint_exposer.def( - "nPoints" - , nPoints_function_value + "ktheta" + , ktheta_function_value , bp::release_gil_policy() - , "This restraint involves three points" ); + , "Return the force constant for the restraint" ); } AngleRestraint_exposer.def( bp::self != bp::self ); + AngleRestraint_exposer.def( bp::self + bp::self ); + AngleRestraint_exposer.def( bp::self + bp::other< SireMM::AngleRestraints >() ); { //::SireMM::AngleRestraint::operator= typedef ::SireMM::AngleRestraint & ( ::SireMM::AngleRestraint::*assign_function_type)( ::SireMM::AngleRestraint const & ) ; @@ -217,78 +94,28 @@ void register_AngleRestraint_class(){ } AngleRestraint_exposer.def( bp::self == bp::self ); - { //::SireMM::AngleRestraint::point - - typedef ::SireFF::Point const & ( ::SireMM::AngleRestraint::*point_function_type)( int ) const; - point_function_type point_function_value( &::SireMM::AngleRestraint::point ); - - AngleRestraint_exposer.def( - "point" - , point_function_value - , ( bp::arg("i") ) - , bp::return_value_policy() - , "Return the ith point" ); - - } - { //::SireMM::AngleRestraint::point0 - - typedef ::SireFF::Point const & ( ::SireMM::AngleRestraint::*point0_function_type)( ) const; - point0_function_type point0_function_value( &::SireMM::AngleRestraint::point0 ); - - AngleRestraint_exposer.def( - "point0" - , point0_function_value - , bp::return_value_policy() - , "Return the first point" ); - - } - { //::SireMM::AngleRestraint::point1 - - typedef ::SireFF::Point const & ( ::SireMM::AngleRestraint::*point1_function_type)( ) const; - point1_function_type point1_function_value( &::SireMM::AngleRestraint::point1 ); - - AngleRestraint_exposer.def( - "point1" - , point1_function_value - , bp::return_value_policy() - , "Return the second point" ); - - } - { //::SireMM::AngleRestraint::point2 - - typedef ::SireFF::Point const & ( ::SireMM::AngleRestraint::*point2_function_type)( ) const; - point2_function_type point2_function_value( &::SireMM::AngleRestraint::point2 ); - - AngleRestraint_exposer.def( - "point2" - , point2_function_value - , bp::return_value_policy() - , "Return the third point" ); - - } - { //::SireMM::AngleRestraint::setSpace + { //::SireMM::AngleRestraint::theta0 - typedef void ( ::SireMM::AngleRestraint::*setSpace_function_type)( ::SireVol::Space const & ) ; - setSpace_function_type setSpace_function_value( &::SireMM::AngleRestraint::setSpace ); + typedef ::SireUnits::Dimension::Angle ( ::SireMM::AngleRestraint::*theta0_function_type)( ) const; + theta0_function_type theta0_function_value( &::SireMM::AngleRestraint::theta0 ); AngleRestraint_exposer.def( - "setSpace" - , setSpace_function_value - , ( bp::arg("space") ) + "theta0" + , theta0_function_value , bp::release_gil_policy() - , "Set the space used to evaluate the energy of this restraint\nThrow: SireVol::incompatible_space\n" ); + , "Return the equilibrium angle for the restraint" ); } - { //::SireMM::AngleRestraint::theta + { //::SireMM::AngleRestraint::toString - typedef ::SireCAS::Symbol const & ( *theta_function_type )( ); - theta_function_type theta_function_value( &::SireMM::AngleRestraint::theta ); + typedef ::QString ( ::SireMM::AngleRestraint::*toString_function_type)( ) const; + toString_function_type toString_function_value( &::SireMM::AngleRestraint::toString ); AngleRestraint_exposer.def( - "theta" - , theta_function_value - , bp::return_value_policy() - , "Return the symbol that represents the angle between the points (theta)" ); + "toString" + , toString_function_value + , bp::release_gil_policy() + , "" ); } { //::SireMM::AngleRestraint::typeName @@ -303,69 +130,26 @@ void register_AngleRestraint_class(){ , "" ); } - { //::SireMM::AngleRestraint::update + { //::SireMM::AngleRestraint::what - typedef void ( ::SireMM::AngleRestraint::*update_function_type)( ::SireMol::MoleculeData const & ) ; - update_function_type update_function_value( &::SireMM::AngleRestraint::update ); + typedef char const * ( ::SireMM::AngleRestraint::*what_function_type)( ) const; + what_function_type what_function_value( &::SireMM::AngleRestraint::what ); AngleRestraint_exposer.def( - "update" - , update_function_value - , ( bp::arg("moldata") ) + "what" + , what_function_value , bp::release_gil_policy() - , "Update the points of this restraint using new molecule data from moldata\nThrow: SireBase::missing_property\nThrow: SireError::invalid_cast\nThrow: SireError::incompatible_error\n" ); - - } - { //::SireMM::AngleRestraint::update - - typedef void ( ::SireMM::AngleRestraint::*update_function_type)( ::SireMol::Molecules const & ) ; - update_function_type update_function_value( &::SireMM::AngleRestraint::update ); - - AngleRestraint_exposer.def( - "update" - , update_function_value - , ( bp::arg("molecules") ) - , bp::release_gil_policy() - , "Update the points of this restraint using new molecule data from molecules\nThrow: SireBase::missing_property\nThrow: SireError::invalid_cast\nThrow: SireError::incompatible_error\n" ); - - } - { //::SireMM::AngleRestraint::usesMoleculesIn - - typedef bool ( ::SireMM::AngleRestraint::*usesMoleculesIn_function_type)( ::SireFF::ForceTable const & ) const; - usesMoleculesIn_function_type usesMoleculesIn_function_value( &::SireMM::AngleRestraint::usesMoleculesIn ); - - AngleRestraint_exposer.def( - "usesMoleculesIn" - , usesMoleculesIn_function_value - , ( bp::arg("forcetable") ) - , bp::release_gil_policy() - , "Return whether or not this restraint involves any of the molecules\nthat are in the forcetable forcetable" ); - - } - { //::SireMM::AngleRestraint::usesMoleculesIn - - typedef bool ( ::SireMM::AngleRestraint::*usesMoleculesIn_function_type)( ::SireMol::Molecules const & ) const; - usesMoleculesIn_function_type usesMoleculesIn_function_value( &::SireMM::AngleRestraint::usesMoleculesIn ); - - AngleRestraint_exposer.def( - "usesMoleculesIn" - , usesMoleculesIn_function_value - , ( bp::arg("molecules") ) - , bp::release_gil_policy() - , "Return whether or not this restraint involves any of the molecules\nin molecules" ); + , "" ); } - AngleRestraint_exposer.staticmethod( "halfHarmonic" ); - AngleRestraint_exposer.staticmethod( "harmonic" ); - AngleRestraint_exposer.staticmethod( "theta" ); AngleRestraint_exposer.staticmethod( "typeName" ); AngleRestraint_exposer.def( "__copy__", &__copy__); AngleRestraint_exposer.def( "__deepcopy__", &__copy__); AngleRestraint_exposer.def( "clone", &__copy__); AngleRestraint_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::AngleRestraint >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); AngleRestraint_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::AngleRestraint >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); AngleRestraint_exposer.def_pickle(sire_pickle_suite< ::SireMM::AngleRestraint >()); AngleRestraint_exposer.def( "__str__", &__str__< ::SireMM::AngleRestraint > ); AngleRestraint_exposer.def( "__repr__", &__str__< ::SireMM::AngleRestraint > ); diff --git a/wrapper/MM/AngleRestraints.pypp.cpp b/wrapper/MM/AngleRestraints.pypp.cpp new file mode 100644 index 000000000..07edf38ca --- /dev/null +++ b/wrapper/MM/AngleRestraints.pypp.cpp @@ -0,0 +1,242 @@ +// This file has been generated by Py++. + +// (C) Christopher Woods, GPL >= 3 License + +#include "boost/python.hpp" +#include "Helpers/clone_const_reference.hpp" +#include "AngleRestraints.pypp.hpp" + +namespace bp = boost::python; + +#include "SireCAS/errors.h" + +#include "SireID/index.h" + +#include "SireStream/datastream.h" + +#include "SireStream/shareddatastream.h" + +#include "SireUnits/units.h" + +#include "anglerestraints.h" + +#include + +#include "anglerestraints.h" + +SireMM::AngleRestraints __copy__(const SireMM::AngleRestraints &other){ return SireMM::AngleRestraints(other); } + +#include "Helpers/copy.hpp" + +#include "Qt/qdatastream.hpp" + +#include "Helpers/str.hpp" + +#include "Helpers/release_gil_policy.hpp" + +#include "Helpers/len.hpp" + +void register_AngleRestraints_class(){ + + { //::SireMM::AngleRestraints + typedef bp::class_< SireMM::AngleRestraints, bp::bases< SireMM::Restraints, SireBase::Property > > AngleRestraints_exposer_t; + AngleRestraints_exposer_t AngleRestraints_exposer = AngleRestraints_exposer_t( "AngleRestraints", "This class represents a collection of angle restraints", bp::init< >("Null constructor") ); + bp::scope AngleRestraints_scope( AngleRestraints_exposer ); + AngleRestraints_exposer.def( bp::init< QString const & >(( bp::arg("name") ), "") ); + AngleRestraints_exposer.def( bp::init< SireMM::AngleRestraint const & >(( bp::arg("restraint") ), "") ); + AngleRestraints_exposer.def( bp::init< QList< SireMM::AngleRestraint > const & >(( bp::arg("restraints") ), "") ); + AngleRestraints_exposer.def( bp::init< QString const &, SireMM::AngleRestraint const & >(( bp::arg("name"), bp::arg("restraint") ), "") ); + AngleRestraints_exposer.def( bp::init< QString const &, QList< SireMM::AngleRestraint > const & >(( bp::arg("name"), bp::arg("restraints") ), "") ); + AngleRestraints_exposer.def( bp::init< SireMM::AngleRestraints const & >(( bp::arg("other") ), "") ); + { //::SireMM::AngleRestraints::add + + typedef void ( ::SireMM::AngleRestraints::*add_function_type)( ::SireMM::AngleRestraint const & ) ; + add_function_type add_function_value( &::SireMM::AngleRestraints::add ); + + AngleRestraints_exposer.def( + "add" + , add_function_value + , ( bp::arg("restraint") ) + , bp::release_gil_policy() + , "Add a restraint onto the list" ); + + } + { //::SireMM::AngleRestraints::add + + typedef void ( ::SireMM::AngleRestraints::*add_function_type)( ::SireMM::AngleRestraints const & ) ; + add_function_type add_function_value( &::SireMM::AngleRestraints::add ); + + AngleRestraints_exposer.def( + "add" + , add_function_value + , ( bp::arg("restraints") ) + , bp::release_gil_policy() + , "Add a restraint onto the list" ); + + } + { //::SireMM::AngleRestraints::at + + typedef ::SireMM::AngleRestraint const & ( ::SireMM::AngleRestraints::*at_function_type)( int ) const; + at_function_type at_function_value( &::SireMM::AngleRestraints::at ); + + AngleRestraints_exposer.def( + "at" + , at_function_value + , ( bp::arg("i") ) + , bp::return_value_policy() + , "Return the ith restraint" ); + + } + { //::SireMM::AngleRestraints::count + + typedef int ( ::SireMM::AngleRestraints::*count_function_type)( ) const; + count_function_type count_function_value( &::SireMM::AngleRestraints::count ); + + AngleRestraints_exposer.def( + "count" + , count_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + { //::SireMM::AngleRestraints::isEmpty + + typedef bool ( ::SireMM::AngleRestraints::*isEmpty_function_type)( ) const; + isEmpty_function_type isEmpty_function_value( &::SireMM::AngleRestraints::isEmpty ); + + AngleRestraints_exposer.def( + "isEmpty" + , isEmpty_function_value + , bp::release_gil_policy() + , "Return whether or not this is empty" ); + + } + { //::SireMM::AngleRestraints::isNull + + typedef bool ( ::SireMM::AngleRestraints::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::AngleRestraints::isNull ); + + AngleRestraints_exposer.def( + "isNull" + , isNull_function_value + , bp::release_gil_policy() + , "Return whether or not this is empty" ); + + } + { //::SireMM::AngleRestraints::nRestraints + + typedef int ( ::SireMM::AngleRestraints::*nRestraints_function_type)( ) const; + nRestraints_function_type nRestraints_function_value( &::SireMM::AngleRestraints::nRestraints ); + + AngleRestraints_exposer.def( + "nRestraints" + , nRestraints_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + AngleRestraints_exposer.def( bp::self != bp::self ); + AngleRestraints_exposer.def( bp::self + bp::other< SireMM::AngleRestraint >() ); + AngleRestraints_exposer.def( bp::self + bp::self ); + { //::SireMM::AngleRestraints::operator= + + typedef ::SireMM::AngleRestraints & ( ::SireMM::AngleRestraints::*assign_function_type)( ::SireMM::AngleRestraints const & ) ; + assign_function_type assign_function_value( &::SireMM::AngleRestraints::operator= ); + + AngleRestraints_exposer.def( + "assign" + , assign_function_value + , ( bp::arg("other") ) + , bp::return_self< >() + , "" ); + + } + AngleRestraints_exposer.def( bp::self == bp::self ); + { //::SireMM::AngleRestraints::operator[] + + typedef ::SireMM::AngleRestraint const & ( ::SireMM::AngleRestraints::*__getitem___function_type)( int ) const; + __getitem___function_type __getitem___function_value( &::SireMM::AngleRestraints::operator[] ); + + AngleRestraints_exposer.def( + "__getitem__" + , __getitem___function_value + , ( bp::arg("i") ) + , bp::return_value_policy() + , "" ); + + } + { //::SireMM::AngleRestraints::restraints + + typedef ::QList< SireMM::AngleRestraint > ( ::SireMM::AngleRestraints::*restraints_function_type)( ) const; + restraints_function_type restraints_function_value( &::SireMM::AngleRestraints::restraints ); + + AngleRestraints_exposer.def( + "restraints" + , restraints_function_value + , bp::release_gil_policy() + , "Return all of the restraints" ); + + } + { //::SireMM::AngleRestraints::size + + typedef int ( ::SireMM::AngleRestraints::*size_function_type)( ) const; + size_function_type size_function_value( &::SireMM::AngleRestraints::size ); + + AngleRestraints_exposer.def( + "size" + , size_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + { //::SireMM::AngleRestraints::toString + + typedef ::QString ( ::SireMM::AngleRestraints::*toString_function_type)( ) const; + toString_function_type toString_function_value( &::SireMM::AngleRestraints::toString ); + + AngleRestraints_exposer.def( + "toString" + , toString_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMM::AngleRestraints::typeName + + typedef char const * ( *typeName_function_type )( ); + typeName_function_type typeName_function_value( &::SireMM::AngleRestraints::typeName ); + + AngleRestraints_exposer.def( + "typeName" + , typeName_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMM::AngleRestraints::what + + typedef char const * ( ::SireMM::AngleRestraints::*what_function_type)( ) const; + what_function_type what_function_value( &::SireMM::AngleRestraints::what ); + + AngleRestraints_exposer.def( + "what" + , what_function_value + , bp::release_gil_policy() + , "" ); + + } + AngleRestraints_exposer.staticmethod( "typeName" ); + AngleRestraints_exposer.def( "__copy__", &__copy__); + AngleRestraints_exposer.def( "__deepcopy__", &__copy__); + AngleRestraints_exposer.def( "clone", &__copy__); + AngleRestraints_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::AngleRestraints >, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + AngleRestraints_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::AngleRestraints >, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + AngleRestraints_exposer.def_pickle(sire_pickle_suite< ::SireMM::AngleRestraints >()); + AngleRestraints_exposer.def( "__str__", &__str__< ::SireMM::AngleRestraints > ); + AngleRestraints_exposer.def( "__repr__", &__str__< ::SireMM::AngleRestraints > ); + AngleRestraints_exposer.def( "__len__", &__len_size< ::SireMM::AngleRestraints > ); + } + +} diff --git a/wrapper/MM/AngleRestraints.pypp.hpp b/wrapper/MM/AngleRestraints.pypp.hpp new file mode 100644 index 000000000..b1292a5cd --- /dev/null +++ b/wrapper/MM/AngleRestraints.pypp.hpp @@ -0,0 +1,10 @@ +// This file has been generated by Py++. + +// (C) Christopher Woods, GPL >= 3 License + +#ifndef AngleRestraints_hpp__pyplusplus_wrapper +#define AngleRestraints_hpp__pyplusplus_wrapper + +void register_AngleRestraints_class(); + +#endif//AngleRestraints_hpp__pyplusplus_wrapper diff --git a/wrapper/MM/CMakeAutogenFile.txt b/wrapper/MM/CMakeAutogenFile.txt index 3eeb5da11..3c3a4985a 100644 --- a/wrapper/MM/CMakeAutogenFile.txt +++ b/wrapper/MM/CMakeAutogenFile.txt @@ -68,6 +68,7 @@ set ( PYPP_SOURCES PositionalRestraints.pypp.cpp IntraGroupFF.pypp.cpp DihedralRestraint.pypp.cpp + DihedralRestraints.pypp.cpp IntraFF.pypp.cpp DihedralComponent.pypp.cpp ScaledChargeParameterNames3D.pypp.cpp @@ -181,6 +182,7 @@ set ( PYPP_SOURCES CoulombNBPairs.pypp.cpp AmberDihedral.pypp.cpp AngleRestraint.pypp.cpp + AngleRestraints.pypp.cpp TwoAtomFunctions.pypp.cpp IntraSoftCLJFFBase.pypp.cpp BondSymbols.pypp.cpp diff --git a/wrapper/MM/DihedralRestraint.pypp.cpp b/wrapper/MM/DihedralRestraint.pypp.cpp index 65ed2feb1..a5d04277e 100644 --- a/wrapper/MM/DihedralRestraint.pypp.cpp +++ b/wrapper/MM/DihedralRestraint.pypp.cpp @@ -3,36 +3,25 @@ // (C) Christopher Woods, GPL >= 3 License #include "boost/python.hpp" -#include "Helpers/clone_const_reference.hpp" #include "DihedralRestraint.pypp.hpp" namespace bp = boost::python; -#include "SireCAS/conditional.h" - #include "SireCAS/errors.h" -#include "SireCAS/power.h" - -#include "SireCAS/symbols.h" - -#include "SireCAS/values.h" - -#include "SireFF/forcetable.h" - #include "SireID/index.h" #include "SireStream/datastream.h" #include "SireStream/shareddatastream.h" -#include "SireUnits/angle.h" - #include "SireUnits/units.h" -#include "dihedralrestraint.h" +#include "dihedralrestraints.h" + +#include -#include "dihedralrestraint.h" +#include "dihedralrestraints.h" SireMM::DihedralRestraint __copy__(const SireMM::DihedralRestraint &other){ return SireMM::DihedralRestraint(other); } @@ -47,162 +36,50 @@ SireMM::DihedralRestraint __copy__(const SireMM::DihedralRestraint &other){ retu void register_DihedralRestraint_class(){ { //::SireMM::DihedralRestraint - typedef bp::class_< SireMM::DihedralRestraint, bp::bases< SireMM::Restraint3D, SireMM::Restraint, SireBase::Property > > DihedralRestraint_exposer_t; - DihedralRestraint_exposer_t DihedralRestraint_exposer = DihedralRestraint_exposer_t( "DihedralRestraint", "This is a restraint that operates on the dihedral angle between\nfour SireMM::Point objects (e.g. four atoms in a molecule)\n\nAuthor: Christopher Woods\n", bp::init< >("Constructor") ); + typedef bp::class_< SireMM::DihedralRestraint, bp::bases< SireBase::Property > > DihedralRestraint_exposer_t; + DihedralRestraint_exposer_t DihedralRestraint_exposer = DihedralRestraint_exposer_t( "DihedralRestraint", "This class represents a single angle restraint between any three\natoms in a system\nAuthor: Christopher Woods\n", bp::init< >("Null constructor") ); bp::scope DihedralRestraint_scope( DihedralRestraint_exposer ); - DihedralRestraint_exposer.def( bp::init< SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireCAS::Expression const & >(( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("point3"), bp::arg("restraint") ), "Construct a restraint that acts on the angle within the\nthree points point0, point1 and point2 (theta == a(012)),\nrestraining the angle within these points using the expression\nrestraint") ); - DihedralRestraint_exposer.def( bp::init< SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireFF::PointRef const &, SireCAS::Expression const &, SireCAS::Values const & >(( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("point3"), bp::arg("restraint"), bp::arg("values") ), "Construct a restraint that acts on the angle within the\nthree points point0, point1 and point2 (theta == a(012)),\nrestraining the angle within these points using the expression\nrestraint") ); + DihedralRestraint_exposer.def( bp::init< QList< long long > const &, SireUnits::Dimension::Angle const &, SireUnits::Dimension::HarmonicAngleConstant const & >(( bp::arg("atoms"), bp::arg("phi0"), bp::arg("kphi") ), "Construct a restraint that acts on the angle within the\nfour atoms atom0, atom1, atom2 atom3 (phi == a(0123)),\nrestraining the angle within these atoms") ); DihedralRestraint_exposer.def( bp::init< SireMM::DihedralRestraint const & >(( bp::arg("other") ), "Copy constructor") ); - { //::SireMM::DihedralRestraint::builtinSymbols + { //::SireMM::DihedralRestraint::atoms - typedef ::SireCAS::Symbols ( ::SireMM::DihedralRestraint::*builtinSymbols_function_type)( ) const; - builtinSymbols_function_type builtinSymbols_function_value( &::SireMM::DihedralRestraint::builtinSymbols ); + typedef ::QVector< long long > ( ::SireMM::DihedralRestraint::*atoms_function_type)( ) const; + atoms_function_type atoms_function_value( &::SireMM::DihedralRestraint::atoms ); DihedralRestraint_exposer.def( - "builtinSymbols" - , builtinSymbols_function_value + "atoms" + , atoms_function_value , bp::release_gil_policy() - , "Return the built-in symbols for this restraint" ); + , "Return the atoms involved in the restraint" ); } - { //::SireMM::DihedralRestraint::builtinValues + { //::SireMM::DihedralRestraint::isNull - typedef ::SireCAS::Values ( ::SireMM::DihedralRestraint::*builtinValues_function_type)( ) const; - builtinValues_function_type builtinValues_function_value( &::SireMM::DihedralRestraint::builtinValues ); + typedef bool ( ::SireMM::DihedralRestraint::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::DihedralRestraint::isNull ); DihedralRestraint_exposer.def( - "builtinValues" - , builtinValues_function_value + "isNull" + , isNull_function_value , bp::release_gil_policy() - , "Return the values of the built-in symbols of this restraint" ); - - } - { //::SireMM::DihedralRestraint::contains - - typedef bool ( ::SireMM::DihedralRestraint::*contains_function_type)( ::SireMol::MolNum ) const; - contains_function_type contains_function_value( &::SireMM::DihedralRestraint::contains ); - - DihedralRestraint_exposer.def( - "contains" - , contains_function_value - , ( bp::arg("molnum") ) - , bp::release_gil_policy() - , "Return whether or not this restraint affects the molecule\nwith number molnum" ); - - } - { //::SireMM::DihedralRestraint::contains - - typedef bool ( ::SireMM::DihedralRestraint::*contains_function_type)( ::SireMol::MolID const & ) const; - contains_function_type contains_function_value( &::SireMM::DihedralRestraint::contains ); - - DihedralRestraint_exposer.def( - "contains" - , contains_function_value - , ( bp::arg("molid") ) - , bp::release_gil_policy() - , "Return whether or not this restraint affects the molecule\nwith ID molid" ); - - } - { //::SireMM::DihedralRestraint::differentialRestraintFunction - - typedef ::SireCAS::Expression const & ( ::SireMM::DihedralRestraint::*differentialRestraintFunction_function_type)( ) const; - differentialRestraintFunction_function_type differentialRestraintFunction_function_value( &::SireMM::DihedralRestraint::differentialRestraintFunction ); - - DihedralRestraint_exposer.def( - "differentialRestraintFunction" - , differentialRestraintFunction_function_value - , bp::return_value_policy< bp::copy_const_reference >() - , "Return the function used to calculate the restraint force" ); - - } - { //::SireMM::DihedralRestraint::differentiate - - typedef ::SireMM::RestraintPtr ( ::SireMM::DihedralRestraint::*differentiate_function_type)( ::SireCAS::Symbol const & ) const; - differentiate_function_type differentiate_function_value( &::SireMM::DihedralRestraint::differentiate ); - - DihedralRestraint_exposer.def( - "differentiate" - , differentiate_function_value - , ( bp::arg("symbol") ) - , bp::release_gil_policy() - , "Return the differential of this restraint with respect to\nthe symbol symbol\nThrow: SireCAS::unavailable_differential\n" ); - - } - { //::SireMM::DihedralRestraint::force - - typedef void ( ::SireMM::DihedralRestraint::*force_function_type)( ::SireFF::MolForceTable &,double ) const; - force_function_type force_function_value( &::SireMM::DihedralRestraint::force ); - - DihedralRestraint_exposer.def( - "force" - , force_function_value - , ( bp::arg("forcetable"), bp::arg("scale_force")=1 ) - , "Calculate the force acting on the molecule in the forcetable forcetable\ncaused by this restraint, and add it on to the forcetable scaled by\nscale_force" ); - - } - { //::SireMM::DihedralRestraint::force - - typedef void ( ::SireMM::DihedralRestraint::*force_function_type)( ::SireFF::ForceTable &,double ) const; - force_function_type force_function_value( &::SireMM::DihedralRestraint::force ); - - DihedralRestraint_exposer.def( - "force" - , force_function_value - , ( bp::arg("forcetable"), bp::arg("scale_force")=1 ) - , "Calculate the force acting on the molecules in the forcetable forcetable\ncaused by this restraint, and add it on to the forcetable scaled by\nscale_force" ); - - } - { //::SireMM::DihedralRestraint::halfHarmonic - - typedef ::SireMM::DihedralRestraint ( *halfHarmonic_function_type )( ::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireUnits::Dimension::Angle const &,::SireMM::HarmonicAngleForceConstant const & ); - halfHarmonic_function_type halfHarmonic_function_value( &::SireMM::DihedralRestraint::halfHarmonic ); - - DihedralRestraint_exposer.def( - "halfHarmonic" - , halfHarmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("point3"), bp::arg("angle"), bp::arg("force_constant") ) - , bp::release_gil_policy() - , "Return a distance restraint that applied a half-harmonic potential\nbetween the points point0 and point1 above a distance distance\nusing a force constant force_constant" ); - - } - { //::SireMM::DihedralRestraint::harmonic - - typedef ::SireMM::DihedralRestraint ( *harmonic_function_type )( ::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireFF::PointRef const &,::SireMM::HarmonicAngleForceConstant const & ); - harmonic_function_type harmonic_function_value( &::SireMM::DihedralRestraint::harmonic ); - - DihedralRestraint_exposer.def( - "harmonic" - , harmonic_function_value - , ( bp::arg("point0"), bp::arg("point1"), bp::arg("point2"), bp::arg("point3"), bp::arg("force_constant") ) - , bp::release_gil_policy() - , "Return a distance restraint that applies a harmonic potential between\nthe points point0 and point1 using a force constant force_constant" ); - - } - { //::SireMM::DihedralRestraint::molecules - - typedef ::SireMol::Molecules ( ::SireMM::DihedralRestraint::*molecules_function_type)( ) const; - molecules_function_type molecules_function_value( &::SireMM::DihedralRestraint::molecules ); - - DihedralRestraint_exposer.def( - "molecules" - , molecules_function_value - , bp::release_gil_policy() - , "Return the molecules used in this restraint" ); + , "" ); } - { //::SireMM::DihedralRestraint::nPoints + { //::SireMM::DihedralRestraint::kphi - typedef int ( ::SireMM::DihedralRestraint::*nPoints_function_type)( ) const; - nPoints_function_type nPoints_function_value( &::SireMM::DihedralRestraint::nPoints ); + typedef ::SireUnits::Dimension::HarmonicAngleConstant ( ::SireMM::DihedralRestraint::*kphi_function_type)( ) const; + kphi_function_type kphi_function_value( &::SireMM::DihedralRestraint::kphi ); DihedralRestraint_exposer.def( - "nPoints" - , nPoints_function_value + "kphi" + , kphi_function_value , bp::release_gil_policy() - , "This restraint involves four points" ); + , "Return the force constant for the restraint" ); } DihedralRestraint_exposer.def( bp::self != bp::self ); + DihedralRestraint_exposer.def( bp::self + bp::self ); + DihedralRestraint_exposer.def( bp::self + bp::other< SireMM::DihedralRestraints >() ); { //::SireMM::DihedralRestraint::operator= typedef ::SireMM::DihedralRestraint & ( ::SireMM::DihedralRestraint::*assign_function_type)( ::SireMM::DihedralRestraint const & ) ; @@ -217,90 +94,28 @@ void register_DihedralRestraint_class(){ } DihedralRestraint_exposer.def( bp::self == bp::self ); - { //::SireMM::DihedralRestraint::phi - - typedef ::SireCAS::Symbol const & ( *phi_function_type )( ); - phi_function_type phi_function_value( &::SireMM::DihedralRestraint::phi ); - - DihedralRestraint_exposer.def( - "phi" - , phi_function_value - , bp::return_value_policy() - , "Return the symbol that represents the dihedral angle between the points (phi)" ); - - } - { //::SireMM::DihedralRestraint::point - - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point_function_type)( int ) const; - point_function_type point_function_value( &::SireMM::DihedralRestraint::point ); - - DihedralRestraint_exposer.def( - "point" - , point_function_value - , ( bp::arg("i") ) - , bp::return_value_policy() - , "Return the ith point" ); - - } - { //::SireMM::DihedralRestraint::point0 - - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point0_function_type)( ) const; - point0_function_type point0_function_value( &::SireMM::DihedralRestraint::point0 ); - - DihedralRestraint_exposer.def( - "point0" - , point0_function_value - , bp::return_value_policy() - , "Return the first point" ); - - } - { //::SireMM::DihedralRestraint::point1 - - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point1_function_type)( ) const; - point1_function_type point1_function_value( &::SireMM::DihedralRestraint::point1 ); - - DihedralRestraint_exposer.def( - "point1" - , point1_function_value - , bp::return_value_policy() - , "Return the second point" ); - - } - { //::SireMM::DihedralRestraint::point2 + { //::SireMM::DihedralRestraint::phi0 - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point2_function_type)( ) const; - point2_function_type point2_function_value( &::SireMM::DihedralRestraint::point2 ); + typedef ::SireUnits::Dimension::Angle ( ::SireMM::DihedralRestraint::*phi0_function_type)( ) const; + phi0_function_type phi0_function_value( &::SireMM::DihedralRestraint::phi0 ); DihedralRestraint_exposer.def( - "point2" - , point2_function_value - , bp::return_value_policy() - , "Return the third point" ); - - } - { //::SireMM::DihedralRestraint::point3 - - typedef ::SireFF::Point const & ( ::SireMM::DihedralRestraint::*point3_function_type)( ) const; - point3_function_type point3_function_value( &::SireMM::DihedralRestraint::point3 ); - - DihedralRestraint_exposer.def( - "point3" - , point3_function_value - , bp::return_value_policy() - , "Return the fourth point" ); + "phi0" + , phi0_function_value + , bp::release_gil_policy() + , "Return the equilibrium angle for the restraint" ); } - { //::SireMM::DihedralRestraint::setSpace + { //::SireMM::DihedralRestraint::toString - typedef void ( ::SireMM::DihedralRestraint::*setSpace_function_type)( ::SireVol::Space const & ) ; - setSpace_function_type setSpace_function_value( &::SireMM::DihedralRestraint::setSpace ); + typedef ::QString ( ::SireMM::DihedralRestraint::*toString_function_type)( ) const; + toString_function_type toString_function_value( &::SireMM::DihedralRestraint::toString ); DihedralRestraint_exposer.def( - "setSpace" - , setSpace_function_value - , ( bp::arg("space") ) + "toString" + , toString_function_value , bp::release_gil_policy() - , "Set the space used to evaluate the energy of this restraint\nThrow: SireVol::incompatible_space\n" ); + , "" ); } { //::SireMM::DihedralRestraint::typeName @@ -315,69 +130,26 @@ void register_DihedralRestraint_class(){ , "" ); } - { //::SireMM::DihedralRestraint::update - - typedef void ( ::SireMM::DihedralRestraint::*update_function_type)( ::SireMol::MoleculeData const & ) ; - update_function_type update_function_value( &::SireMM::DihedralRestraint::update ); - - DihedralRestraint_exposer.def( - "update" - , update_function_value - , ( bp::arg("moldata") ) - , bp::release_gil_policy() - , "Update the points of this restraint using new molecule data from moldata\nThrow: SireBase::missing_property\nThrow: SireError::invalid_cast\nThrow: SireError::incompatible_error\n" ); - - } - { //::SireMM::DihedralRestraint::update + { //::SireMM::DihedralRestraint::what - typedef void ( ::SireMM::DihedralRestraint::*update_function_type)( ::SireMol::Molecules const & ) ; - update_function_type update_function_value( &::SireMM::DihedralRestraint::update ); + typedef char const * ( ::SireMM::DihedralRestraint::*what_function_type)( ) const; + what_function_type what_function_value( &::SireMM::DihedralRestraint::what ); DihedralRestraint_exposer.def( - "update" - , update_function_value - , ( bp::arg("molecules") ) + "what" + , what_function_value , bp::release_gil_policy() - , "Update the points of this restraint using new molecule data from molecules\nThrow: SireBase::missing_property\nThrow: SireError::invalid_cast\nThrow: SireError::incompatible_error\n" ); - - } - { //::SireMM::DihedralRestraint::usesMoleculesIn - - typedef bool ( ::SireMM::DihedralRestraint::*usesMoleculesIn_function_type)( ::SireFF::ForceTable const & ) const; - usesMoleculesIn_function_type usesMoleculesIn_function_value( &::SireMM::DihedralRestraint::usesMoleculesIn ); - - DihedralRestraint_exposer.def( - "usesMoleculesIn" - , usesMoleculesIn_function_value - , ( bp::arg("forcetable") ) - , bp::release_gil_policy() - , "Return whether or not this restraint involves any of the molecules\nthat are in the forcetable forcetable" ); - - } - { //::SireMM::DihedralRestraint::usesMoleculesIn - - typedef bool ( ::SireMM::DihedralRestraint::*usesMoleculesIn_function_type)( ::SireMol::Molecules const & ) const; - usesMoleculesIn_function_type usesMoleculesIn_function_value( &::SireMM::DihedralRestraint::usesMoleculesIn ); - - DihedralRestraint_exposer.def( - "usesMoleculesIn" - , usesMoleculesIn_function_value - , ( bp::arg("molecules") ) - , bp::release_gil_policy() - , "Return whether or not this restraint involves any of the molecules\nin molecules" ); + , "" ); } - DihedralRestraint_exposer.staticmethod( "halfHarmonic" ); - DihedralRestraint_exposer.staticmethod( "harmonic" ); - DihedralRestraint_exposer.staticmethod( "phi" ); DihedralRestraint_exposer.staticmethod( "typeName" ); DihedralRestraint_exposer.def( "__copy__", &__copy__); DihedralRestraint_exposer.def( "__deepcopy__", &__copy__); DihedralRestraint_exposer.def( "clone", &__copy__); DihedralRestraint_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::DihedralRestraint >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); DihedralRestraint_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::DihedralRestraint >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); DihedralRestraint_exposer.def_pickle(sire_pickle_suite< ::SireMM::DihedralRestraint >()); DihedralRestraint_exposer.def( "__str__", &__str__< ::SireMM::DihedralRestraint > ); DihedralRestraint_exposer.def( "__repr__", &__str__< ::SireMM::DihedralRestraint > ); diff --git a/wrapper/MM/DihedralRestraints.pypp.cpp b/wrapper/MM/DihedralRestraints.pypp.cpp new file mode 100644 index 000000000..44d49c86a --- /dev/null +++ b/wrapper/MM/DihedralRestraints.pypp.cpp @@ -0,0 +1,242 @@ +// This file has been generated by Py++. + +// (C) Christopher Woods, GPL >= 3 License + +#include "boost/python.hpp" +#include "Helpers/clone_const_reference.hpp" +#include "DihedralRestraints.pypp.hpp" + +namespace bp = boost::python; + +#include "SireCAS/errors.h" + +#include "SireID/index.h" + +#include "SireStream/datastream.h" + +#include "SireStream/shareddatastream.h" + +#include "SireUnits/units.h" + +#include "dihedralrestraints.h" + +#include + +#include "dihedralrestraints.h" + +SireMM::DihedralRestraints __copy__(const SireMM::DihedralRestraints &other){ return SireMM::DihedralRestraints(other); } + +#include "Helpers/copy.hpp" + +#include "Qt/qdatastream.hpp" + +#include "Helpers/str.hpp" + +#include "Helpers/release_gil_policy.hpp" + +#include "Helpers/len.hpp" + +void register_DihedralRestraints_class(){ + + { //::SireMM::DihedralRestraints + typedef bp::class_< SireMM::DihedralRestraints, bp::bases< SireMM::Restraints, SireBase::Property > > DihedralRestraints_exposer_t; + DihedralRestraints_exposer_t DihedralRestraints_exposer = DihedralRestraints_exposer_t( "DihedralRestraints", "This class represents a collection of angle restraints", bp::init< >("Null constructor") ); + bp::scope DihedralRestraints_scope( DihedralRestraints_exposer ); + DihedralRestraints_exposer.def( bp::init< QString const & >(( bp::arg("name") ), "") ); + DihedralRestraints_exposer.def( bp::init< SireMM::DihedralRestraint const & >(( bp::arg("restraint") ), "") ); + DihedralRestraints_exposer.def( bp::init< QList< SireMM::DihedralRestraint > const & >(( bp::arg("restraints") ), "") ); + DihedralRestraints_exposer.def( bp::init< QString const &, SireMM::DihedralRestraint const & >(( bp::arg("name"), bp::arg("restraint") ), "") ); + DihedralRestraints_exposer.def( bp::init< QString const &, QList< SireMM::DihedralRestraint > const & >(( bp::arg("name"), bp::arg("restraints") ), "") ); + DihedralRestraints_exposer.def( bp::init< SireMM::DihedralRestraints const & >(( bp::arg("other") ), "") ); + { //::SireMM::DihedralRestraints::add + + typedef void ( ::SireMM::DihedralRestraints::*add_function_type)( ::SireMM::DihedralRestraint const & ) ; + add_function_type add_function_value( &::SireMM::DihedralRestraints::add ); + + DihedralRestraints_exposer.def( + "add" + , add_function_value + , ( bp::arg("restraint") ) + , bp::release_gil_policy() + , "Add a restraints onto the list" ); + + } + { //::SireMM::DihedralRestraints::add + + typedef void ( ::SireMM::DihedralRestraints::*add_function_type)( ::SireMM::DihedralRestraints const & ) ; + add_function_type add_function_value( &::SireMM::DihedralRestraints::add ); + + DihedralRestraints_exposer.def( + "add" + , add_function_value + , ( bp::arg("restraints") ) + , bp::release_gil_policy() + , "Add a restraint onto the list" ); + + } + { //::SireMM::DihedralRestraints::at + + typedef ::SireMM::DihedralRestraint const & ( ::SireMM::DihedralRestraints::*at_function_type)( int ) const; + at_function_type at_function_value( &::SireMM::DihedralRestraints::at ); + + DihedralRestraints_exposer.def( + "at" + , at_function_value + , ( bp::arg("i") ) + , bp::return_value_policy() + , "Return the ith restraint" ); + + } + { //::SireMM::DihedralRestraints::count + + typedef int ( ::SireMM::DihedralRestraints::*count_function_type)( ) const; + count_function_type count_function_value( &::SireMM::DihedralRestraints::count ); + + DihedralRestraints_exposer.def( + "count" + , count_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + { //::SireMM::DihedralRestraints::isEmpty + + typedef bool ( ::SireMM::DihedralRestraints::*isEmpty_function_type)( ) const; + isEmpty_function_type isEmpty_function_value( &::SireMM::DihedralRestraints::isEmpty ); + + DihedralRestraints_exposer.def( + "isEmpty" + , isEmpty_function_value + , bp::release_gil_policy() + , "Return whether or not this is empty" ); + + } + { //::SireMM::DihedralRestraints::isNull + + typedef bool ( ::SireMM::DihedralRestraints::*isNull_function_type)( ) const; + isNull_function_type isNull_function_value( &::SireMM::DihedralRestraints::isNull ); + + DihedralRestraints_exposer.def( + "isNull" + , isNull_function_value + , bp::release_gil_policy() + , "Return whether or not this is empty" ); + + } + { //::SireMM::DihedralRestraints::nRestraints + + typedef int ( ::SireMM::DihedralRestraints::*nRestraints_function_type)( ) const; + nRestraints_function_type nRestraints_function_value( &::SireMM::DihedralRestraints::nRestraints ); + + DihedralRestraints_exposer.def( + "nRestraints" + , nRestraints_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + DihedralRestraints_exposer.def( bp::self != bp::self ); + DihedralRestraints_exposer.def( bp::self + bp::other< SireMM::DihedralRestraint >() ); + DihedralRestraints_exposer.def( bp::self + bp::self ); + { //::SireMM::DihedralRestraints::operator= + + typedef ::SireMM::DihedralRestraints & ( ::SireMM::DihedralRestraints::*assign_function_type)( ::SireMM::DihedralRestraints const & ) ; + assign_function_type assign_function_value( &::SireMM::DihedralRestraints::operator= ); + + DihedralRestraints_exposer.def( + "assign" + , assign_function_value + , ( bp::arg("other") ) + , bp::return_self< >() + , "" ); + + } + DihedralRestraints_exposer.def( bp::self == bp::self ); + { //::SireMM::DihedralRestraints::operator[] + + typedef ::SireMM::DihedralRestraint const & ( ::SireMM::DihedralRestraints::*__getitem___function_type)( int ) const; + __getitem___function_type __getitem___function_value( &::SireMM::DihedralRestraints::operator[] ); + + DihedralRestraints_exposer.def( + "__getitem__" + , __getitem___function_value + , ( bp::arg("i") ) + , bp::return_value_policy() + , "" ); + + } + { //::SireMM::DihedralRestraints::restraints + + typedef ::QList< SireMM::DihedralRestraint > ( ::SireMM::DihedralRestraints::*restraints_function_type)( ) const; + restraints_function_type restraints_function_value( &::SireMM::DihedralRestraints::restraints ); + + DihedralRestraints_exposer.def( + "restraints" + , restraints_function_value + , bp::release_gil_policy() + , "Return all of the restraints" ); + + } + { //::SireMM::DihedralRestraints::size + + typedef int ( ::SireMM::DihedralRestraints::*size_function_type)( ) const; + size_function_type size_function_value( &::SireMM::DihedralRestraints::size ); + + DihedralRestraints_exposer.def( + "size" + , size_function_value + , bp::release_gil_policy() + , "Return the number of restraints" ); + + } + { //::SireMM::DihedralRestraints::toString + + typedef ::QString ( ::SireMM::DihedralRestraints::*toString_function_type)( ) const; + toString_function_type toString_function_value( &::SireMM::DihedralRestraints::toString ); + + DihedralRestraints_exposer.def( + "toString" + , toString_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMM::DihedralRestraints::typeName + + typedef char const * ( *typeName_function_type )( ); + typeName_function_type typeName_function_value( &::SireMM::DihedralRestraints::typeName ); + + DihedralRestraints_exposer.def( + "typeName" + , typeName_function_value + , bp::release_gil_policy() + , "" ); + + } + { //::SireMM::DihedralRestraints::what + + typedef char const * ( ::SireMM::DihedralRestraints::*what_function_type)( ) const; + what_function_type what_function_value( &::SireMM::DihedralRestraints::what ); + + DihedralRestraints_exposer.def( + "what" + , what_function_value + , bp::release_gil_policy() + , "" ); + + } + DihedralRestraints_exposer.staticmethod( "typeName" ); + DihedralRestraints_exposer.def( "__copy__", &__copy__); + DihedralRestraints_exposer.def( "__deepcopy__", &__copy__); + DihedralRestraints_exposer.def( "clone", &__copy__); + DihedralRestraints_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::DihedralRestraints >, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + DihedralRestraints_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::DihedralRestraints >, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); + DihedralRestraints_exposer.def_pickle(sire_pickle_suite< ::SireMM::DihedralRestraints >()); + DihedralRestraints_exposer.def( "__str__", &__str__< ::SireMM::DihedralRestraints > ); + DihedralRestraints_exposer.def( "__repr__", &__str__< ::SireMM::DihedralRestraints > ); + DihedralRestraints_exposer.def( "__len__", &__len_size< ::SireMM::DihedralRestraints > ); + } + +} diff --git a/wrapper/MM/DihedralRestraints.pypp.hpp b/wrapper/MM/DihedralRestraints.pypp.hpp new file mode 100644 index 000000000..0aab6519a --- /dev/null +++ b/wrapper/MM/DihedralRestraints.pypp.hpp @@ -0,0 +1,10 @@ +// This file has been generated by Py++. + +// (C) Christopher Woods, GPL >= 3 License + +#ifndef DihedralRestraints_hpp__pyplusplus_wrapper +#define DihedralRestraints_hpp__pyplusplus_wrapper + +void register_DihedralRestraints_class(); + +#endif//DihedralRestraints_hpp__pyplusplus_wrapper diff --git a/wrapper/MM/SireMM_registrars.cpp b/wrapper/MM/SireMM_registrars.cpp index b7c70f668..9dabd1bd9 100644 --- a/wrapper/MM/SireMM_registrars.cpp +++ b/wrapper/MM/SireMM_registrars.cpp @@ -16,7 +16,7 @@ #include "selectordihedral.h" #include "clj14group.h" #include "improper.h" -#include "anglerestraint.h" +#include "anglerestraints.h" #include "bondrestraints.h" #include "twoatomfunctions.h" #include "boreschrestraints.h" @@ -54,7 +54,7 @@ #include "cljboxes.h" #include "internalgroupff.h" #include "internalcomponent.h" -#include "dihedralrestraint.h" +#include "dihedralrestraints.h" #include "selectorimproper.h" #include "intergroupff.h" #include "dihedral.h" @@ -112,6 +112,7 @@ void register_SireMM_objects() ObjectRegistry::registerConverterFor< SireMM::Improper >(); ObjectRegistry::registerConverterFor< SireMol::Mover >(); ObjectRegistry::registerConverterFor< SireMM::AngleRestraint >(); + ObjectRegistry::registerConverterFor< SireMM::AngleRestraints >(); ObjectRegistry::registerConverterFor< SireMM::BondRestraint >(); ObjectRegistry::registerConverterFor< SireMM::BondRestraints >(); ObjectRegistry::registerConverterFor< SireMM::TwoAtomFunctions >(); @@ -200,6 +201,7 @@ void register_SireMM_objects() ObjectRegistry::registerConverterFor< SireMM::Intra14Component >(); ObjectRegistry::registerConverterFor< SireMM::InternalComponent >(); ObjectRegistry::registerConverterFor< SireMM::DihedralRestraint >(); + ObjectRegistry::registerConverterFor< SireMM::DihedralRestraints >(); ObjectRegistry::registerConverterFor< SireMM::SelectorImproper >(); ObjectRegistry::registerConverterFor< SireMol::Mover >(); ObjectRegistry::registerConverterFor< SireMM::InterGroupFF >(); diff --git a/wrapper/MM/_MM.main.cpp b/wrapper/MM/_MM.main.cpp index 5a98d9755..addf5ef01 100644 --- a/wrapper/MM/_MM.main.cpp +++ b/wrapper/MM/_MM.main.cpp @@ -27,6 +27,8 @@ #include "AngleRestraint.pypp.hpp" +#include "AngleRestraints.pypp.hpp" + #include "AngleSymbols.pypp.hpp" #include "AtomFunction.pypp.hpp" @@ -161,6 +163,8 @@ #include "DihedralRestraint.pypp.hpp" +#include "DihedralRestraints.pypp.hpp" + #include "DihedralSymbols.pypp.hpp" #include "DistanceRestraint.pypp.hpp" @@ -552,6 +556,10 @@ BOOST_PYTHON_MODULE(_MM){ register_AngleRestraint_class(); + register_Restraints_class(); + + register_AngleRestraints_class(); + register_InternalSymbolsBase_class(); register_AngleSymbols_class(); @@ -586,8 +594,6 @@ BOOST_PYTHON_MODULE(_MM){ register_BondRestraint_class(); - register_Restraints_class(); - register_BondRestraints_class(); register_BondSymbols_class(); @@ -682,6 +688,8 @@ BOOST_PYTHON_MODULE(_MM){ register_DihedralRestraint_class(); + register_DihedralRestraints_class(); + register_DihedralSymbols_class(); register_DistanceRestraint_class(); diff --git a/wrapper/MM/active_headers.h b/wrapper/MM/active_headers.h index 3cb17a25a..c50bbe8c8 100644 --- a/wrapper/MM/active_headers.h +++ b/wrapper/MM/active_headers.h @@ -5,7 +5,7 @@ #include "amberparams.h" #include "angle.h" -#include "anglerestraint.h" +#include "anglerestraints.h" #include "atomfunctions.h" #include "atomljs.h" #include "bond.h" @@ -30,7 +30,7 @@ #include "cljworkspace.h" #include "coulombpotential.h" #include "dihedral.h" -#include "dihedralrestraint.h" +#include "dihedralrestraints.h" #include "distancerestraint.h" #include "excludedpairs.h" #include "fouratomfunctions.h"