Hi
I struggle for getting OBFF parameters again.
In the following example code, a OBForceField instance (pFF) was
down-casted on a new-defined class (ffpGaff)
making an instance (pP).
This works well in some cases, but an exception is thrown in other cases.
Any idea instead of this would be greatly appreciated.
Thanks in advance.
// *****************************************************
// an example code
#include <openbabel/base.h>
#include <openbabel/mol.h>
#include <openbabel/obconversion.h>
#include <openbabel/forcefield.h>
#include <openbabel/forcefields/forcefieldgaff.h>
using namespace std;
using namespace OpenBabel;
// bond parameter
struct bondKind
{
char aT1[3]; // atom Type 1,
char aT2[3]; // atom Type 2,
double Kb; // force constant (kcal/mol).
double req; // equibrium distance (A).
};
// angle parameter
struct angleKind
{
char aT1[3]; // atom Type 1
char aT2[3]; // atom Type 2
char aT3[3]; // atom Type 3
double Kth; // constant of bending motion (kcal/(mol radian^2))
double theq; // equibrium angle (degree)
};
// torsion parameter
struct torKind
{
char aT1[3];
char aT2[3];
char aT3[3];
char aT4[3];
int m; // number of bond paths. ignoired.
double vn2;
double gamma; // phase
double n; // the periodicity of torsion (multiplicity)
};
// out-of-plane motion parameter
struct oopKind
{
char aT1[3];
char aT2[3];
char aT3[3];
char aT4[3];
double vn2; // equals to Kchi for charmm pot (kcal/mol)
double gamma; // equals to delta for charmm pot (degree)
double n; // the periodicity of torsion
};
// non-bonding (vdW) parameter
struct nbKind
{
char aT[3];
double R;
double eps;
};
// a new-defined class
class ffpGaff : public OBForceFieldGaff
{
public:
vector<bondKind> bkmol;
vector<angleKind> akmol;
vector<torKind> tkmol;
vector<oopKind> okmol;
vector<nbKind> nkmol;
void getBkmol( void )
{
if( &( _bondcalculations ) == NULL) return;
bkmol.clear();
vector<OBFFBondCalculationGaff>::iterator i;
for (i = _bondcalculations.begin(); i != _bondcalculations.end(); ++i)
{
bondKind tempBk;
strcpy( tempBk.aT1, (*i).a->GetType() );
strcpy( tempBk.aT2, (*i).b->GetType() );
tempBk.Kb = (*i).kr / KCAL_TO_KJ ;
tempBk.req = (*i).r0;
bkmol.push_back( tempBk );
}
}
void getAkmol( void )
{
if( &( _anglecalculations ) == NULL ) return;
akmol.clear();
vector<OBFFAngleCalculationGaff>::iterator i;
for (i = _anglecalculations.begin(); i != _anglecalculations.end(); ++i)
{
angleKind tempAk;
strcpy( tempAk.aT1, (*i).a->GetType() );
strcpy( tempAk.aT2, (*i).b->GetType() );
strcpy( tempAk.aT3, (*i).c->GetType() );
tempAk.Kth = (*i).kth / KCAL_TO_KJ;
tempAk.theq = (*i).theta0;
akmol.push_back( tempAk );
}
}
void getTkmol( void )
{
if( &( _torsioncalculations ) == NULL ) return;
tkmol.clear();
vector<OBFFTorsionCalculationGaff>::iterator i;
for (i = _torsioncalculations.begin(); i != _torsioncalculations.end();
++i)
{
torKind tempTk;
strcpy( tempTk.aT1, (*i).a->GetType() );
strcpy( tempTk.aT2, (*i).b->GetType() );
strcpy( tempTk.aT3, (*i).c->GetType() );
strcpy( tempTk.aT4, (*i).d->GetType() );
tempTk.m = 0;
tempTk.vn2 = (*i).vn_half / KCAL_TO_KJ;
tempTk.gamma = (*i).gamma;
tempTk.n = (*i).n;
tkmol.push_back( tempTk );
}
}
void getOkmol( void )
{
if( &( _oopcalculations ) == NULL ) return;
okmol.clear();
vector<OBFFOOPCalculationGaff>::iterator i;
for (i = _oopcalculations.begin(); i != _oopcalculations.end(); ++i)
{
oopKind tempOk;
strcpy( tempOk.aT1, (*i).a->GetType() );
strcpy( tempOk.aT2, (*i).b->GetType() );
strcpy( tempOk.aT3, (*i).c->GetType() );
strcpy( tempOk.aT4, (*i).d->GetType() );
tempOk.vn2 = (*i).vn_half / KCAL_TO_KJ;
tempOk.gamma = (*i).gamma;
tempOk.n = (*i).n;
okmol.push_back( tempOk );
}
}
void getNkmol( void )
{
if( &( _ffvdwparams ) == NULL ) return;
nkmol.clear();
for (int i = 0; i < _ffvdwparams.size() ; i++ )
{
nbKind tempNk;
strcpy( tempNk.aT, _ffvdwparams[i]._a.c_str() );
FOR_ATOMS_OF_MOL( a, _mol )
{
char* aType = a->GetType();
if( aType != NULL && strcmp( aType, tempNk.aT ) == 0)
{
tempNk.R = _ffvdwparams[i]._dpar[0];
tempNk.eps = _ffvdwparams[i]._dpar[1];
nkmol.push_back( tempNk );
}
}
}
}
};
// Execute as >program_name molecular_file_name
int main( int argc, char* argv[] )
{
char *program_name= argv[0];
string basename, filename = "", ff = "GAFF";
basename = filename = argv[1];
size_t extPos = filename.rfind('.');
if (extPos != string::npos)
{ basename = filename.substr(0, extPos); }
// Find Input filetype
OBConversion conv;
OBFormat *format_in = conv.FormatFromExt(filename.c_str());
if (!format_in || !conv.SetInFormat(format_in))
{
cout << program_name << ": cannot read input/output format!" << endl;
return -1;
}
ifstream ifs;
// Read the file
ifs.open(filename.c_str());
if (!ifs)
{
cout << program_name << ": cannot read input file!" << endl;
return -1;
}
OBMol mol;
mol.Clear();
if (!conv.Read(&mol, &ifs))
return -1;
if (mol.Empty())
return -1;
OBForceField* pFF = OBForceField::FindForceField(ff);
if (!pFF)
{
cout << program_name << ": could not find forcefield '" << ff << "'." <<
endl;
return -1;
}
mol.AddHydrogens(false, true); // hydrogens must be added before Setup(mol)
is called
if (!pFF->Setup(mol))
{
cout << program_name << ": could not setup force field." << endl;
return -1;
}
ffpGaff* pP = static_cast<ffpGaff*>(pFF); // down-casting
// ffpGaff* pP = dynamic_cast<ffpGaff*>(pFF); // pP is always NULL for
dynamic_cast
if( pP == NULL ) return -1;
// bond stretching
pP->getBkmol();
if( pP->bkmol.empty() ) return -1;
cout << "GAFF bond stretching parameters" << endl;
for( int i = 0 ; i < pP->bkmol.size(); i++ )
{
cout << pP->bkmol[i].aT1 << " "
<< pP->bkmol[i].aT2 << " "
<< pP->bkmol[i].Kb << " "
<< pP->bkmol[i].req << endl;
}
// angle bending
pP->getAkmol();
if( pP->akmol.empty() ) return -1;
cout << "GAFF angle bending parameters" << endl;
for( int i = 0 ; i < pP->akmol.size(); i++ )
{
cout << pP->akmol[i].aT1 << " "
<< pP->akmol[i].aT2 << " "
<< pP->akmol[i].aT3 << " "
<< pP->akmol[i].Kth << " "
<< pP->akmol[i].theq << endl;
}
// torsional motion
pP->getTkmol();
if( pP->tkmol.empty() ) return -1;
cout << "GAFF torsional parameters" << endl;
for( int i = 0 ; i < pP->tkmol.size(); i++ )
{
cout << pP->tkmol[i].aT1 << " "
<< pP->tkmol[i].aT2 << " "
<< pP->tkmol[i].aT3 << " "
<< pP->tkmol[i].aT4 << " "
<< pP->tkmol[i].vn2 << " "
<< pP->tkmol[i].gamma << " "
<< pP->tkmol[i].n << endl;
}
// out-of-plane motion
pP->getOkmol();
if( pP->okmol.empty() ) return -1;
cout << "GAFF out-of-plane motion parameters" << endl;
for( int i = 0 ; i < pP->okmol.size(); i++ )
{
cout << pP->okmol[i].aT1 << " "
<< pP->okmol[i].aT2 << " "
<< pP->okmol[i].aT3 << " "
<< pP->okmol[i].aT4 << " "
<< pP->okmol[i].vn2 << " "
<< pP->okmol[i].gamma << " "
<< pP->okmol[i].n << endl;
}
// non-bonding interactions
pP->getNkmol();
if( pP->nkmol.empty() ) return -1;
cout << "GAFF non-bonding interaction parameters" << endl;
for( int i = 0 ; i < pP->nkmol.size(); i++ )
{
cout << pP->nkmol[i].aT << " "
<< pP->nkmol[i].R << " "
<< pP->nkmol[i].eps << endl;
}
return 0;
}
// end of an example code
// *****************************************************
An execution result (in a well-worked case) on MSVC for toluene is:
// *****************************************************
GAFF bond stretching parameters
c3 ca 323.5 1.513
ca ca 478.4 1.387
ca ca 478.4 1.387
ca ca 478.4 1.387
ca ca 478.4 1.387
ca ca 478.4 1.387
ca ca 478.4 1.387
c3 hc 337.3 1.092
c3 hc 337.3 1.092
c3 hc 337.3 1.092
ca ha 344.3 1.087
ca ha 344.3 1.087
ca ha 344.3 1.087
ca ha 344.3 1.087
ca ha 344.3 1.087
GAFF angle bending parameters
ca c3 hc 46.96 110.15
ca c3 hc 46.96 110.15
ca c3 hc 46.96 110.15
hc c3 hc 39.43 108.35
hc c3 hc 39.43 108.35
hc c3 hc 39.43 108.35
c3 ca ca 63.84 120.63
c3 ca ca 63.84 120.63
ca ca ca 67.18 119.97
ca ca ca 67.18 119.97
ca ca ha 48.46 120.01
ca ca ha 48.46 120.01
ca ca ca 67.18 119.97
ca ca ha 48.46 120.01
ca ca ha 48.46 120.01
ca ca ca 67.18 119.97
ca ca ha 48.46 120.01
ca ca ha 48.46 120.01
ca ca ca 67.18 119.97
ca ca ha 48.46 120.01
ca ca ha 48.46 120.01
ca ca ca 67.18 119.97
ca ca ha 48.46 120.01
ca ca ha 48.46 120.01
GAFF torsional parameters
hc c3 ca ca 0 0 2
hc c3 ca ca 0 0 2
hc c3 ca ca 0 0 2
hc c3 ca ca 0 0 2
hc c3 ca ca 0 0 2
hc c3 ca ca 0 0 2
c3 ca ca ca 3.625 180 2
c3 ca ca ha 3.625 180 2
ca ca ca ca 3.625 180 2
ca ca ca ha 3.625 180 2
ca ca ca ca 3.625 180 2
ca ca ca ha 3.625 180 2
ha ca ca ca 3.625 180 2
ha ca ca ha 3.625 180 2
ca ca ca ca 3.625 180 2
ca ca ca ha 3.625 180 2
ha ca ca ca 3.625 180 2
ha ca ca ha 3.625 180 2
ca ca ca ca 3.625 180 2
ca ca ca ha 3.625 180 2
ha ca ca ca 3.625 180 2
ha ca ca ha 3.625 180 2
ca ca ca ca 3.625 180 2
ca ca ca ha 3.625 180 2
ha ca ca ca 3.625 180 2
ha ca ca ha 3.625 180 2
c3 ca ca ca 3.625 180 2
c3 ca ca ha 3.625 180 2
ca ca ca ca 3.625 180 2
ca ca ca ha 3.625 180 2
GAFF out-of-plane motion parameters
ca ca ca c3 1.1 180 2
ca ca ca ha 1.1 180 2
ca ca ca ha 1.1 180 2
ca ca ca ha 1.1 180 2
ca ca ca ha 1.1 180 2
ca ca ca ha 1.1 180 2
GAFF non-bonding interaction parameters
ha 1.459 0.015
ha 1.459 0.015
ha 1.459 0.015
ha 1.459 0.015
ha 1.459 0.015
hc 1.487 0.0157
hc 1.487 0.0157
hc 1.487 0.0157
c3 1.908 0.1094
ca 1.908 0.086
ca 1.908 0.086
ca 1.908 0.086
ca 1.908 0.086
ca 1.908 0.086
ca 1.908 0.086
// *****************************************************
Used molecular file (toluene.cml) is:
<molecule>
<atomArray>
<atom id="a1" elementType="C" x3="0.932545" y3="-0.068050" z3="0.029397"/>
<atom id="a2" elementType="C" x3="2.429382" y3="-0.034585"
z3="-0.058760"/>
<atom id="a3" elementType="C" x3="3.122390" y3="1.181313" z3="-0.094235"/>
<atom id="a4" elementType="C" x3="4.515950" y3="1.197955" z3="-0.177582"/>
<atom id="a5" elementType="C" x3="5.227527" y3="0.000593" z3="-0.226551"/>
<atom id="a6" elementType="C" x3="4.546450" y3="-1.214677"
z3="-0.192359"/>
<atom id="a7" elementType="C" x3="3.153262" y3="-1.232774"
z3="-0.108991"/>
<atom id="a8" elementType="H" x3="0.509709" y3="0.941061" z3="0.072209"/>
<atom id="a9" elementType="H" x3="0.619747" y3="-0.601149" z3="0.932968"/>
<atom id="a10" elementType="H" x3="0.514468" y3="-0.570641"
z3="-0.848686"/>
<atom id="a11" elementType="H" x3="2.582472" y3="2.124526"
z3="-0.057180"/>
<atom id="a12" elementType="H" x3="5.046193" y3="2.146205"
z3="-0.204403"/>
<atom id="a13" elementType="H" x3="6.312108" y3="0.014698"
z3="-0.291477"/>
<atom id="a14" elementType="H" x3="5.099906" y3="-2.149094"
z3="-0.230715"/>
<atom id="a15" elementType="H" x3="2.633402" y3="-2.187764"
z3="-0.083315"/>
</atomArray>
<bondArray>
<bond atomRefs2="a1 a2" order="1"/>
<bond atomRefs2="a2 a3" order="2"/>
<bond atomRefs2="a3 a4" order="1"/>
<bond atomRefs2="a4 a5" order="2"/>
<bond atomRefs2="a5 a6" order="1"/>
<bond atomRefs2="a6 a7" order="2"/>
<bond atomRefs2="a2 a7" order="1"/>
<bond atomRefs2="a1 a8" order="1"/>
<bond atomRefs2="a1 a9" order="1"/>
<bond atomRefs2="a1 a10" order="1"/>
<bond atomRefs2="a3 a11" order="1"/>
<bond atomRefs2="a4 a12" order="1"/>
<bond atomRefs2="a5 a13" order="1"/>
<bond atomRefs2="a6 a14" order="1"/>
<bond atomRefs2="a7 a15" order="1"/>
</bondArray>
</molecule>
------------------------------------------------------------------------------
Site24x7 APM Insight: Get Deep Visibility into Application Performance
APM + Mobile APM + RUM: Monitor 3 App instances at just $35/Month
Monitor end-to-end web transactions and take corrective actions now
Troubleshoot faster and improve end-user experience. Signup Now!
http://pubads.g.doubleclick.net/gampad/clk?id=267308311&iu=/4140
_______________________________________________
OpenBabel-Devel mailing list
OpenBabel-Devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/openbabel-devel
