ewaldrealspacesum.cc

Go to the documentation of this file.

/*
 *            Copyright 2009-2026 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
 
// Standard includes
#include <algorithm>
#include <exception>
#include <fstream>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <unordered_map>
#include <vector>
 
// Local VOTCA includes
#include "votca/xtp/ewaldrealspacesum.h"
 
namespace votca {
namespace xtp {
 
namespace {
// BUG FIX (this session): PolarSegment::getPos() (via AtomContainer<T>::
// calcPos()) is a MASS-weighted center of mass. Legacy's own equivalent
// (PolarSeg::CalcPos()) is a plain, UNWEIGHTED arithmetic mean of site
// positions instead -- confirmed by direct comparison of the two
// implementations. Since this position feeds directly into the
// minimum-image PBC wrap below (raw_offset = target - source's own
// representative position), even a small difference between the two
// definitions can flip which periodic image is chosen as the true
// minimum for a given pair -- a discrete, not continuous, effect,
// consistent with the scattered (not uniformly scaled) per-site
// mismatch pattern that motivated this fix, rather than a smooth
// site-by-site drift.
//
// Deliberately NOT fixed by changing calcPos() itself: that's a shared
// AtomContainer<T> method used well beyond this file (md2qmengine,
// segmentmapper), where a genuine physical mass-weighted center of mass
// may be exactly what's wanted. This local helper instead gives THIS
// file its own, legacy-matching definition, without touching shared
// code any other caller depends on.
Eigen::Vector3d UnweightedCentroid(const PolarSegment& seg) {
  Eigen::Vector3d pos = Eigen::Vector3d::Zero();
  Index n = 0;
  for (const PolarSite& site : seg) {
    pos += site.getPos();
    ++n;
  }
  if (n > 0) {
    pos /= double(n);
  }
  return pos;
}
}  // namespace
 
EwaldRealSpaceSum::EwaldRealSpaceSum(
    const Eigen::Matrix3d& box, const EwaldRegistry& registry, double alpha,
    double thole_a, double r_min, double field_tol, double shell_width,
    Index n_max, double screening_factor,
    const std::vector<std::pair<Index, Eigen::Vector3d>>& foreground)
    : real_space_cutoff_(screening_factor / alpha),
      segment_radius_(0.0),
      box_(box),
      registry_(registry),
      interactor_(alpha, thole_a),
      r_min_(r_min),
      field_tol_(field_tol),
      shell_width_(shell_width),
      n_max_(n_max) {
  translations_ = GenerateSortedTranslations();
 
  for (const auto& entry : foreground) {
    foreground_[entry.first].push_back(entry.second);
  }
 
  // Largest site-to-centroid distance anywhere in the registry -- the
  // margin the segment-granular distance cull needs (see AddFieldAt).
  // Computed once here rather than per call: the registry's geometry is
  // fixed for this object's lifetime, the same assumption the neighbor
  // cache already relies on.
  for (Index id : registry_.AllIds()) {
    for (EwaldChargeState state :
         {EwaldChargeState::Neutral, EwaldChargeState::Electron,
          EwaldChargeState::Hole}) {
      if (!registry_.Has(id, state)) {
        continue;
      }
      const PolarSegment& segment = registry_.Get(id, state);
      const Eigen::Vector3d centroid = UnweightedCentroid(segment);
      for (const PolarSite& site : segment) {
        segment_radius_ =
            std::max(segment_radius_, (site.getPos() - centroid).norm());
      }
    }
  }
}
 
std::vector<EwaldRealSpaceSum::Translation>
    EwaldRealSpaceSum::GenerateSortedTranslations() const {
  std::vector<Translation> translations;
  const Eigen::Vector3d a = box_.col(0);
  const Eigen::Vector3d b = box_.col(1);
  const Eigen::Vector3d c = box_.col(2);
 
  for (Index na = -n_max_; na <= n_max_; ++na) {
    for (Index nb = -n_max_; nb <= n_max_; ++nb) {
      for (Index nc = -n_max_; nc <= n_max_; ++nc) {
        Eigen::Vector3d t = double(na) * a + double(nb) * b + double(nc) * c;
        translations.push_back({t, t.norm()});
      }
    }
  }
 
  std::sort(
      translations.begin(), translations.end(),
      [](const Translation& x, const Translation& y) { return x.r < y.r; });
  return translations;
}
 
template <enum Estatic CE>
void EwaldRealSpaceSum::ApplyNeighborList(const NeighborList& list,
                                          PolarSite& target,
                                          bool include_static) const {
  for (const auto& entry : list) {
    // Stored by pointer rather than by id -- see neighbor_cache_'s own
    // declaration for why. The Has()/Get() guard the id-based version
    // needed is gone with it: the cache is keyed on source_state, so
    // every entry in THIS list was registered at THIS state when the
    // list was built, and segment addresses are stable for the
    // registry's lifetime.
    const PolarSegment& source_segment = *std::get<0>(entry);
    const Index translation_idx = std::get<1>(entry);
    const Eigen::Vector3d& baseline_shift = std::get<2>(entry);
    const Eigen::Vector3d t = baseline_shift + translations_[translation_idx].t;
    for (const PolarSite& source_site : source_segment) {
      // Shift passed through rather than applied to a copy of the
      // source site -- see ApplyStaticField's own source_shift note.
      if (include_static) {
        interactor_.ApplyStaticField<PolarSite, CE>(source_site, target, t);
      }
      interactor_.ApplyInducedField<CE>(source_site, target, t);
    }
  }
}
 
void EwaldRealSpaceSum::InsertNeighborList(
    const std::pair<const PolarSite*, EwaldChargeState>& key, NeighborList list,
    const SearchCounts& counts) const {
  culled_entries_ += counts.culled;
  foreground_entries_ += counts.foreground;
  ++cached_targets_;
  cached_entries_ += Index(list.size());
  neighbor_cache_.emplace(key, std::move(list));
}
 
template <enum Estatic CE>
void EwaldRealSpaceSum::AddFieldAt(Index target_segment_id, PolarSite& target,
                                   EwaldChargeState source_state,
                                   bool include_static) const {
  const std::pair<const PolarSite*, EwaldChargeState> cache_key(&target,
                                                                source_state);
  auto cached = neighbor_cache_.find(cache_key);
  if (cached == neighbor_cache_.end()) {
    // First call for this target: run the convergence search on a
    // scratch copy, store the list, then fall through to the same path
    // every later call takes. Applying the list to target adds exactly
    // the contributions the search added to the copy, in the same order
    // and starting from the same V, so the field is bit-identical to
    // accumulating it during the search.
    SearchCounts counts;
    PolarSite scratch = target;
    NeighborList list = SearchNeighbors<CE>(
        target_segment_id, scratch, source_state, include_static, counts);
    InsertNeighborList(cache_key, std::move(list), counts);
    cached = neighbor_cache_.find(cache_key);
  }
  // Fast path: the real, geometry-determined neighbor set for this
  // target (see this class's own documentation for why reusing it across
  // calls for the same target is safe). Just apply every cached (source,
  // translation) pair directly -- no shell-by-shell search, no scan of
  // registry_.AllIds() members that turned out not to matter.
  ApplyNeighborList<CE>(cached->second, target, include_static);
}
 
template <enum Estatic CE>
EwaldRealSpaceSum::NeighborList EwaldRealSpaceSum::SearchNeighbors(
    Index target_segment_id, PolarSite& target, EwaldChargeState source_state,
    bool include_static, SearchCounts& counts) const {
  // The full convergence search. Records every (source, translation)
  // pair actually visited; the field it accumulates into target (a
  // scratch copy, see the declaration) is only there for the
  // shell-convergence test. Nothing in this object is written.
  NeighborList visited_pairs;
 
  // See the cull inside the shell loop below for what this is and why.
  const double cutoff_with_margin = real_space_cutoff_ + 2.0 * segment_radius_;
 
  Index shell_start = 0;
  double shell_edge = 0.0;
  bool converged = false;
 
  while (shell_start < Index(translations_.size())) {
    // Advance shell_edge by shell_width_ until it covers at least the next
    // not-yet-processed translation, then collect every translation up to
    // that edge into this shell. Since translations_ is sorted by
    // distance, this always yields a contiguous, radially-ordered shell.
    shell_edge =
        translations_[shell_start].r +
        (shell_edge > translations_[shell_start].r ? 0.0 : shell_width_);
    Index shell_end = shell_start;
    while (shell_end < Index(translations_.size()) &&
           translations_[shell_end].r <= shell_edge) {
      ++shell_end;
    }
 
    // Field is written into target.V()/V_noE() by ApplyStaticField/
    // ApplyInducedField as a side effect; recover this shell's own
    // contribution by differencing target's accumulator before and after
    // processing the whole shell.
    const Eigen::Vector3d before_shell =
        (CE == Estatic::noE_V) ? target.V_noE() : target.V();
 
    for (Index source_id : registry_.AllIds()) {
      // No per-segment skip: a target's own segment is treated like any
      // other. Its COINCIDENT copy is dropped by the foreground
      // suppression below; its other lattice images are ordinary
      // background molecules and stay. Legacy agrees -- SetupMidground
      // excludes by (segment id, na, nb, nc), not by segment id.
      if (!registry_.Has(source_id, source_state)) {
        continue;
      }
      const PolarSegment& source_segment =
          registry_.Get(source_id, source_state);
 
      // Per-source baseline wrap (computed once per source_id, not per
      // shell-translation index, since it depends only on the target
      // and source_segment's own representative positions, not on
      // which small shell-translation t is currently being applied).
      //
      // BUG FIX (this session, found via a real production comparison
      // against legacy): translations_ is sorted and shell-converged
      // by |t| alone (the raw lattice-translation magnitude), but the
      // quantity that actually determines a pair's true separation is
      // |raw_offset + t|, not |t| by itself. A translation with a
      // LARGE |t| can still be the one that gives the SMALLEST actual
      // distance, if it happens to cancel a large raw target-source
      // offset -- and such a translation is exactly what the
      // shell-convergence criterion (stop once shell radius >= r_min_)
      // can miss entirely, since it never even considers translations
      // whose OWN magnitude exceeds r_min_, regardless of what they'd
      // produce once combined with the raw offset. Confirmed on a real
      // pair: legacy's own minimum-image search (round-based,
      // independent of any shell-magnitude cutoff) found r=13.6 bohr
      // via t=(0,-1,-1)*box (|t|~97 bohr, comfortably past r_min_ in
      // this run's own ~75.6 bohr setting) for a pair this class's own
      // unfixed search reported no closer than r=55.4 bohr on, because
      // it never got there.
      //
      // Fixed by mirroring legacy's own two-step structure exactly:
      // wrap the raw (large, unbounded) target-source separation into
      // its true minimum image FIRST (fractional-coordinate rounding,
      // correct for any box_ shape, not just orthorhombic), THEN run
      // the existing small, bounded shell search of translations_ on
      // top of that already-small baseline -- exactly as legacy's own
      // dr12_pbc + L structure does. The shell search only ever needs
      // to explore the neighborhood right around the true minimum
      // image, not the raw, unbounded separation, so its own
      // shell_radius >= r_min_ stopping criterion is now sound (it was
      // never wrong in isolation -- it was being applied to the wrong
      // starting point).
      const Eigen::Vector3d source_centroid =
          UnweightedCentroid(source_segment);
      const Eigen::Vector3d raw_offset = target.getPos() - source_centroid;
      const Eigen::Vector3d frac = box_.inverse() * raw_offset;
      const Eigen::Vector3d wrapped_frac = frac - frac.array().round().matrix();
      const Eigen::Vector3d min_image_offset = box_ * wrapped_frac;
      const Eigen::Vector3d baseline_shift = raw_offset - min_image_offset;
 
      for (Index idx = shell_start; idx < shell_end; ++idx) {
        // Distance cull. The erfc(alpha*r) screening means a pair's
        // contribution falls off as fast as erfc does: at alpha*r = 6 a
        // single pair contributes ~1e-20 of the total field, and the
        // whole remaining tail sums to far below field_tol_. Without
        // this test every registered segment is evaluated at every
        // accepted translation, however far away it is -- measured on a
        // real 1000-segment system as 14985 (segment, translation)
        // entries per target, of which only ~22% lie inside the cutoff.
        // That factor of ~4.6, not any per-pair cost, was the whole
        // remaining real-space gap against legacy (whose own
        // PolarNbs list is distance-built, so it never had the excess).
        //
        // r_min_ does NOT already do this: it bounds the shell search by
        // TRANSLATION magnitude |t|, which is unrelated to how far a
        // given source segment is once that translation is applied.
        //
        // The true separation for this (segment, translation) pair is
        // |min_image_offset - translations_[idx].t| -- the same quantity
        // the interactor will form, since it evaluates
        // target - (source + t) with t = baseline_shift +
        // translations_[idx].t and baseline_shift = raw_offset -
        // min_image_offset. Compared at segment granularity, so
        // segment_radius_ (the largest site-to-centroid distance in the
        // registry) is added twice as a margin: no individual site pair
        // can then be closer than the cutoff while its segment pair is
        // culled.
        const double pair_distance =
            (min_image_offset - translations_[idx].t).norm();
        if (pair_distance > cutoff_with_margin) {
          ++counts.culled;
          continue;
        }
        const Eigen::Vector3d t = baseline_shift + translations_[idx].t;
 
        // Foreground suppression. This one periodic copy of this segment
        // is handled explicitly elsewhere (a polar region), so it must
        // not also appear in the periodic background -- see the
        // constructor's own foreground documentation. Only the copy
        // sitting at the recorded position is dropped; the segment's
        // other lattice images stay in the sum.
        if (!foreground_.empty()) {
          auto fg = foreground_.find(source_id);
          if (fg != foreground_.end()) {
            const Eigen::Vector3d shifted_centroid = source_centroid + t;
            bool suppressed = false;
            for (const Eigen::Vector3d& fg_pos : fg->second) {
              if ((shifted_centroid - fg_pos).norm() < kForegroundMatchTol) {
                suppressed = true;
                break;
              }
            }
            if (suppressed) {
              ++counts.foreground;
              continue;
            }
          }
        }
 
        // The target's own segment at ZERO translation is the segment
        // itself: intramolecular pairs plus the r = 0 self-pair, never
        // part of this sum (AddIntraSegmentCoupling owns the induced
        // side, the permanent side has its own compensation pass).
        //
        // Only where the segment has no recorded foreground copy. Where
        // it does, the suppression above already decides which copy is
        // carved out, and that need not be the t = 0 one: a foreground
        // segment sitting at a nonzero image has its t = 0 copy as a
        // genuine neighbour. Omitting this guard made the field infinite
        // in every background solve.
        const bool source_has_foreground =
            !foreground_.empty() && foreground_.count(source_id) > 0;
        if (!source_has_foreground && source_id == target_segment_id &&
            t.squaredNorm() < kSelfTranslationTol * kSelfTranslationTol) {
          continue;
        }
 
        visited_pairs.emplace_back(&source_segment, idx, baseline_shift);
        for (const PolarSite& source_site : source_segment) {
          if (include_static) {
            interactor_.ApplyStaticField<PolarSite, CE>(source_site, target, t);
          }
          interactor_.ApplyInducedField<CE>(source_site, target, t);
        }
      }
    }
 
    const Eigen::Vector3d after_shell =
        (CE == Estatic::noE_V) ? target.V_noE() : target.V();
    const Eigen::Vector3d shell_field = after_shell - before_shell;
 
    const double shell_radius =
        translations_[shell_end > shell_start ? shell_end - 1 : shell_start].r;
    if (shell_radius >= r_min_ && shell_field.norm() < field_tol_) {
      converged = true;
      break;
    }
 
    shell_start = shell_end;
  }
 
  if (!converged) {
    std::stringstream message;
    message << "EwaldRealSpaceSum: real-space sum for segment "
            << target_segment_id
            << " did not converge within the n_max lattice-translation "
               "search box; increase n_max or field_tol.";
    throw std::runtime_error(message.str());
  }
 
  return visited_pairs;
}
 
double EwaldRealSpaceSum::CalcStaticEnergyAt(
    const PolarSite& target, EwaldChargeState source_state) const {
  // See this method's own declaration for what is and is not included.
  const std::pair<const PolarSite*, EwaldChargeState> cache_key(&target,
                                                                source_state);
  auto cached = neighbor_cache_.find(cache_key);
  if (cached == neighbor_cache_.end()) {
    throw std::runtime_error(
        "EwaldRealSpaceSum::CalcStaticEnergyAt: no neighbour list for this "
        "target. AddFieldAt or PrepareNeighborCache must run first -- this "
        "method deliberately does not build one, so that an energy query "
        "cannot silently become the expensive shell search.");
  }
 
  double energy = 0.0;
  for (const auto& entry : cached->second) {
    const PolarSegment& source_segment = *std::get<0>(entry);
    const Index translation_idx = std::get<1>(entry);
    const Eigen::Vector3d& baseline_shift = std::get<2>(entry);
    const Eigen::Vector3d t = baseline_shift + translations_[translation_idx].t;
    for (const PolarSite& source_site : source_segment) {
      energy += interactor_.CalcStaticEnergy<PolarSite, PolarSite>(source_site,
                                                                   target, t);
    }
  }
  return energy;
}
 
Eigen::VectorXd EwaldRealSpaceSum::PotentialAtMany(
    Index target_segment_id, const std::vector<Eigen::Vector3d>& points,
    EwaldChargeState source_state) const {
  // Sources flattened once: AllIds() returns by value and the centroids
  // are fixed, so recomputing either per point would cost more than the
  // sum itself on a grid.
  // The sites are flattened too: for a unit probe only position, charge and
  // total (static + induced) dipole matter, and a compact array of those is
  // far cheaper to stream per grid point than the PolarSite objects.
  struct ProbeSource {
    Eigen::Vector3d pos;
    double q;
    Eigen::Vector3d mu;
  };
  std::vector<std::vector<ProbeSource>> segments;
  std::vector<Index> segment_ids;
  std::vector<Eigen::Vector3d> centroids;
  for (Index source_id : registry_.AllIds()) {
    if (!registry_.Has(source_id, source_state)) {
      continue;
    }
    const PolarSegment& segment = registry_.Get(source_id, source_state);
    std::vector<ProbeSource> sites;
    sites.reserve(segment.size());
    for (const PolarSite& site : segment) {
      sites.push_back({site.getPos(), site.getCharge(),
                       site.getStaticDipole() + site.getInducedDipole()});
    }
    segments.push_back(std::move(sites));
    segment_ids.push_back(source_id);
    centroids.push_back(UnweightedCentroid(segment));
  }
 
  const double cutoff_with_margin = real_space_cutoff_ + 2.0 * segment_radius_;
  const Index n_points = Index(points.size());
  const Index n_sources = Index(segments.size());
  Eigen::VectorXd phi = Eigen::VectorXd::Zero(n_points);
 
  // Hoisted out of both loops below: the box is fixed for this object's
  // lifetime, and inverting it per (point, source) is n_points * n_sources
  // 3x3 inversions for one matrix.
  const Eigen::Matrix3d box_inv = box_.inverse();
 
#pragma omp parallel for schedule(dynamic, 8)
  for (Index p = 0; p < n_points; ++p) {
    const Eigen::Vector3d& point = points[std::size_t(p)];
 
    double acc = 0.0;
    for (Index s_i = 0; s_i < n_sources; ++s_i) {
      const std::vector<ProbeSource>& source_segment =
          segments[std::size_t(s_i)];
      const Index source_id = segment_ids[std::size_t(s_i)];
      const Eigen::Vector3d& source_centroid = centroids[std::size_t(s_i)];
 
      // Minimum image first, shell search on top of it -- see AddFieldAt's
      // own account of why the raw separation cannot be handed to a
      // search bounded by |t|.
      const Eigen::Vector3d raw_offset = point - source_centroid;
      const Eigen::Vector3d frac = box_inv * raw_offset;
      const Eigen::Vector3d wrapped_frac = frac - frac.array().round().matrix();
      const Eigen::Vector3d min_image_offset = box_ * wrapped_frac;
      const Eigen::Vector3d baseline_shift = raw_offset - min_image_offset;
 
      const auto fg = foreground_.find(source_id);
      const bool source_has_foreground = fg != foreground_.end();
 
      // STOP, don't skip. translations_ is sorted by |t| ascending, and
      // |min_image_offset - t| >= |t| - |min_image_offset|, so once |t|
      // passes this limit every remaining translation is culled as well.
      // The bound is exact, not a heuristic: nothing inside the cutoff can
      // sit beyond it.
      //
      // This matters far more here than in the site-based paths above,
      // which never walk the whole list -- their shell-convergence search
      // hands them a [shell_start, shell_end) window. Without the break
      // this loop runs (2*n_max+1)^3 translations, 29791 at the default
      // n_max of 15, for EVERY (point, source) pair. On a DFT integration
      // grid against a 1000-segment registry that is of order 1e12
      // distance evaluations and hours of wall time; with it the surviving
      // range is a couple of lattice shells.
      const double t_limit = cutoff_with_margin + min_image_offset.norm();
 
      for (std::size_t idx = 0; idx < translations_.size(); ++idx) {
        if (translations_[idx].r > t_limit) {
          break;
        }
        const double pair_distance =
            (min_image_offset - translations_[idx].t).norm();
        if (pair_distance > cutoff_with_margin) {
          continue;
        }
        const Eigen::Vector3d t = baseline_shift + translations_[idx].t;
 
        // Foreground suppression: this one periodic copy is handled
        // explicitly elsewhere, so it must not also appear here.
        if (source_has_foreground) {
          const Eigen::Vector3d shifted_centroid = source_centroid + t;
          bool suppressed = false;
          for (const Eigen::Vector3d& fg_pos : fg->second) {
            if ((shifted_centroid - fg_pos).norm() < kForegroundMatchTol) {
              suppressed = true;
              break;
            }
          }
          if (suppressed) {
            continue;
          }
        }
 
        if (!source_has_foreground && source_id == target_segment_id &&
            t.squaredNorm() < kSelfTranslationTol * kSelfTranslationTol) {
          continue;
        }
 
        // Unit-charge probe at the point: what CalcStaticEnergy plus the
        // UNDAMPED CalcInducedSourceEnergy return for it. Undamped because
        // the target is a point in space, not a point-polarizable site, so
        // there is no overlap for Thole to correct -- and the rest of the
        // package already treats [induced dipole] x [QM density] undamped
        // (AOMultipole and DFTEngine::ExternalRepulsion apply no Thole).
        const Eigen::Vector3d probe_offset = point - t;
        for (const ProbeSource& source_site : source_segment) {
          acc += interactor_.ScreenedPotential(source_site.q, source_site.mu,
                                               probe_offset - source_site.pos);
        }
      }
    }
    phi[p] = acc;
  }
  return phi;
}
 
double EwaldRealSpaceSum::CalcInducedSourceEnergyAt(
    const PolarSite& target, EwaldChargeState source_state) const {
  // Deliberately a near-copy of CalcStaticEnergyAt above rather than a
  // shared template over the interactor call: the two differ only in
  // which moments they contract, but they are validated separately and
  // against different legacy channels (_pp and half of _pu), and a
  // shared body would let a change to one silently move the other.
  const std::pair<const PolarSite*, EwaldChargeState> cache_key(&target,
                                                                source_state);
  auto cached = neighbor_cache_.find(cache_key);
  if (cached == neighbor_cache_.end()) {
    throw std::runtime_error(
        "EwaldRealSpaceSum::CalcInducedSourceEnergyAt: no neighbour list for "
        "this target. AddFieldAt or PrepareNeighborCache must run first -- "
        "this method deliberately does not build one, so that an energy "
        "query cannot silently become the expensive shell search.");
  }
 
  double energy = 0.0;
  for (const auto& entry : cached->second) {
    const PolarSegment& source_segment = *std::get<0>(entry);
    const Index translation_idx = std::get<1>(entry);
    const Eigen::Vector3d& baseline_shift = std::get<2>(entry);
    const Eigen::Vector3d t = baseline_shift + translations_[translation_idx].t;
    for (const PolarSite& source_site : source_segment) {
      energy += interactor_.CalcInducedSourceEnergy(source_site, target, t);
    }
  }
  return energy;
}
 
void EwaldRealSpaceSum::PrepareNeighborCache(
    const std::vector<std::pair<Index, PolarSite*>>& targets,
    EwaldChargeState source_state, bool converge_static_field) const {
  // The list AddFieldAt<V>(..., converge_static_field) would build.
  PrepareNeighborCacheImpl<Estatic::V>(targets, source_state,
                                       converge_static_field);
}
 
template <enum Estatic CE>
void EwaldRealSpaceSum::PrepareNeighborCacheImpl(
    const std::vector<std::pair<Index, PolarSite*>>& targets,
    EwaldChargeState source_state, bool include_static) const {
  // See PrepareNeighborCache's own declaration for what this is for. The search
  // is the one AddFieldAt runs on a miss (SearchNeighbors), so the two
  // paths cannot drift apart.
  //
  // Serial selection: skip targets that already have a list, and
  // repeated entries, so each missing list is searched exactly once.
  std::vector<std::size_t> todo;
  {
    std::unordered_map<std::pair<const PolarSite*, EwaldChargeState>, bool,
                       PairHash>
        seen;
    for (std::size_t i = 0; i < targets.size(); ++i) {
      const std::pair<const PolarSite*, EwaldChargeState> key(targets[i].second,
                                                              source_state);
      if (neighbor_cache_.count(key) > 0 || !seen.emplace(key, true).second) {
        continue;
      }
      todo.push_back(i);
    }
  }
  if (todo.empty()) {
    return;
  }
 
  // Parallel search. Each iteration reads this object and writes only
  // its own slots below. An exception cannot leave an omp region, so
  // the first one is carried out and rethrown afterwards.
  std::vector<NeighborList> lists(todo.size());
  std::vector<SearchCounts> counts(todo.size());
  std::exception_ptr failure = nullptr;
  const Index n_todo = Index(todo.size());
#pragma omp parallel for schedule(dynamic, 16)
  for (Index k = 0; k < n_todo; ++k) {
    try {
      const auto& entry = targets[todo[std::size_t(k)]];
      PolarSite scratch = *entry.second;
      lists[std::size_t(k)] =
          SearchNeighbors<CE>(entry.first, scratch, source_state,
                              include_static, counts[std::size_t(k)]);
    } catch (...) {
#pragma omp critical(ewald_real_prepare_failure)
      {
        if (!failure) {
          failure = std::current_exception();
        }
      }
    }
  }
  if (failure) {
    std::rethrow_exception(failure);
  }
 
  // Serial insertion, in target order.
  for (std::size_t k = 0; k < todo.size(); ++k) {
    const std::pair<const PolarSite*, EwaldChargeState> key(
        targets[todo[k]].second, source_state);
    InsertNeighborList(key, std::move(lists[k]), counts[k]);
  }
}
 
template <enum Estatic CE>
void EwaldRealSpaceSum::AddFieldAtMany(
    const std::vector<std::pair<Index, PolarSite*>>& targets,
    EwaldChargeState source_state, bool include_static) const {
  // Lists first: afterwards every AddFieldAt below takes the read-only
  // fast path. The list is converged on the same field AddFieldAt's own
  // on-demand search would have used, so results do not depend on
  // whether a target was prepared here or found on demand.
  PrepareNeighborCacheImpl<CE>(targets, source_state, include_static);
  // A site may appear more than once in targets. Its occurrences are
  // grouped and applied by one thread, in input order, so no two threads
  // add to the same site and the sum matches the serial one bit for bit.
  std::unordered_map<const PolarSite*, std::size_t> group_of;
  std::vector<std::vector<std::size_t>> groups;
  group_of.reserve(targets.size());
  groups.reserve(targets.size());
  for (std::size_t i = 0; i < targets.size(); ++i) {
    const auto inserted = group_of.emplace(targets[i].second, groups.size());
    if (inserted.second) {
      groups.emplace_back();
    }
    groups[inserted.first->second].push_back(i);
  }
  const Index n_groups = Index(groups.size());
#pragma omp parallel for schedule(dynamic, 16)
  for (Index g = 0; g < n_groups; ++g) {
    for (const std::size_t i : groups[std::size_t(g)]) {
      PolarSite& target = *targets[i].second;
      const auto cached = neighbor_cache_.find({&target, source_state});
      ApplyNeighborList<CE>(cached->second, target, include_static);
    }
  }
}
 
template void EwaldRealSpaceSum::AddFieldAt<Estatic::V>(Index, PolarSite&,
                                                        EwaldChargeState,
                                                        bool) const;
template void EwaldRealSpaceSum::AddFieldAt<Estatic::noE_V>(Index, PolarSite&,
                                                            EwaldChargeState,
                                                            bool) const;
template void EwaldRealSpaceSum::AddFieldAtMany<Estatic::V>(
    const std::vector<std::pair<Index, PolarSite*>>&, EwaldChargeState,
    bool) const;
template void EwaldRealSpaceSum::AddFieldAtMany<Estatic::noE_V>(
    const std::vector<std::pair<Index, PolarSite*>>&, EwaldChargeState,
    bool) const;
 
}  // namespace xtp
}  // namespace votca