source: src/bondgraph.hpp@ 111f4a

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

  /** 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 or others.
   *
   * \param &Set AtomSetMixin with all particles to consider
   */
  template <class container_type,
            class iterator_type,
            class const_iterator_type>
  double getMaxPossibleBondDistance(
      const AtomSetMixin<container_type,iterator_type,const_iterator_type> &Set) const
  {
    double max_distance = 0.;
    // get all elements
    std::set< const element *> PresentElements;
    for(const_iterator_type AtomRunner = Set.begin(); AtomRunner != Set.end(); ++AtomRunner) {
      PresentElements.insert( (*AtomRunner)->getType() );
    }
    // 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) {
        range<double> MinMaxDistance(0.,0.);
        getMinMaxDistance((*iter),(*otheriter), MinMaxDistance, IsAngstroem);
        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
   * \param &MinMaxDistance Range for interval on return
   * \param IsAngstroem whether units are in angstroem or bohr radii
   */
  void getMinMaxDistance(
      const BondedParticle * const Walker,
      const BondedParticle * const OtherWalker,
      range<double> &MinMaxDistance,
      bool IsAngstroem) 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
   * \param &MinMaxDistance Range for interval on return
   * \param IsAngstroem whether units are in angstroem or bohr radii
   */
  void getMinMaxDistanceSquared(
      const BondedParticle * const Walker,
      const BondedParticle * const OtherWalker,
      range<double> &MinMaxDistance,
      bool IsAngstroem) const;

  /** 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) const
  {
    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.");
      const double max_distance = getMaxPossibleBondDistance(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 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
    unsigned int bondcount = 0;
    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 << ".");
      bond * const Binder = new bond(Walker, OtherWalker, 1, -1);
      Walker->RegisterBond(WorldTime::getTime(), Binder);
      OtherWalker->RegisterBond(WorldTime::getTime(), Binder);
      bondcounter++;
    }
    LOG(1, "STATUS: "<< bondcounter << " bonds have been parsed.");
  }

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

  /** 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) {
      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 the BondLengthMatrix 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) 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 &MinMaxDistance Range for interval on return
   * \param IsAngstroem whether units are in angstroem or bohr radii
   */
  void getMinMaxDistance(
      const element * const Walker,
      const element * const OtherWalker,
      range<double> &MinMaxDistance,
      bool IsAngstroem) 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
   * @param MinMaxDistance reference to range type set on return
   * \param IsAngstroem whether units are in angstroem or bohr radii
   */
  void BondLengthMatrixMinMaxDistance(
      const element * const Walker,
      const element * const OtherWalker,
      range<double> &DistanceInterval,
      bool IsAngstroem) const;

  /** Returns bond criterion for given pair of elements based on covalent radius.
   * \param *Walker first element
   * \param *OtherWalker second element
   * @param MinMaxDistance reference to range type set on return
   * \param IsAngstroem whether units are in angstroem or bohr radii
   */
  void CovalentMinMaxDistance(
      const element * const Walker,
      const element * const OtherWalker,
      range<double> &DistanceInterval,
      bool IsAngstroem) 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) {
      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;
  //!> distance units are angstroem (true), bohr radii (false)
  bool IsAngstroem;
};

#endif /* BONDGRAPH_HPP_ */