From fc5f81ecf42068cdae3f1fb341d04e9f9a7c8039 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Wed, 3 Jan 2024 04:03:23 +0800 Subject: [PATCH 1/3] unify matrix writing function --- source/Makefile.Objects | 1 - source/module_esolver/esolver_ks_lcao.cpp | 13 +- .../module_esolver/esolver_ks_lcao_tddft.cpp | 13 +- .../hamilt_lcaodft/FORCE_gamma.cpp | 10 +- .../hamilt_lcaodft/FORCE_k.cpp | 10 +- source/module_io/CMakeLists.txt | 1 - source/module_io/write_HS.cpp | 1363 ----------------- source/module_io/write_HS.h | 81 +- source/module_io/write_HS.hpp | 486 ++++++ 9 files changed, 525 insertions(+), 1453 deletions(-) delete mode 100644 source/module_io/write_HS.cpp create mode 100644 source/module_io/write_HS.hpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 6178ef5e32..b9768dab5f 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -456,7 +456,6 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\ unk_overlap_lcao.o\ read_wfc_nao.o\ write_wfc_nao.o\ - write_HS.o\ write_HS_sparse.o\ single_R_io.o\ write_HS_R.o\ diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index bcb23c4765..4c48eadcc6 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -686,14 +686,11 @@ namespace ModuleESolver { hamilt::MatrixBlock h_mat, s_mat; this->p_hamilt->matrix(h_mat, s_mat); - ModuleIO::saving_HS(istep, - h_mat.p, - s_mat.p, - bit, - hsolver::HSolverLCAO::out_mat_hs, - "data-" + std::to_string(ik), - this->LOWF.ParaV[0], - 1); // LiuXh, 2017-03-21 + if (hsolver::HSolverLCAO::out_mat_hs) + { + ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, bit, 1, GlobalV::out_app_flag, "H", "data-" + std::to_string(ik), *this->LOWF.ParaV, GlobalV::DRANK); + ModuleIO::save_mat(istep, s_mat.p, GlobalV::NLOCAL, bit, 1, GlobalV::out_app_flag, "S", "data-" + std::to_string(ik), *this->LOWF.ParaV, GlobalV::DRANK); + } } } } diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index a26e97f75e..4f99c016e6 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -258,14 +258,11 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) { hamilt::MatrixBlock> h_mat, s_mat; this->p_hamilt->matrix(h_mat, s_mat); - ModuleIO::saving_HS(istep, - h_mat.p, - s_mat.p, - bit, - hsolver::HSolverLCAO>::out_mat_hs, - "data-" + std::to_string(ik), - this->LOWF.ParaV[0], - 1); // LiuXh, 2017-03-21 + if (hsolver::HSolverLCAO>::out_mat_hs) + { + ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, bit, 1, GlobalV::out_app_flag, "H", "data-" + std::to_string(ik), *this->LOWF.ParaV, GlobalV::DRANK); + ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, bit, 1, GlobalV::out_app_flag, "S", "data-" + std::to_string(ik), *this->LOWF.ParaV, GlobalV::DRANK); + } } } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp index 7c84846bc8..f1ba3c8b69 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp @@ -210,14 +210,8 @@ void Force_LCAO_gamma::allocate_gamma(const Parallel_Orbitals& pv) this->UHM->genH .build_ST_new('S', cal_deri, GlobalC::ucell, this->UHM->genH.LM->Sloc.data(), INPUT.cal_syns, INPUT.dmax); bool bit = false; // LiuXh, 2017-03-21 - ModuleIO::saving_HS(0, - this->UHM->genH.LM->Hloc.data(), - this->UHM->genH.LM->Sloc.data(), - bit, - 1, - "data-" + std::to_string(0), - this->ParaV[0], - 0); // LiuXh, 2017-03-21 + ModuleIO::save_mat(0, this->UHM->genH.LM->Hloc.data(), GlobalV::NLOCAL, bit, 0, GlobalV::out_app_flag, "H", "data-" + std::to_string(0), *this->ParaV, GlobalV::DRANK); + ModuleIO::save_mat(0, this->UHM->genH.LM->Sloc.data(), GlobalV::NLOCAL, bit, 0, GlobalV::out_app_flag, "S", "data-" + std::to_string(0), *this->ParaV, GlobalV::DRANK); } ModuleBase::timer::tick("Force_LCAO_gamma", "allocate_gamma"); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index db82831068..3767105b22 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -243,14 +243,8 @@ void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv, this->UHM->genH.LM->zeros_HSk('S'); this->UHM->genH.LM->folding_fixedH(ik, kvec_d, 1); bool bit = false; // LiuXh, 2017-03-21 - ModuleIO::saving_HS(0, - this->UHM->genH.LM->Hloc2.data(), - this->UHM->genH.LM->Sloc2.data(), - bit, - 1, - "data-" + std::to_string(ik), - this->ParaV[0], - 0); // LiuXh, 2017-03-21 + ModuleIO::save_mat(0, this->UHM->genH.LM->Hloc2.data(), GlobalV::NLOCAL, bit, 0, GlobalV::out_app_flag, "H", "data-" + std::to_string(ik), *this->ParaV, GlobalV::DRANK); + ModuleIO::save_mat(0, this->UHM->genH.LM->Sloc2.data(), GlobalV::NLOCAL, bit, 0, GlobalV::out_app_flag, "S", "data-" + std::to_string(ik), *this->ParaV, GlobalV::DRANK); } } diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index f36fa60c54..cc9682d472 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -53,7 +53,6 @@ if(ENABLE_LCAO) read_dm.cpp read_wfc_nao.cpp write_wfc_nao.cpp - write_HS.cpp write_dm.cpp dos_nao.cpp output_dm.cpp diff --git a/source/module_io/write_HS.cpp b/source/module_io/write_HS.cpp deleted file mode 100644 index 0e677f8f6d..0000000000 --- a/source/module_io/write_HS.cpp +++ /dev/null @@ -1,1363 +0,0 @@ -#include "write_HS.h" - -#include "module_base/parallel_reduce.h" -#include "module_base/timer.h" -#include "module_cell/module_neighbor/sltk_grid_driver.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" - -void ModuleIO::saving_HS(const int istep, - const double* Hloc, - const double* Sloc, - const bool bit, - const int& out_mat_hs, - const std::string& file_name, - const Parallel_Orbitals& pv, - bool tri) -{ - - if (out_mat_hs == 1) - { - if (tri) - { - save_HS_triangle(istep, Hloc, Sloc, bit, file_name, pv); - } - else - { - save_HS_complete(istep, Hloc, Sloc, bit, file_name, pv); - } - } - else if (out_mat_hs == 2) - { - if (tri) - { - save_HS_triangle(istep, Hloc, Sloc, bit, file_name, pv); - } - else - { - save_HS_complete(istep, Hloc, Sloc, bit, file_name, pv); - } - } - else if (out_mat_hs == 3) - { - // please call individually - } - else if (out_mat_hs == 0) - { - // do nothing. - } - else - { - ModuleBase::WARNING("Diago_LCAO_Matrix", "unrecorganized out_mat_hs value."); - } - return; -} - -/* -void ModuleIO::save_HS_ccf(const int &iter, const int &Hnnz, const int *colptr_H, const int *rowind_H, - const double *nzval_H, const double *nzval_S, bool bit) -{ - ModuleBase::TITLE("ModuleIO","save_HS_ccf"); - - if(GlobalV::DRANK!=0)return; - - std::stringstream ssh; - std::stringstream sss; - - if(bit) - { - ssh << GlobalV::global_out_dir << "H_bit.ccf"; - sss << GlobalV::global_out_dir << "S_bit.ccf"; - } - else - { - // mohan update 2021-02-10 - ssh << GlobalV::global_out_dir << "H" << ELEC_scf::iter << "_" << iter+1 << ".ccf"; - sss << GlobalV::global_out_dir << "S" << ELEC_scf::iter << "_" << iter+1 << ".ccf"; - } - - if(bit) - { - FILE *g1 = fopen(ssh.str().c_str(),"wb"); - FILE *g2 = fopen(sss.str().c_str(),"wb"); - - fwrite(&GlobalV::NLOCAL,sizeof(int),1,g1); - fwrite(&Hnnz,sizeof(int),1,g1); - fwrite(&GlobalV::NLOCAL,sizeof(int),1,g2); - fwrite(&Hnnz,sizeof(int),1,g2); - - fclose(g1); - fclose(g2); - } - - - if(!bit) - { - std::ofstream g1(ssh.str().c_str()); - std::ofstream g2(sss.str().c_str()); - - g1 << GlobalV::NLOCAL << " " << Hnnz << std::endl; - g2 << GlobalV::NLOCAL << " " << Hnnz << std::endl; - - for(int i=0; i= 0) - { - // data collection - for (int j = i; j < GlobalV::NLOCAL; j++) - { - ic = pv.global2local_col(j); - if (ic >= 0) - { - int iic; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { - iic = ir + ic * pv.nrow; - } - else - { - iic = ir * pv.ncol + ic; - } - // lineH[j-i] = H[ir*pv.ncol+ic]; - // lineS[j-i] = S[ir*pv.ncol+ic]; - lineH[j - i] = H[iic]; - lineS[j - i] = S[iic]; - } - } - } - else - { - // do nothing - } - - Parallel_Reduce::reduce_all(lineH, GlobalV::NLOCAL - i); - Parallel_Reduce::reduce_all(lineS, GlobalV::NLOCAL - i); - - if (GlobalV::DRANK == 0) - { - for (int j = i; j < GlobalV::NLOCAL; j++) - { - fwrite(&lineH[j - i], sizeof(double), 1, g1); - fwrite(&lineS[j - i], sizeof(double), 1, g2); - } - } - delete[] lineH; - delete[] lineS; - - MPI_Barrier(DIAG_WORLD); - } - - if (GlobalV::DRANK == 0) - { - fclose(g1); - fclose(g2); - } -#else - FILE *g1 = fopen(ssh.str().c_str(), "wb"); - FILE *g2 = fopen(sss.str().c_str(), "wb"); - - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g1); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g2); - - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = i; j < GlobalV::NLOCAL; j++) - { - fwrite(&H[i * GlobalV::NLOCAL + j], sizeof(double), 1, g1); - fwrite(&S[i * GlobalV::NLOCAL + j], sizeof(double), 1, g2); - } - } - fclose(g1); - fclose(g2); -#endif - } // end bit - else - { -#ifdef __MPI - std::ofstream g1; - std::ofstream g2; - - if (GlobalV::DRANK == 0) - { - if (GlobalV::out_app_flag) - { - g1.open(ssh.str().c_str(), std::ofstream::app); - g2.open(sss.str().c_str(), std::ofstream::app); - } - else - { - g1.open(ssh.str().c_str()); - g2.open(sss.str().c_str()); - } - g1 << GlobalV::NLOCAL; - g2 << GlobalV::NLOCAL; - } - - int ir, ic; - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - double *lineH = new double[GlobalV::NLOCAL - i]; - double *lineS = new double[GlobalV::NLOCAL - i]; - ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL - i); - ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL - i); - - ir = pv.global2local_row(i); - if (ir >= 0) - { - // data collection - for (int j = i; j < GlobalV::NLOCAL; j++) - { - ic = pv.global2local_col(j); - if (ic >= 0) - { - int iic; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { - iic = ir + ic * pv.nrow; - } - else - { - iic = ir * pv.ncol + ic; - } - // lineH[j-i] = H[ir*pv.ncol+ic]; - // lineS[j-i] = S[ir*pv.ncol+ic]; - lineH[j - i] = H[iic]; - lineS[j - i] = S[iic]; - } - } - } - else - { - // do nothing - } - - Parallel_Reduce::reduce_all(lineH, GlobalV::NLOCAL - i); - Parallel_Reduce::reduce_all(lineS, GlobalV::NLOCAL - i); - - if (GlobalV::DRANK == 0) - { - for (int j = i; j < GlobalV::NLOCAL; j++) - { - g1 << " " << lineH[j - i]; - g2 << " " << lineS[j - i]; - } - g1 << std::endl; - g2 << std::endl; - } - delete[] lineH; - delete[] lineS; - } - - // if (GlobalV::DRANK==0); - if (GlobalV::DRANK == 0) // Peize Lin delete ; at 2020.01.31 - { - g1.close(); - g2.close(); - } -#else - if (GlobalV::out_app_flag) - { - std::ofstream g1(ssh.str().c_str(), ofstream::app); - std::ofstream g2(sss.str().c_str(), ofstream::app); - } - else - { - std::ofstream g1(ssh.str().c_str()); - std::ofstream g2(sss.str().c_str()); - } - - g1 << GlobalV::NLOCAL; - g2 << GlobalV::NLOCAL; - - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = i; j < GlobalV::NLOCAL; j++) - { - g1 << " " << H[i * GlobalV::NLOCAL + j]; - g2 << " " << S[i * GlobalV::NLOCAL + j]; - } - g1 << std::endl; - g2 << std::endl; - } - g1.close(); - g2.close(); -#endif - } - - ModuleBase::timer::tick("ModuleIO", "save_HS_bit"); - return; -} - -// mohan add 2010/3/20, output H and S matrix, convinence for diagonalization -// test or save the middle information for next start. -void ModuleIO::save_HS_complete(const int istep, - const double* H, - const double* S, - const bool bit, - const std::string& file_name, - const Parallel_Orbitals& pv) -{ - ModuleBase::TITLE("ModuleIO", "save_HS_bit"); - ModuleBase::timer::tick("ModuleIO", "save_HS_bit"); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Dimension of H and S", GlobalV::NLOCAL); - - std::stringstream ssh; - std::stringstream sss; - - if (bit) - { - ssh << GlobalV::global_out_dir << file_name + "-H-bit"; - sss << GlobalV::global_out_dir << file_name + "-S-bit"; - } - else - { - if (GlobalV::out_app_flag) - { - ssh << GlobalV::global_out_dir << file_name + "-H"; - sss << GlobalV::global_out_dir << file_name + "-S"; - } - else - { - ssh << GlobalV::global_out_dir << istep << "_" << file_name + "-H"; - sss << GlobalV::global_out_dir << istep << "_" << file_name + "-S"; - } - } - if (bit) - { -#ifdef __MPI - FILE *g1 = nullptr; - FILE *g2 = nullptr; - - if (GlobalV::DRANK == 0) - { - g1 = fopen(ssh.str().c_str(), "wb"); - g2 = fopen(sss.str().c_str(), "wb"); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g1); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g2); - } - - int ir, ic; - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - double *lineH = new double[GlobalV::NLOCAL]; - double *lineS = new double[GlobalV::NLOCAL]; - ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL); - ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL); - - ir = pv.global2local_row(i); - if (ir >= 0) - { - // data collection - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - ic = pv.global2local_col(j); - if (ic >= 0) - { - int iic; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { - iic = ir + ic * pv.nrow; - } - else - { - iic = ir * pv.ncol + ic; - } - // lineH[j-i] = H[ir*pv.ncol+ic]; - // lineS[j-i] = S[ir*pv.ncol+ic]; - lineH[j] = H[iic]; - lineS[j] = S[iic]; - } - } - } - else - { - // do nothing - } - - Parallel_Reduce::reduce_all(lineH, GlobalV::NLOCAL); - Parallel_Reduce::reduce_all(lineS, GlobalV::NLOCAL); - - if (GlobalV::DRANK == 0) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - fwrite(&lineH[j], sizeof(double), 1, g1); - fwrite(&lineS[j], sizeof(double), 1, g2); - } - } - delete[] lineH; - delete[] lineS; - - MPI_Barrier(DIAG_WORLD); - } - - if (GlobalV::DRANK == 0) - { - fclose(g1); - fclose(g2); - } -#else - FILE *g1 = fopen(ssh.str().c_str(), "wb"); - FILE *g2 = fopen(sss.str().c_str(), "wb"); - - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g1); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g2); - - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - fwrite(&H[i * GlobalV::NLOCAL + j], sizeof(double), 1, g1); - fwrite(&S[i * GlobalV::NLOCAL + j], sizeof(double), 1, g2); - } - } - fclose(g1); - fclose(g2); -#endif - } // end bit - else - { -#ifdef __MPI - std::ofstream g1; - std::ofstream g2; - - if (GlobalV::DRANK == 0) - { - if (GlobalV::out_app_flag) - { - g1.open(ssh.str().c_str(), std::ofstream::app); - g2.open(sss.str().c_str(), std::ofstream::app); - } - else - { - g1.open(ssh.str().c_str()); - g2.open(sss.str().c_str()); - } - g1 << GlobalV::NLOCAL; - g2 << GlobalV::NLOCAL; - } - - int ir, ic; - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - double *lineH = new double[GlobalV::NLOCAL]; - double *lineS = new double[GlobalV::NLOCAL]; - ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL); - ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL); - - ir = pv.global2local_row(i); - if (ir >= 0) - { - // data collection - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - ic = pv.global2local_col(j); - if (ic >= 0) - { - int iic; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { - iic = ir + ic * pv.nrow; - } - else - { - iic = ir * pv.ncol + ic; - } - // lineH[j-i] = H[ir*pv.ncol+ic]; - // lineS[j-i] = S[ir*pv.ncol+ic]; - lineH[j] = H[iic]; - lineS[j] = S[iic]; - } - } - } - else - { - // do nothing - } - - Parallel_Reduce::reduce_all(lineH, GlobalV::NLOCAL); - Parallel_Reduce::reduce_all(lineS, GlobalV::NLOCAL); - - if (GlobalV::DRANK == 0) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - g1 << " " << lineH[j]; - g2 << " " << lineS[j]; - } - g1 << std::endl; - g2 << std::endl; - } - delete[] lineH; - delete[] lineS; - } - - // if (GlobalV::DRANK==0); - if (GlobalV::DRANK == 0) // Peize Lin delete ; at 2020.01.31 - { - g1.close(); - g2.close(); - } -#else - if (GlobalV::out_app_flag) - { - std::ofstream g1(ssh.str().c_str(), ofstream::app); - std::ofstream g2(sss.str().c_str(), ofstream::app); - } - else - { - std::ofstream g1(ssh.str().c_str()); - std::ofstream g2(sss.str().c_str()); - } - - g1 << GlobalV::NLOCAL; - g2 << GlobalV::NLOCAL; - - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - g1 << " " << H[i * GlobalV::NLOCAL + j]; - g2 << " " << S[i * GlobalV::NLOCAL + j]; - } - g1 << std::endl; - g2 << std::endl; - } - g1.close(); - g2.close(); -#endif - } - - ModuleBase::timer::tick("ModuleIO", "save_HS_bit"); - return; -} - -// LiuXh, 2017-03-21 -void ModuleIO::saving_HS(const int istep, - std::complex* Hloc, - std::complex* Sloc, - const bool bit, - const int& out_mat_hs, - const std::string& file_name, - const Parallel_Orbitals& pv, - bool tri) -{ - if (out_mat_hs == 1) - { - if (tri) - { - save_HS_complex_triangle(istep, Hloc, Sloc, bit, file_name, pv); - } - else - { - save_HS_complex_complete(istep, Hloc, Sloc, bit, file_name, pv); - } - } - else if (out_mat_hs == 0) - { - // do nothing. - } - else - { - ModuleBase::WARNING("Diago_LCAO_Matrix", "unrecorganized out_mat_hs value."); - } - return; -} - -// LiuXh, 2017-03-21 -void ModuleIO::save_HS_complex_triangle(const int istep, - std::complex* H, - std::complex* S, - const bool bit, - const std::string& file_name, - const Parallel_Orbitals& pv) -{ - ModuleBase::TITLE("ModuleIO", "save_HS_bit"); - ModuleBase::timer::tick("ModuleIO", "save_HS_bit"); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Dimension of H and S", GlobalV::NLOCAL); - - std::stringstream ssh; - std::stringstream sss; - - if (bit) - { - ssh << GlobalV::global_out_dir << file_name + "-H-bit"; - sss << GlobalV::global_out_dir << file_name + "-S-bit"; - } - else - { - if (GlobalV::out_app_flag) - { - ssh << GlobalV::global_out_dir << file_name + "-H"; - sss << GlobalV::global_out_dir << file_name + "-S"; - } - else - { - ssh << GlobalV::global_out_dir << istep << "_" << file_name + "-H"; - sss << GlobalV::global_out_dir << istep << "_" << file_name + "-S"; - } - } - - if (bit) - { -#ifdef __MPI - FILE *g1 = nullptr; - FILE *g2 = nullptr; - - if (GlobalV::DRANK == 0) - { - g1 = fopen(ssh.str().c_str(), "wb"); - g2 = fopen(sss.str().c_str(), "wb"); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g1); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g2); - } - - int ir, ic; - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - std::complex *lineH = new std::complex[GlobalV::NLOCAL - i]; - std::complex *lineS = new std::complex[GlobalV::NLOCAL - i]; - ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL - i); - ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL - i); - - ir = pv.global2local_row(i); - if (ir >= 0) - { - // data collection - for (int j = i; j < GlobalV::NLOCAL; j++) - { - ic = pv.global2local_col(j); - if (ic >= 0) - { - int iic; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { - iic = ir + ic * pv.nrow; - } - else - { - iic = ir * pv.ncol + ic; - } - // lineH[j-i] = H[ir*pv.ncol+ic]; - // lineS[j-i] = S[ir*pv.ncol+ic]; - lineH[j - i] = H[iic]; - lineS[j - i] = S[iic]; - } - } - } - else - { - // do nothing - } - - Parallel_Reduce::reduce_pool(lineH, GlobalV::NLOCAL - i); - Parallel_Reduce::reduce_pool(lineS, GlobalV::NLOCAL - i); - - if (GlobalV::DRANK == 0) - { - for (int j = i; j < GlobalV::NLOCAL; j++) - { - fwrite(&lineH[j - i], sizeof(std::complex), 1, g1); - fwrite(&lineS[j - i], sizeof(std::complex), 1, g2); - } - } - delete[] lineH; - delete[] lineS; - - MPI_Barrier(DIAG_WORLD); - } - - if (GlobalV::DRANK == 0) - { - fclose(g1); - fclose(g2); - } -#else - FILE *g1 = fopen(ssh.str().c_str(), "wb"); - FILE *g2 = fopen(sss.str().c_str(), "wb"); - - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g1); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g2); - - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = i; j < GlobalV::NLOCAL; j++) - { - fwrite(&H[i * GlobalV::NLOCAL + j], sizeof(std::complex), 1, g1); - fwrite(&S[i * GlobalV::NLOCAL + j], sizeof(std::complex), 1, g2); - } - } - fclose(g1); - fclose(g2); -#endif - } // end bit - else - { -#ifdef __MPI - std::ofstream g1; - std::ofstream g2; - - if (GlobalV::DRANK == 0) - { - if (GlobalV::out_app_flag) - { - g1.open(ssh.str().c_str(), std::ofstream::app); - g2.open(sss.str().c_str(), std::ofstream::app); - } - else - { - g1.open(ssh.str().c_str()); - g2.open(sss.str().c_str()); - } - g1 << GlobalV::NLOCAL; - g2 << GlobalV::NLOCAL; - } - - int ir, ic; - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - std::complex *lineH = new std::complex[GlobalV::NLOCAL - i]; - std::complex *lineS = new std::complex[GlobalV::NLOCAL - i]; - ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL - i); - ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL - i); - - ir = pv.global2local_row(i); - if (ir >= 0) - { - // data collection - for (int j = i; j < GlobalV::NLOCAL; j++) - { - ic = pv.global2local_col(j); - if (ic >= 0) - { - int iic; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { - iic = ir + ic * pv.nrow; - } - else - { - iic = ir * pv.ncol + ic; - } - // lineH[j-i] = H[ir*pv.ncol+ic]; - // lineS[j-i] = S[ir*pv.ncol+ic]; - lineH[j - i] = H[iic]; - lineS[j - i] = S[iic]; - } - } - } - else - { - // do nothing - } - - Parallel_Reduce::reduce_pool(lineH, GlobalV::NLOCAL - i); - Parallel_Reduce::reduce_pool(lineS, GlobalV::NLOCAL - i); - - if (GlobalV::DRANK == 0) - { - for (int j = i; j < GlobalV::NLOCAL; j++) - { - g1 << " " << lineH[j - i]; - g2 << " " << lineS[j - i]; - } - g1 << std::endl; - g2 << std::endl; - } - delete[] lineH; - delete[] lineS; - } - - // if (GlobalV::DRANK==0); - if (GlobalV::DRANK == 0) // Peize Lin delete ; at 2020.01.31 - { - g1.close(); - g2.close(); - } -#else - - if (GlobalV::out_app_flag) - { - std::ofstream g1(ssh.str().c_str(), ofstream::app); - std::ofstream g2(sss.str().c_str(), ofstream::app); - } - else - { - std::ofstream g1(ssh.str().c_str()); - std::ofstream g2(sss.str().c_str()); - } - - g1 << GlobalV::NLOCAL; - g2 << GlobalV::NLOCAL; - - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = i; j < GlobalV::NLOCAL; j++) - { - g1 << " " << H[i * GlobalV::NLOCAL + j]; - g2 << " " << S[i * GlobalV::NLOCAL + j]; - } - g1 << std::endl; - g2 << std::endl; - } - g1.close(); - g2.close(); -#endif - } - - ModuleBase::timer::tick("ModuleIO", "save_HS_bit"); - return; -} - -// LiuXh, 2017-03-21 -void ModuleIO::save_HS_complex_complete(const int istep, - std::complex* H, - std::complex* S, - const bool bit, - const std::string& file_name, - const Parallel_Orbitals& pv) -{ - ModuleBase::TITLE("ModuleIO", "save_HS_bit"); - ModuleBase::timer::tick("ModuleIO", "save_HS_bit"); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Dimension of H and S", GlobalV::NLOCAL); - - std::stringstream ssh; - std::stringstream sss; - - if (bit) - { - ssh << GlobalV::global_out_dir << file_name + "-H-bit"; - sss << GlobalV::global_out_dir << file_name + "-S-bit"; - } - else - { - if (GlobalV::out_app_flag) - { - ssh << GlobalV::global_out_dir << file_name + "-H"; - sss << GlobalV::global_out_dir << file_name + "-S"; - } - else - { - ssh << GlobalV::global_out_dir << istep << "_" << file_name + "-H"; - sss << GlobalV::global_out_dir << istep << "_" << file_name + "-S"; - } - } - - if (bit) - { -#ifdef __MPI - FILE *g1 = nullptr; - FILE *g2 = nullptr; - - if (GlobalV::DRANK == 0) - { - g1 = fopen(ssh.str().c_str(), "wb"); - g2 = fopen(sss.str().c_str(), "wb"); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g1); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g2); - } - - int ir, ic; - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - std::complex *lineH = new std::complex[GlobalV::NLOCAL]; - std::complex *lineS = new std::complex[GlobalV::NLOCAL]; - ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL); - ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL); - - ir = pv.global2local_row(i); - if (ir >= 0) - { - // data collection - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - ic = pv.global2local_col(j); - if (ic >= 0) - { - int iic; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { - iic = ir + ic * pv.nrow; - } - else - { - iic = ir * pv.ncol + ic; - } - // lineH[j-i] = H[ir*pv.ncol+ic]; - // lineS[j-i] = S[ir*pv.ncol+ic]; - lineH[j] = H[iic]; - lineS[j] = S[iic]; - } - } - } - else - { - // do nothing - } - - Parallel_Reduce::reduce_pool(lineH, GlobalV::NLOCAL); - Parallel_Reduce::reduce_pool(lineS, GlobalV::NLOCAL); - - if (GlobalV::DRANK == 0) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - fwrite(&lineH[j], sizeof(std::complex), 1, g1); - fwrite(&lineS[j], sizeof(std::complex), 1, g2); - } - } - delete[] lineH; - delete[] lineS; - - MPI_Barrier(DIAG_WORLD); - } - - if (GlobalV::DRANK == 0) - { - fclose(g1); - fclose(g2); - } -#else - FILE *g1 = fopen(ssh.str().c_str(), "wb"); - FILE *g2 = fopen(sss.str().c_str(), "wb"); - - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g1); - fwrite(&GlobalV::NLOCAL, sizeof(int), 1, g2); - - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - fwrite(&H[i * GlobalV::NLOCAL + j], sizeof(std::complex), 1, g1); - fwrite(&S[i * GlobalV::NLOCAL + j], sizeof(std::complex), 1, g2); - } - } - fclose(g1); - fclose(g2); -#endif - } // end bit - else - { -#ifdef __MPI - std::ofstream g1; - std::ofstream g2; - - if (GlobalV::DRANK == 0) - { - if (GlobalV::out_app_flag) - { - g1.open(ssh.str().c_str(), std::ofstream::app); - g2.open(sss.str().c_str(), std::ofstream::app); - } - else - { - g1.open(ssh.str().c_str()); - g2.open(sss.str().c_str()); - } - g1 << GlobalV::NLOCAL; - g2 << GlobalV::NLOCAL; - } - - int ir, ic; - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - std::complex *lineH = new std::complex[GlobalV::NLOCAL]; - std::complex *lineS = new std::complex[GlobalV::NLOCAL]; - ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL); - ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL); - - ir = pv.global2local_row(i); - if (ir >= 0) - { - // data collection - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - ic = pv.global2local_col(j); - if (ic >= 0) - { - int iic; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { - iic = ir + ic * pv.nrow; - } - else - { - iic = ir * pv.ncol + ic; - } - // lineH[j-i] = H[ir*pv.ncol+ic]; - // lineS[j-i] = S[ir*pv.ncol+ic]; - lineH[j] = H[iic]; - lineS[j] = S[iic]; - } - } - } - else - { - // do nothing - } - - Parallel_Reduce::reduce_pool(lineH, GlobalV::NLOCAL); - Parallel_Reduce::reduce_pool(lineS, GlobalV::NLOCAL); - - if (GlobalV::DRANK == 0) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - g1 << " " << lineH[j]; - g2 << " " << lineS[j]; - } - g1 << std::endl; - g2 << std::endl; - } - delete[] lineH; - delete[] lineS; - } - - // if (GlobalV::DRANK==0); - if (GlobalV::DRANK == 0) // Peize Lin delete ; at 2020.01.31 - { - g1.close(); - g2.close(); - } -#else - if (GlobalV::out_app_flag) - { - std::ofstream g1(ssh.str().c_str(), ofstream::app); - std::ofstream g2(sss.str().c_str(), ofstream::app); - } - else - { - std::ofstream g1(ssh.str().c_str()); - std::ofstream g2(sss.str().c_str()); - } - - g1 << GlobalV::NLOCAL; - g2 << GlobalV::NLOCAL; - - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - g1 << " " << H[i * GlobalV::NLOCAL + j]; - g2 << " " << S[i * GlobalV::NLOCAL + j]; - } - g1 << std::endl; - g2 << std::endl; - } - g1.close(); - g2.close(); -#endif - } - - ModuleBase::timer::tick("ModuleIO", "save_HS_bit"); - return; -} - -// void ModuleIO::save_HSR_tr(const int Rx, const int Ry, const int Rz, const double *H, const double *S) -void ModuleIO::save_HSR_tr(const int current_spin, LCAO_Matrix &lm) -// void ModuleIO::save_HSR_tr(void) -{ - ModuleBase::TITLE("ModuleIO", "save_HSR_tr"); - ModuleBase::timer::tick("ModuleIO", "save_HSR_tr"); - - std::stringstream ssh; - std::stringstream sss; - - ssh << GlobalV::global_out_dir << "data-HR-tr_SPIN" << current_spin; - sss << GlobalV::global_out_dir << "data-SR-tr_SPIN" << current_spin; - // ssh << GlobalV::global_out_dir << "data-HR-tr_SPIN"; - // sss << GlobalV::global_out_dir << "data-SR-tr_SPIN"; - -#ifdef __MPI - std::ofstream g1; - std::ofstream g2; - - if (GlobalV::DRANK == 0) - { - g1.open(ssh.str().c_str()); - g2.open(sss.str().c_str()); - g1 << "Matrix Dimension of H(R): " << GlobalV::NLOCAL << std::endl; - g2 << "Matrix Dimension of S(R): " << GlobalV::NLOCAL << std::endl; - } - - int R_x = GlobalC::GridD.getCellX(); - int R_y = GlobalC::GridD.getCellY(); - int R_z = GlobalC::GridD.getCellZ(); - - // std::cout<<"R_x: "< *lineH_soc = nullptr; - std::complex *lineS_soc = nullptr; - if (GlobalV::NSPIN != 4) - { - lineH = new double[GlobalV::NLOCAL]; - lineS = new double[GlobalV::NLOCAL]; - ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL); - ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL); - } - else - { - lineH_soc = new std::complex[GlobalV::NLOCAL]; - lineS_soc = new std::complex[GlobalV::NLOCAL]; - ModuleBase::GlobalFunc::ZEROS(lineH_soc, GlobalV::NLOCAL); - ModuleBase::GlobalFunc::ZEROS(lineS_soc, GlobalV::NLOCAL); - } - // ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL-i); - // ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL-i); - // ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL); - // ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL); - - ir = lm.ParaV->global2local_row(i); - if (ir >= 0) - { - // for(int j=i; jglobal2local_col(j); - if (ic >= 0) - { - // lineH[j-i] = H[ir*lm.ParaV->ncol+ic]; - // lineS[j-i] = S[ir*lm.ParaV->ncol+ic]; - int iic; - if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) - { - iic = ir + ic * lm.ParaV->nrow; - } - else - { - iic = ir * lm.ParaV->ncol + ic; - } - if (GlobalV::NSPIN != 4) - { - lineH[j] = lm.HR_tr[ix][iy][iz][iic]; - lineS[j] = lm.SlocR_tr[ix][iy][iz][iic]; - } - else - { - lineH_soc[j] = lm.HR_tr_soc[ix][iy][iz][iic]; - lineS_soc[j] = lm.SlocR_tr_soc[ix][iy][iz][iic]; - } - } - } - } - else - { - // do nothing - } - - // Parallel_Reduce::reduce_all(lineH,GlobalV::NLOCAL-i); - // Parallel_Reduce::reduce_all(lineS,GlobalV::NLOCAL-i); - if (GlobalV::NSPIN != 4) - { - Parallel_Reduce::reduce_all(lineH, GlobalV::NLOCAL); - Parallel_Reduce::reduce_all(lineS, GlobalV::NLOCAL); - } - else - { - Parallel_Reduce::reduce_all(lineH_soc, GlobalV::NLOCAL); - Parallel_Reduce::reduce_all(lineS_soc, GlobalV::NLOCAL); - } - - if (GlobalV::DRANK == 0) - { - // for(int j=i; j(0.0, lineH_soc[j].imag()); - if (std::abs(lineH_soc[j].imag()) < 1.0e-12) - lineH_soc[j] = std::complex(lineH_soc[j].real(), 0.0); - if (std::abs(lineS_soc[j].real()) < 1.0e-12) - lineS_soc[j] = std::complex(0.0, lineS_soc[j].imag()); - if (std::abs(lineS_soc[j].imag()) < 1.0e-12) - lineS_soc[j] = std::complex(lineS_soc[j].real(), 0.0); - g1 << " " << lineH_soc[j]; - g2 << " " << lineS_soc[j]; - } - } - g1 << std::endl; - g2 << std::endl; - } - if (GlobalV::NSPIN != 4) - { - delete[] lineH; - delete[] lineS; - } - else - { - delete[] lineH_soc; - delete[] lineS_soc; - } - } - /* - if(GlobalV::DRANK==0); - { - g1.close(); - g2.close(); - } - */ - } - } - } - // if(GlobalV::DRANK==0); - if (GlobalV::DRANK == 0) // Peize Lin delete ; at 2020.01.31 - { - g1.close(); - g2.close(); - } - -#else - std::ofstream g1(ssh.str().c_str()); - std::ofstream g2(sss.str().c_str()); - - g1 << GlobalV::NLOCAL; - g2 << GlobalV::NLOCAL; - - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = i; j < GlobalV::NLOCAL; j++) - { - // not correct for sequential version; need change - // g1 << " " << H[i*GlobalV::NLOCAL+j]; - // g2 << " " << S[i*GlobalV::NLOCAL+j]; - } - } - g1.close(); - g2.close(); -#endif - - ModuleBase::timer::tick("ModuleIO", "save_HSR_tr"); - return; -} diff --git a/source/module_io/write_HS.h b/source/module_io/write_HS.h index 26d1595cea..93365a48bf 100644 --- a/source/module_io/write_HS.h +++ b/source/module_io/write_HS.h @@ -10,65 +10,34 @@ // mohan add this file 2010-09-10 namespace ModuleIO { -void saving_HS(const int istep, - const double* Hloc, - const double* Sloc, - const bool bit, - const int& out_hs, - const std::string& file_name, - const Parallel_Orbitals& pv, - bool tri); -// tri : if set 1 , output triangle matrix ; if set 0 , output complete matrix . - -void save_HS_triangle(const int istep, - const double* H, - const double* S, - const bool bit, - const std::string& file_name, - const Parallel_Orbitals& pv); -void save_HS_complete(const int istep, - const double* H, - const double* S, - const bool bit, - const std::string& file_name, - const Parallel_Orbitals& pv); - -void save_HS_complex(const int istep, - const std::complex* H, - const std::complex* S, - const bool bit, - const std::string& file__name, - const Parallel_Orbitals& pv, - bool tri); - -void save_HSR_tr(const int current_spin, LCAO_Matrix& lm); // LiuXh add 2019-07-15 + /// @brief save a square matrix + /// @param[in] istep : the step of the calculation + /// @param[in] mat : the local matrix + /// @param[in] bit : true for binary, false for decimal + /// @param[in] tri : true for upper triangle, false for full matrix + /// @param[in] app : true for append, false for overwrite + /// @param[in] label : the symbol of the matrix, like "H", "S" + /// @param[in] file_name : the name of the output file + /// @param[in] pv : the 2d-block parallelization information + /// @param[in] drank : the rank of the current process + template + void save_mat(const int istep, + const T* mat, + const int dim, + const bool bit, + const bool tri, + const bool app, + const std::string label, + const std::string& file_name, + const Parallel_2D& pv, + const int drank); + + // comment out this function for not used + // void save_HSR_tr(const int current_spin, LCAO_Matrix& lm); // LiuXh add 2019-07-15 // mohan comment out 2021-02-10 // void save_HS_ccf(const int &iter, const int &Hnnz, const int *colptr_H, const int *rowind_H, // const double *nzval_H, const double *nzval_S, bool bit); - -void saving_HS(const int istep, - std::complex* Hloc, - std::complex* Sloc, - bool bit, - const int& out_hs, - const std::string& file__name, - const Parallel_Orbitals& pv, - bool tri); // LiuXh, 2017-03-21 - -void save_HS_complex_triangle(const int istep, - std::complex* H, - std::complex* S, - const bool bit, - const std::string& file__name, - const Parallel_Orbitals& pv); // LiuXh, 2017-03-21 - -void save_HS_complex_complete(const int istep, - std::complex* H, - std::complex* S, - const bool bit, - const std::string& file__name, - const Parallel_Orbitals& pv); // LiuXh, 2017-03-21 } - +#include "write_HS.hpp" #endif diff --git a/source/module_io/write_HS.hpp b/source/module_io/write_HS.hpp new file mode 100644 index 0000000000..7ab6c1f02c --- /dev/null +++ b/source/module_io/write_HS.hpp @@ -0,0 +1,486 @@ +#include "write_HS.h" + +#include "module_base/parallel_reduce.h" +#include "module_base/timer.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" + + +/* +void ModuleIO::save_HS_ccf(const int &iter, const int &Hnnz, const int *colptr_H, const int *rowind_H, + const double *nzval_H, const double *nzval_S, bool bit) +{ + ModuleBase::TITLE("ModuleIO","save_HS_ccf"); + + if(GlobalV::DRANK!=0)return; + + std::stringstream ssh; + std::stringstream sss; + + if(bit) + { + ssh << GlobalV::global_out_dir << "H_bit.ccf"; + sss << GlobalV::global_out_dir << "S_bit.ccf"; + } + else + { + // mohan update 2021-02-10 + ssh << GlobalV::global_out_dir << "H" << ELEC_scf::iter << "_" << iter+1 << ".ccf"; + sss << GlobalV::global_out_dir << "S" << ELEC_scf::iter << "_" << iter+1 << ".ccf"; + } + + if(bit) + { + FILE *g1 = fopen(ssh.str().c_str(),"wb"); + FILE *g2 = fopen(sss.str().c_str(),"wb"); + + fwrite(&GlobalV::NLOCAL,sizeof(int),1,g1); + fwrite(&Hnnz,sizeof(int),1,g1); + fwrite(&GlobalV::NLOCAL,sizeof(int),1,g2); + fwrite(&Hnnz,sizeof(int),1,g2); + + fclose(g1); + fclose(g2); + } + + + if(!bit) + { + std::ofstream g1(ssh.str().c_str()); + std::ofstream g2(sss.str().c_str()); + + g1 << GlobalV::NLOCAL << " " << Hnnz << std::endl; + g2 << GlobalV::NLOCAL << " " << Hnnz << std::endl; + + for(int i=0; i +void ModuleIO::save_mat(const int istep, + const T* mat, + const int dim, + const bool bit, + const bool tri, + const bool app, + const std::string label, + const std::string& file_name, + const Parallel_2D& pv, + const int drank) +{ + ModuleBase::TITLE("ModuleIO", "save_mat"); + ModuleBase::timer::tick("ModuleIO", "save_mat"); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Dimension of " + label + " : ", dim); + + std::stringstream ss; + + if (bit)ss << GlobalV::global_out_dir << file_name + "-" + label + "-bit"; + else + { + if (app || istep < 0) + ss << GlobalV::global_out_dir << file_name + "-" + label; + else + ss << GlobalV::global_out_dir << istep << "_" << file_name + "-" + label; + } + if (bit) + { +#ifdef __MPI + FILE* g = nullptr; + + if (drank == 0) + { + g = fopen(ss.str().c_str(), "wb"); + fwrite(&dim, sizeof(int), 1, g); + } + + int ir, ic; + for (int i = 0; i < dim; ++i) + { + T* line = new T[tri ? dim - i : dim]; + ModuleBase::GlobalFunc::ZEROS(line, tri ? dim - i : dim); + + ir = pv.global2local_row(i); + if (ir >= 0) + { + // data collection + for (int j = (tri ? i : 0); j < dim; ++j) + { + ic = pv.global2local_col(j); + if (ic >= 0) + { + int iic; + if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) + { + iic = ir + ic * pv.nrow; + } + else + { + iic = ir * pv.ncol + ic; + } + line[tri ? j - i : j] = mat[iic]; + } + } + } + + Parallel_Reduce::reduce_all(line, tri ? dim - i : dim); + + if (drank == 0) + { + for (int j = (tri ? i : 0); j < dim; ++j) + { + fwrite(&line[tri ? j - i : j], sizeof(T), 1, g); + } + } + delete[] line; + + MPI_Barrier(DIAG_WORLD); + } + + if (drank == 0) + fclose(g); +#else + FILE* g = fopen(ss.str().c_str(), "wb"); + + fwrite(&dim, sizeof(int), 1, g); + + for (int i = 0; i < dim; i++) + { + for (int j = (tri ? i : 0); j < dim; j++) + { + fwrite(&mat[i * dim + j], sizeof(T), 1, g); + } + } + fclose(g); +#endif + } // end bit + else + { + std::ofstream g; +#ifdef __MPI + if (drank == 0) + { + if (app) + g.open(ss.str().c_str(), std::ofstream::app); + else + g.open(ss.str().c_str()); + g << dim; + } + + int ir, ic; + for (int i = 0; i < dim; i++) + { + T* line = new T[tri ? dim - i : dim]; + ModuleBase::GlobalFunc::ZEROS(line, tri ? dim - i : dim); + + ir = pv.global2local_row(i); + if (ir >= 0) + { + // data collection + for (int j = (tri ? i : 0); j < dim; ++j) + { + ic = pv.global2local_col(j); + if (ic >= 0) + { + int iic; + if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) + { + iic = ir + ic * pv.nrow; + } + else + { + iic = ir * pv.ncol + ic; + } + line[tri ? j - i : j] = mat[iic]; + } + } + } + + Parallel_Reduce::reduce_all(line, tri ? dim - i : dim); + + if (drank == 0) + { + for (int j = (tri ? i : 0); j < dim; j++) g << " " << line[tri ? j - i : j]; + g << std::endl; + } + delete[] line; + } + + if (drank == 0) // Peize Lin delete ; at 2020.01.31 + g.close(); +#else + if (app) + std::ofstream g(ss.str().c_str(), std::ofstream::app); + else + std::ofstream g(ss.str().c_str()); + + g << dim; + + for (int i = 0; i < dim; i++) + { + for (int j = (tri ? i : 0); j < dim; j++) + { + g << " " << mat[i * dim + j]; + } + g << std::endl; + } + g.close(); +#endif + } + ModuleBase::timer::tick("ModuleIO", "save_mat"); + return; +} + + +/* comment out this function for not used +// void ModuleIO::save_HSR_tr(const int Rx, const int Ry, const int Rz, const double *H, const double *S) +void ModuleIO::save_HSR_tr(const int current_spin, LCAO_Matrix &lm) +// void ModuleIO::save_HSR_tr(void) +{ + ModuleBase::TITLE("ModuleIO", "save_HSR_tr"); + ModuleBase::timer::tick("ModuleIO", "save_HSR_tr"); + + std::stringstream ssh; + std::stringstream sss; + + ssh << GlobalV::global_out_dir << "data-HR-tr_SPIN" << current_spin; + sss << GlobalV::global_out_dir << "data-SR-tr_SPIN" << current_spin; + // ssh << GlobalV::global_out_dir << "data-HR-tr_SPIN"; + // sss << GlobalV::global_out_dir << "data-SR-tr_SPIN"; + +#ifdef __MPI + std::ofstream g1; + std::ofstream g2; + + if (GlobalV::DRANK == 0) + { + g1.open(ssh.str().c_str()); + g2.open(sss.str().c_str()); + g1 << "Matrix Dimension of H(R): " << GlobalV::NLOCAL << std::endl; + g2 << "Matrix Dimension of S(R): " << GlobalV::NLOCAL << std::endl; + } + + int R_x = GlobalC::GridD.getCellX(); + int R_y = GlobalC::GridD.getCellY(); + int R_z = GlobalC::GridD.getCellZ(); + + // std::cout<<"R_x: "< *lineH_soc = nullptr; + std::complex *lineS_soc = nullptr; + if (GlobalV::NSPIN != 4) + { + lineH = new double[GlobalV::NLOCAL]; + lineS = new double[GlobalV::NLOCAL]; + ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL); + ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL); + } + else + { + lineH_soc = new std::complex[GlobalV::NLOCAL]; + lineS_soc = new std::complex[GlobalV::NLOCAL]; + ModuleBase::GlobalFunc::ZEROS(lineH_soc, GlobalV::NLOCAL); + ModuleBase::GlobalFunc::ZEROS(lineS_soc, GlobalV::NLOCAL); + } + // ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL-i); + // ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL-i); + // ModuleBase::GlobalFunc::ZEROS(lineH, GlobalV::NLOCAL); + // ModuleBase::GlobalFunc::ZEROS(lineS, GlobalV::NLOCAL); + + ir = lm.ParaV->global2local_row(i); + if (ir >= 0) + { + // for(int j=i; jglobal2local_col(j); + if (ic >= 0) + { + // lineH[j-i] = H[ir*lm.ParaV->ncol+ic]; + // lineS[j-i] = S[ir*lm.ParaV->ncol+ic]; + int iic; + if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER()) + { + iic = ir + ic * lm.ParaV->nrow; + } + else + { + iic = ir * lm.ParaV->ncol + ic; + } + if (GlobalV::NSPIN != 4) + { + lineH[j] = lm.HR_tr[ix][iy][iz][iic]; + lineS[j] = lm.SlocR_tr[ix][iy][iz][iic]; + } + else + { + lineH_soc[j] = lm.HR_tr_soc[ix][iy][iz][iic]; + lineS_soc[j] = lm.SlocR_tr_soc[ix][iy][iz][iic]; + } + } + } + } + else + { + // do nothing + } + + // Parallel_Reduce::reduce_all(lineH,GlobalV::NLOCAL-i); + // Parallel_Reduce::reduce_all(lineS,GlobalV::NLOCAL-i); + if (GlobalV::NSPIN != 4) + { + Parallel_Reduce::reduce_all(lineH, GlobalV::NLOCAL); + Parallel_Reduce::reduce_all(lineS, GlobalV::NLOCAL); + } + else + { + Parallel_Reduce::reduce_all(lineH_soc, GlobalV::NLOCAL); + Parallel_Reduce::reduce_all(lineS_soc, GlobalV::NLOCAL); + } + + if (GlobalV::DRANK == 0) + { + // for(int j=i; j(0.0, lineH_soc[j].imag()); + if (std::abs(lineH_soc[j].imag()) < 1.0e-12) + lineH_soc[j] = std::complex(lineH_soc[j].real(), 0.0); + if (std::abs(lineS_soc[j].real()) < 1.0e-12) + lineS_soc[j] = std::complex(0.0, lineS_soc[j].imag()); + if (std::abs(lineS_soc[j].imag()) < 1.0e-12) + lineS_soc[j] = std::complex(lineS_soc[j].real(), 0.0); + g1 << " " << lineH_soc[j]; + g2 << " " << lineS_soc[j]; + } + } + g1 << std::endl; + g2 << std::endl; + } + if (GlobalV::NSPIN != 4) + { + delete[] lineH; + delete[] lineS; + } + else + { + delete[] lineH_soc; + delete[] lineS_soc; + } + } + + // if(GlobalV::DRANK==0); + // { + // g1.close(); + // g2.close(); + // } + + } + } + } + // if(GlobalV::DRANK==0); + if (GlobalV::DRANK == 0) // Peize Lin delete ; at 2020.01.31 + { + g1.close(); + g2.close(); + } + +#else + std::ofstream g1(ssh.str().c_str()); + std::ofstream g2(sss.str().c_str()); + + g1 << GlobalV::NLOCAL; + g2 << GlobalV::NLOCAL; + + for (int i = 0; i < GlobalV::NLOCAL; i++) + { + for (int j = i; j < GlobalV::NLOCAL; j++) + { + // not correct for sequential version; need change + // g1 << " " << H[i*GlobalV::NLOCAL+j]; + // g2 << " " << S[i*GlobalV::NLOCAL+j]; + } + } + g1.close(); + g2.close(); +#endif + + ModuleBase::timer::tick("ModuleIO", "save_HSR_tr"); + return; +} +*/ From 9fc2e2fd40827307affd3176be983c3d44de5f7c Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 4 Jan 2024 21:12:06 +0800 Subject: [PATCH 2/3] feature: write Vxc with test --- docs/advanced/input_files/input-main.md | 8 + source/module_base/global_variable.cpp | 2 + source/module_base/global_variable.h | 2 + source/module_esolver/esolver_ks_lcao.cpp | 7 + source/module_io/input.cpp | 6 + source/module_io/input.h | 1 + source/module_io/input_conv.cpp | 1 + source/module_io/parameter_pool.cpp | 4 + source/module_io/test/input_conv_test.cpp | 1 + source/module_io/test/input_test.cpp | 2 + source/module_io/test/input_test_para.cpp | 1 + source/module_io/test/write_input_test.cpp | 1 + source/module_io/write_Vxc.hpp | 210 +++++++++++++++++++++ source/module_io/write_input.cpp | 1 + tests/integrate/207_NO_KP_OXC/INPUT | 30 +++ tests/integrate/207_NO_KP_OXC/KPT | 4 + tests/integrate/207_NO_KP_OXC/STRU | 22 +++ tests/integrate/207_NO_KP_OXC/jd | 1 + tests/integrate/207_NO_KP_OXC/k-1-Vxc.ref | 6 + tests/integrate/207_NO_KP_OXC/result.ref | 4 + tests/integrate/307_NO_GO_OXC/INPUT | 36 ++++ tests/integrate/307_NO_GO_OXC/KPT | 4 + tests/integrate/307_NO_GO_OXC/STRU | 22 +++ tests/integrate/307_NO_GO_OXC/jd | 1 + tests/integrate/307_NO_GO_OXC/k-0-Vxc.ref | 6 + tests/integrate/307_NO_GO_OXC/result.ref | 7 + tests/integrate/CASES | 2 + tests/integrate/tools/catch_properties.sh | 13 ++ 28 files changed, 405 insertions(+) create mode 100644 source/module_io/write_Vxc.hpp create mode 100644 tests/integrate/207_NO_KP_OXC/INPUT create mode 100644 tests/integrate/207_NO_KP_OXC/KPT create mode 100644 tests/integrate/207_NO_KP_OXC/STRU create mode 100644 tests/integrate/207_NO_KP_OXC/jd create mode 100644 tests/integrate/207_NO_KP_OXC/k-1-Vxc.ref create mode 100644 tests/integrate/207_NO_KP_OXC/result.ref create mode 100644 tests/integrate/307_NO_GO_OXC/INPUT create mode 100644 tests/integrate/307_NO_GO_OXC/KPT create mode 100644 tests/integrate/307_NO_GO_OXC/STRU create mode 100644 tests/integrate/307_NO_GO_OXC/jd create mode 100644 tests/integrate/307_NO_GO_OXC/k-0-Vxc.ref create mode 100644 tests/integrate/307_NO_GO_OXC/result.ref diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 5822142af2..3119d96e42 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -141,6 +141,7 @@ - [out\_mat\_hs2](#out_mat_hs2) - [out\_mat\_t](#out_mat_t) - [out\_mat\_dh](#out_mat_dh) + - [out\_mat\_xc](#out_mat_xc) - [out\_app\_flag](#out_app_flag) - [out\_interval](#out_interval) - [out\_element\_info](#out_element_info) @@ -1569,6 +1570,13 @@ These variables are used to control the output of properties. - **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix (in Ry/Bohr). The format will be the same as the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ as mentioned in [out_mat_hs2](#out_mat_hs2). The name of the files will be `data-dHRx-sparse_SPIN0.csr` and so on. Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag). - **Default**: False +### out_mat_xc + +- **Type**: Boolean +- **Availability**: Numerical atomic orbital basis +- **Description**: Whether to print the upper triangular part of the exchange-correlation matrices in **Kohn-Sham orbital representation** (unit: Ry) for each k point into files in the directory `OUT.${suffix}`, which is useful for the subsequent GW calculation . The files are named `k-$k-Vxc`, the meaning of `$k`corresponding to k point and spin is same as [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs). +- **Default**: False + ### out_app_flag - **Type**: Boolean diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp index 09ebde269a..48e2437b4d 100644 --- a/source/module_base/global_variable.cpp +++ b/source/module_base/global_variable.cpp @@ -277,6 +277,8 @@ double nelec = 0; bool out_bandgap = false; // QO added for bandgap printing int out_interval = 1; // convert from out_hsR_interval liuyu 2023-04-18 +bool out_mat_xc = false; // output Vxc in KS-wfc representation for GW calculation + //========================================================== // Deltaspin related //========================================================== diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h index 4b00d92d85..d482038868 100644 --- a/source/module_base/global_variable.h +++ b/source/module_base/global_variable.h @@ -308,6 +308,8 @@ extern double nelec; extern bool out_bandgap; extern int out_interval; +extern bool out_mat_xc; // output Vxc in KS-wfc representation for GW calculation + // Deltaspin related extern bool sc_mag_switch; // 0: no deltaspin; 1: constrain atomic magnetic moments; extern bool decay_grad_switch; // 0: decay grad will be set to zero; 1: with decay grad set for some elements diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 4c48eadcc6..2d318827d4 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -8,6 +8,7 @@ #include "module_io/write_proj_band_lcao.h" #include "module_io/output_log.h" #include "module_io/to_qo.h" +#include "module_io/write_Vxc.hpp" //--------------temporary---------------------------- #include "module_base/global_function.h" @@ -821,6 +822,12 @@ namespace ModuleESolver } } + bool out_exc = true; // tmp, add parameter! + if (GlobalV::out_mat_xc) + ModuleIO::write_Vxc(GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::DRANK, + *this->psi, GlobalC::ucell, this->sf, *this->pw_rho, *this->pw_rhod, GlobalC::ppcell.vloc, + *this->pelec->charge, this->UHM, this->LM, this->LOC, this->kv, this->pelec->wg, GlobalC::GridD); + #ifdef __EXX if (GlobalC::exx_info.info_global.cal_exx) // Peize Lin add if 2022.11.14 { diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 006a4a4f92..4eb177810a 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -337,6 +337,7 @@ void Input::Default(void) out_band = 0; out_proj_band = 0; out_mat_hs = 0; + out_mat_xc = 0; cal_syns = 0; dmax = 0.01; out_mat_hs2 = 0; // LiuXh add 2019-07-15 @@ -1396,6 +1397,10 @@ bool Input::Read(const std::string &fn) { read_bool(ifs, out_mat_dh); } + else if (strcmp("out_mat_xc", word) == 0) + { + read_bool(ifs, out_mat_xc); + } else if (strcmp("out_interval", word) == 0) { read_value(ifs, out_interval); @@ -3279,6 +3284,7 @@ void Input::Bcast() Parallel_Common::bcast_bool(out_mat_hs2); // LiuXh add 2019-07-15 Parallel_Common::bcast_bool(out_mat_t); Parallel_Common::bcast_bool(out_mat_dh); + Parallel_Common::bcast_bool(out_mat_xc); Parallel_Common::bcast_bool(out_mat_r); // jingan add 2019-8-14 Parallel_Common::bcast_int(out_wfc_lcao); Parallel_Common::bcast_bool(out_alllog); diff --git a/source/module_io/input.h b/source/module_io/input.h index d9dd3b982a..e2f72e8357 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -265,6 +265,7 @@ class Input bool out_band; // band calculation pengfei 2014-10-13 bool out_proj_band; // projected band structure calculation jiyy add 2022-05-11 bool out_mat_hs; // output H matrix and S matrix in local basis. + bool out_mat_xc; // output exchange-correlation matrix in KS-orbital representation. bool cal_syns; // calculate asynchronous S matrix to output double dmax; // maximum displacement of all atoms in one step (bohr) bool out_mat_hs2; // LiuXh add 2019-07-16, output H(R) matrix and S(R) matrix in local basis. diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 9871ebe254..c3af6c5c89 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -638,6 +638,7 @@ void Input_Conv::Convert(void) hsolver::HSolverLCAO>::out_mat_hsR = INPUT.out_mat_hs2; // LiuXh add 2019-07-16 hsolver::HSolverLCAO>::out_mat_t = INPUT.out_mat_t; hsolver::HSolverLCAO>::out_mat_dh = INPUT.out_mat_dh; + GlobalV::out_mat_xc = INPUT.out_mat_xc; if (GlobalV::GAMMA_ONLY_LOCAL) { elecstate::ElecStateLCAO::out_wfc_lcao = INPUT.out_wfc_lcao; diff --git a/source/module_io/parameter_pool.cpp b/source/module_io/parameter_pool.cpp index 45d2262f94..d66fba4996 100644 --- a/source/module_io/parameter_pool.cpp +++ b/source/module_io/parameter_pool.cpp @@ -927,6 +927,10 @@ bool input_parameters_set(std::map input_parameters { INPUT.out_mat_hs = *static_cast(input_parameters["out_mat_hs"].get()); } + else if (input_parameters.count("out_mat_xc") != 0) + { + INPUT.out_mat_xc = *static_cast(input_parameters["out_mat_xc"].get()); + } else if (input_parameters.count("cal_syns") != 0) { INPUT.cal_syns = *static_cast(input_parameters["cal_syns"].get()); diff --git a/source/module_io/test/input_conv_test.cpp b/source/module_io/test/input_conv_test.cpp index 6f0493167d..7c1e4b18a8 100644 --- a/source/module_io/test/input_conv_test.cpp +++ b/source/module_io/test/input_conv_test.cpp @@ -153,6 +153,7 @@ TEST_F(InputConvTest, Conv) EXPECT_EQ(hsolver::HSolverLCAO>::out_mat_t, false); EXPECT_EQ(hsolver::HSolverLCAO::out_mat_dh, INPUT.out_mat_dh); EXPECT_EQ(hsolver::HSolverLCAO>::out_mat_dh, INPUT.out_mat_dh); + EXPECT_EQ(GlobalV::out_mat_xc, false); EXPECT_EQ(GlobalV::out_interval, 1); EXPECT_EQ(elecstate::ElecStateLCAO::out_wfc_lcao, false); EXPECT_EQ(berryphase::berry_phase_flag, false); diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index 32cbec3b24..8843d03454 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -180,6 +180,7 @@ TEST_F(InputTest, Default) EXPECT_EQ(INPUT.out_proj_band,0); EXPECT_EQ(INPUT.out_mat_hs,0); EXPECT_EQ(INPUT.out_mat_hs2,0); + EXPECT_EQ(INPUT.out_mat_xc, 0); EXPECT_EQ(INPUT.out_interval,1); EXPECT_EQ(INPUT.out_app_flag,1); EXPECT_EQ(INPUT.out_mat_r,0); @@ -542,6 +543,7 @@ TEST_F(InputTest, Read) EXPECT_EQ(INPUT.out_proj_band,0); EXPECT_EQ(INPUT.out_mat_hs,0); EXPECT_EQ(INPUT.out_mat_hs2,0); + EXPECT_EQ(INPUT.out_mat_xc, 0); EXPECT_EQ(INPUT.out_interval,1); EXPECT_EQ(INPUT.out_app_flag,0); EXPECT_EQ(INPUT.out_mat_r,0); diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp index ace61d068a..dcdf01ab41 100644 --- a/source/module_io/test/input_test_para.cpp +++ b/source/module_io/test/input_test_para.cpp @@ -184,6 +184,7 @@ TEST_F(InputParaTest, Bcast) EXPECT_EQ(INPUT.out_proj_band, 0); EXPECT_EQ(INPUT.out_mat_hs, 0); EXPECT_EQ(INPUT.out_mat_hs2, 0); + EXPECT_EQ(INPUT.out_mat_xc, 0); EXPECT_EQ(INPUT.out_interval, 1); EXPECT_EQ(INPUT.out_app_flag, 1); EXPECT_EQ(INPUT.out_mat_r, 0); diff --git a/source/module_io/test/write_input_test.cpp b/source/module_io/test/write_input_test.cpp index 439d163834..d61133715d 100644 --- a/source/module_io/test/write_input_test.cpp +++ b/source/module_io/test/write_input_test.cpp @@ -323,6 +323,7 @@ TEST_F(write_input, LCAO5) EXPECT_THAT(output, testing::HasSubstr("lcao_rmax 30 #max R for 1D two-center integration table")); EXPECT_THAT(output, testing::HasSubstr("out_mat_hs 0 #output H and S matrix")); + EXPECT_THAT(output, testing::HasSubstr("out_mat_xc 0 #output exchange-correlation matrix in KS-orbital representation")); EXPECT_THAT(output, testing::HasSubstr("out_mat_hs2 0 #output H(R) and S(R) matrix")); EXPECT_THAT(output, testing::HasSubstr("out_mat_dh 0 #output of derivative of H(R) matrix")); EXPECT_THAT( diff --git a/source/module_io/write_Vxc.hpp b/source/module_io/write_Vxc.hpp new file mode 100644 index 0000000000..7bb63329b2 --- /dev/null +++ b/source/module_io/write_Vxc.hpp @@ -0,0 +1,210 @@ +#include "module_psi/psi.h" +#include "module_hamilt_lcao/module_gint/gint_tools.h" +#include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" +#include "module_elecstate/potentials/pot_xc.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" +#include "module_base/scalapack_connector.h" +#include "module_base/parallel_reduce.h" +template struct TGint; +template <> +struct TGint { + using type = Gint_Gamma; +}; +template <> +struct TGint> { + using type = Gint_k; +}; +namespace ModuleIO +{ + inline void gint_vl(Gint_Gamma& gg, Gint_inout& io, LCAO_Matrix& lm) { gg.cal_vlocal(&io, &lm, false); }; + inline void gint_vl(Gint_k& gk, Gint_inout& io, LCAO_Matrix& lm, ModuleBase::matrix& wg) { gk.cal_gint(&io); }; + + void set_para2d_MO(const Parallel_Orbitals& pv, const int nbands, Parallel_2D& p2d) + { + std::ofstream ofs; +#ifdef __MPI + int dsize; + MPI_Comm_size(MPI_COMM_WORLD, &dsize); + p2d.set_block_size(pv.nb); + p2d.set_proc_dim(dsize); + p2d.comm_2D = pv.comm_2D; + p2d.blacs_ctxt = pv.blacs_ctxt; + p2d.set_local2global(nbands, nbands, ofs, ofs); + p2d.set_desc(nbands, nbands, p2d.get_row_size(), false); + p2d.set_global2local(nbands, nbands, true, ofs); +#else + p2d.set_proc_dim(1); + p2d.set_serial(nbands, nbands); + p2d.set_global2local(nbands, nbands, false, ofs); +#endif + } + + std::vector> cVc(std::complex* V, std::complex* c, int nbasis, int nbands, const Parallel_Orbitals& pv, const Parallel_2D& p2d) + { + std::vector> Vc(pv.nloc_wfc, 0.0); + char transa = 'N'; + char transb = 'N'; + const std::complex alpha(1.0, 0.0); + const std::complex beta(0.0, 0.0); +#ifdef __MPI + const int i1 = 1; + pzgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &i1, &i1, pv.desc, c, &i1, &i1, pv.desc_wfc, &beta, Vc.data(), &i1, &i1, pv.desc_wfc); +#else + zgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, c, &nbasis, &beta, Vc.data(), &nbasis); +#endif + std::vector> cVc(p2d.nloc, 0.0); + transa = 'C'; +#ifdef __MPI + pzgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &i1, &i1, pv.desc_wfc, Vc.data(), &i1, &i1, pv.desc_wfc, &beta, cVc.data(), &i1, &i1, p2d.desc); +#else + zgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &nbasis, Vc.data(), &nbasis, &beta, cVc.data(), &nbasis); +#endif + return cVc; + } + + std::vector cVc(double* V, double* c, int nbasis, int nbands, const Parallel_Orbitals& pv, const Parallel_2D& p2d) + { + std::vector Vc(pv.nloc_wfc, 0.0); + char transa = 'N'; + char transb = 'N'; + const double alpha = 1.0; + const double beta = 0.0; +#ifdef __MPI + const int i1 = 1; + pdgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &i1, &i1, pv.desc, c, &i1, &i1, pv.desc_wfc, &beta, Vc.data(), &i1, &i1, pv.desc_wfc); +#else + dgemm_(&transa, &transb, &nbasis, &nbands, &nbasis, &alpha, V, &nbasis, c, &nbasis, &beta, Vc.data(), &nbasis); +#endif + std::vector cVc(p2d.nloc, 0.0); + transa = 'T'; +#ifdef __MPI + pdgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &i1, &i1, pv.desc_wfc, Vc.data(), &i1, &i1, pv.desc_wfc, &beta, cVc.data(), &i1, &i1, p2d.desc); +#else + dgemm_(&transa, &transb, &nbands, &nbands, &nbasis, &alpha, c, &nbasis, Vc.data(), &nbasis, &beta, cVc.data(), &nbasis); +#endif + return cVc; + } + + inline double get_real(const std::complex& c) { return c.real(); } + inline double get_real(const double& d) { return d; } + + template + double all_band_energy(const int ik, const std::vector& mat_mo, const Parallel_2D& p2d, const ModuleBase::matrix& wg) + { + double e = 0.0; + for (int i = 0;i < p2d.get_row_size();++i) + for (int j = 0;j < p2d.get_col_size();++j) + if (p2d.local2global_row(i) == p2d.local2global_col(j)) + e += get_real(mat_mo[j * p2d.get_row_size() + i]) * wg(ik, p2d.local2global_row(i)); + Parallel_Reduce::reduce_all(e); + return e; + } + + template + void set_gint_pointer(LCAO_Hamilt& uhm, typename TGint::type*& gint); + + template <> + void set_gint_pointer(LCAO_Hamilt& uhm, typename TGint::type*& gint) + { + gint = &uhm.GG; + } + template <> + void set_gint_pointer>(LCAO_Hamilt& uhm, typename TGint>::type*& gint) + { + gint = &uhm.GK; + } + + + /// @brief write the Vxc matrix in KS orbital representation, usefull for GW calculation + template + void write_Vxc(int nspin, int nbasis, int drank, const psi::Psi& psi, const UnitCell& ucell, Structure_Factor& sf, + const ModulePW::PW_Basis& rho_basis, const ModulePW::PW_Basis& rhod_basis, const ModuleBase::matrix& vloc, + const Charge& chg, LCAO_Hamilt& uhm, LCAO_Matrix& lm, Local_Orbital_Charge& loc, + const K_Vectors& kv, const ModuleBase::matrix& wg, Grid_Driver& gd) + { + ModuleBase::TITLE("ModuleIO", "write_Vxc"); + const Parallel_Orbitals* pv = lm.ParaV; + int nbands = wg.nc; + // 1. real-space xc potential + // ModuleBase::matrix vr_xc(nspin, chg.nrxx); + double etxc = 0.0; + double vtxc = 0.0; + // elecstate::PotXC* potxc(&rho_basis, &etxc, vtxc, nullptr); + // potxc.cal_v_eff(&chg, &ucell, vr_xc); + elecstate::Potential* potxc = new elecstate::Potential(&rhod_basis, &rho_basis, &ucell, &vloc, &sf, &etxc, &vtxc); + std::vector compnents_list = { "xc" }; + potxc->pot_register(compnents_list); + potxc->update_from_charge(&chg, &ucell); + + // 2. allocate AO-matrix + // R (the number of hR: 1 for nspin=1, 4; 2 for nspin=2) + int nspin0 = (nspin == 2) ? 2 : 1; + std::vector> vxcs_R_ao(nspin0, hamilt::HContainer(pv)); + for (int is = 0;is < nspin0;++is) vxcs_R_ao[is].set_zero(); + // k (size for each k-point) + std::vector vxc_k_ao(pv->nloc); + + // 3. allocate operators and contribute HR + // op (corresponding to hR) + typename TGint::type* gint = nullptr; + set_gint_pointer(uhm, gint); + std::vector>*> vxcs_op_ao(nspin0); + for (int is = 0;is < nspin0;++is) + { + vxcs_op_ao[is] = new hamilt::Veff>(gint, &loc, &lm, kv.kvec_d, potxc, &vxcs_R_ao[is], &vxc_k_ao, &ucell, &gd, pv); + GlobalV::CURRENT_SPIN = is; //caution: Veff::contributeHR depends on GlobalV::CURRENT_SPIN + vxcs_op_ao[is]->contributeHR(); + } +#ifdef __EXX + hamilt::OperatorEXX> vexx_op_ao(&lm, nullptr, &vxc_k_ao, kv); + // ======test======= + // std::vector test_vexxonly_k_ao(pv->nloc); + // hamilt::OperatorEXX> test_vexxonly_op_ao(&lm, nullptr, &test_vexxonly_k_ao, kv); + // ======test======= +#endif + + //3. calculate and write the MO-matrix Exc + Parallel_2D p2d; + set_para2d_MO(*pv, nbands, p2d); + + // ======test======= + // double total_energy = 0.0; + // double exx_energy = 0.0; + // ======test======= + for (int ik = 0;ik < kv.nks;++ik) + { + ModuleBase::GlobalFunc::ZEROS(vxc_k_ao.data(), pv->nloc); + int is = GlobalV::CURRENT_SPIN = kv.isk[ik]; + dynamic_cast*>(vxcs_op_ao[is])->contributeHk(ik); +#ifdef __EXX + if (GlobalC::exx_info.info_global.cal_exx) vexx_op_ao.contributeHk(ik); + // ======test======= + // ModuleBase::GlobalFunc::ZEROS(test_vexxonly_k_ao.data(), pv->nloc); + // if (GlobalC::exx_info.info_global.cal_exx) test_vexxonly_op_ao.contributeHk(ik); + // std::vector test_vexxonly_k_mo = cVc(test_vexxonly_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); + // exx_energy += all_band_energy(ik, test_vexxonly_k_mo, p2d, wg); + // ======test======= +#endif + std::vector vxc_k_mo = cVc(vxc_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); + // write + ModuleIO::save_mat(-1, vxc_k_mo.data(), nbands, + false/*binary*/, true/*triangle*/, false/*append*/, + "Vxc", "k-" + std::to_string(ik), p2d, drank); + // ======test======= + // total_energy += all_band_energy(ik, vxc_k_mo, p2d, wg); + // ======test======= + } + // ======test======= + // total_energy -= 0.5 * exx_energy; + // std::cout << "total energy: " << total_energy << std::endl; + // std::cout << "etxc: " << etxc << std::endl; + // std::cout << "vtxc_cal: " << total_energy - 0.5 * exx_energy << std::endl; + // std::cout << "vtxc_ref: " << vtxc << std::endl; + // std::cout << "exx_energy: " << 0.5 * exx_energy << std::endl; + // ======test======= + delete potxc; + for (int is = 0;is < nspin0;++is) delete vxcs_op_ao[is]; + } +} \ No newline at end of file diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index b01376a011..db7315e7a7 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -225,6 +225,7 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_hs", out_mat_hs, "output H and S matrix"); ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_hs2", out_mat_hs2, "output H(R) and S(R) matrix"); ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_dh", out_mat_dh, "output of derivative of H(R) matrix"); + ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_xc", out_mat_xc, "output exchange-correlation matrix in KS-orbital representation"); ModuleBase::GlobalFunc::OUTP(ofs, "out_interval", out_interval, "interval for printing H(R) and S(R) matrix during MD"); ModuleBase::GlobalFunc::OUTP(ofs, "out_app_flag", out_app_flag, "whether output r(R), H(R), S(R), T(R), and dH(R) matrices in an append manner during MD"); ModuleBase::GlobalFunc::OUTP(ofs, "out_mat_t", out_mat_t, "output T(R) matrix"); diff --git a/tests/integrate/207_NO_KP_OXC/INPUT b/tests/integrate/207_NO_KP_OXC/INPUT new file mode 100644 index 0000000000..a97c26dd21 --- /dev/null +++ b/tests/integrate/207_NO_KP_OXC/INPUT @@ -0,0 +1,30 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation scf + +nbands 6 +symmetry 0 +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB +gamma_only 0 + +#Parameters (2.Iteration) +ecutwfc 20 +scf_thr 1e-8 +scf_nmax 50 + +#Parameters (3.Basis) +basis_type lcao + +#Parameters (4.Smearing) +smearing_method gauss +smearing_sigma 0.002 + +#Parameters (5.Mixing) +mixing_type broyden +mixing_beta 0.7 + +out_mat_xc 1 +ks_solver scalapack_gvx + diff --git a/tests/integrate/207_NO_KP_OXC/KPT b/tests/integrate/207_NO_KP_OXC/KPT new file mode 100644 index 0000000000..f5f7f4ec34 --- /dev/null +++ b/tests/integrate/207_NO_KP_OXC/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/tests/integrate/207_NO_KP_OXC/STRU b/tests/integrate/207_NO_KP_OXC/STRU new file mode 100644 index 0000000000..37d6597727 --- /dev/null +++ b/tests/integrate/207_NO_KP_OXC/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +Si 14 Si_ONCV_PBE-1.0.upf upf201 + +NUMERICAL_ORBITAL +Si_gga_8au_60Ry_2s2p1d.orb + +LATTICE_CONSTANT +10.2 // add lattice constant + +LATTICE_VECTORS +0.5 0.5 0.0 +0.5 0.0 0.5 +0.0 0.5 0.5 + +ATOMIC_POSITIONS +Direct + +Si // Element type +0.0 // magnetism +2 +0.00 0.00 0.00 1 1 1 +0.25 0.25 0.25 1 1 1 diff --git a/tests/integrate/207_NO_KP_OXC/jd b/tests/integrate/207_NO_KP_OXC/jd new file mode 100644 index 0000000000..5f49129742 --- /dev/null +++ b/tests/integrate/207_NO_KP_OXC/jd @@ -0,0 +1 @@ +test for output XC matrix with multipy k-points diff --git a/tests/integrate/207_NO_KP_OXC/k-1-Vxc.ref b/tests/integrate/207_NO_KP_OXC/k-1-Vxc.ref new file mode 100644 index 0000000000..6651d03e42 --- /dev/null +++ b/tests/integrate/207_NO_KP_OXC/k-1-Vxc.ref @@ -0,0 +1,6 @@ +6 (-0.805858,3.20475e-31) (-1.38778e-16,4.41729e-16) (-9.28077e-17,2.16937e-30) (1.03216e-16,1.97215e-31) (7.21645e-16,-6.64683e-16) (-1.04083e-17,-5.91646e-31) + (-0.745176,-2.69938e-30) (-1.38778e-16,3.9443e-31) (2.77556e-17,-2.36658e-30) (-0.100899,-4.89558e-15) (6.93889e-17,7.09975e-30) + (-0.827683,1.12166e-30) (-1.8735e-16,-2.21867e-30) (1.11022e-16,1.97215e-30) (4.996e-16,-1.19047e-16) + (-0.827683,3.57453e-31) (3.1225e-17,-1.72563e-30) (-3.88578e-16,1.19825e-16) + (-0.745849,3.20475e-31) (-5.38225e-17,7.39557e-31) + (-0.739459,3.69779e-31) diff --git a/tests/integrate/207_NO_KP_OXC/result.ref b/tests/integrate/207_NO_KP_OXC/result.ref new file mode 100644 index 0000000000..7d89999bce --- /dev/null +++ b/tests/integrate/207_NO_KP_OXC/result.ref @@ -0,0 +1,4 @@ +etotref -211.4466153062836 +etotperatomref -105.7233076531 +CompareXC_pass 0 +totaltimeref 0.62879 diff --git a/tests/integrate/307_NO_GO_OXC/INPUT b/tests/integrate/307_NO_GO_OXC/INPUT new file mode 100644 index 0000000000..1b7f4fb774 --- /dev/null +++ b/tests/integrate/307_NO_GO_OXC/INPUT @@ -0,0 +1,36 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation scf + +nbands 6 +symmetry 1 +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB +gamma_only 1 + +#Parameters (2.Iteration) +ecutwfc 20 +scf_thr 1e-8 +scf_nmax 100 + + +#Parameters (3.Basis) +basis_type lcao + +#Parameters (4.Smearing) +smearing_method gauss +smearing_sigma 0.002 + +out_mat_xc 1 + +#Parameters (5.Mixing) +mixing_type broyden +mixing_beta 0.7 +mixing_gg0 1.5 + +ks_solver scalapack_gvx + +bx 2 +by 2 +bz 2 \ No newline at end of file diff --git a/tests/integrate/307_NO_GO_OXC/KPT b/tests/integrate/307_NO_GO_OXC/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/integrate/307_NO_GO_OXC/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/integrate/307_NO_GO_OXC/STRU b/tests/integrate/307_NO_GO_OXC/STRU new file mode 100644 index 0000000000..8202fbee08 --- /dev/null +++ b/tests/integrate/307_NO_GO_OXC/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +Si 14 Si_ONCV_PBE-1.0.upf upf201 + +NUMERICAL_ORBITAL +Si_gga_8au_60Ry_2s2p1d.orb + +LATTICE_CONSTANT +20 // add lattice constant + +LATTICE_VECTORS +0.5 0.5 0.0 +0.5 0.0 0.5 +0.0 0.5 0.5 + +ATOMIC_POSITIONS +Direct + +Si // Element type +0.0 // magnetism +2 +0.00 0.00 0.00 1 1 1 +0.25 0.25 0.25 1 1 1 diff --git a/tests/integrate/307_NO_GO_OXC/jd b/tests/integrate/307_NO_GO_OXC/jd new file mode 100644 index 0000000000..13e6c830e3 --- /dev/null +++ b/tests/integrate/307_NO_GO_OXC/jd @@ -0,0 +1 @@ +test for output XC matrix with gamma_only diff --git a/tests/integrate/307_NO_GO_OXC/k-0-Vxc.ref b/tests/integrate/307_NO_GO_OXC/k-0-Vxc.ref new file mode 100644 index 0000000000..3e11a37e42 --- /dev/null +++ b/tests/integrate/307_NO_GO_OXC/k-0-Vxc.ref @@ -0,0 +1,6 @@ +6 -0.684691 -8.83712e-07 1.13515e-16 1.82341e-16 2.58361e-17 3.54064e-17 + -0.71493 -8.00801e-18 2.58871e-17 -5.41766e-17 6.18429e-17 + -0.593224 7.78235e-17 -6.23292e-17 -1.32096e-06 + -0.593224 -8.16693e-17 6.51331e-07 + -0.593224 8.8859e-07 + -0.619565 diff --git a/tests/integrate/307_NO_GO_OXC/result.ref b/tests/integrate/307_NO_GO_OXC/result.ref new file mode 100644 index 0000000000..8b42dbf7ee --- /dev/null +++ b/tests/integrate/307_NO_GO_OXC/result.ref @@ -0,0 +1,7 @@ +etotref -204.1062806122546 +etotperatomref -102.0531403061 +CompareXC_pass 0 +pointgroupref T_d +spacegroupref O_h +nksibzref 1 +totaltimeref 0.52294 diff --git a/tests/integrate/CASES b/tests/integrate/CASES index 27cb8694de..742073013a 100644 --- a/tests/integrate/CASES +++ b/tests/integrate/CASES @@ -150,6 +150,7 @@ 207_NO_KP_OH2 207_NO_KP_OHS_SPIN4 207_NO_KP_OTdH +207_NO_KP_OXC 207_NO_OK 207_NO_KP_OS 208_NO_KP_CF_RE @@ -196,6 +197,7 @@ 304_NO_GO_FM 304_NO_GO_ocp 307_NO_GO_OH +307_NO_GO_OXC 308_NO_GO_CF_RE 308_NO_GO_CS_CR 308_NO_GO_RE_MB diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index 25d18844de..e6936c6b03 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -34,6 +34,7 @@ has_dos=`grep -En '(^|[[:space:]])out_dos($|[[:space:]])' INPUT | awk '{print $2 has_cond=`grep -En '(^|[[:space:]])cal_cond($|[[:space:]])' INPUT | awk '{print $2}'` has_hs=`grep -En '(^|[[:space:]])out_mat_hs($|[[:space:]])' INPUT | awk '{print $2}'` has_hs2=`grep -En '(^|[[:space:]])out_mat_hs2($|[[:space:]])' INPUT | awk '{print $2}'` +has_xc=`grep -En '(^|[[:space:]])out_mat_xc($|[[:space:]])' INPUT | awk '{print $2}'` has_r=`grep -En '(^|[[:space:]])out_mat_r($|[[:space:]])' INPUT | awk '{print $2}'` deepks_out_labels=`grep deepks_out_labels INPUT | awk '{print $2}' | sed s/[[:space:]]//g` deepks_bandgap=`grep deepks_bandgap INPUT | awk '{print $2}' | sed s/[[:space:]]//g` @@ -212,6 +213,18 @@ if ! test -z "$has_hs" && [ $has_hs == 1 ]; then echo "CompareS_pass $?" >>$1 fi +if ! test -z "$has_xc" && [ $has_xc == 1 ]; then + if ! test -z "$gamma_only" && [ $gamma_only == 1 ]; then + xcref=k-0-Vxc.ref + xccal=OUT.autotest/k-0-Vxc + else + xcref=k-1-Vxc.ref + xccal=OUT.autotest/k-1-Vxc + fi + python3 ../tools/CompareFile.py $xcref $xccal 4 + echo "CompareXC_pass $?" >>$1 +fi + #echo $has_hs2 if ! test -z "$has_hs2" && [ $has_hs2 == 1 ]; then #python3 ../tools/CompareFile.py data-HR-sparse_SPIN0.csr.ref OUT.autotest/data-HR-sparse_SPIN0.csr 8 From e605c60c02e5e4ca35166a1dd21b51bba691005d Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Fri, 5 Jan 2024 16:33:34 +0800 Subject: [PATCH 3/3] add DFT+U term --- docs/advanced/input_files/input-main.md | 2 +- source/module_io/write_Vxc.hpp | 10 ++++------ 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 9d76f07bd5..46da3473ae 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -1575,7 +1575,7 @@ These variables are used to control the output of properties. - **Type**: Boolean - **Availability**: Numerical atomic orbital basis -- **Description**: Whether to print the upper triangular part of the exchange-correlation matrices in **Kohn-Sham orbital representation** (unit: Ry): $\braket{\psi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}|\psi_j}$ for each k point into files in the directory `OUT.${suffix}`, which is useful for the subsequent GW calculation . The files are named `k-$k-Vxc`, the meaning of `$k`corresponding to k point and spin is same as [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs). +- **Description**: Whether to print the upper triangular part of the exchange-correlation matrices in **Kohn-Sham orbital representation** (unit: Ry): $\braket{\psi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}+V_\text{DFTU}|\psi_j}$ for each k point into files in the directory `OUT.${suffix}`, which is useful for the subsequent GW calculation. (Note that currently DeePKS term is not included. ) The files are named `k-$k-Vxc`, the meaning of `$k`corresponding to k point and spin is same as [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs). - **Default**: False ### out_app_flag diff --git a/source/module_io/write_Vxc.hpp b/source/module_io/write_Vxc.hpp index 1b351d1a9b..c150825ce1 100644 --- a/source/module_io/write_Vxc.hpp +++ b/source/module_io/write_Vxc.hpp @@ -1,10 +1,6 @@ #include "module_psi/psi.h" -#include "module_hamilt_lcao/module_gint/gint_tools.h" -#include "module_hamilt_lcao/module_gint/gint_gamma.h" -#include "module_hamilt_lcao/module_gint/gint_k.h" -#include "module_elecstate/potentials/pot_xc.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h" -#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h" #include "module_base/scalapack_connector.h" #include "module_base/parallel_reduce.h" template struct TGint; @@ -118,7 +114,7 @@ namespace ModuleIO /// @brief write the Vxc matrix in KS orbital representation, usefull for GW calculation - /// including terms: Vxc(local/semi-local), Vexx + /// including terms: local/semi-local XC, EXX, DFTU template void write_Vxc(int nspin, int nbasis, int drank, const psi::Psi& psi, const UnitCell& ucell, Structure_Factor& sf, const ModulePW::PW_Basis& rho_basis, const ModulePW::PW_Basis& rhod_basis, const ModuleBase::matrix& vloc, @@ -165,6 +161,7 @@ namespace ModuleIO // hamilt::OperatorEXX> test_vexxonly_op_ao(&lm, nullptr, &test_vexxonly_k_ao, kv); // ======test======= #endif + hamilt::OperatorDFTU> vdftu_op_ao(&lm, kv.kvec_d, nullptr, &vxc_k_ao, kv.isk); //4. calculate and write the MO-matrix Exc Parallel_2D p2d; @@ -188,6 +185,7 @@ namespace ModuleIO // exx_energy += all_band_energy(ik, test_vexxonly_k_mo, p2d, wg); // ======test======= #endif + if (GlobalV::dft_plus_u) vdftu_op_ao.contributeHk(ik); std::vector vxc_k_mo = cVc(vxc_k_ao.data(), &psi(ik, 0, 0), nbasis, nbands, *pv, p2d); // write ModuleIO::save_mat(-1, vxc_k_mo.data(), nbands,