Changeset c1ec8e for src/Actions

Timestamp:

May 18, 2016, 10:02:53 PM (9 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, 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_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, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, 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:

91f907

Parents:

01120c

git-author:

Frederik Heber <heber@…> (03/07/16 21:03:37)

git-committer:

Frederik Heber <heber@…> (05/18/16 22:02:53)

Message:

Split FitPartialChargesAction's extensive performCall() into subfunctions.

File:

• src/Actions/PotentialAction/FitPartialChargesAction.cpp (modified) (5 diffs)

src/Actions/PotentialAction/FitPartialChargesAction.cpp

@@ -171 +171 @@
 }
 
-ActionState::ptr PotentialFitPartialChargesAction::performCall()
-{
-  // check for selected atoms
-  const World &world = const_cast<const World &>(World::getInstance());
-  if (world.beginAtomSelection() == world.endAtomSelection()) {
-    STATUS("There are no atoms selected for fitting partial charges to.");
-    return Action::failure;
-  }
-
-  /// obtain possible fragments to each selected atom
-  const AtomFragmentsMap &atomfragments = AtomFragmentsMap::getConstInstance();
-  if (!atomfragments.checkCompleteness()) {
-    STATUS("AtomFragmentsMap failed internal consistency check, missing forcekeysets?");
-    return Action::failure;
-  }
+static
+std::set<KeySet> accumulateKeySetsForAtoms(
+    const AtomFragmentsMap::AtomFragmentsMap_t &_atommap,
+    const std::vector<const atom *> &_selected_atoms)
+{
   std::set<KeySet> fragments;
-  const AtomFragmentsMap::AtomFragmentsMap_t &atommap = atomfragments.getMap();
-  for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
-      iter != world.endAtomSelection(); ++iter) {
+  for (std::vector<const atom *>::const_iterator iter = _selected_atoms.begin();
+      iter != _selected_atoms.end(); ++iter) {
+    const atomId_t atomid = (*iter)->getId();
     const AtomFragmentsMap::AtomFragmentsMap_t::const_iterator atomiter =
-        atommap.find(iter->first);
-    if (iter->second->getElementNo() != 1) {
-      if (atomiter == atommap.end()) {
-        ELOG(2, "There are no fragments associated to " << iter->first << ".");
+        _atommap.find(atomid);
+    if ((*iter)->getElementNo() != 1) {
+      if (atomiter == _atommap.end()) {
+        ELOG(2, "There are no fragments associated to " << atomid << ".");
         continue;
       }
       const AtomFragmentsMap::keysets_t &keysets = atomiter->second;
-      LOG(2, "DEBUG: atom " << iter->first << " has " << keysets.size() << " fragments.");
+      LOG(2, "DEBUG: atom " << atomid << " has " << keysets.size() << " fragments.");
       fragments.insert( keysets.begin(), keysets.end() );
     } else {
-      LOG(3, "DEBUG: Skipping atom " << iter->first << " as it's hydrogen.");
-    }
-  }
-
-  // reduce given fragments to homologous graphs to avoid multiple fittings
-  typedef std::map<KeySet, HomologyGraph> KeysetsToGraph_t;
-  KeysetsToGraph_t keyset_graphs;
-  typedef std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> GraphFittedChargeMap_t;
-  GraphFittedChargeMap_t fittedcharges_per_fragment;
-  HomologyContainer &homologies = World::getInstance().getHomologies();
-  for (std::set<KeySet>::const_iterator fragmentiter = fragments.begin();
-      fragmentiter != fragments.end(); ++fragmentiter) {
+      LOG(3, "DEBUG: Skipping atom " << atomid << " as it's hydrogen.");
+    }
+  }
+  return fragments;
+}
+
+static
+void getKeySetsToGraphMapping(
+    std::map<KeySet, HomologyGraph> &_keyset_graphs,
+    std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> &_fittedcharges_per_fragment,
+    const std::set<KeySet> &_fragments,
+    const AtomFragmentsMap &_atomfragments)
+{
+  for (std::set<KeySet>::const_iterator fragmentiter = _fragments.begin();
+      fragmentiter != _fragments.end(); ++fragmentiter) {
     const KeySet &keyset = *fragmentiter;
-    const AtomFragmentsMap::indices_t &forceindices = atomfragments.getFullKeyset(keyset);
+    const AtomFragmentsMap::indices_t &forceindices = _atomfragments.getFullKeyset(keyset);
     ASSERT( !forceindices.empty(),
-        "PotentialFitPartialChargesAction::performCall() - force keyset to "
-        +toString(keyset)+" is empty.");
+        "getKeySetsToGraphMapping() - force keyset to "+toString(keyset)+" is empty.");
     KeySet forcekeyset;
     forcekeyset.insert(forceindices.begin(), forceindices.end());
@@ -223 +215 @@
     const HomologyGraph graph(forcekeyset);
     LOG(2, "DEBUG: Associating keyset " << forcekeyset << " with graph " << graph);
-    keyset_graphs.insert( std::make_pair(keyset, graph) );
-    fittedcharges_per_fragment.insert( std::make_pair(graph, PartialNucleiChargeFitter::charges_t()) );
-  }
-
-  /// then go through all fragments and get partial charges for each
-  for (GraphFittedChargeMap_t::iterator graphiter = fittedcharges_per_fragment.begin();
-      graphiter != fittedcharges_per_fragment.end(); ++graphiter) {
+    _keyset_graphs.insert( std::make_pair(keyset, graph) );
+    _fittedcharges_per_fragment.insert( std::make_pair(graph, PartialNucleiChargeFitter::charges_t()) );
+  }
+}
+
+static
+bool getPartialChargesForAllGraphs(
+    std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> &_fittedcharges_per_fragment,
+    const HomologyContainer &_homologies,
+    const double _mask_radius,
+    const bool enforceZeroCharge)
+{
+  typedef std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> GraphFittedChargeMap_t;
+  for (GraphFittedChargeMap_t::iterator graphiter = _fittedcharges_per_fragment.begin();
+      graphiter != _fittedcharges_per_fragment.end(); ++graphiter) {
     const HomologyGraph &graph = graphiter->first;
     LOG(2, "DEBUG: Now fitting charges to graph " << graph);
-    const HomologyContainer::range_t range = homologies.getHomologousGraphs(graph);
+    const HomologyContainer::range_t range = _homologies.getHomologousGraphs(graph);
     if (range.first == range.second) {
-      STATUS("HomologyContainer does not contain specified fragment.");
-      return Action::failure;
+      ELOG(0, "HomologyContainer does not contain specified fragment.");
+      return false;
     }
 
@@ -241 +241 @@
     PartialNucleiChargeFitter::charges_t &averaged_charges = graphiter->second;
     const size_t NoFragments =
-        obtainAverageChargesOverFragments(averaged_charges, range, params.radius.get());
+        obtainAverageChargesOverFragments(averaged_charges, range, _mask_radius);
     if ((NoFragments == 0) && (range.first != range.second)) {
-      STATUS("Fitting for fragment "+toString(*graphiter)+" failed.");
-      return Action::failure;
+      ELOG(0, "Fitting for fragment "+toString(*graphiter)+" failed.");
+      return false;
     }
 
     // make sum of charges zero if desired
-    if (params.enforceZeroCharge.get())
+    if (enforceZeroCharge)
       enforceZeroTotalCharge(averaged_charges);
 
@@ -255 +255 @@
         << averaged_charges << ", averaged over " << NoFragments << " fragments.");
   }
-
-  /// obtain average charge for each atom the fitted charges over all its fragments
+  return true;
+}
+
+double fitAverageChargeToAtom(
+    const atom * const _walker,
+    const AtomFragmentsMap &_atomfragments,
+    const std::map<KeySet, HomologyGraph> &_keyset_graphs,
+    const std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> &_fittedcharges_per_fragment)
+{
+  typedef std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> GraphFittedChargeMap_t;
+  const atom * surrogate = _walker;
+  if (surrogate->getElementNo() == 1) {
+    // it's hydrogen, check its bonding and use its bond partner instead to request
+    // keysets
+    const BondList &ListOfBonds = surrogate->getListOfBonds();
+    if ( ListOfBonds.size() != 1) {
+      ELOG(1, "Solitary hydrogen in atom " << surrogate->getId() << " detected, skipping.");
+      return 0.;
+    }
+    surrogate = (*ListOfBonds.begin())->GetOtherAtom(surrogate);
+  }
+  const atomId_t walkerid = surrogate->getId();
+  const AtomFragmentsMap::AtomFragmentsMap_t &atommap = _atomfragments.getMap();
+  const AtomFragmentsMap::AtomFragmentsMap_t::const_iterator keysetsiter =
+      atommap.find(walkerid);
+  ASSERT(keysetsiter != atommap.end(),
+      "fitAverageChargeToAtom() - we checked already that "+toString(walkerid)
+      +" should be present!");
+  const AtomFragmentsMap::keysets_t & keysets = keysetsiter->second;
+
+  double average_charge = 0.;
+  size_t NoFragments = 0;
+  // go over all fragments associated to this atom
+  for (AtomFragmentsMap::keysets_t::const_iterator keysetsiter = keysets.begin();
+      keysetsiter != keysets.end(); ++keysetsiter) {
+    const KeySet &keyset = *keysetsiter;
+
+    const AtomFragmentsMap::indices_t &forcekeyset = _atomfragments.getFullKeyset(keyset);
+    ASSERT( !forcekeyset.empty(),
+        "fitAverageChargeToAtom() - force keyset to "+toString(keyset)+" is empty.");
+
+    // find the associated charge in the charge vector
+    std::map<KeySet, HomologyGraph>::const_iterator keysetgraphiter =
+        _keyset_graphs.find(keyset);
+    ASSERT( keysetgraphiter != _keyset_graphs.end(),
+        "fitAverageChargeToAtom() - keyset "+toString(keyset)
+        +" not contained in keyset_graphs.");
+    const HomologyGraph &graph = keysetgraphiter->second;
+    const GraphFittedChargeMap_t::const_iterator chargesiter =
+        _fittedcharges_per_fragment.find(graph);
+    ASSERT(chargesiter != _fittedcharges_per_fragment.end(),
+        "fitAverageChargeToAtom() - no charge to "+toString(keyset)
+        +" any longer present in fittedcharges_per_fragment?");
+    const PartialNucleiChargeFitter::charges_t &charges = chargesiter->second;
+    ASSERT( charges.size() == forcekeyset.size(),
+        "fitAverageChargeToAtom() - charges "+toString(charges.size())+" and keyset "
+        +toString(forcekeyset.size())+" do not have the same length?");
+    PartialNucleiChargeFitter::charges_t::const_iterator chargeiter =
+        charges.begin();
+    const AtomFragmentsMap::indices_t::const_iterator forcekeysetiter =
+        std::find(forcekeyset.begin(), forcekeyset.end(), _walker->getId());
+    ASSERT( forcekeysetiter != forcekeyset.end(),
+        "fitAverageChargeToAtom() - atom "+toString(_walker->getId())
+        +" not contained in force keyset "+toString(forcekeyset));
+    std::advance(chargeiter, std::distance(forcekeyset.begin(), forcekeysetiter));
+
+    // and add onto charge sum
+    const double & charge_in_fragment = *chargeiter;
+    average_charge += charge_in_fragment;
+    ++NoFragments;
+  }
+  // average to obtain final partial charge for this atom
+  average_charge *= 1./(double)NoFragments;
+
+  return average_charge;
+}
+
+void addToParticleRegistry(
+    const ParticleFactory &factory,
+    const periodentafel &periode,
+    const std::map<atomId_t, double> &_fitted_charges)
+{
   typedef std::map<atomId_t, double> fitted_charges_t;
-  fitted_charges_t fitted_charges;
-  for (World::AtomSelectionConstIterator atomiter = world.beginAtomSelection();
-      atomiter != world.endAtomSelection(); ++atomiter) {
-    const atom * walker = atomiter->second;
-    if (walker->getElementNo() == 1) {
-      // it's hydrogen, check its bonding and use its bond partner instead to request
-      // keysets
-      const BondList &ListOfBonds = walker->getListOfBonds();
-      if ( ListOfBonds.size() != 1) {
-        ELOG(1, "Solitary hydrogen in atom " << atomiter->first << " detected, skipping.");
-        continue;
-      }
-      walker = (*ListOfBonds.begin())->GetOtherAtom(walker);
-    }
-    const atomId_t walkerid = walker->getId();
-    const AtomFragmentsMap::AtomFragmentsMap_t::const_iterator keysetsiter =
-        atommap.find(walkerid);
-    ASSERT(keysetsiter != atommap.end(),
-        "PotentialFitPartialChargesAction::performCall() - we checked already that "
-        +toString(walkerid)+" should be present!");
-    const AtomFragmentsMap::keysets_t & keysets = keysetsiter->second;
-
-    double average_charge = 0.;
-    size_t NoFragments = 0;
-    // go over all fragments associated to this atom
-    for (AtomFragmentsMap::keysets_t::const_iterator keysetsiter = keysets.begin();
-        keysetsiter != keysets.end(); ++keysetsiter) {
-      const KeySet &keyset = *keysetsiter;
-
-      const AtomFragmentsMap::indices_t &forcekeyset = atomfragments.getFullKeyset(keyset);
-      ASSERT( !forcekeyset.empty(),
-          "PotentialFitPartialChargesAction::performCall() - force keyset to "
-          +toString(keyset)+" is empty.");
-
-      // find the associated charge in the charge vector
-      const HomologyGraph &graph = keyset_graphs[keyset];
-      const GraphFittedChargeMap_t::const_iterator chargesiter =
-          fittedcharges_per_fragment.find(graph);
-      ASSERT(chargesiter != fittedcharges_per_fragment.end(),
-          "PotentialFitPartialChargesAction::performCall() - no charge to "+toString(keyset)
-          +" any longer present in fittedcharges_per_fragment?");
-      const PartialNucleiChargeFitter::charges_t &charges = chargesiter->second;
-      ASSERT( charges.size() == forcekeyset.size(),
-          "PotentialFitPartialChargesAction::performCall() - charges "
-          +toString(charges.size())+" and keyset "+toString(forcekeyset.size())
-          +" do not have the same length?");
-      PartialNucleiChargeFitter::charges_t::const_iterator chargeiter =
-          charges.begin();
-      const AtomFragmentsMap::indices_t::const_iterator forcekeysetiter =
-          std::find(forcekeyset.begin(), forcekeyset.end(), atomiter->first);
-      ASSERT( forcekeysetiter != forcekeyset.end(),
-          "PotentialFitPartialChargesAction::performCall() - atom "
-          +toString(atomiter->first)+" not contained in force keyset "+toString(forcekeyset));
-      std::advance(chargeiter, std::distance(forcekeyset.begin(), forcekeysetiter));
-
-      // and add onto charge sum
-      const double & charge_in_fragment = *chargeiter;
-      average_charge += charge_in_fragment;
-      ++NoFragments;
-    }
-    // average to obtain final partial charge for this atom
-    average_charge *= 1./(double)NoFragments;
-    fitted_charges.insert( std::make_pair(atomiter->first, average_charge) );
-
-    LOG(2, "DEBUG: For atom " << atomiter->first << " we have an average charge of "
-        << average_charge);
-  }
-
-  /// place all fitted charges into ParticleRegistry
-  const ParticleFactory &factory =
-      const_cast<const ParticleFactory&>(ParticleFactory::getInstance());
-  const periodentafel &periode = *World::getInstance().getPeriode();
-  for (fitted_charges_t::const_iterator chargeiter = fitted_charges.begin();
-      chargeiter != fitted_charges.end(); ++chargeiter) {
+  for (fitted_charges_t::const_iterator chargeiter = _fitted_charges.begin();
+      chargeiter != _fitted_charges.end(); ++chargeiter) {
     const atom * const walker = World::getInstance().getAtom(AtomById(chargeiter->first));
     ASSERT( walker != NULL,
@@ -343 +351 @@
     factory.createInstance(name, elementno, charge);
   }
+}
+
+ActionState::ptr PotentialFitPartialChargesAction::performCall()
+{
+  // check for selected atoms
+  const World &world = const_cast<const World &>(World::getInstance());
+  if (world.beginAtomSelection() == world.endAtomSelection()) {
+    STATUS("There are no atoms selected for fitting partial charges to.");
+    return Action::failure;
+  }
+
+  /// obtain possible fragments to each selected atom
+  const AtomFragmentsMap &atomfragments = AtomFragmentsMap::getConstInstance();
+  if (!atomfragments.checkCompleteness()) {
+    STATUS("AtomFragmentsMap failed internal consistency check, missing forcekeysets?");
+    return Action::failure;
+  }
+  std::set<KeySet> fragments =
+      accumulateKeySetsForAtoms( atomfragments.getMap(), world.getSelectedAtoms());
+
+  // reduce given fragments to homologous graphs to avoid multiple fittings
+  typedef std::map<KeySet, HomologyGraph> KeysetsToGraph_t;
+  KeysetsToGraph_t keyset_graphs;
+  typedef std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> GraphFittedChargeMap_t;
+  GraphFittedChargeMap_t fittedcharges_per_fragment;
+  getKeySetsToGraphMapping(keyset_graphs, fittedcharges_per_fragment, fragments, atomfragments);
+
+  /// then go through all fragments and get partial charges for each
+  const bool status = getPartialChargesForAllGraphs(
+      fittedcharges_per_fragment,
+      World::getInstance().getHomologies(),
+      params.radius.get(),
+      params.enforceZeroCharge.get());
+  if (!status)
+    return Action::failure;
+
+  /// obtain average charge for each atom the fitted charges over all its fragments
+  typedef std::map<atomId_t, double> fitted_charges_t;
+  fitted_charges_t fitted_charges;
+  for (World::AtomSelectionConstIterator atomiter = world.beginAtomSelection();
+      atomiter != world.endAtomSelection(); ++atomiter) {
+    const double average_charge = fitAverageChargeToAtom(
+        atomiter->second, atomfragments, keyset_graphs, fittedcharges_per_fragment);
+
+    if (average_charge != 0.) {
+      LOG(2, "DEBUG: For atom " << atomiter->first << " we have an average charge of "
+          << average_charge);
+
+      fitted_charges.insert( std::make_pair(atomiter->first, average_charge) );
+    }
+  }
+
+  /// place all fitted charges into ParticleRegistry
+  addToParticleRegistry(
+      ParticleFactory::getConstInstance(),
+      *World::getInstance().getPeriode(),
+      fitted_charges);
 
   return Action::success;