source: src/moleculelist.cpp@ 311137

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/** \file MoleculeListClass.cpp
 *
 * Function implementations for the class MoleculeListClass.
 *
 */

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

#include "CodePatterns/MemDebug.hpp"

#include <iostream>

//#include <gsl/gsl_inline.h>
#include <gsl/gsl_heapsort.h>

#include "MoleculeListClass.hpp"

#include "CodePatterns/Log.hpp"

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Box.hpp"
#include "config.hpp"
#include "Element/element.hpp"
#include "Element/periodentafel.hpp"
#include "Fragmentation/Graph.hpp"
#include "Fragmentation/KeySet.hpp"
#include "Graph/BondGraph.hpp"
#include "Helpers/helpers.hpp"
#include "molecule.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "Parser/FormatParserStorage.hpp"
#include "World.hpp"


/** Constructor for MoleculeListClass.
 */
MoleculeListClass::MoleculeListClass(World *_world) :
  Observable("MoleculeListClass"),
  MaxIndex(1),
  world(_world)
{};

/** Destructor for MoleculeListClass.
 */
MoleculeListClass::~MoleculeListClass()
{
  LOG(4, "Clearing ListOfMolecules.");
  for(MoleculeList::iterator MolRunner = ListOfMolecules.begin(); MolRunner != ListOfMolecules.end(); ++MolRunner)
    (*MolRunner)->signOff(this);
  ListOfMolecules.clear(); // empty list
};

/** Insert a new molecule into the list and set its number.
 * \param *mol molecule to add to list.
 */
void MoleculeListClass::insert(molecule *mol)
{
  OBSERVE;
  mol->IndexNr = MaxIndex++;
  ListOfMolecules.push_back(mol);
  mol->signOn(this);
};

/** Erases a molecule from the list.
 * \param *mol molecule to add to list.
 */
void MoleculeListClass::erase(molecule *mol)
{
  OBSERVE;
  mol->signOff(this);
  ListOfMolecules.remove(mol);
};

/** Comparison function for two values.
 * \param *a
 * \param *b
 * \return <0, \a *a less than \a *b, ==0 if equal, >0 \a *a greater than \a *b
 */
int CompareDoubles (const void * a, const void * b)
{
  if (*(double *)a > *(double *)b)
    return -1;
  else if (*(double *)a < *(double *)b)
    return 1;
  else
    return 0;
};


/** Compare whether two molecules are equal.
 * \param *a molecule one
 * \param *n molecule two
 * \return lexical value (-1, 0, +1)
 */
int MolCompare(const void *a, const void *b)
{
  int *aList = NULL, *bList = NULL;
  int Count, Counter, aCounter, bCounter;
  int flag;

  // sort each atom list and put the numbers into a list, then go through
  //LOG(0, "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << ".");
  // Yes those types are awkward... but check it for yourself it checks out this way
  molecule *const *mol1_ptr= static_cast<molecule *const *>(a);
  molecule *mol1 = *mol1_ptr;
  molecule *const *mol2_ptr= static_cast<molecule *const *>(b);
  molecule *mol2 = *mol2_ptr;
  if (mol1->getAtomCount() < mol2->getAtomCount()) {
    return -1;
  } else {
    if (mol1->getAtomCount() > mol2->getAtomCount())
      return +1;
    else {
      Count = mol1->getAtomCount();
      aList = new int[Count];
      bList = new int[Count];

      // fill the lists
      Counter = 0;
      aCounter = 0;
      bCounter = 0;
      molecule::const_iterator aiter = mol1->begin();
      molecule::const_iterator biter = mol2->begin();
      for (;(aiter != mol1->end()) && (biter != mol2->end());
          ++aiter, ++biter) {
        if ((*aiter)->GetTrueFather() == NULL)
          aList[Counter] = Count + (aCounter++);
        else
          aList[Counter] = (*aiter)->GetTrueFather()->getNr();
        if ((*biter)->GetTrueFather() == NULL)
          bList[Counter] = Count + (bCounter++);
        else
          bList[Counter] = (*biter)->GetTrueFather()->getNr();
        Counter++;
      }
      // check if AtomCount was for real
      flag = 0;
      if ((aiter == mol1->end()) && (biter != mol2->end())) {
        flag = -1;
      } else {
        if ((aiter != mol1->end()) && (biter == mol2->end()))
          flag = 1;
      }
      if (flag == 0) {
        // sort the lists
        gsl_heapsort(aList, Count, sizeof(int), CompareDoubles);
        gsl_heapsort(bList, Count, sizeof(int), CompareDoubles);
        // compare the lists

        flag = 0;
        for (int i = 0; i < Count; i++) {
          if (aList[i] < bList[i]) {
            flag = -1;
          } else {
            if (aList[i] > bList[i])
              flag = 1;
          }
          if (flag != 0)
            break;
        }
      }
      delete[] (aList);
      delete[] (bList);
      return flag;
    }
  }
  return -1;
};

