RDKit
Open-source cheminformatics and machine learning.
ForceField.h
Go to the documentation of this file.
1 //
2 // Copyright (C) 2004-2006 Rational Discovery LLC
3 //
4 // @@ All Rights Reserved @@
5 // This file is part of the RDKit.
6 // The contents are covered by the terms of the BSD license
7 // which is included in the file license.txt, found at the root
8 // of the RDKit source tree.
9 //
10 #include <RDGeneral/export.h>
11 #ifndef __RD_FORCEFIELD_H__
12 #define __RD_FORCEFIELD_H__
13 
14 #include <vector>
15 #include <boost/smart_ptr.hpp>
16 #include <boost/foreach.hpp>
17 #include <Geometry/point.h>
19 
20 namespace RDKit {
21 namespace ForceFieldsHelper {
22  void RDKIT_FORCEFIELD_EXPORT normalizeAngleDeg(double &angleDeg);
23  void RDKIT_FORCEFIELD_EXPORT computeDihedral(const RDGeom::PointPtrVect &pos, unsigned int idx1,
24  unsigned int idx2, unsigned int idx3, unsigned int idx4,
25  double *dihedral = NULL, double *cosPhi = NULL,
26  RDGeom::Point3D r[4] = NULL, RDGeom::Point3D t[2] = NULL,
27  double d[2] = NULL);
28  void RDKIT_FORCEFIELD_EXPORT computeDihedral(const double *pos, unsigned int idx1,
29  unsigned int idx2, unsigned int idx3, unsigned int idx4,
30  double *dihedral = NULL, double *cosPhi = NULL,
31  RDGeom::Point3D r[4] = NULL, RDGeom::Point3D t[2] = NULL,
32  double d[2] = NULL);
34  const RDGeom::Point3D *p3, const RDGeom::Point3D *p4,
35  double *dihedral = NULL, double *cosPhi = NULL,
36  RDGeom::Point3D r[4] = NULL, RDGeom::Point3D t[2] = NULL,
37  double d[2] = NULL);
38 }
39 }
40 
41 namespace ForceFields {
43 typedef std::vector<int> INT_VECT;
44 typedef boost::shared_ptr<const ForceFieldContrib> ContribPtr;
45 typedef std::vector<ContribPtr> ContribPtrVect;
46 
47 //-------------------------------------------------------
48 //! A class to store forcefields and handle minimization
49 /*!
50  A force field is used like this (schematically):
51 
52  \verbatim
53  ForceField ff;
54 
55  // add contributions:
56  for contrib in contribs:
57  ff.contribs().push_back(contrib);
58 
59  // set up the points:
60  for positionPtr in positions:
61  ff.positions().push_back(point);
62 
63  // initialize:
64  ff.initialize()
65 
66  // and minimize:
67  needsMore = ff.minimize();
68 
69  \endverbatim
70 
71  <b>Notes:</b>
72  - The ForceField owns its contributions, which are stored using smart
73  pointers.
74  - Distance calculations are currently lazy; the full distance matrix is
75  never generated. In systems where the distance matrix is not sparse,
76  this is almost certainly inefficient.
77 
78 */
80  public:
81  //! construct with a dimension
82  ForceField(unsigned int dimension = 3)
83  : d_dimension(dimension), df_init(false), d_numPoints(0), dp_distMat(0){};
84 
85  ~ForceField();
86 
87  //! copy ctor, copies contribs.
88  ForceField(const ForceField &other);
89 
90  //! does initialization
91  void initialize();
92 
93  //! calculates and returns the energy (in kcal/mol) based on existing
94  // positions in the forcefield
95  /*!
96 
97  \return the current energy
98 
99  <b>Note:</b>
100  This function is less efficient than calcEnergy with postions passed in as
101  double *
102  the positions need to be converted to double * here
103  */
104  double calcEnergy(std::vector<double> *contribs = NULL) const;
105 
106  // these next two aren't const because they may update our
107  // distance matrix
108 
109  //! calculates and returns the energy of the position passed in
110  /*!
111  \param pos an array of doubles. Should be \c 3*this->numPoints() long.
112 
113  \return the current energy
114 
115  <b>Side effects:</b>
116  - Calling this resets the current distance matrix
117  - The individual contributions may further update the distance matrix
118  */
119  double calcEnergy(double *pos);
120 
121  //! calculates the gradient of the energy at the current position
122  /*!
123 
124  \param forces an array of doubles. Should be \c 3*this->numPoints() long.
125 
126  <b>Note:</b>
127  This function is less efficient than calcGrad with positions passed in
128  the positions need to be converted to double * here
129  */
130  void calcGrad(double *forces) const;
131 
132  //! calculates the gradient of the energy at the provided position
133  /*!
134 
135  \param pos an array of doubles. Should be \c 3*this->numPoints() long.
136  \param forces an array of doubles. Should be \c 3*this->numPoints() long.
137 
138  <b>Side effects:</b>
139  - The individual contributions may modify the distance matrix
140  */
141  void calcGrad(double *pos, double *forces);
142 
143  //! minimizes the energy of the system by following gradients
144  /*!
145  \param maxIts the maximum number of iterations to try
146  \param forceTol the convergence criterion for forces
147  \param energyTol the convergence criterion for energies
148 
149  \return an integer value indicating whether or not the convergence
150  criteria were achieved:
151  - 0: indicates success
152  - 1: the minimization did not converge in \c maxIts iterations.
153  */
154  int minimize(unsigned int snapshotFreq, RDKit::SnapshotVect *snapshotVect,
155  unsigned int maxIts = 200, double forceTol = 1e-4,
156  double energyTol = 1e-6);
157 
158  //! minimizes the energy of the system by following gradients
159  /*!
160  \param maxIts the maximum number of iterations to try
161  \param forceTol the convergence criterion for forces
162  \param energyTol the convergence criterion for energies
163  \param snapshotFreq a snapshot of the minimization trajectory
164  will be stored after as many steps as indicated
165  through this parameter; defaults to 0 (no
166  trajectory stored)
167  \param snapshotVect a pointer to a std::vector<Snapshot> where
168  coordinates and energies will be stored
169 
170  \return an integer value indicating whether or not the convergence
171  criteria were achieved:
172  - 0: indicates success
173  - 1: the minimization did not converge in \c maxIts iterations.
174  */
175  int minimize(unsigned int maxIts = 200, double forceTol = 1e-4,
176  double energyTol = 1e-6);
177 
178  // ---------------------------
179  // setters and getters
180 
181  //! returns a reference to our points (a PointPtrVect)
182  RDGeom::PointPtrVect &positions() { return d_positions; };
183  const RDGeom::PointPtrVect &positions() const { return d_positions; };
184 
185  //! returns a reference to our contribs (a ContribPtrVect)
186  ContribPtrVect &contribs() { return d_contribs; };
187  const ContribPtrVect &contribs() const { return d_contribs; };
188 
189  //! returns the distance between two points
190  /*!
191  \param i point index
192  \param j point index
193  \param pos (optional) If this argument is provided, it will be used
194  to provide the positions of points. \c pos should be
195  \c 3*this->numPoints() long.
196 
197  \return the distance
198 
199  <b>Side effects:</b>
200  - if the distance between i and j has not previously been calculated,
201  our internal distance matrix will be updated.
202  */
203  double distance(unsigned int i, unsigned int j, double *pos = 0);
204 
205  //! returns the distance between two points
206  /*!
207  \param i point index
208  \param j point index
209  \param pos (optional) If this argument is provided, it will be used
210  to provide the positions of points. \c pos should be
211  \c 3*this->numPoints() long.
212 
213  \return the distance
214 
215  <b>Note:</b>
216  The internal distance matrix is not updated in this case
217  */
218  double distance(unsigned int i, unsigned int j, double *pos = 0) const;
219 
220  //! returns the dimension of the forcefield
221  unsigned int dimension() const { return d_dimension; }
222 
223  //! returns the number of points the ForceField is handling
224  unsigned int numPoints() const { return d_numPoints; };
225 
226  INT_VECT &fixedPoints() { return d_fixedPoints; };
227  const INT_VECT &fixedPoints() const { return d_fixedPoints; };
228 
229  protected:
230  unsigned int d_dimension;
231  bool df_init; //!< whether or not we've been initialized
232  unsigned int d_numPoints; //!< the number of active points
233  double *dp_distMat; //!< our internal distance matrix
234  RDGeom::PointPtrVect d_positions; //!< pointers to the points we're using
235  ContribPtrVect d_contribs; //!< contributions to the energy
237  unsigned int d_matSize = 0;
238  //! scatter our positions into an array
239  /*!
240  \param pos should be \c 3*this->numPoints() long;
241  */
242  void scatter(double *pos) const;
243 
244  //! update our positions from an array
245  /*!
246  \param pos should be \c 3*this->numPoints() long;
247  */
248  void gather(double *pos);
249 
250  //! initializes our internal distance matrix
251  void initDistanceMatrix();
252 };
253 } // namespace ForceFields
254 #endif
ForceFields::ContribPtrVect
std::vector< ContribPtr > ContribPtrVect
Definition: ForceField.h:45
point.h
ForceFields::ForceField::contribs
const ContribPtrVect & contribs() const
Definition: ForceField.h:187
ForceFields::ContribPtr
boost::shared_ptr< const ForceFieldContrib > ContribPtr
Definition: ForceField.h:44
ForceFields::ForceField::dimension
unsigned int dimension() const
returns the dimension of the forcefield
Definition: ForceField.h:221
RDGeom::Point3D
Definition: point.h:46
Snapshot.h
ForceFields::ForceField::d_contribs
ContribPtrVect d_contribs
contributions to the energy
Definition: ForceField.h:235
ForceFields::ForceField
A class to store forcefields and handle minimization.
Definition: ForceField.h:79
ForceFields::ForceField::numPoints
unsigned int numPoints() const
returns the number of points the ForceField is handling
Definition: ForceField.h:224
RDKit::ForceFieldsHelper::computeDihedral
void RDKIT_FORCEFIELD_EXPORT computeDihedral(const RDGeom::PointPtrVect &pos, unsigned int idx1, unsigned int idx2, unsigned int idx3, unsigned int idx4, double *dihedral=NULL, double *cosPhi=NULL, RDGeom::Point3D r[4]=NULL, RDGeom::Point3D t[2]=NULL, double d[2]=NULL)
ForceFields::ForceField::d_positions
RDGeom::PointPtrVect d_positions
pointers to the points we're using
Definition: ForceField.h:234
ForceFields::ForceFieldContrib
abstract base class for contributions to ForceFields
Definition: Contrib.h:18
ForceFields::ForceField::fixedPoints
INT_VECT & fixedPoints()
Definition: ForceField.h:226
ForceFields::ForceField::contribs
ContribPtrVect & contribs()
returns a reference to our contribs (a ContribPtrVect)
Definition: ForceField.h:186
ForceFields::ForceField::d_fixedPoints
INT_VECT d_fixedPoints
Definition: ForceField.h:236
ForceFields::ForceField::positions
RDGeom::PointPtrVect & positions()
returns a reference to our points (a PointPtrVect)
Definition: ForceField.h:182
ForceFields::ForceField::dp_distMat
double * dp_distMat
our internal distance matrix
Definition: ForceField.h:233
RDKIT_FORCEFIELD_EXPORT
#define RDKIT_FORCEFIELD_EXPORT
Definition: export.h:255
ForceFields::ForceField::fixedPoints
const INT_VECT & fixedPoints() const
Definition: ForceField.h:227
ForceFields::INT_VECT
std::vector< int > INT_VECT
Definition: ForceField.h:42
RDKit::SnapshotVect
std::vector< Snapshot > SnapshotVect
Definition: Snapshot.h:19
ForceFields::ForceField::ForceField
ForceField(unsigned int dimension=3)
construct with a dimension
Definition: ForceField.h:82
ForceFields
Definition: TorsionAngleM6.h:24
BFGSOpt::minimize
int minimize(unsigned int dim, double *pos, double gradTol, unsigned int &numIters, double &funcVal, EnergyFunctor func, GradientFunctor gradFunc, unsigned int snapshotFreq, RDKit::SnapshotVect *snapshotVect, double funcTol=TOLX, unsigned int maxIts=MAXITS)
Do a BFGS minimization of a function.
Definition: BFGSOpt.h:187
ForceFields::ForceField::df_init
bool df_init
whether or not we've been initialized
Definition: ForceField.h:231
RDKit
Std stuff.
Definition: Atom.h:30
RDKit::ForceFieldsHelper::normalizeAngleDeg
void RDKIT_FORCEFIELD_EXPORT normalizeAngleDeg(double &angleDeg)
ForceFields::ForceField::d_numPoints
unsigned int d_numPoints
the number of active points
Definition: ForceField.h:232
ForceFields::ForceField::positions
const RDGeom::PointPtrVect & positions() const
Definition: ForceField.h:183
RDGeom::PointPtrVect
std::vector< RDGeom::Point * > PointPtrVect
Definition: point.h:499
export.h