RDKit
Open-source cheminformatics and machine learning.
MolDrawing.h
Go to the documentation of this file.
1 // $Id$
2 //
3 // Copyright (C) 2009-2012 Greg Landrum
4 //
5 // @@ All Rights Reserved @@
6 // This file is part of the RDKit.
7 // The contents are covered by the terms of the BSD license
8 // which is included in the file license.txt, found at the root
9 // of the RDKit source tree.
10 //
11 // Includes contributions from Dave Cosgrove (davidacosgroveaz@gmail.com)
12 //
13 #ifndef _RD_MOLDRAWING_H_
14 #define _RD_MOLDRAWING_H_
15 
16 #include <vector>
17 #include <boost/foreach.hpp>
18 #include <boost/lexical_cast.hpp>
19 #include <GraphMol/RDKitBase.h>
21 #include <Geometry/point.h>
22 
23 /***********
24  Return Format: vector of ints
25 
26  RESOLUTION dots_per_angstrom
27  BOUNDS x1 y1 x2 y2
28  LINE width dashed atom1_atnum atom2_atnum x1 y1 x2 y2
29  WEDGE dashed atom1_atnum atom2_atnum x1 y1 x2 y2 x3 y3
30  ATOM idx atnum x y num_chars char1-charx orient
31 
32 
33 
34 *************/
35 
36 namespace RDKit {
37  namespace Drawing {
38  typedef int ElementType;
39 
40  typedef enum {
41  LINE=1,
46  } PrimType;
47  typedef enum {
48  C=0,
49  N,
50  E,
51  S,
52  W
53  } OrientType;
54 
55  namespace detail {
56  // **************************************************************************
57  void drawLine( std::vector<ElementType> &res ,
58  int atnum1 , int atnum2 , int lineWidth , int dashed ,
59  double x1 , double y1 ,
60  double x2 , double y2 ) {
61 
62  res.push_back( LINE );
63  res.push_back( static_cast<ElementType>(lineWidth) );
64  res.push_back(dashed);
65  res.push_back( static_cast<ElementType>(atnum1) );
66  res.push_back( static_cast<ElementType>(atnum2) );
67  res.push_back( static_cast<ElementType>(x1) );
68  res.push_back( static_cast<ElementType>(y1) );
69  res.push_back( static_cast<ElementType>(x2) );
70  res.push_back( static_cast<ElementType>(y2) );
71 
72  }
73  std::pair<std::string,OrientType> getAtomSymbolAndOrientation(const Atom &atom,RDGeom::Point2D nbrSum){
74  std::string symbol="";
75  OrientType orient=C;
76  int isotope = atom.getIsotope();
77  if(atom.getAtomicNum()!=6 ||
78  atom.getFormalCharge()!=0 ||
79  isotope ||
80  atom.getNumRadicalElectrons()!=0 ||
82  atom.getDegree()==0 ){
83  symbol=atom.getSymbol();
84  bool leftToRight=true;
85  if(atom.getDegree()==1 && nbrSum.x>0){
86  leftToRight=false;
87  }
88  if(isotope){
89  symbol = boost::lexical_cast<std::string>(isotope)+symbol;
90  }
92  std::string mapNum;
94  symbol += ":" + mapNum;
95  }
96  int nHs=atom.getTotalNumHs();
97  if(nHs>0){
98  std::string h="H";
99  if(nHs>1) {
100  h += boost::lexical_cast<std::string>(nHs);
101  }
102  if(leftToRight) symbol += h;
103  else symbol = h+symbol;
104  }
105  if( atom.getFormalCharge()!=0 ){
106  int chg=atom.getFormalCharge();
107  std::string sgn="+";
108  if(chg<0){
109  sgn="-";
110  }
111  chg=abs(chg);
112  if(chg>1){
113  sgn += boost::lexical_cast<std::string>(chg);
114  }
115  if(leftToRight) symbol+=sgn;
116  else symbol = sgn+symbol;
117  }
118 
119  if(atom.getDegree()==1){
120  double islope=0;
121  if(fabs(nbrSum.y)>1){
122  islope=nbrSum.x/fabs(nbrSum.y);
123  } else {
124  islope=nbrSum.x;
125  }
126  if(fabs(islope)>.85){
127  if(islope>0){
128  orient=W;
129  } else {
130  orient=E;
131  }
132  } else {
133  if(nbrSum.y>0){
134  orient=N;
135  } else {
136  orient=S;
137  }
138  }
139  }
140  }
141  return std::make_pair(symbol,orient);
142  }
143  } // end of detail namespace
144  // **************************************************************************
145  std::vector<ElementType> DrawMol(const ROMol &mol,int confId=-1,
146  const std::vector<int> *highlightAtoms=0,
147  bool includeAtomCircles=false,
148  unsigned int dotsPerAngstrom=100,
149  double dblBondOffset=0.3,
150  double dblBondLengthFrac=0.8,
151  double angstromsPerChar=0.20
152  ){
153  if(!mol.getRingInfo()->isInitialized()){
154  MolOps::findSSSR(mol);
155  }
156  std::vector<ElementType> res;
157  res.push_back(RESOLUTION);
158  res.push_back(static_cast<ElementType>(dotsPerAngstrom));
159 
160  const Conformer &conf=mol.getConformer(confId);
161  const RDGeom::POINT3D_VECT &locs=conf.getPositions();
162 
163  // get atom symbols and orientations
164  // (we need them for the bounding box calculation)
165  std::vector< std::pair<std::string,OrientType> > atomSymbols;
166  ROMol::VERTEX_ITER bAts,eAts;
167  boost::tie(bAts,eAts)=mol.getVertices();
168  while(bAts!=eAts){
169  ROMol::OEDGE_ITER nbr,endNbrs;
170  RDGeom::Point2D nbrSum(0,0);
171  boost::tie(nbr,endNbrs) = mol.getAtomBonds(mol[*bAts].get());
172  RDGeom::Point2D a1(locs[mol[*bAts]->getIdx()].x,
173  locs[mol[*bAts]->getIdx()].y);
174  while(nbr!=endNbrs){
175  const BOND_SPTR bond=mol[*nbr];
176  ++nbr;
177  int a2Idx=bond->getOtherAtomIdx(mol[*bAts]->getIdx());
178  RDGeom::Point2D a2(locs[a2Idx].x,locs[a2Idx].y);
179  nbrSum+=a2-a1;
180  }
181  atomSymbols.push_back(detail::getAtomSymbolAndOrientation(*mol[*bAts],nbrSum));
182  ++bAts;
183  }
184 
185 
186  //------------
187  // do the bounding box
188  //------------
189  double minx=1e6,miny=1e6,maxx=-1e6,maxy=-1e6;
190  for(unsigned int i=0;i<mol.getNumAtoms();++i){
191  RDGeom::Point3D pt=locs[i];
192  std::string symbol;
193  OrientType orient;
194  boost::tie(symbol,orient)=atomSymbols[i];
195  if(symbol!=""){
196  // account for a possible expansion of the bounding box by the symbol
197  if(pt.x<=minx){
198  switch(orient){
199  case C:
200  case N:
201  case S:
202  case E:
203  minx = pt.x-symbol.size()/2*angstromsPerChar;
204  break;
205  case W:
206  minx = pt.x-symbol.size()*angstromsPerChar;
207  break;
208  }
209  }
210  if(pt.x>=maxx){
211  switch(orient){
212  case C:
213  case N:
214  case S:
215  case W:
216  maxx = pt.x+symbol.size()/2*angstromsPerChar;
217  break;
218  case E:
219  maxx = pt.x+symbol.size()*angstromsPerChar;
220  break;
221  }
222  }
223 
224  if(pt.y<=miny){
225  miny = pt.y-1.5*angstromsPerChar;
226  }
227  if(pt.y>=maxy){
228  maxy = pt.y+angstromsPerChar;
229  }
230  } else {
231  minx=std::min(pt.x,minx);
232  miny=std::min(pt.y,miny);
233  maxx=std::max(pt.x,maxx);
234  maxy=std::max(pt.y,maxy);
235  }
236  }
237  double dimx=(maxx-minx),dimy=(maxy-miny);
238  res.push_back(BOUNDS);
239  res.push_back(static_cast<ElementType>(dotsPerAngstrom*0));
240  res.push_back(static_cast<ElementType>(dotsPerAngstrom*0));
241  res.push_back(static_cast<ElementType>(dotsPerAngstrom*dimx));
242  res.push_back(static_cast<ElementType>(dotsPerAngstrom*dimy));
243 
244  // loop over atoms:
245  boost::tie(bAts,eAts)=mol.getVertices();
246  while(bAts!=eAts){
247  int a1Idx=mol[*bAts]->getIdx();
248  RDGeom::Point2D a1(locs[a1Idx].x-minx,locs[a1Idx].y-miny);
249  ROMol::OEDGE_ITER nbr,endNbrs;
250  RDGeom::Point2D nbrSum(0,0);
251  boost::tie(nbr,endNbrs) = mol.getAtomBonds(mol[*bAts].get());
252  while(nbr!=endNbrs){
253  const BOND_SPTR bond=mol[*nbr];
254  ++nbr;
255  int a2Idx=bond->getOtherAtomIdx(a1Idx);
256  int lineWidth=1;
257  if(highlightAtoms
258  &&
259  std::find(highlightAtoms->begin(),highlightAtoms->end(),a1Idx)
260  != highlightAtoms->end()
261  &&
262  std::find(highlightAtoms->begin(),highlightAtoms->end(),a2Idx)
263  != highlightAtoms->end() ){
264  lineWidth=3;
265  }
266  RDGeom::Point2D a2(locs[a2Idx].x-minx,locs[a2Idx].y-miny);
267  nbrSum+=a2-a1;
268  if(a2Idx<a1Idx) continue;
269 
270  // draw bond from a1 to a2.
271  int atnum1 = mol[*bAts]->getAtomicNum();
272  int atnum2 = mol.getAtomWithIdx(a2Idx)->getAtomicNum();
273 
274  if( !mol.getRingInfo()->numBondRings(bond->getIdx()) && bond->getBondType()!=Bond::AROMATIC ) {
275  // acyclic bonds
276  RDGeom::Point2D obv=a2-a1;
277  RDGeom::Point2D perp=obv;
278  perp.rotate90();
279  perp.normalize();
280 
281  if( bond->getBondType()==Bond::DOUBLE || bond->getBondType()==Bond::TRIPLE) {
282  RDGeom::Point2D startP=a1,endP=a2;
283  if( bond->getBondType()==Bond::TRIPLE){
284  perp *= dblBondOffset;
285  startP+=(obv*(1.-dblBondLengthFrac)/2);
286  endP-=(obv*(1.-dblBondLengthFrac)/2);
287  } else {
288  perp *= 0.5 * dblBondOffset;
289  }
290  detail::drawLine( res , atnum1 , atnum2 , lineWidth , 0 ,
291  dotsPerAngstrom*(startP.x+perp.x) ,
292  dotsPerAngstrom*(startP.y+perp.y) ,
293  dotsPerAngstrom*(endP.x+perp.x) ,
294  dotsPerAngstrom*(endP.y+perp.y) );
295  if(bond->getBondType() != Bond::AROMATIC){
296  detail::drawLine( res , atnum1 , atnum2 , lineWidth , 0 ,
297  dotsPerAngstrom*(startP.x-perp.x) ,
298  dotsPerAngstrom*(startP.y-perp.y) ,
299  dotsPerAngstrom*(endP.x-perp.x) ,
300  dotsPerAngstrom*(endP.y-perp.y) );
301  } else {
302  detail::drawLine( res , atnum1 , atnum2 , lineWidth , 1 ,
303  dotsPerAngstrom*(startP.x-perp.x) ,
304  dotsPerAngstrom*(startP.y-perp.y) ,
305  dotsPerAngstrom*(endP.x-perp.x) ,
306  dotsPerAngstrom*(endP.y-perp.y) );
307 
308  }
309  }
310  if( bond->getBondType()==Bond::SINGLE || bond->getBondType()==Bond::TRIPLE ) {
311  detail::drawLine( res , atnum1 , atnum2 , lineWidth , 0 ,
312  dotsPerAngstrom*(a1.x) ,
313  dotsPerAngstrom*(a1.y) ,
314  dotsPerAngstrom*(a2.x) ,
315  dotsPerAngstrom*(a2.y) );
316  } else if( bond->getBondType()!=Bond::DOUBLE ) {
317  detail::drawLine( res , atnum1 , atnum2 , lineWidth , 2 ,
318  dotsPerAngstrom*(a1.x) ,
319  dotsPerAngstrom*(a1.y) ,
320  dotsPerAngstrom*(a2.x) ,
321  dotsPerAngstrom*(a2.y) );
322  }
323  } else {
324  // cyclic bonds
325  detail::drawLine( res , atnum1 , atnum2 , lineWidth , 0 ,
326  dotsPerAngstrom*a1.x ,
327  dotsPerAngstrom*a1.y ,
328  dotsPerAngstrom*a2.x ,
329  dotsPerAngstrom*a2.y );
330 
331  if(bond->getBondType()==Bond::DOUBLE ||
332  bond->getBondType()==Bond::AROMATIC ||
333  bond->getBondType()==Bond::TRIPLE ){
334  RDGeom::Point2D obv=a2-a1;
335  RDGeom::Point2D perp=obv;
336  perp.rotate90();
337  perp.normalize();
338 
339  if( (bond->getBondType()==Bond::DOUBLE ||
340  bond->getBondType()==Bond::AROMATIC) &&
341  mol.getRingInfo()->numBondRings(bond->getIdx())){
342  // we're in a ring... we might need to flip sides:
343  ROMol::OEDGE_ITER nbr2,endNbrs2;
344  boost::tie(nbr2,endNbrs2) = mol.getAtomBonds(mol[*bAts].get());
345  while(nbr2!=endNbrs2){
346  const BOND_SPTR bond2=mol[*nbr2];
347  ++nbr2;
348  if(bond2->getIdx()==bond->getIdx() ||
349  !mol.getRingInfo()->numBondRings(bond2->getIdx())) continue;
350  bool sharedRing=false;
351  BOOST_FOREACH(const INT_VECT &ring,mol.getRingInfo()->bondRings()){
352  if(std::find(ring.begin(),ring.end(),bond->getIdx())!=ring.end() &&
353  std::find(ring.begin(),ring.end(),bond2->getIdx())!=ring.end()){
354  sharedRing=true;
355  break;
356  }
357  }
358  if(sharedRing){
359  // these two bonds share a ring.
360  int a3Idx=bond2->getOtherAtomIdx(a1Idx);
361  if(a3Idx!=a2Idx){
362  RDGeom::Point2D a3(locs[a3Idx].x-minx,locs[a3Idx].y-miny);
363  RDGeom::Point2D obv2=a3-a1;
364  if(obv2.dotProduct(perp)<0){
365  perp*=-1;
366  }
367  }
368  }
369  }
370  }
371  perp *= dblBondOffset;
372 
373  RDGeom::Point2D offsetStart=a1 + obv*(.5*(1.-dblBondLengthFrac));
374 
375  obv *= dblBondLengthFrac;
376 
377  detail::drawLine( res , atnum1 , atnum2 , lineWidth , (bond->getBondType()==Bond::AROMATIC),
378  dotsPerAngstrom*(offsetStart.x+perp.x) ,
379  dotsPerAngstrom*(offsetStart.y+perp.y) ,
380  dotsPerAngstrom*(offsetStart.x+obv.x+perp.x) ,
381  dotsPerAngstrom*(offsetStart.y+obv.y+perp.y) );
382  }
383  }
384  }
385  std::string symbol;
386  OrientType orient;
387  boost::tie(symbol,orient)=atomSymbols[a1Idx];
388  if(symbol!="" || includeAtomCircles ){
389  res.push_back(ATOM);
390  res.push_back(mol[*bAts]->getAtomicNum());
391  res.push_back(static_cast<ElementType>(dotsPerAngstrom*a1.x));
392  res.push_back(static_cast<ElementType>(dotsPerAngstrom*a1.y));
393  res.push_back(static_cast<ElementType>(symbol.length()));
394  if(symbol.length()){
395  BOOST_FOREACH(char c, symbol){
396  res.push_back(static_cast<ElementType>(c));
397  }
398  }
399  res.push_back(static_cast<ElementType>(orient));
400  }
401 
402  ++bAts;
403  }
404 
405  return res;
406  }
407 
408  std::vector<int> MolToDrawing(const RDKit::ROMol &mol,const std::vector<int> *highlightAtoms=0,
409  bool kekulize=true,bool includeAtomCircles=false){
410  RDKit::RWMol *cp = new RDKit::RWMol(mol);
411  if(kekulize){
412  try{
414  } catch (...) {
415  delete cp;
416  cp = new RDKit::RWMol(mol);
417  }
418  }
419  if(!mol.getNumConformers()) {
421  }
422  std::vector<int> drawing=DrawMol(*cp,-1,highlightAtoms,includeAtomCircles);
423  delete cp;
424  return drawing;
425  }
426 
427  } // end of namespace Drawing
428 } // end of namespace RDKit
429 
430 #endif
boost::shared_ptr< Bond > BOND_SPTR
Definition: ROMol.h:38
void Kekulize(RWMol &mol, bool markAtomsBonds=true, unsigned int maxBackTracks=100)
Kekulizes the molecule.
unsigned int getNumConformers() const
Definition: ROMol.h:315
std::vector< int > MolToDrawing(const RDKit::ROMol &mol, const std::vector< int > *highlightAtoms=0, bool kekulize=true, bool includeAtomCircles=false)
Definition: MolDrawing.h:408
ATOM_ITER_PAIR getVertices()
returns an iterator pair for looping over all Atoms
int findSSSR(const ROMol &mol, std::vector< std::vector< int > > &res)
finds a molecule&#39;s Smallest Set of Smallest Rings
void rotate90()
Definition: point.h:344
double y
Definition: point.h:265
RWMol is a molecule class that is intended to be edited.
Definition: RWMol.h:30
unsigned int getTotalNumHs(bool includeNeighbors=false) const
returns the total number of Hs (implicit and explicit) that this Atom is bound to ...
unsigned int getNumRadicalElectrons() const
returns the number of radical electrons for this Atom
Definition: Atom.h:187
unsigned int getIsotope() const
returns our isotope number
Definition: Atom.h:218
const std::string molAtomMapNumber
unsigned int getNumAtoms(bool onlyExplicit=1) const
returns our number of atoms
std::vector< Point3D > POINT3D_VECT
Definition: point.h:529
void getProp(const char *key, T &res) const
allows retrieval of a particular property value
Definition: Atom.h:350
pulls in the core RDKit functionality
ROMol is a molecule class that is intended to have a fixed topology.
Definition: ROMol.h:105
unsigned int getDegree() const
const Conformer & getConformer(int id=-1) const
Definition: types.h:23
double dotProduct(const Point2D &other) const
Definition: point.h:361
int getAtomicNum() const
returns our atomic number
Definition: Atom.h:113
std::vector< int > INT_VECT
Definition: types.h:146
void normalize()
Definition: point.h:338
Includes a bunch of functionality for handling Atom and Bond queries.
Definition: Atom.h:28
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
const VECT_INT_VECT & bondRings() const
returns our bond-rings vectors
Definition: RingInfo.h:124
const RDGeom::POINT3D_VECT & getPositions() const
Get a const reference to the vector of atom positions.
std::pair< std::string, OrientType > getAtomSymbolAndOrientation(const Atom &atom, RDGeom::Point2D nbrSum)
Definition: MolDrawing.h:73
double y
Definition: point.h:51
RingInfo * getRingInfo() const
Definition: ROMol.h:327
bool hasProp(const char *key) const
returns whether or not we have a property with name key
Definition: Atom.h:383
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.
The class for representing 2D or 3D conformation of a molecule.
Definition: Conformer.h:41
double x
Definition: point.h:265
unsigned int numBondRings(unsigned int idx) const
returns the number of rings bond idx is involved in
std::string getSymbol() const
returns our symbol (determined by our atomic number)
Atom * getAtomWithIdx(unsigned int idx)
returns a pointer to a particular Atom
double x
Definition: point.h:51
int getFormalCharge() const
returns the formal charge of this atom
Definition: Atom.h:192
OBOND_ITER_PAIR getAtomBonds(Atom const *at) const
provides access to all Bond objects connected to an Atom
bool isInitialized() const
checks to see if we&#39;ve been properly initialized
Definition: RingInfo.h:37
The class for representing atoms.
Definition: Atom.h:67
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