ewaldrealspaceinteractor.h

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.
 *
 */
 
#pragma once
#ifndef VOTCA_XTP_EWALDREALSPACEINTERACTOR_H
#define VOTCA_XTP_EWALDREALSPACEINTERACTOR_H
 
// Local VOTCA includes
#include "classicalsegment.h"
#include "eeinteractor.h"
#include "eigen.h"

 
namespace votca {
namespace xtp {
 
class EwaldRealSpaceInteractor {
 public:
  explicit EwaldRealSpaceInteractor(double alpha, double thole_a = 0.39)
      : alpha_(alpha), thole_(thole_a) {};
 
  // Applies the erfc-screened field generated by a permanent multipole
  // (charge + dipole, taken from site1's static moments) at site2, and
  // returns the corresponding screened interaction energy. Undamped by
  // Thole -- site1 is treated as a permanent source regardless of whether
  // it happens to be a PolarSite.
  //
  // CE selects which of site2's field accumulators (V() or V_noE()) is
  // updated, mirroring eeInteractor::ApplyStaticField.
  //
  // source_shift is a periodic-image translation added to site1's own
  // position, WITHOUT modifying site1. It exists purely to avoid
  // materializing a shifted copy of the source site in the caller's hot
  // loop: EwaldRealSpaceSum evaluates every (target site, source site,
  // translation) triple, and constructing a PolarSite copy per triple
  // dominated the real-space cost (a PolarSite carries a Vector9d of
  // multipoles, two 3x3 matrices, several vectors, and a std::string
  // element name -- so the copy is far from free, and there are ~1e8 of
  // them per solver iteration on a 1000-segment system). Legacy's own
  // equivalent takes the shift as an argument for the same reason (see
  // EwdInteractor::FU12_ERFC_At_By's own vec &s overload). Defaults to
  // zero, so existing two-argument calls are unaffected.
  template <class T, enum Estatic CE>
  double ApplyStaticField(
      const T& site1, PolarSite& site2,
      const Eigen::Vector3d& source_shift = Eigen::Vector3d::Zero()) const;
 
  // SUBTRACTS the erf(alpha*r)-screened (not erfc) field of a permanent
  // multipole at site2 for two sites of the same molecule, as legacy
  // EwdInteractor::FP12_ERF_At_By ("compensation term"). The reciprocal-space
  // sum excludes no pair, so it contains the erf(alpha*r)/r part of every
  // intramolecular static-static pair; subtracting it makes the net
  // intramolecular static-static interaction zero, as in legacy and in
  // molecular mechanics (bonded terms already account for it; see
  // EwaldPeriodicDipoleOperator). Same CE/accumulator convention as
  // ApplyStaticField.
  //
  // source_shift (see ApplyStaticField's own note) is last in the
  // parameter list, not next to the sites it belongs with: a defaulted
  // parameter inserted mid-signature silently reinterprets existing
  // positional arguments at call sites that still compile. Keep new
  // parameters at the end.
  template <class T, enum Estatic CE>
  void ApplyErfStaticFieldCorrection(
      const T& site1, PolarSite& site2,
      const Eigen::Vector3d& source_shift = Eigen::Vector3d::Zero()) const;
 
  // The induced-dipole counterpart of ApplyErfStaticFieldCorrection:
  // removes the erf-screened (reciprocal-space) field that site1's
  // INDUCED dipole contributes at site2. Subtracted, like its static
  // sibling.
  //
  // Deliberately NOT Thole-damped. Thole damping is a short-range
  // real-space construct; the reciprocal-space sum this term removes
  // has no damping of its own, so damping the removal would not cancel
  // what was actually added. Legacy agrees -- EwdInteractor::
  // FU12_ERF_At_By uses the bare C1/C2 functions with no l3/l5 factor
  // anywhere, unlike its erfc sibling FU12_ERFC_At_By.
  //
  // Built on the same EvaluateSource path as the static correction
  // rather than re-deriving the field expression, so it inherits that
  // path's already-validated sign convention instead of depending on a
  // fresh derivation of it.
  template <enum Estatic CE>
  void ApplyErfInducedFieldCorrection(
      const PolarSite& site1, PolarSite& site2,
      const Eigen::Vector3d& source_shift = Eigen::Vector3d::Zero()) const;
 
