TransitionChargesCalculator.h

Go to the documentation of this file.

#ifndef SPARROW_TRANSITIONCHARGESCALCULATOR_H
#define SPARROW_TRANSITIONCHARGESCALCULATOR_H

#include <Sparrow/Implementations/TimeDependent/TimeDependentUtils.h>
#include <Utils/DataStructures/AtomsOrbitalsIndexes.h>
#include <Utils/DataStructures/SpinAdaptedMatrix.h>
#include <Utils/Typenames.h>
#include <Eigen/Core>
#include <map>
#include <vector>

namespace Scine {
namespace Utils {
class MolecularOrbitals;
class AtomsOrbitalsIndexes;
namespace LcaoUtils {
class ElectronicOccupation;
} // namespace LcaoUtils
} // namespace Utils
namespace Sparrow {

enum class AngularMomentum { None, S, P, D };
using AtomAndOrbitalShellIndex = std::pair<int, AngularMomentum>;
template<Utils::Reference restrictedness>
using TransitionChargesContainer =
    Utils::SpinAdaptedContainer<restrictedness, std::map<AtomAndOrbitalShellIndex, Eigen::VectorXd>>;

class InvalidOccupationException : public std::exception {
  const char* what() const noexcept final {
    return "Occupation does not respect the Aufbau principle.";
  }
};

class SpinPolarizedTransitionChargesNotImplementedException : public std::exception {
  const char* what() const noexcept final {
    return "Spin polarized version non yet implemented.";
  }
};

class SpinPolarizedOrbitalsNotAvailableException : public std::exception {
  const char* what() const noexcept final {
    return "No unrestricted orbitals present in unrestricted calculation.";
  }
};
class TransitionChargesCalculator {
 public:
  TransitionChargesCalculator(const Utils::MolecularOrbitals& molecularOrbitals, const Eigen::MatrixXd& overlapMatrix,
                              const Utils::AtomsOrbitalsIndexes& aoIndex);

  void fillOverlapProductMatrix();
  std::vector<Eigen::MatrixXd> calculateMORestrictedAtomicChargeMatrices() const;
  std::vector<Eigen::MatrixXd> calculateMOUnrestrictedAtomicChargeMatrices() const;
  template<Utils::Reference restrictedness>
  Eigen::MatrixXd calculateAtomicTransitionChargeMatrices(const Utils::LcaoUtils::ElectronicOccupation& /*occupation*/) const {
    throw std::runtime_error("Wrong specialization in calculateAtomicTransitionChargeMatrices()");
  }
  Eigen::MatrixXd calculateRestrictedTransitionChargeMatrices(const Utils::LcaoUtils::ElectronicOccupation& occupation) const;
  Eigen::MatrixXd calculateUnrestrictedTransitionChargeMatrices(const Utils::LcaoUtils::ElectronicOccupation& occupation) const;

