Index: src/Actions/PotentialAction/FitPartialChargesAction.cpp
===================================================================
--- src/Actions/PotentialAction/FitPartialChargesAction.cpp	(revision 78202bb746daf78c7efca90c8f594001d926da2c)
+++ src/Actions/PotentialAction/FitPartialChargesAction.cpp	(revision d8ed6292b835bf1c1008f7aba846e47d5a612ea4)
@@ -58,4 +58,5 @@
 #include "Element/element.hpp"
 #include "Element/periodentafel.hpp"
+#include "Fragmentation/Homology/AtomFragmentsMap.hpp"
 #include "Fragmentation/Homology/HomologyContainer.hpp"
 #include "Fragmentation/Homology/HomologyGraph.hpp"
@@ -117,38 +118,42 @@
 }
 
-ActionState::ptr PotentialFitPartialChargesAction::performCall() {
-
-  // fragment specifies the homology fragment to use
-  SerializablePotential::ParticleTypes_t fragmentnumbers;
-  {
-    const std::vector<const element *> &fragment = params.fragment.get();
-    std::transform(fragment.begin(), fragment.end(), std::back_inserter(fragmentnumbers),
-        boost::bind(&element::getAtomicNumber, _1));
-  }
-
-  // parse homologies into container
-  HomologyContainer &homologies = World::getInstance().getHomologies();
-  const HomologyGraph graph = getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
-  HomologyContainer::range_t range = homologies.getHomologousGraphs(graph);
-  // for the moment just use the very first fragment
-  if (range.first == range.second) {
-    STATUS("HomologyContainer does not contain specified fragment.");
-    return Action::failure;
-  }
-
-  // average partial charges over all fragments
-  HomologyContainer::const_iterator iter = range.first;
+void enforceZeroTotalCharge(
+    PartialNucleiChargeFitter::charges_t &_averaged_charges)
+{
+  double charge_sum = 0.;
+  charge_sum = std::accumulate(_averaged_charges.begin(), _averaged_charges.end(), charge_sum);
+  if (fabs(charge_sum) > MYEPSILON) {
+    std::transform(
+        _averaged_charges.begin(), _averaged_charges.end(),
+        _averaged_charges.begin(),
+        boost::bind(std::minus<double>(), _1, charge_sum/_averaged_charges.size()));
+  }
+  charge_sum = 0.;
+  charge_sum = std::accumulate(_averaged_charges.begin(), _averaged_charges.end(), charge_sum);
+  ASSERT( fabs(charge_sum) < MYEPSILON,
+      "enforceZeroTotalCharge() - enforcing neutral net charge failed, "
+      +toString(charge_sum)+" is the remaining net charge.");
+}
+
+size_t obtainAverageChargesOverFragments(
+    PartialNucleiChargeFitter::charges_t &_average_charges,
+    const std::pair<
+      HomologyContainer::const_iterator,
+      HomologyContainer::const_iterator> &_range,
+    const double _radius
+    )
+{
+  HomologyContainer::const_iterator iter = _range.first;
   if (!iter->second.containsGrids) {
-    STATUS("This HomologyGraph does not contain sampled grids.");
-    return Action::failure;
-  }
-  PartialNucleiChargeFitter::charges_t averaged_charges;
-  averaged_charges.resize(iter->second.fragment.getCharges().size(), 0.);
+    ELOG(1, "This HomologyGraph does not contain sampled grids.");
+    return 0;
+  }
+  _average_charges.resize(iter->second.fragment.getCharges().size(), 0.);
   size_t NoFragments = 0;
   for (;
-      iter != range.second; ++iter, ++NoFragments) {
+      iter != _range.second; ++iter, ++NoFragments) {
     if (!iter->second.containsGrids) {
       ELOG(2, "This HomologyGraph does not contain sampled grids,\ndid you forget to add '--store-grids 1' to AnalyseFragmentResults.");
-      return Action::failure;
+      return 0;
     }
     const Fragment &fragment = iter->second.fragment;
@@ -161,5 +166,5 @@
             && (potential.begin[2] == potential.end[2]))) {
       ELOG(1, "Sampled grid contains grid made of zero points.");
-      return Action::failure;
+      return 0;
     }
 
