| 1 | /*
 | 
|---|
| 2 |  * Project: MoleCuilder
 | 
|---|
| 3 |  * Description: creates and alters molecular systems
 | 
|---|
| 4 |  * Copyright (C)  2013 Frederik Heber. All rights reserved.
 | 
|---|
| 5 |  * 
 | 
|---|
| 6 |  *
 | 
|---|
| 7 |  *   This file is part of MoleCuilder.
 | 
|---|
| 8 |  *
 | 
|---|
| 9 |  *    MoleCuilder is free software: you can redistribute it and/or modify
 | 
|---|
| 10 |  *    it under the terms of the GNU General Public License as published by
 | 
|---|
| 11 |  *    the Free Software Foundation, either version 2 of the License, or
 | 
|---|
| 12 |  *    (at your option) any later version.
 | 
|---|
| 13 |  *
 | 
|---|
| 14 |  *    MoleCuilder is distributed in the hope that it will be useful,
 | 
|---|
| 15 |  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 | 
|---|
| 16 |  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 | 
|---|
| 17 |  *    GNU General Public License for more details.
 | 
|---|
| 18 |  *
 | 
|---|
| 19 |  *    You should have received a copy of the GNU General Public License
 | 
|---|
| 20 |  *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 | 
|---|
| 21 |  */
 | 
|---|
| 22 | 
 | 
|---|
| 23 | /*
 | 
|---|
| 24 |  * FitPartialChargesAction.cpp
 | 
|---|
| 25 |  *
 | 
|---|
| 26 |  *  Created on: Jul 03, 2013
 | 
|---|
| 27 |  *      Author: heber
 | 
|---|
| 28 |  */
 | 
|---|
| 29 | 
 | 
|---|
| 30 | // include config.h
 | 
|---|
| 31 | #ifdef HAVE_CONFIG_H
 | 
|---|
| 32 | #include <config.h>
 | 
|---|
| 33 | #endif
 | 
|---|
| 34 | 
 | 
|---|
| 35 | // needs to come before MemDebug due to placement new
 | 
|---|
| 36 | #include <boost/archive/text_iarchive.hpp>
 | 
|---|
| 37 | 
 | 
|---|
| 38 | #include "CodePatterns/MemDebug.hpp"
 | 
|---|
| 39 | 
 | 
|---|
| 40 | #include "Atom/atom.hpp"
 | 
|---|
| 41 | #include "CodePatterns/Log.hpp"
 | 
|---|
| 42 | #include "Fragmentation/Exporters/ExportGraph_ToFiles.hpp"
 | 
|---|
| 43 | #include "Fragmentation/Graph.hpp"
 | 
|---|
| 44 | #include "World.hpp"
 | 
|---|
| 45 | 
 | 
|---|
| 46 | #include <boost/bimap.hpp>
 | 
|---|
| 47 | #include <boost/bimap/multiset_of.hpp>
 | 
|---|
| 48 | #include <boost/bind.hpp>
 | 
|---|
| 49 | #include <boost/filesystem.hpp>
 | 
|---|
| 50 | #include <boost/foreach.hpp>
 | 
|---|
| 51 | #include <algorithm>
 | 
|---|
| 52 | #include <functional>
 | 
|---|
| 53 | #include <iostream>
 | 
|---|
| 54 | #include <string>
 | 
|---|
| 55 | 
 | 
|---|
| 56 | #include "Actions/PotentialAction/FitPartialChargesAction.hpp"
 | 
|---|
| 57 | 
 | 
|---|
| 58 | #include "Potentials/PartialNucleiChargeFitter.hpp"
 | 
|---|
| 59 | 
 | 
|---|
| 60 | #include "AtomIdSet.hpp"
 | 
|---|
| 61 | #include "Descriptors/AtomIdDescriptor.hpp"
 | 
|---|
| 62 | #include "Element/element.hpp"
 | 
|---|
| 63 | #include "Element/periodentafel.hpp"
 | 
|---|
| 64 | #include "Fragmentation/Homology/AtomFragmentsMap.hpp"
 | 
|---|
| 65 | #include "Fragmentation/Homology/HomologyContainer.hpp"
 | 
|---|
| 66 | #include "Fragmentation/Homology/HomologyGraph.hpp"
 | 
|---|
| 67 | #include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
 | 
|---|
| 68 | #include "FunctionApproximation/Extractors.hpp"
 | 
|---|
| 69 | #include "Potentials/PartialNucleiChargeFitter.hpp"
 | 
|---|
| 70 | #include "Potentials/Particles/ParticleFactory.hpp"
 | 
|---|
| 71 | #include "Potentials/Particles/ParticleRegistry.hpp"
 | 
|---|
| 72 | #include "Potentials/SerializablePotential.hpp"
 | 
|---|
| 73 | #include "World.hpp"
 | 
|---|
| 74 | 
 | 
|---|
| 75 | using namespace MoleCuilder;
 | 
|---|
| 76 | 
 | 
|---|
| 77 | // and construct the stuff
 | 
|---|
| 78 | #include "FitPartialChargesAction.def"
 | 
|---|
| 79 | #include "Action_impl_pre.hpp"
 | 
|---|
| 80 | /** =========== define the function ====================== */
 | 
|---|
| 81 | 
 | 
