Skip to content
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions source/module_esolver/esolver_ks_lcao.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -183,12 +183,12 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
// initialize 2-center radial tables for EXX-LRI
if (GlobalC::exx_info.info_ri.real_number)
{
this->exx_lri_double->init(MPI_COMM_WORLD, ucell, this->kv, orb_);
this->exd->init(MPI_COMM_WORLD, ucell, this->kv, orb_);
this->exd->exx_before_all_runners(this->kv, ucell, this->pv);
}
else
{
this->exx_lri_complex->init(MPI_COMM_WORLD, ucell, this->kv, orb_);
this->exc->init(MPI_COMM_WORLD, ucell, this->kv, orb_);
this->exc->exx_before_all_runners(this->kv, ucell, this->pv);
}
}
Expand Down
56 changes: 28 additions & 28 deletions source/module_ri/Exx_LRI.h
Original file line number Diff line number Diff line change
Expand Up @@ -21,21 +21,21 @@
#include "module_exx_symmetry/symmetry_rotation.h"

class Parallel_Orbitals;

template<typename T, typename Tdata>
class RPA_LRI;

template<typename T, typename Tdata>
class Exx_LRI_Interface;

namespace LR
{
template<typename T, typename TR>
class ESolver_LR;
namespace LR
{
template<typename T, typename TR>
class ESolver_LR;

template<typename T>
class OperatorLREXX;
}
template<typename T>
class OperatorLREXX;
}

template<typename Tdata>
class Exx_LRI
Expand All @@ -49,37 +49,37 @@ class Exx_LRI
using TatomR = std::array<double,Ndim>; // tmp

public:
Exx_LRI(const Exx_Info::Exx_Info_RI& info_in) :info(info_in) {}
Exx_LRI operator=(const Exx_LRI&) = delete;
Exx_LRI operator=(Exx_LRI&&);
Exx_LRI(const Exx_Info::Exx_Info_RI& info_in) :info(info_in) {}
Exx_LRI operator=(const Exx_LRI&) = delete;
Exx_LRI operator=(Exx_LRI&&);

void reset_Cs(const std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs_in) { this->exx_lri.set_Cs(Cs_in, this->info.C_threshold); }
void reset_Vs(const std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Vs_in) { this->exx_lri.set_Vs(Vs_in, this->info.V_threshold); }
void reset_Cs(const std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs_in) { this->exx_lri.set_Cs(Cs_in, this->info.C_threshold); }
void reset_Vs(const std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Vs_in) { this->exx_lri.set_Vs(Vs_in, this->info.V_threshold); }

void init(const MPI_Comm &mpi_comm_in,
void init(const MPI_Comm &mpi_comm_in,
const UnitCell &ucell,
const K_Vectors &kv_in,
const K_Vectors &kv_in,
const LCAO_Orbitals& orb);
void cal_exx_force(const int& nat);
void cal_exx_stress(const double& omega, const double& lat0);
void cal_exx_stress(const double& omega, const double& lat0);
void cal_exx_ions(const UnitCell& ucell, const bool write_cv = false);
void cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>& Ds,
void cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>& Ds,
const UnitCell& ucell,
const Parallel_Orbitals& pv,
const ModuleSymmetry::Symmetry_rotation* p_symrot = nullptr);
std::vector<std::vector<int>> get_abfs_nchis() const;
const Parallel_Orbitals& pv,
const ModuleSymmetry::Symmetry_rotation* p_symrot = nullptr);
std::vector<std::vector<int>> get_abfs_nchis() const;

std::vector< std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> Hexxs;
double Eexx;
double Eexx;
ModuleBase::matrix force_exx;
ModuleBase::matrix stress_exx;


private:
const Exx_Info::Exx_Info_RI &info;
MPI_Comm mpi_comm;
const K_Vectors *p_kv = nullptr;
std::vector<double> orb_cutoff_;
std::vector<double> orb_cutoff_;

std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> lcaos;
std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs;
Expand All @@ -89,16 +89,16 @@ class Exx_LRI
RI::Exx<TA,Tcell,Ndim,Tdata> exx_lri;

void post_process_Hexx( std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> &Hexxs_io ) const;
double post_process_Eexx(const double& Eexx_in) const;
double post_process_Eexx(const double& Eexx_in) const;

friend class RPA_LRI<double, Tdata>;
friend class RPA_LRI<std::complex<double>, Tdata>;
friend class Exx_LRI_Interface<double, Tdata>;
friend class Exx_LRI_Interface<std::complex<double>, Tdata>;
friend class LR::ESolver_LR<double, double>;
friend class LR::ESolver_LR<std::complex<double>, double>;
friend class LR::OperatorLREXX<double>;
friend class LR::OperatorLREXX<std::complex<double>>;
friend class LR::ESolver_LR<double, double>;
friend class LR::ESolver_LR<std::complex<double>, double>;
friend class LR::OperatorLREXX<double>;
friend class LR::OperatorLREXX<std::complex<double>>;
};

#include "Exx_LRI.hpp"
Expand Down
13 changes: 4 additions & 9 deletions source/module_ri/Exx_LRI.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,9 @@
#include <string>

