TimeDependentUtils.h

Go to the documentation of this file.

#ifndef SPARROW_TIMEDEPENDENTUTILS_H
#define SPARROW_TIMEDEPENDENTUTILS_H

#include <Utils/DataStructures/MolecularOrbitals.h>
#include <Utils/DataStructures/SingleParticleEnergies.h>
#include <Utils/Scf/LcaoUtils/ElectronicOccupation.h>
#include <Utils/TimeDependent/TransitionDipoleCalculator.h>
#include <Eigen/Core>
#include <cereal/cereal.hpp>
#include <map>

namespace Scine {

namespace Sparrow {

struct ExcitedStatesParam {
  double c1 = 1, c2 = 1, d = 1;

  template<class Archive>
  void serialize(Archive& archive) {
    archive(CEREAL_NVP(c1), CEREAL_NVP(c2), CEREAL_NVP(d));
  }
};

class TimeDependentUtils {
 public:
  static Eigen::VectorXd generateEnergyDifferenceVector(int nOccupiedOrbitals, int nVirtualOrbitals,
                                                        const Utils::SingleParticleEnergies& energies);
  template<Utils::Reference restrictedness>
  static std::vector<std::string>
  generateExcitationsLabels(const Utils::SpinAdaptedContainer<restrictedness, std::vector<Utils::Excitation>>& excitations);

  static std::vector<std::string> generateExcitationsLabels(const std::vector<Utils::Excitation>& excitations,
                                                            const Eigen::Matrix<bool, -1, 1>& isBeta = {});

  template<Utils::Reference restrictedness>
  static void generateEnergyDifferenceVector(const Utils::SingleParticleEnergies& energies,
                                             const Utils::LcaoUtils::ElectronicOccupation& occupation,
                                             Utils::SpinAdaptedContainer<restrictedness, Eigen::VectorXd>& energyDifferenceVector);

  static void generateVirtualOrbitalIndices(const std::vector<int>& occupiedOrbitalsIndices,
                                            std::vector<int>& virtualOrbitalsIndices);
  static std::vector<Utils::Excitation> generateExcitations(const std::vector<int>& occupiedOrbitalsIndices,
                                                            const std::vector<int>& virtualOrbitalsIndices);
  template<Utils::Reference restrictedness>
  static Utils::SpinAdaptedContainer<restrictedness, std::vector<Utils::Excitation>>
  generateExcitations(const Utils::MolecularOrbitals& mos, const Utils::LcaoUtils::ElectronicOccupation& occupation);

  template<Utils::Reference restrictedness>
  static std::vector<int>
  generateEnergyOrderMap(const Utils::SpinAdaptedContainer<restrictedness, Eigen::VectorXd>& energyDifferenceVector);

  enum class Direction { To, From };
  template<typename Derived, typename DerivedResult>
  static void transformOrder(const Eigen::MatrixBase<Derived>& toTransform, const Eigen::MatrixBase<DerivedResult>& result,
                             const std::vector<int>& orderVector, Direction direction);
  template<typename Type>
  static void transformOrder(const std::vector<Type>& toTransform, std::vector<Type>& result,
                             const std::vector<int>& orderVector, Direction direction);

  static auto flatten(const Utils::SpinAdaptedContainer<Utils::Reference::Restricted, Eigen::VectorXd>& in) -> Eigen::VectorXd;
  static auto flatten(const Utils::SpinAdaptedContainer<Utils::Reference::Unrestricted, Eigen::VectorXd>& in) -> Eigen::VectorXd;

  template<typename T>
  static auto flatten(Utils::SpinAdaptedContainer<Utils::Reference::Restricted, std::vector<T>> in) -> std::vector<T>;

