GSpinorBasis.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/ScalarBasis/GTOShell.hpp"
#include "Basis/ScalarBasis/LondonGTOShell.hpp"
#include "Basis/ScalarBasis/ScalarBasis.hpp"
#include "Basis/SpinorBasis/RSpinOrbitalBasis.hpp"
#include "Basis/SpinorBasis/SimpleSpinorBasis.hpp"
#include "Basis/SpinorBasis/USpinOrbitalBasis.hpp"
#include "Basis/Transformations/GTransformation.hpp"
#include "Domain/MullikenDomain/GMullikenDomain.hpp"
#include "Molecule/Molecule.hpp"
#include "Operator/FirstQuantized/CoulombRepulsionOperator.hpp"
#include "Operator/FirstQuantized/DiamagneticOperator.hpp"
#include "Operator/FirstQuantized/ElectronicDipoleOperator.hpp"
#include "Operator/FirstQuantized/ElectronicSpinOperator.hpp"
#include "Operator/FirstQuantized/ElectronicSpinSquaredOperator.hpp"
#include "Operator/FirstQuantized/ElectronicSpin_zOperator.hpp"
#include "Operator/FirstQuantized/FQMolecularHamiltonian.hpp"
#include "Operator/FirstQuantized/FQMolecularPauliHamiltonian.hpp"
#include "Operator/FirstQuantized/KineticOperator.hpp"
#include "Operator/FirstQuantized/NuclearAttractionOperator.hpp"
#include "Operator/FirstQuantized/OrbitalZeemanOperator.hpp"
#include "Operator/FirstQuantized/OverlapOperator.hpp"
#include "Operator/FirstQuantized/SpinZeemanOperator.hpp"
#include "Operator/SecondQuantized/GSQOneElectronOperator.hpp"
#include "Operator/SecondQuantized/GSQTwoElectronOperator.hpp"
#include "Operator/SecondQuantized/SQHamiltonian.hpp"
#include "Utilities/type_traits.hpp"
 
 
namespace GQCP {
 
 
template <typename _ExpansionScalar, typename _Shell>
class GSpinorBasis:
    public SimpleSpinorBasis<_ExpansionScalar, GSpinorBasis<_ExpansionScalar, _Shell>> {
public:
    // The scalar type used to represent an expansion coefficient of the spinors in the underlying scalar orbitals: real or complex.
    using ExpansionScalar = _ExpansionScalar;
 
    // The type of shell the underlying scalar bases contain.
    using Shell = _Shell;
 
    // The type of the base spinor basis.
    using Base = SimpleSpinorBasis<_ExpansionScalar, GSpinorBasis<_ExpansionScalar, _Shell>>;
 
    // The type of transformation that is naturally related to a `GSpinorBasis`.
    using Transformation = GTransformation<ExpansionScalar>;
 
    // The type that is used for representing the primitive for a basis function of this spin-orbital basis' underlying AO basis.
    using Primitive = typename Shell::Primitive;
 
    // The type that is used for representing the underlying basis functions of this spin-orbital basis.
    using BasisFunction = typename Shell::BasisFunction;
 
 
private:
    // The scalar bases for the alpha and beta components of the spinors.
    SpinResolved<ScalarBasis<Shell>> scalar_bases;
 
 
public:
    /*
     *  MARK: Constructors
     */
 
    GSpinorBasis(const SpinResolved<ScalarBasis<Shell>>& scalar_bases, const Transformation& C) :
        Base(C),
        scalar_bases {scalar_bases} {
 
        // Check if the dimensions of the given objects are compatible
        const auto K_alpha = scalar_bases.alpha().numberOfBasisFunctions();
        const auto K_beta = scalar_bases.beta().numberOfBasisFunctions();
 
        if (C.numberOfOrbitals() != K_alpha + K_beta) {
            throw std::invalid_argument("GSpinorBasis(const SpinResolved<ScalarBasis<Shell>>&, const Transformation&): The given dimensions of the scalar bases and coefficient matrix are incompatible.");
        }
    }
 
    GSpinorBasis(const ScalarBasis<Shell>& alpha_scalar_basis, const ScalarBasis<Shell>& beta_scalar_basis, const Transformation& C) :
        GSpinorBasis(SpinResolved<ScalarBasis<Shell>> {alpha_scalar_basis, beta_scalar_basis}, C) {}
 
 
    GSpinorBasis(const ScalarBasis<Shell>& scalar_basis, const Transformation& C) :
        GSpinorBasis(scalar_basis, scalar_basis, C) {}
 
 
    GSpinorBasis(const ScalarBasis<Shell>& alpha_scalar_basis, const ScalarBasis<Shell>& beta_scalar_basis) :
        GSpinorBasis(alpha_scalar_basis, beta_scalar_basis,
                     Transformation::Identity(alpha_scalar_basis.numberOfBasisFunctions(), beta_scalar_basis.numberOfBasisFunctions())) {}
 
 
    GSpinorBasis(const ScalarBasis<Shell>& scalar_basis) :
        GSpinorBasis(scalar_basis, scalar_basis) {}
 
 
    GSpinorBasis(const NuclearFramework& nuclear_framework, const std::string& basisset_name) :
        GSpinorBasis(ScalarBasis<Shell>(nuclear_framework, basisset_name)) {}
 
 
    GSpinorBasis(const Molecule& molecule, const std::string& basisset_name) :
        GSpinorBasis(ScalarBasis<Shell>(molecule.nuclearFramework(), basisset_name)) {}
 
 
    template <typename Z = Shell>
    GSpinorBasis(const Molecule& molecule, const std::string& basisset_name, const HomogeneousMagneticField& B,
                 typename std::enable_if<std::is_same<Z, LondonGTOShell>::value>::type* = 0) :
        GSpinorBasis(ScalarBasis<Shell>(molecule.nuclearFramework(), basisset_name, B)) {}
 
 
    GSpinorBasis(const NuclearFramework& nuclear_framework, const std::string& basisset_name_alpha, const std::string& basisset_name_beta) :
        GSpinorBasis(ScalarBasis<Shell>(nuclear_framework, basisset_name_alpha),
                     ScalarBasis<Shell>(nuclear_framework, basisset_name_beta)) {}
 
 
    GSpinorBasis(const Molecule& molecule, const std::string& basisset_name_alpha, const std::string& basisset_name_beta) :
        GSpinorBasis(ScalarBasis<Shell>(molecule.nuclearFramework(), basisset_name_alpha),
                     ScalarBasis<Shell>(molecule.nuclearFramework(), basisset_name_beta)) {}
 
 
    /*
     *  MARK: Named constructors
     */
 
    static GSpinorBasis<ExpansionScalar, Shell> FromRestricted(const RSpinOrbitalBasis<ExpansionScalar, Shell>& r_spinor_basis) {
 
        // Create an USpinOrbitalBasis from the restricted one and use ::FromUnrestricted.
        const auto u_spinor_basis = USpinOrbitalBasis<ExpansionScalar, Shell>::FromRestricted(r_spinor_basis);
 
        return GSpinorBasis<ExpansionScalar, Shell>::FromUnrestricted(u_spinor_basis);
    }
 
 
    static GSpinorBasis<ExpansionScalar, Shell> FromUnrestricted(const USpinOrbitalBasis<ExpansionScalar, Shell>& u_spinor_basis) {
 
        const auto C_general = GTransformation<ExpansionScalar>::FromUnrestricted(u_spinor_basis.expansion());
 
        return GSpinorBasis<ExpansionScalar, Shell>(u_spinor_basis.alpha().scalarBasis(), C_general);  // Assume the alpha- and beta- scalar bases are equal.
    }
 
 
    /*
     *  MARK: Coefficients
     */
 
    size_t numberOfCoefficients(const Spin& sigma) const { return this->scalarBases().component(sigma).numberOfBasisFunctions(); }
 
 
    /*
     *  MARK: General info
     */
 
    size_t numberOfSpinors() const {
 
        const auto K_alpha = this->numberOfCoefficients(Spin::alpha);
        const auto K_beta = this->numberOfCoefficients(Spin::beta);
 
        return K_alpha + K_beta;
    }
 
 
    /*
     *  MARK: Quantization of first-quantized operators (GTOShell)
     */
 
    template <typename FQOneElectronOperator, typename Z = Shell>
    auto quantize(const FQOneElectronOperator& fq_one_op) const -> enable_if_t<std::is_same<Z, GTOShell>::value, GSQOneElectronOperator<product_t<typename FQOneElectronOperator::Scalar, ExpansionScalar>, typename FQOneElectronOperator::Vectorizer>> {
 
        using ResultScalar = product_t<typename FQOneElectronOperator::Scalar, ExpansionScalar>;
        using ResultOperator = GSQOneElectronOperator<ResultScalar, typename FQOneElectronOperator::Vectorizer>;
 
 
        // The strategy for calculating the matrix representation of the one-electron operator in this spinor basis is to:
        //  1. Express the operator in the underlying scalar bases; and
        //  2. Afterwards transform them using the current coefficient matrix.
        const auto K_alpha = this->numberOfCoefficients(Spin::alpha);
        const auto K_beta = this->numberOfCoefficients(Spin::beta);
        const auto M = this->numberOfSpinors();
        SquareMatrix<ResultScalar> f = SquareMatrix<ResultScalar>::Zero(M);  // the total result
 
        // 1. Express the operator in the underlying scalar bases: spin-independent operators only have alpha-alpha and beta-beta blocks.
        const auto F_alpha = IntegralCalculator::calculateLibintIntegrals(fq_one_op, this->scalarBases().alpha());
        const auto F_beta = IntegralCalculator::calculateLibintIntegrals(fq_one_op, this->scalarBases().beta());
 
        f.topLeftCorner(K_alpha, K_alpha) = F_alpha;
        f.bottomRightCorner(K_beta, K_beta) = F_beta;
 
        // 2. Transform using the current coefficient matrix.
        ResultOperator op {{f}};  // op for 'operator'
        op.transform(this->expansion());
        return op;
    }
 
 
    template <typename Z = Shell>
    auto quantize(const ElectronicSpinOperator& fq_one_op) const -> enable_if_t<std::is_same<Z, GTOShell>::value, GSQOneElectronOperator<product_t<ElectronicSpinOperator::Scalar, ExpansionScalar>, ElectronicSpinOperator::Vectorizer>> {
 
        using ResultScalar = product_t<ElectronicSpinOperator::Scalar, ExpansionScalar>;
        using ResultOperator = GSQOneElectronOperator<ResultScalar, ElectronicSpinOperator::Vectorizer>;
 
        const auto K_alpha = this->numberOfCoefficients(Spin::alpha);
        const auto K_beta = this->numberOfCoefficients(Spin::beta);
        const auto M = this->numberOfSpinors();
 
        // The strategy to quantize the spin operator is as follows.
        //  1. First, calculate the necessary overlap integrals over the scalar bases.
        //  2. Then, construct the scalar basis representations of the components of the spin operator by placing the overlaps into the correct blocks.
        //  3. Transform the components (in scalar basis) with the current coefficient matrix to yield the components in spinor basis.
 
        SquareMatrix<ResultScalar> S_x = SquareMatrix<ResultScalar>::Zero(M);
        SquareMatrix<ResultScalar> S_y = SquareMatrix<ResultScalar>::Zero(M);
        SquareMatrix<ResultScalar> S_z = SquareMatrix<ResultScalar>::Zero(M);
 
 
        // 1. Calculate the necessary overlap integrals over the scalar bases.
        const auto S_aa = IntegralCalculator::calculateLibintIntegrals(OverlapOperator(), this->scalarBases().alpha());
        const auto S_ab = IntegralCalculator::calculateLibintIntegrals(OverlapOperator(), this->scalarBases().alpha(), this->scalarBases().beta());
        const auto S_ba = IntegralCalculator::calculateLibintIntegrals(OverlapOperator(), this->scalarBases().beta(), this->scalarBases().alpha());
        const auto S_bb = IntegralCalculator::calculateLibintIntegrals(OverlapOperator(), this->scalarBases().beta());
 
 
        // 2. Place the overlaps into the correct blocks.
        S_x.block(0, K_alpha, K_alpha, K_beta) = 0.5 * S_ab;
        S_x.block(K_alpha, 0, K_beta, K_alpha) = 0.5 * S_ba;
 
        using namespace GQCP::literals;
        S_y.block(0, K_alpha, K_alpha, K_beta) = -0.5 * 1_ii * S_ab;
        S_y.block(K_alpha, 0, K_beta, K_alpha) = 0.5 * 1_ii * S_ba;
 
        S_z.topLeftCorner(K_alpha, K_alpha) = 0.5 * S_aa;
        S_z.bottomRightCorner(K_beta, K_beta) = -0.5 * S_bb;
 
 
        // 3. Transform using the coefficient matrix
        ResultOperator spin_op {std::vector<SquareMatrix<ResultScalar>> {S_x, S_y, S_z}};  // 'op' for operator
        spin_op.transform(this->expansion());
        return spin_op;
    }
 
 
    template <typename Z = Shell>
    auto quantize(const ElectronicSpin_zOperator& fq_one_op) const -> enable_if_t<std::is_same<Z, GTOShell>::value, GSQOneElectronOperator<product_t<ElectronicSpin_zOperator::Scalar, ExpansionScalar>, ElectronicSpin_zOperator::Vectorizer>> {
 
        using ResultScalar = product_t<ElectronicSpin_zOperator::Scalar, ExpansionScalar>;
        using ResultOperator = GSQOneElectronOperator<ResultScalar, ElectronicSpin_zOperator::Vectorizer>;
 
        const auto K_alpha = this->numberOfCoefficients(Spin::alpha);
        const auto K_beta = this->numberOfCoefficients(Spin::beta);
        const auto M = this->numberOfSpinors();
 
        // The strategy to quantize the spin operator is as follows.
        //  1. First, calculate the necessary overlap integrals over the scalar bases.
        //  2. Then, construct the scalar basis representations of the components of the spin operator by placing the overlaps into the correct blocks.
        //  3. Transform the components (in scalar basis) with the current coefficient matrix to yield the components in spinor basis.
        SquareMatrix<ResultScalar> S_z = SquareMatrix<ResultScalar>::Zero(M);
 
        // 1. Calculate the necessary overlap integrals over the scalar bases.
        const auto S_aa = IntegralCalculator::calculateLibintIntegrals(OverlapOperator(), this->scalarBases().alpha());
        const auto S_bb = IntegralCalculator::calculateLibintIntegrals(OverlapOperator(), this->scalarBases().beta());
 
        // 2. Place the overlaps into the correct blocks.
        S_z.topLeftCorner(K_alpha, K_alpha) = 0.5 * S_aa;
        S_z.bottomRightCorner(K_beta, K_beta) = -0.5 * S_bb;
 
        // 3. Transform using the coefficient matrix.
        ResultOperator spin_op {{S_z}};  // 'op' for operator
        spin_op.transform(this->expansion());
        return spin_op;
    }
 
 
    template <typename Z = Shell>
    auto quantize(const CoulombRepulsionOperator& coulomb_op) const -> enable_if_t<std::is_same<Z, GTOShell>::value, GSQTwoElectronOperator<product_t<CoulombRepulsionOperator::Scalar, ExpansionScalar>, CoulombRepulsionOperator::Vectorizer>> {
 
        using ResultScalar = product_t<CoulombRepulsionOperator::Scalar, ExpansionScalar>;
        using ResultOperator = GSQTwoElectronOperator<ResultScalar, CoulombRepulsionOperator::Vectorizer>;
 
        // The strategy for calculating the matrix representation of the two-electron operator in this spinor basis is to:
        //  1. Calculate the Coulomb integrals in the underlying scalar bases;
        //  2. Place the calculated integrals as 'blocks' in the larger representation, so that we can;
        //  3. Transform the operator using the current coefficient matrix.
 
        // 1. Calculate the Coulomb integrals in the underlying scalar bases.
        const auto g_aaaa = IntegralCalculator::calculateLibintIntegrals(CoulombRepulsionOperator(), this->scalarBases().alpha());
        const auto g_aabb = IntegralCalculator::calculateLibintIntegrals(CoulombRepulsionOperator(), this->scalarBases().alpha(), this->scalarBases().beta());
        const auto g_bbaa = IntegralCalculator::calculateLibintIntegrals(CoulombRepulsionOperator(), this->scalarBases().beta(), this->scalarBases().alpha());
        const auto g_bbbb = IntegralCalculator::calculateLibintIntegrals(CoulombRepulsionOperator(), this->scalarBases().beta());
 
 
        // 2. Place the calculated integrals as 'blocks' in the larger representation
        const auto K_alpha = this->numberOfCoefficients(Spin::alpha);
        const auto K_beta = this->numberOfCoefficients(Spin::beta);
 
        const auto M = this->numberOfSpinors();
        auto g_par = SquareRankFourTensor<ResultScalar>::Zero(M);  // 'par' for 'parameters'
 
        // Primed indices are indices in the larger representation, normal ones are those in the smaller tensors.
        for (size_t mu_ = 0; mu_ < M; mu_++) {  // mu 'prime'
            const size_t mu = mu_ % K_alpha;
 
            for (size_t nu_ = 0; nu_ < M; nu_++) {  // nu 'prime'
                const size_t nu = nu_ % K_alpha;
 
                for (size_t rho_ = 0; rho_ < M; rho_++) {  // rho 'prime'
                    const size_t rho = rho_ % K_alpha;
 
                    for (size_t lambda_ = 0; lambda_ < M; lambda_++) {  // lambda 'prime'
                        const size_t lambda = lambda_ % K_alpha;
 
                        if ((mu_ < K_alpha) && (nu_ < K_alpha) && (rho_ < K_alpha) && (lambda_ < K_alpha)) {
                            g_par(mu_, nu_, rho_, lambda_) = g_aaaa(mu, nu, rho, lambda);
                        } else if ((mu_ < K_alpha) && (nu_ < K_alpha) && (rho_ >= K_alpha) && (lambda_ >= K_alpha)) {
                            g_par(mu_, nu_, rho_, lambda_) = g_aabb(mu, nu, rho, lambda);
                        } else if ((mu_ >= K_alpha) && (nu_ >= K_alpha) && (rho_ < K_alpha) && (lambda_ < K_alpha)) {
                            g_par(mu_, nu_, rho_, lambda_) = g_bbaa(mu, nu, rho, lambda);
                        } else if ((mu_ >= K_alpha) && (nu_ >= K_alpha) && (rho_ >= K_alpha) && (lambda_ >= K_alpha)) {
                            g_par(mu_, nu_, rho_, lambda_) = g_bbbb(mu, nu, rho, lambda);
                        }
                    }
                }
            }
        }
 
 
        // 3. Transform the operator using the current coefficient matrix.
        ResultOperator g_op {g_par};  // 'op' for 'operator'
        g_op.transform(this->expansion());
        return g_op;
    }
 
 
    /*
     *  MARK: Quantization of first-quantized operators (LondonGTOShell)
     */
 
    template <typename FQOneElectronOperator, typename Z = Shell>
    auto quantize(const FQOneElectronOperator& fq_one_op) const -> enable_if_t<std::is_same<Z, LondonGTOShell>::value, GSQOneElectronOperator<product_t<typename FQOneElectronOperator::Scalar, ExpansionScalar>, typename FQOneElectronOperator::Vectorizer>> {
 
        using ResultScalar = product_t<typename FQOneElectronOperator::Scalar, ExpansionScalar>;
        using ResultOperator = GSQOneElectronOperator<ResultScalar, typename FQOneElectronOperator::Vectorizer>;
        using Vectorizer = typename FQOneElectronOperator::Vectorizer;
 
        const auto N = FQOneElectronOperator::NumberOfComponents;
        const auto& vectorizer = FQOneElectronOperator::vectorizer;
 
 
        // The strategy for calculating the matrix representation of the one-electron operator in this spinor basis is to:
        //  1. Express the operator in the underlying scalar bases and;
        //  2. Afterwards transform them using the current coefficient matrix.
        const auto K_alpha = this->numberOfCoefficients(Spin::alpha);
        const auto K_beta = this->numberOfCoefficients(Spin::beta);
        const auto M = this->numberOfSpinors();
 
 
        // 1. Express the operator in the underlying scalar bases: spin-independent operators only have alpha-alpha and beta-beta blocks.
        auto engine = GQCP::IntegralEngine::InHouse<GQCP::LondonGTOShell>(fq_one_op);
        const auto F_aa = GQCP::IntegralCalculator::calculate(engine, this->scalarBases().alpha().shellSet(), this->scalarBases().alpha().shellSet());
        const auto F_bb = GQCP::IntegralCalculator::calculate(engine, this->scalarBases().beta().shellSet(), this->scalarBases().beta().shellSet());
 
        // For each of the components of the operator, place the scalar basis representations into the spinor basis representation.
        std::array<SquareMatrix<ResultScalar>, N> fs;
        for (size_t i = 0; i < N; i++) {
            SquareMatrix<ResultScalar> f = SquareMatrix<ResultScalar>::Zero(M);
 
            f.topLeftCorner(K_alpha, K_alpha) = F_aa[i];
            f.bottomRightCorner(K_beta, K_beta) = F_bb[i];
 
            fs[i] = f;
        }
 
 
        // 2. Transform using the current coefficient matrix.
        StorageArray<SquareMatrix<ResultScalar>, Vectorizer> array {fs, vectorizer};
        ResultOperator op {array};  // 'op' for 'operator'.
        op.transform(this->expansion());
        return op;
    }
 
 
    template <typename Z = Shell>
    auto quantize(const ElectronicSpinOperator& fq_one_op) const -> enable_if_t<std::is_same<Z, LondonGTOShell>::value, GSQOneElectronOperator<product_t<ElectronicSpinOperator::Scalar, ExpansionScalar>, ElectronicSpinOperator::Vectorizer>> {
 
        using ResultScalar = product_t<ElectronicSpinOperator::Scalar, ExpansionScalar>;
        using ResultOperator = GSQOneElectronOperator<ResultScalar, ElectronicSpinOperator::Vectorizer>;
 
        const auto K_alpha = this->numberOfCoefficients(Spin::alpha);
        const auto K_beta = this->numberOfCoefficients(Spin::beta);
        const auto M = this->numberOfSpinors();
 
        // The strategy to quantize the spin operator is as follows.
        //  1. First, calculate the necessary overlap integrals over the scalar bases.
        //  2. Then, construct the scalar basis representations of the components of the spin operator by placing the overlaps into the correct blocks.
        //  3. Transform the components (in scalar basis) with the current coefficient matrix to yield the components in spinor basis.
 
        SquareMatrix<ResultScalar> S_x = SquareMatrix<ResultScalar>::Zero(M);
        SquareMatrix<ResultScalar> S_y = SquareMatrix<ResultScalar>::Zero(M);
        SquareMatrix<ResultScalar> S_z = SquareMatrix<ResultScalar>::Zero(M);
 
 
        // 1. Calculate the necessary overlap integrals over the scalar bases.
        auto overlap_engine = GQCP::IntegralEngine::InHouse<GQCP::LondonGTOShell>(OverlapOperator());
        const auto S_aa = GQCP::IntegralCalculator::calculate(overlap_engine, this->scalarBases().alpha().shellSet(), this->scalarBases().alpha().shellSet())[0];
        const auto S_ab = GQCP::IntegralCalculator::calculate(overlap_engine, this->scalarBases().alpha().shellSet(), this->scalarBases().beta().shellSet())[0];
        const auto S_ba = GQCP::IntegralCalculator::calculate(overlap_engine, this->scalarBases().beta().shellSet(), this->scalarBases().alpha().shellSet())[0];
        const auto S_bb = GQCP::IntegralCalculator::calculate(overlap_engine, this->scalarBases().beta().shellSet(), this->scalarBases().beta().shellSet())[0];
 
 
        // 2. Place the overlaps into the correct blocks.
        S_x.block(0, K_alpha, K_alpha, K_beta) = 0.5 * S_ab;
        S_x.block(K_alpha, 0, K_beta, K_alpha) = 0.5 * S_ba;
 
        using namespace GQCP::literals;
        S_y.block(0, K_alpha, K_alpha, K_beta) = -0.5 * 1_ii * S_ab;
        S_y.block(K_alpha, 0, K_beta, K_alpha) = 0.5 * 1_ii * S_ba;
 
        S_z.topLeftCorner(K_alpha, K_alpha) = 0.5 * S_aa;
        S_z.bottomRightCorner(K_beta, K_beta) = -0.5 * S_bb;
 
 
        // 3. Transform using the coefficient matrix.
        ResultOperator spin_op {std::vector<SquareMatrix<ResultScalar>> {S_x, S_y, S_z}};  // 'op' for operator.
        spin_op.transform(this->expansion());
        return spin_op;
    }
 
 
    template <typename Z = Shell>
    auto quantize(const OrbitalZeemanOperator& op) const -> enable_if_t<std::is_same<Z, LondonGTOShell>::value, GSQOneElectronOperator<product_t<OrbitalZeemanOperator::Scalar, ExpansionScalar>, OrbitalZeemanOperator::Vectorizer>> {
 
        // Return the orbital Zeeman operator as a contraction beween the magnetic field and the angular momentum operator.
        const auto L = this->quantize(op.angularMomentum());
        const auto& B = op.magneticField().strength();
        return 0.5 * L.dot(B);
    }
 
 
    template <typename Z = Shell>
    auto quantize(const DiamagneticOperator& op) const -> enable_if_t<std::is_same<Z, LondonGTOShell>::value, GSQOneElectronOperator<product_t<DiamagneticOperator::Scalar, ExpansionScalar>, DiamagneticOperator::Vectorizer>> {
 
        using ResultScalar = product_t<DiamagneticOperator::Scalar, ExpansionScalar>;
        using ResultOperator = GSQOneElectronOperator<ResultScalar, DiamagneticOperator::Vectorizer>;
 
 
        // Return the diamagnetic operator as a contraction beween the magnetic field and the electronic quadrupole operator.
        const auto Q = this->quantize(op.electronicQuadrupole()).allParameters();
        const auto& B = op.magneticField().strength();
 
        // Prepare some variables.
        const auto& B_x = B(CartesianDirection::x);
        const auto& B_y = B(CartesianDirection::y);
        const auto& B_z = B(CartesianDirection::z);
 
        const auto& Q_xx = Q[DyadicCartesianDirection::xx];
        const auto& Q_xy = Q[DyadicCartesianDirection::xy];
        const auto& Q_xz = Q[DyadicCartesianDirection::xz];
        const auto& Q_yy = Q[DyadicCartesianDirection::yy];
        const auto& Q_yz = Q[DyadicCartesianDirection::yz];
        const auto& Q_zz = Q[DyadicCartesianDirection::zz];
 
 
        SquareMatrix<ResultScalar> D_par = 0.125 * ((std::pow(B_y, 2) + std::pow(B_z, 2)) * Q_xx +
                                                    (std::pow(B_x, 2) + std::pow(B_z, 2)) * Q_yy +
                                                    (std::pow(B_x, 2) + std::pow(B_y, 2)) * Q_zz -
                                                    2 * B_x * B_y * Q_xy -
                                                    2 * B_x * B_z * Q_xz -
                                                    2 * B_y * B_z * Q_yz);
 
        return ResultOperator {D_par};
    }
 
 
    template <typename Z = Shell>
    auto quantize(const SpinZeemanOperator& op) const -> enable_if_t<std::is_same<Z, LondonGTOShell>::value, GSQOneElectronOperator<product_t<SpinZeemanOperator::Scalar, ExpansionScalar>, SpinZeemanOperator::Vectorizer>> {
 
        // Return the spin Zeeman operator as a contraction beween the magnetic field and the spin operator.
        const auto S = this->quantize(ElectronicSpinOperator());
        const auto& B = op.magneticField().strength();
        return S.dot(B);
    }
 
 
    template <typename Z = Shell>
    auto quantize(const CoulombRepulsionOperator& coulomb_op) const -> enable_if_t<std::is_same<Z, LondonGTOShell>::value, GSQTwoElectronOperator<product_t<CoulombRepulsionOperator::Scalar, ExpansionScalar>, CoulombRepulsionOperator::Vectorizer>> {
 
        using ResultScalar = product_t<CoulombRepulsionOperator::Scalar, ExpansionScalar>;
        using ResultOperator = GSQTwoElectronOperator<ResultScalar, CoulombRepulsionOperator::Vectorizer>;
 
        // The strategy for calculating the matrix representation of the two-electron operator in this spinor basis is to:
        //  1. Calculate the Coulomb integrals in the underlying scalar bases;
        //  2. Place the calculated integrals as 'blocks' in the larger representation, so that we can;
        //  3. Transform the operator using the current coefficient matrix.
 
        // 1. Calculate the Coulomb integrals in the underlying scalar bases.
        auto coulomb_engine = GQCP::IntegralEngine::InHouse<GQCP::LondonGTOShell>(CoulombRepulsionOperator());
        const auto g = GQCP::IntegralCalculator::calculate(coulomb_engine, this->scalarBases().alpha().shellSet(), this->scalarBases().alpha().shellSet())[0];
 
 
        // 2. Place the calculated integrals as 'blocks' in the larger representation.
        const auto K_alpha = this->numberOfCoefficients(Spin::alpha);
        const auto K_beta = this->numberOfCoefficients(Spin::beta);
 
        const auto M = this->numberOfSpinors();
        auto g_par = SquareRankFourTensor<ResultScalar>::Zero(M);  // 'par' for 'parameters'.
 
        // Primed indices are indices in the larger representation, normal ones are those in the smaller tensors.
        for (size_t mu_ = 0; mu_ < M; mu_++) {  // Mu 'prime'.
            const size_t mu = mu_ % K_alpha;
 
            for (size_t nu_ = 0; nu_ < M; nu_++) {  // Nu 'prime'.
                const size_t nu = nu_ % K_alpha;
 
                for (size_t rho_ = 0; rho_ < M; rho_++) {  // Rho 'prime'.
                    const size_t rho = rho_ % K_alpha;
 
                    for (size_t lambda_ = 0; lambda_ < M; lambda_++) {  // Lambda 'prime'.
                        const size_t lambda = lambda_ % K_alpha;
 
                        if ((mu_ < K_alpha) && (nu_ < K_alpha) && (rho_ < K_alpha) && (lambda_ < K_alpha)) {
                            g_par(mu_, nu_, rho_, lambda_) = g(mu, nu, rho, lambda);
                        } else if ((mu_ < K_alpha) && (nu_ < K_alpha) && (rho_ >= K_alpha) && (lambda_ >= K_alpha)) {
                            g_par(mu_, nu_, rho_, lambda_) = g(mu, nu, rho, lambda);
                        } else if ((mu_ >= K_alpha) && (nu_ >= K_alpha) && (rho_ < K_alpha) && (lambda_ < K_alpha)) {
                            g_par(mu_, nu_, rho_, lambda_) = g(mu, nu, rho, lambda);
                        } else if ((mu_ >= K_alpha) && (nu_ >= K_alpha) && (rho_ >= K_alpha) && (lambda_ >= K_alpha)) {
                            g_par(mu_, nu_, rho_, lambda_) = g(mu, nu, rho, lambda);
                        }
                    }
                }
            }
        }
 
 
        // 3. Transform the operator using the current coefficient matrix.
        ResultOperator g_op {g_par};  // 'op' for 'operator'.
        g_op.transform(this->expansion());
        return g_op;
    }
 
 
    template <typename Z = Shell>
    enable_if_t<std::is_same<Z, LondonGTOShell>::value, GSQHamiltonian<ExpansionScalar>> quantize(const FQMolecularMagneticHamiltonian& fq_hamiltonian) const {
 
        const auto T = this->quantize(fq_hamiltonian.kinetic());
        const auto OZ = this->quantize(fq_hamiltonian.orbitalZeeman());
        const auto D = this->quantize(fq_hamiltonian.diamagnetic());
 
        const auto V = this->quantize(fq_hamiltonian.nuclearAttraction());
 
        const auto g = this->quantize(fq_hamiltonian.coulombRepulsion());
 
        return GSQHamiltonian<ExpansionScalar> {{T, OZ, D, V}, {g}};
    }
 
 
    template <typename Z = Shell>
    enable_if_t<std::is_same<Z, LondonGTOShell>::value, GSQHamiltonian<ExpansionScalar>> quantize(const FQMolecularPauliHamiltonian& fq_hamiltonian) const {
 
        const auto T = this->quantize(fq_hamiltonian.kinetic());
        const auto OZ = this->quantize(fq_hamiltonian.orbitalZeeman());
        const auto D = this->quantize(fq_hamiltonian.diamagnetic());
 
        const auto SZ = this->quantize(fq_hamiltonian.spinZeeman());
 
        const auto V = this->quantize(fq_hamiltonian.nuclearAttraction());
 
        const auto g = this->quantize(fq_hamiltonian.coulombRepulsion());
 
        return GSQHamiltonian<ExpansionScalar> {{T, OZ, D, SZ, V}, {g}};
    }
 
 
    /*
     *  MARK: Quantization of first-quantized operators
     */
 
    GSQHamiltonian<ExpansionScalar> quantize(const FQMolecularHamiltonian& fq_hamiltonian) const {
 
        const auto T = this->quantize(fq_hamiltonian.kinetic());
        const auto V = this->quantize(fq_hamiltonian.nuclearAttraction());
 
        const auto g = this->quantize(fq_hamiltonian.coulombRepulsion());
 
        return GSQHamiltonian<ExpansionScalar> {T + V, g};
    }
 
 
    template <typename S = ExpansionScalar>
    enable_if_t<std::is_same<S, complex>::value, ScalarGSQOneElectronOperatorProduct<complex>> quantize(const ElectronicSpinSquaredOperator& fq_S2_op) const {
 
        // Prepare the 3 vector components of the electronic spin operator S = (S_x, S_y, S_z).
        const auto S_op = this->quantize(ElectronicSpinOperator());
        const auto& S_x_op = S_op(CartesianDirection::x);
        const auto& S_y_op = S_op(CartesianDirection::y);
        const auto& S_z_op = S_op(CartesianDirection::z);
 
 
        // Calculate S^2 = S_x^2 + S_y^2 + S_z^2.
        const auto S_x_2 = S_x_op * S_x_op;
        const auto S_y_2 = S_y_op * S_y_op;
        const auto S_z_2 = S_z_op * S_z_op;
 
        return S_x_2 + S_y_2 + S_z_2;
    }
 
 
    template <typename S = ExpansionScalar>
    enable_if_t<std::is_same<S, double>::value, ScalarGSQOneElectronOperatorProduct<double>> quantize(const ElectronicSpinSquaredOperator& fq_S2_op) const {
        // Define the dimensions and empty matrices.
        const auto K = this->numberOfCoefficients(Spin::alpha);
        const auto M = this->numberOfSpinors();
 
        SquareMatrix<double> S_plus = SquareMatrix<double>::Zero(M);
        SquareMatrix<double> S_minus = SquareMatrix<double>::Zero(M);
        SquareMatrix<double> S_plus_minus = SquareMatrix<double>::Zero(M);
        SquareMatrix<double> S_z_squared = SquareMatrix<double>::Zero(M);
 
        // Calculate the necessary overlap integrals over the scalar bases.
        const auto S_aa = IntegralCalculator::calculateLibintIntegrals(OverlapOperator(), this->scalarBases().alpha());
 
        // Calculate the S_z operator.
        const auto S_z = this->quantize(ElectronicSpin_zOperator()).parameters();
 
        // Fill in the overlap in the correct positions.
        S_plus.topRightCorner(K, K) = S_aa;
        S_minus.bottomLeftCorner(K, K) = S_aa;
        S_plus_minus.bottomRightCorner(K, K) = S_aa;
 
        S_z_squared.topLeftCorner(K, K) = 0.25 * S_aa;
        S_z_squared.bottomRightCorner(K, K) = 0.25 * S_aa;
 
        // Create the one electron part.
        ScalarGSQOneElectronOperator<double> one_electron_part {S_z_squared + S_z + S_plus_minus};
 
        // Determine the tensors of the two-electron contributions.
        SquareRankFourTensor<double> T_S_z = SquareRankFourTensor<double>::Zero(M);
        SquareRankFourTensor<double> T_S_minus_plus = SquareRankFourTensor<double>::Zero(M);
        for (size_t p = 0; p < M; p++) {
            for (size_t q = 0; q < M; q++) {
                for (size_t r = 0; r < M; r++) {
                    for (size_t s = 0; s < M; s++) {
                        T_S_z(p, q, r, s) = 2.0 * S_z(p, q) * S_z(r, s);                  // Include the prefactor '2' because we're going to encapsulate these matrix elements with a `ScalarGSQTwoElectronOperator`, whose matrix elements should not embet the prefactor 0.5 for two-electron operators.
                        T_S_minus_plus(p, q, r, s) = 2.0 * S_minus(p, q) * S_plus(r, s);  // Include the prefactor '2' because we're going to encapsulate these matrix elements with a `ScalarGSQTwoElectronOperator`, whose matrix elements should not embet the prefactor 0.5 for two-electron operators.
                    }
                }
            }
        }
        // Create the two electron part.
        ScalarGSQTwoElectronOperator<double> two_electron_part = ScalarGSQTwoElectronOperator<double> {T_S_z} + ScalarGSQTwoElectronOperator<double> {T_S_minus_plus};
 
        return ScalarGSQOneElectronOperatorProduct<double>(one_electron_part, two_electron_part);
    }
 
 
    /*
     *  MARK: Scalar basis
     */
 
    const SpinResolved<ScalarBasis<Shell>>& scalarBases() const { return this->scalar_bases; }
 
 
    GMullikenDomain<ExpansionScalar> mullikenDomain(const std::function<bool(const BasisFunction&)>& selector) const {
 
        // Assume the underlying scalar bases are equal, and proceed to work with the one for the alpha component.
        const auto ao_indices = this->scalarBases().alpha().basisFunctionIndices(selector);
        return GMullikenDomain<ExpansionScalar> {ao_indices, ao_indices.size()};
    }
 
 
    GMullikenDomain<ExpansionScalar> mullikenDomain(const std::function<bool(const Shell&)>& selector) const {
 
        // Assume the underlying scalar bases are equal, and proceed to work with the one for the alpha component.
        const auto ao_indices = this->scalarBases().alpha().basisFunctionIndices(selector);
        return GMullikenDomain<ExpansionScalar> {ao_indices, ao_indices.size()};
    }
};
 
 
/*
 *  MARK: SpinorBasisTraits
 */
 
template <typename _ExpansionScalar, typename _Shell>
struct SpinorBasisTraits<GSpinorBasis<_ExpansionScalar, _Shell>> {
    // The scalar type used to represent an expansion coefficient of the spinors in the underlying scalar orbitals: real or complex.
    using ExpansionScalar = _ExpansionScalar;
 
    // The type of transformation that is naturally related to a `GSpinorBasis`.
    using Transformation = GTransformation<ExpansionScalar>;
 
    // The second-quantized representation of the overlap operator related to the derived spinor basis.
    using SQOverlapOperator = ScalarGSQOneElectronOperator<ExpansionScalar>;
};
 
 
/*
 *  MARK: BasisTransformableTraits
 */
 
template <typename _ExpansionScalar, typename _Shell>
struct BasisTransformableTraits<GSpinorBasis<_ExpansionScalar, _Shell>> {
 
    // The type of transformation that is naturally related to a `GSpinorBasis`.
    using Transformation = GTransformation<_ExpansionScalar>;
};
 
 
/*
 *  MARK: JacobiRotatableTraits
 */
 
template <typename _ExpansionScalar, typename _Shell>
struct JacobiRotatableTraits<GSpinorBasis<_ExpansionScalar, _Shell>> {
 
    // The type of Jacobi rotation that is naturally related to a `GSpinorBasis`.
    using JacobiRotationType = JacobiRotation;
};
 
 
}  // namespace GQCP