/** Output of a list of all molecules.
 * \param *out output stream
 */
void MoleculeListClass::Enumerate(std::ostream *out)
{
  periodentafel *periode = World::getInstance().getPeriode();
  std::map<atomicNumber_t,unsigned int> counts;
  double size=0;
  Vector Origin;

  // header
  (*out) << "Index\tName\t\tAtoms\tFormula\tCenter\tSize" << endl;
  (*out) << "-----------------------------------------------" << endl;
  if (ListOfMolecules.size() == 0)
    (*out) << "\tNone" << endl;
  else {
    Origin.Zero();
    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
      // count atoms per element and determine size of bounding sphere
      size=0.;
      for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
        counts[(*iter)->getType()->getAtomicNumber()]++;
        if ((*iter)->DistanceSquared(Origin) > size)
          size = (*iter)->DistanceSquared(Origin);
      }
      // output Index, Name, number of atoms, chemical formula
      (*out) << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->getAtomCount() << "\t";

      std::map<atomicNumber_t,unsigned int>::reverse_iterator iter;
      for(iter=counts.rbegin(); iter!=counts.rend();++iter){
        atomicNumber_t Z =(*iter).first;
        (*out) << periode->FindElement(Z)->getSymbol() << (*iter).second;
      }
      // Center and size
      Vector *Center = (*ListRunner)->DetermineCenterOfAll();
      (*out) << "\t" << *Center << "\t" << sqrt(size) << endl;
      delete(Center);
    }
  }
};

/** Returns the molecule with the given index \a index.
 * \param index index of the desired molecule
 * \return pointer to molecule structure, NULL if not found
 */
