source: src/moleculelist.cpp@ d2a294

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

#include "molecules.hpp"

/*********************************** Functions for class MoleculeListClass *************************/

/** Constructor for MoleculeListClass.
 */
MoleculeListClass::MoleculeListClass()
{
};

/** constructor for MoleculeListClass.
 * \param NumMolecules number of molecules to allocate for
 * \param NumAtoms number of atoms to allocate for
 */
MoleculeListClass::MoleculeListClass(int NumMolecules, int NumAtoms)
{
  ListOfMolecules = (molecule **) Malloc(sizeof(molecule *)*NumMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
  for (int i=0;i<NumMolecules;i++)
    ListOfMolecules[i] = NULL;
  NumberOfMolecules = NumMolecules;
  NumberOfTopAtoms = NumAtoms;
};


/** Destructor for MoleculeListClass.
 */
MoleculeListClass::~MoleculeListClass()
{
  cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
  for (int i=0;i<NumberOfMolecules;i++) {
    if (ListOfMolecules[i] != NULL) { // if NULL don't free
      cout << Verbose(4) << "ListOfMolecules: Freeing " << ListOfMolecules[i] << "." << endl;
      delete(ListOfMolecules[i]);
      ListOfMolecules[i] = NULL;
    }
  }
  cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
  Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
};

/** 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;
  atom *aWalker = NULL;
  atom *bWalker = NULL;
  
  // sort each atom list and put the numbers into a list, then go through
  //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
  if ( (**(molecule **)a).AtomCount < (**(molecule **)b).AtomCount ) {
    return -1;
  } else { if ((**(molecule **)a).AtomCount > (**(molecule **)b).AtomCount)
    return +1;
    else {
      Count = (**(molecule **)a).AtomCount;
      aList = (int *) Malloc(sizeof(int)*Count, "MolCompare: *aList");
      bList = (int *) Malloc(sizeof(int)*Count, "MolCompare: *bList");
  
      // fill the lists
      aWalker = (**(molecule **)a).start;
      bWalker = (**(molecule **)b).start;
      Counter = 0;
      aCounter = 0;
      bCounter = 0;
      while ((aWalker->next != (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
        aWalker = aWalker->next;
        bWalker = bWalker->next;
        if (aWalker->GetTrueFather() == NULL)
          aList[Counter] = Count + (aCounter++);
        else
          aList[Counter] = aWalker->GetTrueFather()->nr;
        if (bWalker->GetTrueFather() == NULL)
          bList[Counter] = Count + (bCounter++);
        else
          bList[Counter] = bWalker->GetTrueFather()->nr;
        Counter++;
      }
      // check if AtomCount was for real
      flag = 0;
      if ((aWalker->next == (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
        flag = -1;
      } else {
        if ((aWalker->next != (**(molecule **)a).end) && (bWalker->next == (**(molecule **)b).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;
        }
      }
      Free((void **)&aList, "MolCompare: *aList");
      Free((void **)&bList, "MolCompare: *bList");
      return flag;
    }
  }
  return  -1;
};

/** Simple output of the pointers in ListOfMolecules.
 * \param *out output stream
 */
void MoleculeListClass::Output(ofstream *out)
{
  *out<< Verbose(1) << "MoleculeList: ";
  for (int i=0;i<NumberOfMolecules;i++)
    *out << ListOfMolecules[i] << "\t";
  *out << endl;
};


/** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
 * \param *out output stream for debugging
 * \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(ofstream *out, char *path, int *SortIndex)
{
  bool status = true;
  ofstream ForcesFile;
  stringstream line;
  atom *Walker = NULL;
  element *runner = NULL;

  // open file for the force factors
  *out << Verbose(1) << "Saving  force factors ... ";
  line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
  ForcesFile.open(line.str().c_str(), ios::out);
  if (ForcesFile != NULL) {
    //cout << Verbose(1) << "Final AtomicForcesList: ";
    //output << prefix << "Forces" << endl;
    for(int j=0;j<NumberOfMolecules;j++) {
      //if (TEList[j] != 0) {
      runner = ListOfMolecules[j]->elemente->start;
      while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element
        runner = runner->next;
        if (ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms
          Walker = ListOfMolecules[j]->start;
          while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element
            Walker = Walker->next;
            if (Walker->type->Z == runner->Z) {
              if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
                //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
                ForcesFile <<  SortIndex[Walker->GetTrueFather()->nr] << "\t";
                } else  // otherwise a -1 to indicate an added saturation hydrogen
                  ForcesFile << "-1\t";
              }
            }
          }
      }
      ForcesFile << endl;
    }
    ForcesFile.close();
    *out << Verbose(1) << "done." << endl;
  } else {
    status = false;
    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
  }
  ForcesFile.close();
  
  return status;
};

/** Writes a config file for each molecule in the given \a **FragmentList.
 * \param *out output stream for debugging
 * \param *configuration standard configuration to attach atoms in fragment molecule to.
 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
 * \return true - success (each file was written), false - something went wrong.
 */
bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
{
  ofstream outputFragment;
  char FragmentName[MAXSTRINGSIZE];
  char PathBackup[MAXSTRINGSIZE];
  bool result = true;
  bool intermediateResult = true;
  atom *Walker = NULL;
  vector BoxDimension;
  char *FragmentNumber;
  int FragmentCounter = 0;
  ofstream output;
  
  // store the fragments as config and as xyz
  for(int i=0;i<NumberOfMolecules;i++) {
    // save default path as it is changed for each fragment
    strcpy(PathBackup, configuration->GetDefaultPath());

    // correct periodic
    ListOfMolecules[i]->ScanForPeriodicCorrection(out);

    // output xyz file
    FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
    sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    outputFragment.open(FragmentName, ios::out);
    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as XYZ ...";
    if (intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment))
      *out << " done." << endl;
    else
      *out << " failed." << endl;
    result = result && intermediateResult;
    outputFragment.close();
    outputFragment.clear();

    *out << Verbose(2) << "Contained atoms: ";
    Walker = ListOfMolecules[i]->start;
    while (Walker->next != ListOfMolecules[i]->end) {
      Walker = Walker->next;
      *out << Walker->Name << " ";
    }
    *out << endl;
    // prepare output of config file
    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    outputFragment.open(FragmentName, ios::out);
    strcpy(PathBackup, configuration->configpath);
    sprintf(FragmentName, "%s/%s%s/", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    
    // center on edge
    ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
    ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
    int j = -1;
    for (int k=0;k<3;k++) {
      j += k+1;
      BoxDimension.x[k] = 5.;
      ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
    }
    ListOfMolecules[i]->Translate(&BoxDimension);

    // also calculate necessary orbitals
    ListOfMolecules[i]->CountElements();  // this is a bugfix, atoms should should actually be added correctly to this fragment 
    ListOfMolecules[i]->CalculateOrbitals(*configuration);
    
    // change path in config
    configuration->SetDefaultPath(FragmentName);
    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << " as config ...";
    if (intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i]))
      *out << " done." << endl;
    else
      *out << " failed." << endl;
    // restore old config
    configuration->SetDefaultPath(PathBackup);
    
    result = result && intermediateResult;
    outputFragment.close();
    outputFragment.clear();
    Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
  }
  cout << " done." << endl;
  
  // printing final number
  *out << "Final number of fragments: " << FragmentCounter << "." << endl;
      
  return result;
};

