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