RPASusceptibilityCalculator.h

/* Copyright 2017 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_RPA_SUSCEPTIBILITY_CALCULATOR
#define COM_DAFER45_TBTK_RPA_SUSCEPTIBILITY_CALCULATOR

#include "TBTK/InteractionAmplitude.h"
#include "TBTK/RPA/MomentumSpaceContext.h"
#include "TBTK/SerializableVector.h"
#include "TBTK/RPA/SusceptibilityCalculator.h"
#include "TBTK/RPA/LindhardSusceptibilityCalculator.h"

#include <complex>

//#include <omp.h>

namespace TBTK{

class RPASusceptibilityCalculator{
public:
    RPASusceptibilityCalculator(
        const RPA::MomentumSpaceContext &momentumSpaceContext,
        SusceptibilityCalculator::Algorithm algorithm
            = SusceptibilityCalculator::Algorithm::Lindhard
    );

    ~RPASusceptibilityCalculator();

    RPASusceptibilityCalculator* createSlave();

    void precompute(unsigned int numWorkers = 129);

    const RPA::MomentumSpaceContext& getMomentumSpaceContext() const;

    void setInteractionAmplitudes(
        const std::vector<InteractionAmplitude> &interactionAmplitudes
    );

    enum class EnergyType {Real, Imaginary, Complex};

//  void setSusceptibilityMode(SusceptibilityCalculator::Mode mode);

//  SusceptibilityCalculator::Mode getSusceptibilityMode() const;

    void setEnergyType(EnergyType energyType);

    EnergyType getEnergyType() const;

    void setEnergies(
        const std::vector<std::complex<double>> &energies
    );

    void setEnergiesAreInversionSymmetric(
        bool energiesAreInversionSymmetric
    );

    bool getEnergiesAreInversionSymmetric() const;

    void setSusceptibilityIsSafeFromPoles(
        bool susceptibilityIsSafeFromPoles
    );

    bool getSusceptibilityIsSafeFromPoles() const;

//  void setNumSummationEnergies(unsigned int numSummationEnergies);

    void saveSusceptibilities(const std::string &filename) const;

    void loadSusceptibilities(const std::string &filename);
private:
    IndexedDataTree<SerializableVector<std::complex<double>>> rpaSusceptibilityTree;

    IndexedDataTree<SerializableVector<std::complex<double>>> rpaChargeSusceptibilityTree;

    IndexedDataTree<SerializableVector<std::complex<double>>> rpaSpinSusceptibilityTree;

    EnergyType energyType;

    std::vector<std::complex<double>> energies;
public:
    std::vector<std::complex<double>> calculateRPASusceptibility(
        const DualIndex &kDual,
        const std::vector<int> &orbitalIndices
    );

    std::vector<std::complex<double>> calculateRPASusceptibility(
        const std::vector<double> &k,
        const std::vector<int> &orbitalIndices
    );

    std::vector<std::complex<double>> calculateChargeRPASusceptibility(
        const DualIndex &kDual,
        const std::vector<int> &orbitalIndices
    );

    std::vector<std::complex<double>> calculateChargeRPASusceptibility(
        const std::vector<double> &k,
        const std::vector<int> &orbitalIndices
    );

    std::vector<std::complex<double>> calculateSpinRPASusceptibility(
        const DualIndex &kDual,
        const std::vector<int> &orbitalIndices
    );

    std::vector<std::complex<double>> calculateSpinRPASusceptibility(
        const std::vector<double> &k,
        const std::vector<int> &orbitalIndices
    );

    void setU(std::complex<double> U);

    void setUp(std::complex<double> Up);

    void setJ(std::complex<double> J);

    void setJp(std::complex<double> Jp);
private:
    SusceptibilityCalculator *susceptibilityCalculator;

    std::vector<InteractionAmplitude> interactionAmplitudes;

    std::vector<InteractionAmplitude> interactionAmplitudesCharge;

    std::vector<InteractionAmplitude> interactionAmplitudesSpin;

    bool interactionAmplitudesAreGenerated;

    RPASusceptibilityCalculator(
        SusceptibilityCalculator &susceptibilityCalculator
    );

    Index getSusceptibilityResultIndex(
        const Index &kIndex,
        const std::vector<int> &orbitalIndices
    ) const;

    void invertMatrix(std::complex<double> *matrix, unsigned int dimensions);

    void multiplyMatrices(
        std::complex<double> *matrix1,
        std::complex<double> *matrix2,
        std::complex<double> *result,
        unsigned int dimensions
    );

    std::vector<std::vector<std::vector<std::complex<double>>>> rpaSusceptibilityMainAlgorithm(
        const DualIndex &kDual,
        const std::vector<int> &orbitalIndices,
        const std::vector<InteractionAmplitude> &interactionAmpltiudes
    );

    std::complex<double> U, Up, J, Jp;

    void generateInteractionAmplitudes();
};

inline const RPA::MomentumSpaceContext& RPASusceptibilityCalculator::getMomentumSpaceContext(
) const{
    return susceptibilityCalculator->getMomentumSpaceContext();
}

inline void RPASusceptibilityCalculator::setInteractionAmplitudes(
    const std::vector<InteractionAmplitude> &interactionAmplitudes
){
    this->interactionAmplitudes = interactionAmplitudes;
}

inline Index RPASusceptibilityCalculator::getSusceptibilityResultIndex(
    const Index &kIndex,
    const std::vector<int> &orbitalIndices
) const{
    return Index(
        kIndex,
        {
            orbitalIndices.at(0),
            orbitalIndices.at(1),
            orbitalIndices.at(2),
            orbitalIndices.at(3)
        }
    );
}

/*inline void RPASusceptibilityCalculator::setSusceptibilityMode(
    SusceptibilityCalculator::Mode mode
){
    susceptibilityCalculator->setMode(mode);
}

inline SusceptibilityCalculator::Mode RPASusceptibilityCalculator::getSusceptibilityMode(
) const{
    return susceptibilityCalculator->getMode();
}*/

inline void RPASusceptibilityCalculator::setEnergyType(
    EnergyType energyType
){
    this->energyType = energyType;
    switch(energyType){
    case EnergyType::Real:
        susceptibilityCalculator->setEnergyType(
            SusceptibilityCalculator::EnergyType::Real
        );
        break;
    case EnergyType::Imaginary:
        susceptibilityCalculator->setEnergyType(
            SusceptibilityCalculator::EnergyType::Imaginary
        );
        break;
    case EnergyType::Complex:
        susceptibilityCalculator->setEnergyType(
            SusceptibilityCalculator::EnergyType::Complex
        );
        break;
    default:
        TBTKExit(
            "RPASusceptibilityCalculator::setEnergyType()",
            "Unknown energy type.",
            "This should never happen, contact the developer."
        );
    }
}

inline RPASusceptibilityCalculator::EnergyType RPASusceptibilityCalculator::getEnergyType(
) const{
    return energyType;
}

inline void RPASusceptibilityCalculator::setEnergies(
    const std::vector<std::complex<double>> &energies
){
    susceptibilityCalculator->setEnergies(
        energies
    );

    this->energies = energies;

    rpaSusceptibilityTree.clear();
    rpaChargeSusceptibilityTree.clear();
    rpaSpinSusceptibilityTree.clear();
}

inline void RPASusceptibilityCalculator::setEnergiesAreInversionSymmetric(
    bool energiesAreInversionSymmetric
){
    susceptibilityCalculator->setEnergiesAreInversionSymmetric(
        energiesAreInversionSymmetric
    );
}

inline bool RPASusceptibilityCalculator::getEnergiesAreInversionSymmetric(
) const{
    return susceptibilityCalculator->getEnergiesAreInversionSymmetric();
}

inline void RPASusceptibilityCalculator::setSusceptibilityIsSafeFromPoles(
    bool susceptibilityIsSafeFromPoles
){
    TBTKAssert(
        susceptibilityCalculator->getAlgorithm()
            == SusceptibilityCalculator::Algorithm::Lindhard,
        "RPASusceptibilityCalculator::setSusceptibilityIsSafeFromPoles()",
        "Only valid function call if the underlying algorithm is"
        << " SusceptibilityCalculator::Algorithm::Lindhard.",
        ""
    );

    ((LindhardSusceptibilityCalculator*)susceptibilityCalculator)->setSusceptibilityIsSafeFromPoles(
        susceptibilityIsSafeFromPoles
    );
}

inline bool RPASusceptibilityCalculator::getSusceptibilityIsSafeFromPoles() const{
    TBTKAssert(
        susceptibilityCalculator->getAlgorithm()
            == SusceptibilityCalculator::Algorithm::Lindhard,
        "RPASusceptibilityCalculator::getSusceptibilityIsSafeFromPoles()",
        "Only valid function call if the underlying algorithm is"
        << " SusceptibilityCalculator::Algorithm::Lindhard.",
        ""
    );

    return ((LindhardSusceptibilityCalculator*)susceptibilityCalculator)->getSusceptibilityIsSafeFromPoles();
}

/*inline void RPASusceptibilityCalculator::setNumSummationEnergies(
    unsigned int numSummationEnergies
){
    susceptibilityCalculator->setNumSummationEnergies(
        numSummationEnergies
    );
}*/

inline void RPASusceptibilityCalculator::saveSusceptibilities(
    const std::string &filename
) const{
    susceptibilityCalculator->saveSusceptibilities(filename);
}

inline void RPASusceptibilityCalculator::loadSusceptibilities(
    const std::string &filename
){
    susceptibilityCalculator->loadSusceptibilities(filename);
}

inline std::vector<std::complex<double>> RPASusceptibilityCalculator::calculateRPASusceptibility(
        const std::vector<double> &k,
        const std::vector<int> &orbitalIndices
){
    return calculateRPASusceptibility(
        DualIndex(
            susceptibilityCalculator->getMomentumSpaceContext().getKIndex(k),
            k
        ),
        orbitalIndices
    );
}

inline std::vector<std::complex<double>> RPASusceptibilityCalculator::calculateChargeRPASusceptibility(
        const std::vector<double> &k,
        const std::vector<int> &orbitalIndices
){
    return calculateChargeRPASusceptibility(
        DualIndex(
            susceptibilityCalculator->getMomentumSpaceContext().getKIndex(k),
            k
        ),
        orbitalIndices
    );
}

inline std::vector<std::complex<double>> RPASusceptibilityCalculator::calculateSpinRPASusceptibility(
        const std::vector<double> &k,
        const std::vector<int> &orbitalIndices
){
    return calculateSpinRPASusceptibility(
        DualIndex(
            susceptibilityCalculator->getMomentumSpaceContext().getKIndex(k),
            k
        ),
        orbitalIndices
    );
}

inline void RPASusceptibilityCalculator::setU(std::complex<double> U){
    this->U = U;
    rpaSusceptibilityTree.clear();
    rpaChargeSusceptibilityTree.clear();
    rpaSpinSusceptibilityTree.clear();
    interactionAmplitudesAreGenerated = false;
}

inline void RPASusceptibilityCalculator::setUp(std::complex<double> Up){
    this->Up = Up;
    rpaSusceptibilityTree.clear();
    rpaChargeSusceptibilityTree.clear();
    rpaSpinSusceptibilityTree.clear();
    interactionAmplitudesAreGenerated = false;
}

inline void RPASusceptibilityCalculator::setJ(std::complex<double> J){
    this->J = J;
    rpaSusceptibilityTree.clear();
    rpaChargeSusceptibilityTree.clear();
    rpaSpinSusceptibilityTree.clear();
    interactionAmplitudesAreGenerated = false;
}

inline void RPASusceptibilityCalculator::setJp(std::complex<double> Jp){
    this->Jp = Jp;
    rpaSusceptibilityTree.clear();
    rpaChargeSusceptibilityTree.clear();
    rpaSpinSusceptibilityTree.clear();
    interactionAmplitudesAreGenerated = false;
}

};  //End of namespace TBTK

#endif