/******************************************* Class MoleculeLeafClass ************************************************/

/** Constructor for MoleculeLeafClass root leaf.
 * \param *Up Leaf on upper level
 * \param *PreviousLeaf NULL - We are the first leaf on this level, otherwise points to previous in list
 */
//MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *Up = NULL, MoleculeLeafClass *Previous = NULL)
MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL)
{
//  if (Up != NULL)
//    if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
//      Up->DownLeaf = this;
//  UpLeaf = Up;
//  DownLeaf = NULL;
  Leaf = NULL;
  previous = PreviousLeaf;
  if (previous != NULL) {
    MoleculeLeafClass *Walker = previous->next;
    previous->next = this;
    next = Walker;
  } else {
    next = NULL;
  }
};

/** Destructor for MoleculeLeafClass.
 */
MoleculeLeafClass::~MoleculeLeafClass()
{
//  if (DownLeaf != NULL) {// drop leaves further down
//    MoleculeLeafClass *Walker = DownLeaf;
//    MoleculeLeafClass *Next;
//    do {
//      Next = Walker->NextLeaf;
//      delete(Walker);
//      Walker = Next;
//    } while (Walker != NULL);
//    // Last Walker sets DownLeaf automatically to NULL
//  }
  // remove the leaf itself
  if (Leaf != NULL) {
    delete(Leaf);
    Leaf = NULL;
  }
  // remove this Leaf from level list
  if (previous != NULL)   
    previous->next = next;
//  } else { // we are first in list (connects to UpLeaf->DownLeaf)
//    if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
//      NextLeaf->UpLeaf = UpLeaf;  // either null as we are top level or the upleaf of the first node
//    if (UpLeaf != NULL)
//      UpLeaf->DownLeaf = NextLeaf;  // either null as we are only leaf or NextLeaf if we are just the first
//  }
//  UpLeaf = NULL;
  if (next != NULL) // are we last in list
    next->previous = previous;
  next = NULL;
  previous = NULL;
};

