Low-level CI improvements

[lelemmen]

Unfortunately, the main developer for the low-level CI routines left our group (a while ago). Even though the high-level APIs for doing CI calculations is clear, the actual low-level implementation instantly became legacy code upon his departure. In order to make this library more robust, we could really benefit from extra documentation, refactors, ... an overal maintenance of the low-level CI codes.

• LinearExpansion::basisTransform
• SpinUnresolvedONVBasis::calculateOneElectronCouplings
• DOCI (seniority zero) 1 & 2-DMs
• FCI 1 & 2-DMs
• RSQOperators in a seniority-zero ONV basis
• RSQOperators in a full, spin-resolved ONV basis
• RSQOperators in a selected ONV basis
• USQOperators in a selected ONV basis
• GSQOperators in a spin-unresolved ONV basis
• Hubbard evaluations in a full, spin-resolved ONV basis

[lelemmen]

We have some other code lying around for selected CI evaluations for restricted operators

SelectedCI.hpp

// // 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 "ONVBasis/SpinResolvedSelectedONVBasis.hpp"
// #include "QCMethod/CI/HamiltonianBuilder/HamiltonianBuilder.hpp"


// namespace GQCP {


// /**
//  *  Typedef for a type of function that handles where to 'put' a calculated value (e.g. matrix or vector).
//  *
//  *  @param I        index or address of an SpinUnresolvedONV
//  *  @param J        index or address of SpinUnresolvedONV that couples with SpinUnresolvedONV I
//  *  @param value    value related to the coupling
//  */
// using PassToMethod = std::function<void(size_t I, size_t J, double value)>;


// /**
//  *  SelectedCI builds a Hamiltonian matrix in a spin-resolved selected ONV basis
//  */
// class SelectedCI: public HamiltonianBuilder {
// private:
//     SpinResolvedSelectedONVBasis onv_basis;  // the spin-resolved selected ONV basis


// public:
//     // CONSTRUCTORS

//     /**
//      *  @param onv_basis               the spin-resolved selected ONV basis
//      */
//     explicit SelectedCI(const SpinResolvedSelectedONVBasis& onv_basis);


//     // DESTRUCTOR

//     /**
//      *  The default destructor.
//      */
//     ~SelectedCI() = default;


//     // PUBLIC OVERRIDDEN METHODS

//     /**
//      *  @param sq_hamiltonian               the Hamiltonian expressed in an orthonormal basis
//      *
//      *  @return the diagonal of the matrix representation of the SelectedCI Hamiltonian
//      */
//     VectorX<double> calculateDiagonal(const RSQHamiltonian<double>& sq_hamiltonian) const override;

//     /**
//      *  @param sq_hamiltonian               the Hamiltonian expressed in an orthonormal basis
//      *
//      *  @return the SelectedCI Hamiltonian matrix
//      */
//     SquareMatrix<double> constructHamiltonian(const RSQHamiltonian<double>& sq_hamiltonian) const override;

//     /**
//      *  @param sq_hamiltonian               the Hamiltonian expressed in an orthonormal basis
//      *  @param x                            the vector upon which the SelectedCI Hamiltonian acts
//      *  @param diagonal                     the diagonal of the SelectedCI Hamiltonian matrix
//      *
//      *  @return the action of the SelectedCI Hamiltonian on the coefficient vector
//      */
//     VectorX<double> matrixVectorProduct(const RSQHamiltonian<double>& sq_hamiltonian, const VectorX<double>& x, const VectorX<double>& diagonal) const override;

//     /**
//      *  @return the ONV basis that is associated with this HamiltonianBuilder
//      */
//     const BaseONVBasis* onvBasis() const override { return &onv_basis; }


//     // PUBLIC METHODS

//     /**
//      *  Evaluate all Hamiltonian elements, putting the results in the Hamiltonian matrix or matvec through the `method` function
//      *  This function is used both in `constructHamiltonian()` and `matrixVectorProduct()` to avoid duplicate code.
//      *
//      *  @param sq_hamiltonian           the Hamiltonian expressed in an orthonormal basis
//      *  @param method                   the method depending to how you wish to construct the Hamiltonian
//      */
//     void evaluateHamiltonianElements(const RSQHamiltonian<double>& sq_hamiltonian, const PassToMethod& method) const;
// };


// }  // namespace GQCP

SelectedCI.cpp

