source: src/bondgraph.cpp@ e7350d4

/*
 * 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 "element.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Verbose.hpp"
#include "CodePatterns/Log.hpp"
#include "molecule.hpp"
#include "parser.hpp"
#include "periodentafel.hpp"
#include "LinearAlgebra/Vector.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),
    max_distance(0),
    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) {
    DoLog(1) && (Log() << Verbose(1) << "MatrixContainer for Bond length already present, removing." << endl);
    delete(BondLengthMatrix);
  }
  TempContainer = new MatrixContainer;

  // parse in matrix
  if ((input.good()) && (status = TempContainer->ParseMatrix(input, 0, 1, 0))) {
    DoLog(1) && (Log() << Verbose(1) << "Parsing bond length matrix successful." << endl);
  } else {
    DoeLog(1) && (eLog()<< Verbose(1) << "Parsing bond length matrix failed." << endl);
    status = false;
  }

  // find greatest distance
  max_distance=0;
  if (status) {
    for(int i=0;i<TempContainer->RowCounter[0];i++)
      for(int j=i;j<TempContainer->ColumnCounter[0];j++)
        if (TempContainer->Matrix[0][i][j] > max_distance)
          max_distance = TempContainer->Matrix[0][i][j];
  }
  max_distance += BondThreshold;

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

bool BondGraph::ConstructBondGraph(molecule * const mol)
{
  Info FunctionInfo(__func__);
  bool status = true;

  if (mol->empty()) // only construct if molecule is not empty
    return false;

  SetMaxDistanceToMaxOfCovalentRadii(mol);
  mol->CreateAdjacencyList(max_distance, IsAngstroem, &BondGraph::getMinMaxDistance, this);

  return status;
};

double BondGraph::GetBondLength(int firstZ, int secondZ)
{
  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]);
  }
};

double BondGraph::SetMaxDistanceToMaxOfCovalentRadii(const molecule * const mol)
{
  Info FunctionInfo(__func__);
  max_distance = 0.;

  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    if ((*iter)->getType()->getCovalentRadius() > max_distance)
      max_distance = (*iter)->getType()->getCovalentRadius();
  }
  max_distance *= 2.;
  max_distance += BondThreshold;

  return max_distance;
};

double BondGraph::getMaxDistance() const
{
  return max_distance;
}


void BondGraph::CovalentMinMaxDistance(BondedParticle * const Walker, BondedParticle * const OtherWalker, double &MinDistance, double &MaxDistance, bool IsAngstroem)
{
  MinDistance = OtherWalker->getType()->getCovalentRadius() + Walker->getType()->getCovalentRadius();
  MinDistance *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
  MaxDistance = MinDistance + BondThreshold;
  MinDistance -= BondThreshold;
};

void BondGraph::BondLengthMatrixMinMaxDistance(BondedParticle * const Walker, BondedParticle * const OtherWalker, double &MinDistance, double &MaxDistance, bool IsAngstroem)
{
  ASSERT(BondLengthMatrix != NULL,
      "BondGraph::BondLengthMatrixMinMaxDistance() called without NULL BondLengthMatrix.");
  MinDistance = GetBondLength(Walker->getType()->getAtomicNumber()-1, OtherWalker->getType()->getAtomicNumber()-1);
  MinDistance *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
  MaxDistance = MinDistance + BondThreshold;
  MinDistance -= BondThreshold;
};

void BondGraph::getMinMaxDistance(BondedParticle * const Walker, BondedParticle * const OtherWalker, double &MinDistance, double &MaxDistance, bool IsAngstroem)
{
  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, MinDistance, MaxDistance, IsAngstroem);
  } else {
    LOG(2, "INFO: Using Covalent radii criterion for [min,max] distances.");
    BondLengthMatrixMinMaxDistance(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
  }
}