source: src/Graph/BondGraph.cpp@ 908dc7

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2011 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

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

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

#include "CodePatterns/MemDebug.hpp"

#include <iostream>

#include "atom.hpp"
#include "Bond/bond.hpp"
#include "Graph/BondGraph.hpp"
#include "Box.hpp"
#include "Element/element.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Range.hpp"
#include "CodePatterns/Verbose.hpp"
#include "molecule.hpp"
#include "Element/periodentafel.hpp"
#include "Fragmentation/MatrixContainer.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "World.hpp"
#include "WorldTime.hpp"

const double BondGraph::BondThreshold = 0.4;   //!< CSD threshold in bond check which is the width of the interval whose center is the sum of the covalent radii

BondGraph::BondGraph() :
    BondLengthMatrix(NULL),
    IsAngstroem(true)
{}

BondGraph::BondGraph(bool IsA) :
    BondLengthMatrix(NULL),
    IsAngstroem(IsA)
{}

BondGraph::~BondGraph()
{
  CleanupBondLengthTable();
}

void BondGraph::CleanupBondLengthTable()
{
  if (BondLengthMatrix != NULL) {
    delete(BondLengthMatrix);
  }
}

bool BondGraph::LoadBondLengthTable(
    std::istream &input)
{
  Info FunctionInfo(__func__);
  bool status = true;
  MatrixContainer *TempContainer = NULL;

  // allocate MatrixContainer
  if (BondLengthMatrix != NULL) {
    LOG(1, "MatrixContainer for Bond length already present, removing.");
    delete(BondLengthMatrix);
    BondLengthMatrix = NULL;
  }
  TempContainer = new MatrixContainer;

  // parse in matrix
  if ((input.good()) && (status = TempContainer->ParseMatrix(input, 0, 1, 0))) {
    LOG(1, "Parsing bond length matrix successful.");
  } else {
    ELOG(1, "Parsing bond length matrix failed.");
    status = false;
  }

  if (status) // set to not NULL only if matrix was parsed
    BondLengthMatrix = TempContainer;
  else {
    BondLengthMatrix = NULL;
    delete(TempContainer);
  }
  return status;
}

double BondGraph::GetBondLength(
    int firstZ,
    int secondZ) const
{
  std::cout << "Request for length between " << firstZ << " and " << secondZ << ": ";
  if (BondLengthMatrix == NULL) {
    std::cout << "-1." << std::endl;
    return( -1. );
  } else {
    std::cout << BondLengthMatrix->Matrix[0][firstZ][secondZ] << std::endl;
    return (BondLengthMatrix->Matrix[0][firstZ][secondZ]);
  }
}

range<double> BondGraph::CovalentMinMaxDistance(
    const element * const Walker,
    const element * const OtherWalker) const
{
  range<double> MinMaxDistance(0.,0.);
  MinMaxDistance.first = OtherWalker->getCovalentRadius() + Walker->getCovalentRadius();
  MinMaxDistance.first *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
  MinMaxDistance.last = MinMaxDistance.first + BondThreshold;
  MinMaxDistance.first -= BondThreshold;
  return MinMaxDistance;
}

range<double> BondGraph::BondLengthMatrixMinMaxDistance(
    const element * const Walker,
    const element * const OtherWalker) const
{
  ASSERT(BondLengthMatrix, "BondGraph::BondLengthMatrixMinMaxDistance() called without NULL BondLengthMatrix.");
  ASSERT(Walker, "BondGraph::BondLengthMatrixMinMaxDistance() - illegal element given.");
  ASSERT(OtherWalker, "BondGraph::BondLengthMatrixMinMaxDistance() - illegal other element given.");
  range<double> MinMaxDistance(0.,0.);
  MinMaxDistance.first = GetBondLength(Walker->getAtomicNumber()-1, OtherWalker->getAtomicNumber()-1);
  MinMaxDistance.first *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
  MinMaxDistance.last = MinMaxDistance.first + BondThreshold;
  MinMaxDistance.first -= BondThreshold;
  return MinMaxDistance;
}

range<double> BondGraph::getMinMaxDistance(
    const element * const Walker,
    const element * const OtherWalker) const
{
  range<double> MinMaxDistance(0.,0.);
  if (BondLengthMatrix == NULL) {// safety measure if no matrix has been parsed yet
    LOG(2, "INFO: Using Covalent radii criterion for [min,max) distances.");
    MinMaxDistance = CovalentMinMaxDistance(Walker, OtherWalker);
  } else {
    LOG(2, "INFO: Using Covalent radii criterion for [min,max) distances.");
    MinMaxDistance = BondLengthMatrixMinMaxDistance(Walker, OtherWalker);
  }
  return MinMaxDistance;
}