molecule * MoleculeListClass::ReturnIndex(int index)
{
  for(MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
    if ((*ListRunner)->IndexNr == index)
      return (*ListRunner);
  return NULL;
};


/** Simple output of the pointers in ListOfMolecules.
 * \param *out output stream
 */
void MoleculeListClass::Output(std::ostream *out)
{
  if (DoLog(1)) {
    std::stringstream output;
    output << "MoleculeList: ";
    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
      output << *ListRunner << "\t";
      LOG(1, output.str());
  }
};

/** Returns a string with \a i prefixed with 0s to match order of total number of molecules in digits.
 * \param FragmentNumber total number of fragments to determine necessary number of digits
 * \param digits number to create with 0 prefixed
 * \return allocated(!) char array with number in digits, ten base.
 */
inline char *FixedDigitNumber(const int FragmentNumber, const int digits)
{
  char *returnstring;
  int number = FragmentNumber;
  int order = 0;
  while (number != 0) { // determine number of digits needed
    number = (int)floor(((double)number / 10.));
    order++;
    //LOG(0, "Number is " << number << ", order is " << order << ".");
  }
  // allocate string
  returnstring = new char[order + 2];
  // terminate  and fill string array from end backward
  returnstring[order] = '\0';
  number = digits;
  for (int i=order;i--;) {
    returnstring[i] = '0' + (char)(number % 10);
    number = (int)floor(((double)number / 10.));
  }
  //LOG(0, returnstring);
  return returnstring;
};

/** Calculates necessary hydrogen correction due to unwanted interaction between saturated ones.
 * If for a pair of two hydrogen atoms a and b, at least is a saturated one, and a and b are not
 * bonded to the same atom, then we add for this pair a correction term constructed from a Morse
 * potential function fit to QM calculations with respecting to the interatomic hydrogen distance.
 * \param &path path to file
 */
bool MoleculeListClass::AddHydrogenCorrection(std::string &path)
{
  const bond *Binder = NULL;
  double ***FitConstant = NULL, **correction = NULL;
  int a, b;
  ofstream output;
  ifstream input;
  string line;
  stringstream zeile;
  double distance;
  char ParsedLine[1023];
  double tmp;
  char *FragmentNumber = NULL;

  LOG(1, "Saving hydrogen saturation correction ... ");
  // 0. parse in fit constant files that should have the same dimension as the final energy files
  // 0a. find dimension of matrices with constants
  line = path;
  line += "1";
  line += FITCONSTANTSUFFIX;
  input.open(line.c_str());
  if (input.fail()) {
    LOG(1, endl << "Unable to open " << line << ", is the directory correct?");
    return false;
  }
  a = 0;
  b = -1; // we overcount by one
  while (!input.eof()) {
    input.getline(ParsedLine, 1023);
    zeile.str(ParsedLine);
    int i = 0;
    while (!zeile.eof()) {
      zeile >> distance;
      i++;
    }
    if (i > a)
      a = i;
    b++;
  }
  LOG(0, "I recognized " << a << " columns and " << b << " rows, ");
  input.close();

  // 0b. allocate memory for constants
  FitConstant = new double**[3];
  for (int k = 0; k < 3; k++) {
    FitConstant[k] = new double*[a];
    for (int i = a; i--;) {
      FitConstant[k][i] = new double[b];
      for (int j = b; j--;) {
        FitConstant[k][i][j] = 0.;
      }
    }
  }
  // 0c. parse in constants
  for (int i = 0; i < 3; i++) {
    line = path;
    line.append("/");
    line += FRAGMENTPREFIX;
    sprintf(ParsedLine, "%d", i + 1);
    line += ParsedLine;
    line += FITCONSTANTSUFFIX;
    input.open(line.c_str());
    if (input.fail()) {
      ELOG(0, endl << "Unable to open " << line << ", is the directory correct?");
      performCriticalExit();
      return false;
    }
    int k = 0, l;
    while ((!input.eof()) && (k < b)) {
      input.getline(ParsedLine, 1023);
      //LOG(1, "INFO: Current Line: " << ParsedLine);
      zeile.str(ParsedLine);
      zeile.clear();
      l = 0;
      //std::stringstream output;
      while ((!zeile.eof()) && (l < a)) {
        zeile >> FitConstant[i][l][k];
        //output << FitConstant[i][l][k] << "\t";
        l++;
      }
      //LOG(1, "INFO: fit constant are " << output.str());
      k++;
    }
    input.close();
  }
  if (DoLog(1)) {
    for (int k = 0; k < 3; k++) {
      std::stringstream output;
      output << "Constants " << k << ": ";
      for (int j = 0; j < b; j++) {
        for (int i = 0; i < a; i++) {
          output << FitConstant[k][i][j] << "\t";
        }
        output << std::endl;
      }
      LOG(0, output.str());
    }
  }

  // 0d. allocate final correction matrix
  correction = new double*[a];
  for (int i = a; i--;)
    correction[i] = new double[b];

  // 1a. go through every molecule in the list
  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    // 1b. zero final correction matrix
    for (int k = a; k--;)
      for (int j = b; j--;)
        correction[k][j] = 0.;
    // 2. take every hydrogen that is a saturated one
    for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
      //LOG(1, "(*iter): " << *(*iter) << " with first bond " << *((*iter)->getListOfBonds().begin()) << ".");
      if (((*iter)->getType()->getAtomicNumber() == 1) && (((*iter)->father == NULL)
          || ((*iter)->father->getType()->getAtomicNumber() != 1))) { // if it's a hydrogen
        for (molecule::const_iterator runner = (*ListRunner)->begin(); runner != (*ListRunner)->end(); ++runner) {
          //LOG(2, "Runner: " << *(*runner) << " with first bond " << *((*iter)->getListOfBonds().begin()) << ".");
          // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
          const BondList &bondlist = (*runner)->getListOfBonds();
          Binder = *(bondlist.begin());
          if (((*runner)->getType()->getAtomicNumber() == 1) && ((*runner)->getNr() > (*iter)->getNr()) && (Binder->GetOtherAtom((*runner)) != Binder->GetOtherAtom((*iter)))) { // (hydrogens have only one bonding partner!)
            // 4. evaluate the morse potential for each matrix component and add up
            distance = (*runner)->distance(*(*iter));
            //std::stringstream output;
            //output << "Fragment " << (*ListRunner)->name << ": " << *(*runner) << "<= " << distance << "=>" << *(*iter) << ":";
            for (int k = 0; k < a; k++) {
              for (int j = 0; j < b; j++) {
                switch (k) {
                  case 1:
                  case 7:
                  case 11:
                    tmp = pow(FitConstant[0][k][j] * (1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]))), 2);
                    break;
                  default:
                    tmp = FitConstant[0][k][j] * pow(distance, FitConstant[1][k][j]) + FitConstant[2][k][j];
                    break;
                };
                correction[k][j] -= tmp; // ground state is actually lower (disturbed by additional interaction)
                //output << tmp << "\t";
              }
              //output << endl;
            }
            //LOG(0, output.str());
          }
        }
      }
    }
    // 5. write final matrix to file
    line = path;
    line.append("/");
    line += FRAGMENTPREFIX;
    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr);
    line += FragmentNumber;
    delete[] (FragmentNumber);
    line += HCORRECTIONSUFFIX;
    output.open(line.c_str());
    output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
    for (int j = 0; j < b; j++) {
      for (int i = 0; i < a; i++)
        output << correction[i][j] << "\t";
      output << endl;
    }
    output.close();
  }
  for (int i = a; i--;)
    delete[](correction[i]);
  delete[](correction);

  line = path;
  line.append("/");
  line += HCORRECTIONSUFFIX;
  output.open(line.c_str());
  output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
  for (int j = 0; j < b; j++) {
    for (int i = 0; i < a; i++)
      output << 0 << "\t";
    output << endl;
  }
  output.close();
  // 6. free memory of parsed matrices
  for (int k = 0; k < 3; k++) {
    for (int i = a; i--;) {
      delete[](FitConstant[k][i]);
    }
    delete[](FitConstant[k]);
  }
  delete[](FitConstant);
  LOG(0, "done.");
  return true;
};

/** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
 * \param &path path to file
 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
 * \return true - file written successfully, false - writing failed
 */
bool MoleculeListClass::StoreForcesFile(std::string &path, int *SortIndex)
{
  bool status = true;
  string filename(path);
  filename += FORCESFILE;
  ofstream ForcesFile(filename.c_str());
  periodentafel *periode=World::getInstance().getPeriode();

  // open file for the force factors
  LOG(1, "Saving  force factors ... ");
  if (!ForcesFile.fail()) {
    //output << prefix << "Forces" << endl;
    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
      periodentafel::const_iterator elemIter;
      for(elemIter=periode->begin();elemIter!=periode->end();++elemIter){
        if ((*ListRunner)->hasElement((*elemIter).first)) { // if this element got atoms
          for(molecule::iterator atomIter = (*ListRunner)->begin(); atomIter !=(*ListRunner)->end();++atomIter){
            if ((*atomIter)->getType()->getAtomicNumber() == (*elemIter).first) {
              if (((*atomIter)->GetTrueFather() != NULL) && ((*atomIter)->GetTrueFather() != (*atomIter))) {// if there is a rea
                ForcesFile << SortIndex[(*atomIter)->GetTrueFather()->getNr()] << "\t";
              } else
                // otherwise a -1 to indicate an added saturation hydrogen
                ForcesFile << "-1\t";
            }
          }
        }
      }
      ForcesFile << endl;
    }
    ForcesFile.close();
    LOG(1, "done.");
  } else {
    status = false;
    LOG(1, "failed to open file " << filename << ".");
  }
  ForcesFile.close();

  return status;
};

/** Writes a config file for each molecule in the given \a **FragmentList.
 * \param *out output stream for debugging
 * \param &prefix path and prefix to the fragment config files
 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
 * \param type desired type to store
 * \return true - success (each file was written), false - something went wrong.
 */
bool MoleculeListClass::OutputConfigForListOfFragments(std::string &prefix, int *SortIndex, ParserTypes type)
{
  ofstream outputFragment;
  std::string FragmentName;
  bool result = true;
  bool intermediateResult = true;
  Vector BoxDimension;
  char *FragmentNumber = NULL;
  int FragmentCounter = 0;
  ofstream output;
  RealSpaceMatrix cell_size = World::getInstance().getDomain().getM();
  RealSpaceMatrix cell_size_backup = cell_size;
  int count=0;

  // store the fragments as config and as xyz
  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    // correct periodic
    if ((*ListRunner)->ScanForPeriodicCorrection()) {
      count++;
    }

    {
      // list atoms in fragment for debugging
      std::stringstream output;
      output << "Contained atoms: ";
      for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
        output << (*iter)->getName() << " ";
      }
      LOG(2, output.str());
    }

    {
      // center on edge
      (*ListRunner)->CenterEdge(&BoxDimension);
      for (int k = 0; k < NDIM; k++)  // if one edge is to small, set at least to 1 angstroem
        if (BoxDimension[k] < 1.)
          BoxDimension[k] += 1.;
      (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
      for (int k = 0; k < NDIM; k++) {
        BoxDimension[k] = 2.5 * (World::getInstance().getConfig()->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
        cell_size.at(k,k) = BoxDimension[k] * 2.;
      }
      World::getInstance().setDomain(cell_size);
      (*ListRunner)->Translate(&BoxDimension);

      // output file
      std::vector<atom *> atoms;
      // TODO: Convert iterator to const_iterator when FormatParserStorage::save() has vector<const atom *>
      // We need iterator here because FormatParserStorage::save() need vector<atom *> not const refs.
      for (molecule::iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
        atoms.push_back(*iter);
      }
      FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++);
      FragmentName = prefix + FragmentNumber + "." + FormatParserStorage::getInstance().getSuffixFromType(type);
      outputFragment.open(FragmentName.c_str(), ios::out);
      std::stringstream output;
      output << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as XYZ ... ";
      if ((intermediateResult = FormatParserStorage::getInstance().save(
          outputFragment,
          FormatParserStorage::getInstance().getSuffixFromType(type),
          atoms)))
        output << " done.";
      else
        output << " failed.";
      LOG(3, output.str());
      delete[](FragmentNumber);

      result = result && intermediateResult;
      outputFragment.close();
      outputFragment.clear();
    }
  }
  LOG(0, "STATUS: done.");

  // printing final number
  LOG(2, "INFO: Final number of fragments: " << FragmentCounter << ".");

  // printing final number
  LOG(2, "INFO: For " << count << " fragments periodic correction would have been necessary.");

  // restore cell_size
  World::getInstance().setDomain(cell_size_backup);

  return result;
};

/** Counts the number of molecules with the molecule::ActiveFlag set.
 * \return number of molecules with ActiveFlag set to true.
 */
int MoleculeListClass::NumberOfActiveMolecules()
{
  int count = 0;
  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
    count += ((*ListRunner)->ActiveFlag ? 1 : 0);
  return count;
};

/** Count all atoms in each molecule.
 * \return number of atoms in the MoleculeListClass.
 * TODO: the inner loop should be done by some (double molecule::CountAtom()) function
 */
int MoleculeListClass::CountAllAtoms() const
{
  int AtomNo = 0;
  for (MoleculeList::const_iterator MolWalker = ListOfMolecules.begin(); MolWalker != ListOfMolecules.end(); MolWalker++) {
    AtomNo += (*MolWalker)->size();
  }
  return AtomNo;
}