MolDrawing.h
Go to the documentation of this file.
// $Id$
//
// Copyright (C) 2009-2012 Greg Landrum
//
// @@ 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.
//
// Includes contributions from Dave Cosgrove (davidacosgroveaz@gmail.com)
//
#ifndef _RD_MOLDRAWING_H_
#define _RD_MOLDRAWING_H_
#include <vector>
#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>
#include <GraphMol/RDKitBase.h>
#include <GraphMol/Depictor/RDDepictor.h>
#include <Geometry/point.h>
/***********
Return Format: vector of ints
RESOLUTION dots_per_angstrom
BOUNDS x1 y1 x2 y2
LINE width dashed atom1_atnum atom2_atnum x1 y1 x2 y2
WEDGE dashed atom1_atnum atom2_atnum x1 y1 x2 y2 x3 y3
ATOM idx atnum x y num_chars char1-charx orient
*************/
namespace RDKit {
namespace Drawing {
typedef int ElementType;
typedef enum {
LINE=1,
WEDGE,
ATOM,
BOUNDS,
RESOLUTION
} PrimType;
typedef enum {
C=0,
N,
E,
S,
W
} OrientType;
namespace detail {
// **************************************************************************
void drawLine( std::vector<ElementType> &res ,
int atnum1 , int atnum2 , int lineWidth , int dashed ,
double x1 , double y1 ,
double x2 , double y2 ) {
res.push_back( LINE );
res.push_back( static_cast<ElementType>(lineWidth) );
res.push_back(dashed);
res.push_back( static_cast<ElementType>(atnum1) );
res.push_back( static_cast<ElementType>(atnum2) );
res.push_back( static_cast<ElementType>(x1) );
res.push_back( static_cast<ElementType>(y1) );
res.push_back( static_cast<ElementType>(x2) );
res.push_back( static_cast<ElementType>(y2) );
}
std::pair<std::string,OrientType> getAtomSymbolAndOrientation(const Atom &atom,RDGeom::Point2D nbrSum){
std::string symbol="";
OrientType orient=C;
int isotope = atom.getIsotope();
if(atom.getAtomicNum()!=6 ||
atom.getFormalCharge()!=0 ||
isotope ||
atom.getNumRadicalElectrons()!=0 ||
atom.hasProp(common_properties::molAtomMapNumber) ||
atom.getDegree()==0 ){
symbol=atom.getSymbol();
bool leftToRight=true;
if(atom.getDegree()==1 && nbrSum.x>0){
leftToRight=false;
}
if(isotope){
symbol = boost::lexical_cast<std::string>(isotope)+symbol;
}
if(atom.hasProp(common_properties::molAtomMapNumber)){
std::string mapNum;
atom.getProp(common_properties::molAtomMapNumber,mapNum);
symbol += ":" + mapNum;
}
int nHs=atom.getTotalNumHs();
if(nHs>0){
std::string h="H";
if(nHs>1) {
h += boost::lexical_cast<std::string>(nHs);
}
if(leftToRight) symbol += h;
else symbol = h+symbol;
}
if( atom.getFormalCharge()!=0 ){
int chg=atom.getFormalCharge();
std::string sgn="+";
if(chg<0){
sgn="-";
}
chg=abs(chg);
if(chg>1){
sgn += boost::lexical_cast<std::string>(chg);
}
if(leftToRight) symbol+=sgn;
else symbol = sgn+symbol;
}
if(atom.getDegree()==1){
double islope=0;
if(fabs(nbrSum.y)>1){
islope=nbrSum.x/fabs(nbrSum.y);
} else {
islope=nbrSum.x;
}
if(fabs(islope)>.85){
if(islope>0){
orient=W;
} else {
orient=E;
}
} else {
if(nbrSum.y>0){
orient=N;
} else {
orient=S;
}
}
}
}
return std::make_pair(symbol,orient);
}
} // end of detail namespace
// **************************************************************************
std::vector<ElementType> DrawMol(const ROMol &mol,int confId=-1,
const std::vector<int> *highlightAtoms=0,
bool includeAtomCircles=false,
unsigned int dotsPerAngstrom=100,
double dblBondOffset=0.3,
double dblBondLengthFrac=0.8,
double angstromsPerChar=0.20
){
if(!mol.getRingInfo()->isInitialized()){
MolOps::findSSSR(mol);
}
std::vector<ElementType> res;
res.push_back(RESOLUTION);
res.push_back(static_cast<ElementType>(dotsPerAngstrom));
const Conformer &conf=mol.getConformer(confId);
const RDGeom::POINT3D_VECT &locs=conf.getPositions();
// get atom symbols and orientations
// (we need them for the bounding box calculation)
std::vector< std::pair<std::string,OrientType> > atomSymbols;
ROMol::VERTEX_ITER bAts,eAts;
boost::tie(bAts,eAts)=mol.getVertices();
while(bAts!=eAts){
ROMol::OEDGE_ITER nbr,endNbrs;
RDGeom::Point2D nbrSum(0,0);
boost::tie(nbr,endNbrs) = mol.getAtomBonds(mol[*bAts].get());
RDGeom::Point2D a1(locs[mol[*bAts]->getIdx()].x,
locs[mol[*bAts]->getIdx()].y);
while(nbr!=endNbrs){
const BOND_SPTR bond=mol[*nbr];
++nbr;
int a2Idx=bond->getOtherAtomIdx(mol[*bAts]->getIdx());
RDGeom::Point2D a2(locs[a2Idx].x,locs[a2Idx].y);
nbrSum+=a2-a1;
}
atomSymbols.push_back(detail::getAtomSymbolAndOrientation(*mol[*bAts],nbrSum));
++bAts;
}
//------------
// do the bounding box
//------------
double minx=1e6,miny=1e6,maxx=-1e6,maxy=-1e6;
for(unsigned int i=0;i<mol.getNumAtoms();++i){
RDGeom::Point3D pt=locs[i];
std::string symbol;
OrientType orient;
boost::tie(symbol,orient)=atomSymbols[i];
if(symbol!=""){
// account for a possible expansion of the bounding box by the symbol
if(pt.x<=minx){
switch(orient){
case C:
case N:
case S:
case E:
minx = pt.x-symbol.size()/2*angstromsPerChar;
break;
case W:
minx = pt.x-symbol.size()*angstromsPerChar;
break;
}
}
if(pt.x>=maxx){
switch(orient){
case C:
case N:
case S:
case W:
maxx = pt.x+symbol.size()/2*angstromsPerChar;
break;
case E:
maxx = pt.x+symbol.size()*angstromsPerChar;
break;
}
}
if(pt.y<=miny){
miny = pt.y-1.5*angstromsPerChar;
}
if(pt.y>=maxy){
maxy = pt.y+angstromsPerChar;
}
} else {
minx=std::min(pt.x,minx);
miny=std::min(pt.y,miny);
maxx=std::max(pt.x,maxx);
maxy=std::max(pt.y,maxy);
}
}
double dimx=(maxx-minx),dimy=(maxy-miny);
res.push_back(BOUNDS);
res.push_back(static_cast<ElementType>(dotsPerAngstrom*0));
res.push_back(static_cast<ElementType>(dotsPerAngstrom*0));
res.push_back(static_cast<ElementType>(dotsPerAngstrom*dimx));
res.push_back(static_cast<ElementType>(dotsPerAngstrom*dimy));
// loop over atoms:
boost::tie(bAts,eAts)=mol.getVertices();
while(bAts!=eAts){
int a1Idx=mol[*bAts]->getIdx();
RDGeom::Point2D a1(locs[a1Idx].x-minx,locs[a1Idx].y-miny);
ROMol::OEDGE_ITER nbr,endNbrs;
RDGeom::Point2D nbrSum(0,0);
boost::tie(nbr,endNbrs) = mol.getAtomBonds(mol[*bAts].get());
while(nbr!=endNbrs){
const BOND_SPTR bond=mol[*nbr];
++nbr;
int a2Idx=bond->getOtherAtomIdx(a1Idx);
int lineWidth=1;
if(highlightAtoms
&&
std::find(highlightAtoms->begin(),highlightAtoms->end(),a1Idx)
!= highlightAtoms->end()
&&
std::find(highlightAtoms->begin(),highlightAtoms->end(),a2Idx)
!= highlightAtoms->end() ){
lineWidth=3;
}
RDGeom::Point2D a2(locs[a2Idx].x-minx,locs[a2Idx].y-miny);
nbrSum+=a2-a1;
if(a2Idx<a1Idx) continue;
// draw bond from a1 to a2.
int atnum1 = mol[*bAts]->getAtomicNum();
int atnum2 = mol.getAtomWithIdx(a2Idx)->getAtomicNum();
if( !mol.getRingInfo()->numBondRings(bond->getIdx()) && bond->getBondType()!=Bond::AROMATIC ) {
// acyclic bonds
RDGeom::Point2D obv=a2-a1;
RDGeom::Point2D perp=obv;
perp.rotate90();
perp.normalize();
if( bond->getBondType()==Bond::DOUBLE || bond->getBondType()==Bond::TRIPLE) {
RDGeom::Point2D startP=a1,endP=a2;
if( bond->getBondType()==Bond::TRIPLE){
perp *= dblBondOffset;
startP+=(obv*(1.-dblBondLengthFrac)/2);
endP-=(obv*(1.-dblBondLengthFrac)/2);
} else {
perp *= 0.5 * dblBondOffset;
}
detail::drawLine( res , atnum1 , atnum2 , lineWidth , 0 ,
dotsPerAngstrom*(startP.x+perp.x) ,
dotsPerAngstrom*(startP.y+perp.y) ,
dotsPerAngstrom*(endP.x+perp.x) ,
dotsPerAngstrom*(endP.y+perp.y) );
if(bond->getBondType() != Bond::AROMATIC){
detail::drawLine( res , atnum1 , atnum2 , lineWidth , 0 ,
dotsPerAngstrom*(startP.x-perp.x) ,
dotsPerAngstrom*(startP.y-perp.y) ,
dotsPerAngstrom*(endP.x-perp.x) ,
dotsPerAngstrom*(endP.y-perp.y) );
} else {
detail::drawLine( res , atnum1 , atnum2 , lineWidth , 1 ,
dotsPerAngstrom*(startP.x-perp.x) ,
dotsPerAngstrom*(startP.y-perp.y) ,
dotsPerAngstrom*(endP.x-perp.x) ,
dotsPerAngstrom*(endP.y-perp.y) );
}
}
if( bond->getBondType()==Bond::SINGLE || bond->getBondType()==Bond::TRIPLE ) {
detail::drawLine( res , atnum1 , atnum2 , lineWidth , 0 ,
dotsPerAngstrom*(a1.x) ,
dotsPerAngstrom*(a1.y) ,
dotsPerAngstrom*(a2.x) ,
dotsPerAngstrom*(a2.y) );
} else if( bond->getBondType()!=Bond::DOUBLE ) {
detail::drawLine( res , atnum1 , atnum2 , lineWidth , 2 ,
dotsPerAngstrom*(a1.x) ,
dotsPerAngstrom*(a1.y) ,
dotsPerAngstrom*(a2.x) ,
dotsPerAngstrom*(a2.y) );
}
} else {
// cyclic bonds
detail::drawLine( res , atnum1 , atnum2 , lineWidth , 0 ,
dotsPerAngstrom*a1.x ,
dotsPerAngstrom*a1.y ,
dotsPerAngstrom*a2.x ,
dotsPerAngstrom*a2.y );
if(bond->getBondType()==Bond::DOUBLE ||
bond->getBondType()==Bond::AROMATIC ||
bond->getBondType()==Bond::TRIPLE ){
RDGeom::Point2D obv=a2-a1;
RDGeom::Point2D perp=obv;
perp.rotate90();
perp.normalize();
if( (bond->getBondType()==Bond::DOUBLE ||
bond->getBondType()==Bond::AROMATIC) &&
mol.getRingInfo()->numBondRings(bond->getIdx())){
// we're in a ring... we might need to flip sides:
ROMol::OEDGE_ITER nbr2,endNbrs2;
boost::tie(nbr2,endNbrs2) = mol.getAtomBonds(mol[*bAts].get());
while(nbr2!=endNbrs2){
const BOND_SPTR bond2=mol[*nbr2];
++nbr2;
if(bond2->getIdx()==bond->getIdx() ||
!mol.getRingInfo()->numBondRings(bond2->getIdx())) continue;
bool sharedRing=false;
BOOST_FOREACH(const INT_VECT &ring,mol.getRingInfo()->bondRings()){
if(std::find(ring.begin(),ring.end(),bond->getIdx())!=ring.end() &&
std::find(ring.begin(),ring.end(),bond2->getIdx())!=ring.end()){
sharedRing=true;
break;
}
}
if(sharedRing){
// these two bonds share a ring.
int a3Idx=bond2->getOtherAtomIdx(a1Idx);
if(a3Idx!=a2Idx){
RDGeom::Point2D a3(locs[a3Idx].x-minx,locs[a3Idx].y-miny);
RDGeom::Point2D obv2=a3-a1;
if(obv2.dotProduct(perp)<0){
perp*=-1;
}
}
}
}
}
perp *= dblBondOffset;
RDGeom::Point2D offsetStart=a1 + obv*(.5*(1.-dblBondLengthFrac));
obv *= dblBondLengthFrac;
detail::drawLine( res , atnum1 , atnum2 , lineWidth , (bond->getBondType()==Bond::AROMATIC),
dotsPerAngstrom*(offsetStart.x+perp.x) ,
dotsPerAngstrom*(offsetStart.y+perp.y) ,
dotsPerAngstrom*(offsetStart.x+obv.x+perp.x) ,
dotsPerAngstrom*(offsetStart.y+obv.y+perp.y) );
}
}
}
std::string symbol;
OrientType orient;
boost::tie(symbol,orient)=atomSymbols[a1Idx];
if(symbol!="" || includeAtomCircles ){
res.push_back(ATOM);
res.push_back(mol[*bAts]->getAtomicNum());
res.push_back(static_cast<ElementType>(dotsPerAngstrom*a1.x));
res.push_back(static_cast<ElementType>(dotsPerAngstrom*a1.y));
res.push_back(static_cast<ElementType>(symbol.length()));
if(symbol.length()){
BOOST_FOREACH(char c, symbol){
res.push_back(static_cast<ElementType>(c));
}
}
res.push_back(static_cast<ElementType>(orient));
}
++bAts;
}
return res;
}
std::vector<int> MolToDrawing(const RDKit::ROMol &mol,const std::vector<int> *highlightAtoms=0,
bool kekulize=true,bool includeAtomCircles=false){
RDKit::RWMol *cp = new RDKit::RWMol(mol);
if(kekulize){
try{
RDKit::MolOps::Kekulize(*cp);
} catch (...) {
delete cp;
cp = new RDKit::RWMol(mol);
}
}
if(!mol.getNumConformers()) {
RDDepict::compute2DCoords(*cp);
}
std::vector<int> drawing=DrawMol(*cp,-1,highlightAtoms,includeAtomCircles);
delete cp;
return drawing;
}
} // end of namespace Drawing
} // end of namespace RDKit
#endif
RDKit::BOND_SPTR
boost::shared_ptr< Bond > BOND_SPTR
Definition: ROMol.h:38
RDKit::Drawing::N
Definition: MolDrawing.h:49
RDKit::MolOps::Kekulize
void Kekulize(RWMol &mol, bool markAtomsBonds=true, unsigned int maxBackTracks=100)
Kekulizes the molecule.
RDKit::ROMol::getNumConformers
unsigned int getNumConformers() const
Definition: ROMol.h:315
RDKit::Drawing::MolToDrawing
std::vector< int > MolToDrawing(const RDKit::ROMol &mol, const std::vector< int > *highlightAtoms=0, bool kekulize=true, bool includeAtomCircles=false)
Definition: MolDrawing.h:408
RDKit::ROMol::getVertices
ATOM_ITER_PAIR getVertices()
returns an iterator pair for looping over all Atoms
RDKit::Drawing::WEDGE
Definition: MolDrawing.h:42
RDKit::MolOps::findSSSR
int findSSSR(const ROMol &mol, std::vector< std::vector< int > > &res)
finds a molecule's Smallest Set of Smallest Rings
RDGeom::Point2D::rotate90
void rotate90()
Definition: point.h:344
RDGeom::Point2D::y
double y
Definition: point.h:265
RDKit::RWMol
RWMol is a molecule class that is intended to be edited.
Definition: RWMol.h:30
RDKit::Atom::getTotalNumHs
unsigned int getTotalNumHs(bool includeNeighbors=false) const
returns the total number of Hs (implicit and explicit) that this Atom is bound to ...
RDKit::Atom::getNumRadicalElectrons
unsigned int getNumRadicalElectrons() const
returns the number of radical electrons for this Atom
Definition: Atom.h:187
point.h
RDKit::Atom::getIsotope
unsigned int getIsotope() const
returns our isotope number
Definition: Atom.h:218
RDKit::common_properties::molAtomMapNumber
const std::string molAtomMapNumber
RDKit::Drawing::C
Definition: MolDrawing.h:48
RDGeom::Point2D
Definition: point.h:263
RDKit::Drawing::S
Definition: MolDrawing.h:51
RDKit::ROMol::getNumAtoms
unsigned int getNumAtoms(bool onlyExplicit=1) const
returns our number of atoms
RDGeom::POINT3D_VECT
std::vector< Point3D > POINT3D_VECT
Definition: point.h:529
RDKit::Atom::getProp
void getProp(const char *key, T &res) const
allows retrieval of a particular property value
Definition: Atom.h:350
RDKitBase.h
pulls in the core RDKit functionality
RDKit::ROMol
ROMol is a molecule class that is intended to have a fixed topology.
Definition: ROMol.h:105
RDKit::Atom::getDegree
unsigned int getDegree() const
RDKit::ROMol::getConformer
const Conformer & getConformer(int id=-1) const
RDKit::Drawing::RESOLUTION
Definition: MolDrawing.h:45
RDKit::Bond::DOUBLE
Definition: Bond.h:58
detail
Definition: types.h:23
RDKit::Drawing::LINE
Definition: MolDrawing.h:41
RDDepictor.h
RDGeom::Point2D::dotProduct
double dotProduct(const Point2D &other) const
Definition: point.h:361
RDKit::Atom::getAtomicNum
int getAtomicNum() const
returns our atomic number
Definition: Atom.h:113
RDKit::INT_VECT
std::vector< int > INT_VECT
Definition: types.h:146
RDGeom::Point2D::normalize
void normalize()
Definition: point.h:338
RDKit::Drawing::BOUNDS
Definition: MolDrawing.h:44
RDKit::Drawing::E
Definition: MolDrawing.h:50
RDKit
Includes a bunch of functionality for handling Atom and Bond queries.
Definition: Atom.h:28
RDKit::Drawing::ATOM
Definition: MolDrawing.h:43
RDKit::Drawing::DrawMol
std::vector< ElementType > DrawMol(const ROMol &mol, int confId=-1, const std::vector< int > *highlightAtoms=0, bool includeAtomCircles=false, unsigned int dotsPerAngstrom=100, double dblBondOffset=0.3, double dblBondLengthFrac=0.8, double angstromsPerChar=0.20)
Definition: MolDrawing.h:145
RDKit::Drawing::OrientType
OrientType
Definition: MolDrawing.h:47
RDKit::RingInfo::bondRings
const VECT_INT_VECT & bondRings() const
returns our bond-rings vectors
Definition: RingInfo.h:124
RDKit::Conformer::getPositions
const RDGeom::POINT3D_VECT & getPositions() const
Get a const reference to the vector of atom positions.
RDKit::Drawing::detail::getAtomSymbolAndOrientation
std::pair< std::string, OrientType > getAtomSymbolAndOrientation(const Atom &atom, RDGeom::Point2D nbrSum)
Definition: MolDrawing.h:73
RDGeom::Point3D::y
double y
Definition: point.h:51
RDGeom::Point3D
Definition: point.h:48
RDKit::ROMol::getRingInfo
RingInfo * getRingInfo() const
Definition: ROMol.h:327
RDKit::Drawing::W
Definition: MolDrawing.h:52
RDKit::Bond::TRIPLE
Definition: Bond.h:59
RDKit::Atom::hasProp
bool hasProp(const char *key) const
returns whether or not we have a property with name key
Definition: Atom.h:383
RDDepict::compute2DCoords
unsigned int compute2DCoords(RDKit::ROMol &mol, const RDGeom::INT_POINT2D_MAP *coordMap=0, bool canonOrient=false, bool clearConfs=true, unsigned int nFlipsPerSample=0, unsigned int nSamples=0, int sampleSeed=0, bool permuteDeg4Nodes=false)
Generate 2D coordinates (a depiction) for a molecule.
RDKit::Conformer
The class for representing 2D or 3D conformation of a molecule.
Definition: Conformer.h:41
RDKit::Drawing::PrimType
PrimType
Definition: MolDrawing.h:40
RDKit::Bond::AROMATIC
Definition: Bond.h:68
RDGeom::Point2D::x
double x
Definition: point.h:265
RDKit::RingInfo::numBondRings
unsigned int numBondRings(unsigned int idx) const
returns the number of rings bond idx is involved in
RDKit::Atom::getSymbol
std::string getSymbol() const
returns our symbol (determined by our atomic number)
RDKit::ROMol::getAtomWithIdx
Atom * getAtomWithIdx(unsigned int idx)
returns a pointer to a particular Atom
RDGeom::Point3D::x
double x
Definition: point.h:51
RDKit::Atom::getFormalCharge
int getFormalCharge() const
returns the formal charge of this atom
Definition: Atom.h:192
RDKit::ROMol::getAtomBonds
OBOND_ITER_PAIR getAtomBonds(Atom const *at) const
provides access to all Bond objects connected to an Atom
RDKit::Drawing::ElementType
int ElementType
Definition: MolDrawing.h:38
RDKit::RingInfo::isInitialized
bool isInitialized() const
checks to see if we've been properly initialized
Definition: RingInfo.h:37
RDKit::Bond::SINGLE
Definition: Bond.h:57
RDKit::Atom
The class for representing atoms.
Definition: Atom.h:67
RDKit::Drawing::detail::drawLine
void drawLine(std::vector< ElementType > &res, int atnum1, int atnum2, int lineWidth, int dashed, double x1, double y1, double x2, double y2)
Definition: MolDrawing.h:57