source: src/Fragmentation/HessianMatrix.cpp@ d4d7a1

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

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

  LOG(0, "Parsing hessian indices for " << MatrixCounter << " matrices.");
  Indices.resize(MatrixCounter + 1);
  line << name << FRAGMENTPREFIX << FORCESFILE;
  input.open(line.str().c_str(), ios::in);
  //LOG(0, "Opening " << line.str() << " ... "  << input);
  if (input.fail()) {
    LOG(0, endl << "HessianMatrix::ParseIndices: Unable to open " << line.str() << ", is the directory correct?");
    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].resize(RowCounter[i]);
    FragmentNumber = FixedDigitNumber(MatrixCounter, i);
    //std::stringstream output;
    //output << FRAGMENTPREFIX << FragmentNumber << "[" << RowCounter[i] << "]:";
    delete[](FragmentNumber);
    for(int j=0;(j<RowCounter[i]) && (!line.eof());j++) {
      line >> Indices[i][j];
      //output << " " << Indices[i][j];
    }
    //LOG(0, output.str());
  }
  Indices[MatrixCounter].resize(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]) {
        ELOG(0, "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!");
        performCriticalExit();
        return false;
      }
      if (j != -1) {
        for(int m=0;m<ColumnCounter[ FragmentNr ];m++) {
          int k = Indices[ FragmentNr ][m];
          if (k > ColumnCounter[MatrixCounter]) {
            ELOG(0, "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!");
            performCriticalExit();
            return false;
          }
          if (k != -1) {
            //LOG(0, "Adding " << sign*Fragments.Matrix[ FragmentNr ][l][m] << " from [" << l << "][" << m << "] onto [" << j << "][" << k << "].");
            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(0, "Current Fragment is " << CurrentFragment << "/" << KeySets.OrderSet[Order][CurrentFragment] << ".");
    // 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(0, "Current other fragment is " << j << "/" << KeySets.OrderSet[SubOrder][j] << ".");
          // 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(0, "Current row index is " << k << "/" << m << ".");
            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(0, "Comparing " << m << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l] << ".");
                if (m == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][l]) {
                  m = l;
                  break;
                }
              }
              //LOG(0, "Corresponding row index for " << k << " in CurrentFragment is " << m << ".");
              if (m > RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
                ELOG(0, "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current row index " << m << " is greater than " << RowCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!");
                performCriticalExit();
                return false;
              }

              for(int l=0;l<ColumnCounter[ KeySets.OrderSet[SubOrder][j] ];l++) {
                int n = MatrixValues.Indices[ KeySets.OrderSet[SubOrder][j] ][l];
                //LOG(0, "Current column index is " << l << "/" << n << ".");
                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(0, "Comparing " << n << " with " << MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p] << ".");
                    if (n == MatrixValues.Indices[ KeySets.OrderSet[Order][CurrentFragment] ][p]) {
                      n = p;
                      break;
                    }
                  }
                  //LOG(0, "Corresponding column index for " << l << " in CurrentFragment is " << n << ".");
                  if (n > ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ]) {
                    ELOG(0, "In fragment No. " << KeySets.OrderSet[Order][CurrentFragment]   << " current column index " << n << " is greater than " << ColumnCounter[ KeySets.OrderSet[Order][CurrentFragment] ] << "!");
                    performCriticalExit();
                    return false;
                  }
                  if (Order == SubOrder) { // equal order is always copy from Energies
                    //LOG(0, "Adding " << MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "].");
                    Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l];
                  } else {
                    //LOG(0, "Subtracting " << Matrix[ KeySets.OrderSet[SubOrder][j] ][k][l] << " from [" << k << "][" << l << "] onto [" << m << "][" << n << "].");
                    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(0, "Fragments[ KeySets.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySets.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " <<  Matrix[ KeySets.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1]);
          }
        } else {
          //LOG(0, "Fragment " << KeySets.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySets.OrderSet[Order][CurrentFragment] << ".");
        }
      }
    }
   //LOG(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]);
  }

  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, std::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()) {
      LOG(0, endl << "HessianMatrix::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?");
      return false;
    }
    RowCounter[MatrixCounter] = 0;
    ColumnCounter[MatrixCounter] = 0;
    while (!input.eof()) {
      input.getline(filename, 1023);
      stringstream zeile(filename);
      while (!zeile.eof()) {
        zeile >> nr;
        //LOG(0, "Current index: " << getNr() << ".");
        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
    LOG(0, "Allocating last plus one matrix with " << (RowCounter[MatrixCounter]+1) << " rows and " << ColumnCounter[MatrixCounter] << " columns.");
    if ((int)Matrix[MatrixCounter].size() <= RowCounter[MatrixCounter] + 2)
      Matrix[MatrixCounter].resize(RowCounter[MatrixCounter] + 1);
    for(int j=0;j<=RowCounter[MatrixCounter];j++)
      if ((int)Matrix[MatrixCounter][j].size() <= ColumnCounter[MatrixCounter]+1)
        Matrix[MatrixCounter][j].resize(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;
};