source: src/Graph/BondGraph.hpp@ e4b9d5

/*
 * bondgraph.hpp
 *
 *  Created on: Oct 29, 2009
 *      Author: heber
 */

#ifndef BONDGRAPH_HPP_
#define BONDGRAPH_HPP_

using namespace std;

/*********************************************** includes ***********************************/

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

#include <iosfwd>
#include <set>

#include <boost/serialization/array.hpp>

#include "Atom/AtomSet.hpp"
#include "Bond/bond.hpp"
#include "Box.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Range.hpp"
#include "Element/element.hpp"
#include "Fragmentation/MatrixContainer.hpp"
#include "Helpers/defs.hpp"
#include "LinkedCell/LinkedCell_View.hpp"
#include "LinkedCell/IPointCloud.hpp"
#include "LinkedCell/PointCloudAdaptor.hpp"
#include "WorldTime.hpp"

/****************************************** forward declarations *****************************/

class molecule;
class BondedParticle;
class MatrixContainer;

/********************************************** definitions *********************************/

/********************************************** declarations *******************************/


class BondGraph {
  //!> analysis bonds unit test should be friend to access private parts.
  friend class AnalysisBondsTest;
  //!> own bond graph unit test should be friend to access private parts.
  friend class BondGraphTest;
public:
  /** Constructor of class BondGraph.
   * This classes contains typical bond lengths and thus may be used to construct a bond graph for a given molecule.
   */
  BondGraph(bool IsA);

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

  /** Parses the bond lengths in a given file and puts them int a matrix form.
   * Allocates \a MatrixContainer for BondGraph::BondLengthMatrix, using MatrixContainer::ParseMatrix(),
   * but only if parsing is successful. Otherwise variable is left as NULL.
   * \param &input input stream to parse table from
   * \return true - success in parsing file, false - failed to parse the file
   */
  bool LoadBondLengthTable(std::istream &input);

  /** States whether the bond table has been loaded.
   * \return true - the internal matrix container is allocated, false - no bond table structure allocated
   */
  bool IsBondLengthTableLoaded() const { return BondLengthMatrix != NULL; }

  /** Removes allocated bond length table.
   *
   */
  void CleanupBondLengthTable();

  /** Internal helper to convert a set of atomicNumber_t to element refs.
   *
   * @param Set set of atomicNumber_t
   * @return set of element refs
   */
  std::set< const element *> getElementSetFromNumbers(const std::set<atomicNumber_t> &Set) const;

  /** Determines the maximum of all element::CovalentRadius for elements present in \a &Set.
   *
   * I.e. the function returns a sensible cutoff criteria for bond recognition,
   * e.g. to be used for LinkedCell_deprecated or others.
   *
   * \param &PresentElements set of elements whose maximal pair to find
   */
  double getMaxPossibleBondDistance(
      const std::set< const element *> &PresentElements) const
  {
    double max_distance = 0.;

    // create all element combinations
    for (std::set< const element *>::const_iterator iter = PresentElements.begin();
        iter != PresentElements.end();
        ++iter) {
      for (std::set< const element *>::const_iterator otheriter = iter;
          otheriter != PresentElements.end();
          ++otheriter) {
        const range<double> MinMaxDistance(getMinMaxDistance((*iter),(*otheriter)));
        if (MinMaxDistance.last > max_distance)
          max_distance = MinMaxDistance.last;
      }
    }
    return max_distance;
  }

  /** Determines the maximum of all element::CovalentRadius for elements present in \a &Set.
   *
   * I.e. the function returns a sensible cutoff criteria for bond recognition,
   * e.g. to be used for LinkedCell_deprecated or others.
   *
   * \param Walker element first element in the pair
   * \param &PresentElements set of elements whose maximal pair to find
   */
  double getMaxPossibleBondDistance(
      const element * const Walker,
      const std::set< const element *> &PresentElements) const
  {
    double max_distance = 0.;

    // create all element combinations
    for (std::set< const element *>::const_iterator iter = PresentElements.begin();
        iter != PresentElements.end();
        ++iter) {
      const range<double> MinMaxDistance(getMinMaxDistance((*iter),Walker));
      if (MinMaxDistance.last > max_distance)
        max_distance = MinMaxDistance.last;
    }
    return max_distance;
  }

  /** Returns bond criterion for given pair based on a bond length matrix.
   * This calls element-version of getMinMaxDistance().
   * \param *Walker first BondedParticle
   * \param *OtherWalker second BondedParticle
   * \return Range with bond interval
   */
  range<double> getMinMaxDistance(
      const BondedParticle * const Walker,
      const BondedParticle * const OtherWalker) const;