|---|
| 82 | namespace detail {
 | 
|---|
| 83 |   typedef std::map<KeySet, HomologyGraph> KeysetsToGraph_t;
 | 
|---|
| 84 | 
 | 
|---|
| 85 |   typedef std::map<HomologyGraph, PartialNucleiChargeFitter::charges_t> GraphFittedChargeMap_t;
 | 
|---|
| 86 | 
 | 
|---|
| 87 |   typedef std::map<atomId_t, double> fitted_charges_t;
 | 
|---|
| 88 | 
 | 
|---|
| 89 |   typedef std::map<HomologyGraph, size_t> GraphIndex_t;
 | 
|---|
| 90 | 
 | 
|---|
| 91 |   typedef std::set<size_t> AtomsGraphIndices_t;
 | 
|---|
| 92 |   typedef boost::bimaps::bimap<
 | 
|---|
| 93 |       boost::bimaps::multiset_of<AtomsGraphIndices_t>,
 | 
|---|
| 94 |       atomId_t > GraphIndices_t;
 | 
|---|
| 95 | 
 | 
|---|
| 96 |   typedef std::map<std::set<size_t>, std::map<atomicNumber_t, std::string> > AtomParticleNames_t;
 | 
|---|
| 97 | 
 | 
|---|
| 98 |   typedef std::map<std::set<size_t>, std::string> GraphToNameMap_t;
 | 
|---|
| 99 | };
 | 
|---|
| 100 | 
 | 
|---|
| 101 | static void enforceZeroTotalCharge(
 | 
|---|
| 102 |     PartialNucleiChargeFitter::charges_t &_averaged_charges)
 | 
|---|
| 103 | {
 | 
|---|
| 104 |   double charge_sum = 0.;
 | 
|---|
| 105 |   charge_sum = std::accumulate(_averaged_charges.begin(), _averaged_charges.end(), charge_sum);
 | 
|---|
| 106 |   if (fabs(charge_sum) > MYEPSILON) {
 | 
|---|
| 107 |     std::transform(
 | 
|---|
| 108 |         _averaged_charges.begin(), _averaged_charges.end(),
 | 
|---|
| 109 |         _averaged_charges.begin(),
 | 
|---|
| 110 |         boost::bind(std::minus<double>(), _1, charge_sum/_averaged_charges.size()));
 | 
|---|
| 111 |   }
 | 
|---|
| 112 |   charge_sum = 0.;
 | 
|---|
| 113 |   charge_sum = std::accumulate(_averaged_charges.begin(), _averaged_charges.end(), charge_sum);
 | 
|---|
| 114 |   ASSERT( fabs(charge_sum) < MYEPSILON,
 | 
|---|
| 115 |       "enforceZeroTotalCharge() - enforcing neutral net charge failed, "
 | 
|---|
| 116 |       +toString(charge_sum)+" is the remaining net charge.");
 | 
|---|
| 117 | 
 | 
|---|
| 118 |   LOG(2, "DEBUG: final charges with net zero charge are " << _averaged_charges);
 | 
|---|
| 119 | }
 | 
|---|
| 120 | 
 | 
|---|
| 121 | static size_t obtainAverageChargesOverFragments(
 | 
|---|
| 122 |     PartialNucleiChargeFitter::charges_t &_average_charges,
 | 
|---|
| 123 |     const std::pair<
 | 
|---|
| 124 |       HomologyContainer::const_iterator,
 | 
|---|
| 125 |       HomologyContainer::const_iterator> &_range,
 | 
|---|
| 126 |     const double _radius
 | 
|---|
| 127 |     )
 | 
|---|
| 128 | {
 | 
|---|
| 129 |   HomologyContainer::const_iterator iter = _range.first;
 | 
|---|
| 130 |   if (!iter->second.containsGrids) {
 | 
|---|
| 131 |     ELOG(1, "This HomologyGraph does not contain sampled grids.");
 | 
|---|
| 132 |     return 0;
 | 
|---|
| 133 |   }
 | 
|---|
| 134 |   _average_charges.resize(iter->second.fragment.getCharges().size(), 0.);
 | 
|---|
| 135 |   size_t NoFragments = 0;
 | 
|---|
| 136 |   for (;
 | 
|---|
| 137 |       iter != _range.second; ++iter, ++NoFragments) {
 | 
|---|
| 138 |     if (!iter->second.containsGrids) {
 | 
|---|
| 139 |       ELOG(2, "This HomologyGraph does not contain sampled grids,\ndid you forget to add '--store-grids 1' to AnalyseFragmentResults.");
 | 
|---|
| 140 |       return 0;
 | 
|---|
| 141 |     }
 | 
|---|
| 142 |     const Fragment &fragment = iter->second.fragment;
 | 
|---|
| 143 |   //  const double &energy = iter->second.energy;
 | 
|---|
| 144 |   //  const SamplingGrid &charge = iter->second.charge_distribution;
 | 
|---|
| 145 |     const SamplingGrid &potential = iter->second.potential_distribution;
 | 
|---|
| 146 |     if ((potential.level == 0)
 | 
|---|
| 147 |         || ((potential.begin[0] == potential.end[0])
 | 
|---|
| 148 |             && (potential.begin[1] == potential.end[1])
 | 
|---|
| 149 |             && (potential.begin[2] == potential.end[2]))) {
 | 
|---|
| 150 |       ELOG(1, "Sampled grid contains grid made of zero points.");
 | 
|---|
| 151 |       return 0;
 | 
|---|
| 152 |     }
 | 
|---|
| 153 | 
 | 
|---|
| 154 |     // then we extract positions from fragment
 | 
|---|
| 155 |     PartialNucleiChargeFitter::positions_t positions;
 | 
|---|
| 156 |     Fragment::positions_t fragmentpositions = fragment.getPositions();
 | 
|---|
| 157 |     positions.reserve(fragmentpositions.size());
 | 
|---|
| 158 |     BOOST_FOREACH( Fragment::position_t pos, fragmentpositions) {
 | 
|---|
| 159 |       positions.push_back( Vector(pos[0], pos[1], pos[2]) );
 | 
|---|
| 160 |     }
 | 
|---|
| 161 |     PartialNucleiChargeFitter fitter(potential, positions, _radius);
 | 
|---|
| 162 |     fitter();
 | 
|---|
| 163 |     PartialNucleiChargeFitter::charges_t return_charges = fitter.getSolutionAsCharges_t();
 | 
|---|
| 164 |     LOG(2, "DEBUG: fitted charges are " << return_charges);
 | 
|---|
| 165 |     std::transform(
 | 
|---|
| 166 |         return_charges.begin(), return_charges.end(),
 | 
|---|
| 167 |         _average_charges.begin(),
 | 
|---|
| 168 |         _average_charges.begin(),
 | 
|---|
| 169 |         std::plus<PartialNucleiChargeFitter::charge_t>());
 | 
|---|
| 170 |   }
 | 
|---|
| 171 |   if (NoFragments != 0)
 | 
|---|
| 172 |     std::transform(_average_charges.begin(), _average_charges.end(),
 | 
|---|
| 173 |         _average_charges.begin(),
 | 
|---|
| 174 |         std::bind1st(std::multiplies<PartialNucleiChargeFitter::charge_t>(),1./(double)NoFragments)
 | 
|---|
| 175 |     );
 | 
|---|
| 176 |   LOG(2, "DEBUG: final averaged charges are " << _average_charges);
 | 
|---|
| 177 | 
 | 
|---|
| 178 |   return NoFragments;
 | 
|---|
| 179 | }
 | 
