Skip to content

Commit e84abb4

Browse files
linpeizePeizeLin
linpeize
and
PeizeLin
authored
add exx nscf file check (#6289)
* add exx nscf file check * delete some space --------- Co-authored-by: linpz <linpz@mail.ustc.edu.cn>
1 parent 808af53 commit e84abb4

16 files changed

+198
-137
lines changed

‎source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp

Lines changed: 54 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -82,13 +82,13 @@ template <typename TK, typename TR>
8282
OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
8383
HContainer<TR>*hR_in,
8484
const UnitCell& ucell_in,
85-
const K_Vectors& kv_in,
86-
std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd_in,
87-
std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc_in,
85+
const K_Vectors& kv_in,
86+
std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd_in,
87+
std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc_in,
8888
Add_Hexx_Type add_hexx_type_in,
8989
const int istep,
9090
int* two_level_step_in,
91-
const bool restart_in)
91+
const bool restart_in)
9292
: OperatorLCAO<TK, TR>(hsk_in, kv_in.kvec_d, hR_in),
9393
ucell(ucell_in),
9494
kv(kv_in),
@@ -105,42 +105,75 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
105105

106106
if (PARAM.inp.calculation == "nscf" && GlobalC::exx_info.info_global.cal_exx)
107107
{ // if nscf, read HexxR first and reallocate hR according to the read-in HexxR
108-
const std::string file_name_exx = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(GlobalV::MY_RANK);
109-
bool all_exist = true;
110-
for (int is=0;is<PARAM.inp.nspin;++is)
108+
auto file_name_list_csr = []() -> std::vector<std::string>
111109
{
112-
std::ifstream ifs(file_name_exx + "_" + std::to_string(is) + ".csr");
113-
if (!ifs) { all_exist = false; break; }
114-
}
115-
if (all_exist)
110+
std::vector<std::string> file_name_list;
111+
for (int irank=0; irank<PARAM.globalv.nproc; ++irank) {
112+
for (int is=0;is<PARAM.inp.nspin;++is) {
113+
file_name_list.push_back( PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(irank) + "_" + std::to_string(is) + ".csr" );
114+
} }
115+
return file_name_list;
116+
};
117+
auto file_name_list_cereal = []() -> std::vector<std::string>
118+
{
119+
std::vector<std::string> file_name_list;
120+
for (int irank=0; irank<PARAM.globalv.nproc; ++irank)
121+
{ file_name_list.push_back( "HexxR_" + std::to_string(irank) ); }
122+
return file_name_list;
123+
};
124+
auto check_exist = [](const std::vector<std::string> &file_name_list) -> bool
125+
{
126+
for (const std::string &file_name : file_name_list)
127+
{
128+
std::ifstream ifs(file_name);
129+
if (!ifs.is_open())
130+
{ return false; }
131+
}
132+
return true;
133+
};
134+
135+
std::cout<<" Attention: The number of MPI processes must be strictly identical between SCF and NSCF when computing exact-exchange."<<std::endl;
136+
if (check_exist(file_name_list_csr()))
116137
{
138+
const std::string file_name_exx_csr = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(PARAM.globalv.myrank);
117139
// Read HexxR in CSR format
118140
if (GlobalC::exx_info.info_ri.real_number)
119141
{
120-
ModuleIO::read_Hexxs_csr(file_name_exx, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd);
121-
if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxd, this->hR); }
142+
ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd);
143+
if (this->add_hexx_type == Add_Hexx_Type::R)
144+
{ reallocate_hcontainer(*Hexxd, this->hR); }
122145
}
123146
else
124147
{
125-
ModuleIO::read_Hexxs_csr(file_name_exx, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc);
126-
if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxc, this->hR); }
148+
ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc);
149+
if (this->add_hexx_type == Add_Hexx_Type::R)
150+
{ reallocate_hcontainer(*Hexxc, this->hR); }
127151
}
128152
}
129-
else
153+
else if (check_exist(file_name_list_cereal()))
130154
{
131155
// Read HexxR in binary format (old version)
132-
const std::string file_name_exx_cereal = PARAM.globalv.global_readin_dir + "HexxR_" + std::to_string(GlobalV::MY_RANK);
156+
const std::string file_name_exx_cereal = PARAM.globalv.global_readin_dir + "HexxR_" + std::to_string(PARAM.globalv.myrank);
157+
std::ifstream ifs(file_name_exx_cereal, std::ios::binary);
158+
if (!ifs)
159+
{ ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file < " + file_name_exx_cereal + " >."); }
133160
if (GlobalC::exx_info.info_ri.real_number)
134161
{
135162
ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxd);
136-
if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxd, this->hR); }
163+
if (this->add_hexx_type == Add_Hexx_Type::R)
164+
{ reallocate_hcontainer(*Hexxd, this->hR); }
137165
}
138166
else
139167
{
140168
ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxc);
141-
if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxc, this->hR); }
169+
if (this->add_hexx_type == Add_Hexx_Type::R)
170+
{ reallocate_hcontainer(*Hexxc, this->hR); }
142171
}
143172
}
173+
else
174+
{
175+
ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file in " + PARAM.globalv.global_readin_dir);
176+
}
144177
this->use_cell_nearest = false;
145178
}
146179
else
@@ -207,7 +240,7 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
207240
else if (this->add_hexx_type == Add_Hexx_Type::R)
208241
{
209242
// read in Hexx(R)
210-
const std::string restart_HR_path = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(GlobalV::MY_RANK);
243+
const std::string restart_HR_path = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(PARAM.globalv.myrank);
211244
bool all_exist = true;
212245
for (int is = 0; is < PARAM.inp.nspin; ++is)
213246
{
@@ -227,7 +260,7 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
227260
else
228261
{
229262
// Read HexxR in binary format (old version)
230-
const std::string restart_HR_path_cereal = GlobalC::restart.folder + "HexxR_" + std::to_string(GlobalV::MY_RANK);
263+
const std::string restart_HR_path_cereal = GlobalC::restart.folder + "HexxR_" + std::to_string(PARAM.globalv.myrank);
231264
if (GlobalC::exx_info.info_ri.real_number) {
232265
ModuleIO::read_Hexxs_cereal(restart_HR_path_cereal, *Hexxd);
233266
}

‎source/module_io/restart_exx_csr.hpp

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,12 @@ namespace ModuleIO
2828
{
2929
const std::vector<int>& R = csr.getRCoordinate(iR);
3030
TC dR({ R[0], R[1], R[2] });
31-
Hexxs[is][iat1][{iat2, dR}] = RI::Tensor<Tdata>({ static_cast<size_t>(ucell.atoms[ucell.iat2it[iat1]].nw), static_cast<size_t>(ucell.atoms[ucell.iat2it[iat2]].nw) });
31+
Hexxs[is][iat1][{iat2, dR}] = RI::Tensor<Tdata>(
32+
{
33+
static_cast<size_t>(ucell.atoms[ucell.iat2it[iat1]].nw),
34+
static_cast<size_t>(ucell.atoms[ucell.iat2it[iat2]].nw)
35+
}
36+
);
3237
}
3338
}
3439
}
@@ -44,7 +49,12 @@ namespace ModuleIO
4449
const int& npol = ucell.get_npol();
4550
const int& i = ijv.first.first * npol;
4651
const int& j = ijv.first.second * npol;
47-
Hexxs.at(is).at(ucell.iwt2iat[i]).at({ ucell.iwt2iat[j], { R[0], R[1], R[2] } })(ucell.iwt2iw[i] / npol, ucell.iwt2iw[j] / npol) = ijv.second;
52+
Hexxs.at(is).at(ucell.iwt2iat[i]).at(
53+
{
54+
ucell.iwt2iat[j],
55+
{ R[0], R[1], R[2] }
56+
}
57+
)(ucell.iwt2iw[i] / npol, ucell.iwt2iw[j] / npol) = ijv.second;
4858
}
4959
}
5060
}
@@ -57,6 +67,8 @@ namespace ModuleIO
5767
ModuleBase::TITLE("Exx_LRI", "read_Hexxs_cereal");
5868
ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal");
5969
std::ifstream ifs(file_name, std::ios::binary);
70+
if(!ifs.is_open())
71+
{ ModuleBase::WARNING_QUIT("read_Hexxs_cereal", file_name+" not found."); }
6072
cereal::BinaryInputArchive iar(ifs);
6173
iar(Hexxs);
6274
ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal");