  // Applies the erfc-screened, Thole-damped field generated by site1's
  // induced dipole at site2, and returns the corresponding energy.
  // source_shift: see ApplyStaticField's own note above.
  template <enum Estatic CE>
  double ApplyInducedField(
      const PolarSite& site1, PolarSite& site2,
      const Eigen::Vector3d& source_shift = Eigen::Vector3d::Zero()) const;
 
  // erfc-screened interaction energy between two permanent multipole
  // sites (charge + dipole only), undamped by Thole.
  //
  // source_shift, as elsewhere, is last in the parameter list: a
  // defaulted parameter inserted mid-signature silently reinterprets
  // existing positional arguments at call sites that still compile.
  template <class S1, class S2>
  double CalcStaticEnergy(
      const S1& site1, const S2& site2,
      const Eigen::Vector3d& source_shift = Eigen::Vector3d::Zero()) const;
 
  // The erf-screened counterpart of CalcStaticEnergy: the permanent
  // multipole interaction energy that the RECIPROCAL sum contributes for
  // this pair. Subtracting it removes a foreground copy's energy from
  // the periodic total, exactly as ApplyErfStaticFieldCorrection does
  // for the field.
  //
  // Undamped, for the same reason its field sibling is: the
  // reciprocal-space contribution it removes carries no Thole damping.
  // Takes the analytic r -> 0 limit for coincident sites, which is the
  // common case here rather than an edge case -- every foreground target
  // sits inside a suppressed copy.
  template <class S1, class S2>
  double CalcErfStaticEnergy(
      const S1& site1, const S2& site2,
      const Eigen::Vector3d& source_shift = Eigen::Vector3d::Zero()) const;
 
  // erfc-screened, Thole-damped interaction energy between two induced
  // dipoles.
  double CalcInducedEnergy(const PolarSite& site1,
                           const PolarSite& site2) const;
 
  // erfc-screened, Thole-damped interaction energy between site1's
  // INDUCED dipole (as source) and site2's PERMANENT moments (as
  // target) -- the [background induced] x [foreground permanent] term.
  //
  // This is the half of legacy's _pu channel that nothing else in this
  // code path covers. PolarRegion's own E_polar_ext contracts the
  // FOREGROUND's induced dipoles against the delivered field, which
  // gives [fg induced] x [bg permanent + bg induced]; the permanent
  // energies here give [fg permanent] x [bg permanent]. The remaining
  // corner, [fg permanent] x [bg induced], has no other home, and its
  // absence was measurable: it accounts for the factor ~1.9 by which
  // the neutral job's induced energy fell short of legacy's.
  //
  // Damped, unlike the permanent-multipole energies, because the source
  // is an induced dipole and Thole damping is exactly the short-range
  // correction to induced-dipole interactions. l3 multiplies the r^-3
  // terms (the potential, and the field's mu term), l5 the r^-5 term,
  // matching ApplyInducedField's own combination so the two cannot
  // drift apart.
  //
  // NOTE that this makes the total mildly alpha-dependent, since the
  // reciprocal-space partner carries no damping (a long-range sum has
  // no short-range correction to apply). Legacy has the identical
  // structure -- FU12_ERFC_At_By damps, FU12_ERF_At_By does not -- so
  // this reproduces its convention rather than improving on it. The
  // decomposition IS exactly alpha-independent once damping is switched
  // off, which is how the unit test pins it down.
  //
  // damp = false is for a target that is a POINT rather than a site: a
  // field-evaluation point, as EwaldRealSpaceSum::PotentialAtMany uses.
  // Thole damping corrects the overlap of two POINT-POLARIZABLE sites,
  // and a coordinate in space is not one of those -- so there is nothing
  // to correct, and the undamped screened interaction is the whole
  // answer.
  //
  // This is not a free choice, it is what the rest of the package
  // already does for the same physical interaction. A QM region gets the
  // classical regions' induced dipoles through AOMultipole (which reads
  // getDipole(), permanent + induced) and through DFTEngine::
  // ExternalRepulsion -> eeInteractor::CalcStaticEnergy_site: neither
  // applies Thole. In eeInteractor, Thole lives only in
  // FillTholeInteraction, which the induction solve and the polar-polar
  // energies use and the QM path never touches. Damping here would give
  // a three-region job two different conventions for [induced dipole] x
  // [QM density] in one Hamiltonian, decided by which code path the
  // dipole arrived through.
  //
  // It cannot be expressed by giving the probe no polarizability.
  // ComputeThole would then form au3 = 0, read that as COMPLETE overlap,
  // and damp maximally at every distance (l3 = l5 = 0, c3 = -B1_erf
  // where undamped wants +B1). A PolarSite cannot be built without a
  // polarizability in any case -- its constructor assigns one from the
  // element name -- so the probe unavoidably carries a number that means
  // nothing, and this flag is how the caller says so.
  double CalcInducedSourceEnergy(
      const PolarSite& site1, const PolarSite& site2,
      const Eigen::Vector3d& source_shift = Eigen::Vector3d::Zero(),
      bool damp = true) const;
 
