rdkit/Code/GraphMol/Depictor/DepictUtils.h

//
//  Copyright (C) 2003-2022 Greg Landrum and other RDKit contributors
//
//   @@ All Rights Reserved @@
//  This file is part of the RDKit.
//  The contents are covered by the terms of the BSD license
//  which is included in the file license.txt, found at the root
//  of the RDKit source tree.
//
#include <RDGeneral/export.h>
#ifndef RD_DEPICT_UTILS_H
#define RD_DEPICT_UTILS_H

// REVIEW: remove extra headers here
#include <RDGeneral/types.h>
#include <GraphMol/RDKitBase.h>
#include <GraphMol/RWMol.h>
#include <GraphMol/ROMol.h>
#include <GraphMol/Substruct/SubstructUtils.h>
#include <Geometry/Transform2D.h>
#include <Geometry/Transform3D.h>
#include <Geometry/point.h>
#include <queue>

namespace RDDepict {

RDKIT_DEPICTOR_EXPORT extern double BOND_LEN;
RDKIT_DEPICTOR_EXPORT extern double COLLISION_THRES;
RDKIT_DEPICTOR_EXPORT extern double BOND_THRES;
RDKIT_DEPICTOR_EXPORT extern double ANGLE_OPEN;
RDKIT_DEPICTOR_EXPORT extern unsigned int MAX_COLL_ITERS;
RDKIT_DEPICTOR_EXPORT extern double HETEROATOM_COLL_SCALE;
RDKIT_DEPICTOR_EXPORT extern unsigned int NUM_BONDS_FLIPS;

typedef std::vector<const RDGeom::Point2D *> VECT_C_POINT;

typedef std::pair<int, int> PAIR_I_I;
typedef std::vector<PAIR_I_I> VECT_PII;
struct RDKIT_DEPICTOR_EXPORT gtIIPair {
  bool operator()(const PAIR_I_I &pd1, const PAIR_I_I &pd2) const {
    return pd1.first > pd2.first;
  }
};

typedef std::priority_queue<PAIR_I_I, VECT_PII, gtIIPair> PR_QUEUE;

typedef std::pair<double, PAIR_I_I> PAIR_D_I_I;
typedef std::list<PAIR_D_I_I> LIST_PAIR_DII;

//! Some utility functions used in generating 2D coordinates

//! Embed a ring as a convex polygon in 2D
/*!
  The process here is very straightforward:

  We take the center of the ring to lie at the origin, so put the first
  point at the origin and then sweep
  anti-clockwise by an angle A = 360/n for the next point.

  The length of the arm (l) we want to sweep is easy to compute given the
  bond length (b) we want to use for each bond in the ring (for now
  we will assume that this bond legnth is the same for all bonds in the ring:

  l = b/sqrt(2*(1 - cos(A))

  the above formula derives from the triangle formula, where side 'c' is given
  in terms of sides 'a' and 'b' as:

  c = a^2 + b^2 - 2.a.b.cos(A)

  where A is the angle between a and b
 */
RDKIT_DEPICTOR_EXPORT RDGeom::INT_POINT2D_MAP embedRing(
    const RDKit::INT_VECT &ring);

RDKIT_DEPICTOR_EXPORT void transformPoints(RDGeom::INT_POINT2D_MAP &nringCor,
                                           const RDGeom::Transform2D &trans);

//! Find a point that bisects the angle at rcr
/*!
  The new point lies between nb1 and nb2. The line (rcr, newPt) bisects the
  angle 'ang' at rcr
*/
RDKIT_DEPICTOR_EXPORT RDGeom::Point2D computeBisectPoint(
    const RDGeom::Point2D &rcr, double ang, const RDGeom::Point2D &nb1,
    const RDGeom::Point2D &nb2);

//! Reflect a set of point through the line joining two point
/*!
  ARGUMENTS:
  \param coordMap       a map of <int, point2D> going from atom id to current
                        coordinates of the points that need to be reflected:
                        The coordinates are overwritten
  \param loc1           the first point of the line that is to be used as a
  mirror
  \param loc2           the second point of the line to be used as a mirror
 */
RDKIT_DEPICTOR_EXPORT void reflectPoints(RDGeom::INT_POINT2D_MAP &coordMap,
                                         const RDGeom::Point2D &loc1,
                                         const RDGeom::Point2D &loc2);

RDKIT_DEPICTOR_EXPORT RDGeom::Point2D reflectPoint(const RDGeom::Point2D &point,
                                                   const RDGeom::Point2D &loc1,
                                                   const RDGeom::Point2D &loc2);

//! Set the neighbors yet to added to aid such that the atoms with the most subs
/// fall on opposite sides
/*!
  Ok this needs some explanation
  - Let A, B, C, D be the substituent on the central atom X (given
    by atom index aid)
  - also let A be the atom that is already embedded
  - Now we want the order in which the remaining neighbors B,C,D are
    added to X such that the atoms with atom with largest number of
    substituent fall on opposite sides of X so as to minimize atom
    clashes later in the depiction

  E.g. let say we have the following situation
<pre>
          B
       |  |
       A--X--C
       |  |
        --D--
          |
</pre>
  In this case the number substituent of A, B, C, D are 3, 1, 1,
  4 respectively so want to A and D to go opposite sides and so that
  we draw
<pre>
          B
       |  |  |
       A--X--D--
       |  |  |
          C
</pre>
  And the correct ordering of the neighbors is B,D,C
*/
RDKIT_DEPICTOR_EXPORT RDKit::INT_VECT setNbrOrder(unsigned int aid,
                                                  const RDKit::INT_VECT &nbrs,
                                                  const RDKit::ROMol &mol);

//! \brief From a given set of fused rings find the "core" rings, i.e. the rings
//! that are left after iteratively removing rings that are fused with only one
//! other ring by one or two atoms
/*
  \param fusedRings   list of all the rings in the fused system
  \param coreRingsIds this is where the IDs of the core rings are written
    \param mol          the molecule of interest

  \return list of rings that represent the core
*/
RDKIT_DEPICTOR_EXPORT RDKit::VECT_INT_VECT findCoreRings(
    const RDKit::VECT_INT_VECT &fusedRings, RDKit::INT_VECT &coreRingsIds,
    const RDKit::ROMol &mol);

//! \brief From a given set of rings find the ring the largest common elements
/// with other rings
/*
  Bit of a weird function - this is typically called once we have embedded some
  of the rings in a fused system and we are looking for the ring that must be
  embedded (or merged) next. The heuristic used here is to pick the rings with
  the maximum number of atoms in common with the rings that are already
  embedded.

  \param doneRings    a vector of ring IDs that have been embedded already
  \param fusedRings   list of all the rings in the fused system \param nextId
  this is where the ID for the next ring is written

  \return list of atom ids that are common
*/
RDKIT_DEPICTOR_EXPORT RDKit::INT_VECT findNextRingToEmbed(
    const RDKit::INT_VECT &doneRings, const RDKit::VECT_INT_VECT &fusedRings,
    int &nextId);

typedef std::pair<int, int> INT_PAIR;
typedef std::vector<INT_PAIR> INT_PAIR_VECT;
typedef INT_PAIR_VECT::const_iterator INT_PAIR_VECT_CI;

typedef std::pair<double, INT_PAIR> DOUBLE_INT_PAIR;

//! Sort a list of atoms by their  CIP rank
/*!
  \param mol        molecule of interest
  \param commAtms   atoms that need to be ranked
  \param ascending  sort to an ascending order or a descending order
*/
template <class T>
RDKIT_DEPICTOR_EXPORT T rankAtomsByRank(const RDKit::ROMol &mol,
                                        const T &commAtms,
                                        bool ascending = true);

//! computes a subangle for an atom of given hybridization and degree
/*!
  \param degree the degree of the atom (number of neighbors)
  \param htype  the atom's hybridization

  \return the subangle (in radians)
*/
inline double computeSubAngle(unsigned int degree,
                              RDKit::Atom::HybridizationType htype) {
  double angle = M_PI;
  switch (htype) {
    case RDKit::Atom::UNSPECIFIED:
    case RDKit::Atom::SP3:
      if (degree == 4) {
        angle = M_PI / 2;
      } else {
        angle = 2 * M_PI / 3;
      }
      break;
    case RDKit::Atom::SP2:
      angle = 2 * M_PI / 3;
      break;
    default:
      angle = 2. * M_PI / degree;
  }
  return angle;
}

//! computes the rotation direction between two vectors
/*!

  Let:

     v1 = loc1 - center

     v2 = loc2 - center

  If remaining angle(v1, v2) is < 180 and corss(v1, v2) > 0.0 then the rotation
  dir is +1.0

  else if remAngle(v1, v2) is > 180 and cross(v1, v2) < 0.0 then rotation dir is
  -1.0

  else if remAngle(v1, v2) is < 180 and cross(v1, v2) < 0.0 then rotation dir is
  -1.0

  finally if remAngle(v1, v2) is > 180 and cross(v1, v2) < 0.0 then rotation dir
  is +1.0

  \param center     the common point
  \param loc1       endpoint 1
  \param loc2       endpoint 2
  \param remAngle   the remaining angle about center in radians

  \return the rotation direction (1 or -1)
*/
inline int rotationDir(const RDGeom::Point2D &center,
                       const RDGeom::Point2D &loc1, const RDGeom::Point2D &loc2,
                       double remAngle) {
  auto pt1 = loc1 - center;
  auto pt2 = loc2 - center;
  auto cross = pt1.x * pt2.y - pt1.y * pt2.x;
  auto diffAngle = M_PI - remAngle;
  cross *= diffAngle;
  if (cross >= 0.0) {
    return -1;
  } else {
    return 1;
  }
}

//! computes and return the normal of a vector between two points
/*!
  \param center     the common point
  \param other      the endpoint

  \return the normal
*/
inline RDGeom::Point2D computeNormal(const RDGeom::Point2D &center,
                                     const RDGeom::Point2D &other) {
  auto res = other - center;
  res.normalize();
  std::swap(res.x, res.y);
  res.x *= -1;
  return res;
}

//! computes the rotation angle between two vectors
/*!
  \param center     the common point
  \param loc1       endpoint 1
  \param loc2       endpoint 2

  \return the angle (in radians)
*/
inline double computeAngle(const RDGeom::Point2D &center,
                           const RDGeom::Point2D &loc1,
                           const RDGeom::Point2D &loc2) {
  auto v1 = loc1 - center;
  auto v2 = loc2 - center;
  return v1.angleTo(v2);
}

//! \brief pick the ring to embed first in a fused system
/*!
  \param mol        the molecule of interest
  \param fusedRings the collection of the molecule's fused rings

  \return the index of the ring with the least number of substitutions
*/
RDKIT_DEPICTOR_EXPORT int pickFirstRingToEmbed(
    const RDKit::ROMol &mol, const RDKit::VECT_INT_VECT &fusedRings);

//! \brief find the rotatable bonds on the shortest path between two atoms
//!   we will ignore ring atoms, and double bonds which are marked cis/trans
/*!
  <b>Note</b> that rotatable in this context doesn't connect to the
    standard chemical definition of a rotatable bond; we're just talking
    about bonds than can be flipped in order to clean up the depiction.

  \param mol   the molecule of interest
  \param aid1  index of the first atom
  \param aid2  index of the second atom

  \return a set of the indices of the rotatable bonds
*/
RDKIT_DEPICTOR_EXPORT RDKit::INT_VECT getRotatableBonds(const RDKit::ROMol &mol,
                                                        unsigned int aid1,
                                                        unsigned int aid2);

//! \brief find all the rotatable bonds in a molecule
//!   we will ignore ring atoms, and double bonds which are marked cis/trans
/*!
  <b>Note</b> that rotatable in this context doesn't connect to the
    standard chemical definition of a rotatable bond; we're just talking
    about bonds than can be flipped in order to clean up the depiction.

  \param mol   the molecule of interest

  \return a set of the indices of the rotatable bonds
*/
RDKIT_DEPICTOR_EXPORT RDKit::INT_VECT getAllRotatableBonds(
    const RDKit::ROMol &mol);

//! \brief check if an atom is a spiro center (belongs to exactly 2 rings)
/*!
  \param aid  index of the atom
  \param mol  the molecule of interest

  \return true if the atom is a spiro center
*/
RDKIT_DEPICTOR_EXPORT bool isSpiroCenter(unsigned int aid,
                                          const RDKit::ROMol *mol);

//! \brief find spiro centers on the shortest path between two atoms
/*!
  \param mol   the molecule of interest
  \param aid1  index of the first atom
  \param aid2  index of the second atom

  \return a set of the indices of spiro centers
*/
RDKIT_DEPICTOR_EXPORT RDKit::INT_VECT getSpiroCenters(
    const RDKit::ROMol &mol, unsigned int aid1, unsigned int aid2);

//! Get the ids of the atoms and bonds that are connected to aid
RDKIT_DEPICTOR_EXPORT void getNbrAtomAndBondIds(unsigned int aid,
                                                const RDKit::ROMol *mol,
                                                RDKit::INT_VECT &aids,
                                                RDKit::INT_VECT &bids);

//! Find pairs of bonds that can be permuted at a non-ring degree 4 atom
/*!
  This function will return only those pairs that cannot be
  permuted by flipping a rotatble bond

         D
         |
         b3
         |
    A-b1-B-b2-C
         |
         b4
         |
         E
   For example in the above situation on the pairs (b1, b3) and (b1, b4) will be
  returned
  // All other permutations can be achieved via a rotatable bond flip.

   ARGUMENTS:
   \param center - location of the central atom
   \param nbrBids - a vector (of length 4) containing the ids of the bonds to
  the neighbors
   \param nbrLocs - locations of the neighbors
*/
RDKIT_DEPICTOR_EXPORT INT_PAIR_VECT findBondsPairsToPermuteDeg4(
    const RDGeom::Point2D &center, const RDKit::INT_VECT &nbrBids,
    const VECT_C_POINT &nbrLocs);

//! returns the rank of the atom for determining draw order
inline int getAtomDepictRank(const RDKit::Atom *at) {
  const int maxAtNum = 1000;
  const int maxDeg = 100;
  int anum = at->getAtomicNum();
  anum = anum == 1 ? maxAtNum : anum;  // favor non-hydrogen atoms
  int deg = at->getDegree();
  return maxDeg * anum + deg;
}

RDKIT_DEPICTOR_EXPORT bool hasTerminalRGroupOrQueryHydrogen(
    const RDKit::ROMol &query);
RDKIT_DEPICTOR_EXPORT std::unique_ptr<RDKit::RWMol> prepareTemplateForRGroups(
    RDKit::RWMol &templateMol);
RDKIT_DEPICTOR_EXPORT void reducedToFullMatches(
    const RDKit::RWMol &reducedQuery, const RDKit::RWMol &molHs,
    std::vector<RDKit::MatchVectType> &matches);
RDKIT_DEPICTOR_EXPORT bool invertWedgingIfMolHasFlipped(
    RDKit::ROMol &mol, const RDGeom::Transform3D &trans);
}  // namespace RDDepict

#endif