| [c4a323] | 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 | /*
 | 
|---|
| [50d49d] | 24 |  * FitPartialChargesAction.cpp
 | 
|---|
| [c4a323] | 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 | 
 | 
|---|
| [9eb71b3] | 38 | //#include "CodePatterns/MemDebug.hpp"
 | 
|---|
| [c4a323] | 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 | 
 | 
|---|
| [01120c] | 46 | #include <boost/bimap.hpp>
 | 
|---|
| [91f907] | 47 | #include <boost/bimap/multiset_of.hpp>
 | 
|---|
| [a9099d] | 48 | #include <boost/bind.hpp>
 | 
|---|
| [c4a323] | 49 | #include <boost/filesystem.hpp>
 | 
|---|
 | 50 | #include <boost/foreach.hpp>
 | 
|---|
 | 51 | #include <algorithm>
 | 
|---|
| [9b68fc] | 52 | #include <functional>
 | 
|---|
| [c4a323] | 53 | #include <iostream>
 | 
|---|
 | 54 | #include <string>
 | 
|---|
 | 55 | 
 | 
|---|
| [50d49d] | 56 | #include "Actions/PotentialAction/FitPartialChargesAction.hpp"
 | 
|---|
| [c4a323] | 57 | 
 | 
|---|
 | 58 | #include "Potentials/PartialNucleiChargeFitter.hpp"
 | 
|---|
 | 59 | 
 | 
|---|
| [01120c] | 60 | #include "AtomIdSet.hpp"
 | 
|---|
 | 61 | #include "Descriptors/AtomIdDescriptor.hpp"
 | 
|---|
| [c4a323] | 62 | #include "Element/element.hpp"
 | 
|---|
| [2082637] | 63 | #include "Element/periodentafel.hpp"
 | 
|---|
| [d8ed62] | 64 | #include "Fragmentation/Homology/AtomFragmentsMap.hpp"
 | 
|---|
| [c4a323] | 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"
 | 
|---|
| [708ec1] | 70 | #include "Potentials/Particles/ParticleFactory.hpp"
 | 
|---|
 | 71 | #include "Potentials/Particles/ParticleRegistry.hpp"
 | 
|---|
| [c4a323] | 72 | #include "Potentials/SerializablePotential.hpp"
 | 
|---|
 | 73 | #include "World.hpp"
 | 
|---|
 | 74 | 
 | 
|---|
 | 75 | using namespace MoleCuilder;
 | 
|---|
 | 76 | 
 | 
|---|
 | 77 | // and construct the stuff
 | 
|---|
| [50d49d] | 78 | #include "FitPartialChargesAction.def"
 | 
|---|
| [c4a323] | 79 | #include "Action_impl_pre.hpp"
 | 
|---|
 | 80 | /** =========== define the function ====================== */
 | 
|---|
 | 81 | 
 | 
|---|
| [91f907] | 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<
 | 
|---|
| [ed0f88] | 93 |       boost::bimaps::multiset_of<AtomsGraphIndices_t>,
 | 
|---|
 | 94 |       atomId_t > GraphIndices_t;
 | 
|---|
| [91f907] | 95 | 
 | 
|---|
| [ed0f88] | 96 |   typedef std::map<std::set<size_t>, std::map<atomicNumber_t, std::string> > AtomParticleNames_t;
 | 
|---|
| [91f907] | 97 | 
 | 
|---|
 | 98 |   typedef std::map<std::set<size_t>, std::string> GraphToNameMap_t;
 | 
|---|
 | 99 | };
 | 
|---|
 | 100 | 
 | 
|---|
| [01120c] | 101 | static void enforceZeroTotalCharge(
 | 
|---|
| [d8ed62] | 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()));
 | 
|---|
| [c4a323] | 111 |   }
 | 
|---|
| [d8ed62] | 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.");
 | 
|---|
| [01120c] | 117 | 
 | 
|---|
 | 118 |   LOG(2, "DEBUG: final charges with net zero charge are " << _averaged_charges);
 | 
|---|
| [d8ed62] | 119 | }
 | 
