ewaldrealspaceinteractor.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 <cmath>
 
// Local VOTCA includes
#include "votca/xtp/ewaldrealspaceinteractor.h"
 
namespace votca {
namespace xtp {
 
namespace {
constexpr double kRSqrtPi = 0.5641895835477563;  // 1/sqrt(pi)
}  // namespace
 
EwaldRealSpaceInteractor::BFunctions EwaldRealSpaceInteractor::ComputeB(
    double r) const {
  const double rr1 = 1.0 / r;
  const double rr2 = rr1 * rr1;
  const double a1 = alpha_;
  const double a3 = a1 * a1 * a1;
  const double expTerm = kRSqrtPi * std::exp(-a1 * a1 * r * r);
 
  BFunctions b;
  b.B0 = std::erfc(a1 * r) * rr1;
  b.B1 = rr2 * (b.B0 + 2.0 * a1 * expTerm);
  b.B2 = rr2 * (3.0 * b.B1 + 4.0 * a3 * expTerm);
  return b;
}
 
EwaldRealSpaceInteractor::BFunctions EwaldRealSpaceInteractor::ComputeErfB(
    double r) const {
  const BFunctions erfc_b = ComputeB(r);
  const double rr1 = 1.0 / r;
  const double rr3 = rr1 * rr1 * rr1;
  BFunctions b;
  b.B0 = rr1 - erfc_b.B0;
  b.B1 = rr3 - erfc_b.B1;
  b.B2 = 3.0 * rr3 * rr1 * rr1 - erfc_b.B2;
  return b;
}
 
double EwaldRealSpaceInteractor::ScreenedPotential(
    double q, const Eigen::Vector3d& mu, const Eigen::Vector3d& r_vec) const {
  const BFunctions b = ComputeB(r_vec.norm());
  return q * b.B0 + mu.dot(r_vec) * b.B1;
}
 
double EwaldRealSpaceInteractor::ErfPotential(
    double q, const Eigen::Vector3d& mu, const Eigen::Vector3d& r_vec) const {
  const double r = r_vec.norm();
  if (r < kCoincidenceTol) {
    // limit of q erf(alpha r)/r; the dipole's potential vanishes there
    return q * 2.0 * alpha_ / std::sqrt(votca::tools::conv::Pi);
  }
  const BFunctions b = ComputeErfB(r);
  return q * b.B0 + mu.dot(r_vec) * b.B1;
}
 
EwaldRealSpaceInteractor::TholeFactors EwaldRealSpaceInteractor::ComputeThole(
    double r, const PolarSite& site1, const PolarSite& site2) const {
  // Reuses the exponential Thole model already implemented in
  // eeInteractor::FillTholeInteraction, so that the induced-induced
  // damping used here is identical to the one PolarRegion's own
  // finite-region induction solve already uses.
  const double au3 = thole_ * r * r * r * site1.getSqrtInvEigenDamp() *
                     site2.getSqrtInvEigenDamp();
  TholeFactors t{1.0, 1.0};
  if (au3 < 40.0) {
    const double expUa = std::exp(-au3);
    t.l3 = 1.0 - expUa;
    t.l5 = 1.0 - (1.0 + au3) * expUa;
  }
  return t;
}
 
namespace {
 
// Potential and field generated by a screened point charge q and dipole mu
// (source), evaluated at a point separated from the source by r_vec
// (pointing from source to target). Field is what a target site
// experiences; the target's own energy contribution is computed by the
// caller from this and the target's own moments.
struct ScreenedPotentialField {
  double phi;
  Eigen::Vector3d field;
};
 
ScreenedPotentialField EvaluateSource(
    double q, const Eigen::Vector3d& mu, const Eigen::Vector3d& r_vec, double r,
    const EwaldRealSpaceInteractor::BFunctions& b) {
  ScreenedPotentialField result;
  const double mu_dot_r = mu.dot(r_vec);
 
  // phi(r) = q*B0 + (mu . r_vec)*B1
  result.phi = q * b.B0 + mu_dot_r * b.B1;
 
  // E(r) = -grad(phi) = q*r_vec*B1 + [(mu.r_vec)*r_vec*B2 - mu*B1]
  // Note: B2 (from ComputeB's recursion) already carries the factor of
  // (2l-1)=3 for l=2, so the r_vec*r_vec dyadic term is (mu.r_vec)*B2,
  // not 3*(mu.r_vec)*B2 -- an earlier version of this file double-counted
  // that factor of 3 and was caught by the dipole-dipole cross-check in
  // test_ewaldrealspaceinteractor.cc.
  result.field = q * b.B1 * r_vec + mu_dot_r * b.B2 * r_vec - b.B1 * mu;
  (void)r;
  return result;
}
 
}  // namespace
 
template <class T, enum Estatic CE>
double EwaldRealSpaceInteractor::ApplyStaticField(
    const T& site1, PolarSite& site2,
    const Eigen::Vector3d& source_shift) const {
  const Eigen::Vector3d r_vec =
      site2.getPos() - (site1.getPos() + source_shift);
  const double r = r_vec.norm();
  const BFunctions b = ComputeB(r);
 
  const double q1 = site1.getCharge();
  const Eigen::Vector3d mu1 = site1.getStaticDipole();
  const ScreenedPotentialField src = EvaluateSource(q1, mu1, r_vec, r, b);
 
  if (CE == Estatic::noE_V) {
    site2.V_noE() += src.field;
  } else {
    site2.V() += src.field;
  }
 
  const double q2 = site2.getCharge();
  const Eigen::Vector3d mu2 = site2.getStaticDipole();
  // U = q2*phi1(r2) - mu2.E1(r2)
  return q2 * src.phi - mu2.dot(src.field);
}
 
template <class T, enum Estatic CE>
void EwaldRealSpaceInteractor::ApplyErfStaticFieldCorrection(
    const T& site1, PolarSite& site2,
    const Eigen::Vector3d& source_shift) const {
  const Eigen::Vector3d r_vec =
      site2.getPos() - (site1.getPos() + source_shift);
  const double r = r_vec.norm();
 
  // Coincident sites. The erf-screened functions are singular at r = 0
  // in form only: erf(ar)/r tends to 2a/sqrt(pi), and the dipole field
  // tends to the analytic Ewald self-term below. This case is reached
  // whenever the correction is applied to a foreground copy that
  // contains the target itself, which happens for every target -- so it
  // is the common path, not an edge case. Legacy guards it identically
  // (EwdInteractor::FP12_ERF_At_By / FU12_ERF_At_By, R1 < 1e-2 branch),
  // with the same 4/3 * alpha^3 / sqrt(pi) coefficient that
  // EwaldReciprocalSpaceSum::SelfFieldMatrix returns.
  //
  // A charge contributes no field to itself (no direction), so only the
  // dipole term survives.
  if (r < kCoincidenceTol) {
    const double self_coeff = (4.0 / 3.0) * alpha_ * alpha_ * alpha_ /
                              std::sqrt(votca::tools::conv::Pi);
    // SIGN: this branch must be the r -> 0 limit of src.field below, so
    // that the "-=" means the same on both sides of kCoincidenceTol.
    // That limit is NEGATIVE -- EvaluateSource's dipole term is -B1*mu
    // and B1_erf -> +4/3 alpha^3/sqrt(pi) -- so the field tends to
    // -self_coeff*mu. The other sign makes the correction jump by
    // 2*self_coeff*mu across the threshold, which no continuous
    // function does.
    const Eigen::Vector3d self_field = -self_coeff * site1.getStaticDipole();
    if (CE == Estatic::noE_V) {
      site2.V_noE() -= self_field;
    } else {
      site2.V() -= self_field;
    }
    return;
  }
  const BFunctions b = ComputeErfB(r);
 
  const double q1 = site1.getCharge();
  const Eigen::Vector3d mu1 = site1.getStaticDipole();
  const ScreenedPotentialField src = EvaluateSource(q1, mu1, r_vec, r, b);
 
  // Subtracted, not added -- see this method's own header documentation
  // (and legacy EwdInteractor::FP12_ERF_At_By's own explicit "Note the
  // (-)" comment) for why: this removes the reciprocal-space leak for
  // this intramolecular pair, it does not add a genuine field
  // contribution of its own.
  if (CE == Estatic::noE_V) {
    site2.V_noE() -= src.field;
  } else {
    site2.V() -= src.field;
  }
}
 
template <enum Estatic CE>
void EwaldRealSpaceInteractor::ApplyErfInducedFieldCorrection(
    const PolarSite& site1, PolarSite& site2,
    const Eigen::Vector3d& source_shift) const {
  const Eigen::Vector3d r_vec =
      site2.getPos() - (site1.getPos() + source_shift);
  const double r = r_vec.norm();
 
  // Coincident sites. The erf-screened functions are singular at r = 0
  // in form only: erf(ar)/r tends to 2a/sqrt(pi), and the dipole field
  // tends to the analytic Ewald self-term below. This case is reached
  // whenever the correction is applied to a foreground copy that
  // contains the target itself, which happens for every target -- so it
  // is the common path, not an edge case. Legacy guards it identically
  // (EwdInteractor::FP12_ERF_At_By / FU12_ERF_At_By, R1 < 1e-2 branch),
  // with the same 4/3 * alpha^3 / sqrt(pi) coefficient that
  // EwaldReciprocalSpaceSum::SelfFieldMatrix returns.
  //
  // A charge contributes no field to itself (no direction), so only the
  // dipole term survives.
  if (r < kCoincidenceTol) {
    const double self_coeff = (4.0 / 3.0) * alpha_ * alpha_ * alpha_ /
                              std::sqrt(votca::tools::conv::Pi);
    // SIGN: negative, as in ApplyErfStaticFieldCorrection's branch.
    //
    // This is the one of the four coincidence branches that is live for
    // a rank-0 system -- the other three contract against a STATIC
    // dipole, zero when the .mps carry only charges -- and it was the
    // whole of an alpha^3 drift in the delivered induced field, measured
    // at alpha^2.9 over alpha = 1...4 1/nm on an 18-segment job and flat
    // once the background's induced dipoles were switched off.
    // EwaldPeriodicDipoleOperator::multiply has this sign right (it
    // forms V + S*mu), which is why the background solve was already
    // alpha-independent while this path was not.
    const Eigen::Vector3d self_field = -self_coeff * site1.getInducedDipole();
    if (CE == Estatic::noE_V) {
      site2.V_noE() -= self_field;
    } else {
      site2.V() -= self_field;
    }
    return;
  }
  const BFunctions b = ComputeErfB(r);
 
  // Charge deliberately zero: this removes the field of site1's INDUCED
  // dipole only. Its permanent multipoles are handled by
  // ApplyErfStaticFieldCorrection, and passing them here too would
  // remove them twice.
  //
  // No ComputeThole call -- see this method's own header documentation.
  const ScreenedPotentialField src =
      EvaluateSource(0.0, site1.getInducedDipole(), r_vec, r, b);
 
  if (CE == Estatic::noE_V) {
    site2.V_noE() -= src.field;
  } else {
    site2.V() -= src.field;
  }
}
 
template <enum Estatic CE>
double EwaldRealSpaceInteractor::ApplyInducedField(
    const PolarSite& site1, PolarSite& site2,
    const Eigen::Vector3d& source_shift) const {
  const Eigen::Vector3d r_vec =
      site2.getPos() - (site1.getPos() + source_shift);
  const double r = r_vec.norm();
  const BFunctions b = ComputeB(r);
  const BFunctions berf = ComputeErfB(r);
  const TholeFactors t = ComputeThole(r, site1, site2);
 
  const Eigen::Vector3d mu1 = site1.getInducedDipole();
  const double mu_dot_r = mu1.dot(r_vec);
 
  // Thole damping combines as l3*B1 + l5*B2 (legacy's rule). B2 already
  // carries its (2l-1)=3 recursion factor -- see EvaluateSource.
  //
  // WHAT IS DAMPED: the BARE interaction, not the erfc-screened one.
  // Thole damping is physics (overlap of two smeared densities, so the
  // real interaction is l*T_bare); the Ewald split is bookkeeping, with
  // alpha carrying no physical content. So l*T_bare is the object to be
  // split, and since the reciprocal sum contributes an UNDAMPED B_erf,
  // this half must supply
  //
  //     l*B_bare - B_erf  ==  l*B_erfc + (l-1)*B_erf
  //
  // l*B_erfc alone sums to l*T_bare + (1-l)*erf(ar)/r, which depends on
  // alpha -- not a split at all. Written in the second form so l = 1 is
  // visibly unchanged and no large bare terms cancel at small r.
  const double c3 = t.l3 * b.B1 + (t.l3 - 1.0) * berf.B1;
  const double c5 = t.l5 * b.B2 + (t.l5 - 1.0) * berf.B2;
  const Eigen::Vector3d field = mu_dot_r * c5 * r_vec - c3 * mu1;
 
  if (CE == Estatic::noE_V) {
    site2.V_noE() += field;
  } else {
    site2.V() += field;
  }
 
  const Eigen::Vector3d mu2 = site2.getInducedDipole();
  return -mu2.dot(field);
}
 
template <class S1, class S2>
double EwaldRealSpaceInteractor::CalcStaticEnergy(
    const S1& site1, const S2& site2,
    const Eigen::Vector3d& source_shift) const {
  const Eigen::Vector3d r_vec =
      site2.getPos() - (site1.getPos() + source_shift);
  const double r = r_vec.norm();
  const BFunctions b = ComputeB(r);
 
  const double q1 = site1.getCharge();
  const Eigen::Vector3d mu1 = site1.getStaticDipole();
  const ScreenedPotentialField src = EvaluateSource(q1, mu1, r_vec, r, b);
 
  const double q2 = site2.getCharge();
  const Eigen::Vector3d mu2 = site2.getStaticDipole();
  return q2 * src.phi - mu2.dot(src.field);
}
 
template <class S1, class S2>
double EwaldRealSpaceInteractor::CalcErfStaticEnergy(
    const S1& site1, const S2& site2,
    const Eigen::Vector3d& source_shift) const {
  const Eigen::Vector3d r_vec =
      site2.getPos() - (site1.getPos() + source_shift);
  const double r = r_vec.norm();
 
  const double q1 = site1.getCharge();
  const Eigen::Vector3d mu1 = site1.getStaticDipole();
  const double q2 = site2.getCharge();
  const Eigen::Vector3d mu2 = site2.getStaticDipole();
 
  // Coincident sites. As in the field corrections, the erf-screened
  // functions are singular in form only at r = 0: erf(ar)/r tends to
  // 2a/sqrt(pi), and the dipole-dipole term to the analytic Ewald
  // self-value. Reached for every foreground target, since each sits
  // inside a copy that is being removed.
  if (r < kCoincidenceTol) {
    const double sqrt_pi = std::sqrt(votca::tools::conv::Pi);
    const double phi_self = 2.0 * alpha_ / sqrt_pi;
    const double dip_self = (4.0 / 3.0) * alpha_ * alpha_ * alpha_ / sqrt_pi;
    // Charge-charge through the self-potential, dipole-dipole through
    // the self-field. The cross terms vanish: a charge produces no field
    // at its own position, and a dipole no potential.
    //
    // SIGN: the general branch returns q2*src.phi - mu2.src.field, and
    // src.field tends to -dip_self*mu1 (see
    // ApplyErfStaticFieldCorrection's coincidence branch), so the
    // dipole-dipole piece enters with a PLUS here. Inert while the .mps
    // files are rank 0, since mu1 and mu2 are then both zero.
    return q2 * q1 * phi_self + dip_self * mu2.dot(mu1);
  }
 
  const BFunctions b = ComputeErfB(r);
  const ScreenedPotentialField src = EvaluateSource(q1, mu1, r_vec, r, b);
  return q2 * src.phi - mu2.dot(src.field);
}
 
double EwaldRealSpaceInteractor::CalcErfInducedSourceEnergy(
    const PolarSite& site1, const PolarSite& site2,
    const Eigen::Vector3d& source_shift) const {
  // See this method's own declaration. Mirrors CalcErfStaticEnergy, with
  // the source's INDUCED dipole in place of its permanent moments and no
  // Thole damping.
  const Eigen::Vector3d r_vec =
      site2.getPos() - (site1.getPos() + source_shift);
  const double r = r_vec.norm();
 
  const Eigen::Vector3d mu1 = site1.getInducedDipole();
  const double q2 = site2.getCharge();
  const Eigen::Vector3d mu2 = site2.getStaticDipole();
 
  if (r < kCoincidenceTol) {
    const double sqrt_pi = std::sqrt(votca::tools::conv::Pi);
    const double dip_self = (4.0 / 3.0) * alpha_ * alpha_ * alpha_ / sqrt_pi;
    // A dipole produces no potential at its own position, so the
    // charge-dipole cross term drops; only dipole-dipole survives.
    //
    // SIGN: plus, matching the r -> 0 limit of the general branch below
    // (see CalcErfStaticEnergy's own coincidence branch). Inert at rank
    // 0, since mu2 is the TARGET's static dipole.
    return dip_self * mu2.dot(mu1);
  }
 
  const BFunctions b = ComputeErfB(r);
  const ScreenedPotentialField src = EvaluateSource(0.0, mu1, r_vec, r, b);
  return q2 * src.phi - mu2.dot(src.field);
}
 
double EwaldRealSpaceInteractor::CalcInducedEnergy(
    const PolarSite& site1, const PolarSite& site2) const {
  const Eigen::Vector3d r_vec = site2.getPos() - site1.getPos();
  const double r = r_vec.norm();
  const BFunctions b = ComputeB(r);
  const BFunctions berf = ComputeErfB(r);
  const TholeFactors t = ComputeThole(r, site1, site2);
 
  const Eigen::Vector3d mu1 = site1.getInducedDipole();
  const Eigen::Vector3d mu2 = site2.getInducedDipole();
  const double mu1_dot_r = mu1.dot(r_vec);
 
  // Same damped-bare-minus-erf combination as ApplyInducedField; see the
  // long note there for why l multiplies the bare interaction and not the
  // erfc-screened one. This method has no production caller at present
  // (only tests), and is converted with the others so that a future
  // caller does not inherit a convention the rest of the class has left
  // behind.
  const double c3 = t.l3 * b.B1 + (t.l3 - 1.0) * berf.B1;
  const double c5 = t.l5 * b.B2 + (t.l5 - 1.0) * berf.B2;
  const Eigen::Vector3d field1 = mu1_dot_r * c5 * r_vec - c3 * mu1;
  return -mu2.dot(field1);
}
 
double EwaldRealSpaceInteractor::CalcInducedSourceEnergy(
    const PolarSite& site1, const PolarSite& site2,
    const Eigen::Vector3d& source_shift, bool damp) const {
  // See this method's own declaration for what this term is, why it
  // needs a home of its own, why it is damped by default, and why a
  // field-evaluation point asks for damp = false rather than being
  // handed a zero polarizability.
  const Eigen::Vector3d r_vec =
      site2.getPos() - (site1.getPos() + source_shift);
  const double r = r_vec.norm();
  const BFunctions b = ComputeB(r);
  const BFunctions berf = ComputeErfB(r);
  const TholeFactors t =
      damp ? ComputeThole(r, site1, site2) : TholeFactors{1.0, 1.0};
 
  const Eigen::Vector3d mu1 = site1.getInducedDipole();
  const double mu_dot_r = mu1.dot(r_vec);
 
  // Potential and field of a screened, damped point dipole. The field
  // is character-for-character ApplyInducedField's own expression -- see
  // the long note there for why the damping multiplies the BARE
  // interaction and the undamped erf piece is subtracted back off. The
  // potential is its r^-3 partner, phi = (mu . r_vec) * B1, carrying the
  // same c3 that the field's mu term does. Writing them together here
  // keeps the two from drifting apart.
  const double c3 = t.l3 * b.B1 + (t.l3 - 1.0) * berf.B1;
  const double c5 = t.l5 * b.B2 + (t.l5 - 1.0) * berf.B2;
  const double phi = c3 * mu_dot_r;
  const Eigen::Vector3d field = mu_dot_r * c5 * r_vec - c3 * mu1;
 
  // Target's PERMANENT moments only: its induced dipole is PolarRegion's
  // business, via the field this code separately delivers.
  const double q2 = site2.getCharge();
  const Eigen::Vector3d mu2 = site2.getStaticDipole();
  return q2 * phi - mu2.dot(field);
}
 
// Explicit instantiations for the source types actually used.
template double
    EwaldRealSpaceInteractor::ApplyStaticField<StaticSite, Estatic::V>(
        const StaticSite&, PolarSite&, const Eigen::Vector3d&) const;
template double
    EwaldRealSpaceInteractor::ApplyStaticField<StaticSite, Estatic::noE_V>(
        const StaticSite&, PolarSite&, const Eigen::Vector3d&) const;
template double
    EwaldRealSpaceInteractor::ApplyStaticField<PolarSite, Estatic::V>(
        const PolarSite&, PolarSite&, const Eigen::Vector3d&) const;
template double
    EwaldRealSpaceInteractor::ApplyStaticField<PolarSite, Estatic::noE_V>(
        const PolarSite&, PolarSite&, const Eigen::Vector3d&) const;
 
template void EwaldRealSpaceInteractor::ApplyErfStaticFieldCorrection<
    StaticSite, Estatic::V>(const StaticSite&, PolarSite&,
                            const Eigen::Vector3d&) const;
template void EwaldRealSpaceInteractor::ApplyErfStaticFieldCorrection<
    StaticSite, Estatic::noE_V>(const StaticSite&, PolarSite&,
                                const Eigen::Vector3d&) const;
template void EwaldRealSpaceInteractor::ApplyErfStaticFieldCorrection<
    PolarSite, Estatic::V>(const PolarSite&, PolarSite&,
                           const Eigen::Vector3d&) const;
template void EwaldRealSpaceInteractor::ApplyErfStaticFieldCorrection<
    PolarSite, Estatic::noE_V>(const PolarSite&, PolarSite&,
                               const Eigen::Vector3d&) const;
 
template void EwaldRealSpaceInteractor::ApplyErfInducedFieldCorrection<
    Estatic::V>(const PolarSite&, PolarSite&, const Eigen::Vector3d&) const;
template void EwaldRealSpaceInteractor::ApplyErfInducedFieldCorrection<
    Estatic::noE_V>(const PolarSite&, PolarSite&, const Eigen::Vector3d&) const;
 
template double EwaldRealSpaceInteractor::ApplyInducedField<Estatic::V>(
    const PolarSite&, PolarSite&, const Eigen::Vector3d&) const;
template double EwaldRealSpaceInteractor::ApplyInducedField<Estatic::noE_V>(
    const PolarSite&, PolarSite&, const Eigen::Vector3d&) const;
 
template double
    EwaldRealSpaceInteractor::CalcStaticEnergy<StaticSite, StaticSite>(
        const StaticSite&, const StaticSite&, const Eigen::Vector3d&) const;
template double
    EwaldRealSpaceInteractor::CalcStaticEnergy<StaticSite, PolarSite>(
        const StaticSite&, const PolarSite&, const Eigen::Vector3d&) const;
template double
    EwaldRealSpaceInteractor::CalcStaticEnergy<PolarSite, StaticSite>(
        const PolarSite&, const StaticSite&, const Eigen::Vector3d&) const;
template double
    EwaldRealSpaceInteractor::CalcStaticEnergy<PolarSite, PolarSite>(
        const PolarSite&, const PolarSite&, const Eigen::Vector3d&) const;
 
template double
    EwaldRealSpaceInteractor::CalcErfStaticEnergy<StaticSite, StaticSite>(
        const StaticSite&, const StaticSite&, const Eigen::Vector3d&) const;
template double
    EwaldRealSpaceInteractor::CalcErfStaticEnergy<StaticSite, PolarSite>(
        const StaticSite&, const PolarSite&, const Eigen::Vector3d&) const;
template double
    EwaldRealSpaceInteractor::CalcErfStaticEnergy<PolarSite, StaticSite>(
        const PolarSite&, const StaticSite&, const Eigen::Vector3d&) const;
template double
    EwaldRealSpaceInteractor::CalcErfStaticEnergy<PolarSite, PolarSite>(
        const PolarSite&, const PolarSite&, const Eigen::Vector3d&) const;
 
}  // namespace xtp
}  // namespace votca