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 
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 
72  void transformPoints(RDGeom::INT_POINT2D_MAP &nringCor, 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 angle
77  'ang' at rcr
78  */
80  double ang, const RDGeom::Point2D &nb1,
81  const RDGeom::Point2D &nb2);
82 
83  //! Reflect a set of point through a the line joining two point
84  /*!
85  ARGUMENTS:
86  \param coordMap a map of <int, point2D> going from atom id to current
87  coordinates of the points that need to be reflected:
88  The coordinates are overwritten
89  \param loc1 the first point of the line that is to be used as a mirror
90  \param loc2 the second point of the line to be used as a mirror
91  */
92  void reflectPoints(RDGeom::INT_POINT2D_MAP &coordMap, const RDGeom::Point2D &loc1,
93  const RDGeom::Point2D &loc2);
94 
96  const RDGeom::Point2D &loc2);
97 
98  //! Set the neighbors yet to added to aid such that the atoms with the most subs fall on opposite sides
99  /*!
100  Ok this needs some explanation
101  - Let A, B, C, D be the substituent on the central atom X (given
102  by atom index aid)
103  - also let A be the atom that is already embedded
104  - Now we want the order in which the remaining neighbors B,C,D are
105  added to X such that the atoms with atom with largest number of
106  substituent fall on opposite sides of X so as to minimize atom
107  clashes later in the depiction
108 
109  E.g. let say we have the following situation
110 <pre>
111  B
112  | |
113  A--X--C
114  | |
115  --D--
116  |
117 </pre>
118  In this case the the number substituent of A, B, C, D are 3, 1, 1,
119  4 respectively so want to A and D to go opposite sides and so that
120  we draw
121 <pre>
122  B
123  | | |
124  A--X--D--
125  | | |
126  C
127 </pre>
128  And the correct ordering of the neighbors is B,D,C
129  */
130  RDKit::INT_VECT setNbrOrder(unsigned int aid, const RDKit::INT_VECT &nbrs, const RDKit::ROMol &mol);
131 
132  //! \brief From a given set of rings find the ring the largest common elements with other rings
133  /*
134  Bit of a weird function - this is typically called once we have embedded some of the
135  rings in a fused system and we are looking for the ring that must be embedded (or merged)
136  next. The heuristic used here is to pick the rings with the maximum number of atoms
137  in common with the rings that are already embedded.
138 
139  \param doneRings a vertor of ring IDs that have been embedded already
140  \param fusedRings list of all the rings in the the fused system
141  \param nextId this is where the ID for the next ring is written
142 
143  \return list of atom ids that are common
144  */
146  const RDKit::VECT_INT_VECT &fusedRings,
147  int &nextId);
148 
149  typedef std::pair<int,int> INT_PAIR;
150  typedef std::vector<INT_PAIR> INT_PAIR_VECT;
151  typedef INT_PAIR_VECT::const_iterator INT_PAIR_VECT_CI;
152 
153  typedef std::pair<double, INT_PAIR> DOUBLE_INT_PAIR;
154 
155  //! Sort a list of atoms by their CIP rank
156  /*!
157  \param mol molecule of interest
158  \param commAtms atoms that need to be ranked
159  \param ascending sort to an ascending order or a descending order
160  */
161  template<class T> T rankAtomsByRank(const RDKit::ROMol &mol, const T &commAtms,
162  bool ascending=true);
163 
164 
165  //! computes a subangle for an atom of given hybridization and degree
166  /*!
167  \param degree the degree of the atom (number of neighbors)
168  \param htype the atom's hybridization
169 
170  \return the subangle (in radians)
171  */
172  inline double computeSubAngle(unsigned int degree, RDKit::Atom::HybridizationType htype) {
173  double angle = M_PI;
174  switch (htype) {
176  case RDKit::Atom::SP3 :
177  if (degree == 4) {
178  angle = M_PI/2;
179  } else {
180  angle = 2*M_PI/3;
181  }
182  break;
183  case RDKit::Atom::SP2 :
184  angle = 2*M_PI/3;
185  break;
186  default:
187  angle = 2.*M_PI/degree;
188  }
189  return angle;
190  }
191 
192  //! computes the rotation direction between two vectors
193  /*!
194 
195  Let:
196 
197  v1 = loc1 - center
198 
199  v2 = loc2 - center
200 
201  If remaining angle(v1, v2) is < 180 and corss(v1, v2) > 0.0 then the rotation dir is +1.0
202 
203  else if remAngle(v1, v2) is > 180 and cross(v1, v2) < 0.0 then rotation dir is -1.0
204 
205  else if remAngle(v1, v2) is < 180 and cross(v1, v2) < 0.0 then rotation dir is -1.0
206 
207  finally if remAngle(v1, v2) is > 180 and cross(v1, v2) < 0.0 then rotation dir is +1.0
208 
209  \param center the common point
210  \param loc1 endpoint 1
211  \param loc2 endpoint 2
212  \param remAngle the remaining angle about center in radians
213 
214  \return the rotation direction (1 or -1)
215  */
216  inline int rotationDir(const RDGeom::Point2D &center, const RDGeom::Point2D &loc1,
217  const RDGeom::Point2D &loc2, double remAngle) {
218  RDGeom::Point2D pt1 = loc1 - center;
219  RDGeom::Point2D pt2 = loc2 - center;
220  double cross = pt1.x*pt2.y - pt1.y*pt2.x;
221  double diffAngle = M_PI - remAngle;
222  cross *= diffAngle;
223  if (cross >= 0.0) {
224  return -1;
225  } else {
226  return 1;
227  }
228  }
229 
230  //! computes and return the normal of a vector between two points
231  /*!
232  \param center the common point
233  \param other the endpoint
234 
235  \return the normal
236  */
238  const RDGeom::Point2D &other) {
239  RDGeom::Point2D res = other - center;
240  res.normalize();
241  double tmp = res.x;
242  res.x = -res.y;
243  res.y = tmp;
244  return res;
245  }
246 
247  //! computes the rotation angle between two vectors
248  /*!
249  \param center the common point
250  \param loc1 endpoint 1
251  \param loc2 endpoint 2
252 
253  \return the angle (in radians)
254  */
255  inline double computeAngle(const RDGeom::Point2D &center,
256  const RDGeom::Point2D &loc1,
257  const RDGeom::Point2D &loc2) {
258  RDGeom::Point2D v1 = loc1 - center;
259  RDGeom::Point2D v2 = loc2 - center;
260  return v1.angleTo(v2);
261  }
262 
263  //! \brief pick the ring to embed first in a fused system
264  /*!
265  \param mol the molecule of interest
266  \param fusedRings the collection of the molecule's fused rings
267 
268  \return the index of the ring with the least number of substitutions
269  */
270  int pickFirstRingToEmbed(const RDKit::ROMol &mol, const RDKit::VECT_INT_VECT &fusedRings);
271 
272  //! \brief find the rotatable bonds on the shortest path between two atoms
273  //! we will ignore ring atoms, and double bonds which are marked cis/trans
274  /*!
275  <b>Note</b> that rotatable in this context doesn't connect to the
276  standard chemical definition of a rotatable bond; we're just talking
277  about bonds than can be flipped in order to clean up the depiction.
278 
279  \param mol the molecule of interest
280  \param aid1 index of the first atom
281  \param aid2 index of the second atom
282 
283  \return a set of the indices of the rotatable bonds
284  */
285  RDKit::INT_VECT getRotatableBonds(const RDKit::ROMol &mol, unsigned int aid1,
286  unsigned int aid2);
287 
288  //! \brief find all the rotatable bonds in a molecule
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 
297  \return a set of the indices of the rotatable bonds
298  */
300 
301  //! Get the ids of the atoms and bonds that are connected to aid
302  void getNbrAtomAndBondIds(unsigned int aid, const RDKit::ROMol *mol,
303  RDKit::INT_VECT &aids, RDKit::INT_VECT &bids);
304 
305  //! Find pairs of bonds that can be permuted at a non-ring degree 4 atom
306  /*!
307  This function will return only those pairs that cannot be
308  permuted by flipping a rotatble bond
309 
310  D
311  |
312  b3
313  |
314  A-b1-B-b2-C
315  |
316  b4
317  |
318  E
319  For example in teh above situation on the pairs (b1, b3) and (b1, b4) will be returned
320  All other permutations can be achieved via a rotatable bond flip.
321 
322  ARGUMENTS:
323  \param center - location of the central atom
324  \param nbrBids - a vector (of length 4) containing the ids of the bonds to the neighbors
325  \param nbrLocs - locations of the neighbors
326  */
327  INT_PAIR_VECT findBondsPairsToPermuteDeg4(const RDGeom::Point2D &center, const RDKit::INT_VECT &nbrBids,
328  const VECT_C_POINT &nbrLocs);
329 }
330 
331 #endif
unsigned int MAX_COLL_ITERS
std::pair< int, int > PAIR_I_I
Definition: DepictUtils.h:33
std::pair< int, int > INT_PAIR
Definition: DepictUtils.h:149
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:265
double computeSubAngle(unsigned int degree, RDKit::Atom::HybridizationType htype)
computes a subangle for an atom of given hybridization and degree
Definition: DepictUtils.h:172
hybridization that hasn&#39;t been specified
Definition: Atom.h:80
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:255
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:105
std::vector< INT_VECT > VECT_INT_VECT
Definition: types.h:160
std::vector< const RDGeom::Point2D * > VECT_C_POINT
Definition: DepictUtils.h:31
std::map< int, Point2D > INT_POINT2D_MAP
Definition: point.h:533
std::vector< INT_PAIR > INT_PAIR_VECT
Definition: DepictUtils.h:150
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
std::vector< int > INT_VECT
Definition: types.h:146
void normalize()
Definition: point.h:338
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 fall on opposite sides...
void transformPoints(RDGeom::INT_POINT2D_MAP &nringCor, const RDGeom::Transform2D &trans)
std::pair< double, INT_PAIR > DOUBLE_INT_PAIR
Definition: DepictUtils.h:153
INT_PAIR_VECT::const_iterator INT_PAIR_VECT_CI
Definition: DepictUtils.h:151
double angleTo(const Point2D &other) const
Definition: point.h:366
HybridizationType
store hybridization
Definition: Atom.h:79
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:237
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:265
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 with other rings.
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:216
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...
#define M_PI
Definition: MMFF/Params.h:25
Defines the editable molecule class RWMol.