// // 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/>.

// #include "QCMethod/CI/HamiltonianBuilder/SelectedCI.hpp"


// namespace GQCP {


// /*
//  *  CONSTRUCTORS
//  */

// /**
//  *  @param onv_basis       the full spin-resolved ONV basis
//  */
// SelectedCI::SelectedCI(const SpinResolvedSelectedONVBasis& onv_basis) :
//     HamiltonianBuilder(),
//     onv_basis {onv_basis} {}


// /*
//  *  PUBLIC OVERRIDDEN METHODS
//  */

// /**
//  *  @param sq_hamiltonian               the SelectedCI Hamiltonian parameters in an orthonormal orbital basis
//  *
//  *  @return the diagonal of the matrix representation of the SelectedCI Hamiltonian
//  */
// // VectorX<double> SelectedCI::calculateDiagonal(const RSQHamiltonian<double>& sq_hamiltonian) const {
// //     return this->onv_basis.evaluateOperatorDiagonal(sq_hamiltonian);
// // }


// /**
//  *  @param sq_hamiltonian           the SelectedCI Hamiltonian parameters in an orthonormal orbital basis
//  *
//  *  @return the SelectedCI Hamiltonian matrix
//  */
// SquareMatrix<double> SelectedCI::constructHamiltonian(const RSQHamiltonian<double>& sq_hamiltonian) const {
//     auto K = sq_hamiltonian.core().numberOfOrbitals();
//     if (K != this->onv_basis.numberOfOrbitals()) {
//         throw std::invalid_argument("SelectedCI::constructHamiltonian(RSQHamiltonian<double>): Basis functions of the ONV basis and sq_hamiltonian are incompatible.");
//     }

//     auto dim = onv_basis.dimension();

//     SquareMatrix<double> result_matrix = SquareMatrix<double>::Zero(dim);
//     result_matrix += this->calculateDiagonal(sq_hamiltonian).asDiagonal();

//     // We should put the calculated elements inside the result matrix
//     PassToMethod addToMatrix = [&result_matrix](size_t I, size_t J, double value) { result_matrix(I, J) += value; };

//     this->evaluateHamiltonianElements(sq_hamiltonian, addToMatrix);
//     return result_matrix;
// }


// /**
//  *  @param sq_hamiltonian               the SelectedCI Hamiltonian parameters in an orthonormal orbital basis
//  *  @param x                            the vector upon which the SelectedCI Hamiltonian acts
//  *  @param diagonal                     the diagonal of the SelectedCI Hamiltonian matrix
//  *
//  *  @return the action of the SelectedCI Hamiltonian on the coefficient vector
//  */
// VectorX<double> SelectedCI::matrixVectorProduct(const RSQHamiltonian<double>& sq_hamiltonian, const VectorX<double>& x, const VectorX<double>& diagonal) const {

//     auto K = sq_hamiltonian.core().numberOfOrbitals();
//     if (K != this->onv_basis.numberOfOrbitals()) {
//         throw std::invalid_argument("SelectedCI::matrixVectorProduct(RSQHamiltonian<double>, VectorX<double>, VectorX<double>): Basis functions of the ONV basis and sq_hamiltonian are incompatible.");
//     }

//     VectorX<double> matvec = diagonal.cwiseProduct(x);

//     // We should pass the calculated elements to the resulting vector and perform the product
//     PassToMethod addToMatvec = [&matvec, &x](size_t I, size_t J, double value) { matvec(I) += value * x(J); };

//     this->evaluateHamiltonianElements(sq_hamiltonian, addToMatvec);

//     return matvec;
// }


// /*
//  *  PUBLIC METHODS
//  */

// /**
//  *  Evaluate all Hamiltonian elements, putting the results in the Hamiltonian matrix or matvec through the `method` function
//  *  This function is used both in `constructHamiltonian()` and `matrixVectorProduct()` to avoid duplicate code.
//  *
//  *  @param sq_hamiltonian           the Hamiltonian expressed in an orthonormal basis
//  *  @param method                   the method depending to how you wish to construct the Hamiltonian
//  */
// void SelectedCI::evaluateHamiltonianElements(const RSQHamiltonian<double>& sq_hamiltonian, const PassToMethod& method) const {

//     const size_t dim = onv_basis.dimension();
//     const size_t K = onv_basis.numberOfOrbitals();

//     const auto& h = sq_hamiltonian.core().parameters();
//     const auto& g = sq_hamiltonian.twoElectron().parameters();

//     for (size_t I = 0; I < dim; I++) {  // loop over all addresses (1)
//         SpinResolvedONV configuration_I = this->onv_basis.onvWithIndex(I);
//         SpinUnresolvedONV alpha_I = configuration_I.onv(Spin::alpha);
//         SpinUnresolvedONV beta_I = configuration_I.onv(Spin::beta);

//         // Calculate the off-diagonal elements, by going over all other ONVs
//         for (size_t J = I + 1; J < dim; J++) {

//             SpinResolvedONV configuration_J = this->onv_basis.onvWithIndex(J);
//             SpinUnresolvedONV alpha_J = configuration_J.onv(Spin::alpha);
//             SpinUnresolvedONV beta_J = configuration_J.onv(Spin::beta);

//             if ((alpha_I.countNumberOfDifferences(alpha_J) == 2) && (beta_I.countNumberOfDifferences(beta_J) == 0)) {

//                 // Find the orbitals that are occupied in one string, and aren't in the other
//                 size_t p = alpha_I.findDifferentOccupations(alpha_J)[0];  // we're sure that there is only 1 element in the std::vector<size_t>
//                 size_t q = alpha_J.findDifferentOccupations(alpha_I)[0];  // we're sure that there is only 1 element in the std::vector<size_t>

//                 // Calculate the total sign
//                 int sign = alpha_I.operatorPhaseFactor(p) * alpha_J.operatorPhaseFactor(q);

//                 double value = h(p, q);

//                 method(I, J, sign * value);
//                 method(J, I, sign * value);

//                 for (size_t r = 0; r < K; r++) {  // r loops over spatial orbitals

//                     if (alpha_I.isOccupied(r) && alpha_J.isOccupied(r)) {  // r must be occupied on the left and on the right
//                         if ((p != r) && (q != r)) {                        // can't create or annihilate the same orbital

//                             double value = 0.5 * (g(p, q, r, r) - g(r, q, p, r) - g(p, r, r, q) + g(r, r, p, q));

//                             method(I, J, sign * value);
//                             method(J, I, sign * value);
//                         }
//                     }

//                     if (beta_I.isOccupied(r)) {  // beta_I == beta_J from the previous if-branch

//                         double value = 0.5 * (g(p, q, r, r) + g(r, r, p, q));

//                         method(I, J, sign * value);
//                         method(J, I, sign * value);
//                     }
//                 }
//             }

//             // 0 electron excitations in alpha, 1 in beta
//             if ((alpha_I.countNumberOfDifferences(alpha_J) == 0) && (beta_I.countNumberOfDifferences(beta_J) == 2)) {


//                 // Find the orbitals that are occupied in one string, and aren't in the other
//                 size_t p = beta_I.findDifferentOccupations(beta_J)[0];  // we're sure that there is only 1 element in the std::vector<size_t>
//                 size_t q = beta_J.findDifferentOccupations(beta_I)[0];  // we're sure that there is only 1 element in the std::vector<size_t>

//                 // Calculate the total sign
//                 int sign = beta_I.operatorPhaseFactor(p) * beta_J.operatorPhaseFactor(q);

//                 double value = h(p, q);

//                 method(I, J, sign * value);
//                 method(J, I, sign * value);

//                 for (size_t r = 0; r < K; r++) {  // r loops over spatial orbitals

//                     if (beta_I.isOccupied(r) && beta_J.isOccupied(r)) {  // r must be occupied on the left and on the right
//                         if ((p != r) && (q != r)) {                      // can't create or annihilate the same orbital
//                             double value = 0.5 * (g(p, q, r, r) - g(r, q, p, r) - g(p, r, r, q) + g(r, r, p, q));

//                             method(I, J, sign * value);
//                             method(J, I, sign * value);
//                         }
//                     }

//                     if (alpha_I.isOccupied(r)) {  // alpha_I == alpha_J from the previous if-branch

//                         double value = 0.5 * (g(p, q, r, r) + g(r, r, p, q));

//                         method(I, J, sign * value);
//                         method(J, I, sign * value);
//                     }
//                 }
//             }

