Context Navigation

-              rc1ec8e
+              r91f907
 #include <boost/bimap.hpp>
+#include <boost/bimap/multiset_of.hpp>
 #include <boost/bind.hpp>
 #include <boost/filesystem.hpp>
 …
 #include "Action_impl_pre.hpp"
 /** =========== define the function ====================== */
+namespace detail {
+  typedef std::map<KeySet, HomologyGraph> KeysetsToGraph_t;
+  typedef std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> GraphFittedChargeMap_t;
+  typedef std::map<atomId_t, double> fitted_charges_t;
+  typedef std::map<HomologyGraph, size_t> GraphIndex_t;
+  typedef std::set<size_t> AtomsGraphIndices_t;
+  typedef boost::bimaps::bimap<
+      AtomsGraphIndices_t,
+      boost::bimaps::multiset_of<atomId_t> > GraphIndices_t;
+  typedef std::map<atomId_t, std::string> AtomParticleNames_t;
+  typedef std::map<std::set<size_t>, std::string> GraphToNameMap_t;
+};
 static void enforceZeroTotalCharge(
 …
 static
 void getKeySetsToGraphMapping(
     std::map<KeySet, HomologyGraph> &_keyset_graphs,
     std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> &_fittedcharges_per_fragment,
+    detail::KeysetsToGraph_t &_keyset_graphs,
+    detail::GraphFittedChargeMap_t &_fittedcharges_per_fragment,
     const std::set<KeySet> &_fragments,
     const AtomFragmentsMap &_atomfragments)
 …
 static
 bool getPartialChargesForAllGraphs(
     std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> &_fittedcharges_per_fragment,
+    detail::GraphFittedChargeMap_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();
+  for (detail::GraphFittedChargeMap_t::iterator graphiter = _fittedcharges_per_fragment.begin();
       graphiter != _fittedcharges_per_fragment.end(); ++graphiter) {
     const HomologyGraph &graph = graphiter->first;
 …
     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 detail::KeysetsToGraph_t &_keyset_graphs,
+    const detail::GraphFittedChargeMap_t &_fittedcharges_per_fragment)
+{
   const atom * surrogate = _walker;
   if (surrogate->getElementNo() == 1) {
 …
     // find the associated charge in the charge vector
     std::map<KeySet, HomologyGraph>::const_iterator keysetgraphiter =
+    const std::map<KeySet, HomologyGraph>::const_iterator keysetgraphiter =
         _keyset_graphs.find(keyset);
     ASSERT( keysetgraphiter != _keyset_graphs.end(),
 …
         +" not contained in keyset_graphs.");
     const HomologyGraph &graph = keysetgraphiter->second;
     const GraphFittedChargeMap_t::const_iterator chargesiter =
+    const detail::GraphFittedChargeMap_t::const_iterator chargesiter =
         _fittedcharges_per_fragment.find(graph);
     ASSERT(chargesiter != _fittedcharges_per_fragment.end(),
 …
     const ParticleFactory &factory,
     const periodentafel &periode,
+    const std::map<atomId_t, double> &_fitted_charges)
+{
+  typedef std::map<atomId_t, double> fitted_charges_t;
+  for (fitted_charges_t::const_iterator chargeiter = _fitted_charges.begin();
+    const detail::fitted_charges_t &_fitted_charges,
+    const detail::GraphIndices_t &_GraphIndices,
+    detail::AtomParticleNames_t &_atom_particlenames)
+{
+  detail::GraphToNameMap_t GraphToNameMap;
+  for (detail::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);
+    const atomId_t &atomid = chargeiter->first;
+    const detail::GraphIndices_t::right_const_iterator graphiter = _GraphIndices.right.find(atomid);
+    const detail::GraphToNameMap_t::const_iterator nameiter = GraphToNameMap.find(graphiter->second);
+    std::string name;
+    if (nameiter != GraphToNameMap.end()) {
+      name = nameiter->second;
+    } else {
+      const atom * const walker = World::getInstance().getAtom(AtomById(atomid));
+      ASSERT( walker != NULL,
+          "addToParticleRegistry() - atom "+toString(atomid)+" not present in the World?");
+      const double &charge = chargeiter->second;
+      const atomicNumber_t elementno = walker->getElementNo();
+      name = Particle::findFreeName(periode, elementno);
+      GraphToNameMap[graphiter->second] = name;
+      LOG(1, "INFO: Adding particle " << name << " for atom "
+          << *walker << " with element " << elementno << ", charge " << charge);
+      factory.createInstance(name, elementno, charge);
+    }
+    _atom_particlenames.insert( std::make_pair(atomid, name) );
+  }
+}
 …
   // 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;
+  detail::KeysetsToGraph_t keyset_graphs;
+  detail::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 HomologyContainer &homologies = World::getInstance().getHomologies();
   const bool status = getPartialChargesForAllGraphs(
       fittedcharges_per_fragment,
       World::getInstance().getHomologies(),
+      homologies,
       params.radius.get(),
       params.enforceZeroCharge.get());
 …
   /// 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;
+  detail::fitted_charges_t fitted_charges;
   for (World::AtomSelectionConstIterator atomiter = world.beginAtomSelection();
       atomiter != world.endAtomSelection(); ++atomiter) {
 …
+  }
+  /// make Particles be used for every atom that was fitted on the same number of graphs
+  detail::GraphIndex_t GraphIndex;
+  size_t index = 0;
+  for (HomologyContainer::const_key_iterator iter = homologies.key_begin();
+      iter != homologies.key_end(); iter = homologies.getNextKey(iter)) {
+    GraphIndex.insert( std::make_pair( *iter, index++));
+  }
+  LOG(2, "DEBUG: There are " << index << " unique graphs in the homology container.");
+  // go through every atom, get all graphs, convert to GraphIndex and store
+  detail::GraphIndices_t GraphIndices;
+  const AtomFragmentsMap::AtomFragmentsMap_t &atommap = atomfragments.getMap();
+  for (World::AtomSelectionConstIterator atomiter = world.beginAtomSelection();
+      atomiter != world.endAtomSelection(); ++atomiter) {
+    const atomId_t walkerid = atomiter->first;
+    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;
+    detail::AtomsGraphIndices_t AtomsGraphIndices;
+    // 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 std::map<KeySet, HomologyGraph>::const_iterator keysetgraphiter =
+          keyset_graphs.find(keyset);
+      ASSERT( keysetgraphiter != keyset_graphs.end(),
+          "PotentialFitPartialChargesAction::performCall() - keyset "+toString(keyset)
+          +" not contained in keyset_graphs.");
+      const HomologyGraph &graph = keysetgraphiter->second;
+      const detail::GraphIndex_t::const_iterator indexiter = GraphIndex.find(graph);
+      ASSERT( indexiter != GraphIndex.end(),
+          "PotentialFitPartialChargesAction::performCall() - graph "+toString(graph)
+          +" not contained in GraphIndex.");
+      AtomsGraphIndices.insert( indexiter->second );
+    }
+    GraphIndices.left.insert( std::make_pair(AtomsGraphIndices, walkerid) );
+    LOG(2, "DEBUG: Atom " << *atomiter->second << " has graph indices " << AtomsGraphIndices);
+  }
   /// place all fitted charges into ParticleRegistry
+  detail::AtomParticleNames_t atom_particlenames;
   addToParticleRegistry(
       ParticleFactory::getConstInstance(),
       *World::getInstance().getPeriode(),
+      fitted_charges);
+      fitted_charges,
+      GraphIndices,
+      atom_particlenames);
+  LOG(1, "INFO: The atoms received the following particles " << atom_particlenames);
   return Action::success;

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

Context Navigation

Changeset 91f907 for src/Actions

Legend:

src/Actions/PotentialAction/FitPartialChargesAction.cpp

Download in other formats: