source: src/Fragmentation/Exporters/SaturatedFragment.hpp@ 20fc6f

/*
 * SaturatedFragment.hpp
 *
 *  Created on: Mar 3, 2013
 *      Author: heber
 */

#ifndef SATURATEDFRAGMENT_HPP_
#define SATURATEDFRAGMENT_HPP_

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iosfwd>
#include <map>
#include <set>
#include <string>

#include "Atom/atom_bondedparticleinfo.hpp"
#include "Bond/bond.hpp"
#include "Fragmentation/EdgesPerFragment.hpp"
#include "Fragmentation/KeySet.hpp"
#include "Fragmentation/HydrogenSaturation_enum.hpp"
#include "Parser/FormatParserStorage.hpp"

#include "LinearAlgebra/Vector.hpp"

class atom;
class HydrogenPool;

/** The SaturatedFragment class acts as a wrapper to a KeySet by adding a list
 * of saturation hydrogens.
 *
 * This SaturatedFragment along with a currently leased set of hydrogens from the
 * HydrogenPool is all that is required to create a fully storable molecular
 * fragment from a given Keyset.
 *
 * The instance notes down its existence in an external container.
 *
 */
class SaturatedFragment
{
public:
  //!> typedef to a container to mark keysets that are in use
  typedef std::set<KeySet> KeySetsInUse_t;

  //!> List of points giving saturation positions for a single bond neighbor
  typedef std::list<Vector> SaturationsPositions_t;
  //!> map for one atom, containing the saturation points for all its neighbors
  typedef std::map<int, SaturationsPositions_t> SaturationsPositionsPerNeighbor_t;
  //!> containing the saturation points over all desired atoms required
  typedef std::map<int, SaturationsPositionsPerNeighbor_t> GlobalSaturationPositions_t;

  /** Constructor of SaturatedFragment requires \a set which we are tightly
   * associated.
   *
   * \param _set KeySet which this instance is associated with
   * \param _container container to add KeySet as in-use
   * \param _hydrogens pool with hydrogens for saturation
   * \param _globalsaturationpositions saturation positions to be used
   */
  SaturatedFragment(
      const KeySet &_set,
      KeySetsInUse_t &_container,
      HydrogenPool &_hydrogens,
      const enum HydrogenTreatment _treatment,
      const enum HydrogenSaturation saturation,
      const GlobalSaturationPositions_t &_globalsaturationpositions);

  /** Destructor of class SaturatedFragment.
   *
   */
  ~SaturatedFragment();

  /** Getter for the KeySet this instance is associated with.
   *
   * \return const ref to KeySet
   */
  const KeySet & getKeySet() const
  {
    return set;
  }

  /** Getter for the FullMolecule this instance is associated with.
   *
   * \return const ref to FullMolecule
   */
  const KeySet & getFullMolecule() const
  {
    return FullMolecule;
  }

  /** Getter for the SaturationHydrogens this instance is associated with.
   *
   * \return const ref to SaturationHydrogens
   */
  const KeySet & getSaturationHydrogens() const
  {
    return SaturationHydrogens;
  }

  /** Determines the bounding box of the fragment without performing the hydrogen
   * saturation but by simply using the atoms whose bonds are cut.
   *
   * \return pair of Vector with min and max coordinates
   */
  std::pair<Vector, Vector> getRoughBoundingBox() const;

  /** Returns the edges of the bond graph of this fragment.
   *
   * \return set of edges
   */
  FragmentationEdges::edges_t getEdges() const;

  /** Prints the config of the fragment of \a _type to \a out.
   *
   * \param out output stream to write to
   * \param _type parser type to write config
   */
  bool OutputConfig(
      std::ostream &out,
      const ParserTypes _type) const;

  //!> typedef for the replaced atom (due to cut bonds) and their replacement hydrogens
  typedef std::multimap<atomId_t, atomId_t> replaced_atoms_t;

private:
  /** Helper function to lease and bring in place saturation hydrogens.
   *
   * Here, we use local information to calculate saturation positions.
   *
   */
  void saturate();

  /** Helper function to lease and bring in place saturation hydrogens.
   *
   * Here, saturation positions have to be calculated before and are fully
   * stored in \a _globalsaturationpositions.
   *
   * \param_globalsaturationpositions
   */
  void saturate(const GlobalSaturationPositions_t &_globalsaturationpositions);

