// 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/Transformations/BasisTransformable.hpp"
#include "Mathematical/Representation/SquareRankFourTensor.hpp"
#include "Molecule/Molecule.hpp"
#include "Operator/FirstQuantized/AngularMomentumOperator.hpp"
#include "Operator/FirstQuantized/CoulombRepulsionOperator.hpp"
#include "Operator/FirstQuantized/DiamagneticOperator.hpp"
#include "Operator/FirstQuantized/ElectronicQuadrupoleOperator.hpp"
#include "Operator/FirstQuantized/FQMolecularHamiltonian.hpp"
#include "Operator/FirstQuantized/FQMolecularMagneticHamiltonian.hpp"
#include "Operator/FirstQuantized/KineticOperator.hpp"
#include "Operator/FirstQuantized/NuclearAttractionOperator.hpp"
#include "Operator/FirstQuantized/OrbitalZeemanOperator.hpp"
#include "Operator/FirstQuantized/OverlapOperator.hpp"
#include "Operator/SecondQuantized/OperatorTraits.hpp"
#include "Partition/DomainPartition.hpp"
#include "Partition/SimplePartition.hpp"
#include "Physical/HomogeneousMagneticField.hpp"
#include "QuantumChemical/SpinResolvedBase.hpp"
#include "gqcpy/include/utilities.hpp"

#include <pybind11/eigen/matrix.h>
#include <pybind11/eigen/tensor.h>
#include <pybind11/functional.h>
#include <pybind11/iostream.h>
#include <pybind11/operators.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>


namespace gqcpy {


// Provide some shortcuts for frequent namespaces.
namespace py = pybind11;
using namespace GQCP;


/**
 *  This header contains methods that bind the public interface of some C++ classes/interfaces to Pybind11 classes.
 */

/*
 *  MARK: Abstract interfaces
 */

/**
 *  Add Python bindings for the `BasisTransformable` interface.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `BasisTransformable` APIs.
 */
template <typename Class>
void bindBasisTransformableInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Transformation = typename BasisTransformableTraits<Type>::Transformation;

    py_class
        .def(
            "rotate",
            [](Type& transformable, const Transformation& U) {
                transformable.rotate(U);
            },
            py::arg("U"),
            "In-place apply the basis rotation.")

        .def(
            "rotated",
            [](const Type& transformable, const Transformation& U) {
                return transformable.rotated(U);
            },
            py::arg("U"),
            "Apply the basis rotation and return the result.")

        .def(
            "transform",
            [](Type& transformable, const Transformation& T) {
                transformable.transform(T);
            },
            py::arg("T"),
            "In-place apply the basis transformation.")

        .def(
            "transformed",
            [](const Type& transformable, const Transformation& T) {
                return transformable.transformed(T);
            },
            py::arg("T"),
            "Apply the basis transformation and return the result.");
}


/**
 *  Add Python bindings for the `VectorSpaceArithmetic` C++ interface to a given `pybind11::class_`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain `VectorSpaceArithmetic`'s APIs.
 */
template <typename Class>
void bindVectorSpaceArithmeticInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;

    py_class
        .def(py::self += py::self)
        .def(py::self + py::self)

        .def(py::self -= py::self)
        .def(py::self - py::self)

        .def(py::self *= Scalar())
        .def(Scalar() * py::self)
        .def(py::self * Scalar())

        .def(py::self / Scalar())

        .def(-py::self);
}


/**
 *  Add Python bindings for the `SpinResolvedBase` C++ interface to a given `pybind11::class_`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain `SpinResolvedBase`'s APIs.
 */
template <typename Class>
void bindSpinResolvedBaseInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using ComponentType = typename Type::ComponentType;

    py_class

        .def(py::init<const ComponentType&, const ComponentType&>())

        .def_property(
            "alpha",
            [](const Type& spin_resolved_object) {
                return spin_resolved_object.alpha();
            },
            [](Type& spin_resolved_object, const ComponentType& new_component) {
                spin_resolved_object.alpha() = new_component;
            },
            "A read-writable reference to the alpha object.")

        .def_property(
            "beta",
            [](const Type& spin_resolved_object) {
                return spin_resolved_object.beta();
            },
            [](Type& spin_resolved_object, const ComponentType& new_component) {
                spin_resolved_object.beta() = new_component;
            },
            "A read-writable reference to the beta object.")

        .def_property(
            "component",
            [](const Type& spin_resolved_object, const Spin sigma) {
                return spin_resolved_object.component(sigma);
            },
            [](Type& spin_resolved_object, const ComponentType& new_component, const Spin sigma) {
                spin_resolved_object.component(sigma) = new_component;
            },
            "A read-writable reference to the alpha orbeta object.");
}


