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