RDKit
Open-source cheminformatics and machine learning.
FFConvenience.h
Go to the documentation of this file.
1 //
2 // Copyright (C) 2019 Paolo Tosco
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_FFCONVENIENCE_H
12 #define RD_FFCONVENIENCE_H
13 #include <ForceField/ForceField.h>
14 #include <RDGeneral/RDThreads.h>
15 
16 namespace RDKit {
17 class ROMol;
18 namespace ForceFieldsHelper {
19 namespace detail {
20 #ifdef RDK_BUILD_THREADSAFE_SSS
21 void OptimizeMoleculeConfsHelper_(ForceFields::ForceField ff, ROMol *mol,
22  std::vector<std::pair<int, double>> *res,
23  unsigned int threadIdx,
24  unsigned int numThreads, int maxIters) {
25  PRECONDITION(mol, "mol must not be nullptr");
26  PRECONDITION(res, "res must not be nullptr");
27  PRECONDITION(res->size() >= mol->getNumConformers(),
28  "res->size() must be >= mol->getNumConformers()");
29  unsigned int i = 0;
30  ff.positions().resize(mol->getNumAtoms());
31  for (ROMol::ConformerIterator cit = mol->beginConformers();
32  cit != mol->endConformers(); ++cit, ++i) {
33  if (i % numThreads != threadIdx) {
34  continue;
35  }
36  for (unsigned int aidx = 0; aidx < mol->getNumAtoms(); ++aidx) {
37  ff.positions()[aidx] = &(*cit)->getAtomPos(aidx);
38  }
39  ff.initialize();
40  int needsMore = ff.minimize(maxIters);
41  double e = ff.calcEnergy();
42  (*res)[i] = std::make_pair(needsMore, e);
43  }
44 }
45 
46 void OptimizeMoleculeConfsMT(ROMol &mol, const ForceFields::ForceField &ff,
47  std::vector<std::pair<int, double>> &res,
48  int numThreads, int maxIters) {
49  std::vector<std::thread> tg;
50  for (int ti = 0; ti < numThreads; ++ti) {
51  tg.emplace_back(std::thread(detail::OptimizeMoleculeConfsHelper_, ff, &mol,
52  &res, ti, numThreads, maxIters));
53  }
54  for (auto &thread : tg) {
55  if (thread.joinable()) {
56  thread.join();
57  }
58  }
59 }
60 #endif
61 
63  std::vector<std::pair<int, double>> &res,
64  int maxIters) {
65  PRECONDITION(res.size() >= mol.getNumConformers(),
66  "res.size() must be >= mol.getNumConformers()");
67  unsigned int i = 0;
68  for (ROMol::ConformerIterator cit = mol.beginConformers();
69  cit != mol.endConformers(); ++cit, ++i) {
70  for (unsigned int aidx = 0; aidx < mol.getNumAtoms(); ++aidx) {
71  ff.positions()[aidx] = &(*cit)->getAtomPos(aidx);
72  }
73  ff.initialize();
74  int needsMore = ff.minimize(maxIters);
75  double e = ff.calcEnergy();
76  res[i] = std::make_pair(needsMore, e);
77  }
78 }
79 } // namespace detail
80 
81 //! Convenience function for optimizing a molecule using a pre-generated
82 //! force-field
83 /*
84  \param ff the force-field
85  \param res vector of (needsMore,energy) pairs
86  \param maxIters the maximum number of force-field iterations
87 
88  \return a pair with:
89  first: -1 if parameters were missing, 0 if the optimization converged, 1 if
90  more iterations are required.
91  second: the energy
92 */
93 std::pair<int, double> OptimizeMolecule(ForceFields::ForceField &ff,
94  int maxIters = 1000) {
95  ff.initialize();
96  int res = ff.minimize(maxIters);
97  double e = ff.calcEnergy();
98  return std::make_pair(res, e);
99 }
100 
101 //! Convenience function for optimizing all of a molecule's conformations using
102 /// a pre-generated force-field
103 /*
104  \param mol the molecule to use
105  \param ff the force-field
106  \param res vector of (needsMore,energy) pairs
107  \param numThreads the number of simultaneous threads to use (only has an
108  effect if the RDKit is compiled with thread support).
109  If set to zero, the max supported by the system will be
110  used.
111  \param maxIters the maximum number of force-field iterations
112 
113 */
115  std::vector<std::pair<int, double>> &res,
116  int numThreads = 1, int maxIters = 1000) {
117  res.resize(mol.getNumConformers());
118  numThreads = getNumThreadsToUse(numThreads);
119  if (numThreads == 1) {
120  detail::OptimizeMoleculeConfsST(mol, ff, res, maxIters);
121  }
122 #ifdef RDK_BUILD_THREADSAFE_SSS
123  else {
124  detail::OptimizeMoleculeConfsMT(mol, ff, res, numThreads, maxIters);
125  }
126 #endif
127 }
128 } // end of namespace ForceFieldsHelper
129 } // end of namespace RDKit
130 #endif
#define PRECONDITION(expr, mess)
Definition: Invariant.h:109
A class to store forcefields and handle minimization.
Definition: ForceField.h:79
double calcEnergy(std::vector< double > *contribs=nullptr) const
void initialize()
does initialization
RDGeom::PointPtrVect & positions()
returns a reference to our points (a PointPtrVect)
Definition: ForceField.h:181
int minimize(unsigned int snapshotFreq, RDKit::SnapshotVect *snapshotVect, unsigned int maxIts=200, double forceTol=1e-4, double energyTol=1e-6)
minimizes the energy of the system by following gradients
unsigned int getNumConformers() const
Definition: ROMol.h:560
unsigned int getNumAtoms() const
returns our number of atoms
Definition: ROMol.h:413
ConformerIterator beginConformers()
Definition: ROMol.h:732
ConformerIterator endConformers()
Definition: ROMol.h:734
void OptimizeMoleculeConfsST(ROMol &mol, ForceFields::ForceField &ff, std::vector< std::pair< int, double >> &res, int maxIters)
Definition: FFConvenience.h:62
std::pair< int, double > OptimizeMolecule(ForceFields::ForceField &ff, int maxIters=1000)
Definition: FFConvenience.h:93
void OptimizeMoleculeConfs(ROMol &mol, ForceFields::ForceField &ff, std::vector< std::pair< int, double >> &res, int numThreads=1, int maxIters=1000)
Std stuff.
Definition: Abbreviations.h:19
unsigned int getNumThreadsToUse(int target)
Definition: RDThreads.h:37