|---|
| 180 | 
 | 
|---|
| 181 | inline SerializablePotential::ParticleTypes_t
 | 
|---|
| 182 | getParticleTypesForAtomIdSet(const AtomIdSet &_atoms)
 | 
|---|
| 183 | {
 | 
|---|
| 184 |   SerializablePotential::ParticleTypes_t particletypes;
 | 
|---|
| 185 |   particletypes.resize(_atoms.size());
 | 
|---|
| 186 |   std::transform(
 | 
|---|
| 187 |       _atoms.begin(), _atoms.end(),
 | 
|---|
| 188 |       particletypes.begin(),
 | 
|---|
| 189 |       boost::bind(&atom::getElementNo, _1));
 | 
|---|
| 190 |   return particletypes;
 | 
|---|
| 191 | }
 | 
|---|
| 192 | 
 | 
|---|
| 193 | static
 | 
|---|
| 194 | std::set<KeySet> accumulateKeySetsForAtoms(
 | 
|---|
| 195 |     const AtomFragmentsMap::AtomFragmentsMap_t &_atommap,
 | 
|---|
| 196 |     const std::vector<const atom *> &_selected_atoms)
 | 
|---|
| 197 | {
 | 
|---|
| 198 |   std::set<KeySet> fragments;
 | 
|---|
| 199 |   for (std::vector<const atom *>::const_iterator iter = _selected_atoms.begin();
 | 
|---|
| 200 |       iter != _selected_atoms.end(); ++iter) {
 | 
|---|
| 201 |     const atomId_t atomid = (*iter)->getId();
 | 
|---|
| 202 |     const AtomFragmentsMap::AtomFragmentsMap_t::const_iterator atomiter =
 | 
|---|
| 203 |         _atommap.find(atomid);
 | 
|---|
| 204 |     if ((*iter)->getElementNo() != 1) {
 | 
|---|
| 205 |       if (atomiter == _atommap.end()) {
 | 
|---|
| 206 |         ELOG(2, "There are no fragments associated to " << atomid << ".");
 | 
|---|
| 207 |         continue;
 | 
|---|
| 208 |       }
 | 
|---|
| 209 |       const AtomFragmentsMap::keysets_t &keysets = atomiter->second;
 | 
|---|
| 210 |       LOG(2, "DEBUG: atom " << atomid << " has " << keysets.size() << " fragments.");
 | 
|---|
| 211 |       fragments.insert( keysets.begin(), keysets.end() );
 | 
|---|
| 212 |     } else {
 | 
|---|
| 213 |       LOG(3, "DEBUG: Skipping atom " << atomid << " as it's hydrogen.");
 | 
|---|
| 214 |     }
 | 
|---|
| 215 |   }
 | 
|---|
| 216 |   return fragments;
 | 
|---|
| 217 | }
 | 
|---|
| 218 | 
 | 
|---|
| 219 | static
 | 
|---|
| 220 | void getKeySetsToGraphMapping(
 | 
|---|
| 221 |     detail::KeysetsToGraph_t &_keyset_graphs,
 | 
|---|
| 222 |     detail::GraphFittedChargeMap_t &_fittedcharges_per_fragment,
 | 
|---|
| 223 |     const std::set<KeySet> &_fragments,
 | 
|---|
| 224 |     const AtomFragmentsMap &_atomfragments)
 | 
