source: src/Fragmentation/ForceMatrix.cpp

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. 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/>.
 */

/*
 * ForceMatrix.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 "ForceMatrix.hpp"

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

  LOG(0, "Parsing force 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 << "ForceMatrix::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;
};

/** Parsing force Indices of each fragment
 * \param *name directory with \a ForcesFile
 * \return parsing successful
 */
bool ForceMatrix::ParseIndices(const KeySetsContainer &ForceKeySet)
{
  Indices = ForceKeySet.KeySets;
  return true;
}

/** Sums the forces and puts into last element of \a ForceMatrix::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 ForceMatrix::SumSubForces(class ForceMatrix &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 force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!");
        performCriticalExit();
        return false;
      }
      if (j != -1) {
        //if (j == 0) LOG(0, "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << ".");
        for(int k=2;k<ColumnCounter[MatrixCounter];k++)
          Matrix[MatrixCounter][j][k] += sign*Fragments.Matrix[ FragmentNr ][l][k];
      }
    }
  }
  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 ForceMatrix::ParseFragmentMatrix(const char *name, const char *prefix, std::string suffix, int skiplines, int skipcolumns)
{
  char filename[1023];
  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 << "ForceMatrix::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?");
      return false;
    }
    RowCounter[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;
      }
    }
    RowCounter[MatrixCounter]++;    // Nr start at 0, count starts at 1
    input.close();

    ColumnCounter[MatrixCounter] = 0;
    for(int j=0; j < MatrixCounter;j++) { // (energy matrix might be bigger than number of atoms in terms of rows)
      if (ColumnCounter[j] > ColumnCounter[MatrixCounter])  // take maximum of all for last matrix
        ColumnCounter[MatrixCounter] = ColumnCounter[j];
    }

    // allocate last plus one matrix
    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;
};