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
|
//
// Copyright (C) 2020 Brian P. Kelley
//
// @@ 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.
//
#ifdef RDK_HAS_EIGEN3
#include "BCUT.h"
#include "Crippen.h"
#include <Eigen/Dense>
#include <GraphMol/RDKitBase.h>
#include "GraphMol/PartialCharges/GasteigerCharges.h"
#include "GraphMol/PartialCharges/GasteigerParams.h"
#include <RDGeneral/types.h>
namespace RDKit {
namespace Descriptors {
constexpr int NUM_BCUTS = 8;
// diagonal elements are a property (atomic num, charge, etc)
// off diagonal are 1/sqrt(bond_order)
// Original burden matrix was .1, .2, .3, .15 for single,double,triple or
// aromatic all other elements are .001
namespace {
std::unique_ptr<Eigen::MatrixXd> make_burden(const ROMol &m) {
auto num_atoms = m.getNumAtoms();
std::unique_ptr<Eigen::MatrixXd> burden(
new Eigen::MatrixXd(num_atoms, num_atoms));
for (unsigned int i = 0; i < num_atoms; ++i) {
for (unsigned int j = 0; j < num_atoms; ++j) {
(*burden)(i, j) = (*burden)(j, i) = 0.001;
}
}
for (auto &bond : m.bonds()) {
unsigned int i = bond->getBeginAtomIdx();
unsigned int j = bond->getEndAtomIdx();
double score = 0.0;
switch (bond->getBondType()) {
case Bond::AROMATIC:
// score = 0.15; orig burden
score = 0.8164965809277261; // 1/sqrt(1.5)
break;
case Bond::SINGLE:
// score = 0.1;
score = 1.0; // 1/sqrt(1.0)
break;
case Bond::DOUBLE:
// score = 0.2;
score = 0.7071067811865475; // 1/sqrt(2.0)
break;
case Bond::TRIPLE:
// score = 0.3;
score = 0.5773502691896258; // 1/sqrt(3);
break;
default:
CHECK_INVARIANT(
0, "Bond order must be Single, Double, Triple or Aromatic");
}
(*burden)(i, j) = (*burden)(j, i) = score;
}
return burden;
}
std::pair<double, double> BCUT2D(std::unique_ptr<Eigen::MatrixXd> &burden,
const std::vector<double> &atom_props) {
if (atom_props.empty()) {
return std::pair<double, double>(0.0, 0.0);
}
for (unsigned int i = 0; i < atom_props.size(); ++i) {
(*burden)(i, i) = atom_props[i];
}
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(*burden);
auto eivals = es.eigenvalues();
double lowest = eivals(0);
double highest = eivals(atom_props.size() - 1);
return std::pair<double, double>(highest, lowest);
}
} // namespace
std::pair<double, double> BCUT2D(const ROMol &m,
const std::vector<double> &atom_props) {
unsigned int num_atoms = m.getNumAtoms();
PRECONDITION(atom_props.size() == num_atoms,
"Number of atom props not equal to number of atoms");
if (num_atoms == 0) {
return std::pair<double, double>(0, 0);
}
auto burden = make_burden(m);
return BCUT2D(burden, atom_props);
}
std::pair<double, double> BCUT2D(const ROMol &m,
const std::string &atom_double_prop) {
std::vector<double> props;
props.reserve(m.getNumAtoms());
for (auto &atom : m.atoms()) {
props.push_back(atom->getProp<double>(atom_double_prop));
}
return BCUT2D(m, props);
}
std::vector<double> BCUT2D(const ROMol &m) {
if (!m.getNumAtoms()) {
return std::vector<double>(NUM_BCUTS, 0.0);
}
std::unique_ptr<ROMol> mol(MolOps::removeAllHs(m));
unsigned int numAtoms = mol->getNumAtoms();
std::vector<double> masses;
std::vector<double> charges;
masses.reserve(numAtoms);
charges.reserve(numAtoms);
RDKit::computeGasteigerCharges(*mol, 12, true);
for (auto &atom : mol->atoms()) {
masses.push_back(atom->getMass());
charges.push_back(
atom->getProp<double>(common_properties::_GasteigerCharge));
}
std::vector<double> slogp(numAtoms, 0.0);
std::vector<double> cmr(numAtoms, 0.0);
getCrippenAtomContribs(*mol, slogp, cmr);
// polarizability? - need model
// slogp? sasa?
auto burden = make_burden(*mol);
auto atom_bcut = BCUT2D(burden, masses);
auto gasteiger = BCUT2D(burden, charges);
auto logp = BCUT2D(burden, slogp);
auto mr = BCUT2D(burden, cmr);
std::vector<double> res = {
atom_bcut.first, atom_bcut.second, gasteiger.first, gasteiger.second,
logp.first, logp.second, mr.first, mr.second};
return res;
}
} // namespace Descriptors
} // namespace RDKit
#endif
|
