forked from OpenChemistry/avogadrolibs
-
Notifications
You must be signed in to change notification settings - Fork 2
/
obenergy.cpp
185 lines (155 loc) · 4.99 KB
/
obenergy.cpp
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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
/******************************************************************************
This source file is part of the Avogadro project.
This source code is released under the 3-Clause BSD License, (see "LICENSE").
It links to the Open Babel library, which is released under the GNU GPL v2.
******************************************************************************/
#include "obenergy.h"
#include <avogadro/core/molecule.h>
#include <openbabel/babelconfig.h>
#include <openbabel/atom.h>
#include <openbabel/base.h>
#include <openbabel/forcefield.h>
#include <openbabel/math/vector3.h>
#include <openbabel/mol.h>
#include <openbabel/obconversion.h>
#include <openbabel/obiter.h>
#include <QDebug>
#include <QDir>
using namespace OpenBabel;
namespace Avogadro::QtPlugins {
class OBEnergy::Private
{
public:
// OBMol and OBForceField are owned by this class
OBMol* m_obmol = nullptr;
OBForceField* m_forceField = nullptr;
~Private()
{
delete m_obmol;
delete m_forceField;
}
};
OBEnergy::OBEnergy(const std::string& method)
: m_identifier(method), m_name(method), m_molecule(nullptr)
{
d = new Private;
// Ensure the plugins are loaded
OBPlugin::LoadAllPlugins();
d->m_forceField = static_cast<OBForceField*>(
OBPlugin::GetPlugin("forcefields", method.c_str()));
if (method == "UFF") {
m_description = tr("Universal Force Field");
m_elements.reset();
for (unsigned int i = 1; i < 102; ++i)
m_elements.set(i);
} else if (method == "GAFF") {
m_description = tr("Generalized Amber Force Field");
// H, C, N, O, F, P, S, Cl, Br, and I
m_elements.set(1);
m_elements.set(6);
m_elements.set(7);
m_elements.set(8);
m_elements.set(9);
m_elements.set(15);
m_elements.set(16);
m_elements.set(17);
m_elements.set(35);
m_elements.set(53);
} else if (method == "MMFF94") {
m_description = tr("Merck Molecular Force Field 94");
m_elements.reset();
// H, C, N, O, F, Si, P, S, Cl, Br, and I
m_elements.set(1);
m_elements.set(6);
m_elements.set(7);
m_elements.set(8);
m_elements.set(9);
m_elements.set(14);
m_elements.set(15);
m_elements.set(16);
m_elements.set(17);
m_elements.set(35);
m_elements.set(53);
}
}
OBEnergy::~OBEnergy() {}
Calc::EnergyCalculator* OBEnergy::newInstance() const
{
return new OBEnergy(m_name);
}
void OBEnergy::setMolecule(Core::Molecule* mol)
{
m_molecule = mol;
if (mol == nullptr || mol->atomCount() == 0) {
return; // nothing to do
}
// set up our internal OBMol
d->m_obmol = new OBMol;
// copy the atoms, bonds, and coordinates
d->m_obmol->BeginModify();
for (size_t i = 0; i < mol->atomCount(); ++i) {
const Core::Atom& atom = mol->atom(i);
OBAtom* obAtom = d->m_obmol->NewAtom();
obAtom->SetAtomicNum(atom.atomicNumber());
auto pos = atom.position3d().cast<double>();
obAtom->SetVector(pos.x(), pos.y(), pos.z());
}
for (size_t i = 0; i < mol->bondCount(); ++i) {
const Core::Bond& bond = mol->bond(i);
d->m_obmol->AddBond(bond.atom1().index() + 1, bond.atom2().index() + 1,
bond.order());
}
d->m_obmol->EndModify();
// make sure we can set up the force field
if (d->m_forceField != nullptr) {
d->m_forceField->Setup(*d->m_obmol);
} else {
d->m_forceField = static_cast<OBForceField*>(
OBPlugin::GetPlugin("forcefields", m_identifier.c_str()));
if (d->m_forceField != nullptr) {
d->m_forceField->Setup(*d->m_obmol);
}
}
}
Real OBEnergy::value(const Eigen::VectorXd& x)
{
if (m_molecule == nullptr || m_molecule->atomCount() == 0)
return 0.0; // nothing to do
// update coordinates in our private OBMol
for (size_t i = 0; i < m_molecule->atomCount(); ++i) {
Eigen::Vector3d pos(x[i * 3], x[i * 3 + 1], x[i * 3 + 2]);
d->m_obmol->GetAtom(i + 1)->SetVector(pos.x(), pos.y(), pos.z());
}
double energy = 0.0;
if (d->m_forceField != nullptr) {
d->m_forceField->SetCoordinates(*d->m_obmol);
energy = d->m_forceField->Energy(false);
}
return energy;
}
void OBEnergy::gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad)
{
if (m_molecule == nullptr || m_molecule->atomCount() == 0)
return;
// update coordinates in our private OBMol
for (size_t i = 0; i < m_molecule->atomCount(); ++i) {
Eigen::Vector3d pos(x[i * 3], x[i * 3 + 1], x[i * 3 + 2]);
d->m_obmol->GetAtom(i + 1)->SetVector(pos.x(), pos.y(), pos.z());
}
if (d->m_forceField != nullptr) {
d->m_forceField->SetCoordinates(*d->m_obmol);
// make sure gradients are calculated
double energy = d->m_forceField->Energy(true);
for (size_t i = 0; i < m_molecule->atomCount(); ++i) {
OBAtom* atom = d->m_obmol->GetAtom(i + 1);
OpenBabel::vector3 obGrad = d->m_forceField->GetGradient(atom);
grad[3 * i] = obGrad.x();
grad[3 * i + 1] = obGrad.y();
grad[3 * i + 2] = obGrad.z();
}
grad *= -1; // OpenBabel outputs forces, not grads
cleanGradients(grad);
}
}
} // namespace Avogadro::QtPlugins
#include "obenergy.moc"