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
236  INT_VECT d_fixedPoints;
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
INT_VECT & fixedPoints()
Definition: ForceField.h:226
unsigned int numPoints() const
returns the number of points the ForceField is handling
Definition: ForceField.h:224
const ContribPtrVect & contribs() const
Definition: ForceField.h:187
boost::shared_ptr< const ForceFieldContrib > ContribPtr
Definition: ForceField.h:44
ContribPtrVect & contribs()
returns a reference to our contribs (a ContribPtrVect)
Definition: ForceField.h:186
RDGeom::PointPtrVect d_positions
pointers to the points we&#39;re using
Definition: ForceField.h:234
unsigned int dimension() const
returns the dimension of the forcefield
Definition: ForceField.h:221
std::vector< RDGeom::Point * > PointPtrVect
Definition: point.h:499
ContribPtrVect d_contribs
contributions to the energy
Definition: ForceField.h:235
RDGeom::PointPtrVect & positions()
returns a reference to our points (a PointPtrVect)
Definition: ForceField.h:182
abstract base class for contributions to ForceFields
Definition: Contrib.h:18
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)
ForceField(unsigned int dimension=3)
construct with a dimension
Definition: ForceField.h:82
#define RDKIT_FORCEFIELD_EXPORT
Definition: export.h:255
std::vector< int > INT_VECT
Definition: types.h:254
Std stuff.
Definition: Atom.h:30
std::vector< Snapshot > SnapshotVect
Definition: Snapshot.h:19
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
bool df_init
whether or not we&#39;ve been initialized
Definition: ForceField.h:231
const RDGeom::PointPtrVect & positions() const
Definition: ForceField.h:183
double * dp_distMat
our internal distance matrix
Definition: ForceField.h:233
A class to store forcefields and handle minimization.
Definition: ForceField.h:79
const INT_VECT & fixedPoints() const
Definition: ForceField.h:227
void RDKIT_FORCEFIELD_EXPORT normalizeAngleDeg(double &angleDeg)
unsigned int d_numPoints
the number of active points
Definition: ForceField.h:232
std::vector< ContribPtr > ContribPtrVect
Definition: ForceField.h:45