RDKit
Open-source cheminformatics and machine learning.
EmbeddedFrag.h
Go to the documentation of this file.
1 //
2 // Copyright (C) 2003-2017 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_EMBEDDED_FRAG_H_
11 #define _RD_EMBEDDED_FRAG_H_
12 
13 #include <RDGeneral/types.h>
14 #include <Geometry/Transform2D.h>
15 #include <Geometry/point.h>
16 #include "DepictUtils.h"
17 #include <boost/smart_ptr.hpp>
18 
19 namespace RDKit {
20 class ROMol;
21 class Bond;
22 }
23 
24 namespace RDDepict {
25 typedef boost::shared_array<double> DOUBLE_SMART_PTR;
26 
27 //! Class that contains the data for an atoms that has alredy been embedded
28 class EmbeddedAtom {
29  public:
30  typedef enum { UNSPECIFIED = 0, CISTRANS, RING } EAtomType;
31 
33  : aid(0),
34  angle(-1.0),
35  nbr1(-1),
36  nbr2(-1),
37  CisTransNbr(-1),
38  ccw(true),
39  rotDir(0),
40  d_density(-1.0),
41  df_fixed(false) {
42  neighs.clear();
43  }
44 
45  EmbeddedAtom(unsigned int aid, const RDGeom::Point2D &pos)
46  : aid(aid),
47  angle(-1.0),
48  nbr1(-1),
49  nbr2(-1),
50  CisTransNbr(-1),
51  ccw(true),
52  rotDir(0),
53  d_density(-1.0),
54  df_fixed(false) {
55  loc = pos;
56  }
57 
59  loc = other.loc;
60  angle = other.angle;
61  nbr1 = other.nbr1;
62  nbr2 = other.nbr2;
63  CisTransNbr = other.CisTransNbr;
64  rotDir = other.rotDir;
65  normal = other.normal;
66  ccw = other.ccw;
67  neighs = other.neighs;
68  d_density = other.d_density;
69  df_fixed = other.df_fixed;
70  return *this;
71  }
72 
73  void Transform(const RDGeom::Transform2D &trans) {
74  RDGeom::Point2D temp = loc + normal;
75  trans.TransformPoint(loc);
76  trans.TransformPoint(temp);
77  normal = temp - loc;
78  }
79 
80  void Reflect(const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2) {
81  RDGeom::Point2D temp = loc + normal;
82  loc = reflectPoint(loc, loc1, loc2);
83  temp = reflectPoint(temp, loc1, loc2);
84  normal = temp - loc;
85  ccw = (!ccw);
86  }
87 
88  unsigned int aid; // the id of the atom
89 
90  //! the angle that is already takes at this atom, so any new atom attaching to
91  // this
92  //! atom with have to fall in the available part
93  double angle;
94 
95  //! the first neighbor of this atom that form the 'angle'
96  int nbr1;
97 
98  //! the second neighbor of atom that from the 'angle'
99  int nbr2;
100 
101  //! is this is a cis/trans atom the neighbor of this atom that is involved in
102  // the
103  //! cis/trans system - defaults to -1
105 
106  //! which direction do we rotate this normal to add the next bond
107  //! if ccw is true we rotate counter cloack wise, otherwise rotate clock wise,
108  // by an angle that is
109  //! <= PI/2
110  bool ccw;
111 
112  //! rotation direction around this atom when adding new atoms,
113  //! we determine this for the first neighbor and stick to this direction after
114  // that
115  //! useful only on atoms that are degree >= 4
116  int rotDir;
117 
118  RDGeom::Point2D loc; // the current location of this atom
119  //! this is a normal vector to one of the bonds that added this atom
120  //! it provides the side on which we want to add a new bond to this atom,
121  //! this is only relevant when we are dealing with non ring atoms. We would
122  // like to draw chains in
123  //! a zig-zag manner
125 
126  //! and these are the atom IDs of the neighbors that still need to be embedded
128 
129  // density of the atoms around this atoms
130  // - this is sum of inverse of the square of distances to other atoms from
131  // this atom
132  // Used in the collision removal code - initialized to -1.0
133  double d_density;
134 
135  //! if set this atom is fixed: further operations on the fragment may not
136  //! move it.
137  bool df_fixed;
138 };
139 
140 typedef std::map<unsigned int, EmbeddedAtom> INT_EATOM_MAP;
141 typedef INT_EATOM_MAP::iterator INT_EATOM_MAP_I;
142 typedef INT_EATOM_MAP::const_iterator INT_EATOM_MAP_CI;
143 
144 //! Class containing a fragment of a molecule that has already been embedded
145 /*
146  Here is how this class is designed to be used
147  - find a set of fused rings and compte the coordinates for the atoms in those
148  ring
149  - them grow this sytem either by adding non ring neighbors
150  - or by adding other embedded fragment
151  - so at the end of the process the whole molecule end up being one these
152  embedded frag objects
153 */
155  // REVIEW: think about moving member functions up to global level and just
156  // using
157  // this class as a container
158 
159  public:
160  //! Default constructor
161  EmbeddedFrag() : d_done(false), dp_mol(0) {
162  d_eatoms.clear();
163  d_attachPts.clear();
164  };
165 
166  //! Intializer from a single atom id
167  /*!
168  A single Embedded Atom with this atom ID is added and placed at the origin
169  */
170  EmbeddedFrag(unsigned int aid, const RDKit::ROMol *mol);
171 
172  //! Constructor when the coordinates have been specified for a set of atoms
173  /*!
174  This simply initialized a set of EmbeddedAtom to have the same coordinates
175  as the
176  one's specified. No testing is done to verify any kind of ocrrectlness.
177  Also
178  this fragment is less ready (to expand and add new neighbors) than when
179  using other
180  constructors. This is because:
181  - the user may have specified coords for only a part of the atoms in a
182  fused ring systems
183  in which case we need to find these atoms and merge these ring systems to
184  this fragment
185  - The atoms are not yet aware of their neighbor (what is left to add etc.)
186  this again
187  depends on atoms properly so that new
188  neighbors can be added to them
189  */
190  EmbeddedFrag(const RDKit::ROMol *mol,
191  const RDGeom::INT_POINT2D_MAP &coordMap);
192 
193  //! Initializer from a set of fused rings
194  /*!
195  ARGUMENTS:
196  \param mol the molecule of interest
197  \param fusedRings a vector of rings, each ring is a list of atom ids
198  */
199  EmbeddedFrag(const RDKit::ROMol *mol, const RDKit::VECT_INT_VECT &fusedRings);
200 
201  //! Initializer for a cis/trans system using the double bond
202  /*!
203  ARGUMENTS:
204  \param dblBond the double bond that is involed in the cis/trans
205  configuration
206  */
207  explicit EmbeddedFrag(const RDKit::Bond *dblBond);
208 
209  //! Expand this embedded system by adding negihboring atoms or other embedded
210  // systems
211  /*!
212 
213  Note that both nratms and efrags are modified in this function
214  as we start merging them with the current fragment
215 
216  */
217  void expandEfrag(RDKit::INT_LIST &nratms, std::list<EmbeddedFrag> &efrags);
218 
219  //! Add a new non-ring atom to this object
220  /*
221  ARGUMENTS:
222  \param aid ID of the atom to be added
223  \param toAid ID of the atom that is already in this object to which this
224  atom is added
225  */
226  void addNonRingAtom(unsigned int aid, unsigned int toAid);
227 
228  //! Merge this embedded object with another embedded fragment
229  /*!
230 
231  The transformation (rotation + translation required to attached
232  the passed in object will be computed and applied. The
233  coordinates of the atoms in this object will remain fixed We
234  will assume that there are no common atoms between the two
235  fragments to start with
236 
237  ARGUMENTS:
238  \param embObj another EmbeddedFrag object to be merged with this object
239  \param toAid the atom in this embedded fragment to which the new object
240  will be attached
241  \param nbrAid the atom in the other fragment to attach to
242  */
243  void mergeNoCommon(EmbeddedFrag &embObj, unsigned int toAid,
244  unsigned int nbrAid);
245 
246  //! Merge this embedded object with another embedded fragment
247  /*!
248 
249  The transformation (rotation + translation required to attached
250  the passed in object will be computed and applied. The
251  coordinates of the atoms in this object will remain fixed This
252  already know there are a atoms in common and we will use them to
253  merge things
254 
255  ARGUMENTS:
256  \param embObj another EmbeddedFrag object to be merged with this object
257  \param commAtms a vector of ids of the common atoms
258 
259  */
260  void mergeWithCommon(EmbeddedFrag &embObj, RDKit::INT_VECT &commAtms);
261 
262  void mergeFragsWithComm(std::list<EmbeddedFrag> &efrags);
263 
264  //! Mark this fragment to be done for final embedding
265  void markDone() { d_done = true; }
266 
267  //! If this fragment done for the final embedding
268  bool isDone() { return d_done; }
269 
270  //! Get the molecule that this embedded fragmetn blongs to
271  const RDKit::ROMol *getMol() const { return dp_mol; }
272 
273  //! Find the common atom ids between this fragment and a second one
274  RDKit::INT_VECT findCommonAtoms(const EmbeddedFrag &efrag2);
275 
276  //! Find a neighbor to a non-ring atom among the already embedded atoms
277  /*!
278  ARGUMENTS:
279  \param aid the atom id of interest
280 
281  RETURNS:
282  \return the id of the atom if we found a neighbor
283  -1 otherwise
284 
285  NOTE: by definition we can have only one neighbor in the embdded system.
286  */
287  int findNeighbor(unsigned int aid);
288 
289  //! Tranform this object to a new coordinates system
290  /*!
291  ARGUMENTS:
292  \param trans : the transformation that need to be applied to the atoms in
293  this object
294  */
295  void Transform(const RDGeom::Transform2D &trans);
296 
297  void Reflect(const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2);
298 
299  const INT_EATOM_MAP &GetEmbeddedAtoms() const { return d_eatoms; }
300 
301  void Translate(const RDGeom::Point2D &shift) {
302  INT_EATOM_MAP_I eari;
303  for (eari = d_eatoms.begin(); eari != d_eatoms.end(); eari++) {
304  eari->second.loc += shift;
305  }
306  }
307 
308  EmbeddedAtom GetEmbeddedAtom(unsigned int aid) const {
309  INT_EATOM_MAP_CI posi = d_eatoms.find(aid);
310  if (posi == d_eatoms.end()) {
311  PRECONDITION(0, "Embedded atom does not contain embedded atom specified");
312  }
313  return posi->second;
314  }
315 
316  //! the number of atoms in the embedded system
317  int Size() const { return d_eatoms.size(); }
318 
319  //! \brief compute a box that encloses the fragment
320  void computeBox();
321 
322  //! \brief Flip atoms on one side of a bond - used in removing colissions
323  /*!
324  ARGUMENTS:
325  \param bondId - the bond used as the mirror to flip
326  \param flipEnd - flip the atoms at the end of the bond
327 
328  */
329  void flipAboutBond(unsigned int bondId, bool flipEnd = true);
330 
331  void openAngles(const double *dmat, unsigned int aid1, unsigned int aid2);
332 
333  std::vector<PAIR_I_I> findCollisions(const double *dmat,
334  bool includeBonds = 1);
335 
336  void computeDistMat(DOUBLE_SMART_PTR &dmat);
337 
338  double mimicDistMatAndDensityCostFunc(const DOUBLE_SMART_PTR *dmat,
339  double mimicDmatWt);
340 
341  void permuteBonds(unsigned int aid, unsigned int aid1, unsigned int aid2);
342 
343  void randomSampleFlipsAndPermutations(unsigned int nBondsPerSample = 3,
344  unsigned int nSamples = 100,
345  int seed = 100,
346  const DOUBLE_SMART_PTR *dmat = 0,
347  double mimicDmatWt = 0.0,
348  bool permuteDeg4Nodes = false);
349 
350  //! Remove collisions in a structure by flipping rotable bonds
351  //! along the shortest path between two colliding atoms
352  void removeCollisionsBondFlip();
353 
354  //! Remove collision by opening angles at the offending atoms
355  void removeCollisionsOpenAngles();
356 
357  //! Remove collisions by shortening bonds along the shortest path between the
358  // atoms
359  void removeCollisionsShortenBonds();
360 
361  //! helpers funtions to
362 
363  //! \brief make list of neighbors for each atom in the embedded system that
364  //! still need to be embedded
365  void setupNewNeighs();
366 
367  //! update the unembedded neighbor atom list for a specified atom
368  void updateNewNeighs(unsigned int aid);
369 
370  //! \brief Find all atoms in this embedded system that are
371  //! within a specified distant of a point
372  int findNumNeigh(const RDGeom::Point2D &pt, double radius);
373 
374  inline double getBoxPx() { return d_px; }
375  inline double getBoxNx() { return d_nx; }
376  inline double getBoxPy() { return d_py; }
377  inline double getBoxNy() { return d_ny; }
378 
379  void canonicalizeOrientation();
380 
381  private:
382  double totalDensity();
383 
384  void embedFusedRings(const RDKit::VECT_INT_VECT &fusedRings);
385 
386  //! \brief Find a transform to join a ring to the current embedded frag when
387  // we
388  //! have only on common atom
389  /*!
390  So this is the state of affairs assumed here:
391  - we already have some rings in the fused system embeded and the
392  coordinates for the atoms
393  - the coordinates for the atoms in the new ring (with the center
394  of rings at the origin) are available nringCors. we want to
395  translate and rotate this ring to join with the already
396  embeded rings.
397  - only one atom is common between this new ring and the atoms
398  that are already embedded
399  - so we need to compute a transform that includes a translation
400  so that the common atom overlaps and the rotation to minimize
401  overalp with other atoms.
402 
403  Here's what is done:
404  - we bisect the remaining sweep angle at the common atom and
405  attach the new ring such that the center of the new ring falls
406  on this bisecting line
407 
408  NOTE: It is assumed here that the original coordinates for the
409  new ring are such that the center is at the origin (this is the
410  way rings come out of embedRing)
411  */
412  RDGeom::Transform2D computeOneAtomTrans(unsigned int commAid,
413  const EmbeddedFrag &other);
414 
415  RDGeom::Transform2D computeTwoAtomTrans(
416  unsigned int aid1, unsigned int aid2,
417  const RDGeom::INT_POINT2D_MAP &nringCor);
418 
419  //! Merge a ring with already embedded atoms
420  /*!
421  It is assumed that the new rings has already been oriented
422  correctly etc. This function just update all the relevant data,
423  like the neighbor information and the sweep angle
424  */
425  void mergeRing(const EmbeddedFrag &embRing, unsigned int nCommon,
426  const RDKit::INT_VECT &pinAtoms);
427 
428  //! Reflect a fragment if necessary through a line connecting two atoms
429  /*!
430 
431  We want add the new fragment such that, most of its atoms fall
432  on the side opoiste to where the atoms already embedded are aid1
433  and aid2 give the atoms that were used to align the new ring to
434  the embedded atoms and we will assume that that process has
435  already taken place (i.e. transformRing has been called)
436 
437  */
438  void reflectIfNecessaryDensity(EmbeddedFrag &embFrag, unsigned int aid1,
439  unsigned int aid2);
440 
441  //! Reflect a fragment if necessary based on the cis/trans specification
442  /*!
443 
444  we want to add the new fragment such that the cis/trans
445  specification on bond between aid1 and aid2 is not violated. We
446  will assume that aid1 and aid2 from this fragments as well as
447  embFrag are already aligned to each other.
448 
449  \param embFrag the fragment that will be reflected if necessary
450  \param ctCase which fragment if the cis/trans dbl bond
451  - 1 means embFrag is the cis/trans fragment
452  - 2 mean "this" is the cis/trans fragment
453  \param aid1 first atom that forms the plane (line) of reflection
454  \param aid2 seconf atom that forms the plane of reflection
455  */
456  void reflectIfNecessaryCisTrans(EmbeddedFrag &embFrag, unsigned int ctCase,
457  unsigned int aid1, unsigned int aid2);
458 
459  //! Reflect a fragment if necessary based on a thrid common point
460  /*!
461 
462  we want add the new fragment such that the thrid point falls on
463  the same side of aid1 and aid2. We will assume that aid1 and
464  aid2 from this fragments as well as embFrag are already aligned
465  to each other.
466 
467  */
468  void reflectIfNecessaryThirdPt(EmbeddedFrag &embFrag, unsigned int aid1,
469  unsigned int aid2, unsigned int aid3);
470 
471  //! \brief Initialize this fragment from a ring and coordinates for its atoms
472  /*!
473  ARGUMENTS:
474  /param ring a vector of atom ids in the ring; it is assumed that there
475  in
476  clockwise or anti-clockwise order
477  /param nringMap a map of atomId to coordinate map for the atoms in the ring
478  */
479  void initFromRingCoords(const RDKit::INT_VECT &ring,
480  const RDGeom::INT_POINT2D_MAP &nringMap);
481 
482  //! Helper function to addNonRingAtom to a specified atoms in the fragment
483  /*
484  Add an atom to this embedded fragment when the fragment already
485  has a atleast two previously added neighbors to 'toAid'. In this
486  case we have to choose where the the new neighbor goes based on
487  the angle that is already taken around the atom.
488 
489  ARGUMENTS:
490  \param aid ID of the atom to be added
491  \param toAid ID of the atom that is already in this object to which this
492  atom is added
493  */
494  void addAtomToAtomWithAng(unsigned int aid, unsigned int toAid);
495 
496  //! Helper function to addNonRingAtom to a specified atoms in the fragment
497  /*!
498 
499  Add an atom (aid) to an atom (toAid) in this embedded fragment
500  when 'toAid' has one or no neighbors previously added. In this
501  case where the new atom should fall is determined by the degree
502  of 'toAid' and the congestion around it.
503 
504  ARGUMENTS:
505  \param aid ID of the atom to be added
506  \param toAid ID of the atom that is already in this object to which this
507  atom is added
508  \param mol the molecule we are dealing with
509  */
510  void addAtomToAtomWithNoAng(
511  unsigned int aid,
512  unsigned int toAid); //, const RDKit::ROMol *mol);
513 
514  //! Helper funtion to contructor that takes predefined coordinates
515  /*!
516 
517  Given an atom with more than 2 neighbors all embedded in this
518  fragment this function tries to determine
519 
520  - how much of an angle if left for any new neighbors yet to be
521  added
522  - which atom should we rotate when we add a new neighbor and in
523  which direction (clockwise or anticlockwise
524 
525  This is how it works
526  - find the pair of nbrs that have the largest angle
527  - this will most likely be the angle that is available - unless
528  we have fused rings and we found on of the ring angle !!!! -
529  in this cae we find the next best
530  - find the smallest anngle that contains one of these nbrs -
531  this determined which
532  - way we want to rotate
533 
534  ARGUMENTS:
535  \param aid the atom id where we are centered right now
536  \param doneNbrs list of neighbors that are already embedded around aid
537  */
538  void computeNbrsAndAng(unsigned int aid, const RDKit::INT_VECT &doneNbrs);
539  // const RDKit::ROMol *mol);
540 
541  //! are we embedded with the final (molecule) coordinates
542  bool d_done;
543  double d_px, d_nx, d_py, d_ny;
544 
545  //! a map that takes one from teh atom id to the embeddedatom object for that
546  // atom.
547  INT_EATOM_MAP d_eatoms;
548 
549  // RDKit::INT_DEQUE d_attachPts;
550  RDKit::INT_LIST d_attachPts;
551 
552  // pointer to the owning molecule
553  const RDKit::ROMol *dp_mol;
554 };
555 }
556 
557 #endif
std::list< int > INT_LIST
Definition: types.h:196
EmbeddedAtom & operator=(const EmbeddedAtom &other)
Definition: EmbeddedFrag.h:58
const RDKit::ROMol * getMol() const
Get the molecule that this embedded fragmetn blongs to.
Definition: EmbeddedFrag.h:271
EmbeddedAtom(unsigned int aid, const RDGeom::Point2D &pos)
Definition: EmbeddedFrag.h:45
RDGeom::Point2D loc
Definition: EmbeddedFrag.h:118
RDGeom::Point2D normal
a zig-zag manner
Definition: EmbeddedFrag.h:124
double angle
the angle that is already takes at this atom, so any new atom attaching to
Definition: EmbeddedFrag.h:93
void Transform(const RDGeom::Transform2D &trans)
Definition: EmbeddedFrag.h:73
void Translate(const RDGeom::Point2D &shift)
Definition: EmbeddedFrag.h:301
ROMol is a molecule class that is intended to have a fixed topology.
Definition: ROMol.h:106
std::vector< INT_VECT > VECT_INT_VECT
Definition: types.h:204
std::map< int, Point2D > INT_POINT2D_MAP
Definition: point.h:511
boost::shared_array< double > DOUBLE_SMART_PTR
Definition: EmbeddedFrag.h:25
const INT_EATOM_MAP & GetEmbeddedAtoms() const
Definition: EmbeddedFrag.h:299
std::vector< int > INT_VECT
Definition: types.h:190
Class containing a fragment of a molecule that has already been embedded.
Definition: EmbeddedFrag.h:154
Includes a bunch of functionality for handling Atom and Bond queries.
Definition: Atom.h:29
Class that contains the data for an atoms that has alredy been embedded.
Definition: EmbeddedFrag.h:28
int Size() const
the number of atoms in the embedded system
Definition: EmbeddedFrag.h:317
int rotDir
useful only on atoms that are degree >= 4
Definition: EmbeddedFrag.h:116
int nbr2
the second neighbor of atom that from the &#39;angle&#39;
Definition: EmbeddedFrag.h:99
class for representing a bond
Definition: Bond.h:47
std::map< unsigned int, EmbeddedAtom > INT_EATOM_MAP
Definition: EmbeddedFrag.h:140
RDGeom::Point2D reflectPoint(const RDGeom::Point2D &point, const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2)
EmbeddedFrag()
Default constructor.
Definition: EmbeddedFrag.h:161
INT_EATOM_MAP::const_iterator INT_EATOM_MAP_CI
Definition: EmbeddedFrag.h:142
void Reflect(const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2)
Definition: EmbeddedFrag.h:80
void markDone()
Mark this fragment to be done for final embedding.
Definition: EmbeddedFrag.h:265
bool isDone()
If this fragment done for the final embedding.
Definition: EmbeddedFrag.h:268
#define PRECONDITION(expr, mess)
Definition: Invariant.h:107
EmbeddedAtom GetEmbeddedAtom(unsigned int aid) const
Definition: EmbeddedFrag.h:308
INT_EATOM_MAP::iterator INT_EATOM_MAP_I
Definition: EmbeddedFrag.h:141
RDKit::INT_VECT neighs
and these are the atom IDs of the neighbors that still need to be embedded
Definition: EmbeddedFrag.h:127
int CisTransNbr
is this is a cis/trans atom the neighbor of this atom that is involved in
Definition: EmbeddedFrag.h:104
void TransformPoint(Point2D &pt) const
int nbr1
the first neighbor of this atom that form the &#39;angle&#39;
Definition: EmbeddedFrag.h:96