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