environmentscreening.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_ENVIRONMENTSCREENING_H
#define VOTCA_XTP_ENVIRONMENTSCREENING_H
 
// Standard includes
#include <string>
#include <vector>
 
// Local VOTCA includes
#include "aobasis.h"
#include "classicalsegment.h"
#include "eigen.h"
#include "qmmolecule.h"

 
namespace votca {
namespace xtp {

struct ScreeningEnvironment {
  std::vector<PolarSegment> explicit_segments;
  double exp_damp = 0.39;
  std::vector<PolarSegment> shell_segments;
  double shell_dielectric = 4.0;
  bool include_kreac = true;
  // Width of the polar sites as seen by the QM charge density, in units of
  // alpha^(1/3): each site responds to the QM field averaged over a
  // unit-charge Gaussian exp(-r^2/R^2)/(pi^(3/2) R^3), R = site_width *
  // alpha_iso^(1/3). 0 gives point sites. See SiteWidths.
  double site_width = 0.5;
 
  bool empty() const {
    return explicit_segments.empty() && shell_segments.empty();
  }
};

struct ScreeningCheck {
  // Lowest eigenvalue of R = T^T B T. Must be > -1.
  double lowest = 0.0;
  // Number of eigenvalues at or below -1.
  Index n_unstable = 0;
  bool ok() const { return lowest > -1.0; }
  // Multi-line report: spectrum, the geometry of the closest contacts, and
  // for the lowest modes where their charge sits on the QM side and which
  // polar sites carry their reaction. Always filled.
  std::string report;
};
 
class EnvironmentScreening {
 public:
  static Eigen::MatrixXd AuxFieldAtPoints(
      const AOBasis& auxbasis, const std::vector<Eigen::Vector3d>& points,
      const std::vector<double>& widths = {});
  static Eigen::MatrixXd AuxPotentialAtPoints(
      const AOBasis& auxbasis, const std::vector<Eigen::Vector3d>& points,
      const std::vector<double>& widths);
  // widths (bohr, one per point, or empty for all zero): a point with
  // width R > 0 is a unit-charge Gaussian exp(-r^2/R^2)/(pi^(3/2) R^3)
  // rather than a point, and F is minus the gradient of the Coulomb
  // interaction (chi_Q | g_R) with respect to its centre -- a libint2
  // two-centre Coulomb integral, same stencil. Outside the extent of the
  // auxiliary functions a spherical charge acts as a point (Newton), so
  // this changes F only where a site sits inside their tails -- which is
  // exactly where point sites over-respond (see Check). R = 0 is the
  // nuclear-attraction path above, unchanged.

  static std::vector<double> SiteWidths(
      const std::vector<PolarSegment>& segments, double site_width);

  static Eigen::MatrixXd ReactionFieldKernel(
      const Eigen::MatrixXd& F, const std::vector<PolarSegment>& segments,
      double exp_damp);
  static Eigen::MatrixXd ChargeResponseKernel(
      const Eigen::MatrixXd& F, const Eigen::MatrixXd& Phi,
      const std::vector<PolarSegment>& segments, double exp_damp);

  static Eigen::MatrixXd ShellKernel(const Eigen::MatrixXd& F,
                                     const std::vector<PolarSegment>& segments,
                                     double epsilon);

  static Eigen::MatrixXd SymmetrizedReactionField(const Eigen::MatrixXd& B,
                                                  const Eigen::MatrixXd& T);

  static Eigen::MatrixXd DressingMatrix(const Eigen::MatrixXd& R);

  static Eigen::MatrixXd Kernel(const AOBasis& auxbasis,
                                const ScreeningEnvironment& env);

  static Eigen::MatrixXd Metric(const AOBasis& auxbasis);

  static ScreeningCheck Check(const AOBasis& auxbasis, const QMMolecule& atoms,
                              const ScreeningEnvironment& env,
                              const Eigen::MatrixXd& T, Index n_modes = 3);
 
 private:
  // Net charge of every auxiliary function, integral chi_Q d^3r, in the
  // order and normalization of the three-centre integrals. Defined, like
  // AuxFieldAtPoints, in libint2_calls.cc, the translation unit that
  // includes libint2.
  static Eigen::VectorXd AuxCharges(const AOBasis& auxbasis);
};
 
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_ENVIRONMENTSCREENING_H