Changeset fa2c89

Timestamp:

Mar 24, 2017, 10:12:25 AM (8 years ago)

Author:

Frederik Heber <heber@…>

Branches:

FitPartialCharges_GlobalError

Children:

9afdef

Parents:

e6a2ef

git-author:

Frederik Heber <heber@…> (10/10/16 12:35:24)

git-committer:

Frederik Heber <heber@…> (03/24/17 10:12:25)

Message:

Added calculation of l2 error of fitted charges against each full_potential level.

File:

• src/Actions/PotentialAction/FitFragmentPartialChargesAction.cpp (modified) (2 diffs)

src/Actions/PotentialAction/FitFragmentPartialChargesAction.cpp

@@ -330 +330 @@
 }
 
+static bool getFullPositions(
+    PartialNucleiChargeFitter::positions_t &positions,
+    PartialNucleiChargeFitter::charges_t &charges)
+{
+  /// \note we cannot use the summed up Fragment here, as the saturation hydrogens
+  /// are in the way and cannot be sorted out properly/in a simple fashion.
+  const World &world = const_cast<const World &>(World::getInstance());
+  const ParticleRegistry &registry = const_cast<const ParticleRegistry &>(ParticleRegistry::getInstance());
+  {
+    const World::ConstAtomComposite &atoms = world.getAllAtoms();
+    positions.reserve(atoms.size());
+    charges.reserve(atoms.size());
+    for (World::ConstAtomComposite::const_iterator iter = atoms.begin();
+        iter != atoms.end(); ++iter) {
+      // set position for this atom
+      const Vector &pos = (*iter)->getPosition();
+
+      // use the fitted charges
+      const atomId_t atomid = (*iter)->getId();
+      std::string particlename = (*iter)->getParticleName();
+      if (particlename.empty())
+        particlename = (*iter)->getElement().getSymbol();
+      if (registry.isPresentByName(particlename)) {
+        const Particle * const particle = registry.getByName(particlename);
+        LOG(3, "DEBUG: Using implicit charge " << particle->charge << " of particle "
+            << particle->getName() << " for atom " << atomid);
+        positions.push_back((1./AtomicLengthToAngstroem)*pos);
+        charges.push_back(particle->charge);
+      } else {
+        ELOG(1, "Could not find fitted charge for " << particlename);
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+
 static bool isNotHydrogen(const atom * const _atom)
 {
@@ -588 +626 @@
   }
 
+  // compute the global error using the full solutions
+  {
+    FragmentationResultContainer::longrangedata_t longrangeData = resultscontainer.getLongRangeResults();
+    FragmentationResultContainer::longrangedata_t::iterator iter = longrangeData.end();
+    const size_t MaxLevel = std::max_element(fragments.begin(), fragments.end(), keyset_comparator)->size();
+    std::advance(iter, -MaxLevel);
+    FragmentationResultContainer::longrangedata_t::iterator remove_iter = iter;
+    std::vector<VMGData> fullsolutionData;
+    for (; iter != longrangeData.end(); ++iter)
+      fullsolutionData.push_back(iter->second);
+    if (longrangeData.size() > 1) // when there's just a single fragment, it corresponds to full solution
+      longrangeData.erase(remove_iter, longrangeData.end());
+
+    // extraction of all positions is not easy as we cannot use the summed up
+    // fragments due to saturated hydrogens. Also, we cannot easily removed the
+    // non-selected atoms from the set. Hence, the error is always calculated
+    // against all atoms.
+    PartialNucleiChargeFitter::positions_t fullpositions;
+    PartialNucleiChargeFitter::charges_t fullcharges;
+    getFullPositions(fullpositions, fullcharges);
+    PartialNucleiChargeFitter::charges_t zerocharges(fullpositions.size(), 0.);
+    LOG(4, "DEBUG: Full positions are " << fullpositions);
+    LOG(4, "DEBUG: Full (fitted) charges are " << fullcharges);
+
+    for (size_t i=0;i<fullsolutionData.size();++i) {
+      // get positions and potential
+      const SamplingGrid &fullpotential = fullsolutionData[i].both_sampled_potential;
+      if ((fullpotential.level == 0)
+          || ((fullpotential.begin[0] == fullpotential.end[0])
+              && (fullpotential.begin[1] == fullpotential.end[1])
+              && (fullpotential.begin[2] == fullpotential.end[2]))) {
+        ELOG(1, "Sampled grid contains grid made of zero points.");
+        return Action::failure;
+      }
+
+      // set up matrix and calculate residual
+      PartialNucleiChargeFitter fitter(fullpotential, fullpositions, params.radius.get());
+      if (!fitter.setCharges(zerocharges)) {
+        ELOG(1, "fullcharges and fullpositions do not coincide in size, internal error.");
+        return Action::failure;
+      }
+      VectorContent potential_norm = fitter.calculateResiduum();
+      if (!fitter.setCharges(fullcharges)) {
+        ELOG(1, "fullcharges and fullpositions do not coincide in size, internal error.");
+        return Action::failure;
+      }
+      VectorContent residuum = fitter.calculateResiduum();
+
+      LOG(1, "INFO: At order " << i+1 << " we have an error of " << residuum.Norm()
+          << " against the potential norm of " << potential_norm.Norm());
+    }
+  }
+
   return Action::success;
 #else