  template<typename T>
  static auto flatten(Utils::SpinAdaptedContainer<Utils::Reference::Unrestricted, std::vector<T>> in) -> std::vector<T>;
};

inline auto TimeDependentUtils::flatten(const Utils::SpinAdaptedContainer<Utils::Reference::Restricted, Eigen::VectorXd>& in)
    -> Eigen::VectorXd {
  return in.restricted;
}
inline auto TimeDependentUtils::flatten(const Utils::SpinAdaptedContainer<Utils::Reference::Unrestricted, Eigen::VectorXd>& in)
    -> Eigen::VectorXd {
  Eigen::VectorXd out(in.alpha.size() + in.beta.size());
  out << in.alpha, in.beta;
  return out;
}

template<typename T>
inline auto TimeDependentUtils::flatten(Utils::SpinAdaptedContainer<Utils::Reference::Restricted, std::vector<T>> in)
    -> std::vector<T> {
  return in.restricted;
}
template<typename T>
inline auto TimeDependentUtils::flatten(Utils::SpinAdaptedContainer<Utils::Reference::Unrestricted, std::vector<T>> in)
    -> std::vector<T> {
  std::vector<T> out;
  out.reserve(in.alpha.size() + in.beta.size());
  std::copy(in.alpha.begin(), in.alpha.end(), std::back_inserter(out));
  std::copy(in.beta.begin(), in.beta.end(), std::back_inserter(out));
  return out;
}

inline Eigen::VectorXd TimeDependentUtils::generateEnergyDifferenceVector(int nOccupiedOrbitals, int nVirtualOrbitals,
                                                                          const Utils::SingleParticleEnergies& energies) {
  auto const& energyLevels = energies.getRestrictedEnergies();

  Eigen::VectorXd energyDifferenceVector(nOccupiedOrbitals * nVirtualOrbitals);
  for (int occ = 0; occ < nOccupiedOrbitals; ++occ) {
    for (int vir = 0; vir < nVirtualOrbitals; ++vir) {
      energyDifferenceVector(nVirtualOrbitals * occ + vir) = energyLevels[nOccupiedOrbitals + vir] - energyLevels[occ];
    }
  }
  return energyDifferenceVector;
}

template<>
inline void TimeDependentUtils::generateEnergyDifferenceVector<Utils::Reference::Unrestricted>(
    const Utils::SingleParticleEnergies& energies, const Utils::LcaoUtils::ElectronicOccupation& occupation,
    Utils::SpinAdaptedContainer<Utils::Reference::Unrestricted, Eigen::VectorXd>& energyDifferenceVector) {
  const auto& alphaEnergies = energies.getAlphaEnergies();
  const auto& betaEnergies = energies.getBetaEnergies();
  auto const nAlphaMOs = alphaEnergies.size();
  auto const nBetaMOs = betaEnergies.size();
  const std::vector<int>& occupiedAlphaOrbitals = occupation.getFilledAlphaOrbitals();
  const std::vector<int>& occupiedBetaOrbitals = occupation.getFilledBetaOrbitals();
  auto const nOccAlpha = occupiedAlphaOrbitals.size();
  auto const nOccBeta = occupiedBetaOrbitals.size();
  auto const nVirAlpha = nAlphaMOs - nOccAlpha;
  auto const nVirBeta = nBetaMOs - nOccBeta;

  energyDifferenceVector.alpha.resize(nOccAlpha * nVirAlpha);
  energyDifferenceVector.beta.resize(nOccBeta * nVirBeta);

  int alphaIter = 0;
  for (auto const& occ : occupiedAlphaOrbitals) {
    for (unsigned int mo = 0; mo < nAlphaMOs; ++mo) {
      if ((std::find(occupiedAlphaOrbitals.begin(), occupiedAlphaOrbitals.end(), mo) == occupiedAlphaOrbitals.end())) {
        energyDifferenceVector.alpha[alphaIter] = alphaEnergies[mo] - alphaEnergies[occ];
        alphaIter++;
      }
    }
  }
  int betaIter = 0;
  for (auto const& occ : occupiedBetaOrbitals) {
    for (unsigned int mo = 0; mo < nBetaMOs; ++mo) {
      if ((std::find(occupiedBetaOrbitals.begin(), occupiedBetaOrbitals.end(), mo) == occupiedBetaOrbitals.end())) {
        energyDifferenceVector.beta[betaIter] = betaEnergies[mo] - betaEnergies[occ];
        betaIter++;
      }
    }
  }
}

template<>
inline void TimeDependentUtils::generateEnergyDifferenceVector<Utils::Reference::Restricted>(
    const Utils::SingleParticleEnergies& energies, const Utils::LcaoUtils::ElectronicOccupation& occupation,
    Utils::SpinAdaptedContainer<Utils::Reference::Restricted, Eigen::VectorXd>& energyDifferenceVector) {
  std::vector<int> occupiedOrbitals = occupation.getFilledRestrictedOrbitals();
  const auto& energiesLevels = energies.getRestrictedEnergies();
  auto const nMOs = energiesLevels.size();
  auto const nOcc = occupiedOrbitals.size();
  auto const nVir = nMOs - nOcc;

  energyDifferenceVector.restricted.resize(nOcc * nVir);

  int iter = 0;
  for (int occ : occupiedOrbitals) {
    for (unsigned int mo = 0; mo < nMOs; mo++) {
      if ((std::find(occupiedOrbitals.begin(), occupiedOrbitals.end(), mo) == occupiedOrbitals.end())) {
        energyDifferenceVector.restricted[iter] = energiesLevels[mo] - energiesLevels[occ];
        iter++;
      }
    }
  }
}

inline void TimeDependentUtils::generateVirtualOrbitalIndices(const std::vector<int>& occupiedOrbitalsIndices,
                                                              std::vector<int>& virtualOrbitalsIndices) {
  unsigned int nMOs = occupiedOrbitalsIndices.size() + virtualOrbitalsIndices.size();
  unsigned int iterEmpty = 0;
  unsigned int iterFilled = 0;
  for (int i = 0; i < static_cast<int>(nMOs); ++i) {
    if (iterFilled < occupiedOrbitalsIndices.size() && i == occupiedOrbitalsIndices[iterFilled]) {
      ++iterFilled;
    }
    else {
      virtualOrbitalsIndices[iterEmpty] = i;
      ++iterEmpty;
    }
  }
}

inline std::vector<Utils::Excitation> TimeDependentUtils::generateExcitations(const std::vector<int>& occupiedOrbitalsIndices,
                                                                              const std::vector<int>& virtualOrbitalsIndices) {
  std::vector<Utils::Excitation> excitations(occupiedOrbitalsIndices.size() * virtualOrbitalsIndices.size());

  int iter = 0;
  for (int occ : occupiedOrbitalsIndices) {
    for (int vir : virtualOrbitalsIndices) {
      excitations[iter].occ = occ;
      excitations[iter].vir = vir;
      iter++;
    }
  }
  return excitations;
}
template<>
inline Utils::SpinAdaptedContainer<Utils::Reference::Restricted, std::vector<Utils::Excitation>>
TimeDependentUtils::generateExcitations<Utils::Reference::Restricted>(const Utils::MolecularOrbitals& mos,
                                                                      const Utils::LcaoUtils::ElectronicOccupation& occupation) {
  int nMOs = mos.restrictedMatrix().cols();
  const auto& occupiedOrbitals = occupation.getFilledRestrictedOrbitals();
  std::vector<int> virtualOrbitals(nMOs - occupiedOrbitals.size());

  generateVirtualOrbitalIndices(occupiedOrbitals, virtualOrbitals);
  Utils::SpinAdaptedContainer<Utils::Reference::Restricted, std::vector<Utils::Excitation>> excitations;
  excitations.restricted.resize(occupiedOrbitals.size() * virtualOrbitals.size());

  int iter = 0;
  for (int occ : occupiedOrbitals) {
    for (int vir : virtualOrbitals) {
      excitations.restricted[iter].occ = occ;
      excitations.restricted[iter].vir = vir;
      iter++;
    }
  }
  return excitations;
}

template<>
inline Utils::SpinAdaptedContainer<Utils::Reference::Unrestricted, std::vector<Utils::Excitation>>
TimeDependentUtils::generateExcitations<Utils::Reference::Unrestricted>(const Utils::MolecularOrbitals& mos,
                                                                        const Utils::LcaoUtils::ElectronicOccupation& occupation) {
  int nMOs = mos.alphaMatrix().cols();
  const auto& occupiedAlphaOrbitals = occupation.getFilledAlphaOrbitals();
  const auto& occupiedBetaOrbitals = occupation.getFilledBetaOrbitals();
  std::vector<int> virtualAlphaOrbitals(nMOs - occupiedAlphaOrbitals.size());
  std::vector<int> virtualBetaOrbitals(nMOs - occupiedBetaOrbitals.size());

  generateVirtualOrbitalIndices(occupiedAlphaOrbitals, virtualAlphaOrbitals);
  generateVirtualOrbitalIndices(occupiedBetaOrbitals, virtualBetaOrbitals);
  Utils::SpinAdaptedContainer<Utils::Reference::Unrestricted, std::vector<Utils::Excitation>> excitations;
  excitations.alpha.resize(occupiedAlphaOrbitals.size() * virtualAlphaOrbitals.size());
  excitations.beta.resize(occupiedBetaOrbitals.size() * virtualBetaOrbitals.size());

  int iter = 0;
  for (int occ : occupiedAlphaOrbitals) {
    for (int vir : virtualAlphaOrbitals) {
      excitations.alpha[iter].occ = occ;
      excitations.alpha[iter].vir = vir;
      iter++;
    }
  }
  iter = 0;
  for (int occ : occupiedBetaOrbitals) {
    for (int vir : virtualBetaOrbitals) {
      excitations.beta[iter].occ = occ;
      excitations.beta[iter].vir = vir;
      iter++;
    }
  }
  return excitations;
}

namespace detail {
struct BetaExcitation {
  bool value = true;
};

inline auto occToVirLabel(int occupiedIndex, int virtualIndex) -> std::string {
  return std::to_string(occupiedIndex) + " -> " + std::to_string(virtualIndex);
}

inline auto occToVirLabel(int occupiedIndex, int virtualIndex, BetaExcitation isBeta) -> std::string {
  std::string spinLabel = isBeta.value ? "b" : "a";
  return std::to_string(occupiedIndex) + spinLabel + " -> " + std::to_string(virtualIndex) + spinLabel;
}
} // namespace detail

template<>
inline std::vector<std::string> TimeDependentUtils::generateExcitationsLabels<Utils::Reference::Restricted>(
    const Utils::SpinAdaptedContainer<Utils::Reference::Restricted, std::vector<Utils::Excitation>>& excitations) {
  int nExc = excitations.restricted.size();
  std::vector<std::string> labels;

  labels.reserve(nExc);

  for (int i = 0; i < nExc; ++i) {
    int indexOccupied = excitations.restricted[i].occ;
    int indexVirtual = excitations.restricted[i].vir;
    labels.push_back(detail::occToVirLabel(indexOccupied, indexVirtual));
  }
  return labels;
}

template<>
inline std::vector<std::string> TimeDependentUtils::generateExcitationsLabels<Utils::Reference::Unrestricted>(
    const Utils::SpinAdaptedContainer<Utils::Reference::Unrestricted, std::vector<Utils::Excitation>>& excitations) {
  int nAlphaExc = excitations.alpha.size();
  int nBetaExc = excitations.beta.size();

  std::vector<std::string> labels;
  labels.reserve(nAlphaExc + nBetaExc);

  for (int i = 0; i < nAlphaExc; ++i) {
    int indexOccupied = excitations.alpha[i].occ;
    int indexVirtual = excitations.alpha[i].vir;
    labels.push_back(detail::occToVirLabel(indexOccupied, indexVirtual, detail::BetaExcitation{false}));
  }

  for (int i = 0; i < nBetaExc; ++i) {
    int indexOccupied = excitations.beta[i].occ;
    int indexVirtual = excitations.beta[i].vir;
    labels.push_back(detail::occToVirLabel(indexOccupied, indexVirtual, detail::BetaExcitation{true}));
  }
  return labels;
}

inline std::vector<std::string> TimeDependentUtils::generateExcitationsLabels(const std::vector<Utils::Excitation>& excitations,
                                                                              const Eigen::Matrix<bool, -1, 1>& isBeta) {
  assert(static_cast<Eigen::Index>(excitations.size()) == isBeta.size() || isBeta.size() == 0);
  std::vector<std::string> labels;
  labels.reserve(excitations.size());

  for (unsigned int i = 0; i < excitations.size(); ++i) {
    int occ = excitations[i].occ;
    int vir = excitations[i].vir;
    if (isBeta.size() != 0) {
      labels.push_back(detail::occToVirLabel(occ, vir, detail::BetaExcitation{isBeta(i)}));
    }
    else {
      labels.push_back(detail::occToVirLabel(occ, vir));
    }
  }
  return labels;
}

template<>
inline std::vector<int> TimeDependentUtils::generateEnergyOrderMap<Utils::Reference::Restricted>(
    const Utils::SpinAdaptedContainer<Utils::Reference::Restricted, Eigen::VectorXd>& energyDifferenceVector) {
  std::multimap<double, int> energyOrderMap;
  std::vector<int> orderedIndices;
  const Eigen::VectorXd& energies = energyDifferenceVector.restricted;
  orderedIndices.reserve(energies.size());
  for (int i = 0; i < energies.size(); ++i) {
    energyOrderMap.insert({energies(i), i});
  }
  for (const auto& element : energyOrderMap) {
    orderedIndices.push_back(element.second);
  }
  return orderedIndices;
}

template<>
inline std::vector<int> TimeDependentUtils::generateEnergyOrderMap<Utils::Reference::Unrestricted>(
    const Utils::SpinAdaptedContainer<Utils::Reference::Unrestricted, Eigen::VectorXd>& energyDifferenceVector) {
  std::multimap<double, int> energyOrderMap;
  std::vector<int> orderedIndices;
  const Eigen::VectorXd& alpha = energyDifferenceVector.alpha;
  const Eigen::VectorXd& beta = energyDifferenceVector.beta;
  Eigen::VectorXd concatenated(alpha.size() + beta.size());
  concatenated << alpha, beta;
  orderedIndices.reserve(concatenated.size());
  for (int i = 0; i < concatenated.size(); ++i) {
    energyOrderMap.insert({concatenated(i), i});
  }
  for (const auto& element : energyOrderMap) {
    orderedIndices.push_back(element.second);
  }
  return orderedIndices;
}

template<typename Derived, typename DerivedResult>
inline void TimeDependentUtils::transformOrder(const Eigen::MatrixBase<Derived>& toTransform,
                                               const Eigen::MatrixBase<DerivedResult>& result,
                                               const std::vector<int>& orderVector, Direction direction) {
  assert(direction == Direction::To || direction == Direction::From);
  assert(static_cast<Eigen::Index>(orderVector.size()) == toTransform.rows());
  auto& castResult = const_cast<Eigen::MatrixBase<DerivedResult>&>(result).derived();
  castResult.resize(toTransform.rows(), toTransform.cols());
  if (direction == Direction::To) {
    for (unsigned int i = 0; i < orderVector.size(); ++i) {
      castResult.row(i) = toTransform.row(orderVector[i]);
    }
  }
  else { // Direction::From
    for (unsigned int i = 0; i < orderVector.size(); ++i) {
      castResult.row(orderVector[i]) = toTransform.row(i);
    }
  }
}

template<typename Type>
inline void TimeDependentUtils::transformOrder(const std::vector<Type>& toTransform, std::vector<Type>& result,
                                               const std::vector<int>& orderVector, Direction direction) {
  assert(direction == Direction::To || direction == Direction::From);
  assert(static_cast<Eigen::Index>(orderVector.size()) == static_cast<Eigen::Index>(toTransform.size()));
  result.resize(toTransform.size());
  if (direction == Direction::To) {
    for (unsigned int i = 0; i < orderVector.size(); ++i) {
      result[i] = toTransform[orderVector[i]];
    }
  }
  else { // Direction::From
    for (unsigned int i = 0; i < orderVector.size(); ++i) {
      result[orderVector[i]] = toTransform[i];
    }
  }
}

} // namespace Sparrow
} // namespace Scine

#endif // SPARROW_TIMEDEPENDENTUTILS_H