//             // 1 electron excitation in alpha, 1 in beta
//             if ((alpha_I.countNumberOfDifferences(alpha_J) == 2) && (beta_I.countNumberOfDifferences(beta_J) == 2)) {

//                 // Find the orbitals that are occupied in one string, and aren't in the other
//                 size_t p = alpha_I.findDifferentOccupations(alpha_J)[0];  // we're sure that there is only 1 element in the std::vector<size_t>
//                 size_t q = alpha_J.findDifferentOccupations(alpha_I)[0];  // we're sure that there is only 1 element in the std::vector<size_t>

//                 size_t r = beta_I.findDifferentOccupations(beta_J)[0];  // we're sure that there is only 1 element in the std::vector<size_t>
//                 size_t s = beta_J.findDifferentOccupations(beta_I)[0];  // we're sure that there is only 1 element in the std::vector<size_t>

//                 int sign = alpha_I.operatorPhaseFactor(p) * alpha_J.operatorPhaseFactor(q) * beta_I.operatorPhaseFactor(r) * beta_J.operatorPhaseFactor(s);
//                 double value = 0.5 * (g(p, q, r, s) + g(r, s, p, q));

//                 method(I, J, sign * value);
//                 method(J, I, sign * value);
//             }

//             // 2 electron excitations in alpha, 0 in beta
//             if ((alpha_I.countNumberOfDifferences(alpha_J) == 4) && (beta_I.countNumberOfDifferences(beta_J) == 0)) {

//                 // Find the orbitals that are occupied in one string, and aren't in the other
//                 std::vector<size_t> occupied_indices_I = alpha_I.findDifferentOccupations(alpha_J);  // we're sure this has two elements
//                 size_t p = occupied_indices_I[0];
//                 size_t r = occupied_indices_I[1];

//                 std::vector<size_t> occupied_indices_J = alpha_J.findDifferentOccupations(alpha_I);  // we're sure this has two elements
//                 size_t q = occupied_indices_J[0];
//                 size_t s = occupied_indices_J[1];

//                 int sign = alpha_I.operatorPhaseFactor(p) * alpha_I.operatorPhaseFactor(r) * alpha_J.operatorPhaseFactor(q) * alpha_J.operatorPhaseFactor(s);

//                 double value = 0.5 * (g(p, q, r, s) - g(p, s, r, q) - g(r, q, p, s) + g(r, s, p, q));

//                 method(I, J, sign * value);
//                 method(J, I, sign * value);
//             }

//             // 0 electron excitations in alpha, 2 in beta
//             if ((alpha_I.countNumberOfDifferences(alpha_J) == 0) && (beta_I.countNumberOfDifferences(beta_J) == 4)) {

//                 // Find the orbitals that are occupied in one string, and aren't in the other
//                 std::vector<size_t> occupied_indices_I = beta_I.findDifferentOccupations(beta_J);  // we're sure this has two elements
//                 size_t p = occupied_indices_I[0];
//                 size_t r = occupied_indices_I[1];

//                 std::vector<size_t> occupied_indices_J = beta_J.findDifferentOccupations(beta_I);  // we're sure this has two elements
//                 size_t q = occupied_indices_J[0];
//                 size_t s = occupied_indices_J[1];

//                 int sign = beta_I.operatorPhaseFactor(p) * beta_I.operatorPhaseFactor(r) * beta_J.operatorPhaseFactor(q) * beta_J.operatorPhaseFactor(s);

//                 double value = 0.5 * (g(p, q, r, s) - g(p, s, r, q) - g(r, q, p, s) + g(r, s, p, q));

//                 method(I, J, sign * value);
//                 method(J, I, sign * value);
//             }
//         }  // loop over addresses J > I
//     }      // loop over addresses I
// }


// }  // namespace GQCP

[lelemmen]

Here are some Hubbard implementations we had

Hubbard.hpp

// // 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 "ONVBasis/SpinResolvedONVBasis.hpp"
// #include "Operator/SecondQuantized/ModelHamiltonian/HubbardHamiltonian.hpp"


// namespace GQCP {


// /**
//  *  Hubbard builds a Hubbard Hamiltonian matrix in the full spin-resolved ONV basis
//  *
//  *  Hubbard distinguishes itself from FCI by explicitly implementing simplified matrix construction and matrix-vector implementations.
//  */
// class Hubbard {
// private:
//     SpinResolvedONVBasis onv_basis;  // the full spin-resolved ONV basis


