charge_mixing.h

Go to the documentation of this file.

#ifndef CHARGE_MIXING_H
#define CHARGE_MIXING_H
#include "charge.h"
#include "source_estate/module_dm/density_matrix.h"
#include "source_base/module_mixing/mixing.h"
#include "source_base/module_mixing/plain_mixing.h"
 
class Charge_Mixing
{
 
  public:
    Charge_Mixing();
    ~Charge_Mixing();
 
    void set_mixing(const std::string& mixing_mode_in,
                    const double& mixing_beta_in,
                    const int& mixing_ndim_in,
                    const double& mixing_gg0_in,
                    const bool& mixing_tau_in,
                    const double& mixing_beta_mag_in,
                    const double& mixing_gg0_mag_in,
                    const double& mixing_gg0_min_in,
                    const double& mixing_angle_in,
                    const bool& mixing_dmr_in,
                    double& omega_in,
                    double& tpiba_in);
 
    void close_kerker_gg0() { mixing_gg0 = 0.0; mixing_gg0_mag = 0.0; }
    void init_mixing();
 
    void allocate_mixing_dmr(const int nnr);
 
    void mix_rho(Charge* chr);
 
    void mix_dmr(elecstate::DensityMatrix<double, double>* DM);
    void mix_dmr(elecstate::DensityMatrix<std::complex<double>, double>* DM);
    
    double get_drho(Charge* chr, const double nelec);
    double get_dkin(Charge* chr, const double nelec);
 
    void mix_reset();
    
    void set_rhopw(ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis* rhodpw_in);
 
    // extracting parameters normally these parameters will not be used outside charge mixing
    // while Exx is using them as well as some other places
    const std::string& get_mixing_mode() const {return mixing_mode;}
    double get_mixing_beta() const {return mixing_beta;}
    int get_mixing_ndim() const {return mixing_ndim;}
    double get_mixing_gg0() const {return mixing_gg0;}
    Base_Mixing::Mixing* get_mixing() const {return mixing;}
 
    // for mixing restart
    int mixing_restart_step = 0; //which step to restart mixing during SCF, always equal to scf_namx except for the mixing restart
    int mixing_restart_count = 0; // the number of restart mixing during SCF. Do not set mixing_restart_count as bool since I want to keep some flexibility in the future
    int mixing_restart_last = 0; // the label of mixing restart step, store the step number of the last mixing restart
 
    // to calculate the slope of drho curve during SCF, which is used to determine if SCF oscillate
    bool if_scf_oscillate(const int iteration, const double drho, const int iternum_used, const double threshold);
    
  private:
  
    // mixing_data
    Base_Mixing::Mixing* mixing = nullptr; 
    Base_Mixing::Mixing_Data rho_mdata;    
    Base_Mixing::Mixing_Data tau_mdata;    
    Base_Mixing::Mixing_Data nhat_mdata;   
    Base_Mixing::Mixing_Data dmr_mdata;    
    Base_Mixing::Plain_Mixing* mixing_highf = nullptr; 
 
    //======================================
    // private mixing parameters
    //======================================
    std::string mixing_mode = "broyden"; 
    double mixing_beta = 0.8;            
    double mixing_beta_mag = 1.6;        
    int mixing_ndim = 8;                 
    double mixing_gg0 = 0.0;             
    bool mixing_tau = false;             
    double mixing_gg0_mag = 0.0;         
    double mixing_gg0_min = 0.1;         
    double mixing_angle = 0.0;           
    bool mixing_dmr = false;             
    double* omega = nullptr;                  
    double* tpiba = nullptr;                  
    double* tpiba2 = nullptr;                 
    std::vector<double> _drho_history; 
    
    bool new_e_iteration = true;
 
    ModulePW::PW_Basis* rhopw = nullptr;  
    ModulePW::PW_Basis* rhodpw = nullptr; 
 
    void mix_rho_recip(Charge* chr);
 
    void mix_rho_real(Charge* chr);
 
    void Kerker_screen_recip(std::complex<double>* rhog);
 
    void Kerker_screen_real(double* rho);
 
    double inner_product_recip_rho(std::complex<double>* rho1, std::complex<double>* rho2);
    double inner_product_recip_simple(std::complex<double>* rho1, std::complex<double>* rho2);
    double inner_product_recip_hartree(std::complex<double>* rho1, std::complex<double>* rho2);
 
    double inner_product_real(double* rho1, double* rho2);
 
    void divide_data(std::complex<double>* data_d, std::complex<double>*& data_s, std::complex<double>*& data_hf);
    void combine_data(std::complex<double>* data_d, std::complex<double>*& data_s, std::complex<double>*& data_hf);
    void clean_data(std::complex<double>*& data_s, std::complex<double>*& data_hf);
};
 
#endif