template<typename Tdata>
void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
const UnitCell &ucell,
const K_Vectors &kv_in,
const K_Vectors &kv_in,
const LCAO_Orbitals& orb)
{
ModuleBase::TITLE("Exx_LRI","init");
Expand Down Expand Up @@ -130,7 +130,7 @@ void Exx_LRI<Tdata>::cal_exx_ions(const UnitCell& ucell,
this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period);

// std::max(3) for gamma_only, list_A2 should contain cell {-1,0,1}. In the future distribute will be neighbour.
const std::array<Tcell,Ndim> period_Vs = LRI_CV_Tools::cal_latvec_range<Tcell>(1+this->info.ccp_rmesh_times, ucell, orb_cutoff_);
const std::array<Tcell,Ndim> period_Vs = LRI_CV_Tools::cal_latvec_range<Tcell>(1+this->info.ccp_rmesh_times, ucell, orb_cutoff_);
const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false);

Expand Down Expand Up @@ -237,7 +237,7 @@ void Exx_LRI<Tdata>::cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, R
}
this->Eexx = post_process_Eexx(this->Eexx);
this->exx_lri.set_symmetry(false, {});
ModuleBase::timer::tick("Exx_LRI", "cal_exx_elec");
ModuleBase::timer::tick("Exx_LRI", "cal_exx_elec");
}

template<typename Tdata>
Expand Down Expand Up @@ -283,11 +283,6 @@ void Exx_LRI<Tdata>::cal_exx_force(const int& nat)
ModuleBase::TITLE("Exx_LRI","cal_exx_force");
ModuleBase::timer::tick("Exx_LRI", "cal_exx_force");

if (!this->exx_lri.flag_finish.D)
{
ModuleBase::WARNING_QUIT("Force_Stress_LCAO", "Cannot calculate EXX force when the first PBE loop is not converged.");
}

this->force_exx.create(nat, Ndim);
for(int is=0; is<PARAM.inp.nspin; ++is)
{
Expand Down
44 changes: 34 additions & 10 deletions source/module_ri/Exx_LRI_interface.h
Original file line number Diff line number Diff line change
Expand Up @@ -44,12 +44,19 @@ class Exx_LRI_Interface
void read_Hexxs_cereal(const std::string& file_name);

std::vector<std::map<int, std::map<TAC, RI::Tensor<Tdata>>>>& get_Hexxs() const { return this->exx_ptr->Hexxs; }

double& get_Eexx() const { return this->exx_ptr->Eexx; }

// Processes in ESolver_KS_LCAO
/// @brief Exx_LRI::init()
void init(const MPI_Comm &mpi_comm,
const UnitCell &ucell,
const K_Vectors &kv,
const LCAO_Orbitals& orb);

// Processes in ESolver_KS_LCAO
/// @brief in before_all_runners: set symmetry according to irreducible k-points
/// since k-points are not reduced again after the variation of the cell and exx-symmetry must be consistent with k-points.
/// since k-points are not reduced again after the variation of the cell and exx-symmetry must be consistent with k-points.
/// In the future, we will reduce k-points again during cell-relax, then this setting can be moved to `exx_beforescf`.
void exx_before_all_runners(const K_Vectors& kv, const UnitCell& ucell, const Parallel_2D& pv);

Expand All @@ -58,23 +65,23 @@ class Exx_LRI_Interface

/// @brief in eachiterinit: do DM mixing and calculate Hexx when entering 2nd SCF
void exx_eachiterinit(const int istep,
const UnitCell& ucell,
const UnitCell& ucell,
const elecstate::DensityMatrix<T, double>& dm/**< double should be Tdata if complex-PBE-DM is supported*/,
const K_Vectors& kv,
const K_Vectors& kv,
const int& iter);

/// @brief in hamilt2density: calculate Hexx and Eexx
void exx_hamilt2density(elecstate::ElecState& elec, const Parallel_Orbitals& pv, const int iter);

/// @brief in iter_finish: write Hexx, do something according to whether SCF is converged
void exx_iter_finish(const K_Vectors& kv,
void exx_iter_finish(const K_Vectors& kv,
const UnitCell& ucell,
hamilt::Hamilt<T>& hamilt,
elecstate::ElecState& elec,
hamilt::Hamilt<T>& hamilt,
elecstate::ElecState& elec,
Charge_Mixing& chgmix,
const double& scf_ene_thr,
int& iter,
const int istep,
const double& scf_ene_thr,
int& iter,
const int istep,
bool& conv_esolver);
/// @brief: in do_after_converge: add exx operators; do DM mixing if seperate loop
bool exx_after_converge(const UnitCell& ucell,
Expand All @@ -86,6 +93,13 @@ class Exx_LRI_Interface
const int& istep,
const double& etot,
const double& scf_ene_thr);

/// @brief: in cal_exx_force: Exx_LRI::cal_exx_force()
void cal_exx_force(const double& omega, const double& lat0);

/// @brief: in cal_exx_stress: Exx_LRI::cal_exx_stress()
void cal_exx_stress(const double& omega, const double& lat0);

int two_level_step = 0;
double etot_last_outer_loop = 0.0;
elecstate::DensityMatrix<T, double>* dm_last_step;
Expand All @@ -95,6 +109,16 @@ class Exx_LRI_Interface

bool exx_spacegroup_symmetry = false;
ModuleSymmetry::Symmetry_rotation symrot_;

struct Flag_Finish
{
bool init = false;
bool ions = false;
bool elec = false;
bool force = false;
bool stress = false;
};
Flag_Finish flag_finish;
};

#include "Exx_LRI_interface.hpp"
Expand Down
Loading
Loading