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