Changeset 3690e4 for src/Fragmentation/Interfragmenter.cpp

Timestamp:

Feb 2, 2016, 5:50:29 PM (9 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

9ae11c

Parents:

d1831e (diff), 62d092 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'Enhancing_Interdistance' into Candidate_v1.5.1

File:

• src/Fragmentation/Interfragmenter.cpp (modified) (7 diffs)

src/Fragmentation/Interfragmenter.cpp

@@ -37 +37 @@
 #include "Interfragmenter.hpp"
 
-#include <list>
-#include <map>
-
 #include "CodePatterns/Assert.hpp"
 #include "CodePatterns/Log.hpp"
@@ -45 +42 @@
 #include "LinearAlgebra/Vector.hpp"
 
-#include "AtomIdSet.hpp"
+#include "Descriptors/AtomDescriptor.hpp"
 #include "Element/element.hpp"
 #include "Fragmentation/Graph.hpp"
@@ -53 +50 @@
 #include "World.hpp"
 
-void Interfragmenter::operator()(
-    size_t MaxOrder,
-    double Rcut,
-    const enum HydrogenTreatment treatment)
+
+Interfragmenter::atomkeyset_t Interfragmenter::getAtomKeySetMap(
+    size_t _MaxOrder
+    ) const
 {
   /// 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.");
@@ -69 +64 @@
     const AtomIdSet atoms(keyset);
     const size_t atoms_size = atoms.getAtomIds().size();
-    if ((atoms_size > MaxOrder) || (atoms_size == 0))
+    if ((atoms_size > _MaxOrder) || (atoms_size == 0))
       continue;
     for (AtomIdSet::const_iterator atomiter = atoms.begin();
@@ -87 +82 @@
   LOG(2, "DEBUG: There are " << atomkeyset.size() << " entries in lookup.");
 
+  return atomkeyset;
+}
+
+static Vector getAtomIdSetCenter(
+    const AtomIdSet &_atoms)
+{
+  const molecule * const _mol = (*_atoms.begin())->getMolecule();
+  const size_t atoms_size = _atoms.getAtomIds().size();
+  Vector center;
+  for (AtomIdSet::const_iterator iter = _atoms.begin();
+      iter != _atoms.end(); ++iter) {
+    center += (*iter)->getPosition();
+    ASSERT ( _mol == (*iter)->getMolecule(),
+        "getAtomIdSetCenter() - ids in same keyset belong to different molecule.");
+  }
+  center *= 1./(double)atoms_size;
+
+  return center;
+}
+
+Interfragmenter::candidates_t Interfragmenter::getNeighborsOutsideMolecule(
+    const AtomIdSet &_atoms,
+    double _Rcut,
+    const enum HydrogenTreatment _treatment) const
+{
+  /// go through linked cell and get all neighboring atoms up to Rcut
+  const LinkedCell::LinkedCell_View view = World::getInstance().getLinkedCell(_Rcut);
+  const Vector center = getAtomIdSetCenter(_atoms);
+  const 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
+  const molecule * const _mol = (*_atoms.begin())->getMolecule();
+  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::getNeighborsOutsideMolecule() - 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.");
+
+  return candidates;
+}
+
+Interfragmenter::atomkeyset_t Interfragmenter::getCandidatesSpecificKeySetMap(
+    const candidates_t &_candidates,
+    const atomkeyset_t &_atomkeyset) const
+{
+  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::getAtomSpecificKeySetMap() - 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.");
+
+  return fragmentmap;
+}
+
+void Interfragmenter::combineFragments(
+    const size_t _MaxOrder,
+    const candidates_t &_candidates,
+    atomkeyset_t &_fragmentmap,
+    const KeySet &_keyset,
+    Graph &_InterFragments,
+    int &_counter)
+{
+  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::combineFragments() - could not find atom "
+        +toString(_atom->getId())+" in fragment specific lookup.");
+    keysets_t &othersets = finditer->second;
+    ASSERT( !othersets.empty(),
+        "Interfragmenter::combineFragments() - keysets to "+toString(_atom->getId())+
+        "is empty.");
+    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;
+      }
+      // only add if combined they don't exceed the desired maxorder
+      if (otherset.size() + _keyset.size() > _MaxOrder) {
+        LOG(3, "INFO: Rejecting " << otherset << " as in sum their orders exceed " << _MaxOrder);
+        ++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);
+    }
+  }
+}
+
+void Interfragmenter::operator()(
+    const size_t MaxOrder,
+    const double Rcut,
+    const enum HydrogenTreatment treatment)
+{
+  atomkeyset_t atomkeyset = getAtomKeySetMap(MaxOrder);
+
   Graph InterFragments;
   int counter = TotalGraph.size();
@@ -101 +220 @@
       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.");
+    // get neighboring atoms outside the current molecule
+    candidates_t candidates = getNeighborsOutsideMolecule(atoms, Rcut, treatment);
 
     // 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.");
+    atomkeyset_t fragmentmap = getCandidatesSpecificKeySetMap(candidates, atomkeyset);
 
     /// 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);
-      }
-    }
+    /// if keyset is less (to prevent addding same inter-fragment twice)
+    combineFragments(MaxOrder, candidates, fragmentmap, keyset, InterFragments, counter);
   }
 
@@ -181 +235 @@
   TotalGraph.InsertGraph(InterFragments, counter);
 }
+
+double Interfragmenter::findLargestCutoff(
+    const size_t _MaxOrder,
+    const double _upperbound,
+    const enum HydrogenTreatment _treatment) const
+{
+  double Rcut = _upperbound*_upperbound;
+  std::pair<atomId_t, atomId_t> ClosestPair;
+  // place all atoms into LC grid with some upper bound
+  const LinkedCell::LinkedCell_View view = World::getInstance().getLinkedCell(_upperbound);
+
+  atomkeyset_t atomkeyset = getAtomKeySetMap(_MaxOrder);
+
+  // go through each atom and find closest atom not in the same keyset
+  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);
+
+    // get neighboring atoms outside the current molecule
+    const candidates_t candidates = getNeighborsOutsideMolecule(atoms, _upperbound, _treatment);
+    const Vector center = getAtomIdSetCenter(atoms);
+
+    for (candidates_t::const_iterator candidateiter = candidates.begin();
+        candidateiter != candidates.end(); ++candidateiter) {
+      atom const * const Walker = *candidateiter;
+      // go through each atom in set and pick minimal distance
+      for (AtomIdSet::const_iterator setiter = atoms.begin(); setiter != atoms.end(); ++setiter) {
+        const double distanceSquared = Walker->getPosition().DistanceSquared(center);
+        // pick the smallest compared to current Rcut if smaller
+        if (distanceSquared < Rcut) {
+          Rcut = distanceSquared;
+          ClosestPair.first = (*setiter)->getId();
+          ClosestPair.second = Walker->getId();
+          LOG(2, "DEBUG: Found new pair " << ClosestPair << " with distance " << sqrt(Rcut));
+        }
+      }
+    }
+  }
+  const double largest_distance = sqrt(Rcut);
+  LOG(1, "INFO: Largest inter-fragment distance to cause no additional fragments: "
+      << largest_distance);
+
+  return largest_distance;
+}