/**
 *  Add some APIs related to operations on `BasisTransformable` objects.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `BasisTransformable` operations.
 */
template <typename Class>
void bindBasisTransformableOperationsInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class
        .def(
            "inverse",
            &Type::inverse,
            "Return the inverse transformation of this transformation matrix.")

        .def(
            "isUnitary",
            &Type::isUnitary,
            py::arg("threshold") = 1.0e-12,
            "Return if this transformation is unitary, within the given threshold");
}


/**
 *  Add Python bindings for the `DoublySpinResolvedBase` C++ interface to a given `pybind11::class_`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain `DoublySpinResolvedBase`'s APIs.
 */
template <typename Class>
void bindDoublySpinResolvedBaseInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using PureComponentType = typename Type::PureComponentType;
    using MixedComponentType = typename Type::MixedComponentType;

    py_class
        .def(py::init<const PureComponentType&, const MixedComponentType&, const MixedComponentType&, const PureComponentType&>())

        .def(
            "alphaAlpha",
            [](const Type& spin_resolved_object) {
                return spin_resolved_object.alphaAlpha();
            },
            "A read-only reference to the alpha-alpha object.")

        .def(
            "alphaAlpha",
            [](Type& spin_resolved_object) {
                return spin_resolved_object.alphaAlpha();
            },
            "A writable reference to the alpha-alpha object.")

        .def(
            "alphaBeta",
            [](const Type& spin_resolved_object) {
                return spin_resolved_object.alphaBeta();
            },
            "A read-only reference to the alpha-beta object.")

        .def(
            "alphaBeta",
            [](Type& spin_resolved_object) {
                return spin_resolved_object.alphaBeta();
            },
            "A writable reference to the alpha-beta object.")

        .def(
            "betaAlpha",
            [](const Type& spin_resolved_object) {
                return spin_resolved_object.betaAlpha();
            },
            "A read-only reference to the beta-alpha object.")

        .def(
            "betaAlpha",
            [](Type& spin_resolved_object) {
                return spin_resolved_object.betaAlpha();
            },
            "A writable reference to the beta-alpha object.")

        .def(
            "betaBeta",
            [](const Type& spin_resolved_object) {
                return spin_resolved_object.betaBeta();
            },
            "A read-only reference to the beta-beta object.")

        .def(
            "betaBeta",
            [](Type& spin_resolved_object) {
                return spin_resolved_object.betaBeta();
            },
            "A writable reference to the beta-beta object.");
}


/*
 *  MARK: Basis - SpinorBasis
 */

/**
 *  Add Python bindings for some APIs related to spinor bases.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to spinor bases.
 */
template <typename Class>
void bindSpinorBasisInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class

        /*
         *  MARK: Coefficient access
         */

        .def(
            "expansion",
            &Type::expansion,
            "A read-only reference to the transformation that relates the current set of spinors with the atomic spinors.")


        /*
         *  MARK: Orthonormality
         */

        .def(
            "overlap",
            &Type::overlap,
            "The overlap (one-electron) operator of this spinor basis.")

        .def(
            "isOrthonormal",
            &Type::isOrthonormal,
            py::arg("precision") = 1.0e-08,
            "Return if this spinor basis is orthonormal within the given precision")

        .def(
            "lowdinOrthonormalization",
            &Type::lowdinOrthonormalization,
            "The transformation to the Löwdin basis: T = S_current^{-1/2}.")

        .def(
            "lowdinOrthonormalize",
            &Type::lowdinOrthonormalize,
            "Transform the spinor basis to the 'Löwdin basis', which is the orthonormal basis that we transform to with T = S^{-1/2}, where S is the current overlap matrix.");


    /*
     *  MARK: Conforming to `BasisTransformable`
     */
    bindBasisTransformableInterface(py_class);
}


/**
 *  Add Python bindings for the Mulliken domains for spinor bases.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to Mulliken domain.
 */
template <typename Class>
void bindSpinorBasisMullikenInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class.def(
        "mullikenDomain",
        [](const Type& spinor_basis, const std::function<bool(const GTOShell&)>& selector) {
            return spinor_basis.mullikenDomain(selector);
        },
        py::arg("selector"),
        "A Mulliken domain of the AOs selected by the supplied selector function.");
}


