pub_soft/rdkit
3
0
Fork 0
mirror of https://github.com/rdkit/rdkit.git synced 2026-06-06 22:39:55 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
AlternateValence_21Nov2012
rdkit/Code/GraphMol/ForceFieldHelpers/UFF/Builder.cpp
Greg Landrum f3fbef45c5 update copyright statements
2010-09-26 17:04:37 +00:00

555 lines
23 KiB
C++
Raw Permalink Blame History

// $Id$
//
// Copyright (C) 2004-2010 Greg Landrum and Rational Discovery LLC
//
// @@ 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.
//
#include <iostream>
#include <cmath>
#include <RDGeneral/Invariant.h>
#include <GraphMol/RDKitBase.h>
#include <GraphMol/SmilesParse/SmilesParse.h>
#include <GraphMol/Substruct/SubstructMatch.h>
#include <ForceField/ForceField.h>
#include <ForceField/UFF/Params.h>
#include <ForceField/UFF/Contribs.h>
#include "AtomTyper.h"
#include "Builder.h"
namespace RDKit {
namespace UFF {
using namespace ForceFields::UFF;
namespace Tools {
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
void addBonds(const ROMol &mol,const AtomicParamVect &params,
ForceFields::ForceField *field){
PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
PRECONDITION(field,"bad forcefield");
for (ROMol::ConstBondIterator bi=mol.beginBonds();
bi != mol.endBonds();
bi++) {
int idx1=(*bi)->getBeginAtomIdx();
int idx2=(*bi)->getEndAtomIdx();
// FIX: recognize amide bonds here.
if(params[idx1]&&params[idx2]){
BondStretchContrib *contrib;
contrib = new BondStretchContrib(field,idx1,idx2,
(*bi)->getBondTypeAsDouble(),
params[idx1],params[idx2]);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
}
}
// ------------------------------------------------------------------------
//
// the matrix returned by this contains:
// -1: if atoms i and j are directly connected
// idx: if atoms i and j are connected via atom idx
// -2: otherwise
//
// NOTE: the caller is responsible for calling delete []
// on the result
//
// ------------------------------------------------------------------------
boost::shared_array<int> buildNeighborMatrix(const ROMol &mol){
unsigned int nAtoms = mol.getNumAtoms();
boost::shared_array<int> res(new int[nAtoms*nAtoms]);
for(unsigned int i=0;i<nAtoms;i++){
unsigned int iTab=i*nAtoms;
for(unsigned int j=i;j<nAtoms;j++){
res[iTab+j] = -2;
res[i+j*nAtoms] = -2;
}
}
for(unsigned int i=0;i<mol.getNumBonds();i++){
const Bond *bondi=mol.getBondWithIdx(i);
res[bondi->getBeginAtomIdx()*nAtoms+bondi->getEndAtomIdx()] = -1;
res[bondi->getEndAtomIdx()*nAtoms+bondi->getBeginAtomIdx()] = -1;
for(unsigned int j=i+1;j<mol.getNumBonds();j++){
const Bond *bondj=mol.getBondWithIdx(j);
int idx1=-1,idx2=-1,idx3=-1;
if(bondi->getBeginAtomIdx()==bondj->getBeginAtomIdx()){
idx1 = bondi->getEndAtomIdx();
idx2 = bondi->getBeginAtomIdx();
idx3 = bondj->getEndAtomIdx();
} else if(bondi->getBeginAtomIdx()==bondj->getEndAtomIdx()){
idx1 = bondi->getEndAtomIdx();
idx2 = bondi->getBeginAtomIdx();
idx3 = bondj->getBeginAtomIdx();
} else if(bondi->getEndAtomIdx()==bondj->getBeginAtomIdx()){
idx1 = bondi->getBeginAtomIdx();
idx2 = bondi->getEndAtomIdx();
idx3 = bondj->getEndAtomIdx();
} else if(bondi->getEndAtomIdx()==bondj->getEndAtomIdx()){
idx1 = bondi->getBeginAtomIdx();
idx2 = bondi->getEndAtomIdx();
idx3 = bondj->getBeginAtomIdx();
}
if(idx1>-1){
res[idx1*nAtoms+idx3] = idx2;
res[idx3*nAtoms+idx1] = idx2;
}
}
}
return res;
}
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
void addAngles(const ROMol &mol,const AtomicParamVect &params,
ForceFields::ForceField *field,boost::shared_array<int> neighborMatrix){
PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
PRECONDITION(field,"bad forcefield");
unsigned int nAtoms=mol.getNumAtoms();
for(unsigned int i=0;i<nAtoms;i++){
if(!params[i]) continue;
for(unsigned int j=i+1;j<nAtoms;j++){
if(!params[j]) continue;
if(neighborMatrix[i*nAtoms+j]>-1){
int k = neighborMatrix[i*nAtoms+j];
if(!params[k]) continue;
const Atom *atomK = mol.getAtomWithIdx(k);
// skip special cases:
if( !(atomK->getHybridization()==Atom::SP3D && atomK->getDegree()==5) ){
const Bond *b1 =mol.getBondBetweenAtoms(i,k);
const Bond *b2 =mol.getBondBetweenAtoms(j,k);
// FIX: recognize amide bonds here.
AngleBendContrib *contrib;
int order=0;
switch(atomK->getHybridization()){
case Atom::SP:
order=2;
break;
case Atom::SP3D2:
order=4;
break;
default:
order=0;
break;
}
contrib = new AngleBendContrib(field,i,k,j,
b1->getBondTypeAsDouble(),
b2->getBondTypeAsDouble(),
params[i],params[k],params[j],order);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
}
}
}
}
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
void addTrigonalBipyramidAngles(const Atom *atom,const ROMol &mol, int confId,
const AtomicParamVect &params,
ForceFields::ForceField *field){
PRECONDITION(atom,"bad atom");
PRECONDITION(atom->getHybridization()==Atom::SP3D,"bad hybridization");
PRECONDITION(atom->getDegree()==5,"bad degree");
PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
PRECONDITION(field,"bad forcefield");
const Bond *ax1=0,*ax2=0;
const Bond *eq1=0,*eq2=0,*eq3=0;
const Conformer &conf = mol.getConformer(confId);
//------------------------------------------------------------
// identify the axial and equatorial bonds:
double mostNeg=100.0;
ROMol::OEDGE_ITER beg1,end1;
boost::tie(beg1,end1) = mol.getAtomBonds(atom);
unsigned int aid = atom->getIdx();
while(beg1!=end1){
const Bond *bond1=mol[*beg1].get();
unsigned int oaid = bond1->getOtherAtomIdx(aid);
RDGeom::Point3D v1=conf.getAtomPos(aid).directionVector(conf.getAtomPos(oaid));
ROMol::OEDGE_ITER beg2,end2;
boost::tie(beg2,end2) = mol.getAtomBonds(atom);
while(beg2 != end2){
const Bond *bond2=mol[*beg2].get();
if(bond2->getIdx() > bond1->getIdx()){
unsigned int oaid2 = bond2->getOtherAtomIdx(aid);
RDGeom::Point3D v2=conf.getAtomPos(aid).directionVector(conf.getAtomPos(oaid2));
double dot=v1.dotProduct(v2);
if(dot<mostNeg){
mostNeg = dot;
ax1 = bond1;
ax2 = bond2;
}
}
++beg2;
}
++beg1;
}
CHECK_INVARIANT(ax1,"axial bond not found");
CHECK_INVARIANT(ax2,"axial bond not found");
boost::tie(beg1,end1) = mol.getAtomBonds(atom);
while(beg1!=end1){
const Bond *bond=mol[*beg1].get();
++beg1;
if(bond==ax1 || bond==ax2) continue;
if(!eq1) eq1=bond;
else if(!eq2) eq2=bond;
else if(!eq3) eq3=bond;
}
CHECK_INVARIANT(eq1,"equatorial bond not found");
CHECK_INVARIANT(eq2,"equatorial bond not found");
CHECK_INVARIANT(eq3,"equatorial bond not found");
//------------------------------------------------------------
// alright, add the angles:
AngleBendContrib *contrib;
int atomIdx=atom->getIdx();
int i,j;
// Axial-Axial
i=ax1->getOtherAtomIdx(atomIdx);
j=ax2->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
ax1->getBondTypeAsDouble(),
ax2->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j],2);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
// Equatorial-Equatorial
i=eq1->getOtherAtomIdx(atomIdx);
j=eq2->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
eq1->getBondTypeAsDouble(),
eq2->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j],3);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
i=eq1->getOtherAtomIdx(atomIdx);
j=eq3->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
eq1->getBondTypeAsDouble(),
eq3->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j],3);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
i=eq2->getOtherAtomIdx(atomIdx);
j=eq3->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
eq2->getBondTypeAsDouble(),
eq3->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j],3);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
// Axial-Equatorial
i=ax1->getOtherAtomIdx(atomIdx);
j=eq1->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
ax1->getBondTypeAsDouble(),
eq1->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j]);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
i=ax1->getOtherAtomIdx(atomIdx);
j=eq2->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
ax1->getBondTypeAsDouble(),
eq2->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j]);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
i=ax1->getOtherAtomIdx(atomIdx);
j=eq3->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
ax1->getBondTypeAsDouble(),
eq3->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j]);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
i=ax2->getOtherAtomIdx(atomIdx);
j=eq1->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
ax2->getBondTypeAsDouble(),
eq1->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j]);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
i=ax2->getOtherAtomIdx(atomIdx);
j=eq2->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
ax2->getBondTypeAsDouble(),
eq2->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j]);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
i=ax2->getOtherAtomIdx(atomIdx);
j=eq3->getOtherAtomIdx(atomIdx);
if(params[i]&&params[j]){
contrib = new AngleBendContrib(field,i,atomIdx,j,
ax2->getBondTypeAsDouble(),
eq3->getBondTypeAsDouble(),
params[i],params[atomIdx],params[j]);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
}
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
void addAngleSpecialCases(const ROMol &mol, int confId, const AtomicParamVect &params,
ForceFields::ForceField *field){
PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
PRECONDITION(field,"bad forcefield");
unsigned int nAtoms=mol.getNumAtoms();
for(unsigned int i=0;i<nAtoms;i++){
const Atom *atom = mol.getAtomWithIdx(i);
// trigonal bipyramidal:
if( (atom->getHybridization()==Atom::SP3D && atom->getDegree()==5) ){
addTrigonalBipyramidAngles(atom,mol,confId, params,field);
}
}
}
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
void addNonbonded(const ROMol &mol,int confId,const AtomicParamVect &params,
ForceFields::ForceField *field,boost::shared_array<int> neighborMatrix,
double vdwThresh,bool ignoreInterfragInteractions){
PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
PRECONDITION(field,"bad forcefield");
INT_VECT fragMapping;
if(ignoreInterfragInteractions){
std::vector<ROMOL_SPTR> molFrags=MolOps::getMolFrags(mol,true,&fragMapping);
}
unsigned int nAtoms=mol.getNumAtoms();
const Conformer &conf = mol.getConformer(confId);
for(unsigned int i=0;i<nAtoms;i++){
if(!params[i]) continue;
for(unsigned int j=i+1;j<nAtoms;j++){
if(!params[j] || (ignoreInterfragInteractions && fragMapping[i]!=fragMapping[j])){
continue;
}
if(neighborMatrix[i*nAtoms+j]==-2){
double dist=(conf.getAtomPos(i) - conf.getAtomPos(j)).length();
if(dist <
vdwThresh*Utils::calcNonbondedMinimum(params[i],params[j])){
vdWContrib *contrib;
contrib = new vdWContrib(field,i,j,params[i],params[j]);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
}
}
}
}
}
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
bool okToIncludeTorsion(const ROMol &mol,const Bond *bond,
int idx1,int idx2,int idx3,int idx4){
bool res=true;
RingInfo *rings=mol.getRingInfo();
// having torsions in small rings makes the solver unstable
// and tends to yield poor-quality geometries, so filter those out:
if(rings->isBondInRingOfSize(bond->getIdx(),3)){
res = false;
}// else if(rings->isBondInRingOfSize(bond->getIdx(),4)){
// res = false;
//}
return res;
}
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
void addTorsions(const ROMol &mol,const AtomicParamVect &params,
ForceFields::ForceField *field,
std::string torsionBondSmarts){
PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
PRECONDITION(field,"bad forcefield");
// find all of the torsion bonds:
std::vector<MatchVectType> matchVect;
ROMol *query=SmartsToMol(torsionBondSmarts);
TEST_ASSERT(query);
unsigned int nHits=SubstructMatch(mol,*query,matchVect);
delete query;
for(unsigned int i=0; i<nHits; i++){
MatchVectType match=matchVect[i];
TEST_ASSERT(match.size()==2);
int idx1=match[0].second;
int idx2=match[1].second;
if(!params[idx1]||!params[idx2]) continue;
const Bond *bond=mol.getBondBetweenAtoms(idx1,idx2);
std::vector<TorsionAngleContrib *> contribsHere;
TEST_ASSERT(bond);
const Atom *atom1=mol.getAtomWithIdx(idx1);
const Atom *atom2=mol.getAtomWithIdx(idx2);
if( (atom1->getHybridization()==Atom::SP2||atom1->getHybridization()==Atom::SP3) &&
(atom2->getHybridization()==Atom::SP2||atom2->getHybridization()==Atom::SP3) ){
ROMol::OEDGE_ITER beg1,end1;
boost::tie(beg1,end1) = mol.getAtomBonds(atom1);
while(beg1!=end1){
const Bond *tBond1=mol[*beg1].get();
if(tBond1!=bond){
int bIdx = tBond1->getOtherAtomIdx(idx1);
ROMol::OEDGE_ITER beg2,end2;
boost::tie(beg2,end2) = mol.getAtomBonds(atom2);
while(beg2 != end2){
const Bond *tBond2=mol[*beg2].get();
if(tBond2!=bond && tBond2!=tBond1){
int eIdx=tBond2->getOtherAtomIdx(idx2);
// make sure this isn't a three-membered ring:
if(eIdx != bIdx){
// we now have a torsion involving atoms (bonds):
// bIdx - (tBond1) - idx1 - (bond) - idx2 - (tBond2) - eIdx
TorsionAngleContrib *contrib;
// if either of the end atoms is SP2 hybridized, set a flag
// here.
bool hasSP2=false;
if(mol.getAtomWithIdx(bIdx)->getHybridization()==Atom::SP2 ||
mol.getAtomWithIdx(bIdx)->getHybridization()==Atom::SP2) {
hasSP2 = true;
}
//std::cout << "Torsion: " << bIdx << "-" << idx1 << "-" << idx2 << "-" << eIdx << std::endl;
if(okToIncludeTorsion(mol,bond,bIdx,idx1,idx2,eIdx)){
//std::cout << " INCLUDED" << std::endl;
contrib = new TorsionAngleContrib(field,bIdx,idx1,idx2,eIdx,
bond->getBondTypeAsDouble(),
atom1->getAtomicNum(),
atom2->getAtomicNum(),
atom1->getHybridization(),
atom2->getHybridization(),
params[idx1],params[idx2],
hasSP2);
field->contribs().push_back(ForceFields::ContribPtr(contrib));
contribsHere.push_back(contrib);
}
}
}
beg2++;
}
}
beg1++;
}
}
// now divide the force constant for each contribution to the torsion energy
// about this bond by the number of contribs about this bond:
for(std::vector<TorsionAngleContrib *>::iterator chI=contribsHere.begin();
chI!=contribsHere.end();++chI){
(*chI)->scaleForceConstant(contribsHere.size());
}
}
}
#if 0
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
void addInversions(const ROMol &mol,const AtomicParamVect &params,
ForceFields::ForceField *field){
PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
PRECONDITION(field,"bad forcefield");
unsigned int nAtoms=mol.getNumAtoms();
}
#endif
} // end of namespace Tools
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
ForceFields::ForceField *constructForceField(ROMol &mol,
const AtomicParamVect &params,
double vdwThresh, int confId,
bool ignoreInterfragInteractions){
PRECONDITION(mol.getNumAtoms()==params.size(),"bad parameters");
ForceFields::ForceField *res=new ForceFields::ForceField();
// add the atomic positions:
Conformer &conf = mol.getConformer(confId);
for(unsigned int i=0;i<mol.getNumAtoms();i++){
res->positions().push_back(&conf.getAtomPos(i));
}
Tools::addBonds(mol,params,res);
boost::shared_array<int> neighborMat = Tools::buildNeighborMatrix(mol);
Tools::addAngles(mol,params,res,neighborMat);
Tools::addAngleSpecialCases(mol,confId,params,res);
Tools::addNonbonded(mol,confId,params,res,neighborMat,vdwThresh,ignoreInterfragInteractions);
Tools::addTorsions(mol,params,res);
//Tools::addInversions(mol,params,res);
return res;
}
// ------------------------------------------------------------------------
//
//
//
// ------------------------------------------------------------------------
ForceFields::ForceField *constructForceField(ROMol &mol,double vdwThresh, int confId,
bool ignoreInterfragInteractions){
bool foundAll;
AtomicParamVect params;
boost::tie(params,foundAll)=getAtomTypes(mol);
return constructForceField(mol,params,vdwThresh, confId,ignoreInterfragInteractions);
}
}
}
Reference in New Issue View Git Blame Copy Permalink