source: src/molecule_graph.cpp@ c44322

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * molecule_graph.cpp
 *
 *  Created on: Oct 5, 2009
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include <stack>

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Box.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "config.hpp"
#include "Graph/DepthFirstSearchAnalysis.hpp"
#include "Element/element.hpp"
#include "Graph/BondGraph.hpp"
#include "Helpers/defs.hpp"
#include "Helpers/helpers.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "LinkedCell/linkedcell.hpp"
#include "LinkedCell/PointCloudAdaptor.hpp"
#include "molecule.hpp"
#include "World.hpp"
#include "WorldTime.hpp"


/** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule.
 * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL.
 * \param *reference reference molecule with the bond structure to be copied
 * \param **&ListOfLocalAtoms Lookup table for this subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled
 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 * \return true - success, false - failure
 */
bool molecule::FillBondStructureFromReference(const molecule * const reference, atom **&ListOfLocalAtoms, bool FreeList)
{
  bool status = true;

  LOG(1, "Begin of FillBondStructureFromReference.");
  // fill ListOfLocalAtoms if NULL was given
  if (!FillListOfLocalAtoms(ListOfLocalAtoms, reference->getAtomCount())) {
    LOG(1, "Filling of ListOfLocalAtoms failed.");
    return false;
  }

  if (status) {
    LOG(1, "Creating adjacency list for molecule " << getName() << ".");
    // remove every bond from the list
    for_each(begin(), end(),
        boost::bind(&BondedParticle::ClearBondsAtStep,_1,WorldTime::getTime()));


    for(molecule::iterator iter = begin(); iter != end(); ++iter) {
      const atom * const Father = (*iter)->GetTrueFather();
      //const int AtomNo = Father->getNr(); // global id of the current walker
      const BondList& ListOfBonds = Father->getListOfBonds();
      for (BondList::const_iterator Runner = ListOfBonds.begin();
          Runner != ListOfBonds.end();
          ++Runner) {
        atom * const OtherAtom = (*Runner)->GetOtherAtom((*iter)->GetTrueFather());
        atom * const OtherWalker = ListOfLocalAtoms[OtherAtom->getNr()]; // local copy of current bond partner of walker
        if (OtherWalker != NULL) {
          if (OtherWalker->getNr() > (*iter)->getNr())
            AddBond((*iter), OtherWalker, (*Runner)->BondDegree);
        } else {
          LOG(1, "OtherWalker = ListOfLocalAtoms[" << OtherAtom->getNr() << "] is NULL!");
          status = false;
        }
      }
    }
  }

  if ((FreeList) && (ListOfLocalAtoms != NULL)) {
    // free the index lookup list
    delete[](ListOfLocalAtoms);
  }
  LOG(1, "End of FillBondStructureFromReference.");
  return status;
};

/** Checks for presence of bonds within atom list.
 * TODO: more sophisticated check for bond structure (e.g. connected subgraph, ...)
 * \return true - bonds present, false - no bonds
 */
bool molecule::hasBondStructure() const
{
  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
    //LOG(0, "molecule::hasBondStructure() - checking bond list of atom " << (*AtomRunner)->getId() << ".");
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    if (!ListOfBonds.empty())
      return true;
  }
  return false;
}

/** Prints a list of all bonds to \a *out.
 */
void molecule::OutputBondsList() const
{
  if (DoLog(1)) {
    std::stringstream output;
    output << std::endl << "From contents of bond chain list:";
    for(molecule::const_iterator AtomRunner = molecule::begin(); AtomRunner != molecule::end(); ++AtomRunner) {
      const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
      for(BondList::const_iterator BondRunner = ListOfBonds.begin();
          BondRunner != ListOfBonds.end();
          ++BondRunner)
        if ((*BondRunner)->leftatom == *AtomRunner) {
          output << *(*BondRunner) << "\t";
        }
    }
    LOG(1, output.str());
  }
}


/** Storing the bond structure of a molecule to file.
 * Simply stores Atom::Nr and then the Atom::Nr of all bond partners per line.
 * \param &filename name of file
 * \param path path to file, defaults to empty
 * \return true - file written successfully, false - writing failed
 */
