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:42
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.
unsigned int getNumConformers() const
Definition: ROMol.h:363
std::vector< int > MolToDrawing(const RDKit::ROMol &mol, const std::vector< int > *highlightAtoms=0, bool kekulize=true, bool includeAtomCircles=false)
Definition: MolDrawing.h:396
void getProp(const std::string &key, T &res) const
allows retrieval of a particular property value
Definition: RDProps.h:95
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:103
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:188
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:133
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:374
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:114
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:45