@@ -171,51 +176,154 @@
       positions.push_back( Vector(pos[0], pos[1], pos[2]) );
     }
-    PartialNucleiChargeFitter fitter(potential, positions, params.radius.get());
+    PartialNucleiChargeFitter fitter(potential, positions, _radius);
     fitter();
     PartialNucleiChargeFitter::charges_t return_charges = fitter.getSolutionAsCharges_t();
+    LOG(2, "DEBUG: fitted charges are " << return_charges);
     std::transform(
         return_charges.begin(), return_charges.end(),
-        averaged_charges.begin(),
-        averaged_charges.begin(),
+        _average_charges.begin(),
+        _average_charges.begin(),
         std::plus<PartialNucleiChargeFitter::charge_t>());
   }
-  std::transform(averaged_charges.begin(),averaged_charges.end(),
-      averaged_charges.begin(),
-      std::bind1st(std::multiplies<PartialNucleiChargeFitter::charge_t>(),1./NoFragments)
-  );
-
-  // make sum of charges zero if desired
-  if (params.enforceZeroCharge.get()) {
-    double charge_sum = 0.;
-    charge_sum = std::accumulate(averaged_charges.begin(), averaged_charges.end(), charge_sum);
-    if (fabs(charge_sum) > MYEPSILON) {
-      std::transform(
-          averaged_charges.begin(), averaged_charges.end(),
-          averaged_charges.begin(),
-          boost::bind(std::minus<double>(), _1, charge_sum/averaged_charges.size()));
-    }
-    charge_sum = 0.;
-    charge_sum = std::accumulate(averaged_charges.begin(), averaged_charges.end(), charge_sum);
-    ASSERT( fabs(charge_sum) < MYEPSILON,
-        "PotentialFitPartialChargesAction::performCall() - enforcing neutral net charge failed, "
-        +toString(charge_sum)+" is the remaining net charge.");
-  }
-
-  // place all fitted charges into ParticleRegistry
+  if (NoFragments != 0)
+    std::transform(_average_charges.begin(), _average_charges.end(),
+        _average_charges.begin(),
+        std::bind1st(std::multiplies<PartialNucleiChargeFitter::charge_t>(),1./(double)NoFragments)
+    );
+  LOG(2, "DEBUG: final averaged charges are " << _average_charges);
+
+  return NoFragments;
+}
+
+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
+  AtomFragmentsMap::keysets_t fragments;
+  const AtomFragmentsMap::AtomFragmentsMap_t &atommap =
+      AtomFragmentsMap::getConstInstance().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() );
+  }
+
+  /// 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);
+    const HomologyContainer::range_t range = homologies.getHomologousGraphs(graph);
+    if (range.first == range.second) {
+      STATUS("HomologyContainer does not contain specified fragment.");
+      return Action::failure;
+    }
+
+    // fit and average partial charges over all homologous fragments
+    PartialNucleiChargeFitter::charges_t averaged_charges;
+    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.");
+      return Action::failure;
+    }
+
+    // make sum of charges zero if desired
+    if (params.enforceZeroCharge.get())
+      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 "
+        << averaged_charges << ", averaged over " << NoFragments << " fragments.");
+  }
+
+  /// 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 AtomFragmentsMap::AtomFragmentsMap_t::const_iterator keysetsiter =
+        atommap.find(atomiter->first);
+    ASSERT(keysetsiter != atommap.end(),
+        "PotentialFitPartialChargesAction::performCall() - we checked already that "
+        +toString(atomiter->first)+" 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;
+
+      // find the associated charge in the charge vector
+      const FragmentFittedChargeMap_t::const_iterator chargesiter =
+          fittedcharges_per_fragment.find(keyset);
+      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;
+      PartialNucleiChargeFitter::charges_t::const_iterator chargeiter =
+          charges.begin();
+      const KeySet::const_iterator keysetiter = keyset.find(atomiter->first);
+      ASSERT( keysetiter != keyset.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?");
+
+      // 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./(size_t)NoFragments;
+    fitted_charges.insert( std::make_pair(atomiter->first, average_charge) );
+  }
+
+  /// place all fitted charges into ParticleRegistry
   const ParticleFactory &factory =
       const_cast<const ParticleFactory&>(ParticleFactory::getInstance());
   const periodentafel &periode = *World::getInstance().getPeriode();
-  ASSERT(averaged_charges.size() == fragmentnumbers.size(),
-      "PotentialFitPartialChargesAction::performCall() - charges and elements length mismatch.");
-  for (size_t i=0;i<averaged_charges.size(); ++i) {
-    const std::string name = Particle::findFreeName(periode, fragmentnumbers[i]);
-    LOG(2, "INFO: Adding particle " << name << " for element "
-        << fragmentnumbers[i] << ", charge " << averaged_charges[i]);
-    factory.createInstance(name, fragmentnumbers[i], averaged_charges[i]);
-  }
-
-  // output fitted charges
-  LOG(0, "STATUS: We have fitted the following charges " << averaged_charges
-      << ", averaged over " << NoFragments << " fragments.");
+  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,
+        "PotentialFitPartialChargesAction::performCall() - atom "
+        +toString(chargeiter->first)+" not present in the World?");
+    const double &charge = chargeiter->second;
+    const atomicNumber_t elementno = walker->getElementNo();
+    const std::string name = Particle::findFreeName(periode, elementno);
+    LOG(1, "INFO: Adding particle " << name << " for atom "
+        << *walker << " with element " << elementno << ", charge " << charge);
+    factory.createInstance(name, elementno, charge);
+  }
 
   return Action::success;