bool molecule::StoreAdjacencyToFile(std::string filename, std::string path)
{
  ofstream AdjacencyFile;
  string line;
  bool status = true;

  if (path != "")
    line = path + "/" + filename;
  else
    line = filename;
  AdjacencyFile.open(line.c_str(), ios::out);
  LOG(1, "Saving adjacency list ... ");
  if (AdjacencyFile.good()) {
    AdjacencyFile << "m\tn" << endl;
    for_each(begin(),end(),bind2nd(mem_fun(&atom::OutputAdjacency),&AdjacencyFile));
    AdjacencyFile.close();
    LOG(1, "\t... done.");
  } else {
    LOG(1, "\t... failed to open file " << line << ".");
    status = false;
  }

  return status;
}
;

/** Storing the bond structure of a molecule to file.
 * Simply stores Atom::Nr and then the Atom::Nr of all bond partners, one per line.
 * \param &filename name of file
 * \param path path to file, defaults to empty
 * \return true - file written successfully, false - writing failed
 */
bool molecule::StoreBondsToFile(std::string filename, std::string path)
{
  ofstream BondFile;
  string line;
  bool status = true;

  if (path != "")
    line = path + "/" + filename;
  else
    line = filename;
  BondFile.open(line.c_str(), ios::out);
  LOG(1, "Saving adjacency list ... ");
  if (BondFile.good()) {
    BondFile << "m\tn" << endl;
    for_each(begin(),end(),bind2nd(mem_fun(&atom::OutputBonds),&BondFile));
    BondFile.close();
    LOG(1, "\t... done.");
  } else {
    LOG(1, "\t... failed to open file " << line << ".");
    status = false;
  }

  return status;
}
;

/** Adds a bond as a copy to a given one
 * \param *left leftatom of new bond
 * \param *right rightatom of new bond
 * \param *CopyBond rest of fields in bond are copied from this
 * \return pointer to new bond
 */
bond * molecule::CopyBond(atom *left, atom *right, bond *CopyBond)
{
  bond *Binder = AddBond(left, right, CopyBond->BondDegree);
  Binder->Cyclic = CopyBond->Cyclic;
  Binder->Type = CopyBond->Type;
  return Binder;
}
;

/** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
 * \param **&ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
 * \param GlobalAtomCount number of atoms in the complete molecule
 * \return true - success, false - failure (ListOfLocalAtoms != NULL)
 */
bool molecule::FillListOfLocalAtoms(atom **&ListOfLocalAtoms, const int GlobalAtomCount)
{
  bool status = true;

  if (ListOfLocalAtoms == NULL) { // allocate and fill list of this fragment/subgraph
    status = status && CreateFatherLookupTable(ListOfLocalAtoms, GlobalAtomCount);
  } else
    return false;

  return status;
}


/** Creates a lookup table for true father's Atom::Nr -> atom ptr.
 * \param *start begin of list (STL iterator, i.e. first item)
 * \paran *end end of list (STL iterator, i.e. one past last item)
 * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
 * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
 * \return true - success, false - failure
 */
bool molecule::CreateFatherLookupTable(atom **&LookupTable, int count)
{
  bool status = true;
  int AtomNo;

  if (LookupTable != NULL) {
    ELOG(1, "Pointer for Lookup table is not NULL! Aborting ...");
    return false;
  }

  // count them
  if (count == 0) {
    for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::Nr -> atom) used as a marker table lateron
      count = (count < (*iter)->GetTrueFather()->getNr()) ? (*iter)->GetTrueFather()->getNr() : count;
    }
  }
  if (count <= 0) {
    ELOG(1, "Count of lookup list is 0 or less.");
    return false;
  }

  // allocate and fill
  LookupTable = new atom *[count+1];
  if (LookupTable == NULL) {
    ELOG(0, "LookupTable memory allocation failed!");
    performCriticalExit();
    status = false;
  } else {
    for (int i=0;i<=count;i++)
      LookupTable[i] = NULL;
    for (molecule::iterator iter = begin(); iter != end(); ++iter) {
      AtomNo = (*iter)->GetTrueFather()->getNr();
      if ((AtomNo >= 0) && (AtomNo <= count)) {
        LOG(3, "DEBUG: Setting LookupTable[" << AtomNo << "] to " << *(*iter));
        LookupTable[AtomNo] = (*iter);
      } else {
        ELOG(1, "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << "].");
        status = false;
        break;
      }
    }
  }

  return status;
};



