GQCG / GQCP

The Ghent Quantum Chemistry Package for electronic structure calculations
https://gqcg.github.io/GQCP/
GNU Lesser General Public License v3.0
32 stars 10 forks source link

Enable the evaluation of `GSQOperators` in a spin-unresolved frozen-core ONV basis #729

Closed lelemmen closed 3 years ago

lelemmen commented 3 years ago

In recent refactors (#688), we had to temporarily disable some of the CI functionality. This issue tracks the re-enabling of the API to evaluate generalized operators in a spin-unresolved frozen ONV basis.

lelemmen commented 3 years ago

Here's the code that we used to support:

Header 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 "ONVBasis/BaseFrozenCoreONVBasis.hpp"
#include "ONVBasis/BaseONVBasis.hpp"
#include "ONVBasis/ONVManipulator.hpp"
#include "ONVBasis/SpinUnresolvedONVBasis.hpp"

namespace GQCP {

/**
 *  A spin-unresolved frozen ONV basis: this is a subset of the N-electron spin-unresolved ONV basis in which the first X orbitals are always occupied.
 */
class SpinUnresolvedFrozenONVBasis: public BaseFrozenCoreONVBasis, public ONVManipulator<SpinUnresolvedFrozenONVBasis> {
protected:
    size_t X;                                 // number of frozen orbitals/electrons
    SpinUnresolvedONVBasis active_onv_basis;  // active (non-frozen) spin-unresolved ONV basis containing only the active electrons (N-X) and orbitals (K-X)

public:
    // CONSTRUCTORS
    /**
     *  @param K        the total number of orbitals
     *  @param N        the total number of electrons
     *  @param X        the number of frozen orbitals and electrons
     */
    SpinUnresolvedFrozenONVBasis(const size_t K, const size_t N, const size_t X);

    /**
     *  @param onv_basis        (to be frozen) full spin-unresolved ONV basis
     *  @param X                the number of frozen orbitals and electrons
     */
    SpinUnresolvedFrozenONVBasis(const SpinUnresolvedONVBasis& onv_basis, const size_t X);

    // DESTRUCTOR

    /**
     *  The default destructor.
     */
    ~SpinUnresolvedFrozenONVBasis() override = default;

    // OVERRIDEN PUBLIC METHODS

    /**
     *  @param representation      a representation of a spin-unresolved ONV
     *
     *  @return the address (i.e. the ordering number) of the given spin-unresolved ONV
     */
    size_t addressOf(const size_t representation) const override;

    /**
     *  @param onv       the spin-unresolved ONV
     *
     *  @return the number of ONVs (with a larger address) this spin-unresolved ONV would couple with given a one electron operator
     */
    size_t countOneElectronCouplings(const SpinUnresolvedONV& onv) const override;

    /**
     *  @return the number of non-zero (non-diagonal) couplings of a one electron coupling scheme in the spin-unresolved ONV basis
     */
    size_t countTotalOneElectronCouplings() const override { return this->active_onv_basis.countTotalOneElectronCouplings(); }

    /**
     *  @return the number of non-zero (non-diagonal) couplings of a two electron coupling scheme in the spin-unresolved ONV basis
     */
    size_t countTotalTwoElectronCouplings() const override { return this->active_onv_basis.countTotalTwoElectronCouplings(); }

    /**
     *  @param onv       the spin-unresolved ONV
     *
     *  @return the number of ONVs (with a larger address) this spin-unresolved ONV would couple with given a two electron operator
     */
    size_t countTwoElectronCouplings(const SpinUnresolvedONV& onv) const override;

    /**
     *  @param representation       a representation of a spin-unresolved ONV
     *
     *  @return the next bitstring permutation in the frozen spin-unresolved ONV basis
     *
     *      Examples (first orbital is frozen):
     *          0101 -> 1001
     *         01101 -> 10011
     */
    size_t nextPermutationOf(const size_t representation) const override;

    /**
      *  Calculate unsigned representation for a given address
      *
      *  @param address                 the address of the representation is calculated
      *
      *  @return unsigned representation of the address
      */
    size_t representationOf(const size_t address) const override;

    // PUBLIC METHODS

    /**
     *  @return this' ONV basis that is considered active
     */
    const SpinUnresolvedONVBasis& activeONVBasis() const { return this->active_onv_basis; }

    using ONVManipulator<SpinUnresolvedFrozenONVBasis>::addressOf;

    /**
     *  @return the number of orbitals that are considered frozen in this ONV basis
     */
    size_t numberOfFrozenOrbitals() const { return this->X; }
};

}  // namespace GQCP

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

#include "ONVBasis/SpinUnresolvedFrozenONVBasis.hpp"

namespace GQCP {

/*
 *  CONSTRUCTORS
 */

/**
 *  @param M        the total number of orbitals
 *  @param N        the total number of electrons
 *  @param X        the number of frozen orbitals and electrons
 */
SpinUnresolvedFrozenONVBasis::SpinUnresolvedFrozenONVBasis(const size_t M, const size_t N, const size_t X) :
    BaseFrozenCoreONVBasis(std::make_shared<SpinUnresolvedONVBasis>(SpinUnresolvedONVBasis(M - X, N - X)), X),
    ONVManipulator(N),
    active_onv_basis {M - X, N - X},
    X {X} {}

/**
 *  @param fock_space       (to be frozen) full spin-unresolved ONV basis
 *  @param X                the number of frozen orbitals and electrons
 */
SpinUnresolvedFrozenONVBasis::SpinUnresolvedFrozenONVBasis(const SpinUnresolvedONVBasis& fock_space, const size_t X) :
    SpinUnresolvedFrozenONVBasis(fock_space.numberOfOrbitals(), fock_space.numberOfElectrons(), X) {}

/*
 *  OVERRIDDEN PUBLIC METHODS
 */

/**
 *  @param representation      a representation of an ONV
 *
 *  @return the address (i.e. the ordering number) of the given ONV
 */
size_t SpinUnresolvedFrozenONVBasis::addressOf(const size_t representation) const {

    // Transform the representation to the sub space, by bitshifting left X number of times to remove the frozen orbital indices
    // Address of the total ONV in this frozen ONV basis is identical to that of the sub ONV in the sub ONV basis.
    return this->active_onv_basis.addressOf(representation >> this->X);
};

/**
 *  @param onv       the ONV
 *
 *  @return the number of ONVs (with a larger address) the given ONV would couple with given a one electron operator
 */
size_t SpinUnresolvedFrozenONVBasis::countOneElectronCouplings(const SpinUnresolvedONV& onv) const {

    size_t V = this->M - N;  // number of virtual orbitals
    size_t coupling_count = 0;

    for (size_t e1 = this->X; e1 < this->N; e1++) {  // start from X to ignore the frozen electrons
        size_t p = onv.occupationIndexOf(e1);
        coupling_count += (V + e1 - p);  // number of virtuals with an index larger than p
    }

    return coupling_count;
}

/**
 *  @param onv       the ONV
 *
 *  @return the number of ONVs (with a larger address) the given ONV would couple with given a two electron operator
 */
size_t SpinUnresolvedFrozenONVBasis::countTwoElectronCouplings(const SpinUnresolvedONV& onv) const {

    size_t V = this->M - N;  // number of virtual orbitals
    size_t coupling_count = 0;

    for (size_t e1 = this->X; e1 < this->N; e1++) {  // start from X to ignore the frozen electrons

        size_t p = onv.occupationIndexOf(e1);
        coupling_count += (V + e1 - p);  //  one electron part

        for (size_t e2 = e1 + 1; e2 < this->N; e2++) {

            size_t q = onv.occupationIndexOf(e2);
            size_t coupling_count2 = (V + e2 - q);
            coupling_count += (V - coupling_count2) * coupling_count2;

            if (coupling_count2 > 1) {
                coupling_count += SpinUnresolvedONVBasis::calculateDimension(coupling_count2, 2);
            }
        }
    }

    return coupling_count;
}

/**
 *  @param representation       a representation of an ONV
 *
 *  @return the next bitstring permutation in this ONV basis
 *
 *      Examples (first orbital is frozen):
 *          0101 -> 1001
 *         01101 -> 10011
 */
size_t SpinUnresolvedFrozenONVBasis::nextPermutationOf(size_t representation) const {

    // Generate the permutation from the active space, bitshift left X number of times to remove the frozen orbital indices before passing it to the active space
    size_t sub_permutation = this->active_onv_basis.nextPermutationOf(representation >> this->X);

    // Transform the permutation to the frozen core space, by bitshifting right X number of times and filling the new 0 bits with 1's (by adding 2^X - 1)
    sub_permutation <<= X;
    sub_permutation += static_cast<size_t>(pow(2, X) - 1);
    return sub_permutation;
};

/**
  *  Calculate unsigned representation for a given address
  *
  *  @param address                 the address of the representation is calculated
  *
  *  @return unsigned representation of the address
  */
size_t SpinUnresolvedFrozenONVBasis::representationOf(size_t address) const {
    // generate the representation in the active space
    size_t representation = this->active_onv_basis.representationOf(address);

    // transform the permutation to the frozen core space, by bitshifting right X number of times and filling the new 0 bits with 1's (by adding 2^X - 1)
    representation <<= X;
    representation += static_cast<size_t>(pow(2, X) - 1);
    return representation;
};

}  // namespace GQCP
lelemmen commented 3 years ago

BaseFrozenCoreONVBasis header

// 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/RSpinOrbitalBasis.hpp"
#include "ONVBasis/BaseONVBasis.hpp"

namespace GQCP {

/**
 *  A small struct to accommodate the return type of "freezeOperator(const ScalarSQTwoElectronOperator<double>&)", as it should return both a one- and two-electron operator.
 */
struct FrozenOperators {
    ScalarSQOneElectronOperator<double> one_op;
    ScalarSQTwoElectronOperator<double> two_op;
};

/**
 *  A base class for a "frozen" ONV basis: this represents an ONV basis in which the first X occupation numbers are always 1
 */
class BaseFrozenCoreONVBasis: public BaseONVBasis {
protected:
    size_t X;                                         // the number of frozen orbitals/electrons
    std::shared_ptr<BaseONVBasis> active_fock_space;  // active (non-frozen) ONV basis containing only the active electrons (N-X) and orbitals (K-X)

public:
    // CONSTRUCTORS

    /**
     *  @param fock_space                    shared pointer to active (non-frozen core) ONV basis
     *  @param X                             the number of frozen orbitals
     */
    BaseFrozenCoreONVBasis(const std::shared_ptr<BaseONVBasis> fock_space, const size_t X);

    // STATIC PUBLIC METHODS

    /**
     *  @param one_op       the one-electron operator in an orthonormal orbital basis
     *  @param X            the number of frozen orbitals
     *
     *  @return 'frozen' one-electron operator which cover evaluations from the active and inactive orbitals
     */
    static ScalarSQOneElectronOperator<double> freezeOperator(const ScalarSQOneElectronOperator<double>& one_op, const size_t X);

    /**
     *  @param two_op       the two-electron operator in an orthonormal orbital basis
     *  @param X            the number of frozen orbitals
     *
     *  @return 'frozen' two-electron operators as a struct of a one- and two-electron operator which cover evaluations from the active and inactive orbitals
     */
    static FrozenOperators freezeOperator(const ScalarSQTwoElectronOperator<double>& two_op, const size_t X);

    /**
     *  @param sq_hamiltonian       the Hamiltonian expressed in an orthonormal basis
     *  @param X                    the number of frozen orbitals
     *
     *  @return a 'frozen' Hamiltonian which cover two-electron integral evaluations from the active and inactive orbitals
     *  (see https://drive.google.com/file/d/1Fnhv2XyNO9Xw9YDoJOXU21_6_x2llntI/view?usp=sharing)
     */
    static RSQHamiltonian<double> freezeOperator(const RSQHamiltonian<double>& sq_hamiltonian, const size_t X);

    /**
     *  @param usq_hamiltonian      the Hamiltonian expressed in an unrestricted orthonormal basis
     *  @param X                    the number of frozen orbitals
     *
     *  @return a 'frozen' Hamiltonian which cover two-electron integral evaluations from the active and inactive orbitals
     *  (see https://drive.google.com/file/d/1Fnhv2XyNO9Xw9YDoJOXU21_6_x2llntI/view?usp=sharing)
     */
    static USQHamiltonian<double> freezeOperator(const USQHamiltonian<double>& usq_hamiltonian, const size_t X);

    /**
     *  @param one_op       the one-electron operator in an orthonormal orbital basis
     *  @param X            the number of frozen orbitals
     *
     *  @return the operator diagonal from strictly evaluating the frozen orbitals in the ONV basis
     */
    static VectorX<double> frozenCoreDiagonal(const ScalarSQOneElectronOperator<double>& one_op, const size_t X, const size_t dimension);

    /**
     *  @param two_op       the two-electron operator in an orthonormal orbital basis
     *  @param X            the number of frozen orbitals
     *
     *  @return the operator diagonal from strictly evaluating the frozen orbitals in the ONV basis
     */
    static VectorX<double> frozenCoreDiagonal(const ScalarSQTwoElectronOperator<double>& two_op, const size_t X, const size_t dimension);

    /**
     *  @param sq_hamiltonian       the Hamiltonian expressed in an orthonormal basis
     *  @param X                    the number of frozen orbitals
     *  @param dimension            the dimension of the diagonal
     *
     *  @return the Hamiltonian diagonal from strictly evaluating the frozen orbitals in the ONV basis
     */
    static VectorX<double> frozenCoreDiagonal(const RSQHamiltonian<double>& sq_hamiltonian, const size_t X, const size_t dimension);

    /**
     *  @param usq_hamiltonian      the Hamiltonian expressed in an unrestricted orthonormal basis
     *  @param X                    the number of frozen orbitals
     *  @param dimension            the dimension of the diagonal
     *
     *  @return the Hamiltonian diagonal from strictly evaluating the frozen orbitals in the ONV basis
     */
    static VectorX<double> frozenCoreDiagonal(const USQHamiltonian<double>& usq_hamiltonian, const size_t X, const size_t dimension);

    // OVERRIDDEN PUBLIC METHODS

    /**
     *  Evaluate the operator in a dense matrix
     *
     *  @param one_op               the one-electron operator in an orthonormal orbital basis to be evaluated in the ONV basis
     *  @param diagonal_values      bool to indicate if diagonal values will be calculated
     *
     *  @return the operator's evaluation in a dense matrix with the dimensions of the ONV basis
     */
    SquareMatrix<double> evaluateOperatorDense(const ScalarSQOneElectronOperator<double>& one_op, const bool diagonal_values) const override;

    /**
     *  Evaluate the operator in a dense matrix
     *
     *  @param two_op               the two-electron operator in an orthonormal orbital basis to be evaluated in the ONV basis
     *  @param diagonal_values      bool to indicate if diagonal values will be calculated
     *
     *  @return the operator's evaluation in a dense matrix with the dimensions of the ONV basis
     */
    SquareMatrix<double> evaluateOperatorDense(const ScalarSQTwoElectronOperator<double>& two_op, const bool diagonal_values) const override;

    /**
     *  Evaluate the Hamiltonian in a dense matrix
     *
     *  @param sq_hamiltonian           the Hamiltonian expressed in an orthonormal basis
     *  @param diagonal_values          bool to indicate if diagonal values will be calculated
     *
     *  @return the Hamiltonian's evaluation in a dense matrix with the dimensions of the ONV basis
     */
    SquareMatrix<double> evaluateOperatorDense(const RSQHamiltonian<double>& sq_hamiltonian, const bool diagonal_values) const override;

    /**
     *  Evaluate the diagonal of the operator
     *
     *  @param one_op               the one-electron operator in an orthonormal orbital basis to be evaluated in the ONV basis
     *
     *  @return the operator's diagonal evaluation in a vector with the dimension of the ONV basis
     */
    VectorX<double> evaluateOperatorDiagonal(const ScalarSQOneElectronOperator<double>& one_op) const override;

    /**
     *  Evaluate the diagonal of the operator
     *
     *  @param two_op               the two-electron operator in an orthonormal orbital basis to be evaluated in the ONV basis
     *
     *  @return the operator's diagonal evaluation in a vector with the dimension of the ONV basis
     */
    VectorX<double> evaluateOperatorDiagonal(const ScalarSQTwoElectronOperator<double>& two_op) const override;

    /**
     *  Evaluate the diagonal of the Hamiltonian
     *
     *  @param sq_hamiltonian              the Hamiltonian expressed in an orthonormal basis
     *
     *  @return the Hamiltonian's diagonal evaluation in a vector with the dimension of the ONV basis
     */
    VectorX<double> evaluateOperatorDiagonal(const RSQHamiltonian<double>& sq_hamiltonian) const override;

    /**
     *  Evaluate the operator in a sparse matrix
     *
     *  @param one_op               the one-electron operator in an orthonormal orbital basis to be evaluated in the ONV basis
     *  @param diagonal_values      bool to indicate if diagonal values will be calculated
     *
     *  @return the operator's evaluation in a sparse matrix with the dimensions of the ONV basis
     */
    Eigen::SparseMatrix<double> evaluateOperatorSparse(const ScalarSQOneElectronOperator<double>& one_op, const bool diagonal_values) const override;

    /**
     *  Evaluate the operator in a sparse matrix
     *
     *  @param two_op               the two-electron operator in an orthonormal orbital basis to be evaluated in the ONV basis
     *  @param diagonal_values      bool to indicate if diagonal values will be calculated
     *
     *  @return the operator's evaluation in a sparse matrix with the dimensions of the ONV basis
     */
    Eigen::SparseMatrix<double> evaluateOperatorSparse(const ScalarSQTwoElectronOperator<double>& two_op, const bool diagonal_values) const override;

    /**
     *  Evaluate the Hamiltonian in a sparse matrix
     *
     *  @param sq_hamiltonian           the Hamiltonian expressed in an orthonormal basis
     *  @param diagonal_values          bool to indicate if diagonal values will be calculated
     *
     *  @return the Hamiltonian's evaluation in a sparse matrix with the dimensions of the ONV basis
     */
    Eigen::SparseMatrix<double> evaluateOperatorSparse(const SQHamiltonian<double>& sq_hamiltonian, const bool diagonal_values) const override;
};

}  // namespace GQCP

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

#include "ONVBasis/BaseFrozenCoreONVBasis.hpp"

namespace GQCP {

/*
 *  CONSTRUCTORS
 */

/**
 *  @param onv_basis                     shared pointer to active (non-frozen core) frozen ONV basis
 *  @param X                             the number of frozen orbitals
 */
BaseFrozenCoreONVBasis::BaseFrozenCoreONVBasis(const std::shared_ptr<GQCP::BaseONVBasis> onv_basis, const size_t X) :
    BaseONVBasis(onv_basis->numberOfOrbitals() + X, onv_basis->dimension()),
    active_fock_space {std::move(onv_basis)},
    X {X} {}

/*
 *  STATIC PUBLIC METHODS
 */

/**
 *  @param one_op       the one-electron operator in an orthonormal orbital basis
 *  @param X            the number of frozen orbitals
 *
 *  @return 'frozen' one-electron operator which cover evaluations from the active and inactive orbitals
 */
ScalarSQOneElectronOperator<double> BaseFrozenCoreONVBasis::freezeOperator(const ScalarSQOneElectronOperator<double>& one_op, const size_t X) {

    size_t K_active = one_op.numberOfOrbitals() - X;
    return ScalarSQOneElectronOperator<double> {one_op.parameters().block(X, X, K_active, K_active)};
}

/**
 *  @param two_op       the two-electron operator in an orthonormal orbital basis
 *  @param X            the number of frozen orbitals
 *
 *  @return 'frozen' two-electron operators as a struct of a one- and two-electron operator which cover evaluations from the active and inactive orbitals
 */
FrozenOperators BaseFrozenCoreONVBasis::freezeOperator(const ScalarSQTwoElectronOperator<double>& two_op, const size_t X) {

    size_t K_active = two_op.numberOfOrbitals() - X;
    SquareMatrix<double> frozen_one_op_par = SquareMatrix<double>::Zero(K_active);
    const auto& two_op_par = two_op.parameters();
    const auto frozen_two_op_par = SquareRankFourTensor<double>::FromBlock(two_op_par, X, X, X, X);

    // Frozen two-electron integrals can be rewritten partially as one electron integrals.
    for (size_t i = 0; i < K_active; i++) {  // iterate over the active orbitals
        size_t q = i + X;                    // map active orbitals indexes to total orbital indexes (those including the frozen orbitals)

        for (size_t l = 0; l < X; l++) {  // iterate over the frozen orbitals
            frozen_one_op_par(i, i) += two_op_par(q, q, l, l) + two_op_par(l, l, q, q) - two_op_par(q, l, l, q) / 2 - two_op_par(l, q, q, l) / 2;
        }

        for (size_t j = i + 1; j < K_active; j++) {  // iterate over the active orbitals
            size_t p = j + X;                        // map active orbitals indexes to total orbital indexes (those including the frozen orbitals)

            for (size_t l = 0; l < X; l++) {  // iterate over the frozen orbitals
                frozen_one_op_par(i, j) += two_op_par(q, p, l, l) + two_op_par(l, l, q, p) - two_op_par(q, l, l, p) / 2 - two_op_par(l, p, q, l) / 2;
                frozen_one_op_par(j, i) += two_op_par(p, q, l, l) + two_op_par(l, l, p, q) - two_op_par(p, l, l, q) / 2 - two_op_par(l, q, p, l) / 2;
            }
        }
    }

    return {ScalarSQOneElectronOperator<double> {frozen_one_op_par}, ScalarSQTwoElectronOperator<double> {frozen_two_op_par}};
}

/**
 *  @param sq_hamiltonian       the Hamiltonian expressed in an orthonormal basis
 *  @param spinor_basis         the spinor basis
 *  @param X                    the number of frozen orbitals
 *
 *  @return a 'frozen' Hamiltonian which cover two-electron integral evaluations from the active and inactive orbitals
 *  (see https://drive.google.com/file/d/1Fnhv2XyNO9Xw9YDoJOXU21_6_x2llntI/view?usp=sharing)
 */
RSQHamiltonian<double> BaseFrozenCoreONVBasis::freezeOperator(const RSQHamiltonian<double>& sq_hamiltonian, const size_t X) {

    size_t K_active = sq_hamiltonian.numberOfOrbitals() - X;  // number of non-frozen orbitals

    const auto frozen_components_g = BaseFrozenCoreONVBasis::freezeOperator(sq_hamiltonian.twoElectron(), X);

    ScalarSQOneElectronOperator<double> h = BaseFrozenCoreONVBasis::freezeOperator(sq_hamiltonian.core(), X) + frozen_components_g.one_op;  // active
    ScalarSQTwoElectronOperator<double> g = frozen_components_g.two_op;

    return RSQHamiltonian<double>(h, g);
}

/**
 *  @param usq_hamiltonian      the Hamiltonian expressed in an unrestricted orthonormal basis
 *  @param X                    the number of frozen orbitals
 *
 *  @return a 'frozen' Hamiltonian which cover two-electron integral evaluations from the active and inactive orbitals
 *  (see https://drive.google.com/file/d/1Fnhv2XyNO9Xw9YDoJOXU21_6_x2llntI/view?usp=sharing)
 */
USQHamiltonian<double> BaseFrozenCoreONVBasis::freezeOperator(const USQHamiltonian<double>& usq_hamiltonian, size_t X) {

    size_t K_active = usq_hamiltonian.numberOfOrbitals() / 2 - X;  // number of non-frozen orbitals

    SquareMatrix<double> frozen_one_op_par_alpha = usq_hamiltonian.spinHamiltonian(Spin::alpha).core().parameters().block(X, X, K_active, K_active);
    SquareMatrix<double> frozen_one_op_par_beta = usq_hamiltonian.spinHamiltonian(Spin::beta).core().parameters().block(X, X, K_active, K_active);

    const auto& two_op_par_alpha = usq_hamiltonian.spinHamiltonian(Spin::alpha).twoElectron().parameters();
    const auto& two_op_par_beta = usq_hamiltonian.spinHamiltonian(Spin::beta).twoElectron().parameters();
    const auto& two_op_par_mixed = usq_hamiltonian.twoElectronMixed().parameters();

    SquareRankFourTensor<double> frozen_two_op_par_alpha = SquareRankFourTensor<double>::FromBlock(usq_hamiltonian.spinHamiltonian(Spin::alpha).twoElectron().parameters(), X, X, X, X);
    const auto frozen_two_op_par_beta = SquareRankFourTensor<double>::FromBlock(usq_hamiltonian.spinHamiltonian(Spin::beta).twoElectron().parameters(), X, X, X, X);
    const auto frozen_two_op_par_mixed = SquareRankFourTensor<double>::FromBlock(usq_hamiltonian.twoElectronMixed().parameters(), X, X, X, X);

    // Frozen two-electron integrals can be rewritten partially as one electron integrals
    for (size_t i = 0; i < K_active; i++) {  // iterate over the active orbitals
        size_t q = i + X;                    // map active orbitals indexes to total orbital indexes (those including the frozen orbitals)

        for (size_t l = 0; l < X; l++) {  // iterate over the frozen orbitals
            frozen_one_op_par_alpha(i, i) += two_op_par_mixed(q, q, l, l) + two_op_par_alpha(l, l, q, q) - two_op_par_alpha(q, l, l, q) / 2 - two_op_par_alpha(l, q, q, l) / 2;
            // There's a small catch here, the variable 'two_op_par_mixed' is represented as g_aabb, the formulas dictate we need g_bbaa for the beta component, hence we switched the indices
            frozen_one_op_par_beta(i, i) += two_op_par_mixed(l, l, q, q) + two_op_par_beta(l, l, q, q) - two_op_par_beta(q, l, l, q) / 2 - two_op_par_beta(l, q, q, l) / 2;
        }

        for (size_t j = i + 1; j < K_active; j++) {  // iterate over the active orbitals
            size_t p = j + X;                        // map active orbitals indexes to total orbital indexes (those including the frozen orbitals)

            for (size_t l = 0; l < X; l++) {  // iterate over the frozen orbitals
                frozen_one_op_par_alpha(i, j) += two_op_par_mixed(q, p, l, l) + two_op_par_alpha(l, l, q, p) - two_op_par_alpha(q, l, l, p) / 2 - two_op_par_alpha(l, p, q, l) / 2;
                frozen_one_op_par_beta(i, j) += two_op_par_mixed(l, l, q, p) + two_op_par_beta(l, l, q, p) - two_op_par_beta(q, l, l, p) / 2 - two_op_par_beta(l, p, q, l) / 2;

                frozen_one_op_par_alpha(j, i) += two_op_par_mixed(p, q, l, l) + two_op_par_alpha(l, l, p, q) - two_op_par_alpha(p, l, l, q) / 2 - two_op_par_alpha(l, q, p, l) / 2;
                frozen_one_op_par_beta(j, i) += two_op_par_mixed(l, l, p, q) + two_op_par_beta(l, l, p, q) - two_op_par_beta(p, l, l, q) / 2 - two_op_par_beta(l, q, p, l) / 2;
            }
        }
    }

    return USQHamiltonian<double>(RSQHamiltonian<double>(ScalarSQOneElectronOperator<double> {frozen_one_op_par_alpha}, ScalarSQTwoElectronOperator<double> {frozen_two_op_par_alpha}),
                                  RSQHamiltonian<double>(ScalarSQOneElectronOperator<double> {frozen_one_op_par_beta}, ScalarSQTwoElectronOperator<double> {frozen_two_op_par_beta}),
                                  ScalarSQTwoElectronOperator<double>(frozen_two_op_par_mixed));
}

/**
 *  @param one_op       the one-electron operator in an orthonormal orbital basis
 *  @param X            the number of frozen orbitals
 *
 *  @return the operator diagonal from strictly evaluating the frozen orbitals in the frozen ONV basis
 */
VectorX<double> BaseFrozenCoreONVBasis::frozenCoreDiagonal(const ScalarSQOneElectronOperator<double>& one_op, const size_t X, const size_t dimension) {

    const auto& one_op_par = one_op.parameters();

    // The diagonal value for the frozen orbitals is the same for each frozen ONV
    double value = 0;
    for (size_t i = 0; i < X; i++) {
        value += 2 * one_op_par(i, i);
    }

    return VectorX<double>::Constant(dimension, value);
}

/**
 *  @param two_op       the two-electron operator in an orthonormal orbital basis
 *  @param X            the number of frozen orbitals
 *
 *  @return the operator diagonal from strictly evaluating the frozen orbitals in the frozen ONV basis
 */
VectorX<double> BaseFrozenCoreONVBasis::frozenCoreDiagonal(const ScalarSQTwoElectronOperator<double>& two_op, const size_t X, const size_t dimension) {

    const auto& two_op_par = two_op.parameters();

    // The diagonal value for the frozen orbitals is the same for each ONV
    double value = 0;
    for (size_t i = 0; i < X; i++) {
        value += two_op_par(i, i, i, i);

        for (size_t j = i + 1; j < X; j++) {
            value += 2 * two_op_par(i, i, j, j) + 2 * two_op_par(j, j, i, i) - two_op_par(j, i, i, j) - two_op_par(i, j, j, i);
        }
    }

    return VectorX<double>::Constant(dimension, value);
}

/**
 *  @param sq_hamiltonian       the Hamiltonian expressed in an orthonormal basis
 *  @param X                    the number of frozen orbitals
 *
 *  @return the Hamiltonian diagonal from strictly evaluating the frozen orbitals in the frozen ONV basis
 */
VectorX<double> BaseFrozenCoreONVBasis::frozenCoreDiagonal(const RSQHamiltonian<double>& sq_hamiltonian, const size_t X, const size_t dimension) {

    return BaseFrozenCoreONVBasis::frozenCoreDiagonal(sq_hamiltonian.core(), X, dimension) + BaseFrozenCoreONVBasis::frozenCoreDiagonal(sq_hamiltonian.twoElectron(), X, dimension);
}

/**
 *  @param usq_hamiltonian      the Hamiltonian expressed in an unrestricted orthonormal basis
 *  @param X                    the number of frozen orbitals
 *  @param dimension            the dimension of the diagonal
 *
 *  @return the Hamiltonian diagonal from strictly evaluating the frozen orbitals in a (any) frozen ONV basis
 */
VectorX<double> BaseFrozenCoreONVBasis::frozenCoreDiagonal(const USQHamiltonian<double>& usq_hamiltonian, const size_t X, const size_t dimension) {

    SquareMatrix<double> one_op_par_alpha = usq_hamiltonian.spinHamiltonian(Spin::alpha).core().parameters();
    SquareMatrix<double> one_op_par_beta = usq_hamiltonian.spinHamiltonian(Spin::beta).core().parameters();

    const auto& two_op_par_mixed = usq_hamiltonian.twoElectronMixed().parameters();
    const auto& two_op_par_a = usq_hamiltonian.spinHamiltonian(Spin::alpha).twoElectron().parameters();
    const auto& two_op_par_b = usq_hamiltonian.spinHamiltonian(Spin::beta).twoElectron().parameters();

    // The diagonal value for the frozen orbitals is the same for each ONV
    double value = 0;
    for (size_t i = 0; i < X; i++) {
        value += two_op_par_mixed(i, i, i, i);
        value += one_op_par_alpha(i, i) + one_op_par_beta(i, i);
        for (size_t j = i + 1; j < X; j++) {
            value += two_op_par_mixed(i, i, j, j) + two_op_par_mixed(j, j, i, i) +
                     two_op_par_a(i, i, j, j) / 2 + two_op_par_a(j, j, i, i) / 2 - two_op_par_a(j, i, i, j) / 2 - two_op_par_a(i, j, j, i) / 2 +
                     two_op_par_b(i, i, j, j) / 2 + two_op_par_b(j, j, i, i) / 2 - two_op_par_b(j, i, i, j) / 2 - two_op_par_b(i, j, j, i) / 2;
        }
    }

    return VectorX<double>::Constant(dimension, value);
}

/*
 *  PUBLIC OVERRIDDEN METHODS
 */

/**
 *  Evaluate the operator in a dense matrix
 *
 *  @param one_op               the one-electron operator in an orthonormal orbital basis to be evaluated in the frozen ONV basis
 *  @param diagonal_values      bool to indicate if diagonal values will be calculated
 *
 *  @return the operator's evaluation in a dense matrix with the dimensions of the frozen ONV basis
 */
SquareMatrix<double> BaseFrozenCoreONVBasis::evaluateOperatorDense(const ScalarSQOneElectronOperator<double>& one_op, const bool diagonal_values) const {

    // Freeze the Hamiltonian
    ScalarSQOneElectronOperator<double> frozen_one_op = BaseFrozenCoreONVBasis::freezeOperator(one_op, this->X);

    // Evaluate the frozen operator in the active space
    auto evaluation = this->active_fock_space->evaluateOperatorDense(frozen_one_op, diagonal_values);

    if (diagonal_values) {
        // Diagonal correction
        const auto frozen_core_diagonal = BaseFrozenCoreONVBasis::frozenCoreDiagonal(one_op, this->X, this->active_fock_space->dimension());
        evaluation += frozen_core_diagonal.asDiagonal();
    }

    return evaluation;
}

/**
 *  Evaluate the operator in a dense matrix
 *
 *  @param two_op               the two-electron operator in an orthonormal orbital basis to be evaluated in the frozen ONV basis
 *  @param diagonal_values      bool to indicate if diagonal values will be calculated
 *
 *  @return the operator's evaluation in a dense matrix with the dimensions of the frozen ONV basis
 */
SquareMatrix<double> BaseFrozenCoreONVBasis::evaluateOperatorDense(const ScalarSQTwoElectronOperator<double>& two_op, const bool diagonal_values) const {

    // Freeze the operators
    const auto frozen_ops = BaseFrozenCoreONVBasis::freezeOperator(two_op, this->X);

    // Evaluate the frozen operator in the active space
    auto evaluation = this->active_fock_space->evaluateOperatorDense(frozen_ops.one_op, diagonal_values);
    evaluation += this->active_fock_space->evaluateOperatorDense(frozen_ops.two_op, diagonal_values);

    if (diagonal_values) {
        // Diagonal correction
        const auto frozen_core_diagonal = BaseFrozenCoreONVBasis::frozenCoreDiagonal(two_op, this->X, this->active_fock_space->dimension());
        evaluation += frozen_core_diagonal.asDiagonal();
    }

    return evaluation;
}

/**
 *  Evaluate the Hamiltonian in a dense matrix
 *
 *  @param sq_hamiltonian               the Hamiltonian expressed in an orthonormal basis
 *  @param diagonal_values              bool to indicate if diagonal values will be calculated
 *
 *  @return the Hamiltonian's evaluation in a dense matrix with the dimensions of the frozen ONV basis
 */
SquareMatrix<double> BaseFrozenCoreONVBasis::evaluateOperatorDense(const RSQHamiltonian<double>& sq_hamiltonian, const bool diagonal_values) const {
    // Freeze the operators
    const auto frozen_ham_par = BaseFrozenCoreONVBasis::freezeOperator(sq_hamiltonian, this->X);

    // Evaluate the frozen operator in the active space
    auto evaluation = this->active_fock_space->evaluateOperatorDense(frozen_ham_par, diagonal_values);

    if (diagonal_values) {
        // Diagonal correction
        const auto frozen_core_diagonal = BaseFrozenCoreONVBasis::frozenCoreDiagonal(sq_hamiltonian, this->X, this->active_fock_space->dimension());
        evaluation += frozen_core_diagonal.asDiagonal();
    }

    return evaluation;
}

/**
 *  Evaluate the diagonal of the operator
 *
 *  @param one_op               the one-electron operator in an orthonormal orbital basis to be evaluated in the frozen ONV basis
 *
 *  @return the operator's diagonal evaluation in a vector with the dimension of the frozen ONV basis
 */
VectorX<double> BaseFrozenCoreONVBasis::evaluateOperatorDiagonal(const ScalarSQOneElectronOperator<double>& one_op) const {

    const auto frozen_op = BaseFrozenCoreONVBasis::freezeOperator(one_op, this->X);

    // Calculate diagonal in the active space with the "frozen" Hamiltonian
    const auto diagonal = this->active_fock_space->evaluateOperatorDiagonal(frozen_op);

    // Calculate diagonal for the frozen orbitals
    const auto frozen_core_diagonal = BaseFrozenCoreONVBasis::frozenCoreDiagonal(one_op, this->X, this->active_fock_space->dimension());

    return diagonal + frozen_core_diagonal;
};

/**
 *  Evaluate the diagonal of the operator
 *
 *  @param two_op               the two-electron operator in an orthonormal orbital basis to be evaluated in the frozen ONV basis
 *
 *  @return the operator's diagonal evaluation in a vector with the dimension of the frozen ONV basis
 */
VectorX<double> BaseFrozenCoreONVBasis::evaluateOperatorDiagonal(const ScalarSQTwoElectronOperator<double>& two_op) const {

    const auto frozen_ops = BaseFrozenCoreONVBasis::freezeOperator(two_op, this->X);

    // Calculate diagonal in the active space with the "frozen" Hamiltonian
    auto diagonal = this->active_fock_space->evaluateOperatorDiagonal(frozen_ops.two_op);
    diagonal += this->active_fock_space->evaluateOperatorDiagonal(frozen_ops.one_op);

    // Calculate diagonal for the frozen orbitals
    const auto frozen_core_diagonal = BaseFrozenCoreONVBasis::frozenCoreDiagonal(two_op, this->X, this->active_fock_space->dimension());

    return diagonal + frozen_core_diagonal;
}

/**
 *  Evaluate the diagonal of the Hamiltonian
 *
 *  @param sq_hamiltonian               the Hamiltonian expressed in an orthonormal basis
 *
 *  @return the Hamiltonian's diagonal evaluation in a vector with the dimension of the frozen ONV basis
 */
VectorX<double> BaseFrozenCoreONVBasis::evaluateOperatorDiagonal(const RSQHamiltonian<double>& sq_hamiltonian) const {

    const auto frozen_ham_par = BaseFrozenCoreONVBasis::freezeOperator(sq_hamiltonian, this->X);

    // Calculate diagonal in the active space with the "frozen" Hamiltonian
    const auto diagonal = this->active_fock_space->evaluateOperatorDiagonal(frozen_ham_par);

    // Calculate diagonal for the frozen orbitals
    const auto frozen_core_diagonal = BaseFrozenCoreONVBasis::frozenCoreDiagonal(sq_hamiltonian, this->X, this->active_fock_space->dimension());

    return diagonal + frozen_core_diagonal;
}

/**
 *  Evaluate the operator in a sparse matrix
 *
 *  @param one_op               the one-electron operator in an orthonormal orbital basis to be evaluated in the frozen ONV basis
 *  @param diagonal_values      bool to indicate if diagonal values will be calculated
 *
 *  @return the operator's evaluation in a sparse matrix with the dimensions of the frozen ONV basis
 */
Eigen::SparseMatrix<double> BaseFrozenCoreONVBasis::evaluateOperatorSparse(const ScalarSQOneElectronOperator<double>& one_op, const bool diagonal_values) const {

    // Freeze the operator
    ScalarSQOneElectronOperator<double> frozen_one_op = BaseFrozenCoreONVBasis::freezeOperator(one_op, this->X);

    // Evaluate the frozen operator in the active space
    auto evaluation = this->active_fock_space->evaluateOperatorSparse(frozen_one_op, diagonal_values);

    if (diagonal_values) {
        // Diagonal correction
        const auto frozen_core_diagonal = BaseFrozenCoreONVBasis::frozenCoreDiagonal(one_op, this->X, this->active_fock_space->dimension());
        evaluation += frozen_core_diagonal.asDiagonal();
    }

    return evaluation;
}

/**
 *  Evaluate the operator in a sparse matrix
 *
 *  @param two_op               the two-electron operator in an orthonormal orbital basis to be evaluated in the frozen ONV basis
 *  @param diagonal_values      bool to indicate if diagonal values will be calculated
 *
 *  @return the operator's evaluation in a sparse matrix with the dimensions of the frozen ONV basis
 */
Eigen::SparseMatrix<double> BaseFrozenCoreONVBasis::evaluateOperatorSparse(const ScalarSQTwoElectronOperator<double>& two_op, const bool diagonal_values) const {

    // Freeze the operators
    const auto frozen_ops = BaseFrozenCoreONVBasis::freezeOperator(two_op, this->X);

    // Evaluate the frozen operator in the active space
    auto evaluation = this->active_fock_space->evaluateOperatorSparse(frozen_ops.one_op, diagonal_values);
    evaluation += this->active_fock_space->evaluateOperatorSparse(frozen_ops.two_op, diagonal_values);

    if (diagonal_values) {
        // Diagonal correction
        const auto frozen_core_diagonal = BaseFrozenCoreONVBasis::frozenCoreDiagonal(two_op, this->X, this->active_fock_space->dimension());
        evaluation += frozen_core_diagonal.asDiagonal();
    }

    return evaluation;
}

/**
 *  Evaluate the Hamiltonian in a sparse matrix
 *
 *  @param sq_hamiltonian               the Hamiltonian expressed in an orthonormal basis
 *  @param diagonal_values              bool to indicate if diagonal values will be calculated
 *
 *  @return the Hamiltonian's evaluation in a sparse matrix with the dimensions of the frozen ONV basis
 */
Eigen::SparseMatrix<double> BaseFrozenCoreONVBasis::evaluateOperatorSparse(const RSQHamiltonian<double>& sq_hamiltonian, const bool diagonal_values) const {

    // Freeze the operators
    const auto frozen_ham_par = BaseFrozenCoreONVBasis::freezeOperator(sq_hamiltonian, this->X);

    // Evaluate the frozen operator in the active space
    auto evaluation = this->active_fock_space->evaluateOperatorSparse(frozen_ham_par, diagonal_values);

    if (diagonal_values) {
        // Diagonal correction
        const auto frozen_core_diagonal = BaseFrozenCoreONVBasis::frozenCoreDiagonal(sq_hamiltonian, this->X, this->active_fock_space->dimension());
        evaluation += frozen_core_diagonal.asDiagonal();
    }

    return evaluation;
}

}  // namespace GQCP
lelemmen commented 3 years ago

A 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 "SpinUnresolvedFrozenONVBasis"

#include <boost/test/unit_test.hpp>

#include "ONVBasis/SpinResolvedFrozenONVBasis.hpp"
#include "ONVBasis/SpinUnresolvedFrozenONVBasis.hpp"

BOOST_AUTO_TEST_CASE(FrozenONVBasis_constructor) {

    BOOST_CHECK_NO_THROW(GQCP::SpinUnresolvedFrozenONVBasis(10, 5, 1));
    BOOST_CHECK_NO_THROW(GQCP::SpinUnresolvedFrozenONVBasis(10, 5, 2));
    BOOST_CHECK_NO_THROW(GQCP::SpinUnresolvedFrozenONVBasis(10, 5, 3));
    BOOST_CHECK_NO_THROW(GQCP::SpinUnresolvedFrozenONVBasis(10, 5, 4));
    BOOST_CHECK_NO_THROW(GQCP::SpinUnresolvedFrozenONVBasis(10, 5, 5));
}

BOOST_AUTO_TEST_CASE(coupling_count) {

    GQCP::SpinUnresolvedFrozenONVBasis fock_space {7, 5, 2};
    GQCP::SpinUnresolvedONV onv = fock_space.constructONVFromAddress(3);

    // We only count couplings with larger addresses

    BOOST_CHECK(fock_space.countOneElectronCouplings(onv) == 3);
    BOOST_CHECK(fock_space.countTwoElectronCouplings(onv) == 3 + 3);

    onv = fock_space.constructONVFromAddress(0);

    BOOST_CHECK(fock_space.countOneElectronCouplings(onv) == 6);
    BOOST_CHECK(fock_space.countTwoElectronCouplings(onv) == 6 + 3);

    // Test whether the total count matches that of individual counts of all ONVs in the ONV basis
    GQCP::SpinUnresolvedFrozenONVBasis fock_space2 {18, 10, 2};

    size_t coupling_count1 = 0;
    size_t coupling_count2 = 0;
    onv = fock_space2.constructONVFromAddress(0);  // spin string with address 0

    for (size_t I = 0; I < fock_space2.dimension(); I++) {  // I_alpha loops over all addresses of alpha spin strings
        if (I > 0) {
            fock_space2.transformONVToNextPermutation(onv);
        }
        coupling_count1 += fock_space2.countOneElectronCouplings(onv);
        coupling_count2 += fock_space2.countTwoElectronCouplings(onv);
    }

    BOOST_CHECK(2 * coupling_count1 == fock_space2.countTotalOneElectronCouplings());
    BOOST_CHECK(2 * coupling_count2 == fock_space2.countTotalTwoElectronCouplings());
}

BOOST_AUTO_TEST_CASE(ulongNextPermutation_getAddress_calculateRepresentation) {

    size_t K = 5;
    size_t N = 3;
    size_t X = 1;  // number of frozen orbitals/electrons

    GQCP::SpinUnresolvedFrozenONVBasis frozen_fock_space {K, N, X};

    // "00111", "01011", "01101", "10011", "10101", "11001"
    size_t bit_frozen_fock_space[6] = {7, 11, 13, 19, 21, 25};

    size_t bit_onv = bit_frozen_fock_space[0];
    BOOST_CHECK(frozen_fock_space.addressOf(bit_onv) == 0);
    for (size_t i = 1; i < frozen_fock_space.dimension(); i++) {
        bit_onv = frozen_fock_space.nextPermutationOf(bit_onv);
        BOOST_CHECK(bit_onv == bit_frozen_fock_space[i]);
        BOOST_CHECK(frozen_fock_space.addressOf(bit_onv) == i);
        BOOST_CHECK(frozen_fock_space.representationOf(i) == bit_onv);
    }
}

BOOST_AUTO_TEST_CASE(address_setNext_frozen_space) {

    // Here we will test a set of permutations through a frozen ONV basis of K = 15, N = 5, X = 2
    GQCP::SpinUnresolvedFrozenONVBasis fock_space {15, 5, 2};

    // Retrieve the first ONV of the ONV basis
    GQCP::SpinUnresolvedONV onv_test = fock_space.constructONVFromAddress(0);

    const size_t permutations = 285;
    bool is_correct = true;  // variable that is updated to false if an unexpected result occurs

    // Iterate through the ONV basis in reverse lexicographical order and test whether address matches
    for (size_t i = 0; i < permutations; i++) {

        // Tests address
        if (i != fock_space.addressOf(onv_test)) {
            is_correct = false;
        }

        // transforms the given ONV to the next ONV in the ONV basis
        if (i < permutations - 1) {
            fock_space.transformONVToNextPermutation(onv_test);
        }
    }

    // Checks if no unexpected results occurred in a full iteration
    BOOST_CHECK(is_correct);
}

// BOOST_AUTO_TEST_CASE(ONVBasis_EvaluateOperator_diagonal_vs_no_diagonal) {

//     GQCP::Molecule hchain = GQCP::Molecule::HChain(6, 0.742, 2);
//     GQCP::RSpinOrbitalBasis<double, GQCP::GTOShell> spinor_basis {hchain, "STO-3G"};
//     spinor_basis.lowdinOrthonormalize();
//     auto sq_hamiltonian = GQCP::RSQHamiltonian<double>::Molecular(spinor_basis, hchain);  // in the Löwdin basis

//     GQCP::SpinUnresolvedFrozenONVBasis fock_space {6, 4, 2};

//     GQCP::SquareMatrix<double> hamiltonian = fock_space.evaluateOperatorDense(sq_hamiltonian, true);
//     GQCP::SquareMatrix<double> hamiltonian_no_diagonal = fock_space.evaluateOperatorDense(sq_hamiltonian, false);
//     GQCP::VectorX<double> hamiltonian_diagonal = fock_space.evaluateOperatorDiagonal(sq_hamiltonian);

//     // Test if non-diagonal evaluation and diagonal evaluations are correct
//     BOOST_CHECK(hamiltonian.isApprox(hamiltonian_no_diagonal + GQCP::SquareMatrix<double>(hamiltonian_diagonal.asDiagonal())));
// }
lelemmen commented 3 years ago

ONVManipulator header

// 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/SpinUnresolvedONV.hpp"
#include "Utilities/CRTP.hpp"

namespace GQCP {

/**
 *  An interface for ONV bases that generate, permutate and can calculate the address of a spin-unresolved ONV given its unsigned representation.
 * 
 *  @tparam _DerivedManipulator             the type of the derived class
 */
template <typename _DerivedManipulator>
class ONVManipulator:
    public CRTP<_DerivedManipulator> {

protected:
    size_t N;  // number of electrons

public:
    /*
     *  CONSTRUCTORS
     */

    /**
     *  Default constructor setting everything to zero
     */
    ONVManipulator() :
        ONVManipulator(0) {}

    /**
     *  @param N            the number of electrons
     */
    ONVManipulator(const size_t N) :
        N {N} {}

    // DESTRUCTOR
    virtual ~ONVManipulator() = default;

    // VIRTUAL PUBLIC METHODS

    /**
     *  @param representation           a representation of a spin-unresolved ONV
     *
     *  @return the address (i.e. the ordering number) of the given spin-unresolved ONV's representation
     */
    virtual size_t addressOf(const size_t representation) const = 0;

    /**
     *  @param onv              the spin-unresolved ONV
     *
     *  @return the number of ONVs (with a larger address) the given spin-unresolved ONV would couple with when evaluating a one electron operator
     */
    virtual size_t countOneElectronCouplings(const SpinUnresolvedONV& onv) const = 0;

    /**
     *  @return the number non-zero (non-diagonal) couplings of a one electron coupling scheme
     */
    virtual size_t countTotalOneElectronCouplings() const = 0;

    /**
     *  @return the number non-zero (non-diagonal) couplings of a two electron coupling scheme
     */
    virtual size_t countTotalTwoElectronCouplings() const = 0;

    /**
     *  @param onv       the spin-resolved ONV
     *
     *  @return the number of ONVs (with a larger address) the given spin-resolved ONV would couple with when evaluating a two electron operator
     */
    virtual size_t countTwoElectronCouplings(const SpinUnresolvedONV& onv) const = 0;

    /**
      *  Calculate unsigned representation for a given address
      *
      *  @param address                 the address of the representation is calculated
      *
      *  @return unsigned representation of the address
      */
    virtual size_t representationOf(const size_t address) const = 0;

    /**
     *  @param representation           a representation of a spin-unresolved ONV
     *
     *  @return the next bitstring permutation in the spin-unresolved ONV basis
     */
    virtual size_t nextPermutationOf(const size_t representation) const = 0;

    // PUBLIC METHODS

    /**
     *  @param onv      the spin-unresolved ONV
     *
     *  @return the address (i.e. the ordering number) of the given spin-unresolved ONV
     */
    size_t addressOf(const SpinUnresolvedONV& onv) const {

        const auto& fock_space = this->derived();
        return fock_space.addressOf(onv.unsignedRepresentation());
    };

    /**
     *  @param address          the address (i.e. the ordening number) of a spin-unresolved ONV
     *
     *  @return the spin-unresolved ONV with the corresponding address
     */
    SpinUnresolvedONV constructONVFromAddress(const size_t address) const {

        const auto& fock_space = this->derived();

        SpinUnresolvedONV onv {fock_space.numberOfOrbitals(), this->N};
        fock_space.transformONVCorrespondingToAddress(onv, address);
        return onv;
    }

    /**
     *  @return the number of electrons that are associated to this ONV manipulator
     */
    size_t numberOfElectrons() const { return this->N; }

    /**
     *  Set the current SpinUnresolvedONV to the next SpinUnresolvedONV: performs nextPermutationOf() and updates the corresponding occupation indices of the SpinUnresolvedONV occupation array
     *
     *  @param onv      the current SpinUnresolvedONV
     */
    void transformONVToNextPermutation(SpinUnresolvedONV& onv) const {

        const auto& fock_space = this->derived();
        onv.replaceRepresentationWith(fock_space.nextPermutationOf(onv.unsignedRepresentation()));
    }

    /**
     *  Transform a spin-unresolved ONV to one corresponding to the given address
     *
     *  @param onv          the spin-unresolved ONV
     *  @param address      the address to which the spin-unresolved ONV will be set
     */
    void transformONVCorrespondingToAddress(SpinUnresolvedONV& onv, const size_t address) const {

        const auto& fock_space = this->derived();
        onv.replaceRepresentationWith((fock_space.representationOf(address)));
    }
};

}  // namespace GQCP