// public:
//     // CONSTRUCTORS

//     /**
//      *  @param onv_basis       the full spin-resolved ONV basis
//      */
//     explicit Hubbard(const SpinResolvedONVBasis& onv_basis);


//     // DESTRUCTOR
//     ~Hubbard() = default;


//     // PUBLIC METHODS

//     /**
//      *  @param hubbard_hamiltonian              the Hubbard model Hamiltonian
//      *
//      *  @return the diagonal of the matrix representation of the Hubbard model Hamiltonian
//      */
//     VectorX<double> calculateDiagonal(const HubbardHamiltonian<double>& hubbard_hamiltonian) const;

//     /**
//      *  @param hubbard_hamiltonian              the Hubbard model Hamiltonian
//      *
//      *  @return the Hubbard Hamiltonian matrix
//      */
//     SquareMatrix<double> constructHamiltonian(const HubbardHamiltonian<double>& hubbard_hamiltonian) const;

//     /**
//      *  @param hubbard_hamiltonian              the Hubbard model Hamiltonian
//      *  @param x                                the vector upon which the Hubbard Hamiltonian acts
//      *  @param diagonal                         the diagonal of the Hubbard Hamiltonian matrix
//      *
//      *  @return the action of the Hubbard model Hamiltonian on the coefficient vector
//      */
//     VectorX<double> matrixVectorProduct(const HubbardHamiltonian<double>& hubbard_hamiltonian, const VectorX<double>& x, const VectorX<double>& diagonal) const;
// };


// }  // namespace GQCP

Hubbard.cpp

// // 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/>.

// #include "QCMethod/CI/HamiltonianBuilder/Hubbard.hpp"


// namespace GQCP {


// /*
//  *  CONSTRUCTORS
//  */

// /**
//  *  @param onv_basis       the full spin-resolved ONV basis
//  */
// Hubbard::Hubbard(const SpinResolvedONVBasis& onv_basis) :
//     onv_basis {onv_basis} {}


// /*
//  *  PUBLIC METHODS
//  */

// /**
//  *  @param hubbard_hamiltonian              the Hubbard model Hamiltonian
//  *
//  *  @return the diagonal of the matrix representation of the Hubbard model Hamiltonian
//  */
// VectorX<double> Hubbard::calculateDiagonal(const HubbardHamiltonian<double>& hubbard_hamiltonian) const {

//     const auto K = hubbard_hamiltonian.numberOfLatticeSites();
//     if (K != this->onv_basis.numberOfOrbitals()) {
//         throw std::invalid_argument("Hubbard::constructHamiltonian(const HubbardHamiltonian<double>&): The number of spatial orbitals of this ONV basis and the number of lattice sites for the Hubbard Hamiltonian are incompatible.");
//     }


//     // Set up ONV bases for alpha and beta
//     const SpinUnresolvedONVBasis onv_basis_alpha = onv_basis.onvBasisAlpha();
//     const auto dim_alpha = onv_basis_alpha.dimension();

//     SpinUnresolvedONVBasis onv_basis_beta = onv_basis.onvBasisBeta();
//     const auto dim_beta = onv_basis_beta.dimension();

//     const auto dim = onv_basis.dimension();


//     // Calculate the diagonal contributions resulting from the two-electron on-site interactions by iterating over all ONVs.
//     VectorX<double> diagonal = VectorX<double>::Zero(dim);

//     const auto& H = hubbard_hamiltonian.hoppingMatrix();


//     SpinUnresolvedONV onv_alpha = onv_basis_alpha.constructONVFromAddress(0);
//     SpinUnresolvedONV onv_beta = onv_basis_beta.constructONVFromAddress(0);
//     for (size_t Ia = 0; Ia < dim_alpha; Ia++) {  // Ia loops over the addresses of alpha ONVs
//         onv_basis_beta.transformONVCorrespondingToAddress(onv_beta, 0);

//         for (size_t Ib = 0; Ib < dim_beta; Ib++) {      // Ib loops over addresses of beta ONVs
//             const size_t address = Ia * dim_beta + Ib;  // compound address in the spin-resolved ONV basis

//             // There is a contribution for all orbital indices p that are occupied both in the alpha- and beta ONV.
//             std::vector<size_t> occupations = onv_alpha.findMatchingOccupations(onv_beta);
//             for (const auto& p : occupations) {
//                 diagonal(address) += H(p, p);  // the two-electron (on-site repulsion) contributions are on the diagonal of the hopping matrix
//             }