  // The erf-screened counterpart of CalcInducedSourceEnergy, used to
  // remove a coincident foreground copy from the reciprocal side, just
  // as CalcErfStaticEnergy does for the permanent channel.
  //
  // UNDAMPED, matching ApplyErfInducedFieldCorrection: the
  // reciprocal-space contribution this removes carries no Thole damping,
  // so damping the removal would not cancel what was actually added.
  double CalcErfInducedSourceEnergy(
      const PolarSite& site1, const PolarSite& site2,
      const Eigen::Vector3d& source_shift = Eigen::Vector3d::Zero()) const;
 
  // The screened-Coulomb derivative functions B0, B1, B2 for separation r:
  //   B0 = erfc(alpha*r) / r
  //   B_l = [ (2l-1)*B_{l-1} + (2*alpha)^(2l-1) / sqrt(pi) * exp(-alpha^2 r^2)
  //   ] / r^2
  // B0 supports charge-charge, B1 additionally supports charge-dipole,
  // B2 additionally supports dipole-dipole. Cross-checked against
  // EwdInteractor::UpdateAllBls (xtp/ewald/ewaldactor.h) for correctness
  // -- that file was deleted with the legacy code; the check is in the
  // git history, not in the tree.
  // Public (rather than private) so the free helper functions in
  // ewaldrealspaceinteractor.cc's anonymous namespace can take it by
  // value/reference as a plain parameter type.
  // Potential at a point, r_vec = point - source position, of a source with
  // charge q and dipole mu (static plus induced). For a unit-charge probe
  // ScreenedPotential equals CalcStaticEnergy + CalcInducedSourceEnergy with
  // damp=false, and ErfPotential equals CalcErfStaticEnergy +
  // CalcErfInducedSourceEnergy, coincidence limit included -- with one erfc
  // and one exp instead of three of each. Used on DFT integration grids.
  double ScreenedPotential(double q, const Eigen::Vector3d& mu,
                           const Eigen::Vector3d& r_vec) const;
  double ErfPotential(double q, const Eigen::Vector3d& mu,
                      const Eigen::Vector3d& r_vec) const;
 
  struct BFunctions {
    double B0;
    double B1;
    double B2;
  };
 
  // Thole damping factors l3, l5 (multiplying the r^-3 and r^-5 field
  // terms respectively) for a pair of induced-dipole sites, using the
  // exponential Thole model already implemented in
  // eeInteractor::FillTholeInteraction. Returns {1.0, 1.0} (no damping)
  // once the pair is far enough apart that the exponential term is
  // numerically negligible, matching FillTholeInteraction's own cutoff.
  struct TholeFactors {
    double l3;
    double l5;
  };
 
  BFunctions ComputeB(double r) const;
 
  // The erf(alpha*r)-screened complement to ComputeB's own erfc-screened
  // B-functions -- same B0/B1/B2 roles, but for erf(alpha*r)/r rather
  // than erfc(alpha*r)/r. Computed as the bare (undamped, alpha-
  // independent) Coulomb B-functions minus ComputeB's own result: since
  // erf(x)+erfc(x)=1, and the B-function recursion is linear in B0, this
  // decomposition is exact, not an approximation (confirmed both
  // algebraically and against a finite-difference derivative of
  // erf(alpha*r)/r directly, matching to the precision the finite
  // difference itself allows). B0_bare=1/r, B1_bare=1/r^3, B2_bare=3/r^5
  // are the standard undamped point-multipole tensor coefficients.
  BFunctions ComputeErfB(double r) const;
  TholeFactors ComputeThole(double r, const PolarSite& site1,
                            const PolarSite& site2) const;
 
 private:
  // Below this separation two sites are treated as coincident by the
  // erf-screened corrections, which take their analytic r -> 0 limit
  // there. Matches legacy's own 1e-2 threshold in EwdInteractor.
  static constexpr double kCoincidenceTol = 1e-2;
 
  double alpha_;
  double thole_;
};
 
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_EWALDREALSPACEINTERACTOR_H