range<double> BondGraph::getMinMaxDistance(
    const BondedParticle * const Walker,
    const BondedParticle * const OtherWalker) const
{
  return getMinMaxDistance(Walker->getType(), OtherWalker->getType());
}

range<double> BondGraph::getMinMaxDistanceSquared(
    const BondedParticle * const Walker,
    const BondedParticle * const OtherWalker) const
{
  // use non-squared version
  range<double> MinMaxDistance(getMinMaxDistance(Walker, OtherWalker));
  // and square
  MinMaxDistance.first *= MinMaxDistance.first;
  MinMaxDistance.last *= MinMaxDistance.last;
  return MinMaxDistance;
}

void BondGraph::CreateAdjacency(LinkedCell &LC) const
{
  atom *Walker = NULL;
  atom *OtherWalker = NULL;
  int n[NDIM];
  Box &domain = World::getInstance().getDomain();
  size_t CurrentTime = WorldTime::getTime();

  unsigned int BondCount = 0;
  // 3a. go through every cell
  LOG(3, "INFO: Celling ... ");
  for (LC.n[0] = 0; LC.n[0] < LC.N[0]; LC.n[0]++)
    for (LC.n[1] = 0; LC.n[1] < LC.N[1]; LC.n[1]++)
      for (LC.n[2] = 0; LC.n[2] < LC.N[2]; LC.n[2]++) {
        const TesselPointSTLList *List = LC.GetCurrentCell();
        LOG(2, "Current cell is " << LC.n[0] << ", " << LC.n[1] << ", " << LC.n[2] << " with No. " << LC.index << " containing " << List->size() << " points.");
        if (List != NULL) {
          for (TesselPointSTLList::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
            Walker = dynamic_cast<atom*>(*Runner);
            ASSERT(Walker != NULL,
                "BondGraph::CreateAdjacency() - Tesselpoint that was not an atom retrieved from LinkedNode");
            LOG(2, "INFO: Current Atom is " << *Walker << ".");
            // 3c. check for possible bond between each atom in this and every one in the 27 cells
            for (n[0] = -1; n[0] <= 1; n[0]++)
              for (n[1] = -1; n[1] <= 1; n[1]++)
                for (n[2] = -1; n[2] <= 1; n[2]++) {
                  const TesselPointSTLList *OtherList = LC.GetRelativeToCurrentCell(n);
                  if (OtherList != NULL) {
                    LOG(3, "INFO: Current relative cell is " << n[0] << ", " << n[1] << ", " << n[2] << " with No. " << LC.index << " containing " << List->size() << " points.");
                    for (TesselPointSTLList::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
                      if ((*OtherRunner) > Walker) {  // just to not add bonds from both sides
                        OtherWalker = dynamic_cast<atom*>(*OtherRunner);
                        ASSERT(OtherWalker != NULL,
                            "BondGraph::CreateAdjacency() - TesselPoint that was not an atom retrieved from LinkedNode");
                        LOG(3, "INFO: Current other atom is " << *OtherWalker << ".");
                        if (OtherWalker->father > Walker->father ) { // just to not add bonds from both sides
                          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(1, "ACCEPT: Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << ".");
                            //const bond * Binder =
                              Walker->father->addBond(CurrentTime, OtherWalker->father);
                            BondCount++;
                          } else {
                            LOG(2, "REJECT: Squared distance "
                                << distance << " is out of squared covalent bounds "
                                << MinMaxDistanceSquared << ".");
                          }
                        } else {
                          LOG(5, "REJECT: Not Adding: Wrong order of father's.");
                        }
                      } else {
                        LOG(5, "REJECT: Not Adding: Wrong order.");
                      }
                    }
                  }
                }
          }
        }
      }
    LOG(1, "I detected " << BondCount << " bonds in the molecule.");
}

bool BondGraph::operator==(const BondGraph &other) const
{
  if (IsAngstroem != other.IsAngstroem)
    return false;
  if (((BondLengthMatrix != NULL) && (other.BondLengthMatrix == NULL))
      || ((BondLengthMatrix == NULL) && (other.BondLengthMatrix != NULL)))
    return false;
  if ((BondLengthMatrix != NULL) && (other.BondLengthMatrix != NULL))
    if (*BondLengthMatrix != *other.BondLengthMatrix)
      return false;
  return true;
}