Skip to content
Permalink
96353c3a8b
Switch branches/tags
Go to file
 
 
Cannot retrieve contributors at this time
134 lines (122 sloc) 3.99 KB
//
// 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 {
// 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) {
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);
}
}
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) {
std::unique_ptr<ROMol> mol(MolOps::removeHs(m));
std::vector<double> masses;
std::vector<double> charges;
unsigned int num_atoms = mol->getNumAtoms();
masses.reserve(num_atoms);
charges.reserve(num_atoms);
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(num_atoms, 0.0);
std::vector<double> cmr(num_atoms, 0.0);
getCrippenAtomContribs(*mol, slogp, cmr);
// polarizability? - need model
// slogp? sasa?
auto burden = make_burden(m);
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;
}
}
}
#endif