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 ForceFields {
22 typedef std::vector<int> INT_VECT;
23 typedef boost::shared_ptr<const ForceFieldContrib> ContribPtr;
24 typedef std::vector<ContribPtr> ContribPtrVect;
25 
26 //-------------------------------------------------------
27 //! A class to store forcefields and handle minimization
28 /*!
29  A force field is used like this (schematically):
30 
31  \verbatim
32  ForceField ff;
33 
34  // add contributions:
35  for contrib in contribs:
36  ff.contribs().push_back(contrib);
37 
38  // set up the points:
39  for positionPtr in positions:
40  ff.positions().push_back(point);
41 
42  // initialize:
43  ff.initialize()
44 
45  // and minimize:
46  needsMore = ff.minimize();
47 
48  \endverbatim
49 
50  <b>Notes:</b>
51  - The ForceField owns its contributions, which are stored using smart
52  pointers.
53  - Distance calculations are currently lazy; the full distance matrix is
54  never generated. In systems where the distance matrix is not sparse,
55  this is almost certainly inefficient.
56 
57 */
59  public:
60  //! construct with a dimension
61  ForceField(unsigned int dimension = 3)
62  : d_dimension(dimension), df_init(false), d_numPoints(0), dp_distMat(0){};
63 
64  ~ForceField();
65 
66  //! copy ctor, copies contribs.
67  ForceField(const ForceField &other);
68 
69  //! does initialization
70  void initialize();
71 
72  //! calculates and returns the energy (in kcal/mol) based on existing
73  // positions in the forcefield
74  /*!
75 
76  \return the current energy
77 
78  <b>Note:</b>
79  This function is less efficient than calcEnergy with postions passed in as
80  double *
81  the positions need to be converted to double * here
82  */
83  double calcEnergy(std::vector<double> *contribs = NULL) const;
84 
85  // these next two aren't const because they may update our
86  // distance matrix
87 
88  //! calculates and returns the energy of the position passed in
89  /*!
90  \param pos an array of doubles. Should be \c 3*this->numPoints() long.
91 
92  \return the current energy
93 
94  <b>Side effects:</b>
95  - Calling this resets the current distance matrix
96  - The individual contributions may further update the distance matrix
97  */
98  double calcEnergy(double *pos);
99 
100  //! calculates the gradient of the energy at the current position
101  /*!
102 
103  \param forces an array of doubles. Should be \c 3*this->numPoints() long.
104 
105  <b>Note:</b>
106  This function is less efficient than calcGrad with positions passed in
107  the positions need to be converted to double * here
108  */
109  void calcGrad(double *forces) const;
110 
111  //! calculates the gradient of the energy at the provided position
112  /*!
113 
114  \param pos an array of doubles. Should be \c 3*this->numPoints() long.
115  \param forces an array of doubles. Should be \c 3*this->numPoints() long.
116 
117  <b>Side effects:</b>
118  - The individual contributions may modify the distance matrix
119  */
120  void calcGrad(double *pos, double *forces);
121 
122  //! minimizes the energy of the system by following gradients
123  /*!
124  \param maxIts the maximum number of iterations to try
125  \param forceTol the convergence criterion for forces
126  \param energyTol the convergence criterion for energies
127 
128  \return an integer value indicating whether or not the convergence
129  criteria were achieved:
130  - 0: indicates success
131  - 1: the minimization did not converge in \c maxIts iterations.
132  */
133  int minimize(unsigned int snapshotFreq, RDKit::SnapshotVect *snapshotVect,
134  unsigned int maxIts = 200, double forceTol = 1e-4,
135  double energyTol = 1e-6);
136 
137  //! minimizes the energy of the system by following gradients
138  /*!
139  \param maxIts the maximum number of iterations to try
140  \param forceTol the convergence criterion for forces
141  \param energyTol the convergence criterion for energies
142  \param snapshotFreq a snapshot of the minimization trajectory
143  will be stored after as many steps as indicated
144  through this parameter; defaults to 0 (no
145  trajectory stored)
146  \param snapshotVect a pointer to a std::vector<Snapshot> where
147  coordinates and energies will be stored
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 maxIts = 200, double forceTol = 1e-4,
155  double energyTol = 1e-6);
156 
157  // ---------------------------
158  // setters and getters
159 
160  //! returns a reference to our points (a PointPtrVect)
161  RDGeom::PointPtrVect &positions() { return d_positions; };
162  const RDGeom::PointPtrVect &positions() const { return d_positions; };
163 
164  //! returns a reference to our contribs (a ContribPtrVect)
165  ContribPtrVect &contribs() { return d_contribs; };
166  const ContribPtrVect &contribs() const { return d_contribs; };
167 
168  //! returns the distance between two points
169  /*!
170  \param i point index
171  \param j point index
172  \param pos (optional) If this argument is provided, it will be used
173  to provide the positions of points. \c pos should be
174  \c 3*this->numPoints() long.
175 
176  \return the distance
177 
178  <b>Side effects:</b>
179  - if the distance between i and j has not previously been calculated,
180  our internal distance matrix will be updated.
181  */
182  double distance(unsigned int i, unsigned int j, double *pos = 0);
183 
184  //! returns the distance between two points
185  /*!
186  \param i point index
187  \param j point index
188  \param pos (optional) If this argument is provided, it will be used
189  to provide the positions of points. \c pos should be
190  \c 3*this->numPoints() long.
191 
192  \return the distance
193 
194  <b>Note:</b>
195  The internal distance matrix is not updated in this case
196  */
197  double distance(unsigned int i, unsigned int j, double *pos = 0) const;
198 
199  //! returns the dimension of the forcefield
200  unsigned int dimension() const { return d_dimension; }
201 
202  //! returns the number of points the ForceField is handling
203  unsigned int numPoints() const { return d_numPoints; };
204 
205  INT_VECT &fixedPoints() { return d_fixedPoints; };
206  const INT_VECT &fixedPoints() const { return d_fixedPoints; };
207 
208  protected:
209  unsigned int d_dimension;
210  bool df_init; //!< whether or not we've been initialized
211  unsigned int d_numPoints; //!< the number of active points
212  double *dp_distMat; //!< our internal distance matrix
213  RDGeom::PointPtrVect d_positions; //!< pointers to the points we're using
214  ContribPtrVect d_contribs; //!< contributions to the energy
215  INT_VECT d_fixedPoints;
216  unsigned int d_matSize;
217  //! scatter our positions into an array
218  /*!
219  \param pos should be \c 3*this->numPoints() long;
220  */
221  void scatter(double *pos) const;
222 
223  //! update our positions from an array
224  /*!
225  \param pos should be \c 3*this->numPoints() long;
226  */
227  void gather(double *pos);
228 
229  //! initializes our internal distance matrix
230  void initDistanceMatrix();
231 };
232 }
233 #endif
INT_VECT & fixedPoints()
Definition: ForceField.h:205
unsigned int numPoints() const
returns the number of points the ForceField is handling
Definition: ForceField.h:203
const ContribPtrVect & contribs() const
Definition: ForceField.h:166
boost::shared_ptr< const ForceFieldContrib > ContribPtr
Definition: ForceField.h:23
ContribPtrVect & contribs()
returns a reference to our contribs (a ContribPtrVect)
Definition: ForceField.h:165
RDGeom::PointPtrVect d_positions
pointers to the points we&#39;re using
Definition: ForceField.h:213
unsigned int dimension() const
returns the dimension of the forcefield
Definition: ForceField.h:200
std::vector< RDGeom::Point * > PointPtrVect
Definition: point.h:493
ContribPtrVect d_contribs
contributions to the energy
Definition: ForceField.h:214
RDGeom::PointPtrVect & positions()
returns a reference to our points (a PointPtrVect)
Definition: ForceField.h:161
abstract base class for contributions to ForceFields
Definition: Contrib.h:18
unsigned int d_matSize
Definition: ForceField.h:216
std::vector< int > INT_VECT
Definition: ForceField.h:21
ForceField(unsigned int dimension=3)
construct with a dimension
Definition: ForceField.h:61
#define RDKIT_FORCEFIELD_EXPORT
Definition: export.h:242
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:210
const RDGeom::PointPtrVect & positions() const
Definition: ForceField.h:162
double * dp_distMat
our internal distance matrix
Definition: ForceField.h:212
A class to store forcefields and handle minimization.
Definition: ForceField.h:58
const INT_VECT & fixedPoints() const
Definition: ForceField.h:206
unsigned int d_numPoints
the number of active points
Definition: ForceField.h:211
std::vector< ContribPtr > ContribPtrVect
Definition: ForceField.h:24