SimpleTransformation.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/Transformations/BasisTransformable.hpp"
#include "Basis/Transformations/JacobiRotatable.hpp"
#include "Basis/Transformations/OrbitalRotationGeneratorTraits.hpp"
#include "Mathematical/Representation/SquareMatrix.hpp"
 
#include <unsupported/Eigen/MatrixFunctions>
 
 
namespace GQCP {
 
 
/*
 *  MARK: SimpleTransformation implementation
 */
 
template <typename _Scalar, typename _DerivedTransformation>
class SimpleTransformation:
    public BasisTransformable<_DerivedTransformation>,
    public JacobiRotatable<_DerivedTransformation> {
 
public:
    // The scalar type used for a transformation coefficient: real or complex.
    using Scalar = _Scalar;
 
    // The type of the transformation that derives from this class, enabling CRTP and compile-time polymorphism.
    using DerivedTransformation = _DerivedTransformation;
 
    // The type of 'this'.
    using Self = SimpleTransformation<Scalar, DerivedTransformation>;
 
    // The type of Jacobi rotation for which the Jacobi rotation should be defined.
    using JacobiRotationType = JacobiRotation;
 
    // The type of the orbital rotation generators that is naturally associated with the derived transformation.
    using OrbitalRotationGeneratorType = typename OrbitalRotationGeneratorTraits<DerivedTransformation>::OrbitalRotationGenerators;
 
 
protected:
    // The transformation matrix that collects the expansion coefficients of the new basis (vectors) in the old basis as columns.
    SquareMatrix<Scalar> T;
 
 
public:
    /*
     *  MARK: Constructors
     */
 
    SimpleTransformation(const SquareMatrix<Scalar>& T) :
        T {T} {}
 
 
    SimpleTransformation(const OrbitalRotationGeneratorType& kappa) :
        T {(-kappa.asMatrix()).exp()} {}
 
 
    /*
     *  MARK: Named constructors
     */
 
    static DerivedTransformation FromJacobi(const JacobiRotation& jacobi_rotation, const size_t dim) {
 
        // Create an identity transformation matrix and apply a Jacobi rotation.
        DerivedTransformation J = SquareMatrix<Scalar>::Identity(dim);
 
        return J.rotated(jacobi_rotation);
    }
 
 
    static DerivedTransformation Identity(const size_t dim) { return DerivedTransformation {SquareMatrix<Scalar>::Identity(dim)}; }
 
    static DerivedTransformation Random(const size_t dim) { return DerivedTransformation {SquareMatrix<Scalar>::Random(dim)}; }
 
    static DerivedTransformation RandomUnitary(const size_t dim) { return DerivedTransformation {SquareMatrix<Scalar>::RandomUnitary(dim)}; }
 
    static DerivedTransformation Zero(const size_t dim) { return DerivedTransformation {SquareMatrix<Scalar>::Zero(dim)}; }
 
 
    /*
     *  MARK: General information
     */
 
    size_t numberOfOrbitals() const { return this->matrix().dimension(); }
 
    size_t dimension() const { return this->numberOfOrbitals(); }
 
 
    /*
     *  MARK: Transformation matrix
     */
 
    const SquareMatrix<Scalar>& matrix() const { return this->T; }
 
 
    /*
     *  MARK: Linear algebra
     */
 
    DerivedTransformation adjoint() const { return DerivedTransformation {this->matrix().adjoint()}; }
 
    DerivedTransformation inverse() const { return DerivedTransformation {this->matrix().inverse()}; }
 
    bool isUnitary(const double threshold = 1.0e-12) const { return this->matrix().isUnitary(threshold); }
 
 
    /*
     *  MARK: Conforming to `BasisTransformable`
     */
 
    DerivedTransformation transformed(const DerivedTransformation& T) const override { return DerivedTransformation {this->matrix() * T.matrix()}; }
 
    // Allow the `rotate` method from `BasisTransformable`, since there's also a `rotate` from `JacobiRotatable`.
    using BasisTransformable<DerivedTransformation>::rotate;
 
    // Allow the `rotated` method from `BasisTransformable`, since there's also a `rotated` from `JacobiRotatable`.
    using BasisTransformable<DerivedTransformation>::rotated;
 
 
    /*
     *  MARK: Conforming to JacobiRotatable
     */
 
    DerivedTransformation rotated(const JacobiRotationType& jacobi_rotation) const override {
 
        const auto p = jacobi_rotation.p();
        const auto q = jacobi_rotation.q();
 
        // Create the matrix representation of a Jacobi rotation using Eigen's APIs.
        // We're applying the Jacobi rotation as J = I * jacobi_rotation (cfr. B' = B T).
        const auto eigen_jacobi_rotation = jacobi_rotation.Eigen();
        auto result = this->matrix();
        result.applyOnTheRight(p, q, eigen_jacobi_rotation);
 
        return DerivedTransformation {result};
    }
 
    // Allow the `rotate` method from `JacobiRotatable`, since there's also a `rotate` from `BasisTransformable`.
    using JacobiRotatable<DerivedTransformation>::rotate;
};
 
 
/*
 *  MARK: BasisTransformableTraits
 */
 
template <typename Scalar, typename DerivedTransformation>
struct BasisTransformableTraits<SimpleTransformation<Scalar, DerivedTransformation>> {
 
    // The type of the transformation for which the basis transformation should be defined. A transformation matrix should naturally be transformable with itself.
    using Transformation = DerivedTransformation;
};
 
 
/*
 *  MARK: JacobiRotatableTraits
 */
 
template <typename Scalar, typename DerivedTransformation>
struct JacobiRotatableTraits<SimpleTransformation<Scalar, DerivedTransformation>> {
 
    // The type of Jacobi rotation for which the Jacobi rotation should be defined.
    using JacobiRotationType = JacobiRotation;
};
 
 
}  // namespace GQCP