‎source/module_ri/ABFs_Construct-PCA.cpp

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -32,14 +32,14 @@ namespace PCA
3232
dsyev_(&jobz, &uplo, &nr, a.ptr(), &nc, w, work.data(), &lwork, &info);
3333
}
3434

35-
RI::Tensor<double> get_sub_matrix(
35+
RI::Tensor<double> get_sub_matrix(
3636
const RI::Tensor<double> & m, // size: (lcaos, lcaos, abfs)
3737
const std::size_t & T,
3838
const std::size_t & L,
3939
const ModuleBase::Element_Basis_Index::Range & range,
4040
const ModuleBase::Element_Basis_Index::IndexLNM & index )
4141
{
42-
ModuleBase::TITLE("ABFs_Construct::PCA::get_sub_matrix");
42+
ModuleBase::TITLE("ABFs_Construct::PCA::get_sub_matrix");
4343
assert(m.shape.size() == 3);
4444
RI::Tensor<double> m_sub({ m.shape[0], m.shape[1], range[T][L].N });
4545
for (std::size_t ir=0; ir!=m.shape[0]; ++ir) {
@@ -74,12 +74,12 @@ namespace PCA
7474
std::vector<std::vector<std::pair<std::vector<double>, RI::Tensor<double>>>> cal_PCA(
7575
const UnitCell &ucell,
7676
const LCAO_Orbitals& orb,
77-
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &lcaos,
77+
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &lcaos,
7878
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs,
7979
const double kmesh_times )
8080
{
8181
ModuleBase::TITLE("ABFs_Construct::PCA::cal_PCA");
82-
82+
8383
const ModuleBase::Element_Basis_Index::Range
8484
range_lcaos = Exx_Abfs::Abfs_Index::construct_range( lcaos );
8585
const ModuleBase::Element_Basis_Index::IndexLNM
@@ -107,27 +107,27 @@ namespace PCA
107107
m_abfslcaos_lcaos.init_radial_table(delta_R);
108108

109109
GlobalC::exx_info.info_ri.abfs_Lmax = Lmax_bak;
110-
110+
111111
std::vector<std::vector<std::pair<std::vector<double>,RI::Tensor<double>>>> eig(abfs.size());
112112
for( std::size_t T=0; T!=abfs.size(); ++T )
113113
{
114-
const RI::Tensor<double> A = m_abfslcaos_lcaos.cal_overlap_matrix<double>(
115-
T,
116-
T,
114+
const RI::Tensor<double> A = m_abfslcaos_lcaos.cal_overlap_matrix<double>(
115+
T,
116+
T,
117+
ModuleBase::Vector3<double>{0,0,0},
117118
ModuleBase::Vector3<double>{0,0,0},
118-
ModuleBase::Vector3<double>{0,0,0},
119-
index_abfs,
119+
index_abfs,
120120
index_lcaos,
121121
index_lcaos,
122122
Matrix_Orbs21::Matrix_Order::A2BA1);
123-
123+
124124
eig[T].resize(abfs[T].size());
125125
for( std::size_t L=0; L!=abfs[T].size(); ++L )
126126
{
127127
const RI::Tensor<double> A_sub = get_sub_matrix( A, T, L, range_abfs, index_abfs );
128128
RI::Tensor<double> mm = A_sub.transpose() * A_sub;
129129
std::vector<double> eig_value(mm.shape[0]);
130-
130+
131131
int info=1;
132132

133133
tensor_dsyev('V', 'L', mm, eig_value.data(), info);
@@ -158,7 +158,7 @@ namespace PCA
158158
eig[T][L] = std::make_pair( eig_value, mm );
159159
}
160160
}
161-
161+
162162
return eig;
163163
}
164164

‎source/module_ri/ABFs_Construct-PCA.h

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -15,10 +15,10 @@ namespace ABFs_Construct
1515
{
1616
namespace PCA
1717
{
18-
extern std::vector<std::vector<std::pair<std::vector<double>,RI::Tensor<double>>>> cal_PCA(
18+
extern std::vector<std::vector<std::pair<std::vector<double>,RI::Tensor<double>>>> cal_PCA(
1919
const UnitCell& ucell,
2020
const LCAO_Orbitals &orb,
21-
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &lcaos,
21+
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &lcaos,
2222
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs, // abfs must be orthonormal
2323
const double kmesh_times );
2424
}

‎source/module_ri/Exx_LRI_interface.h

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -42,9 +42,9 @@ class Exx_LRI_Interface
4242
}
4343
Exx_LRI_Interface() = delete;
4444

45-
/// read and write Hexxs using cereal
46-
void write_Hexxs_cereal(const std::string& file_name) const;
47-
void read_Hexxs_cereal(const std::string& file_name);
45+
///// read and write Hexxs using cereal
46+
//void write_Hexxs_cereal(const std::string& file_name) const;
47+
//void read_Hexxs_cereal(const std::string& file_name);
4848

4949
std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>& get_Hexxs() const { return this->exx_ptr->Hexxs; }
5050
double &get_Eexx() const { return this->exx_ptr->Eexx; }

‎source/module_ri/Exx_LRI_interface.hpp

Lines changed: 26 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
#include <stdexcept>
1919
#include <string>
2020

21+
/*
2122
template<typename T, typename Tdata>
2223
void Exx_LRI_Interface<T, Tdata>::write_Hexxs_cereal(const std::string& file_name) const
2324
{
@@ -32,13 +33,17 @@ void Exx_LRI_Interface<T, Tdata>::write_Hexxs_cereal(const std::string& file_nam
3233
template<typename T, typename Tdata>
3334
void Exx_LRI_Interface<T, Tdata>::read_Hexxs_cereal(const std::string& file_name)
3435
{
35-
ModuleBase::TITLE("Exx_LRI", "read_Hexxs_cereal");
36-
ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal");
37-
std::ifstream ifs(file_name + "_" + std::to_string(GlobalV::MY_RANK), std::ofstream::binary);
36+
ModuleBase::TITLE("Exx_LRI_Interface", "read_Hexxs_cereal");
37+
ModuleBase::timer::tick("Exx_LRI_Interface", "read_Hexxs_cereal");
38+
const std::string file_name_rank = file_name + "_" + std::to_string(GlobalV::MY_RANK);
39+
std::ifstream ifs(file_name_rank, std::ofstream::binary);
40+
if(!ifs.is_open())
41+
{ ModuleBase::WARNING_QUIT("Exx_LRI_Interface", file_name_rank+" not found."); }
3842
cereal::BinaryInputArchive iar(ifs);
3943
iar(this->exx_ptr->Hexxs);
4044
ModuleBase::timer::tick("Exx_LRI", "read_Hexxs_cereal");
4145
}
46+
*/
4247

4348
template<typename T, typename Tdata>
4449
void Exx_LRI_Interface<T, Tdata>::init(const MPI_Comm &mpi_comm,
@@ -115,7 +120,8 @@ void Exx_LRI_Interface<T, Tdata>::exx_before_all_runners(const K_Vectors& kv, co
115120
if (this->exx_spacegroup_symmetry)
116121
{
117122
const std::array<int, 3>& period = RI_Util::get_Born_vonKarmen_period(kv);
118-
this->symrot_.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st,
123+
this->symrot_.find_irreducible_sector(
124+
ucell.symm, ucell.atoms, ucell.st,
119125
RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat);
120126
// this->symrot_.set_Cs_rotation(this->exx_ptr->get_abfs_nchis());
121127
this->symrot_.cal_Ms(kv, ucell, pv);
@@ -209,8 +215,19 @@ void Exx_LRI_Interface<T, Tdata>::exx_eachiterinit(const int istep,
209215
{ this->mix_DMk_2D.mix(dm_in.get_DMK_vector(), flag_restart); }
210216
const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>
211217
Ds = PARAM.globalv.gamma_only_local
212-
? RI_2D_Comm::split_m2D_ktoR<Tdata>(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm_in.get_paraV_pointer(), PARAM.inp.nspin)
213-
: RI_2D_Comm::split_m2D_ktoR<Tdata>(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm_in.get_paraV_pointer(), PARAM.inp.nspin, this->exx_spacegroup_symmetry);
218+
? RI_2D_Comm::split_m2D_ktoR<Tdata>(
219+
ucell,
220+
*this->exx_ptr->p_kv,
221+
this->mix_DMk_2D.get_DMk_gamma_out(),
222+
*dm_in.get_paraV_pointer(),
223+
PARAM.inp.nspin)
224+
: RI_2D_Comm::split_m2D_ktoR<Tdata>(
225+
ucell,
226+
*this->exx_ptr->p_kv,
227+
this->mix_DMk_2D.get_DMk_k_out(),
228+
*dm_in.get_paraV_pointer(),
229+
PARAM.inp.nspin,
230+
this->exx_spacegroup_symmetry);
214231

215232
if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace)
216233
{ this->cal_exx_elec(Ds, ucell,*dm_in.get_paraV_pointer(), &this->symrot_); }
@@ -240,11 +257,10 @@ void Exx_LRI_Interface<T, Tdata>::exx_hamilt2density(elecstate::ElecState& elec,
240257
{
241258
if (GlobalV::MY_RANK == 0)
242259
{
243-
try { GlobalC::restart.load_disk("Eexx", 0, 1, &this->exx_ptr->Eexx); }
260+
try
261+
{ GlobalC::restart.load_disk("Eexx", 0, 1, &this->exx_ptr->Eexx); }
244262
catch (const std::exception& e)
245-
{
246-
std::cout << "WARNING: Cannot read Eexx from disk, the energy of the 1st loop will be wrong, sbut it does not influence the subsequent loops." << std::endl;
247-
}
263+
{ std::cout << "WARNING: Cannot read Eexx from disk, the energy of the 1st loop will be wrong, sbut it does not influence the subsequent loops." << std::endl; }
248264
}
249265
Parallel_Common::bcast_double(this->exx_ptr->Eexx);
250266
this->exx_ptr->Eexx /= GlobalC::exx_info.info_global.hybrid_alpha;

‎source/module_ri/exx_abfs-abfs_index.cpp

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,9 @@ ModuleBase::Element_Basis_Index::Range
2929
range[T].resize( orb[T].size() );
3030
for( size_t L=0; L!=range[T].size(); ++L )
3131
{
32-
range[T][L].N = orb[T][L].size();
33-
range[T][L].M = 2*L+1;
34-
}
32+
range[T][L].N = orb[T][L].size();
33+
range[T][L].M = 2*L+1;
34+
}
3535
}
3636
return range;
3737
}

‎source/module_ri/exx_abfs-abfs_index.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ class LCAO_Orbitals;
1212
class Exx_Abfs::Abfs_Index
1313
{
1414
public:
15-
static ModuleBase::Element_Basis_Index::Range construct_range( const LCAO_Orbitals &orb );
15+
static ModuleBase::Element_Basis_Index::Range construct_range( const LCAO_Orbitals &orb );
1616
static ModuleBase::Element_Basis_Index::Range construct_range( const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orb );
1717
};
1818

0 commit comments

Comments
 (0)