|---|
| 225 | {
 | 
|---|
| 226 |   for (std::set<KeySet>::const_iterator fragmentiter = _fragments.begin();
 | 
|---|
| 227 |       fragmentiter != _fragments.end(); ++fragmentiter) {
 | 
|---|
| 228 |     const KeySet &keyset = *fragmentiter;
 | 
|---|
| 229 |     const AtomFragmentsMap::indices_t &forceindices = _atomfragments.getFullKeyset(keyset);
 | 
|---|
| 230 |     ASSERT( !forceindices.empty(),
 | 
|---|
| 231 |         "getKeySetsToGraphMapping() - force keyset to "+toString(keyset)+" is empty.");
 | 
|---|
| 232 |     KeySet forcekeyset;
 | 
|---|
| 233 |     forcekeyset.insert(forceindices.begin(), forceindices.end());
 | 
|---|
| 234 |     forcekeyset.erase(-1);
 | 
|---|
| 235 |     const HomologyGraph graph(forcekeyset);
 | 
|---|
| 236 |     LOG(2, "DEBUG: Associating keyset " << forcekeyset << " with graph " << graph);
 | 
|---|
| 237 |     _keyset_graphs.insert( std::make_pair(keyset, graph) );
 | 
|---|
| 238 |     _fittedcharges_per_fragment.insert( std::make_pair(graph, PartialNucleiChargeFitter::charges_t()) );
 | 
|---|
| 239 |   }
 | 
|---|
| 240 | }
 | 
|---|
| 241 | 
 | 
|---|
| 242 | static
 | 
|---|
| 243 | bool getPartialChargesForAllGraphs(
 | 
|---|
| 244 |     detail::GraphFittedChargeMap_t &_fittedcharges_per_fragment,
 | 
|---|
| 245 |     const HomologyContainer &_homologies,
 | 
|---|
| 246 |     const double _mask_radius,
 | 
|---|
| 247 |     const bool enforceZeroCharge)
 | 
|---|
| 248 | {
 | 
|---|
| 249 |   for (detail::GraphFittedChargeMap_t::iterator graphiter = _fittedcharges_per_fragment.begin();
 | 
|---|
| 250 |       graphiter != _fittedcharges_per_fragment.end(); ++graphiter) {
 | 
|---|
| 251 |     const HomologyGraph &graph = graphiter->first;
 | 
|---|
| 252 |     LOG(2, "DEBUG: Now fitting charges to graph " << graph);
 | 
|---|
| 253 |     const HomologyContainer::range_t range = _homologies.getHomologousGraphs(graph);
 | 
|---|
| 254 |     if (range.first == range.second) {
 | 
|---|
| 255 |       ELOG(0, "HomologyContainer does not contain specified fragment.");
 | 
|---|
| 256 |       return false;
 | 
|---|
| 257 |     }
 | 
|---|
| 258 | 
 | 
|---|
| 259 |     // fit and average partial charges over all homologous fragments
 | 
|---|
| 260 |     PartialNucleiChargeFitter::charges_t &averaged_charges = graphiter->second;
 | 
|---|
| 261 |     const size_t NoFragments =
 | 
|---|
| 262 |         obtainAverageChargesOverFragments(averaged_charges, range, _mask_radius);
 | 
|---|
| 263 |     if ((NoFragments == 0) && (range.first != range.second)) {
 | 
|---|
| 264 |       ELOG(0, "Fitting for fragment "+toString(*graphiter)+" failed.");
 | 
|---|
| 265 |       return false;
 | 
|---|
| 266 |     }
 | 
|---|
| 267 | 
 | 
|---|
| 268 |     // make sum of charges zero if desired
 | 
|---|
| 269 |     if (enforceZeroCharge)
 | 
|---|
| 270 |       enforceZeroTotalCharge(averaged_charges);
 | 
|---|
| 271 | 
 | 
|---|
| 272 |     // output status info fitted charges
 | 
|---|
| 273 |     LOG(2, "DEBUG: For fragment " << *graphiter << " we have fitted the following charges "
 | 
|---|
| 274 |         << averaged_charges << ", averaged over " << NoFragments << " fragments.");
 | 
|---|
| 275 |   }
 | 
|---|
| 276 |   return true;
 | 
|---|
| 277 | }
 | 
|---|
| 278 | 
 | 
|---|
| 279 | const atom * getNonHydrogenSurrogate(const atom * const _walker)
 | 
|---|
| 280 | {
 | 
|---|
| 281 |   const atom * surrogate = _walker;
 | 
|---|
| 282 |   if (surrogate->getElementNo() == 1) {
 | 
|---|
| 283 |     // it's hydrogen, check its bonding and use its bond partner instead to request
 | 
|---|
| 284 |     // keysets
 | 
|---|
| 285 |     const BondList &ListOfBonds = surrogate->getListOfBonds();
 | 
|---|
| 286 |     if ( ListOfBonds.size() != 1) {
 | 
|---|
| 287 |       ELOG(1, "Solitary hydrogen in atom " << surrogate->getId() << " detected.");
 | 
|---|
| 288 |       return _walker;
 | 
|---|
| 289 |     }
 | 
|---|
| 290 |     surrogate = (*ListOfBonds.begin())->GetOtherAtom(surrogate);
 | 
|---|
| 291 |   }
 | 
|---|
| 292 |   return surrogate;
 | 
|---|
| 293 | }
 | 
|---|
| 294 | 
 | 
|---|
| 295 | double fitAverageChargeToAtom(
 | 
|---|
| 296 |     const atom * const _walker,
 | 
|---|
| 297 |     const AtomFragmentsMap &_atomfragments,
 | 
|---|
| 298 |     const detail::KeysetsToGraph_t &_keyset_graphs,
 | 
|---|
| 299 |     const detail::GraphFittedChargeMap_t &_fittedcharges_per_fragment)
 | 
