source: src/Fragmentation/HessianMatrix.cpp@ a9b86d

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

/*
 * HessianMatrix.cpp
 *
 *  Created on: Sep 15, 2011
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include <cstring>
#include <fstream>

#include "CodePatterns/Log.hpp"
#include "KeySetsContainer.hpp"

#include "Fragmentation/helpers.hpp"
#include "Helpers/defs.hpp"
#include "Helpers/helpers.hpp"

#include "HessianMatrix.hpp"

/** Parsing force Indices of each fragment
 * \param *name directory with \a ForcesFile
 * \return parsing successful
 */
bool HessianMatrix::ParseIndices(char *name)
{
  std::ifstream input;
  char *FragmentNumber = NULL;
  char filename[1023];
  stringstream line;

  DoLog(0) && (Log() << Verbose(0) << "Parsing hessian indices for " << MatrixCounter << " matrices." << endl);
  Indices = new int*[MatrixCounter + 1];
  line << name << FRAGMENTPREFIX << FORCESFILE;
  input.open(line.str().c_str(), ios::in);
  //Log() << Verbose(0) << "Opening " << line.str() << " ... "  << input << endl;
  if (input.fail()) {
    DoLog(0) && (Log() << Verbose(0) << endl << "HessianMatrix::ParseIndices: Unable to open " << line.str() << ", is the directory correct?" << endl);
    return false;
  }
  for (int i=0;(i<MatrixCounter) && (!input.eof());i++) {
    // get the number of atoms for this fragment
    input.getline(filename, 1023);
    line.str(filename);
    // parse the values
    Indices[i] = new int[RowCounter[i]];
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    //Log() << Verbose(0) << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    delete[](FragmentNumber);
    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
      line >> Indices[i][j];
      //Log() << Verbose(0) << " " << Indices[i][j];
    }
    //Log() << Verbose(0) << endl;
  }
  Indices[MatrixCounter] = new int[RowCounter[MatrixCounter]];
  for(int j=RowCounter[MatrixCounter];j--;) {
    Indices[MatrixCounter][j] = j;
  }
  input.close();
  return true;
};


/** Sums the hessian entries and puts into last element of \a HessianMatrix::Matrix.
 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
 * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 *  \param sign +1 or -1
 * \return true if summing was successful
 */
bool HessianMatrix::SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySets, int Order, double sign)
{
  int FragmentNr;
  // sum forces
  for(int i=0;i<KeySets.FragmentsPerOrder[Order];i++) {
    FragmentNr = KeySets.OrderSet[Order][i];
    for(int l=0;l<RowCounter[ FragmentNr ];l++) {
      int j = Indices[ FragmentNr ][l];
      if (j > RowCounter[MatrixCounter]) {
        DoeLog(0) && (eLog()<< Verbose(0) << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl);
        performCriticalExit();
        return false;
      }
      if (j != -1) {
        for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
          int k = Indices[ FragmentNr ][m];
          if (k > ColumnCounter[MatrixCounter]) {
            DoeLog(0) && (eLog()<< Verbose(0) << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl);
            performCriticalExit();
            return false;
          }
          if (k != -1) {
            //Log() << Verbose(0) << "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "]." << endl;
            Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][m];
          }
        }
      }
    }
  }
  return true;
};

/** Constructor for class HessianMatrix.
 */
HessianMatrix::HessianMatrix() :
  MatrixContainer(),
  IsSymmetric(true)
{}

/** Sums the hessian entries with each factor and put into last element of \a ***Matrix.
 * Sums over "E"-terms to create the "F"-terms
 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies.
 * \param KeySets KeySetContainer with bond Order and association mapping of each fragment to an order
 * \param Order bond order
 * \return true if summing was successful
 */
bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySets, int Order)
{
  // go through each order
  for (int CurrentFragment=0;CurrentFragment<KeySets.FragmentsPerOrder[Order];CurrentFragment++) {
    //Log() << Verbose(0) << "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
    // then go per order through each suborder and pick together all the terms that contain this fragment
    for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order
      for (int j=0;j<KeySets.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder
        if (KeySets.Contains(KeySets.OrderSet[Order][CurrentFragment], KeySets.OrderSet[SubOrder][j])) {
          //Log() << Verbose(0) << "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << "." << endl;
          // if the fragment's indices are all in the current fragment
          for(int k=0;k<RowCounter[ KeySets.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment
            int m = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][k];
            //Log() << Verbose(0) << "Current row index is " << k << "/" << m << "." << endl;
            if (m != -1) { // if it's not an added hydrogen
              for (int l=0;l<RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment
                //Log() << Verbose(0) << "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << "." << endl;
                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
                  m = l;
                  break;
                }
              }
              //Log() << Verbose(0) << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl;
              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
                DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
                performCriticalExit();
                return false;
              }

              for(int l=0;l<ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l++) {
                int n = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][l];
                //Log() << Verbose(0) << "Current column index is " << l << "/" << n << "." << endl;
                if (n != -1) { // if it's not an added hydrogen
                  for (int p=0;p<ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment
                    //Log() << Verbose(0) << "Comparing " << n << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p] << "." << endl;
                    if (n == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p]) {
                      n = p;
                      break;
                    }
                  }
                  //Log() << Verbose(0) << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl;
                  if (n > ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
                    DoeLog(0) && (eLog()<< Verbose(0) << "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!" << endl);
                    performCriticalExit();
                    return false;
                  }
                  if (Order == SubOrder) { // equal order is always copy from Energies
                    //Log() << Verbose(0) << "Adding " << MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                    Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
                  } else {
                    //Log() << Verbose(0) << "Subtracting " << Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "]." << endl;
                    Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
                  }
                }
              }
            }
            //if ((ColumnCounter[ KeySets.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1))
             //Log() << Verbose(0) << "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl;
          }
        } else {
          //Log() << Verbose(0) << "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << "." << endl;
        }
      }
    }
   //Log() << Verbose(0) << "Final Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << KeySets.AtomCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][KeySets.AtomCounter[0]-1][1] << endl;
  }

  return true;
};

/** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix.
 * \param *name directory with files
 * \param *prefix prefix of each matrix file
 * \param *suffix suffix of each matrix file
 * \param skiplines number of inital lines to skip
 * \param skiplines number of inital columns to skip
 * \return parsing successful
 */
bool HessianMatrix::ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns)
{
  char filename[1023];
  std::ifstream input;
  stringstream file;
  int nr;
  bool status = MatrixContainer::ParseFragmentMatrix(name, prefix, suffix, skiplines, skipcolumns);

  if (status) {
    // count number of atoms for last plus one matrix
    file << name << FRAGMENTPREFIX << KEYSETFILE;
    input.open(file.str().c_str(), ios::in);
    if (input.fail()) {
      DoLog(0) && (Log() << Verbose(0) << endl << "HessianMatrix::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?" << endl);
      return false;
    }
    RowCounter[MatrixCounter] = 0;
    ColumnCounter[MatrixCounter] = 0;
    while (!input.eof()) {
      input.getline(filename, 1023);
      stringstream zeile(filename);
      while (!zeile.eof()) {
        zeile >> nr;
        //Log() << Verbose(0) << "Current index: " << getNr() << "." << endl;
        if (nr > RowCounter[MatrixCounter]) {
          RowCounter[MatrixCounter] = nr;
          ColumnCounter[MatrixCounter] = nr;
        }
      }
    }
    RowCounter[MatrixCounter]++;    // Nr start at 0, count starts at 1
    ColumnCounter[MatrixCounter]++;    // Nr start at 0, count starts at 1
    input.close();

    // allocate last plus one matrix
    DoLog(0) && (Log() << Verbose(0) << "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns." << endl);
    Matrix[MatrixCounter] = new double*[RowCounter[MatrixCounter] + 1];
    for(int j=0;j<=RowCounter[MatrixCounter];j++)
      Matrix[MatrixCounter][j] = new double[ColumnCounter[MatrixCounter]];

    // try independently to parse global forcesuffix file if present
    strncpy(filename, name, 1023);
    strncat(filename, prefix, 1023-strlen(filename));
    strncat(filename, suffix.c_str(), 1023-strlen(filename));
    std::ifstream input(filename);
    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
    input.close();
  }


  return status;
};