//             if (Ib < dim_beta - 1) {  // prevent the last permutation from occurring
//                 onv_basis_beta.transformONVToNextPermutation(onv_beta);
//             }
//         }  // beta address (Ib) loop

//         if (Ia < dim_alpha - 1) {  // prevent the last permutation from occurring
//             onv_basis_alpha.transformONVToNextPermutation(onv_alpha);
//         }
//     }  // alpha address (Ia) loop
//     return diagonal;
// }


// /**
//  *  @param hubbard_hamiltonian              the Hubbard model Hamiltonian
//  *
//  *  @return the Hubbard Hamiltonian matrix
//  */
// // SquareMatrix<double> Hubbard::constructHamiltonian(const HubbardHamiltonian<double>& hubbard_hamiltonian) const {

// //     const auto K = hubbard_hamiltonian.numberOfLatticeSites();
// //     if (K != this->onv_basis.numberOfOrbitals()) {
// //         throw std::invalid_argument("Hubbard::constructHamiltonian(const HubbardHamiltonian<double>&): The number of spatial orbitals of this ONV basis and the number of lattice sites for the Hubbard Hamiltonian are incompatible.");
// //     }


// //     return this->onv_basis.evaluateOperatorDense(hubbard_hamiltonian.core(), false) + SquareMatrix<double>(this->calculateDiagonal(hubbard_hamiltonian).asDiagonal());
// // }


// // /**
// //  *  @param hubbard_hamiltonian              the Hubbard model Hamiltonian
// //  *  @param x                                the vector upon which the Hubbard Hamiltonian acts
// //  *  @param diagonal                         the diagonal of the Hubbard Hamiltonian matrix
// //  *
// //  *  @return the action of the Hubbard model Hamiltonian on the coefficient vector
// //  */
// // VectorX<double> Hubbard::matrixVectorProduct(const HubbardHamiltonian<double>& hubbard_hamiltonian, const VectorX<double>& x, const VectorX<double>& diagonal) const {

// //     const auto K = hubbard_hamiltonian.numberOfLatticeSites();
// //     if (K != this->onv_basis.numberOfOrbitals()) {
// //         throw std::invalid_argument("Hubbard::constructHamiltonian(const HubbardHamiltonian<double>&): The number of spatial orbitals of this ONV basis and the number of lattice sites for the Hubbard Hamiltonian are incompatible.");
// //     }


// //     // Set up ONV bases for alpha and beta
// //     const SpinUnresolvedONVBasis onv_basis_alpha = onv_basis.onvBasisAlpha();
// //     const auto dim_alpha = onv_basis_alpha.dimension();

// //     const SpinUnresolvedONVBasis onv_basis_beta = onv_basis.onvBasisBeta();
// //     const auto dim_beta = onv_basis_beta.dimension();


// //     // Calculate the Hubbard matrix-vector product, which is the sum of the alpha- and beta one-electron contributions plus the diagonal (which contains the two-electron contributions).
// //     VectorX<double> matvec = diagonal.cwiseProduct(x);

// //     Eigen::Map<Eigen::MatrixXd> matvecmap {matvec.data(), static_cast<long>(dim_beta), static_cast<long>(dim_alpha)};
// //     Eigen::Map<const Eigen::MatrixXd> xmap {x.data(), static_cast<long>(dim_beta), static_cast<long>(dim_alpha)};

// //     const auto H_alpha = onv_basis_alpha.evaluateOperatorSparse(hubbard_hamiltonian.core(), false);  // false: no diagonal contributions
// //     const auto H_beta = onv_basis_beta.evaluateOperatorSparse(hubbard_hamiltonian.core(), false);    // false: no diagonal contributions

// //     matvecmap += xmap * H_alpha + H_beta * xmap;
// //     return matvec;
// // }


// }  // namespace GQCP

Test case

// // 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/>.

// #define BOOST_TEST_MODULE "Hubbard"

// #include <boost/test/unit_test.hpp>


// /**
//  *  Check if a correct constructor does not throw.
//  */
// BOOST_AUTO_TEST_CASE(Hubbard_constructor_throws) {