|---|
| [9b68fc] | 120 | 
 | 
|---|
| [01120c] | 121 | static size_t obtainAverageChargesOverFragments(
 | 
|---|
| [d8ed62] | 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;
 | 
|---|
| [c4a323] | 130 |   if (!iter->second.containsGrids) {
 | 
|---|
| [d8ed62] | 131 |     ELOG(1, "This HomologyGraph does not contain sampled grids.");
 | 
|---|
 | 132 |     return 0;
 | 
|---|
| [c4a323] | 133 |   }
 | 
|---|
| [d8ed62] | 134 |   _average_charges.resize(iter->second.fragment.getCharges().size(), 0.);
 | 
|---|
| [9b68fc] | 135 |   size_t NoFragments = 0;
 | 
|---|
 | 136 |   for (;
 | 
|---|
| [d8ed62] | 137 |       iter != _range.second; ++iter, ++NoFragments) {
 | 
|---|
| [9b68fc] | 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.");
 | 
|---|
| [d8ed62] | 140 |       return 0;
 | 
|---|
| [9b68fc] | 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.");
 | 
|---|
| [d8ed62] | 151 |       return 0;
 | 
|---|
| [9b68fc] | 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 |     }
 | 
|---|
| [d8ed62] | 161 |     PartialNucleiChargeFitter fitter(potential, positions, _radius);
 | 
|---|
| [9b68fc] | 162 |     fitter();
 | 
|---|
 | 163 |     PartialNucleiChargeFitter::charges_t return_charges = fitter.getSolutionAsCharges_t();
 | 
|---|
| [d8ed62] | 164 |     LOG(2, "DEBUG: fitted charges are " << return_charges);
 | 
|---|
| [9b68fc] | 165 |     std::transform(
 | 
|---|
 | 166 |         return_charges.begin(), return_charges.end(),
 | 
|---|
| [d8ed62] | 167 |         _average_charges.begin(),
 | 
|---|
 | 168 |         _average_charges.begin(),
 | 
|---|
| [9b68fc] | 169 |         std::plus<PartialNucleiChargeFitter::charge_t>());
 | 
|---|
| [c4a323] | 170 |   }
 | 
|---|
| [d8ed62] | 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 | 
 | 
|---|
| [01120c] | 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 | 
 | 
|---|
| [c1ec8e] | 193 | static
 | 
|---|
 | 194 | std::set<KeySet> accumulateKeySetsForAtoms(
 | 
|---|
 | 195 |     const AtomFragmentsMap::AtomFragmentsMap_t &_atommap,
 | 
|---|
 | 196 |     const std::vector<const atom *> &_selected_atoms)
 | 
|---|
| [d8ed62] | 197 | {
 | 
|---|
| [01120c] | 198 |   std::set<KeySet> fragments;
 | 
|---|
| [c1ec8e] | 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();
 | 
|---|
| [d8ed62] | 202 |     const AtomFragmentsMap::AtomFragmentsMap_t::const_iterator atomiter =
 | 
|---|
| [c1ec8e] | 203 |         _atommap.find(atomid);
 | 
|---|
 | 204 |     if ((*iter)->getElementNo() != 1) {
 | 
|---|
 | 205 |       if (atomiter == _atommap.end()) {
 | 
|---|
 | 206 |         ELOG(2, "There are no fragments associated to " << atomid << ".");
 | 
|---|
| [01120c] | 207 |         continue;
 | 
|---|
 | 208 |       }
 | 
|---|
 | 209 |       const AtomFragmentsMap::keysets_t &keysets = atomiter->second;
 | 
|---|
| [c1ec8e] | 210 |       LOG(2, "DEBUG: atom " << atomid << " has " << keysets.size() << " fragments.");
 | 
|---|
| [01120c] | 211 |       fragments.insert( keysets.begin(), keysets.end() );
 | 
|---|
 | 212 |     } else {
 | 
|---|
| [c1ec8e] | 213 |       LOG(3, "DEBUG: Skipping atom " << atomid << " as it's hydrogen.");
 | 
|---|
| [a9099d] | 214 |     }
 | 
|---|
 | 215 |   }
 | 
