Context Navigation

-              rd8ed62
+              r01120c
 #include "World.hpp"
+#include <boost/bimap.hpp>
 #include <boost/bind.hpp>
 #include <boost/filesystem.hpp>
 …
 #include "Potentials/PartialNucleiChargeFitter.hpp"
+#include "AtomIdSet.hpp"
+#include "Descriptors/AtomIdDescriptor.hpp"
 #include "Element/element.hpp"
 #include "Element/periodentafel.hpp"
 …
 /** =========== 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)
+{
 …
       "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<
 …
+}
+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()
+{
 …
   /// 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) {
 …
     // 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;
+    }
 …
       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.");
+  }
 …
   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;
 …
       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
 …
+    }
     // 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);
+  }

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

Context Navigation

Changeset 01120c for src/Actions

Legend:

src/Actions/PotentialAction/FitPartialChargesAction.cpp

Download in other formats: