source: src/bondgraph.cpp@ 111f4a

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 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.hpp"
#include "bondgraph.hpp"
#include "Box.hpp"
#include "element.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Range.hpp"
#include "CodePatterns/Verbose.hpp"
#include "molecule.hpp"
#include "parser.hpp"
#include "periodentafel.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(bool IsA) :
    BondLengthMatrix(NULL),
    IsAngstroem(IsA)
{}

BondGraph::~BondGraph()
{
  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);
  }
  TempContainer = new MatrixContainer;

  // parse in matrix
  if ((input.good()) && (status = TempContainer->ParseMatrix(input, 0, 1, 0))) {
    LOG(1, "Parsing bond length matrix successful.");
  } else {
    DoeLog(1) && (eLog()<< Verbose(1) << "Parsing bond length matrix failed." << endl);
    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]);
  }
}

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

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

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

void BondGraph::getMinMaxDistance(
    const BondedParticle * const Walker,
    const BondedParticle * const OtherWalker,
    range<double> &MinMaxDistance,
    bool IsAngstroem) const
{
  getMinMaxDistance(Walker->getType(), OtherWalker->getType(),MinMaxDistance, IsAngstroem);
}

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

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

  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 " << LC.n[0] << ", " << LC.n[1] << ", " << LC.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");
                        range<double> MinMaxDistanceSquared(0.,0.);
                        getMinMaxDistanceSquared(Walker, OtherWalker, MinMaxDistanceSquared, IsAngstroem);
                        const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
                        LOG(2, "INFO: Checking squared distance " << distance << " against typical bond length of " << MinMaxDistanceSquared << ".");
                        const bool status = MinMaxDistanceSquared.isInRange(distance);
                        if (OtherWalker->father > Walker->father ) { // just to not add bonds from both sides
                          if (status) { // create bond if distance is smaller
                            LOG(1, "ACCEPT: Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << ".");
                            bond * Binder = new bond(Walker->father, OtherWalker->father, 1);
                            Walker->father->RegisterBond(WorldTime::getTime(), Binder);
                            OtherWalker->father->RegisterBond(WorldTime::getTime(), Binder);
                            BondCount++;
                          } else {
                            LOG(1, "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.");
}