|---|
| [c1ec8e] | 216 |   return fragments;
 | 
|---|
 | 217 | }
 | 
|---|
| [a9099d] | 218 | 
 | 
|---|
| [c1ec8e] | 219 | static
 | 
|---|
 | 220 | void getKeySetsToGraphMapping(
 | 
|---|
| [91f907] | 221 |     detail::KeysetsToGraph_t &_keyset_graphs,
 | 
|---|
 | 222 |     detail::GraphFittedChargeMap_t &_fittedcharges_per_fragment,
 | 
|---|
| [c1ec8e] | 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) {
 | 
|---|
| [01120c] | 228 |     const KeySet &keyset = *fragmentiter;
 | 
|---|
| [c1ec8e] | 229 |     const AtomFragmentsMap::indices_t &forceindices = _atomfragments.getFullKeyset(keyset);
 | 
|---|
| [01120c] | 230 |     ASSERT( !forceindices.empty(),
 | 
|---|
| [c1ec8e] | 231 |         "getKeySetsToGraphMapping() - force keyset to "+toString(keyset)+" is empty.");
 | 
|---|
| [01120c] | 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);
 | 
|---|
| [c1ec8e] | 237 |     _keyset_graphs.insert( std::make_pair(keyset, graph) );
 | 
|---|
 | 238 |     _fittedcharges_per_fragment.insert( std::make_pair(graph, PartialNucleiChargeFitter::charges_t()) );
 | 
|---|
| [01120c] | 239 |   }
 | 
|---|
| [c1ec8e] | 240 | }
 | 
|---|
| [d8ed62] | 241 | 
 | 
|---|
| [c1ec8e] | 242 | static
 | 
|---|
 | 243 | bool getPartialChargesForAllGraphs(
 | 
|---|
| [91f907] | 244 |     detail::GraphFittedChargeMap_t &_fittedcharges_per_fragment,
 | 
|---|
| [c1ec8e] | 245 |     const HomologyContainer &_homologies,
 | 
|---|
 | 246 |     const double _mask_radius,
 | 
|---|
 | 247 |     const bool enforceZeroCharge)
 | 
|---|
 | 248 | {
 | 
|---|
| [91f907] | 249 |   for (detail::GraphFittedChargeMap_t::iterator graphiter = _fittedcharges_per_fragment.begin();
 | 
|---|
| [c1ec8e] | 250 |       graphiter != _fittedcharges_per_fragment.end(); ++graphiter) {
 | 
|---|
| [01120c] | 251 |     const HomologyGraph &graph = graphiter->first;
 | 
|---|
 | 252 |     LOG(2, "DEBUG: Now fitting charges to graph " << graph);
 | 
|---|
| [c1ec8e] | 253 |     const HomologyContainer::range_t range = _homologies.getHomologousGraphs(graph);
 | 
|---|
| [d8ed62] | 254 |     if (range.first == range.second) {
 | 
|---|
| [c1ec8e] | 255 |       ELOG(0, "HomologyContainer does not contain specified fragment.");
 | 
|---|
 | 256 |       return false;
 | 
|---|
| [d8ed62] | 257 |     }
 | 
|---|
 | 258 | 
 | 
|---|
 | 259 |     // fit and average partial charges over all homologous fragments
 | 
|---|
| [01120c] | 260 |     PartialNucleiChargeFitter::charges_t &averaged_charges = graphiter->second;
 | 
|---|
| [d8ed62] | 261 |     const size_t NoFragments =
 | 
|---|
| [c1ec8e] | 262 |         obtainAverageChargesOverFragments(averaged_charges, range, _mask_radius);
 | 
|---|
| [d8ed62] | 263 |     if ((NoFragments == 0) && (range.first != range.second)) {
 | 
|---|
| [c1ec8e] | 264 |       ELOG(0, "Fitting for fragment "+toString(*graphiter)+" failed.");
 | 
|---|
 | 265 |       return false;
 | 
|---|
| [d8ed62] | 266 |     }
 | 
|---|
 | 267 | 
 | 
|---|
 | 268 |     // make sum of charges zero if desired
 | 
|---|
| [c1ec8e] | 269 |     if (enforceZeroCharge)
 | 
|---|
| [d8ed62] | 270 |       enforceZeroTotalCharge(averaged_charges);
 | 
|---|
 | 271 | 
 | 
|---|
 | 272 |     // output status info fitted charges
 | 
|---|
| [01120c] | 273 |     LOG(2, "DEBUG: For fragment " << *graphiter << " we have fitted the following charges "
 | 
|---|
| [d8ed62] | 274 |         << averaged_charges << ", averaged over " << NoFragments << " fragments.");
 | 
|---|
 | 275 |   }
 | 
|---|
| [c1ec8e] | 276 |   return true;
 | 
|---|
 | 277 | }
 | 
