Index: src/Actions/PotentialAction/FitPartialChargesAction.cpp
===================================================================
--- src/Actions/PotentialAction/FitPartialChargesAction.cpp	(revision d8ed6292b835bf1c1008f7aba846e47d5a612ea4)
+++ src/Actions/PotentialAction/FitPartialChargesAction.cpp	(revision 01120c6d716ee08fe2e1ee35a1c1daa426350b9d)
@@ -44,4 +44,5 @@
 #include "World.hpp"
 
+#include <boost/bimap.hpp>
 #include <boost/bind.hpp>
 #include <boost/filesystem.hpp>
@@ -56,4 +57,6 @@
 #include "Potentials/PartialNucleiChargeFitter.hpp"
 
+#include "AtomIdSet.hpp"
+#include "Descriptors/AtomIdDescriptor.hpp"
 #include "Element/element.hpp"
 #include "Element/periodentafel.hpp"
@@ -76,47 +79,5 @@
 /** =========== define the function ====================== */
 
-inline
-HomologyGraph getFirstGraphwithSpecifiedElements(
-    const HomologyContainer &homologies,
-    const SerializablePotential::ParticleTypes_t &types)
-{
-  ASSERT( !types.empty(),
-      "getFirstGraphwithSpecifiedElements() - charges is empty?");
-  // create charges
-  Fragment::charges_t charges;
-  charges.resize(types.size());
-  std::transform(types.begin(), types.end(),
-      charges.begin(), boost::lambda::_1);
-  // convert into count map
-  Extractors::elementcounts_t counts_per_charge =
-      Extractors::_detail::getElementCounts(charges);
-  ASSERT( !counts_per_charge.empty(),
-      "getFirstGraphwithSpecifiedElements() - charge counts are empty?");
-  LOG(2, "DEBUG: counts_per_charge is " << counts_per_charge << ".");
-  // we want to check each (unique) key only once
-  HomologyContainer::const_key_iterator olditer = homologies.key_end();
-  for (HomologyContainer::const_key_iterator iter =
-      homologies.key_begin(); iter != homologies.key_end();
-      iter = homologies.getNextKey(iter)) {
-    // if it's the same as the old one, skip it
-    if (olditer == iter)
-      continue;
-    else
-      olditer = iter;
-    // if it's a new key, check if every element has the right number of counts
-    Extractors::elementcounts_t::const_iterator countiter = counts_per_charge.begin();
-    for (; countiter != counts_per_charge.end(); ++countiter)
-      if (!(*iter).hasTimesAtomicNumber(
-          static_cast<size_t>(countiter->first),
-          static_cast<size_t>(countiter->second))
-          )
-        break;
-    if( countiter == counts_per_charge.end())
-      return *iter;
-  }
-  return HomologyGraph();
-}
-
-void enforceZeroTotalCharge(
+static void enforceZeroTotalCharge(
     PartialNucleiChargeFitter::charges_t &_averaged_charges)
 {
@@ -134,7 +95,9 @@
       "enforceZeroTotalCharge() - enforcing neutral net charge failed, "
       +toString(charge_sum)+" is the remaining net charge.");
-}
-
-size_t obtainAverageChargesOverFragments(
+
+  LOG(2, "DEBUG: final charges with net zero charge are " << _averaged_charges);
+}
+
+static size_t obtainAverageChargesOverFragments(
     PartialNucleiChargeFitter::charges_t &_average_charges,
     const std::pair<
@@ -196,4 +159,16 @@
 }
 