  /** Returns SQUARED bond criterion for given pair based on a bond length matrix.
   * This calls element-version of getMinMaxDistance() and squares the values
   * of either interval end.
   * \param *Walker first BondedParticle
   * \param *OtherWalker second BondedParticle
   * \return Range with bond interval
   */
  range<double> getMinMaxDistanceSquared(
      const BondedParticle * const Walker,
      const BondedParticle * const OtherWalker) const;

  /** Creates an adjacency list of the molecule.
   * Generally, we use the CSD approach to bond recognition, that is the the distance
   * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
   * a threshold t = 0.4 Angstroem.
   * To make it O(N log N) the function uses the linked-cell technique as follows:
   * The procedure is step-wise:
   *  -# Remove every bond in list
   *  -# go through every atom in given \a set, check the atoms therein against all possible bond partners, add bond if true
   *  -# correct the bond degree iteratively (single->double->triple bond)
   *  -# finally print the bond list to \a *out if desired
   * \param &set Container with all atoms to create adjacency for
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  void CreateAdjacency(
      AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
  {
    LOG(1, "STATUS: Removing all present bonds.");
    cleanAdjacencyList(Set);

    // gather set of all present elements
    std::set<atomicNumber_t> elements;
    for(typename AtomSetMixin<container_type,iterator_type,const_iterator_type>::iterator iter = Set.begin();
        iter != Set.end(); ++iter) {
      const atom * const Walker = dynamic_cast<const atom *>(*iter);
      ASSERT(Walker != NULL,
          "BondGraph::CreateAdjacency() - TesselPoint "
          +(*iter)->getName()+" that was not an atom retrieved from given set");
      elements.insert( Walker->getElementNo() );
    }
    // get all elements
    const std::set< const element *> PresentElements = getElementSetFromNumbers(elements);

    // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
    const unsigned int counter = Set.size();
    if (counter > 1) {
      LOG(1, "STATUS: Setting max bond distance.");
      LOG(1, "STATUS: Creating LinkedCell structure for given " << counter << " atoms.");
      LinkedCell::LinkedCell_View LC = getLinkedCell(getMaxPossibleBondDistance(PresentElements));

      LOG(1, "STATUS: Evaluating distance criterion for each atom.");

      const Box &domain = getDomain();
      const unsigned int CurrentTime = getTime();

      unsigned int BondCount = 0;
      // go through every atom in the set (observed cause we change its bonds)
      for(typename AtomSetMixin<container_type,iterator_type,const_iterator_type>::iterator iter = Set.begin();
          iter != Set.end(); ++iter) {
        const atom * const Walker = dynamic_cast<const atom *>(*iter);
        ASSERT(Walker != NULL,
            "BondGraph::CreateAdjacency() - TesselPoint "
            +(*iter)->getName()+" that was not an atom retrieved from given set");
        LOG(2, "INFO: Current Atom is " << *Walker << ".");

        // obtain all possible neighbors
        LinkedCell::LinkedList ListOfNeighbors = LC.getAllNeighbors(
            getMaxPossibleBondDistance(Walker->getType(), PresentElements),
            Walker->getPosition());
        if (!ListOfNeighbors.empty()) {
          // we have some possible candidates, go through each
          for (LinkedCell::LinkedList::const_iterator neighboriter = ListOfNeighbors.begin();
              neighboriter != ListOfNeighbors.end();
              ++neighboriter) {
            const atom * const OtherWalker = dynamic_cast<const atom *>(*neighboriter);
            ASSERT(OtherWalker != NULL,
                "BondGraph::CreateAdjacency() - TesselPoint "
                +(*neighboriter)->getName()+" that was not an atom retrieved from LinkedList");
            if (OtherWalker->getId() > Walker->getId()) {  // just to not add bonds from both sides
              LOG(4, "INFO: Current other atom is " << *OtherWalker << ".");

              const range<double> MinMaxDistanceSquared(
                  getMinMaxDistanceSquared(Walker, OtherWalker));
              const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
              LOG(3, "INFO: Checking squared distance " << distance << " against typical bond length of " << MinMaxDistanceSquared << ".");
              const bool status = MinMaxDistanceSquared.isInRange(distance);
              if (status) { // create bond if distance is smaller
                LOG(2, "ACCEPT: Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << ".");
                // directly use iter to avoid const_cast'ing Walker, too
                //const bond::ptr Binder =
                const_cast<atom *>(Walker)->addBond(CurrentTime, const_cast<atom *>(OtherWalker));
                ++BondCount;
              } else {
                LOG(3, "REJECT: Squared distance "
                    << distance << " is out of squared covalent bounds "
                    << MinMaxDistanceSquared << ".");
              }
            } else {
              LOG(4, "REJECT: Not Adding: Wrong order.");
            }
          }
        }
      }
      LOG(1, "I detected " << BondCount << " bonds in the molecule.");

      // correct bond degree by comparing valence and bond degree
      LOG(1, "STATUS: Correcting bond degree.");
      CorrectBondDegree(Set);

      // output bonds for debugging (if bond chain list was correctly installed)
      LOG(3, "STATUS: Printing list of created bonds.");
      if (DoLog(3)) {
        std::stringstream output;
        for(const_iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
          (*AtomRunner)->OutputBondOfAtom(output);
          output << std::endl << "\t\t";
        }
        LOG(3, output.str());
      }
    } else {
      LOG(1, "REJECT: AtomCount is " << counter << ", thus no bonds, no connections.");
    }
  }

  /** Creates an adjacency list of the given \a Set of atoms.
   *
   * Note that the input stream is required to refer to the same number of
   * atoms also contained in \a Set.
   *
   * \param &Set container with atoms
   * \param *input input stream to parse
 * \param skiplines how many header lines to skip
 * \param id_offset is base id compared to World startin at 0
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  void CreateAdjacencyListFromDbondFile(
      AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set,
      ifstream *input,
      unsigned int skiplines,
      int id_offset) const
  {
    char line[MAXSTRINGSIZE];

    // check input stream
    if (input->fail()) {
      ELOG(0, "Opening of bond file failed \n");
      return;
    };
    // skip headers
    for (unsigned int i=0;i<skiplines;i++)
      input->getline(line,MAXSTRINGSIZE);

    // create lookup map
    LOG(1, "STATUS: Creating lookup map.");
    std::map< unsigned int, atom *> AtomLookup;
    unsigned int counter = id_offset; // if ids do not start at 0
    for (iterator_type iter = Set.begin(); iter != Set.end(); ++iter) {
      AtomLookup.insert( make_pair( counter++, *iter) );
    }
    LOG(2, "INFO: There are " << counter << " atoms in the given set.");

    LOG(1, "STATUS: Scanning file.");
    unsigned int atom1, atom2;
    unsigned int bondcounter = 0;
    while (!input->eof()) // Check whether we read everything already
    {
      input->getline(line,MAXSTRINGSIZE);
      stringstream zeile(line);
      if (zeile.str().empty())
        continue;
      zeile >> atom1;
      zeile >> atom2;

      LOG(4, "INFO: Looking for atoms " << atom1 << " and " << atom2 << ".");
      if (atom2 < atom1) //Sort indices of atoms in order
        std::swap(atom1, atom2);
      ASSERT(atom2 < counter,
          "BondGraph::CreateAdjacencyListFromDbondFile() - ID "
          +toString(atom2)+" exceeds number of present atoms "+toString(counter)+".");
      ASSERT(AtomLookup.count(atom1),
          "BondGraph::CreateAdjacencyListFromDbondFile() - Could not find an atom with the ID given in dbond file");
      ASSERT(AtomLookup.count(atom2),
          "BondGraph::CreateAdjacencyListFromDbondFile() - Could not find an atom with the ID given in dbond file");
      atom * const Walker = AtomLookup[atom1];
      atom * const OtherWalker = AtomLookup[atom2];

      LOG(3, "INFO: Creating bond between atoms " << atom1 << " and " << atom2 << ".");
      //const bond::ptr Binder =
          Walker->addBond(WorldTime::getTime(), OtherWalker);
      bondcounter++;
    }
    LOG(1, "STATUS: "<< bondcounter << " bonds have been parsed.");
  }

  /** Removes all bonds within the given set of iterable atoms.
   *
   * @param Set Range with atoms
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  void cleanAdjacencyList(
      AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
  {
    // remove every bond from the list
    for(iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
      (*AtomRunner)->removeAllBonds(WorldTime::getTime());
    }
  }

  /** Checks whether the bond degree for each atom on the set matches with its valency.
   *
   * @param Set Range with atoms
   * @return true - bond degrees match valency, false - bond degree correction is needed
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  bool checkBondDegree(
      AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
  {
    std::set<atom *> allatoms;
    for(iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner)
      allatoms.insert(*AtomRunner);
    return checkBondDegree(allatoms);
  }

  /** correct bond degree by comparing valence and bond degree.
   * correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
   * iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
   * preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
   * double bonds as was expected.
   * @param Set Range with atoms
   * \return number of bonds that could not be corrected
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  int CorrectBondDegree(
      AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
  {
    // reset
    resetBondDegree(Set);
    // re-calculate
    return calculateBondDegree(Set);
  }

  /** Equality comparator for class BondGraph.
   *
   * @param other other instance to compare to
   * @return true - if equal in every member variable, except static
   *         \a BondGraph::BondThreshold.
   */
  bool operator==(const BondGraph &other) const;

