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
const Conformer & getConformer(int id=-1) const
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.
std::vector< int > MolToDrawing(const RDKit::ROMol &mol, const std::vector< int > *highlightAtoms=0, bool kekulize=true, bool includeAtomCircles=false)
Definition: MolDrawing.h:399
ATOM_ITER_PAIR getVertices()
returns an iterator pair for looping over all Atoms
unsigned int numBondRings(unsigned int idx) const
returns the number of rings bond idx is involved in
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 getNumAtoms(bool onlyExplicit=1) const
returns our number of atoms
const VECT_INT_VECT & bondRings() const
returns our bond-rings vectors
Definition: RingInfo.h:126
const std::string molAtomMapNumber
unsigned int getTotalNumHs(bool includeNeighbors=false) const
returns the total number of Hs (implicit and explicit) that this Atom is bound to ...
unsigned int getNumConformers() const
Definition: ROMol.h:366
unsigned int getNumRadicalElectrons() const
returns the number of radical electrons for this Atom
Definition: Atom.h:197
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:103
RingInfo * getRingInfo() const
Definition: ROMol.h:377
const RDGeom::POINT3D_VECT & getPositions() const
Get a const reference to the vector of atom positions.
unsigned int getIsotope() const
returns our isotope number
Definition: Atom.h:229
double dotProduct(const Point2D &other) const
Definition: point.h:347
std::vector< int > INT_VECT
Definition: types.h:191
void normalize()
Definition: point.h:324
int getFormalCharge() const
returns the formal charge of this atom
Definition: Atom.h:203
Std stuff.
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
int getAtomicNum() const
returns our atomic number
Definition: Atom.h:116
bool hasProp(const std::string &key) const
Definition: RDProps.h:116
void getProp(const std::string &key, T &res) const
allows retrieval of a particular property value
Definition: RDProps.h:97
OBOND_ITER_PAIR getAtomBonds(Atom const *at) const
provides access to all Bond objects connected to an Atom
double y
Definition: point.h:47
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 getDegree() const
Atom * getAtomWithIdx(unsigned int idx)
returns a pointer to a particular Atom
double x
Definition: point.h:47
std::string getSymbol() const
returns our symbol (determined by our atomic number)
The class for representing atoms.
Definition: Atom.h:68
bool isInitialized() const
checks to see if we&#39;ve been properly initialized
Definition: RingInfo.h:39
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