source: src/bondgraph.hpp@ 3738f0

/*
 * 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 "AtomSet.hpp"
#include "bond.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Range.hpp"
#include "CodePatterns/Verbose.hpp"
#include "element.hpp"
#include "linkedcell.hpp"
#include "IPointCloud.hpp"
#include "PointCloudAdaptor.hpp"

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

class molecule;
class BondedParticle;
class MatrixContainer;

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

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


class BondGraph {
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);

  /** Returns the entry for a given index pair.
   * \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);

  /** Determines the maximum of all element::CovalentRadius for elements present in \a &Set.
   *
   * Sets BondGraph::max_distance
   *
   * \param &Set PointCloud with all particles
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  double SetMaxDistanceToMaxOfCovalentRadii(const AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set)
  {
    max_distance = 0.;

    for(const_iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
      const double radius = (*AtomRunner)->getType()->getCovalentRadius();
      if (radius > max_distance)
        max_distance = radius;
    }
    max_distance *= 2.;
    max_distance += BondThreshold;

    return max_distance;
  }

  /** Returns the upper limit on possible bond distances.
   *
   * Is set by SetMaxDistanceToMaxOfCovalentRadii(), is needed e.g. as cutoff
   * for linked-cell.
   *
   * \return maximum possible bond distance
   */
  double getMaxDistance() 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
   * \param &MinDistance lower bond bound on return
   * \param &MaxDistance upper bond bound on return
   * \param IsAngstroem whether units are in angstroem or bohr radii
   */
  void getMinMaxDistance(const BondedParticle * const Walker, const BondedParticle * const OtherWalker, double &MinDistance, double &MaxDistance, bool IsAngstroem);

  /** Creates the adjacency list for a given \a Range of iterable atoms.
   *
   * @param Set Range with begin and end iterator
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  void CreateAdjacency(AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set)
  {
    LOG(1, "STATUS: Removing all present bonds.");
    cleanAdjacencyList(Set);

    // 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.");
      SetMaxDistanceToMaxOfCovalentRadii(Set);

      LOG(1, "STATUS: Creating LinkedCell structure for given " << counter << " atoms.");
      PointCloudAdaptor< AtomSetMixin<container_type,iterator_type> > cloud(&Set, "SetOfAtoms");
      LinkedCell *LC = new LinkedCell(cloud, max_distance);

      CreateAdjacency(*LC);
      delete (LC);

      // 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(2, "STATUS: Printing list of created bonds.");
      std::stringstream output;
      for(const_iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
        (*AtomRunner)->OutputBondOfAtom(output);
        output << std::endl << "\t\t";
      }
      LOG(2, output.str());
    } else {
      LOG(1, "REJECT: AtomCount is " << counter << ", thus no bonds, no connections.");
    }
  }

  /** 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
   *  -# Count the atoms in the molecule with CountAtoms()
   *  -# partition cell into smaller linked cells of size \a bonddistance
   *  -# put each atom into its corresponding cell
   *  -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
   *  -# correct the bond degree iteratively (single->double->triple bond)
   *  -# finally print the bond list to \a *out if desired
   * \param &LC Linked Cell Container with all atoms
   */
  void CreateAdjacency(LinkedCell &LC);

  /** 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)
  {
    // remove every bond from the list
    for(iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
      BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
      for(BondList::iterator BondRunner = ListOfBonds.begin();
          !ListOfBonds.empty();
          BondRunner = ListOfBonds.begin()) {
        ASSERT((*BondRunner)->Contains(*AtomRunner),
            "BondGraph::cleanAdjacencyList() - "+
            toString(*BondRunner)+" does not contain "+
            toString(*AtomRunner)+".");
        delete((*BondRunner));
      }
    }
  }

  /** 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);
  }

private:
  static const double BondThreshold;

  /** Returns bond criterion for given pair 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 BondedParticle
   * \param *OtherWalker second BondedParticle
   * \param &MinDistance lower bond bound on return
   * \param &MaxDistance upper bond bound on return
   * \param IsAngstroem whether units are in angstroem or bohr radii
   */
  void BondLengthMatrixMinMaxDistance(const BondedParticle * const Walker, const BondedParticle * const OtherWalker, double &MinDistance, double &MaxDistance, bool IsAngstroem);

  /** Returns bond criterion for given pair based on covalent radius.
   * \param *Walker first BondedParticle
   * \param *OtherWalker second BondedParticle
   * \param &MinDistance lower bond bound on return
   * \param &MaxDistance upper bond bound on return
   * \param IsAngstroem whether units are in angstroem or bohr radii
   */
  void CovalentMinMaxDistance(const BondedParticle * const Walker, const BondedParticle * const OtherWalker, double &MinDistance, double &MaxDistance, bool IsAngstroem);


  /** 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) {
      BondList &ListOfBonds = (*AtomRunner)->getListOfBonds();
      for (BondList::iterator BondIter = ListOfBonds.begin();
          BondIter != ListOfBonds.end();
          ++BondIter)
        (*BondIter)->BondDegree = 1;
    }
  }

  /** 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
  {
    //DoLog(1) && (Log() << Verbose(1) << "Correcting Bond degree of each bond ... " << endl);
    int No = 0, OldNo = -1;
    do {
      OldNo = No;
      No=0;
      for(iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
        No+=(*AtomRunner)->CorrectBondDegree();
      }
    } while (OldNo != No);
    //DoLog(0) && (Log() << Verbose(0) << " done." << endl);
    return No;
  }

  //!> Matrix with bond lenth per two elements
  MatrixContainer *BondLengthMatrix;
  //!> maximum distance over all bonds possible
  double max_distance;
  //!> distance units are angstroem (true), bohr radii (false)
  bool IsAngstroem;
};

#endif /* BONDGRAPH_HPP_ */