  /** Replaces all cut bonds with respect to the given atom by hydrogens.
   *
   * \param _atom atom whose cut bonds to saturate
   * \param _cutbonds list of cut bonds for \a _atom
   * \return true - bonds saturated, false - something went wrong
   */
  bool saturateAtom(atom * const _atom, const BondList &_cutbonds);

  /** Helper function to get a hydrogen replacement for a given \a replacement.
   *
   * \param replacement atom to "replace" with hydrogen in a fragment.
   * \return ref to leased hydrogen atom
   */
  atom * const getHydrogenReplacement(atom * const replacement);

  /** Sets a saturation hydrogen at the given position in place of \a _father.
   *
   * \param _OwnerAtom atom "owning" the hydrogen (i.e. the one getting saturated)
   * \param _position new position relative to \a _OwnerAtom
   * \param _distance scale factor to the distance (default 1.)
   * \param _father bond partner of \a _OwnerAtom that is replaced
   */
  const atom& setHydrogenReplacement(
      const atom * const _OwnerAtom,
      const Vector &_position,
      const double _distance,
      atom * const _father);

  /** Leases and adds a Hydrogen atom in replacement for the given atom \a *partner in bond with a *origin.
   * Here, we have to distinguish between single, double or triple bonds as stated by \a BondDegree, that each demand
   * a different scheme when adding \a *replacement atom for the given one.
   * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one
   * -# Double Bond: Here, we need the **BondList of the \a *origin atom, by scanning for the other bonds instead of
   *    *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalvector().
   *    The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two
   *    replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the
   *    element of *origin and form the sin/cos admixture of both plane vectors for the new coordinates of the two
   *    hydrogens forming this angle with *origin.
   * -# Triple Bond: The idea is to set up a tetraoid (C1-H1-H2-H3) (however the lengths \f$b\f$ of the sides of the base
   *    triangle formed by the to be added hydrogens are not equal to the typical bond distance \f$l\f$ but have to be
   *    determined from the typical angle \f$\alpha\f$ for a hydrogen triple connected to the element of *origin):
   *    We have the height \f$d\f$ as the vector in *Bond direction (from triangle C1-H1-H2).
   *    \f[ h = l \cdot \cos{\left (\frac{\alpha}{2} \right )} \qquad b = 2l \cdot \sin{\left (\frac{\alpha}{2} \right)} \quad \rightarrow \quad d = l \cdot \sqrt{\cos^2{\left (\frac{\alpha}{2} \right)}-\frac{1}{3}\cdot\sin^2{\left (\frac{\alpha}{2}\right )}}
   *    \f]
   *    vector::GetNormalvector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalvector creates
   *    the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above.
   *    The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
   *    the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
   *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
   *    \f]
   *    as the coordination of all three atoms in the coordinate system of these three vectors:
   *    \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
   *
   * \param TopBond pointer to bond between \a *origin and \a *replacement
   * \param Origin atom that is actually contained in the fragment
   * \param Replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
   * \param isAngstroem whether the coordination of the given atoms is in AtomicLength (false) or Angstrom(true)
   * \return number of atoms added, if < bond::BondDegree then something went wrong
   */
  bool AddHydrogenReplacementAtom(
      bond::ptr TopBond,
      atom *Origin,
      atom *Replacement,
      bool IsAngstroem);

private:
  //!> container to mark ourselves RAII-style
  KeySetsInUse_t &container;
  //!> key set this fragment is associated with.
  const KeySet &set;
  //!> pool with saturation hydrogens
  HydrogenPool &hydrogens;
  //!> key set containing all atoms used for e.g. storing this to file
  KeySet FullMolecule;
  //!> key set containing the ids of all hydrogens added for saturation
  KeySet SaturationHydrogens;
  //!> whether hydrogens are generally contained in set or not
  const enum HydrogenTreatment treatment;
  //!> whether to actually saturate or not
  const enum HydrogenSaturation saturation;
  //!> list of saturated hydrogens to their original atoms
  replaced_atoms_t replaced_atoms;
};

#endif /* SATURATEDFRAGMENT_HPP_ */