/**
 *  Add Python bindings for the quantization of first-quantized operators.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to the quantization of first-quantized operators.
 */
template <typename Class>
void bindSpinorBasisQuantizationInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class
        .def(
            "quantize",
            [](const Type& spinor_basis, const CoulombRepulsionOperator& op) {
                return spinor_basis.quantize(op);
            },
            "Return the Coulomb operator expressed in this spinor basis.")

        .def(
            "quantize",
            [](const Type& spinor_basis, const KineticOperator& op) {
                return spinor_basis.quantize(op);
            },
            "Return the kinetic energy operator expressed in this spinor basis.")

        .def(
            "quantize",
            [](const Type& spinor_basis, const NuclearAttractionOperator& op) {
                return spinor_basis.quantize(op);
            },
            "Return the nuclear attraction operator expressed in this spinor basis.")

        .def(
            "quantize",
            [](const Type& spinor_basis, const OverlapOperator& op) {
                return spinor_basis.quantize(op);
            },
            "Return the overlap operator expressed in this spinor basis.")

        .def(
            "quantize",
            [](const Type& spinor_basis, const FQMolecularHamiltonian& op) {
                return spinor_basis.quantize(op);
            },
            "Return the molecular Hamiltonian expressed in this spinor basis.");
}


/**
 *  Add Python bindings for the quantization of first-quantized operators in a London spinor basis.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to the quantization of first-quantized operators in a London spinor basis
 */
template <typename Class>
void bindLondonSpinorBasisQuantizationInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class
        .def(py::init<const Molecule&, const std::string&, const HomogeneousMagneticField&>(),
             py::arg("molecule"),
             py::arg("basisset_name"),
             py::arg("B"))

        .def(
            "quantize",
            [](const Type& spinor_basis, const AngularMomentumOperator& op) {
                return spinor_basis.quantize(op);
            },
            "Return the angular momentum operator expressed in this spinor basis.")

        .def(
            "quantize",
            [](const Type& spinor_basis, const ElectronicQuadrupoleOperator& op) {
                return spinor_basis.quantize(op);
            },
            "Return the electronic quadrupole operator expressed in this spinor basis.")


        .def(
            "quantize",
            [](const Type& spinor_basis, const OrbitalZeemanOperator& op) {
                return spinor_basis.quantize(op);
            },
            "Return the orbital Zeeman operator expressed in this spinor basis.")

        .def(
            "quantize",
            [](const Type& spinor_basis, const DiamagneticOperator& op) {
                return spinor_basis.quantize(op);
            },
            "Return the diamagnetic operator expressed in this spinor basis.")

        .def(
            "quantize",
            [](const Type& spinor_basis, const FQMolecularMagneticHamiltonian& op) {
                return spinor_basis.quantize(op);
            },
            "Return the molecular magnetic Hamiltonian expressed in this spinor basis.");
}


/*
 * MARK: Basis - Transformations
 */

/**
 *  Add Python bindings for some APIs related to `SimpleTransformation`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `SimpleTransformation`.
 */
template <typename Class>
void bindSimpleTransformationInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;


    py_class

        /*
         *  MARK: Constructors
         */

        .def(py::init<>([](const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>& T) {
                 return Type {T};
             }),
             py::arg("T"))


        /*
         *  MARK: Transformation matrix
         */

        .def("matrix",
             &Type::matrix,
             "Return the transformation matrix that collects the expansion coefficients of the new basis (vectors) in the old basis as columns.");


    // Expose the `BasisTransformable` APIs.
    bindBasisTransformableInterface(py_class);

    // Add some APIs related to operations on `BasisTransformable` objects.
    bindBasisTransformableOperationsInterface(py_class);
}


/*
 *  MARK: Domain - SimpleDomain
 */

/**
 *  Add Python bindings to a type that represents a Simple domain.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to Simple domain.
 */
template <typename Class>
void bindSimpleDomainInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class
        .def(
            "__eq__",
            [](const Type& domain, const Type& other) {
                return domain == other;
            },
            py::arg("other"),
            "Returns whether this domain is equal to another domain.")

        .def(
            "__neq__",
            [](const Type& domain, const Type& other) {
                return domain != other;
            },
            py::arg("other"),
            "Returns whether this domain is not equal to another domain.")

        .def(
            "__repr__",
            &Type::asString)

        .def(
            "asString",
            &Type::asString)

        .def(
            "inDomain",
            &Type::inDomain);
}


