ewaldperiodicdipoleoperator.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_EWALDPERIODICDIPOLEOPERATOR_H
#define VOTCA_XTP_EWALDPERIODICDIPOLEOPERATOR_H
 
// Standard includes
#include <string>
#include <vector>
 
// Local VOTCA includes
#include "eeinteractor.h"
#include "eigen.h"
#include "ewaldrealspacesum.h"
#include "ewaldreciprocalspacesum.h"
#include "ewaldregistry.h"
#include "ewaldshapecorrection.h"

 
namespace votca {
namespace xtp {
class EwaldPeriodicDipoleOperator;
}
}  // namespace votca
 
namespace Eigen {
namespace internal {
// EwaldPeriodicDipoleOperator's traits must be specialized before the
// class itself inherits from Eigen::EigenBase<EwaldPeriodicDipoleOperator>
// below (Eigen::EigenBase's own definition instantiates
// internal::traits<Derived> immediately) -- matching the exact ordering
// DipoleDipoleInteraction's own header (dipoledipoleinteraction.h) uses
// for the same reason.
template <>
struct traits<votca::xtp::EwaldPeriodicDipoleOperator>
    : public Eigen::internal::traits<Eigen::MatrixXd> {};
}  // namespace internal
}  // namespace Eigen
 
namespace votca {
namespace xtp {
 
class EwaldPeriodicDipoleOperator
    : public Eigen::EigenBase<EwaldPeriodicDipoleOperator> {
 public:
  // Required typedefs, constants, and method:
  using Scalar = double;
  using RealScalar = double;
  using StorageIndex = votca::Index;
  enum {
    ColsAtCompileTime = Eigen::Dynamic,
    MaxColsAtCompileTime = Eigen::Dynamic,
    IsRowMajor = false
  };
 
  // ids: the segments this operator solves for, in the order that
  //   defines the layout of any vector this operator is applied to (see
  //   class documentation for the exact, variable-width layout). Every
  //   id must be registered in registry at EwaldChargeState::Neutral.
  // alpha_ewald: Ewald splitting parameter, shared with real_sum/
  //   recip_sum's own construction, needed here to build a matching
  //   erfc-screened intramolecular interactor.
  // thole_a: the SAME Thole damping parameter real_sum's own
  //   construction uses (polarmethod.aDamp in ewdbgpol.xml terms). This
  //   used to be an unused placeholder here (0.39, hardcoded, since
  //   intra_interactor_ was originally only ever used for its own public
  //   ComputeB): AddIntraSegmentCoupling did not apply Thole damping to
  //   the intramolecular term at all, on the -- since corrected --
  //   assumption that legacy's own treatment was undamped there too (see
  //   AddIntraSegmentCoupling's own documentation for the fuller
  //   history). It now genuinely calls ComputeThole, so this parameter
  //   is genuinely used and must be the real value, not a placeholder.
  // shape: the SAME EwaldShapeCorrection instance EwaldBackground's own
  //   permanent-field generation uses. This class previously had NO
  //   shape-correction mechanism at all for the induced case -- a real,
  //   previously-unnoticed gap, found only once EwdInteractor's own
  //   FU12_ShapeField_At_By was traced directly (called unconditionally
  //   inside PolarBackground's own induction-iteration loop, mirroring
  //   FP12_ShapeField_At_By's own permanent-field treatment exactly, but
  //   using each site's own CURRENT induced dipole rather than its
  //   static one). EwaldShapeCorrection::TotalDipoleMoment already sums
  //   charge + static dipole + induced dipole together (see its own
  //   documentation), so no new shape-correction class or formula is
  //   needed here -- only a new call site, added inside RawMultiply
  //   itself (see there for why this is genuinely linear in v and does
  //   not disturb baseline_'s own established v=0 subtraction pattern).
  EwaldPeriodicDipoleOperator(EwaldRegistry& registry,
                              const EwaldRealSpaceSum& real_sum,
                              const EwaldReciprocalSpaceSum& recip_sum,
                              const EwaldShapeCorrection& shape,
                              std::vector<Index> ids, double alpha_ewald,
                              double thole_a);
 
  // Debug/experimental. Accumulated wall-clock spent inside RawMultiply,
  // split by phase, across every call for the lifetime of this object.
  // Added to attribute the per-iteration cost of the solve without a
  // sampling profiler: the phases below are each a single loop, so a
  // deterministic accumulator answers "where does the matvec time go"
  // exactly, where sampling an inlined hot function tends to report
  // little beyond "in AddFieldAt". Times are in seconds; n_calls is the
  // number of RawMultiply invocations they are summed over (one per
  // solver iteration, plus one for the baseline_ construction).
  struct RawMultiplyTimings {
    double setup = 0.0;       // dipole set + Reset over every site
    double real_space = 0.0;  // EwaldRealSpaceSum::AddFieldAt
    double reciprocal = 0.0;  // EwaldReciprocalSpaceSum::AddFieldAtMany
    double shape = 0.0;       // EwaldShapeCorrection::AddFieldAt
    double assemble = 0.0;    // P^-1*v - V - M*v result assembly
    double intra = 0.0;       // AddIntraSegmentCoupling
    Index n_calls = 0;
    double total() const {
      return setup + real_space + reciprocal + shape + assemble + intra;
    }
  };
  const RawMultiplyTimings& Timings() const { return timings_; }
 
  class InnerIterator {
   public:
    InnerIterator(const EwaldPeriodicDipoleOperator& xpr, const Index& id)
        : xpr_(xpr), id_(id) {};
 
    InnerIterator& operator++() {
      row_++;
      return *this;
    }
    operator bool() const { return row_ < xpr_.size_; }
    double value() const { return xpr_(row_, id_); }
    Index row() const { return row_; }
    Index col() const { return id_; }
    Index index() const { return row(); }
 
   private:
    const EwaldPeriodicDipoleOperator& xpr_;
    const Index id_;
    Index row_ = 0;
  };
 
  Index rows() const { return size_; }
  Index cols() const { return size_; }
  Index outerSize() const { return size_; }
 
  template <typename Vtype>
  Eigen::Product<EwaldPeriodicDipoleOperator, Vtype, Eigen::AliasFreeProduct>
      operator*(const Eigen::MatrixBase<Vtype>& x) const {
    return Eigen::Product<EwaldPeriodicDipoleOperator, Vtype,
                          Eigen::AliasFreeProduct>(*this, x.derived());
  }
 
  // See class documentation: same-site entries return the real getPInv()
  // block; cross-site entries return 0.0 unconditionally (never actually
  // read by DiagonalPreconditioner, and expensive to compute honestly
  // for a periodic pair).
  double operator()(Index i, Index j) const;
 
  // The actual matrix-vector product; see class documentation for what
  // this computes.
  Eigen::VectorXd multiply(const Eigen::VectorXd& v) const;
 
 private:
  // Computes getPInv()*v + (periodic field from every OTHER registered
  // segment, including those not in ids_) at every ids_ site. Not itself
  // linear in v (it carries the constant leak from non-ids_ sources) --
  // multiply() below fixes that by subtracting baseline_ = RawMultiply(0).
  Eigen::VectorXd RawMultiply(const Eigen::VectorXd& v) const;
 
  // Adds the intramolecular erfc-screened, Thole-damped contribution
  // (see class documentation for the fuller history) into result in
  // place, for every segment in ids_ with more than one site. This term
  // is linear in v and contributes 0 at v=0, so it is folded directly
  // into RawMultiply's own result rather than needing any separate
  // baseline_ handling.
  void AddIntraSegmentCoupling(const Eigen::VectorXd& v,
                               Eigen::VectorXd& result) const;
 
  // Global vector index i's (segment index within ids_, site index
  // within that segment) -- the reverse of the offset table below,
  // found by binary search (std::upper_bound) on offsets_. Used only by
  // operator()(i,j), which is not performance-critical (see class
  // documentation).
  std::pair<Index, Index> LocateSite(Index i) const;
 
  EwaldRegistry& registry_;
  const EwaldRealSpaceSum& real_sum_;
  const EwaldReciprocalSpaceSum& recip_sum_;
  const EwaldShapeCorrection& shape_;
  std::vector<Index> ids_;
  // Used for its own public ComputeB (erfc-screened B-functions) AND
  // ComputeThole (Thole damping factors), constructed with the REAL
  // thole_a passed to this class's own constructor -- unlike an earlier
  // version of this class (see AddIntraSegmentCoupling's own
  // documentation), ComputeThole genuinely is called from this class
  // now, so thole_a here must be the real value, not a placeholder.
  EwaldRealSpaceInteractor intra_interactor_;
  // offsets_[n] is the global vector index (in units of a single scalar,
  // not a 3-vector) where segment ids_[n]'s own block starts;
  // offsets_[n+1]-offsets_[n] = 3 * (that segment's own site count).
  // offsets_.size() == ids_.size()+1, with offsets_.back() == size_.
  std::vector<Index> offsets_;
  Index size_;
  Eigen::VectorXd baseline_;
  mutable RawMultiplyTimings timings_;
  // The position-independent 3x3 matrix M such that a site's own trial
  // dipole v_i produces a spurious reciprocal-space self-field -M*v_i at
  // its own position (see EwaldReciprocalSpaceSum::SelfFieldMatrix's own
  // documentation for the derivation). Legacy explicitly removes this
  // same leak for the induced-dipole case (EwdInteractor::
  // FU12_ERF_At_By, called with the SAME site passed as both arguments
  // -- confirmed by tracing PolarBackground's own SOR iteration
  // directly -- with its own explicit "Note the (-): This is a
  // compensation term" comment); this class did not, until this was
  // traced down as the explanation for a large, carbon-specific
  // discrepancy that survived fixing the (separate) intramolecular
  // static-static leak. Computed once here (position-independent, same
  // for every site, depends only on the k-vector set/alpha/volume) and
  // added back into RawMultiply's own result for every site, cancelling
  // that site's own -M*v_i leak already present via site.V() (recip_sum_
  // never excludes anything, self included -- see its own class
  // documentation).
  Eigen::Matrix3d self_field_matrix_;
};
 
}  // namespace xtp
}  // namespace votca
 
namespace Eigen {
namespace internal {
template <typename Vtype>
struct generic_product_impl<votca::xtp::EwaldPeriodicDipoleOperator, Vtype,
                            DenseShape, DenseShape, GemvProduct>
    : generic_product_impl_base<
          votca::xtp::EwaldPeriodicDipoleOperator, Vtype,
          generic_product_impl<votca::xtp::EwaldPeriodicDipoleOperator,
                               Vtype>> {
  typedef
      typename Product<votca::xtp::EwaldPeriodicDipoleOperator, Vtype>::Scalar
          Scalar;
 
  template <typename Dest>
  static void scaleAndAddTo(Dest& dst,
                            const votca::xtp::EwaldPeriodicDipoleOperator& op,
                            const Vtype& v, const Scalar& alpha) {
    assert(alpha == Scalar(1) && "scaling is not implemented");
    EIGEN_ONLY_USED_FOR_DEBUG(alpha);
    Eigen::VectorXd temp = op.multiply(v);
    dst = temp.cast<Scalar>();
  }
};
}  // namespace internal
}  // namespace Eigen
 
#endif  // VOTCA_XTP_EWALDPERIODICDIPOLEOPERATOR_H