threecenter.h

Go to the documentation of this file.

/*
 *            Copyright 2009-2020 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_THREECENTER_H
#define VOTCA_XTP_THREECENTER_H
 
// Local VOTCA includes
#include "aobasis.h"
#include "eigen.h"
#include "symmetric_matrix.h"

 
namespace votca {
namespace xtp {
 
// due to different requirements for the data format for DFT and GW we have two
// different classes TCMatrix_gwbse and TCMatrix_dft which inherit from TCMatrix
class TCMatrix {
 
 public:
  TCMatrix() = default;
  virtual ~TCMatrix() = default;
  // the virtual destructor would otherwise suppress moving: moving the DFT
  // tensor (tens of GB) must not fall back to a copy
  TCMatrix(const TCMatrix&) = default;
  TCMatrix(TCMatrix&&) = default;
  TCMatrix& operator=(const TCMatrix&) = default;
  TCMatrix& operator=(TCMatrix&&) = default;
  Index Removedfunctions() const { return removedfunctions_; }
  enum class SpinChannel { Alpha, Beta };
 
 protected:
  Index removedfunctions_ = 0;
  Eigen::MatrixXd inv_sqrt_;
};

class TCMatrix_dft final : public TCMatrix {
 public:
  void Fill(const AOBasis& auxbasis, const AOBasis& dftbasis,
            double pair_threshold = 0.0, Index build_columns = kBuildColumns);

  Index size() const { return Index(data_.cols()); }
  Index basissize() const { return basissize_; }

  Eigen::MatrixXd FullMatrix(Index P) const;

  void FillFullMatrix(Index P, Eigen::MatrixXd& full) const;

  Eigen::VectorXd PackWeighted(const Eigen::MatrixXd& full) const;

  Eigen::MatrixXd Unpack(const Eigen::VectorXd& packed) const;

  const Eigen::MatrixXd& Data() const { return data_; }

  Index StoredPairs() const { return Index(data_.rows()); }
  Index AllPairs() const { return basissize_ * (basissize_ + 1) / 2; }
 
  double MetricSeconds() const { return metric_seconds_; }
 
 private:
  struct Run {
    Index col;
    Index row;
    Index length;
    Index offset;
  };

  std::vector<std::vector<Index>> SignificantPairs(
      const AOBasis& dftbasis, const Eigen::MatrixXd& metric,
      double pair_threshold) const;
  void SetupLayout(const AOBasis& dftbasis,
                   const std::vector<std::vector<Index>>& kept);
 
  Eigen::MatrixXd data_;
  std::vector<Run> runs_;
  // column mu starts at col_offset_[mu] in a data column
  std::vector<Index> col_offset_;
  Index basissize_ = 0;
  double metric_seconds_ = 0.0;

  static constexpr Index kBuildColumns = 512;
};
 
class TCMatrix_gwbse final : public TCMatrix {
 public:
  // Eigenvalue tolerance of the metric (Pseudo_InvSqrt_GWBSE) folded into
  // the stored integrals. Shared with EnvironmentScreening::Metric, which
  // must build the same T.
  static constexpr double metric_tolerance = 5e-7;
 
  // returns one level as a constant reference
  const Eigen::MatrixXd& operator[](Index i) const { return matrix_[i]; }
 
  // returns one level as a reference
  Eigen::MatrixXd& operator[](Index i) { return matrix_[i]; }
  // returns auxbasissize
  Index auxsize() const { return auxbasissize_; }
 
  Index get_mmin() const { return mmin_; }
 
  Index get_mmax() const { return mmax_; }
 
  Index get_nmin() const { return nmin_; }
 
  Index get_nmax() const { return nmax_; }
 
  Index msize() const { return mtotal_; }
 
  Index nsize() const { return ntotal_; }

  static Eigen::MatrixXd CoreExchange(const AOBasis& auxbasis,
                                      const AOBasis& dftbasis,
                                      const Eigen::MatrixXd& mos, Index ncore,
                                      Index qpmin, Index qpmax);
 
  void Initialize(Index basissize, Index mmin, Index mmax, Index nmin,
                  Index nmax);
 
  void Fill(const AOBasis& auxbasis, const AOBasis& dftbasis,
            const Eigen::MatrixXd& dft_orbitals);
  // Rebuilds ThreeCenterIntegrals, only works if the original basisobjects
  // still exist
  void Rebuild() { Fill(*auxbasis_, *dftbasis_, *dft_orbitals_); }
 
  void MultiplyRightWithAuxMatrix(const Eigen::MatrixXd& matrix);

  void MultiplyRightWithAuxMatrixLeading(const Eigen::MatrixXd& U,
                                         Index full_slices, Index lead);
  bool PartiallyRotated() const { return partially_rotated_; }
 
  // The metric that was folded into the stored integrals at Fill time:
  // T = Pseudo_InvSqrt_GWBSE, with T T^T = V^-1 (pseudo-inverse) and
  // T^T V T the projector onto the retained auxiliary functions. After
  // Fill, operator[] holds M = (mn|Q) T, so the bare Coulomb operator
  // exists only implicitly, as M M^T.
  //
  // Kept because anything that must act on the Coulomb interaction in
  // the same space as M needs it. The environment reaction field is the
  // case in point: a kernel B between auxiliary functions, taken as charge
  // densities, enters as M (T^T B T) M^T -- see EnvironmentScreening.
  //
  // IN THE FILL-TIME AUXILIARY BASIS. MultiplyRightWithAuxMatrix, which the
  // PPM and the BSE both call with an eigenvector matrix U, rotates the
  // auxiliary index of M. This matrix is not rotated with it: a kernel
  // built from it has to be brought into the current frame with
  // ToCurrentAuxFrame before it meets operator[]. Refilled, and so reset,
  // by every Fill and Rebuild.
  const Eigen::MatrixXd& InvSqrt() const { return inv_sqrt_; }
 
  // The auxiliary frame operator[] is currently in. Every
  // MultiplyRightWithAuxMatrix since the last Fill or Rebuild is recorded
  // here, so that operator[] holds M_fill * AuxFrame() with M_fill the
  // integrals as Fill left them. Empty means the identity: nothing has
  // rotated them yet.
  //
  // The bare Coulomb interaction does not care -- it is the identity in
  // every orthonormal frame, which is why nothing tracked this before. A
  // kernel K between the columns of M_fill does care: M_fill K M_fill^T is
  // M K' M^T only with K' = ToCurrentAuxFrame(K).
  const Eigen::MatrixXd& AuxFrame() const { return aux_frame_; }
 
  // For code that snapshots slices through operator[] and later writes
  // them back, which bypasses the bookkeeping above: restore the frame
  // that was current when the snapshot was taken, alongside the slices.
  void RestoreAuxFrame(const Eigen::MatrixXd& frame) {
    aux_frame_ = frame;
    aux_frame_known_ = true;
  }
 
  // K, a kernel between the columns of M_fill, as it must be applied to
  // operator[] now: F^-1 K F^-T with F = AuxFrame(), which is F^T K F for
  // the orthogonal eigenvector matrices the PPM and the BSE rotate with.
  // Throws if the frame is not known -- after a non-square
  // MultiplyRightWithAuxMatrix, which is not a change of frame.
  Eigen::MatrixXd ToCurrentAuxFrame(const Eigen::MatrixXd& K) const;
 
  // Whether AuxFrame() is orthogonal, so that the bare interaction M M^T is
  // what it was at Fill time. False once the index has been dressed.
  bool AuxFrameIsOrthogonal() const;

  void DressAuxIndex(const Eigen::MatrixXd& S);
  // Back to M_fill (times whatever rotations happened since), bare v.
  void UndressAuxIndex();
  bool Dressed() const { return dressing_.size() > 0; }

  void RotateOrbitals(const Eigen::MatrixXd& U, Index first, Index last);
 
 private:
  // store vector of matrices
  std::vector<Eigen::MatrixXd> matrix_;
 
  // band summation indices
  Index mmin_;
  Index mmax_;
  Index nmin_;
  Index nmax_;
  Index ntotal_;
  Index mtotal_;
  Index auxbasissize_;
 
  Eigen::MatrixXd aux_frame_;  // empty: identity
  bool aux_frame_known_ = true;
  bool partially_rotated_ = false;
  Eigen::MatrixXd dressing_;  // S of DressAuxIndex; empty: bare
 
  // M_fill -> M_fill A, whatever frame the integrals are in.
  void MultiplyInFillFrame(const Eigen::MatrixXd& A);
 
  const AOBasis* auxbasis_ = nullptr;
  const AOBasis* dftbasis_ = nullptr;
  const Eigen::MatrixXd* dft_orbitals_ = nullptr;
 
  void Fill3cMO(const AOBasis& auxbasis, const AOBasis& dftbasis,
                const Eigen::MatrixXd& dft_orbitals);
};
 
struct TCMatrix_gwbse_spin {
  TCMatrix_gwbse alpha;
  TCMatrix_gwbse beta;
 
  TCMatrix_gwbse& operator[](TCMatrix::SpinChannel spin) {
    return (spin == TCMatrix::SpinChannel::Alpha) ? alpha : beta;
  }
 
  const TCMatrix_gwbse& operator[](TCMatrix::SpinChannel spin) const {
    return (spin == TCMatrix::SpinChannel::Alpha) ? alpha : beta;
  }
};
 
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_THREECENTER_H