T1Amplitudes.hpp

Go to the documentation of this file.

// This file is part of GQCG-GQCP.
//
// Copyright (C) 2017-2020  the GQCG developers
//
// GQCG-GQCP is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// GQCG-GQCP is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with GQCG-GQCP.  If not, see <http://www.gnu.org/licenses/>.
 
#pragma once
 
 
#include "Basis/SpinorBasis/OrbitalSpace.hpp"
#include "Mathematical/Functions/VectorSpaceArithmetic.hpp"
#include "Mathematical/Representation/ImplicitMatrixSlice.hpp"
#include "Operator/SecondQuantized/SQHamiltonian.hpp"
#include "Utilities/type_traits.hpp"
 
 
namespace GQCP {
 
 
template <typename _Scalar>
class T1Amplitudes:
    public VectorSpaceArithmetic<T1Amplitudes<_Scalar>, _Scalar> {
 
public:
    // The scalar type that represents one of the amplitudes.
    using Scalar = _Scalar;
 
    // The type of 'this'.
    using Self = T1Amplitudes<Scalar>;
 
 
private:
    // The orbital space which encapsulates the occupied-virtual separation.
    OrbitalSpace orbital_space;
 
    // The T1-amplitudes as an implicit matrix, implementing intuitive operator(i,a) calls.
    ImplicitMatrixSlice<Scalar> t;
 
 
public:
    /*
     *  MARK: Constructors
     */
 
    T1Amplitudes(const ImplicitMatrixSlice<Scalar>& t, const OrbitalSpace& orbital_space) :
        orbital_space {orbital_space},
        t {t} {}
 
 
    T1Amplitudes(const ImplicitMatrixSlice<Scalar>& t, const size_t N, const size_t M) :
        T1Amplitudes(t, OrbitalSpace::Implicit({{OccupationType::k_occupied, N}, {OccupationType::k_virtual, M - N}})) {}
 
 
    /*
     *  MARK: Named constructors
     */
 
    static T1Amplitudes<Scalar> Perturbative(const SquareMatrix<Scalar>& f, const OrbitalSpace& orbital_space) {
 
        // Zero-initialize a matrix representation for the (occupied-virtual) T1-amplitudes t_i^a.
        auto t1 = orbital_space.initializeRepresentableObjectFor<Scalar>(OccupationType::k_occupied, OccupationType::k_virtual);
 
        // Provide the perturbative T1-amplitudes, by setting the RHS amplitudes in Stanton1991, equation (1) to zero.
        for (const auto& i : orbital_space.indices(OccupationType::k_occupied)) {
            for (const auto& a : orbital_space.indices(OccupationType::k_virtual)) {
                const auto denominator = f(i, i) - f(a, a);
 
                t1(i, a) = f(i, a) / denominator;
            }
        }
 
        return T1Amplitudes<Scalar>(t1, orbital_space);
    }
 
 
    /*
     *  MARK: Access.
     */
 
    Scalar operator()(const size_t i, const size_t a) const { return this->t(i, a); }
 
    Scalar& operator()(const size_t i, const size_t a) { return this->t(i, a); }
 
    const ImplicitMatrixSlice<Scalar>& asImplicitMatrixSlice() const { return this->t; }
 
    const OrbitalSpace& orbitalSpace() const { return this->orbital_space; }
 
 
    /*
     *  MARK: Linear algebra
     */
 
    Scalar norm() const { return this->asImplicitMatrixSlice().asMatrix().norm(); }
 
 
    /*
     *  MARK: Vector space arithmetic
     */
 
    Self& operator+=(const Self& rhs) override {
 
        // Prepare some variables.
        const auto& row_map = this->asImplicitMatrixSlice().rowIndexMap();
        const auto& col_map = this->asImplicitMatrixSlice().columnIndexMap();
 
        // Sum the matrix representations.
        const MatrixX<Scalar> t_sum_dense = this->asImplicitMatrixSlice().asMatrix() + rhs.asImplicitMatrixSlice().asMatrix();
        const ImplicitMatrixSlice<Scalar> t_sum_slice {row_map, col_map, t_sum_dense};
 
        this->t = t_sum_slice;
 
        return *this;
    }
 
 
    Self& operator*=(const Scalar& a) override {
 
        // Prepare some variables.
        const auto& row_map = this->asImplicitMatrixSlice().rowIndexMap();
        const auto& col_map = this->asImplicitMatrixSlice().columnIndexMap();
 
        // Multiply the matrix representation.
        const MatrixX<Scalar> t_multiplied_dense = a * this->asImplicitMatrixSlice().asMatrix();
        const ImplicitMatrixSlice<Scalar> t_multiplied_slice {row_map, col_map, t_multiplied_dense};
 
        this->t = t_multiplied_slice;
 
        return *this;
    }
};
 
}  // namespace GQCP