|---|
| [d8ed62] | 278 | 
 | 
|---|
| [ed0f88] | 279 | const atom * getNonHydrogenSurrogate(const atom * const _walker)
 | 
|---|
| [c1ec8e] | 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) {
 | 
|---|
| [ed0f88] | 287 |       ELOG(1, "Solitary hydrogen in atom " << surrogate->getId() << " detected.");
 | 
|---|
 | 288 |       return _walker;
 | 
|---|
| [d8ed62] | 289 |     }
 | 
|---|
| [c1ec8e] | 290 |     surrogate = (*ListOfBonds.begin())->GetOtherAtom(surrogate);
 | 
|---|
 | 291 |   }
 | 
|---|
| [ed0f88] | 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);
 | 
|---|
| [c1ec8e] | 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
 | 
|---|
| [91f907] | 323 |     const std::map<KeySet, HomologyGraph>::const_iterator keysetgraphiter =
 | 
|---|
| [c1ec8e] | 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;
 | 
|---|
| [91f907] | 329 |     const detail::GraphFittedChargeMap_t::const_iterator chargesiter =
 | 
|---|
| [c1ec8e] | 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;
 | 
|---|
| [d8ed62] | 351 |   }
 | 
|---|
| [c1ec8e] | 352 |   // average to obtain final partial charge for this atom
 | 
|---|
 | 353 |   average_charge *= 1./(double)NoFragments;
 | 
|---|
| [d8ed62] | 354 | 
 | 
|---|
| [c1ec8e] | 355 |   return average_charge;
 | 
|---|
 | 356 | }
 | 
|---|
 | 357 | 
 | 
|---|
 | 358 | void addToParticleRegistry(
 | 
|---|
 | 359 |     const ParticleFactory &factory,
 | 
|---|
 | 360 |     const periodentafel &periode,
 | 
|---|
| [91f907] | 361 |     const detail::fitted_charges_t &_fitted_charges,
 | 
|---|
 | 362 |     const detail::GraphIndices_t &_GraphIndices,
 | 
|---|
 | 363 |     detail::AtomParticleNames_t &_atom_particlenames)
 | 
|---|
| [c1ec8e] | 364 | {
 | 
|---|
| [91f907] | 365 |   for (detail::fitted_charges_t::const_iterator chargeiter = _fitted_charges.begin();
 | 
|---|
| [c1ec8e] | 366 |       chargeiter != _fitted_charges.end(); ++chargeiter) {
 | 
|---|
| [91f907] | 367 |     const atomId_t &atomid = chargeiter->first;
 | 
|---|
| [ed0f88] | 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();
 | 
|---|
| [91f907] | 380 |     std::string name;
 | 
|---|
| [ed0f88] | 381 |     if ((nameiter != _atom_particlenames.end()) && (nameiter->second.count(elementno))) {
 | 
|---|
 | 382 |         name = (nameiter->second)[elementno];
 | 
|---|
| [91f907] | 383 |     } else {
 | 
|---|
| [ed0f88] | 384 |       if (nameiter == _atom_particlenames.end())
 | 
|---|
 | 385 |         _atom_particlenames.insert(
 | 
|---|
 | 386 |             std::make_pair(graphiter->second, std::map<atomicNumber_t, std::string>()) );
 | 
|---|
| [91f907] | 387 |       const double &charge = chargeiter->second;
 | 
|---|
 | 388 |       name = Particle::findFreeName(periode, elementno);
 | 
|---|
| [ed0f88] | 389 |       _atom_particlenames[graphiter->second][elementno] = name;
 | 
|---|
| [91f907] | 390 |       LOG(1, "INFO: Adding particle " << name << " for atom "
 | 
|---|
 | 391 |           << *walker << " with element " << elementno << ", charge " << charge);
 | 
|---|
 | 392 |       factory.createInstance(name, elementno, charge);
 | 
|---|
 | 393 |     }
 | 
|---|
| [708ec1] | 394 |   }
 | 
|---|
| [c1ec8e] | 395 | }
 | 
