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_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) continue;
34  for (unsigned int aidx = 0; aidx < mol->getNumAtoms(); ++aidx) {
35  ff.positions()[aidx] = &(*cit)->getAtomPos(aidx);
36  }
37  ff.initialize();
38  int needsMore = ff.minimize(maxIters);
39  double e = ff.calcEnergy();
40  (*res)[i] = std::make_pair(needsMore, e);
41  }
42 }
43 
44 void OptimizeMoleculeConfsMT(ROMol &mol, const ForceFields::ForceField &ff,
45  std::vector<std::pair<int, double>> &res,
46  int numThreads, int maxIters) {
47  std::vector<std::thread> tg;
48  for (int ti = 0; ti < numThreads; ++ti) {
49  tg.emplace_back(std::thread(detail::OptimizeMoleculeConfsHelper_, ff, &mol,
50  &res, ti, numThreads, maxIters));
51  }
52  for (auto &thread : tg) {
53  if (thread.joinable()) thread.join();
54  }
55 }
56 #endif
57 
59  std::vector<std::pair<int, double>> &res,
60  int maxIters) {
61  PRECONDITION(res.size() >= mol.getNumConformers(),
62  "res.size() must be >= mol.getNumConformers()");
63  unsigned int i = 0;
64  for (ROMol::ConformerIterator cit = mol.beginConformers();
65  cit != mol.endConformers(); ++cit, ++i) {
66  for (unsigned int aidx = 0; aidx < mol.getNumAtoms(); ++aidx) {
67  ff.positions()[aidx] = &(*cit)->getAtomPos(aidx);
68  }
69  ff.initialize();
70  int needsMore = ff.minimize(maxIters);
71  double e = ff.calcEnergy();
72  res[i] = std::make_pair(needsMore, e);
73  }
74 }
75 } // namespace detail
76 
77 //! Convenience function for optimizing a molecule using a pre-generated
78 //! force-field
79 /*
80  \param ff the force-field
81  \param res vector of (needsMore,energy) pairs
82  \param maxIters the maximum number of force-field iterations
83 
84  \return a pair with:
85  first: -1 if parameters were missing, 0 if the optimization converged, 1 if
86  more iterations are required.
87  second: the energy
88 */
89 std::pair<int, double> OptimizeMolecule(ForceFields::ForceField &ff,
90  int maxIters = 1000) {
91  ff.initialize();
92  int res = ff.minimize(maxIters);
93  double e = ff.calcEnergy();
94  return std::make_pair(res, e);
95 }
96 
97 //! Convenience function for optimizing all of a molecule's conformations using
98 /// a pre-generated force-field
99 /*
100  \param mol the molecule to use
101  \param ff the force-field
102  \param res vector of (needsMore,energy) pairs
103  \param numThreads the number of simultaneous threads to use (only has an
104  effect if the RDKit is compiled with thread support).
105  If set to zero, the max supported by the system will be
106  used.
107  \param maxIters the maximum number of force-field iterations
108 
109 */
111  std::vector<std::pair<int, double>> &res,
112  int numThreads = 1, int maxIters = 1000) {
113  res.resize(mol.getNumConformers());
114  numThreads = getNumThreadsToUse(numThreads);
115  if (numThreads == 1) {
116  detail::OptimizeMoleculeConfsST(mol, ff, res, maxIters);
117  }
118 #ifdef RDK_THREADSAFE_SSS
119  else {
120  detail::OptimizeMoleculeConfsMT(mol, ff, res, numThreads, maxIters);
121  }
122 #endif
123 }
124 } // end of namespace ForceFieldsHelper
125 } // end of namespace RDKit
126 #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:479
unsigned int getNumAtoms() const
returns our number of atoms
Definition: ROMol.h:332
ConformerIterator beginConformers()
Definition: ROMol.h:651
ConformerIterator endConformers()
Definition: ROMol.h:653
void OptimizeMoleculeConfsST(ROMol &mol, ForceFields::ForceField &ff, std::vector< std::pair< int, double >> &res, int maxIters)
Definition: FFConvenience.h:58
std::pair< int, double > OptimizeMolecule(ForceFields::ForceField &ff, int maxIters=1000)
Definition: FFConvenience.h:89
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:18
unsigned int getNumThreadsToUse(int target)
Definition: RDThreads.h:37