/*
 *  MARK: Domain - DiscreteDomain
 */

/**
 *  Add Python bindings for some APIs related to the construction of `DiscreteDomain`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to the construction of `DiscreteDomain`.
 */
template <typename Class>
void bindDiscreteDomainConstructorInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class

        /*
         *  MARK: Constructors
         */

        .def(py::init<std::vector<size_t>&, const size_t>(),
             py::arg("domain_indices"),
             py::arg("M"),
             "Create a discrete domain of size M containing the given domain indices.")

        .def(py::init<const size_t, const size_t>(),
             py::arg("unsigned_representation"),
             py::arg("M"),
             "Create a discrete domain of size M from a given unsigned integer.");
}


/**
 *  Add Python bindings to a type that represents a Discrete domain.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to Discrete domain.
 */
template <typename Class>
void bindDiscreteDomainInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class
        .def("addElement",
             &Type::addElement,
             py::arg("i"),
             "Add an element to the domain at position i.")

        .def("asString",
             &Type::asString,
             "Return the domain string as a string.")

        .def("dimension",
             &Type::dimension,
             "Return the dimension of this discrete domain, i.e. the maximum number of elements that the discrete domain can contain.")

        .def("domainIndices",
             &Type::domainIndices,
             "Return the set of indices that correspond to the elements that are included in the domain.")

        .def("numberOfElements",
             &Type::numberOfElements,
             "Return the number of elements the domain contains.")

        .def(
            "overlapWith",
            [](const Type& domain, const Type& other) {
                return domain.overlapWith(other);
            },
            py::arg("domain"),
            "Return the number of domain elements that are equal between `this` and the other discrete domain.")

        .def("removeElement",
             &Type::removeElement,
             py::arg("i"),
             "Remove an element from the domain at position i.")

        .def("unsignedRepresentation",
             &Type::unsignedRepresentation,
             "Return the unsigned representation of this domain.");

    // Expose the `SimpleDomain` APIs.
    bindSimpleDomainInterface(py_class);
}


/*
 *  MARK: Domain - MullikenDomain
 */


/**
 *  Add Python bindings for the Mulliken matrices to a type that represents a Mulliken domain.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to Mulliken domain.
 */
template <typename Class>
void bindMullikenDomainMatricesInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class
        .def("partitionMatrix",
             &Type::partitionMatrix,
             py::arg("C"),
             "Return the partition matrix 'P_A' related to this Mulliken domain.")

        .def("projectionMatrix",
             &Type::projectionMatrix,
             py::arg("C"),
             "Return the Mulliken projection, defined as C^{-1} P_A C, where C is the transformation matrix and P_A is the partition matrix.");
}


/*
 *  MARK: Partition - SimplePartition
 */

/**
 *  Add Python bindings to a type that represents a Simple partition.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `SimplePartition`.
 */
template <typename Class>
void bindSimplePartitionInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class

        /*
         *  MARK: Constructors
         */

        .def(py::init<std::vector<typename Type::ElementType>&>(),
             py::arg("partition_elements"),
             "Create a partition from a vector of `element`s.")

        /*
         *  MARK: General info
         */

        .def(
            "__eq__",
            [](const Type& domain_partition, const Type& other) {
                return domain_partition == other;
            },
            py::arg("other"),
            "Returns whether this partition is equal to the other partition.")

        .def(
            "__ne__",
            [](const Type& domain_partition, const Type& other) {
                return domain_partition != other;
            },
            py::arg("other"),
            "Returns whether this partition is not equal to the other partition.")

        .def(
            "__repr__",
            &Type::asString)

        .def(
            "__call__",
            &Type::operator())

        .def(
            "dimension",
            &Type::dimension)

        .def(
            "partitionElements",
            &Type::partitionElements);
}


/*
 *  MARK: Partition - DomainPartition
 */

/**
 *  Add Python bindings to a type that represents a domain partition.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `DomainPartition`.
 */
template <typename Class>
void bindDomainPartitionInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class

        .def(
            "asString",
            &Type::asString,
            "Return the domain string as a string.");

    // Expose the `SimplePartition` APIs.
    bindSimplePartitionInterface(py_class);
}


/*
 *  MARK: Operator - SecondQuantized
 */


/**
 *  Add Python bindings for some Mulliken domain APIs to a given `pybind11::class_` that represents a one-electron operator.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` (that represents a one-electron operator) that should obtain some Mulliken domain APIs.
 */
