1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
//
// Copyright (C) 2019 Paolo Tosco
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#include <RDGeneral/export.h>
#ifndef RD_FFCONVENIENCE_H
#define RD_FFCONVENIENCE_H
#include <ForceField/ForceField.h>
#include <RDGeneral/RDThreads.h>
namespace RDKit {
class ROMol;
namespace ForceFieldsHelper {
namespace detail {
#ifdef RDK_BUILD_THREADSAFE_SSS
void OptimizeMoleculeConfsHelper_(ForceFields::ForceField ff, ROMol *mol,
std::vector<std::pair<int, double>> *res,
unsigned int threadIdx,
unsigned int numThreads, int maxIters) {
PRECONDITION(mol, "mol must not be nullptr");
PRECONDITION(res, "res must not be nullptr");
PRECONDITION(res->size() >= mol->getNumConformers(),
"res->size() must be >= mol->getNumConformers()");
unsigned int i = 0;
ff.positions().resize(mol->getNumAtoms());
for (ROMol::ConformerIterator cit = mol->beginConformers();
cit != mol->endConformers(); ++cit, ++i) {
if (i % numThreads != threadIdx) {
continue;
}
for (unsigned int aidx = 0; aidx < mol->getNumAtoms(); ++aidx) {
ff.positions()[aidx] = &(*cit)->getAtomPos(aidx);
}
ff.initialize();
int needsMore = ff.minimize(maxIters);
double e = ff.calcEnergy();
(*res)[i] = std::make_pair(needsMore, e);
}
}
void OptimizeMoleculeConfsMT(ROMol &mol, const ForceFields::ForceField &ff,
std::vector<std::pair<int, double>> &res,
int numThreads, int maxIters) {
std::vector<std::thread> tg;
for (int ti = 0; ti < numThreads; ++ti) {
tg.emplace_back(std::thread(detail::OptimizeMoleculeConfsHelper_, ff, &mol,
&res, ti, numThreads, maxIters));
}
for (auto &thread : tg) {
if (thread.joinable()) {
thread.join();
}
}
}
#endif
void OptimizeMoleculeConfsST(ROMol &mol, ForceFields::ForceField &ff,
std::vector<std::pair<int, double>> &res,
int maxIters) {
PRECONDITION(res.size() >= mol.getNumConformers(),
"res.size() must be >= mol.getNumConformers()");
unsigned int i = 0;
for (ROMol::ConformerIterator cit = mol.beginConformers();
cit != mol.endConformers(); ++cit, ++i) {
for (unsigned int aidx = 0; aidx < mol.getNumAtoms(); ++aidx) {
ff.positions()[aidx] = &(*cit)->getAtomPos(aidx);
}
ff.initialize();
int needsMore = ff.minimize(maxIters);
double e = ff.calcEnergy();
res[i] = std::make_pair(needsMore, e);
}
}
} // namespace detail
//! Convenience function for optimizing a molecule using a pre-generated
//! force-field
/*
\param ff the force-field
\param res vector of (needsMore,energy) pairs
\param maxIters the maximum number of force-field iterations
\return a pair with:
first: -1 if parameters were missing, 0 if the optimization converged, 1 if
more iterations are required.
second: the energy
*/
std::pair<int, double> OptimizeMolecule(ForceFields::ForceField &ff,
int maxIters = 1000) {
ff.initialize();
int res = ff.minimize(maxIters);
double e = ff.calcEnergy();
return std::make_pair(res, e);
}
//! Convenience function for optimizing all of a molecule's conformations using
/// a pre-generated force-field
/*
\param mol the molecule to use
\param ff the force-field
\param res vector of (needsMore,energy) pairs
\param numThreads the number of simultaneous threads to use (only has an
effect if the RDKit is compiled with thread support).
If set to zero, the max supported by the system will be
used.
\param maxIters the maximum number of force-field iterations
*/
void OptimizeMoleculeConfs(ROMol &mol, ForceFields::ForceField &ff,
std::vector<std::pair<int, double>> &res,
int numThreads = 1, int maxIters = 1000) {
res.resize(mol.getNumConformers());
numThreads = getNumThreadsToUse(numThreads);
if (numThreads == 1) {
detail::OptimizeMoleculeConfsST(mol, ff, res, maxIters);
}
#ifdef RDK_BUILD_THREADSAFE_SSS
else {
detail::OptimizeMoleculeConfsMT(mol, ff, res, numThreads, maxIters);
}
#endif
}
} // end of namespace ForceFieldsHelper
} // end of namespace RDKit
#endif
|
