Context Navigation

-              rd9dbef
+              rd8ed62
 #include "Element/element.hpp"
 #include "Element/periodentafel.hpp"
+#include "Fragmentation/Homology/AtomFragmentsMap.hpp"
 #include "Fragmentation/Homology/HomologyContainer.hpp"
 #include "Fragmentation/Homology/HomologyGraph.hpp"
 …
+}
+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;
 …
             && (potential.begin[2] == potential.end[2]))) {
       ELOG(1, "Sampled grid contains grid made of zero points.");
       return Action::failure;
+      return 0;
+    }
 …
       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;

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset d8ed62 for src/Actions/PotentialAction/FitPartialChargesAction.cpp

Legend:

src/Actions/PotentialAction/FitPartialChargesAction.cpp

Download in other formats: