FLEX.h

/* Copyright 2018 Kristofer Björnson
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */


#ifndef COM_DAFER45_TBTK_SOLVER_FLEX
#define COM_DAFER45_TBTK_SOLVER_FLEX

#include "TBTK/MomentumSpaceContext.h"
#include "TBTK/Property/GreensFunction.h"
#include "TBTK/Property/Susceptibility.h"
#include "TBTK/Property/InteractionVertex.h"
#include "TBTK/Property/SelfEnergy.h"
#include "TBTK/Solver/Solver.h"

namespace TBTK{
namespace Solver{

class FLEX : public Solver{
public:
    enum class State {
        NotYetStarted,
        GreensFunctionCalculated,
        BareSusceptibilityCalculated,
        RPASusceptibilitiesCalculated,
        InteractionVertexCalculated,
        SelfEnergyCalculated
    };

    enum class Norm{Max, L2};

    FLEX(const MomentumSpaceContext &momentumSpaceContext);

    const MomentumSpaceContext& getMomentumSpaceContext() const;

    const Property::GreensFunction& getGreensFunction() const;

    const Property::Susceptibility& getBareSusceptibility() const;

    const Property::Susceptibility& getRPASpinSusceptibility() const;

    const Property::Susceptibility& getRPAChargeSusceptibility() const;

    const Property::InteractionVertex& getInteractionVertex() const;

    void setSelfEnergy(const Property::SelfEnergy &selfEnergy);

    const Property::SelfEnergy& getSelfEnergy() const;

    void setSelfEnergyMixingParameter(double selfEnergyMixingParameter);

    void setEnergyWindow(
        int lowerFermionicMatsubaraEnergyIndex,
        int upperFermionicMatsubaraEnergyIndex,
        int lowerBosonicMatsubaraEnergyIndex,
        int upperBosonicMatsubaraEnergyIndex
    );

    void setU(double U);

    void setJ(double J);

    void setUp(double Up);

    void setJp(double Jp);

    void setNumOrbitals(unsigned int numOrbitals);

    double getDensity() const;

    void setTargetDensity(double targetDensity);

    void setDensityTolerance(double densityTolerance);

    void setChemicalPotentialStepSize(double chemicalPotentialStepSize);

    State getState() const;

    void setMaxIterations(unsigned int maxIterations);

    void setNorm(Norm norm);

    void setTolerance(double tolerance);

    void setCallback(void (*callback)(FLEX &solver));

    void setNumSlices(unsigned int numSlices);

    void run();
private:
    const MomentumSpaceContext momentumSpaceContext;

    Property::GreensFunction greensFunction0;

    Property::GreensFunction greensFunction;

    Property::GreensFunction oldGreensFunction;

    Property::Susceptibility bareSusceptibility;

    Property::Susceptibility rpaSpinSusceptibility;

    Property::Susceptibility rpaChargeSusceptibility;

    Property::InteractionVertex interactionVertex;

    Property::SelfEnergy selfEnergy;

    Property::SelfEnergy previousSelfEnergy;

    double selfEnergyMixingParameter;

    double density;

    double targetDensity;

    double densityTolerance;

    double chemicalPotentialStepSize;

    int lowerFermionicMatsubaraEnergyIndex;

    int upperFermionicMatsubaraEnergyIndex;

    int lowerBosonicMatsubaraEnergyIndex;

    int upperBosonicMatsubaraEnergyIndex;

    double U;

    double J;

    double Up;

    double Jp;

    unsigned int numOrbitals;

    State state;

    unsigned int maxIterations;

    Norm norm;

    double tolerance;

    unsigned int numSlices;

    double convergenceParameter;

    void (*callback)(FLEX &solver);

    void calculateBareGreensFunction();

    void calculateBareSusceptibility(unsigned int slice);

    void calculateRPASusceptibilities();

    void calculateInteractionVertex();

    void calculateSelfEnergy(unsigned int slice);

    void calculateGreensFunction();

    void convertSelfEnergyIndexStructure();

    void calculateConvergenceParameter();

    std::vector<InteractionAmplitude>
        generateRPAChargeSusceptibilityInteractionAmplitudes();

    std::vector<InteractionAmplitude>
        generateRPASpinSusceptibilityInteractionAmplitudes();

    std::vector<InteractionAmplitude> generateU1();

    std::vector<InteractionAmplitude> generateU2();

    std::vector<InteractionAmplitude> generateU3();

    std::vector<InteractionAmplitude> generateU4();

    Property::InteractionVertex generateHartreeFockTerm();

    void calculateDensity();

    int getLowerBosonicMatsubaraEnergyIndex(unsigned int slice);

