RDKit
Open-source cheminformatics and machine learning.
DepictUtils.h
Go to the documentation of this file.
1 //
2 // Copyright (C) 2003-2010 Greg Landrum and Rational Discovery LLC
3 //
4 // @@ All Rights Reserved @@
5 // This file is part of the RDKit.
6 // The contents are covered by the terms of the BSD license
7 // which is included in the file license.txt, found at the root
8 // of the RDKit source tree.
9 //
10 #include <RDGeneral/export.h>
11 #ifndef _RD_DEPICT_UTILS_H_
12 #define _RD_DEPICT_UTILS_H_
13 
14 // REVIEW: remove extra headers here
15 #include <RDGeneral/types.h>
16 #include <GraphMol/RDKitBase.h>
17 #include <GraphMol/RWMol.h>
18 #include <GraphMol/ROMol.h>
19 #include <Geometry/Transform2D.h>
20 #include <Geometry/point.h>
21 #include <queue>
22 
23 namespace RDDepict {
24 RDKIT_DEPICTOR_EXPORT extern double BOND_LEN;
26 RDKIT_DEPICTOR_EXPORT extern double BOND_THRES;
27 RDKIT_DEPICTOR_EXPORT extern double ANGLE_OPEN;
28 RDKIT_DEPICTOR_EXPORT extern unsigned int MAX_COLL_ITERS;
30 RDKIT_DEPICTOR_EXPORT extern unsigned int NUM_BONDS_FLIPS;
31 
32 typedef std::vector<const RDGeom::Point2D *> VECT_C_POINT;
33 
34 typedef std::pair<int, int> PAIR_I_I;
35 typedef std::vector<PAIR_I_I> VECT_PII;
37  bool operator()(const PAIR_I_I &pd1, const PAIR_I_I &pd2) const {
38  return pd1.first > pd2.first;
39  }
40 };
41 
42 typedef std::priority_queue<PAIR_I_I, VECT_PII, gtIIPair> PR_QUEUE;
43 
44 typedef std::pair<double, PAIR_I_I> PAIR_D_I_I;
45 typedef std::list<PAIR_D_I_I> LIST_PAIR_DII;
46 
47 //! Some utility functions used in generating 2D coordinates
48 
49 //! Embed a ring as a convex polygon in 2D
50 /*!
51  The process here is very straightforward:
52 
53  We take the center of the ring to lie at the origin, so put the first
54  point at the origin and then sweep
55  anti-clockwise by an angle A = 360/n for the next point.
56 
57  The length of the arm (l) we want to sweep is easy to compute given the
58  bond length (b) we want to use for each bond in the ring (for now
59  we will assume that this bond legnth is the same for all bonds in the ring:
60 
61  l = b/sqrt(2*(1 - cos(A))
62 
63  the above formula derives from the triangle formula, where side 'c' is given
64  in terms of sides 'a' and 'b' as:
65 
66  c = a^2 + b^2 - 2.a.b.cos(A)
67 
68  where A is the angle between a and b
69  */
71 
73  const RDGeom::Transform2D &trans);
74 
75 //! Find a point that bisects the angle at rcr
76 /*!
77  The new point lies between nb1 and nb2. The line (rcr, newPt) bisects the
78  angle
79  'ang' at rcr
80 */
82  const RDGeom::Point2D &nb1,
83  const RDGeom::Point2D &nb2);
84 
85 //! Reflect a set of point through a the line joining two point
86 /*!
87  ARGUMENTS:
88  \param coordMap a map of <int, point2D> going from atom id to current
89  coordinates of the points that need to be reflected:
90  The coordinates are overwritten
91  \param loc1 the first point of the line that is to be used as a
92  mirror
93  \param loc2 the second point of the line to be used as a mirror
94  */
96  const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2);
97 
99  const RDGeom::Point2D &loc1,
100  const RDGeom::Point2D &loc2);
101 
102 //! Set the neighbors yet to added to aid such that the atoms with the most subs
103 // fall on opposite sides
104 /*!
105  Ok this needs some explanation
106  - Let A, B, C, D be the substituent on the central atom X (given
107  by atom index aid)
108  - also let A be the atom that is already embedded
109  - Now we want the order in which the remaining neighbors B,C,D are
110  added to X such that the atoms with atom with largest number of
111  substituent fall on opposite sides of X so as to minimize atom
112  clashes later in the depiction
113 
114  E.g. let say we have the following situation
115 <pre>
116  B
117  | |
118  A--X--C
119  | |
120  --D--
121  |
122 </pre>
123  In this case the the number substituent of A, B, C, D are 3, 1, 1,
124  4 respectively so want to A and D to go opposite sides and so that
125  we draw
126 <pre>
127  B
128  | | |
129  A--X--D--
130  | | |
131  C
132 </pre>
133  And the correct ordering of the neighbors is B,D,C
134 */
136  const RDKit::ROMol &mol);
137 
138 //! \brief From a given set of rings find the ring the largest common elements
139 // with other rings
140 /*
141  Bit of a weird function - this is typically called once we have embedded some
142  of the
143  rings in a fused system and we are looking for the ring that must be embedded
144  (or merged)
145  next. The heuristic used here is to pick the rings with the maximum number of
146  atoms
147  in common with the rings that are already embedded.
148 
149  \param doneRings a vertor of ring IDs that have been embedded already
150  \param fusedRings list of all the rings in the the fused system
151  \param nextId this is where the ID for the next ring is written
152 
153  \return list of atom ids that are common
154 */
156  const RDKit::VECT_INT_VECT &fusedRings,
157  int &nextId);
158 
159 typedef std::pair<int, int> INT_PAIR;
160 typedef std::vector<INT_PAIR> INT_PAIR_VECT;
161 typedef INT_PAIR_VECT::const_iterator INT_PAIR_VECT_CI;
162 
163 typedef std::pair<double, INT_PAIR> DOUBLE_INT_PAIR;
164 
165 //! Sort a list of atoms by their CIP rank
166 /*!
167  \param mol molecule of interest
168  \param commAtms atoms that need to be ranked
169  \param ascending sort to an ascending order or a descending order
170 */
171 template <class T>
172 RDKIT_DEPICTOR_EXPORT T rankAtomsByRank(const RDKit::ROMol &mol, const T &commAtms,
173  bool ascending = true);
174 
175 //! computes a subangle for an atom of given hybridization and degree
176 /*!
177  \param degree the degree of the atom (number of neighbors)
178  \param htype the atom's hybridization
179 
180  \return the subangle (in radians)
181 */
182 inline double computeSubAngle(unsigned int degree,
184  double angle = M_PI;
185  switch (htype) {
187  case RDKit::Atom::SP3:
188  if (degree == 4) {
189  angle = M_PI / 2;
190  } else {
191  angle = 2 * M_PI / 3;
192  }
193  break;
194  case RDKit::Atom::SP2:
195  angle = 2 * M_PI / 3;
196  break;
197  default:
198  angle = 2. * M_PI / degree;
199  }
200  return angle;
201 }
202 
203 //! computes the rotation direction between two vectors
204 /*!
205 
206  Let:
207 
208  v1 = loc1 - center
209 
210  v2 = loc2 - center
211 
212  If remaining angle(v1, v2) is < 180 and corss(v1, v2) > 0.0 then the rotation
213  dir is +1.0
214 
215  else if remAngle(v1, v2) is > 180 and cross(v1, v2) < 0.0 then rotation dir is
216  -1.0
217 
218  else if remAngle(v1, v2) is < 180 and cross(v1, v2) < 0.0 then rotation dir is
219  -1.0
220 
221  finally if remAngle(v1, v2) is > 180 and cross(v1, v2) < 0.0 then rotation dir
222  is +1.0
223 
224  \param center the common point
225  \param loc1 endpoint 1
226  \param loc2 endpoint 2
227  \param remAngle the remaining angle about center in radians
228 
229  \return the rotation direction (1 or -1)
230 */
231 inline int rotationDir(const RDGeom::Point2D &center,
232  const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2,
233  double remAngle) {
234  RDGeom::Point2D pt1 = loc1 - center;
235  RDGeom::Point2D pt2 = loc2 - center;
236  double cross = pt1.x * pt2.y - pt1.y * pt2.x;
237  double diffAngle = M_PI - remAngle;
238  cross *= diffAngle;
239  if (cross >= 0.0) {
240  return -1;
241  } else {
242  return 1;
243  }
244 }
245 
246 //! computes and return the normal of a vector between two points
247 /*!
248  \param center the common point
249  \param other the endpoint
250 
251  \return the normal
252 */
254  const RDGeom::Point2D &other) {
255  RDGeom::Point2D res = other - center;
256  res.normalize();
257  double tmp = res.x;
258  res.x = -res.y;
259  res.y = tmp;
260  return res;
261 }
262 
263 //! computes the rotation angle between two vectors
264 /*!
265  \param center the common point
266  \param loc1 endpoint 1
267  \param loc2 endpoint 2
268 
269  \return the angle (in radians)
270 */
271 inline double computeAngle(const RDGeom::Point2D &center,
272  const RDGeom::Point2D &loc1,
273  const RDGeom::Point2D &loc2) {
274  RDGeom::Point2D v1 = loc1 - center;
275  RDGeom::Point2D v2 = loc2 - center;
276  return v1.angleTo(v2);
277 }
278 
279 //! \brief pick the ring to embed first in a fused system
280 /*!
281  \param mol the molecule of interest
282  \param fusedRings the collection of the molecule's fused rings
283 
284  \return the index of the ring with the least number of substitutions
285 */
287  const RDKit::VECT_INT_VECT &fusedRings);
288 
289 //! \brief find the rotatable bonds on the shortest path between two atoms
290 //! we will ignore ring atoms, and double bonds which are marked cis/trans
291 /*!
292  <b>Note</b> that rotatable in this context doesn't connect to the
293  standard chemical definition of a rotatable bond; we're just talking
294  about bonds than can be flipped in order to clean up the depiction.
295 
296  \param mol the molecule of interest
297  \param aid1 index of the first atom
298  \param aid2 index of the second atom
299 
300  \return a set of the indices of the rotatable bonds
301 */
303  unsigned int aid2);
304 
305 //! \brief find all the rotatable bonds in a molecule
306 //! we will ignore ring atoms, and double bonds which are marked cis/trans
307 /*!
308  <b>Note</b> that rotatable in this context doesn't connect to the
309  standard chemical definition of a rotatable bond; we're just talking
310  about bonds than can be flipped in order to clean up the depiction.
311 
312  \param mol the molecule of interest
313 
314  \return a set of the indices of the rotatable bonds
315 */
317 
318 //! Get the ids of the atoms and bonds that are connected to aid
319 RDKIT_DEPICTOR_EXPORT void getNbrAtomAndBondIds(unsigned int aid, const RDKit::ROMol *mol,
320  RDKit::INT_VECT &aids, RDKit::INT_VECT &bids);
321 
322 //! Find pairs of bonds that can be permuted at a non-ring degree 4 atom
323 /*!
324  This function will return only those pairs that cannot be
325  permuted by flipping a rotatble bond
326 
327  D
328  |
329  b3
330  |
331  A-b1-B-b2-C
332  |
333  b4
334  |
335  E
336  For example in teh above situation on the pairs (b1, b3) and (b1, b4) will be
337  returned
338  All other permutations can be achieved via a rotatable bond flip.
339 
340  ARGUMENTS:
341  \param center - location of the central atom
342  \param nbrBids - a vector (of length 4) containing the ids of the bonds to
343  the neighbors
344  \param nbrLocs - locations of the neighbors
345 */
347  const RDKit::INT_VECT &nbrBids,
348  const VECT_C_POINT &nbrLocs);
349 
350 //! returns the rank of the atom for determining draw order
351 inline int getAtomDepictRank(const RDKit::Atom *at) {
352  const int maxAtNum = 1000;
353  const int maxDeg = 100;
354  int anum = at->getAtomicNum();
355  anum = anum == 1 ? maxAtNum : anum; // favor non-hydrogen atoms
356  int deg = at->getDegree();
357  return maxDeg * anum + deg;
358 }
359 }
360 
361 #endif
RDKIT_DEPICTOR_EXPORT void getNbrAtomAndBondIds(unsigned int aid, const RDKit::ROMol *mol, RDKit::INT_VECT &aids, RDKit::INT_VECT &bids)
Get the ids of the atoms and bonds that are connected to aid.
std::pair< int, int > PAIR_I_I
Definition: DepictUtils.h:34
RDKIT_DEPICTOR_EXPORT RDKit::INT_VECT getRotatableBonds(const RDKit::ROMol &mol, unsigned int aid1, unsigned int aid2)
find the rotatable bonds on the shortest path between two atoms we will ignore ring atoms...
double y
Definition: point.h:257
std::pair< int, int > INT_PAIR
Definition: DepictUtils.h:159
double computeSubAngle(unsigned int degree, RDKit::Atom::HybridizationType htype)
computes a subangle for an atom of given hybridization and degree
Definition: DepictUtils.h:182
hybridization that hasn&#39;t been specified
Definition: Atom.h:80
double angleTo(const Point2D &other) const
Definition: point.h:353
RDKIT_DEPICTOR_EXPORT RDKit::INT_VECT setNbrOrder(unsigned int aid, const RDKit::INT_VECT &nbrs, const RDKit::ROMol &mol)
Set the neighbors yet to added to aid such that the atoms with the most subs.
RDKIT_DEPICTOR_EXPORT RDKit::INT_VECT getAllRotatableBonds(const RDKit::ROMol &mol)
find all the rotatable bonds in a molecule we will ignore ring atoms, and double bonds which are mark...
Defines the primary molecule class ROMol as well as associated typedefs.
std::pair< double, PAIR_I_I > PAIR_D_I_I
Definition: DepictUtils.h:44
double computeAngle(const RDGeom::Point2D &center, const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2)
computes the rotation angle between two vectors
Definition: DepictUtils.h:271
pulls in the core RDKit functionality
RDKIT_DEPICTOR_EXPORT double BOND_LEN
std::vector< INT_VECT > VECT_INT_VECT
Definition: types.h:261
RDKIT_DEPICTOR_EXPORT int pickFirstRingToEmbed(const RDKit::ROMol &mol, const RDKit::VECT_INT_VECT &fusedRings)
pick the ring to embed first in a fused system
std::vector< const RDGeom::Point2D * > VECT_C_POINT
Definition: DepictUtils.h:32
RDKIT_DEPICTOR_EXPORT RDGeom::Point2D reflectPoint(const RDGeom::Point2D &point, const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2)
RDKIT_DEPICTOR_EXPORT void transformPoints(RDGeom::INT_POINT2D_MAP &nringCor, const RDGeom::Transform2D &trans)
RDKIT_DEPICTOR_EXPORT RDGeom::INT_POINT2D_MAP embedRing(const RDKit::INT_VECT &ring)
Some utility functions used in generating 2D coordinates.
std::map< int, Point2D > INT_POINT2D_MAP
Definition: point.h:512
std::vector< INT_PAIR > INT_PAIR_VECT
Definition: DepictUtils.h:160
std::priority_queue< PAIR_I_I, VECT_PII, gtIIPair > PR_QUEUE
Definition: DepictUtils.h:42
RDKIT_DEPICTOR_EXPORT unsigned int NUM_BONDS_FLIPS
std::vector< int > INT_VECT
Definition: types.h:247
void normalize()
Definition: point.h:325
std::pair< double, INT_PAIR > DOUBLE_INT_PAIR
Definition: DepictUtils.h:163
INT_PAIR_VECT::const_iterator INT_PAIR_VECT_CI
Definition: DepictUtils.h:161
int getAtomDepictRank(const RDKit::Atom *at)
returns the rank of the atom for determining draw order
Definition: DepictUtils.h:351
RDKIT_DEPICTOR_EXPORT void reflectPoints(RDGeom::INT_POINT2D_MAP &coordMap, const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2)
Reflect a set of point through a the line joining two point.
RDKIT_DEPICTOR_EXPORT double COLLISION_THRES
HybridizationType
store hybridization
Definition: Atom.h:79
int getAtomicNum() const
returns our atomic number
Definition: Atom.h:115
RDKIT_DEPICTOR_EXPORT double HETEROATOM_COLL_SCALE
std::vector< PAIR_I_I > VECT_PII
Definition: DepictUtils.h:35
RDGeom::Point2D computeNormal(const RDGeom::Point2D &center, const RDGeom::Point2D &other)
computes and return the normal of a vector between two points
Definition: DepictUtils.h:253
RDKIT_DEPICTOR_EXPORT INT_PAIR_VECT findBondsPairsToPermuteDeg4(const RDGeom::Point2D &center, const RDKit::INT_VECT &nbrBids, const VECT_C_POINT &nbrLocs)
Find pairs of bonds that can be permuted at a non-ring degree 4 atom.
double x
Definition: point.h:257
bool operator()(const PAIR_I_I &pd1, const PAIR_I_I &pd2) const
Definition: DepictUtils.h:37
unsigned int getDegree() const
RDKIT_DEPICTOR_EXPORT RDGeom::Point2D computeBisectPoint(const RDGeom::Point2D &rcr, double ang, const RDGeom::Point2D &nb1, const RDGeom::Point2D &nb2)
Find a point that bisects the angle at rcr.
int rotationDir(const RDGeom::Point2D &center, const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2, double remAngle)
computes the rotation direction between two vectors
Definition: DepictUtils.h:231
RDKIT_DEPICTOR_EXPORT RDKit::INT_VECT findNextRingToEmbed(const RDKit::INT_VECT &doneRings, const RDKit::VECT_INT_VECT &fusedRings, int &nextId)
From a given set of rings find the ring the largest common elements.
RDKIT_DEPICTOR_EXPORT double ANGLE_OPEN
RDKIT_DEPICTOR_EXPORT T rankAtomsByRank(const RDKit::ROMol &mol, const T &commAtms, bool ascending=true)
Sort a list of atoms by their CIP rank.
std::list< PAIR_D_I_I > LIST_PAIR_DII
Definition: DepictUtils.h:45
RDKIT_DEPICTOR_EXPORT unsigned int MAX_COLL_ITERS
The class for representing atoms.
Definition: Atom.h:69
#define M_PI
Definition: MMFF/Params.h:26
Defines the editable molecule class RWMol.
#define RDKIT_DEPICTOR_EXPORT
Definition: export.h:125
RDKIT_DEPICTOR_EXPORT double BOND_THRES