/** Corrects the nuclei position if the fragment was created over the cell borders.
 * Scans all bonds, checks the distance, if greater than typical, we have a candidate for the correction.
 * We remove the bond whereafter the graph probably separates. Then, we translate the one component periodically
 * and re-add the bond. Looping on the distance check.
 * \param *out ofstream for debugging messages
 */
bool molecule::ScanForPeriodicCorrection()
{
  bond *Binder = NULL;
  //bond *OtherBinder = NULL;
  atom *Walker = NULL;
  atom *OtherWalker = NULL;
  RealSpaceMatrix matrix = World::getInstance().getDomain().getM();
  enum GraphEdge::Shading *ColorList = NULL;
  double tmp;
  //bool LastBond = true; // only needed to due list construct
  Vector Translationvector;
  //std::deque<atom *> *CompStack = NULL;
  std::deque<atom *> *AtomStack = new std::deque<atom *>; // (getAtomCount());
  bool flag = true;
  BondGraph *BG = World::getInstance().getBondGraph();

  LOG(2, "Begin of ScanForPeriodicCorrection.");

  ColorList = new enum GraphEdge::Shading[getAtomCount()];
  for (int i=0;i<getAtomCount();i++)
    ColorList[i] = (enum GraphEdge::Shading)0;
  if (flag) {
    // remove bonds that are beyond bonddistance
    Translationvector.Zero();
    // scan all bonds
    flag = false;
    for(molecule::iterator AtomRunner = begin(); (!flag) && (AtomRunner != end()); ++AtomRunner) {
      const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
      for(BondList::const_iterator BondRunner = ListOfBonds.begin();
          (!flag) && (BondRunner != ListOfBonds.end());
          ++BondRunner) {
        Binder = (*BondRunner);
        for (int i=NDIM;i--;) {
          tmp = fabs(Binder->leftatom->at(i) - Binder->rightatom->at(i));
          //LOG(3, "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << ".");
          const range<double> MinMaxDistance(
              BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom));
          if (!MinMaxDistance.isInRange(tmp)) {
            LOG(2, "Correcting at bond " << *Binder << ".");
            flag = true;
            break;
          }
        }
      }
    }
    //if (flag) {
    if (0) {
      // create translation vector from their periodically modified distance
      for (int i=NDIM;i--;) {
        tmp = Binder->leftatom->at(i) - Binder->rightatom->at(i);
        const range<double> MinMaxDistance(
            BG->getMinMaxDistance(Binder->leftatom, Binder->rightatom));
        if (fabs(tmp) > MinMaxDistance.last)  // check against Min is not useful for components
          Translationvector[i] = (tmp < 0) ? +1. : -1.;
      }
      Translationvector *= matrix;
      LOG(3, "INFO: Translation vector is " << Translationvector << ".");
      // apply to all atoms of first component via BFS
      for (int i=getAtomCount();i--;)
        ColorList[i] = GraphEdge::white;
      AtomStack->push_front(Binder->leftatom);
      while (!AtomStack->empty()) {
        Walker = AtomStack->front();
        AtomStack->pop_front();
        //LOG(3, "INFO: Current Walker is: " << *Walker << ".");
        ColorList[Walker->getNr()] = GraphEdge::black;    // mark as explored
        *Walker += Translationvector; // translate
        const BondList& ListOfBonds = Walker->getListOfBonds();
        for (BondList::const_iterator Runner = ListOfBonds.begin();
            Runner != ListOfBonds.end();
            ++Runner) {
          if ((*Runner) != Binder) {
            OtherWalker = (*Runner)->GetOtherAtom(Walker);
            if (ColorList[OtherWalker->getNr()] == GraphEdge::white) {
              AtomStack->push_front(OtherWalker); // push if yet unexplored
            }
          }
        }
      }
//      // re-add bond
//      if (OtherBinder == NULL) { // is the only bond?
//        //Do nothing
//      } else {
//        if (!LastBond) {
//          link(Binder, OtherBinder); // no more implemented bond::previous ...
//        } else {
//          link(OtherBinder, Binder); // no more implemented bond::previous ...
//        }
//      }
    } else {
      LOG(3, "No corrections for this fragment.");
    }
    //delete(CompStack);
  }
  // free allocated space from ReturnFullMatrixforSymmetric()
  delete(AtomStack);
  delete[](ColorList);
  LOG(2, "End of ScanForPeriodicCorrection.");

  return flag;
};