    int getUpperBosonicMatsubaraEnergyIndex(unsigned int slice);
};

inline const MomentumSpaceContext& FLEX::getMomentumSpaceContext() const{
    return momentumSpaceContext;
}

inline const Property::GreensFunction& FLEX::getGreensFunction() const{
    return greensFunction;
}

inline const Property::Susceptibility& FLEX::getBareSusceptibility() const{
    return bareSusceptibility;
}

inline const Property::Susceptibility& FLEX::getRPASpinSusceptibility() const{
    return rpaSpinSusceptibility;
}

inline const Property::Susceptibility& FLEX::getRPAChargeSusceptibility(
) const{
    return rpaChargeSusceptibility;
}

inline const Property::InteractionVertex& FLEX::getInteractionVertex() const{
    return interactionVertex;
}

inline void FLEX::setSelfEnergy(const Property::SelfEnergy &selfEnergy){
    this->selfEnergy = selfEnergy;
}

inline const Property::SelfEnergy& FLEX::getSelfEnergy() const{
    return selfEnergy;
}

inline void FLEX::setSelfEnergyMixingParameter(
    double selfEnergyMixingParameter
){
    this->selfEnergyMixingParameter = selfEnergyMixingParameter;
}

inline double FLEX::getDensity() const{
    TBTKAssert(
        targetDensity >= 0,
        "Solver::FLEX::getDensity()",
        "This function is only available for non-negative target"
        << " densities.",
        "Use Solver::FLEX::setTargetDensity() to set the target"
        << " density."
    );

    return density;
}

inline void FLEX::setTargetDensity(double targetDensity){
    this->targetDensity = targetDensity;
}

inline void FLEX::setDensityTolerance(double densityTolerance){
    this->densityTolerance = densityTolerance;
}

inline void FLEX::setChemicalPotentialStepSize(
    double chemicalPotentialStepSize
){
    this->chemicalPotentialStepSize = chemicalPotentialStepSize;
}

inline void FLEX::setEnergyWindow(
    int lowerFermionicMatsubaraEnergyIndex,
    int upperFermionicMatsubaraEnergyIndex,
    int lowerBosonicMatsubaraEnergyIndex,
    int upperBosonicMatsubaraEnergyIndex
){
    TBTKAssert(
        lowerFermionicMatsubaraEnergyIndex
            <= upperFermionicMatsubaraEnergyIndex,
        "Solver::FLEX::setEnergyWindow()",
        "The 'lowerFermionicMatsubaraEnergyIndex="
        << lowerFermionicMatsubaraEnergyIndex << "' must be less or"
        << " equal to the 'upperFermionicMatsubaraEnergyIndex="
        << upperFermionicMatsubaraEnergyIndex << "'.",
        ""
    );
    TBTKAssert(
        lowerBosonicMatsubaraEnergyIndex
            <= upperBosonicMatsubaraEnergyIndex,
        "Solver::FLEX::setEnergyWindow()",
        "The 'lowerBosonicMatsubaraEnergyIndex="
        << lowerBosonicMatsubaraEnergyIndex << "' must be less or"
        << " equal to the 'upperBosonicMatsubaraEnergyIndex="
        << upperBosonicMatsubaraEnergyIndex << "'.",
        ""
    );

    this->lowerFermionicMatsubaraEnergyIndex
        = lowerFermionicMatsubaraEnergyIndex;
    this->upperFermionicMatsubaraEnergyIndex
        = upperFermionicMatsubaraEnergyIndex;
    this->lowerBosonicMatsubaraEnergyIndex
        = lowerBosonicMatsubaraEnergyIndex;
    this->upperBosonicMatsubaraEnergyIndex
        = upperBosonicMatsubaraEnergyIndex;
}

inline void FLEX::setU(double U){
    this->U = U;
}

inline void FLEX::setJ(double J){
    this->J = J;
}

inline void FLEX::setUp(double Up){
    this->Up = Up;
}

inline void FLEX::setJp(double Jp){
    this->Jp = Jp;
}

inline void FLEX::setNumOrbitals(unsigned int numOrbitals){
    this->numOrbitals = numOrbitals;
}

inline FLEX::State FLEX::getState() const{
    return state;
}

inline void FLEX::setMaxIterations(unsigned int maxIterations){
    TBTKAssert(
        maxIterations > 0,
        "Solver::FLEX::setMaxIterations()",
        "'maxIterations=" << maxIterations << "' must be larger than"
        << " zero.",
        ""
    );

    this->maxIterations = maxIterations;
}

inline void FLEX::setNorm(Norm norm){
    this->norm = norm;
}

inline void FLEX::setTolerance(double tolerance){
    this->tolerance = tolerance;
}

inline void FLEX::setCallback(void (*callback)(FLEX &solver)){
    this->callback = callback;
}

inline void FLEX::setNumSlices(unsigned int numSlices){
    TBTKAssert(
        numSlices != 0,
        "Solver::FLEX::setNumSlices()",
        "'numSlices' must be non-zero.",
        ""
    );

    this->numSlices = numSlices;
}

};  //End of namespace Solver
};  //End of namespace TBTK

#endif