source: src/Graph/BondGraph.cpp@ 7b0494

/*
 * 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/>.
 */

/*
 * 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/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
{
  double return_length;
  if ((firstZ < 0) || (firstZ >= (int)BondLengthMatrix->Matrix[0].size()))
    return -1.;
  if ((secondZ < 0) || (secondZ >= (int)BondLengthMatrix->Matrix[0][firstZ].size()))
    return -1.;
  if (BondLengthMatrix == NULL) {
    return_length = -1.;
  } else {
    return_length = BondLengthMatrix->Matrix[0][firstZ][secondZ];
  }
  LOG(4, "INFO: Request for length between " << firstZ << " and "
      << secondZ << ": " << return_length << ".");
  return return_length;
}

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(5, "DEBUG: Using Covalent radii criterion for [min,max) distances.");
    MinMaxDistance = CovalentMinMaxDistance(Walker, OtherWalker);
  } else {
    LOG(5, "DEBUG: Using tabulated bond table 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;
}

LinkedCell::LinkedCell_View BondGraph::getLinkedCell(const double max_distance) const
{
  return World::getInstance().getLinkedCell(max_distance);
}

std::set< const element *> BondGraph::getElementSetFromNumbers(const std::set<atomicNumber_t> &Set) const
{
  std::set< const element *> PresentElements;
  for(std::set<atomicNumber_t>::const_iterator iter = Set.begin(); iter != Set.end(); ++iter)
    PresentElements.insert( World::getInstance().getPeriode()->FindElement(*iter) );
  return PresentElements;
}

Box &BondGraph::getDomain() const
{
  return World::getInstance().getDomain();
}

unsigned int BondGraph::getTime() const
{
  return WorldTime::getTime();
}

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;
}