source: src/Fragmentation/Interfragmenter.cpp@ c6719d

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * 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/>.
 */

/*
 * Interfragmenter.cpp
 *
 *  Created on: Jul 5, 2013
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "Interfragmenter.hpp"

#include <list>
#include <map>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"

#include "LinearAlgebra/Vector.hpp"

#include "AtomIdSet.hpp"
#include "Element/element.hpp"
#include "Fragmentation/Graph.hpp"
#include "Fragmentation/KeySet.hpp"
#include "LinkedCell/LinkedCell_View.hpp"
#include "LinkedCell/types.hpp"
#include "World.hpp"

void Interfragmenter::operator()(
    size_t MaxOrder,
    double Rcut,
    const enum HydrogenTreatment treatment)
{
  /// create a map of atom to keyset (below equal MaxOrder)
  typedef std::list<const KeySet *> keysets_t;
  typedef std::map<const atom *, keysets_t > atomkeyset_t;
  atomkeyset_t atomkeyset;
  LOG(1, "INFO: Placing all atoms and their keysets into a map.");
  for (Graph::const_iterator keysetiter = TotalGraph.begin();
      keysetiter != TotalGraph.end(); ++keysetiter) {
    const KeySet &keyset = keysetiter->first;
    LOG(2, "DEBUG: Current keyset is " << keyset);
    const AtomIdSet atoms(keyset);
    const size_t atoms_size = atoms.getAtomIds().size();
    if ((atoms_size > MaxOrder) || (atoms_size == 0))
      continue;
    for (AtomIdSet::const_iterator atomiter = atoms.begin();
        atomiter != atoms.end(); ++atomiter) {
      // either create new list ...
      std::pair<atomkeyset_t::iterator, bool> inserter =
          atomkeyset.insert( std::make_pair(*atomiter, keysets_t(1, &keyset) ));
      // ... or push to end
      if (inserter.second) {
        LOG(3, "DEBUG: Created new entry in map.");
      } else {
        LOG(3, "DEBUG: Added keyset to present entry.");
        inserter.first->second.push_back(&keyset);
      }
    }
  }
  LOG(2, "DEBUG: There are " << atomkeyset.size() << " entries in lookup.");

  Graph InterFragments;
  int counter = TotalGraph.size();

  /// go through all fragments up to MaxOrder
  LOG(1, "INFO: Creating inter-fragments.");
  for (Graph::const_iterator keysetiter = TotalGraph.begin();
      keysetiter != TotalGraph.end(); ++keysetiter) {
    const KeySet &keyset = keysetiter->first;
    LOG(2, "DEBUG: Current keyset is " << keyset);
    const AtomIdSet atoms(keyset);
    const size_t atoms_size = atoms.getAtomIds().size();
    if ((atoms_size > MaxOrder) || (atoms_size == 0))
      continue;

    /// go through linked cell and get all neighboring atoms up to Rcut
    Vector center;
    const molecule *_mol = (*atoms.begin())->getMolecule();
    for (AtomIdSet::const_iterator iter = atoms.begin();
        iter != atoms.end(); ++iter) {
      center += (*iter)->getPosition();
      ASSERT ( _mol == (*iter)->getMolecule(),
          "Interfragmenter::operator() - ids in same keyset belong to different molecule.");
    }
    center *= 1./(double)atoms_size;
    LinkedCell::LinkedCell_View view = World::getInstance().getLinkedCell(Rcut);
    LinkedCell::LinkedList neighbors = view.getAllNeighbors(Rcut, center);
    LOG(4, "DEBUG: Obtained " << neighbors.size() << " neighbors in distance of "
        << Rcut << " around " << center << ".");

    /// remove all atoms that belong to same molecule as does one of the
    /// fragment's atoms
    typedef std::vector<const atom *> candidates_t;
    candidates_t candidates;
    candidates.reserve(neighbors.size());
    for (LinkedCell::LinkedList::const_iterator iter = neighbors.begin();
        iter != neighbors.end(); ++iter) {
      const atom * const _atom = static_cast<const atom * const >(*iter);
      ASSERT( _atom != NULL,
          "Interfragmenter::operator() - a neighbor is not actually an atom?");
      if ((_atom->getMolecule() != _mol)
          && (_atom->getPosition().DistanceSquared(center) < Rcut*Rcut)
          && ((treatment == IncludeHydrogen) || (_atom->getType()->getAtomicNumber() != 1))) {
        candidates.push_back(_atom);
      }
    }
    LOG(3, "DEBUG: There remain " << candidates.size() << " candidates.");

    // create a lookup specific to this fragment
    atomkeyset_t fragmentmap;
    for (candidates_t::const_iterator candidateiter = candidates.begin();
        candidateiter != candidates.end(); ++candidateiter) {
      const atom * _atom = *candidateiter;
      atomkeyset_t::const_iterator iter = atomkeyset.find(_atom);
      ASSERT( iter != atomkeyset.end(),
          "Interfragmenter::operator() - could not find atom "
          +toString(_atom->getId())+" in lookup.");
      fragmentmap.insert( std::make_pair( _atom, iter->second) );
    }
    LOG(4, "DEBUG: Copied part of lookup map contains " << fragmentmap.size() << " keys.");

    /// combine each remaining fragment with current fragment to a new fragment
    /// if keyset is less
    for (candidates_t::const_iterator candidateiter = candidates.begin();
        candidateiter != candidates.end(); ++candidateiter) {
      const atom *_atom = *candidateiter;
      LOG(3, "DEBUG: Current candidate is " << *_atom << ".");
      atomkeyset_t::iterator finditer = fragmentmap.find(_atom);
      ASSERT( finditer != fragmentmap.end(),
          "Interfragmenter::operator() - could not find atom "
          +toString(_atom->getId())+" in fragment specific lookup.");
      keysets_t &othersets = finditer->second;
      keysets_t::iterator otheriter = othersets.begin();
      while (otheriter != othersets.end()) {
        const KeySet &otherset = **otheriter;
        LOG(3, "DEBUG: Current keyset is " << otherset << ".");
        // only add them one way round and not the other
        if (otherset < keyset) {
          ++otheriter;
          continue;
        }
        KeySet newset(otherset);
        newset.insert(keyset.begin(), keyset.end());
        LOG(3, "DEBUG: Inserting new combined set " << newset << ".");
        InterFragments.insert( std::make_pair(newset, std::make_pair(++counter, 1.)));
        // finally, remove the set such that no other combination exists
        otheriter = othersets.erase(otheriter);
      }
    }
  }

  /// eventually, add all new fragments to the Graph
  counter = TotalGraph.size();
  TotalGraph.InsertGraph(InterFragments, counter);
}