template <typename Class>
void bindSQOneElectronOperatorMullikenInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    py_class.def("partitioned",
                 &Type::partitioned,
                 py::arg("P"),
                 "Partition this one-electron operator according to the supplied Mulliken partitioning scheme.");
}


/**
 *  Add Python bindings for APIs related to one-electron operators to a given `pybind11::class_`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain one-electron operator APIs.
 */
template <typename Class>
void bindSQOneElectronOperatorInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    using OneDM = typename OperatorTraits<Type>::OneDM;


    py_class

        /*
         *  MARK: Calculations
         */

        .def(
            "calculateExpectationValue",
            [](const Type& op, const OneDM& D) {
                return op.calculateExpectationValue(D).elements();
            },
            "Calculate the expectation value of this one-electron operator.");


    // Expose the `VectorSpaceArithmetic` APIs.
    bindVectorSpaceArithmeticInterface(py_class);

    // Expose the `BasisTransformable` APIs.
    bindBasisTransformableInterface(py_class);

    // Expose the Mulliken-domain APIs.
    bindSQOneElectronOperatorMullikenInterface(py_class);
}


/**
 *  Add Python bindings for the `SimpleSQOneElectronOperator` C++ interface to a given `pybind11::class_`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain `SimpleSQOneElectronOperator`'s APIs.
 */
template <typename Class>
void bindSimpleSQOneElectronOperatorInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;


    py_class

        /*
         *  MARK: Parameter access
         */

        .def(
            "allParameters",
            [](const Type& op) {
                return op.allParameters();
            },
            "A vector of read-only matrix representions of the parameters/matrix elements/integrals of this operator.");


    // Expose general one-electron operator APIs.
    bindSQOneElectronOperatorInterface(py_class);
}


/**
 *  Add Python bindings for the parameter access for scalar one-electron operators.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain access bindings like a scalar one-electron operator.
 */
template <typename Class>
void bindScalarSQOneElectronOperatorParameterInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;

    py_class

        /*
         *  MARK: Constructors
         */

        .def(py::init<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>())


        .def(
            "parameters",
            [](const Type& op) {
                return op.parameters();
            },
            "A read-only matrix representation of the parameters/matrix elements/integrals of one of the tensor components of this operator.");
}


/**
 *  Add Python bindings for the some of the two-electron operator interface to a given `pybind11::class_`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain some two-electron operator APIs.
 */
template <typename Class>
void bindSQTwoElectronOperatorInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    using TwoDM = typename OperatorTraits<Type>::TwoDM;


    py_class

        /*
         *  MARK: Calculations
         */

        .def(
            "calculateExpectationValue",
            [](const Type& op, const TwoDM& D) {
                return op.calculateExpectationValue(D).elements();
            },
            "Calculate the expectation value of this one-electron operator.");


    // Expose the `VectorSpaceArithmetic` APIs.
    bindVectorSpaceArithmeticInterface(py_class);
}


/**
 *  Add Python bindings for the `SimpleSQTwoElectronOperator` C++ interface to a given `pybind11::class_`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain `SimpleSQTwoElectronOperator`'s APIs.
 */
template <typename Class>
void bindSimpleSQTwoElectronOperatorInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;


    // Add some two-electron operator APIs.
    bindSQTwoElectronOperatorInterface(py_class);

    // Expose the `BasisTransformable` APIs.
    bindBasisTransformableInterface(py_class);
}


/**
 *  Add Python bindings for the parameter access for scalar two-electron operators.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain access bindings like a scalar two-electron operator.
 */
template <typename Class>
void bindScalarSQTwoElectronOperatorParameterInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    using Scalar = typename Type::Scalar;


    py_class

        /*
         *  MARK: Constructors
         */

        .def(py::init<>([](const Eigen::Tensor<Scalar, 4>& g) {
                 return Type(SquareRankFourTensor<Scalar>(g));
             }),
             py::arg("g"))


        /*
         *  MARK: Parameter access
         */

        .def(
            "parameters",
            [](const Type& op) {
                return asNumpyArray(op.parameters().Eigen());
            },
            "A read-only matrix representation of the parameters/matrix elements/integrals of one of the tensor components of this operator.");
}


/*
 *  MARK: DensityMatrix
 */

/**
 *  Add Python bindings for some APIs related to `Simple1DM`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `Simple1DM`.
 */
