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