Index: src/Actions/PotentialAction/FitPartialChargesAction.def
===================================================================
--- src/Actions/PotentialAction/FitPartialChargesAction.def	(revision 78202bb746daf78c7efca90c8f594001d926da2c)
+++ src/Actions/PotentialAction/FitPartialChargesAction.def	(revision d8ed6292b835bf1c1008f7aba846e47d5a612ea4)
@@ -19,11 +19,10 @@
 // ValueStorage by the token "Z" -> first column: int, Z, "Z"
 // "undefine" if no parameters are required, use (NOPARAM_DEFAULT) for each (undefined) default value
-#define paramtypes (std::vector<const element *>)(double)(bool)
-#define paramtokens ("fragment-charges")("radius")("enforce-net-zero-charge")
-#define paramdescriptions ("charges specifying the fragment")("radius of sphere around nuclei where potential does not need to match")("whether to make the sum of fitted charged become zero")
-#define paramdefaults (NOPARAM_DEFAULT)(PARAM_DEFAULT(0))(PARAM_DEFAULT(1))
-#define paramreferences (fragment)(radius)(enforceZeroCharge)
+#define paramtypes (double)(bool)
+#define paramtokens ("radius")("enforce-net-zero-charge")
+#define paramdescriptions ("radius of sphere around nuclei where potential does not need to match")("whether to make the sum of fitted charged become zero")
+#define paramdefaults (PARAM_DEFAULT(0))(PARAM_DEFAULT(1))
+#define paramreferences (radius)(enforceZeroCharge)
 #define paramvalids \
-(STLVectorValidator< std::vector<const element *> >(1,99, ElementValidator())) \
 (NonNegativeValidator<double>()) \
 (DummyValidator<bool>())
@@ -41,4 +40,4 @@
 
 // finally the information stored in the ActionTrait specialization
-#define DESCRIPTION "fits partial nuclear charges to present particle types from loaded homologies containing sampled grids."
+#define DESCRIPTION "fits partial nuclear charges to selected atoms from averages over homologous fragments containing sampled charge grids."
 #undef SHORTFORM