template <typename Class>
void bindSimple1DMInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    using Scalar = typename Type::Scalar;


    // Add specific APIs for `Simple1DM`.
    py_class
        .def(py::init<>([](const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>& D) {
            return Type {D};
        }))

        .def(
            "matrix",
            [](const Type& D) {
                return D.matrix();
            },
            "Return a read-only reference to the matrix representation of the 1-DM.");


    // Add common APIs.
    bindVectorSpaceArithmeticInterface(py_class);
}


/**
 *  Add Python bindings for some APIs related to `Simple2DM`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `Simple2DM`.
 */
template <typename Class>
void bindSimple2DMInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;

    using Scalar = typename Type::Scalar;


    // Add specific APIs for `Simple2DM`.
    py_class
        .def(py::init<>([](const Eigen::Tensor<Scalar, 4>& d) {
            return Type(SquareRankFourTensor<Scalar>(d));
        }))

        .def(
            "tensor",
            [](const Type& d) {
                return asNumpyArray(d.tensor().Eigen());
            },
            "Return a read-only reference to the matrix representation of this 2-DM.")

        .def(
            "reduce",
            &Type::reduce,
            "Return a partial contraction of the 2-DM, where D(p,q) = d(p,q,r,r).")

        .def(
            "trace",
            &Type::trace,
            "Return the trace of the 2-DM, i.e. d(p,p,q,q).");


    // Add common APIs.
    bindVectorSpaceArithmeticInterface(py_class);
}


/*
 *  MARK: QCMethod - Hartree-Fock
 */

/**
 *  Add Python bindings for some APIs related to `QCMethod::R/U/GHF`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `QCMethod::R/U/GHF`.
 */
template <typename Class>
void bindQCMethodHartreeFockSCFEnvironmentInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;


    py_class

        /*
         *  MARK: Read-write members & properties
         */

        .def_readwrite("N", &Type::N)

        .def_readwrite("electronic_energies", &Type::electronic_energies)

        .def_readwrite("orbital_energies", &Type::orbital_energies)

        /*
         *  MARK: Read-only 'getters'
         */

        .def_readonly(
            "coefficient_matrices",
            &Type::coefficient_matrices)

        .def_readonly(
            "density_matrices",
            &Type::density_matrices)

        .def_readonly(
            "fock_matrices",
            &Type::fock_matrices)

        .def_readonly(
            "error_vectors",
            &Type::error_vectors)

        .def_readonly(
            "sq_hamiltonian",
            &Type::sq_hamiltonian);
}


/**
 *  Add Python bindings for some APIs related to R/U/GHF SCF solvers.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to R/U/GHF SCF solvers.
 */
template <typename Class>
void bindHartreeFockSCFSolverInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;


    py_class
        .def_static(
            "DIIS",
            &Type::DIIS,
            py::arg("minimum_subspace_dimension") = 6,
            py::arg("maximum_subspace_dimension") = 6,
            py::arg("threshold") = 1.0e-08,
            py::arg("maximum_number_of_iterations") = 128,
            "Return a DIIS HF SCF solver that uses the combination of norm of the difference of two consecutive density matrices as a convergence criterion.")

        .def_static(
            "Plain",
            &Type::Plain,
            py::arg("threshold") = 1.0e-08,
            py::arg("maximum_number_of_iterations") = 128,
            "Return a plain GHF SCF solver that uses the norm of the difference of two consecutive density matrices as a convergence criterion.");
}


/*
 *  MARK: QCModel - Hartree-Fock
 */

/**
 *  Add Python bindings for some APIs related to `QCModel::R/U/GHF`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `QCModel::R/U/GHF`.
 */
