aoshell.h

Go to the documentation of this file.

/*
 *            Copyright 2009-2021 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_AOSHELL_H
#define VOTCA_XTP_AOSHELL_H
 
// Third party includes
#include <boost/math/constants/constants.hpp>
 
// VOTCA includes
#include <votca/tools/constants.h>
 
// Local VOTCA includes
#include "basisset.h"
#include "eigen.h"
#include "qmatom.h"
// include libint last otherwise it overrides eigen
#include <libint2/shell.h>
 
namespace votca {
namespace xtp {
 
class AOBasis;
class AOShell;
class SetupCptTable;
 
class AOGaussianPrimitive {
 
 public:
  AOGaussianPrimitive(const GaussianPrimitive& gaussian);
  friend class AOShell;
 
  struct data {
    Index atomid;
    Index l;
    Index startindex;
    double decay;
    double contraction;
    double x;
    double y;
    double z;
    double scale;
  };
 
  AOGaussianPrimitive(const AOGaussianPrimitive::data& d) {
    decay_ = d.decay;
    contraction_ = d.contraction;
    powfactor_ = CalcPowFactor(decay_);
  }
 
  static void SetupCptTable(CptTable& table);
 
  void WriteData(data& d, const AOShell& s) const;
 
  double getPowfactor() const { return powfactor_; }
  double getDecay() const { return decay_; }
  double getContraction() const { return contraction_; }
 
 private:
  static double CalcPowFactor(double decay) {
    return std::pow(2.0 * decay / boost::math::constants::pi<double>(), 0.75);
  }
  double decay_;
  double contraction_;
  double powfactor_;  // used in evalspace to speed up DFT
};
 
/*
 * shells in a Gaussian-basis expansion
 */
class AOShell {
  friend AOBasis;
 
 public:
  AOShell(const Shell& shell, const QMAtom& atom, Index startIndex);
 
  AOShell(const AOGaussianPrimitive::data& d) {
    l_ = static_cast<L>(d.l);
    startIndex_ = d.startindex;
    atomindex_ = d.atomid;
    pos_ = Eigen::Vector3d(d.x, d.y, d.z);
    gaussians_.push_back(AOGaussianPrimitive(d));
  }
 
  L getL() const { return l_; }
  Index getNumFunc() const { return NumFuncShell(l_); };
  Index getCartesianNumFunc() const { return NumFuncShell_cartesian(l_); };
  Index getStartIndex() const { return startIndex_; }
  Index getOffset() const { return OffsetFuncShell(l_); }
  Index getCartesianOffset() const { return OffsetFuncShell_cartesian(l_); }
  Index getAtomIndex() const { return atomindex_; }
  Index getSize() const { return gaussians_.size(); }
 
  libint2::Shell LibintShell() const;
 
  const Eigen::Vector3d& getPos() const { return pos_; }
 
  void CalcMinDecay() {
    mindecay_ = std::numeric_limits<double>::max();
    for (auto& gaussian : gaussians_) {
      mindecay_ = std::min(mindecay_, gaussian.getDecay());
    }
  }
 
  double getMinDecay() const { return mindecay_; }
 
  struct AOValues {
 
    AOValues(Index size) {
      values = Eigen::VectorXd::Zero(size);
      derivatives = Eigen::MatrixX3d::Zero(size, 3);
    }
    Eigen::VectorXd values;
    Eigen::MatrixX3d derivatives;
  };
 
  AOValues EvalAOspace(const Eigen::Vector3d& grid_pos) const;

  using AODerivativeBlock = Eigen::Ref<Eigen::Matrix<double, Eigen::Dynamic, 3>,
                                       0, Eigen::OuterStride<>>;
  void EvalAOspace(const Eigen::Vector3d& grid_pos,
                   Eigen::Ref<Eigen::VectorXd> values,
                   AODerivativeBlock derivatives) const;
 
  // ===========================================================================
  // STATUS: written but NOT yet compiled/tested (see aoshell.cc for the
  // detailed derivation). Added for GGA XC-gradient support -- the
  // existing AOValues/EvalAOspace (first derivatives only) is untouched;
  // this is a new, separate struct/function reusing the SAME per-shell
  // polynomial structure, extended with second r-derivatives (Hessians),
  // needed to differentiate sigma=|grad(rho)|^2 for GGA functionals.
  //
  // Supports S, P, D, F, G shells (l=0..4), matching EvalAOspace's own
  // supported range exactly -- not a new limitation, same shells this
  // codebase's AO evaluator has always supported for values/gradients.
  //
  // General formula (verified symbolically/exactly against direct
  // double differentiation for representative polynomials spanning
  // every supported shell, via sympy, before writing any of this):
  // for chi(r) = exp(-alpha*d^2) * P(center), d^2=|center|^2,
  //   H_ij = prefactor * d2P/dxi_dxj
  //          - 2*alpha*delta_ij*AOvalue
  //          - 2*alpha*xi*grad_j - 2*alpha*xj*grad_i
  //          - 4*alpha^2*xi*xj*AOvalue
  // where prefactor, AOvalue, grad_i are exactly the quantities
  // EvalAOspace already computes for the SAME shell -- only the new
  // d2P/dxi_dxj terms (the polynomial P's own second derivatives,
  // exact per shell, also sympy-verified) are new per-shell content.
  // ===========================================================================
  struct AOValuesHessian {
    AOValuesHessian(Index size) {
      values = Eigen::VectorXd::Zero(size);
      derivatives = Eigen::MatrixX3d::Zero(size, 3);
      hessians = std::vector<Eigen::Matrix3d>(size, Eigen::Matrix3d::Zero());
    }
    Eigen::VectorXd values;
    Eigen::MatrixX3d derivatives;
    std::vector<Eigen::Matrix3d> hessians;
  };
 
  AOValuesHessian EvalAOspaceHessian(const Eigen::Vector3d& grid_pos) const;
 
  // iterator over pairs (decay constant; contraction coefficient)
  using GaussianIterator = std::vector<AOGaussianPrimitive>::const_iterator;
  GaussianIterator begin() const { return gaussians_.begin(); }
  GaussianIterator end() const { return gaussians_.end(); }
 
  // adds a Gaussian
  void addGaussian(const GaussianPrimitive& gaussian) {
    gaussians_.push_back(AOGaussianPrimitive(gaussian));
    return;
  }
 
  void normalizeContraction();
 
  friend std::ostream& operator<<(std::ostream& out, const AOShell& shell);
 
 private:
  L l_;
  // scaling factor
  // number of functions in shell
  double mindecay_;
  Index startIndex_;
  Eigen::Vector3d pos_;
  Index atomindex_;
 
  // vector of pairs of decay constants and contraction coefficients
  std::vector<AOGaussianPrimitive> gaussians_;
};
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_AOSHELL_H