|---|
 | 396 | 
 | 
|---|
| [f33ef9] | 397 | bool isNotHydrogen(const atom * const _atom)
 | 
|---|
 | 398 | {
 | 
|---|
 | 399 |   return (_atom->getElementNo() != (atomicNumber_t) 1);
 | 
|---|
 | 400 | }
 | 
|---|
 | 401 | 
 | 
|---|
| [c1ec8e] | 402 | ActionState::ptr PotentialFitPartialChargesAction::performCall()
 | 
|---|
 | 403 | {
 | 
|---|
 | 404 |   // check for selected atoms
 | 
|---|
| [f20b4b] | 405 |   const World &world = World::getConstInstance();
 | 
|---|
| [f33ef9] | 406 |   const std::vector<const atom *> selected_atoms = world.getSelectedAtoms();
 | 
|---|
 | 407 |   if (selected_atoms.empty()) {
 | 
|---|
| [c1ec8e] | 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()) {
 | 
|---|
| [f33ef9] | 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?");
 | 
|---|
| [c1ec8e] | 424 |     return Action::failure;
 | 
|---|
 | 425 |   }
 | 
|---|
 | 426 | 
 | 
|---|
 | 427 |   // reduce given fragments to homologous graphs to avoid multiple fittings
 | 
|---|
| [91f907] | 428 |   detail::KeysetsToGraph_t keyset_graphs;
 | 
|---|
 | 429 |   detail::GraphFittedChargeMap_t fittedcharges_per_fragment;
 | 
|---|
| [c1ec8e] | 430 |   getKeySetsToGraphMapping(keyset_graphs, fittedcharges_per_fragment, fragments, atomfragments);
 | 
|---|
 | 431 | 
 | 
|---|
 | 432 |   /// then go through all fragments and get partial charges for each
 | 
|---|
| [91f907] | 433 |   const HomologyContainer &homologies = World::getInstance().getHomologies();
 | 
|---|
| [c1ec8e] | 434 |   const bool status = getPartialChargesForAllGraphs(
 | 
|---|
 | 435 |       fittedcharges_per_fragment,
 | 
|---|
| [91f907] | 436 |       homologies,
 | 
|---|
| [c1ec8e] | 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
 | 
|---|
| [91f907] | 443 |   detail::fitted_charges_t fitted_charges;
 | 
|---|
| [f33ef9] | 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();
 | 
|---|
| [c1ec8e] | 447 |     const double average_charge = fitAverageChargeToAtom(
 | 
|---|
| [f33ef9] | 448 |         *atomiter, atomfragments, keyset_graphs, fittedcharges_per_fragment);
 | 
|---|
| [c1ec8e] | 449 | 
 | 
|---|
 | 450 |     if (average_charge != 0.) {
 | 
|---|
| [f5dbea] | 451 |       LOG(2, "DEBUG: For atom " << **atomiter << " we have an average charge of "
 | 
|---|
| [c1ec8e] | 452 |           << average_charge);
 | 
|---|
 | 453 | 
 | 
|---|
| [f33ef9] | 454 |       fitted_charges.insert( std::make_pair(walkerid, average_charge) );
 | 
|---|
| [c1ec8e] | 455 |     }
 | 
|---|
 | 456 |   }
 | 
|---|
 | 457 | 
 | 
|---|
| [91f907] | 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 | 
 | 
|---|
| [ed0f88] | 467 |   // go through every non-hydrogen atom, get all graphs, convert to GraphIndex and store
 | 
|---|
| [91f907] | 468 |   detail::GraphIndices_t GraphIndices;
 | 
|---|
 | 469 |   const AtomFragmentsMap::AtomFragmentsMap_t &atommap = atomfragments.getMap();
 | 
|---|
| [f33ef9] | 470 |   for (std::vector<const atom *>::const_iterator atomiter = selected_atoms.begin();
 | 
|---|
 | 471 |       atomiter != selected_atoms.end(); ++atomiter) {
 | 
|---|
| [ed0f88] | 472 |     // use the non-hydrogen here
 | 
|---|
| [f33ef9] | 473 |     const atomId_t walkerid = (*atomiter)->getId();
 | 
|---|
 | 474 |     const atomId_t surrogateid = getNonHydrogenSurrogate(*atomiter)->getId();
 | 
|---|
 | 475 |     if (surrogateid != walkerid)
 | 
|---|
| [ed0f88] | 476 |       continue;
 | 
|---|
| [91f907] | 477 |     const AtomFragmentsMap::AtomFragmentsMap_t::const_iterator keysetsiter =
 | 
|---|
| [f33ef9] | 478 |         atommap.find(walkerid);
 | 
|---|
| [91f907] | 479 |     ASSERT(keysetsiter != atommap.end(),
 | 
|---|
 | 480 |         "PotentialFitPartialChargesAction::performCall() - we checked already that "
 | 
|---|
| [f33ef9] | 481 |         +toString(surrogateid)+" should be present!");
 | 
|---|
| [91f907] | 482 |     const AtomFragmentsMap::keysets_t & keysets = keysetsiter->second;
 | 
|---|
 | 483 | 
 | 
|---|
 | 484 |     // go over all fragments associated to this atom
 | 
|---|
| [ed0f88] | 485 |     detail::AtomsGraphIndices_t AtomsGraphIndices;
 | 
|---|
| [91f907] | 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 |     }
 | 
|---|
| [ed0f88] | 501 | 
 | 
|---|
| [f33ef9] | 502 |     GraphIndices.insert( detail::GraphIndices_t::value_type(AtomsGraphIndices, walkerid) );
 | 
|---|
| [ed0f88] | 503 | 
 | 
|---|
| [f5dbea] | 504 |     LOG(2, "DEBUG: Atom #" << walkerid << "," << **atomiter << ". has graph indices "
 | 
|---|
| [f33ef9] | 505 |         << AtomsGraphIndices);
 | 
|---|
| [ed0f88] | 506 |   }
 | 
