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