|---|
| 300 | {
 | 
|---|
| 301 |   const atom * const surrogate = getNonHydrogenSurrogate(_walker);
 | 
|---|
| 302 |   const atomId_t walkerid = surrogate->getId();
 | 
|---|
| 303 |   const AtomFragmentsMap::AtomFragmentsMap_t &atommap = _atomfragments.getMap();
 | 
|---|
| 304 |   const AtomFragmentsMap::AtomFragmentsMap_t::const_iterator keysetsiter =
 | 
|---|
| 305 |       atommap.find(walkerid);
 | 
|---|
| 306 |   ASSERT(keysetsiter != atommap.end(),
 | 
|---|
| 307 |       "fitAverageChargeToAtom() - we checked already that "+toString(walkerid)
 | 
|---|
| 308 |       +" should be present!");
 | 
|---|
| 309 |   const AtomFragmentsMap::keysets_t & keysets = keysetsiter->second;
 | 
|---|
| 310 | 
 | 
|---|
| 311 |   double average_charge = 0.;
 | 
|---|
| 312 |   size_t NoFragments = 0;
 | 
|---|
| 313 |   // go over all fragments associated to this atom
 | 
|---|
| 314 |   for (AtomFragmentsMap::keysets_t::const_iterator keysetsiter = keysets.begin();
 | 
|---|
| 315 |       keysetsiter != keysets.end(); ++keysetsiter) {
 | 
|---|
| 316 |     const KeySet &keyset = *keysetsiter;
 | 
|---|
| 317 | 
 | 
|---|
| 318 |     const AtomFragmentsMap::indices_t &forcekeyset = _atomfragments.getFullKeyset(keyset);
 | 
|---|
| 319 |     ASSERT( !forcekeyset.empty(),
 | 
|---|
| 320 |         "fitAverageChargeToAtom() - force keyset to "+toString(keyset)+" is empty.");
 | 
|---|
| 321 | 
 | 
|---|
| 322 |     // find the associated charge in the charge vector
 | 
|---|
| 323 |     const std::map<KeySet, HomologyGraph>::const_iterator keysetgraphiter =
 | 
|---|
| 324 |         _keyset_graphs.find(keyset);
 | 
|---|
| 325 |     ASSERT( keysetgraphiter != _keyset_graphs.end(),
 | 
|---|
| 326 |         "fitAverageChargeToAtom() - keyset "+toString(keyset)
 | 
|---|
| 327 |         +" not contained in keyset_graphs.");
 | 
|---|
| 328 |     const HomologyGraph &graph = keysetgraphiter->second;
 | 
|---|
| 329 |     const detail::GraphFittedChargeMap_t::const_iterator chargesiter =
 | 
|---|
| 330 |         _fittedcharges_per_fragment.find(graph);
 | 
|---|
| 331 |     ASSERT(chargesiter != _fittedcharges_per_fragment.end(),
 | 
|---|
| 332 |         "fitAverageChargeToAtom() - no charge to "+toString(keyset)
 | 
|---|
| 333 |         +" any longer present in fittedcharges_per_fragment?");
 | 
|---|
| 334 |     const PartialNucleiChargeFitter::charges_t &charges = chargesiter->second;
 | 
|---|
| 335 |     ASSERT( charges.size() == forcekeyset.size(),
 | 
|---|
| 336 |         "fitAverageChargeToAtom() - charges "+toString(charges.size())+" and keyset "
 | 
|---|
| 337 |         +toString(forcekeyset.size())+" do not have the same length?");
 | 
|---|
| 338 |     PartialNucleiChargeFitter::charges_t::const_iterator chargeiter =
 | 
|---|
| 339 |         charges.begin();
 | 
|---|
| 340 |     const AtomFragmentsMap::indices_t::const_iterator forcekeysetiter =
 | 
|---|
| 341 |         std::find(forcekeyset.begin(), forcekeyset.end(), _walker->getId());
 | 
|---|
| 342 |     ASSERT( forcekeysetiter != forcekeyset.end(),
 | 
|---|
| 343 |         "fitAverageChargeToAtom() - atom "+toString(_walker->getId())
 | 
|---|
| 344 |         +" not contained in force keyset "+toString(forcekeyset));
 | 
|---|
| 345 |     std::advance(chargeiter, std::distance(forcekeyset.begin(), forcekeysetiter));
 | 
|---|
| 346 | 
 | 
|---|
| 347 |     // and add onto charge sum
 | 
|---|
| 348 |     const double & charge_in_fragment = *chargeiter;
 | 
|---|
| 349 |     average_charge += charge_in_fragment;
 | 
|---|
| 350 |     ++NoFragments;
 | 
|---|
| 351 |   }
 | 
|---|
| 352 |   // average to obtain final partial charge for this atom
 | 
|---|
| 353 |   average_charge *= 1./(double)NoFragments;
 | 
|---|
| 354 | 
 | 
|---|
| 355 |   return average_charge;
 | 
|---|
| 356 | }
 | 
|---|
| 357 | 
 | 
|---|
| 358 | void addToParticleRegistry(
 | 
|---|
| 359 |     const ParticleFactory &factory,
 | 
|---|
| 360 |     const periodentafel &periode,
 | 
|---|
| 361 |     const detail::fitted_charges_t &_fitted_charges,
 | 
|---|
| 362 |     const detail::GraphIndices_t &_GraphIndices,
 | 
|---|
| 363 |     detail::AtomParticleNames_t &_atom_particlenames)
 | 
