Global2c2eMatrix.h

Go to the documentation of this file.

#ifndef SPARROW_GLOBAL2C2EMATRIX_H
#define SPARROW_GLOBAL2C2EMATRIX_H

#include "Global2c2eTerms.h"
#include "Local2c2eMatrix.h"
#include "multipoleTypes.h"
#include <Utils/Math/AutomaticDifferentiation/MethodsTypesHelper.h>
#include <Eigen/Core>

namespace Scine {
namespace Sparrow {

namespace nddo {

namespace multipole {
class ChargeSeparationParameter;
class KlopmanParameter;

class Global2c2eMatrix {
 public:
  using orb_index_t = int;
  using orbPair_index_t = int;

  Global2c2eMatrix(int l1, int l2, const ChargeSeparationParameter& D1, const ChargeSeparationParameter& D2,
                   const KlopmanParameter& r1, const KlopmanParameter& r2);

  void setSymmetric(bool sym) {
    localZero_.setSymmetric(sym);
    localOne_.setSymmetric(sym);
    sameElement_ = sym;
  }

  template<Utils::DerivativeOrder O>
  void calculate(const Eigen::Vector3d& Rab);

  template<Utils::Derivative O>
  Utils::AutomaticDifferentiation::DerivativeType<O> getDerivative(orb_index_t o1, orb_index_t o2, orb_index_t o3,
                                                                   orb_index_t o4) const;
  template<Utils::Derivative O>
  Utils::AutomaticDifferentiation::DerivativeType<O> getDerivative(orbPair_index_t op1, orbPair_index_t op2) const;
  double get(orb_index_t o1, orb_index_t o2, orb_index_t o3, orb_index_t o4) const;
  double get(orbPair_index_t op1, orbPair_index_t op2) const;
  orbPair_index_t getPairIndex(orb_index_t o1, orb_index_t o2) const;
  const Eigen::MatrixXd& getGlobalMatrix() const;
  void output() const;

 private:
  Eigen::RowVector3d getFirstDerivative(orbPair_index_t op1, orbPair_index_t op2) const;
  Utils::AutomaticDifferentiation::Second3D getSecondDerivative(orbPair_index_t op1, orbPair_index_t op2) const;

  template<Utils::DerivativeOrder O>
  void evaluate();

  template<Utils::DerivativeOrder O>
  Utils::AutomaticDifferentiation::Value3DType<O> evaluateMatrixElement(orbPair_index_t op1, orbPair_index_t op2);

  const int d1_, d2_;
  Global2c2eTerms terms_;
  const ChargeSeparationParameter &dist1, &dist2;
  const KlopmanParameter &rho1, &rho2;
  Local2c2eMatrix<Utils::DerivativeOrder::Zero> localZero_;
  Local2c2eMatrix<Utils::DerivativeOrder::One> localOne_;
  Local2c2eMatrix<Utils::DerivativeOrder::Two> localTwo_;
  OrbitalRotation<Utils::DerivativeOrder::Zero> rotationsZero_;
  OrbitalRotation<Utils::DerivativeOrder::One> rotationsOne_;
  OrbitalRotation<Utils::DerivativeOrder::Two> rotationsTwo_;
  double R_, RNormX_, RNormY_, RNormZ_;
  Utils::AutomaticDifferentiation::First3D nullDeriv, oneDeriv;
  Eigen::MatrixXd globalMatrix_;
  Eigen::Matrix<Utils::AutomaticDifferentiation::First3D, Eigen::Dynamic, Eigen::Dynamic> globalMatrixOne_;
  Eigen::Matrix<Utils::AutomaticDifferentiation::Second3D, Eigen::Dynamic, Eigen::Dynamic> globalMatrixTwo_;
  bool sameElement_;
  TwoElectronIntegralIndexes pairIndexes_;
};

inline Eigen::RowVector3d Global2c2eMatrix::getFirstDerivative(orbPair_index_t op1, orbPair_index_t op2) const {
  if (op1 == 100 || op2 == 100)
    return Utils::Gradient::Zero();

  return globalMatrixOne_(op1, op2).derivatives();
}

inline Utils::AutomaticDifferentiation::Second3D Global2c2eMatrix::getSecondDerivative(orbPair_index_t op1,
                                                                                       orbPair_index_t op2) const {
  if (op1 == 100 || op2 == 100)
    return {};

  return globalMatrixTwo_(op1, op2);
}

template<Utils::Derivative O>
Utils::AutomaticDifferentiation::DerivativeType<O> Global2c2eMatrix::getDerivative(orb_index_t o1, orb_index_t o2,
                                                                                   orb_index_t o3, orb_index_t o4) const {
  return getDerivative<O>(getPairIndex(o1, o2), getPairIndex(o3, o4));
}

template<>
inline Utils::AutomaticDifferentiation::DerivativeType<Utils::Derivative::First>
Global2c2eMatrix::getDerivative<Utils::Derivative::First>(orbPair_index_t op1, orbPair_index_t op2) const {
  return getFirstDerivative(op1, op2);
}

inline const Eigen::MatrixXd& Global2c2eMatrix::getGlobalMatrix() const {
  return globalMatrix_;
}

template<>
inline Utils::AutomaticDifferentiation::DerivativeType<Utils::Derivative::SecondAtomic>
Global2c2eMatrix::getDerivative<Utils::Derivative::SecondAtomic>(orbPair_index_t op1, orbPair_index_t op2) const {
  return getSecondDerivative(op1, op2);
}

template<>
inline Utils::AutomaticDifferentiation::DerivativeType<Utils::Derivative::SecondFull>
Global2c2eMatrix::getDerivative<Utils::Derivative::SecondFull>(orbPair_index_t op1, orbPair_index_t op2) const {
  return getSecondDerivative(op1, op2);
}

} // namespace multipole

} // namespace nddo

} // namespace Sparrow
} // namespace Scine
#endif // SPARROW_GLOBAL2C2EMATRIX_H