USQTwoElectronOperator.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/Transformations/DoublySpinResolvedBasisTransformable.hpp"
#include "Basis/Transformations/DoublySpinResolvedJacobiRotatable.hpp"
#include "Basis/Transformations/UJacobiRotation.hpp"
#include "Basis/Transformations/UTransformation.hpp"
#include "DensityMatrix/SpinResolved1DM.hpp"
#include "DensityMatrix/SpinResolved2DM.hpp"
#include "Mathematical/Functions/VectorSpaceArithmetic.hpp"
#include "Operator/SecondQuantized/MixedUSQTwoElectronOperatorComponent.hpp"
#include "Operator/SecondQuantized/PureUSQTwoElectronOperatorComponent.hpp"
#include "Operator/SecondQuantized/RSQTwoElectronOperator.hpp"
#include "Operator/SecondQuantized/USQOneElectronOperator.hpp"
#include "QuantumChemical/DoublySpinResolvedBase.hpp"
#include "QuantumChemical/spinor_tags.hpp"
 
 
namespace GQCP {
 
 
template <typename _Scalar, typename _Vectorizer>
class USQTwoElectronOperator:
    public DoublySpinResolvedBase<PureUSQTwoElectronOperatorComponent<_Scalar, _Vectorizer>, MixedUSQTwoElectronOperatorComponent<_Scalar, _Vectorizer>, USQTwoElectronOperator<_Scalar, _Vectorizer>>,
    public DoublySpinResolvedBasisTransformable<USQTwoElectronOperator<_Scalar, _Vectorizer>>,
    public DoublySpinResolvedJacobiRotatable<USQTwoElectronOperator<_Scalar, _Vectorizer>>,
    public VectorSpaceArithmetic<USQTwoElectronOperator<_Scalar, _Vectorizer>, _Scalar> {
public:
    // The scalar type used for a single parameter: real or complex.
    using Scalar = _Scalar;
 
    // The type of the vectorizer that relates a one-dimensional storage of matrices to the tensor structure of two-electron operators. This distinction is carried over from `USQOneElectronOperator`.
    using Vectorizer = _Vectorizer;
 
    // The type of 'this'.
    using Self = USQTwoElectronOperator<Scalar, Vectorizer>;
 
    // The type of the transformation that is naturally related to an unrestricted two-electron operator.
    using Transformation = UTransformation<Scalar>;
 
    // The type of Jacobi rotation that is naturally related to an unrestricted two-electron operator.
    using JacobiRotationType = UJacobiRotation;
 
    // The spinor tag corresponding to a `USQTwoElectronOperator`.
    using SpinorTag = UnrestrictedSpinOrbitalTag;
 
    // The type components this spin resolved object is made of.
    using PureComponentType = typename DoublySpinResolvedBase<PureUSQTwoElectronOperatorComponent<Scalar, Vectorizer>, MixedUSQTwoElectronOperatorComponent<Scalar, Vectorizer>, USQTwoElectronOperator<Scalar, Vectorizer>>::Pure;
    using MixedComponentType = typename DoublySpinResolvedBase<PureUSQTwoElectronOperatorComponent<Scalar, Vectorizer>, MixedUSQTwoElectronOperatorComponent<Scalar, Vectorizer>, USQTwoElectronOperator<Scalar, Vectorizer>>::Mixed;
 
 
public:
    /*
     *  MARK: Constructors
     */
 
    // Inherit `DoublySpinResolvedBase`'s constructors.
    using DoublySpinResolvedBase<PureUSQTwoElectronOperatorComponent<Scalar, Vectorizer>, MixedUSQTwoElectronOperatorComponent<Scalar, Vectorizer>, USQTwoElectronOperator<Scalar, Vectorizer>>::DoublySpinResolvedBase;
 
 
    template <size_t N>
    USQTwoElectronOperator(const std::array<SquareRankFourTensor<Scalar>, N>& gs_aa, const std::array<SquareRankFourTensor<Scalar>, N>& gs_ab, const std::array<SquareRankFourTensor<Scalar>, N>& gs_ba, const std::array<SquareRankFourTensor<Scalar>, N>& gs_bb, const Vectorizer& vectorizer) :
 
        // Encapsulate the array of matrix representations in the different spin-components components, and put them together to form the `USQTwoElectronOperator`.
        DoublySpinResolvedBase<PureUSQTwoElectronOperatorComponent<Scalar, Vectorizer>, MixedUSQTwoElectronOperatorComponent<Scalar, Vectorizer>, USQTwoElectronOperator<Scalar, Vectorizer>>(
            PureUSQTwoElectronOperatorComponent<Scalar, Vectorizer> {StorageArray<SquareRankFourTensor<Scalar>, Vectorizer> {gs_aa, vectorizer}},
            MixedUSQTwoElectronOperatorComponent<Scalar, Vectorizer> {StorageArray<SquareRankFourTensor<Scalar>, Vectorizer> {gs_ab, vectorizer}},
            MixedUSQTwoElectronOperatorComponent<Scalar, Vectorizer> {StorageArray<SquareRankFourTensor<Scalar>, Vectorizer> {gs_ba, vectorizer}},
            PureUSQTwoElectronOperatorComponent<Scalar, Vectorizer> {StorageArray<SquareRankFourTensor<Scalar>, Vectorizer> {gs_bb, vectorizer}}) {
 
        // Check if the given tensor representations have the same dimensions, for each spin part.
        const auto dimension_of_first_aa = gs_aa[0].dimension();
        const auto dimension_of_first_ab = gs_ab[0].dimension();
        const auto dimension_of_first_ba = gs_ba[0].dimension();
        const auto dimension_of_first_bb = gs_bb[0].dimension();
 
        for (size_t i = 1; i < N; i++) {
 
            const auto dimension_of_ith_aa = gs_aa[i].dimension();
            const auto dimension_of_ith_ab = gs_ab[i].dimension();
            const auto dimension_of_ith_ba = gs_ba[i].dimension();
            const auto dimension_of_ith_bb = gs_bb[i].dimension();
 
            if ((dimension_of_first_aa != dimension_of_ith_aa) || (dimension_of_first_ab != dimension_of_ith_ab) || (dimension_of_first_ba != dimension_of_ith_ba) || (dimension_of_first_bb != dimension_of_ith_bb)) {
                throw std::invalid_argument("USQTwoElectronOperator(const std::array<SquareRankFourTensor<Scalar>, N>&, const std::array<SquareRankFourTensor<Scalar>, N>&, const std::array<SquareRankFourTensor<Scalar>, N>&, const std::array<SquareRankFourTensor<Scalar>, N>&): The given tensor representations do not have the same dimensions for either the alpha, beta or one of the mixed components.");
            }
        }
    }
 
 
    template <typename Z = Vectorizer>
    USQTwoElectronOperator(const SquareRankFourTensor<Scalar>& g_aa, const SquareRankFourTensor<Scalar>& g_ab, const SquareRankFourTensor<Scalar>& g_ba, const SquareRankFourTensor<Scalar>& g_bb,
                           typename std::enable_if<std::is_same<Z, ScalarVectorizer>::value>::type* = 0) :
        USQTwoElectronOperator(std::array<SquareRankFourTensor<Scalar>, 1> {g_aa}, std::array<SquareRankFourTensor<Scalar>, 1> {g_ab}, std::array<SquareRankFourTensor<Scalar>, 1> {g_ba}, std::array<SquareRankFourTensor<Scalar>, 1> {g_bb}, ScalarVectorizer()) {}
 
 
    USQTwoElectronOperator() :
        USQTwoElectronOperator(USQTwoElectronOperator<Scalar, Vectorizer>::Zero(0)) {}
 
 
    /*
     *  MARK: Named constructors
     */
 
    static USQTwoElectronOperator<Scalar, Vectorizer> Zero(const size_t dim) {
 
        const auto zero_component_pure = PureUSQTwoElectronOperatorComponent<Scalar, Vectorizer>::Zero(dim);
        const auto zero_component_mixed = MixedUSQTwoElectronOperatorComponent<Scalar, Vectorizer>::Zero(dim);
 
        return USQTwoElectronOperator<Scalar, Vectorizer> {zero_component_pure, zero_component_mixed, zero_component_mixed, zero_component_pure};
    }
 
 
    static USQTwoElectronOperator<Scalar, Vectorizer> FromRestricted(const RSQTwoElectronOperator<Scalar, Vectorizer>& g_restricted) {
 
        return USQTwoElectronOperator<Scalar, Vectorizer> {g_restricted.alphaAlpha(), g_restricted.alphaBeta(), g_restricted.betaAlpha(), g_restricted.betaBeta()};
    }
 
 
    /*
     *  MARK: Calculations
     */
 
    StorageArray<Scalar, Vectorizer> calculateExpectationValue(const SpinResolved2DM<Scalar>& d) const {
 
        // Calculate the sum of all spin contributions.
        // Unfortunately, we can't give `StorageArray` out-of-the-box vector-space arithmetic (since the default scalar type for scalar multiplication is unknown), so we'll have to do the summation ourselves.
        const auto expectation_value_elements_aa = this->alphaAlpha().calculateExpectationValue(d.alphaAlpha()).elements();
        const auto expectation_value_elements_ab = this->alphaBeta().calculateExpectationValue(d.alphaBeta()).elements();
        const auto expectation_value_elements_ba = this->betaAlpha().calculateExpectationValue(d.betaAlpha()).elements();
        const auto expectation_value_elements_bb = this->betaBeta().calculateExpectationValue(d.betaBeta()).elements();
 
        auto result_elements = expectation_value_elements_aa;
        const auto vectorizer = this->alphaAlpha().vectorizer();
        for (size_t i = 0; i < vectorizer.numberOfElements(); i++) {
            result_elements[i] += expectation_value_elements_ab[i] + expectation_value_elements_ba[i] + expectation_value_elements_bb[i];
        }
 
        return StorageArray<Scalar, Vectorizer> {result_elements, vectorizer};
    }
 
 
    /*
     *  MARK: General information
     */
 
    size_t numberOfOrbitals(const Spin sigma, const Spin tau) const {
 
        if (sigma == Spin::alpha && tau == Spin::beta) {
            return this->alphaAlpha().numberOfOrbitals();
        } else if (sigma == Spin::alpha && tau == Spin::beta) {
            return this->alphaBeta().numberOfOrbitals();
        } else if (sigma == Spin::beta && tau == Spin::alpha) {
            return this->betaAlpha().numberOfOrbitals();
        } else {
            return this->betaBeta().numberOfOrbitals();
        }
    }
 
 
    /*
     *  @return The number of orbital related to the alpha-alpha part of the unrestricted two-electron operator.
     *
     *  @note It is advised to only use this API when it is known that all spin-components of the two-electron operator are equal.
     */
    size_t numberOfOrbitals() const { return this->alphaAlpha().numberOfOrbitals(); }
 
 
    // Since `rotate` and `rotated` are both defined in `DoublySpinResolvedBasisTransformable` and `DoublySpinResolvedJacobiRotatable`, we have to explicitly enable these methods here.
 
    // Allow the `rotate` method from `DoublySpinResolvedBasisTransformable`, since there's also a `rotate` from `DoublySpinResolvedJacobiRotatable`.
    using DoublySpinResolvedBasisTransformable<Self>::rotate;
 
    // Allow the `rotated` method from `DoublySpinResolvedBasisTransformable`, since there's also a `rotated` from `DoublySpinResolvedJacobiRotatable`.
    using DoublySpinResolvedBasisTransformable<Self>::rotated;
 
    // Allow the `rotate` method from `DoublySpinResolvedJacobiRotatable`, since there's also a `rotate` from `DoublySpinResolvedBasisTransformable`.
    using DoublySpinResolvedJacobiRotatable<Self>::rotate;
 
    // Allow the `rotated` method from `DoublySpinResolvedJacobiRotatable`, since there's also a `rotated` from `DoublySpinResolvedBasisTransformable`.
    using DoublySpinResolvedJacobiRotatable<Self>::rotated;
 
 
    /*
     *  MARK: Conforming to VectorSpaceArithmetic
     */
 
    Self& operator+=(const Self& rhs) override {
 
        // Add the spin-components.
        this->alphaAlpha() += rhs.alphaAlpha();
        this->alphaBeta() += rhs.alphaBeta();
        this->betaAlpha() += rhs.betaAlpha();
        this->betaBeta() += rhs.betaBeta();
 
        return *this;
    }
 
 
    Self& operator*=(const Scalar& a) override {
 
        // Multiply the spin-components with the scalar.
        this->alphaAlpha() *= a;
        this->alphaBeta() *= a;
        this->betaAlpha() *= a;
        this->betaBeta() *= a;
 
        return *this;
    }
};
 
 
/*
 *  MARK: Convenience aliases
 */
 
// A scalar-like USQTwoElectronOperator, i.e. with scalar-like access.
template <typename Scalar>
using ScalarUSQTwoElectronOperator = USQTwoElectronOperator<Scalar, ScalarVectorizer>;
 
// A vector-like USQTwoElectronOperator, i.e. with vector-like access.
template <typename Scalar>
using VectorUSQTwoElectronOperator = USQTwoElectronOperator<Scalar, VectorVectorizer>;
 
// A matrix-like USQTwoElectronOperator, i.e. with matrix-like access.
template <typename Scalar>
using MatrixUSQTwoElectronOperator = USQTwoElectronOperator<Scalar, MatrixVectorizer>;
 
// A tensor-like USQTwoElectronOperator, i.e. with tensor-like access.
template <typename Scalar, size_t N>
using TensorUSQTwoElectronOperator = USQTwoElectronOperator<Scalar, TensorVectorizer<N>>;
 
 
template <typename Scalar, typename Vectorizer>
struct OperatorTraits<USQTwoElectronOperator<Scalar, Vectorizer>> {
 