template <typename Class>
void bindQCModelHartreeFockInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;

    py_class
        .def(
            "orbitalSpace",
            [](const Type& parameters) {
                return parameters.orbitalSpace();
            },
            "Return the implicit (i.e. with ascending and contiguous orbital indices) occupied-virtual orbital space that corresponds to these RHF model parameters.")

        .def(
            "calculateOrthonormalBasis1DM",
            [](const Type& parameters) {
                return parameters.calculateOrthonormalBasis1DM();
            },
            "Return the 1-DM expressed in an orthonormal basis related to these optimal HF parameters.")

        .def(
            "calculateOrthonormalBasis2DM",
            [](const Type& parameters) {
                return parameters.calculateOrthonormalBasis2DM();
            },
            "Return the 2-DM expressed in an orthonormal basis related to these optimal HF parameters.")

        .def(
            "calculateScalarBasis1DM",
            [](const Type& parameters) {
                return parameters.calculateScalarBasis1DM();
            },
            "Return the HF 1-DM in the scalar/AO basis related to these optimal parameters.")

        .def(
            "calculateScalarBasis2DM",
            [](const Type& parameters) {
                return parameters.calculateScalarBasis2DM();
            },
            "Return the HF 2-DM in the scalar/AO basis related to these optimal parameters.")

        .def(
            "calculateStabilityMatrices",
            &Type::calculateStabilityMatrices,
            py::arg("hamiltonian"),
            "Return the HF stability matrices parameters.")

        .def("expansion",
             &Type::expansion,
             "Return the coefficient matrix that expresses every spatial orbital (as a column) in its underlying scalar basis.")

        .def(
            "numberOfElectrons",
            [](const Type& parameters) {
                return parameters.numberOfElectrons();
            },
            "Return the number of electrons that these HF model parameters describe, i.e. the number of occupied orbitals.")

        .def("occupiedOrbitalEnergies",
             &Type::occupiedOrbitalEnergies,
             "Return the orbital energies of the occupied orbitals.")

        .def("orbitalEnergies",
             &Type::orbitalEnergies,
             "Return all of the orbital energies.")

        .def("virtualOrbitalEnergies",
             &Type::virtualOrbitalEnergies,
             "Return the orbital energies of the virtual orbitals.");
}


/**
 *  Add Python bindings for some APIs related to `complex R/U/GHF Stability matrices`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `R/U/GHFStabilityMatrices`.
 */
template <typename Class>
void bindQCModelHartreeFockComplexStabilityInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;

    py_class
        .def(
            "internal",
            [](const Type& stability_matrices) {
                return stability_matrices.internal();
            },
            "Return the internal HF stability matrix.")

        .def("isInternallyStable",
             &Type::isInternallyStable,
             py::arg("threshold") = -1.0e-05,
             "Return a boolean, telling us if the real or complex valued internal stability matrix belongs to a stable or unstable set of parameters.")

        .def(
            "printStabilityDescription",
            [](const Type& stability_matrices) {
                py::scoped_ostream_redirect stream {std::cout, py::module::import("sys").attr("stdout")};
                stability_matrices.printStabilityDescription();
            },
            "Print the description of the stability properties of a HF wavefunction.");
}


/**
 *  Add Python bindings for some APIs related to `Real R/U/GHF Stability matrices`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `R/U/GHFStabilityMatrices`.
 */
template <typename Class>
void bindQCModelHartreeFockRealStabilityInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;

    py_class
        .def(
            "realComplex",
            [](const Type& stability_matrices) {
                return stability_matrices.realComplex();
            },
            "Return the real->complex external HF stability matrix.")

        .def(
            "isExternallyStable",
            [](const Type& stability_matrices, const double threshold) {
                return stability_matrices.isExternallyStable(threshold);
            },
            py::arg("threshold") = -1.0e-05,
            "Return a boolean, telling us if the real valued external stability matrix belongs to a stable or unstable set of parameters.");

    // Add the APIs from the complex interface, as every real method has those as well.
    bindQCModelHartreeFockComplexStabilityInterface(py_class);
}


/**
 *  Add Python bindings for some APIs related to `R/U/GNonOrthogonalStateBases`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `R/U/GNonOrthogonalStateBases`.
 */
