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