Dear Paul,
I'm afraid currently there is no public method to get r0 and kb from
ForceField contribs. I have attached a C++ sample program which shows
how to get the information you need from a ROMol. Please feel free to
contact me if you need further help.
Cheers,
p.
On 09/29/2014 07:36 PM, Paul Emsley wrote:
Dear Paolo (I presume :-),
I'd like to get the ro values for the bonds in an MMFF. I have the
following so far:
void mmff_bonds(RDKit::ROMol mol_in) {
RDKit::MMFF::MMFFMolProperties *mmffMolProperties = new
RDKit::MMFF::MMFFMolProperties(mol_in);
ForceFields::ForceField *field = new ForceFields::ForceField();
RDKit::MMFF::Tools::addBonds( mol_in, mmffMolProperties, field);
ForceFields::ContribPtrVect::const_iterator contrib;
for(contrib=field-contribs().begin();
contrib!=field-contribs().end();
contrib++) {
ForceFields::MMFF::MMFFBond *mmffBondParams;
// something xxx!!!
double d = something(contrib???)-calcBondRestLength(mmffBondParams);
}
}
As you can see, I got lost between contrib and bond parameters. Please
can you tell me what I am missing? I hope it's only a few lines - but I
can't determine which ones :)
Thanks,
Paul.
--
Slashdot TV. Videos for Nerds. Stuff that Matters.
http://pubads.g.doubleclick.net/gampad/clk?id=160591471iu=/4140/ostg.clktrk
___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
--
==
Paolo Tosco, Ph.D.
Department of Drug Science and Technology
Via Pietro Giuria, 9 - 10125 Torino (Italy)
Tel: +39 011 670 7680 | Mob: +39 348 5537206
Fax: +39 011 670 7687 | E-mail: paolo.to...@unito.it
http://open3dqsar.org | http://open3dalign.org
==
/*
Build with:
g++ -Wall -o get_r0_kb get_r0_kb.cpp \
-I$BOOST_INCLUDE -I$RDBASE/Code -L$BOOST_LIB -L$RDBASE/lib \
-lFileParsers -lRDGeneral -lForceFieldHelpers
*/
#include GraphMol/RDKitBase.h
#include GraphMol/FileParsers/FileParsers.h
#include ForceField/MMFF/Params.h
#include ForceField/MMFF/Contribs.h
#include ForceField/MMFF/BondStretch.h
#include GraphMol/ForceFieldHelpers/MMFF/AtomTyper.h
#include GraphMol/ForceFieldHelpers/MMFF/Builder.h
#include ForceField/ForceField.h
#include string
#include iostream
using namespace RDKit;
int main(int argc, char **argv)
{
std::stringstream ss;
if (argc 2) {
std::cerr Usage: get_r0_kb molFile std::endl;
return 0;
}
std::string sdf(argv[1]);
ROMol *mol = static_castROMol *(MolFileToMol(sdf, true, false));
MMFF::MMFFMolProperties *mp = new MMFF::MMFFMolProperties(*mol);
if (!(mp-isValid())) {
std::cerr Error setting up force-field std::endl;
}
else {
ForceFields::MMFF::MMFFBondCollection *mmffBond =
ForceFields::MMFF::MMFFBondCollection::getMMFFBond();
std::cout std::setw(8) idx1 std::setw(8) idx2
std::setw(8) r0 std::setw(8) kb std::endl;
for (ROMol::BondIterator bi = mol-beginBonds();
bi != mol-endBonds(); ++bi) {
unsigned int idx1 = (*bi)-getBeginAtomIdx();
unsigned int idx2 = (*bi)-getEndAtomIdx();
unsigned int iAtomType = mp-getMMFFAtomType(idx1);
unsigned int jAtomType = mp-getMMFFAtomType(idx2);
unsigned int bondType = mp-getMMFFBondType(*bi);
bool areMMFFBondParamsEmpirical = false;
const ForceFields::MMFF::MMFFBond *mmffBondParams =
(*mmffBond)(bondType, iAtomType, jAtomType);
if (!mmffBondParams) {
mmffBondParams = mp-getMMFFBondStretchEmpiricalRuleParams(*mol, *bi);
areMMFFBondParamsEmpirical = true;
}
double r0 = ForceFields::MMFF::Utils::calcBondRestLength(mmffBondParams);
double kb = ForceFields::MMFF::Utils::calcBondForceConstant(mmffBondParams);
if (areMMFFBondParamsEmpirical) {
delete mmffBondParams;
}
std::cout std::setw(8) idx1 + 1 std::setw(8) idx2 + 1
std::setw(8) std::fixed std::setprecision(4) r0
std::setw(8) std::fixed std::setprecision(4) kb
std::endl;
}
}
if (mp) delete mp;
if (mol) delete mol;
return 0;
}
--
Slashdot TV. Videos for Nerds. Stuff that Matters.
http://pubads.g.doubleclick.net/gampad/clk?id=160591471iu=/4140/ostg.clktrk___
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
