thank you Craig and Geoff for your answers to my first question.
I'm still having problems though, possibly related to the difference
between implicit and explicit hydrogens and/or to the less usual atoms that
I'm using ?
In my code below I create three molecules [Si][Si], [Si]=[Si] and [Si]#[Si]
first by hand and then the same three by reading from smiles. For all six
molecules I output smiles, inchikey and energy. The following worries me:
- The third smilescode is output as [SiH3]#[SiH3] instead of
[SiH]#[SiH], i.e. too many implicit hydrogens are shown. The molecule seems
correct though, since the Inchi-key is identical to the one created for
[SiH]#[SiH] in another package (CDK).
- The inchi-keys generated for the second set of molecules seem
incorrect and their generation throws warnings. The energies seem ok, i.e.
are identical to those calculated for the first set.
- If I substitute Silicium by Carbon (i.e. creating ethane, ethene,
ethyn), all smiles and inchikeys seem consistent, but the energies
calculated for the second set are ~10x higher.
Btw my target is to create the molecules by hand and calculate the energies
as an estimate for the enthalpy of formation. I know this may be
problematic but a discussion on calculating heat of formation is probably a
different post ...
regards,
Rafel
#include <iostream>
#include <openbabel/obconversion.h>
#include <openbabel/mol.h>
#include <openbabel/forcefield.h>
#include <openbabel/conformersearch.h>
#include <openbabel/builder.h>
OpenBabel::OBMol obMol;
OpenBabel::OBBuilder builder;
OpenBabel::OBForceField *ff;
OpenBabel::OBConversion conv, convToInchi;
void calcEnergy_and_output() {
conv.Write(&obMol);
convToInchi.Write(&obMol);
builder.Build(obMol);
if (!ff->Setup(obMol)) std::cout << "error setup forcefield \n";
std::cout << ff->Energy() << " " << ff->GetUnit() << std::endl;
}
int main() {
// prepare ernergy calculation, conversion to inchi-key and to smiles
conv.SetOutStream(&std::cout);
conv.SetInAndOutFormats("SMI","smi");
convToInchi.SetOutStream(&std::cout);
convToInchi.SetOutFormat("inchi");
convToInchi.AddOption("K",OpenBabel::OBConversion::OUTOPTIONS);
ff = OpenBabel::OBForceField::FindType("UFF");
if (!ff) std::cout << "error finding forcefield \n";
// create a di-Silicium molecule "by Hand"
obMol.NewAtom(1);
obMol.NewAtom(2);
obMol.AddBond(1,2,1);
obMol.GetAtomById(1)->SetAtomicNum(14);
obMol.GetAtomById(2)->SetAtomicNum(14);
calcEnergy_and_output();
// increase bond type (2x) and output each
obMol.GetBond(1,2)->SetBondOrder(2);
calcEnergy_and_output();
obMol.GetBond(1,2)->SetBondOrder(3);
calcEnergy_and_output();
// read same molecules from smiles
conv.ReadString(&obMol,"[Si][Si]");
calcEnergy_and_output();
conv.ReadString(&obMol,"[Si]=[Si]");
calcEnergy_and_output();
conv.ReadString(&obMol,"[Si]#[Si]");
calcEnergy_and_output();
}
output:
%./test
[SiH3][SiH3] PZPGRFITIJYNEJ-UHFFFAOYSA-N 0.131888 kJ/mol
[SiH2]=[SiH2] IKQCSDAPCHQIFR-UHFFFAOYSA-N 0.0511706 kJ/mol
[SiH3]#[SiH3] NWEREQWWNYLPTF-UHFFFAOYSA-N 0.639372 kJ/mol
[Si][Si]
==============================
*** Open Babel Warning in InChI code
#0 :Accepted unusual valence(s): Si(1)
NTQGILPNLZZOJH-UHFFFAOYSA-N0.131888 kJ/mol
[Si]=[Si]
NTQGILPNLZZOJH-UHFFFAOYSA-N
0.0511706 kJ/mol
[Si]#[Si]
==============================
*** Open Babel Warning in InChI code
#0 :Accepted unusual valence(s): Si(3)
NTQGILPNLZZOJH-UHFFFAOYSA-N
0.639372 kJ/mol
------------------------------------------------------------------------------
LIMITED TIME SALE - Full Year of Microsoft Training For Just $49.99!
1,500+ hours of tutorials including VisualStudio 2012, Windows 8, SharePoint
2013, SQL 2012, MVC 4, more. BEST VALUE: New Multi-Library Power Pack includes
Mobile, Cloud, Java, and UX Design. Lowest price ever! Ends 9/20/13.
http://pubads.g.doubleclick.net/gampad/clk?id=58041151&iu=/4140/ostg.clktrk
_______________________________________________
OpenBabel-discuss mailing list
OpenBabel-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/openbabel-discuss