/** Adds \a molecule leaf to the tree.
 * \param *ptr ptr to molecule to be added
 * \param *Previous previous MoleculeLeafClass referencing level and which on the level
 * \return true - success, false - something went wrong
 */
bool MoleculeLeafClass::AddLeaf(molecule *ptr, MoleculeLeafClass *Previous)
{
  return false;
};

/** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule.
 * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL.
 * \param *out output stream for debugging
 * \param *reference reference molecule with the bond structure to be copied
 * \param &FragmentCounter Counter needed to address \a **ListOfLocalAtoms
 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled
 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 * \return true - success, false - faoilure
 */
bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
{
  atom *Walker = NULL, *OtherWalker = NULL;
  bond *Binder = NULL;
  bool status = true;
  int AtomNo;

  // fill ListOfLocalAtoms if NULL was given
  if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
    return false;
  }
  
  if (status) {
    *out << Verbose(1) << "Creating adjacency list for subgraph " << this << "." << endl;
    Walker = Leaf->start;
    while (Walker->next != Leaf->end) {
      Walker = Walker->next;
      AtomNo = Walker->father->nr;  // global id of the current walker
      for(int i=0;i<reference->NumberOfBondsPerAtom[AtomNo];i++) { // go through father's bonds and copy them all
        Binder = reference->ListOfBondsPerAtom[AtomNo][i];
        OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->father)->nr];    // local copy of current bond partner of walker
        if (OtherWalker != NULL) {
          if (OtherWalker->nr > Walker->nr)
          Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
        } else {
          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->father)->nr << "] is NULL!" << endl;
          status = false;
        }
      }
    }
    Leaf->CreateListOfBondsPerAtom(out);
    FragmentCounter++;
    if (next != NULL)
      status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
  }
  
  if (FreeList) {
    // free the index lookup list
    Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
    if (ListOfLocalAtoms[FragmentCounter] == NULL)
      Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
  }
  FragmentCounter--;
  return status;
};

/** Fills the root stack for sites to be used as root in fragmentation depending on order or adaptivity criteria
 * Again, as in \sa FillBondStructureFromReference steps recursively through each Leaf in this chain list of molecule's.
 * \param *out output stream for debugging
 * \param *&RootStack stack to be filled
 * \param Order desired bond order for all sites 
 * \param &FragmentCounter counts through the fragments in this MoleculeLeafClass
 * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
 */
bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, int Order, int &FragmentCounter)
{
  atom *Walker = NULL;

  if (RootStack != NULL) {
    if (Order == -1) {  // means we want to increase order adaptively
      cerr << "Adaptive chosing of sites not yet implemented!" << endl;
    } else {  // means we just want to increase it at every site by one
      // find first root candidates 
      if (&(RootStack[FragmentCounter]) != NULL) {
        RootStack[FragmentCounter].clear();  
        Walker = Leaf->start;
        while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
          Walker = Walker->next;
      #ifdef ADDHYDROGEN
          if (Walker->type->Z != 1) // skip hydrogen
      #endif
            if (Walker->GetTrueFather()->AdaptiveOrder < Order) // only if Order is still greater
              RootStack[FragmentCounter].push_front(Walker->nr);
        }
        if (next != NULL)
          next->FillRootStackForSubgraphs(out, RootStack, Order, ++FragmentCounter);
      }  else {
        *out << Verbose(1) << "Rootstack[" << FragmentCounter  << "] is NULL." << endl;
        return false;
      }
    }
    FragmentCounter--;
    return true;
  } else {
    *out << Verbose(1) << "Rootstack is NULL." << endl;
    return false;
  }
};

/** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
 * \param *out output stream fro debugging
 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
 * \param &FragmentCounter counts the fragments as we move along the list
 * \param GlobalAtomCount number of atoms in the complete molecule
 * \param &FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 * \return true - succes, false - failure
 */
bool MoleculeLeafClass::FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, int &FragmentCounter, int GlobalAtomCount, bool &FreeList)
{
  bool status = true;
  
  int Counter = Count();
  if (ListOfLocalAtoms == NULL) { // allocated initial pointer
    // allocate and set each field to NULL
    ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **)*Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
    if (ListOfLocalAtoms != NULL) {
      for (int i=0;i<Counter;i++)
        ListOfLocalAtoms[i] = NULL;
      FreeList = FreeList && true;
    } else
      status = false;
  }
  
  if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
    status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
    FreeList = FreeList && true;
  }
  
  return status;
};

/** The indices per keyset are compared to the respective father's Atom::nr in each subgraph and thus put into \a **&FragmentList.
 * \param *out output stream fro debugging
 * \param *reference reference molecule with the bond structure to be copied
 * \param *KeySetList list with all keysets
 * \param ***ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
 * \param **&FragmentList list to be allocated and returned
 * \param &FragmentCounter counts the fragments as we move along the list
 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 * \retuen true - success, false - failure
 */
bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList)
{
  bool status = true;
  int KeySetCounter = 0;
  
  // fill ListOfLocalAtoms if NULL was given
  if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
    return false;
  }

  // allocate fragment list
  if (FragmentList == NULL) {
    KeySetCounter = Count();
    FragmentList = (Graph **) Malloc(sizeof(Graph *)*KeySetCounter, "molecule::FragmentMolecule - **BondFragments");
    for(int i=0;i<KeySetCounter;i++)
      FragmentList[i] = NULL;
    KeySetCounter = 0;
  }
  
  if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
    // assign scanned keysets
    if (FragmentList[FragmentCounter] == NULL)
      FragmentList[FragmentCounter] = new Graph;
    KeySet *TempSet = new KeySet;
    for(Graph::iterator runner = KeySetList->begin();runner != KeySetList->end(); runner++) { // key sets contain global numbers!
      if ( ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr]->nr != -1) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
        // translate keyset to local numbers
        for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
          TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
        // insert into FragmentList
        FragmentList[FragmentCounter]->insert(GraphPair (*TempSet, pair<int,double>(KeySetCounter++, (*runner).second.second)));
      }
      TempSet->clear();
    }
    delete(TempSet);
    if (KeySetCounter == 0) {// if there are no keysets, delete the list
      *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
      delete(FragmentList[FragmentCounter]);
    } else
      *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
    FragmentCounter++;
    if (next != NULL)
      next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
    FragmentCounter--;
  } else
    *out << Verbose(1) << "KeySetList is NULL or empty." << endl;
  
  return status;
};

/** Simply counts the number of items in the list, from given MoleculeLeafClass.
 * \return number of items
 */
int MoleculeLeafClass::Count() const
{
  if (next != NULL)
    return next->Count()+1;
  else
    return 1; 
};