//     const GQCP::SpinResolvedONVBasis onv_basis {15, 3, 3};
//     BOOST_CHECK_NO_THROW(GQCP::Hubbard hubbard {onv_basis});
// }


// /**
//  *  Check if correct arguments for the Hubbard HamiltonianBuilder API throw (and don't) as expected.
//  */
// BOOST_AUTO_TEST_CASE(Hubbard_public_methods) {

//     // Create a Hubbard model Hamiltonian.
//     const size_t K = 7;  // number of lattice sites
//     const auto H = GQCP::HoppingMatrix<double>::Random(K);
//     const GQCP::HubbardHamiltonian<double> hubbard_hamiltonian {H};


//     // Create a compatible ONV basis.
//     const GQCP::SpinResolvedONVBasis onv_basis_valid {K, 3, 3};
//     const auto dim = onv_basis_valid.dimension();
//     const GQCP::VectorX<double> x = GQCP::VectorX<double>::Random(dim);

//     const GQCP::Hubbard hubbard_builder_valid(onv_basis_valid);
//     const auto diagonal = hubbard_builder_valid.calculateDiagonal(hubbard_hamiltonian);


//     BOOST_CHECK_NO_THROW(hubbard_builder_valid.constructHamiltonian(hubbard_hamiltonian));
//     BOOST_CHECK_NO_THROW(hubbard_builder_valid.matrixVectorProduct(hubbard_hamiltonian, x, diagonal));


//     // Create an incompatible ONV basis.
//     const GQCP::SpinResolvedONVBasis onv_basis_invalid {K + 1, 3, 3};

//     const GQCP::Hubbard hubbard_builder_invalid(onv_basis_invalid);
//     BOOST_CHECK_THROW(hubbard_builder_invalid.constructHamiltonian(hubbard_hamiltonian), std::invalid_argument);
//     BOOST_CHECK_THROW(hubbard_builder_invalid.matrixVectorProduct(hubbard_hamiltonian, x, diagonal), std::invalid_argument);
// }


// /**
//  *  Check if the specialized Hubbard Hamiltonian matrix construction algorithm produces the same result as the unspecialized algorithm.
//  */
// BOOST_AUTO_TEST_CASE(Hubbard_matrix_specialized_vs_unspecialized) {

//     // Create the Hubbard model Hamiltonian and an appropriate ONV basis.
//     const auto K = 4;      // number of lattice sites
//     const size_t N_P = 2;  // number of electron pairs

//     const auto H = GQCP::HoppingMatrix<double>::Random(K);
//     const GQCP::HubbardHamiltonian<double> hubbard_hamiltonian {H};

//     GQCP::SpinResolvedONVBasis onv_basis {K, N_P, N_P};


//     // Create an identical, but 'unspecialized' second-quantized Hamiltonian.
//     const auto sq_hamiltonian = GQCP::RSQHamiltonian<double>::FromHubbard(hubbard_hamiltonian);


//     // Check if the matrix construction produces identical results.
//     const GQCP::Hubbard specialized_builder {onv_basis};
//     const auto H_specialized = specialized_builder.constructHamiltonian(hubbard_hamiltonian);

//     const GQCP::FCI unspecialized_builder {onv_basis};
//     const auto H_unspecialized = unspecialized_builder.constructHamiltonian(sq_hamiltonian);

//     BOOST_CHECK(H_specialized.isApprox(H_unspecialized));
// }


// /**
//  *  Check if the specialized Hubbard Hamiltonian matrix construction algorithm produces the same result as the unspecialized algorithm, for a larger number of lattice sites.
//  */
// BOOST_AUTO_TEST_CASE(Hubbard_matrix_specialized_vs_unspecialized_large) {

//     // Create the Hubbard model Hamiltonian and an appropriate ONV basis.
//     const auto K = 6;      // number of lattice sites
//     const size_t N_P = 3;  // number of electron pairs

//     const auto H = GQCP::HoppingMatrix<double>::Random(K);
//     const GQCP::HubbardHamiltonian<double> hubbard_hamiltonian {H};

//     GQCP::SpinResolvedONVBasis onv_basis {K, N_P, N_P};


//     // Create an identical, but 'unspecialized' second-quantized Hamiltonian.
//     const auto sq_hamiltonian = GQCP::RSQHamiltonian<double>::FromHubbard(hubbard_hamiltonian);