|---|
| 364 | {
 | 
|---|
| 365 |   for (detail::fitted_charges_t::const_iterator chargeiter = _fitted_charges.begin();
 | 
|---|
| 366 |       chargeiter != _fitted_charges.end(); ++chargeiter) {
 | 
|---|
| 367 |     const atomId_t &atomid = chargeiter->first;
 | 
|---|
| 368 |     const atom * const walker = World::getInstance().getAtom(AtomById(atomid));
 | 
|---|
| 369 |     ASSERT( walker != NULL,
 | 
|---|
| 370 |         "addToParticleRegistry() - atom "+toString(atomid)
 | 
|---|
| 371 |         +" not present in the World?");
 | 
|---|
| 372 |     const detail::GraphIndices_t::right_const_iterator graphiter =
 | 
|---|
| 373 |         _GraphIndices.right.find(atomid);
 | 
|---|
| 374 |     ASSERT(graphiter != _GraphIndices.right.end(),
 | 
|---|
| 375 |         "addToParticleRegistry() - atom #"+toString(atomid)
 | 
|---|
| 376 |         +" not contained in GraphIndices.");
 | 
|---|
| 377 |     const detail::AtomParticleNames_t::iterator nameiter =
 | 
|---|
| 378 |         _atom_particlenames.find(graphiter->second);
 | 
|---|
| 379 |     const atomicNumber_t elementno = walker->getElementNo();
 | 
|---|
| 380 |     std::string name;
 | 
|---|
| 381 |     if ((nameiter != _atom_particlenames.end()) && (nameiter->second.count(elementno))) {
 | 
|---|
| 382 |         name = (nameiter->second)[elementno];
 | 
|---|
| 383 |     } else {
 | 
|---|
| 384 |       if (nameiter == _atom_particlenames.end())
 | 
|---|
| 385 |         _atom_particlenames.insert(
 | 
|---|
| 386 |             std::make_pair(graphiter->second, std::map<atomicNumber_t, std::string>()) );
 | 
|---|
| 387 |       const double &charge = chargeiter->second;
 | 
|---|
| 388 |       name = Particle::findFreeName(periode, elementno);
 | 
|---|
| 389 |       _atom_particlenames[graphiter->second][elementno] = name;
 | 
|---|
| 390 |       LOG(1, "INFO: Adding particle " << name << " for atom "
 | 
|---|
| 391 |           << *walker << " with element " << elementno << ", charge " << charge);
 | 
|---|
| 392 |       factory.createInstance(name, elementno, charge);
 | 
|---|
| 393 |     }
 | 
|---|
| 394 |   }
 | 
|---|
| 395 | }
 | 
|---|
| 396 | 
 | 
|---|
| 397 | bool isNotHydrogen(const atom * const _atom)
 | 
|---|
| 398 | {
 | 
|---|
| 399 |   return (_atom->getElementNo() != (atomicNumber_t) 1);
 | 
|---|
| 400 | }
 | 
|---|
| 401 | 
 | 
|---|
| 402 | ActionState::ptr PotentialFitPartialChargesAction::performCall()
 | 
