From e15e8f532a692d1aa37e49430c0ad9898419ce98 Mon Sep 17 00:00:00 2001 From: pxlxingliang Date: Fri, 28 Jun 2024 10:52:55 +0800 Subject: [PATCH 1/8] refactor: replace the interface of out_dm by using DensityMatrix.DMK instead of loc.DM --- source/Makefile.Objects | 4 +- source/module_elecstate/elecstate_lcao.cpp | 28 -- source/module_esolver/esolver_ks_lcao.cpp | 32 +- source/module_esolver/esolver_ks_lcao.h | 4 - .../module_esolver/esolver_ks_lcao_elec.cpp | 1 - .../module_esolver/esolver_ks_lcao_tddft.cpp | 1 - source/module_io/CMakeLists.txt | 4 +- source/module_io/dm_io.h | 44 --- source/module_io/io_dmk.cpp | 323 ++++++++++++++++++ source/module_io/io_dmk.h | 70 ++++ source/module_io/output_dm.cpp | 55 --- source/module_io/output_dm.h | 43 --- source/module_io/read_dm.cpp | 173 ---------- source/module_io/test_serial/CMakeLists.txt | 4 +- .../{dm_io_test.cpp => io_dmk_test.cpp} | 34 +- source/module_io/write_dm.cpp | 195 ----------- tests/integrate/314_NO_GO_dm_out/SPIN1_DM.ref | 35 ++ tests/integrate/314_NO_GO_dm_out/result.ref | 1 + tests/integrate/tools/catch_properties.sh | 25 +- 19 files changed, 473 insertions(+), 603 deletions(-) delete mode 100644 source/module_io/dm_io.h create mode 100644 source/module_io/io_dmk.cpp create mode 100644 source/module_io/io_dmk.h delete mode 100644 source/module_io/output_dm.cpp delete mode 100644 source/module_io/output_dm.h delete mode 100644 source/module_io/read_dm.cpp rename source/module_io/test_serial/{dm_io_test.cpp => io_dmk_test.cpp} (77%) delete mode 100644 source/module_io/write_dm.cpp create mode 100644 tests/integrate/314_NO_GO_dm_out/SPIN1_DM.ref diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 635b842c46..d3265d2d10 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -479,7 +479,7 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\ nscf_fermi_surf.o\ istate_charge.o\ istate_envelope.o\ - read_dm.o\ + io_dmk.o\ unk_overlap_lcao.o\ read_wfc_nao.o\ read_wfc_lcao.o\ @@ -487,8 +487,6 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\ write_HS_sparse.o\ single_R_io.o\ write_HS_R.o\ - write_dm.o\ - output_dm.o\ write_dmr.o\ sparse_matrix.o\ output_mulliken.o\ diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index 9b2a294dd4..1156f67dea 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -106,24 +106,6 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) // cal_dm(this->loc->ParaV, this->wg, psi, this->loc->dm_gamma); elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM)); this->DM->cal_DMR(); - if (this->loc->out_dm) // keep interface for old Output_DM until new one is ready - { - this->loc->dm_gamma.resize(GlobalV::NSPIN); - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - this->loc->set_dm_gamma(is, this->DM->get_DMK_pointer(is)); - } - } - ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d"); - for (int ik = 0; ik < psi.get_nk(); ++ik) - { - // for gamma_only case, no convertion occured, just for print. - // old 2D-to-Grid conversion has been replaced by new Gint Refactor 2023/09/25 - if (this->loc->out_dm) // keep interface for old Output_DM until new one is ready - { - this->loc->cal_dk_gamma_from_2D_pub(); - } - } } for (int is = 0; is < GlobalV::NSPIN; is++) @@ -191,16 +173,6 @@ void ElecStateLCAO::dmToRho(std::vector pexsi_DM, std::vectorloc->out_dm) // keep interface for old Output_DM until new one is ready - { - for (int is = 0; is < nspin; ++is) - { - this->loc->set_dm_gamma(is, pexsi_DM[is]); - } - this->loc->cal_dk_gamma_from_2D_pub(); - } - this->get_DM()->pexsi_EDM = pexsi_EDM; for (int is = 0; is < nspin; is++) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 0f522ab7a8..fed5d1e1d7 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -46,6 +46,8 @@ #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" #include "module_io/write_dmr.h" #include "module_io/write_wfc_nao.h" +#include "module_io/io_dmk.h" + namespace ModuleESolver { @@ -1000,7 +1002,6 @@ void ESolver_KS_LCAO::iter_finish(int iter) for (int is = 0; is < GlobalV::NSPIN; is++) { this->create_Output_Rho(is, iter, "tmp_").write(); - this->create_Output_DM(is, iter).write(); if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { this->create_Output_Kin(is, iter, "tmp_").write(); @@ -1081,12 +1082,19 @@ void ESolver_KS_LCAO::after_scf(const int istep) } // 4) write density matrix - if (this->LOC.out_dm) + if (INPUT.out_dm) { - for (int is = 0; is < GlobalV::NSPIN; is++) + std::vector efermis(GlobalV::NSPIN == 2 ? 2 : 1); + for (int ispin = 0; ispin < efermis.size(); ispin++) { - this->create_Output_DM(is, istep).write(); + efermis[ispin] = this->pelec->eferm.get_efval(ispin); } + const int precision = 3; + ModuleIO::write_dmk(dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(), + precision, + efermis, + &(GlobalC::ucell), + this->orb_con.ParaV); } // 5) write Vxc @@ -1295,22 +1303,6 @@ bool ESolver_KS_LCAO::do_after_converge(int& iter) //! the 16th function of ESolver_KS_LCAO: create_Output_DM //! mohan add 2024-05-11 //------------------------------------------------------------------------------ -template -ModuleIO::Output_DM ESolver_KS_LCAO::create_Output_DM(int is, int iter) -{ - const int precision = 3; - - return ModuleIO::Output_DM(this->GridT, - is, - iter, - precision, - this->LOC.out_dm, - this->LOC.DM, - this->pelec->eferm.get_efval(is), - &(GlobalC::ucell), - GlobalV::global_out_dir, - GlobalV::GAMMA_ONLY_LOCAL); -} //------------------------------------------------------------------------------ //! the 17th function of ESolver_KS_LCAO: create_Output_DM1 diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 2ecdb67de9..f362d28cd4 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -12,7 +12,6 @@ #include "module_ri/Mix_DMk_2D.h" #endif #include "module_basis/module_nao/two_center_bundle.h" -#include "module_io/output_dm.h" #include "module_io/output_mat_sparse.h" #include @@ -99,9 +98,6 @@ class ESolver_KS_LCAO : public ESolver_KS void beforesolver(const int istep); //---------------------------------------------------------------------- - /// @brief create ModuleIO::Output_DM object to output density matrix - ModuleIO::Output_DM create_Output_DM(int is, int iter); - /// @brief create ModuleIO::Output_Mat_Sparse object to output sparse density matrix of H, S, T, r ModuleIO::Output_Mat_Sparse create_Output_Mat_Sparse(int istep); diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 41426c335c..bd0a510e3c 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -24,7 +24,6 @@ #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" -#include "module_io/dm_io.h" #include "module_io/rho_io.h" #include "module_io/write_pot.h" diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index e05512f81c..cc035fbaf6 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -2,7 +2,6 @@ #include "module_io/cal_r_overlap_R.h" #include "module_io/dipole_io.h" -#include "module_io/dm_io.h" #include "module_io/rho_io.h" #include "module_io/td_current_io.h" #include "module_io/write_HS.h" diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index 47a4e30316..1a5337ff6f 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -53,14 +53,12 @@ if(ENABLE_LCAO) nscf_fermi_surf.cpp istate_charge.cpp istate_envelope.cpp - read_dm.cpp read_wfc_nao.cpp read_wfc_lcao.cpp write_wfc_nao.cpp - write_dm.cpp + io_dmk.cpp write_dmr.cpp dos_nao.cpp - output_dm.cpp sparse_matrix.cpp file_reader.cpp csr_reader.cpp diff --git a/source/module_io/dm_io.h b/source/module_io/dm_io.h deleted file mode 100644 index 346763b2a6..0000000000 --- a/source/module_io/dm_io.h +++ /dev/null @@ -1,44 +0,0 @@ -#ifndef DM_IO_H -#define DM_IO_H - -#include -#include "module_cell/unitcell.h" - -namespace ModuleIO -{ - -void read_dm( -#ifdef __MPI - const int nnrg, - const int* trace_lo, -#endif - const bool gamma_only_local, - const int nlocal, - const int nspin, - const int &is, - const std::string &fn, - double*** DM, - double** DM_R, - double& ef, - const UnitCell* ucell); - -void write_dm( -#ifdef __MPI - const int* trace_lo, -#endif - const int &is, - const int &iter, - const std::string &fn, - int precision, - int out_dm, - double*** DM, - const double& ef, - const UnitCell* ucell, - const int my_rank, - const int nspin, - const int nlocal); - -} - -#endif - diff --git a/source/module_io/io_dmk.cpp b/source/module_io/io_dmk.cpp new file mode 100644 index 0000000000..fc29f86bd6 --- /dev/null +++ b/source/module_io/io_dmk.cpp @@ -0,0 +1,323 @@ +#include "module_base/parallel_common.h" +#include "module_base/timer.h" +#include "module_io/io_dmk.h" +#include "module_base/scalapack_connector.h" + + +void ModuleIO::read_dmk( +#ifdef __MPI + const int nnrg, + const int* trace_lo, +#endif + const bool gamma_only_local, + const int nlocal, + const int nspin, + const int& is, + const std::string& fn, + double*** DM, + double** DM_R, + double& ef, + const UnitCell* ucell) +{ + ModuleBase::TITLE("ModuleIO", "read_dmk"); + ModuleBase::timer::tick("ModuleIO", "read_dmk"); + + GlobalV::ofs_running << "\n processor 0 is reading density matrix from file < " << fn << " > " << std::endl; + // xiaohui modify 2015-03-25 + // bool quit_mesia = false; + bool quit_abacus = false; + + std::ifstream ifs; + if (GlobalV::MY_RANK == 0) + { + ifs.open(fn.c_str()); + if (!ifs) + { + // xiaohui modify 2015-03-25 + // quit_mesia = true; + quit_abacus = true; + } + else + { + // if the number is not match, + // quit the program or not. + bool quit = false; + + std::string name; + ifs >> name; + + // check lattice constant, unit is Angstrom + ModuleBase::CHECK_DOUBLE(ifs, ucell->lat0 * ModuleBase::BOHR_TO_A, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e11, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e12, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e13, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e21, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e22, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e23, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e31, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e32, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e33, quit); + + for (int it = 0; it < ucell->ntype; it++) + { + ModuleBase::CHECK_STRING(ifs, ucell->atoms[it].label, quit); + } + + for (int it = 0; it < ucell->ntype; it++) + { + ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].na, quit); + } + + std::string coordinate; + ifs >> coordinate; + + for (int it = 0; it < ucell->ntype; it++) + { + for (int ia = 0; ia < ucell->atoms[it].na; ia++) + { + ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].x, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].y, quit); + ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].z, quit); + } + } + + ModuleBase::CHECK_INT(ifs, nspin); + ModuleBase::GlobalFunc::READ_VALUE(ifs, ef); + ModuleBase::CHECK_INT(ifs, nlocal); + ModuleBase::CHECK_INT(ifs, nlocal); + } // If file exist, read in data. + } // Finish reading the first part of density matrix. + +#ifndef __MPI + GlobalV::ofs_running << " Read SPIN = " << is + 1 << " density matrix now." << std::endl; + + if (gamma_only_local) + { + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) + { + ifs >> DM[is][i][j]; + } + } + } + else + { +#ifdef __MPI + ModuleBase::WARNING_QUIT("ModuleIO::read_dmk", "The nnrg should not be update"); + ModuleBase::CHECK_INT(ifs, nnrg); + + for (int i = 0; i < nnrg; ++i) + { + ifs >> DM_R[is][i]; + } +#endif + } +#else + + // distribution of necessary data + // xiaohui modify 2015-03-25 + // Parallel_Common::bcast_bool(quit_mesia); + Parallel_Common::bcast_bool(quit_abacus); + // xiaohui modify 2015-03-25 + // if(quit_mesia) + if (quit_abacus) + { + ModuleBase::WARNING_QUIT("ModuleIO::read_dmk", "Can not find the density matrix file."); + } + + Parallel_Common::bcast_double(ef); + + if (gamma_only_local) + { + + double* tmp = new double[nlocal]; + for (int i = 0; i < nlocal; ++i) + { + // GlobalV::ofs_running << " i=" << i << std::endl; + ModuleBase::GlobalFunc::ZEROS(tmp, nlocal); + if (GlobalV::MY_RANK == 0) + { + for (int j = 0; j < nlocal; ++j) + { + ifs >> tmp[j]; + } + } + Parallel_Common::bcast_double(tmp, nlocal); + + const int mu = trace_lo[i]; + if (mu >= 0) + { + for (int j = 0; j < nlocal; ++j) + { + const int nu = trace_lo[j]; + if (nu >= 0) + { + DM[is][mu][nu] = tmp[j]; + } + } + } + } // i + delete[] tmp; + } + else + { + ModuleBase::WARNING_QUIT("ModuleIO::read_dmk", "not ready to readin DM_R"); + } +#endif + if (GlobalV::MY_RANK == 0) + ifs.close(); + + GlobalV::ofs_running << " Finish reading density matrix." << std::endl; + + ModuleBase::timer::tick("ModuleIO", "read_dmk"); + return; +} + +std::string ModuleIO::dmk_gen_fname(const bool gamma_only, const int ispin, const int ik) +{ + if(gamma_only) + { + return "SPIN" + std::to_string(ispin + 1) + "_DM"; + } + else + { + // this case is not implemented now. + ModuleBase::WARNING_QUIT("dmk_gen_fname","Not implemented for non-gamma_only case."); + } +} + +template +void ModuleIO::write_dmk( + const std::vector> dmk, + const int precision, + const std::vector efs, + const UnitCell* ucell, + const Parallel_2D& pv + ) +{ + ModuleBase::TITLE("ModuleIO","write_dmk"); + ModuleBase::timer::tick("ModuleIO","write_dmk"); + + int my_rank = 0; +#ifdef __MPI + MPI_Comm_rank(pv.comm_2D, &my_rank); +#endif + + bool gamma_only = std::is_same::value; + int nlocal = pv.get_global_row_size(); + int nspin = efs.size(); + int nk = dmk.size() / nspin; + if (nk * nspin != dmk.size()) + { + ModuleBase::WARNING_QUIT("write_dmk","The size of dmk is not consistent with nspin and nk."); + } + Parallel_2D pv_glb; + + // when nspin == 2, assume the order of K in dmk is K1_up, K2_up, ..., K1_down, K2_down, ... + for (int ispin = 0; ispin < nspin; ispin++) + { + for (int ik = 0; ik < nk; ik++) + { + // gather dmk[ik] to dmk_global + std::vector dmk_global(my_rank == 0 ? nlocal * nlocal : 0); +#ifdef __MPI + pv_glb.set(nlocal, nlocal, nlocal, pv.comm_2D, pv.blacs_ctxt); + Cpxgemr2d(nlocal, + nlocal, + const_cast(dmk[ik + nk * ispin].data()), + 1, + 1, + const_cast(pv.desc), + dmk_global.data(), + 1, + 1, + pv_glb.desc, + pv_glb.blacs_ctxt); +#else + dmk_global = dmk[ik + nk * ispin]; +#endif + + if (my_rank == 0) + { + std::string fn = GlobalV::global_out_dir + dmk_gen_fname(gamma_only, ispin, ik); + std::ofstream ofs(fn.c_str()); + + if (!ofs) + { + ModuleBase::WARNING("ModuleIO::write_dmk", "Can't create DENSITY MATRIX File!"); + } + + // write the UnitCell information + ofs << ucell->latName << std::endl; + ofs << " " << ucell->lat0 * ModuleBase::BOHR_TO_A << std::endl; + ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " << ucell->latvec.e13 << std::endl; + ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " << ucell->latvec.e23 << std::endl; + ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " << ucell->latvec.e33 << std::endl; + for(int it=0; itntype; it++) + { + ofs << " " << ucell->atoms[it].label; + } + ofs << std::endl; + for(int it=0; itntype; it++) + { + ofs << " " << ucell->atoms[it].na; + } + ofs << std::endl; + ofs << "Direct" << std::endl; + for(int it=0; itntype; it++) + { + Atom* atom = &ucell->atoms[it]; + ofs << std::setprecision(15); + for(int ia=0; iaatoms[it].na; ia++) + { + ofs << " " << atom->taud[ia].x + << " " << atom->taud[ia].y + << " " << atom->taud[ia].z << std::endl; + } + } + ofs << "\n " << dmk.size(); // nspin + ofs << "\n " << efs[ispin] << " (fermi energy)"; + ofs << "\n " << nlocal << " " << nlocal << std::endl; + + ofs << std::setprecision(precision); + ofs << std::scientific; + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) + { + if (j % 8 == 0) + { + ofs << "\n"; + } + if (std::is_same::value) + { + ofs << " " << dmk_global[i * nlocal + j]; + } + else if (std::is_same, T>::value) + { + ofs << " (" << std::real(dmk_global[i * nlocal + j]) << "," << std::imag(dmk_global[i * nlocal + j]) << ")"; + } + } + } + ofs.close(); + } // rank0 + } // ik + } // ispin + + ModuleBase::timer::tick("ModuleIO","write_dmk"); +} + +template void ModuleIO::write_dmk( + const std::vector> dmk, + const int precision, + const std::vector efs, + const UnitCell* ucell, + const Parallel_2D& pv); + +template void ModuleIO::write_dmk>( + const std::vector>> dmk, + const int precision, + const std::vector efs, + const UnitCell* ucell, + const Parallel_2D& pv); \ No newline at end of file diff --git a/source/module_io/io_dmk.h b/source/module_io/io_dmk.h new file mode 100644 index 0000000000..8ac084cf6f --- /dev/null +++ b/source/module_io/io_dmk.h @@ -0,0 +1,70 @@ +#ifndef DM_IO_H +#define DM_IO_H + +#include +#include "module_cell/unitcell.h" +#include "module_basis/module_ao/parallel_2d.h" + +namespace ModuleIO +{ +/** + * @brief Reads the DMK file and stores the data in the provided arrays. + * + * @param nnrg The number of energy points. + * @param trace_lo An array containing the lower trace indices. + * @param gamma_only_local A boolean indicating whether to read only the gamma-only part of the DMK file. + * @param nlocal The number of local orbitals. + * @param nspin The number of spin components. + * @param is The index of the spin component. + * @param fn The filename of the DMK file. + * @param DM A 3D array to store the DMK data. + * @param DM_R A 2D array to store the real part of the DMK data. + * @param ef The Fermi energy. + * @param ucell A pointer to the UnitCell object. + */ +void read_dmk( +#ifdef __MPI + const int nnrg, + const int* trace_lo, +#endif + const bool gamma_only_local, + const int nlocal, + const int nspin, + const int &is, + const std::string &fn, + double*** DM, + double** DM_R, + double& ef, + const UnitCell* ucell); + +/** + * @brief Generates the filename for the DMK file based on the given parameters. + * + * @param gamma_only A boolean indicating whether to generate the filename for the gamma-only part. + * @param ispin The index of the spin component. + * @param ik The index of the k-point. + * @return The generated filename. + */ +std::string dmk_gen_fname(const bool gamma_only, const int ispin, const int ik); + +/** + * @brief Writes the DMK data to a file. + * + * @tparam T The type of the DMK data. + * @param dmk A vector containing the DMK data. The first dimension is nspin*nk, and the second dimension is nlocal*nlocal. DMK is parallel in 2d-block type if using MPI. + * @param precision The precision of the output of DMK. + * @param efs A vector containing the Fermi energies, and should have the same size as the number of SPIN. + * @param ucell A pointer to the UnitCell object. + * @param pv The Parallel_2D object. The 2d-block parallel information of DMK. + */ +template +void write_dmk( + const std::vector> dmk, + const int precision, + const std::vector efs, + const UnitCell* ucell, + const Parallel_2D& pv + ); +} +#endif + diff --git a/source/module_io/output_dm.cpp b/source/module_io/output_dm.cpp deleted file mode 100644 index 2b08c85a7a..0000000000 --- a/source/module_io/output_dm.cpp +++ /dev/null @@ -1,55 +0,0 @@ -#include "module_io/output_dm.h" -#include "module_io/dm_io.h" - -namespace ModuleIO -{ - -Output_DM::Output_DM(const Grid_Technique& GridT, - int is, - int iter, - int precision, - int out_dm, - double*** DM, - const double& ef, - const UnitCell* ucell, - const std::string& directory, - bool gamma_only_local) - : _GridT(GridT), - _is(is), - _iter(iter), - _precision(precision), - _out_dm(out_dm), - _DM(DM), - _ef(ef), - _ucell(ucell), - _directory(directory), - _gamma_only_local(gamma_only_local) -{ - if (gamma_only_local) - { - this->_fn = this->_directory + "/SPIN" + std::to_string(this->_is + 1) + "_DM"; - } - else - { - this->_fn = this->_directory + "/SPIN" + std::to_string(this->_is + 1) + "_DM_R"; - } -} -void Output_DM::write() -{ - ModuleIO::write_dm( -#ifdef __MPI - _GridT.trace_lo.data(), -#endif - _is, - _iter, - _fn, - _precision, - _out_dm, - _DM, - _ef, - _ucell, - GlobalV::MY_RANK, - GlobalV::NSPIN, - GlobalV::NLOCAL); -} -} // namespace ModuleIO diff --git a/source/module_io/output_dm.h b/source/module_io/output_dm.h deleted file mode 100644 index b1415f4c75..0000000000 --- a/source/module_io/output_dm.h +++ /dev/null @@ -1,43 +0,0 @@ -#ifndef OUTPUT_DM_H -#define OUTPUT_DM_H - -#include "module_cell/unitcell.h" -#include "module_hamilt_lcao/module_gint/grid_technique.h" - -#include - -namespace ModuleIO -{ - -/// @brief the output interface to write the density matrix -class Output_DM -{ - public: - Output_DM(const Grid_Technique& GridT, - int is, - int iter, - int precision, - int out_dm, - double*** DM, - const double& ef, - const UnitCell* ucell, - const std::string& directory, - bool gamma_only_local); - void write(); - - private: - const Grid_Technique& _GridT; - int _is; - int _iter; - std::string _fn; - int _precision; - int _out_dm; - double*** _DM; - const double& _ef; - const UnitCell* _ucell; - const std::string& _directory; - bool _gamma_only_local; -}; -} // namespace ModuleIO - -#endif // OUTPUT_DM_H \ No newline at end of file diff --git a/source/module_io/read_dm.cpp b/source/module_io/read_dm.cpp deleted file mode 100644 index 2332866cb3..0000000000 --- a/source/module_io/read_dm.cpp +++ /dev/null @@ -1,173 +0,0 @@ -#include "module_base/parallel_common.h" -#include "module_base/timer.h" -#include "module_io/dm_io.h" - -void ModuleIO::read_dm( -#ifdef __MPI - const int nnrg, - const int* trace_lo, -#endif - const bool gamma_only_local, - const int nlocal, - const int nspin, - const int& is, - const std::string& fn, - double*** DM, - double** DM_R, - double& ef, - const UnitCell* ucell) -{ - ModuleBase::TITLE("ModuleIO", "read_dm"); - ModuleBase::timer::tick("ModuleIO", "read_dm"); - - GlobalV::ofs_running << "\n processor 0 is reading density matrix from file < " << fn << " > " << std::endl; - // xiaohui modify 2015-03-25 - // bool quit_mesia = false; - bool quit_abacus = false; - - std::ifstream ifs; - if (GlobalV::MY_RANK == 0) - { - ifs.open(fn.c_str()); - if (!ifs) - { - // xiaohui modify 2015-03-25 - // quit_mesia = true; - quit_abacus = true; - } - else - { - // if the number is not match, - // quit the program or not. - bool quit = false; - - std::string name; - ifs >> name; - - // check lattice constant, unit is Angstrom - ModuleBase::CHECK_DOUBLE(ifs, ucell->lat0 * ModuleBase::BOHR_TO_A, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e11, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e12, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e13, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e21, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e22, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e23, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e31, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e32, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e33, quit); - - for (int it = 0; it < ucell->ntype; it++) - { - ModuleBase::CHECK_STRING(ifs, ucell->atoms[it].label, quit); - } - - for (int it = 0; it < ucell->ntype; it++) - { - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].na, quit); - } - - std::string coordinate; - ifs >> coordinate; - - for (int it = 0; it < ucell->ntype; it++) - { - for (int ia = 0; ia < ucell->atoms[it].na; ia++) - { - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].x, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].y, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].z, quit); - } - } - - ModuleBase::CHECK_INT(ifs, nspin); - ModuleBase::GlobalFunc::READ_VALUE(ifs, ef); - ModuleBase::CHECK_INT(ifs, nlocal); - ModuleBase::CHECK_INT(ifs, nlocal); - } // If file exist, read in data. - } // Finish reading the first part of density matrix. - -#ifndef __MPI - GlobalV::ofs_running << " Read SPIN = " << is + 1 << " density matrix now." << std::endl; - - if (gamma_only_local) - { - for (int i = 0; i < nlocal; ++i) - { - for (int j = 0; j < nlocal; ++j) - { - ifs >> DM[is][i][j]; - } - } - } - else - { -#ifdef __MPI - ModuleBase::WARNING_QUIT("ModuleIO::read_dm", "The nnrg should not be update"); - ModuleBase::CHECK_INT(ifs, nnrg); - - for (int i = 0; i < nnrg; ++i) - { - ifs >> DM_R[is][i]; - } -#endif - } -#else - - // distribution of necessary data - // xiaohui modify 2015-03-25 - // Parallel_Common::bcast_bool(quit_mesia); - Parallel_Common::bcast_bool(quit_abacus); - // xiaohui modify 2015-03-25 - // if(quit_mesia) - if (quit_abacus) - { - ModuleBase::WARNING_QUIT("ModuleIO::read_dm", "Can not find the density matrix file."); - } - - Parallel_Common::bcast_double(ef); - - if (gamma_only_local) - { - - double* tmp = new double[nlocal]; - for (int i = 0; i < nlocal; ++i) - { - // GlobalV::ofs_running << " i=" << i << std::endl; - ModuleBase::GlobalFunc::ZEROS(tmp, nlocal); - if (GlobalV::MY_RANK == 0) - { - for (int j = 0; j < nlocal; ++j) - { - ifs >> tmp[j]; - } - } - Parallel_Common::bcast_double(tmp, nlocal); - - const int mu = trace_lo[i]; - if (mu >= 0) - { - for (int j = 0; j < nlocal; ++j) - { - const int nu = trace_lo[j]; - if (nu >= 0) - { - DM[is][mu][nu] = tmp[j]; - } - } - } - } // i - delete[] tmp; - } - else - { - ModuleBase::WARNING_QUIT("ModuleIO::read_dm", "not ready to readin DM_R"); - } -#endif - if (GlobalV::MY_RANK == 0) - ifs.close(); - - GlobalV::ofs_running << " Finish reading density matrix." << std::endl; - - ModuleBase::timer::tick("ModuleIO", "read_dm"); - return; -} diff --git a/source/module_io/test_serial/CMakeLists.txt b/source/module_io/test_serial/CMakeLists.txt index 1e6aebe635..8a643d6308 100644 --- a/source/module_io/test_serial/CMakeLists.txt +++ b/source/module_io/test_serial/CMakeLists.txt @@ -12,9 +12,9 @@ AddTest( ) AddTest( - TARGET io_dm_io + TARGET io_dmk_io LIBS ${math_libs} base device cell_info - SOURCES dm_io_test.cpp ../read_dm.cpp ../write_dm.cpp ../output.cpp + SOURCES io_dmk_test.cpp ../io_dmk.cpp ../output.cpp ../../module_basis/module_ao/parallel_2d.cpp ) AddTest( diff --git a/source/module_io/test_serial/dm_io_test.cpp b/source/module_io/test_serial/io_dmk_test.cpp similarity index 77% rename from source/module_io/test_serial/dm_io_test.cpp rename to source/module_io/test_serial/io_dmk_test.cpp index 63ddf38a53..8fda1123e2 100644 --- a/source/module_io/test_serial/dm_io_test.cpp +++ b/source/module_io/test_serial/io_dmk_test.cpp @@ -1,6 +1,6 @@ #include "gtest/gtest.h" #include "gmock/gmock.h" -#include "module_io/dm_io.h" +#include "module_io/io_dmk.h" #include "module_base/global_variable.h" #include "prepare_unitcell.h" @@ -22,16 +22,16 @@ Magnetism::~Magnetism() } /************************************************ - * unit test of read_dm and write_dm + * unit test of read_dmk and write_dmk ***********************************************/ /** * - Tested Functions: - * - read_dm() - * - the function to read density matrix from file + * - read_dmk() + * - the function to read density matrix K from file * - the serial version without MPI - * - write_dm() - * - the function to write density matrix to file + * - write_dmk() + * - the function to write density matrix K to file * - the serial version without MPI */ @@ -86,7 +86,7 @@ TEST_F(DMIOTest,Read) double ef; UcellTestPrepare utp = UcellTestLib["Si"]; ucell = utp.SetUcellInfo(); - ModuleIO::read_dm(GlobalV::GAMMA_ONLY_LOCAL, GlobalV::NLOCAL, + ModuleIO::read_dmk(GlobalV::GAMMA_ONLY_LOCAL, GlobalV::NLOCAL, GlobalV::NSPIN,is,fn,DM,DM_R,ef,ucell); EXPECT_DOUBLE_EQ(ef,0.570336288802337); EXPECT_NEAR(DM[0][0][0],3.904e-01,1e-6); @@ -104,16 +104,28 @@ TEST_F(DMIOTest,Write) double ef; UcellTestPrepare utp = UcellTestLib["Si"]; ucell = utp.SetUcellInfo(); - ModuleIO::read_dm(GlobalV::GAMMA_ONLY_LOCAL, GlobalV::NLOCAL, + ModuleIO::read_dmk(GlobalV::GAMMA_ONLY_LOCAL, GlobalV::NLOCAL, GlobalV::NSPIN,is,fn,DM,DM_R,ef,ucell); EXPECT_DOUBLE_EQ(ef,0.570336288802337); EXPECT_NEAR(DM[0][0][0],3.904e-01,1e-6); EXPECT_NEAR(DM[0][25][25],3.445e-02,1e-6); //then write - std::string ssd = "SPIN1_DM"; int precision = 3; - int out_dm = 1; - ModuleIO::write_dm(is,0,ssd,precision,out_dm,DM,ef,ucell,GlobalV::MY_RANK,GlobalV::NSPIN,GlobalV::NLOCAL); + std::vector> dmk(nspin,std::vector(lgd*lgd,0.0)); + for(int is=0; is(nspin,ef),ucell,pv); std::ifstream ifs; ifs.open("SPIN1_DM"); std::string str((std::istreambuf_iterator(ifs)),std::istreambuf_iterator()); diff --git a/source/module_io/write_dm.cpp b/source/module_io/write_dm.cpp deleted file mode 100644 index ebfd7532b6..0000000000 --- a/source/module_io/write_dm.cpp +++ /dev/null @@ -1,195 +0,0 @@ -#include "module_base/parallel_reduce.h" -#include "module_base/timer.h" -#include "module_io/dm_io.h" - -//------------------------------------------------- -// NOTE for ModuleIO::write_dm -// I will give an example here, suppose we have a 4*4 -// density matrix (symmetry) which is -// 1.1 2.3 3.6 4.2 -// 2.3 5.2 7.1 8.9 -// 3.6 7.1 1.3 3.2 -// 4.2 8.9 3.2 2.4 -// we use two processors, each one has 3 orbitals -// processor 1 has orbital index 1, 2, 4: -// ('no' means no value on this processor) -// 1.1 2.3 no 4.2 -// 2.3 5.2 no 8.9 -// no no no no -// 4.2 8.9 no 2.4 -// processor 2 has orbital index 1, 3, 4; -// 1.1 no 3.6 4.2 -// no no no no -// 3.6 no 1.3 3.2 -// 4.2 no 3.2 2.4 -// now we want to reduce them and print out, -// we plan to reduce them one row by one row, -// then for the first row, we need to set the -// temparary array to 4 (GlobalV::NLOCAL in code), -// then we reduce first row, it will become -// 2.2 2.3 3.6 8.4, -// we notice that the first element and fourth -// element is doubled, that's because the density -// may have overlap, so we need to first count -// for each element, how many times it is duplicated -// on other processors, this is why there is -// a 'count' integer array in the code. -// UPDATED BY MOHAN 2014-05-18 -void ModuleIO::write_dm( -#ifdef __MPI - const int* trace_lo, -#endif - const int &is, - const int &iter, - const std::string &fn, - int precision, - int out_dm, - double*** DM, - const double& ef, - const UnitCell* ucell, - const int my_rank, - const int nspin, - const int nlocal) -{ - ModuleBase::TITLE("ModuleIO","write_dm"); - - if (out_dm==0) - { - return; - } - - ModuleBase::timer::tick("ModuleIO","write_dm"); - - time_t start, end; - std::ofstream ofs; - - if(my_rank==0) - { - start = time(NULL); - - ofs.open(fn.c_str()); - if (!ofs) - { - ModuleBase::WARNING("ModuleIO::write_dm","Can't create DENSITY MATRIX File!"); - } - - //GlobalV::ofs_running << "\n Output charge file." << std::endl; - - ofs << ucell->latName << std::endl;//1 - ofs << " " << ucell->lat0 * ModuleBase::BOHR_TO_A << std::endl; - ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " << ucell->latvec.e13 << std::endl; - ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " << ucell->latvec.e23 << std::endl; - ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " << ucell->latvec.e33 << std::endl; - for(int it=0; itntype; it++) - { - ofs << " " << ucell->atoms[it].label; - } - ofs << std::endl; - for(int it=0; itntype; it++) - { - ofs << " " << ucell->atoms[it].na; - } - ofs << std::endl; - ofs << "Direct" << std::endl; - - for(int it=0; itntype; it++) - { - Atom* atom = &ucell->atoms[it]; - ofs << std::setprecision(15); - for(int ia=0; iaatoms[it].na; ia++) - { - ofs << " " << atom->taud[ia].x - << " " << atom->taud[ia].y - << " " << atom->taud[ia].z << std::endl; - } - } - - ofs << "\n " << nspin; - ofs << "\n " << ef << " (fermi energy)"; - - ofs << "\n " << nlocal << " " << nlocal << std::endl; - - ofs << std::setprecision(precision); - ofs << std::scientific; - - } - - //ofs << "\n " << GlobalV::GAMMA_ONLY_LOCAL << " (GAMMA ONLY LOCAL)" << std::endl; -#ifndef __MPI - - for(int i=0; i tmp(nlocal,0); - std::vector count(nlocal,0); - for (int i=0; i= 0) - { - for (int j=0; j= 0) - { - count[j]=1; - } - } - } - Parallel_Reduce::reduce_all(count.data(), nlocal); - - // reduce the density matrix for 'i' line. - tmp.assign(tmp.size(),0); - if (mu >= 0) - { - for (int j=0; j=0) - { - tmp[j] = DM[is][mu][nu]; - //GlobalV::ofs_running << " dmi=" << i << " j=" << j << " " << DM[is][mu][nu] << std::endl; - } - } - } - Parallel_Reduce::reduce_all(tmp.data(), nlocal); - - if(my_rank==0) - { - for (int j=0; j0) - { - ofs << " " << tmp[j]/(double)count[j]; - } - else - { - ofs << " 0"; - } - } - } - } - std::vector().swap(tmp); - std::vector().swap(count); -#endif - if(my_rank==0) - { - end = time(NULL); - ModuleBase::GlobalFunc::OUT_TIME("write_dm",start,end); - ofs.close(); - } - ModuleBase::timer::tick("ModuleIO","write_dm"); - - return; -} diff --git a/tests/integrate/314_NO_GO_dm_out/SPIN1_DM.ref b/tests/integrate/314_NO_GO_dm_out/SPIN1_DM.ref new file mode 100644 index 0000000000..cdc3ba9105 --- /dev/null +++ b/tests/integrate/314_NO_GO_dm_out/SPIN1_DM.ref @@ -0,0 +1,35 @@ +sc + 5.29177 + 1 0 0 + 0 1 0 + 0 0 1 + H + 2 +Direct + 0 0 0.933859999999186 + 0 0 0.0661400000008143 + + 1 + -0.0883978533958687 (fermi energy) + 10 10 + + 5.773e-01 3.902e-02 1.661e-02 4.797e-17 -2.255e-17 5.773e-01 3.902e-02 -1.661e-02 + -1.461e-17 -4.414e-17 + 3.902e-02 2.637e-03 1.122e-03 3.242e-18 -1.524e-18 3.902e-02 2.637e-03 -1.122e-03 + -9.878e-19 -2.984e-18 + 1.661e-02 1.122e-03 4.777e-04 1.380e-18 -6.487e-19 1.661e-02 1.122e-03 -4.777e-04 + -4.204e-19 -1.270e-18 + 4.797e-17 3.242e-18 1.380e-18 3.986e-33 -1.874e-33 4.797e-17 3.242e-18 -1.380e-18 + -1.214e-33 -3.668e-33 + -2.255e-17 -1.524e-18 -6.487e-19 -1.874e-33 8.809e-34 -2.255e-17 -1.524e-18 6.487e-19 + 5.709e-34 1.724e-33 + 5.773e-01 3.902e-02 1.661e-02 4.797e-17 -2.255e-17 5.773e-01 3.902e-02 -1.661e-02 + -1.461e-17 -4.414e-17 + 3.902e-02 2.637e-03 1.122e-03 3.242e-18 -1.524e-18 3.902e-02 2.637e-03 -1.122e-03 + -9.878e-19 -2.984e-18 + -1.661e-02 -1.122e-03 -4.777e-04 -1.380e-18 6.487e-19 -1.661e-02 -1.122e-03 4.777e-04 + 4.204e-19 1.270e-18 + -1.461e-17 -9.878e-19 -4.204e-19 -1.214e-33 5.709e-34 -1.461e-17 -9.878e-19 4.204e-19 + 3.700e-34 1.118e-33 + -4.414e-17 -2.984e-18 -1.270e-18 -3.668e-33 1.724e-33 -4.414e-17 -2.984e-18 1.270e-18 + 1.118e-33 3.376e-33 \ No newline at end of file diff --git a/tests/integrate/314_NO_GO_dm_out/result.ref b/tests/integrate/314_NO_GO_dm_out/result.ref index 16627e485b..670d9daff7 100644 --- a/tests/integrate/314_NO_GO_dm_out/result.ref +++ b/tests/integrate/314_NO_GO_dm_out/result.ref @@ -1,3 +1,4 @@ etotref -31.58774083965451 etotperatomref -15.7938704198 +DM_different 0 totaltimeref 0.18486 diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index 8f5aa446e7..4fdf6c16ac 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -42,7 +42,7 @@ has_lowf=`grep out_wfc_lcao INPUT | awk '{print $2}' | sed s/[[:space:]]//g` out_app_flag=`grep out_app_flag INPUT | awk '{print $2}' | sed s/[[:space:]]//g` has_wfc_r=`grep out_wfc_r INPUT | awk '{print $2}' | sed s/[[:space:]]//g` has_wfc_pw=`grep out_wfc_pw INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_dm=`grep "out_dm " INPUT | awk '{print $2}' | sed s/[[:space:]]//g` +out_dm=`grep "out_dm" INPUT | awk '{print $2}' | sed s/[[:space:]]//g` out_mul=`grep out_mul INPUT | awk '{print $2}' | sed s/[[:space:]]//g` gamma_only=`grep gamma_only INPUT | awk '{print $2}' | sed s/[[:space:]]//g` imp_sol=`grep imp_sol INPUT | awk '{print $2}' | sed s/[[:space:]]//g` @@ -341,29 +341,14 @@ if ! test -z "$has_lowf" && [ $has_lowf == 1 ]; then fi if ! test -z "$out_dm" && [ $out_dm == 1 ]; then - dmfile=`ls OUT.autotest/ | grep "^SPIN1_DM"` + dmfile=OUT.autotest/SPIN1_DM + dmref=SPIN1_DM.ref if test -z "$dmfile"; then echo "Can't find DM files" exit 1 else - for dm in $dmfile; - do - if ! test -f OUT.autotest/$dm; then - echo "Irregular DM file found" - exit 1 - else - total_dm=$(awk 'BEGIN {sum=0.0;startline=999} - { - if(NR==7){startline=$1+14;} - else if(NR>=startline) - { - for(i=1;i<=NF;i++){sum+=sqrt($i*$i)} - } - } - END {printf"%.6f",sum}' OUT.autotest/$dm) - echo "$dm $total_dm" >>$1 - fi - done + python3 ../tools/CompareFile.py $dmref $dm 5 + echo "DM_different $?" >>$1 fi fi From a5ff3d25674abed1016afe220158b03c616eb9c6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Fri, 28 Jun 2024 03:30:56 +0000 Subject: [PATCH 2/8] [pre-commit.ci lite] apply automatic fixes --- source/module_esolver/esolver_ks_lcao.cpp | 3 +- source/module_io/io_dmk.cpp | 114 +++++----- source/module_io/io_dmk.h | 51 +++-- source/module_io/test_serial/io_dmk_test.cpp | 206 ++++++++++--------- 4 files changed, 188 insertions(+), 186 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index fed5d1e1d7..bc886ebad9 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -44,10 +44,9 @@ //--------------------------------------------------- #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" +#include "module_io/io_dmk.h" #include "module_io/write_dmr.h" #include "module_io/write_wfc_nao.h" -#include "module_io/io_dmk.h" - namespace ModuleESolver { diff --git a/source/module_io/io_dmk.cpp b/source/module_io/io_dmk.cpp index fc29f86bd6..5fe919b3d2 100644 --- a/source/module_io/io_dmk.cpp +++ b/source/module_io/io_dmk.cpp @@ -1,8 +1,8 @@ -#include "module_base/parallel_common.h" -#include "module_base/timer.h" #include "module_io/io_dmk.h" -#include "module_base/scalapack_connector.h" +#include "module_base/parallel_common.h" +#include "module_base/scalapack_connector.h" +#include "module_base/timer.h" void ModuleIO::read_dmk( #ifdef __MPI @@ -176,28 +176,26 @@ void ModuleIO::read_dmk( std::string ModuleIO::dmk_gen_fname(const bool gamma_only, const int ispin, const int ik) { - if(gamma_only) + if (gamma_only) { return "SPIN" + std::to_string(ispin + 1) + "_DM"; } else { // this case is not implemented now. - ModuleBase::WARNING_QUIT("dmk_gen_fname","Not implemented for non-gamma_only case."); + ModuleBase::WARNING_QUIT("dmk_gen_fname", "Not implemented for non-gamma_only case."); } } template -void ModuleIO::write_dmk( - const std::vector> dmk, - const int precision, - const std::vector efs, - const UnitCell* ucell, - const Parallel_2D& pv - ) +void ModuleIO::write_dmk(const std::vector> dmk, + const int precision, + const std::vector efs, + const UnitCell* ucell, + const Parallel_2D& pv) { - ModuleBase::TITLE("ModuleIO","write_dmk"); - ModuleBase::timer::tick("ModuleIO","write_dmk"); + ModuleBase::TITLE("ModuleIO", "write_dmk"); + ModuleBase::timer::tick("ModuleIO", "write_dmk"); int my_rank = 0; #ifdef __MPI @@ -210,10 +208,10 @@ void ModuleIO::write_dmk( int nk = dmk.size() / nspin; if (nk * nspin != dmk.size()) { - ModuleBase::WARNING_QUIT("write_dmk","The size of dmk is not consistent with nspin and nk."); + ModuleBase::WARNING_QUIT("write_dmk", "The size of dmk is not consistent with nspin and nk."); } Parallel_2D pv_glb; - + // when nspin == 2, assume the order of K in dmk is K1_up, K2_up, ..., K1_down, K2_down, ... for (int ispin = 0; ispin < nspin; ispin++) { @@ -250,32 +248,31 @@ void ModuleIO::write_dmk( // write the UnitCell information ofs << ucell->latName << std::endl; - ofs << " " << ucell->lat0 * ModuleBase::BOHR_TO_A << std::endl; - ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " << ucell->latvec.e13 << std::endl; - ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " << ucell->latvec.e23 << std::endl; - ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " << ucell->latvec.e33 << std::endl; - for(int it=0; itntype; it++) - { - ofs << " " << ucell->atoms[it].label; - } - ofs << std::endl; - for(int it=0; itntype; it++) - { - ofs << " " << ucell->atoms[it].na; - } - ofs << std::endl; - ofs << "Direct" << std::endl; - for(int it=0; itntype; it++) - { - Atom* atom = &ucell->atoms[it]; - ofs << std::setprecision(15); - for(int ia=0; iaatoms[it].na; ia++) - { - ofs << " " << atom->taud[ia].x - << " " << atom->taud[ia].y - << " " << atom->taud[ia].z << std::endl; - } - } + ofs << " " << ucell->lat0 * ModuleBase::BOHR_TO_A << std::endl; + ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " << ucell->latvec.e13 << std::endl; + ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " << ucell->latvec.e23 << std::endl; + ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " << ucell->latvec.e33 << std::endl; + for (int it = 0; it < ucell->ntype; it++) + { + ofs << " " << ucell->atoms[it].label; + } + ofs << std::endl; + for (int it = 0; it < ucell->ntype; it++) + { + ofs << " " << ucell->atoms[it].na; + } + ofs << std::endl; + ofs << "Direct" << std::endl; + for (int it = 0; it < ucell->ntype; it++) + { + Atom* atom = &ucell->atoms[it]; + ofs << std::setprecision(15); + for (int ia = 0; ia < ucell->atoms[it].na; ia++) + { + ofs << " " << atom->taud[ia].x << " " << atom->taud[ia].y << " " << atom->taud[ia].z + << std::endl; + } + } ofs << "\n " << dmk.size(); // nspin ofs << "\n " << efs[ispin] << " (fermi energy)"; ofs << "\n " << nlocal << " " << nlocal << std::endl; @@ -296,28 +293,27 @@ void ModuleIO::write_dmk( } else if (std::is_same, T>::value) { - ofs << " (" << std::real(dmk_global[i * nlocal + j]) << "," << std::imag(dmk_global[i * nlocal + j]) << ")"; + ofs << " (" << std::real(dmk_global[i * nlocal + j]) << "," + << std::imag(dmk_global[i * nlocal + j]) << ")"; } } } ofs.close(); } // rank0 - } // ik - } // ispin + } // ik + } // ispin - ModuleBase::timer::tick("ModuleIO","write_dmk"); + ModuleBase::timer::tick("ModuleIO", "write_dmk"); } -template void ModuleIO::write_dmk( - const std::vector> dmk, - const int precision, - const std::vector efs, - const UnitCell* ucell, - const Parallel_2D& pv); - -template void ModuleIO::write_dmk>( - const std::vector>> dmk, - const int precision, - const std::vector efs, - const UnitCell* ucell, - const Parallel_2D& pv); \ No newline at end of file +template void ModuleIO::write_dmk(const std::vector> dmk, + const int precision, + const std::vector efs, + const UnitCell* ucell, + const Parallel_2D& pv); + +template void ModuleIO::write_dmk>(const std::vector>> dmk, + const int precision, + const std::vector efs, + const UnitCell* ucell, + const Parallel_2D& pv); \ No newline at end of file diff --git a/source/module_io/io_dmk.h b/source/module_io/io_dmk.h index 8ac084cf6f..91afd2ee07 100644 --- a/source/module_io/io_dmk.h +++ b/source/module_io/io_dmk.h @@ -1,15 +1,16 @@ #ifndef DM_IO_H #define DM_IO_H -#include -#include "module_cell/unitcell.h" #include "module_basis/module_ao/parallel_2d.h" +#include "module_cell/unitcell.h" + +#include namespace ModuleIO { /** * @brief Reads the DMK file and stores the data in the provided arrays. - * + * * @param nnrg The number of energy points. * @param trace_lo An array containing the lower trace indices. * @param gamma_only_local A boolean indicating whether to read only the gamma-only part of the DMK file. @@ -24,22 +25,22 @@ namespace ModuleIO */ void read_dmk( #ifdef __MPI - const int nnrg, - const int* trace_lo, + const int nnrg, + const int* trace_lo, #endif - const bool gamma_only_local, - const int nlocal, - const int nspin, - const int &is, - const std::string &fn, - double*** DM, - double** DM_R, - double& ef, - const UnitCell* ucell); + const bool gamma_only_local, + const int nlocal, + const int nspin, + const int& is, + const std::string& fn, + double*** DM, + double** DM_R, + double& ef, + const UnitCell* ucell); /** * @brief Generates the filename for the DMK file based on the given parameters. - * + * * @param gamma_only A boolean indicating whether to generate the filename for the gamma-only part. * @param ispin The index of the spin component. * @param ik The index of the k-point. @@ -49,22 +50,20 @@ std::string dmk_gen_fname(const bool gamma_only, const int ispin, const int ik); /** * @brief Writes the DMK data to a file. - * + * * @tparam T The type of the DMK data. - * @param dmk A vector containing the DMK data. The first dimension is nspin*nk, and the second dimension is nlocal*nlocal. DMK is parallel in 2d-block type if using MPI. + * @param dmk A vector containing the DMK data. The first dimension is nspin*nk, and the second dimension is + * nlocal*nlocal. DMK is parallel in 2d-block type if using MPI. * @param precision The precision of the output of DMK. * @param efs A vector containing the Fermi energies, and should have the same size as the number of SPIN. * @param ucell A pointer to the UnitCell object. * @param pv The Parallel_2D object. The 2d-block parallel information of DMK. */ template -void write_dmk( - const std::vector> dmk, - const int precision, - const std::vector efs, - const UnitCell* ucell, - const Parallel_2D& pv - ); -} +void write_dmk(const std::vector> dmk, + const int precision, + const std::vector efs, + const UnitCell* ucell, + const Parallel_2D& pv); +} // namespace ModuleIO #endif - diff --git a/source/module_io/test_serial/io_dmk_test.cpp b/source/module_io/test_serial/io_dmk_test.cpp index 8fda1123e2..b1e3853f06 100644 --- a/source/module_io/test_serial/io_dmk_test.cpp +++ b/source/module_io/test_serial/io_dmk_test.cpp @@ -1,24 +1,34 @@ -#include "gtest/gtest.h" -#include "gmock/gmock.h" #include "module_io/io_dmk.h" + #include "module_base/global_variable.h" #include "prepare_unitcell.h" +#include "gmock/gmock.h" +#include "gtest/gtest.h" + #ifdef __LCAO -InfoNonlocal::InfoNonlocal(){} -InfoNonlocal::~InfoNonlocal(){} -LCAO_Orbitals::LCAO_Orbitals(){} -LCAO_Orbitals::~LCAO_Orbitals(){} +InfoNonlocal::InfoNonlocal() +{ +} +InfoNonlocal::~InfoNonlocal() +{ +} +LCAO_Orbitals::LCAO_Orbitals() +{ +} +LCAO_Orbitals::~LCAO_Orbitals() +{ +} #endif Magnetism::Magnetism() { - this->tot_magnetization = 0.0; - this->abs_magnetization = 0.0; - this->start_magnetization = nullptr; + this->tot_magnetization = 0.0; + this->abs_magnetization = 0.0; + this->start_magnetization = nullptr; } Magnetism::~Magnetism() { - delete[] this->start_magnetization; + delete[] this->start_magnetization; } /************************************************ @@ -37,99 +47,97 @@ Magnetism::~Magnetism() class DMIOTest : public ::testing::Test { -protected: - int nspin = 1; - int lgd = 26; - int nnrg = 26*26; - double*** DM; - double** DM_R; - UnitCell* ucell; - void SetUp() - { - DM = new double**[nspin]; - DM_R = new double*[nspin]; - ucell = new UnitCell; - for(int is=0; is> dmk(nspin,std::vector(lgd*lgd,0.0)); - for(int is=0; is(nspin,ef),ucell,pv); - std::ifstream ifs; - ifs.open("SPIN1_DM"); - std::string str((std::istreambuf_iterator(ifs)),std::istreambuf_iterator()); - EXPECT_THAT(str, testing::HasSubstr("0.570336288802337 (fermi energy)")); - ifs.close(); - //remove("SPIN1_DM"); + // first read + GlobalV::MY_RANK = 0; + GlobalV::NLOCAL = lgd; + GlobalV::GAMMA_ONLY_LOCAL = true; + int is = 0; + std::string fn = "./support/SPIN1_DM"; + double ef; + UcellTestPrepare utp = UcellTestLib["Si"]; + ucell = utp.SetUcellInfo(); + ModuleIO::read_dmk(GlobalV::GAMMA_ONLY_LOCAL, GlobalV::NLOCAL, GlobalV::NSPIN, is, fn, DM, DM_R, ef, ucell); + EXPECT_DOUBLE_EQ(ef, 0.570336288802337); + EXPECT_NEAR(DM[0][0][0], 3.904e-01, 1e-6); + EXPECT_NEAR(DM[0][25][25], 3.445e-02, 1e-6); + // then write + int precision = 3; + std::vector> dmk(nspin, std::vector(lgd * lgd, 0.0)); + for (int is = 0; is < nspin; ++is) + { + for (int ig = 0; ig < lgd; ++ig) + { + for (int jg = 0; jg < lgd; ++jg) + { + dmk[is][ig * lgd + jg] = DM[is][ig][jg]; + } + } + } + Parallel_2D pv; + pv.nrow = lgd; + pv.ncol = lgd; + ModuleIO::write_dmk(dmk, precision, std::vector(nspin, ef), ucell, pv); + std::ifstream ifs; + ifs.open("SPIN1_DM"); + std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); + EXPECT_THAT(str, testing::HasSubstr("0.570336288802337 (fermi energy)")); + ifs.close(); + // remove("SPIN1_DM"); } From cae3aac7f6bf21ee75ad4d4d8d0c0153b3d8664d Mon Sep 17 00:00:00 2001 From: pxlxingliang Date: Fri, 28 Jun 2024 15:19:29 +0800 Subject: [PATCH 3/8] use awk to get key value of INPUT in catch_properties.sh --- tests/integrate/tools/catch_properties.sh | 73 +++++++++++++---------- 1 file changed, 40 insertions(+), 33 deletions(-) diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index 4fdf6c16ac..cf992207d6 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -21,46 +21,53 @@ sum_file(){ #echo $answer #exit 0 +get_input_key_value(){ + key=$1 + inputf=$2 + value=$(awk -v key=$key '{if($1==key) a=$2} END {print a}' $inputf) + echo $value +} + file=$1 #echo $1 calculation=`grep calculation INPUT | awk '{print $2}' | sed s/[[:space:]]//g` running_path=`echo "OUT.autotest/running_$calculation"".log"` natom=`grep -En '(^|[[:space:]])TOTAL ATOM NUMBER($|[[:space:]])' $running_path | tail -1 | awk '{print $6}'` -has_force=`grep -En '(^|[[:space:]])cal_force($|[[:space:]])' INPUT | awk '{print $2}'` -has_stress=`grep -En '(^|[[:space:]])cal_stress($|[[:space:]])' INPUT | awk '{print $2}'` -has_dftu=`grep -En '(^|[[:space:]])dft_plus_u($|[[:space:]])' INPUT | awk '{print $2}'` -has_band=`grep -En '(^|[[:space:]])out_band($|[[:space:]])' INPUT | awk '{print $2}'` -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` -has_lowf=`grep out_wfc_lcao INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_app_flag=`grep out_app_flag INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_wfc_r=`grep out_wfc_r INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_wfc_pw=`grep out_wfc_pw INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_dm=`grep "out_dm" INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_mul=`grep out_mul INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -gamma_only=`grep gamma_only INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -imp_sol=`grep imp_sol INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -run_rpa=`grep rpa INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_pot=`grep out_pot INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_dm1=`grep out_dm1 INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -get_s=`grep calculation INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_pband=`grep out_proj_band INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -toW90=`grep towannier90 INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_mat_r=`grep out_mat_r INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_mat_t=`grep out_mat_t INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_mat_dh=`grep out_mat_dh INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_scan=`grep dft_functional INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_chg=`grep out_chg INPUT | awk '{print $2}' | sed s/[[:space:]]//g` +has_force=$(get_input_key_value "cal_force" "INPUT") +has_stress=$(get_input_key_value "cal_stress" "INPUT") +has_dftu=$(get_input_key_value "dft_plus_u" "INPUT") +has_band=$(get_input_key_value "out_band" "INPUT") +has_dos=$(get_input_key_value "out_dos" "INPUT") +has_cond=$(get_input_key_value "cal_cond" "INPUT") +has_hs=$(get_input_key_value "out_mat_hs" "INPUT") +has_hs2=$(get_input_key_value "out_mat_hs2" "INPUT") +has_xc=$(get_input_key_value "out_mat_xc" "INPUT") +has_r=$(get_input_key_value "out_mat_r" "INPUT") +deepks_out_labels=$(get_input_key_value "deepks_out_labels" "INPUT") +deepks_bandgap=$(get_input_key_value "deepks_bandgap" "INPUT") +has_lowf=$(get_input_key_value "out_wfc_lcao" "INPUT") +out_app_flag=$(get_input_key_value "out_app_flag" "INPUT") +has_wfc_r=$(get_input_key_value "out_wfc_r" "INPUT") +has_wfc_pw=$(get_input_key_value "out_wfc_pw" "INPUT") +out_dm=$(get_input_key_value "out_dm" "INPUT") +out_mul=$(get_input_key_value "out_mul" "INPUT") +gamma_only=$(get_input_key_value "gamma_only" "INPUT") +imp_sol=$(get_input_key_value "imp_sol" "INPUT") +run_rpa=$(get_input_key_value "rpa" "INPUT") +out_pot=$(get_input_key_value "out_pot" "INPUT") +out_dm1=$(get_input_key_value "out_dm1" "INPUT") +get_s=$(get_input_key_value "calculation" "INPUT") +out_pband=$(get_input_key_value "out_proj_band" "INPUT") +toW90=$(get_input_key_value "towannier90" "INPUT") +has_mat_r=$(get_input_key_value "out_mat_r" "INPUT") +has_mat_t=$(get_input_key_value "out_mat_t" "INPUT") +has_mat_dh=$(get_input_key_value "out_mat_dh" "INPUT") +has_scan=$(get_input_key_value "dft_functional" "INPUT") +out_chg=$(get_input_key_value "out_chg" "INPUT") #echo $running_path -base=`grep -En '(^|[[:space:]])basis_type($|[[:space:]])' INPUT | awk '{print $2}' | sed s/[[:space:]]//g` +base=$(get_input_key_value "basis_type" "INPUT") word="driver_line" -symmetry=`grep -w "symmetry" INPUT | awk '{print $2}' | sed s/[[:space:]]//g` +symmetry=$(get_input_key_value "symmetry" "INPUT") test -e $1 && rm $1 #-------------------------------------------- # if NOT non-self-consistent calculations From 2ef59ab5ef571257b5878cc6413d75888bb56548 Mon Sep 17 00:00:00 2001 From: pxlxingliang Date: Fri, 28 Jun 2024 17:10:34 +0800 Subject: [PATCH 4/8] fix bug in catch_properties.sh --- tests/integrate/tools/catch_properties.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index cf992207d6..56ac8fb216 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -354,7 +354,7 @@ if ! test -z "$out_dm" && [ $out_dm == 1 ]; then echo "Can't find DM files" exit 1 else - python3 ../tools/CompareFile.py $dmref $dm 5 + python3 ../tools/CompareFile.py $dmref $dmfile 5 echo "DM_different $?" >>$1 fi fi From 34a3bbb25fbecbb5aa9be9769057363036a91b78 Mon Sep 17 00:00:00 2001 From: pxlxingliang Date: Mon, 1 Jul 2024 09:38:13 +0800 Subject: [PATCH 5/8] replace pointer to vector in read_cmk --- source/module_io/io_dmk.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/source/module_io/io_dmk.cpp b/source/module_io/io_dmk.cpp index 5fe919b3d2..5340eeffae 100644 --- a/source/module_io/io_dmk.cpp +++ b/source/module_io/io_dmk.cpp @@ -130,8 +130,7 @@ void ModuleIO::read_dmk( if (gamma_only_local) { - - double* tmp = new double[nlocal]; + std::vector tmp(nlocal); for (int i = 0; i < nlocal; ++i) { // GlobalV::ofs_running << " i=" << i << std::endl; @@ -143,7 +142,7 @@ void ModuleIO::read_dmk( ifs >> tmp[j]; } } - Parallel_Common::bcast_double(tmp, nlocal); + Parallel_Common::bcast_double(tmp.data(), nlocal); const int mu = trace_lo[i]; if (mu >= 0) @@ -158,7 +157,6 @@ void ModuleIO::read_dmk( } } } // i - delete[] tmp; } else { From 8988f8e6aed84bb26c466e538b0a87611f50aa1e Mon Sep 17 00:00:00 2001 From: pxlxingliang Date: Mon, 1 Jul 2024 09:44:08 +0800 Subject: [PATCH 6/8] fix --- source/module_io/io_dmk.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_io/io_dmk.cpp b/source/module_io/io_dmk.cpp index 5340eeffae..31499f2458 100644 --- a/source/module_io/io_dmk.cpp +++ b/source/module_io/io_dmk.cpp @@ -134,7 +134,7 @@ void ModuleIO::read_dmk( for (int i = 0; i < nlocal; ++i) { // GlobalV::ofs_running << " i=" << i << std::endl; - ModuleBase::GlobalFunc::ZEROS(tmp, nlocal); + ModuleBase::GlobalFunc::ZEROS(tmp.data(), nlocal); if (GlobalV::MY_RANK == 0) { for (int j = 0; j < nlocal; ++j) From e9bdf83be8f74cb42806ca48102206fa7fa8524c Mon Sep 17 00:00:00 2001 From: pxlxingliang Date: Tue, 2 Jul 2024 16:42:48 +0800 Subject: [PATCH 7/8] refactor the read_dmk interface --- source/module_io/io_dmk.cpp | 437 +++++++++++-------- source/module_io/io_dmk.h | 83 ++-- source/module_io/test_serial/io_dmk_test.cpp | 148 +++---- 3 files changed, 374 insertions(+), 294 deletions(-) diff --git a/source/module_io/io_dmk.cpp b/source/module_io/io_dmk.cpp index 31499f2458..24fd730a4c 100644 --- a/source/module_io/io_dmk.cpp +++ b/source/module_io/io_dmk.cpp @@ -4,191 +4,269 @@ #include "module_base/scalapack_connector.h" #include "module_base/timer.h" -void ModuleIO::read_dmk( -#ifdef __MPI - const int nnrg, - const int* trace_lo, -#endif - const bool gamma_only_local, - const int nlocal, - const int nspin, - const int& is, - const std::string& fn, - double*** DM, - double** DM_R, - double& ef, - const UnitCell* ucell) -{ - ModuleBase::TITLE("ModuleIO", "read_dmk"); - ModuleBase::timer::tick("ModuleIO", "read_dmk"); +/* +The format of the DMK file is as follows: +''' + + + + + + ... + ... +Direct + + +... + + +... + + + (fermi energy) + + + +... +''' + + +Example: +''' +sc + 5.29177 + 1 0 0 + 0 1 0 + 0 0 1 + H + 2 +Direct + 0 0 0.933859999999186 + 0 0 0.0661400000008143 + + 1 + -0.0883978533958687 (fermi energy) + 10 10 + + 5.773e-01 3.902e-02 1.661e-02 4.797e-17 -2.255e-17 5.773e-01 3.902e-02 -1.661e-02 + -1.461e-17 -4.414e-17 + ... + ''' + */ - GlobalV::ofs_running << "\n processor 0 is reading density matrix from file < " << fn << " > " << std::endl; - // xiaohui modify 2015-03-25 - // bool quit_mesia = false; - bool quit_abacus = false; - std::ifstream ifs; - if (GlobalV::MY_RANK == 0) +std::string ModuleIO::dmk_gen_fname(const bool gamma_only, const int ispin, const int ik) +{ + if (gamma_only) { - ifs.open(fn.c_str()); - if (!ifs) + return "SPIN" + std::to_string(ispin + 1) + "_DM"; + } + else + { + // this case is not implemented now. + ModuleBase::WARNING_QUIT("dmk_gen_fname", "Not implemented for non-gamma_only case."); + } +} + +void ModuleIO::dmk_write_ucell(std::ofstream& ofs, const UnitCell* ucell) +{ + // write the UnitCell information + ofs << ucell->latName << std::endl; + ofs << " " << ucell->lat0 * ModuleBase::BOHR_TO_A << std::endl; + ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " << ucell->latvec.e13 << std::endl; + ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " << ucell->latvec.e23 << std::endl; + ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " << ucell->latvec.e33 << std::endl; + for (int it = 0; it < ucell->ntype; it++) + { + ofs << " " << ucell->atoms[it].label; + } + ofs << std::endl; + for (int it = 0; it < ucell->ntype; it++) + { + ofs << " " << ucell->atoms[it].na; + } + ofs << std::endl; + ofs << "Direct" << std::endl; + for (int it = 0; it < ucell->ntype; it++) + { + Atom* atom = &ucell->atoms[it]; + ofs << std::setprecision(15); + for (int ia = 0; ia < ucell->atoms[it].na; ia++) { - // xiaohui modify 2015-03-25 - // quit_mesia = true; - quit_abacus = true; + ofs << " " << atom->taud[ia].x << " " << atom->taud[ia].y << " " << atom->taud[ia].z + << std::endl; } - else - { - // if the number is not match, - // quit the program or not. - bool quit = false; - - std::string name; - ifs >> name; - - // check lattice constant, unit is Angstrom - ModuleBase::CHECK_DOUBLE(ifs, ucell->lat0 * ModuleBase::BOHR_TO_A, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e11, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e12, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e13, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e21, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e22, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e23, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e31, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e32, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e33, quit); - - for (int it = 0; it < ucell->ntype; it++) - { - ModuleBase::CHECK_STRING(ifs, ucell->atoms[it].label, quit); - } + } +} - for (int it = 0; it < ucell->ntype; it++) - { - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].na, quit); - } +void ModuleIO::dmk_read_ucell(std::ifstream& ifs) +{ + std::string tmp; + for (int i=0;i<6;i++) + { + std::getline(ifs, tmp); // latName + lat0 + latvec + atom label + } + std::getline(ifs, tmp); // atom number of each type - std::string coordinate; - ifs >> coordinate; + std::istringstream iss(tmp); + int natom = 0; + int total_natom = 0; + while (iss >> natom) + { + total_natom += natom; + } + for (int i=0;intype; it++) - { - for (int ia = 0; ia < ucell->atoms[it].na; ia++) - { - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].x, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].y, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].z, quit); - } - } +void ModuleIO::dmk_readData(std::ifstream& ifs, double& data) +{ + ifs >> data; +} - ModuleBase::CHECK_INT(ifs, nspin); - ModuleBase::GlobalFunc::READ_VALUE(ifs, ef); - ModuleBase::CHECK_INT(ifs, nlocal); - ModuleBase::CHECK_INT(ifs, nlocal); - } // If file exist, read in data. - } // Finish reading the first part of density matrix. +void ModuleIO::dmk_readData(std::ifstream& ifs, std::complex& data) +{ + double real, imag; + ifs >> real; + ifs >> imag; + data = std::complex(real, imag); +} -#ifndef __MPI - GlobalV::ofs_running << " Read SPIN = " << is + 1 << " density matrix now." << std::endl; +template +bool ModuleIO::read_dmk( + const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>& dmk +) +{ + ModuleBase::TITLE("ModuleIO", "read_dmk"); + ModuleBase::timer::tick("ModuleIO", "read_dmk"); - if (gamma_only_local) - { - for (int i = 0; i < nlocal; ++i) - { - for (int j = 0; j < nlocal; ++j) - { - ifs >> DM[is][i][j]; - } - } - } - else - { + int my_rank = 0; #ifdef __MPI - ModuleBase::WARNING_QUIT("ModuleIO::read_dmk", "The nnrg should not be update"); - ModuleBase::CHECK_INT(ifs, nnrg); - - for (int i = 0; i < nnrg; ++i) + MPI_Comm_rank(pv.comm_2D, &my_rank); +#endif + + int nlocal = pv.get_global_row_size(); + bool gamma_only = std::is_same::value; + std::vector> dmk_global; + + // write a lambda function to check the consistency of the data + auto check_consistency = [&](const std::string& fn, const std::string& name, const std::string& value, const int& target) { + if (std::stoi(value) != target) { - ifs >> DM_R[is][i]; + ModuleBase::WARNING("ModuleIO::read_dmk", name + " is not consistent in file < " + fn + " >."); + std::cout << name << " = " << target << ", " << name << " in file = " << value << std::endl; + return false; } -#endif - } -#else + return true; + }; - // distribution of necessary data - // xiaohui modify 2015-03-25 - // Parallel_Common::bcast_bool(quit_mesia); - Parallel_Common::bcast_bool(quit_abacus); - // xiaohui modify 2015-03-25 - // if(quit_mesia) - if (quit_abacus) + bool read_success = true; + std::string tmp; + if (my_rank == 0) { - ModuleBase::WARNING_QUIT("ModuleIO::read_dmk", "Can not find the density matrix file."); - } + dmk_global.resize(nspin * nk, std::vector(nlocal*nlocal)); - Parallel_Common::bcast_double(ef); - - if (gamma_only_local) - { - std::vector tmp(nlocal); - for (int i = 0; i < nlocal; ++i) + for (int ispin = 0; ispin < nspin; ispin++) { - // GlobalV::ofs_running << " i=" << i << std::endl; - ModuleBase::GlobalFunc::ZEROS(tmp.data(), nlocal); - if (GlobalV::MY_RANK == 0) + for (int ik = 0; ik < nk; ik++) { - for (int j = 0; j < nlocal; ++j) + std::string fn = dmk_dir + dmk_gen_fname(gamma_only, ispin, ik); + std::ifstream ifs(fn.c_str()); + + if (!ifs) { - ifs >> tmp[j]; + ModuleBase::WARNING("ModuleIO::read_dmk", "Can't open DENSITY MATRIX File < " + fn + " >."); + read_success = false; + break; } - } - Parallel_Common::bcast_double(tmp.data(), nlocal); - const int mu = trace_lo[i]; - if (mu >= 0) - { - for (int j = 0; j < nlocal; ++j) + // read the UnitCell + dmk_read_ucell(ifs); + + ifs >> tmp; // nspin + if (!check_consistency(fn, "ispin", tmp, ispin+1)) + { + read_success = false; + ifs.close(); + break; + } + ifs >> tmp; ifs >> tmp; ifs >> tmp;// fermi energy + ifs >> tmp; // nlocal + if (!check_consistency(fn, "nlocal", tmp, nlocal)) { - const int nu = trace_lo[j]; - if (nu >= 0) + read_success = false; + ifs.close(); + break; + } + ifs >> tmp; // nlocal + if (!check_consistency(fn, "nlocal", tmp, nlocal)) + { + read_success = false; + ifs.close(); + break; + } + + // read the DMK data + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) { - DM[is][mu][nu] = tmp[j]; + dmk_readData(ifs, dmk_global[ik + nk * ispin][i * nlocal + j]); } } + ifs.close(); + } // ik + if (!read_success) + { + break; } - } // i - } - else - { - ModuleBase::WARNING_QUIT("ModuleIO::read_dmk", "not ready to readin DM_R"); - } -#endif - if (GlobalV::MY_RANK == 0) - ifs.close(); - - GlobalV::ofs_running << " Finish reading density matrix." << std::endl; + } // ispin + } // rank0 - ModuleBase::timer::tick("ModuleIO", "read_dmk"); - return; -} +#ifdef __MPI + MPI_Bcast(&read_success, 1, MPI_C_BOOL, 0, pv.comm_2D); +#endif -std::string ModuleIO::dmk_gen_fname(const bool gamma_only, const int ispin, const int ik) -{ - if (gamma_only) - { - return "SPIN" + std::to_string(ispin + 1) + "_DM"; - } - else + if (read_success) { - // this case is not implemented now. - ModuleBase::WARNING_QUIT("dmk_gen_fname", "Not implemented for non-gamma_only case."); +#ifdef __MPI + // seperate dmk data to each processor with 2D block distribution + dmk.resize(nspin * nk, std::vector(pv.get_row_size()*pv.get_col_size())); + Parallel_2D pv_glb; + pv_glb.set(nlocal, nlocal, nlocal, pv.comm_2D, pv.blacs_ctxt); + for(int ik=0;ik < nspin * nk; ik++) + { + Cpxgemr2d(nlocal, + nlocal, + dmk_global[ik].data(), + 1, + 1, + pv_glb.desc, + dmk[ik].data(), + 1, + 1, + const_cast(pv.desc), + pv_glb.blacs_ctxt); + } +#else + dmk = dmk_global; +#endif } + ModuleBase::timer::tick("ModuleIO", "read_dmk"); + return read_success; } + template -void ModuleIO::write_dmk(const std::vector> dmk, +void ModuleIO::write_dmk(const std::vector>& dmk, const int precision, - const std::vector efs, + const std::vector& efs, const UnitCell* ucell, const Parallel_2D& pv) { @@ -241,38 +319,16 @@ void ModuleIO::write_dmk(const std::vector> dmk, if (!ofs) { - ModuleBase::WARNING("ModuleIO::write_dmk", "Can't create DENSITY MATRIX File!"); + ModuleBase::WARNING("ModuleIO::write_dmk", "Can't create DENSITY MATRIX File < " + fn + " >."); + continue; } // write the UnitCell information - ofs << ucell->latName << std::endl; - ofs << " " << ucell->lat0 * ModuleBase::BOHR_TO_A << std::endl; - ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " << ucell->latvec.e13 << std::endl; - ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " << ucell->latvec.e23 << std::endl; - ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " << ucell->latvec.e33 << std::endl; - for (int it = 0; it < ucell->ntype; it++) - { - ofs << " " << ucell->atoms[it].label; - } - ofs << std::endl; - for (int it = 0; it < ucell->ntype; it++) - { - ofs << " " << ucell->atoms[it].na; - } - ofs << std::endl; - ofs << "Direct" << std::endl; - for (int it = 0; it < ucell->ntype; it++) - { - Atom* atom = &ucell->atoms[it]; - ofs << std::setprecision(15); - for (int ia = 0; ia < ucell->atoms[it].na; ia++) - { - ofs << " " << atom->taud[ia].x << " " << atom->taud[ia].y << " " << atom->taud[ia].z - << std::endl; - } - } + dmk_write_ucell(ofs, ucell); + + ofs << "\n " << dmk.size(); // nspin - ofs << "\n " << efs[ispin] << " (fermi energy)"; + ofs << "\n " << std::fixed << std::setprecision(5) << efs[ispin] << " (fermi energy)"; ofs << "\n " << nlocal << " " << nlocal << std::endl; ofs << std::setprecision(precision); @@ -285,14 +341,13 @@ void ModuleIO::write_dmk(const std::vector> dmk, { ofs << "\n"; } - if (std::is_same::value) + if(std::is_same::value) { ofs << " " << dmk_global[i * nlocal + j]; } else if (std::is_same, T>::value) { - ofs << " (" << std::real(dmk_global[i * nlocal + j]) << "," - << std::imag(dmk_global[i * nlocal + j]) << ")"; + ofs << " (" << std::real(dmk_global[i * nlocal + j]) << "," << std::imag(dmk_global[i * nlocal + j]) << ")"; } } } @@ -304,14 +359,26 @@ void ModuleIO::write_dmk(const std::vector> dmk, ModuleBase::timer::tick("ModuleIO", "write_dmk"); } -template void ModuleIO::write_dmk(const std::vector> dmk, +template bool ModuleIO::read_dmk(const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>& dmk); + +template bool ModuleIO::read_dmk>(const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>>& dmk); + +template void ModuleIO::write_dmk(const std::vector>& dmk, const int precision, - const std::vector efs, + const std::vector& efs, const UnitCell* ucell, const Parallel_2D& pv); -template void ModuleIO::write_dmk>(const std::vector>> dmk, +template void ModuleIO::write_dmk>(const std::vector>>& dmk, const int precision, - const std::vector efs, + const std::vector& efs, const UnitCell* ucell, const Parallel_2D& pv); \ No newline at end of file diff --git a/source/module_io/io_dmk.h b/source/module_io/io_dmk.h index 91afd2ee07..9af9ef1551 100644 --- a/source/module_io/io_dmk.h +++ b/source/module_io/io_dmk.h @@ -5,49 +5,64 @@ #include "module_cell/unitcell.h" #include +#include namespace ModuleIO { -/** - * @brief Reads the DMK file and stores the data in the provided arrays. - * - * @param nnrg The number of energy points. - * @param trace_lo An array containing the lower trace indices. - * @param gamma_only_local A boolean indicating whether to read only the gamma-only part of the DMK file. - * @param nlocal The number of local orbitals. - * @param nspin The number of spin components. - * @param is The index of the spin component. - * @param fn The filename of the DMK file. - * @param DM A 3D array to store the DMK data. - * @param DM_R A 2D array to store the real part of the DMK data. - * @param ef The Fermi energy. - * @param ucell A pointer to the UnitCell object. - */ -void read_dmk( -#ifdef __MPI - const int nnrg, - const int* trace_lo, -#endif - const bool gamma_only_local, - const int nlocal, - const int nspin, - const int& is, - const std::string& fn, - double*** DM, - double** DM_R, - double& ef, - const UnitCell* ucell); /** * @brief Generates the filename for the DMK file based on the given parameters. * - * @param gamma_only A boolean indicating whether to generate the filename for the gamma-only part. + * @param gamma_only Whether the calculation is gamma_only. * @param ispin The index of the spin component. * @param ik The index of the k-point. * @return The generated filename. */ std::string dmk_gen_fname(const bool gamma_only, const int ispin, const int ik); +/** + * @brief Writes the unit cell information to a DMK file. + * + * @param ofs The output file stream. + * @param ucell A pointer to the UnitCell object. + */ +void dmk_write_ucell(std::ofstream& ofs, const UnitCell* ucell); + +/** + * @brief Reads the unit cell information lines in a DMK file. + * + * @param ifs The input file stream. + */ +void dmk_read_ucell(std::ifstream& ifs); + +/** + * @brief Read one double from a file. + */ +void dmk_readData(std::ifstream& ifs, double& data); + +/** + * @brief Read one complex double from a file. + */ +void dmk_readData(std::ifstream& ifs, std::complex& data); + +/** + * @brief Reads the DMK data from a file. + * + * @tparam T The type of the DMK data. + * @param nspin The number of spin components. + * @param nk The number of k-points. + * @param pv The Parallel_2D object. Will get the global size and local size from it, and seperate the data into different processors accordingly. + * @param dmk_dir The directory path of the DMK file. + * @param dmk A vector to store the DMK data. If use MPI parallel, the data will be seperated into different processors based on the Parallel_2D object. + * @return True if the DMK data is successfully read, false otherwise. + */ +template +bool read_dmk(const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>& dmk); + /** * @brief Writes the DMK data to a file. * @@ -60,10 +75,12 @@ std::string dmk_gen_fname(const bool gamma_only, const int ispin, const int ik); * @param pv The Parallel_2D object. The 2d-block parallel information of DMK. */ template -void write_dmk(const std::vector> dmk, +void write_dmk(const std::vector>& dmk, const int precision, - const std::vector efs, + const std::vector& efs, const UnitCell* ucell, const Parallel_2D& pv); + } // namespace ModuleIO -#endif + +#endif // IO_DMK_H diff --git a/source/module_io/test_serial/io_dmk_test.cpp b/source/module_io/test_serial/io_dmk_test.cpp index b1e3853f06..f8c5ff11cf 100644 --- a/source/module_io/test_serial/io_dmk_test.cpp +++ b/source/module_io/test_serial/io_dmk_test.cpp @@ -45,99 +45,95 @@ Magnetism::~Magnetism() * - the serial version without MPI */ -class DMIOTest : public ::testing::Test +TEST(DMKTest,GenFileName) +{ + std::string fname = ModuleIO::dmk_gen_fname(true, 0, 0); + EXPECT_EQ(fname, "SPIN1_DM"); + fname = ModuleIO::dmk_gen_fname(true, 1, 1); + EXPECT_EQ(fname, "SPIN2_DM"); + + // do not support non-gamma-only case now + std::string output; + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleIO::dmk_gen_fname(false, 2, 0), ::testing::ExitedWithCode(0), ""); + output = testing::internal::GetCapturedStdout(); +}; + +class DMKIOTest : public ::testing::Test { protected: - int nspin = 1; - int lgd = 26; - int nnrg = 26 * 26; - double*** DM; - double** DM_R; - UnitCell* ucell; + int nspin = 2; + int nk = 1; + int nlocal = 20; + std::vector> dmk; + Parallel_2D pv; + std::vector efs; + void SetUp() { - DM = new double**[nspin]; - DM_R = new double*[nspin]; - ucell = new UnitCell; - for (int is = 0; is < nspin; ++is) + dmk.resize(nspin*nk, std::vector(nlocal * nlocal, 0.0)); + for (int i=0;i> dmk(nspin, std::vector(lgd * lgd, 0.0)); - for (int is = 0; is < nspin; ++is) - { - for (int ig = 0; ig < lgd; ++ig) - { - for (int jg = 0; jg < lgd; ++jg) - { - dmk[is][ig * lgd + jg] = DM[is][ig][jg]; - } - } - } - Parallel_2D pv; - pv.nrow = lgd; - pv.ncol = lgd; - ModuleIO::write_dmk(dmk, precision, std::vector(nspin, ef), ucell, pv); + ModuleIO::write_dmk(dmk, 3, efs, ucell, pv); std::ifstream ifs; - ifs.open("SPIN1_DM"); + + std::string fn = "SPIN1_DM"; + ifs.open(fn); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); - EXPECT_THAT(str, testing::HasSubstr("0.570336288802337 (fermi energy)")); + EXPECT_THAT(str, testing::HasSubstr("0.00000 (fermi energy)")); + EXPECT_THAT(str, testing::HasSubstr("20 20")); + EXPECT_THAT(str, testing::HasSubstr("0.000e+00 1.000e-01 2.000e-01 3.000e-01 4.000e-01 5.000e-01 6.000e-01 7.000e-01\n")); + EXPECT_THAT(str, testing::HasSubstr("8.000e-01 9.000e-01 1.000e+00 1.100e+00 1.200e+00 1.300e+00 1.400e+00 1.500e+00\n")); + EXPECT_THAT(str, testing::HasSubstr("1.600e+00 1.700e+00 1.800e+00 1.900e+00\n")); + ifs.close(); + + fn = "SPIN2_DM"; + ifs.open(fn); + str = std::string((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); + EXPECT_THAT(str, testing::HasSubstr("0.10000 (fermi energy)")); + EXPECT_THAT(str, testing::HasSubstr("20 20")); + EXPECT_THAT(str, testing::HasSubstr("1.000e+00 1.100e+00 1.200e+00 1.300e+00 1.400e+00 1.500e+00 1.600e+00 1.700e+00\n")); + EXPECT_THAT(str, testing::HasSubstr("1.800e+00 1.900e+00 2.000e+00 2.100e+00 2.200e+00 2.300e+00 2.400e+00 2.500e+00\n")); + EXPECT_THAT(str, testing::HasSubstr("2.600e+00 2.700e+00 2.800e+00 2.900e+00\n")); ifs.close(); - // remove("SPIN1_DM"); + + delete ucell; + // remove the generated files + remove("SPIN1_DM"); + remove("SPIN2_DM"); +}; + +TEST_F(DMKIOTest,ReadDMK) +{ + pv.nrow = 26; + pv.ncol = 26; + EXPECT_TRUE(ModuleIO::read_dmk(1, 1, pv, "./support/", dmk)); + EXPECT_EQ(dmk.size(), 1); + EXPECT_EQ(dmk[0].size(), 26*26); + EXPECT_NEAR(dmk[0][0], 3.904e-01, 1e-6); + EXPECT_NEAR(dmk[0][25*26+25], 3.445e-02, 1e-6); } From 2c086c2ad1c078b67db497a95bdf94aeb9edfece Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Tue, 2 Jul 2024 15:26:50 +0000 Subject: [PATCH 8/8] [pre-commit.ci lite] apply automatic fixes --- source/module_elecstate/elecstate_lcao.cpp | 144 +-- source/module_esolver/esolver_ks_lcao.cpp | 883 +++++++++--------- source/module_esolver/esolver_ks_lcao.h | 22 +- .../module_esolver/esolver_ks_lcao_elec.cpp | 421 +++++---- .../module_esolver/esolver_ks_lcao_tddft.cpp | 334 +++---- source/module_io/io_dmk.cpp | 258 +++-- source/module_io/io_dmk.h | 17 +- source/module_io/test_serial/io_dmk_test.cpp | 96 +- 8 files changed, 1107 insertions(+), 1068 deletions(-) diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index 1156f67dea..85a6791f30 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -10,46 +10,51 @@ #include -namespace elecstate -{ +namespace elecstate { // multi-k case template <> -void ElecStateLCAO>::psiToRho(const psi::Psi>& psi) -{ +void ElecStateLCAO>::psiToRho( + const psi::Psi>& psi) { ModuleBase::TITLE("ElecStateLCAO", "psiToRho"); ModuleBase::timer::tick("ElecStateLCAO", "psiToRho"); this->calculate_weights(); - // the calculations of dm, and dm -> rho are, technically, two separate functionalities, as we cannot - // rule out the possibility that we may have a dm from other sources, such as read from file. - // However, since we are not separating them now, I opt to add a flag to control how dm is obtained as of now - if (!GlobalV::dm_to_rho) - { + // the calculations of dm, and dm -> rho are, technically, two separate + // functionalities, as we cannot rule out the possibility that we may have a + // dm from other sources, such as read from file. However, since we are not + // separating them now, I opt to add a flag to control how dm is obtained as + // of now + if (!GlobalV::dm_to_rho) { this->calEBand(); ModuleBase::GlobalFunc::NOTE("Calculate the density matrix."); - // this part for calculating DMK in 2d-block format, not used for charge now + // this part for calculating DMK in 2d-block format, not used for charge + // now // psi::Psi> dm_k_2d(); - if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx" || GlobalV::KS_SOLVER == "lapack" - || GlobalV::KS_SOLVER == "cusolver" || GlobalV::KS_SOLVER == "cusolvermp" + if (GlobalV::KS_SOLVER == "genelpa" + || GlobalV::KS_SOLVER == "scalapack_gvx" + || GlobalV::KS_SOLVER == "lapack" + || GlobalV::KS_SOLVER == "cusolver" + || GlobalV::KS_SOLVER == "cusolvermp" || GlobalV::KS_SOLVER == "cg_in_lcao") // Peize Lin test 2019-05-15 { // cal_dm(this->loc->ParaV, this->wg, psi, this->loc->dm_k); - elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM)); + elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), + this->wg, + psi, + *(this->DM)); this->DM->cal_DMR(); // interface for RI-related calculation, which needs loc.dm_k #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) - { + if (GlobalC::exx_info.info_global.cal_exx) { const K_Vectors* kv = this->DM->get_kv_pointer(); this->loc->dm_k.resize(kv->get_nks()); - for (int ik = 0; ik < kv->get_nks(); ++ik) - { + for (int ik = 0; ik < kv->get_nks(); ++ik) { this->loc->set_dm_k(ik, this->DM->get_DMK_pointer(ik)); } } @@ -57,9 +62,9 @@ void ElecStateLCAO>::psiToRho(const psi::Psicharge->rho[is], this->charge->nrxx); // mohan 2009-11-10 + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], + this->charge->nrxx); // mohan 2009-11-10 } //------------------------------------------------------------ @@ -67,15 +72,16 @@ void ElecStateLCAO>::psiToRho(const psi::Psigint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + this->gint_k->transfer_DM2DtoGrid( + this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho); this->gint_k->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx); + if (XC_Functional::get_func_type() == 3 + || XC_Functional::get_func_type() == 5) { + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], + this->charge->nrxx); } Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau); this->gint_k->cal_gint(&inout1); @@ -89,28 +95,31 @@ void ElecStateLCAO>::psiToRho(const psi::Psi -void ElecStateLCAO::psiToRho(const psi::Psi& psi) -{ +void ElecStateLCAO::psiToRho(const psi::Psi& psi) { ModuleBase::TITLE("ElecStateLCAO", "psiToRho"); ModuleBase::timer::tick("ElecStateLCAO", "psiToRho"); this->calculate_weights(); this->calEBand(); - if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx" || GlobalV::KS_SOLVER == "lapack" - || GlobalV::KS_SOLVER == "cusolver" || GlobalV::KS_SOLVER == "cusolvermp" || GlobalV::KS_SOLVER == "cg_in_lcao") - { + if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx" + || GlobalV::KS_SOLVER == "lapack" || GlobalV::KS_SOLVER == "cusolver" + || GlobalV::KS_SOLVER == "cusolvermp" + || GlobalV::KS_SOLVER == "cg_in_lcao") { ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d"); // get DMK in 2d-block format // cal_dm(this->loc->ParaV, this->wg, psi, this->loc->dm_gamma); - elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM)); + elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), + this->wg, + psi, + *(this->DM)); this->DM->cal_DMR(); } - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx); // mohan 2009-11-10 + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], + this->charge->nrxx); // mohan 2009-11-10 } //------------------------------------------------------------ @@ -118,17 +127,18 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) //------------------------------------------------------------ ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!"); - this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + this->gint_gamma->transfer_DM2DtoGrid( + this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho); this->gint_gamma->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx); + if (XC_Functional::get_func_type() == 3 + || XC_Functional::get_func_type() == 5) { + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], + this->charge->nrxx); } Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau); this->gint_gamma->cal_gint(&inout1); @@ -141,21 +151,21 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) } template -void ElecStateLCAO::init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin) -{ +void ElecStateLCAO::init_DM(const K_Vectors* kv, + const Parallel_Orbitals* paraV, + const int nspin) { this->DM = new DensityMatrix(kv, paraV, nspin); } template <> -double ElecStateLCAO::get_spin_constrain_energy() -{ - SpinConstrain& sc = SpinConstrain::getScInstance(); +double ElecStateLCAO::get_spin_constrain_energy() { + SpinConstrain& sc + = SpinConstrain::getScInstance(); return sc.cal_escon(); } template <> -double ElecStateLCAO>::get_spin_constrain_energy() -{ +double ElecStateLCAO>::get_spin_constrain_energy() { SpinConstrain, base_device::DEVICE_CPU>& sc = SpinConstrain>::getScInstance(); return sc.cal_escon(); @@ -163,38 +173,37 @@ double ElecStateLCAO>::get_spin_constrain_energy() #ifdef __PEXSI template <> -void ElecStateLCAO::dmToRho(std::vector pexsi_DM, std::vector pexsi_EDM) -{ +void ElecStateLCAO::dmToRho(std::vector pexsi_DM, + std::vector pexsi_EDM) { ModuleBase::timer::tick("ElecStateLCAO", "dmToRho"); int nspin = GlobalV::NSPIN; - if (GlobalV::NSPIN == 4) - { + if (GlobalV::NSPIN == 4) { nspin = 1; } this->get_DM()->pexsi_EDM = pexsi_EDM; - for (int is = 0; is < nspin; is++) - { + for (int is = 0; is < nspin; is++) { this->DM->set_DMK_pointer(is, pexsi_DM[is]); } DM->cal_DMR(); - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx); // mohan 2009-11-10 + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], + this->charge->nrxx); // mohan 2009-11-10 } ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!"); - this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + this->gint_gamma->transfer_DM2DtoGrid( + this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho); this->gint_gamma->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0], this->charge->nrxx); + if (XC_Functional::get_func_type() == 3 + || XC_Functional::get_func_type() == 5) { + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0], + this->charge->nrxx); } Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau); this->gint_gamma->cal_gint(&inout1); @@ -207,10 +216,11 @@ void ElecStateLCAO::dmToRho(std::vector pexsi_DM, std::vector -void ElecStateLCAO>::dmToRho(std::vector*> pexsi_DM, - std::vector*> pexsi_EDM) -{ - ModuleBase::WARNING_QUIT("ElecStateLCAO", "pexsi is not completed for multi-k case"); +void ElecStateLCAO>::dmToRho( + std::vector*> pexsi_DM, + std::vector*> pexsi_EDM) { + ModuleBase::WARNING_QUIT("ElecStateLCAO", + "pexsi is not completed for multi-k case"); } #endif diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 32106a9fbd..0a8ac41f9e 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -48,31 +48,31 @@ #include "module_io/write_dmr.h" #include "module_io/write_wfc_nao.h" -namespace ModuleESolver -{ +namespace ModuleESolver { //------------------------------------------------------------------------------ //! the 1st function of ESolver_KS_LCAO: constructor //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -ESolver_KS_LCAO::ESolver_KS_LCAO() -{ +ESolver_KS_LCAO::ESolver_KS_LCAO() { this->classname = "ESolver_KS_LCAO"; this->basisname = "LCAO"; // the following EXX code should be removed to other places, mohan 2024/05/11 #ifdef __EXX - if (GlobalC::exx_info.info_ri.real_number) - { - this->exx_lri_double = std::make_shared>(GlobalC::exx_info.info_ri); - this->exd = std::make_shared>(this->exx_lri_double); + if (GlobalC::exx_info.info_ri.real_number) { + this->exx_lri_double + = std::make_shared>(GlobalC::exx_info.info_ri); + this->exd = std::make_shared>( + this->exx_lri_double); this->LM.Hexxd = &this->exd->get_Hexxs(); - } - else - { - this->exx_lri_complex = std::make_shared>>(GlobalC::exx_info.info_ri); - this->exc = std::make_shared>>(this->exx_lri_complex); + } else { + this->exx_lri_complex = std::make_shared>>( + GlobalC::exx_info.info_ri); + this->exc + = std::make_shared>>( + this->exx_lri_complex); this->LM.Hexxc = &this->exc->get_Hexxs(); } #endif @@ -83,9 +83,7 @@ ESolver_KS_LCAO::ESolver_KS_LCAO() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -ESolver_KS_LCAO::~ESolver_KS_LCAO() -{ -} +ESolver_KS_LCAO::~ESolver_KS_LCAO() {} //------------------------------------------------------------------------------ //! the 3rd function of ESolver_KS_LCAO: init @@ -106,50 +104,53 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() //! 14) set occupations? //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) -{ +void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS_LCAO", "before_all_runners"); ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); // 1) calculate overlap matrix S - if (GlobalV::CALCULATION == "get_S") - { + if (GlobalV::CALCULATION == "get_S") { // 1.1) read pseudopotentials ucell.read_pseudo(GlobalV::ofs_running); // 1.2) symmetrize things - if (ModuleSymmetry::Symmetry::symm_flag == 1) - { - ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running); + if (ModuleSymmetry::Symmetry::symm_flag == 1) { + ucell.symm.analy_sys(ucell.lat, + ucell.st, + ucell.atoms, + GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } // 1.3) Setup k-points according to symmetry. - this->kv - .set(ucell.symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec, GlobalV::ofs_running); + this->kv.set(ucell.symm, + GlobalV::global_kpoint_card, + GlobalV::NSPIN, + ucell.G, + ucell.latvec, + GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); Print_Info::setup_parameters(ucell, this->kv); - } - else - { + } else { // 1) else, call before_all_runners() in ESolver_KS ESolver_KS::before_all_runners(inp, ucell); } // end ifnot get_S // 2) init ElecState - // autoset nbands in ElecState, it should before basis_init (for Psi 2d divid) - if (this->pelec == nullptr) - { + // autoset nbands in ElecState, it should before basis_init (for Psi 2d + // divid) + if (this->pelec == nullptr) { // TK stands for double and complex? - this->pelec = new elecstate::ElecStateLCAO(&(this->chr), // use which parameter? - &(this->kv), - this->kv.get_nks(), - &(this->LOC), // use which parameter? - &(this->GG), // mohan add 2024-04-01 - &(this->GK), // mohan add 2024-04-01 - this->pw_rho, - this->pw_big); + this->pelec = new elecstate::ElecStateLCAO( + &(this->chr), // use which parameter? + &(this->kv), + this->kv.get_nks(), + &(this->LOC), // use which parameter? + &(this->GG), // mohan add 2024-04-01 + &(this->GK), // mohan add 2024-04-01 + this->pw_rho, + this->pw_big); } // 3) init LCAO basis @@ -167,8 +168,7 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) ->init_DM(&this->kv, &(this->orb_con.ParaV), GlobalV::NSPIN); // this function should be removed outside of the function - if (GlobalV::CALCULATION == "get_S") - { + if (GlobalV::CALCULATION == "get_S") { ModuleBase::timer::tick("ESolver_KS_LCAO", "init"); return; } @@ -176,36 +176,33 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) // 6) initialize Hamilt in LCAO // * allocate H and S matrices according to computational resources // * set the 'trace' between local H/S and global H/S - this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, orb_con.ParaV, this->kv.get_nks()); + this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, + orb_con.ParaV, + this->kv.get_nks()); #ifdef __EXX // 7) initialize exx // PLEASE simplify the Exx_Global interface - if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "relax" || GlobalV::CALCULATION == "cell-relax" - || GlobalV::CALCULATION == "md") - { - if (GlobalC::exx_info.info_global.cal_exx) - { + if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "relax" + || GlobalV::CALCULATION == "cell-relax" + || GlobalV::CALCULATION == "md") { + if (GlobalC::exx_info.info_global.cal_exx) { /* In the special "two-level" calculation case, - first scf iteration only calculate the functional without exact exchange. - but in "nscf" calculation, there is no need of "two-level" method. */ - if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "PBE0" - || ucell.atoms[0].ncpp.xc_func == "HSE") - { + first scf iteration only calculate the functional without exact + exchange. but in "nscf" calculation, there is no need of "two-level" + method. */ + if (ucell.atoms[0].ncpp.xc_func == "HF" + || ucell.atoms[0].ncpp.xc_func == "PBE0" + || ucell.atoms[0].ncpp.xc_func == "HSE") { XC_Functional::set_xc_type("pbe"); - } - else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") - { + } else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") { XC_Functional::set_xc_type("scan"); } // GlobalC::exx_lcao.init(); - if (GlobalC::exx_info.info_ri.real_number) - { + if (GlobalC::exx_info.info_ri.real_number) { this->exx_lri_double->init(MPI_COMM_WORLD, this->kv); - } - else - { + } else { this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv); } } @@ -213,8 +210,7 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) #endif // 8) initialize DFT+U - if (GlobalV::dft_plus_u) - { + if (GlobalV::dft_plus_u) { GlobalC::dftu.init(ucell, this->LM, this->kv.get_nks()); } @@ -222,8 +218,7 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho); // 10) initialize the HSolver - if (this->phsol == nullptr) - { + if (this->phsol == nullptr) { this->phsol = new hsolver::HSolverLCAO(&(this->orb_con.ParaV)); this->phsol->method = GlobalV::KS_SOLVER; } @@ -233,8 +228,7 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) this->pelec->omega = GlobalC::ucell.omega; // 12) initialize the potential - if (this->pelec->pot == nullptr) - { + if (this->pelec->pot == nullptr) { this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho, &GlobalC::ucell, @@ -246,17 +240,17 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) #ifdef __DEEPKS // 13) initialize deepks - if (GlobalV::deepks_scf) - { + if (GlobalV::deepks_scf) { // load the DeePKS model from deep neural network GlobalC::ld.load_model(INPUT.deepks_model); } #endif // 14) set occupations - if (GlobalV::ocp) - { - this->pelec->fixed_weights(GlobalV::ocp_kb, GlobalV::NBANDS, GlobalV::nelec); + if (GlobalV::ocp) { + this->pelec->fixed_weights(GlobalV::ocp_kb, + GlobalV::NBANDS, + GlobalV::nelec); } ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); @@ -268,25 +262,25 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) -{ +void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS_LCAO", "init_after_vc"); ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc"); ESolver_KS::init_after_vc(inp, ucell); - if (GlobalV::md_prec_level == 2) - { + if (GlobalV::md_prec_level == 2) { delete this->pelec; - this->pelec = new elecstate::ElecStateLCAO(&(this->chr), - &(this->kv), - this->kv.get_nks(), - &(this->LOC), - &(this->GG), // mohan add 2024-04-01 - &(this->GK), // mohan add 2024-04-01 - this->pw_rho, - this->pw_big); - - dynamic_cast*>(this->pelec)->init_DM(&this->kv, this->LM.ParaV, GlobalV::NSPIN); + this->pelec = new elecstate::ElecStateLCAO( + &(this->chr), + &(this->kv), + this->kv.get_nks(), + &(this->LOC), + &(this->GG), // mohan add 2024-04-01 + &(this->GK), // mohan add 2024-04-01 + this->pw_rho, + this->pw_big); + + dynamic_cast*>(this->pelec) + ->init_DM(&this->kv, this->LM.ParaV, GlobalV::NSPIN); GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho); @@ -294,15 +288,15 @@ void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) this->pelec->omega = GlobalC::ucell.omega; // Initialize the potential. - if (this->pelec->pot == nullptr) - { - this->pelec->pot = new elecstate::Potential(this->pw_rhod, - this->pw_rho, - &GlobalC::ucell, - &(GlobalC::ppcell.vloc), - &(this->sf), - &(this->pelec->f_en.etxc), - &(this->pelec->f_en.vtxc)); + if (this->pelec->pot == nullptr) { + this->pelec->pot + = new elecstate::Potential(this->pw_rhod, + this->pw_rho, + &GlobalC::ucell, + &(GlobalC::ppcell.vloc), + &(this->sf), + &(this->pelec->f_en.etxc), + &(this->pelec->f_en.vtxc)); } } @@ -315,8 +309,7 @@ void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -double ESolver_KS_LCAO::cal_energy() -{ +double ESolver_KS_LCAO::cal_energy() { ModuleBase::TITLE("ESolver_KS_LCAO", "cal_energy"); return this->pelec->f_en.etot; @@ -327,8 +320,7 @@ double ESolver_KS_LCAO::cal_energy() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) -{ +void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) { ModuleBase::TITLE("ESolver_KS_LCAO", "cal_force"); ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_force"); @@ -370,13 +362,11 @@ void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) -{ +void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) { ModuleBase::TITLE("ESolver_KS_LCAO", "cal_stress"); ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress"); - if (!this->have_force) - { + if (!this->have_force) { ModuleBase::matrix fcs; this->cal_force(fcs); } @@ -391,45 +381,68 @@ void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::after_all_runners() -{ +void ESolver_KS_LCAO::after_all_runners() { ModuleBase::TITLE("ESolver_KS_LCAO", "after_all_runners"); ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); - GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; + GlobalV::ofs_running << "\n\n --------------------------------------------" + << std::endl; GlobalV::ofs_running << std::setprecision(16); - GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl; - GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl; - - if (INPUT.out_dos != 0 || INPUT.out_band[0] != 0 || INPUT.out_proj_band != 0) - { + GlobalV::ofs_running << " !FINAL_ETOT_IS " + << this->pelec->f_en.etot * ModuleBase::Ry_to_eV + << " eV" << std::endl; + GlobalV::ofs_running << " --------------------------------------------\n\n" + << std::endl; + + if (INPUT.out_dos != 0 || INPUT.out_band[0] != 0 + || INPUT.out_proj_band != 0) { GlobalV::ofs_running << "\n\n\n\n"; - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - GlobalV::ofs_running << " | |" << std::endl; - GlobalV::ofs_running << " | Post-processing of data: |" << std::endl; - GlobalV::ofs_running << " | DOS (density of states) and bands will be output here. |" << std::endl; - GlobalV::ofs_running << " | If atomic orbitals are used, Mulliken charge analysis can be done. |" << std::endl; - GlobalV::ofs_running << " | Also the .bxsf file containing fermi surface information can be |" << std::endl; - GlobalV::ofs_running << " | done here. |" << std::endl; - GlobalV::ofs_running << " | |" << std::endl; - GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" + ">>>>>>>>>>>>>>>>>>>>>>>>>" + << std::endl; + GlobalV::ofs_running << " | " + " |" + << std::endl; + GlobalV::ofs_running << " | Post-processing of data: " + " |" + << std::endl; + GlobalV::ofs_running << " | DOS (density of states) and bands will be " + "output here. |" + << std::endl; + GlobalV::ofs_running << " | If atomic orbitals are used, Mulliken " + "charge analysis can be done. |" + << std::endl; + GlobalV::ofs_running << " | Also the .bxsf file containing fermi " + "surface information can be |" + << std::endl; + GlobalV::ofs_running << " | done here. " + " |" + << std::endl; + GlobalV::ofs_running << " | " + " |" + << std::endl; + GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" + "<<<<<<<<<<<<<<<<<<<<<<<<<" + << std::endl; GlobalV::ofs_running << "\n\n\n\n"; } // qianrui modify 2020-10-18 - if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "md" || GlobalV::CALCULATION == "relax") - { - ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv, &(GlobalC::Pkpoints)); + if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "md" + || GlobalV::CALCULATION == "relax") { + ModuleIO::write_istate_info(this->pelec->ekb, + this->pelec->wg, + this->kv, + &(GlobalC::Pkpoints)); } const int nspin0 = (GlobalV::NSPIN == 2) ? 2 : 1; - if (INPUT.out_band[0]) - { - for (int is = 0; is < nspin0; is++) - { + if (INPUT.out_band[0]) { + for (int is = 0; is < nspin0; is++) { std::stringstream ss2; ss2 << GlobalV::global_out_dir << "BANDS_" << is + 1 << ".dat"; - GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() << std::endl; + GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() + << std::endl; ModuleIO::nscf_band(is, ss2.str(), GlobalV::NBANDS, @@ -443,11 +456,15 @@ void ESolver_KS_LCAO::after_all_runners() if (INPUT.out_proj_band) // Projeced band structure added by jiyy-2022-4-20 { - ModuleIO::write_proj_band_lcao(this->psi, this->LM, this->pelec, this->kv, GlobalC::ucell, this->p_hamilt); + ModuleIO::write_proj_band_lcao(this->psi, + this->LM, + this->pelec, + this->kv, + GlobalC::ucell, + this->p_hamilt); } - if (INPUT.out_dos) - { + if (INPUT.out_dos) { ModuleIO::out_dos_nao(this->psi, this->LM, this->orb_con.ParaV, @@ -464,8 +481,7 @@ void ESolver_KS_LCAO::after_all_runners() this->p_hamilt); } - if (INPUT.out_mat_xc) - { + if (INPUT.out_mat_xc) { ModuleIO::write_Vxc(GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::DRANK, @@ -492,23 +508,20 @@ void ESolver_KS_LCAO::after_all_runners() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, Input& inp, UnitCell& ucell) -{ +void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, + Input& inp, + UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS_LCAO", "init_basis_lcao"); // autoset NB2D first - if (GlobalV::NB2D == 0) - { - if (GlobalV::NLOCAL > 0) - { + if (GlobalV::NB2D == 0) { + if (GlobalV::NLOCAL > 0) { GlobalV::NB2D = (GlobalV::NSPIN == 4) ? 2 : 1; } - if (GlobalV::NLOCAL > 500) - { + if (GlobalV::NLOCAL > 500) { GlobalV::NB2D = 32; } - if (GlobalV::NLOCAL > 1000) - { + if (GlobalV::NLOCAL > 1000) { GlobalV::NB2D = 64; } } @@ -530,15 +543,24 @@ void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, Input& inp, // * construct the interpolation tables. two_center_bundle_.build_orb(ucell.ntype, ucell.orbital_fn); - two_center_bundle_.build_alpha(GlobalV::deepks_setorb, &ucell.descriptor_file); + two_center_bundle_.build_alpha(GlobalV::deepks_setorb, + &ucell.descriptor_file); two_center_bundle_.build_orb_onsite(ucell.ntype, GlobalV::onsite_radius); - // currently deepks only use one descriptor file, so cast bool to int is fine + // currently deepks only use one descriptor file, so cast bool to int is + // fine // TODO Due to the omnipresence of GlobalC::ORB, we still have to rely // on the old interface for now. - two_center_bundle_.to_LCAO_Orbitals(GlobalC::ORB, inp.lcao_ecut, inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax); + two_center_bundle_.to_LCAO_Orbitals(GlobalC::ORB, + inp.lcao_ecut, + inp.lcao_dk, + inp.lcao_dr, + inp.lcao_rmax); - ucell.infoNL.setupNonlocal(ucell.ntype, ucell.atoms, GlobalV::ofs_running, GlobalC::ORB); + ucell.infoNL.setupNonlocal(ucell.ntype, + ucell.atoms, + GlobalV::ofs_running, + GlobalC::ORB); two_center_bundle_.build_beta(ucell.ntype, ucell.infoNL.Beta); @@ -550,13 +572,18 @@ void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, Input& inp, #ifdef USE_NEW_TWO_CENTER two_center_bundle_.tabulate(); #else - two_center_bundle_.tabulate(inp.lcao_ecut, inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax); + two_center_bundle_.tabulate(inp.lcao_ecut, + inp.lcao_dk, + inp.lcao_dr, + inp.lcao_rmax); #endif - if (this->orb_con.setup_2d) - { - this->orb_con.setup_2d_division(GlobalV::ofs_running, GlobalV::ofs_warning); - this->orb_con.ParaV.set_atomic_trace(GlobalC::ucell.get_iat2iwt(), GlobalC::ucell.nat, GlobalV::NLOCAL); + if (this->orb_con.setup_2d) { + this->orb_con.setup_2d_division(GlobalV::ofs_running, + GlobalV::ofs_warning); + this->orb_con.ParaV.set_atomic_trace(GlobalC::ucell.get_iat2iwt(), + GlobalC::ucell.nat, + GlobalV::NLOCAL); } return; @@ -567,45 +594,39 @@ void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, Input& inp, //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::iter_init(const int istep, const int iter) -{ +void ESolver_KS_LCAO::iter_init(const int istep, const int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_init"); - if (iter == 1) - { + if (iter == 1) { this->p_chgmix->init_mixing(); // init mixing this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1; this->p_chgmix->mixing_restart_count = 0; // this output will be removed once the feeature is stable - if (GlobalC::dftu.uramping > 0.01) - { + if (GlobalC::dftu.uramping > 0.01) { std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < GlobalC::dftu.U0.size(); i++) - { + for (int i = 0; i < GlobalC::dftu.U0.size(); i++) { std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; } std::cout << " eV " << std::endl; } } // for mixing restart - if (iter == this->p_chgmix->mixing_restart_step && GlobalV::MIXING_RESTART > 0.0) - { + if (iter == this->p_chgmix->mixing_restart_step + && GlobalV::MIXING_RESTART > 0.0) { this->p_chgmix->init_mixing(); this->p_chgmix->mixing_restart_count++; - if (GlobalV::dft_plus_u) - { - GlobalC::dftu.uramping_update(); // update U by uramping if uramping > 0.01 - if (GlobalC::dftu.uramping > 0.01) - { + if (GlobalV::dft_plus_u) { + GlobalC::dftu + .uramping_update(); // update U by uramping if uramping > 0.01 + if (GlobalC::dftu.uramping > 0.01) { std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < GlobalC::dftu.U0.size(); i++) - { - std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; + for (int i = 0; i < GlobalC::dftu.U0.size(); i++) { + std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV + << " "; } std::cout << " eV " << std::endl; } - if (GlobalC::dftu.uramping > 0.01 && !GlobalC::dftu.u_converged()) - { + if (GlobalC::dftu.uramping > 0.01 && !GlobalC::dftu.u_converged()) { this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1; } } @@ -613,7 +634,8 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) { // allocate memory for dmr_mdata const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); + = dynamic_cast*>(this->pelec) + ->get_DM(); int nnr_tmp = dm->get_DMR_pointer(1)->get_nnr(); this->p_chgmix->allocate_mixing_dmr(nnr_tmp); } @@ -625,10 +647,8 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) // mohan move it outside 2011-01-13 // first need to calculate the weight according to // electrons number. - if (istep == 0 && this->wf.init_wfc == "file") - { - if (iter == 1) - { + if (istep == 0 && this->wf.init_wfc == "file") { + if (iter == 1) { std::cout << " WAVEFUN -> CHARGE " << std::endl; // calculate the density matrix using read in wave functions @@ -660,13 +680,13 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) // rho1 and rho2 are the same rho. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if (GlobalV::NSPIN == 4) - { + if (GlobalV::NSPIN == 4) { GlobalC::ucell.cal_ux(); } //! update the potentials by using new electron charge density - this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->pelec->charge, + &GlobalC::ucell); //! compute the correction energy for metals this->pelec->f_en.descf = this->pelec->cal_delta_escf(); @@ -675,24 +695,28 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) #ifdef __EXX // calculate exact-exchange - if (GlobalC::exx_info.info_ri.real_number) - { - this->exd->exx_eachiterinit(*dynamic_cast*>(this->pelec)->get_DM(), iter); - } - else - { - this->exc->exx_eachiterinit(*dynamic_cast*>(this->pelec)->get_DM(), iter); + if (GlobalC::exx_info.info_ri.real_number) { + this->exd->exx_eachiterinit( + *dynamic_cast*>(this->pelec) + ->get_DM(), + iter); + } else { + this->exc->exx_eachiterinit( + *dynamic_cast*>(this->pelec) + ->get_DM(), + iter); } #endif - if (GlobalV::dft_plus_u) - { - if (istep != 0 || iter != 1) - { - GlobalC::dftu.set_dmr(dynamic_cast*>(this->pelec)->get_DM()); + if (GlobalV::dft_plus_u) { + if (istep != 0 || iter != 1) { + GlobalC::dftu.set_dmr( + dynamic_cast*>(this->pelec) + ->get_DM()); } // Calculate U and J if Yukawa potential is used - GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx); + GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, + this->pw_rho->nrxx); } #ifdef __DEEPKS @@ -700,27 +724,25 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) GlobalC::ld.set_hr_cal(true); // HR in HamiltLCAO should be recalculate - if (GlobalV::deepks_scf) - { + if (GlobalV::deepks_scf) { this->p_hamilt->refresh(); } #endif - if (GlobalV::VL_IN_H) - { + if (GlobalV::VL_IN_H) { // update Gint_K - if (!GlobalV::GAMMA_ONLY_LOCAL) - { + if (!GlobalV::GAMMA_ONLY_LOCAL) { this->GK.renew(); } // update real space Hamiltonian this->p_hamilt->refresh(); } - // run the inner lambda loop to contrain atomic moments with the DeltaSpin method - if (GlobalV::sc_mag_switch && iter > GlobalV::sc_scf_nmin) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + // run the inner lambda loop to contrain atomic moments with the DeltaSpin + // method + if (GlobalV::sc_mag_switch && iter > GlobalV::sc_scf_nmin) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.run_lambda_loop(iter - 1); } } @@ -742,75 +764,69 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) //! 12) calculate delta energy //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) -{ +void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) { ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density"); // 1) save input rho this->pelec->charge->save_rho_before_sum_band(); // 2) save density matrix DMR for mixing - if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR && this->p_chgmix->mixing_restart_count > 0) - { - elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); + if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR + && this->p_chgmix->mixing_restart_count > 0) { + elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec) + ->get_DM(); dm->save_DMR(); } // 3) solve the Hamiltonian and output band gap - if (this->phsol != nullptr) - { + if (this->phsol != nullptr) { // reset energy this->pelec->f_en.eband = 0.0; this->pelec->f_en.demet = 0.0; - this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec, GlobalV::KS_SOLVER); + this->phsol->solve(this->p_hamilt, + this->psi[0], + this->pelec, + GlobalV::KS_SOLVER); - if (GlobalV::out_bandgap) - { - if (!GlobalV::TWO_EFERMI) - { + if (GlobalV::out_bandgap) { + if (!GlobalV::TWO_EFERMI) { this->pelec->cal_bandgap(); - } - else - { + } else { this->pelec->cal_bandgap_updw(); } } - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!"); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", + "HSolver has not been initialed!"); } // 4) print bands for each k-point and each band - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { + for (int ik = 0; ik < this->kv.get_nks(); ++ik) { this->pelec->print_band(ik, INPUT.printe, iter); } // 5) what's the exd used for? #ifdef __EXX - if (GlobalC::exx_info.info_ri.real_number) - { + if (GlobalC::exx_info.info_ri.real_number) { this->exd->exx_hamilt2density(*this->pelec, this->orb_con.ParaV, iter); - } - else - { + } else { this->exc->exx_hamilt2density(*this->pelec, this->orb_con.ParaV, iter); } #endif - // 6) calculate the local occupation number matrix and energy correction in DFT+U - if (GlobalV::dft_plus_u) - { + // 6) calculate the local occupation number matrix and energy correction in + // DFT+U + if (GlobalV::dft_plus_u) { // only old DFT+U method should calculated energy correction in esolver, // new DFT+U method will calculate energy in calculating Hamiltonian - if (GlobalV::dft_plus_u == 2) - { - if (GlobalC::dftu.omc != 2) - { + if (GlobalV::dft_plus_u == 2) { + if (GlobalC::dftu.omc != 2) { const std::vector>& tmp_dm - = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); + = dynamic_cast*>(this->pelec) + ->get_DM() + ->get_DMK_vector(); this->dftu_cal_occup_m(iter, tmp_dm); } GlobalC::dftu.cal_energy_correction(istep); @@ -820,19 +836,20 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) // (7) for deepks, calculate delta_e #ifdef __DEEPKS - if (GlobalV::deepks_scf) - { + if (GlobalV::deepks_scf) { const std::vector>& dm - = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); + = dynamic_cast*>(this->pelec) + ->get_DM() + ->get_DMK_vector(); this->dpks_cal_e_delta_band(dm); } #endif // 8) for delta spin - if (GlobalV::sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + if (GlobalV::sc_mag_switch) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.cal_MW(iter, &(this->LM)); } @@ -841,9 +858,12 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) // 10) symmetrize the charge density Symmetry_rho srho; - for (int is = 0; is < GlobalV::NSPIN; is++) - { - srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm); + for (int is = 0; is < GlobalV::NSPIN; is++) { + srho.begin(is, + *(this->pelec->charge), + this->pw_rho, + GlobalC::Pgrid, + GlobalC::ucell.symm); } // 11) compute magnetization, only for spin==2 @@ -864,31 +884,26 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) //! 3) print potential //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::update_pot(const int istep, const int iter) -{ +void ESolver_KS_LCAO::update_pot(const int istep, const int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "update_pot"); // 1) print Hamiltonian and Overlap matrix - if (this->conv_elec || iter == GlobalV::SCF_NMAX) - { - if (!GlobalV::GAMMA_ONLY_LOCAL && hsolver::HSolverLCAO::out_mat_hs[0]) - { + if (this->conv_elec || iter == GlobalV::SCF_NMAX) { + if (!GlobalV::GAMMA_ONLY_LOCAL + && hsolver::HSolverLCAO::out_mat_hs[0]) { this->GK.renew(true); } - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { - if (hsolver::HSolverLCAO::out_mat_hs[0]) - { + for (int ik = 0; ik < this->kv.get_nks(); ++ik) { + if (hsolver::HSolverLCAO::out_mat_hs[0]) { this->p_hamilt->updateHk(ik); } bool bit = false; // LiuXh, 2017-03-21 - // if set bit = true, there would be error in soc-multi-core calculation, noted by zhengdy-soc - if (this->psi != nullptr && (istep % GlobalV::out_interval == 0)) - { + // if set bit = true, there would be error in soc-multi-core + // calculation, noted by zhengdy-soc + if (this->psi != nullptr && (istep % GlobalV::out_interval == 0)) { hamilt::MatrixBlock h_mat, s_mat; this->p_hamilt->matrix(h_mat, s_mat); - if (hsolver::HSolverLCAO::out_mat_hs[0]) - { + if (hsolver::HSolverLCAO::out_mat_hs[0]) { ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, @@ -917,9 +932,9 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) } // 2) print wavefunctions - if (elecstate::ElecStateLCAO::out_wfc_lcao && (this->conv_elec || iter == GlobalV::SCF_NMAX) - && (istep % GlobalV::out_interval == 0)) - { + if (elecstate::ElecStateLCAO::out_wfc_lcao + && (this->conv_elec || iter == GlobalV::SCF_NMAX) + && (istep % GlobalV::out_interval == 0)) { ModuleIO::write_wfc_nao(elecstate::ElecStateLCAO::out_wfc_lcao, this->psi[0], this->pelec->ekb, @@ -930,25 +945,21 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) } // 3) print potential - if (this->conv_elec || iter == GlobalV::SCF_NMAX) - { + if (this->conv_elec || iter == GlobalV::SCF_NMAX) { if (GlobalV::out_pot < 0) // mohan add 2011-10-10 { GlobalV::out_pot = -2; } } - if (!this->conv_elec) - { - if (GlobalV::NSPIN == 4) - { + if (!this->conv_elec) { + if (GlobalV::NSPIN == 4) { GlobalC::ucell.cal_ux(); } - this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->pelec->charge, + &GlobalC::ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); - } - else - { + } else { this->pelec->cal_converged(); } } @@ -964,47 +975,58 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) //! 6) calculate the total energy? //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::iter_finish(int iter) -{ +void ESolver_KS_LCAO::iter_finish(int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish"); - // 1) mix density matrix if mixing_restart + mixing_dmr + not first mixing_restart at every iter - if (GlobalV::MIXING_RESTART > 0 && this->p_chgmix->mixing_restart_count > 0 && GlobalV::MIXING_DMR) - { - elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); + // 1) mix density matrix if mixing_restart + mixing_dmr + not first + // mixing_restart at every iter + if (GlobalV::MIXING_RESTART > 0 && this->p_chgmix->mixing_restart_count > 0 + && GlobalV::MIXING_DMR) { + elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec) + ->get_DM(); this->p_chgmix->mix_dmr(dm); } // 2) save charge density // Peize Lin add 2020.04.04 - if (GlobalC::restart.info_save.save_charge) - { - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - GlobalC::restart.save_disk("charge", is, this->pelec->charge->nrxx, this->pelec->charge->rho[is]); + if (GlobalC::restart.info_save.save_charge) { + for (int is = 0; is < GlobalV::NSPIN; ++is) { + GlobalC::restart.save_disk("charge", + is, + this->pelec->charge->nrxx, + this->pelec->charge->rho[is]); } } #ifdef __EXX // 3) save exx matrix - int two_level_step = GlobalC::exx_info.info_ri.real_number ? this->exd->two_level_step : this->exc->two_level_step; + int two_level_step = GlobalC::exx_info.info_ri.real_number + ? this->exd->two_level_step + : this->exc->two_level_step; if (GlobalC::restart.info_save.save_H && two_level_step > 0 - && (!GlobalC::exx_info.info_global.separate_loop || iter == 1)) // to avoid saving the same value repeatedly + && (!GlobalC::exx_info.info_global.separate_loop + || iter == 1)) // to avoid saving the same value repeatedly { std::vector Hexxk_save(this->orb_con.ParaV.get_local_size()); - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { + for (int ik = 0; ik < this->kv.get_nks(); ++ik) { ModuleBase::GlobalFunc::ZEROS(Hexxk_save.data(), Hexxk_save.size()); - hamilt::OperatorEXX> opexx_save(&this->LM, nullptr, &Hexxk_save, this->kv); + hamilt::OperatorEXX> opexx_save( + &this->LM, + nullptr, + &Hexxk_save, + this->kv); opexx_save.contributeHk(ik); - GlobalC::restart.save_disk("Hexx", ik, this->orb_con.ParaV.get_local_size(), Hexxk_save.data()); + GlobalC::restart.save_disk("Hexx", + ik, + this->orb_con.ParaV.get_local_size(), + Hexxk_save.data()); } - if (GlobalV::MY_RANK == 0) - { + if (GlobalV::MY_RANK == 0) { GlobalC::restart.save_disk("Eexx", 0, 1, &this->pelec->f_en.exx); } } @@ -1012,28 +1034,26 @@ void ESolver_KS_LCAO::iter_finish(int iter) // 4) output charge density and density matrix bool print = false; - if (this->out_freq_elec && iter % this->out_freq_elec == 0) - { + if (this->out_freq_elec && iter % this->out_freq_elec == 0) { print = true; } - if (print) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { + if (print) { + for (int is = 0; is < GlobalV::NSPIN; is++) { this->create_Output_Rho(is, iter, "tmp_").write(); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { + if (XC_Functional::get_func_type() == 3 + || XC_Functional::get_func_type() == 5) { this->create_Output_Kin(is, iter, "tmp_").write(); } } } // 5) cal_MW? - // escon: energy of spin constraint depends on Mi, so cal_energies should be called after cal_MW - if (GlobalV::sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + // escon: energy of spin constraint depends on Mi, so cal_energies should be + // called after cal_MW + if (GlobalV::sc_mag_switch) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.cal_MW(iter, &(this->LM)); } @@ -1064,13 +1084,11 @@ void ESolver_KS_LCAO::iter_finish(int iter) //! 18) write quasi-orbitals //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::after_scf(const int istep) -{ +void ESolver_KS_LCAO::after_scf(const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf"); // 1) write charge difference into files for charge extrapolation - if (GlobalV::CALCULATION != "scf") - { + if (GlobalV::CALCULATION != "scf") { this->CE.save_files(istep, GlobalC::ucell, #ifdef __MPI @@ -1081,53 +1099,53 @@ void ESolver_KS_LCAO::after_scf(const int istep) } // 2) write density matrix for sparse matrix - ModuleIO::write_dmr(dynamic_cast*>(this->pelec)->get_DM()->get_DMR_vector(), - this->orb_con.ParaV, - INPUT.out_dm1, - false, - GlobalV::out_app_flag, - istep); + ModuleIO::write_dmr( + dynamic_cast*>(this->pelec) + ->get_DM() + ->get_DMR_vector(), + this->orb_con.ParaV, + INPUT.out_dm1, + false, + GlobalV::out_app_flag, + istep); // 3) write charge density - if (GlobalV::out_chg) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { + if (GlobalV::out_chg) { + for (int is = 0; is < GlobalV::NSPIN; is++) { this->create_Output_Rho(is, istep).write(); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { + if (XC_Functional::get_func_type() == 3 + || XC_Functional::get_func_type() == 5) { this->create_Output_Kin(is, istep).write(); } } } // 4) write density matrix - if (INPUT.out_dm) - { + if (INPUT.out_dm) { std::vector efermis(GlobalV::NSPIN == 2 ? 2 : 1); - for (int ispin = 0; ispin < efermis.size(); ispin++) - { + for (int ispin = 0; ispin < efermis.size(); ispin++) { efermis[ispin] = this->pelec->eferm.get_efval(ispin); } const int precision = 3; - ModuleIO::write_dmk(dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(), - precision, - efermis, - &(GlobalC::ucell), - this->orb_con.ParaV); + ModuleIO::write_dmk( + dynamic_cast*>(this->pelec) + ->get_DM() + ->get_DMK_vector(), + precision, + efermis, + &(GlobalC::ucell), + this->orb_con.ParaV); } #ifdef __EXX // 5) write Exx matrix if (GlobalC::exx_info.info_global.cal_exx) // Peize Lin add if 2022.11.14 { - const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); - if (GlobalC::exx_info.info_ri.real_number) - { + const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + + std::to_string(GlobalV::MY_RANK); + if (GlobalC::exx_info.info_ri.real_number) { this->exd->write_Hexxs_csr(file_name_exx, GlobalC::ucell); - } - else - { + } else { this->exc->write_Hexxs_csr(file_name_exx, GlobalC::ucell); } } @@ -1137,116 +1155,124 @@ void ESolver_KS_LCAO::after_scf(const int istep) this->create_Output_Potential(istep).write(); // 7) write convergence - ModuleIO::output_convergence_after_scf(this->conv_elec, this->pelec->f_en.etot); + ModuleIO::output_convergence_after_scf(this->conv_elec, + this->pelec->f_en.etot); // 8) write fermi energy ModuleIO::output_efermi(this->conv_elec, this->pelec->eferm.ef); // 9) write eigenvalues - if (GlobalV::OUT_LEVEL != "m") - { + if (GlobalV::OUT_LEVEL != "m") { this->pelec->print_eigenvalue(GlobalV::ofs_running); } // 10) write deepks information #ifdef __DEEPKS - std::shared_ptr ld_shared_ptr(&GlobalC::ld, [](LCAO_Deepks*) {}); + std::shared_ptr ld_shared_ptr(&GlobalC::ld, + [](LCAO_Deepks*) {}); LCAO_Deepks_Interface LDI = LCAO_Deepks_Interface(ld_shared_ptr); ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels"); - LDI.out_deepks_labels(this->pelec->f_en.etot, - this->pelec->klist->get_nks(), - GlobalC::ucell.nat, - this->pelec->ekb, - this->pelec->klist->kvec_d, - GlobalC::ucell, - GlobalC::ORB, - GlobalC::GridD, - &(this->orb_con.ParaV), - *(this->psi), - dynamic_cast*>(this->pelec)->get_DM()); + LDI.out_deepks_labels( + this->pelec->f_en.etot, + this->pelec->klist->get_nks(), + GlobalC::ucell.nat, + this->pelec->ekb, + this->pelec->klist->kvec_d, + GlobalC::ucell, + GlobalC::ORB, + GlobalC::GridD, + &(this->orb_con.ParaV), + *(this->psi), + dynamic_cast*>(this->pelec) + ->get_DM()); ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels"); #endif #ifdef __EXX // 11) write rpa information - if (INPUT.rpa) - { - // ModuleRPA::DFT_RPA_interface rpa_interface(GlobalC::exx_info.info_global); + if (INPUT.rpa) { + // ModuleRPA::DFT_RPA_interface + // rpa_interface(GlobalC::exx_info.info_global); // rpa_interface.rpa_exx_lcao().info.files_abfs = GlobalV::rpa_orbitals; RPA_LRI rpa_lri_double(GlobalC::exx_info.info_ri); - rpa_lri_double.cal_postSCF_exx(*dynamic_cast*>(this->pelec)->get_DM(), - MPI_COMM_WORLD, - this->kv); + rpa_lri_double.cal_postSCF_exx( + *dynamic_cast*>(this->pelec) + ->get_DM(), + MPI_COMM_WORLD, + this->kv); rpa_lri_double.init(MPI_COMM_WORLD, this->kv); - rpa_lri_double.out_for_RPA(this->orb_con.ParaV, *(this->psi), this->pelec); + rpa_lri_double.out_for_RPA(this->orb_con.ParaV, + *(this->psi), + this->pelec); } #endif // 12) write HR in npz format - if (GlobalV::out_hr_npz) - { + if (GlobalV::out_hr_npz) { this->p_hamilt->updateHk(0); // first k point, up spin hamilt::HamiltLCAO, double>* p_ham_lcao - = dynamic_cast, double>*>(this->p_hamilt); + = dynamic_cast, double>*>( + this->p_hamilt); std::string zipname = "output_HR0.npz"; this->output_mat_npz(zipname, *(p_ham_lcao->getHR())); - if (GlobalV::NSPIN == 2) - { - this->p_hamilt->updateHk(this->kv.get_nks() / 2); // the other half of k points, down spin + if (GlobalV::NSPIN == 2) { + this->p_hamilt->updateHk( + this->kv.get_nks() + / 2); // the other half of k points, down spin hamilt::HamiltLCAO, double>* p_ham_lcao - = dynamic_cast, double>*>(this->p_hamilt); + = dynamic_cast< + hamilt::HamiltLCAO, double>*>( + this->p_hamilt); zipname = "output_HR1.npz"; this->output_mat_npz(zipname, *(p_ham_lcao->getHR())); } } // 13) write dm in npz format - if (GlobalV::out_dm_npz) - { + if (GlobalV::out_dm_npz) { const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); + = dynamic_cast*>(this->pelec) + ->get_DM(); std::string zipname = "output_DM0.npz"; this->output_mat_npz(zipname, *(dm->get_DMR_pointer(1))); - if (GlobalV::NSPIN == 2) - { + if (GlobalV::NSPIN == 2) { zipname = "output_DM1.npz"; this->output_mat_npz(zipname, *(dm->get_DMR_pointer(2))); } } // 14) write md related - if (!md_skip_out(GlobalV::CALCULATION, istep, GlobalV::out_interval)) - { + if (!md_skip_out(GlobalV::CALCULATION, istep, GlobalV::out_interval)) { this->create_Output_Mat_Sparse(istep).write(); // mulliken charge analysis - if (GlobalV::out_mul) - { + if (GlobalV::out_mul) { this->cal_mag(istep, true); } } // 15) write spin constrian MW? // spin constrain calculations, added by Tianqi Zhao. - if (GlobalV::sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + if (GlobalV::sc_mag_switch) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.cal_MW(istep, &(this->LM), true); sc.print_Mag_Force(); } // 16) delete grid - if (!GlobalV::CAL_FORCE && !GlobalV::CAL_STRESS) - { + if (!GlobalV::CAL_FORCE && !GlobalV::CAL_STRESS) { RA.delete_grid(); } // 17) write quasi-orbitals, added by Yike Huang. - if (GlobalV::qo_switch) - { - toQO tqo(GlobalV::qo_basis, GlobalV::qo_strategy, GlobalV::qo_thr, GlobalV::qo_screening_coeff); + if (GlobalV::qo_switch) { + toQO tqo(GlobalV::qo_basis, + GlobalV::qo_strategy, + GlobalV::qo_thr, + GlobalV::qo_screening_coeff); tqo.initialize(GlobalV::global_out_dir, GlobalV::global_pseudo_dir, GlobalV::global_orbital_dir, @@ -1264,31 +1290,30 @@ void ESolver_KS_LCAO::after_scf(const int istep) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -bool ESolver_KS_LCAO::do_after_converge(int& iter) -{ +bool ESolver_KS_LCAO::do_after_converge(int& iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "do_after_converge"); #ifdef __EXX - if (GlobalC::exx_info.info_ri.real_number) - { - return this->exd->exx_after_converge(*this->p_hamilt, - this->LM, - *dynamic_cast*>(this->pelec)->get_DM(), - this->kv, - iter); - } - else - { - return this->exc->exx_after_converge(*this->p_hamilt, - this->LM, - *dynamic_cast*>(this->pelec)->get_DM(), - this->kv, - iter); + if (GlobalC::exx_info.info_ri.real_number) { + return this->exd->exx_after_converge( + *this->p_hamilt, + this->LM, + *dynamic_cast*>(this->pelec) + ->get_DM(), + this->kv, + iter); + } else { + return this->exc->exx_after_converge( + *this->p_hamilt, + this->LM, + *dynamic_cast*>(this->pelec) + ->get_DM(), + this->kv, + iter); } #endif // __EXX - if (GlobalV::dft_plus_u) - { + if (GlobalV::dft_plus_u) { // use the converged occupation matrix for next MD/Relax SCF calculation GlobalC::dftu.initialed_locale = true; } @@ -1311,21 +1336,22 @@ bool ESolver_KS_LCAO::do_after_converge(int& iter) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Sparse(int istep) -{ - return ModuleIO::Output_Mat_Sparse(hsolver::HSolverLCAO::out_mat_hsR, - hsolver::HSolverLCAO::out_mat_dh, - hsolver::HSolverLCAO::out_mat_t, - INPUT.out_mat_r, - istep, - this->pelec->pot->get_effective_v(), - this->orb_con.ParaV, - this->GK, // mohan add 2024-04-01 - two_center_bundle_, - this->LM, - GlobalC::GridD, // mohan add 2024-04-06 - this->kv, - this->p_hamilt); +ModuleIO::Output_Mat_Sparse + ESolver_KS_LCAO::create_Output_Mat_Sparse(int istep) { + return ModuleIO::Output_Mat_Sparse( + hsolver::HSolverLCAO::out_mat_hsR, + hsolver::HSolverLCAO::out_mat_dh, + hsolver::HSolverLCAO::out_mat_t, + INPUT.out_mat_r, + istep, + this->pelec->pot->get_effective_v(), + this->orb_con.ParaV, + this->GK, // mohan add 2024-04-01 + two_center_bundle_, + this->LM, + GlobalC::GridD, // mohan add 2024-04-06 + this->kv, + this->p_hamilt); } //------------------------------------------------------------------------------ @@ -1333,12 +1359,11 @@ ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Spars //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -bool ESolver_KS_LCAO::md_skip_out(std::string calculation, int istep, int interval) -{ - if (calculation == "md") - { - if (istep % interval != 0) - { +bool ESolver_KS_LCAO::md_skip_out(std::string calculation, + int istep, + int interval) { + if (calculation == "md") { + if (istep % interval != 0) { return true; } } @@ -1346,8 +1371,7 @@ bool ESolver_KS_LCAO::md_skip_out(std::string calculation, int istep, in } template -void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) -{ +void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) { auto cell_index = CellIndex(GlobalC::ucell.get_atomLabels(), GlobalC::ucell.get_atomCounts(), GlobalC::ucell.get_lnchiCounts(), @@ -1357,10 +1381,12 @@ void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) &(this->orb_con.ParaV), GlobalV::NSPIN, this->kv.get_nks()); - auto out_dmk = ModuleIO::Output_DMK(dynamic_cast*>(this->pelec)->get_DM(), - &(this->orb_con.ParaV), - GlobalV::NSPIN, - this->kv.get_nks()); + auto out_dmk = ModuleIO::Output_DMK( + dynamic_cast*>(this->pelec) + ->get_DM(), + &(this->orb_con.ParaV), + GlobalV::NSPIN, + this->kv.get_nks()); auto mulp = ModuleIO::Output_Mulliken(&(out_sk), &(out_dmk), &(this->orb_con.ParaV), @@ -1370,8 +1396,7 @@ void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) auto atom_chg = mulp.get_atom_chg(); /// used in updating mag info in STRU file GlobalC::ucell.atom_mulliken = mulp.get_atom_mulliken(atom_chg); - if (print && GlobalV::MY_RANK == 0) - { + if (print && GlobalV::MY_RANK == 0) { /// write the Orbital file cell_index.write_orb_info(GlobalV::global_out_dir); /// write mulliken.txt diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index f362d28cd4..d7aa653f3f 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -16,11 +16,9 @@ #include -namespace ModuleESolver -{ +namespace ModuleESolver { template -class ESolver_KS_LCAO : public ESolver_KS -{ +class ESolver_KS_LCAO : public ESolver_KS { public: ESolver_KS_LCAO(); ~ESolver_KS_LCAO(); @@ -48,7 +46,9 @@ class ESolver_KS_LCAO : public ESolver_KS virtual void iter_init(const int istep, const int iter) override; - virtual void hamilt2density(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density(const int istep, + const int iter, + const double ethr) override; virtual void update_pot(const int istep, const int iter) override; @@ -98,11 +98,13 @@ class ESolver_KS_LCAO : public ESolver_KS void beforesolver(const int istep); //---------------------------------------------------------------------- - /// @brief create ModuleIO::Output_Mat_Sparse object to output sparse density matrix of H, S, T, r + /// @brief create ModuleIO::Output_Mat_Sparse object to output sparse + /// density matrix of H, S, T, r ModuleIO::Output_Mat_Sparse create_Output_Mat_Sparse(int istep); void read_mat_npz(std::string& zipname, hamilt::HContainer& hR); - void output_mat_npz(std::string& zipname, const hamilt::HContainer& hR); + void output_mat_npz(std::string& zipname, + const hamilt::HContainer& hR); /// @brief check if skip the corresponding output in md calculation bool md_skip_out(std::string calculation, int istep, int interval); @@ -116,12 +118,14 @@ class ESolver_KS_LCAO : public ESolver_KS private: // tmp interfaces before sub-modules are refactored - void dftu_cal_occup_m(const int& iter, const std::vector>& dm) const; + void dftu_cal_occup_m(const int& iter, + const std::vector>& dm) const; #ifdef __DEEPKS void dpks_cal_e_delta_band(const std::vector>& dm) const; - void dpks_cal_projected_DM(const elecstate::DensityMatrix* dm) const; + void dpks_cal_projected_DM( + const elecstate::DensityMatrix* dm) const; #endif }; } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index bd0a510e3c..ff54dbfd7d 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -27,21 +27,20 @@ #include "module_io/rho_io.h" #include "module_io/write_pot.h" -namespace ModuleESolver -{ +namespace ModuleESolver { template -void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) -{ +void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) { ModuleBase::TITLE("ESolver_KS_LCAO", "set_matrix_grid"); ModuleBase::timer::tick("ESolver_KS_LCAO", "set_matrix_grid"); // (1) Find adjacent atoms for each atom. - GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, - GlobalV::OUT_LEVEL, - GlobalC::ORB.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), - GlobalV::GAMMA_ONLY_LOCAL); + GlobalV::SEARCH_RADIUS + = atom_arrange::set_sr_NL(GlobalV::ofs_running, + GlobalV::OUT_LEVEL, + GlobalC::ORB.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + GlobalV::GAMMA_ONLY_LOCAL); atom_arrange::search(GlobalV::SEARCH_PBC, GlobalV::ofs_running, @@ -50,7 +49,8 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) GlobalV::SEARCH_RADIUS, GlobalV::test_atom_input); - // ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"SEARCH ADJACENT ATOMS"); + // ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"SEARCH ADJACENT + // ATOMS"); // (3) Periodic condition search for each grid. double dr_uniform = 0.001; @@ -59,7 +59,12 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) std::vector> dpsi_u; std::vector> d2psi_u; - Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, psi_u, dpsi_u, d2psi_u); + Gint_Tools::init_orb(dr_uniform, + rcuts, + GlobalC::ucell, + psi_u, + dpsi_u, + d2psi_u); this->GridT.set_pbc_grid(this->pw_rho->nx, this->pw_rho->ny, @@ -94,8 +99,7 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) // If k point is used here, allocate HlocR after atom_arrange. Parallel_Orbitals* pv = this->LM.ParaV; ra.for_2d(*pv, GlobalV::GAMMA_ONLY_LOCAL); - if (!GlobalV::GAMMA_ONLY_LOCAL) - { + if (!GlobalV::GAMMA_ONLY_LOCAL) { // need to first calculae lgd. // using GridT.init. this->GridT.cal_nnrg(pv); @@ -106,8 +110,7 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) } template -void ESolver_KS_LCAO::beforesolver(const int istep) -{ +void ESolver_KS_LCAO::beforesolver(const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "beforesolver"); ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); @@ -122,22 +125,20 @@ void ESolver_KS_LCAO::beforesolver(const int istep) // the force. // init psi - if (this->psi == nullptr) - { + if (this->psi == nullptr) { int nsk = 0; int ncol = 0; - if (GlobalV::GAMMA_ONLY_LOCAL) - { + if (GlobalV::GAMMA_ONLY_LOCAL) { nsk = GlobalV::NSPIN; ncol = this->orb_con.ParaV.ncol_bands; - if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "lapack_gvx" || GlobalV::KS_SOLVER == "pexsi" - || GlobalV::KS_SOLVER == "cusolver" || GlobalV::KS_SOLVER == "cusolvermp") - { + if (GlobalV::KS_SOLVER == "genelpa" + || GlobalV::KS_SOLVER == "lapack_gvx" + || GlobalV::KS_SOLVER == "pexsi" + || GlobalV::KS_SOLVER == "cusolver" + || GlobalV::KS_SOLVER == "cusolvermp") { ncol = this->orb_con.ParaV.ncol; } - } - else - { + } else { nsk = this->kv.get_nks(); #ifdef __MPI ncol = this->orb_con.ParaV.ncol_bands; @@ -145,21 +146,26 @@ void ESolver_KS_LCAO::beforesolver(const int istep) ncol = GlobalV::NBANDS; #endif } - this->psi = new psi::Psi(nsk, ncol, this->orb_con.ParaV.nrow, nullptr); + this->psi + = new psi::Psi(nsk, ncol, this->orb_con.ParaV.nrow, nullptr); } // prepare grid in Gint - LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, *this->pw_rho, *this->pw_big); + LCAO_domain::grid_prepare(this->GridT, + this->GG, + this->GK, + *this->pw_rho, + *this->pw_big); // init Hamiltonian - if (this->p_hamilt != nullptr) - { + if (this->p_hamilt != nullptr) { delete this->p_hamilt; this->p_hamilt = nullptr; } - if (this->p_hamilt == nullptr) - { - elecstate::DensityMatrix* DM = dynamic_cast*>(this->pelec)->get_DM(); + if (this->p_hamilt == nullptr) { + elecstate::DensityMatrix* DM + = dynamic_cast*>(this->pelec) + ->get_DM(); this->p_hamilt = new hamilt::HamiltLCAO( GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr, GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK), @@ -170,19 +176,23 @@ void ESolver_KS_LCAO::beforesolver(const int istep) two_center_bundle_, #ifdef __EXX DM, - GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step); + GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step + : &this->exc->two_level_step); #else DM); #endif } // init density kernel and wave functions. - this->LOC.allocate_dm_wfc(this->GridT, this->pelec, this->psi, this->kv, istep); + this->LOC.allocate_dm_wfc(this->GridT, + this->pelec, + this->psi, + this->kv, + istep); #ifdef __DEEPKS // for each ionic step, the overlap must be rebuilt // since it depends on ionic positions - if (GlobalV::deepks_setorb) - { + if (GlobalV::deepks_setorb) { const Parallel_Orbitals* pv = this->LM.ParaV; // build and save at beginning GlobalC::ld.build_psialpha(GlobalV::CAL_FORCE, @@ -191,15 +201,17 @@ void ESolver_KS_LCAO::beforesolver(const int istep) GlobalC::GridD, *(two_center_bundle_.overlap_orb_alpha)); - if (GlobalV::deepks_out_unittest) - { - GlobalC::ld.check_psialpha(GlobalV::CAL_FORCE, GlobalC::ucell, GlobalC::ORB, GlobalC::GridD); + if (GlobalV::deepks_out_unittest) { + GlobalC::ld.check_psialpha(GlobalV::CAL_FORCE, + GlobalC::ucell, + GlobalC::ORB, + GlobalC::GridD); } } #endif - if (GlobalV::sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + if (GlobalV::sc_mag_switch) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.init_sc(GlobalV::sc_thr, GlobalV::nsc, GlobalV::nsc_min, @@ -223,25 +235,21 @@ void ESolver_KS_LCAO::beforesolver(const int istep) // cal_ux should be called before init_scf because // the direction of ux is used in noncoline_rho //========================================================= - if (GlobalV::NSPIN == 4 && GlobalV::DOMAG) - { + if (GlobalV::NSPIN == 4 && GlobalV::DOMAG) { GlobalC::ucell.cal_ux(); } ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); } template -void ESolver_KS_LCAO::before_scf(const int istep) -{ +void ESolver_KS_LCAO::before_scf(const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf"); ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); - if (GlobalC::ucell.cell_parameter_updated) - { + if (GlobalC::ucell.cell_parameter_updated) { this->init_after_vc(INPUT, GlobalC::ucell); } - if (GlobalC::ucell.ionic_position_updated) - { + if (GlobalC::ucell.ionic_position_updated) { this->CE.update_all_dis(GlobalC::ucell); this->CE.extrapolate_charge( #ifdef __MPI @@ -256,8 +264,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) // about vdw, jiyy add vdwd3 and linpz add vdwd2 //---------------------------------------------------------- auto vdw_solver = vdw::make_vdw(GlobalC::ucell, INPUT); - if (vdw_solver != nullptr) - { + if (vdw_solver != nullptr) { this->pelec->f_en.evdw = vdw_solver->get_energy(); } @@ -265,23 +272,19 @@ void ESolver_KS_LCAO::before_scf(const int istep) // Peize Lin add 2016-12-03 #ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf - if (GlobalC::exx_info.info_ri.real_number) - { + if (GlobalC::exx_info.info_ri.real_number) { this->exd->exx_beforescf(this->kv, *this->p_chgmix); - } - else - { + } else { this->exc->exx_beforescf(this->kv, *this->p_chgmix); } #endif // __EXX this->pelec->init_scf(istep, this->sf.strucFac); - if (GlobalV::out_chg == 2) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { + if (GlobalV::out_chg == 2) { + for (int is = 0; is < GlobalV::NSPIN; is++) { std::stringstream ss; - ss << GlobalV::global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube"; + ss << GlobalV::global_out_dir << "SPIN" << is + 1 + << "_CHG_INI.cube"; ModuleIO::write_rho( #ifdef __MPI this->pw_big->bz, // bz first, then nbz @@ -321,16 +324,16 @@ void ESolver_KS_LCAO::before_scf(const int istep) // DMR should be same size with Hamiltonian(R) dynamic_cast*>(this->pelec) ->get_DM() - ->init_DMR(*(dynamic_cast*>(this->p_hamilt)->getHR())); + ->init_DMR(*(dynamic_cast*>(this->p_hamilt) + ->getHR())); - if (GlobalV::dm_to_rho) - { + if (GlobalV::dm_to_rho) { std::string zipname = "output_DM0.npz"; elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); + = dynamic_cast*>(this->pelec) + ->get_DM(); this->read_mat_npz(zipname, *(dm->get_DMR_pointer(1))); - if (GlobalV::NSPIN == 2) - { + if (GlobalV::NSPIN == 2) { zipname = "output_DM1.npz"; this->read_mat_npz(zipname, *(dm->get_DMR_pointer(2))); } @@ -338,8 +341,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) this->pelec->psiToRho(*this->psi); this->create_Output_Rho(0, istep).write(); - if (GlobalV::NSPIN == 2) - { + if (GlobalV::NSPIN == 2) { this->create_Output_Rho(1, istep).write(); } @@ -349,16 +351,21 @@ void ESolver_KS_LCAO::before_scf(const int istep) // the electron charge density should be symmetrized, // here is the initialization Symmetry_rho srho; - for (int is = 0; is < GlobalV::NSPIN; is++) - { - srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm); + for (int is = 0; is < GlobalV::NSPIN; is++) { + srho.begin(is, + *(this->pelec->charge), + this->pw_rho, + GlobalC::Pgrid, + GlobalC::ucell.symm); } // 1. calculate ewald energy. // mohan update 2021-02-25 - if (!GlobalV::test_skip_ewald) - { - this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(GlobalC::ucell, this->pw_rho, this->sf.strucFac); + if (!GlobalV::test_skip_ewald) { + this->pelec->f_en.ewald_energy + = H_Ewald_pw::compute_ewald(GlobalC::ucell, + this->pw_rho, + this->sf.strucFac); } this->p_hamilt->non_first_scf = istep; @@ -368,15 +375,13 @@ void ESolver_KS_LCAO::before_scf(const int istep) } template -void ESolver_KS_LCAO::others(const int istep) -{ +void ESolver_KS_LCAO::others(const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "others"); ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); const std::string cal_type = GlobalV::CALCULATION; - if (cal_type == "get_S") - { + if (cal_type == "get_S") { std::cout << "\n * * * * * *" << std::endl; std::cout << " << Start writing the overlap matrix." << std::endl; this->get_S(); @@ -385,22 +390,19 @@ void ESolver_KS_LCAO::others(const int istep) ModuleBase::QUIT(); - // return; // use 'return' will cause segmentation fault. by mohan 2024-06-09 - } - else if (cal_type == "test_memory") - { + // return; // use 'return' will cause segmentation fault. by mohan + // 2024-06-09 + } else if (cal_type == "test_memory") { Cal_Test::test_memory(this->pw_rho, this->pw_wfc, this->p_chgmix->get_mixing_mode(), this->p_chgmix->get_mixing_ndim()); return; - } - else if (cal_type == "test_neighbour") - { + } else if (cal_type == "test_neighbour") { // test_search_neighbor(); - if (GlobalV::SEARCH_RADIUS < 0) - { - std::cout << " SEARCH_RADIUS : " << GlobalV::SEARCH_RADIUS << std::endl; + if (GlobalV::SEARCH_RADIUS < 0) { + std::cout << " SEARCH_RADIUS : " << GlobalV::SEARCH_RADIUS + << std::endl; std::cout << " please make sure search_radius > 0" << std::endl; } @@ -418,12 +420,9 @@ void ESolver_KS_LCAO::others(const int istep) // pelec should be initialized before these calculations this->pelec->init_scf(istep, this->sf.strucFac); // self consistent calculations for electronic ground state - if (GlobalV::CALCULATION == "nscf") - { + if (GlobalV::CALCULATION == "nscf") { this->nscf(); - } - else if (cal_type == "get_pchg") - { + } else if (cal_type == "get_pchg") { IState_Charge ISC(this->psi, &(this->orb_con.ParaV)); ISC.begin(this->GG, this->pelec->charge->rho, @@ -449,12 +448,9 @@ void ESolver_KS_LCAO::others(const int istep) GlobalV::ofs_warning, &GlobalC::ucell, &GlobalC::GridD); - } - else if (cal_type == "get_wf") - { + } else if (cal_type == "get_wf") { IState_Envelope IEP(this->pelec); - if (GlobalV::GAMMA_ONLY_LOCAL) - { + if (GlobalV::GAMMA_ONLY_LOCAL) { IEP.begin(this->psi, this->pw_rho, this->pw_wfc, @@ -471,9 +467,7 @@ void ESolver_KS_LCAO::others(const int istep) GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::global_out_dir); - } - else - { + } else { IEP.begin(this->psi, this->pw_rho, this->pw_wfc, @@ -491,10 +485,9 @@ void ESolver_KS_LCAO::others(const int istep) GlobalV::NLOCAL, GlobalV::global_out_dir); } - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", "CALCULATION type not supported"); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", + "CALCULATION type not supported"); } ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); @@ -502,22 +495,22 @@ void ESolver_KS_LCAO::others(const int istep) } template <> -void ESolver_KS_LCAO::get_S(void) -{ +void ESolver_KS_LCAO::get_S(void) { ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::get_S", "not implemented for"); + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::get_S", + "not implemented for"); } template <> -void ESolver_KS_LCAO, double>::get_S(void) -{ +void ESolver_KS_LCAO, double>::get_S(void) { ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); // (1) Find adjacent atoms for each atom. - GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, - GlobalV::OUT_LEVEL, - GlobalC::ORB.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), - GlobalV::GAMMA_ONLY_LOCAL); + GlobalV::SEARCH_RADIUS + = atom_arrange::set_sr_NL(GlobalV::ofs_running, + GlobalV::OUT_LEVEL, + GlobalC::ORB.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + GlobalV::GAMMA_ONLY_LOCAL); atom_arrange::search(GlobalV::SEARCH_PBC, GlobalV::ofs_running, @@ -530,12 +523,14 @@ void ESolver_KS_LCAO, double>::get_S(void) this->LM.ParaV = &this->orb_con.ParaV; - if (this->p_hamilt == nullptr) - { - this->p_hamilt = new hamilt::HamiltLCAO, double>(&this->LM, - this->kv, - *(two_center_bundle_.overlap_orb)); - dynamic_cast, double>*>(this->p_hamilt->ops)->contributeHR(); + if (this->p_hamilt == nullptr) { + this->p_hamilt = new hamilt::HamiltLCAO, double>( + &this->LM, + this->kv, + *(two_center_bundle_.overlap_orb)); + dynamic_cast, double>*>( + this->p_hamilt->ops) + ->contributeHR(); } // mohan add 2024-06-09 @@ -543,21 +538,25 @@ void ESolver_KS_LCAO, double>::get_S(void) std::cout << " The file is saved in " << fn << std::endl; - ModuleIO::output_SR(orb_con.ParaV, this->LM, GlobalC::GridD, this->p_hamilt, fn); + ModuleIO::output_SR(orb_con.ParaV, + this->LM, + GlobalC::GridD, + this->p_hamilt, + fn); return; } template <> -void ESolver_KS_LCAO, std::complex>::get_S(void) -{ +void ESolver_KS_LCAO, std::complex>::get_S(void) { ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); // (1) Find adjacent atoms for each atom. - GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, - GlobalV::OUT_LEVEL, - GlobalC::ORB.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), - GlobalV::GAMMA_ONLY_LOCAL); + GlobalV::SEARCH_RADIUS + = atom_arrange::set_sr_NL(GlobalV::ofs_running, + GlobalV::OUT_LEVEL, + GlobalC::ORB.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + GlobalV::GAMMA_ONLY_LOCAL); atom_arrange::search(GlobalV::SEARCH_PBC, GlobalV::ofs_running, @@ -568,13 +567,15 @@ void ESolver_KS_LCAO, std::complex>::get_S(void) this->RA.for_2d(this->orb_con.ParaV, GlobalV::GAMMA_ONLY_LOCAL); this->LM.ParaV = &this->orb_con.ParaV; - if (this->p_hamilt == nullptr) - { - this->p_hamilt - = new hamilt::HamiltLCAO, std::complex>(&this->LM, - this->kv, - *(two_center_bundle_.overlap_orb)); - dynamic_cast, std::complex>*>(this->p_hamilt->ops) + if (this->p_hamilt == nullptr) { + this->p_hamilt = new hamilt::HamiltLCAO, + std::complex>( + &this->LM, + this->kv, + *(two_center_bundle_.overlap_orb)); + dynamic_cast< + hamilt::OperatorLCAO, std::complex>*>( + this->p_hamilt->ops) ->contributeHR(); } @@ -583,14 +584,17 @@ void ESolver_KS_LCAO, std::complex>::get_S(void) std::cout << " The file is saved in " << fn << std::endl; - ModuleIO::output_SR(orb_con.ParaV, this->LM, GlobalC::GridD, this->p_hamilt, fn); + ModuleIO::output_SR(orb_con.ParaV, + this->LM, + GlobalC::GridD, + this->p_hamilt, + fn); return; } template -void ESolver_KS_LCAO::nscf(void) -{ +void ESolver_KS_LCAO::nscf(void) { ModuleBase::TITLE("ESolver_KS_LCAO", "nscf"); std::cout << " NON-SELF CONSISTENT CALCULATIONS" << std::endl; @@ -600,15 +604,12 @@ void ESolver_KS_LCAO::nscf(void) #ifdef __EXX #ifdef __MPI // Peize Lin add 2018-08-14 - if (GlobalC::exx_info.info_global.cal_exx) - { - const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); - if (GlobalC::exx_info.info_ri.real_number) - { + if (GlobalC::exx_info.info_global.cal_exx) { + const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + + std::to_string(GlobalV::MY_RANK); + if (GlobalC::exx_info.info_ri.real_number) { this->exd->read_Hexxs_csr(file_name_exx, GlobalC::ucell); - } - else - { + } else { this->exc->read_Hexxs_csr(file_name_exx, GlobalC::ucell); } } @@ -620,72 +621,73 @@ void ESolver_KS_LCAO::nscf(void) // then when the istep is a variable of scf or nscf, // istep becomes istep-1, this should be fixed in future int istep = 0; - if (this->phsol != nullptr) - { - this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec, GlobalV::KS_SOLVER, true); - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!"); + if (this->phsol != nullptr) { + this->phsol->solve(this->p_hamilt, + this->psi[0], + this->pelec, + GlobalV::KS_SOLVER, + true); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", + "HSolver has not been initialed!"); } time_t time_finish = std::time(NULL); ModuleBase::GlobalFunc::OUT_TIME("cal_bands", time_start, time_finish); GlobalV::ofs_running << " end of band structure calculation " << std::endl; - GlobalV::ofs_running << " band eigenvalue in this processor (eV) :" << std::endl; + GlobalV::ofs_running << " band eigenvalue in this processor (eV) :" + << std::endl; const int nspin = GlobalV::NSPIN; const int nbands = GlobalV::NBANDS; - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { - if (nspin == 2) - { - if (ik == 0) - { + for (int ik = 0; ik < this->kv.get_nks(); ++ik) { + if (nspin == 2) { + if (ik == 0) { GlobalV::ofs_running << " spin up :" << std::endl; } - if (ik == (this->kv.get_nks() / 2)) - { + if (ik == (this->kv.get_nks() / 2)) { GlobalV::ofs_running << " spin down :" << std::endl; } } - GlobalV::ofs_running << " k-points" << ik + 1 << "(" << this->kv.get_nkstot() << "): " << this->kv.kvec_c[ik].x - << " " << this->kv.kvec_c[ik].y << " " << this->kv.kvec_c[ik].z << std::endl; + GlobalV::ofs_running + << " k-points" << ik + 1 << "(" << this->kv.get_nkstot() + << "): " << this->kv.kvec_c[ik].x << " " << this->kv.kvec_c[ik].y + << " " << this->kv.kvec_c[ik].z << std::endl; - for (int ib = 0; ib < nbands; ++ib) - { - GlobalV::ofs_running << " spin" << this->kv.isk[ik] + 1 << "final_state " << ib + 1 << " " - << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV << " " - << this->pelec->wg(ik, ib) * this->kv.get_nks() << std::endl; + for (int ib = 0; ib < nbands; ++ib) { + GlobalV::ofs_running + << " spin" << this->kv.isk[ik] + 1 << "final_state " << ib + 1 + << " " << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV << " " + << this->pelec->wg(ik, ib) * this->kv.get_nks() << std::endl; } GlobalV::ofs_running << std::endl; } - if (GlobalV::out_bandgap) - { - if (!GlobalV::TWO_EFERMI) - { + if (GlobalV::out_bandgap) { + if (!GlobalV::TWO_EFERMI) { this->pelec->cal_bandgap(); - GlobalV::ofs_running << " E_bandgap " << this->pelec->bandgap * ModuleBase::Ry_to_eV << " eV" << std::endl; - } - else - { + GlobalV::ofs_running << " E_bandgap " + << this->pelec->bandgap * ModuleBase::Ry_to_eV + << " eV" << std::endl; + } else { this->pelec->cal_bandgap_updw(); - GlobalV::ofs_running << " E_bandgap_up " << this->pelec->bandgap_up * ModuleBase::Ry_to_eV << " eV" - << std::endl; - GlobalV::ofs_running << " E_bandgap_dw " << this->pelec->bandgap_dw * ModuleBase::Ry_to_eV << " eV" - << std::endl; + GlobalV::ofs_running + << " E_bandgap_up " + << this->pelec->bandgap_up * ModuleBase::Ry_to_eV << " eV" + << std::endl; + GlobalV::ofs_running + << " E_bandgap_dw " + << this->pelec->bandgap_dw * ModuleBase::Ry_to_eV << " eV" + << std::endl; } } // add by jingan in 2018.11.7 - if (GlobalV::CALCULATION == "nscf" && INPUT.towannier90) - { + if (GlobalV::CALCULATION == "nscf" && INPUT.towannier90) { #ifdef __LCAO - if (INPUT.wannier_method == 1) - { + if (INPUT.wannier_method == 1) { toWannier90_LCAO_IN_PW myWannier(INPUT.out_wannier_mmn, INPUT.out_wannier_amn, INPUT.out_wannier_unk, @@ -701,9 +703,7 @@ void ESolver_KS_LCAO::nscf(void) this->kv, this->psi, &(this->orb_con.ParaV)); - } - else if (INPUT.wannier_method == 2) - { + } else if (INPUT.wannier_method == 2) { toWannier90_LCAO myWannier(INPUT.out_wannier_mmn, INPUT.out_wannier_amn, INPUT.out_wannier_unk, @@ -712,27 +712,35 @@ void ESolver_KS_LCAO::nscf(void) INPUT.nnkpfile, INPUT.wannier_spin); - myWannier.calculate(this->pelec->ekb, this->kv, *(this->psi), &(this->orb_con.ParaV)); + myWannier.calculate(this->pelec->ekb, + this->kv, + *(this->psi), + &(this->orb_con.ParaV)); } #endif } // add by jingan - if (berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1) - { + if (berryphase::berry_phase_flag + && ModuleSymmetry::Symmetry::symm_flag != 1) { berryphase bp(this->LOC); - bp.lcao_init( - this->kv, - this->GridT); // additional step before calling macroscopic_polarization (why capitalize the function name?) - bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); + bp.lcao_init(this->kv, + this->GridT); // additional step before calling + // macroscopic_polarization (why capitalize + // the function name?) + bp.Macroscopic_polarization(this->pw_wfc->npwk_max, + this->psi, + this->pw_rho, + this->pw_wfc, + this->kv); } // below is for DeePKS NSCF calculation #ifdef __DEEPKS - if (GlobalV::deepks_out_labels || GlobalV::deepks_scf) - { + if (GlobalV::deepks_out_labels || GlobalV::deepks_scf) { const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); + = dynamic_cast*>(this->pelec) + ->get_DM(); this->dpks_cal_projected_DM(dm); GlobalC::ld.cal_descriptor(GlobalC::ucell.nat); // final descriptor GlobalC::ld.cal_gedm(GlobalC::ucell.nat); @@ -742,12 +750,15 @@ void ESolver_KS_LCAO::nscf(void) this->create_Output_Mat_Sparse(0).write(); // mulliken charge analysis - if (GlobalV::out_mul) - { - elecstate::ElecStateLCAO* pelec_lcao = dynamic_cast*>(this->pelec); + if (GlobalV::out_mul) { + elecstate::ElecStateLCAO* pelec_lcao + = dynamic_cast*>(this->pelec); this->pelec->calculate_weights(); this->pelec->calEBand(); - elecstate::cal_dm_psi(&(this->orb_con.ParaV), pelec_lcao->wg, *(this->psi), *(pelec_lcao->get_DM())); + elecstate::cal_dm_psi(&(this->orb_con.ParaV), + pelec_lcao->wg, + *(this->psi), + *(pelec_lcao->get_DM())); this->cal_mag(istep, true); } diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index cc035fbaf6..628e76c374 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -31,38 +31,30 @@ //--------------------------------------------------- -namespace ModuleESolver -{ +namespace ModuleESolver { -ESolver_KS_LCAO_TDDFT::ESolver_KS_LCAO_TDDFT() -{ +ESolver_KS_LCAO_TDDFT::ESolver_KS_LCAO_TDDFT() { classname = "ESolver_KS_LCAO_TDDFT"; basisname = "LCAO"; } -ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() -{ +ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() { delete psi_laststep; - if (Hk_laststep != nullptr) - { - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + if (Hk_laststep != nullptr) { + for (int ik = 0; ik < kv.get_nks(); ++ik) { delete[] Hk_laststep[ik]; } delete[] Hk_laststep; } - if (Sk_laststep != nullptr) - { - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + if (Sk_laststep != nullptr) { + for (int ik = 0; ik < kv.get_nks(); ++ik) { delete[] Sk_laststep[ik]; } delete[] Sk_laststep; } } -void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) -{ +void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) { // 1) run "before_all_runners" in ESolver_KS ESolver_KS::before_all_runners(inp, ucell); @@ -70,22 +62,24 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, pw_rho); // 3) initialize the electronic states for TDDFT - if (this->pelec == nullptr) - { - this->pelec = new elecstate::ElecStateLCAO_TDDFT(&this->chr, - &kv, - kv.get_nks(), - &this->LOC, - &this->GK, // mohan add 2024-04-01 - this->pw_rho, - pw_big); + if (this->pelec == nullptr) { + this->pelec = new elecstate::ElecStateLCAO_TDDFT( + &this->chr, + &kv, + kv.get_nks(), + &this->LOC, + &this->GK, // mohan add 2024-04-01 + this->pw_rho, + pw_big); } // 4) read the local orbitals and construct the interpolation tables. this->init_basis_lcao(this->orb_con, inp, ucell); // 5) allocate H and S matrices according to computational resources - this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, orb_con.ParaV, kv.get_nks()); + this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, + orb_con.ParaV, + kv.get_nks()); // this part will be updated soon // pass Hamilt-pointer to Operator @@ -96,9 +90,9 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) ->init_DM(&kv, this->LM.ParaV, GlobalV::NSPIN); // 7) initialize Hsolver - if (this->phsol == nullptr) - { - this->phsol = new hsolver::HSolverLCAO>(this->LM.ParaV); + if (this->phsol == nullptr) { + this->phsol + = new hsolver::HSolverLCAO>(this->LM.ParaV); this->phsol->method = GlobalV::KS_SOLVER; } @@ -119,14 +113,13 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) this->pelec_td = dynamic_cast(this->pelec); } -void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, const double ethr) -{ +void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, + const int iter, + const double ethr) { pelec->charge->save_rho_before_sum_band(); - if (wf.init_wfc == "file") - { - if (istep >= 1) - { + if (wf.init_wfc == "file") { + if (istep >= 1) { module_tddft::Evolve_elec::solve_psi(istep, GlobalV::NBANDS, GlobalV::NLOCAL, @@ -143,9 +136,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons this->pelec_td->psiToRho_td(this->psi[0]); } this->pelec_td->psiToRho_td(this->psi[0]); - } - else if (istep >= 2) - { + } else if (istep >= 2) { module_tddft::Evolve_elec::solve_psi(istep, GlobalV::NBANDS, GlobalV::NLOCAL, @@ -160,49 +151,46 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons INPUT.propagator, kv.get_nks()); this->pelec_td->psiToRho_td(this->psi[0]); - } - else if (this->phsol != nullptr) - { + } else if (this->phsol != nullptr) { // reset energy this->pelec->f_en.eband = 0.0; this->pelec->f_en.demet = 0.0; - if (this->psi != nullptr) - { - this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec_td, GlobalV::KS_SOLVER); + if (this->psi != nullptr) { + this->phsol->solve(this->p_hamilt, + this->psi[0], + this->pelec_td, + GlobalV::KS_SOLVER); } - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "HSolver has not been initialed!"); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", + "HSolver has not been initialed!"); } // print occupation of each band - if (iter == 1 && istep <= 2) - { + if (iter == 1 && istep <= 2) { GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" + << "---------------------------------------------------------------" + "---------------------------------" << std::endl; GlobalV::ofs_running << "occupation : " << std::endl; GlobalV::ofs_running << "ik iband occ " << std::endl; GlobalV::ofs_running << std::setprecision(6); GlobalV::ofs_running << std::setiosflags(std::ios::showpoint); - for (int ik = 0; ik < kv.get_nks(); ik++) - { - for (int ib = 0; ib < GlobalV::NBANDS; ib++) - { + for (int ik = 0; ik < kv.get_nks(); ik++) { + for (int ib = 0; ib < GlobalV::NBANDS; ib++) { std::setprecision(6); - GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " << this->pelec_td->wg(ik, ib) - << std::endl; + GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " + << this->pelec_td->wg(ik, ib) << std::endl; } } GlobalV::ofs_running << std::endl; GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" + << "---------------------------------------------------------------" + "---------------------------------" << std::endl; } - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { this->pelec_td->print_band(ik, INPUT.printe, iter); } @@ -210,12 +198,14 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons this->pelec->cal_energies(1); // symmetrize the charge density only for ground state - if (istep <= 1) - { + if (istep <= 1) { Symmetry_rho srho; - for (int is = 0; is < GlobalV::NSPIN; is++) - { - srho.begin(is, *(pelec->charge), pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm); + for (int is = 0; is < GlobalV::NSPIN; is++) { + srho.begin(is, + *(pelec->charge), + pw_rho, + GlobalC::Pgrid, + GlobalC::ucell.symm); } } @@ -229,34 +219,29 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons this->pelec->f_en.deband = this->pelec->cal_delta_eband(); } -void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) -{ +void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) { // print Hamiltonian and Overlap matrix - if (this->conv_elec) - { - if (!GlobalV::GAMMA_ONLY_LOCAL) - { + if (this->conv_elec) { + if (!GlobalV::GAMMA_ONLY_LOCAL) { this->GK.renew(true); } - for (int ik = 0; ik < kv.get_nks(); ++ik) - { - if (hsolver::HSolverLCAO>::out_mat_hs[0]) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { + if (hsolver::HSolverLCAO>::out_mat_hs[0]) { this->p_hamilt->updateHk(ik); } bool bit = false; // LiuXh, 2017-03-21 - // if set bit = true, there would be error in soc-multi-core calculation, noted by zhengdy-soc - if (this->psi != nullptr && (istep % GlobalV::out_interval == 0)) - { + // if set bit = true, there would be error in soc-multi-core + // calculation, noted by zhengdy-soc + if (this->psi != nullptr && (istep % GlobalV::out_interval == 0)) { hamilt::MatrixBlock> h_mat, s_mat; this->p_hamilt->matrix(h_mat, s_mat); - if (hsolver::HSolverLCAO>::out_mat_hs[0]) - { + if (hsolver::HSolverLCAO>::out_mat_hs[0]) { ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, bit, - hsolver::HSolverLCAO>::out_mat_hs[1], + hsolver::HSolverLCAO< + std::complex>::out_mat_hs[1], 1, GlobalV::out_app_flag, "H", @@ -268,7 +253,8 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) s_mat.p, GlobalV::NLOCAL, bit, - hsolver::HSolverLCAO>::out_mat_hs[1], + hsolver::HSolverLCAO< + std::complex>::out_mat_hs[1], 1, GlobalV::out_app_flag, "S", @@ -280,30 +266,28 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) } } - if (elecstate::ElecStateLCAO>::out_wfc_lcao && (this->conv_elec || iter == GlobalV::SCF_NMAX) - && (istep % GlobalV::out_interval == 0)) - { - ModuleIO::write_wfc_nao(elecstate::ElecStateLCAO>::out_wfc_lcao, - this->psi[0], - this->pelec->ekb, - this->pelec->wg, - this->pelec->klist->kvec_c, - this->orb_con.ParaV, - istep); + if (elecstate::ElecStateLCAO>::out_wfc_lcao + && (this->conv_elec || iter == GlobalV::SCF_NMAX) + && (istep % GlobalV::out_interval == 0)) { + ModuleIO::write_wfc_nao( + elecstate::ElecStateLCAO>::out_wfc_lcao, + this->psi[0], + this->pelec->ekb, + this->pelec->wg, + this->pelec->klist->kvec_c, + this->orb_con.ParaV, + istep); } // Calculate new potential according to new Charge Density - if (!this->conv_elec) - { - if (GlobalV::NSPIN == 4) - { + if (!this->conv_elec) { + if (GlobalV::NSPIN == 4) { GlobalC::ucell.cal_ux(); } - this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->pelec->charge, + &GlobalC::ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); - } - else - { + } else { this->pelec->cal_converged(); } @@ -314,52 +298,51 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) const int nlocal = GlobalV::NLOCAL; // store wfc and Hk laststep - if (istep >= (wf.init_wfc == "file" ? 0 : 1) && this->conv_elec) - { - if (this->psi_laststep == nullptr) - { + if (istep >= (wf.init_wfc == "file" ? 0 : 1) && this->conv_elec) { + if (this->psi_laststep == nullptr) { #ifdef __MPI - this->psi_laststep = new psi::Psi>(kv.get_nks(), ncol_nbands, nrow, nullptr); + this->psi_laststep + = new psi::Psi>(kv.get_nks(), + ncol_nbands, + nrow, + nullptr); #else - this->psi_laststep = new psi::Psi>(kv.get_nks(), nbands, nlocal, nullptr); + this->psi_laststep + = new psi::Psi>(kv.get_nks(), + nbands, + nlocal, + nullptr); #endif } - if (td_htype == 1) - { - if (this->Hk_laststep == nullptr) - { + if (td_htype == 1) { + if (this->Hk_laststep == nullptr) { this->Hk_laststep = new std::complex*[kv.get_nks()]; - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { this->Hk_laststep[ik] = new std::complex[nloc]; ModuleBase::GlobalFunc::ZEROS(Hk_laststep[ik], nloc); } } - if (this->Sk_laststep == nullptr) - { + if (this->Sk_laststep == nullptr) { this->Sk_laststep = new std::complex*[kv.get_nks()]; - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { this->Sk_laststep[ik] = new std::complex[nloc]; ModuleBase::GlobalFunc::ZEROS(Sk_laststep[ik], nloc); } } } - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { this->psi->fix_k(ik); this->psi_laststep->fix_k(ik); int size0 = psi->get_nbands() * psi->get_nbasis(); - for (int index = 0; index < size0; ++index) - { - psi_laststep[0].get_pointer()[index] = psi[0].get_pointer()[index]; + for (int index = 0; index < size0; ++index) { + psi_laststep[0].get_pointer()[index] + = psi[0].get_pointer()[index]; } // store Hamiltonian - if (td_htype == 1) - { + if (td_htype == 1) { this->p_hamilt->updateHk(ik); hamilt::MatrixBlock> h_mat, s_mat; this->p_hamilt->matrix(h_mat, s_mat); @@ -369,53 +352,57 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) } // calculate energy density matrix for tddft - if (istep >= (wf.init_wfc == "file" ? 0 : 2) && module_tddft::Evolve_elec::td_edm == 0) - { + if (istep >= (wf.init_wfc == "file" ? 0 : 2) + && module_tddft::Evolve_elec::td_edm == 0) { this->cal_edm_tddft(); } } // print "eigen value" for tddft - if (this->conv_elec) - { + if (this->conv_elec) { GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" + << "---------------------------------------------------------------" + "---------------------------------" << std::endl; GlobalV::ofs_running << "Eii : " << std::endl; GlobalV::ofs_running << "ik iband Eii (eV)" << std::endl; GlobalV::ofs_running << std::setprecision(6); GlobalV::ofs_running << std::setiosflags(std::ios::showpoint); - for (int ik = 0; ik < kv.get_nks(); ik++) - { - for (int ib = 0; ib < GlobalV::NBANDS; ib++) - { - GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " - << this->pelec_td->ekb(ik, ib) * ModuleBase::Ry_to_eV << std::endl; + for (int ik = 0; ik < kv.get_nks(); ik++) { + for (int ib = 0; ib < GlobalV::NBANDS; ib++) { + GlobalV::ofs_running + << ik + 1 << " " << ib + 1 << " " + << this->pelec_td->ekb(ik, ib) * ModuleBase::Ry_to_eV + << std::endl; } } GlobalV::ofs_running << std::endl; GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" + << "---------------------------------------------------------------" + "---------------------------------" << std::endl; } } -void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) -{ - for (int is = 0; is < GlobalV::NSPIN; is++) - { - if (module_tddft::Evolve_elec::out_dipole == 1) - { +void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) { + for (int is = 0; is < GlobalV::NSPIN; is++) { + if (module_tddft::Evolve_elec::out_dipole == 1) { std::stringstream ss_dipole; - ss_dipole << GlobalV::global_out_dir << "SPIN" << is + 1 << "_DIPOLE"; - ModuleIO::write_dipole(pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str()); + ss_dipole << GlobalV::global_out_dir << "SPIN" << is + 1 + << "_DIPOLE"; + ModuleIO::write_dipole(pelec->charge->rho_save[is], + pelec->charge->rhopw, + is, + istep, + ss_dipole.str()); } } - if (module_tddft::Evolve_elec::out_current == 1) - { + if (module_tddft::Evolve_elec::out_current == 1) { elecstate::DensityMatrix, double>* tmp_DM - = dynamic_cast>*>(this->pelec)->get_DM(); + = dynamic_cast>*>( + this->pelec) + ->get_DM(); ModuleIO::write_current(istep, this->psi, pelec, @@ -428,25 +415,35 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) ESolver_KS_LCAO, double>::after_scf(istep); } -// use the original formula (Hamiltonian matrix) to calculate energy density matrix -void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() -{ +// use the original formula (Hamiltonian matrix) to calculate energy density +// matrix +void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() { // mohan add 2024-03-27 const int nlocal = GlobalV::NLOCAL; assert(nlocal >= 0); // this->LOC.edm_k_tddft.resize(kv.get_nks()); - dynamic_cast>*>(this->pelec)->get_DM()->EDMK.resize(kv.get_nks()); - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + dynamic_cast>*>(this->pelec) + ->get_DM() + ->EDMK.resize(kv.get_nks()); + for (int ik = 0; ik < kv.get_nks(); ++ik) { std::complex* tmp_dmk - = dynamic_cast>*>(this->pelec)->get_DM()->get_DMK_pointer(ik); + = dynamic_cast>*>( + this->pelec) + ->get_DM() + ->get_DMK_pointer(ik); ModuleBase::ComplexMatrix& tmp_edmk - = dynamic_cast>*>(this->pelec)->get_DM()->EDMK[ik]; + = dynamic_cast>*>( + this->pelec) + ->get_DM() + ->EDMK[ik]; const Parallel_Orbitals* tmp_pv - = dynamic_cast>*>(this->pelec)->get_DM()->get_paraV_pointer(); + = dynamic_cast>*>( + this->pelec) + ->get_DM() + ->get_paraV_pointer(); #ifdef __MPI @@ -485,7 +482,14 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() int info = 0; const int one_int = 1; - pzgetrf_(&nlocal, &nlocal, Sinv, &one_int, &one_int, this->orb_con.ParaV.desc, ipiv.data(), &info); + pzgetrf_(&nlocal, + &nlocal, + Sinv, + &one_int, + &one_int, + this->orb_con.ParaV.desc, + ipiv.data(), + &info); int lwork = -1; int liwork = -1; @@ -644,10 +648,8 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() p_hamilt->matrix(h_mat, s_mat); // cout<<"hmat "<* + // I just use ModuleBase::ComplexMatrix temporarily, and will change it + // to complex* ModuleBase::ComplexMatrix tmp_dmk_base(nlocal, nlocal); - for (int i = 0; i < nlocal; i++) - { - for (int j = 0; j < nlocal; j++) - { + for (int i = 0; i < nlocal; i++) { + for (int j = 0; j < nlocal; j++) { tmp_dmk_base(i, j) = tmp_dmk[i * nlocal + j]; } } - tmp_edmk = 0.5 * (Sinv * Htmp * tmp_dmk_base + tmp_dmk_base * Htmp * Sinv); + tmp_edmk + = 0.5 * (Sinv * Htmp * tmp_dmk_base + tmp_dmk_base * Htmp * Sinv); delete[] work; #endif } diff --git a/source/module_io/io_dmk.cpp b/source/module_io/io_dmk.cpp index 24fd730a4c..23776ac3d4 100644 --- a/source/module_io/io_dmk.cpp +++ b/source/module_io/io_dmk.cpp @@ -48,62 +48,56 @@ Direct -0.0883978533958687 (fermi energy) 10 10 - 5.773e-01 3.902e-02 1.661e-02 4.797e-17 -2.255e-17 5.773e-01 3.902e-02 -1.661e-02 - -1.461e-17 -4.414e-17 + 5.773e-01 3.902e-02 1.661e-02 4.797e-17 -2.255e-17 5.773e-01 3.902e-02 +-1.661e-02 -1.461e-17 -4.414e-17 ... ''' */ - -std::string ModuleIO::dmk_gen_fname(const bool gamma_only, const int ispin, const int ik) -{ - if (gamma_only) - { +std::string ModuleIO::dmk_gen_fname(const bool gamma_only, + const int ispin, + const int ik) { + if (gamma_only) { return "SPIN" + std::to_string(ispin + 1) + "_DM"; - } - else - { + } else { // this case is not implemented now. - ModuleBase::WARNING_QUIT("dmk_gen_fname", "Not implemented for non-gamma_only case."); + ModuleBase::WARNING_QUIT("dmk_gen_fname", + "Not implemented for non-gamma_only case."); } } -void ModuleIO::dmk_write_ucell(std::ofstream& ofs, const UnitCell* ucell) -{ +void ModuleIO::dmk_write_ucell(std::ofstream& ofs, const UnitCell* ucell) { // write the UnitCell information ofs << ucell->latName << std::endl; ofs << " " << ucell->lat0 * ModuleBase::BOHR_TO_A << std::endl; - ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " << ucell->latvec.e13 << std::endl; - ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " << ucell->latvec.e23 << std::endl; - ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " << ucell->latvec.e33 << std::endl; - for (int it = 0; it < ucell->ntype; it++) - { + ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " + << ucell->latvec.e13 << std::endl; + ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " + << ucell->latvec.e23 << std::endl; + ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " + << ucell->latvec.e33 << std::endl; + for (int it = 0; it < ucell->ntype; it++) { ofs << " " << ucell->atoms[it].label; } ofs << std::endl; - for (int it = 0; it < ucell->ntype; it++) - { + for (int it = 0; it < ucell->ntype; it++) { ofs << " " << ucell->atoms[it].na; } ofs << std::endl; ofs << "Direct" << std::endl; - for (int it = 0; it < ucell->ntype; it++) - { + for (int it = 0; it < ucell->ntype; it++) { Atom* atom = &ucell->atoms[it]; ofs << std::setprecision(15); - for (int ia = 0; ia < ucell->atoms[it].na; ia++) - { - ofs << " " << atom->taud[ia].x << " " << atom->taud[ia].y << " " << atom->taud[ia].z - << std::endl; + for (int ia = 0; ia < ucell->atoms[it].na; ia++) { + ofs << " " << atom->taud[ia].x << " " << atom->taud[ia].y << " " + << atom->taud[ia].z << std::endl; } } } -void ModuleIO::dmk_read_ucell(std::ifstream& ifs) -{ +void ModuleIO::dmk_read_ucell(std::ifstream& ifs) { std::string tmp; - for (int i=0;i<6;i++) - { + for (int i = 0; i < 6; i++) { std::getline(ifs, tmp); // latName + lat0 + latvec + atom label } std::getline(ifs, tmp); // atom number of each type @@ -111,23 +105,17 @@ void ModuleIO::dmk_read_ucell(std::ifstream& ifs) std::istringstream iss(tmp); int natom = 0; int total_natom = 0; - while (iss >> natom) - { + while (iss >> natom) { total_natom += natom; } - for (int i=0;i> data; -} +void ModuleIO::dmk_readData(std::ifstream& ifs, double& data) { ifs >> data; } -void ModuleIO::dmk_readData(std::ifstream& ifs, std::complex& data) -{ +void ModuleIO::dmk_readData(std::ifstream& ifs, std::complex& data) { double real, imag; ifs >> real; ifs >> imag; @@ -135,14 +123,11 @@ void ModuleIO::dmk_readData(std::ifstream& ifs, std::complex& data) } template -bool ModuleIO::read_dmk( - const int nspin, - const int nk, - const Parallel_2D& pv, - const std::string& dmk_dir, - std::vector>& dmk -) -{ +bool ModuleIO::read_dmk(const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>& dmk) { ModuleBase::TITLE("ModuleIO", "read_dmk"); ModuleBase::timer::tick("ModuleIO", "read_dmk"); @@ -150,17 +135,22 @@ bool ModuleIO::read_dmk( #ifdef __MPI MPI_Comm_rank(pv.comm_2D, &my_rank); #endif - + int nlocal = pv.get_global_row_size(); bool gamma_only = std::is_same::value; std::vector> dmk_global; - + // write a lambda function to check the consistency of the data - auto check_consistency = [&](const std::string& fn, const std::string& name, const std::string& value, const int& target) { - if (std::stoi(value) != target) - { - ModuleBase::WARNING("ModuleIO::read_dmk", name + " is not consistent in file < " + fn + " >."); - std::cout << name << " = " << target << ", " << name << " in file = " << value << std::endl; + auto check_consistency = [&](const std::string& fn, + const std::string& name, + const std::string& value, + const int& target) { + if (std::stoi(value) != target) { + ModuleBase::WARNING("ModuleIO::read_dmk", + name + " is not consistent in file < " + fn + + " >."); + std::cout << name << " = " << target << ", " << name + << " in file = " << value << std::endl; return false; } return true; @@ -168,20 +158,18 @@ bool ModuleIO::read_dmk( bool read_success = true; std::string tmp; - if (my_rank == 0) - { - dmk_global.resize(nspin * nk, std::vector(nlocal*nlocal)); - - for (int ispin = 0; ispin < nspin; ispin++) - { - for (int ik = 0; ik < nk; ik++) - { + if (my_rank == 0) { + dmk_global.resize(nspin * nk, std::vector(nlocal * nlocal)); + + for (int ispin = 0; ispin < nspin; ispin++) { + for (int ik = 0; ik < nk; ik++) { std::string fn = dmk_dir + dmk_gen_fname(gamma_only, ispin, ik); std::ifstream ifs(fn.c_str()); - if (!ifs) - { - ModuleBase::WARNING("ModuleIO::read_dmk", "Can't open DENSITY MATRIX File < " + fn + " >."); + if (!ifs) { + ModuleBase::WARNING("ModuleIO::read_dmk", + "Can't open DENSITY MATRIX File < " + fn + + " >."); read_success = false; break; } @@ -190,58 +178,55 @@ bool ModuleIO::read_dmk( dmk_read_ucell(ifs); ifs >> tmp; // nspin - if (!check_consistency(fn, "ispin", tmp, ispin+1)) - { + if (!check_consistency(fn, "ispin", tmp, ispin + 1)) { read_success = false; ifs.close(); break; } - ifs >> tmp; ifs >> tmp; ifs >> tmp;// fermi energy + ifs >> tmp; + ifs >> tmp; + ifs >> tmp; // fermi energy ifs >> tmp; // nlocal - if (!check_consistency(fn, "nlocal", tmp, nlocal)) - { + if (!check_consistency(fn, "nlocal", tmp, nlocal)) { read_success = false; ifs.close(); break; } ifs >> tmp; // nlocal - if (!check_consistency(fn, "nlocal", tmp, nlocal)) - { + if (!check_consistency(fn, "nlocal", tmp, nlocal)) { read_success = false; ifs.close(); break; } // read the DMK data - for (int i = 0; i < nlocal; ++i) - { - for (int j = 0; j < nlocal; ++j) - { - dmk_readData(ifs, dmk_global[ik + nk * ispin][i * nlocal + j]); + for (int i = 0; i < nlocal; ++i) { + for (int j = 0; j < nlocal; ++j) { + dmk_readData( + ifs, + dmk_global[ik + nk * ispin][i * nlocal + j]); } } ifs.close(); } // ik - if (!read_success) - { + if (!read_success) { break; } - } // ispin - } // rank0 + } // ispin + } // rank0 #ifdef __MPI MPI_Bcast(&read_success, 1, MPI_C_BOOL, 0, pv.comm_2D); #endif - if (read_success) - { -#ifdef __MPI + if (read_success) { +#ifdef __MPI // seperate dmk data to each processor with 2D block distribution - dmk.resize(nspin * nk, std::vector(pv.get_row_size()*pv.get_col_size())); + dmk.resize(nspin * nk, + std::vector(pv.get_row_size() * pv.get_col_size())); Parallel_2D pv_glb; pv_glb.set(nlocal, nlocal, nlocal, pv.comm_2D, pv.blacs_ctxt); - for(int ik=0;ik < nspin * nk; ik++) - { + for (int ik = 0; ik < nspin * nk; ik++) { Cpxgemr2d(nlocal, nlocal, dmk_global[ik].data(), @@ -262,14 +247,12 @@ bool ModuleIO::read_dmk( return read_success; } - template void ModuleIO::write_dmk(const std::vector>& dmk, const int precision, const std::vector& efs, const UnitCell* ucell, - const Parallel_2D& pv) -{ + const Parallel_2D& pv) { ModuleBase::TITLE("ModuleIO", "write_dmk"); ModuleBase::timer::tick("ModuleIO", "write_dmk"); @@ -282,17 +265,17 @@ void ModuleIO::write_dmk(const std::vector>& dmk, int nlocal = pv.get_global_row_size(); int nspin = efs.size(); int nk = dmk.size() / nspin; - if (nk * nspin != dmk.size()) - { - ModuleBase::WARNING_QUIT("write_dmk", "The size of dmk is not consistent with nspin and nk."); + if (nk * nspin != dmk.size()) { + ModuleBase::WARNING_QUIT( + "write_dmk", + "The size of dmk is not consistent with nspin and nk."); } Parallel_2D pv_glb; - // when nspin == 2, assume the order of K in dmk is K1_up, K2_up, ..., K1_down, K2_down, ... - for (int ispin = 0; ispin < nspin; ispin++) - { - for (int ik = 0; ik < nk; ik++) - { + // when nspin == 2, assume the order of K in dmk is K1_up, K2_up, ..., + // K1_down, K2_down, ... + for (int ispin = 0; ispin < nspin; ispin++) { + for (int ik = 0; ik < nk; ik++) { // gather dmk[ik] to dmk_global std::vector dmk_global(my_rank == 0 ? nlocal * nlocal : 0); #ifdef __MPI @@ -312,42 +295,40 @@ void ModuleIO::write_dmk(const std::vector>& dmk, dmk_global = dmk[ik + nk * ispin]; #endif - if (my_rank == 0) - { - std::string fn = GlobalV::global_out_dir + dmk_gen_fname(gamma_only, ispin, ik); + if (my_rank == 0) { + std::string fn = GlobalV::global_out_dir + + dmk_gen_fname(gamma_only, ispin, ik); std::ofstream ofs(fn.c_str()); - if (!ofs) - { - ModuleBase::WARNING("ModuleIO::write_dmk", "Can't create DENSITY MATRIX File < " + fn + " >."); + if (!ofs) { + ModuleBase::WARNING("ModuleIO::write_dmk", + "Can't create DENSITY MATRIX File < " + + fn + " >."); continue; } // write the UnitCell information dmk_write_ucell(ofs, ucell); - - + ofs << "\n " << dmk.size(); // nspin - ofs << "\n " << std::fixed << std::setprecision(5) << efs[ispin] << " (fermi energy)"; + ofs << "\n " << std::fixed << std::setprecision(5) << efs[ispin] + << " (fermi energy)"; ofs << "\n " << nlocal << " " << nlocal << std::endl; ofs << std::setprecision(precision); ofs << std::scientific; - for (int i = 0; i < nlocal; ++i) - { - for (int j = 0; j < nlocal; ++j) - { - if (j % 8 == 0) - { + for (int i = 0; i < nlocal; ++i) { + for (int j = 0; j < nlocal; ++j) { + if (j % 8 == 0) { ofs << "\n"; } - if(std::is_same::value) - { + if (std::is_same::value) { ofs << " " << dmk_global[i * nlocal + j]; - } - else if (std::is_same, T>::value) - { - ofs << " (" << std::real(dmk_global[i * nlocal + j]) << "," << std::imag(dmk_global[i * nlocal + j]) << ")"; + } else if (std::is_same, + T>::value) { + ofs << " (" << std::real(dmk_global[i * nlocal + j]) + << "," << std::imag(dmk_global[i * nlocal + j]) + << ")"; } } } @@ -365,20 +346,23 @@ template bool ModuleIO::read_dmk(const int nspin, const std::string& dmk_dir, std::vector>& dmk); -template bool ModuleIO::read_dmk>(const int nspin, - const int nk, - const Parallel_2D& pv, - const std::string& dmk_dir, - std::vector>>& dmk); - -template void ModuleIO::write_dmk(const std::vector>& dmk, - const int precision, - const std::vector& efs, - const UnitCell* ucell, - const Parallel_2D& pv); - -template void ModuleIO::write_dmk>(const std::vector>>& dmk, - const int precision, - const std::vector& efs, - const UnitCell* ucell, - const Parallel_2D& pv); \ No newline at end of file +template bool ModuleIO::read_dmk>( + const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>>& dmk); + +template void + ModuleIO::write_dmk(const std::vector>& dmk, + const int precision, + const std::vector& efs, + const UnitCell* ucell, + const Parallel_2D& pv); + +template void ModuleIO::write_dmk>( + const std::vector>>& dmk, + const int precision, + const std::vector& efs, + const UnitCell* ucell, + const Parallel_2D& pv); \ No newline at end of file diff --git a/source/module_io/io_dmk.h b/source/module_io/io_dmk.h index 9af9ef1551..bbd35946d8 100644 --- a/source/module_io/io_dmk.h +++ b/source/module_io/io_dmk.h @@ -7,8 +7,7 @@ #include #include -namespace ModuleIO -{ +namespace ModuleIO { /** * @brief Generates the filename for the DMK file based on the given parameters. @@ -51,9 +50,11 @@ void dmk_readData(std::ifstream& ifs, std::complex& data); * @tparam T The type of the DMK data. * @param nspin The number of spin components. * @param nk The number of k-points. - * @param pv The Parallel_2D object. Will get the global size and local size from it, and seperate the data into different processors accordingly. + * @param pv The Parallel_2D object. Will get the global size and local size + * from it, and seperate the data into different processors accordingly. * @param dmk_dir The directory path of the DMK file. - * @param dmk A vector to store the DMK data. If use MPI parallel, the data will be seperated into different processors based on the Parallel_2D object. + * @param dmk A vector to store the DMK data. If use MPI parallel, the data will + * be seperated into different processors based on the Parallel_2D object. * @return True if the DMK data is successfully read, false otherwise. */ template @@ -67,10 +68,12 @@ bool read_dmk(const int nspin, * @brief Writes the DMK data to a file. * * @tparam T The type of the DMK data. - * @param dmk A vector containing the DMK data. The first dimension is nspin*nk, and the second dimension is - * nlocal*nlocal. DMK is parallel in 2d-block type if using MPI. + * @param dmk A vector containing the DMK data. The first dimension is nspin*nk, + * and the second dimension is nlocal*nlocal. DMK is parallel in 2d-block type + * if using MPI. * @param precision The precision of the output of DMK. - * @param efs A vector containing the Fermi energies, and should have the same size as the number of SPIN. + * @param efs A vector containing the Fermi energies, and should have the same + * size as the number of SPIN. * @param ucell A pointer to the UnitCell object. * @param pv The Parallel_2D object. The 2d-block parallel information of DMK. */ diff --git a/source/module_io/test_serial/io_dmk_test.cpp b/source/module_io/test_serial/io_dmk_test.cpp index f8c5ff11cf..9d4baca6e5 100644 --- a/source/module_io/test_serial/io_dmk_test.cpp +++ b/source/module_io/test_serial/io_dmk_test.cpp @@ -7,29 +7,17 @@ #include "gtest/gtest.h" #ifdef __LCAO -InfoNonlocal::InfoNonlocal() -{ -} -InfoNonlocal::~InfoNonlocal() -{ -} -LCAO_Orbitals::LCAO_Orbitals() -{ -} -LCAO_Orbitals::~LCAO_Orbitals() -{ -} +InfoNonlocal::InfoNonlocal() {} +InfoNonlocal::~InfoNonlocal() {} +LCAO_Orbitals::LCAO_Orbitals() {} +LCAO_Orbitals::~LCAO_Orbitals() {} #endif -Magnetism::Magnetism() -{ +Magnetism::Magnetism() { this->tot_magnetization = 0.0; this->abs_magnetization = 0.0; this->start_magnetization = nullptr; } -Magnetism::~Magnetism() -{ - delete[] this->start_magnetization; -} +Magnetism::~Magnetism() { delete[] this->start_magnetization; } /************************************************ * unit test of read_dmk and write_dmk @@ -45,8 +33,7 @@ Magnetism::~Magnetism() * - the serial version without MPI */ -TEST(DMKTest,GenFileName) -{ +TEST(DMKTest, GenFileName) { std::string fname = ModuleIO::dmk_gen_fname(true, 0, 0); EXPECT_EQ(fname, "SPIN1_DM"); fname = ModuleIO::dmk_gen_fname(true, 1, 1); @@ -55,12 +42,13 @@ TEST(DMKTest,GenFileName) // do not support non-gamma-only case now std::string output; testing::internal::CaptureStdout(); - EXPECT_EXIT(ModuleIO::dmk_gen_fname(false, 2, 0), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(ModuleIO::dmk_gen_fname(false, 2, 0), + ::testing::ExitedWithCode(0), + ""); output = testing::internal::GetCapturedStdout(); }; -class DMKIOTest : public ::testing::Test -{ +class DMKIOTest : public ::testing::Test { protected: int nspin = 2; int nk = 1; @@ -69,21 +57,17 @@ class DMKIOTest : public ::testing::Test Parallel_2D pv; std::vector efs; - void SetUp() - { - dmk.resize(nspin*nk, std::vector(nlocal * nlocal, 0.0)); - for (int i=0;i(nlocal * nlocal, 0.0)); + for (int i = 0; i < nspin * nk; i++) { + for (int j = 0; j < nlocal * nlocal; j++) { + dmk[i][j] = 1.0 * i + 0.1 * j; } } efs.resize(nspin, 0.0); - for (int i=0;i(ifs)), std::istreambuf_iterator()); + std::string str((std::istreambuf_iterator(ifs)), + std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("0.00000 (fermi energy)")); EXPECT_THAT(str, testing::HasSubstr("20 20")); - EXPECT_THAT(str, testing::HasSubstr("0.000e+00 1.000e-01 2.000e-01 3.000e-01 4.000e-01 5.000e-01 6.000e-01 7.000e-01\n")); - EXPECT_THAT(str, testing::HasSubstr("8.000e-01 9.000e-01 1.000e+00 1.100e+00 1.200e+00 1.300e+00 1.400e+00 1.500e+00\n")); - EXPECT_THAT(str, testing::HasSubstr("1.600e+00 1.700e+00 1.800e+00 1.900e+00\n")); + EXPECT_THAT( + str, + testing::HasSubstr("0.000e+00 1.000e-01 2.000e-01 3.000e-01 4.000e-01 " + "5.000e-01 6.000e-01 7.000e-01\n")); + EXPECT_THAT( + str, + testing::HasSubstr("8.000e-01 9.000e-01 1.000e+00 1.100e+00 1.200e+00 " + "1.300e+00 1.400e+00 1.500e+00\n")); + EXPECT_THAT( + str, + testing::HasSubstr("1.600e+00 1.700e+00 1.800e+00 1.900e+00\n")); ifs.close(); fn = "SPIN2_DM"; ifs.open(fn); - str = std::string((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); + str = std::string((std::istreambuf_iterator(ifs)), + std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("0.10000 (fermi energy)")); EXPECT_THAT(str, testing::HasSubstr("20 20")); - EXPECT_THAT(str, testing::HasSubstr("1.000e+00 1.100e+00 1.200e+00 1.300e+00 1.400e+00 1.500e+00 1.600e+00 1.700e+00\n")); - EXPECT_THAT(str, testing::HasSubstr("1.800e+00 1.900e+00 2.000e+00 2.100e+00 2.200e+00 2.300e+00 2.400e+00 2.500e+00\n")); - EXPECT_THAT(str, testing::HasSubstr("2.600e+00 2.700e+00 2.800e+00 2.900e+00\n")); + EXPECT_THAT( + str, + testing::HasSubstr("1.000e+00 1.100e+00 1.200e+00 1.300e+00 1.400e+00 " + "1.500e+00 1.600e+00 1.700e+00\n")); + EXPECT_THAT( + str, + testing::HasSubstr("1.800e+00 1.900e+00 2.000e+00 2.100e+00 2.200e+00 " + "2.300e+00 2.400e+00 2.500e+00\n")); + EXPECT_THAT( + str, + testing::HasSubstr("2.600e+00 2.700e+00 2.800e+00 2.900e+00\n")); ifs.close(); delete ucell; @@ -127,13 +128,12 @@ TEST_F(DMKIOTest,WriteDMK) remove("SPIN2_DM"); }; -TEST_F(DMKIOTest,ReadDMK) -{ +TEST_F(DMKIOTest, ReadDMK) { pv.nrow = 26; pv.ncol = 26; EXPECT_TRUE(ModuleIO::read_dmk(1, 1, pv, "./support/", dmk)); EXPECT_EQ(dmk.size(), 1); - EXPECT_EQ(dmk[0].size(), 26*26); + EXPECT_EQ(dmk[0].size(), 26 * 26); EXPECT_NEAR(dmk[0][0], 3.904e-01, 1e-6); - EXPECT_NEAR(dmk[0][25*26+25], 3.445e-02, 1e-6); + EXPECT_NEAR(dmk[0][25 * 26 + 25], 3.445e-02, 1e-6); }