//     // Check if the matrix construction produces identical results.
//     const GQCP::Hubbard specialized_builder {onv_basis};
//     const auto H_specialized = specialized_builder.constructHamiltonian(hubbard_hamiltonian);

//     const GQCP::FCI unspecialized_builder {onv_basis};
//     const auto H_unspecialized = unspecialized_builder.constructHamiltonian(sq_hamiltonian);

//     BOOST_CHECK(H_specialized.isApprox(H_unspecialized));
// }


// /**
//  *  Check if the specialized Hubbard matrix-vector product algorithm produces the same result as the unspecialized algorithm.
//  */
// BOOST_AUTO_TEST_CASE(Hubbard_matvec_specialized_vs_unspecialized) {

//     // Create the Hubbard model Hamiltonian and an appropriate ONV basis.
//     const auto K = 4;      // number of lattice sites
//     const size_t N_P = 2;  // number of electron pairs

//     const auto H = GQCP::HoppingMatrix<double>::Random(K);
//     const GQCP::HubbardHamiltonian<double> hubbard_hamiltonian {H};

//     GQCP::SpinResolvedONVBasis onv_basis {K, N_P, N_P};
//     const auto dim = onv_basis.dimension();
//     const GQCP::VectorX<double> x = GQCP::VectorX<double>::Random(dim);

//     // Create an identical, but 'unspecialized' second-quantized Hamiltonian.
//     const auto sq_hamiltonian = GQCP::RSQHamiltonian<double>::FromHubbard(hubbard_hamiltonian);


//     // Check if the matrix-vector product produces identical results.
//     const GQCP::Hubbard specialized_builder {onv_basis};
//     const auto specialized_diagonal = specialized_builder.calculateDiagonal(hubbard_hamiltonian);
//     const auto spezialized_matvec = specialized_builder.matrixVectorProduct(hubbard_hamiltonian, x, specialized_diagonal);

//     const GQCP::FCI unspecialized_builder {onv_basis};
//     const auto unspecialized_diagonal = unspecialized_builder.calculateDiagonal(sq_hamiltonian);
//     const auto unspecialized_matvec = unspecialized_builder.matrixVectorProduct(sq_hamiltonian, x, unspecialized_diagonal);

//     BOOST_CHECK(specialized_diagonal.isApprox(unspecialized_diagonal));
//     BOOST_CHECK(spezialized_matvec.isApprox(unspecialized_matvec));
// }


// /**
//  *  Check if the specialized Hubbard matrix-vector product algorithm produces the same result as the unspecialized algorithm, for a larger number of lattice sites.
//  */
// BOOST_AUTO_TEST_CASE(Hubbard_matvec_specialized_vs_unspecialized_large) {

//     // Create the Hubbard model Hamiltonian and an appropriate ONV basis.
//     const auto K = 6;      // number of lattice sites
//     const size_t N_P = 3;  // number of electron pairs

//     const auto H = GQCP::HoppingMatrix<double>::Random(K);
//     const GQCP::HubbardHamiltonian<double> hubbard_hamiltonian {H};

//     GQCP::SpinResolvedONVBasis onv_basis {K, N_P, N_P};
//     const auto dim = onv_basis.dimension();
//     const GQCP::VectorX<double> x = GQCP::VectorX<double>::Random(dim);

//     // Create an identical, but 'unspecialized' second-quantized Hamiltonian.
//     const auto sq_hamiltonian = GQCP::RSQHamiltonian<double>::FromHubbard(hubbard_hamiltonian);


//     // Check if the matrix-vector product produces identical results.
//     const GQCP::Hubbard specialized_builder {onv_basis};
//     const auto specialized_diagonal = specialized_builder.calculateDiagonal(hubbard_hamiltonian);
//     const auto spezialized_matvec = specialized_builder.matrixVectorProduct(hubbard_hamiltonian, x, specialized_diagonal);

//     const GQCP::FCI unspecialized_builder {onv_basis};
//     const auto unspecialized_diagonal = unspecialized_builder.calculateDiagonal(sq_hamiltonian);
//     const auto unspecialized_matvec = unspecialized_builder.matrixVectorProduct(sq_hamiltonian, x, unspecialized_diagonal);

//     BOOST_CHECK(specialized_diagonal.isApprox(unspecialized_diagonal));
//     BOOST_CHECK(spezialized_matvec.isApprox(unspecialized_matvec));
// }