|---|
| 403 | {
 | 
|---|
| 404 |   // check for selected atoms
 | 
|---|
| 405 |   const World &world = World::getConstInstance();
 | 
|---|
| 406 |   const std::vector<const atom *> selected_atoms = world.getSelectedAtoms();
 | 
|---|
| 407 |   if (selected_atoms.empty()) {
 | 
|---|
| 408 |     STATUS("There are no atoms selected for fitting partial charges to.");
 | 
|---|
| 409 |     return Action::failure;
 | 
|---|
| 410 |   }
 | 
|---|
| 411 | 
 | 
|---|
| 412 |   /// obtain possible fragments to each selected atom
 | 
|---|
| 413 |   const AtomFragmentsMap &atomfragments = AtomFragmentsMap::getConstInstance();
 | 
|---|
| 414 |   if (!atomfragments.checkCompleteness()) {
 | 
|---|
| 415 |     ELOG(0, "AtomFragmentsMap failed internal consistency check, missing forcekeysets?");
 | 
|---|
| 416 |     return Action::failure;
 | 
|---|
| 417 |   }
 | 
|---|
| 418 |   const std::set<KeySet> fragments =
 | 
|---|
| 419 |       accumulateKeySetsForAtoms( atomfragments.getMap(), selected_atoms);
 | 
|---|
| 420 |   const size_t NoNonHydrogens =
 | 
|---|
| 421 |       std::count_if(selected_atoms.begin(), selected_atoms.end(), isNotHydrogen);
 | 
|---|
| 422 |   if (fragments.size() < NoNonHydrogens) {
 | 
|---|
| 423 |     ELOG(0, "Obtained fewer fragments than there are atoms, has AtomFragments been loaded?");
 | 
|---|
| 424 |     return Action::failure;
 | 
|---|
| 425 |   }
 | 
|---|
| 426 | 
 | 
|---|
| 427 |   // reduce given fragments to homologous graphs to avoid multiple fittings
 | 
|---|
| 428 |   detail::KeysetsToGraph_t keyset_graphs;
 | 
|---|
| 429 |   detail::GraphFittedChargeMap_t fittedcharges_per_fragment;
 | 
|---|
| 430 |   getKeySetsToGraphMapping(keyset_graphs, fittedcharges_per_fragment, fragments, atomfragments);
 | 
|---|
| 431 | 
 | 
|---|
| 432 |   /// then go through all fragments and get partial charges for each
 | 
|---|
| 433 |   const HomologyContainer &homologies = World::getInstance().getHomologies();
 | 
|---|
| 434 |   const bool status = getPartialChargesForAllGraphs(
 | 
|---|
| 435 |       fittedcharges_per_fragment,
 | 
|---|
| 436 |       homologies,
 | 
|---|
| 437 |       params.radius.get(),
 | 
|---|
| 438 |       params.enforceZeroCharge.get());
 | 
|---|
| 439 |   if (!status)
 | 
|---|
| 440 |     return Action::failure;
 | 
|---|
| 441 | 
 | 
|---|
| 442 |   /// obtain average charge for each atom the fitted charges over all its fragments
 | 
|---|
| 443 |   detail::fitted_charges_t fitted_charges;
 | 
|---|
| 444 |   for (std::vector<const atom *>::const_iterator atomiter = selected_atoms.begin();
 | 
|---|
| 445 |       atomiter != selected_atoms.end(); ++atomiter) {
 | 
|---|
| 446 |     const atomId_t walkerid = (*atomiter)->getId();
 | 
|---|
| 447 |     const double average_charge = fitAverageChargeToAtom(
 | 
|---|
| 448 |         *atomiter, atomfragments, keyset_graphs, fittedcharges_per_fragment);
 | 
|---|
| 449 | 
 | 
|---|
| 450 |     if (average_charge != 0.) {
 | 
|---|
| 451 |       LOG(2, "DEBUG: For atom " << *atomiter << " we have an average charge of "
 | 
|---|
| 452 |           << average_charge);
 | 
|---|
| 453 | 
 | 
|---|
| 454 |       fitted_charges.insert( std::make_pair(walkerid, average_charge) );
 | 
|---|
| 455 |     }
 | 
|---|
| 456 |   }
 | 
|---|
| 457 | 
 | 
|---|
| 458 |   /// make Particles be used for every atom that was fitted on the same number of graphs
 | 
|---|
| 459 |   detail::GraphIndex_t GraphIndex;
 | 
|---|
| 460 |   size_t index = 0;
 | 
|---|
| 461 |   for (HomologyContainer::const_key_iterator iter = homologies.key_begin();
 | 
|---|
| 462 |       iter != homologies.key_end(); iter = homologies.getNextKey(iter)) {
 | 
|---|
| 463 |     GraphIndex.insert( std::make_pair( *iter, index++));
 | 
|---|
| 464 |   }
 | 
|---|
| 465 |   LOG(2, "DEBUG: There are " << index << " unique graphs in the homology container.");
 | 
|---|
| 466 | 
 | 
|---|
| 467 |   // go through every non-hydrogen atom, get all graphs, convert to GraphIndex and store
 | 
|---|
| 468 |   detail::GraphIndices_t GraphIndices;
 | 
|---|
| 469 |   const AtomFragmentsMap::AtomFragmentsMap_t &atommap = atomfragments.getMap();
 | 
|---|
| 470 |   for (std::vector<const atom *>::const_iterator atomiter = selected_atoms.begin();
 | 
|---|
| 471 |       atomiter != selected_atoms.end(); ++atomiter) {
 | 
|---|
| 472 |     // use the non-hydrogen here
 | 
|---|
| 473 |     const atomId_t walkerid = (*atomiter)->getId();
 | 
|---|
| 474 |     const atomId_t surrogateid = getNonHydrogenSurrogate(*atomiter)->getId();
 | 
|---|
| 475 |     if (surrogateid != walkerid)
 | 
|---|
| 476 |       continue;
 | 
|---|
| 477 |     const AtomFragmentsMap::AtomFragmentsMap_t::const_iterator keysetsiter =
 | 
|---|
| 478 |         atommap.find(walkerid);
 | 
|---|
| 479 |     ASSERT(keysetsiter != atommap.end(),
 | 
|---|
| 480 |         "PotentialFitPartialChargesAction::performCall() - we checked already that "
 | 
|---|
| 481 |         +toString(surrogateid)+" should be present!");
 | 
|---|
| 482 |     const AtomFragmentsMap::keysets_t & keysets = keysetsiter->second;
 | 
|---|
| 483 | 
 | 
|---|
| 484 |     // go over all fragments associated to this atom
 | 
|---|
| 485 |     detail::AtomsGraphIndices_t AtomsGraphIndices;
 | 
|---|
| 486 |     for (AtomFragmentsMap::keysets_t::const_iterator keysetsiter = keysets.begin();
 | 
|---|
| 487 |         keysetsiter != keysets.end(); ++keysetsiter) {
 | 
|---|
| 488 |       const KeySet &keyset = *keysetsiter;
 | 
|---|
| 489 |       const std::map<KeySet, HomologyGraph>::const_iterator keysetgraphiter =
 | 
|---|
| 490 |           keyset_graphs.find(keyset);
 | 
|---|
| 491 |       ASSERT( keysetgraphiter != keyset_graphs.end(),
 | 
|---|
| 492 |           "PotentialFitPartialChargesAction::performCall() - keyset "+toString(keyset)
 | 
|---|
| 493 |           +" not contained in keyset_graphs.");
 | 
|---|
| 494 |       const HomologyGraph &graph = keysetgraphiter->second;
 | 
|---|
| 495 |       const detail::GraphIndex_t::const_iterator indexiter = GraphIndex.find(graph);
 | 
|---|
| 496 |       ASSERT( indexiter != GraphIndex.end(),
 | 
|---|
| 497 |           "PotentialFitPartialChargesAction::performCall() - graph "+toString(graph)
 | 
|---|
| 498 |           +" not contained in GraphIndex.");
 | 
|---|
| 499 |       AtomsGraphIndices.insert( indexiter->second );
 | 
|---|
| 500 |     }
 | 
|---|
| 501 | 
 | 
|---|
| 502 |     GraphIndices.insert( detail::GraphIndices_t::value_type(AtomsGraphIndices, walkerid) );
 | 
|---|
| 503 | 
 | 
|---|
| 504 |     LOG(2, "DEBUG: Atom #" << walkerid << "," << *atomiter << ". has graph indices "
 | 
|---|
| 505 |         << AtomsGraphIndices);
 | 
|---|
| 506 |   }
 | 
|---|
| 507 |   // then graphs from non-hydrogen bond partner for all hydrogens
 | 
|---|
| 508 |   for (std::vector<const atom *>::const_iterator atomiter = selected_atoms.begin();
 | 
|---|
| 509 |       atomiter != selected_atoms.end(); ++atomiter) {
 | 
|---|
| 510 |     // use the non-hydrogen here
 | 
|---|
| 511 |     const atomId_t walkerid = (*atomiter)->getId();
 | 
|---|
| 512 |     const atomId_t surrogateid = getNonHydrogenSurrogate((*atomiter))->getId();
 | 
|---|
| 513 |     if (surrogateid == walkerid)
 | 
|---|
| 514 |       continue;
 | 
|---|
| 515 |     detail::GraphIndices_t::right_const_iterator graphiter = GraphIndices.right.find(surrogateid);
 | 
|---|
| 516 |     ASSERT( graphiter != GraphIndices.right.end(),
 | 
|---|
| 517 |         "PotentialFitPartialChargesAction::performCall() - atom #"+toString(surrogateid)
 | 
|---|
| 518 |         +" not contained in GraphIndices.");
 | 
|---|
| 519 |     const detail::AtomsGraphIndices_t &AtomsGraphIndices = graphiter->second;
 | 
|---|
| 520 |     GraphIndices.insert( detail::GraphIndices_t::value_type(AtomsGraphIndices, walkerid) );
 | 
|---|
| 521 |     LOG(2, "DEBUG: Hydrogen #" << walkerid << ", " << *atomiter
 | 
|---|
| 522 |         << ", has graph indices " << AtomsGraphIndices);
 | 
|---|
| 523 |   }
 | 
|---|
| 524 | 
 | 
|---|
| 525 |   /// place all fitted charges into ParticleRegistry
 | 
|---|
| 526 |   detail::AtomParticleNames_t atom_particlenames;
 | 
|---|
| 527 |   addToParticleRegistry(
 | 
|---|
| 528 |       ParticleFactory::getConstInstance(),
 | 
|---|
| 529 |       *World::getInstance().getPeriode(),
 | 
|---|
| 530 |       fitted_charges,
 | 
|---|
| 531 |       GraphIndices,
 | 
|---|
| 532 |       atom_particlenames);
 | 
|---|
| 533 |   for (World::AtomSelectionIterator atomiter = World::getInstance().beginAtomSelection();
 | 
|---|
| 534 |       atomiter != World::getInstance().endAtomSelection(); ++atomiter) {
 | 
|---|
| 535 |     atom * const walker = atomiter->second;
 | 
|---|
| 536 |     const atomId_t walkerid = atomiter->first;
 | 
|---|
| 537 |     const detail::GraphIndices_t::right_const_iterator graphiter =
 | 
|---|
| 538 |         GraphIndices.right.find(walkerid);
 | 
|---|
| 539 |     ASSERT( graphiter != GraphIndices.right.end(),
 | 
|---|
| 540 |         "PotentialFitPartialChargesAction::performCall() - cannot find "
 | 
|---|
| 541 |         +toString(walkerid)+" in GraphIndices.");
 | 
|---|
| 542 |     const detail::AtomsGraphIndices_t &graphindex = graphiter->second;
 | 
|---|
| 543 |     const detail::AtomParticleNames_t::const_iterator particlesetiter =
 | 
|---|
| 544 |         atom_particlenames.find(graphindex);
 | 
|---|
| 545 |     ASSERT( particlesetiter != atom_particlenames.end(),
 | 
|---|
| 546 |         "PotentialFitPartialChargesAction::performCall() - cannot find "
 | 
|---|
| 547 |         +toString(graphindex)+" in atom_particlenames.");
 | 
|---|
| 548 |     const std::map<atomicNumber_t, std::string>::const_iterator nameiter =
 | 
|---|
| 549 |         particlesetiter->second.find(walker->getElementNo());
 | 
|---|
| 550 |     ASSERT( nameiter != particlesetiter->second.end(),
 | 
|---|
| 551 |         "PotentialFitPartialChargesAction::performCall() - ");
 | 
|---|
| 552 |     walker->setParticleName(nameiter->second);
 | 
|---|
| 553 |     LOG(1, "INFO: atom " << *walker << " received the following particle "
 | 
|---|
| 554 |         << walker->getParticleName());
 | 
|---|
| 555 |   }
 | 
|---|
| 556 | 
 | 
|---|
| 557 |   return Action::success;
 | 
|---|
| 558 | }
 | 
|---|
| 559 | 
 | 
|---|
| 560 | ActionState::ptr PotentialFitPartialChargesAction::performUndo(ActionState::ptr _state) {
 | 
|---|
| 561 |   return Action::success;
 | 
|---|
| 562 | }
 | 
|---|
| 563 | 
 | 
|---|
| 564 | ActionState::ptr PotentialFitPartialChargesAction::performRedo(ActionState::ptr _state){
 | 
|---|
| 565 |   return Action::success;
 | 
|---|
| 566 | }
 | 
|---|
| 567 | 
 | 
|---|
| 568 | bool PotentialFitPartialChargesAction::canUndo() {
 | 
|---|
| 569 |   return false;
 | 
|---|
| 570 | }
 | 
|---|
| 571 | 
 | 
|---|
| 572 | bool PotentialFitPartialChargesAction::shouldUndo() {
 | 
|---|
| 573 |   return false;
 | 
|---|
| 574 | }
 | 
|---|
| 575 | /** =========== end of function ====================== */
 | 
|---|