  /** Unequality comparator for class BondGraph.
   *
   * @param other other instance to compare to
   * @return false - if equal in every member variable, except static
   *         \a BondGraph::BondThreshold.
   */
  bool operator!=(const BondGraph &other) const {
    return !(*this == other);
  }

  /** Returns the BondLengthMatrix entry for a given index pair.
   *
   * \note This will return a bond length of 1. if we cannot find an entry.
   *
   * \param firstelement index/atom number of first element (row index)
   * \param secondelement index/atom number of second element (column index)
   * \note matrix is of course symmetric.
   */
  double GetBondLength(
      int firstelement,
      int secondelement) const;

private:
  /** Convenience function to place access to World::getLinkedCell() into source module.
   *
   * @return ref to LinkedCell_View
   */
  LinkedCell::LinkedCell_View getLinkedCell(const double max_distance) const;

  /** Convenience function to place access to World::getDomain() into source module.
   *
   * @return ref to Box
   */
  Box &getDomain() const;

  /** Convenience function to place access to WorldTime::getTime() into source module.
   *
   * @return current time step
   */
  unsigned int getTime() const;

  /** Returns bond criterion for given pair based on a bond length matrix.
   * This calls either the covalent or the bond matrix criterion.
   * \param *Walker first BondedParticle
   * \param *OtherWalker second BondedParticle
   * \return Range with bond interval
   */
  range<double> getMinMaxDistance(
      const element * const Walker,
      const element * const OtherWalker) const;

