SQHamiltonian.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 "Basis/Transformations/BasisTransformable.hpp"
#include "Basis/Transformations/JacobiRotatable.hpp"
#include "Operator/SecondQuantized/GSQOneElectronOperator.hpp"
#include "Operator/SecondQuantized/GSQTwoElectronOperator.hpp"
#include "Operator/SecondQuantized/ModelHamiltonian/HubbardHamiltonian.hpp"
#include "Operator/SecondQuantized/RSQOneElectronOperator.hpp"
#include "Operator/SecondQuantized/RSQTwoElectronOperator.hpp"
#include "Operator/SecondQuantized/USQOneElectronOperator.hpp"
#include "Operator/SecondQuantized/USQTwoElectronOperator.hpp"
#include "QuantumChemical/spinor_tags.hpp"
#include "Utilities/type_traits.hpp"
 
#include <numeric>
#include <vector>
 
 
namespace GQCP {
 
 
/*
 *  MARK: SQHamiltonian implementation
 */
 
template <typename _ScalarSQOneElectronOperator, typename _ScalarSQTwoElectronOperator>
class SQHamiltonian:
    public BasisTransformable<SQHamiltonian<_ScalarSQOneElectronOperator, _ScalarSQTwoElectronOperator>>,
    public JacobiRotatable<SQHamiltonian<_ScalarSQOneElectronOperator, _ScalarSQTwoElectronOperator>> {
public:
    // The type of second-quantized one-electron operator underlying this Hamiltonian.
    using ScalarSQOneElectronOperator = _ScalarSQOneElectronOperator;
 
    // The type of second-quantized two-electron operator underlying this Hamiltonian.
    using ScalarSQTwoElectronOperator = _ScalarSQTwoElectronOperator;
 
    // Check if the spinor tags of the one- and two-electron operators match.
    static_assert(std::is_same<typename ScalarSQOneElectronOperator::SpinorTag, typename ScalarSQTwoElectronOperator::SpinorTag>::value, "The spinor tags of the one- and two-electron operators do not match.");
 
    // The spinor tag associated to this Hamiltonian.
    using SpinorTag = typename ScalarSQOneElectronOperator::SpinorTag;
 
    // Check if the scalar type of the parameters/matrix elements of the one-electron and two-electron operators are equal.
    static_assert(std::is_same<typename ScalarSQOneElectronOperator::Scalar, typename ScalarSQTwoElectronOperator::Scalar>::value, "The scalar type of the one- and two-electron parameters/matrix elements do not match.");
 
    // The scalar type used for a single parameter/matrix element: real or complex.
    using Scalar = typename ScalarSQOneElectronOperator::Scalar;
 
    // The type of 'this'
    using Self = SQHamiltonian<ScalarSQOneElectronOperator, ScalarSQTwoElectronOperator>;
 
    // The type of the one-particle density matrix that is naturally associated to the second-quantized Hamiltonian.
    using OneDM = typename OperatorTraits<ScalarSQOneElectronOperator>::OneDM;
 
    // The type of the two-particle density matrix that is naturally associated to the second-quantized Hamiltonian.
    using TwoDM = typename OperatorTraits<ScalarSQOneElectronOperator>::TwoDM;
 
    // The type of transformation that is naturally associated to the Hamiltonian.
    using Transformation = typename OperatorTraits<ScalarSQOneElectronOperator>::Transformation;
 
    // The type of Jacobi rotation for which the Jacobi rotation should be defined.
    using JacobiRotationType = typename JacobiRotatableTraits<ScalarSQOneElectronOperator>::JacobiRotationType;
 
 
private:
    // The total one-electron interaction operator, i.e. the core Hamiltonian.
    ScalarSQOneElectronOperator h;
 
    // The total two-electron interaction operator.
    ScalarSQTwoElectronOperator g;
 
    // The contributions to the total one-electron interaction operator.
    std::vector<ScalarSQOneElectronOperator> h_contributions;
 
    // The contributions to the total two-electron interaction operator.
    std::vector<ScalarSQTwoElectronOperator> g_contributions;
 
 
public:
    /*
     *  MARK: Constructors
     */
 
    SQHamiltonian(const std::vector<ScalarSQOneElectronOperator>& h_contributions, const std::vector<ScalarSQTwoElectronOperator>& g_contributions) :
        h_contributions {h_contributions},
        g_contributions {g_contributions} {
 
        // Check if the dimensions of the one- and two-electron operators are compatible.
        const std::invalid_argument dimension_error {"SQHamiltonian::SQHamiltonian(const std::vector<ScalarSQOneElectronOperator<Scalar>& h_contributions, const std::vector<ScalarSQTwoElectronOperator<Scalar>& g_contributions: The dimensions of the contributing operators are incompatible"};
 
        const auto dim = h_contributions[0].numberOfOrbitals();
        for (const auto& h : this->h_contributions) {
            if (h.numberOfOrbitals() != dim) {
                throw dimension_error;
            }
        }
 
        for (const auto& g : this->g_contributions) {
            if (g.numberOfOrbitals() != dim) {
                throw dimension_error;
            }
        }
 
 
        // Calculate the total one- and two-electron operators.
        this->h = std::accumulate(h_contributions.begin(), h_contributions.end(), ScalarSQOneElectronOperator::Zero(dim));
        this->g = std::accumulate(g_contributions.begin(), g_contributions.end(), ScalarSQTwoElectronOperator::Zero(dim));
    }
 
 
    SQHamiltonian(const ScalarSQOneElectronOperator& h, const ScalarSQTwoElectronOperator& g) :
        SQHamiltonian(std::vector<ScalarSQOneElectronOperator> {h},
                      std::vector<ScalarSQTwoElectronOperator> {g}) {}
 
 
    template <typename Z1 = Scalar, typename Z2 = SpinorTag>
    static enable_if_t<std::is_same<Z1, double>::value && std::is_same<Z2, RestrictedSpinorTag>::value, SQHamiltonian<ScalarSQOneElectronOperator, ScalarSQTwoElectronOperator>> FromHubbard(const HubbardHamiltonian<double>& hubbard_hamiltonian) {
 
        const auto h = hubbard_hamiltonian.core();
        const auto g = hubbard_hamiltonian.twoElectron();
 
        return Self {h, g};
    }
 
 
    template <typename Z = SpinorTag>
    static enable_if_t<std::is_same<Z, UnrestrictedSpinorTag>::value, SQHamiltonian<ScalarSQOneElectronOperator, ScalarSQTwoElectronOperator>> FromRestricted(const SQHamiltonian<ScalarRSQOneElectronOperator<Scalar>, ScalarRSQTwoElectronOperator<Scalar>>& r_hamiltonian) {
 
        const auto h_u = ScalarUSQOneElectronOperator<Scalar>::FromRestricted(r_hamiltonian.core());
        const auto g_u = ScalarUSQTwoElectronOperator<Scalar>::FromRestricted(r_hamiltonian.twoElectron());
 
        return Self {h_u, g_u};
    }
 
 
    template <typename Z = Scalar>
    static enable_if_t<std::is_same<Z, double>::value, Self> Random(const size_t dim) {
 
        const auto h = ScalarSQOneElectronOperator::Random(dim);
        const auto g = ScalarSQTwoElectronOperator::Random(dim);
 
        return Self {h, g};
    }
 
 
    template <typename Z1 = Scalar, typename Z2 = SpinorTag>
    static enable_if_t<std::is_same<Z1, double>::value && std::is_same<Z2, RestrictedSpinOrbitalTag>::value, SQHamiltonian<ScalarSQOneElectronOperator, ScalarSQTwoElectronOperator>> FromFCIDUMP(const std::string& fcidump_filename) {
 
        std::ifstream input_file_stream = validateAndOpen(fcidump_filename, "FCIDUMP");
 
 
        // Do the actual parsing
 
        //  Get the number of orbitals to check if it's a valid FCIDUMP file
        std::string start_line;  // first line contains orbitals and electron count
        std::getline(input_file_stream, start_line);
        std::stringstream linestream {start_line};
 
        size_t K = 0;
        char iter;
 
        while (linestream >> iter) {
            if (iter == '=') {
                linestream >> K;  // right here we have the number of orbitals
                break;            // we can finish reading the linestream after we found K
            }
        }
 
        if (K == 0) {
            throw std::invalid_argument("SQHamiltonian::FromFCIDUMP(std::string): The .FCIDUMP-file is invalid: could not read a number of orbitals.");
        }
 
 
        SquareMatrix<double> h_core = SquareMatrix<double>::Zero(K);
        SquareRankFourTensor<double> g {K};
        g.setZero();
 
        //  Skip 3 lines
        for (size_t counter = 0; counter < 3; counter++) {
            std::getline(input_file_stream, start_line);
        }
 
 
        //  Start reading in the one- and two-electron integrals
        double x;
        size_t i, j, a, b;
 
        std::string line;
        while (std::getline(input_file_stream, line)) {
            std::istringstream iss {line};
 
            // Based on what the values of the indices are, we can read one-electron integrals, two-electron integrals and the internuclear repulsion energy
            //  See also (http://hande.readthedocs.io/en/latest/manual/integrals.html)
            //  I think the documentation is a bit unclear for the two-electron integrals, but we can rest assured that FCIDUMP files give the two-electron integrals in CHEMIST's notation.
            iss >> x >> i >> a >> j >> b;
 
            //  Single-particle eigenvalues (skipped)
            if ((a == 0) && (j == 0) && (b == 0)) {
            }
 
            //  One-electron integrals (h_core)
            else if ((j == 0) && (b == 0)) {
                size_t p = i - 1;
                size_t q = a - 1;
                h_core(p, q) = x;
 
                // Apply the permutational symmetry for real orbitals
                h_core(q, p) = x;
            }
 
            //  Two-electron integrals are given in CHEMIST'S NOTATION, so just copy them over
            else if ((i > 0) && (a > 0) && (j > 0) && (b > 0)) {
                size_t p = i - 1;
                size_t q = a - 1;
                size_t r = j - 1;
                size_t s = b - 1;
                g(p, q, r, s) = x;
 
                // Apply the permutational symmetries for real orbitals
                g(p, q, s, r) = x;
                g(q, p, r, s) = x;
                g(q, p, s, r) = x;
 
                g(r, s, p, q) = x;
                g(s, r, p, q) = x;
                g(r, s, q, p) = x;
                g(s, r, q, p) = x;
            }
        }  // while loop
 
 
        return SQHamiltonian(ScalarSQOneElectronOperator(h_core), ScalarSQTwoElectronOperator(g));
    }
 
 
    /*
     *  MARK: Access
     */
 
    const ScalarSQOneElectronOperator& core() const { return this->h; }
 
    ScalarSQOneElectronOperator& core() { return this->h; }
 
    const std::vector<ScalarSQOneElectronOperator>& coreContributions() const { return this->h_contributions; }
 
    std::vector<ScalarSQOneElectronOperator>& coreContributions() { return this->h_contributions; }
 
    const ScalarSQTwoElectronOperator& twoElectron() const { return this->g; }
 
    ScalarSQTwoElectronOperator& twoElectron() { return this->g; }
 
    const std::vector<ScalarSQTwoElectronOperator>& twoElectronContributions() const { return this->g_contributions; }
 
    std::vector<ScalarSQTwoElectronOperator>& twoElectronContributions() { return this->g_contributions; }
 
 
    /*
     *  MARK: General information
     */
 
    size_t numberOfOrbitals() const { return this->core().numberOfOrbitals(); }
 
 
    /*
     *  MARK: Calculations
     */
 
    template <typename Z = SpinorTag>
    enable_if_t<std::is_same<Z, RestrictedSpinOrbitalTag>::value || std::is_same<Z, GeneralSpinorTag>::value, double> calculateEdmistonRuedenbergLocalizationIndex(const OrbitalSpace orbital_space) const {
 
        const auto& g_total_par = this->twoElectron().parameters();
 
        // TODO: When Eigen releases TensorTrace, use it here.
        double localization_index = 0.0;
        for (const auto& i : orbital_space.indices(OccupationType::k_occupied)) {
            localization_index += g_total_par(i, i, i, i);
        }
 
        return localization_index;
    }
 
 
    template <typename Z = SpinorTag>
    enable_if_t<std::is_same<Z, RestrictedSpinOrbitalTag>::value || std::is_same<Z, GeneralSpinorTag>::value, ScalarSQOneElectronOperator> calculateEffectiveOneElectronIntegrals() const {
 
        return this->core() + this->twoElectron().effectiveOneElectronPartition();
    }
 
 
    Scalar calculateExpectationValue(const OneDM& D, const TwoDM& d) const {
 
        // An SQHamiltonian contains ScalarSQOperators, so we access their expectation values with ().
        return this->core().calculateExpectationValue(D)() + this->twoElectron().calculateExpectationValue(d)();
    }
 
 
    template <typename Z = SpinorTag>
    enable_if_t<std::is_same<Z, RestrictedSpinOrbitalTag>::value || std::is_same<Z, GeneralSpinorTag>::value, SquareMatrix<Scalar>> calculateFockianMatrix(const OneDM& D, const TwoDM& d) const {
 
        // An SQHamiltonian contains ScalarSQOperators, so we access their Fockian matrices with (0).
        return this->core().calculateFockianMatrix(D, d)() + this->twoElectron().calculateFockianMatrix(D, d)();
    }
 
 
    template <typename Z = SpinorTag>
    enable_if_t<std::is_same<Z, RestrictedSpinOrbitalTag>::value || std::is_same<Z, GeneralSpinorTag>::value, SquareRankFourTensor<Scalar>> calculateSuperFockianMatrix(const OneDM& D, const TwoDM& d) const {
 
        // An SQHamiltonian contains ScalarSQOperators, so we access their Fockian matrices with (0).
        return this->core().calculateSuperFockianMatrix(D, d)().Eigen() + this->twoElectron().calculateSuperFockianMatrix(D, d)().Eigen();  // We have to call .Eigen() because operator+ isn't fully enabled on SquareRankFourTensor.
    }
 
 
    template <typename Z = SpinorTag>
    enable_if_t<std::is_same<Z, GeneralSpinorTag>::value, ScalarSQOneElectronOperator> calculateInactiveFockian(const OrbitalSpace orbital_space) const {
 
        const auto& h_par = this->core().parameters();
        const auto& g_par = this->twoElectron().parameters();
 
        // A KISS implementation of the calculation of the (general) inactive Fockian matrix
        auto F_par = h_par;  // one-electron part
 
        // Two-electron part
        for (const auto& p : orbital_space.indices()) {
            for (const auto& q : orbital_space.indices()) {
 
                for (const auto& i : orbital_space.indices(OccupationType::k_occupied)) {
                    F_par(p, q) += g_par(p, q, i, i) - g_par(p, i, i, q);
                }
            }
        }  // F elements loop
 
        return ScalarSQOneElectronOperator {F_par};
    }
 
 
    template <typename Z = SpinorTag>
    enable_if_t<std::is_same<Z, RestrictedSpinOrbitalTag>::value, ScalarSQOneElectronOperator> calculateInactiveFockian(const OrbitalSpace orbital_space) const {
        const auto& h_par = this->core().parameters();
        const auto& g_par = this->twoElectron().parameters();
 
 
        // A KISS implementation of the calculation of the (restricted) inactive Fockian matrix
        auto F_par = h_par;  // one-electron part
 
        // Two-electron part
        for (const auto& p : orbital_space.indices()) {
            for (const auto& q : orbital_space.indices()) {
 
                for (const auto& i : orbital_space.indices(OccupationType::k_occupied)) {
                    F_par(p, q) += 2 * g_par(p, q, i, i) - g_par(p, i, i, q);
                }
            }
        }  // F elements loop
 
        return ScalarSQOneElectronOperator {F_par};
    }
 
 
    /*
     *  MARK: Conforming to BasisTransformable
     */
 
    Self transformed(const Transformation& T) const override {
 
        auto result = *this;
 
        // Transform the one and two-electron contributions.
        for (auto& h : result.coreContributions()) {
            h.transform(T);
        }
 
        for (auto& g : result.twoElectronContributions()) {
            g.transform(T);
        }
 
        // Transform the total one- and two-electron interactions.
        result.core().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 JacobiRotationType& jacobi_rotation) const override {
 
        auto result = *this;
 
        // Transform the one and two-electron contributions.
        for (auto& h : result.coreContributions()) {
            h.rotate(jacobi_rotation);
        }
 
        for (auto& g : result.twoElectronContributions()) {
            g.rotate(jacobi_rotation);
        }
 
        // Transform the total one- and two-electron interactions.
        result.core().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: Operations related to one-electron operators
     */
 
    Self& operator+=(const ScalarSQOneElectronOperator& sq_one_op) {
 
        // Update the one-electron contributions and the total one-electron operator.
        this->h_contributions.push_back(sq_one_op);
        this->h += sq_one_op;
 
        return *this;
    }
 
 
    friend Self operator+(Self lhs, const ScalarSQOneElectronOperator& rhs) {
        lhs += rhs;
        return lhs;
    }
 
 
    friend Self operator+(const ScalarSQOneElectronOperator& lhs, Self rhs) {
        return rhs += lhs;
    }
 
 
    Self& operator-=(const ScalarSQOneElectronOperator& sq_one_op) {
 
        *this += -sq_one_op;
        return *this;
    }
 
 
    friend Self operator-(Self lhs, const ScalarSQOneElectronOperator& rhs) {
        lhs -= rhs;
        return lhs;
    }
 
 
    /*
     *  MARK: Operations related to two-electron operators
     */
 
    Self& operator+=(const ScalarSQTwoElectronOperator& sq_two_op) {
 
        // Update the one-electron contributions and the total one-electron operator.
        this->g_contributions.push_back(sq_two_op);
        this->g += sq_two_op;
 
        return *this;
    }
 
 
    friend Self operator+(Self lhs, const ScalarSQTwoElectronOperator& rhs) {
        lhs += rhs;
        return lhs;
    }
 
 
    friend Self operator+(const ScalarSQTwoElectronOperator& lhs, Self rhs) {
        return rhs += lhs;
    }
 
 
    Self& operator-=(const ScalarSQTwoElectronOperator& sq_two_op) {
 
        *this += -sq_two_op;
        return *this;
    }
 
 
    friend Self operator-(Self lhs, const ScalarSQTwoElectronOperator& rhs) {
        lhs -= rhs;
        return lhs;
    }
};
 
 
/*
 *  MARK: Operator traits
 */
 
template <typename ScalarSQOneElectronOperator, typename ScalarSQTwoElectronOperator>
struct OperatorTraits<SQHamiltonian<ScalarSQOneElectronOperator, ScalarSQTwoElectronOperator>> {
 
    // The type of transformation that is naturally associated to the second-quantized Hamiltonian.
    using Transformation = typename OperatorTraits<ScalarSQOneElectronOperator>::Transformation;
 
    // The type of the one-particle density matrix that is naturally associated to the second-quantized Hamiltonian.
    using OneDM = typename OperatorTraits<ScalarSQOneElectronOperator>::OneDM;
 
    // The type of the two-particle density matrix that is naturally associated to the second-quantized Hamiltonian.
    using TwoDM = typename OperatorTraits<ScalarSQOneElectronOperator>::TwoDM;
};
 
 
/*
 *  MARK: BasisTransformableTraits
 */
 
template <typename ScalarSQOneElectronOperator, typename ScalarSQTwoElectronOperator>
struct BasisTransformableTraits<SQHamiltonian<ScalarSQOneElectronOperator, ScalarSQTwoElectronOperator>> {
 
    // The type of the transformation for which the basis transformation should be defined.
    using Transformation = typename OperatorTraits<ScalarSQOneElectronOperator>::Transformation;
};
 
 
/*
 *  MARK: JacobiRotatableTraits
 */
 
template <typename ScalarSQOneElectronOperator, typename ScalarSQTwoElectronOperator>
struct JacobiRotatableTraits<SQHamiltonian<ScalarSQOneElectronOperator, ScalarSQTwoElectronOperator>> {
 
    // The type of Jacobi rotation for which the Jacobi rotation should be defined.
    using JacobiRotationType = typename JacobiRotatableTraits<ScalarSQOneElectronOperator>::JacobiRotationType;
};
 
 
/*
 *  MARK: Convenience aliases
 */
 
 
// An `SQHamiltonian` related to restricted spin-orbitals. See `RestrictedSpinOrbitalTag`.
template <typename Scalar>
using RSQHamiltonian = SQHamiltonian<ScalarRSQOneElectronOperator<Scalar>, ScalarRSQTwoElectronOperator<Scalar>>;
 
 
// An `SQHamiltonian` related to unrestricted spin-orbitals. See `UnrestrictedSpinOrbitalTag`.
template <typename Scalar>
using USQHamiltonian = SQHamiltonian<ScalarUSQOneElectronOperator<Scalar>, ScalarUSQTwoElectronOperator<Scalar>>;
 
 
// An `SQHamiltonian` related to general spinors. See `GeneralSpinorTag`.
template <typename Scalar>
using GSQHamiltonian = SQHamiltonian<ScalarGSQOneElectronOperator<Scalar>, ScalarGSQTwoElectronOperator<Scalar>>;
 
 
}  // namespace GQCP