 private:
  const Utils::MolecularOrbitals& mos_;
  const Eigen::MatrixXd& overlapMatrix_;
  const Utils::AtomsOrbitalsIndexes& aoIndex_;
  Utils::SpinAdaptedMatrix overlapProductMatrix_;
};

template<>
inline Eigen::MatrixXd TransitionChargesCalculator::calculateAtomicTransitionChargeMatrices<Utils::Reference::Restricted>(
    const Utils::LcaoUtils::ElectronicOccupation& occupation) const {
  if (!occupation.isFilledUpFromTheBottom()) {
    throw InvalidOccupationException();
  }

  auto nOccupied = occupation.numberOccupiedRestrictedOrbitals();
  auto nVirtual = mos_.numberOrbitals() - nOccupied;
  Eigen::MatrixXd transitionAtomicCharges = Eigen::MatrixXd::Zero(nOccupied * nVirtual, aoIndex_.getNAtoms());

  // Get correct coefficient matrix
  const Eigen::MatrixXd& occupiedBlock = mos_.restrictedMatrix().leftCols(nOccupied);
  const Eigen::MatrixXd& virtualBlock = mos_.restrictedMatrix().rightCols(nVirtual);
  const Eigen::MatrixXd& SCvir = overlapProductMatrix_.restrictedMatrix().rightCols(nVirtual);
  const Eigen::MatrixXd& SCocc = overlapProductMatrix_.restrictedMatrix().leftCols(nOccupied);
  for (int atom = 0; atom < aoIndex_.getNAtoms(); ++atom) {
    auto firstIndex = aoIndex_.getFirstOrbitalIndex(atom);
    auto nAOsOnAtom = aoIndex_.getNOrbitals(atom);
    Eigen::MatrixXd frontTransition =
        SCvir.middleRows(firstIndex, nAOsOnAtom).transpose() * occupiedBlock.middleRows(firstIndex, nAOsOnAtom);
    Eigen::MatrixXd backTransition =
        virtualBlock.middleRows(firstIndex, nAOsOnAtom).transpose() * SCocc.middleRows(firstIndex, nAOsOnAtom);

    transitionAtomicCharges.col(atom) = 0.5 * (Eigen::Map<Eigen::VectorXd>(frontTransition.data(), frontTransition.size()) +
                                               Eigen::Map<Eigen::VectorXd>(backTransition.data(), backTransition.size()));
  }
  return transitionAtomicCharges;
}

template<>
inline Eigen::MatrixXd TransitionChargesCalculator::calculateAtomicTransitionChargeMatrices<Utils::Reference::Unrestricted>(
    const Utils::LcaoUtils::ElectronicOccupation& occupation) const {
  if (!occupation.isFilledUpFromTheBottom()) {
    throw InvalidOccupationException();
  }
  if (mos_.alphaMatrix().cols() == 0 || overlapProductMatrix_.alphaMatrix().cols() == 0) {
    throw SpinPolarizedOrbitalsNotAvailableException();
  }

  // Assumes Aufbau principle construction
  auto nOccupiedAlpha = occupation.numberAlphaElectrons();
  auto nVirtualAlpha = mos_.alphaMatrix().cols() - nOccupiedAlpha;
  auto nOccupiedBeta = occupation.numberBetaElectrons();
  auto nVirtualBeta = mos_.betaMatrix().cols() - nOccupiedBeta;
  Eigen::MatrixXd transitionChargeMatrix(nOccupiedAlpha * nVirtualAlpha + nOccupiedBeta * nVirtualBeta, aoIndex_.getNAtoms());
  const Eigen::MatrixXd& occupiedBlockAlpha = mos_.alphaMatrix().leftCols(nOccupiedAlpha);
  const Eigen::MatrixXd& virtualBlockAlpha = mos_.alphaMatrix().rightCols(nVirtualAlpha);
  const Eigen::MatrixXd& occupiedBlockBeta = mos_.betaMatrix().leftCols(nOccupiedBeta);
  const Eigen::MatrixXd& virtualBlockBeta = mos_.betaMatrix().rightCols(nVirtualBeta);
  const Eigen::MatrixXd& SCAlphavir = overlapProductMatrix_.alphaMatrix().rightCols(nVirtualAlpha);
  const Eigen::MatrixXd& SCAlphaocc = overlapProductMatrix_.alphaMatrix().leftCols(nOccupiedAlpha);
  const Eigen::MatrixXd& SCBetavir = overlapProductMatrix_.betaMatrix().rightCols(nVirtualBeta);
  const Eigen::MatrixXd& SCBetaocc = overlapProductMatrix_.betaMatrix().leftCols(nOccupiedBeta);

  for (int atom = 0; atom < aoIndex_.getNAtoms(); ++atom) {
    int firstIndex = aoIndex_.getFirstOrbitalIndex(atom);
    int nAOsOnAtom = aoIndex_.getNOrbitals(atom);
    Eigen::MatrixXd transitionChargesAlpha =
        0.5 *
        (SCAlphavir.middleRows(firstIndex, nAOsOnAtom).transpose() * occupiedBlockAlpha.middleRows(firstIndex, nAOsOnAtom) +
         virtualBlockAlpha.middleRows(firstIndex, nAOsOnAtom).transpose() * SCAlphaocc.middleRows(firstIndex, nAOsOnAtom));
    Eigen::MatrixXd transitionChargesBeta =
        0.5 *
        (SCBetavir.middleRows(firstIndex, nAOsOnAtom).transpose() * occupiedBlockBeta.middleRows(firstIndex, nAOsOnAtom) +
         virtualBlockBeta.middleRows(firstIndex, nAOsOnAtom).transpose() * SCBetaocc.middleRows(firstIndex, nAOsOnAtom));
    transitionChargeMatrix.col(atom)
        << Eigen::Map<Eigen::VectorXd>(transitionChargesAlpha.data(), transitionChargesAlpha.size()),
        Eigen::Map<Eigen::VectorXd>(transitionChargesBeta.data(), transitionChargesBeta.size());
  }
  return transitionChargeMatrix;
}

} // namespace Sparrow
} // namespace Scine

#endif // SPARROW_TRANSITIONCHARGESCALCULATOR_H