source: src/Actions/FragmentationAction/CreateAdjacencyAction.cpp@ 3738f0

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

/*
 * CreateAdjacencyAction.cpp
 *
 *  Created on: May 9, 2010
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "Descriptors/AtomIdDescriptor.hpp"
#include "Descriptors/MoleculeDescriptor.hpp"

#include "atom.hpp"
#include "bond.hpp"
#include "bondgraph.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "config.hpp"
#include "linkedcell.hpp"
#include "molecule.hpp"
#include "PointCloudAdaptor.hpp"
#include "World.hpp"
#include "WorldTime.hpp"

#include <iostream>
#include <list>
#include <string>

typedef std::map< moleculeId_t, std::vector<atomId_t> > MolAtomList;

using namespace std;

#include "Actions/FragmentationAction/CreateAdjacencyAction.hpp"

// and construct the stuff
#include "CreateAdjacencyAction.def"
#include "Action_impl_pre.hpp"
/** =========== define the function ====================== */
Action::state_ptr FragmentationCreateAdjacencyAction::performCall() {
  // obtain information
  getParametersfromValueStorage();

  DoLog(1) && (Log() << Verbose(1) << "Constructing bond graph for selected atoms ... " << endl);

  config *configuration = World::getInstance().getConfig();
  BondGraph *BG = World::getInstance().getBondGraph();
  ASSERT(BG != NULL, "FragmentationCreateAdjacencyAction: BondGraph is NULL.");
  BG->SetMaxDistanceToMaxOfCovalentRadii(AtomSetMixin<std::vector<atom *> >(World::getInstance().getSelectedAtoms()));
  double BondDistance = BG->getMaxDistance();
  bool IsAngstroem = configuration->GetIsAngstroem();

  atom *Walker = NULL;
  atom *OtherWalker = NULL;
  int n[NDIM];
  double MinDistance, MaxDistance;
  LinkedCell *LC = NULL;
  Box &domain = World::getInstance().getDomain();

  // remove every bond from the selected atoms' list
  int AtomCount = 0;
  for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection();
      AtomRunner != World::getInstance().endAtomSelection();
      ++AtomRunner) {
    AtomCount++;
    BondList& ListOfBonds = (AtomRunner->second)->getListOfBonds();
    for(BondList::iterator BondRunner = ListOfBonds.begin();
        !ListOfBonds.empty();
        BondRunner = ListOfBonds.begin())
      if ((*BondRunner)->leftatom == AtomRunner->second)
        delete((*BondRunner));
  }
  int BondCount = 0;

  // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
  DoLog(1) && (Log() << Verbose(1) << "AtomCount " << AtomCount << " and bonddistance is " << BondDistance << "." << endl);

  if ((AtomCount > 1) && (BondDistance > 1.)) {
    DoLog(2) && (Log() << Verbose(2) << "Creating Linked Cell structure ... " << endl);
    TesselPointSTLList list;
    for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection();
        AtomRunner != World::getInstance().endAtomSelection();
        ++AtomRunner) {
      list.push_back(AtomRunner->second);
    }
    PointCloudAdaptor< TesselPointSTLList > cloud(&list, "AtomSelection");
    LC = new LinkedCell(cloud, BondDistance);

    // create a list to map Tesselpoint::Nr to atom *
    DoLog(2) && (Log() << Verbose(2) << "Creating TesselPoint to atom map ... " << endl);

    // set numbers for atoms that can later be used
    std::map<TesselPoint *, int> AtomIds;
    int i=0;
    for (World::AtomSelectionIterator AtomRunner = World::getInstance().beginAtomSelection();
        AtomRunner != World::getInstance().endAtomSelection();
        ++AtomRunner) {
      AtomIds.insert(pair<TesselPoint *, int> (AtomRunner->second, i++) );
    }

    // 3a. go through every cell
    DoLog(2) && (Log() << Verbose(2) << "Celling ... " << endl);
    for (LC->n[0] = 0; LC->n[0] < LC->N[0]; LC->n[0]++)
      for (LC->n[1] = 0; LC->n[1] < LC->N[1]; LC->n[1]++)
        for (LC->n[2] = 0; LC->n[2] < LC->N[2]; LC->n[2]++) {
          const TesselPointSTLList *List = LC->GetCurrentCell();
//          Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
          if (List != NULL) {
            for (TesselPointSTLList::const_iterator Runner = List->begin();
                Runner != List->end();
                Runner++) {
              Walker = dynamic_cast<atom*>(*Runner);
              ASSERT(Walker,"Tesselpoint that was not an atom retrieved from LinkedNode");
              //Log() << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
              // 3c. check for possible bond between each atom in this and every one in the 27 cells
              for (n[0] = -1; n[0] <= 1; n[0]++)
                for (n[1] = -1; n[1] <= 1; n[1]++)
                  for (n[2] = -1; n[2] <= 1; n[2]++) {
                    const TesselPointSTLList *OtherList = LC->GetRelativeToCurrentCell(n);
//                    Log() << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
                    if (OtherList != NULL) {
                      for (TesselPointSTLList::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
                        if (AtomIds.find(*OtherRunner)->second > AtomIds.find(Walker)->second) {
                          OtherWalker = dynamic_cast<atom*>(*OtherRunner);
                          ASSERT(OtherWalker,"TesselPoint that was not an atom retrieved from LinkedNode");
                          //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
                          BG->getMinMaxDistance(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
                          const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
                          const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
//                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
                          if (AtomIds[OtherWalker->father] > AtomIds[Walker->father]) {
                            if (status) { // create bond if distance is smaller
//                              Log() << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
                              bond * const Binder = new bond(Walker->father, OtherWalker->father, 1, BondCount++);
                              Walker->father->RegisterBond(WorldTime::getTime(),Binder);
                              OtherWalker->father->RegisterBond(WorldTime::getTime(),Binder);
                            } else {
//                              Log() << Verbose(1) << "Not Adding: distance too great." << endl;
                            }
                          } else {
//                            Log() << Verbose(1) << "Not Adding: Wrong order of labels." << endl;
                          }
                        }
                      }
                    }
                  }
            }
          }
        }
    delete (LC);
    DoLog(1) && (Log() << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl);

    // correct bond degree by comparing valence and bond degree
    DoLog(2) && (Log() << Verbose(2) << "Correcting bond degree ... " << endl);
    //CorrectBondDegree();

  } else
    DoLog(1) && (Log() << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl);
  DoLog(0) && (Log() << Verbose(0) << "End of CreateAdjacencyList." << endl);

  return Action::success;
}

Action::state_ptr FragmentationCreateAdjacencyAction::performUndo(Action::state_ptr _state) {
//  FragmentationCreateAdjacencyState *state = assert_cast<FragmentationCreateAdjacencyState*>(_state.get());

  return Action::success;
}

Action::state_ptr FragmentationCreateAdjacencyAction::performRedo(Action::state_ptr _state){
  return Action::success;
}

bool FragmentationCreateAdjacencyAction::canUndo() {
  return false;
}

bool FragmentationCreateAdjacencyAction::shouldUndo() {
  return false;
}
/** =========== end of function ====================== */