template <typename Class>
void bindNonOrthogonalStateBasisInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;

    py_class

        /**
         * MARK: Access
         */

        .def("basisState",
             &Type::basisState,
             py::arg("index"),
             "Return the i'th basis state in this non-orthogonal basis.")

        .def("basisStates",
             &Type::basisStates,
             "Return the number of basis states in this non-orthogonal basis.")

        .def("basisStateDimension",
             &Type::basisStateDimension,
             "Return the dimension of the basis states in this non-orthogonal basis.")

        .def("metric",
             &Type::metric,
             "Return the overlap metric of this non-orthogonal basis as an operator.")

        .def("numberOfBasisStates",
             &Type::numberOfBasisStates,
             "Return the number of basis states in this non-orthogonal basis.")

        .def("numberOfOccupiedOrbitals",
             &Type::numberOfOccupiedOrbitals,
             "Return the number of orbitals in each state that should be considered occupied.")

        .def("threshold",
             &Type::threshold,
             "Return the threshold against which values are compared in oreder to check whether they are zero or not.")

        /**
         * MARK: Operator evaluations
         */

        .def("evaluateHamiltonianOperator",
             &Type::evaluateHamiltonianOperator,
             py::arg("hamiltonian"),
             py::arg("nuclear_repulsion_operator"),
             "Return the matrix representation of any second quantized hamiltonian evaluated in this non-orthogonal basis.")

        .def("evaluateOneElectronOperator",
             &Type::evaluateOneElectronOperator,
             py::arg("one_electron_operator"),
             "Return the matrix representation of any second quantized one-electron operator evaluated in this non-orthogonal basis.")

        .def("evaluateOverlapOperator",
             &Type::evaluateOverlapOperator,
             "Return the matrix representation of the overlap metric evaluated in this non-orthogonal basis.")

        .def("evaluateTwoElectronOperator",
             &Type::evaluateTwoElectronOperator,
             py::arg("two_electron_operator"),
             "Return the matrix representation of any second quantized two-electron operator evaluated in this non-orthogonal basis.");

    /**
     *  MARK: Conforming to `BasisTransformable`
     */

    bindBasisTransformableInterface(py_class);
}


/**
 *  Add Python bindings for some APIs related to `R/U/GLowdinPairingBasis`.
 *
 *  @tparam Class               The type of the Pybind11 `class_` (generated by the compiler).
 *
 *  @param py_class             The Pybind11 `class_` that should obtain APIs related to `R/U/GLowdinPairingBasis`.
 */
template <typename Class>
void bindLowdinPairingBasisInterface(Class& py_class) {

    // The C++ type corresponding to the Python class.
    using Type = typename Class::type;
    using Scalar = typename Type::Scalar;

    py_class

        /**
         * MARK: Access
         */

        .def(
            "biorthogonalBraExpansion",
            [](const Type& lp_basis) {
                return lp_basis.biorthogonalBraExpansion();
            },
            "Return the biorthogonalized expansion coefficients of the bra.")

        .def(
            "biorthogonalKetExpansion",
            [](const Type& lp_basis) {
                return lp_basis.biorthogonalKetExpansion();
            },
            "Return the biorthogonalized expansion coefficients of the ket.")

        /**
         * MARK: Overlap metrics
         */

        .def(
            "biorthogonalOverlaps",
            [](const Type& lp_basis) {
                return lp_basis.biorthogonalOverlaps();
            },
            "Return the overlap values of the biorthogonalized expansion coefficients.")

        .def(
            "numberOfZeroOverlaps",
            [](const Type& lp_basis) {
                return lp_basis.numberOfZeroOverlaps();
            },
            "Determine the number of zero overlaps in the biorthogonal overlap vector.")

        .def(
            "reducedOverlap",
            [](const Type& lp_basis) {
                return lp_basis.reducedOverlap();
            },
            "Calculate the reduced overlap. I.e. the biorthogonal overlaps with the zero values get removed and the total remaining overlap is calculated.")


        .def(
            "totalOverlap",
            [](const Type& lp_basis) {
                return lp_basis.totalOverlap();
            },
            "Calculate and return the total overlap value of the biorthogonal coefficients.")


        .def(
            "zeroOverlapIndices",
            [](const Type& lp_basis) {
                return lp_basis.zeroOverlapIndices();
            },
            "Determine the indices of the zero overlap values in the biorthogonal overlap vector.")


        /**
         * MARK: Density matrices
         */

        .def("coDensity",
             &Type::coDensity,
             py::arg("k"),
             "Calculate and return the co-density matrix for the given occupied orbital index.")

        .def("coDensitySum",
             &Type::coDensitySum,
             "Calculate and return the sum of the zero overlap co-density matrix, the weighted co-density matrix and the zero overlap co-density matrix of the biorthogonalized basis of the bra with itself.")

        .def("transition1DM",
             &Type::transition1DM,
             "Calculate and return the transition one-electron density matrix.")

        .def(
            "weightedCoDensity",
            [](const Type& lp_basis) {
                return lp_basis.weightedCoDensity();
            },
            "Calculate and return the weighted co-density matrix. It calculates the co-density matrix contributions corresponding to the non-zero overlap values and divides them by that overlap value.")

        .def("zeroOverlapCoDensity",
             &Type::zeroOverlapCoDensity,
             "Calculate and return the zero-overlap co-density matrix. It takes the indices of the zero overlaps and calculates the co-density matrix at each of these indices. The sum of those co-density matrices is returned.");
}


}  // namespace gqcpy