  /** Returns bond criterion for given pair of elements based on a bond length matrix.
   * The matrix should be contained in \a this BondGraph and contain an element-
   * to-element length.
   * \param *Walker first element
   * \param *OtherWalker second element
   * \return Range with bond interval
   */
  range<double> BondLengthMatrixMinMaxDistance(
      const element * const Walker,
      const element * const OtherWalker) const;

  /** Returns bond criterion for given pair of elements based on covalent radius.
   * \param *Walker first element
   * \param *OtherWalker second element
   * \return Range with bond interval
   */
  range<double> CovalentMinMaxDistance(
      const element * const Walker,
      const element * const OtherWalker) const;


  /** Resets the bond::BondDegree of all atoms in the set to 1.
   *
   * @param Set Range with atoms
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  void resetBondDegree (
      AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
  {
    // reset bond degrees
    for(iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
      resetBondDegree(*AtomRunner);
    }
  }

  bool checkBondDegree(const std::set<atom *> &allatoms) const;

  /** Calculates the bond degree for each atom on the set.
   *
   * @param Set Range with atoms
   * @return number of non-matching bonds
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  int calculateBondDegree(
      AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
  {
    std::set<atom *> allatoms;
    for(iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner)
      allatoms.insert(*AtomRunner);
    return calculateBondDegree(allatoms);
  }

  int calculateBondDegree(const std::set<atom *> &allatoms) const;

  bool operator==(const periodentafel &other) const;

  bool operator!=(const periodentafel &other) const {
    return !(*this == other);
  }

  std::set<bond::ptr> getBackEdges() const;
  std::set<bond::ptr> getTreeEdges() const;

  int CorrectBondDegree(atom *_atom, const std::set<bond::ptr>& skipedges) const;

  void resetBondDegree(atom *_atom) const;

private:
  // default constructor for serialization
  BondGraph();

  friend class boost::serialization::access;
  // serialization
  template<class Archive>
  void serialize(Archive & ar, const unsigned int version)
  {
    //ar & const_cast<double &>(BondThreshold);
    ar & BondLengthMatrix;
    ar & IsAngstroem;
  }

  //!> half width of the interval for allowed bond distances
  static const double BondThreshold;
  //!> Matrix with bond lenth per two elements
  MatrixContainer *BondLengthMatrix;
  //!> distance units are angstroem (true), bohr radii (false)
  bool IsAngstroem;
};

#endif /* BONDGRAPH_HPP_ */