|---|
 | 507 |   // then graphs from non-hydrogen bond partner for all hydrogens
 | 
|---|
| [f33ef9] | 508 |   for (std::vector<const atom *>::const_iterator atomiter = selected_atoms.begin();
 | 
|---|
 | 509 |       atomiter != selected_atoms.end(); ++atomiter) {
 | 
|---|
| [ed0f88] | 510 |     // use the non-hydrogen here
 | 
|---|
| [f33ef9] | 511 |     const atomId_t walkerid = (*atomiter)->getId();
 | 
|---|
 | 512 |     const atomId_t surrogateid = getNonHydrogenSurrogate((*atomiter))->getId();
 | 
|---|
 | 513 |     if (surrogateid == walkerid)
 | 
|---|
| [ed0f88] | 514 |       continue;
 | 
|---|
| [f33ef9] | 515 |     detail::GraphIndices_t::right_const_iterator graphiter = GraphIndices.right.find(surrogateid);
 | 
|---|
| [ed0f88] | 516 |     ASSERT( graphiter != GraphIndices.right.end(),
 | 
|---|
| [f33ef9] | 517 |         "PotentialFitPartialChargesAction::performCall() - atom #"+toString(surrogateid)
 | 
|---|
| [ed0f88] | 518 |         +" not contained in GraphIndices.");
 | 
|---|
 | 519 |     const detail::AtomsGraphIndices_t &AtomsGraphIndices = graphiter->second;
 | 
|---|
| [f33ef9] | 520 |     GraphIndices.insert( detail::GraphIndices_t::value_type(AtomsGraphIndices, walkerid) );
 | 
|---|
 | 521 |     LOG(2, "DEBUG: Hydrogen #" << walkerid << ", " << *atomiter
 | 
|---|
| [ed0f88] | 522 |         << ", has graph indices " << AtomsGraphIndices);
 | 
|---|
| [91f907] | 523 |   }
 | 
