source: src/Fragmentation/ForceMatrix.cpp@ 50e4e5

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

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

  DoLog(0) && (Log() << Verbose(0) << "Parsing force indices for " << MatrixCounter << " matrices." << endl);
  Indices.resize(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 << "ForceMatrix::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].resize(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].resize(RowCounter[MatrixCounter]);
  for(int j=RowCounter[MatrixCounter];j--;) {
    Indices[MatrixCounter][j] = j;
  }
  input.close();
  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]) {
        DoeLog(0) && (eLog()<< Verbose(0) << "Current force index " << j << " is greater than " << RowCounter[MatrixCounter] << "!" << endl);
        performCriticalExit();
        return false;
      }
      if (j != -1) {
        //if (j == 0) Log() << Verbose(0) << "Summing onto ion 0, type 0 from fragment " << FragmentNr << ", ion " << l << "." << endl;
        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, 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()) {
      DoLog(0) && (Log() << Verbose(0) << endl << "ForceMatrix::ParseFragmentMatrix: Unable to open " << file.str() << ", is the directory correct?" << endl);
      return false;
    }
    RowCounter[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;
      }
    }
    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 (Matrix[MatrixCounter].size() <= RowCounter[MatrixCounter] + 2)
      Matrix[MatrixCounter].resize(RowCounter[MatrixCounter] + 1);
    for(int j=0;j<=RowCounter[MatrixCounter];j++)
      if (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;
};