From 7354d2360b12b3fefe9fdf596d411d8bc87eb883 Mon Sep 17 00:00:00 2001 From: Qianrui Liu <76200646+Qianruipku@users.noreply.github.com> Date: Fri, 14 Jun 2024 09:22:21 +0800 Subject: [PATCH] Refactor: remove Kubo-Greenwood and DOS functions from esolver (#4355) * Fix: makefile cannot compile tests in module_pw * Refactor: remove some functions from esolver mv Kubo-Greenwood and DOS functions to specific classes * fix when merge * fix wrong use of sKG * change pointer to vector * [pre-commit.ci lite] apply automatic fixes * update results * merge --------- Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> --- source/Makefile.Objects | 6 +- source/module_esolver/CMakeLists.txt | 2 - source/module_esolver/esolver_ks_pw.cpp | 789 ++++----- source/module_esolver/esolver_ks_pw.h | 151 +- source/module_esolver/esolver_sdft_pw.cpp | 185 +-- source/module_esolver/esolver_sdft_pw.h | 87 - .../module_esolver/esolver_sdft_pw_tool.cpp | 1433 ----------------- .../hamilt_pwdft/CMakeLists.txt | 1 + .../hamilt_pwdft/elecond.cpp} | 204 +-- .../module_hamilt_pw/hamilt_pwdft/elecond.h | 74 + .../hamilt_stodft/CMakeLists.txt | 7 + .../hamilt_stodft/sto_dos.cpp | 238 +++ .../module_hamilt_pw/hamilt_stodft/sto_dos.h | 51 + .../hamilt_stodft/sto_elecond.cpp | 814 ++++++++++ .../hamilt_stodft/sto_elecond.h | 72 + .../hamilt_stodft/sto_tool.cpp | 125 ++ .../module_hamilt_pw/hamilt_stodft/sto_tool.h | 108 ++ .../hamilt_stodft/test/CMakeLists.txt | 7 + .../hamilt_stodft/test/test_sto_tool.cpp | 75 + .../186_PW_SKG_10D10S/refOnsager.txt | 70 +- .../186_PW_SNLKG_10D10S/refOnsager.txt | 64 +- 21 files changed, 2087 insertions(+), 2476 deletions(-) delete mode 100644 source/module_esolver/esolver_sdft_pw_tool.cpp rename source/{module_esolver/esolver_ks_pw_tool.cpp => module_hamilt_pw/hamilt_pwdft/elecond.cpp} (53%) create mode 100644 source/module_hamilt_pw/hamilt_pwdft/elecond.h create mode 100644 source/module_hamilt_pw/hamilt_stodft/sto_dos.cpp create mode 100644 source/module_hamilt_pw/hamilt_stodft/sto_dos.h create mode 100644 source/module_hamilt_pw/hamilt_stodft/sto_elecond.cpp create mode 100644 source/module_hamilt_pw/hamilt_stodft/sto_elecond.h create mode 100644 source/module_hamilt_pw/hamilt_stodft/sto_tool.cpp create mode 100644 source/module_hamilt_pw/hamilt_stodft/sto_tool.h create mode 100644 source/module_hamilt_pw/hamilt_stodft/test/CMakeLists.txt create mode 100644 source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index c515dc426b..7e9d8f2884 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -210,9 +210,7 @@ OBJS_ESOLVER=esolver.o\ esolver_ks.o\ esolver_fp.o\ esolver_ks_pw.o\ - esolver_ks_pw_tool.o\ esolver_sdft_pw.o\ - esolver_sdft_pw_tool.o\ esolver_lj.o\ esolver_dp.o\ esolver_of.o\ @@ -595,6 +593,10 @@ OBJS_SRCPW=H_Ewald_pw.o\ symmetry_rhog.o\ wavefunc.o\ wf_atomic.o\ + elecond.o\ + sto_tool.o\ + sto_elecond.o\ + sto_dos.o OBJS_VDW=vdw.o\ vdwd2_parameters.o\ diff --git a/source/module_esolver/CMakeLists.txt b/source/module_esolver/CMakeLists.txt index 44fe0eabe1..2ea3a9ce78 100644 --- a/source/module_esolver/CMakeLists.txt +++ b/source/module_esolver/CMakeLists.txt @@ -3,9 +3,7 @@ list(APPEND objects esolver_ks.cpp esolver_fp.cpp esolver_ks_pw.cpp - esolver_ks_pw_tool.cpp esolver_sdft_pw.cpp - esolver_sdft_pw_tool.cpp esolver_lj.cpp esolver_dp.cpp esolver_of.cpp diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index d43739c4b0..4daaf5fbbf 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -1,13 +1,14 @@ #include "esolver_ks_pw.h" -#include - +#include "module_hamilt_pw/hamilt_pwdft/elecond.h" +#include "module_io/input_conv.h" #include "module_io/nscf_band.h" +#include "module_io/output_log.h" #include "module_io/write_dos_pw.h" #include "module_io/write_istate_info.h" #include "module_io/write_wfc_pw.h" -#include "module_io/output_log.h" -#include "module_io/input_conv.h" + +#include //--------------temporary---------------------------- #include "module_elecstate/module_charge/symmetry_rho.h" @@ -21,6 +22,7 @@ #include "module_hamilt_pw/hamilt_pwdft/stress_pw.h" //--------------------------------------------------- #include "module_base/memory.h" +#include "module_base/module_device/device.h" #include "module_elecstate/elecstate_pw.h" #include "module_hamilt_general/module_vdw/vdw.h" #include "module_hamilt_pw/hamilt_pwdft/hamilt_pw.h" @@ -32,17 +34,16 @@ #include "module_io/numerical_basis.h" #include "module_io/numerical_descriptor.h" #include "module_io/rho_io.h" -#include "module_io/write_pot.h" #include "module_io/to_wannier90_pw.h" #include "module_io/winput.h" +#include "module_io/write_pot.h" #include "module_io/write_wfc_r.h" -#include "module_base/module_device/device.h" //--------------------------------------------------- #include "module_psi/psi_initializer_atomic.h" -#include "module_psi/psi_initializer_nao.h" -#include "module_psi/psi_initializer_random.h" #include "module_psi/psi_initializer_atomic_random.h" +#include "module_psi/psi_initializer_nao.h" #include "module_psi/psi_initializer_nao_random.h" +#include "module_psi/psi_initializer_random.h" //--------------------------------------------------- #ifdef USE_PAW #include "module_cell/module_paw/paw_cell.h" @@ -107,19 +108,15 @@ ESolver_KS_PW::~ESolver_KS_PW() delete this->psi; } - template -void ESolver_KS_PW::Init_GlobalC( - Input &inp, - UnitCell &ucell, - pseudopot_cell_vnl &ppcell) +void ESolver_KS_PW::Init_GlobalC(Input& inp, UnitCell& ucell, pseudopot_cell_vnl& ppcell) { // GlobalC is a historically left-over namespace, it is used to store global classes, // including: // pseudopot_cell_vnl: pseudopotential in cell, V non-local // UnitCell: cell information with atomic properties - // Grid_Driver: - // Parallel_Grid: + // Grid_Driver: + // Parallel_Grid: // Parallel_Kpoints: // Restart: // Exx_Info: @@ -158,9 +155,10 @@ void ESolver_KS_PW::Init_GlobalC( else // old method { // old method explicitly requires variables such as total number of kpoints, number - // of bands, number of G-vectors, and so on. Comparatively in new method, these + // of bands, number of G-vectors, and so on. Comparatively in new method, these // variables are imported in function called initialize. - this->psi = this->wf.allocate(this->kv.get_nkstot(), this->kv.get_nks(), this->kv.ngk.data(), this->pw_wfc->npwk_max); + this->psi + = this->wf.allocate(this->kv.get_nkstot(), this->kv.get_nks(), this->kv.ngk.data(), this->pw_wfc->npwk_max); } // --------------------------------------------------------------------------------- @@ -192,9 +190,7 @@ void ESolver_KS_PW::Init_GlobalC( } // --------------------------------------------------------------------------------- - - this->kspw_psi = GlobalV::device_flag == "gpu" - || GlobalV::precision_flag == "single" + this->kspw_psi = GlobalV::device_flag == "gpu" || GlobalV::precision_flag == "single" ? new psi::Psi(this->psi[0]) : reinterpret_cast*>(this->psi); @@ -206,14 +202,12 @@ void ESolver_KS_PW::Init_GlobalC( ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS"); } - template void ESolver_KS_PW::before_all_runners(Input& inp, UnitCell& ucell) { // 1) call before_all_runners() of ESolver_KS ESolver_KS::before_all_runners(inp, ucell); - // 2) initialize HSolver if (this->phsol == nullptr) { @@ -223,14 +217,9 @@ void ESolver_KS_PW::before_all_runners(Input& inp, UnitCell& ucell) // 3) initialize ElecState, if (this->pelec == nullptr) { - this->pelec = new elecstate::ElecStatePW(this->pw_wfc, - &(this->chr), - &(this->kv), - &ucell, - &GlobalC::ppcell, - this->pw_rhod, - this->pw_rho, - this->pw_big); + this->pelec + = new elecstate::ElecStatePW(this->pw_wfc, &(this->chr), &(this->kv), &ucell, &GlobalC::ppcell, + this->pw_rhod, this->pw_rho, this->pw_big); } //! 4) inititlize the charge density. @@ -242,27 +231,19 @@ void ESolver_KS_PW::before_all_runners(Input& inp, UnitCell& ucell) //! 6) initialize the potential. if (this->pelec->pot == nullptr) { - this->pelec->pot = new elecstate::Potential(this->pw_rhod, - this->pw_rho, - &ucell, - &GlobalC::ppcell.vloc, - &(this->sf), - &(this->pelec->f_en.etxc), - &(this->pelec->f_en.vtxc)); + this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho, &ucell, &GlobalC::ppcell.vloc, + &(this->sf), &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc)); } - //! 7) initialize the electronic wave functions if (GlobalV::psi_initializer) { this->allocate_psi_init(); } - //! 8) setup global classes this->Init_GlobalC(inp, ucell, GlobalC::ppcell); - //! 9) setup occupations if (GlobalV::ocp) { @@ -270,7 +251,6 @@ void ESolver_KS_PW::before_all_runners(Input& inp, UnitCell& ucell) } } - template void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) { @@ -281,25 +261,16 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) if (GlobalV::md_prec_level == 2) { - this->pw_wfc->initgrids( - ucell.lat0, - ucell.latvec, - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz); + this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz); - this->pw_wfc->initparameters( - false, - inp.ecutwfc, - this->kv.get_nks(), - this->kv.kvec_d.data()); + this->pw_wfc->initparameters(false, inp.ecutwfc, this->kv.get_nks(), this->kv.kvec_d.data()); #ifdef __MPI - if (INPUT.pw_seed > 0) - { - MPI_Allreduce(MPI_IN_PLACE, &this->pw_wfc->ggecut, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - } - // qianrui add 2021-8-13 to make different kpar parameters can get the same results + if (INPUT.pw_seed > 0) + { + MPI_Allreduce(MPI_IN_PLACE, &this->pw_wfc->ggecut, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + } + // qianrui add 2021-8-13 to make different kpar parameters can get the same results #endif this->pw_wfc->setuptransform(); @@ -315,14 +286,9 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) this->phsol = new hsolver::HSolverPW(this->pw_wfc, &this->wf); delete this->pelec; - this->pelec = new elecstate::ElecStatePW(this->pw_wfc, - &(this->chr), - (K_Vectors*)(&(this->kv)), - &ucell, - &GlobalC::ppcell, - this->pw_rhod, - this->pw_rho, - this->pw_big); + this->pelec + = new elecstate::ElecStatePW(this->pw_wfc, &(this->chr), (K_Vectors*)(&(this->kv)), &ucell, + &GlobalC::ppcell, this->pw_rhod, this->pw_rho, this->pw_big); this->pelec->charge->allocate(GlobalV::NSPIN); @@ -331,13 +297,8 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) delete this->pelec->pot; - this->pelec->pot = new elecstate::Potential(this->pw_rhod, - this->pw_rho, - &ucell, - &GlobalC::ppcell.vloc, - &(this->sf), - &(this->pelec->f_en.etxc), - &(this->pelec->f_en.vtxc)); + this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho, &ucell, &GlobalC::ppcell.vloc, + &(this->sf), &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc)); // temporary this->Init_GlobalC(inp, ucell, GlobalC::ppcell); @@ -347,28 +308,24 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) GlobalC::ppcell.init_vnl(ucell, this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); - this->pw_wfc->initgrids(ucell.lat0, - ucell.latvec, - this->pw_wfc->nx, - this->pw_wfc->ny, - this->pw_wfc->nz); + this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz); this->pw_wfc->initparameters(false, INPUT.ecutwfc, this->kv.get_nks(), this->kv.kvec_d.data()); this->pw_wfc->collect_local_pw(inp.erf_ecut, inp.erf_height, inp.erf_sigma); - - if(GlobalV::psi_initializer) // new initialization method, used in KSDFT and LCAO_IN_PW calculation + + if (GlobalV::psi_initializer) // new initialization method, used in KSDFT and LCAO_IN_PW calculation { // re-tabulate because GlobalV::DQ may change due to the change of atomic positions and cell parameters // for nao, we recalculate the overlap matrix between flz and jlq // for atomic, we recalculate the overlap matrix between pswfc and jlq - // for psig is not read-only, its value will be overwritten in initialize_psi(), dont need delete and reallocate - if((GlobalV::init_wfc.substr(0, 3) == "nao") - ||(GlobalV::init_wfc.substr(0, 6) == "atomic")) - { - this->psi_init->tabulate(); - } - } + // for psig is not read-only, its value will be overwritten in initialize_psi(), dont need delete and + // reallocate + if ((GlobalV::init_wfc.substr(0, 3) == "nao") || (GlobalV::init_wfc.substr(0, 6) == "atomic")) + { + this->psi_init->tabulate(); + } + } else // old initialization method, used in EXX calculation { this->wf.init_after_vc(this->kv.get_nks()); // reallocate wanf2, the planewave expansion of lcao @@ -377,18 +334,17 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) } #ifdef USE_PAW - if(GlobalV::use_paw) + if (GlobalV::use_paw) { - GlobalC::paw_cell.set_libpaw_ecut(INPUT.ecutwfc/2.0,INPUT.ecutwfc/2.0); //in Hartree - GlobalC::paw_cell.set_libpaw_fft(this->pw_wfc->nx,this->pw_wfc->ny,this->pw_wfc->nz, - this->pw_wfc->nx,this->pw_wfc->ny,this->pw_wfc->nz, - this->pw_wfc->startz,this->pw_wfc->numz); + GlobalC::paw_cell.set_libpaw_ecut(INPUT.ecutwfc / 2.0, INPUT.ecutwfc / 2.0); // in Hartree + GlobalC::paw_cell.set_libpaw_fft(this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz, this->pw_wfc->nx, + this->pw_wfc->ny, this->pw_wfc->nz, this->pw_wfc->startz, this->pw_wfc->numz); #ifdef __MPI - if(GlobalV::RANK_IN_POOL == 0) - { - GlobalC::paw_cell.prepare_paw(); - } + if (GlobalV::RANK_IN_POOL == 0) + { + GlobalC::paw_cell.prepare_paw(); + } #else GlobalC::paw_cell.prepare_paw(); #endif @@ -398,30 +354,24 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) std::vector> rhoijselect; std::vector nrhoijsel; #ifdef __MPI - if(GlobalV::RANK_IN_POOL == 0) + if (GlobalV::RANK_IN_POOL == 0) { GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - for(int iat = 0; iat < ucell.nat; iat ++) + for (int iat = 0; iat < ucell.nat; iat++) { - GlobalC::paw_cell.set_rhoij(iat, - nrhoijsel[iat], - rhoijselect[iat].size(), - rhoijselect[iat].data(), - rhoijp[iat].data()); - } + GlobalC::paw_cell.set_rhoij(iat, nrhoijsel[iat], rhoijselect[iat].size(), rhoijselect[iat].data(), + rhoijp[iat].data()); + } } #else GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - for(int iat = 0; iat < ucell.nat; iat ++) + for (int iat = 0; iat < ucell.nat; iat++) { - GlobalC::paw_cell.set_rhoij(iat, - nrhoijsel[iat], - rhoijselect[iat].size(), - rhoijselect[iat].data(), - rhoijp[iat].data()); - } + GlobalC::paw_cell.set_rhoij(iat, nrhoijsel[iat], rhoijselect[iat].size(), rhoijselect[iat].data(), + rhoijp[iat].data()); + } #endif } #endif @@ -429,7 +379,6 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) ModuleBase::timer::tick("ESolver_KS_PW", "init_after_vc"); } - template void ESolver_KS_PW::before_scf(const int istep) { @@ -446,9 +395,7 @@ void ESolver_KS_PW::before_scf(const int istep) #ifdef __MPI &(GlobalC::Pgrid), #endif - GlobalC::ucell, - this->pelec->charge, - &this->sf); + GlobalC::ucell, this->pelec->charge, &this->sf); } // init Hamilt, this should be allocated before each scf loop @@ -463,10 +410,7 @@ void ESolver_KS_PW::before_scf(const int istep) // allocate HamiltPW if (this->p_hamilt == nullptr) { - this->p_hamilt = new hamilt::HamiltPW( - this->pelec->pot, - this->pw_wfc, - &this->kv); + this->p_hamilt = new hamilt::HamiltPW(this->pelec->pot, this->pw_wfc, &this->kv); } //---------------------------------------------------------- @@ -481,77 +425,47 @@ void ESolver_KS_PW::before_scf(const int istep) // calculate ewald energy if (!GlobalV::test_skip_ewald) { - this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald( - GlobalC::ucell, - this->pw_rhod, - this->sf.strucFac); - } + this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(GlobalC::ucell, this->pw_rhod, this->sf.strucFac); + } //! cal_ux should be called before init_scf because //! the direction of ux is used in noncoline_rho - if(GlobalV::NSPIN == 4 && GlobalV::DOMAG) - { - GlobalC::ucell.cal_ux(); - } + if (GlobalV::NSPIN == 4 && GlobalV::DOMAG) + { + GlobalC::ucell.cal_ux(); + } //! calculate the total local pseudopotential in real space this->pelec->init_scf(istep, this->sf.strucFac); - - - if(GlobalV::out_chg == 2) + if (GlobalV::out_chg == 2) { - for(int is = 0; is < GlobalV::NSPIN; is++) + 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, - this->pw_big->nbz, - this->pw_rho->nplane, - this->pw_rho->startz_current, + this->pw_big->bz, this->pw_big->nbz, this->pw_rho->nplane, this->pw_rho->startz_current, #endif - this->pelec->charge->rho[is], - is, - GlobalV::NSPIN, - 0, - ss.str(), - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pelec->eferm.ef, - &(GlobalC::ucell), - 11); + this->pelec->charge->rho[is], is, GlobalV::NSPIN, 0, ss.str(), this->pw_rho->nx, this->pw_rho->ny, + this->pw_rho->nz, this->pelec->eferm.ef, &(GlobalC::ucell), 11); } } - ModuleIO::write_pot( - GlobalV::out_pot, - GlobalV::NSPIN, - GlobalV::global_out_dir, + ModuleIO::write_pot(GlobalV::out_pot, GlobalV::NSPIN, GlobalV::global_out_dir, #ifdef __MPI - this->pw_big->bz, - this->pw_big->nbz, - this->pw_rho->nplane, - this->pw_rho->startz_current, + this->pw_big->bz, this->pw_big->nbz, this->pw_rho->nplane, this->pw_rho->startz_current, #endif - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pelec->pot->get_effective_v()); + this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz, this->pelec->pot->get_effective_v()); //! Symmetry_rho should behind init_scf, because charge should be initialized first. //! liuyu comment: Symmetry_rho should be located between init_rho and v_of_rho? Symmetry_rho srho; for (int is = 0; is < GlobalV::NSPIN; is++) { - srho.begin(is, - *(this->pelec->charge), - this->pw_rhod, - GlobalC::Pgrid, - GlobalC::ucell.symm); - } + srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::Pgrid, GlobalC::ucell.symm); + } // liuyu move here 2023-10-09 // D in uspp need vloc, thus behind init_scf() @@ -562,22 +476,20 @@ void ESolver_KS_PW::before_scf(const int istep) // after init_rho (in pelec->init_scf), we have rho now. // before hamilt2density, we update Hk and initialize psi - if(GlobalV::psi_initializer) + if (GlobalV::psi_initializer) { // before_scf function will be called everytime before scf. However, once atomic coordinates changed, // structure factor will change, therefore all atomwise properties will change. So we need to reinitialize // psi every time before scf. But for random wavefunction, we dont, because random wavefunction is not // related to atomic coordinates. // What the old strategy does is only to initialize for once... - if(((GlobalV::init_wfc == "random")&&(istep == 0)) - ||(GlobalV::init_wfc != "random")) - { - this->initialize_psi(); - } + if (((GlobalV::init_wfc == "random") && (istep == 0)) || (GlobalV::init_wfc != "random")) + { + this->initialize_psi(); + } } } - template void ESolver_KS_PW::others(const int istep) { @@ -587,19 +499,14 @@ void ESolver_KS_PW::others(const int istep) if (cal_type == "test_memory") { - Cal_Test::test_memory(this->pw_rho, - this->pw_wfc, - this->p_chgmix->get_mixing_mode(), + Cal_Test::test_memory(this->pw_rho, this->pw_wfc, this->p_chgmix->get_mixing_mode(), this->p_chgmix->get_mixing_ndim()); } else if (cal_type == "gen_bessel") { Numerical_Descriptor nc; - nc.output_descriptor(this->psi[0], - INPUT.bessel_descriptor_lmax, - INPUT.bessel_descriptor_rcut, - INPUT.bessel_descriptor_tolerence, - this->kv.get_nks()); + nc.output_descriptor(this->psi[0], INPUT.bessel_descriptor_lmax, INPUT.bessel_descriptor_rcut, + INPUT.bessel_descriptor_tolerence, this->kv.get_nks()); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "GENERATE DESCRIPTOR FOR DEEPKS"); } else if (cal_type == "nscf") @@ -639,77 +546,56 @@ void ESolver_KS_PW::iter_init(const int istep, const int iter) } } - template void ESolver_KS_PW::allocate_psi_init() { // under restriction of C++11, std::unique_ptr can not be allocate via std::make_unique // use new instead, but will cause asymmetric allocation and deallocation, in literal aspect ModuleBase::timer::tick("ESolver_KS_PW", "allocate_psi_init"); - if((GlobalV::init_wfc.substr(0, 6) == "atomic")&&(GlobalC::ucell.natomwfc == 0)) - { - this->psi_init = std::unique_ptr>( - new psi_initializer_random()); - } - else if(GlobalV::init_wfc == "atomic") - { - this->psi_init = std::unique_ptr>( - new psi_initializer_atomic()); - } - else if(GlobalV::init_wfc == "random") - { - this->psi_init = std::unique_ptr>( - new psi_initializer_random()); - } - else if(GlobalV::init_wfc == "nao") - { - this->psi_init = std::unique_ptr>( - new psi_initializer_nao()); - } - else if(GlobalV::init_wfc == "atomic+random") - { - this->psi_init = std::unique_ptr>( - new psi_initializer_atomic_random()); - } - else if(GlobalV::init_wfc == "nao+random") - { - this->psi_init = std::unique_ptr>( - new psi_initializer_nao_random()); - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW::allocate_psi_init", - "for new psi initializer, init_wfc type not supported"); - } - - //! function polymorphism is moved from constructor to function initialize. + if ((GlobalV::init_wfc.substr(0, 6) == "atomic") && (GlobalC::ucell.natomwfc == 0)) + { + this->psi_init = std::unique_ptr>(new psi_initializer_random()); + } + else if (GlobalV::init_wfc == "atomic") + { + this->psi_init = std::unique_ptr>(new psi_initializer_atomic()); + } + else if (GlobalV::init_wfc == "random") + { + this->psi_init = std::unique_ptr>(new psi_initializer_random()); + } + else if (GlobalV::init_wfc == "nao") + { + this->psi_init = std::unique_ptr>(new psi_initializer_nao()); + } + else if (GlobalV::init_wfc == "atomic+random") + { + this->psi_init = std::unique_ptr>(new psi_initializer_atomic_random()); + } + else if (GlobalV::init_wfc == "nao+random") + { + this->psi_init = std::unique_ptr>(new psi_initializer_nao_random()); + } + else + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::allocate_psi_init", + "for new psi initializer, init_wfc type not supported"); + } + + //! function polymorphism is moved from constructor to function initialize. //! Two slightly different implementation are for MPI and serial case, respectively. #ifdef __MPI - this->psi_init->initialize( - &this->sf, - this->pw_wfc, - &GlobalC::ucell, - &GlobalC::Pkpoints, - 1, - &GlobalC::ppcell, - GlobalV::MY_RANK); + this->psi_init->initialize(&this->sf, this->pw_wfc, &GlobalC::ucell, &GlobalC::Pkpoints, 1, &GlobalC::ppcell, + GlobalV::MY_RANK); #else - this->psi_init->initialize( - &this->sf, - this->pw_wfc, - &GlobalC::ucell, - 1, - &GlobalC::ppcell); + this->psi_init->initialize(&this->sf, this->pw_wfc, &GlobalC::ucell, 1, &GlobalC::ppcell); #endif // always new->initialize->tabulate->allocate->proj_ao_onkG this->psi_init->tabulate(); ModuleBase::timer::tick("ESolver_KS_PW", "allocate_psi_init"); - } - - //! Although ESolver_KS_PW supports template, but in this function it has no relationship with //! heterogeneous calculation, so all templates function are specialized to double template @@ -736,15 +622,15 @@ void ESolver_KS_PW::initialize_psi(void) //! this way, we can avoid memory leak and undefined behavior std::weak_ptr> psig = this->psi_init->share_psig(); - if(psig.expired()) - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW::initialize_psi", "psig lifetime is expired"); - } + if (psig.expired()) + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::initialize_psi", "psig lifetime is expired"); + } //! to use psig, we need to lock it to get a shared pointer version, - //! then switch kpoint of psig to the given one - auto psig_ = psig.lock(); - psig_->fix_k(ik); + //! then switch kpoint of psig to the given one + auto psig_ = psig.lock(); + psig_->fix_k(ik); std::vector etatom(psig_->get_nbands(), 0.0); @@ -753,31 +639,21 @@ void ESolver_KS_PW::initialize_psi(void) if (this->psi_init->method() != "random") { // lcao_in_pw and pw share the same esolver. In the future, we will have different esolver - if ( - ((GlobalV::KS_SOLVER == "cg")||(GlobalV::KS_SOLVER == "lapack")) - &&(GlobalV::BASIS_TYPE == "pw") - ) + if (((GlobalV::KS_SOLVER == "cg") || (GlobalV::KS_SOLVER == "lapack")) && (GlobalV::BASIS_TYPE == "pw")) { // the following function is only run serially, to be improved - hsolver::DiagoIterAssist::diagH_subspace_init( - this->p_hamilt, - psig_->get_pointer(), - psig_->get_nbands(), - psig_->get_nbasis(), - *(this->kspw_psi), - etatom.data() - ); + hsolver::DiagoIterAssist::diagH_subspace_init(this->p_hamilt, psig_->get_pointer(), + psig_->get_nbands(), psig_->get_nbasis(), + *(this->kspw_psi), etatom.data()); continue; } - else if ((GlobalV::KS_SOLVER == "lapack") - && (GlobalV::BASIS_TYPE == "lcao_in_pw")) - { - if(ik == 0) - { - GlobalV::ofs_running - << " START WAVEFUNCTION: LCAO_IN_PW, psi initialization skipped " - << std::endl; - } + else if ((GlobalV::KS_SOLVER == "lapack") && (GlobalV::BASIS_TYPE == "lcao_in_pw")) + { + if (ik == 0) + { + GlobalV::ofs_running << " START WAVEFUNCTION: LCAO_IN_PW, psi initialization skipped " + << std::endl; + } continue; } // else the case is davidson @@ -786,12 +662,8 @@ void ESolver_KS_PW::initialize_psi(void) { if (GlobalV::KS_SOLVER == "cg") { - hsolver::DiagoIterAssist::diagH_subspace( - this->p_hamilt, - *(psig_), - *(this->kspw_psi), - etatom.data() - ); + hsolver::DiagoIterAssist::diagH_subspace(this->p_hamilt, *(psig_), *(this->kspw_psi), + etatom.data()); continue; } // else the case is davidson @@ -805,21 +677,17 @@ void ESolver_KS_PW::initialize_psi(void) (*(this->kspw_psi))(iband, ibasis) = (*psig_)(iband, ibasis); } } - }// end k-point loop + } // end k-point loop this->psi_init->set_initialized(true); - } // end GlobalV::psi_initializer + } // end GlobalV::psi_initializer ModuleBase::timer::tick("ESolver_KS_PW", "initialize_psi"); } - // Temporary, it should be replaced by hsolver later. template -void ESolver_KS_PW::hamilt2density( - const int istep, - const int iter, - const double ethr) +void ESolver_KS_PW::hamilt2density(const int istep, const int iter, const double ethr) { ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density"); @@ -831,14 +699,14 @@ void ESolver_KS_PW::hamilt2density( // choose if psi should be diag in subspace // be careful that istep start from 0 and iter start from 1 // if (iter == 1) - hsolver::DiagoIterAssist::need_subspace = ((istep == 0 || istep == 1) && iter == 1)? false : true; + hsolver::DiagoIterAssist::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true; hsolver::DiagoIterAssist::SCF_ITER = iter; hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; hsolver::DiagoIterAssist::PW_DIAG_NMAX = GlobalV::PW_DIAG_NMAX; // after init_rho (in pelec->init_scf), we have rho now. // before hamilt2density, we update Hk and initialize psi - if(GlobalV::psi_initializer) + if (GlobalV::psi_initializer) { // before_scf function will be called everytime before scf. However, once atomic coordinates changed, // structure factor will change, therefore all atomwise properties will change. So we need to reinitialize @@ -847,14 +715,13 @@ void ESolver_KS_PW::hamilt2density( // What the old strategy does is only to initialize for once... we also initialize only once here because // this can save a lot of time. But if cell and ion change significantly, re-initialization psi will be - // more efficient. Or an extrapolation strategy can be used. - if((istep == 0)&&(iter == 1) - &&!(this->psi_init->initialized())) - { - this->initialize_psi(); - } + // more efficient. Or an extrapolation strategy can be used. + if ((istep == 0) && (iter == 1) && !(this->psi_init->initialized())) + { + this->initialize_psi(); + } } - if(GlobalV::BASIS_TYPE != "lcao_in_pw") + if (GlobalV::BASIS_TYPE != "lcao_in_pw") { // from HSolverPW this->phsol->solve(this->p_hamilt, // hamilt::Hamilt* pHamilt, @@ -864,33 +731,33 @@ void ESolver_KS_PW::hamilt2density( } else { - // It is not a good choice to overload another solve function here, this will spoil the concept of + // It is not a good choice to overload another solve function here, this will spoil the concept of // multiple inheritance and polymorphism. But for now, we just do it in this way. // In the future, there will be a series of class ESolver_KS_LCAO_PW, HSolver_LCAO_PW and so on. std::weak_ptr> psig = this->psi_init->share_psig(); - if(psig.expired()) - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW::hamilt2density", "psig lifetime is expired"); - } + if (psig.expired()) + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::hamilt2density", "psig lifetime is expired"); + } - // from HSolverPW - this->phsol->solve(this->p_hamilt, // hamilt::Hamilt* pHamilt, - this->kspw_psi[0], // psi::Psi& psi, - this->pelec, // elecstate::ElecState* pelec, - psig.lock().get()[0]); // psi::Psi& transform, + // from HSolverPW + this->phsol->solve(this->p_hamilt, // hamilt::Hamilt* pHamilt, + this->kspw_psi[0], // psi::Psi& psi, + this->pelec, // elecstate::ElecState* pelec, + psig.lock().get()[0]); // psi::Psi& transform, } if (GlobalV::out_bandgap) { - if (!GlobalV::TWO_EFERMI) - { - this->pelec->cal_bandgap(); - } - else - { - this->pelec->cal_bandgap_updw(); - } - } + if (!GlobalV::TWO_EFERMI) + { + this->pelec->cal_bandgap(); + } + else + { + this->pelec->cal_bandgap_updw(); + } + } } else { @@ -910,19 +777,13 @@ void ESolver_KS_PW::hamilt2density( Symmetry_rho srho; for (int is = 0; is < GlobalV::NSPIN; is++) - { - srho.begin(is, - *(this->pelec->charge), - this->pw_rhod, - GlobalC::Pgrid, - GlobalC::ucell.symm); - } + { + srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::Pgrid, GlobalC::ucell.symm); + } // compute magnetization, only for LSDA(spin==2) - GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, - this->pelec->charge->nxyz, - this->pelec->charge->rho, - this->pelec->nelec_spin.data()); + GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, this->pelec->charge->nxyz, + this->pelec->charge->rho, this->pelec->nelec_spin.data()); // deband is calculated from "output" charge density calculated // in sum_band @@ -932,18 +793,17 @@ void ESolver_KS_PW::hamilt2density( ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density"); } - // Temporary, it should be rewritten with Hamilt class. template void ESolver_KS_PW::update_pot(const int istep, const int iter) { if (!this->conv_elec) { - if (GlobalV::NSPIN == 4) - { - GlobalC::ucell.cal_ux(); - } - this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + if (GlobalV::NSPIN == 4) + { + GlobalC::ucell.cal_ux(); + } + this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); } else @@ -952,7 +812,6 @@ void ESolver_KS_PW::update_pot(const int istep, const int iter) } } - template void ESolver_KS_PW::iter_finish(const int iter) { @@ -965,14 +824,13 @@ void ESolver_KS_PW::iter_finish(const int iter) GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, GlobalC::ucell); } - // 1 means Harris-Foulkes functional + // 1 means Harris-Foulkes functional // 2 means Kohn-Sham functional - const int energy_type = 2; + const int energy_type = 2; this->pelec->cal_energies(2); 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; } @@ -1003,7 +861,6 @@ void ESolver_KS_PW::iter_finish(const int iter) } } - template void ESolver_KS_PW::after_scf(const int istep) { @@ -1012,13 +869,11 @@ void ESolver_KS_PW::after_scf(const int istep) // save charge difference into files for charge extrapolation if (GlobalV::CALCULATION != "scf") { - this->CE.save_files(istep, - GlobalC::ucell, + this->CE.save_files(istep, GlobalC::ucell, #ifdef __MPI this->pw_big, #endif - this->pelec->charge, - &this->sf); + this->pelec->charge, &this->sf); } if (GlobalV::out_chg) @@ -1033,9 +888,8 @@ void ESolver_KS_PW::after_scf(const int istep) } } - if (this->wf.out_wfc_pw == 1 - || this->wf.out_wfc_pw == 2) - { + if (this->wf.out_wfc_pw == 1 || this->wf.out_wfc_pw == 2) + { std::stringstream ssw; ssw << GlobalV::global_out_dir << "WAVEFUNC"; ModuleIO::write_wfc_pw(ssw.str(), this->psi[0], this->kv, this->pw_wfc); @@ -1043,7 +897,7 @@ void ESolver_KS_PW::after_scf(const int istep) ModuleIO::output_convergence_after_scf(this->conv_elec, this->pelec->f_en.etot); - ModuleIO::output_efermi(this->conv_elec, this->pelec->eferm.ef); + ModuleIO::output_efermi(this->conv_elec, this->pelec->eferm.ef); if (GlobalV::OUT_LEVEL != "m") { @@ -1052,16 +906,14 @@ void ESolver_KS_PW::after_scf(const int istep) if (this->device == base_device::GpuDevice) { - castmem_2d_d2h_op()(this->psi[0].get_device(), - this->kspw_psi[0].get_device(), + castmem_2d_d2h_op()(this->psi[0].get_device(), this->kspw_psi[0].get_device(), this->psi[0].get_pointer() - this->psi[0].get_psi_bias(), - this->kspw_psi[0].get_pointer() - this->kspw_psi[0].get_psi_bias(), - this->psi[0].size()); + this->kspw_psi[0].get_pointer() - this->kspw_psi[0].get_psi_bias(), this->psi[0].size()); } // Get bands_to_print through public function of INPUT (returns a const pointer to string) std::string bands_to_print = *INPUT.get_bands_to_print(); - if(!bands_to_print.empty()) + if (!bands_to_print.empty()) { std::vector out_band_kb; Input_Conv::parse_expression(bands_to_print, out_band_kb); @@ -1109,7 +961,7 @@ void ESolver_KS_PW::after_scf(const int istep) if (!bands_picked[ib]) { continue; - } + } for (int i = 0; i < this->pw_rho->nxyz; i++) { @@ -1135,46 +987,32 @@ void ESolver_KS_PW::after_scf(const int istep) ModuleIO::write_rho( #ifdef __MPI - this->pw_big->bz, - this->pw_big->nbz, - this->pw_big->nplane, - this->pw_big->startz_current, + this->pw_big->bz, this->pw_big->nbz, this->pw_big->nplane, this->pw_big->startz_current, #endif - rho_band, - 0, - GlobalV::NSPIN, - 0, - ssc.str(), - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - 0.0, - &(GlobalC::ucell), - 11); + rho_band, 0, GlobalV::NSPIN, 0, ssc.str(), this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz, 0.0, + &(GlobalC::ucell), 11); } delete[] wfcr; delete[] rho_band; } } - template double ESolver_KS_PW::cal_energy() { return this->pelec->f_en.etot; } - template void ESolver_KS_PW::cal_force(ModuleBase::matrix& force) { Forces ff(GlobalC::ucell.nat); - if (this->__kspw_psi != nullptr) - { - this->__kspw_psi = nullptr; - } + if (this->__kspw_psi != nullptr) + { + this->__kspw_psi = nullptr; + } - if (this->__kspw_psi == nullptr) + if (this->__kspw_psi == nullptr) { this->__kspw_psi = GlobalV::precision_flag == "single" ? new psi::Psi, Device>(this->kspw_psi[0]) @@ -1182,17 +1020,10 @@ void ESolver_KS_PW::cal_force(ModuleBase::matrix& force) } //! Calculate forces - ff.cal_force(force, - *this->pelec, - this->pw_rhod, - &GlobalC::ucell.symm, - &this->sf, - &this->kv, - this->pw_wfc, - this->__kspw_psi); + ff.cal_force(force, *this->pelec, this->pw_rhod, &GlobalC::ucell.symm, &this->sf, &this->kv, this->pw_wfc, + this->__kspw_psi); } - template void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) { @@ -1204,24 +1035,17 @@ void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) if (this->__kspw_psi == nullptr) { - this->__kspw_psi = GlobalV::precision_flag == "single" - ? new psi::Psi, Device>(this->kspw_psi[0]) - : reinterpret_cast, Device>*>(this->kspw_psi); - } - ss.cal_stress(stress, - GlobalC::ucell, - this->pw_rhod, - &GlobalC::ucell.symm, - &this->sf, - &this->kv, - this->pw_wfc, - this->psi, - this->__kspw_psi); - - // external stress - double unit_transform = 0.0; - unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; - double external_stress[3] = {GlobalV::PRESS1, GlobalV::PRESS2, GlobalV::PRESS3}; + this->__kspw_psi = GlobalV::precision_flag == "single" + ? new psi::Psi, Device>(this->kspw_psi[0]) + : reinterpret_cast, Device>*>(this->kspw_psi); + } + ss.cal_stress(stress, GlobalC::ucell, this->pw_rhod, &GlobalC::ucell.symm, &this->sf, &this->kv, this->pw_wfc, + this->psi, this->__kspw_psi); + + // external stress + double unit_transform = 0.0; + unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; + double external_stress[3] = {GlobalV::PRESS1, GlobalV::PRESS2, GlobalV::PRESS3}; for (int i = 0; i < 3; i++) { stress(i, i) -= external_stress[i] / unit_transform; @@ -1229,7 +1053,6 @@ void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) GlobalV::PRESSURE = (stress(0, 0) + stress(1, 1) + stress(2, 2)) / 3; } - template void ESolver_KS_PW::after_all_runners(void) { @@ -1247,28 +1070,24 @@ void ESolver_KS_PW::after_all_runners(void) 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"; - } - - int nspin0 = 1; - if (GlobalV::NSPIN == 2) - { - nspin0 = 2; - } - //! print occupation in istate.info - ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv, &(GlobalC::Pkpoints)); - - //! compute density of states + GlobalV::ofs_running << " | done here. |" << std::endl; + GlobalV::ofs_running << " | |" << std::endl; + GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + GlobalV::ofs_running << "\n\n\n\n"; + } + + int nspin0 = 1; + if (GlobalV::NSPIN == 2) + { + nspin0 = 2; + } + //! print occupation in istate.info + ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv, &(GlobalC::Pkpoints)); + + //! compute density of states if (INPUT.out_dos) { - ModuleIO::write_dos_pw(this->pelec->ekb, - this->pelec->wg, - this->kv, - INPUT.dos_edelta_ev, - INPUT.dos_scale, + ModuleIO::write_dos_pw(this->pelec->ekb, this->pelec->wg, this->kv, INPUT.dos_edelta_ev, INPUT.dos_scale, INPUT.dos_sigma); if (nspin0 == 1) @@ -1276,19 +1095,14 @@ void ESolver_KS_PW::after_all_runners(void) GlobalV::ofs_running << " Fermi energy is " << this->pelec->eferm.ef << " Rydberg" << std::endl; } else if (nspin0 == 2) - { - GlobalV::ofs_running << " Fermi energy (spin = 1) is " - << this->pelec->eferm.ef_up - << " Rydberg" - << std::endl; - GlobalV::ofs_running << " Fermi energy (spin = 2) is " - << this->pelec->eferm.ef_dw - << " Rydberg" - << std::endl; - } + { + GlobalV::ofs_running << " Fermi energy (spin = 1) is " << this->pelec->eferm.ef_up << " Rydberg" + << std::endl; + GlobalV::ofs_running << " Fermi energy (spin = 2) is " << this->pelec->eferm.ef_dw << " Rydberg" + << std::endl; + } } - if (INPUT.out_band[0]) // pengfei 2014-10-13 { for (int is = 0; is < nspin0; is++) @@ -1296,18 +1110,11 @@ void ESolver_KS_PW::after_all_runners(void) std::stringstream ss2; ss2 << GlobalV::global_out_dir << "BANDS_" << is + 1 << ".dat"; GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() << std::endl; - ModuleIO::nscf_band(is, - ss2.str(), - GlobalV::NBANDS, - 0.0, - INPUT.out_band[1], - this->pelec->ekb, - this->kv, + ModuleIO::nscf_band(is, ss2.str(), GlobalV::NBANDS, 0.0, INPUT.out_band[1], this->pelec->ekb, this->kv, &(GlobalC::Pkpoints)); } } - if (GlobalV::BASIS_TYPE == "pw" && winput::out_spillage) // xiaohui add 2013-09-01 { // calculate spillage value. @@ -1315,12 +1122,12 @@ void ESolver_KS_PW::after_all_runners(void) // ! Print out overlap before spillage optimization to generate atomic orbitals if (winput::out_spillage <= 2) { - for(int i = 0; i < INPUT.bessel_nao_rcuts.size(); i++) + for (int i = 0; i < INPUT.bessel_nao_rcuts.size(); i++) { - if(GlobalV::MY_RANK == 0) + if (GlobalV::MY_RANK == 0) { - std::cout << "update value: bessel_nao_rcut <- " - << std::fixed << INPUT.bessel_nao_rcuts[i] << " a.u." << std::endl; + std::cout << "update value: bessel_nao_rcut <- " << std::fixed << INPUT.bessel_nao_rcuts[i] + << " a.u." << std::endl; } INPUT.bessel_nao_rcut = INPUT.bessel_nao_rcuts[i]; Numerical_Basis numerical_basis; @@ -1333,25 +1140,16 @@ void ESolver_KS_PW::after_all_runners(void) //! Print out wave functions in real space if (this->wf.out_wfc_r == 1) // Peize Lin add 2021.11.21 { - ModuleIO::write_psi_r_1( - this->psi[0], - this->pw_wfc, - "wfc_realspace", - true, - this->kv); + ModuleIO::write_psi_r_1(this->psi[0], this->pw_wfc, "wfc_realspace", true, this->kv); } //! Use Kubo-Greenwood method to compute conductivities if (INPUT.cal_cond) { - this->KG( - INPUT.cond_smear, - INPUT.cond_fwhm, - INPUT.cond_wcut, - INPUT.cond_dw, - INPUT.cond_dt, - this->pelec->wg); - } + EleCond elec_cond(&GlobalC::ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &GlobalC::ppcell); + elec_cond.KG(INPUT.cond_smear, INPUT.cond_fwhm, INPUT.cond_wcut, INPUT.cond_dw, INPUT.cond_dt, + INPUT.cond_nonlocal, this->pelec->wg); + } } template @@ -1402,32 +1200,24 @@ void ESolver_KS_PW::nscf(void) if (nspin == 2) { if (ik == 0) - { - GlobalV::ofs_running << " spin up :" << std::endl; - } - if (ik == (this->kv.get_nks() / 2)) - { - GlobalV::ofs_running << " spin down :" << std::endl; - } + { + GlobalV::ofs_running << " spin up :" << std::endl; + } + 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_band " - << 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 << " spin" << this->kv.isk[ik] + 1 << "_final_band " << 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; } @@ -1437,58 +1227,33 @@ void ESolver_KS_PW::nscf(void) if (!GlobalV::TWO_EFERMI) { this->pelec->cal_bandgap(); - GlobalV::ofs_running << " E_bandgap " - << this->pelec->bandgap * ModuleBase::Ry_to_eV - << " eV" << std::endl; - } + 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; - } + 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; + } } //! 6) calculate Wannier functions if (INPUT.towannier90) { - toWannier90_PW wan( - INPUT.out_wannier_mmn, - INPUT.out_wannier_amn, - INPUT.out_wannier_unk, - INPUT.out_wannier_eig, - INPUT.out_wannier_wvfn_formatted, - INPUT.nnkpfile, - INPUT.wannier_spin - ); - - wan.calculate( - this->pelec->ekb, - this->pw_wfc, - this->pw_big, - this->kv, - this->psi); - } - - - //! 7) calculate Berry phase polarization - if (berryphase::berry_phase_flag - && ModuleSymmetry::Symmetry::symm_flag != 1) + toWannier90_PW wan(INPUT.out_wannier_mmn, INPUT.out_wannier_amn, INPUT.out_wannier_unk, INPUT.out_wannier_eig, + INPUT.out_wannier_wvfn_formatted, INPUT.nnkpfile, INPUT.wannier_spin); + + wan.calculate(this->pelec->ekb, this->pw_wfc, this->pw_big, this->kv, this->psi); + } + + //! 7) calculate Berry phase polarization + if (berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1) { berryphase bp; - bp.Macroscopic_polarization( - this->pw_wfc->npwk_max, - this->psi, - this->pw_rho, - this->pw_wfc, - this->kv); - } + bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); + } ModuleBase::timer::tick("ESolver_KS_PW", "nscf"); return; diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index 235d242117..8a65dc78e7 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -3,145 +3,86 @@ #include "./esolver_ks.h" #include "module_hamilt_pw/hamilt_pwdft/operator_pw/velocity_pw.h" #include "module_psi/psi_initializer.h" + #include #include - namespace ModuleESolver { template class ESolver_KS_PW : public ESolver_KS { - private: - - using Real = typename GetTypeReal::type; - - public: - - ESolver_KS_PW(); + private: + using Real = typename GetTypeReal::type; - ~ESolver_KS_PW(); + public: + ESolver_KS_PW(); - void before_all_runners(Input& inp, UnitCell& cell) override; + ~ESolver_KS_PW(); - void init_after_vc(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; - double cal_energy() override; + void init_after_vc(Input& inp, UnitCell& cell) override; - void cal_force(ModuleBase::matrix& force) override; + double cal_energy() override; - void cal_stress(ModuleBase::matrix& stress) override; + void cal_force(ModuleBase::matrix& force) override; - virtual void hamilt2density(const int istep, const int iter, const double ethr) override; + void cal_stress(ModuleBase::matrix& stress) override; - virtual void hamilt2estates(const double ethr) override; + virtual void hamilt2density(const int istep, const int iter, const double ethr) override; - virtual void nscf() override; + virtual void hamilt2estates(const double ethr) override; - void after_all_runners() override; + virtual void nscf() override; - /** - * @brief calculate Onsager coefficients Lmn(\omega) and conductivities with Kubo-Greenwood formula - * - * @param fwhmin FWHM for delta function - * @param smear_type 1: Gaussian, 2: Lorentzian - * @param wcut cutoff \omega for Lmn(\omega) - * @param dw_in \omega step - * @param dt_in time step - * @param wg wg(ik,ib) occupation for the ib-th band in the ik-th kpoint - */ - void KG(const int& smear_type, - const double fwhmin, - const double wcut, - const double dw_in, - const double dt_in, - ModuleBase::matrix& wg); + void after_all_runners() override; - /** - * @brief calculate the response function Cmn(t) for currents - * - * @param ik k point - * @param nt number of steps of time - * @param dt time step - * @param decut ignore dE which is larger than decut - * @param wg wg(ik,ib) occupation for the ib-th band in the ik-th kpoint - * @param velop velocity operator - * @param ct11 C11(t) - * @param ct12 C12(t) - * @param ct22 C22(t) - */ - void jjcorr_ks(const int ik, - const int nt, - const double dt, - const double decut, - ModuleBase::matrix& wg, - hamilt::Velocity& velop, - double* ct11, - double* ct12, - double* ct22); + protected: + virtual void before_scf(const int istep) override; - protected: + virtual void iter_init(const int istep, const int iter) override; - virtual void before_scf(const int istep) override; + virtual void update_pot(const int istep, const int iter) override; - virtual void iter_init(const int istep, const int iter) override; + virtual void iter_finish(const int iter) override; - virtual void update_pot(const int istep, const int iter) override; + virtual void after_scf(const int istep) override; - virtual void iter_finish(const int iter) override; + virtual void others(const int istep) override; - virtual void after_scf(const int istep) override; + // temporary, this will be removed in the future; + // Init Global class + void Init_GlobalC(Input& inp, UnitCell& ucell, pseudopot_cell_vnl& ppcell); - virtual void others(const int istep)override; + /// @brief allocate psi_init the new psi_initializer + void allocate_psi_init(); - //temporary, this will be removed in the future; - //Init Global class - void Init_GlobalC( - Input &inp, - UnitCell &ucell, - pseudopot_cell_vnl &ppcell); + /// @brief initialize psi + void initialize_psi(); - /// @brief calculate conductivities from j-j correlation function - void calcondw(const int nt, - const double dt, - const int& smear_type, - const double fwhmin, - const double wcut, - const double dw_in, - double* ct11, - double* ct12, - double* ct22); + protected: + //! hide the psi in ESolver_KS for tmp use + psi::Psi, base_device::DEVICE_CPU>* psi = nullptr; - /// @brief allocate psi_init the new psi_initializer - void allocate_psi_init(); + private: + // psi_initializer* psi_init = nullptr; + // change to use smart pointer to manage the memory, and avoid memory leak + // while the std::make_unique() is not supported till C++14, + // so use the new and std::unique_ptr to manage the memory, but this makes new-delete not symmetric + std::unique_ptr> psi_init; - /// @brief initialize psi - void initialize_psi(); + Device* ctx = {}; - protected: + base_device::AbacusDevice_t device = {}; - //! hide the psi in ESolver_KS for tmp use - psi::Psi, base_device::DEVICE_CPU>* psi = nullptr; + psi::Psi* kspw_psi = nullptr; - private: + psi::Psi, Device>* __kspw_psi = nullptr; - // psi_initializer* psi_init = nullptr; - // change to use smart pointer to manage the memory, and avoid memory leak - // while the std::make_unique() is not supported till C++14, - // so use the new and std::unique_ptr to manage the memory, but this makes new-delete not symmetric - std::unique_ptr> psi_init; - - Device * ctx = {}; - - base_device::AbacusDevice_t device = {}; - - psi::Psi* kspw_psi = nullptr; - - psi::Psi, Device>* __kspw_psi = nullptr; - - using castmem_2d_d2h_op - = base_device::memory::cast_memory_op, T, base_device::DEVICE_CPU, Device>; - }; -} // namespace ModuleESolver + using castmem_2d_d2h_op + = base_device::memory::cast_memory_op, T, base_device::DEVICE_CPU, Device>; +}; +} // namespace ModuleESolver #endif diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index 175a40186c..7901f22568 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -1,17 +1,19 @@ #include "esolver_sdft_pw.h" -#include -#include - #include "module_base/memory.h" #include "module_base/timer.h" #include "module_elecstate/elecstate_pw_sdft.h" +#include "module_hamilt_pw/hamilt_stodft/sto_dos.h" +#include "module_hamilt_pw/hamilt_stodft/sto_elecond.h" #include "module_hsolver/diago_iter_assist.h" #include "module_hsolver/hsolver_pw_sdft.h" #include "module_io/output_log.h" #include "module_io/rho_io.h" #include "module_io/write_istate_info.h" +#include +#include + //-------------------Temporary------------------ #include "module_base/global_variable.h" #include "module_elecstate/module_charge/symmetry_rho.h" @@ -46,14 +48,8 @@ void ESolver_SDFT_PW::before_all_runners(Input& inp, UnitCell& ucell) ESolver_KS::before_all_runners(inp, ucell); // 3) initialize the pointer for electronic states of SDFT - this->pelec = new elecstate::ElecStatePW_SDFT(pw_wfc, - &(chr), - (K_Vectors*)(&(kv)), - &ucell, - &(GlobalC::ppcell), - this->pw_rhod, - this->pw_rho, - pw_big); + this->pelec = new elecstate::ElecStatePW_SDFT(pw_wfc, &(chr), (K_Vectors*)(&(kv)), &ucell, &(GlobalC::ppcell), + this->pw_rhod, this->pw_rho, pw_big); // 4) inititlize the charge density. this->pelec->charge->allocate(GlobalV::NSPIN); @@ -62,13 +58,8 @@ void ESolver_SDFT_PW::before_all_runners(Input& inp, UnitCell& ucell) // 5) initialize the potential. if (this->pelec->pot == nullptr) { - this->pelec->pot = new elecstate::Potential(pw_rhod, - pw_rho, - &ucell, - &(GlobalC::ppcell.vloc), - &(sf), - &(this->pelec->f_en.etxc), - &(this->pelec->f_en.vtxc)); + this->pelec->pot = new elecstate::Potential(pw_rhod, pw_rho, &ucell, &(GlobalC::ppcell.vloc), &(sf), + &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc)); GlobalTemp::veff = &(this->pelec->pot->get_effective_v()); } @@ -95,10 +86,10 @@ void ESolver_SDFT_PW::before_all_runners(Input& inp, UnitCell& ucell) Init_Sto_Orbitals_Ecut(this->stowf, inp.seed_sto, kv, *pw_wfc, inp.initsto_ecut); } } - else - { - Init_Com_Orbitals(this->stowf); - } + else + { + Init_Com_Orbitals(this->stowf); + } size_t size = stowf.chi0->size(); @@ -108,7 +99,8 @@ void ESolver_SDFT_PW::before_all_runners(Input& inp, UnitCell& ucell) if (GlobalV::NBANDS > 0) { - this->stowf.chiortho = new psi::Psi>(kv.get_nks(), stowf.nchip_max, wf.npwx, kv.ngk.data()); + this->stowf.chiortho + = new psi::Psi>(kv.get_nks(), stowf.nchip_max, wf.npwx, kv.ngk.data()); ModuleBase::Memory::record("SDFT::chiortho", size * sizeof(std::complex)); } @@ -118,14 +110,13 @@ void ESolver_SDFT_PW::before_all_runners(Input& inp, UnitCell& ucell) return; } - void ESolver_SDFT_PW::before_scf(const int istep) { ESolver_KS_PW::before_scf(istep); - if (istep > 0 && INPUT.nbands_sto != 0 && INPUT.initsto_freq > 0 && istep % INPUT.initsto_freq == 0) - { - Update_Sto_Orbitals(this->stowf, INPUT.seed_sto); - } + if (istep > 0 && INPUT.nbands_sto != 0 && INPUT.initsto_freq > 0 && istep % INPUT.initsto_freq == 0) + { + Update_Sto_Orbitals(this->stowf, INPUT.seed_sto); + } } void ESolver_SDFT_PW::iter_finish(int iter) @@ -134,19 +125,16 @@ void ESolver_SDFT_PW::iter_finish(int iter) this->pelec->cal_energies(2); } - void ESolver_SDFT_PW::after_scf(const int istep) { // save charge difference into files for charge extrapolation if (GlobalV::CALCULATION != "scf") { - this->CE.save_files(istep, - GlobalC::ucell, + this->CE.save_files(istep, GlobalC::ucell, #ifdef __MPI this->pw_big, #endif - this->pelec->charge, - &this->sf); + this->pelec->charge, &this->sf); } if (GlobalV::out_chg > 0) @@ -158,21 +146,10 @@ void ESolver_SDFT_PW::after_scf(const int istep) const double ef_tmp = this->pelec->eferm.get_efval(is); ModuleIO::write_rho( #ifdef __MPI - pw_big->bz, - pw_big->nbz, - pw_rho->nplane, - pw_rho->startz_current, + pw_big->bz, pw_big->nbz, pw_rho->nplane, pw_rho->startz_current, #endif - pelec->charge->rho_save[is], - is, - GlobalV::NSPIN, - 0, - ssc.str(), - pw_rho->nx, - pw_rho->ny, - pw_rho->nz, - ef_tmp, - &(GlobalC::ucell)); + pelec->charge->rho_save[is], is, GlobalV::NSPIN, 0, ssc.str(), pw_rho->nx, pw_rho->ny, pw_rho->nz, + ef_tmp, &(GlobalC::ucell)); } } @@ -199,15 +176,7 @@ void ESolver_SDFT_PW::hamilt2density(int istep, int iter, double ethr) hsolver::DiagoIterAssist>::PW_DIAG_NMAX = GlobalV::PW_DIAG_NMAX; - this->phsol->solve( - this->p_hamilt, - this->psi[0], - this->pelec, - pw_wfc, - this->stowf, - istep, - iter, - GlobalV::KS_SOLVER); + this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec, pw_wfc, this->stowf, istep, iter, GlobalV::KS_SOLVER); if (GlobalV::MY_STOGROUP == 0) { @@ -221,10 +190,10 @@ void ESolver_SDFT_PW::hamilt2density(int istep, int iter, double ethr) else { #ifdef __MPI - if (ModuleSymmetry::Symmetry::symm_flag == 1) - { - MPI_Barrier(MPI_COMM_WORLD); - } + if (ModuleSymmetry::Symmetry::symm_flag == 1) + { + MPI_Barrier(MPI_COMM_WORLD); + } #endif } } @@ -234,40 +203,20 @@ double ESolver_SDFT_PW::cal_energy() return this->pelec->f_en.etot; } - void ESolver_SDFT_PW::cal_force(ModuleBase::matrix& force) { Sto_Forces ff(GlobalC::ucell.nat); - ff.cal_stoforce(force, - *this->pelec, - pw_rho, - &GlobalC::ucell.symm, - &sf, - &kv, - pw_wfc, - this->psi, - this->stowf); + ff.cal_stoforce(force, *this->pelec, pw_rho, &GlobalC::ucell.symm, &sf, &kv, pw_wfc, this->psi, this->stowf); } - void ESolver_SDFT_PW::cal_stress(ModuleBase::matrix& stress) { Sto_Stress_PW ss; - ss.cal_stress( - stress, - *this->pelec, - pw_rho, - &GlobalC::ucell.symm, - &sf, - &kv, - pw_wfc, - this->psi, - this->stowf, - pelec->charge); + ss.cal_stress(stress, *this->pelec, pw_rho, &GlobalC::ucell.symm, &sf, &kv, pw_wfc, this->psi, this->stowf, + pelec->charge); } - void ESolver_SDFT_PW::after_all_runners(void) { @@ -279,74 +228,27 @@ void ESolver_SDFT_PW::after_all_runners(void) ((hsolver::HSolverPW_SDFT*)phsol)->stoiter.cleanchiallorder(); // release lots of memories - int nche_test = this->nche_sto; - - if (INPUT.out_dos) - { - nche_test = std::max(nche_test, INPUT.dos_nche); - } - - int cond_nche = 0; - if (INPUT.cal_cond) - { - cond_nche = set_cond_nche(INPUT.cond_dt, INPUT.cond_dtbatch, 1e-8, nche_test, INPUT.emin_sto, INPUT.emax_sto); - } - else - { - check_che(nche_test, INPUT.emin_sto, INPUT.emax_sto); - } - if (INPUT.out_dos) { - double emax=0.0; - double emin=0.0; - if (INPUT.dos_setemax) - { - emax = INPUT.dos_emax_ev; - } - else - { - emax = ((hsolver::HSolverPW_SDFT*)phsol)->stoiter.stohchi.Emax * ModuleBase::Ry_to_eV; - } - if (INPUT.dos_setemin) - { - emin = INPUT.dos_emin_ev; - } - else - { - emin = ((hsolver::HSolverPW_SDFT*)phsol)->stoiter.stohchi.Emin * ModuleBase::Ry_to_eV; - } - - if (!INPUT.dos_setemax && !INPUT.dos_setemin) - { - double delta = (emax - emin) * INPUT.dos_scale; - emax = emax + delta / 2.0; - emin = emin - delta / 2.0; - } - this->caldos( - INPUT.dos_nche, - INPUT.dos_sigma, - emin, - emax, - INPUT.dos_edelta_ev, - INPUT.npart_sto); + Sto_DOS sto_dos(this->pw_wfc, &this->kv, this->pelec, this->psi, this->p_hamilt, + (hsolver::HSolverPW_SDFT*)phsol, &stowf); + sto_dos.decide_param(INPUT.dos_nche, INPUT.emin_sto, INPUT.emax_sto, INPUT.dos_setemin, INPUT.dos_setemax, + INPUT.dos_emin_ev, INPUT.dos_emax_ev, INPUT.dos_scale); + sto_dos.caldos(INPUT.dos_sigma, INPUT.dos_edelta_ev, INPUT.npart_sto); } // sKG cost memory, and it should be placed at the end of the program if (INPUT.cal_cond) { - this->sKG(cond_nche, - INPUT.cond_smear, - INPUT.cond_fwhm, - INPUT.cond_wcut, - INPUT.cond_dw, - INPUT.cond_dt, - INPUT.cond_dtbatch, - INPUT.npart_sto); + Sto_EleCond sto_elecond(&GlobalC::ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &GlobalC::ppcell, + this->p_hamilt, (hsolver::HSolverPW_SDFT*)phsol, &stowf); + sto_elecond.decide_nche(INPUT.cond_dt, INPUT.cond_dtbatch, 1e-8, this->nche_sto, INPUT.emin_sto, + INPUT.emax_sto); + sto_elecond.sKG(INPUT.cond_smear, INPUT.cond_fwhm, INPUT.cond_wcut, INPUT.cond_dw, INPUT.cond_dt, + INPUT.cond_nonlocal, INPUT.cond_dtbatch, INPUT.npart_sto); } } - void ESolver_SDFT_PW::others(const int istep) { ModuleBase::TITLE("ESolver_SDFT_PW", "others"); @@ -363,7 +265,6 @@ void ESolver_SDFT_PW::others(const int istep) return; } - void ESolver_SDFT_PW::nscf(void) { ModuleBase::TITLE("ESolver_SDFT_PW", "nscf"); diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h index c5692a7303..f63b3c0a05 100644 --- a/source/module_esolver/esolver_sdft_pw.h +++ b/source/module_esolver/esolver_sdft_pw.h @@ -2,7 +2,6 @@ #define ESOLVER_SDFT_PW_H #include "esolver_ks_pw.h" -#include "module_hamilt_pw/hamilt_pwdft/operator_pw/velocity_pw.h" #include "module_hamilt_pw/hamilt_stodft/sto_hchi.h" #include "module_hamilt_pw/hamilt_stodft/sto_iter.h" #include "module_hamilt_pw/hamilt_stodft/sto_wf.h" @@ -25,7 +24,6 @@ class ESolver_SDFT_PW : public ESolver_KS_PW> void cal_stress(ModuleBase::matrix& stress) override; public: - Stochastic_WF stowf; protected: @@ -43,94 +41,9 @@ class ESolver_SDFT_PW : public ESolver_KS_PW> virtual void after_all_runners() override; - public: - /** - * @brief calculate Stochastic Kubo-Greenwood - * - * @param nche_KG Number Chebyshev orders - * @param fwhmin FWHM - * @param smear_type 1: Gaussian, 2: Lorentzian - * @param wcut cutoff omega - * @param dw_in omega step - * @param dt_in t step - * @param nbatch t step batch - * @param npart_sto number stochastic wavefunctions parts to evalution simultaneously - */ - void sKG(const int nche_KG, - const int& smear_type, - const double fwhmin, - const double wcut, - const double dw_in, - const double dt_in, - const int nbatch, - const int npart_sto); - // calculate DOS - void caldos(const int nche_dos, - const double sigmain, - const double emin, - const double emax, - const double de, - const int npart); - private: int nche_sto; ///< norder of Chebyshev int method_sto; ///< method of SDFT - - /** - * @brief Check if Emin and Emax are converged - * - * @param nche_in N order of Chebyshev expansion - * @param try_emin trial Emin - * @param try_emax trial Emax - */ - void check_che(const int nche_in, const double try_emin, const double try_emax); - - /** - * @brief Set the N order of Chebyshev expansion for conductivities - * - * @param dt t step - * @param nbatch number of t batch - * @param cond_thr threshold of errors for conductivities - * @param nche_min minimum N order of Chebyshev - * @param try_emin trial Emin - * @param try_emax trial Emax - * @return N order of Chebyshev - */ - int set_cond_nche(const double dt, - int& nbatch, - const double cond_thr, - const int& nche_min, - double try_emin, - double try_emax); - - /** - * @brief calculate Jmatrix - * - */ - void cal_jmatrix(const psi::Psi>& kspsi_all, - const psi::Psi>& vkspsi, - const double* en, - const double* en_all, - std::complex* leftfact, - std::complex* rightfact, - const psi::Psi>& leftchi, - psi::Psi>& rightchi, - psi::Psi>& left_hchi, - psi::Psi>& batch_vchi, - psi::Psi>& batch_vhchi, -#ifdef __MPI - psi::Psi>& chi_all, - psi::Psi>& hchi_all, - void* gatherinfo_ks, - void* gatherinfo_sto, -#endif - const int& bsize_psi, - std::vector>& j1, - std::vector>& j2, - hamilt::Velocity& velop, - const int& ik, - const std::complex& factor, - const int bandinfo[6]); }; } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_sdft_pw_tool.cpp b/source/module_esolver/esolver_sdft_pw_tool.cpp deleted file mode 100644 index 257dd7f2e5..0000000000 --- a/source/module_esolver/esolver_sdft_pw_tool.cpp +++ /dev/null @@ -1,1433 +0,0 @@ -#include - -#include "esolver_sdft_pw.h" -#include "module_base/complexmatrix.h" -#include "module_base/constants.h" -#include "module_base/global_function.h" -#include "module_base/global_variable.h" -#include "module_base/memory.h" -#include "module_base/timer.h" -#include "module_base/vector3.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_hsolver/hsolver_pw_sdft.h" - -#define TWOSQRT2LN2 2.354820045030949 // FWHM = 2sqrt(2ln2) * \sigma -#define FACTOR 1.839939223835727e7 - -namespace ModuleESolver -{ -struct parallel_distribution -{ - parallel_distribution(const int& num_all, const int& np, const int myrank) - { - int num_per = num_all / np; - int st_per = num_per * myrank; - int re = num_all % np; - if (myrank < re) - { - ++num_per; - st_per += myrank; - } - else - { - st_per += re; - } - this->start = st_per; - this->num_per = num_per; - } - int start; - int num_per; -}; - -#ifdef __MPI -struct info_gatherv -{ - info_gatherv(const int& ngroup_per, const int& np, const int& num_in_group, MPI_Comm comm_world) - { - nrecv = new int[np]; - displs = new int[np]; - MPI_Allgather(&ngroup_per, 1, MPI_INT, nrecv, 1, MPI_INT, comm_world); - displs[0] = 0; - for (int i = 1; i < np; ++i) - { - displs[i] = displs[i - 1] + nrecv[i - 1]; - } - for (int i = 0; i < np; ++i) - { - nrecv[i] *= num_in_group; - displs[i] *= num_in_group; - } - } - ~info_gatherv() - { - delete[] nrecv; - delete[] displs; - } - int* nrecv = nullptr; - int* displs = nullptr; -}; -#endif - - -void convert_psi(const psi::Psi>& psi_in, psi::Psi>& psi_out) -{ - psi_in.fix_k(0); - psi_out.fix_k(0); - for (int i = 0; i < psi_in.size(); ++i) - { - psi_out.get_pointer()[i] = static_cast>(psi_in.get_pointer()[i]); - } - return; -} - - -psi::Psi>* gatherchi(psi::Psi>& chi, - psi::Psi>& chi_all, - const int& npwx, - int* nrecv_sto, - int* displs_sto, - const int perbands_sto) -{ - psi::Psi>* p_chi; - p_chi = χ -#ifdef __MPI - if (GlobalV::NSTOGROUP > 1) - { - p_chi = &chi_all; - ModuleBase::timer::tick("sKG", "bands_gather"); - MPI_Allgatherv(chi.get_pointer(), - perbands_sto * npwx, - MPI_COMPLEX, - chi_all.get_pointer(), - nrecv_sto, - displs_sto, - MPI_COMPLEX, - PARAPW_WORLD); - ModuleBase::timer::tick("sKG", "bands_gather"); - } -#endif - return p_chi; -} - - -void ESolver_SDFT_PW::check_che(const int nche_in, const double try_emin, const double try_emax) -{ - //------------------------------ - // Convergence test - //------------------------------ - bool change = false; - const int nk = kv.get_nks(); - ModuleBase::Chebyshev chetest(nche_in); - Stochastic_Iter& stoiter = ((hsolver::HSolverPW_SDFT*)phsol)->stoiter; - Stochastic_hchi& stohchi = stoiter.stohchi; - int ntest0 = 5; - stohchi.Emax = try_emax; - stohchi.Emin = try_emin; - // if (GlobalV::NBANDS > 0) - // { - // double tmpemin = 1e10; - // for (int ik = 0; ik < nk; ++ik) - // { - // tmpemin = std::min(tmpemin, this->pelec->ekb(ik, GlobalV::NBANDS - 1)); - // } - // stohchi.Emin = tmpemin; - // } - // else - // { - // stohchi.Emin = 0; - // } - for (int ik = 0; ik < nk; ++ik) - { - this->p_hamilt->updateHk(ik); - stoiter.stohchi.current_ik = ik; - const int npw = kv.ngk[ik]; - std::complex* pchi = nullptr; - int ntest = std::min(ntest0, stowf.nchip[ik]); - for (int i = 0; i < ntest; ++i) - { - if (INPUT.nbands_sto == 0) - { - pchi = new std::complex[npw]; - for (int ig = 0; ig < npw; ++ig) - { - double rr = std::rand() / double(RAND_MAX); - double arg = std::rand() / double(RAND_MAX); - pchi[ig] = std::complex(rr * cos(arg), rr * sin(arg)); - } - } - else if (GlobalV::NBANDS > 0) - { - pchi = &stowf.chiortho[0](ik, i, 0); - } - else - { - pchi = &stowf.chi0[0](ik, i, 0); - } - while (1) - { - bool converge; - converge = chetest.checkconverge(&stohchi, - &Stochastic_hchi::hchi_norm, - pchi, - npw, - stohchi.Emax, - stohchi.Emin, - 2.0); - - if (!converge) - { - change = true; - } - else - { - break; - } - } - if (INPUT.nbands_sto == 0) - { - delete[] pchi; - } - } - - if (ik == nk - 1) - { -#ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, &stohchi.Emax, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, &stohchi.Emin, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); -#endif - stoiter.stofunc.Emax = stohchi.Emax; - stoiter.stofunc.Emin = stohchi.Emin; - GlobalV::ofs_running << "New Emax " << stohchi.Emax << " Ry; new Emin " << stohchi.Emin <<" Ry" << std::endl; - change = false; - } - } -} - -int ESolver_SDFT_PW::set_cond_nche(const double dt, - int& nbatch, - const double cond_thr, - const int& nche_min, - double try_emin, - double try_emax) -{ - int nche_guess = 1000; - ModuleBase::Chebyshev chet(nche_guess); - Stochastic_Iter& stoiter = ((hsolver::HSolverPW_SDFT*)phsol)->stoiter; - const double mu = this->pelec->eferm.ef; - stoiter.stofunc.mu = mu; - // try to find nbatch - if(nbatch == 0) - { - for (int test_nbatch = 128; test_nbatch >= 1; test_nbatch /= 2) - { - nbatch = test_nbatch; - stoiter.stofunc.t = 0.5 * dt * nbatch; - chet.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::n_sin); - double minerror = std::abs(chet.coef_complex[nche_guess - 1] / chet.coef_complex[0]); - if (minerror < cond_thr) - { - for (int i = 1; i < nche_guess; ++i) - { - double error = std::abs(chet.coef_complex[i] / chet.coef_complex[0]); - if (error < cond_thr) - { - //To make nche to around 100 - nbatch = ceil(double(test_nbatch) / i * 100.0); - std::cout << "set cond_dtbatch to " << nbatch << std::endl; - break; - } - } - break; - } - } - } - - // first try to find nche - stoiter.stofunc.t = 0.5 * dt * nbatch; - auto getnche = [&](int& nche) - { - chet.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::n_sin); - for (int i = 1; i < nche_guess; ++i) - { - double error = std::abs(chet.coef_complex[i] / chet.coef_complex[0]); - if(error < cond_thr) - { - nche = i + 1; - break; - } - } - }; - int nche_old = 0; - getnche(nche_old); - - int nche_new = 0; -loop: - // re-set Emin & Emax both in stohchi & stofunc - check_che(std::max(nche_old * 2, nche_min), try_emin, try_emax); - - // second try to find nche with new Emin & Emax - getnche(nche_new); - - if(nche_new > nche_old * 2) - { - nche_old = nche_new; - try_emin = stoiter.stohchi.Emin; - try_emax = stoiter.stohchi.Emax; - goto loop; - } - - std::cout << "set N order of Chebyshev for KG as " << nche_new << std::endl; - std::ofstream cheofs("Chebycoef"); - for (int i = 1; i < nche_guess; ++i) - { - double error = std::abs(chet.coef_complex[i] / chet.coef_complex[0]); - cheofs << std::setw(5) << i << std::setw(20) << error << std::endl; - } - cheofs.close(); - - if (nche_new >= 1000) - { - ModuleBase::WARNING_QUIT("ESolver_SDFT_PW", "N order of Chebyshev for KG will be larger than 1000!"); - } - - return nche_new; -} - -void ESolver_SDFT_PW::cal_jmatrix(const psi::Psi>& kspsi_all, - const psi::Psi>& vkspsi, - const double* en, - const double* en_all, - std::complex* leftfact, - std::complex* rightfact, - const psi::Psi>& leftchi, - psi::Psi>& rightchi, - psi::Psi>& left_hchi, - psi::Psi>& batch_vchi, - psi::Psi>& batch_vhchi, -#ifdef __MPI - psi::Psi>& chi_all, - psi::Psi>& hchi_all, - void* gatherinfo_ks, - void* gatherinfo_sto, -#endif - const int& bsize_psi, - std::vector>& j1, - std::vector>& j2, - hamilt::Velocity& velop, - const int& ik, - const std::complex& factor, - const int bandinfo[6]) -{ - ModuleBase::timer::tick(this->classname, "cal_jmatrix"); - const char transa = 'C'; - const char transb = 'N'; - const std::complex float_factor = static_cast>(factor); - const std::complex conjfactor = std::conj(float_factor); - const float mu = static_cast(this->pelec->eferm.ef); - const std::complex zero = 0.0; - const int npwx = wf.npwx; - const int npw = kv.ngk[ik]; - const int ndim = 3; - const int perbands_ks = bandinfo[0]; - const int perbands_sto = bandinfo[1]; - const int perbands = bandinfo[2]; - const int allbands_ks = bandinfo[3]; - const int allbands_sto = bandinfo[4]; - const int allbands = bandinfo[5]; - const int dim_jmatrix = perbands_ks * allbands_sto + perbands_sto * allbands; - Stochastic_Iter& stoiter = ((hsolver::HSolverPW_SDFT*)phsol)->stoiter; - - psi::Psi> right_hchi(1, perbands_sto, npwx, kv.ngk.data()); - psi::Psi> f_rightchi(1, perbands_sto, npwx, kv.ngk.data()); - psi::Psi> f_right_hchi(1, perbands_sto, npwx, kv.ngk.data()); - - stoiter.stohchi.hchi(leftchi.get_pointer(), left_hchi.get_pointer(), perbands_sto); - stoiter.stohchi.hchi(rightchi.get_pointer(), right_hchi.get_pointer(), perbands_sto); - convert_psi(rightchi, f_rightchi); - convert_psi(right_hchi, f_right_hchi); - right_hchi.resize(1, 1, 1); - - psi::Psi>* rightchi_all = &f_rightchi; - psi::Psi>* righthchi_all = &f_right_hchi; - std::complex*tmprightf_all = nullptr, *rightf_all = rightfact; -#ifdef __MPI - info_gatherv* ks_fact = static_cast(gatherinfo_ks); - info_gatherv* sto_npwx = static_cast(gatherinfo_sto); - rightchi_all = gatherchi(f_rightchi, chi_all, npwx, sto_npwx->nrecv, sto_npwx->displs, perbands_sto); - righthchi_all = gatherchi(f_right_hchi, hchi_all, npwx, sto_npwx->nrecv, sto_npwx->displs, perbands_sto); - if (GlobalV::NSTOGROUP > 1 && rightfact != nullptr) - { - tmprightf_all = new std::complex[allbands_ks]; - rightf_all = tmprightf_all; - MPI_Allgatherv(rightfact, - perbands_ks, - MPI_DOUBLE_COMPLEX, - rightf_all, - ks_fact->nrecv, - ks_fact->displs, - MPI_DOUBLE_COMPLEX, - PARAPW_WORLD); - } -#endif - - psi::Psi> f_batch_vchi(1, bsize_psi * ndim, npwx, kv.ngk.data()); - psi::Psi> f_batch_vhchi(1, bsize_psi * ndim, npwx, kv.ngk.data()); - std::vector> tmpj(ndim * allbands_sto * perbands_sto); - - // 1. (<\psi|J|\chi>)^T - // (allbands_sto, perbands_ks) - if (perbands_ks > 0) - { - for (int id = 0; id < ndim; ++id) - { - const int idnb = id * perbands_ks; - const int jbais = 0; - std::complex* j1mat = &j1[id * dim_jmatrix]; - std::complex* j2mat = &j2[id * dim_jmatrix]; - //(<\psi|v|\chi>)^T - cgemm_(&transa, - &transb, - &allbands_sto, - &perbands_ks, - &npw, - &conjfactor, - rightchi_all->get_pointer(), - &npwx, - &vkspsi(idnb, 0), - &npwx, - &zero, - j1mat, - &allbands_sto); - - //(<\psi|Hv|\chi>)^T - // for(int i = 0 ; i < perbands_ks ; ++i) - // { - // double* evmat = &j2(id, 0 + i * allbands_sto); - // double* vmat = &j1(id, 0 + i * allbands_sto); - // double ei = en[i]; - // for(int j = 0 ; j < allbands_sto ; ++j) - // { - // evmat[j] = vmat[j] * ei; - // } - // } - - //(<\psi|vH|\chi>)^T - cgemm_(&transa, - &transb, - &allbands_sto, - &perbands_ks, - &npw, - &conjfactor, - righthchi_all->get_pointer(), - &npwx, - &vkspsi(idnb, 0), - &npwx, - &zero, - j2mat, - &allbands_sto); - } - } - - int remain = perbands_sto; - int startnb = 0; - while (remain > 0) - { - int tmpnb = std::min(remain, bsize_psi); - // v|\chi> - velop.act(&leftchi, tmpnb, &leftchi(0, startnb, 0), batch_vchi.get_pointer()); - convert_psi(batch_vchi, f_batch_vchi); - // v|H\chi> - velop.act(&leftchi, tmpnb, &left_hchi(0, startnb, 0), batch_vhchi.get_pointer()); - convert_psi(batch_vhchi, f_batch_vhchi); - // 2. <\chi|J|\psi> - // (perbands_sto, allbands_ks) - if (allbands_ks > 0) - { - for (int id = 0; id < ndim; ++id) - { - const int idnb = id * tmpnb; - const int jbais = perbands_ks * allbands_sto + startnb; - std::complex* j1mat = &j1[id * dim_jmatrix + jbais]; - std::complex* j2mat = &j2[id * dim_jmatrix + jbais]; - //<\chi|v|\psi> - cgemm_(&transa, - &transb, - &tmpnb, - &allbands_ks, - &npw, - &float_factor, - &f_batch_vchi(idnb, 0), - &npwx, - kspsi_all.get_pointer(), - &npwx, - &zero, - j1mat, - &perbands_sto); - - //<\chi|vH|\psi> = \epsilon * <\chi|v|\psi> - // for(int i = 0 ; i < allbands_ks ; ++i) - // { - // double* evmat = &j2(id, jbais + i * allbands_sto); - // double* vmat = &j1(id, jbais + i * allbands_sto); - // double ei = en_all[i]; - // for(int j = 0 ; j < tmpnb ; ++j) - // { - // evmat[j] = vmat[j] * en_all[i]; - // } - // } - - //<\chi|Hv|\psi> - cgemm_(&transa, - &transb, - &tmpnb, - &allbands_ks, - &npw, - &float_factor, - &f_batch_vhchi(idnb, 0), - &npwx, - kspsi_all.get_pointer(), - &npwx, - &zero, - j2mat, - &perbands_sto); - } - } - - // 3. <\chi|J|\chi> - // (perbands_sto, allbands_sto) - for (int id = 0; id < ndim; ++id) - { - const int idnb = id * tmpnb; - const int jbais = perbands_ks * allbands_sto + perbands_sto * allbands_ks + startnb; - std::complex* j1mat = &j1[id * dim_jmatrix + jbais]; - std::complex* j2mat = &j2[id * dim_jmatrix + jbais]; - std::complex* tmpjmat = &tmpj[id * allbands_sto * perbands_sto + startnb]; - //<\chi|v|\chi> - cgemm_(&transa, - &transb, - &tmpnb, - &allbands_sto, - &npw, - &float_factor, - &f_batch_vchi(idnb, 0), - &npwx, - rightchi_all->get_pointer(), - &npwx, - &zero, - j1mat, - &perbands_sto); - - //<\chi|Hv|\chi> - cgemm_(&transa, - &transb, - &tmpnb, - &allbands_sto, - &npw, - &float_factor, - &f_batch_vhchi(idnb, 0), - &npwx, - rightchi_all->get_pointer(), - &npwx, - &zero, - j2mat, - &perbands_sto); - - //<\chi|vH|\chi> - cgemm_(&transa, - &transb, - &tmpnb, - &allbands_sto, - &npw, - &float_factor, - &f_batch_vchi(idnb, 0), - &npwx, - righthchi_all->get_pointer(), - &npwx, - &zero, - tmpjmat, - &perbands_sto); - } - - remain -= tmpnb; - startnb += tmpnb; - if (remain == 0) - break; - } - - for (int id = 0; id < ndim; ++id) - { - for (int i = 0; i < perbands_ks; ++i) - { - const float ei = static_cast(en[i]); - const int jst = i * allbands_sto; - std::complex* j2mat = j2.data() + id * dim_jmatrix + jst; - std::complex* j1mat = j1.data() + id * dim_jmatrix + jst; - if (leftfact == nullptr) - { - for (int j = 0; j < allbands_sto; ++j) - { - j2mat[j] = 0.5f * j2mat[j] + (0.5f * ei - mu) * j1mat[j]; - } - } - else - { - const std::complex jfac = static_cast>(leftfact[i]); - for (int j = 0; j < allbands_sto; ++j) - { - j2mat[j] = jfac * (0.5f * j2mat[j] + (0.5f * ei - mu) * j1mat[j]); - j1mat[j] *= jfac; - } - } - } - - for (int i = 0; i < allbands_ks; ++i) - { - const float ei = static_cast(en_all[i]); - const int jst = perbands_ks * allbands_sto + i * perbands_sto; - std::complex* j2mat = j2.data() + id * dim_jmatrix + jst; - std::complex* j1mat = j1.data() + id * dim_jmatrix + jst; - if (rightfact == nullptr) - { - for (int j = 0; j < perbands_sto; ++j) - { - j2mat[j] = 0.5f * j2mat[j] + (0.5f * ei - mu) * j1mat[j]; - } - } - else - { - const std::complex jfac = static_cast>(rightf_all[i]); - for (int j = 0; j < perbands_sto; ++j) - { - j2mat[j] = jfac * (0.5f * j2mat[j] + (0.5f * ei - mu) * j1mat[j]); - j1mat[j] *= jfac; - } - } - } - - const int jst = perbands_ks * allbands_sto + perbands_sto * allbands_ks; - const int ed = dim_jmatrix - jst; - std::complex* j2mat = j2.data() + id * dim_jmatrix + jst; - std::complex* j1mat = j1.data() + id * dim_jmatrix + jst; - std::complex* tmpjmat = tmpj.data() + id * allbands_sto * perbands_sto; - - for (int j = 0; j < ed; ++j) - { - j2mat[j] = 0.5f * (j2mat[j] + tmpjmat[j]) - mu * j1mat[j]; - } - } - -#ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, j1.data(), ndim * dim_jmatrix, MPI_COMPLEX, MPI_SUM, POOL_WORLD); - MPI_Allreduce(MPI_IN_PLACE, j2.data(), ndim * dim_jmatrix, MPI_COMPLEX, MPI_SUM, POOL_WORLD); -#endif - delete[] tmprightf_all; - - ModuleBase::timer::tick(this->classname, "cal_jmatrix"); - - return; -} - -void ESolver_SDFT_PW::sKG(const int nche_KG, - const int& smear_type, - const double fwhmin, - const double wcut, - const double dw_in, - const double dt_in, - const int nbatch, - const int npart_sto) -{ - ModuleBase::TITLE(this->classname, "sKG"); - ModuleBase::timer::tick(this->classname, "sKG"); - std::cout << "Calculating conductivity...." << std::endl; - // if (GlobalV::NSTOGROUP > 1) - // { - // ModuleBase::WARNING_QUIT("ESolver_SDFT_PW", "sKG is not supported in parallel!"); - // } - - //------------------------------------------------------------------ - // Init - //------------------------------------------------------------------ - // Parameters - int nw = ceil(wcut / dw_in); - double dw = dw_in / ModuleBase::Ry_to_eV; // converge unit in eV to Ry - double sigma = fwhmin / TWOSQRT2LN2 / ModuleBase::Ry_to_eV; - double gamma = fwhmin / 2.0 / ModuleBase::Ry_to_eV; - double dt = dt_in; // unit in a.u., 1 a.u. = 4.837771834548454e-17 s - const double expfactor = 18.42; // exp(-18.42) = 1e-8 - int nt=0; // set nt empirically - if(smear_type == 1) - { - nt = ceil(sqrt(2 * expfactor) / sigma / dt); - } - else if (smear_type == 2) - { - nt = ceil(expfactor / gamma / dt); - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW::calcondw", "smear_type should be 0 or 1"); - } - std::cout << "nw: " << nw << " ; dw: " << dw * ModuleBase::Ry_to_eV << " eV" << std::endl; - std::cout << "nt: " << nt << " ; dt: " << dt << " a.u.(ry^-1)" << std::endl; - assert(nw >= 1); - assert(nt >= 1); - const int ndim = 3; - const int nk = kv.get_nks(); - const int npwx = wf.npwx; - const double tpiba = GlobalC::ucell.tpiba; - psi::Psi>* stopsi; - if (GlobalV::NBANDS > 0) - { - stopsi = stowf.chiortho; - // clean memories //Note shchi is different from \sqrt(fH_here)|chi>, since veffs are different - stowf.shchi->resize(1, 1, 1); - stowf.chi0->resize(1, 1, 1); // clean memories - } - else - { - stopsi = stowf.chi0; - stowf.shchi->resize(1, 1, 1); // clean memories - } - const double dEcut = (wcut + fwhmin) / ModuleBase::Ry_to_eV; - - // response funtion - double* ct11 = new double[nt]; - double* ct12 = new double[nt]; - double* ct22 = new double[nt]; - ModuleBase::GlobalFunc::ZEROS(ct11, nt); - ModuleBase::GlobalFunc::ZEROS(ct12, nt); - ModuleBase::GlobalFunc::ZEROS(ct22, nt); - - // Init Chebyshev - ModuleBase::Chebyshev che(this->nche_sto); - ModuleBase::Chebyshev chet(nche_KG); - ModuleBase::Chebyshev chemt(nche_KG); - Stochastic_Iter& stoiter = ((hsolver::HSolverPW_SDFT*)phsol)->stoiter; - Stochastic_hchi& stohchi = stoiter.stohchi; - - //------------------------------------------------------------------ - // Calculate - //------------------------------------------------------------------ - - // Prepare Chebyshev coefficients for exp(-i H/\hbar t) - const double mu = this->pelec->eferm.ef; - stoiter.stofunc.mu = mu; - stoiter.stofunc.t = 0.5 * dt * nbatch; - chet.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::nsin); - chemt.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::n_sin); - std::complex*batchcoef = nullptr, *batchmcoef = nullptr; - if (nbatch > 1) - { - batchcoef = new std::complex[nche_KG * nbatch]; - std::complex* tmpcoef = batchcoef + (nbatch - 1) * nche_KG; - batchmcoef = new std::complex[nche_KG * nbatch]; - std::complex* tmpmcoef = batchmcoef + (nbatch - 1) * nche_KG; - for (int i = 0; i < nche_KG; ++i) - { - tmpcoef[i] = chet.coef_complex[i]; - tmpmcoef[i] = chemt.coef_complex[i]; - } - for (int ib = 0; ib < nbatch - 1; ++ib) - { - tmpcoef = batchcoef + ib * nche_KG; - tmpmcoef = batchmcoef + ib * nche_KG; - stoiter.stofunc.t = 0.5 * dt * (ib + 1); - chet.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::nsin); - chemt.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::n_sin); - for (int i = 0; i < nche_KG; ++i) - { - tmpcoef[i] = chet.coef_complex[i]; - tmpmcoef[i] = chemt.coef_complex[i]; - } - } - stoiter.stofunc.t = 0.5 * dt * nbatch; - } - - // ik loop - ModuleBase::timer::tick(this->classname, "kloop"); - hamilt::Velocity velop(pw_wfc, kv.isk.data(), &GlobalC::ppcell, &GlobalC::ucell, INPUT.cond_nonlocal); - for (int ik = 0; ik < nk; ++ik) - { - velop.init(ik); - stopsi->fix_k(ik); - psi->fix_k(ik); - if (nk > 1) - { - this->p_hamilt->updateHk(ik); - } - stoiter.stohchi.current_ik = ik; - const int npw = kv.ngk[ik]; - - // get allbands_ks - int cutib0 = 0; - if (GlobalV::NBANDS > 0) - { - double Emax_KS = std::max(stoiter.stofunc.Emin, this->pelec->ekb(ik, GlobalV::NBANDS - 1)); - for (cutib0 = GlobalV::NBANDS - 1; cutib0 >= 0; --cutib0) - { - if (Emax_KS - this->pelec->ekb(ik, cutib0) > dEcut) - { - break; - } - } - ++cutib0; - double Emin_KS = (cutib0 < GlobalV::NBANDS) ? this->pelec->ekb(ik, cutib0) : stoiter.stofunc.Emin; - double dE = stoiter.stofunc.Emax - Emin_KS + wcut / ModuleBase::Ry_to_eV; - std::cout << "Emin_KS(" << cutib0+1 << "): " << Emin_KS * ModuleBase::Ry_to_eV - << " eV; Emax: " << stoiter.stofunc.Emax * ModuleBase::Ry_to_eV - << " eV; Recommended max dt: " << 2 * M_PI / dE << " a.u." << std::endl; - } - else - { - double dE = stoiter.stofunc.Emax - stoiter.stofunc.Emin + wcut / ModuleBase::Ry_to_eV; - std::cout << "Emin: " << stoiter.stofunc.Emin * ModuleBase::Ry_to_eV - << " eV; Emax: " << stoiter.stofunc.Emax * ModuleBase::Ry_to_eV - << " eV; Recommended max dt: " << 2 * M_PI / dE << " a.u." << std::endl; - } - // Parallel for bands - int allbands_ks = GlobalV::NBANDS - cutib0; - parallel_distribution paraks(allbands_ks, GlobalV::NSTOGROUP, GlobalV::MY_STOGROUP); - int perbands_ks = paraks.num_per; - int ib0_ks = paraks.start; - ib0_ks += GlobalV::NBANDS - allbands_ks; - int perbands_sto = this->stowf.nchip[ik]; - int perbands = perbands_sto + perbands_ks; - int allbands_sto = perbands_sto; - int allbands = perbands; -#ifdef __MPI - MPI_Allreduce(&perbands, &allbands, 1, MPI_INT, MPI_SUM, PARAPW_WORLD); - allbands_sto = allbands - allbands_ks; - info_gatherv ks_fact(perbands_ks, GlobalV::NSTOGROUP, 1, PARAPW_WORLD); - info_gatherv sto_npwx(perbands_sto, GlobalV::NSTOGROUP, npwx, PARAPW_WORLD); -#endif - const int bandsinfo[6]{perbands_ks, perbands_sto, perbands, allbands_ks, allbands_sto, allbands}; - double *en = nullptr, *en_all = nullptr; - if (allbands_ks > 0) - { - en_all = &(this->pelec->ekb(ik, GlobalV::NBANDS - allbands_ks)); - } - if (perbands_ks > 0) - { - en = new double[perbands_ks]; - for (int ib = 0; ib < perbands_ks; ++ib) - { - en[ib] = this->pelec->ekb(ik, ib0_ks + ib); - } - } - - //----------------------------------------------------------- - // ks conductivity - //----------------------------------------------------------- - if (GlobalV::MY_STOGROUP == 0 && allbands_ks > 0) - { - jjcorr_ks(ik, nt, dt, dEcut, this->pelec->wg, velop, ct11, ct12, ct22); - } - - //----------------------------------------------------------- - // sto conductivity - //----------------------------------------------------------- - //------------------- allocate ------------------------- - size_t ks_memory_cost = perbands_ks * npwx * sizeof(std::complex); - psi::Psi> kspsi(1, perbands_ks, npwx, kv.ngk.data()); - psi::Psi> vkspsi(1, perbands_ks * ndim, npwx, kv.ngk.data()); - std::vector> expmtmf_fact(perbands_ks), expmtf_fact(perbands_ks); - psi::Psi> f_kspsi(1, perbands_ks, npwx, kv.ngk.data()); - ModuleBase::Memory::record("SDFT::kspsi", ks_memory_cost); - psi::Psi> f_vkspsi(1, perbands_ks * ndim, npwx, kv.ngk.data()); - ModuleBase::Memory::record("SDFT::vkspsi", ks_memory_cost); - psi::Psi>* kspsi_all = &f_kspsi; - - size_t sto_memory_cost = perbands_sto * npwx * sizeof(std::complex); - psi::Psi> sfchi(1, perbands_sto, npwx, kv.ngk.data()); - ModuleBase::Memory::record("SDFT::sfchi", sto_memory_cost); - psi::Psi> smfchi(1, perbands_sto, npwx, kv.ngk.data()); - ModuleBase::Memory::record("SDFT::smfchi", sto_memory_cost); -#ifdef __MPI - psi::Psi> chi_all, hchi_all, psi_all; - if (GlobalV::NSTOGROUP > 1) - { - chi_all.resize(1, allbands_sto, npwx); - hchi_all.resize(1, allbands_sto, npwx); - ModuleBase::Memory::record("SDFT::chi_all", allbands_sto * npwx * sizeof(std::complex)); - ModuleBase::Memory::record("SDFT::hchi_all", allbands_sto * npwx * sizeof(std::complex)); - psi_all.resize(1, allbands_ks, npwx); - ModuleBase::Memory::record("SDFT::kspsi_all", allbands_ks * npwx * sizeof(std::complex)); - for (int ib = 0; ib < allbands_ks; ++ib) - { - for (int ig = 0; ig < npw; ++ig) - { - psi_all(0, ib, ig) - = static_cast>(psi[0](GlobalV::NBANDS - allbands_ks + ib, ig)); - } - } - kspsi_all = &psi_all; - f_kspsi.resize(1, 1, 1); - } -#endif - - const int nbatch_psi = npart_sto; - const int bsize_psi = ceil(double(perbands_sto) / nbatch_psi); - psi::Psi> batch_vchi(1, bsize_psi * ndim, npwx, kv.ngk.data()); - psi::Psi> batch_vhchi(1, bsize_psi * ndim, npwx, kv.ngk.data()); - ModuleBase::Memory::record("SDFT::batchjpsi", 3 * bsize_psi * ndim * npwx * sizeof(std::complex)); - - //------------------- sqrt(f)|psi> sqrt(1-f)|psi> --------------- - if(perbands_ks > 0) - { - for (int ib = 0; ib < perbands_ks; ++ib) - { - for (int ig = 0; ig < npw; ++ig) - { - kspsi(0, ib, ig) = psi[0](ib0_ks + ib, ig); - } - double fi = stoiter.stofunc.fd(en[ib]); - expmtmf_fact[ib] = 1 - fi; - expmtf_fact[ib] = fi; - } - // v|\psi> - velop.act(&kspsi, perbands_ks, kspsi.get_pointer(), vkspsi.get_pointer()); - // convert to complex - if (GlobalV::NSTOGROUP == 1) - { - convert_psi(kspsi, f_kspsi); - } - convert_psi(vkspsi, f_vkspsi); - kspsi.resize(1, 1, 1); - vkspsi.resize(1, 1, 1); - } - - che.calcoef_real(&stoiter.stofunc, &Sto_Func::nroot_fd); - che.calfinalvec_real(&stohchi, - &Stochastic_hchi::hchi_norm, - stopsi->get_pointer(), - sfchi.get_pointer(), - npw, - npwx, - perbands_sto); - - che.calcoef_real(&stoiter.stofunc, &Sto_Func::nroot_mfd); - - che.calfinalvec_real(&stohchi, - &Stochastic_hchi::hchi_norm, - stopsi->get_pointer(), - smfchi.get_pointer(), - npw, - npwx, - perbands_sto); - - //------------------------ allocate ------------------------ - psi::Psi>& expmtsfchi = sfchi; - psi::Psi>& expmtsmfchi = smfchi; - psi::Psi> exptsfchi = expmtsfchi; - ModuleBase::Memory::record("SDFT::exptsfchi", sto_memory_cost); - psi::Psi> exptsmfchi = expmtsmfchi; - ModuleBase::Memory::record("SDFT::exptsmfchi", sto_memory_cost); - psi::Psi> poly_expmtsfchi, poly_expmtsmfchi; - psi::Psi> poly_exptsfchi, poly_exptsmfchi; - if (nbatch > 1) - { - poly_exptsfchi.resize(nche_KG, perbands_sto, npwx); - ModuleBase::Memory::record("SDFT::poly_exptsfchi", - sizeof(std::complex) * nche_KG * perbands_sto * npwx); - - poly_exptsmfchi.resize(nche_KG, perbands_sto, npwx); - ModuleBase::Memory::record("SDFT::poly_exptsmfchi", - sizeof(std::complex) * nche_KG * perbands_sto * npwx); - - poly_expmtsfchi.resize(nche_KG, perbands_sto, npwx); - ModuleBase::Memory::record("SDFT::poly_expmtsfchi", - sizeof(std::complex) * nche_KG * perbands_sto * npwx); - - poly_expmtsmfchi.resize(nche_KG, perbands_sto, npwx); - ModuleBase::Memory::record("SDFT::poly_expmtsmfchi", - sizeof(std::complex) * nche_KG * perbands_sto * npwx); - } - - const int dim_jmatrix = perbands_ks * allbands_sto + perbands_sto * allbands; - parallel_distribution parajmat(ndim * dim_jmatrix, GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL); - std::vector> j1l(ndim * dim_jmatrix), j2l(ndim * dim_jmatrix); - ModuleBase::Memory::record("SDFT::j1l", sizeof(std::complex) * ndim * dim_jmatrix); - ModuleBase::Memory::record("SDFT::j2l", sizeof(std::complex) * ndim * dim_jmatrix); - std::vector> j1r(ndim * dim_jmatrix), j2r(ndim * dim_jmatrix); - ModuleBase::Memory::record("SDFT::j1r", sizeof(std::complex) * ndim * dim_jmatrix); - ModuleBase::Memory::record("SDFT::j2r", sizeof(std::complex) * ndim * dim_jmatrix); - psi::Psi> tmphchil(1, perbands_sto, npwx, kv.ngk.data()); - ModuleBase::Memory::record("SDFT::tmphchil/r", sto_memory_cost * 2); - - //------------------------ t loop -------------------------- - std::cout << "ik=" << ik << ": "; - auto start = std::chrono::high_resolution_clock::now(); - const int print_step = ceil(20.0 / nbatch) * nbatch; - for (int it = 1; it < nt; ++it) - { - // evaluate time cost - if (it - 1 == print_step) - { - auto end = std::chrono::high_resolution_clock::now(); - std::chrono::duration duration = end - start; - double timeTaken = duration.count(); - std::cout << "(Time left " << timeTaken * (double(nt - 1) / print_step * (nk - ik) - 1) << " s) " - << std::endl; - std::cout << "nt: "< 1) - { - std::cout < - // KS - ModuleBase::timer::tick(this->classname, "evolution"); - for (int ib = 0; ib < perbands_ks; ++ib) - { - double eigen = en[ib]; - const std::complex expmfactor = exp(ModuleBase::NEG_IMAG_UNIT * eigen * dt); - expmtf_fact[ib] *= expmfactor; - expmtmf_fact[ib] *= expmfactor; - } - // Sto - if (nbatch == 1) - { - chemt.calfinalvec_complex(&stohchi, - &Stochastic_hchi::hchi_norm, - expmtsfchi.get_pointer(), - expmtsfchi.get_pointer(), - npw, - npwx, - perbands_sto); - chemt.calfinalvec_complex(&stohchi, - &Stochastic_hchi::hchi_norm, - expmtsmfchi.get_pointer(), - expmtsmfchi.get_pointer(), - npw, - npwx, - perbands_sto); - chet.calfinalvec_complex(&stohchi, - &Stochastic_hchi::hchi_norm, - exptsfchi.get_pointer(), - exptsfchi.get_pointer(), - npw, - npwx, - perbands_sto); - chet.calfinalvec_complex(&stohchi, - &Stochastic_hchi::hchi_norm, - exptsmfchi.get_pointer(), - exptsmfchi.get_pointer(), - npw, - npwx, - perbands_sto); - } - else - { - std::complex* tmppolyexpmtsfchi = poly_expmtsfchi.get_pointer(); - std::complex* tmppolyexpmtsmfchi = poly_expmtsmfchi.get_pointer(); - std::complex* tmppolyexptsfchi = poly_exptsfchi.get_pointer(); - std::complex* tmppolyexptsmfchi = poly_exptsmfchi.get_pointer(); - std::complex* stoexpmtsfchi = expmtsfchi.get_pointer(); - std::complex* stoexpmtsmfchi = expmtsmfchi.get_pointer(); - std::complex* stoexptsfchi = exptsfchi.get_pointer(); - std::complex* stoexptsmfchi = exptsmfchi.get_pointer(); - if ((it - 1) % nbatch == 0) - { - chet.calpolyvec_complex(&stohchi, - &Stochastic_hchi::hchi_norm, - stoexptsfchi, - tmppolyexptsfchi, - npw, - npwx, - perbands_sto); - chet.calpolyvec_complex(&stohchi, - &Stochastic_hchi::hchi_norm, - stoexptsmfchi, - tmppolyexptsmfchi, - npw, - npwx, - perbands_sto); - chemt.calpolyvec_complex(&stohchi, - &Stochastic_hchi::hchi_norm, - stoexpmtsfchi, - tmppolyexpmtsfchi, - npw, - npwx, - perbands_sto); - chemt.calpolyvec_complex(&stohchi, - &Stochastic_hchi::hchi_norm, - stoexpmtsmfchi, - tmppolyexpmtsmfchi, - npw, - npwx, - perbands_sto); - } - - std::complex* tmpcoef = batchcoef + (it - 1) % nbatch * nche_KG; - std::complex* tmpmcoef = batchmcoef + (it - 1) % nbatch * nche_KG; - const char transa = 'N'; - const int LDA = perbands_sto * npwx; - const int M = perbands_sto * npwx; - const int N = nche_KG; - const int inc = 1; - zgemv_(&transa, - &M, - &N, - &ModuleBase::ONE, - tmppolyexptsfchi, - &LDA, - tmpcoef, - &inc, - &ModuleBase::ZERO, - stoexptsfchi, - &inc); - zgemv_(&transa, - &M, - &N, - &ModuleBase::ONE, - tmppolyexptsmfchi, - &LDA, - tmpcoef, - &inc, - &ModuleBase::ZERO, - stoexptsmfchi, - &inc); - zgemv_(&transa, - &M, - &N, - &ModuleBase::ONE, - tmppolyexpmtsfchi, - &LDA, - tmpmcoef, - &inc, - &ModuleBase::ZERO, - stoexpmtsfchi, - &inc); - zgemv_(&transa, - &M, - &N, - &ModuleBase::ONE, - tmppolyexpmtsmfchi, - &LDA, - tmpmcoef, - &inc, - &ModuleBase::ZERO, - stoexpmtsmfchi, - &inc); - } - ModuleBase::timer::tick(this->classname, "evolution"); - - // calculate i<\psi|sqrt(f) exp(-iHt/2)*J*exp(iHt/2) sqrt(1-f)|\psi>^+ - // = i<\psi|sqrt(1-f) exp(-iHt/2)*J*exp(iHt/2) sqrt(f)|\psi> - cal_jmatrix(*kspsi_all, - f_vkspsi, - en, - en_all, - nullptr, - nullptr, - exptsmfchi, - exptsfchi, - tmphchil, - batch_vchi, - batch_vhchi, -#ifdef __MPI - chi_all, - hchi_all, - (void*)&ks_fact, - (void*)&sto_npwx, -#endif - bsize_psi, - j1l, - j2l, - velop, - ik, - ModuleBase::IMAG_UNIT, - bandsinfo); - - // calculate <\psi|sqrt(1-f) exp(iHt/2)*J*exp(-iHt/2) sqrt(f)|\psi> - cal_jmatrix(*kspsi_all, - f_vkspsi, - en, - en_all, - expmtmf_fact.data(), - expmtf_fact.data(), - expmtsmfchi, - expmtsfchi, - tmphchil, - batch_vchi, - batch_vhchi, -#ifdef __MPI - chi_all, - hchi_all, - (void*)&ks_fact, - (void*)&sto_npwx, -#endif - bsize_psi, - j1r, - j2r, - velop, - ik, - ModuleBase::ONE, - bandsinfo); - - // prepare for parallel - int num_per = parajmat.num_per; - int st_per = parajmat.start; - // Re(i) - // Im(l_ij*r_ji) = Re(-il_ij * r_ji) = Re( ((il)^+_ji)^* * r_ji)=Re(((il)^+_i)^* * r^+_i) - // ddot_real = real(A_i^* * B_i) - ModuleBase::timer::tick(this->classname, "ddot_real"); - ct11[it] += static_cast( - ModuleBase::GlobalFunc::ddot_real(num_per, j1l.data() + st_per, j1r.data() + st_per, false) * kv.wk[ik] - / 2.0); - double tmp12 = static_cast( - ModuleBase::GlobalFunc::ddot_real(num_per, j1l.data() + st_per, j2r.data() + st_per, false)); - - double tmp21 = static_cast( - ModuleBase::GlobalFunc::ddot_real(num_per, j2l.data() + st_per, j1r.data() + st_per, false)); - - ct12[it] -= 0.5 * (tmp12 + tmp21) * kv.wk[ik] / 2.0; - - ct22[it] += static_cast( - ModuleBase::GlobalFunc::ddot_real(num_per, j2l.data() + st_per, j2r.data() + st_per, false) * kv.wk[ik] - / 2.0); - - ModuleBase::timer::tick(this->classname, "ddot_real"); - } - std::cout << std::endl; - delete[] en; - } // ik loop - ModuleBase::timer::tick(this->classname, "kloop"); - delete[] batchcoef; - delete[] batchmcoef; - -#ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, ct11, nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, ct12, nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, ct22, nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); -#endif - - //------------------------------------------------------------------ - // Output - //------------------------------------------------------------------ - if (GlobalV::MY_RANK == 0) - { - calcondw(nt, dt, smear_type, fwhmin, wcut, dw_in, ct11, ct12, ct22); - } - delete[] ct11; - delete[] ct12; - delete[] ct22; - ModuleBase::timer::tick(this->classname, "sKG"); -} - -void ESolver_SDFT_PW::caldos(const int nche_dos, - const double sigmain, - const double emin, - const double emax, - const double de, - const int npart) -{ - ModuleBase::TITLE(this->classname, "caldos"); - ModuleBase::timer::tick(this->classname, "caldos"); - std::cout << "=========================" << std::endl; - std::cout << "###Calculating Dos....###" << std::endl; - std::cout << "=========================" << std::endl; - ModuleBase::Chebyshev che(nche_dos); - const int nk = kv.get_nks(); - Stochastic_Iter& stoiter = ((hsolver::HSolverPW_SDFT*)phsol)->stoiter; - Stochastic_hchi& stohchi = stoiter.stohchi; - const int npwx = wf.npwx; - - double* spolyv = nullptr; - std::complex* allorderchi = nullptr; - if (stoiter.method == 1) - { - spolyv = new double[nche_dos]; - ModuleBase::GlobalFunc::ZEROS(spolyv, nche_dos); - } - else - { - spolyv = new double[nche_dos * nche_dos]; - ModuleBase::GlobalFunc::ZEROS(spolyv, nche_dos * nche_dos); - int nchip_new = ceil((double)this->stowf.nchip_max / npart); - allorderchi = new std::complex[nchip_new * npwx * nche_dos]; - } - ModuleBase::timer::tick(this->classname, "Tracepoly"); - std::cout << "1. TracepolyA:" << std::endl; - for (int ik = 0; ik < nk; ik++) - { - std::cout << "ik: " << ik + 1 << std::endl; - if (nk > 1) - { - this->p_hamilt->updateHk(ik); - } - stohchi.current_ik = ik; - const int npw = kv.ngk[ik]; - const int nchipk = this->stowf.nchip[ik]; - - std::complex* pchi; - if (GlobalV::NBANDS > 0) - { - stowf.chiortho->fix_k(ik); - pchi = stowf.chiortho->get_pointer(); - } - else - { - stowf.chi0->fix_k(ik); - pchi = stowf.chi0->get_pointer(); - } - if (stoiter.method == 1) - { - che.tracepolyA(&stohchi, &Stochastic_hchi::hchi_norm, pchi, npw, npwx, nchipk); - for (int i = 0; i < nche_dos; ++i) - { - spolyv[i] += che.polytrace[i] * kv.wk[ik] / 2; - } - } - else - { - int N = nche_dos; - double kweight = kv.wk[ik] / 2; - char trans = 'T'; - char normal = 'N'; - double one = 1; - for (int ipart = 0; ipart < npart; ++ipart) - { - int nchipk_new = nchipk / npart; - int start_nchipk = ipart * nchipk_new + nchipk % npart; - if (ipart < nchipk % npart) - { - nchipk_new++; - start_nchipk = ipart * nchipk_new; - } - ModuleBase::GlobalFunc::ZEROS(allorderchi, nchipk_new * npwx * nche_dos); - std::complex* tmpchi = pchi + start_nchipk * npwx; - che.calpolyvec_complex(&stohchi, - &Stochastic_hchi::hchi_norm, - tmpchi, - allorderchi, - npw, - npwx, - nchipk_new); - double* vec_all = (double*)allorderchi; - int LDA = npwx * nchipk_new * 2; - int M = npwx * nchipk_new * 2; - dgemm_(&trans, &normal, &N, &N, &M, &kweight, vec_all, &LDA, vec_all, &LDA, &one, spolyv, &N); - } - } - } - - if (stoiter.method == 2) - { - delete[] allorderchi; - } - - std::ofstream ofsdos; - int ndos = int((emax - emin) / de) + 1; - stoiter.stofunc.sigma = sigmain / ModuleBase::Ry_to_eV; - ModuleBase::timer::tick(this->classname, "Tracepoly"); - - std::cout << "2. Dos:" << std::endl; - ModuleBase::timer::tick(this->classname, "DOS Loop"); - int n10 = ndos / 10; - int percent = 10; - double* sto_dos = new double[ndos]; - double* ks_dos = new double[ndos]; - double* error = new double[ndos]; - for (int ie = 0; ie < ndos; ++ie) - { - double tmpks = 0; - double tmpsto = 0; - stoiter.stofunc.targ_e = (emin + ie * de) / ModuleBase::Ry_to_eV; - if (stoiter.method == 1) - { - che.calcoef_real(&stoiter.stofunc, &Sto_Func::ngauss); - tmpsto = BlasConnector::dot(nche_dos, che.coef_real, 1, spolyv, 1); - } - else - { - che.calcoef_real(&stoiter.stofunc, &Sto_Func::nroot_gauss); - tmpsto = stoiter.vTMv(che.coef_real, spolyv, nche_dos); - } - if (GlobalV::NBANDS > 0) - { - for (int ik = 0; ik < nk; ++ik) - { - double* en = &(this->pelec->ekb(ik, 0)); - for (int ib = 0; ib < GlobalV::NBANDS; ++ib) - { - tmpks += stoiter.stofunc.gauss(en[ib]) * kv.wk[ik] / 2; - } - } - } - tmpks /= GlobalV::NPROC_IN_POOL; - - double tmperror = 0; - if (stoiter.method == 1) - { - tmperror = che.coef_real[nche_dos - 1] * spolyv[nche_dos - 1]; - } - else - { - const int norder = nche_dos; - double last_coef = che.coef_real[norder - 1]; - double last_spolyv = spolyv[norder * norder - 1]; - tmperror = last_coef - * (BlasConnector::dot(norder, che.coef_real, 1, spolyv + norder * (norder - 1), 1) - + BlasConnector::dot(norder, che.coef_real, 1, spolyv + norder - 1, norder) - - last_coef * last_spolyv); - } - - if (ie % n10 == n10 - 1) - { - std::cout << percent << "%" - << " "; - percent += 10; - } - sto_dos[ie] = tmpsto; - ks_dos[ie] = tmpks; - error[ie] = tmperror; - } -#ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, ks_dos, ndos, MPI_DOUBLE, MPI_SUM, STO_WORLD); - MPI_Allreduce(MPI_IN_PLACE, sto_dos, ndos, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, error, ndos, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); -#endif - if (GlobalV::MY_RANK == 0) - { - std::string dosfile = GlobalV::global_out_dir + "DOS1_smearing.dat"; - ofsdos.open(dosfile.c_str()); - double maxerror = 0; - double sum = 0; - ofsdos << std::setw(8) << "## E(eV) " << std::setw(20) << "dos(eV^-1)" << std::setw(20) << "sum" - << std::setw(20) << "Error(eV^-1)" << std::endl; - for (int ie = 0; ie < ndos; ++ie) - { - double tmperror = 2.0 * std::abs(error[ie]); - if (maxerror < tmperror) - maxerror = tmperror; - double dos = 2.0 * (ks_dos[ie] + sto_dos[ie]) / ModuleBase::Ry_to_eV; - sum += dos; - ofsdos << std::setw(8) << emin + ie * de << std::setw(20) << dos << std::setw(20) << sum * de - << std::setw(20) << tmperror << std::endl; - } - std::cout << std::endl; - std::cout << "Finish DOS" << std::endl; - std::cout << std::scientific << "DOS max absolute Chebyshev Error: " << maxerror << std::endl; - ofsdos.close(); - } - delete[] sto_dos; - delete[] ks_dos; - delete[] error; - delete[] spolyv; - ModuleBase::timer::tick(this->classname, "DOS Loop"); - ModuleBase::timer::tick(this->classname, "caldos"); - return; -} - -} // namespace ModuleESolver - -namespace GlobalTemp -{ - const ModuleBase::matrix* veff; -} diff --git a/source/module_hamilt_pw/hamilt_pwdft/CMakeLists.txt b/source/module_hamilt_pw/hamilt_pwdft/CMakeLists.txt index 977105660e..3ccd067684 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/CMakeLists.txt +++ b/source/module_hamilt_pw/hamilt_pwdft/CMakeLists.txt @@ -30,6 +30,7 @@ list(APPEND objects soc.cpp global.cpp parallel_grid.cpp + elecond.cpp ) add_library( diff --git a/source/module_esolver/esolver_ks_pw_tool.cpp b/source/module_hamilt_pw/hamilt_pwdft/elecond.cpp similarity index 53% rename from source/module_esolver/esolver_ks_pw_tool.cpp rename to source/module_hamilt_pw/hamilt_pwdft/elecond.cpp index 4695311e16..ebb8769e66 100644 --- a/source/module_esolver/esolver_ks_pw_tool.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/elecond.cpp @@ -1,13 +1,11 @@ -#include "esolver_ks_pw.h" +#include "elecond.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_elecstate/occupy.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/binstream.h" -namespace ModuleESolver -{ +#include //------------------------------------------------------------------ // hbar = 6.62607015e-34/2pi @@ -24,13 +22,21 @@ namespace ModuleESolver //------------------------------------------------------------------ #define TWOSQRT2LN2 2.354820045030949 // FWHM = 2sqrt(2ln2) * \sigma #define FACTOR 1.839939223835727e7 -template -void ESolver_KS_PW::KG(const int& smear_type, - const double fwhmin, - const double wcut, - const double dw_in, - const double dt_in, - ModuleBase::matrix& wg) + +EleCond::EleCond(UnitCell* p_ucell_in, K_Vectors* p_kv_in, elecstate::ElecState* p_elec_in, + ModulePW::PW_Basis_K* p_wfcpw_in, psi::Psi>* p_psi_in, + pseudopot_cell_vnl* p_ppcell_in) +{ + this->p_ppcell = p_ppcell_in; + this->p_ucell = p_ucell_in; + this->p_wfcpw = p_wfcpw_in; + this->p_kv = p_kv_in; + this->p_elec = p_elec_in; + this->p_psi = p_psi_in; +} + +void EleCond::KG(const int& smear_type, const double& fwhmin, const double& wcut, const double& dw_in, + const double& dt_in, const bool& nonlocal, ModuleBase::matrix& wg) { //----------------------------------------------------------- // KS conductivity @@ -40,16 +46,16 @@ void ESolver_KS_PW::KG(const int& smear_type, double dw = dw_in / ModuleBase::Ry_to_eV; // converge unit in eV to Ry double sigma = fwhmin / TWOSQRT2LN2 / ModuleBase::Ry_to_eV; const double gamma = fwhmin / 2.0 / ModuleBase::Ry_to_eV; - double dt = dt_in; // unit in a.u., 1 a.u. = 4.837771834548454e-17 s - const double expfactor = 23; //exp(-23) = 1e-10 - int nt; //set nt empirically - if(smear_type == 1) + double dt = dt_in; // unit in a.u., 1 a.u. = 4.837771834548454e-17 s + const double expfactor = 23; // exp(-23) = 1e-10 + int nt; // set nt empirically + if (smear_type == 1) { - nt = ceil(sqrt(2*expfactor)/sigma/dt); + nt = ceil(sqrt(2 * expfactor) / sigma / dt); } - else if(smear_type == 2) + else if (smear_type == 2) { - nt = ceil(expfactor/gamma/dt); + nt = ceil(expfactor / gamma / dt); } else { @@ -59,87 +65,61 @@ void ESolver_KS_PW::KG(const int& smear_type, std::cout << "nt: " << nt << " ; dt: " << dt << " a.u.(ry^-1)" << std::endl; assert(nw >= 1); assert(nt >= 1); - const int nk = this->kv.get_nks(); + const int nk = this->p_kv->get_nks(); - double* ct11 = new double[nt]; - double* ct12 = new double[nt]; - double* ct22 = new double[nt]; - ModuleBase::GlobalFunc::ZEROS(ct11, nt); - ModuleBase::GlobalFunc::ZEROS(ct12, nt); - ModuleBase::GlobalFunc::ZEROS(ct22, nt); + std::vector ct11(nt, 0); + std::vector ct12(nt, 0); + std::vector ct22(nt, 0); - hamilt::Velocity velop(this->pw_wfc, this->kv.isk.data(), &GlobalC::ppcell, &GlobalC::ucell, INPUT.cond_nonlocal); - double decut = (wcut + fwhmin) / ModuleBase::Ry_to_eV; - std::cout<<"Recommended dt: "<<0.25*M_PI/decut<<" a.u."<p_wfcpw, this->p_kv->isk.data(), this->p_ppcell, this->p_ucell, nonlocal); + double decut = (wcut + fwhmin) / ModuleBase::Ry_to_eV; + std::cout << "Recommended dt: " << 0.25 * M_PI / decut << " a.u." << std::endl; for (int ik = 0; ik < nk; ++ik) { velop.init(ik); - jjcorr_ks(ik, nt, dt, decut, wg, velop, ct11, ct12, ct22); + jjresponse_ks(ik, nt, dt, decut, wg, velop, ct11.data(), ct12.data(), ct22.data()); } #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, ct11, nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, ct12, nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, ct22, nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, ct11.data(), nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, ct12.data(), nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, ct22.data(), nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #endif //------------------------------------------------------------------ // Output //------------------------------------------------------------------ if (GlobalV::MY_RANK == 0) { - calcondw(nt, dt, smear_type, fwhmin, wcut, dw_in, ct11, ct12, ct22); + calcondw(nt, dt, smear_type, fwhmin, wcut, dw_in, ct11.data(), ct12.data(), ct22.data()); } - delete[] ct11; - delete[] ct12; - delete[] ct22; } -template -void ESolver_KS_PW::jjcorr_ks(const int ik, - const int nt, - const double dt, - const double decut, - ModuleBase::matrix& wg, - hamilt::Velocity& velop, - double* ct11, - double* ct12, - double* ct22) +void EleCond::jjresponse_ks(const int ik, const int nt, const double dt, const double decut, ModuleBase::matrix& wg, + hamilt::Velocity& velop, double* ct11, double* ct12, double* ct22) { const int nbands = GlobalV::NBANDS; - if(wg(ik, 0) - wg(ik, nbands - 1) < 1e-8 || nbands == 0) + if (wg(ik, 0) - wg(ik, nbands - 1) < 1e-8 || nbands == 0) return; - char transn = 'N'; - char transc = 'C'; + const char transn = 'N'; + const char transc = 'C'; const int ndim = 3; - const int npwx = this->wf.npwx; - const double ef = this->pelec->eferm.ef; - const int npw = this->kv.ngk[ik]; + const int npwx = this->p_wfcpw->npwk_max; + const double ef = this->p_elec->eferm.ef; + const int npw = this->p_kv->ngk[ik]; const int reducenb2 = (nbands - 1) * nbands / 2; - bool gamma_only = false; // ABACUS do not support gamma_only yet. - std::complex* levc = &(this->psi[0](ik, 0, 0)); - std::complex* prevc = new std::complex[3 * npwx * nbands]; - std::complex* pij = new std::complex[nbands * nbands]; - double* pij2 = new double[reducenb2]; - ModuleBase::GlobalFunc::ZEROS(pij2, reducenb2); + const bool gamma_only = false; // ABACUS do not support gamma_only yet. + std::complex* levc = &(this->p_psi[0](ik, 0, 0)); + std::vector> prevc(ndim * npwx * nbands); + std::vector> pij(nbands * nbands); + std::vector pij2(reducenb2, 0); // px|right> - velop.act(this->psi, nbands * GlobalV::NPOL, levc, prevc); + velop.act(this->p_psi, nbands * GlobalV::NPOL, levc, prevc.data()); for (int id = 0; id < ndim; ++id) { - zgemm_(&transc, - &transn, - &nbands, - &nbands, - &npw, - &ModuleBase::ONE, - levc, - &npwx, - prevc + id * npwx * nbands, - &npwx, - &ModuleBase::ZERO, - pij, - &nbands); + zgemm_(&transc, &transn, &nbands, &nbands, &npw, &ModuleBase::ONE, levc, &npwx, + prevc.data() + id * npwx * nbands, &npwx, &ModuleBase::ZERO, pij.data(), &nbands); #ifdef __MPI - MPI_Allreduce(MPI_IN_PLACE, pij, nbands * nbands, MPI_DOUBLE_COMPLEX, MPI_SUM, POOL_WORLD); + MPI_Allreduce(MPI_IN_PLACE, pij.data(), nbands * nbands, MPI_DOUBLE_COMPLEX, MPI_SUM, POOL_WORLD); #endif if (!gamma_only) for (int ib = 0, ijb = 0; ib < nbands; ++ib) @@ -153,13 +133,13 @@ void ESolver_KS_PW::jjcorr_ks(const int ik, if (GlobalV::RANK_IN_POOL == 0) { - int nkstot = this->kv.get_nkstot(); - int ikglobal = this->kv.getik_global(ik); + int nkstot = this->p_kv->get_nkstot(); + int ikglobal = this->p_kv->getik_global(ik); std::stringstream ss; ss << GlobalV::global_out_dir << "vmatrix" << ikglobal + 1 << ".dat"; Binstream binpij(ss.str(), "w"); binpij << 8 * reducenb2; - binpij.write(pij2, reducenb2); + binpij.write(pij2.data(), reducenb2); binpij << 8 * reducenb2; } @@ -180,7 +160,7 @@ void ESolver_KS_PW::jjcorr_ks(const int ik, double tmct11 = 0; double tmct12 = 0; double tmct22 = 0; - double* enb = &(this->pelec->ekb(ik, 0)); + double* enb = &(this->p_elec->ekb(ik, 0)); for (int ib = 0, ijb = 0; ib < nbands; ++ib) { double ei = enb[ib]; @@ -188,7 +168,7 @@ void ESolver_KS_PW::jjcorr_ks(const int ik, for (int jb = ib + 1; jb < nbands; ++jb, ++ijb) { double ej = enb[jb]; - if (ej - ei > decut ) + if (ej - ei > decut) continue; double fj = wg(ik, jb); double tmct = sin((ej - ei) * (it)*dt) * (fi - fj) * pij2[ijb]; @@ -201,22 +181,11 @@ void ESolver_KS_PW::jjcorr_ks(const int ik, ct12[it] += tmct12 / 2.0; ct22[it] += tmct22 / 2.0; } - delete[] pij; - delete[] prevc; - delete[] pij2; return; } -template -void ESolver_KS_PW::calcondw(const int nt, - const double dt, - const int& smear_type, - const double fwhmin, - const double wcut, - const double dw_in, - double* ct11, - double* ct12, - double* ct22) +void EleCond::calcondw(const int nt, const double dt, const int& smear_type, const double fwhmin, const double wcut, + const double dw_in, double* ct11, double* ct12, double* ct22) { double factor = FACTOR; const int ndim = 3; @@ -224,40 +193,36 @@ void ESolver_KS_PW::calcondw(const int nt, double dw = dw_in / ModuleBase::Ry_to_eV; // converge unit in eV to Ry const double sigma = fwhmin / TWOSQRT2LN2 / ModuleBase::Ry_to_eV; const double gamma = fwhmin / 2.0 / ModuleBase::Ry_to_eV; - double* winfunc = new double[nt]; - //1: Gaussian, 2: Lorentzian - if(smear_type == 1) + std::vector winfunc(nt); + // 1: Gaussian, 2: Lorentzian + if (smear_type == 1) { for (int it = 0; it < nt; ++it) { winfunc[it] = exp(-double(1) / 2 * sigma * sigma * pow((it)*dt, 2)); } } - else if(smear_type == 2) + else if (smear_type == 2) { for (int it = 0; it < nt; ++it) { winfunc[it] = exp(-gamma * (it)*dt); } } - + std::ofstream ofscond("je-je.txt"); ofscond << std::setw(8) << "#t(a.u.)" << std::setw(15) << "c11(t)" << std::setw(15) << "c12(t)" << std::setw(15) << "c22(t)" << std::setw(15) << "decay" << std::endl; for (int it = 0; it < nt; ++it) { ofscond << std::setw(8) << (it)*dt << std::setw(15) << -2 * ct11[it] << std::setw(15) << -2 * ct12[it] - << std::setw(15) << -2 * ct22[it] << std::setw(15) - << winfunc[it] << std::endl; + << std::setw(15) << -2 * ct22[it] << std::setw(15) << winfunc[it] << std::endl; } ofscond.close(); - double* cw11 = new double[nw]; - double* cw12 = new double[nw]; - double* cw22 = new double[nw]; - double* kappa = new double[nw]; - ModuleBase::GlobalFunc::ZEROS(cw11, nw); - ModuleBase::GlobalFunc::ZEROS(cw12, nw); - ModuleBase::GlobalFunc::ZEROS(cw22, nw); + std::vector cw11(nw, 0); + std::vector cw12(nw, 0); + std::vector cw22(nw, 0); + std::vector kappa(nw); for (int iw = 0; iw < nw; ++iw) { for (int it = 0; it < nt; ++it) @@ -272,10 +237,10 @@ void ESolver_KS_PW::calcondw(const int nt, << std::setw(20) << "L12/e(Am^-1)" << std::setw(20) << "L22/e^2(Wm^-1)" << std::endl; for (int iw = 0; iw < nw; ++iw) { - cw11[iw] *= double(2) / ndim / GlobalC::ucell.omega * factor; // unit in Sm^-1 + cw11[iw] *= double(2) / ndim / this->p_ucell->omega * factor; // unit in Sm^-1 cw12[iw] - *= double(2) / ndim / GlobalC::ucell.omega * factor * 2.17987092759e-18 / 1.6021766208e-19; // unit in Am^-1 - cw22[iw] *= double(2) / ndim / GlobalC::ucell.omega * factor + *= double(2) / ndim / this->p_ucell->omega * factor * 2.17987092759e-18 / 1.6021766208e-19; // unit in Am^-1 + cw22[iw] *= double(2) / ndim / this->p_ucell->omega * factor * pow(2.17987092759e-18 / 1.6021766208e-19, 2); // unit in Wm^-1 kappa[iw] = (cw22[iw] - pow(cw12[iw], 2) / cw11[iw]) / Occupy::gaussian_parameter / ModuleBase::Ry_to_eV / 11604.518026; @@ -283,24 +248,11 @@ void ESolver_KS_PW::calcondw(const int nt, << kappa[iw] << std::setw(20) << cw12[iw] << std::setw(20) << cw22[iw] << std::endl; } double sigma0 = cw11[0] - (cw11[1] - cw11[0]) * 0.5; - double kappa0 = kappa[0] - (kappa[1] - kappa[0]) * 0.5; - double Lorent0 = kappa0 / sigma0 / Occupy::gaussian_parameter / ModuleBase::Ry_to_eV / 11604.518026 * pow(1.6021766208e-19/1.3806505e-23, 2); + double kappa0 = kappa[0] - (kappa[1] - kappa[0]) * 0.5; + double Lorent0 = kappa0 / sigma0 / Occupy::gaussian_parameter / ModuleBase::Ry_to_eV / 11604.518026 + * pow(1.6021766208e-19 / 1.3806505e-23, 2); std::cout << std::setprecision(6) << "DC electrical conductivity: " << sigma0 << " Sm^-1" << std::endl; std::cout << std::setprecision(6) << "Thermal conductivity: " << kappa0 << " W(mK)^-1" << std::endl; - std::cout << std::setprecision(6) << "Lorenz number: "<, base_device::DEVICE_CPU>; -template class ESolver_KS_PW, base_device::DEVICE_CPU>; -#if ((defined __CUDA) || (defined __ROCM)) -template class ESolver_KS_PW, base_device::DEVICE_GPU>; -template class ESolver_KS_PW, base_device::DEVICE_GPU>; -#endif -} // namespace ModuleESolver \ No newline at end of file +} \ No newline at end of file diff --git a/source/module_hamilt_pw/hamilt_pwdft/elecond.h b/source/module_hamilt_pw/hamilt_pwdft/elecond.h new file mode 100644 index 0000000000..a8c99b22dc --- /dev/null +++ b/source/module_hamilt_pw/hamilt_pwdft/elecond.h @@ -0,0 +1,74 @@ +#ifndef ELECOND_H +#define ELECOND_H + +#include "module_base/matrix.h" +#include "module_basis/module_pw/pw_basis_k.h" +#include "module_cell/klist.h" +#include "module_cell/unitcell.h" +#include "module_elecstate/elecstate.h" +#include "module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h" +#include "module_hamilt_pw/hamilt_pwdft/operator_pw/velocity_pw.h" + +class EleCond +{ + public: + EleCond(UnitCell* p_ucell_in, K_Vectors* p_kv_in, elecstate::ElecState* p_elec_in, ModulePW::PW_Basis_K* p_wfcpw_in, + psi::Psi>* p_psi_in, pseudopot_cell_vnl* p_ppcell_in); + ~EleCond(){}; + + /** + * @brief calculate Onsager coefficients Lmn(\omega) and conductivities with Kubo-Greenwood formula + * + * @param fwhmin FWHM for delta function + * @param smear_type 1: Gaussian, 2: Lorentzian + * @param wcut cutoff \omega for Lmn(\omega) + * @param dw_in \omega step + * @param dt_in time step + * @param nonlocal whether to include the nonlocal potential corrections for velocity operator + * @param wg wg(ik,ib) occupation for the ib-th band in the ik-th kpoint + */ + void KG(const int& smear_type, const double& fwhmin, const double& wcut, const double& dw_in, const double& dt_in, + const bool& nonlocal, ModuleBase::matrix& wg); + + protected: + pseudopot_cell_vnl* p_ppcell = nullptr; ///< pointer to the pseudopotential + UnitCell* p_ucell = nullptr; ///< pointer to the unit cell + ModulePW::PW_Basis_K* p_wfcpw = nullptr; ///< pointer to the plane wave basis + K_Vectors* p_kv = nullptr; ///< pointer to the k vectors + elecstate::ElecState* p_elec = nullptr; ///< pointer to the electronic state + psi::Psi>* p_psi = nullptr; ///< pointer to the wavefunction + + protected: + /** + * @brief calculate the response function Cmn(t) for currents + * + * @param ik k point + * @param nt number of steps of time + * @param dt time step + * @param decut ignore dE which is larger than decut + * @param wg wg(ik,ib) occupation for the ib-th band in the ik-th kpoint + * @param velop velocity operator + * @param ct11 C11(t) + * @param ct12 C12(t) + * @param ct22 C22(t) + */ + void jjresponse_ks(const int ik, const int nt, const double dt, const double decut, ModuleBase::matrix& wg, + hamilt::Velocity& velop, double* ct11, double* ct12, double* ct22); + /** + * @brief Calculate the conductivity using the response function + * + * @param nt number of time steps + * @param dt time step + * @param smear_type smearing type 1: gaussian, 2: lorentzian + * @param fwhmin full width at half maximum of the smearing function + * @param wcut cutoff frequency + * @param dw_in frequency step + * @param ct11 C11 component of the response function + * @param ct12 C12 component of the response function + * @param ct22 C22 component of the response function + */ + void calcondw(const int nt, const double dt, const int& smear_type, const double fwhmin, const double wcut, + const double dw_in, double* ct11, double* ct12, double* ct22); +}; + +#endif // ELECOND_H \ No newline at end of file diff --git a/source/module_hamilt_pw/hamilt_stodft/CMakeLists.txt b/source/module_hamilt_pw/hamilt_stodft/CMakeLists.txt index 793f382dfb..9882a4d868 100644 --- a/source/module_hamilt_pw/hamilt_stodft/CMakeLists.txt +++ b/source/module_hamilt_pw/hamilt_stodft/CMakeLists.txt @@ -5,6 +5,9 @@ list(APPEND hamilt_stodft_srcs sto_func.cpp sto_forces.cpp sto_stress_pw.cpp + sto_tool.cpp + sto_elecond.cpp + sto_dos.cpp ) add_library( @@ -16,3 +19,7 @@ add_library( if(ENABLE_COVERAGE) add_coverage(hamilt_stodft) endif() + +if(BUILD_TESTING) + add_subdirectory(test) +endif() \ No newline at end of file diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_dos.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_dos.cpp new file mode 100644 index 0000000000..a125f5cef8 --- /dev/null +++ b/source/module_hamilt_pw/hamilt_stodft/sto_dos.cpp @@ -0,0 +1,238 @@ +#include "sto_dos.h" + +#include "module_base/timer.h" +#include "module_base/tool_title.h" +#include "sto_tool.h" + +Sto_DOS::Sto_DOS(ModulePW::PW_Basis_K* p_wfcpw_in, K_Vectors* p_kv_in, elecstate::ElecState* p_elec_in, + psi::Psi>* p_psi_in, hamilt::Hamilt>* p_hamilt_in, + hsolver::HSolverPW_SDFT* p_hsol_in, Stochastic_WF* p_stowf_in) +{ + this->p_wfcpw = p_wfcpw_in; + this->p_kv = p_kv_in; + this->p_elec = p_elec_in; + this->p_psi = p_psi_in; + this->p_hamilt = p_hamilt_in; + this->p_hsol = p_hsol_in; + this->p_stowf = p_stowf_in; + this->nbands_ks = p_psi_in->get_nbands(); + this->nbands_sto = p_stowf_in->nchi; +} +void Sto_DOS::decide_param(const int& dos_nche, const double& emin_sto, const double& emax_sto, const bool& dos_setemin, + const bool& dos_setemax, const double& dos_emin_ev, const double& dos_emax_ev, + const double& dos_scale) +{ + this->dos_nche = dos_nche; + check_che(this->dos_nche, emin_sto, emax_sto, this->nbands_sto, this->p_kv, this->p_stowf, this->p_hamilt, + this->p_hsol); + if (dos_setemax) + { + this->emax = dos_emax_ev; + } + else + { + this->emax = p_hsol->stoiter.stohchi.Emax * ModuleBase::Ry_to_eV; + } + if (dos_setemin) + { + this->emin = dos_emin_ev; + } + else + { + this->emin = p_hsol->stoiter.stohchi.Emin * ModuleBase::Ry_to_eV; + } + + if (!dos_setemax && !dos_setemin) + { + double delta = (emax - emin) * dos_scale; + this->emax = emax + delta / 2.0; + this->emin = emin - delta / 2.0; + } +} + +void Sto_DOS::caldos(const double sigmain, const double de, const int npart) +{ + ModuleBase::TITLE("Sto_DOS", "caldos"); + ModuleBase::timer::tick("Sto_DOS", "caldos"); + std::cout << "=========================" << std::endl; + std::cout << "###Calculating Dos....###" << std::endl; + std::cout << "=========================" << std::endl; + ModuleBase::Chebyshev che(dos_nche); + const int nk = p_kv->get_nks(); + Stochastic_Iter& stoiter = p_hsol->stoiter; + Stochastic_hchi& stohchi = stoiter.stohchi; + const int npwx = p_wfcpw->npwk_max; + + std::vector spolyv; + std::vector> allorderchi; + if (stoiter.method == 1) + { + spolyv.resize(dos_nche, 0); + } + else + { + spolyv.resize(dos_nche * dos_nche, 0); + int nchip_new = ceil((double)this->p_stowf->nchip_max / npart); + allorderchi.resize(nchip_new * npwx * dos_nche); + } + ModuleBase::timer::tick("Sto_DOS", "Tracepoly"); + std::cout << "1. TracepolyA:" << std::endl; + for (int ik = 0; ik < nk; ik++) + { + std::cout << "ik: " << ik + 1 << std::endl; + if (nk > 1) + { + this->p_hamilt->updateHk(ik); + } + stohchi.current_ik = ik; + const int npw = p_kv->ngk[ik]; + const int nchipk = this->p_stowf->nchip[ik]; + + std::complex* pchi; + if (GlobalV::NBANDS > 0) + { + p_stowf->chiortho->fix_k(ik); + pchi = p_stowf->chiortho->get_pointer(); + } + else + { + p_stowf->chi0->fix_k(ik); + pchi = p_stowf->chi0->get_pointer(); + } + if (stoiter.method == 1) + { + che.tracepolyA(&stohchi, &Stochastic_hchi::hchi_norm, pchi, npw, npwx, nchipk); + for (int i = 0; i < dos_nche; ++i) + { + spolyv[i] += che.polytrace[i] * p_kv->wk[ik] / 2; + } + } + else + { + int N = dos_nche; + double kweight = p_kv->wk[ik] / 2; + char trans = 'T'; + char normal = 'N'; + double one = 1; + for (int ipart = 0; ipart < npart; ++ipart) + { + int nchipk_new = nchipk / npart; + int start_nchipk = ipart * nchipk_new + nchipk % npart; + if (ipart < nchipk % npart) + { + nchipk_new++; + start_nchipk = ipart * nchipk_new; + } + ModuleBase::GlobalFunc::ZEROS(allorderchi.data(), nchipk_new * npwx * dos_nche); + std::complex* tmpchi = pchi + start_nchipk * npwx; + che.calpolyvec_complex(&stohchi, &Stochastic_hchi::hchi_norm, tmpchi, allorderchi.data(), npw, npwx, + nchipk_new); + double* vec_all = (double*)allorderchi.data(); + int LDA = npwx * nchipk_new * 2; + int M = npwx * nchipk_new * 2; + dgemm_(&trans, &normal, &N, &N, &M, &kweight, vec_all, &LDA, vec_all, &LDA, &one, spolyv.data(), &N); + } + } + } + + allorderchi.clear(); + allorderchi.shrink_to_fit(); + + std::ofstream ofsdos; + int ndos = int((emax - emin) / de) + 1; + stoiter.stofunc.sigma = sigmain / ModuleBase::Ry_to_eV; + ModuleBase::timer::tick("Sto_DOS", "Tracepoly"); + + std::cout << "2. Dos:" << std::endl; + ModuleBase::timer::tick("Sto_DOS", "DOS Loop"); + int n10 = ndos / 10; + int percent = 10; + std::vector sto_dos(ndos); + std::vector ks_dos(ndos); + std::vector error(ndos); + for (int ie = 0; ie < ndos; ++ie) + { + double tmpks = 0; + double tmpsto = 0; + stoiter.stofunc.targ_e = (emin + ie * de) / ModuleBase::Ry_to_eV; + if (stoiter.method == 1) + { + che.calcoef_real(&stoiter.stofunc, &Sto_Func::ngauss); + tmpsto = BlasConnector::dot(dos_nche, che.coef_real, 1, spolyv.data(), 1); + } + else + { + che.calcoef_real(&stoiter.stofunc, &Sto_Func::nroot_gauss); + tmpsto = stoiter.vTMv(che.coef_real, spolyv.data(), dos_nche); + } + if (GlobalV::NBANDS > 0) + { + for (int ik = 0; ik < nk; ++ik) + { + double* en = &(this->p_elec->ekb(ik, 0)); + for (int ib = 0; ib < GlobalV::NBANDS; ++ib) + { + tmpks += stoiter.stofunc.gauss(en[ib]) * p_kv->wk[ik] / 2; + } + } + } + tmpks /= GlobalV::NPROC_IN_POOL; + + double tmperror = 0; + if (stoiter.method == 1) + { + tmperror = che.coef_real[dos_nche - 1] * spolyv[dos_nche - 1]; + } + else + { + const int norder = dos_nche; + double last_coef = che.coef_real[norder - 1]; + double last_spolyv = spolyv[norder * norder - 1]; + tmperror = last_coef + * (BlasConnector::dot(norder, che.coef_real, 1, spolyv.data() + norder * (norder - 1), 1) + + BlasConnector::dot(norder, che.coef_real, 1, spolyv.data() + norder - 1, norder) + - last_coef * last_spolyv); + } + + if (ie % n10 == n10 - 1) + { + std::cout << percent << "%" + << " "; + percent += 10; + } + sto_dos[ie] = tmpsto; + ks_dos[ie] = tmpks; + error[ie] = tmperror; + } +#ifdef __MPI + MPI_Allreduce(MPI_IN_PLACE, ks_dos.data(), ndos, MPI_DOUBLE, MPI_SUM, STO_WORLD); + MPI_Allreduce(MPI_IN_PLACE, sto_dos.data(), ndos, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, error.data(), ndos, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#endif + if (GlobalV::MY_RANK == 0) + { + std::string dosfile = GlobalV::global_out_dir + "DOS1_smearing.dat"; + ofsdos.open(dosfile.c_str()); + double maxerror = 0; + double sum = 0; + ofsdos << std::setw(8) << "## E(eV) " << std::setw(20) << "dos(eV^-1)" << std::setw(20) << "sum" + << std::setw(20) << "Error(eV^-1)" << std::endl; + for (int ie = 0; ie < ndos; ++ie) + { + double tmperror = 2.0 * std::abs(error[ie]); + if (maxerror < tmperror) + maxerror = tmperror; + double dos = 2.0 * (ks_dos[ie] + sto_dos[ie]) / ModuleBase::Ry_to_eV; + sum += dos; + ofsdos << std::setw(8) << emin + ie * de << std::setw(20) << dos << std::setw(20) << sum * de + << std::setw(20) << tmperror << std::endl; + } + std::cout << std::endl; + std::cout << "Finish DOS" << std::endl; + std::cout << std::scientific << "DOS max absolute Chebyshev Error: " << maxerror << std::endl; + ofsdos.close(); + } + ModuleBase::timer::tick("Sto_DOS", "DOS Loop"); + ModuleBase::timer::tick("Sto_DOS", "caldos"); + return; +} diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_dos.h b/source/module_hamilt_pw/hamilt_stodft/sto_dos.h new file mode 100644 index 0000000000..407ec3cad1 --- /dev/null +++ b/source/module_hamilt_pw/hamilt_stodft/sto_dos.h @@ -0,0 +1,51 @@ +#ifndef STO_DOS +#define STO_DOS +#include "module_hamilt_pw/hamilt_stodft/sto_wf.h" +#include "module_hsolver/hsolver_pw_sdft.h" + +class Sto_DOS +{ + public: + Sto_DOS(ModulePW::PW_Basis_K* p_wfcpw_in, K_Vectors* p_kv_in, elecstate::ElecState* p_elec_in, + psi::Psi>* p_psi_in, hamilt::Hamilt>* p_hamilt_in, + hsolver::HSolverPW_SDFT* p_hsol_in, Stochastic_WF* p_stowf_in); + /** + * @brief decide the parameters for the DOS calculation + * + * @param dos_nche Number of Chebyshev orders + * @param emin_sto Emin input for sdft + * @param emax_sto Emax input for sdft + * @param dos_setemin whether to set the minimum energy + * @param dos_setemax whether to set the maximum energy + * @param dos_emin_ev Emin input for DOS + * @param dos_emax_ev Emax input for DOS + * @param dos_scale dos_scale input for DOS + */ + void decide_param(const int& dos_nche, const double& emin_sto, const double& emax_sto, const bool& dos_setemin, + const bool& dos_setemax, const double& dos_emin_ev, const double& dos_emax_ev, + const double& dos_scale); + /** + * @brief Calculate DOS using stochastic wavefunctions + * + * @param sigmain sigma for the gaussian broadening + * @param de energy step + * @param npart number of parts to reduce the memory usage + */ + void caldos(const double sigmain, const double de, const int npart); + + protected: + int nbands_ks = 0; ///< number of KS bands + int nbands_sto = 0; ///< number of stochastic bands + int dos_nche = 0; ///< number of Chebyshev orders for DOS + double emax = 0.0; ///< maximum energy + double emin = 0.0; ///< minimum energy + ModulePW::PW_Basis_K* p_wfcpw = nullptr; ///< pointer to the plane wave basis + K_Vectors* p_kv = nullptr; ///< pointer to the k vectors + elecstate::ElecState* p_elec = nullptr; ///< pointer to the electronic state + psi::Psi>* p_psi = nullptr; ///< pointer to the wavefunction + hamilt::Hamilt>* p_hamilt; ///< pointer to the Hamiltonian + hsolver::HSolverPW_SDFT* p_hsol = nullptr; ///< pointer to the Hamiltonian solver + Stochastic_WF* p_stowf = nullptr; ///< pointer to the stochastic wavefunctions +}; + +#endif // STO_DOS \ No newline at end of file diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_elecond.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_elecond.cpp new file mode 100644 index 0000000000..aa4f15e252 --- /dev/null +++ b/source/module_hamilt_pw/hamilt_stodft/sto_elecond.cpp @@ -0,0 +1,814 @@ +#include "sto_elecond.h" + +#include "module_base/complexmatrix.h" +#include "module_base/constants.h" +#include "module_base/memory.h" +#include "module_base/timer.h" +#include "module_base/vector3.h" +#include "sto_tool.h" + +#include + +#define TWOSQRT2LN2 2.354820045030949 // FWHM = 2sqrt(2ln2) * \sigma +#define FACTOR 1.839939223835727e7 + +Sto_EleCond::Sto_EleCond(UnitCell* p_ucell_in, K_Vectors* p_kv_in, elecstate::ElecState* p_elec_in, + ModulePW::PW_Basis_K* p_wfcpw_in, psi::Psi>* p_psi_in, + pseudopot_cell_vnl* p_ppcell_in, hamilt::Hamilt>* p_hamilt_in, + hsolver::HSolverPW_SDFT* p_hsol_in, Stochastic_WF* p_stowf_in) + : EleCond(p_ucell_in, p_kv_in, p_elec_in, p_wfcpw_in, p_psi_in, p_ppcell_in) +{ + this->p_hamilt = p_hamilt_in; + this->p_hsol = p_hsol_in; + this->p_stowf = p_stowf_in; + this->nbands_ks = p_psi_in->get_nbands(); + this->nbands_sto = p_stowf_in->nchi; +} + +void Sto_EleCond::decide_nche(const double dt, int& nbatch, const double cond_thr, const int& fd_nche, double try_emin, + double try_emax) +{ + int nche_guess = 1000; + ModuleBase::Chebyshev chet(nche_guess); + Stochastic_Iter& stoiter = p_hsol->stoiter; + const double mu = this->p_elec->eferm.ef; + stoiter.stofunc.mu = mu; + // try to find nbatch + if (nbatch == 0) + { + for (int test_nbatch = 128; test_nbatch >= 1; test_nbatch /= 2) + { + nbatch = test_nbatch; + stoiter.stofunc.t = 0.5 * dt * nbatch; + chet.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::n_sin); + double minerror = std::abs(chet.coef_complex[nche_guess - 1] / chet.coef_complex[0]); + if (minerror < cond_thr) + { + for (int i = 1; i < nche_guess; ++i) + { + double error = std::abs(chet.coef_complex[i] / chet.coef_complex[0]); + if (error < cond_thr) + { + // To make nche to around 100 + nbatch = ceil(double(test_nbatch) / i * 100.0); + std::cout << "set cond_dtbatch to " << nbatch << std::endl; + break; + } + } + break; + } + } + } + + // first try to find nche + stoiter.stofunc.t = 0.5 * dt * nbatch; + auto getnche = [&](int& nche) { + chet.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::n_sin); + for (int i = 1; i < nche_guess; ++i) + { + double error = std::abs(chet.coef_complex[i] / chet.coef_complex[0]); + if (error < cond_thr) + { + nche = i + 1; + break; + } + } + }; + int nche_old = 0; + getnche(nche_old); + + int nche_new = 0; +loop: + // re-set Emin & Emax both in stohchi & stofunc + check_che(std::max(nche_old * 2, fd_nche), try_emin, try_emax, this->nbands_sto, this->p_kv, this->p_stowf, + this->p_hamilt, this->p_hsol); + + // second try to find nche with new Emin & Emax + getnche(nche_new); + + if (nche_new > nche_old * 2) + { + nche_old = nche_new; + try_emin = stoiter.stohchi.Emin; + try_emax = stoiter.stohchi.Emax; + goto loop; + } + + std::cout << "set N order of Chebyshev for KG as " << nche_new << std::endl; + std::ofstream cheofs("Chebycoef"); + for (int i = 1; i < nche_guess; ++i) + { + double error = std::abs(chet.coef_complex[i] / chet.coef_complex[0]); + cheofs << std::setw(5) << i << std::setw(20) << error << std::endl; + } + cheofs.close(); + + if (nche_new >= 1000) + { + ModuleBase::WARNING_QUIT("ESolver_SDFT_PW", "N order of Chebyshev for KG will be larger than 1000!"); + } + + this->cond_nche = nche_new; + this->fd_nche = fd_nche; +} + +void Sto_EleCond::cal_jmatrix(const psi::Psi>& kspsi_all, + const psi::Psi>& vkspsi, const double* en, const double* en_all, + std::complex* leftfact, std::complex* rightfact, + const psi::Psi>& leftchi, psi::Psi>& rightchi, + psi::Psi>& left_hchi, psi::Psi>& batch_vchi, + psi::Psi>& batch_vhchi, +#ifdef __MPI + psi::Psi>& chi_all, psi::Psi>& hchi_all, + void* gatherinfo_ks, void* gatherinfo_sto, +#endif + const int& bsize_psi, std::vector>& j1, + std::vector>& j2, hamilt::Velocity& velop, const int& ik, + const std::complex& factor, const int bandinfo[6]) +{ + ModuleBase::timer::tick("Sto_EleCond", "cal_jmatrix"); + const char transa = 'C'; + const char transb = 'N'; + const std::complex float_factor = static_cast>(factor); + const std::complex conjfactor = std::conj(float_factor); + const float mu = static_cast(this->p_elec->eferm.ef); + const std::complex zero = 0.0; + const int npwx = this->p_wfcpw->npwk_max; + const int npw = this->p_wfcpw->npwk[ik]; + const int ndim = 3; + const int perbands_ks = bandinfo[0]; + const int perbands_sto = bandinfo[1]; + const int perbands = bandinfo[2]; + const int allbands_ks = bandinfo[3]; + const int allbands_sto = bandinfo[4]; + const int allbands = bandinfo[5]; + const int dim_jmatrix = perbands_ks * allbands_sto + perbands_sto * allbands; + Stochastic_Iter& stoiter = p_hsol->stoiter; + + psi::Psi> right_hchi(1, perbands_sto, npwx, p_kv->ngk.data()); + psi::Psi> f_rightchi(1, perbands_sto, npwx, p_kv->ngk.data()); + psi::Psi> f_right_hchi(1, perbands_sto, npwx, p_kv->ngk.data()); + + stoiter.stohchi.hchi(leftchi.get_pointer(), left_hchi.get_pointer(), perbands_sto); + stoiter.stohchi.hchi(rightchi.get_pointer(), right_hchi.get_pointer(), perbands_sto); + convert_psi(rightchi, f_rightchi); + convert_psi(right_hchi, f_right_hchi); + right_hchi.resize(1, 1, 1); + + psi::Psi>* rightchi_all = &f_rightchi; + psi::Psi>* righthchi_all = &f_right_hchi; + std::vector> vec_rightf_all; + std::complex* rightf_all = rightfact; +#ifdef __MPI + info_gatherv* ks_fact = static_cast(gatherinfo_ks); + info_gatherv* sto_npwx = static_cast(gatherinfo_sto); + rightchi_all = gatherchi(f_rightchi, chi_all, npwx, sto_npwx->nrecv, sto_npwx->displs, perbands_sto); + righthchi_all = gatherchi(f_right_hchi, hchi_all, npwx, sto_npwx->nrecv, sto_npwx->displs, perbands_sto); + if (GlobalV::NSTOGROUP > 1 && rightfact != nullptr) + { + vec_rightf_all.resize(allbands_ks); + rightf_all = vec_rightf_all.data(); + MPI_Allgatherv(rightfact, perbands_ks, MPI_DOUBLE_COMPLEX, rightf_all, ks_fact->nrecv, ks_fact->displs, + MPI_DOUBLE_COMPLEX, PARAPW_WORLD); + } +#endif + + psi::Psi> f_batch_vchi(1, bsize_psi * ndim, npwx, p_kv->ngk.data()); + psi::Psi> f_batch_vhchi(1, bsize_psi * ndim, npwx, p_kv->ngk.data()); + std::vector> tmpj(ndim * allbands_sto * perbands_sto); + + // 1. (<\psi|J|\chi>)^T + // (allbands_sto, perbands_ks) + if (perbands_ks > 0) + { + for (int id = 0; id < ndim; ++id) + { + const int idnb = id * perbands_ks; + const int jbais = 0; + std::complex* j1mat = &j1[id * dim_jmatrix]; + std::complex* j2mat = &j2[id * dim_jmatrix]; + //(<\psi|v|\chi>)^T + cgemm_(&transa, &transb, &allbands_sto, &perbands_ks, &npw, &conjfactor, rightchi_all->get_pointer(), &npwx, + &vkspsi(idnb, 0), &npwx, &zero, j1mat, &allbands_sto); + + //(<\psi|Hv|\chi>)^T + // for(int i = 0 ; i < perbands_ks ; ++i) + // { + // double* evmat = &j2(id, 0 + i * allbands_sto); + // double* vmat = &j1(id, 0 + i * allbands_sto); + // double ei = en[i]; + // for(int j = 0 ; j < allbands_sto ; ++j) + // { + // evmat[j] = vmat[j] * ei; + // } + // } + + //(<\psi|vH|\chi>)^T + cgemm_(&transa, &transb, &allbands_sto, &perbands_ks, &npw, &conjfactor, righthchi_all->get_pointer(), + &npwx, &vkspsi(idnb, 0), &npwx, &zero, j2mat, &allbands_sto); + } + } + + int remain = perbands_sto; + int startnb = 0; + while (remain > 0) + { + int tmpnb = std::min(remain, bsize_psi); + // v|\chi> + velop.act(&leftchi, tmpnb, &leftchi(0, startnb, 0), batch_vchi.get_pointer()); + convert_psi(batch_vchi, f_batch_vchi); + // v|H\chi> + velop.act(&leftchi, tmpnb, &left_hchi(0, startnb, 0), batch_vhchi.get_pointer()); + convert_psi(batch_vhchi, f_batch_vhchi); + // 2. <\chi|J|\psi> + // (perbands_sto, allbands_ks) + if (allbands_ks > 0) + { + for (int id = 0; id < ndim; ++id) + { + const int idnb = id * tmpnb; + const int jbais = perbands_ks * allbands_sto + startnb; + std::complex* j1mat = &j1[id * dim_jmatrix + jbais]; + std::complex* j2mat = &j2[id * dim_jmatrix + jbais]; + //<\chi|v|\psi> + cgemm_(&transa, &transb, &tmpnb, &allbands_ks, &npw, &float_factor, &f_batch_vchi(idnb, 0), &npwx, + kspsi_all.get_pointer(), &npwx, &zero, j1mat, &perbands_sto); + + //<\chi|vH|\psi> = \epsilon * <\chi|v|\psi> + // for(int i = 0 ; i < allbands_ks ; ++i) + // { + // double* evmat = &j2(id, jbais + i * allbands_sto); + // double* vmat = &j1(id, jbais + i * allbands_sto); + // double ei = en_all[i]; + // for(int j = 0 ; j < tmpnb ; ++j) + // { + // evmat[j] = vmat[j] * en_all[i]; + // } + // } + + //<\chi|Hv|\psi> + cgemm_(&transa, &transb, &tmpnb, &allbands_ks, &npw, &float_factor, &f_batch_vhchi(idnb, 0), &npwx, + kspsi_all.get_pointer(), &npwx, &zero, j2mat, &perbands_sto); + } + } + + // 3. <\chi|J|\chi> + // (perbands_sto, allbands_sto) + for (int id = 0; id < ndim; ++id) + { + const int idnb = id * tmpnb; + const int jbais = perbands_ks * allbands_sto + perbands_sto * allbands_ks + startnb; + std::complex* j1mat = &j1[id * dim_jmatrix + jbais]; + std::complex* j2mat = &j2[id * dim_jmatrix + jbais]; + std::complex* tmpjmat = &tmpj[id * allbands_sto * perbands_sto + startnb]; + //<\chi|v|\chi> + cgemm_(&transa, &transb, &tmpnb, &allbands_sto, &npw, &float_factor, &f_batch_vchi(idnb, 0), &npwx, + rightchi_all->get_pointer(), &npwx, &zero, j1mat, &perbands_sto); + + //<\chi|Hv|\chi> + cgemm_(&transa, &transb, &tmpnb, &allbands_sto, &npw, &float_factor, &f_batch_vhchi(idnb, 0), &npwx, + rightchi_all->get_pointer(), &npwx, &zero, j2mat, &perbands_sto); + + //<\chi|vH|\chi> + cgemm_(&transa, &transb, &tmpnb, &allbands_sto, &npw, &float_factor, &f_batch_vchi(idnb, 0), &npwx, + righthchi_all->get_pointer(), &npwx, &zero, tmpjmat, &perbands_sto); + } + + remain -= tmpnb; + startnb += tmpnb; + if (remain == 0) + break; + } + + for (int id = 0; id < ndim; ++id) + { + for (int i = 0; i < perbands_ks; ++i) + { + const float ei = static_cast(en[i]); + const int jst = i * allbands_sto; + std::complex* j2mat = j2.data() + id * dim_jmatrix + jst; + std::complex* j1mat = j1.data() + id * dim_jmatrix + jst; + if (leftfact == nullptr) + { + for (int j = 0; j < allbands_sto; ++j) + { + j2mat[j] = 0.5f * j2mat[j] + (0.5f * ei - mu) * j1mat[j]; + } + } + else + { + const std::complex jfac = static_cast>(leftfact[i]); + for (int j = 0; j < allbands_sto; ++j) + { + j2mat[j] = jfac * (0.5f * j2mat[j] + (0.5f * ei - mu) * j1mat[j]); + j1mat[j] *= jfac; + } + } + } + + for (int i = 0; i < allbands_ks; ++i) + { + const float ei = static_cast(en_all[i]); + const int jst = perbands_ks * allbands_sto + i * perbands_sto; + std::complex* j2mat = j2.data() + id * dim_jmatrix + jst; + std::complex* j1mat = j1.data() + id * dim_jmatrix + jst; + if (rightfact == nullptr) + { + for (int j = 0; j < perbands_sto; ++j) + { + j2mat[j] = 0.5f * j2mat[j] + (0.5f * ei - mu) * j1mat[j]; + } + } + else + { + const std::complex jfac = static_cast>(rightf_all[i]); + for (int j = 0; j < perbands_sto; ++j) + { + j2mat[j] = jfac * (0.5f * j2mat[j] + (0.5f * ei - mu) * j1mat[j]); + j1mat[j] *= jfac; + } + } + } + + const int jst = perbands_ks * allbands_sto + perbands_sto * allbands_ks; + const int ed = dim_jmatrix - jst; + std::complex* j2mat = j2.data() + id * dim_jmatrix + jst; + std::complex* j1mat = j1.data() + id * dim_jmatrix + jst; + std::complex* tmpjmat = tmpj.data() + id * allbands_sto * perbands_sto; + + for (int j = 0; j < ed; ++j) + { + j2mat[j] = 0.5f * (j2mat[j] + tmpjmat[j]) - mu * j1mat[j]; + } + } + +#ifdef __MPI + MPI_Allreduce(MPI_IN_PLACE, j1.data(), ndim * dim_jmatrix, MPI_COMPLEX, MPI_SUM, POOL_WORLD); + MPI_Allreduce(MPI_IN_PLACE, j2.data(), ndim * dim_jmatrix, MPI_COMPLEX, MPI_SUM, POOL_WORLD); +#endif + + ModuleBase::timer::tick("Sto_EleCond", "cal_jmatrix"); + + return; +} + +void Sto_EleCond::sKG(const int& smear_type, const double& fwhmin, const double& wcut, const double& dw_in, + const double& dt_in, const bool& nonlocal, const int& nbatch, const int& npart_sto) +{ + ModuleBase::TITLE("Sto_EleCond", "sKG"); + ModuleBase::timer::tick("Sto_EleCond", "sKG"); + std::cout << "Calculating conductivity...." << std::endl; + // if (GlobalV::NSTOGROUP > 1) + // { + // ModuleBase::WARNING_QUIT("ESolver_SDFT_PW", "sKG is not supported in parallel!"); + // } + + //------------------------------------------------------------------ + // Init + //------------------------------------------------------------------ + // Parameters + int nw = ceil(wcut / dw_in); + double dw = dw_in / ModuleBase::Ry_to_eV; // converge unit in eV to Ry + double sigma = fwhmin / TWOSQRT2LN2 / ModuleBase::Ry_to_eV; + double gamma = fwhmin / 2.0 / ModuleBase::Ry_to_eV; + double dt = dt_in; // unit in a.u., 1 a.u. = 4.837771834548454e-17 s + const double expfactor = 18.42; // exp(-18.42) = 1e-8 + int nt = 0; // set nt empirically + if (smear_type == 1) + { + nt = ceil(sqrt(2 * expfactor) / sigma / dt); + } + else if (smear_type == 2) + { + nt = ceil(expfactor / gamma / dt); + } + else + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::calcondw", "smear_type should be 0 or 1"); + } + std::cout << "nw: " << nw << " ; dw: " << dw * ModuleBase::Ry_to_eV << " eV" << std::endl; + std::cout << "nt: " << nt << " ; dt: " << dt << " a.u.(ry^-1)" << std::endl; + assert(nw >= 1); + assert(nt >= 1); + const int ndim = 3; + const int nk = p_kv->get_nks(); + const int npwx = p_wfcpw->npwk_max; + const double tpiba = p_wfcpw->tpiba; + psi::Psi>* stopsi; + if (this->nbands_ks > 0) + { + stopsi = p_stowf->chiortho; + // clean memories //Note shchi is different from \sqrt(fH_here)|chi>, since veffs are different + p_stowf->shchi->resize(1, 1, 1); + p_stowf->chi0->resize(1, 1, 1); // clean memories + } + else + { + stopsi = p_stowf->chi0; + p_stowf->shchi->resize(1, 1, 1); // clean memories + } + const double dEcut = (wcut + fwhmin) / ModuleBase::Ry_to_eV; + + // response funtion + std::vector ct11(nt, 0); + std::vector ct12(nt, 0); + std::vector ct22(nt, 0); + + // Init Chebyshev + ModuleBase::Chebyshev che(fd_nche); + ModuleBase::Chebyshev chet(cond_nche); + ModuleBase::Chebyshev chemt(cond_nche); + Stochastic_Iter& stoiter = p_hsol->stoiter; + Stochastic_hchi& stohchi = stoiter.stohchi; + + //------------------------------------------------------------------ + // Calculate + //------------------------------------------------------------------ + + // Prepare Chebyshev coefficients for exp(-i H/\hbar t) + const double mu = this->p_elec->eferm.ef; + stoiter.stofunc.mu = mu; + stoiter.stofunc.t = 0.5 * dt * nbatch; + chet.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::nsin); + chemt.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::n_sin); + std::vector> batchcoef, batchmcoef; + if (nbatch > 1) + { + batchcoef.resize(cond_nche * nbatch); + std::complex* tmpcoef = batchcoef.data() + (nbatch - 1) * cond_nche; + batchmcoef.resize(cond_nche * nbatch); + std::complex* tmpmcoef = batchmcoef.data() + (nbatch - 1) * cond_nche; + for (int i = 0; i < cond_nche; ++i) + { + tmpcoef[i] = chet.coef_complex[i]; + tmpmcoef[i] = chemt.coef_complex[i]; + } + for (int ib = 0; ib < nbatch - 1; ++ib) + { + tmpcoef = batchcoef.data() + ib * cond_nche; + tmpmcoef = batchmcoef.data() + ib * cond_nche; + stoiter.stofunc.t = 0.5 * dt * (ib + 1); + chet.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::nsin); + chemt.calcoef_pair(&stoiter.stofunc, &Sto_Func::ncos, &Sto_Func::n_sin); + for (int i = 0; i < cond_nche; ++i) + { + tmpcoef[i] = chet.coef_complex[i]; + tmpmcoef[i] = chemt.coef_complex[i]; + } + } + stoiter.stofunc.t = 0.5 * dt * nbatch; + } + + // ik loop + ModuleBase::timer::tick("Sto_EleCond", "kloop"); + hamilt::Velocity velop(p_wfcpw, p_kv->isk.data(), p_ppcell, p_ucell, nonlocal); + for (int ik = 0; ik < nk; ++ik) + { + velop.init(ik); + stopsi->fix_k(ik); + p_psi->fix_k(ik); + if (nk > 1) + { + this->p_hamilt->updateHk(ik); + } + stoiter.stohchi.current_ik = ik; + const int npw = p_kv->ngk[ik]; + + // get allbands_ks + int cutib0 = 0; + if (this->nbands_ks > 0) + { + double Emax_KS = std::max(stoiter.stofunc.Emin, this->p_elec->ekb(ik, this->nbands_ks - 1)); + for (cutib0 = this->nbands_ks - 1; cutib0 >= 0; --cutib0) + { + if (Emax_KS - this->p_elec->ekb(ik, cutib0) > dEcut) + { + break; + } + } + ++cutib0; + double Emin_KS = (cutib0 < this->nbands_ks) ? this->p_elec->ekb(ik, cutib0) : stoiter.stofunc.Emin; + double dE = stoiter.stofunc.Emax - Emin_KS + wcut / ModuleBase::Ry_to_eV; + std::cout << "Emin_KS(" << cutib0 + 1 << "): " << Emin_KS * ModuleBase::Ry_to_eV + << " eV; Emax: " << stoiter.stofunc.Emax * ModuleBase::Ry_to_eV + << " eV; Recommended max dt: " << 2 * M_PI / dE << " a.u." << std::endl; + } + else + { + double dE = stoiter.stofunc.Emax - stoiter.stofunc.Emin + wcut / ModuleBase::Ry_to_eV; + std::cout << "Emin: " << stoiter.stofunc.Emin * ModuleBase::Ry_to_eV + << " eV; Emax: " << stoiter.stofunc.Emax * ModuleBase::Ry_to_eV + << " eV; Recommended max dt: " << 2 * M_PI / dE << " a.u." << std::endl; + } + // Parallel for bands + int allbands_ks = this->nbands_ks - cutib0; + parallel_distribution paraks(allbands_ks, GlobalV::NSTOGROUP, GlobalV::MY_STOGROUP); + int perbands_ks = paraks.num_per; + int ib0_ks = paraks.start; + ib0_ks += this->nbands_ks - allbands_ks; + int perbands_sto = this->p_stowf->nchip[ik]; + int perbands = perbands_sto + perbands_ks; + int allbands_sto = perbands_sto; + int allbands = perbands; +#ifdef __MPI + MPI_Allreduce(&perbands, &allbands, 1, MPI_INT, MPI_SUM, PARAPW_WORLD); + allbands_sto = allbands - allbands_ks; + info_gatherv ks_fact(perbands_ks, GlobalV::NSTOGROUP, 1, PARAPW_WORLD); + info_gatherv sto_npwx(perbands_sto, GlobalV::NSTOGROUP, npwx, PARAPW_WORLD); +#endif + const int bandsinfo[6]{perbands_ks, perbands_sto, perbands, allbands_ks, allbands_sto, allbands}; + double* en_all = nullptr; + std::vector en; + if (allbands_ks > 0) + { + en_all = &(this->p_elec->ekb(ik, this->nbands_ks - allbands_ks)); + } + if (perbands_ks > 0) + { + en.resize(perbands_ks); + for (int ib = 0; ib < perbands_ks; ++ib) + { + en[ib] = this->p_elec->ekb(ik, ib0_ks + ib); + } + } + + //----------------------------------------------------------- + // ks conductivity + //----------------------------------------------------------- + if (GlobalV::MY_STOGROUP == 0 && allbands_ks > 0) + { + jjresponse_ks(ik, nt, dt, dEcut, this->p_elec->wg, velop, ct11.data(), ct12.data(), ct22.data()); + } + + //----------------------------------------------------------- + // sto conductivity + //----------------------------------------------------------- + //------------------- allocate ------------------------- + size_t ks_memory_cost = perbands_ks * npwx * sizeof(std::complex); + psi::Psi> kspsi(1, perbands_ks, npwx, p_kv->ngk.data()); + psi::Psi> vkspsi(1, perbands_ks * ndim, npwx, p_kv->ngk.data()); + std::vector> expmtmf_fact(perbands_ks), expmtf_fact(perbands_ks); + psi::Psi> f_kspsi(1, perbands_ks, npwx, p_kv->ngk.data()); + ModuleBase::Memory::record("SDFT::kspsi", ks_memory_cost); + psi::Psi> f_vkspsi(1, perbands_ks * ndim, npwx, p_kv->ngk.data()); + ModuleBase::Memory::record("SDFT::vkspsi", ks_memory_cost); + psi::Psi>* kspsi_all = &f_kspsi; + + size_t sto_memory_cost = perbands_sto * npwx * sizeof(std::complex); + psi::Psi> sfchi(1, perbands_sto, npwx, p_kv->ngk.data()); + ModuleBase::Memory::record("SDFT::sfchi", sto_memory_cost); + psi::Psi> smfchi(1, perbands_sto, npwx, p_kv->ngk.data()); + ModuleBase::Memory::record("SDFT::smfchi", sto_memory_cost); +#ifdef __MPI + psi::Psi> chi_all, hchi_all, psi_all; + if (GlobalV::NSTOGROUP > 1) + { + chi_all.resize(1, allbands_sto, npwx); + hchi_all.resize(1, allbands_sto, npwx); + ModuleBase::Memory::record("SDFT::chi_all", allbands_sto * npwx * sizeof(std::complex)); + ModuleBase::Memory::record("SDFT::hchi_all", allbands_sto * npwx * sizeof(std::complex)); + psi_all.resize(1, allbands_ks, npwx); + ModuleBase::Memory::record("SDFT::kspsi_all", allbands_ks * npwx * sizeof(std::complex)); + for (int ib = 0; ib < allbands_ks; ++ib) + { + for (int ig = 0; ig < npw; ++ig) + { + psi_all(0, ib, ig) + = static_cast>(p_psi[0](this->nbands_ks - allbands_ks + ib, ig)); + } + } + kspsi_all = &psi_all; + f_kspsi.resize(1, 1, 1); + } +#endif + + const int nbatch_psi = npart_sto; + const int bsize_psi = ceil(double(perbands_sto) / nbatch_psi); + psi::Psi> batch_vchi(1, bsize_psi * ndim, npwx, p_kv->ngk.data()); + psi::Psi> batch_vhchi(1, bsize_psi * ndim, npwx, p_kv->ngk.data()); + ModuleBase::Memory::record("SDFT::batchjpsi", 3 * bsize_psi * ndim * npwx * sizeof(std::complex)); + + //------------------- sqrt(f)|psi> sqrt(1-f)|psi> --------------- + if (perbands_ks > 0) + { + for (int ib = 0; ib < perbands_ks; ++ib) + { + for (int ig = 0; ig < npw; ++ig) + { + kspsi(0, ib, ig) = p_psi[0](ib0_ks + ib, ig); + } + double fi = stoiter.stofunc.fd(en[ib]); + expmtmf_fact[ib] = 1 - fi; + expmtf_fact[ib] = fi; + } + // v|\psi> + velop.act(&kspsi, perbands_ks, kspsi.get_pointer(), vkspsi.get_pointer()); + // convert to complex + if (GlobalV::NSTOGROUP == 1) + { + convert_psi(kspsi, f_kspsi); + } + convert_psi(vkspsi, f_vkspsi); + kspsi.resize(1, 1, 1); + vkspsi.resize(1, 1, 1); + } + + che.calcoef_real(&stoiter.stofunc, &Sto_Func::nroot_fd); + che.calfinalvec_real(&stohchi, &Stochastic_hchi::hchi_norm, stopsi->get_pointer(), sfchi.get_pointer(), npw, + npwx, perbands_sto); + + che.calcoef_real(&stoiter.stofunc, &Sto_Func::nroot_mfd); + + che.calfinalvec_real(&stohchi, &Stochastic_hchi::hchi_norm, stopsi->get_pointer(), smfchi.get_pointer(), npw, + npwx, perbands_sto); + + //------------------------ allocate ------------------------ + psi::Psi>& expmtsfchi = sfchi; + psi::Psi>& expmtsmfchi = smfchi; + psi::Psi> exptsfchi = expmtsfchi; + ModuleBase::Memory::record("SDFT::exptsfchi", sto_memory_cost); + psi::Psi> exptsmfchi = expmtsmfchi; + ModuleBase::Memory::record("SDFT::exptsmfchi", sto_memory_cost); + psi::Psi> poly_expmtsfchi, poly_expmtsmfchi; + psi::Psi> poly_exptsfchi, poly_exptsmfchi; + if (nbatch > 1) + { + poly_exptsfchi.resize(cond_nche, perbands_sto, npwx); + ModuleBase::Memory::record("SDFT::poly_exptsfchi", + sizeof(std::complex) * cond_nche * perbands_sto * npwx); + + poly_exptsmfchi.resize(cond_nche, perbands_sto, npwx); + ModuleBase::Memory::record("SDFT::poly_exptsmfchi", + sizeof(std::complex) * cond_nche * perbands_sto * npwx); + + poly_expmtsfchi.resize(cond_nche, perbands_sto, npwx); + ModuleBase::Memory::record("SDFT::poly_expmtsfchi", + sizeof(std::complex) * cond_nche * perbands_sto * npwx); + + poly_expmtsmfchi.resize(cond_nche, perbands_sto, npwx); + ModuleBase::Memory::record("SDFT::poly_expmtsmfchi", + sizeof(std::complex) * cond_nche * perbands_sto * npwx); + } + + const int dim_jmatrix = perbands_ks * allbands_sto + perbands_sto * allbands; + parallel_distribution parajmat(ndim * dim_jmatrix, GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL); + std::vector> j1l(ndim * dim_jmatrix), j2l(ndim * dim_jmatrix); + ModuleBase::Memory::record("SDFT::j1l", sizeof(std::complex) * ndim * dim_jmatrix); + ModuleBase::Memory::record("SDFT::j2l", sizeof(std::complex) * ndim * dim_jmatrix); + std::vector> j1r(ndim * dim_jmatrix), j2r(ndim * dim_jmatrix); + ModuleBase::Memory::record("SDFT::j1r", sizeof(std::complex) * ndim * dim_jmatrix); + ModuleBase::Memory::record("SDFT::j2r", sizeof(std::complex) * ndim * dim_jmatrix); + psi::Psi> tmphchil(1, perbands_sto, npwx, p_kv->ngk.data()); + ModuleBase::Memory::record("SDFT::tmphchil/r", sto_memory_cost * 2); + + //------------------------ t loop -------------------------- + std::cout << "ik=" << ik << ": "; + auto start = std::chrono::high_resolution_clock::now(); + const int print_step = ceil(20.0 / nbatch) * nbatch; + for (int it = 1; it < nt; ++it) + { + // evaluate time cost + if (it - 1 == print_step) + { + auto end = std::chrono::high_resolution_clock::now(); + std::chrono::duration duration = end - start; + double timeTaken = duration.count(); + std::cout << "(Time left " << timeTaken * (double(nt - 1) / print_step * (nk - ik) - 1) << " s) " + << std::endl; + std::cout << "nt: " << std::endl; + } + if ((it - 1) % print_step == 0 && it > 1) + { + std::cout << std::setw(8) << it - 1; + if ((it - 1) % (print_step * 10) == 0) + { + std::cout << std::endl; + } + } + + // time evolution exp(-iHt)|\psi_ks> + // KS + ModuleBase::timer::tick("Sto_EleCond", "evolution"); + for (int ib = 0; ib < perbands_ks; ++ib) + { + double eigen = en[ib]; + const std::complex expmfactor = exp(ModuleBase::NEG_IMAG_UNIT * eigen * dt); + expmtf_fact[ib] *= expmfactor; + expmtmf_fact[ib] *= expmfactor; + } + // Sto + if (nbatch == 1) + { + chemt.calfinalvec_complex(&stohchi, &Stochastic_hchi::hchi_norm, expmtsfchi.get_pointer(), + expmtsfchi.get_pointer(), npw, npwx, perbands_sto); + chemt.calfinalvec_complex(&stohchi, &Stochastic_hchi::hchi_norm, expmtsmfchi.get_pointer(), + expmtsmfchi.get_pointer(), npw, npwx, perbands_sto); + chet.calfinalvec_complex(&stohchi, &Stochastic_hchi::hchi_norm, exptsfchi.get_pointer(), + exptsfchi.get_pointer(), npw, npwx, perbands_sto); + chet.calfinalvec_complex(&stohchi, &Stochastic_hchi::hchi_norm, exptsmfchi.get_pointer(), + exptsmfchi.get_pointer(), npw, npwx, perbands_sto); + } + else + { + std::complex* tmppolyexpmtsfchi = poly_expmtsfchi.get_pointer(); + std::complex* tmppolyexpmtsmfchi = poly_expmtsmfchi.get_pointer(); + std::complex* tmppolyexptsfchi = poly_exptsfchi.get_pointer(); + std::complex* tmppolyexptsmfchi = poly_exptsmfchi.get_pointer(); + std::complex* stoexpmtsfchi = expmtsfchi.get_pointer(); + std::complex* stoexpmtsmfchi = expmtsmfchi.get_pointer(); + std::complex* stoexptsfchi = exptsfchi.get_pointer(); + std::complex* stoexptsmfchi = exptsmfchi.get_pointer(); + if ((it - 1) % nbatch == 0) + { + chet.calpolyvec_complex(&stohchi, &Stochastic_hchi::hchi_norm, stoexptsfchi, tmppolyexptsfchi, npw, + npwx, perbands_sto); + chet.calpolyvec_complex(&stohchi, &Stochastic_hchi::hchi_norm, stoexptsmfchi, tmppolyexptsmfchi, + npw, npwx, perbands_sto); + chemt.calpolyvec_complex(&stohchi, &Stochastic_hchi::hchi_norm, stoexpmtsfchi, tmppolyexpmtsfchi, + npw, npwx, perbands_sto); + chemt.calpolyvec_complex(&stohchi, &Stochastic_hchi::hchi_norm, stoexpmtsmfchi, tmppolyexpmtsmfchi, + npw, npwx, perbands_sto); + } + + std::complex* tmpcoef = batchcoef.data() + (it - 1) % nbatch * cond_nche; + std::complex* tmpmcoef = batchmcoef.data() + (it - 1) % nbatch * cond_nche; + const char transa = 'N'; + const int LDA = perbands_sto * npwx; + const int M = perbands_sto * npwx; + const int N = cond_nche; + const int inc = 1; + zgemv_(&transa, &M, &N, &ModuleBase::ONE, tmppolyexptsfchi, &LDA, tmpcoef, &inc, &ModuleBase::ZERO, + stoexptsfchi, &inc); + zgemv_(&transa, &M, &N, &ModuleBase::ONE, tmppolyexptsmfchi, &LDA, tmpcoef, &inc, &ModuleBase::ZERO, + stoexptsmfchi, &inc); + zgemv_(&transa, &M, &N, &ModuleBase::ONE, tmppolyexpmtsfchi, &LDA, tmpmcoef, &inc, &ModuleBase::ZERO, + stoexpmtsfchi, &inc); + zgemv_(&transa, &M, &N, &ModuleBase::ONE, tmppolyexpmtsmfchi, &LDA, tmpmcoef, &inc, &ModuleBase::ZERO, + stoexpmtsmfchi, &inc); + } + ModuleBase::timer::tick("Sto_EleCond", "evolution"); + + // calculate i<\psi|sqrt(f) exp(-iHt/2)*J*exp(iHt/2) sqrt(1-f)|\psi>^+ + // = i<\psi|sqrt(1-f) exp(-iHt/2)*J*exp(iHt/2) sqrt(f)|\psi> + cal_jmatrix(*kspsi_all, f_vkspsi, en.data(), en_all, nullptr, nullptr, exptsmfchi, exptsfchi, tmphchil, + batch_vchi, batch_vhchi, +#ifdef __MPI + chi_all, hchi_all, (void*)&ks_fact, (void*)&sto_npwx, +#endif + bsize_psi, j1l, j2l, velop, ik, ModuleBase::IMAG_UNIT, bandsinfo); + + // calculate <\psi|sqrt(1-f) exp(iHt/2)*J*exp(-iHt/2) sqrt(f)|\psi> + cal_jmatrix(*kspsi_all, f_vkspsi, en.data(), en_all, expmtmf_fact.data(), expmtf_fact.data(), expmtsmfchi, + expmtsfchi, tmphchil, batch_vchi, batch_vhchi, +#ifdef __MPI + chi_all, hchi_all, (void*)&ks_fact, (void*)&sto_npwx, +#endif + bsize_psi, j1r, j2r, velop, ik, ModuleBase::ONE, bandsinfo); + + // prepare for parallel + int num_per = parajmat.num_per; + int st_per = parajmat.start; + // Re(i) + // Im(l_ij*r_ji) = Re(-il_ij * r_ji) = Re( ((il)^+_ji)^* * r_ji)=Re(((il)^+_i)^* * r^+_i) + // ddot_real = real(A_i^* * B_i) + ModuleBase::timer::tick("Sto_EleCond", "ddot_real"); + ct11[it] += static_cast( + ModuleBase::GlobalFunc::ddot_real(num_per, j1l.data() + st_per, j1r.data() + st_per, false) + * p_kv->wk[ik] / 2.0); + double tmp12 = static_cast( + ModuleBase::GlobalFunc::ddot_real(num_per, j1l.data() + st_per, j2r.data() + st_per, false)); + + double tmp21 = static_cast( + ModuleBase::GlobalFunc::ddot_real(num_per, j2l.data() + st_per, j1r.data() + st_per, false)); + + ct12[it] -= 0.5 * (tmp12 + tmp21) * p_kv->wk[ik] / 2.0; + + ct22[it] += static_cast( + ModuleBase::GlobalFunc::ddot_real(num_per, j2l.data() + st_per, j2r.data() + st_per, false) + * p_kv->wk[ik] / 2.0); + + ModuleBase::timer::tick("Sto_EleCond", "ddot_real"); + } + std::cout << std::endl; + } // ik loop + ModuleBase::timer::tick("Sto_EleCond", "kloop"); +#ifdef __MPI + MPI_Allreduce(MPI_IN_PLACE, ct11.data(), nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, ct12.data(), nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, ct22.data(), nt, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#endif + + //------------------------------------------------------------------ + // Output + //------------------------------------------------------------------ + if (GlobalV::MY_RANK == 0) + { + calcondw(nt, dt, smear_type, fwhmin, wcut, dw_in, ct11.data(), ct12.data(), ct22.data()); + } + ModuleBase::timer::tick("Sto_EleCond", "sKG"); +} + +namespace GlobalTemp +{ +const ModuleBase::matrix* veff; +} diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_elecond.h b/source/module_hamilt_pw/hamilt_stodft/sto_elecond.h new file mode 100644 index 0000000000..8a2bf518cb --- /dev/null +++ b/source/module_hamilt_pw/hamilt_stodft/sto_elecond.h @@ -0,0 +1,72 @@ +#ifndef STOELECOND_H +#define STOELECOND_H + +#include "module_hamilt_general/hamilt.h" +#include "module_hamilt_pw/hamilt_pwdft/elecond.h" +#include "module_hamilt_pw/hamilt_stodft/sto_wf.h" +#include "module_hsolver/hsolver_pw_sdft.h" + +class Sto_EleCond : protected EleCond +{ + public: + Sto_EleCond(UnitCell* p_ucell_in, K_Vectors* p_kv_in, elecstate::ElecState* p_elec_in, + ModulePW::PW_Basis_K* p_wfcpw_in, psi::Psi>* p_psi_in, + pseudopot_cell_vnl* p_ppcell_in, hamilt::Hamilt>* p_hamilt_in, + hsolver::HSolverPW_SDFT* p_hsol_in, Stochastic_WF* p_stowf_in); + ~Sto_EleCond(){}; + /** + * @brief Set the N order of Chebyshev expansion for conductivities + * It will change class member : fd_nche, cond_nche + * + * @param dt t step + * @param nbatch number of t batch + * @param cond_thr threshold of errors for conductivities + * @param fd_nche N order of Chebyshev for Fermi-Dirac function + * @param try_emin trial Emin + * @param try_emax trial Emax + * + */ + void decide_nche(const double dt, int& nbatch, const double cond_thr, const int& fd_nche, double try_emin, + double try_emax); + /** + * @brief calculate Stochastic Kubo-Greenwood + * + * @param fwhmin FWHM + * @param smear_type 1: Gaussian, 2: Lorentzian + * @param wcut cutoff omega + * @param dw_in omega step + * @param dt_in t step + * @param nonlocal whether to include the nonlocal potential corrections for velocity operator + * @param nbatch t step batch + * @param npart_sto number stochastic wavefunctions parts to evalution simultaneously + */ + void sKG(const int& smear_type, const double& fwhmin, const double& wcut, const double& dw_in, const double& dt_in, + const bool& nonlocal, const int& nbatch, const int& npart_sto); + + protected: + int nbands_ks = 0; ///< number of KS bands + int nbands_sto = 0; ///< number of stochastic bands + int cond_nche = 0; ///< number of Chebyshev orders for conductivities + int fd_nche = 0; ///< number of Chebyshev orders for Fermi-Dirac function + hamilt::Hamilt>* p_hamilt; ///< pointer to the Hamiltonian + hsolver::HSolverPW_SDFT* p_hsol = nullptr; ///< pointer to the Hamiltonian solver + Stochastic_WF* p_stowf = nullptr; ///< pointer to the stochastic wavefunctions + + protected: + /** + * @brief calculate Jmatrix + * + */ + void cal_jmatrix(const psi::Psi>& kspsi_all, const psi::Psi>& vkspsi, + const double* en, const double* en_all, std::complex* leftfact, + std::complex* rightfact, const psi::Psi>& leftchi, + psi::Psi>& rightchi, psi::Psi>& left_hchi, + psi::Psi>& batch_vchi, psi::Psi>& batch_vhchi, +#ifdef __MPI + psi::Psi>& chi_all, psi::Psi>& hchi_all, + void* gatherinfo_ks, void* gatherinfo_sto, +#endif + const int& bsize_psi, std::vector>& j1, std::vector>& j2, + hamilt::Velocity& velop, const int& ik, const std::complex& factor, const int bandinfo[6]); +}; +#endif // ELECOND_H \ No newline at end of file diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_tool.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_tool.cpp new file mode 100644 index 0000000000..091092e953 --- /dev/null +++ b/source/module_hamilt_pw/hamilt_stodft/sto_tool.cpp @@ -0,0 +1,125 @@ +#include "sto_tool.h" + +#include "module_base/timer.h" +#ifdef __MPI +#include "mpi.h" +#endif +#include + +void check_che(const int& nche_in, const double& try_emin, const double& try_emax, const int& nbands_sto, + K_Vectors* p_kv, Stochastic_WF* p_stowf, hamilt::Hamilt>* p_hamilt, + hsolver::HSolverPW_SDFT* p_hsol) +{ + //------------------------------ + // Convergence test + //------------------------------ + bool change = false; + const int nk = p_kv->get_nks(); + ModuleBase::Chebyshev chetest(nche_in); + Stochastic_Iter& stoiter = p_hsol->stoiter; + Stochastic_hchi& stohchi = stoiter.stohchi; + int ntest0 = 5; + stohchi.Emax = try_emax; + stohchi.Emin = try_emin; + // if (GlobalV::NBANDS > 0) + // { + // double tmpemin = 1e10; + // for (int ik = 0; ik < nk; ++ik) + // { + // tmpemin = std::min(tmpemin, this->pelec->ekb(ik, GlobalV::NBANDS - 1)); + // } + // stohchi.Emin = tmpemin; + // } + // else + // { + // stohchi.Emin = 0; + // } + for (int ik = 0; ik < nk; ++ik) + { + p_hamilt->updateHk(ik); + stoiter.stohchi.current_ik = ik; + const int npw = p_kv->ngk[ik]; + std::complex* pchi = nullptr; + std::vector> randchi; + int ntest = std::min(ntest0, p_stowf->nchip[ik]); + for (int i = 0; i < ntest; ++i) + { + if (nbands_sto == 0) + { + randchi.resize(npw); + pchi = &randchi[0]; + for (int ig = 0; ig < npw; ++ig) + { + double rr = std::rand() / double(RAND_MAX); + double arg = std::rand() / double(RAND_MAX); + pchi[ig] = std::complex(rr * cos(arg), rr * sin(arg)); + } + } + else if (GlobalV::NBANDS > 0) + { + pchi = &p_stowf->chiortho[0](ik, i, 0); + } + else + { + pchi = &p_stowf->chi0[0](ik, i, 0); + } + while (1) + { + bool converge; + converge = chetest.checkconverge(&stohchi, &Stochastic_hchi::hchi_norm, pchi, npw, stohchi.Emax, + stohchi.Emin, 2.0); + + if (!converge) + { + change = true; + } + else + { + break; + } + } + } + + if (ik == nk - 1) + { +#ifdef __MPI + MPI_Allreduce(MPI_IN_PLACE, &stohchi.Emax, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &stohchi.Emin, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); +#endif + stoiter.stofunc.Emax = stohchi.Emax; + stoiter.stofunc.Emin = stohchi.Emin; + GlobalV::ofs_running << "New Emax " << stohchi.Emax << " Ry; new Emin " << stohchi.Emin << " Ry" + << std::endl; + change = false; + } + } +} + +void convert_psi(const psi::Psi>& psi_in, psi::Psi>& psi_out) +{ + psi_in.fix_k(0); + psi_out.fix_k(0); + for (int i = 0; i < psi_in.size(); ++i) + { + psi_out.get_pointer()[i] = static_cast>(psi_in.get_pointer()[i]); + } + return; +} + +psi::Psi>* gatherchi(psi::Psi>& chi, psi::Psi>& chi_all, + const int& npwx, int* nrecv_sto, int* displs_sto, const int perbands_sto) +{ + psi::Psi>* p_chi; + p_chi = χ +#ifdef __MPI + if (GlobalV::NSTOGROUP > 1) + { + p_chi = &chi_all; + ModuleBase::timer::tick("sKG", "bands_gather"); + MPI_Allgatherv(chi.get_pointer(), perbands_sto * npwx, MPI_COMPLEX, chi_all.get_pointer(), nrecv_sto, + displs_sto, MPI_COMPLEX, PARAPW_WORLD); + ModuleBase::timer::tick("sKG", "bands_gather"); + } +#endif + return p_chi; +} \ No newline at end of file diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_tool.h b/source/module_hamilt_pw/hamilt_stodft/sto_tool.h new file mode 100644 index 0000000000..90fa78a818 --- /dev/null +++ b/source/module_hamilt_pw/hamilt_stodft/sto_tool.h @@ -0,0 +1,108 @@ +#include "module_cell/klist.h" +#include "module_hamilt_general/hamilt.h" +#include "module_hamilt_pw/hamilt_stodft/sto_wf.h" +#include "module_hsolver/hsolver_pw_sdft.h" +#include "module_psi/psi.h" +/** + * @brief Check if Emin and Emax are converged + * + * @param nche_in N order of Chebyshev expansion + * @param try_emin trial Emin + * @param try_emax trial Emax + * @param nbands_sto number of stochastic bands + */ +void check_che(const int& nche_in, const double& try_emin, const double& try_emax, const int& nbands_sto, + K_Vectors* p_kv, Stochastic_WF* p_stowf, hamilt::Hamilt>* p_hamilt, + hsolver::HSolverPW_SDFT* p_hsol_in); + +#ifndef PARALLEL_DISTRIBUTION +#define PARALLEL_DISTRIBUTION +/** + * @brief structure to distribute calculation among processors + * + * @param start start index of this processor + * @param num_per number of elements for this processor + * + */ +struct parallel_distribution +{ + parallel_distribution(const int& num_all, const int& np, const int myrank) + { + int num_per = num_all / np; + int st_per = num_per * myrank; + int re = num_all % np; + if (myrank < re) + { + ++num_per; + st_per += myrank; + } + else + { + st_per += re; + } + this->start = st_per; + this->num_per = num_per; + } + int start; + int num_per; +}; +#endif + +#ifdef __MPI +#ifndef INFO_GATHERV +#define INFO_GATHERV +/** + * @brief gather information from all processors + * + */ +struct info_gatherv +{ + info_gatherv(const int& ngroup_per, const int& np, const int& num_in_group, MPI_Comm comm_world) + { + nrecv = new int[np]; + displs = new int[np]; + MPI_Allgather(&ngroup_per, 1, MPI_INT, nrecv, 1, MPI_INT, comm_world); + displs[0] = 0; + for (int i = 1; i < np; ++i) + { + displs[i] = displs[i - 1] + nrecv[i - 1]; + } + for (int i = 0; i < np; ++i) + { + nrecv[i] *= num_in_group; + displs[i] *= num_in_group; + } + } + ~info_gatherv() + { + delete[] nrecv; + delete[] displs; + } + int* nrecv = nullptr; + int* displs = nullptr; +}; +#endif +#endif + +/** + * @brief convert psi from double to float + * + * @param psi_in input psi of double + * @param psi_out output psi of float + */ +void convert_psi(const psi::Psi>& psi_in, psi::Psi>& psi_out); + +/** + * @brief gather chi from all processors + * + * @param chi stochasitc wave function of this processor + * @param chi_all gathered stochastic wave function + * @param npwx maximum number of plane waves on all processors + * @param nrecv_sto number of stochastic orbitals on each processor + * @param displs_sto displacement of stochastic orbitals on each processor + * @param perbands_sto number of stochastic bands of this processor + * @return psi::Psi> pointer to gathered stochastic wave function + * + */ +psi::Psi>* gatherchi(psi::Psi>& chi, psi::Psi>& chi_all, + const int& npwx, int* nrecv_sto, int* displs_sto, const int perbands_sto); diff --git a/source/module_hamilt_pw/hamilt_stodft/test/CMakeLists.txt b/source/module_hamilt_pw/hamilt_stodft/test/CMakeLists.txt new file mode 100644 index 0000000000..c9bb089dae --- /dev/null +++ b/source/module_hamilt_pw/hamilt_stodft/test/CMakeLists.txt @@ -0,0 +1,7 @@ +remove_definitions(-D__MPI) + +AddTest( + TARGET Sto_Tool_UTs + LIBS ${math_libs} psi base device + SOURCES ../sto_tool.cpp test_sto_tool.cpp +) \ No newline at end of file diff --git a/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp b/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp new file mode 100644 index 0000000000..e57b74e4da --- /dev/null +++ b/source/module_hamilt_pw/hamilt_stodft/test/test_sto_tool.cpp @@ -0,0 +1,75 @@ +#include "../sto_tool.h" +#include "mpi.h" + +#include + +/************************************************ + * unit test of sto_tool.cpp + ***********************************************/ + +void Stochastic_hchi::hchi_norm(std::complex* chig, std::complex* hchig, const int m) +{ + return; +} + +/** + * - Tested Functions: + * - struct parallel_distribution + * - void convert_psi(psi_in, psi_out) + * - psi::Psi>* gatherchi(chi, chi_all, npwx, nrecv_sto, displs_sto, perbands_sto) + */ +class TestStoTool : public ::testing::Test +{ +}; + +TEST_F(TestStoTool, parallel_distribution) +{ + int num_all = 10; + int np = 4; + int myrank = 0; + parallel_distribution pd(num_all, np, myrank); + EXPECT_EQ(pd.start, 0); + EXPECT_EQ(pd.num_per, 3); + + myrank = 1; + parallel_distribution pd1(num_all, np, myrank); + EXPECT_EQ(pd1.start, 3); + EXPECT_EQ(pd1.num_per, 3); + + myrank = 2; + parallel_distribution pd2(num_all, np, myrank); + EXPECT_EQ(pd2.start, 6); + EXPECT_EQ(pd2.num_per, 2); + + myrank = 3; + parallel_distribution pd3(num_all, np, myrank); + EXPECT_EQ(pd3.start, 8); + EXPECT_EQ(pd3.num_per, 2); +} + +TEST_F(TestStoTool, convert_psi) +{ + psi::Psi> psi_in(1, 1, 10); + psi::Psi> psi_out(1, 1, 10); + for (int i = 0; i < 10; ++i) + { + psi_in.get_pointer()[i] = std::complex(i, i); + } + convert_psi(psi_in, psi_out); + for (int i = 0; i < 10; ++i) + { + EXPECT_EQ(psi_out.get_pointer()[i], std::complex(i, i)); + } +} + +TEST_F(TestStoTool, gatherchi) +{ + psi::Psi> chi(1, 1, 10); + psi::Psi> chi_all(1, 1, 10); + int npwx = 10; + int nrecv_sto[4] = {1, 2, 3, 4}; + int displs_sto[4] = {0, 1, 3, 6}; + int perbands_sto = 1; + psi::Psi>* p_chi = gatherchi(chi, chi_all, npwx, nrecv_sto, displs_sto, perbands_sto); + EXPECT_EQ(p_chi, &chi); +} diff --git a/tests/integrate/186_PW_SKG_10D10S/refOnsager.txt b/tests/integrate/186_PW_SKG_10D10S/refOnsager.txt index 5f42b90b71..f4cd8c5530 100644 --- a/tests/integrate/186_PW_SKG_10D10S/refOnsager.txt +++ b/tests/integrate/186_PW_SKG_10D10S/refOnsager.txt @@ -30,7 +30,7 @@ 0.57 530099 798.393 -1.07346e+07 2.93011e+08 0.59 529977 798.11 -1.0733e+07 2.9297e+08 0.61 529852 797.817 -1.07314e+07 2.92929e+08 - 0.63 529722 797.514 -1.07297e+07 2.92885e+08 + 0.63 529722 797.514 -1.07297e+07 2.92886e+08 0.65 529588 797.202 -1.0728e+07 2.92841e+08 0.67 529450 796.88 -1.07262e+07 2.92795e+08 0.69 529308 796.548 -1.07244e+07 2.92748e+08 @@ -191,7 +191,7 @@ 3.79 458881 645.563 -9.73586e+06 2.67717e+08 3.81 458136 644.092 -9.72458e+06 2.67434e+08 3.83 457387 642.616 -9.71324e+06 2.6715e+08 - 3.85 456635 641.137 -9.70182e+06 2.66865e+08 + 3.85 456635 641.136 -9.70182e+06 2.66865e+08 3.87 455880 639.653 -9.69034e+06 2.66577e+08 3.89 455121 638.165 -9.67879e+06 2.66288e+08 3.91 454360 636.673 -9.66718e+06 2.65997e+08 @@ -223,7 +223,7 @@ 4.43 433437 596.62 -9.3416e+06 2.57853e+08 4.45 432591 595.036 -9.32817e+06 2.57517e+08 4.47 431742 593.45 -9.31467e+06 2.57179e+08 - 4.49 430890 591.861 -9.3011e+06 2.5684e+08 + 4.49 430890 591.86 -9.3011e+06 2.5684e+08 4.51 430035 590.269 -9.28747e+06 2.56499e+08 4.53 429178 588.674 -9.27377e+06 2.56157e+08 4.55 428317 587.077 -9.26001e+06 2.55812e+08 @@ -455,10 +455,10 @@ 9.07 197918 231.952 -4.89512e+06 1.43045e+08 9.09 196946 230.68 -4.87402e+06 1.42475e+08 9.11 195977 229.412 -4.85295e+06 1.41906e+08 - 9.13 195009 228.15 -4.83191e+06 1.41338e+08 + 9.13 195009 228.149 -4.83191e+06 1.41338e+08 9.15 194044 226.891 -4.8109e+06 1.4077e+08 9.17 193082 225.638 -4.78993e+06 1.40202e+08 - 9.19 192122 224.389 -4.76898e+06 1.39636e+08 + 9.19 192122 224.388 -4.76898e+06 1.39636e+08 9.21 191164 223.144 -4.74806e+06 1.39069e+08 9.23 190209 221.904 -4.72718e+06 1.38504e+08 9.25 189257 220.67 -4.70633e+06 1.37939e+08 @@ -536,7 +536,7 @@ 10.69 128494 145.035 -3.32862e+06 9.99669e+07 10.71 127773 144.179 -3.31164e+06 9.94905e+07 10.73 127055 143.328 -3.29474e+06 9.90159e+07 - 10.75 126340 142.483 -3.2779e+06 9.8543e+07 + 10.75 126340 142.482 -3.2779e+06 9.8543e+07 10.77 125630 141.643 -3.26114e+06 9.80719e+07 10.79 124923 140.808 -3.24444e+06 9.76025e+07 10.81 124219 139.98 -3.22782e+06 9.71348e+07 @@ -600,9 +600,9 @@ 11.97 89773.5 101.284 -2.39125e+06 7.32896e+07 11.99 89288.5 100.776 -2.3791e+06 7.29381e+07 12.01 88807.2 100.272 -2.36703e+06 7.25887e+07 - 12.03 88329.4 99.7741 -2.35503e+06 7.22413e+07 + 12.03 88329.4 99.774 -2.35503e+06 7.22413e+07 12.05 87855.3 99.281 -2.34311e+06 7.18961e+07 - 12.07 87384.7 98.7932 -2.33127e+06 7.15528e+07 + 12.07 87384.7 98.7931 -2.33127e+06 7.15528e+07 12.09 86917.7 98.3104 -2.3195e+06 7.12117e+07 12.11 86454.3 97.8328 -2.30781e+06 7.08727e+07 12.13 85994.4 97.3603 -2.2962e+06 7.05357e+07 @@ -614,7 +614,7 @@ 12.25 83309.4 94.6319 -2.22814e+06 6.85573e+07 12.27 82874.2 94.1948 -2.21707e+06 6.82348e+07 12.29 82442.5 93.7626 -2.20607e+06 6.79143e+07 - 12.31 82014.3 93.3355 -2.19515e+06 6.75959e+07 + 12.31 82014.3 93.3354 -2.19515e+06 6.75959e+07 12.33 81589.5 92.9132 -2.1843e+06 6.72796e+07 12.35 81168.1 92.496 -2.17353e+06 6.69654e+07 12.37 80750.2 92.0836 -2.16283e+06 6.66532e+07 @@ -626,13 +626,13 @@ 12.49 78314.6 89.7117 -2.10024e+06 6.48229e+07 12.51 77920.5 89.3333 -2.09007e+06 6.4525e+07 12.53 77529.8 88.9596 -2.07998e+06 6.42291e+07 - 12.55 77142.4 88.5907 -2.06995e+06 6.39353e+07 - 12.57 76758.3 88.2265 -2.06001e+06 6.36434e+07 + 12.55 77142.4 88.5906 -2.06995e+06 6.39353e+07 + 12.57 76758.3 88.2264 -2.06001e+06 6.36434e+07 12.59 76377.5 87.867 -2.05013e+06 6.33536e+07 - 12.61 76000 87.5122 -2.04033e+06 6.30658e+07 + 12.61 76000 87.5121 -2.04033e+06 6.30658e+07 12.63 75625.9 87.162 -2.0306e+06 6.278e+07 12.65 75255 86.8165 -2.02095e+06 6.24962e+07 - 12.67 74887.3 86.4757 -2.01137e+06 6.22144e+07 + 12.67 74887.3 86.4756 -2.01137e+06 6.22144e+07 12.69 74523 86.1394 -2.00186e+06 6.19347e+07 12.71 74161.8 85.8077 -1.99242e+06 6.16568e+07 12.73 73803.9 85.4806 -1.98306e+06 6.1381e+07 @@ -644,7 +644,7 @@ 12.85 71723.3 83.6122 -1.92837e+06 5.97675e+07 12.87 71387.6 83.3163 -1.91951e+06 5.95054e+07 12.89 71054.9 83.0248 -1.91071e+06 5.92453e+07 - 12.91 70725.4 82.7377 -1.90198e+06 5.89871e+07 + 12.91 70725.4 82.7376 -1.90198e+06 5.89871e+07 12.93 70398.9 82.4548 -1.89333e+06 5.87308e+07 12.95 70075.5 82.1763 -1.88474e+06 5.84764e+07 12.97 69755.2 81.902 -1.87622e+06 5.8224e+07 @@ -655,19 +655,19 @@ 13.07 68198.8 80.5938 -1.83466e+06 5.69903e+07 13.09 67896.5 80.3447 -1.82655e+06 5.67492e+07 13.11 67597.1 80.0996 -1.81851e+06 5.65099e+07 - 13.13 67300.6 79.8586 -1.81053e+06 5.62726e+07 + 13.13 67300.6 79.8585 -1.81053e+06 5.62726e+07 13.15 67007 79.6216 -1.80263e+06 5.6037e+07 13.17 66716.4 79.3886 -1.79478e+06 5.58034e+07 13.19 66428.6 79.1596 -1.78701e+06 5.55715e+07 13.21 66143.7 78.9345 -1.7793e+06 5.53415e+07 13.23 65861.6 78.7134 -1.77165e+06 5.51133e+07 13.25 65582.4 78.4962 -1.76407e+06 5.48869e+07 - 13.27 65306 78.2829 -1.75655e+06 5.46624e+07 + 13.27 65306 78.2828 -1.75655e+06 5.46624e+07 13.29 65032.4 78.0733 -1.7491e+06 5.44396e+07 13.31 64761.5 77.8676 -1.74171e+06 5.42186e+07 13.33 64493.4 77.6657 -1.73439e+06 5.39994e+07 13.35 64228.1 77.4676 -1.72713e+06 5.37819e+07 - 13.37 63965.5 77.2732 -1.71993e+06 5.35663e+07 + 13.37 63965.5 77.2731 -1.71993e+06 5.35663e+07 13.39 63705.7 77.0824 -1.71279e+06 5.33523e+07 13.41 63448.5 76.8954 -1.70572e+06 5.31402e+07 13.43 63194 76.7119 -1.69871e+06 5.29297e+07 @@ -793,17 +793,17 @@ 15.83 46921.9 71.9123 -1.19665e+06 3.73308e+07 15.85 46863.8 71.9424 -1.1944e+06 3.72564e+07 15.87 46806.4 71.9725 -1.19216e+06 3.71825e+07 - 15.89 46749.6 72.0028 -1.18994e+06 3.71091e+07 + 15.89 46749.6 72.0028 -1.18994e+06 3.71092e+07 15.91 46693.5 72.0331 -1.18774e+06 3.70363e+07 15.93 46638 72.0634 -1.18556e+06 3.6964e+07 15.95 46583.2 72.0938 -1.18339e+06 3.68923e+07 - 15.97 46529 72.1241 -1.18124e+06 3.6821e+07 + 15.97 46529 72.1242 -1.18124e+06 3.6821e+07 15.99 46475.3 72.1545 -1.17911e+06 3.67502e+07 16.01 46422.3 72.1848 -1.177e+06 3.66799e+07 16.03 46369.8 72.2151 -1.1749e+06 3.66101e+07 16.05 46317.9 72.2453 -1.17281e+06 3.65408e+07 16.07 46266.5 72.2754 -1.17075e+06 3.6472e+07 - 16.09 46215.7 72.3053 -1.1687e+06 3.64036e+07 + 16.09 46215.7 72.3054 -1.1687e+06 3.64036e+07 16.11 46165.5 72.3352 -1.16666e+06 3.63357e+07 16.13 46115.7 72.3649 -1.16464e+06 3.62682e+07 16.15 46066.5 72.3945 -1.16264e+06 3.62012e+07 @@ -854,7 +854,7 @@ 17.05 44226.1 73.3138 -1.08481e+06 3.35542e+07 17.07 44190.6 73.3207 -1.08329e+06 3.35019e+07 17.09 44155.2 73.3268 -1.08179e+06 3.34498e+07 - 17.11 44119.9 73.3321 -1.08029e+06 3.33979e+07 + 17.11 44119.9 73.3322 -1.08029e+06 3.33979e+07 17.13 44084.8 73.3368 -1.07879e+06 3.33463e+07 17.15 44049.7 73.3406 -1.0773e+06 3.32949e+07 17.17 44014.8 73.3437 -1.07582e+06 3.32436e+07 @@ -863,8 +863,8 @@ 17.23 43910.5 73.3482 -1.07142e+06 3.30912e+07 17.25 43875.9 73.3481 -1.06997e+06 3.30408e+07 17.27 43841.4 73.3472 -1.06852e+06 3.29906e+07 - 17.29 43807 73.3455 -1.06707e+06 3.29406e+07 - 17.31 43772.6 73.343 -1.06564e+06 3.28908e+07 + 17.29 43807 73.3456 -1.06707e+06 3.29406e+07 + 17.31 43772.6 73.3431 -1.06564e+06 3.28908e+07 17.33 43738.3 73.3397 -1.06421e+06 3.28412e+07 17.35 43704.1 73.3356 -1.06278e+06 3.27917e+07 17.37 43669.9 73.3306 -1.06136e+06 3.27425e+07 @@ -888,12 +888,12 @@ 17.73 43059.2 73.0948 -1.03674e+06 3.18863e+07 17.75 43025.3 73.0734 -1.03542e+06 3.18404e+07 17.77 42991.4 73.051 -1.03411e+06 3.17945e+07 - 17.79 42957.5 73.0278 -1.0328e+06 3.17489e+07 + 17.79 42957.5 73.0279 -1.0328e+06 3.17489e+07 17.81 42923.5 73.0038 -1.03149e+06 3.17034e+07 17.83 42889.5 72.9788 -1.03019e+06 3.1658e+07 17.85 42855.5 72.9529 -1.02889e+06 3.16128e+07 17.87 42821.5 72.9261 -1.02759e+06 3.15677e+07 - 17.89 42787.4 72.8985 -1.0263e+06 3.15228e+07 + 17.89 42787.4 72.8985 -1.0263e+06 3.15229e+07 17.91 42753.4 72.8699 -1.02502e+06 3.14781e+07 17.93 42719.2 72.8404 -1.02374e+06 3.14335e+07 17.95 42685.1 72.8101 -1.02246e+06 3.13891e+07 @@ -911,7 +911,7 @@ 18.19 42271.3 72.3762 -1.00746e+06 3.08675e+07 18.21 42236.5 72.3343 -1.00624e+06 3.0825e+07 18.23 42201.6 72.2914 -1.00502e+06 3.07826e+07 - 18.25 42166.7 72.2477 -1.0038e+06 3.07404e+07 + 18.25 42166.7 72.2478 -1.0038e+06 3.07404e+07 18.27 42131.7 72.2032 -1.00259e+06 3.06983e+07 18.29 42096.6 72.1577 -1.00138e+06 3.06564e+07 18.31 42061.4 72.1114 -1.00018e+06 3.06147e+07 @@ -926,16 +926,16 @@ 18.49 41742 71.6557 -989522 3.02454e+07 18.51 41706.2 71.6008 -988357 3.02051e+07 18.53 41670.3 71.545 -987196 3.0165e+07 - 18.55 41634.3 71.4884 -986039 3.0125e+07 + 18.55 41634.3 71.4885 -986039 3.0125e+07 18.57 41598.2 71.431 -984886 3.00852e+07 18.59 41562.1 71.3728 -983736 3.00455e+07 18.61 41525.8 71.3137 -982591 3.00059e+07 18.63 41489.6 71.2539 -981449 2.99666e+07 - 18.65 41453.2 71.1931 -980311 2.99273e+07 + 18.65 41453.2 71.1932 -980311 2.99273e+07 18.67 41416.8 71.1316 -979177 2.98882e+07 18.69 41380.2 71.0693 -978047 2.98493e+07 18.71 41343.7 71.0062 -976921 2.98106e+07 - 18.73 41307 70.9422 -975799 2.97719e+07 + 18.73 41307 70.9423 -975799 2.97719e+07 18.75 41270.3 70.8775 -974681 2.97335e+07 18.77 41233.4 70.812 -973567 2.96952e+07 18.79 41196.6 70.7457 -972456 2.9657e+07 @@ -943,9 +943,9 @@ 18.83 41122.6 70.6108 -970247 2.95812e+07 18.85 41085.4 70.5421 -969149 2.95435e+07 18.87 41048.3 70.4727 -968054 2.9506e+07 - 18.89 41011 70.4025 -966963 2.94686e+07 + 18.89 41011 70.4026 -966963 2.94686e+07 18.91 40973.7 70.3316 -965877 2.94314e+07 - 18.93 40936.3 70.2599 -964794 2.93944e+07 + 18.93 40936.3 70.26 -964794 2.93944e+07 18.95 40898.8 70.1875 -963715 2.93575e+07 18.97 40861.2 70.1144 -962640 2.93207e+07 18.99 40823.6 70.0405 -961570 2.92842e+07 @@ -959,9 +959,9 @@ 19.15 40520 69.4234 -953150 2.89975e+07 19.17 40481.8 69.3431 -952116 2.89625e+07 19.19 40443.5 69.2621 -951086 2.89275e+07 - 19.21 40405.1 69.1804 -950060 2.88928e+07 + 19.21 40405.1 69.1805 -950060 2.88928e+07 19.23 40366.6 69.0981 -949039 2.88582e+07 - 19.25 40328.1 69.0151 -948021 2.88238e+07 + 19.25 40328.1 69.0152 -948021 2.88238e+07 19.27 40289.5 68.9315 -947008 2.87896e+07 19.29 40250.9 68.8472 -945999 2.87555e+07 19.31 40212.1 68.7623 -944994 2.87216e+07 @@ -975,7 +975,7 @@ 19.47 39900.2 68.0602 -937107 2.84566e+07 19.49 39861 67.9698 -936140 2.84243e+07 19.51 39821.7 67.8787 -935178 2.83921e+07 - 19.53 39782.3 67.787 -934219 2.83602e+07 + 19.53 39782.3 67.7871 -934219 2.83602e+07 19.55 39742.9 67.6948 -933266 2.83284e+07 19.57 39703.4 67.602 -932316 2.82968e+07 19.59 39663.9 67.5087 -931371 2.82653e+07 @@ -984,7 +984,7 @@ 19.65 39545.1 67.2254 -928561 2.81721e+07 19.67 39505.3 67.1298 -927634 2.81413e+07 19.69 39465.6 67.0338 -926710 2.81108e+07 - 19.71 39425.8 66.9372 -925791 2.80804e+07 + 19.71 39425.8 66.9372 -925791 2.80805e+07 19.73 39385.9 66.8401 -924876 2.80503e+07 19.75 39346 66.7425 -923966 2.80203e+07 19.77 39306 66.6444 -923060 2.79905e+07 diff --git a/tests/integrate/186_PW_SNLKG_10D10S/refOnsager.txt b/tests/integrate/186_PW_SNLKG_10D10S/refOnsager.txt index 014f243860..586485603c 100644 --- a/tests/integrate/186_PW_SNLKG_10D10S/refOnsager.txt +++ b/tests/integrate/186_PW_SNLKG_10D10S/refOnsager.txt @@ -134,7 +134,7 @@ 2.65 328061 524.36 -6.61919e+06 1.83227e+08 2.67 327679 523.504 -6.6144e+06 1.83108e+08 2.69 327295 522.643 -6.60958e+06 1.82989e+08 - 2.71 326909 521.777 -6.60471e+06 1.82868e+08 + 2.71 326909 521.778 -6.60471e+06 1.82868e+08 2.73 326520 520.907 -6.59981e+06 1.82746e+08 2.75 326128 520.032 -6.59487e+06 1.82623e+08 2.77 325734 519.152 -6.58989e+06 1.825e+08 @@ -144,7 +144,7 @@ 2.85 324132 515.587 -6.56957e+06 1.81996e+08 2.87 323726 514.684 -6.56439e+06 1.81868e+08 2.89 323316 513.777 -6.55918e+06 1.81739e+08 - 2.91 322905 512.865 -6.55392e+06 1.81608e+08 + 2.91 322905 512.866 -6.55392e+06 1.81608e+08 2.93 322491 511.949 -6.54862e+06 1.81477e+08 2.95 322074 511.029 -6.54329e+06 1.81345e+08 2.97 321655 510.104 -6.53791e+06 1.81212e+08 @@ -163,12 +163,12 @@ 3.23 315987 497.701 -6.46439e+06 1.79395e+08 3.25 315535 496.718 -6.45845e+06 1.79249e+08 3.27 315080 495.732 -6.45247e+06 1.79101e+08 - 3.29 314622 494.742 -6.44646e+06 1.78953e+08 + 3.29 314622 494.743 -6.44646e+06 1.78953e+08 3.31 314163 493.749 -6.4404e+06 1.78804e+08 3.33 313701 492.752 -6.4343e+06 1.78653e+08 3.35 313236 491.751 -6.42816e+06 1.78502e+08 3.37 312770 490.746 -6.42198e+06 1.7835e+08 - 3.39 312301 489.737 -6.41576e+06 1.78196e+08 + 3.39 312301 489.738 -6.41576e+06 1.78196e+08 3.41 311830 488.726 -6.4095e+06 1.78042e+08 3.43 311356 487.71 -6.40319e+06 1.77887e+08 3.45 310881 486.691 -6.39685e+06 1.77731e+08 @@ -195,7 +195,7 @@ 3.87 300381 464.542 -6.25417e+06 1.74224e+08 3.89 299857 463.454 -6.24692e+06 1.74046e+08 3.91 299332 462.364 -6.23963e+06 1.73867e+08 - 3.93 298804 461.271 -6.2323e+06 1.73687e+08 + 3.93 298804 461.272 -6.2323e+06 1.73687e+08 3.95 298274 460.176 -6.22492e+06 1.73506e+08 3.97 297743 459.078 -6.21751e+06 1.73325e+08 3.99 297209 457.977 -6.21005e+06 1.73142e+08 @@ -237,10 +237,10 @@ 4.71 276719 416.949 -5.91396e+06 1.6589e+08 4.73 276117 415.778 -5.90497e+06 1.6567e+08 4.75 275513 414.607 -5.89594e+06 1.65449e+08 - 4.77 274908 413.433 -5.88687e+06 1.65227e+08 + 4.77 274908 413.434 -5.88687e+06 1.65227e+08 4.79 274301 412.259 -5.87775e+06 1.65004e+08 4.81 273692 411.084 -5.8686e+06 1.6478e+08 - 4.83 273082 409.907 -5.8594e+06 1.64554e+08 + 4.83 273082 409.908 -5.8594e+06 1.64554e+08 4.85 272470 408.73 -5.85017e+06 1.64328e+08 4.87 271856 407.551 -5.84089e+06 1.64101e+08 4.89 271241 406.372 -5.83157e+06 1.63873e+08 @@ -263,7 +263,7 @@ 5.23 260556 386.197 -5.667e+06 1.59841e+08 5.25 259914 385.005 -5.65696e+06 1.59595e+08 5.27 259271 383.812 -5.64689e+06 1.59348e+08 - 5.29 258627 382.619 -5.63677e+06 1.591e+08 + 5.29 258627 382.62 -5.63677e+06 1.591e+08 5.31 257982 381.426 -5.62661e+06 1.5885e+08 5.33 257335 380.233 -5.61642e+06 1.586e+08 5.35 256687 379.039 -5.60618e+06 1.58349e+08 @@ -339,7 +339,7 @@ 6.75 208874 296.491 -4.80188e+06 1.3848e+08 6.77 208168 295.348 -4.7893e+06 1.38166e+08 6.79 207461 294.207 -4.77669e+06 1.37852e+08 - 6.81 206755 293.067 -4.76406e+06 1.37537e+08 + 6.81 206755 293.068 -4.76406e+06 1.37537e+08 6.83 206048 291.93 -4.75141e+06 1.37221e+08 6.85 205341 290.794 -4.73873e+06 1.36905e+08 6.87 204634 289.659 -4.72603e+06 1.36588e+08 @@ -373,7 +373,7 @@ 7.43 184805 258.65 -4.36201e+06 1.27461e+08 7.45 184098 257.573 -4.34876e+06 1.27127e+08 7.47 183392 256.497 -4.33549e+06 1.26792e+08 - 7.49 182686 255.423 -4.32222e+06 1.26458e+08 + 7.49 182686 255.424 -4.32222e+06 1.26458e+08 7.51 181980 254.352 -4.30892e+06 1.26122e+08 7.53 181274 253.283 -4.29562e+06 1.25787e+08 7.55 180569 252.217 -4.2823e+06 1.25451e+08 @@ -402,7 +402,7 @@ 8.01 164482 228.344 -3.9735e+06 1.17622e+08 8.03 163790 227.336 -3.96001e+06 1.17279e+08 8.05 163098 226.33 -3.94651e+06 1.16935e+08 - 8.07 162408 225.327 -3.93301e+06 1.16591e+08 + 8.07 162408 225.328 -3.93301e+06 1.16591e+08 8.09 161718 224.327 -3.9195e+06 1.16247e+08 8.11 161029 223.33 -3.906e+06 1.15902e+08 8.13 160341 222.335 -3.89249e+06 1.15558e+08 @@ -709,7 +709,7 @@ 14.15 39903 57.215 -1.26479e+06 4.55096e+07 14.17 39822.4 57.0994 -1.26272e+06 4.54484e+07 14.19 39743.3 56.9862 -1.26068e+06 4.53881e+07 - 14.21 39665.7 56.8755 -1.25868e+06 4.53287e+07 + 14.21 39665.7 56.8756 -1.25868e+06 4.53287e+07 14.23 39589.7 56.7673 -1.25671e+06 4.52702e+07 14.25 39515.1 56.6615 -1.25478e+06 4.52127e+07 14.27 39441.9 56.5581 -1.25288e+06 4.5156e+07 @@ -747,7 +747,7 @@ 14.91 37813.8 54.438 -1.20844e+06 4.37762e+07 14.93 37783.2 54.4059 -1.20752e+06 4.37454e+07 14.95 37753.7 54.3757 -1.20662e+06 4.37152e+07 - 14.97 37725.4 54.3473 -1.20575e+06 4.36857e+07 + 14.97 37725.4 54.3473 -1.20575e+06 4.36858e+07 14.99 37698.2 54.3208 -1.2049e+06 4.36569e+07 15.01 37672 54.2961 -1.20408e+06 4.36287e+07 15.03 37646.9 54.2732 -1.20328e+06 4.36012e+07 @@ -757,7 +757,7 @@ 15.11 37557.1 54.1992 -1.20033e+06 4.34973e+07 15.13 37537.3 54.185 -1.19965e+06 4.34728e+07 15.15 37518.5 54.1726 -1.199e+06 4.3449e+07 - 15.17 37500.6 54.1617 -1.19837e+06 4.34257e+07 + 15.17 37500.6 54.1618 -1.19837e+06 4.34257e+07 15.19 37483.8 54.1526 -1.19776e+06 4.34031e+07 15.21 37468 54.1451 -1.19717e+06 4.3381e+07 15.23 37453.1 54.1391 -1.1966e+06 4.33595e+07 @@ -766,14 +766,14 @@ 15.29 37414.4 54.1309 -1.19504e+06 4.32985e+07 15.31 37403.4 54.1312 -1.19456e+06 4.32792e+07 15.33 37393.3 54.1331 -1.19411e+06 4.32605e+07 - 15.35 37384.1 54.1364 -1.19367e+06 4.32423e+07 + 15.35 37384.1 54.1364 -1.19367e+06 4.32424e+07 15.37 37375.9 54.1413 -1.19326e+06 4.32247e+07 15.39 37368.5 54.1476 -1.19286e+06 4.32076e+07 15.41 37362.1 54.1553 -1.19249e+06 4.3191e+07 15.43 37356.5 54.1645 -1.19213e+06 4.3175e+07 15.45 37351.9 54.175 -1.1918e+06 4.31594e+07 - 15.47 37348.1 54.1869 -1.19149e+06 4.31443e+07 - 15.49 37345.1 54.2002 -1.19119e+06 4.31298e+07 + 15.47 37348.1 54.187 -1.19149e+06 4.31443e+07 + 15.49 37345.1 54.2003 -1.19119e+06 4.31298e+07 15.51 37343 54.2149 -1.19092e+06 4.31157e+07 15.53 37341.8 54.2309 -1.19066e+06 4.31021e+07 15.55 37341.4 54.2481 -1.19042e+06 4.3089e+07 @@ -791,7 +791,7 @@ 15.79 37398.5 54.5502 -1.18896e+06 4.29666e+07 15.81 37408.2 54.5827 -1.18895e+06 4.29592e+07 15.83 37418.6 54.6162 -1.18895e+06 4.29521e+07 - 15.85 37429.7 54.6507 -1.18897e+06 4.29455e+07 + 15.85 37429.7 54.6508 -1.18897e+06 4.29455e+07 15.87 37441.5 54.6863 -1.18901e+06 4.29392e+07 15.89 37454 54.7228 -1.18906e+06 4.29333e+07 15.91 37467.2 54.7603 -1.18912e+06 4.29278e+07 @@ -824,15 +824,15 @@ 16.45 38055 56.0725 -1.196e+06 4.29002e+07 16.47 38084.3 56.1299 -1.19642e+06 4.2903e+07 16.49 38114.2 56.1877 -1.19685e+06 4.29061e+07 - 16.51 38144.5 56.2461 -1.19729e+06 4.29093e+07 + 16.51 38144.5 56.2461 -1.19729e+06 4.29094e+07 16.53 38175.3 56.3048 -1.19774e+06 4.29128e+07 16.55 38206.6 56.364 -1.1982e+06 4.29165e+07 - 16.57 38238.3 56.4236 -1.19867e+06 4.29204e+07 + 16.57 38238.3 56.4236 -1.19867e+06 4.29205e+07 16.59 38270.5 56.4836 -1.19915e+06 4.29246e+07 16.61 38303.1 56.5439 -1.19964e+06 4.29289e+07 16.63 38336.1 56.6047 -1.20014e+06 4.29334e+07 16.65 38369.5 56.6657 -1.20064e+06 4.29381e+07 - 16.67 38403.4 56.7271 -1.20116e+06 4.2943e+07 + 16.67 38403.4 56.7272 -1.20116e+06 4.2943e+07 16.69 38437.7 56.7889 -1.20168e+06 4.29481e+07 16.71 38472.4 56.8509 -1.20222e+06 4.29534e+07 16.73 38507.5 56.9133 -1.20276e+06 4.29588e+07 @@ -857,7 +857,7 @@ 17.11 39243.4 58.1357 -1.21445e+06 4.30907e+07 17.13 39285.4 58.2011 -1.21514e+06 4.30989e+07 17.15 39327.6 58.2666 -1.21582e+06 4.31072e+07 - 17.17 39370.2 58.3321 -1.21652e+06 4.31156e+07 + 17.17 39370.1 58.3321 -1.21652e+06 4.31156e+07 17.19 39413 58.3976 -1.21721e+06 4.3124e+07 17.21 39456 58.4631 -1.21792e+06 4.31326e+07 17.23 39499.4 58.5286 -1.21862e+06 4.31413e+07 @@ -869,7 +869,7 @@ 17.35 39764.8 58.9205 -1.22298e+06 4.31948e+07 17.37 39809.9 58.9855 -1.22372e+06 4.3204e+07 17.39 39855.2 59.0505 -1.22447e+06 4.32132e+07 - 17.41 39900.7 59.1153 -1.22522e+06 4.32225e+07 + 17.41 39900.7 59.1154 -1.22522e+06 4.32225e+07 17.43 39946.4 59.1801 -1.22597e+06 4.32318e+07 17.45 39992.4 59.2447 -1.22673e+06 4.32412e+07 17.47 40038.5 59.3092 -1.22749e+06 4.32507e+07 @@ -883,7 +883,7 @@ 17.63 40414.7 59.8189 -1.23372e+06 4.33277e+07 17.65 40462.5 59.8817 -1.23451e+06 4.33375e+07 17.67 40510.5 59.9443 -1.2353e+06 4.33473e+07 - 17.69 40558.6 60.0066 -1.2361e+06 4.33572e+07 + 17.69 40558.6 60.0067 -1.2361e+06 4.33572e+07 17.71 40606.9 60.0688 -1.2369e+06 4.3367e+07 17.73 40655.4 60.1307 -1.23771e+06 4.33769e+07 17.75 40703.9 60.1923 -1.23851e+06 4.33867e+07 @@ -899,7 +899,7 @@ 17.95 41196.8 60.7927 -1.24668e+06 4.34857e+07 17.97 41246.7 60.851 -1.24751e+06 4.34956e+07 17.99 41296.7 60.909 -1.24834e+06 4.35054e+07 - 18.01 41346.8 60.9666 -1.24917e+06 4.35153e+07 + 18.01 41346.8 60.9667 -1.24917e+06 4.35153e+07 18.03 41396.9 61.0239 -1.25e+06 4.35251e+07 18.05 41447.2 61.0808 -1.25083e+06 4.35349e+07 18.07 41497.5 61.1374 -1.25166e+06 4.35447e+07 @@ -931,13 +931,13 @@ 18.59 42820.3 62.4625 -1.27338e+06 4.37849e+07 18.61 42871.2 62.5073 -1.27421e+06 4.37934e+07 18.63 42922.1 62.5516 -1.27504e+06 4.38018e+07 - 18.65 42973 62.5954 -1.27586e+06 4.38101e+07 + 18.65 42973 62.5955 -1.27586e+06 4.38101e+07 18.67 43023.8 62.6388 -1.27669e+06 4.38184e+07 18.69 43074.6 62.6816 -1.27751e+06 4.38265e+07 18.71 43125.3 62.7239 -1.27833e+06 4.38346e+07 18.73 43176 62.7657 -1.27915e+06 4.38426e+07 18.75 43226.7 62.8071 -1.27997e+06 4.38506e+07 - 18.77 43277.2 62.8478 -1.28079e+06 4.38584e+07 + 18.77 43277.2 62.8479 -1.28079e+06 4.38584e+07 18.79 43327.8 62.8881 -1.2816e+06 4.38662e+07 18.81 43378.2 62.9279 -1.28241e+06 4.38739e+07 18.83 43428.6 62.9671 -1.28322e+06 4.38814e+07 @@ -948,7 +948,7 @@ 18.93 43679.3 63.1555 -1.28724e+06 4.3918e+07 18.95 43729.2 63.1916 -1.28803e+06 4.3925e+07 18.97 43779 63.2271 -1.28883e+06 4.3932e+07 - 18.99 43828.7 63.2621 -1.28962e+06 4.39388e+07 + 18.99 43828.7 63.2622 -1.28962e+06 4.39388e+07 19.01 43878.3 63.2966 -1.29041e+06 4.39455e+07 19.03 43927.8 63.3306 -1.29119e+06 4.39522e+07 19.05 43977.2 63.364 -1.29198e+06 4.39587e+07 @@ -971,18 +971,18 @@ 19.39 44795.9 63.8477 -1.30477e+06 4.40525e+07 19.41 44842.6 63.8712 -1.30549e+06 4.40569e+07 19.43 44889.2 63.8941 -1.3062e+06 4.40612e+07 - 19.45 44935.6 63.9164 -1.30691e+06 4.40653e+07 + 19.45 44935.6 63.9164 -1.30691e+06 4.40654e+07 19.47 44981.7 63.9382 -1.30762e+06 4.40694e+07 19.49 45027.7 63.9594 -1.30832e+06 4.40733e+07 19.51 45073.5 63.9801 -1.30901e+06 4.40771e+07 19.53 45119.1 64.0002 -1.30971e+06 4.40807e+07 - 19.55 45164.5 64.0197 -1.31039e+06 4.40842e+07 + 19.55 45164.5 64.0198 -1.31039e+06 4.40842e+07 19.57 45209.6 64.0387 -1.31107e+06 4.40876e+07 19.59 45254.6 64.0572 -1.31175e+06 4.40908e+07 19.61 45299.3 64.075 -1.31242e+06 4.40939e+07 19.63 45343.8 64.0923 -1.31309e+06 4.40968e+07 19.65 45388.1 64.1091 -1.31375e+06 4.40996e+07 - 19.67 45432.2 64.1252 -1.31441e+06 4.41022e+07 + 19.67 45432.2 64.1253 -1.31441e+06 4.41022e+07 19.69 45476 64.1409 -1.31506e+06 4.41047e+07 19.71 45519.6 64.1559 -1.31571e+06 4.41071e+07 19.73 45563 64.1704 -1.31635e+06 4.41093e+07 @@ -998,4 +998,4 @@ 19.93 45982.5 64.2848 -1.32244e+06 4.41229e+07 19.95 46023 64.2932 -1.32302e+06 4.41234e+07 19.97 46063.1 64.3011 -1.32359e+06 4.41237e+07 - 19.99 46103 64.3083 -1.32415e+06 4.41239e+07 + 19.99 46103 64.3084 -1.32415e+06 4.41239e+07