|---|
 | 524 | 
 | 
|---|
| [c1ec8e] | 525 |   /// place all fitted charges into ParticleRegistry
 | 
|---|
| [91f907] | 526 |   detail::AtomParticleNames_t atom_particlenames;
 | 
|---|
| [c1ec8e] | 527 |   addToParticleRegistry(
 | 
|---|
 | 528 |       ParticleFactory::getConstInstance(),
 | 
|---|
 | 529 |       *World::getInstance().getPeriode(),
 | 
|---|
| [91f907] | 530 |       fitted_charges,
 | 
|---|
 | 531 |       GraphIndices,
 | 
|---|
 | 532 |       atom_particlenames);
 | 
|---|
| [f20b4b] | 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;
 | 
|---|
| [c4aeda] | 537 |     const detail::GraphIndices_t::right_const_iterator graphiter =
 | 
|---|
| [f33ef9] | 538 |         GraphIndices.right.find(walkerid);
 | 
|---|
| [c4aeda] | 539 |     ASSERT( graphiter != GraphIndices.right.end(),
 | 
|---|
| [f33ef9] | 540 |         "PotentialFitPartialChargesAction::performCall() - cannot find "
 | 
|---|
 | 541 |         +toString(walkerid)+" in GraphIndices.");
 | 
|---|
| [c4aeda] | 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(),
 | 
|---|
| [f33ef9] | 546 |         "PotentialFitPartialChargesAction::performCall() - cannot find "
 | 
|---|
 | 547 |         +toString(graphindex)+" in atom_particlenames.");
 | 
|---|
| [c4aeda] | 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() - ");
 | 
|---|
| [f20b4b] | 552 |     walker->setParticleName(nameiter->second);
 | 
|---|
| [f33ef9] | 553 |     LOG(1, "INFO: atom " << *walker << " received the following particle "
 | 
|---|
| [f20b4b] | 554 |         << walker->getParticleName());
 | 
|---|
| [c4aeda] | 555 |   }
 | 
|---|
| [c4a323] | 556 | 
 | 
|---|
 | 557 |   return Action::success;
 | 
|---|
 | 558 | }
 | 
|---|
 | 559 | 
 | 
|---|
| [50d49d] | 560 | ActionState::ptr PotentialFitPartialChargesAction::performUndo(ActionState::ptr _state) {
 | 
|---|
| [c4a323] | 561 |   return Action::success;
 | 
|---|
 | 562 | }
 | 
|---|
 | 563 | 
 | 
|---|
| [50d49d] | 564 | ActionState::ptr PotentialFitPartialChargesAction::performRedo(ActionState::ptr _state){
 | 
|---|
| [c4a323] | 565 |   return Action::success;
 | 
|---|
 | 566 | }
 | 
|---|
 | 567 | 
 | 
|---|
| [50d49d] | 568 | bool PotentialFitPartialChargesAction::canUndo() {
 | 
|---|
| [c4a323] | 569 |   return false;
 | 
|---|
 | 570 | }
 | 
|---|
 | 571 | 
 | 
|---|
| [50d49d] | 572 | bool PotentialFitPartialChargesAction::shouldUndo() {
 | 
|---|
| [c4a323] | 573 |   return false;
 | 
|---|
 | 574 | }
 | 
|---|
 | 575 | /** =========== end of function ====================== */
 | 
|---|