    // A type that corresponds to the scalar version of the associated unrestricted one-electron operator type.
    using ScalarOperator = ScalarUSQTwoElectronOperator<Scalar>;
 
    // The type of transformation that is naturally associated to an unrestricted two-electron operator.
    using Transformation = UTransformation<Scalar>;
 
    // The type of the one-particle density matrix that is naturally associated an unrestricted two-electron operator.
    using OneDM = SpinResolved1DM<Scalar>;
 
    // The type of the two-particle density matrix that is naturally associated an unrestricted two-electron operator.
    using TwoDM = SpinResolved2DM<Scalar>;
};
 
 
/*
 *  MARK: BasisTransformableTraits
 */
 
template <typename Scalar, typename Vectorizer>
struct BasisTransformableTraits<USQTwoElectronOperator<Scalar, Vectorizer>> {
 
    // The type of transformation that is naturally related to a `USQTwoElectronOperator`.
    using Transformation = UTransformation<Scalar>;
};
 
 
/*
 *  MARK: JacobiRotatableTraits
 */
 
template <typename Scalar, typename Vectorizer>
struct JacobiRotatableTraits<USQTwoElectronOperator<Scalar, Vectorizer>> {
 
    // The type of Jacobi rotation for which the Jacobi rotation should be defined.
    using JacobiRotationType = UJacobiRotation;
};
 
 
/*
 *  MARK: One-electron operator products
 */
 
template <typename _Scalar>
class ScalarUSQOneElectronOperatorProduct:
    public VectorSpaceArithmetic<ScalarUSQOneElectronOperatorProduct<_Scalar>, _Scalar>,
    public BasisTransformable<ScalarUSQOneElectronOperatorProduct<_Scalar>>,
    public JacobiRotatable<ScalarUSQOneElectronOperatorProduct<_Scalar>> {
public:
    // The scalar type of one of the matrix elements: real or complex.
    using Scalar = _Scalar;
 
    // The type of 'this'.
    using Self = ScalarUSQOneElectronOperatorProduct<Scalar>;
 
 
private:
    // The one-electron part of the product.
    ScalarUSQOneElectronOperator<Scalar> o;
 
    // The two-electron part of the product. Since it is expressed as a `SQTwoElectronOperator`.
    ScalarUSQTwoElectronOperator<Scalar> t;
 
 
public:
    /*
     *  MARK: Constructors
     */
 
    ScalarUSQOneElectronOperatorProduct(const ScalarUSQOneElectronOperator<Scalar>& o, const ScalarUSQTwoElectronOperator<Scalar>& t) :
        o {o},
        t {t} {}
 
 
    /*
     *  MARK: Access
     */
 
    const ScalarUSQOneElectronOperator<Scalar>& oneElectron() const { return this->o; }
 
    const ScalarUSQTwoElectronOperator<Scalar>& twoElectron() const { return this->t; }
 
    ScalarUSQOneElectronOperator<Scalar>& oneElectron() { return this->o; }
 
    ScalarUSQTwoElectronOperator<Scalar>& twoElectron() { return this->t; }
 
 
    /*
     *  MARK: Conforming to `VectorSpaceArithmetic`
     */
 
    ScalarUSQOneElectronOperatorProduct& operator+=(const ScalarUSQOneElectronOperatorProduct& rhs) override {
 
        // Add the one- and two-electron operator parts.
        this->oneElectron() += rhs.oneElectron();
        this->twoElectron() += rhs.twoElectron();
 
        return *this;
    }
 
 
    ScalarUSQOneElectronOperatorProduct& operator*=(const Scalar& a) override {
 
        // Multiply the one- and two-electron operator parts with the scalar.
        this->oneElectron() *= a;
        this->twoElectron() *= a;
 
        return *this;
    }
 
 
    /*
     *  MARK: Conforming to BasisTransformable
     */
 
    Self transformed(const UTransformation<Scalar>& T) const override {
 
        auto result = *this;
 
        // Transform the one- and two-electron contributions.
        result.oneElectron().transform(T);
        result.twoElectron().transform(T);
 
        return result;
    }
 
 
    // Allow the `rotate` method from `BasisTransformable`, since there's also a `rotate` from `JacobiRotatable`.
    using BasisTransformable<Self>::rotate;
 
    // Allow the `rotated` method from `BasisTransformable`, since there's also a `rotated` from `JacobiRotatable`.
    using BasisTransformable<Self>::rotated;
 
 
    /*
     *  MARK: Conforming to JacobiRotatable
     */
 
    Self rotated(const UJacobiRotation& jacobi_rotation) const override {
 
        auto result = *this;
 
        // Transform the total one- and two-electron contributions.
        result.oneElectron().rotate(jacobi_rotation);
        result.twoElectron().rotate(jacobi_rotation);
 
        return result;
    }
 
    // Allow the `rotate` method from `JacobiRotatable`, since there's also a `rotate` from `BasisTransformable`.
    using JacobiRotatable<Self>::rotate;
 
 
    /*
     *  MARK: Expectation value
     */
 
    Scalar calculateExpectationValue(const SpinResolved1DM<Scalar>& D, const SpinResolved2DM<Scalar>& d) const {
 
        // We can access the expectation value of a scalar operator using an empty call.
        return this->oneElectron().calculateExpectationValue(D)() + this->twoElectron().calculateExpectationValue(d)();
    }
};
 
 
/*
 *  MARK: `BasisTransformableTraits` for `ScalarUSQOneElectronOperatorProduct`.
 */
 
template <typename Scalar>
struct BasisTransformableTraits<ScalarUSQOneElectronOperatorProduct<Scalar>> {
 
    // The type of transformation matrix that is naturally associated to a `ScalarUSQOneElectronOperatorProduct`.
    using Transformation = UTransformation<Scalar>;
};
 
 
/*
 *  MARK: `JacobiRotatableTraits` for `ScalarUSQOneElectronOperatorProduct`.
 */
 
template <typename Scalar>
struct JacobiRotatableTraits<ScalarUSQOneElectronOperatorProduct<Scalar>> {
 
    // The type of Jacobi rotation for which the Jacobi rotation should be defined.
    using JacobiRotationType = UJacobiRotation;
};
 
 
}  // namespace GQCP