Skip to content

Commit

Permalink
Merge branch 'master' of git://github.com/ptosco/rdkit into ptosco-ma…
Browse files Browse the repository at this point in the history 
…ster
  • Loading branch information
@greglandrum
greglandrum committed Sep 26, 2013
2 parents 0b83ef0 + 3646270 commit dad1e1d
Show file tree
Hide file tree
Showing 4 changed files with 190 additions and 85 deletions.
188 changes: 104 additions & 84 deletions Code/ForceField/UFF/AngleBend.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,32 +21,32 @@ namespace ForceFields {

namespace Utils {
double calcAngleForceConstant(double theta0,
double bondOrder12,double bondOrder23,
const AtomicParams *at1Params,
const AtomicParams *at2Params,
const AtomicParams *at3Params){
double cosTheta0=cos(theta0);
double r12 = calcBondRestLength(bondOrder12,at1Params,at2Params);
double r23 = calcBondRestLength(bondOrder23,at2Params,at3Params);
double r13 = sqrt(r12*r12 + r23*r23 - 2.*r12*r23*cosTheta0);
double beta = 2.*Params::G/(r12*r23);

double preFactor = beta*at1Params->Z1*at3Params->Z1 / int_pow<5>(r13);
double rTerm = r12*r23;
double innerBit = 3.*rTerm*(1.-cosTheta0*cosTheta0) - r13*r13*cosTheta0;
double bondOrder12,double bondOrder23,
const AtomicParams *at1Params,
const AtomicParams *at2Params,
const AtomicParams *at3Params){
double cosTheta0=cos(theta0);
double r12 = calcBondRestLength(bondOrder12,at1Params,at2Params);
double r23 = calcBondRestLength(bondOrder23,at2Params,at3Params);
double r13 = sqrt(r12*r12 + r23*r23 - 2.*r12*r23*cosTheta0);
double beta = 2.*Params::G/(r12*r23);

double preFactor = beta*at1Params->Z1*at3Params->Z1 / int_pow<5>(r13);
double rTerm = r12*r23;
double innerBit = 3.*rTerm*(1.-cosTheta0*cosTheta0) - r13*r13*cosTheta0;

double res=preFactor*rTerm*innerBit;
return res;
double res=preFactor*rTerm*innerBit;
return res;
}
} // end of namespace Utils

AngleBendContrib::AngleBendContrib(ForceField *owner,
unsigned int idx1,unsigned int idx2,unsigned int idx3,
double bondOrder12,double bondOrder23,
const AtomicParams *at1Params,
const AtomicParams *at2Params,
const AtomicParams *at3Params,
unsigned int order){
unsigned int idx1,unsigned int idx2,unsigned int idx3,
double bondOrder12,double bondOrder23,
const AtomicParams *at1Params,
const AtomicParams *at2Params,
const AtomicParams *at3Params,
unsigned int order){
PRECONDITION(owner,"bad owner");
PRECONDITION(at1Params,"bad params pointer");
PRECONDITION(at2Params,"bad params pointer");
Expand All @@ -55,20 +55,40 @@ namespace ForceFields {
RANGE_CHECK(0,idx1,owner->positions().size()-1);
RANGE_CHECK(0,idx2,owner->positions().size()-1);
RANGE_CHECK(0,idx3,owner->positions().size()-1);
// the following is a hack to get decent geometries
// with 3- and 4-membered rings incorporating sp2 atoms
double theta0 = at2Params->theta0;
if (order >= 30) {
switch (order) {
case 30:
theta0 = 150.0 / 180.0 * M_PI;
break;
case 35:
theta0 = 60.0 / 180.0 * M_PI;
break;
case 40:
theta0 = 135.0 / 180.0 * M_PI;
break;
case 45:
theta0 = 90.0 / 180.0 * M_PI;
break;
}
order = 0;
}
dp_forceField = owner;
d_at1Idx = idx1;
d_at2Idx = idx2;
d_at3Idx = idx3;
d_order = order;
this->d_forceConstant = Utils::calcAngleForceConstant(at2Params->theta0,
bondOrder12,bondOrder23,
at1Params,at2Params,at3Params);
this->d_forceConstant = Utils::calcAngleForceConstant
(theta0, bondOrder12,bondOrder23,
at1Params,at2Params,at3Params);
if(order==0){
double sinTheta0=sin(at2Params->theta0);
double cosTheta0=cos(at2Params->theta0);
this->d_C2 = 1./(4.*std::max(sinTheta0*sinTheta0,1e-8));
this->d_C1 = -4.*this->d_C2*cosTheta0;
this->d_C0 = this->d_C2*(2.*cosTheta0*cosTheta0 + 1.);
double sinTheta0=sin(theta0);
double cosTheta0=cos(theta0);
this->d_C2 = 1./(4.*std::max(sinTheta0*sinTheta0,1e-8));
this->d_C1 = -4.*this->d_C2*cosTheta0;
this->d_C0 = this->d_C2*(2.*cosTheta0*cosTheta0 + 1.);
}
}

Expand All @@ -80,14 +100,14 @@ namespace ForceFields {
double dist2=this->dp_forceField->distance(this->d_at2Idx,this->d_at3Idx,pos);

RDGeom::Point3D p1(pos[3*this->d_at1Idx],
pos[3*this->d_at1Idx+1],
pos[3*this->d_at1Idx+2]);
pos[3*this->d_at1Idx+1],
pos[3*this->d_at1Idx+2]);
RDGeom::Point3D p2(pos[3*this->d_at2Idx],
pos[3*this->d_at2Idx+1],
pos[3*this->d_at2Idx+2]);
pos[3*this->d_at2Idx+1],
pos[3*this->d_at2Idx+2]);
RDGeom::Point3D p3(pos[3*this->d_at3Idx],
pos[3*this->d_at3Idx+1],
pos[3*this->d_at3Idx+2]);
pos[3*this->d_at3Idx+1],
pos[3*this->d_at3Idx+2]);
RDGeom::Point3D p12=p1-p2;
RDGeom::Point3D p32=p3-p2;
double cosTheta = p12.dotProduct(p32)/(dist1*dist2);
Expand All @@ -111,14 +131,14 @@ namespace ForceFields {
//std::cout << "\tAngle("<<this->d_at1Idx<<","<<this->d_at2Idx<<","<<this->d_at3Idx<<") " << dist1 << " " << dist2 << std::endl;

RDGeom::Point3D p1(pos[3*this->d_at1Idx],
pos[3*this->d_at1Idx+1],
pos[3*this->d_at1Idx+2]);
pos[3*this->d_at1Idx+1],
pos[3*this->d_at1Idx+2]);
RDGeom::Point3D p2(pos[3*this->d_at2Idx],
pos[3*this->d_at2Idx+1],
pos[3*this->d_at2Idx+2]);
pos[3*this->d_at2Idx+1],
pos[3*this->d_at2Idx+2]);
RDGeom::Point3D p3(pos[3*this->d_at3Idx],
pos[3*this->d_at3Idx+1],
pos[3*this->d_at3Idx+2]);
pos[3*this->d_at3Idx+1],
pos[3*this->d_at3Idx+2]);
double *g1=&(grad[3*this->d_at1Idx]);
double *g2=&(grad[3*this->d_at2Idx]);
double *g3=&(grad[3*this->d_at3Idx]);
Expand Down Expand Up @@ -170,26 +190,26 @@ namespace ForceFields {

double res=0.0;
if(this->d_order==0){
res=this->d_C0 + this->d_C1*cosTheta + this->d_C2*cos2Theta;
res=this->d_C0 + this->d_C1*cosTheta + this->d_C2*cos2Theta;
} else {
switch(this->d_order){
case 1:
res=-cosTheta;
break;
case 2:
res=cos2Theta;
break;
case 3:
// cos(3x) = cos^3(x) - 3*cos(x)*sin^2(x)
res = cosTheta*(cosTheta*cosTheta-3.*sinThetaSq);
break;
case 4:
// cos(4x) = cos^4(x) - 6*cos^2(x)*sin^2(x)+sin^4(x)
res = int_pow<4>(cosTheta) - 6.*cosTheta*cosTheta*sinThetaSq + sinThetaSq*sinThetaSq;
break;
}
res = 1.-res;
res /= (double)(this->d_order*this->d_order);
switch(this->d_order){
case 1:
res=-cosTheta;
break;
case 2:
res=cos2Theta;
break;
case 3:
// cos(3x) = cos^3(x) - 3*cos(x)*sin^2(x)
res = cosTheta*(cosTheta*cosTheta-3.*sinThetaSq);
break;
case 4:
// cos(4x) = cos^4(x) - 6*cos^2(x)*sin^2(x)+sin^4(x)
res = int_pow<4>(cosTheta) - 6.*cosTheta*cosTheta*sinThetaSq + sinThetaSq*sinThetaSq;
break;
}
res = 1.-res;
res /= (double)(this->d_order*this->d_order);
}
return res;
}
Expand All @@ -203,33 +223,33 @@ namespace ForceFields {
double sin2Theta = 2.*sinTheta*cosTheta;

if(this->d_order==0){
dE_dTheta = -1.*this->d_forceConstant*(this->d_C1*sinTheta +
2.*this->d_C2*sin2Theta);
dE_dTheta = -1.*this->d_forceConstant*(this->d_C1*sinTheta +
2.*this->d_C2*sin2Theta);
} else {
// E = k/n^2 [1-cos(n theta)]
// dE = - k/n^2 * d cos(n theta)
// E = k/n^2 [1-cos(n theta)]
// dE = - k/n^2 * d cos(n theta)

// these all use:
// d cos(ax) = -a sin(ax)
// these all use:
// d cos(ax) = -a sin(ax)

switch(this->d_order){
case 1:
dE_dTheta = -sinTheta;
break;
case 2:
// sin(2*x) = 2*cos(x)*sin(x)
dE_dTheta = sin2Theta;
break;
case 3:
// sin(3*x) = 3*sin(x) - 4*sin^3(x)
dE_dTheta = sinTheta*(3.-4.*sinTheta*sinTheta);
break;
case 4:
// sin(4*x) = cos(x)*(4*sin(x) - 8*sin^3(x))
dE_dTheta = cosTheta*sinTheta*(4.-8.*sinTheta*sinTheta);
break;
}
dE_dTheta *= this->d_forceConstant/(double)(this->d_order);
switch(this->d_order){
case 1:
dE_dTheta = -sinTheta;
break;
case 2:
// sin(2*x) = 2*cos(x)*sin(x)
dE_dTheta = sin2Theta;
break;
case 3:
// sin(3*x) = 3*sin(x) - 4*sin^3(x)
dE_dTheta = sinTheta*(3.-4.*sinTheta*sinTheta);
break;
case 4:
// sin(4*x) = cos(x)*(4*sin(x) - 8*sin^3(x))
dE_dTheta = cosTheta*sinTheta*(4.-8.*sinTheta*sinTheta);
break;
}
dE_dTheta *= this->d_forceConstant/(double)(this->d_order);
}
return dE_dTheta;
}
Expand Down
31 changes: 31 additions & 0 deletions Code/GraphMol/ForceFieldHelpers/UFF/Builder.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,7 @@ namespace RDKit {
ROMol::ADJ_ITER end1Nbrs;
ROMol::ADJ_ITER nbr2Idx;
ROMol::ADJ_ITER end2Nbrs;
RingInfo *rings=mol.getRingInfo();

unsigned int nAtoms=mol.getNumAtoms();
for(unsigned int j=0;j<nAtoms;j++){
Expand Down Expand Up @@ -180,8 +181,38 @@ namespace RDKit {
case Atom::SP:
order=1;
break;
// the following is a hack to get decent geometries
// with 3- and 4-membered rings incorporating sp2 atoms
case Atom::SP2:
order=3;
// if the central atom is in a ring of size 3
if (rings->isAtomInRingOfSize(j, 3)) {
// if the central atom and one of the bonded atoms, but not the
// other one are inside a ring, then this angle is between a
// ring substituent and a ring edge
if ((rings->isAtomInRingOfSize(i, 3) && !rings->isAtomInRingOfSize(k, 3))
|| (!rings->isAtomInRingOfSize(i, 3) && rings->isAtomInRingOfSize(k, 3))) {
order = 30;
}
// if all atoms are inside the ring, then this is one of ring angles
else if (rings->isAtomInRingOfSize(i, 3) && rings->isAtomInRingOfSize(k, 3)) {
order = 35;
}
}
// if the central atom is in a ring of size 4
else if (rings->isAtomInRingOfSize(j, 4)) {
// if the central atom and one of the bonded atoms, but not the
// other one are inside a ring, then this angle is between a
// ring substituent and a ring edge
if ((rings->isAtomInRingOfSize(i, 4) && !rings->isAtomInRingOfSize(k, 4))
|| (!rings->isAtomInRingOfSize(i, 4) && rings->isAtomInRingOfSize(k, 4))) {
order = 40;
}
// if all atoms are inside the ring, then this is one of ring angles
else if (rings->isAtomInRingOfSize(i, 4) && rings->isAtomInRingOfSize(k, 4)) {
order = 45;
}
}
break;
case Atom::SP3D2:
order=4;
Expand Down
39 changes: 38 additions & 1 deletion Code/GraphMol/ForceFieldHelpers/UFF/testUFFHelpers.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include <GraphMol/SmilesParse/SmilesWrite.h>
#include <GraphMol/FileParsers/FileParsers.h>
#include <GraphMol/FileParsers/MolSupplier.h>
#include <GraphMol/FileParsers/MolWriters.h>

#include <GraphMol/ForceFieldHelpers/UFF/AtomTyper.h>
#include <GraphMol/ForceFieldHelpers/UFF/Builder.h>
Expand Down Expand Up @@ -846,6 +847,42 @@ void testSFIssue3009337(){
}
BOOST_LOG(rdErrorLog) << " done" << std::endl;
}
void testGitHubIssue62() {
BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
BOOST_LOG(rdErrorLog) << " Testing GitHubIssue62." << std::endl;

std::string pathName=getenv("RDBASE");
pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
{
double energyValues[] =
{ 38.687, 174.698, 337.986, 115.248,
2.482, 1.918, 10.165, 98.469, 39.078,
267.236, 15.747, 202.121, 205.539,
20.044, 218.986, 79.627 };
SmilesMolSupplier smiSupplier(pathName + "/Issue62.smi");
SDWriter *sdfWriter = new SDWriter(pathName + "/Issue62.sdf");
for (unsigned int i = 0; i < smiSupplier.length(); ++i) {
ROMol *mol = MolOps::addHs(*(smiSupplier[i]));
TEST_ASSERT(mol);
std::string molName = "";
if (mol->hasProp("_Name")) {
mol->getProp("_Name", molName);
}
DGeomHelpers::EmbedMolecule(*mol);
ForceFields::ForceField *field = UFF::constructForceField(*mol);
TEST_ASSERT(field);
field->initialize();
int needMore = field->minimize(200, 1.e-6, 1.e-3);
TEST_ASSERT(!needMore);
sdfWriter->write(*mol);
double e = field->calcEnergy();
BOOST_LOG(rdErrorLog) << molName << " " << e << std::endl;
TEST_ASSERT(fabs(e - energyValues[i]) < 1.);
}
sdfWriter->close();
BOOST_LOG(rdErrorLog) << " done" << std::endl;
}
}

//-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//
Expand All @@ -866,5 +903,5 @@ int main(){
#endif
testMissingParams();
testSFIssue3009337();

testGitHubIssue62();
}
17 changes: 17 additions & 0 deletions Code/GraphMol/ForceFieldHelpers/UFF/test_data/Issue62.smi
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
SMILES Name
C1C(=O)NC1 beta-lactam
O=C1CC1 cyclopropanone
N1CC1 aziridine
O=C1NC1 aziridinone
C(C)(C)=C isobutene
O=C(C)C acetone
O=C1CCCC1 cyclopentanone
C1C=C1 cyclopropene
C1CC=C1 cyclobutene
C1CC1 cyclopropane
O=C1NCCC1 gamma-lactam
O=C1[C@@H]([C@H]1c1ccccc1)c1ccccc1 trans-2,3-diphenylcyclopropan-1-one
O=C1[C@H]([C@H]1c1ccccc1)c1ccccc1 cis-2,3-diphenylcyclopropan-1-one
O=C1C(=C1c1ccccc1)c1ccccc1 2,3-diphenylcyclopropa-2-en-1-one
C1CC1=C1CCC1 cyclopropylidencyclobutane
C1CCC1=C(c1ccccc1)c1ccccc1 diphenylmethylidencyclobutane

0 comments on commit dad1e1d

Please sign in to comment.