+inline SerializablePotential::ParticleTypes_t
+getParticleTypesForAtomIdSet(const AtomIdSet &_atoms)
+{
+  SerializablePotential::ParticleTypes_t particletypes;
+  particletypes.resize(_atoms.size());
+  std::transform(
+      _atoms.begin(), _atoms.end(),
+      particletypes.begin(),
+      boost::bind(&atom::getElementNo, _1));
+  return particletypes;
+}
+
 ActionState::ptr PotentialFitPartialChargesAction::performCall()
 {
@@ -206,32 +181,55 @@
 
   /// obtain possible fragments to each selected atom
-  AtomFragmentsMap::keysets_t fragments;
-  const AtomFragmentsMap::AtomFragmentsMap_t &atommap =
-      AtomFragmentsMap::getConstInstance().getMap();
+  const AtomFragmentsMap &atomfragments = AtomFragmentsMap::getConstInstance();
+  if (!atomfragments.checkCompleteness()) {
+    STATUS("AtomFragmentsMap failed internal consistency check, missing forcekeysets?");
+    return Action::failure;
+  }
+  std::set<KeySet> fragments;
+  const AtomFragmentsMap::AtomFragmentsMap_t &atommap = atomfragments.getMap();
   for (World::AtomSelectionConstIterator iter = world.beginAtomSelection();
       iter != world.endAtomSelection(); ++iter) {
     const AtomFragmentsMap::AtomFragmentsMap_t::const_iterator atomiter =
         atommap.find(iter->first);
-    if (atomiter == atommap.end()) {
-      ELOG(2, "There are no fragments associated to " << iter->first << ".");
-      continue;
-    }
-    const AtomFragmentsMap::keysets_t &keysets = atomiter->second;
-    fragments.insert( fragments.end(), keysets.begin(), keysets.end() );
+    if (iter->second->getElementNo() != 1) {
+      if (atomiter == atommap.end()) {
+        ELOG(2, "There are no fragments associated to " << iter->first << ".");
+        continue;
+      }
+      const AtomFragmentsMap::keysets_t &keysets = atomiter->second;
+      LOG(2, "DEBUG: atom " << iter->first << " 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) {
+    const KeySet &keyset = *fragmentiter;
+    const AtomFragmentsMap::indices_t &forceindices = atomfragments.getFullKeyset(keyset);
+    ASSERT( !forceindices.empty(),
+        "PotentialFitPartialChargesAction::performCall() - force keyset to "
+        +toString(keyset)+" is empty.");
+    KeySet forcekeyset;
+    forcekeyset.insert(forceindices.begin(), forceindices.end());
+    forcekeyset.erase(-1);
+    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
-  typedef std::map<
-      KeySet,
-      PartialNucleiChargeFitter::charges_t> FragmentFittedChargeMap_t;
-  FragmentFittedChargeMap_t fittedcharges_per_fragment;
-  for (AtomFragmentsMap::keysets_t::const_iterator fragmentiter = fragments.begin();
-      fragmentiter != fragments.end(); ++fragmentiter) {
-    // fragment specifies the homology fragment to use
-    SerializablePotential::ParticleTypes_t fragmentnumbers(fragmentiter->begin(), fragmentiter->end());
-
-    // obtain range of homogolous fragments from container
-    HomologyContainer &homologies = World::getInstance().getHomologies();
-    const HomologyGraph graph = getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
+  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);
     if (range.first == range.second) {
@@ -241,9 +239,9 @@
 
     // fit and average partial charges over all homologous fragments
-    PartialNucleiChargeFitter::charges_t averaged_charges;
+    PartialNucleiChargeFitter::charges_t &averaged_charges = graphiter->second;
     const size_t NoFragments =
         obtainAverageChargesOverFragments(averaged_charges, range, params.radius.get());
     if ((NoFragments == 0) && (range.first != range.second)) {
-      STATUS("Fitting for fragment "+toString(*fragmentiter)+" failed.");
+      STATUS("Fitting for fragment "+toString(*graphiter)+" failed.");
       return Action::failure;
     }
@@ -253,9 +251,6 @@
       enforceZeroTotalCharge(averaged_charges);
 
-    // place into map for later retrieval
-    fittedcharges_per_fragment.insert( std::make_pair(*fragmentiter, averaged_charges));
-
     // output status info fitted charges
-    LOG(2, "DEBUG: For fragment " << *fragmentiter << " we have fitted the following charges "
+    LOG(2, "DEBUG: For fragment " << *graphiter << " we have fitted the following charges "
         << averaged_charges << ", averaged over " << NoFragments << " fragments.");
   }
@@ -266,9 +261,21 @@
   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(atomiter->first);
+        atommap.find(walkerid);
     ASSERT(keysetsiter != atommap.end(),
         "PotentialFitPartialChargesAction::performCall() - we checked already that "
-        +toString(atomiter->first)+" should be present!");
+        +toString(walkerid)+" should be present!");
     const AtomFragmentsMap::keysets_t & keysets = keysetsiter->second;
 
@@ -280,22 +287,29 @@
       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 FragmentFittedChargeMap_t::const_iterator chargesiter =
-          fittedcharges_per_fragment.find(keyset);
+      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 KeySet::const_iterator keysetiter = keyset.find(atomiter->first);
-      ASSERT( keysetiter != keyset.end(),
+      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 keyset "+toString(keyset));
-      std::advance(chargeiter, std::distance(keyset.begin(), keysetiter));
-      ASSERT( chargeiter != charges.end(),
-          "PotentialFitPartialChargesAction::performCall() - charges "
-          +toString(charges.size())+" and keyset "+toString(keyset.size())
-          +" do not have the same length?");
+          +toString(atomiter->first)+" not contained in force keyset "+toString(forcekeyset));
+      std::advance(chargeiter, std::distance(forcekeyset.begin(), forcekeysetiter));
 
       // and add onto charge sum
@@ -305,6 +319,9 @@
     }
     // average to obtain final partial charge for this atom
-    average_charge = 1./(size_t)NoFragments;
+    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);
   }
 
