From 923b69739f3abc8c921996ae6ff987d26925e393 Mon Sep 17 00:00:00 2001 From: Yike Huang <67682086+kirk0830@users.noreply.github.com> Date: Wed, 4 Oct 2023 08:03:06 +0800 Subject: [PATCH] Fix: refactor of ESolver_KS_PW::init_after_vc with new psi initializer (#3031) * add re-initialization after vcrelax * append annotations in esolver_ks_pw * correct esolver_ks_pw.cpp oops, I deleted two lines by mistake, sorry! * improve line coverage of psi_initializer * remove newly added unittests, will re-design later * refactor psi_initializer, add unit tests, add integrated tests, update documents * update docs * eliminate warning about forgetting to add override keyword * fix compile bug from rename pseudo_nc to pseudo * mock the full UnitCell, undefined reference will not see anymore * delete redefinition of parameter and functions * Update input-main.md according to suggestion from Tianqi --- docs/advanced/input_files/input-main.md | 17 + source/module_esolver/esolver_ks_pw.cpp | 157 +++-- source/module_esolver/esolver_ks_pw.h | 4 +- source/module_psi/CMakeLists.txt | 2 - source/module_psi/psi_initializer.cpp | 163 ++--- source/module_psi/psi_initializer.h | 129 +++- source/module_psi/psi_initializer_atomic.cpp | 132 ++-- source/module_psi/psi_initializer_atomic.h | 36 +- .../psi_initializer_atomic_random.cpp | 18 +- .../psi_initializer_atomic_random.h | 28 +- source/module_psi/psi_initializer_nao.cpp | 171 ++--- source/module_psi/psi_initializer_nao.h | 33 +- .../module_psi/psi_initializer_nao_random.cpp | 18 +- .../module_psi/psi_initializer_nao_random.h | 25 +- source/module_psi/psi_initializer_random.cpp | 22 +- source/module_psi/psi_initializer_random.h | 29 +- source/module_psi/test/CMakeLists.txt | 11 +- source/module_psi/test/mock_unitcell.cpp | 107 +++ .../test/psi_initializer_integrated_test.cpp | 230 ------- .../test/psi_initializer_unit_test.cpp | 633 ++++++++++++++++++ tests/integrate/102_PW_PINT_RKS/INPUT | 22 + tests/integrate/102_PW_PINT_RKS/KPT | 4 + tests/integrate/102_PW_PINT_RKS/STRU | 24 + tests/integrate/102_PW_PINT_RKS/jd | 1 + tests/integrate/102_PW_PINT_RKS/result.ref | 3 + tests/integrate/102_PW_PINT_UKS/INPUT | 25 + tests/integrate/102_PW_PINT_UKS/KPT | 4 + tests/integrate/102_PW_PINT_UKS/STRU | 29 + tests/integrate/102_PW_PINT_UKS/jd | 1 + tests/integrate/102_PW_PINT_UKS/result.ref | 3 + tests/integrate/108_PW_RE_PINT_RKS/INPUT | 27 + tests/integrate/108_PW_RE_PINT_RKS/KPT | 4 + tests/integrate/108_PW_RE_PINT_RKS/STRU | 24 + tests/integrate/108_PW_RE_PINT_RKS/jd | 1 + tests/integrate/108_PW_RE_PINT_RKS/result.ref | 5 + tests/integrate/Autotest.sh | 3 +- tests/integrate/CASES | 3 + 37 files changed, 1549 insertions(+), 599 deletions(-) create mode 100644 source/module_psi/test/mock_unitcell.cpp delete mode 100644 source/module_psi/test/psi_initializer_integrated_test.cpp create mode 100644 source/module_psi/test/psi_initializer_unit_test.cpp create mode 100644 tests/integrate/102_PW_PINT_RKS/INPUT create mode 100644 tests/integrate/102_PW_PINT_RKS/KPT create mode 100644 tests/integrate/102_PW_PINT_RKS/STRU create mode 100644 tests/integrate/102_PW_PINT_RKS/jd create mode 100644 tests/integrate/102_PW_PINT_RKS/result.ref create mode 100644 tests/integrate/102_PW_PINT_UKS/INPUT create mode 100644 tests/integrate/102_PW_PINT_UKS/KPT create mode 100644 tests/integrate/102_PW_PINT_UKS/STRU create mode 100644 tests/integrate/102_PW_PINT_UKS/jd create mode 100644 tests/integrate/102_PW_PINT_UKS/result.ref create mode 100644 tests/integrate/108_PW_RE_PINT_RKS/INPUT create mode 100644 tests/integrate/108_PW_RE_PINT_RKS/KPT create mode 100644 tests/integrate/108_PW_RE_PINT_RKS/STRU create mode 100644 tests/integrate/108_PW_RE_PINT_RKS/jd create mode 100644 tests/integrate/108_PW_RE_PINT_RKS/result.ref diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 1d23a726fb..85b88e2378 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -12,6 +12,7 @@ - [kpar](#kpar) - [bndpar](#bndpar) - [latname](#latname) + - [psi\_initializer](#psi_initializer) - [init\_wfc](#init_wfc) - [init\_chg](#init_chg) - [init\_vel](#init_vel) @@ -463,16 +464,32 @@ These variables are used to control general system parameters. - triclinic: triclinic (14) - **Default**: none +### psi_initializer + +- **Type**: Integer +- **Description**: enable the experimental feature psi_initializer, to support use numerical atomic orbitals initialize wavefunction (`basis_type pw` case). + + NOTE: this feature is not well-implemented for `nspin 4` case (closed presently), and cannot use with `calculation nscf`/`esolver_type sdft` cases. + Available options are: + - 0: disable psi_initializer + - 1: enable psi_initializer +- **Default**: 0 + ### init_wfc - **Type**: String - **Description**: Only useful for plane wave basis only now. It is the name of the starting wave functions. In the future. we should also make this variable available for localized orbitals set. Available options are: + - atomic: from atomic pseudo wave functions. If they are not enough, other wave functions are initialized with random numbers. - atomic+random: add small random numbers on atomic pseudo-wavefunctions - file: from file - random: random numbers + + with `psi_initializer 1`, two more options are supported: + - nao: from numerical atomic orbitals. If they are not enough, other wave functions are initialized with random numbers. + - nao+random: add small random numbers on numerical atomic orbitals - **Default**: atomic ### init_chg diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index b99e79cdd6..ae07c26bf3 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -211,43 +211,7 @@ void ESolver_KS_PW::Init(Input& inp, UnitCell& ucell) } if (GlobalV::psi_initializer) { - if(GlobalV::init_wfc == "atomic") - { - this->psi_init = new psi_initializer_atomic(&(this->sf), this->pw_wfc); - // there are things only need to calculate once - //this->psi_init->set_pseudopot_files(GlobalC::ucell.pseudo_fn); - // not parallelized function, but we have GlobalC now, - // in the future once GlobalC is removed, we will parallelize this function - this->psi_init->cal_ovlp_pswfcjlq(); - } - else if(GlobalV::init_wfc == "random") - { - this->psi_init = new psi_initializer_random(&(this->sf), this->pw_wfc); - } - else if(GlobalV::init_wfc == "nao") - { - this->psi_init = new psi_initializer_nao(&(this->sf), this->pw_wfc); - // there are things only need to calculate once - this->psi_init->set_orbital_files(GlobalC::ucell.orbital_fn); - this->psi_init->cal_ovlp_flzjlq(); - } - else if(GlobalV::init_wfc == "atomic+random") - { - this->psi_init = new psi_initializer_atomic_random(&(this->sf), this->pw_wfc); - // there are things only need to calculate once - //this->psi_init->set_pseudopot_files(GlobalC::ucell.pseudo_fn); - // not parallelized function, but we have GlobalC now, - // in the future once GlobalC is removed, we will parallelize this function - this->psi_init->cal_ovlp_pswfcjlq(); - } - else if(GlobalV::init_wfc == "nao+random") - { - this->psi_init = new psi_initializer_nao_random(&(this->sf), this->pw_wfc); - // there are things only need to calculate once - this->psi_init->set_orbital_files(GlobalC::ucell.orbital_fn); - this->psi_init->cal_ovlp_flzjlq(); - } - else ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::init", "for new psi initializer, init_wfc type not supported"); + this->allocate_psi_init(); } // temporary this->Init_GlobalC(inp, ucell); @@ -316,8 +280,23 @@ void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) this->pw_wfc->nz); this->pw_wfc->initparameters(false, INPUT.ecutwfc, this->kv.nks, this->kv.kvec_d.data()); this->pw_wfc->collect_local_pw(inp.erf_ecut, inp.erf_height, inp.erf_sigma); - this->wf.init_after_vc(this->kv.nks); - this->wf.init_at_1(&this->sf); + if(GlobalV::psi_initializer) + { + if(GlobalV::init_wfc.substr(0, 3) == "nao") + { + this->psi_init->cal_ovlp_flzjlq(); // for nao, we recalculate the overlap matrix between flz and jlq + } + else if(GlobalV::init_wfc.substr(0, 6) == "atomic") + { + this->psi_init->cal_ovlp_pswfcjlq(); // 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 + } + else + { + this->wf.init_after_vc(this->kv.nks); // reallocate wanf2, the planewave expansion of lcao + this->wf.init_at_1(&this->sf); // re-calculate tab_at, the overlap matrix between atomic pswfc and jlq + } } #ifdef USE_PAW @@ -401,7 +380,22 @@ void ESolver_KS_PW::beforescf(int istep) */ if(GlobalV::psi_initializer) { - this->initialize_psi(); + /* + beforescf 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") + { + if(istep == 0) this->initialize_psi(); + } + else + { + this->initialize_psi(); + } } } @@ -465,6 +459,89 @@ void ESolver_KS_PW::eachiterinit(const int istep, const int iter) } } +template +void ESolver_KS_PW::allocate_psi_init() +{ + if(this->psi_init != nullptr) + { + delete this->psi_init; + this->psi_init = nullptr; + } + if((GlobalV::init_wfc.substr(0, 6) == "atomic")&&(GlobalC::ucell.natomwfc == 0)) + { + GlobalV::init_wfc = "random"; + std::cout << " WARNING: atomic pseudowavefunction is required but there is NOT ANY, set to random automatically." << std::endl; + #ifdef __MPI + this->psi_init = new psi_initializer_random(&(this->sf), this->pw_wfc, &(GlobalC::ucell), &(GlobalC::Pkpoints), INPUT.pw_seed); + #else + this->psi_init = new psi_initializer_random(&(this->sf), this->pw_wfc, &(GlobalC::ucell), INPUT.pw_seed); + #endif + this->psi_init->initialize_only_once(); + } + else + { + if(GlobalV::init_wfc == "atomic") + { + #ifdef __MPI + this->psi_init = new psi_initializer_atomic(&(this->sf), this->pw_wfc, &(GlobalC::ucell), &(GlobalC::Pkpoints), INPUT.pw_seed); + #else + this->psi_init = new psi_initializer_atomic(&(this->sf), this->pw_wfc, &(GlobalC::ucell), INPUT.pw_seed); + #endif + this->psi_init->initialize_only_once(&(GlobalC::ppcell)); + this->psi_init->cal_ovlp_pswfcjlq(); + } + else if(GlobalV::init_wfc == "random") + { + #ifdef __MPI + this->psi_init = new psi_initializer_random(&(this->sf), this->pw_wfc, &(GlobalC::ucell), &(GlobalC::Pkpoints), INPUT.pw_seed); + #else + this->psi_init = new psi_initializer_random(&(this->sf), this->pw_wfc, &(GlobalC::ucell), INPUT.pw_seed); + #endif + this->psi_init->initialize_only_once(); + } + else if(GlobalV::init_wfc == "nao") + { + if(GlobalV::NSPIN == 4) + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::allocate_psi_init", "for nao, soc this not safely implemented yet. To use it now, comment out this line."); + } + #ifdef __MPI + this->psi_init = new psi_initializer_nao(&(this->sf), this->pw_wfc, &(GlobalC::ucell), &(GlobalC::Pkpoints), INPUT.pw_seed); + #else + this->psi_init = new psi_initializer_nao(&(this->sf), this->pw_wfc, &(GlobalC::ucell), INPUT.pw_seed); + #endif + this->psi_init->set_orbital_files(GlobalC::ucell.orbital_fn); + this->psi_init->initialize_only_once(); + this->psi_init->cal_ovlp_flzjlq(); + } + else if(GlobalV::init_wfc == "atomic+random") + { + #ifdef __MPI + this->psi_init = new psi_initializer_atomic_random(&(this->sf), this->pw_wfc, &(GlobalC::ucell), &(GlobalC::Pkpoints), INPUT.pw_seed); + #else + this->psi_init = new psi_initializer_atomic_random(&(this->sf), this->pw_wfc, &(GlobalC::ucell), INPUT.pw_seed); + #endif + this->psi_init->initialize_only_once(&(GlobalC::ppcell)); + this->psi_init->cal_ovlp_pswfcjlq(); + } + else if(GlobalV::init_wfc == "nao+random") + { + if(GlobalV::NSPIN == 4) + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::allocate_psi_init", "for nao, soc this not safely implemented yet. To use it now, comment out this line."); + } + #ifdef __MPI + this->psi_init = new psi_initializer_nao_random(&(this->sf), this->pw_wfc, &(GlobalC::ucell), &(GlobalC::Pkpoints), INPUT.pw_seed); + #else + this->psi_init = new psi_initializer_nao_random(&(this->sf), this->pw_wfc, &(GlobalC::ucell), INPUT.pw_seed); + #endif + this->psi_init->set_orbital_files(GlobalC::ucell.orbital_fn); + this->psi_init->initialize_only_once(); + this->psi_init->cal_ovlp_flzjlq(); + } + else ModuleBase::WARNING_QUIT("ESolver_KS_PW::allocate_psi_init", "for new psi initializer, init_wfc type not supported"); + } +} /* 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 diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index 9f26407fa6..78441086da 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -83,7 +83,9 @@ namespace ModuleESolver void Init_GlobalC(Input& inp, UnitCell& cell); //calculate conductivities from j-j correlation function void calcondw(const int nt,const double dt, const double fwhmin, const double wcut, const double dw_in, double *ct11, double *ct12, double *ct22); - + /// @brief allocate psi_init the new psi_initializer + void allocate_psi_init(); + /// @brief initialize psi void initialize_psi(); private: psi_initializer* psi_init = nullptr; diff --git a/source/module_psi/CMakeLists.txt b/source/module_psi/CMakeLists.txt index 280fa5c272..cf9ebff49c 100644 --- a/source/module_psi/CMakeLists.txt +++ b/source/module_psi/CMakeLists.txt @@ -23,6 +23,4 @@ endif() if (BUILD_TESTING) add_subdirectory(kernels/test) add_subdirectory(test) - - add_compile_definitions(PSI_INITIALIZER_TEST) endif() \ No newline at end of file diff --git a/source/module_psi/psi_initializer.cpp b/source/module_psi/psi_initializer.cpp index 5a45ba1e81..48edb87813 100644 --- a/source/module_psi/psi_initializer.cpp +++ b/source/module_psi/psi_initializer.cpp @@ -5,10 +5,16 @@ #include "module_base/timer.h" // three global variables definition #include "module_base/global_variable.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -psi_initializer::psi_initializer(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in): sf(sf_in), pw_wfc(pw_wfc_in) +#ifdef __MPI +psi_initializer::psi_initializer(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in) + : sf(sf_in), pw_wfc(pw_wfc_in), p_ucell(p_ucell_in), p_parakpts(p_parakpts_in), random_seed(random_seed_in) +#else +psi_initializer::psi_initializer(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in) + : sf(sf_in), pw_wfc(pw_wfc_in), p_ucell(p_ucell_in), random_seed(random_seed_in) +#endif { + if(this->p_ucell == nullptr) ModuleBase::WARNING_QUIT("psi_initializer", "interface to UnitCell is not valid, quit!"); this->ixy2is = new int[this->pw_wfc->fftnxy]; this->pw_wfc->getfftixy2is(this->ixy2is); } @@ -26,9 +32,10 @@ psi::Psi>* psi_initializer::allocate() /* WARNING: when basis_type = "pw", the variable GlobalV::NLOCAL will also be set, in this case, it is set to 9 = 1 + 3 + 5, which is the maximal number of orbitals spd, I don't think it is reasonable - The way of calculating GlobalC::ucell.natomwfc is, for each atom, read pswfc and for s, it is 1, for p, it is 3 + The way of calculating this->p_ucell->natomwfc is, for each atom, read pswfc and for s, it is 1, for p, it is 3 , then multiplied by the number of atoms, and then add them together. */ + if(this->psig != nullptr) delete this->psig; int prefactor = 1; int nbands_actual = 0; @@ -41,15 +48,15 @@ psi::Psi>* psi_initializer::allocate() { if(GlobalV::init_wfc.substr(0, 6) == "atomic") { - if(GlobalC::ucell.natomwfc >= GlobalV::NBANDS) + if(this->p_ucell->natomwfc >= GlobalV::NBANDS) { - nbands_actual = GlobalC::ucell.natomwfc; + nbands_actual = this->p_ucell->natomwfc; this->nbands_complem = 0; } else { nbands_actual = GlobalV::NBANDS; - this->nbands_complem = GlobalV::NBANDS - GlobalC::ucell.natomwfc; + this->nbands_complem = GlobalV::NBANDS - this->p_ucell->natomwfc; } } else if(GlobalV::init_wfc.substr(0, 3) == "nao") @@ -58,15 +65,26 @@ psi::Psi>* psi_initializer::allocate() previously GlobalV::NLOCAL is used here, however it is wrong. GlobalV::NLOCAL is fixed to 9*nat. */ int nbands_local = 0; - for(int it = 0; it < GlobalC::ucell.ntype; it++) + for(int it = 0; it < this->p_ucell->ntype; it++) { - for(int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + for(int ia = 0; ia < this->p_ucell->atoms[it].na; ia++) { /* FOR EVERY ATOM */ - for(int l = 0; l < GlobalC::ucell.atoms[it].nwl + 1; l++) + for(int l = 0; l < this->p_ucell->atoms[it].nwl + 1; l++) { - /* EVERY ZETA FOR (2l+1) ORBS, for NSPIN = 4, DOUBLE */ - nbands_local += GlobalC::ucell.atoms[it].l_nchi[l]*(2*l+1) * GlobalV::NPOL; + /* EVERY ZETA FOR (2l+1) ORBS, for NSPIN = 4, 4 times (because change to rotate base) */ + nbands_local += this->p_ucell->atoms[it].l_nchi[l]*(2*l+1) * GlobalV::NPOL; + /* the following is for rotate base. However it cannot yield correct result because many zeros yielded and cause diag failure */ + /* + if(l == 0) + { + nbands_local += this->p_ucell->atoms[it].l_nchi[l]*(2*l+1) * GlobalV::NPOL; + } + else + { + nbands_local += this->p_ucell->atoms[it].l_nchi[l]*(2*l+1) * GlobalV::NPOL * GlobalV::NPOL; + } + */ } } } @@ -82,7 +100,6 @@ psi::Psi>* psi_initializer::allocate() } } } - int nkpts_actual = (GlobalV::CALCULATION == "nscf" && this->mem_saver == 1)? 1 : this->pw_wfc->nks; int nbasis_actual = this->pw_wfc->npwk_max * GlobalV::NPOL; @@ -97,97 +114,35 @@ psi::Psi>* psi_initializer::allocate() nbands_actual, nbasis_actual, this->pw_wfc->npwk); - const size_t memory_cost = - nkpts_actual* + const size_t memory_cost_psi = nkpts_actual* - this->pw_wfc->npwk_max * GlobalV::NPOL* + GlobalV::NBANDS * this->pw_wfc->npwk_max * GlobalV::NPOL* sizeof(std::complex); - std::cout << " MEMORY FOR PSI (MB) : " << double(memory_cost)/1024.0/1024.0 << std::endl; - ModuleBase::Memory::record("Psi_PW", memory_cost); + std::cout << " MEMORY FOR PSI PER PROCESSOR (MB) : " << double(memory_cost_psi)/1024.0/1024.0 << std::endl; + const size_t memory_cost_psig = + nkpts_actual* + nbands_actual * this->pw_wfc->npwk_max * GlobalV::NPOL* + sizeof(std::complex); + std::cout << " MEMORY FOR AUXILLARY PSI PER PROCESSOR (MB) : " << double(memory_cost_psig)/1024.0/1024.0 << std::endl; + ModuleBase::Memory::record("Psi_PW", memory_cost_psi); + ModuleBase::Memory::record("PsiG_PW", memory_cost_psig); ModuleBase::timer::tick("psi_initializer", "allocate"); return psi_out; } -void psi_initializer::write_psig() const -{ - for(int ik = 0; ik < this->psig->get_nk(); ik++) - { - std::string filename = "psig_"+std::to_string(ik); - std::ofstream ofs_psig; - ofs_psig.open(filename+"_kpt.out"); - ofs_psig << "N.B.: output data is complex, therefore every data will be enclosed by parenthesis." << std::endl; - ofs_psig << "psig information" << std::endl; - ofs_psig << "number of kpoints: " << this->psig->get_nk() << std::endl; - ofs_psig << "number of bands: " << this->psig->get_nbands() << std::endl; - ofs_psig << "number of planewaves: " << this->psig->get_nbasis() << std::endl; - ofs_psig << "Calculation information" << std::endl; - ofs_psig << "method of psi initialization: " << GlobalV::init_wfc << std::endl; - ofs_psig << "method of diagonalization: " << GlobalV::KS_SOLVER << std::endl; - this->psig->fix_k(ik); - ofs_psig << "k point No. " << ik << std::endl; - for(int iband = 0; iband < this->psig->get_nbands(); iband++) - { - ofs_psig << "energy band No. " << iband << std::endl; - for(int ibasis = 0; ibasis < this->psig->get_nbasis(); ibasis++) - { - ofs_psig << std::setprecision(10) << std::fixed << (*(this->psig))(iband, ibasis) << " "; - } - ofs_psig << std::endl; - } - ofs_psig << std::endl; - ofs_psig.close(); - } -} - -void psi_initializer::write_psig(int ik) const -{ - std::string filename = "psig_"+std::to_string(ik); - std::ofstream ofs_psig; - ofs_psig.open(filename+"_kpt.out"); - ofs_psig << "N.B.: output data is complex, therefore every data will be enclosed by parenthesis." << std::endl; - ofs_psig << "psig information" << std::endl; - ofs_psig << "number of kpoints: " << this->psig->get_nk() << std::endl; - ofs_psig << "number of bands: " << this->psig->get_nbands() << std::endl; - ofs_psig << "number of planewaves: " << this->psig->get_nbasis() << std::endl; - ofs_psig << "Calculation information" << std::endl; - ofs_psig << "method of psi initialization: " << GlobalV::init_wfc << std::endl; - ofs_psig << "method of diagonalization: " << GlobalV::KS_SOLVER << std::endl; - this->psig->fix_k(ik); - ofs_psig << "k point No. " << ik << std::endl; - for(int iband = 0; iband < this->psig->get_nbands(); iband++) - { - ofs_psig << "energy band No. " << iband << std::endl; - for(int ibasis = 0; ibasis < this->psig->get_nbasis(); ibasis++) - { - ofs_psig << std::setprecision(10) << std::fixed << (*(this->psig))(iband, ibasis) << " "; - } - ofs_psig << std::endl; - } - ofs_psig << std::endl; - ofs_psig.close(); -} - -void psi_initializer::print_status(psi::Psi>& psi) const -{ - std::cout << "Current method: " << this->method << std::endl; - std::cout << "Psi status:" << std::endl; - std::cout << " number of kpoints: " << psi.get_nk() << std::endl; - std::cout << " number of bands: " << psi.get_nbands() << std::endl; - std::cout << " number of planewaves: " << psi.get_nbasis() << std::endl; -} -void psi_initializer::random_t(std::complex* psi, const int iw_start, const int iw_end, const int ik, const ModulePW::PW_Basis_K* wfc_basis) +void psi_initializer::random_t(std::complex* psi, const int iw_start, const int iw_end, const int ik) { ModuleBase::timer::tick("psi_initializer", "random_t"); assert(iw_start >= 0); - const int ng = wfc_basis->npwk[ik]; + const int ng = this->pw_wfc->npwk[ik]; #ifdef __MPI - if (INPUT.pw_seed > 0) // qianrui add 2021-8-13 + if (this->random_seed > 0) // qianrui add 2021-8-13 { - srand(unsigned(INPUT.pw_seed + GlobalC::Pkpoints.startk_pool[GlobalV::MY_POOL] + ik)); - const int nxy = wfc_basis->fftnxy; - const int nz = wfc_basis->nz; - const int nstnz = wfc_basis->nst*nz; + srand(unsigned(this->random_seed + this->p_parakpts->startk_pool[GlobalV::MY_POOL] + ik)); + const int nxy = this->pw_wfc->fftnxy; + const int nz = this->pw_wfc->nz; + const int nstnz = this->pw_wfc->nst*nz; double *stickrr = new double[nz]; double *stickarg = new double[nz]; @@ -203,7 +158,7 @@ void psi_initializer::random_t(std::complex* psi, const int iw_start, co for(int ir=0; ir < nxy; ir++) { - if(wfc_basis->fftixy2ip[ir] < 0) continue; + if(this->pw_wfc->fftixy2ip[ir] < 0) continue; if(GlobalV::RANK_IN_POOL==0) { for(int iz=0; iz* psi, const int iw_start, co stickarg[ iz ] = std::rand()/double(RAND_MAX); } } - stick_to_pool(stickrr, ir, tmprr, wfc_basis); - stick_to_pool(stickarg, ir, tmparg, wfc_basis); + stick_to_pool(stickrr, ir, tmprr); + stick_to_pool(stickarg, ir, tmparg); } for (int ig = 0;ig < ng;ig++) { - const double rr = tmprr[wfc_basis->getigl2isz(ik,ig)]; - const double arg= ModuleBase::TWO_PI * tmparg[wfc_basis->getigl2isz(ik,ig)]; - const double gk2 = wfc_basis->getgk2(ik,ig); + const double rr = tmprr[this->pw_wfc->getigl2isz(ik,ig)]; + const double arg= ModuleBase::TWO_PI * tmparg[this->pw_wfc->getigl2isz(ik,ig)]; + const double gk2 = this->pw_wfc->getgk2(ik,ig); psi_slice[ig+startig] = std::complex(rr * cos(arg), rr * sin(arg)) / double(gk2 + 1.0); } startig += this->pw_wfc->npwk_max; @@ -234,9 +189,9 @@ void psi_initializer::random_t(std::complex* psi, const int iw_start, co else { #else // !__MPI - if (INPUT.pw_seed > 0) // qianrui add 2021-8-13 + if (this->random_seed > 0) // qianrui add 2021-8-13 { - srand(unsigned(INPUT.pw_seed + ik)); + srand(unsigned(this->random_seed + ik)); } #endif for (int iw = iw_start ;iw < iw_end; iw++) @@ -246,7 +201,7 @@ void psi_initializer::random_t(std::complex* psi, const int iw_start, co { const double rr = std::rand()/double(RAND_MAX); //qianrui add RAND_MAX const double arg= ModuleBase::TWO_PI * std::rand()/double(RAND_MAX); - const double gk2 = wfc_basis->getgk2(ik,ig); + const double gk2 = this->pw_wfc->getgk2(ik,ig); psi_slice[ig] = std::complex(rr * cos(arg), rr * sin(arg)) / double(gk2 + 1.0); } if(GlobalV::NPOL==2) @@ -255,7 +210,7 @@ void psi_initializer::random_t(std::complex* psi, const int iw_start, co { const double rr = std::rand()/double(RAND_MAX); const double arg= ModuleBase::TWO_PI * std::rand()/double(RAND_MAX); - const double gk2 = wfc_basis->getgk2(ik,ig-this->pw_wfc->npwk_max); + const double gk2 = this->pw_wfc->getgk2(ik,ig-this->pw_wfc->npwk_max); psi_slice[ig] = std::complex(rr * cos(arg), rr * sin(arg)) / double(gk2 + 1.0); } } @@ -267,13 +222,13 @@ void psi_initializer::random_t(std::complex* psi, const int iw_start, co } #ifdef __MPI -void psi_initializer::stick_to_pool(double* stick, const int& ir, double* out, const ModulePW::PW_Basis_K* wfc_basis) const +void psi_initializer::stick_to_pool(double* stick, const int& ir, double* out) const { ModuleBase::timer::tick("psi_initializer", "stick_to_pool"); MPI_Status ierror; const int is = this->ixy2is[ir]; - const int ip = wfc_basis->fftixy2ip[ir]; - const int nz = wfc_basis->nz; + const int ip = this->pw_wfc->fftixy2ip[ir]; + const int nz = this->pw_wfc->nz; if(ip == 0 && GlobalV::RANK_IN_POOL ==0) { diff --git a/source/module_psi/psi_initializer.h b/source/module_psi/psi_initializer.h index 9d8feb7262..241fd801d8 100644 --- a/source/module_psi/psi_initializer.h +++ b/source/module_psi/psi_initializer.h @@ -4,6 +4,7 @@ #include "module_psi/psi.h" // for psi data structure #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" #include "module_basis/module_pw/pw_basis_k.h" // for kpoint related data structure +#include "module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h" // numerical algorithm support #include "module_base/spherical_bessel_transformer.h" // for spherical bessel transform #ifdef __MPI @@ -33,23 +34,31 @@ class psi_initializer { public: /// @brief default constructor of psi initializer - psi_initializer() : sf(nullptr), pw_wfc(nullptr) { }; - /// @brief parameterized constructor of psi initializer - /// @param sf_in structure factor pointer - /// @param pw_wfc_in pw_basis_k pointer - psi_initializer(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in); + psi_initializer() { }; + #ifdef __MPI + /// @brief parameterized constructor of psi initializer (with MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param p_parakpts_in interface, link with Parallel_Kpoints GlobalC::Pkpoints + /// @param random_seed_in random seed + psi_initializer(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in = 1); + #else + /// @brief parameterized constructor of psi initializer (without MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param random_seed_in random seed + psi_initializer(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in = 1); + #endif /// @brief destructor virtual ~psi_initializer(); + // shared methods /// @brief allocate memory for psi /// @return pointer to psi, memory allocated psi::Psi>* allocate(); - /// @brief calculate psi in planewave representation - /// @param psi psi - /// @param ik index of kpoint - virtual psi::Psi>* cal_psig(int ik) = 0; - /// @brief get method of initializing psi /// @return the method std::string get_method() const { return this->method; } @@ -57,19 +66,12 @@ class psi_initializer /// @return nbands_complem int get_nbands_complem() const { return this->nbands_complem; } - void print_status(psi::Psi>& psi) const; - /// @brief set method manually /// @param method_in initialization method void set_method(std::string method_in) { this->method = method_in; } /// @brief set number of complementary bands /// @param nbands_in nbands_complem void set_nbands_complem(int nbands_in) { this->nbands_complem = nbands_in; } - /// @brief output psig to file, given number of kpoints, bands and basis, diagonalization method. Note: because it is complex number, therefore every number has format (real, imag) - void write_psig() const; - /// @brief output psig at given kpoint to file - void write_psig(int ik) const; - // virtual functions, will be implemented in derived classes // random to complement bands not initialized by pswfc or nao, therefore it is a basic function, or psi_initializer_random will be inherented by all other methods. /// @brief kernel to generate and assign random number for psi @@ -77,26 +79,35 @@ class psi_initializer /// @param iw_start index of wavefunction (start), the index of first band to initialize /// @param iw_end index of wavefunction (end) /// @param ik index of kpoint - /// @param wfc_basis because do not want to change contents in it, this parameter is reserved. should give this->pw_wfc - void random_t(std::complex* psi, const int iw_start, const int iw_end, const int ik, const ModulePW::PW_Basis_K* wfc_basis); - + void random_t(std::complex* psi, const int iw_start, const int iw_end, const int ik); // random /// @brief wrapper of random_t /// @param psi for psi::Psi, first use psi.fix(ik), then psi.get_pointer() to pass to parameter list /// @param iw_start index of wavefunction (start), the index of first band to initialize /// @param iw_end index of wavefunction (end) /// @param ik index of kpoint - /// @param wfc_basis because do not want to change contents in it, this parameter is reserved. should give this->pw_wfc virtual void random(std::complex* psi, const int iw_start, const int iw_end, - const int ik, const ModulePW::PW_Basis_K* wfc_basis) { ModuleBase::WARNING_QUIT("psi_initializer::random", "Polymorphism error"); } + const int ik) { ModuleBase::WARNING_QUIT("psi_initializer::random", "Polymorphism error"); } #ifdef __MPI /// @brief (about planewaves distribution) from stick mapping to pool /// @param stick /// @param ir /// @param out - /// @param wfc_basis because do not want to change contents in it, this parameter is reserved. should give this->pw_wfc - void stick_to_pool(double* stick, const int& ir, double* out, const ModulePW::PW_Basis_K* wfc_basis) const; + void stick_to_pool(double* stick, const int& ir, double* out) const; #endif + + // mutual methods, virtual, will be implemented differently in derived classes + /// @brief create table for storing calculated overlap between pseudowavefunction/numerical orbitals with spherical bessel function + virtual void create_ovlp_Xjlq() { ModuleBase::WARNING_QUIT("psi_initializer::create_ovlp_table", "Polymorphism error"); } + /// @brief calculate psi in planewave representation + /// @param psi psi + /// @param ik index of kpoint + virtual psi::Psi>* cal_psig(int ik) = 0; + + /// @brief for variables can be only initialized for once. + /// @param p_pspot_nl_in (for atomic) interfaces to pseudopot_cell_vnl object, in GlobalC now + /// @attention if one variable is necessary for all methods, initialize it in constructor, not here. + virtual void initialize_only_once(pseudopot_cell_vnl* p_pspot_nl_in = nullptr) = 0; // atomic /// @brief setter of pseudopotential files, useful when init_wfc = atomic virtual void set_pseudopot_files(std::string* pseudopot_files) { ModuleBase::WARNING_QUIT("psi_initializer::set_pseudopot_files", "Polymorphism error"); } @@ -110,6 +121,7 @@ class psi_initializer /// @param mode if 1, return cos(arg), 0, return cos(arg)+isin(arg), -1, return sin(arg) /// @return it depends virtual std::complex phase_factor(double arg, int mode = 0) { ModuleBase::WARNING_QUIT("psi_initializer::phase_factor", "Polymorphism error"); return std::complex(0.0,0.0);} + /// @brief calculate overlap table between pseudowavefunction and spherical bessel function virtual void cal_ovlp_pswfcjlq() { ModuleBase::WARNING_QUIT("psi_initializer::calc_ovlp_pswfcjlq", "Polymorphism error"); } // nao @@ -130,16 +142,77 @@ class psi_initializer int* get_ixy2is() const { return this->ixy2is; } /// @brief setter of ixy2is, the mapping from fftixy to stick index void set_ixy2is(int* ixy2is_in) { this->ixy2is = ixy2is_in; } + + /// @brief setter of p_ucell + /// @param p_ucell_in UnitCell pointer + void set_interface_ucell(UnitCell* p_ucell_in) { this->p_ucell = p_ucell_in; } + /// @brief getter of p_ucell + /// @return this->p_ucell + UnitCell* get_interface_ucell() const { return this->p_ucell; } + #ifdef __MPI + /// @brief setter of p_parakpts + /// @param p_parakpts_in Parallel_Kpoints pointer + void set_interface_parakpts(Parallel_Kpoints* p_parakpts_in) { this->p_parakpts = p_parakpts_in; } + /// @brief getter of p_parakpts + /// @return this->p_parakpts + Parallel_Kpoints* get_interface_parakpts() const { return this->p_parakpts; } + #endif + /// @brief setter of p_pspot_nl + /// @param p_pspot_nl_in pseudopot_cell_vnl pointer + void set_interface_pspot_nl(pseudopot_cell_vnl* p_pspot_nl_in) { this->p_pspot_nl = p_pspot_nl_in; } + /// @brief getter of p_pspot_nl + /// @return this->p_pspot_nl + pseudopot_cell_vnl* get_interface_pspot_nl() const { return this->p_pspot_nl; } + /// @brief setter of sf + /// @param sf_in Structure_Factor pointer + void set_interface_sf(Structure_Factor* sf_in) { this->sf = sf_in; } + /// @brief getter of sf + /// @return this->sf + Structure_Factor* get_interface_sf() const { return this->sf; } + /// @brief setter of pw_wfc + /// @param pw_wfc_in ModulePW::PW_Basis_K pointer + void set_interface_pw_wfc(ModulePW::PW_Basis_K* pw_wfc_in) { this->pw_wfc = pw_wfc_in; } + /// @brief getter of pw_wfc + /// @return this->pw_wfc + ModulePW::PW_Basis_K* get_interface_pw_wfc() const { return this->pw_wfc; } + /// @brief setter of random_seed + /// @param random_seed_in new value of random_seed + void set_random_seed(const int random_seed_in) { this->random_seed = random_seed_in; } + /// @brief getter of random_seed + /// @return this->random_seed + int get_random_seed() const { return this->random_seed; } + /// @brief setter of mem_saver + /// @param mem_saver_in new value of mem_saver + void set_mem_saver(const int mem_saver_in) { this->mem_saver = mem_saver_in; } + /// @brief getter of mem_saver + /// @return this->mem_saver + int get_mem_saver() const { return this->mem_saver; } // member variables /// @brief interface to the psi::Psi data structure class - psi::Psi>* psig; + psi::Psi>* psig = nullptr; + + protected: + // interfaces + // ATTENTION: DO NOT USE DELETE ON THESE POINTERS + // normal interfaces /// @brief interface to the Structure_Factor method class - Structure_Factor* sf; + Structure_Factor* sf = nullptr; /// @brief interface to the PW_Basis_K data structure class - ModulePW::PW_Basis_K* pw_wfc; + ModulePW::PW_Basis_K* pw_wfc = nullptr; + // interfaces designed to get rid of dependence troubles of GlobalC in unittest + /// @brief interface to the UnitCell data carrier class, used in all methods + UnitCell* p_ucell = nullptr; + #ifdef __MPI + /// @brief interface to the Parallel_Kpoints method class, used in all methods + Parallel_Kpoints* p_parakpts = nullptr; + #endif + /// @brief interface to the pseudopot_cell_vnl data carrier class, used in atomic + pseudopot_cell_vnl* p_pspot_nl = nullptr; + + int random_seed = 1; // random seed + // tool interfaces /// @brief method of Spherical Bessel Transformation ModuleBase::SphericalBesselTransformer sbt; // useful for atomic-like methods - private: // basic properties int mem_saver = 0; // will deprecated this variable soon diff --git a/source/module_psi/psi_initializer_atomic.cpp b/source/module_psi/psi_initializer_atomic.cpp index e8c93f803d..ad10918455 100644 --- a/source/module_psi/psi_initializer_atomic.cpp +++ b/source/module_psi/psi_initializer_atomic.cpp @@ -10,41 +10,51 @@ #include "module_base/timer.h" // three global variables definition #include "module_base/global_variable.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -psi_initializer_atomic::psi_initializer_atomic(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in) : psi_initializer(sf_in, pw_wfc_in) +#ifdef __MPI +psi_initializer_atomic::psi_initializer_atomic(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in) + : psi_initializer(sf_in, pw_wfc_in, p_ucell_in, p_parakpts_in, random_seed_in) +#else +psi_initializer_atomic::psi_initializer_atomic(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in) + : psi_initializer(sf_in, pw_wfc_in, p_ucell_in, random_seed_in) +#endif { this->set_method("atomic"); - // find correct dimension for ovlp_flzjlq - int dim1 = GlobalC::ucell.ntype; - int dim2 = 0; // dim2 should be the maximum number of pseudo atomic orbitals - for (int it = 0; it < GlobalC::ucell.ntype; it++) - { - dim2 = (GlobalC::ucell.atoms[it].ncpp.nchi > dim2) ? GlobalC::ucell.atoms[it].ncpp.nchi : dim2; - } - if (dim2 == 0) - { - ModuleBase::WARNING_QUIT("psi_initializer_atomic::psi_initializer_atomic", "there is not ANY pseudo atomic orbital read in present system, recommand other methods, quit."); - } - int dim3 = GlobalV::NQX; - // allocate memory for ovlp_flzjlq - this->ovlp_pswfcjlq.create(dim1, dim2, dim3); - this->ovlp_pswfcjlq.zero_out(); } - psi_initializer_atomic::~psi_initializer_atomic() {} +/* I leave this function here for deprecation of UnitCell in the future */ +/* void psi_initializer_atomic::set_pseudopot_files(std::string* pseudopot_files) { ModuleBase::timer::tick("psi_initializer_atomic", "set_pseudopot_files"); - for (int itype = 0; itype < GlobalC::ucell.ntype; itype++) + for (int itype = 0; itype < this->p_ucell->ntype; itype++) { this->pseudopot_files.push_back(pseudopot_files[itype]); } ModuleBase::timer::tick("psi_initializer_atomic", "set_pseudopot_files"); } +*/ +void psi_initializer_atomic::create_ovlp_Xjlq() +{ + // find correct dimension for ovlp_flzjlq + int dim1 = this->p_ucell->ntype; + int dim2 = 0; // dim2 should be the maximum number of pseudo atomic orbitals + for (int it = 0; it < this->p_ucell->ntype; it++) + { + dim2 = (this->p_ucell->atoms[it].ncpp.nchi > dim2) ? this->p_ucell->atoms[it].ncpp.nchi : dim2; + } + if (dim2 == 0) + { + ModuleBase::WARNING_QUIT("psi_initializer_atomic::create_ovlp_Xjlq", "there is not ANY pseudo atomic orbital read in present system, recommand other methods, quit."); + } + int dim3 = GlobalV::NQX; + // allocate memory for ovlp_flzjlq + this->ovlp_pswfcjlq.create(dim1, dim2, dim3); + this->ovlp_pswfcjlq.zero_out(); +} void psi_initializer_atomic::normalize_pswfc(int n_rgrid, double* pswfc, double* rab) { @@ -63,6 +73,19 @@ void psi_initializer_atomic::normalize_pswfc(int n_rgrid, double* pswfc, double* ModuleBase::timer::tick("psi_initializer_atomic", "normalize_pswfc"); } +void psi_initializer_atomic::initialize_only_once(pseudopot_cell_vnl* p_pspot_nl_in) +{ + ModuleBase::timer::tick("psi_initializer_atomic", "initialize_only_once"); + if(p_pspot_nl_in == nullptr) + { + ModuleBase::WARNING_QUIT("psi_initializer_atomic::initialize_only_once", "pseudopot_cell_vnl object cannot be mullptr for atomic, quit."); + } + this->p_pspot_nl = p_pspot_nl_in; + this->create_ovlp_Xjlq(); + //this->set_pseudopot_files(); + //this->cal_ovlp_pswfcjlq(); //because GlobalV::NQX will change during vcrelax, so it should be called in both init and init_after_vc + ModuleBase::timer::tick("psi_initializer_atomic", "initialize_only_once"); +} void psi_initializer_atomic::cal_ovlp_pswfcjlq() { @@ -73,20 +96,20 @@ void psi_initializer_atomic::cal_ovlp_pswfcjlq() { qgrid[iq] = GlobalV::DQ * iq; } - for (int it=0; itp_ucell->ntype; it++) { - maxn_rgrid = (GlobalC::ucell.atoms[it].ncpp.msh > maxn_rgrid) ? GlobalC::ucell.atoms[it].ncpp.msh : maxn_rgrid; + maxn_rgrid = (this->p_ucell->atoms[it].ncpp.msh > maxn_rgrid) ? this->p_ucell->atoms[it].ncpp.msh : maxn_rgrid; } ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"max mesh points in Pseudopotential",maxn_rgrid); - const double pref = ModuleBase::FOUR_PI / sqrt(GlobalC::ucell.omega); + const double pref = ModuleBase::FOUR_PI / sqrt(this->p_ucell->omega); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"dq(describe PAO in reciprocal space)",GlobalV::DQ); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"max q",GlobalV::NQX); - for (int it=0; itp_ucell->ntype; it++) { - Atom* atom = &GlobalC::ucell.atoms[it]; + Atom* atom = &this->p_ucell->atoms[it]; GlobalV::ofs_running<<"\n number of pseudo atomic orbitals for "<label<<" is "<< atom->ncpp.nchi << std::endl; @@ -131,11 +154,11 @@ std::complex psi_initializer_atomic::phase_factor(double arg, int mode) psi::Psi>* psi_initializer_atomic::cal_psig(int ik) { - ModuleBase::timer::tick("psi_initializer_atomic", "initialize"); + ModuleBase::timer::tick("psi_initializer_atomic", "cal_psig"); this->psig->fix_k(ik); //this->print_status(psi); const int npw = this->pw_wfc->npwk[ik]; - int lmax = GlobalC::ucell.lmax_ppwf; + int lmax = this->p_ucell->lmax_ppwf; const int total_lm = (lmax + 1) * (lmax + 1); ModuleBase::matrix ylm(total_lm, npw); std::complex *aux = new std::complex[npw]; @@ -149,34 +172,34 @@ psi::Psi>* psi_initializer_atomic::cal_psig(int ik) ModuleBase::YlmReal::Ylm_Real(total_lm, npw, gk, ylm); int index = 0; double *ovlp_pswfcjlg = new double[npw]; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < this->p_ucell->ntype; it++) { - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + for (int ia = 0; ia < this->p_ucell->atoms[it].na; ia++) { /* FOR EVERY ATOM */ std::complex *sk = this->sf->get_sk(ik, it, ia, this->pw_wfc); - for (int ipswfc = 0; ipswfc < GlobalC::ucell.atoms[it].ncpp.nchi; ipswfc++) + for (int ipswfc = 0; ipswfc < this->p_ucell->atoms[it].ncpp.nchi; ipswfc++) { /* FOR EVERY PSWFC OF ATOM */ - if (GlobalC::ucell.atoms[it].ncpp.oc[ipswfc] >= 0.0) + if (this->p_ucell->atoms[it].ncpp.oc[ipswfc] >= 0.0) { /* IF IS OCCUPIED, GET L */ - const int l = GlobalC::ucell.atoms[it].ncpp.lchi[ipswfc]; + const int l = this->p_ucell->atoms[it].ncpp.lchi[ipswfc]; std::complex lphase = pow(ModuleBase::NEG_IMAG_UNIT, l); for (int ig=0; igovlp_pswfcjlq, it, ipswfc, - GlobalV::NQX, GlobalV::DQ, gk[ig].norm() * GlobalC::ucell.tpiba ); + GlobalV::NQX, GlobalV::DQ, gk[ig].norm() * this->p_ucell->tpiba ); } /* NSPIN == 4 */ if(GlobalV::NSPIN == 4) { - if(GlobalC::ucell.atoms[it].ncpp.has_so) + if(this->p_ucell->atoms[it].ncpp.has_so) { Soc soc; soc.rot_ylm(l + 1); - const double j = GlobalC::ucell.atoms[it].ncpp.jchi[ipswfc]; + const double j = this->p_ucell->atoms[it].ncpp.jchi[ipswfc]; /* NOT NONCOLINEAR CASE, rotation matrix become identity */ if (!(GlobalV::DOMAG||GlobalV::DOMAG_Z)) { @@ -192,7 +215,7 @@ psi::Psi>* psi_initializer_atomic::cal_psig(int ik) if(fabs(cg_coeffs[is]) > 1e-8) { /* GET COMPLEX SPHERICAL HARMONIC FUNCTION */ - const int ind = GlobalC::ppcell.lmaxkb + soc.sph_ind(l,j,m,is); // ind can be l+m, l+m+1, l+m-1 + const int ind = this->p_pspot_nl->lmaxkb + soc.sph_ind(l,j,m,is); // ind can be l+m, l+m+1, l+m-1 ModuleBase::GlobalFunc::ZEROS(aux, npw); for(int n1 = 0; n1 < 2*l+1; n1++) { @@ -244,11 +267,11 @@ psi::Psi>* psi_initializer_atomic::cal_psig(int ik) else { /* L != 0, scan pswfcs that have the same L and satisfy J(pswfc) = L - 0.5 */ - for(int jpsiwfc = 0; jpsiwfc < GlobalC::ucell.atoms[it].ncpp.nchi; jpsiwfc++) + for(int jpsiwfc = 0; jpsiwfc < this->p_ucell->atoms[it].ncpp.nchi; jpsiwfc++) { if( - (GlobalC::ucell.atoms[it].ncpp.lchi[jpsiwfc] == l) - &&(fabs(GlobalC::ucell.atoms[it].ncpp.jchi[jpsiwfc] - l + 0.5) < 1e-4)) + (this->p_ucell->atoms[it].ncpp.lchi[jpsiwfc] == l) + &&(fabs(this->p_ucell->atoms[it].ncpp.jchi[jpsiwfc] - l + 0.5) < 1e-4)) { ipswfc_noncolin_soc = jpsiwfc; break; @@ -260,7 +283,7 @@ psi::Psi>* psi_initializer_atomic::cal_psig(int ik) chiaux[ig] = l * ModuleBase::PolyInt::Polynomial_Interpolation( this->ovlp_pswfcjlq, it, ipswfc_noncolin_soc, - GlobalV::NQX, GlobalV::DQ, gk[ig].norm() * GlobalC::ucell.tpiba); + GlobalV::NQX, GlobalV::DQ, gk[ig].norm() * this->p_ucell->tpiba); chiaux[ig] += ovlp_pswfcjlg[ig] * (l + 1.0) ; chiaux[ig] *= 1/(2.0*l+1.0); } @@ -268,13 +291,17 @@ psi::Psi>* psi_initializer_atomic::cal_psig(int ik) /* ROTATE ACCORDING TO NONCOLINEAR */ double alpha, gamma; std::complex fup, fdw; - alpha = GlobalC::ucell.atoms[it].angle1[ia]; - gamma = -1 * GlobalC::ucell.atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; + alpha = this->p_ucell->atoms[it].angle1[ia]; + gamma = -1 * this->p_ucell->atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; for(int m = 0; m < 2*l+1; m++) { const int lm = l*l +m; - if(index+2*l+1>GlobalC::ucell.natomwfc) ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()","error: too many wfcs"); + if(index+2*l+1 > this->p_ucell->natomwfc) + { + std::cout<<__FILE__<<__LINE__<<" "<p_ucell->natomwfc<>* psi_initializer_atomic::cal_psig(int ik) {//atomic_wfc_nc double alpha, gamman; std::complex fup, fdown; - //alpha = GlobalC::ucell.magnet.angle1_[it]; - //gamman = -GlobalC::ucell.magnet.angle2_[it] + 0.5*ModuleBase::PI; - alpha = GlobalC::ucell.atoms[it].angle1[ia]; - gamman = -1 * GlobalC::ucell.atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; + //alpha = this->p_ucell->magnet.angle1_[it]; + //gamman = -this->p_ucell->magnet.angle2_[it] + 0.5*ModuleBase::PI; + alpha = this->p_ucell->atoms[it].angle1[ia]; + gamman = -1 * this->p_ucell->atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; for(int m = 0; m < 2*l+1; m++) { const int lm = l*l +m; - if(index+2*l+1>GlobalC::ucell.natomwfc) ModuleBase::WARNING_QUIT("GlobalC::wf.atomic_wfc()","error: too many wfcs"); + if(index+2*l+1 > this->p_ucell->natomwfc) + { + std::cout<<__FILE__<<__LINE__<<" "<p_ucell->natomwfc<>* psi_initializer_atomic::cal_psig(int ik) /* complement the rest of bands if there are */ if(this->get_nbands_complem() > 0) { - this->random_t(this->psig->get_pointer(), index, this->psig->get_nbands(), ik, this->pw_wfc); + this->random_t(this->psig->get_pointer(), index, this->psig->get_nbands(), ik); } - ModuleBase::timer::tick("psi_initializer_atomic", "initialize"); -#ifdef PSI_INITIALIZER_TEST - this->write_psig(); -#endif + ModuleBase::timer::tick("psi_initializer_atomic", "cal_psig"); return this->psig; } \ No newline at end of file diff --git a/source/module_psi/psi_initializer_atomic.h b/source/module_psi/psi_initializer_atomic.h index dfb9b8b5d6..82733d0667 100644 --- a/source/module_psi/psi_initializer_atomic.h +++ b/source/module_psi/psi_initializer_atomic.h @@ -9,11 +9,22 @@ Psi (planewave based wavefunction) initializer: atomic class psi_initializer_atomic : public psi_initializer { public: - - /// @brief constructor of psi initializer (atomic) - /// @param sf_in Structure factor interface, link ESolver - /// @param pw_wfc_in ModulePW::PW_Basis_K interface, link ESolver - psi_initializer_atomic(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in); + #ifdef __MPI + /// @brief parameterized constructor of psi initializer (with MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param p_parakpts_in interface, link with Parallel_Kpoints GlobalC::Pkpoints + /// @param random_seed_in random seed + psi_initializer_atomic(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in = 1); + #else + /// @brief parameterized constructor of psi initializer (without MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param random_seed_in random seed + psi_initializer_atomic(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in = 1); + #endif /// @brief default destructor ~psi_initializer_atomic(); @@ -23,25 +34,32 @@ class psi_initializer_atomic : public psi_initializer /// @return initialized planewave wavefunction (psi::Psi>*) psi::Psi>* cal_psig(int ik) override; + /// @brief initialize only once, for atomic, it should be, create ovlp_pswfcjlq, calculate ovlp_pswfcjlq + /// @param p_pspot_nl_in (for atomic) interfaces to pseudopot_cell_vnl object, in GlobalC now + /// @attention if one variable is necessary for all methods, initialize it in constructor, not here. + void initialize_only_once(pseudopot_cell_vnl* p_pspot_nl_in) override; // setters + /* I leave this function here for deprecation of UnitCell in the future */ /// @brief setter of pseudpotential filenames /// @param pseudopot_files pseudpotential filenames organized in an array - void set_pseudopot_files(std::string* pseudopot_files); + //void set_pseudopot_files(std::string* pseudopot_files); // I wont write a function to set ovlp_pswfcjlq, it is totally useless + /// @brief allocate memory for ovlp_pswfcjlq and initialize all elements to 0 + void create_ovlp_Xjlq() override; /// @brief specialized normalization of wfc function /// @param n_rgrid number of grid points in realspace /// @param pswfc pseudo wavefunction in pseudopotential files /// @param rgrid realspace grid points, r1, r2, ... - void normalize_pswfc(int n_rgrid, double* pswfc, double* rgrid); + void normalize_pswfc(int n_rgrid, double* pswfc, double* rgrid) override; /// @brief simple unitary phase factor /// @param arg the argument of the phase factor /// @param mode +1 for real part, -1 for imaginary part, 0 for the whole /// @return the phase factor - std::complex phase_factor(double arg, int mode = 0); + std::complex phase_factor(double arg, int mode = 0) override; /// @brief calculate the overlap between pseudo atomic wavefunctions and planewave basis - void cal_ovlp_pswfcjlq(); + void cal_ovlp_pswfcjlq() override; // historically left functions // getters diff --git a/source/module_psi/psi_initializer_atomic_random.cpp b/source/module_psi/psi_initializer_atomic_random.cpp index 03137c02ab..27ed823668 100644 --- a/source/module_psi/psi_initializer_atomic_random.cpp +++ b/source/module_psi/psi_initializer_atomic_random.cpp @@ -1,6 +1,12 @@ #include "psi_initializer_atomic_random.h" -psi_initializer_atomic_random::psi_initializer_atomic_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in) : psi_initializer_atomic(sf_in, pw_wfc_in) +#ifdef __MPI +psi_initializer_atomic_random::psi_initializer_atomic_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in) + : psi_initializer_atomic(sf_in, pw_wfc_in, p_ucell_in, p_parakpts_in, random_seed_in) +#else +psi_initializer_atomic_random::psi_initializer_atomic_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in) + : psi_initializer_atomic(sf_in, pw_wfc_in, p_ucell_in, random_seed_in) +#endif { this->set_method("atomic+random"); this->set_random_mix(0.05); @@ -8,6 +14,11 @@ psi_initializer_atomic_random::psi_initializer_atomic_random(Structure_Factor* s psi_initializer_atomic_random::~psi_initializer_atomic_random() {} +void psi_initializer_atomic_random::initialize_only_once(pseudopot_cell_vnl* p_pspot_nl_in) +{ + psi_initializer_atomic::initialize_only_once(p_pspot_nl_in); +} + psi::Psi>* psi_initializer_atomic_random::cal_psig(int ik) { double rm = this->get_random_mix(); @@ -15,7 +26,7 @@ psi::Psi>* psi_initializer_atomic_random::cal_psig(int ik) this->psig = psi_initializer_atomic::cal_psig(ik); psi::Psi> psi_random(1, this->psig->get_nbands(), this->psig->get_nbasis(), nullptr); psi_random.fix_k(0); - this->random_t(psi_random.get_pointer(), 0, psi_random.get_nbands(), ik, this->pw_wfc); + this->random_t(psi_random.get_pointer(), 0, psi_random.get_nbands(), ik); for(int iband = 0; iband < this->psig->get_nbands(); iband++) { for(int ibasis = 0; ibasis < this->psig->get_nbasis(); ibasis++) @@ -23,8 +34,5 @@ psi::Psi>* psi_initializer_atomic_random::cal_psig(int ik) (*(this->psig))(iband, ibasis) = (1-rm)*(*(this->psig))(iband, ibasis) + rm*psi_random(iband, ibasis); } } -#ifdef PSI_INITIALIZER_TEST - this->write_psig(); -#endif return this->psig; } \ No newline at end of file diff --git a/source/module_psi/psi_initializer_atomic_random.h b/source/module_psi/psi_initializer_atomic_random.h index e09261d273..99156e246e 100644 --- a/source/module_psi/psi_initializer_atomic_random.h +++ b/source/module_psi/psi_initializer_atomic_random.h @@ -1,22 +1,38 @@ #ifndef PSI_INITIALIZER_ATOMIC_RANDOM_H #define PSI_INITIALIZER_ATOMIC_RANDOM_H #include "psi_initializer_atomic.h" - +#include "module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h" +#include "module_cell/parallel_kpoints.h" /* Psi (planewave based wavefunction) initializer: atomic+random */ class psi_initializer_atomic_random : public psi_initializer_atomic { public: - - /// @brief constructor of psi initializer (atomic+random) - /// @param sf_in Structure factor interface, link ESolver - /// @param pw_wfc_in ModulePW::PW_Basis_K interface, link ESolver - psi_initializer_atomic_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in); + #ifdef __MPI + /// @brief parameterized constructor of psi initializer (with MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param p_parakpts_in interface, link with Parallel_Kpoints GlobalC::Pkpoints + /// @param random_seed_in random seed + psi_initializer_atomic_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in = 1); + #else + /// @brief parameterized constructor of psi initializer (without MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param random_seed_in random seed + psi_initializer_atomic_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in = 1); + #endif /// @brief default destructor ~psi_initializer_atomic_random(); + /// @brief for variables can be only initialized for once. + /// @param p_pspot_nl_in (for atomic) interfaces to pseudopot_cell_vnl object, in GlobalC now + /// @attention if one variable is necessary for all methods, initialize it in constructor, not here. + void initialize_only_once(pseudopot_cell_vnl* p_pspot_nl_in = nullptr) override; // methods /// @brief calculate and output planewave wavefunction /// @param ik kpoint index diff --git a/source/module_psi/psi_initializer_nao.cpp b/source/module_psi/psi_initializer_nao.cpp index 20b072c4a8..e09f53a1d5 100644 --- a/source/module_psi/psi_initializer_nao.cpp +++ b/source/module_psi/psi_initializer_nao.cpp @@ -10,35 +10,20 @@ #include "module_base/timer.h" // three global variables definition #include "module_base/global_variable.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" // parallel communication #ifdef __MPI #include "module_base/parallel_common.h" #endif -psi_initializer_nao::psi_initializer_nao(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in) : psi_initializer(sf_in, pw_wfc_in) +#ifdef __MPI +psi_initializer_nao::psi_initializer_nao(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in) + : psi_initializer(sf_in, pw_wfc_in, p_ucell_in, p_parakpts_in, random_seed_in) +#else +psi_initializer_nao::psi_initializer_nao(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in) + : psi_initializer(sf_in, pw_wfc_in, p_ucell_in, random_seed_in) +#endif { this->set_method("nao"); - // find correct dimension for ovlp_flzjlq - int dim1 = GlobalC::ucell.ntype; - int dim2 = 0; // dim2 should be the maximum number of zeta for each atomtype - for (int it = 0; it < GlobalC::ucell.ntype; it++) - { - int nzeta = 0; - for (int l = 0; l < GlobalC::ucell.atoms[it].nwl+1; l++) - { - nzeta += GlobalC::ucell.atoms[it].l_nchi[l]; - } - dim2 = (nzeta > dim2) ? nzeta : dim2; - } - if (dim2 == 0) - { - ModuleBase::WARNING_QUIT("psi_initializer_nao::psi_initializer_nao", "there is not ANY numerical atomic orbital read in present system, quit."); - } - int dim3 = GlobalV::NQX; - // allocate memory for ovlp_flzjlq - this->ovlp_flzjlq.create(dim1, dim2, dim3); - this->ovlp_flzjlq.zero_out(); } psi_initializer_nao::~psi_initializer_nao() {} @@ -54,7 +39,7 @@ void psi_initializer_nao::set_orbital_files(std::string* orbital_files) if(GlobalV::MY_RANK == 0) { #endif - for (int itype = 0; itype < GlobalC::ucell.ntype; itype++) + for (int itype = 0; itype < this->p_ucell->ntype; itype++) { this->orbital_files.push_back(orbital_files[itype]); } @@ -62,13 +47,37 @@ void psi_initializer_nao::set_orbital_files(std::string* orbital_files) } else { - this->orbital_files.resize(GlobalC::ucell.ntype); + this->orbital_files.resize(this->p_ucell->ntype); } - Parallel_Common::bcast_string(this->orbital_files.data(), GlobalC::ucell.ntype); + Parallel_Common::bcast_string(this->orbital_files.data(), this->p_ucell->ntype); #endif ModuleBase::timer::tick("psi_initializer_nao", "set_orbital_files"); } +void psi_initializer_nao::create_ovlp_Xjlq() +{ + // find correct dimension for ovlp_flzjlq + int dim1 = this->p_ucell->ntype; + int dim2 = 0; // dim2 should be the maximum number of zeta for each atomtype + for (int it = 0; it < this->p_ucell->ntype; it++) + { + int nzeta = 0; + for (int l = 0; l < this->p_ucell->atoms[it].nwl+1; l++) + { + nzeta += this->p_ucell->atoms[it].l_nchi[l]; + } + dim2 = (nzeta > dim2) ? nzeta : dim2; + } + if (dim2 == 0) + { + ModuleBase::WARNING_QUIT("psi_initializer_nao::psi_initializer_nao", "there is not ANY numerical atomic orbital read in present system, quit."); + } + int dim3 = GlobalV::NQX; + // allocate memory for ovlp_flzjlq + this->ovlp_flzjlq.create(dim1, dim2, dim3); + this->ovlp_flzjlq.zero_out(); +} + void psi_initializer_nao::read_orbital_files() { ModuleBase::timer::tick("psi_initializer_nao", "read_orbital_files"); @@ -76,13 +85,13 @@ void psi_initializer_nao::read_orbital_files() if(GlobalV::MY_RANK==0) { #endif - for(int it = 0; it < GlobalC::ucell.ntype; it++) + for(int it = 0; it < this->p_ucell->ntype; it++) { // number of chi per atomtype int nchi = 0; - for(int l = 0; l <= GlobalC::ucell.atoms[it].nwl; l++) + for(int l = 0; l <= this->p_ucell->atoms[it].nwl; l++) { - nchi += GlobalC::ucell.atoms[it].l_nchi[l]; + nchi += this->p_ucell->atoms[it].l_nchi[l]; } std::vector n_rgrid_it; @@ -106,16 +115,16 @@ void psi_initializer_nao::read_orbital_files() int ichi_overall = 0; // check nwl and nchi for each rank - for(int l = 0; l <= GlobalC::ucell.atoms[it].nwl; l++) + for(int l = 0; l <= this->p_ucell->atoms[it].nwl; l++) { - for(int ichi = 0; ichi < GlobalC::ucell.atoms[it].l_nchi[l]; ichi++) + for(int ichi = 0; ichi < this->p_ucell->atoms[it].l_nchi[l]; ichi++) { int n_rgrid_ichi; std::vector rgrid_ichi; std::vector flz_ichi; GlobalV::ofs_running<<"-------------------------------------- "<p_ucell->atoms[it].label<p_ucell->atoms[it].label<create_ovlp_Xjlq(); this->read_orbital_files(); - + //this->cal_ovlp_flzjlq(); //because GlobalV::NQX will change during vcrelax, so it should be called in both init and init_after_vc + ModuleBase::timer::tick("psi_initializer_nao", "initialize_only_once"); +} + +void psi_initializer_nao::cal_ovlp_flzjlq() +{ ModuleBase::timer::tick("psi_initializer_nao", "cal_ovlp_flzjlq"); + //this->read_orbital_files(); this->ovlp_flzjlq.zero_out(); - for(int it=0; itp_ucell->ntype; it++) { int ic=0; - for(int l=0; lp_ucell->atoms[it].nwl+1; l++) { - for(int izeta=0; izetap_ucell->atoms[it].l_nchi[l]; izeta++) { double* ovlp_flzjlq_q = new double[GlobalV::NQX]; double* qgrid = new double[GlobalV::NQX]; @@ -266,13 +282,6 @@ void psi_initializer_nao::cal_ovlp_flzjlq() { qgrid[iq] = iq*GlobalV::DQ; } - /* - this->sbt.direct(l, - GlobalC::ucell.atoms[it].n_rgrid[ic], - GlobalC::ucell.atoms[it].rgrid[ic], - GlobalC::ucell.atoms[it].flz[ic], - GlobalV::NQX, qgrid, ovlp_flzjlq_q); - */ this->sbt.direct(l, this->n_rgrid[it][ic], this->rgrid[it][ic].data(), @@ -290,19 +299,19 @@ void psi_initializer_nao::cal_ovlp_flzjlq() if(GlobalV::MY_RANK==0) { - for(int it = 0; it < GlobalC::ucell.ntype; it++) + for(int it = 0; it < this->p_ucell->ntype; it++) { std::stringstream ss; - ss<p_ucell->atoms[it].label<< "/LOCAL_G.dat"; std::ofstream ofs(ss.str().c_str()); for(int iq = 0; iq < GlobalV::NQX; iq++) { int ic=0; double energy_q = pow((double)iq*GlobalV::DQ, 2); - ofs<p_ucell->tpiba2; + for(int l = 0; lp_ucell->atoms[it].nwl + 1; l++) { - for(int N=0; Np_ucell->atoms[it].l_nchi[l]; N++) { ofs<<" "<>* psi_initializer_nao::cal_psig(int ik) assert(ik>=0); this->psig->fix_k(ik); const int npw = this->pw_wfc->npwk[ik]; - const int total_lm = ( GlobalC::ucell.lmax + 1) * ( GlobalC::ucell.lmax + 1); + const int total_lm = ( this->p_ucell->lmax + 1) * ( this->p_ucell->lmax + 1); ModuleBase::matrix ylm(total_lm, npw); std::complex *aux = new std::complex[npw]; double *chiaux = nullptr; @@ -338,40 +347,39 @@ psi::Psi>* psi_initializer_nao::cal_psig(int ik) //int index = 0; double *ovlp_flzjlg = new double[npw]; int ibasis=0; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < this->p_ucell->ntype; it++) { - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + for (int ia = 0; ia < this->p_ucell->atoms[it].na; ia++) { /* HERE LOOP OVER ALL ATOMIS */ std::complex* sk = this->sf->get_sk(ik, it, ia, this->pw_wfc); int ic = 0; // ic is a flatten index of chi, therefore it is defined here. - for(int L = 0; L < GlobalC::ucell.atoms[it].nwl+1; L++) + for(int L = 0; L < this->p_ucell->atoms[it].nwl+1; L++) { std::complex lphase = pow(ModuleBase::NEG_IMAG_UNIT, L); //mohan 2010-04-19 - for(int N=0; N < GlobalC::ucell.atoms[it].l_nchi[L]; N++) + for(int N=0; N < this->p_ucell->atoms[it].l_nchi[L]; N++) { /* HERE LOOP OVER ALL NAOS */ for(int ig=0; igovlp_flzjlq, - it, ic, GlobalV::NQX, GlobalV::DQ, gk[ig].norm() * GlobalC::ucell.tpiba ); + it, ic, GlobalV::NQX, GlobalV::DQ, gk[ig].norm() * this->p_ucell->tpiba ); } /* FOR EVERY NAO IN EACH ATOM */ if(GlobalV::NSPIN==4) { /* FOR EACH SPIN CHANNEL */ - for(int is_N = 0; is_N < 2; is_N++) + for(int is_N = 0; is_N < 2; is_N++) // rotate base //for(int is_N = 0; is_N < 1; is_N++) { if(L==0 && is_N==1) continue; - if(GlobalC::ucell.atoms[it].ncpp.has_so) + if(this->p_ucell->atoms[it].ncpp.has_so) { const double j = std::abs(double(L+is_N) - 0.5); if (!(GlobalV::DOMAG||GlobalV::DOMAG_Z)) {//atomic_wfc_so for(int m=0; m<2*L+1; m++) { - std::cout<<"ibasis: "<>* psi_initializer_nao::cal_psig(int ik) {//Average the two functions chiaux[ig] = L * ModuleBase::PolyInt::Polynomial_Interpolation(this->ovlp_flzjlq, - it, ic, GlobalV::NQX, GlobalV::DQ, gk[ig].norm() * GlobalC::ucell.tpiba ); + it, ic, GlobalV::NQX, GlobalV::DQ, gk[ig].norm() * this->p_ucell->tpiba ); chiaux[ig] += ovlp_flzjlg[ig] * (L+1.0) ; chiaux[ig] *= 1/(2.0*L+1.0); } } - alpha = GlobalC::ucell.atoms[it].angle1[ia]; - gamma = -1 * GlobalC::ucell.atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; + alpha = this->p_ucell->atoms[it].angle1[ia]; + gamma = -1 * this->p_ucell->atoms[it].angle2[ia] + 0.5 * ModuleBase::PI; for(int m = 0;m<2*L+1;m++) { const int lm = L*L +m; @@ -440,15 +448,15 @@ psi::Psi>* psi_initializer_nao::cal_psig(int ik) } ibasis += 2*L +1; } // end else INPUT.starting_spin_angle || !GlobalV::DOMAG - } // end if GlobalC::ucell.atoms[it].has_so + } // end if this->p_ucell->atoms[it].has_so else {//atomic_wfc_nc double alpha, gamman; std::complex fup, fdown; - //alpha = GlobalC::ucell.magnet.angle1_[it]; - //gamman = -GlobalC::ucell.magnet.angle2_[it] + 0.5*ModuleBase::PI; - alpha = GlobalC::ucell.atoms[it].angle1[ia]; - gamman = -GlobalC::ucell.atoms[it].angle2[ia] + 0.5*ModuleBase::PI; + //alpha = this->p_ucell->magnet.angle1_[it]; + //gamman = -this->p_ucell->magnet.angle2_[it] + 0.5*ModuleBase::PI; + alpha = this->p_ucell->atoms[it].angle1[ia]; + gamman = -this->p_ucell->atoms[it].angle2[ia] + 0.5*ModuleBase::PI; for(int m = 0; m < 2*L+1; m++) { const int lm = L*L +m; @@ -477,7 +485,7 @@ psi::Psi>* psi_initializer_nao::cal_psig(int ik) ibasis++; } // end m ibasis += 2*L+1; - } // end else GlobalC::ucell.atoms[it].has_so + } // end else this->p_ucell->atoms[it].has_so } // end for is_N } // end if GlobalV::NONCOLIN else{//LSDA and nomagnet case @@ -506,12 +514,9 @@ psi::Psi>* psi_initializer_nao::cal_psig(int ik) /* complement the rest of bands if there are */ if(this->get_nbands_complem() > 0) { - this->random_t(this->psig->get_pointer(), ibasis, this->psig->get_nbands(), ik, this->pw_wfc); + this->random_t(this->psig->get_pointer(), ibasis, this->psig->get_nbands(), ik); } ModuleBase::timer::tick("psi_initializer_nao", "initialize"); -#ifdef PSI_INITIALIZER_TEST - this->write_psig(); -#endif return this->psig; } \ No newline at end of file diff --git a/source/module_psi/psi_initializer_nao.h b/source/module_psi/psi_initializer_nao.h index 91b3c72eac..062ac83c88 100644 --- a/source/module_psi/psi_initializer_nao.h +++ b/source/module_psi/psi_initializer_nao.h @@ -9,10 +9,22 @@ Psi (planewave based wavefunction) initializer: numerical atomic orbital method class psi_initializer_nao : public psi_initializer { public: - /// @brief constructor of psi initializer (nao) - /// @param sf_in Structure factor interface, link ESolver - /// @param pw_wfc_in ModulePW::PW_Basis_K interface, link ESolver - psi_initializer_nao(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in); + #ifdef __MPI + /// @brief parameterized constructor of psi initializer (with MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param p_parakpts_in interface, link with Parallel_Kpoints GlobalC::Pkpoints + /// @param random_seed_in random seed + psi_initializer_nao(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in = 1); + #else + /// @brief parameterized constructor of psi initializer (without MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param random_seed_in random seed + psi_initializer_nao(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in = 1); + #endif /// @brief default destructor ~psi_initializer_nao(); @@ -23,18 +35,25 @@ class psi_initializer_nao : public psi_initializer /// @return initialized planewave wavefunction (psi::Psi>*) psi::Psi>* cal_psig(int ik) override; + /// @brief initialize only once, for nao, it should be, read numerical orbitals, create ovlp_Xjlq(, calculate ovlp_flzjlq) + /// @param p_pspot_nl_in (for atomic) interfaces to pseudopot_cell_vnl object, in GlobalC now + /// @attention if one variable is necessary for all methods, initialize it in constructor, not here. + void initialize_only_once(pseudopot_cell_vnl* p_pspot_nl_in = nullptr) override; // setters /// @brief setter of numerical orbital files /// @param orbital_files array storing numerical orbital files - void set_orbital_files(std::string* orbital_files); + void set_orbital_files(std::string* orbital_files) override; // I wont write a function to set ovlp_flzjlq, it is totally useless - + + /// @brief allocate memory for ovlp_flzjlq and initialize all elements to 0 + /// @attention warning! p_ucell must be set in advance! + void create_ovlp_Xjlq() override; /// @brief before refactor and reorganization of UnitCell class, it is temporary to write this function here. /// In future version, it will be moved into UnitCell class. void read_orbital_files(); /// @brief calculate overlap integral between f_{l\\zeta} the radial numerical orbital and spherical Bessel function - void cal_ovlp_flzjlq(); + void cal_ovlp_flzjlq() override; // getters diff --git a/source/module_psi/psi_initializer_nao_random.cpp b/source/module_psi/psi_initializer_nao_random.cpp index 45cddd27fd..0e094d096e 100644 --- a/source/module_psi/psi_initializer_nao_random.cpp +++ b/source/module_psi/psi_initializer_nao_random.cpp @@ -1,6 +1,12 @@ #include "psi_initializer_nao_random.h" -psi_initializer_nao_random::psi_initializer_nao_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in) : psi_initializer_nao(sf_in, pw_wfc_in) +#ifdef __MPI +psi_initializer_nao_random::psi_initializer_nao_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in) + : psi_initializer_nao(sf_in, pw_wfc_in, p_ucell_in, p_parakpts_in, random_seed_in) +#else +psi_initializer_nao_random::psi_initializer_nao_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in) + : psi_initializer_nao(sf_in, pw_wfc_in, p_ucell_in, random_seed_in) +#endif { this->set_random_mix(0.05); this->set_method("nao+random"); @@ -8,6 +14,11 @@ psi_initializer_nao_random::psi_initializer_nao_random(Structure_Factor* sf_in, psi_initializer_nao_random::~psi_initializer_nao_random() {} +void psi_initializer_nao_random::initialize_only_once(pseudopot_cell_vnl* p_pspot_nl_in) +{ + psi_initializer_nao::initialize_only_once(); +} + psi::Psi>* psi_initializer_nao_random::cal_psig(int ik) { double rm = this->get_random_mix(); @@ -15,7 +26,7 @@ psi::Psi>* psi_initializer_nao_random::cal_psig(int ik) this->psig = psi_initializer_nao::cal_psig(ik); psi::Psi> psi_random(1, this->psig->get_nbands(), this->psig->get_nbasis(), nullptr); psi_random.fix_k(0); - this->random_t(psi_random.get_pointer(), 0, psi_random.get_nbands(), ik, this->pw_wfc); + this->random_t(psi_random.get_pointer(), 0, psi_random.get_nbands(), ik); for(int iband = 0; iband < this->psig->get_nbands(); iband++) { for(int ibasis = 0; ibasis < this->psig->get_nbasis(); ibasis++) @@ -23,8 +34,5 @@ psi::Psi>* psi_initializer_nao_random::cal_psig(int ik) (*(this->psig))(iband, ibasis) = (1-rm)*(*(this->psig))(iband, ibasis) + rm*psi_random(iband, ibasis); } } -#ifdef PSI_INITIALIZER_TEST - this->write_psig(); -#endif return this->psig; } \ No newline at end of file diff --git a/source/module_psi/psi_initializer_nao_random.h b/source/module_psi/psi_initializer_nao_random.h index 82774f3040..418b29983a 100644 --- a/source/module_psi/psi_initializer_nao_random.h +++ b/source/module_psi/psi_initializer_nao_random.h @@ -1,21 +1,34 @@ #ifndef PSI_INITIALIZER_NAO_RANDOM_H #define PSI_INITIALIZER_NAO_RANDOM_H #include "psi_initializer_nao.h" - +#include "module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h" +#include "module_cell/parallel_kpoints.h" /* Psi (planewave based wavefunction) initializer: numerical atomic orbital + random method */ class psi_initializer_nao_random : public psi_initializer_nao { public: - - /// @brief constructor of psi initializer (nao+random) - /// @param sf_in Structure factor interface, link ESolver - /// @param pw_wfc_in ModulePW::PW_Basis_K interface, link ESolver - psi_initializer_nao_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in); + #ifdef __MPI + /// @brief parameterized constructor of psi initializer (with MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param p_parakpts_in interface, link with Parallel_Kpoints GlobalC::Pkpoints + /// @param random_seed_in random seed + psi_initializer_nao_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in = 1); + #else + /// @brief parameterized constructor of psi initializer (without MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param random_seed_in random seed + psi_initializer_nao_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in = 1); + #endif /// @brief default destructor ~psi_initializer_nao_random(); + void initialize_only_once(pseudopot_cell_vnl* p_pspot_nl_in = nullptr) override; /// @brief calculate and output planewave wavefunction /// @param ik kpoint index /// @return initialized planewave wavefunction (psi::Psi>*) diff --git a/source/module_psi/psi_initializer_random.cpp b/source/module_psi/psi_initializer_random.cpp index 06309ddfc6..c06c33c2bd 100644 --- a/source/module_psi/psi_initializer_random.cpp +++ b/source/module_psi/psi_initializer_random.cpp @@ -7,23 +7,26 @@ // basic functions support #include "module_base/timer.h" -psi_initializer_random::psi_initializer_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in) : psi_initializer(sf_in, pw_wfc_in) +#ifdef __MPI +psi_initializer_random::psi_initializer_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in) + : psi_initializer(sf_in, pw_wfc_in, p_ucell_in, p_parakpts_in, random_seed_in) +#else +psi_initializer_random::psi_initializer_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in) + : psi_initializer(sf_in, pw_wfc_in, p_ucell_in, random_seed_in) +#endif { this->set_method("random"); } - psi_initializer_random::~psi_initializer_random() {} - void psi_initializer_random::random(std::complex* psi, const int iw_start, const int iw_end, - const int ik, - const ModulePW::PW_Basis_K* wfc_basis) + const int ik) { ModuleBase::timer::tick("psi_initializer_random", "random"); - this->random_t(psi, iw_start, iw_end, ik, wfc_basis); + this->random_t(psi, iw_start, iw_end, ik); ModuleBase::timer::tick("psi_initializer_random", "random"); } @@ -32,12 +35,9 @@ psi::Psi>* psi_initializer_random::cal_psig(int ik) ModuleBase::timer::tick("psi_initializer_random", "initialize"); //this->print_status(psi); this->psig->fix_k(ik); - this->random(this->psig->get_pointer(), 0, this->psig->get_nbands(), ik, this->pw_wfc); + this->random(this->psig->get_pointer(), 0, this->psig->get_nbands(), ik); // we still need to diagonalize the obtained psi from hsolver::DiagoIterAssist::diagH_subspace // will do it in HSolver function... ModuleBase::timer::tick("psi_initializer_random", "initialize"); -#ifdef PSI_INITIALIZER_TEST - this->write_psig(); -#endif return this->psig; -} \ No newline at end of file +} diff --git a/source/module_psi/psi_initializer_random.h b/source/module_psi/psi_initializer_random.h index 1b870d7627..9bcc0d922e 100644 --- a/source/module_psi/psi_initializer_random.h +++ b/source/module_psi/psi_initializer_random.h @@ -2,16 +2,29 @@ #define PSI_INITIALIZER_RANDOM_H #include "psi_initializer.h" +#include "module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h" /* Psi (planewave based wavefunction) initializer: random method */ class psi_initializer_random : public psi_initializer { public: - /// @brief constructor of psi initializer (random) - /// @param sf_in Structure factor interface, link ESolver - /// @param pw_wfc_in ModulePW::PW_Basis_K interface, link ESolver - psi_initializer_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in); + #ifdef __MPI + /// @brief parameterized constructor of psi initializer (with MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param p_parakpts_in interface, link with Parallel_Kpoints GlobalC::Pkpoints + /// @param random_seed_in random seed + psi_initializer_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, Parallel_Kpoints* p_parakpts_in, int random_seed_in = 1); + #else + /// @brief parameterized constructor of psi initializer (without MPI support) + /// @param sf_in interface, link with Structure_Factor ESolver_FP::sf + /// @param pw_wfc_in interface, link with ModulePW::PW_Basis_K* ESolver_FP::pw_wfc + /// @param p_ucell_in interface, link with UnitCell GlobalC::ucell + /// @param random_seed_in random seed + psi_initializer_random(Structure_Factor* sf_in, ModulePW::PW_Basis_K* pw_wfc_in, UnitCell* p_ucell_in, int random_seed_in = 1); + #endif /// @brief default destructor ~psi_initializer_random(); @@ -22,15 +35,17 @@ class psi_initializer_random : public psi_initializer /// @param iw_start the starting band index to fill random value /// @param iw_end the end band index /// @param ik kpoint index - /// @param wfc_basis ModulePW::PW_Basis_K interface, will be automatically linked void random(std::complex* psi, const int iw_start, const int iw_end, - const int ik, - const ModulePW::PW_Basis_K* wfc_basis); + const int ik) override; /// @brief calculate and output planewave wavefunction /// @param ik kpoint index /// @return initialized planewave wavefunction (psi::Psi>*) psi::Psi>* cal_psig(int ik) override; + /// @brief for variables can be only initialized for once. + /// @param p_pspot_nl_in (for atomic) interfaces to pseudopot_cell_vnl object, in GlobalC now + /// @attention if one variable is necessary for all methods, initialize it in constructor, not here. + void initialize_only_once(pseudopot_cell_vnl* p_pspot_nl_in = nullptr) override {}; }; #endif \ No newline at end of file diff --git a/source/module_psi/test/CMakeLists.txt b/source/module_psi/test/CMakeLists.txt index f9696f5c4c..10bc51be83 100644 --- a/source/module_psi/test/CMakeLists.txt +++ b/source/module_psi/test/CMakeLists.txt @@ -8,8 +8,15 @@ AddTest( if(ENABLE_LCAO) AddTest( - TARGET psi_initializer_integrate_test + TARGET psi_initializer_unit_test + LIBS ${math_libs} base device psi psi_initializer planewave SOURCES - psi_initializer_integrated_test.cpp + psi_initializer_unit_test.cpp + ../../module_hamilt_pw/hamilt_pwdft/soc.cpp + ../../module_cell/atom_spec.cpp + ../../module_cell/parallel_kpoints.cpp + mock_unitcell.cpp ) endif() + +install(DIRECTORY support DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) \ No newline at end of file diff --git a/source/module_psi/test/mock_unitcell.cpp b/source/module_psi/test/mock_unitcell.cpp new file mode 100644 index 0000000000..7015f38055 --- /dev/null +++ b/source/module_psi/test/mock_unitcell.cpp @@ -0,0 +1,107 @@ +#include "module_cell/unitcell.h" +/* + README: + This file supports idea like "I dont need any functions of UnitCell, I want to avoid using UnitCell functions because there is GLobalC, which will bring + endless compile troubles like undefined behavior" +*/ +void UnitCell::cal_ux() {} +bool UnitCell::judge_parallel(double a[3], ModuleBase::Vector3 b) {return true;} +void UnitCell::set_iat2iwt(const int& npol_in) {} +UnitCell::UnitCell() +{ + if (GlobalV::test_unitcell) + ModuleBase::TITLE("unitcell", "Constructor"); + Coordinate = "Direct"; + latName = "none"; + lat0 = 0.0; + lat0_angstrom = 0.0; + + ntype = 0; + nat = 0; + namax = 0; + nwmax = 0; + + iat2it = nullptr; + iat2ia = nullptr; + iwt2iat = nullptr; + iwt2iw = nullptr; + + itia2iat.create(1, 1); + lc = new int[3]; + + latvec = ModuleBase::Matrix3(); + latvec_supercell = ModuleBase::Matrix3(); + G = ModuleBase::Matrix3(); + GT = ModuleBase::Matrix3(); + GGT = ModuleBase::Matrix3(); + invGGT = ModuleBase::Matrix3(); + + tpiba = 0.0; + tpiba2 = 0.0; + omega = 0.0; + + atom_label = new std::string[1]; + atom_mass = nullptr; + pseudo_fn = new std::string[1]; + pseudo_type = new std::string[1]; + orbital_fn = new std::string[1]; + + set_atom_flag = false; +} +UnitCell::~UnitCell() +{ + delete[] atom_label; + delete[] atom_mass; + delete[] pseudo_fn; + delete[] pseudo_type; + delete[] orbital_fn; + delete[] iat2it; + delete[] iat2ia; + delete[] iwt2iat; + delete[] iwt2iw; + delete[] lc; + if(set_atom_flag) + { + delete[] atoms; + } +} +void UnitCell::print_cell(std::ofstream& ofs) const {} +void UnitCell::print_cell_xyz(const std::string& fn) const {} +void UnitCell::print_cell_cif(const std::string& fn) const {} +int UnitCell::read_atom_species(std::ifstream &ifa, std::ofstream &ofs_running) {return 0;} +void UnitCell::read_cell_pseudopots(const std::string &pp_dir, std::ofstream &log) {} +bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_running, std::ofstream &ofs_warning) {return true;} +void UnitCell::update_pos_tau(const double* pos) {} +void UnitCell::update_pos_taud(double* posd_in) {} +void UnitCell::update_pos_taud(const ModuleBase::Vector3* posd_in) {} +void UnitCell::update_vel(const ModuleBase::Vector3* vel_in) {} +void UnitCell::periodic_boundary_adjustment() {} +void UnitCell::bcast_atoms_tau() {} +bool UnitCell::judge_big_cell() const {return true;} +void UnitCell::update_stress(ModuleBase::matrix &scs) {} +void UnitCell::update_force(ModuleBase::matrix &fcs) {} +#ifdef __MPI +void UnitCell::bcast_unitcell() {} +void UnitCell::bcast_unitcell2() {} +#endif +void UnitCell::set_iat2itia(void) {} +void UnitCell::setup_cell(const std::string &fn, std::ofstream &log) {} +void UnitCell::read_orb_file(int it, std::string &orb_file, std::ofstream &ofs_running, Atom *atom) {} +void UnitCell::read_pseudo(std::ofstream &ofs) {} +int UnitCell::find_type(const std::string &label) {return 0;} +void UnitCell::print_tau(void) const {} +void UnitCell::print_stru_file(const std::string &fn, const int &type, const int &level)const {} +void UnitCell::check_dtau(void) {} +void UnitCell::setup_cell_after_vc(std::ofstream &log) {} +void UnitCell::remake_cell() {} +void UnitCell::cal_nwfc(std::ofstream &log) {} +void UnitCell::cal_meshx() {} +void UnitCell::cal_natomwfc(std::ofstream &log) {} +void UnitCell::print_unitcell_pseudo(const std::string &fn) {} +bool UnitCell::check_tau(void)const {return true;} +bool UnitCell::if_atoms_can_move()const {return true;} +bool UnitCell::if_cell_can_change()const {return true;} +void UnitCell::setup(const std::string &latname_in, const int &ntype_in, const int &lmaxmax_in, const bool &init_vel_in, const std::string &fixed_axes_in) {} +void UnitCell::check_structure(double factor) {} +void UnitCell::cal_nelec(double& nelec) {} +void UnitCell::compare_atom_labels(std::string label1, std::string label2) {} \ No newline at end of file diff --git a/source/module_psi/test/psi_initializer_integrated_test.cpp b/source/module_psi/test/psi_initializer_integrated_test.cpp deleted file mode 100644 index 09ea124f79..0000000000 --- a/source/module_psi/test/psi_initializer_integrated_test.cpp +++ /dev/null @@ -1,230 +0,0 @@ -#include -#include -#include -#include - -#include - -/* -==================== -psi initializer test -==================== -- Instruction - For psi initialization needs several modules, e.g., hsolver (for diagonalizing band representation Hamiltonian matrix), - hamilt (for imposing operators onto psi), ElecState, Potential and Charge (for initializing operators of hamilt), PW_Basis - and PW_Basis_K (for initiailizing Structure_Factor class and psi::Psi), Structure_Factor (for Fourier transform on sum of atom-wise properties), - UnitCell (for calculating Structure_Factor), therefore the tests are designed for random, atomic and nao for two kind of purpose: - random: directly comparing between psi values with the old implemented version, - one for old version, one for this new version. - atomic and nao: new SphericalBesselTransformer implemented, therefore for nearly all spherical Bessel functions, accurancies are improved, - only prepare reference value here, to ensure if other modules changed, the initialization will not be impacted too much and lead to bad initial guess. -- Tested functions - - cal_psig - - including implementation in atomic, random, nao cases: - - psi_initializer_random::cal_psig - - psi_initializer_atomic::cal_psig - - psi_initializer_nao::cal_psig -- need input file - - INPUT: random_old, random_new, atomic_new, nao_new - - KPT: KPT - - STRU: STRU - - UPF: Si_NCSR_ONCVPSP_v0.5_dojo.upf - - ORB: Si_gga_8au_60Ry_2s2p1d.orb -*/ - -class IntegratedInitializerTest : public ::testing::Test { -protected: - int error = 0; - std::ifstream ifs_new; - std::ifstream ifs_old; - int n_match_max; - int n_match; - - std::regex pattern; - virtual void SetUp() { - std::cout << "IntegratedInitializerTest SetUp" << std::endl; - this->n_match = 0; - this->n_match_max = 10; - this->pattern = std::regex("\\((-?[0-9]+\\.[0-9]*),(-?[0-9]+\\.[0-9]*)\\)"); - this->error = std::system("cp ../../../../source/module_psi/test/support/KPT ./KPT"); - this->error = std::system("cp ../../../../source/module_psi/test/support/STRU ./STRU"); - this->error = std::system("cp ../../../../source/module_psi/test/support/Si_NCSR_ONCVPSP_v0.5_dojo.upf ./Si_NCSR_ONCVPSP_v0.5_dojo.upf"); - this->error = std::system("cp ../../../../source/module_psi/test/support/Si_gga_8au_60Ry_2s2p1d.orb ./Si_gga_8au_60Ry_2s2p1d.orb"); - } - - virtual void TearDown() { - std::cout << "IntegratedInitializerTest TearDown" << std::endl; - this->error = std::system("rm ./psig_0_kpt.out"); - // finalize - this->error = std::system("rm ./Si_NCSR_ONCVPSP_v0.5_dojo.upf"); - this->error = std::system("rm ./Si_gga_8au_60Ry_2s2p1d.orb"); - this->error = std::system("rm ./KPT"); - this->error = std::system("rm ./STRU"); - this->error = std::system("rm -rf ./OUT*"); - this->error = std::system("rm ./time.json"); - } -}; - -TEST_F(IntegratedInitializerTest, CalPsiGRandom) { - - int error = 0; - std::string line_new; - std::smatch match_new; - // compare wavefunctions - // init_wfc = "random" - error = std::system("cp ../../../../source/module_psi/test/support/random_new ./INPUT"); - error = std::system("../../../abacus"); - error = std::system("rm ./INPUT"); - // compare wavefunctions with reference value - std::vector reals_ref = { - 0.0028527200, 0.0021489486, 0.0070330784, -0.0052732148, 0.0008623089, - -0.0194856855, 0.1288743928, -0.6618769676, 0.0586698846, 0.0226938509 - }; - std::vector imags_ref = { - 0.0011941671, -0.0004010668, -0.0095634837, -0.0013454707, -0.0133724755, - -0.0342739193, 0.0324688840, 0.5175270812, -0.1867769238, 0.0663531911 - }; - double* reals_read = new double[reals_ref.size()]; - double* imags_read = new double[imags_ref.size()]; - this->ifs_new.open("./psig_0_kpt.out"); - - bool find_psi = false; - while (ifs_new >> line_new) { - if (line_new == "random") { - std::cout << "Have find correct psi" << std::endl; - find_psi = true; - break; - }; - } - ASSERT_TRUE(find_psi); - while (ifs_new >> line_new) { - if (regex_search(line_new, match_new, this->pattern)) { //check if the line contains a complex number tuple - double real = stod(match_new[1]); //extract the real part of the complex number - double imag = stod(match_new[2]); //extract the imaginary part of the complex number - reals_read[this->n_match] = real; - imags_read[this->n_match] = imag; - this->n_match++; - if (this->n_match == this->n_match_max) break; - } - } - this->ifs_new.close(); - for (int i = 0; i < reals_ref.size(); i++) { - EXPECT_NEAR(reals_read[i], reals_ref[i], 1e-5); - EXPECT_NEAR(imags_read[i], imags_ref[i], 1e-5); - } - delete[] reals_read; - delete[] imags_read; -} - -TEST_F(IntegratedInitializerTest, CalPsiGAtomic) { - - int error = 0; - std::string line_new; - std::smatch match_new; - // compare wavefunctions - // init_wfc = "atomic" - error = std::system("cp ../../../../source/module_psi/test/support/atomic_new ./INPUT"); - error = std::system("../../../abacus"); - error = std::system("rm ./INPUT"); - // compare wavefunctions with reference value - std::vector reals_ref = { - -0.0000060004, 0.0000170688, 0.0000537252, -0.0006497418, -0.0050094394, - -0.0129273778, 0.1125669791, 0.9193169677, 0.1125669791, -0.0129273778 - }; - std::vector imags_ref = { - 0.0, 0.0, 0.0, 0.0, 0.0, - 0.0, 0.0, -0.0, -0.0, -0.0 - }; - double* reals_read = new double[reals_ref.size()]; - double* imags_read = new double[imags_ref.size()]; - this->ifs_new.open("./psig_0_kpt.out"); - - bool find_psi = false; - while (ifs_new >> line_new) { - if (line_new == "atomic") { - std::cout << "Have find correct psi" << std::endl; - find_psi = true; - break; - }; - } - ASSERT_TRUE(find_psi); - while (ifs_new >> line_new) { - if (regex_search(line_new, match_new, this->pattern)) { //check if the line contains a complex number tuple - double real = stod(match_new[1]); //extract the real part of the complex number - double imag = stod(match_new[2]); //extract the imaginary part of the complex number - reals_read[this->n_match] = real; - imags_read[this->n_match] = imag; - this->n_match++; - if (this->n_match == this->n_match_max) break; - } - } - this->ifs_new.close(); - for (int i = 0; i < reals_ref.size(); i++) { - EXPECT_NEAR(reals_read[i], reals_ref[i], 1e-5); - EXPECT_NEAR(imags_read[i], imags_ref[i], 1e-5); - } - delete[] reals_read; - delete[] imags_read; -} - -TEST_F(IntegratedInitializerTest, CalPsiGNao) { - - int error = 0; - std::string line_new; - std::smatch match_new; - // compare wavefunctions - // init_wfc = "nao" - error = std::system("cp ../../../../source/module_psi/test/support/nao_new ./INPUT"); - error = std::system("../../../abacus"); - error = std::system("rm ./INPUT"); - std::vector reals_ref = { - -1.2553e-06, -7.1728e-05, 0.000112783, -0.000599359, -0.00603285, - -0.0164273, 0.14776, 1.16658, 0.14776, -0.0164273 - }; - std::vector imags_ref = { - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 - }; - double* reals_read = new double[reals_ref.size()]; - double* imags_read = new double[imags_ref.size()]; - // compare wavefunctions with reference value - this->ifs_new.open("./psig_0_kpt.out"); - - bool find_psi = false; - while (ifs_new >> line_new) { - if (line_new == "nao") { - std::cout << "Have find correct psi" << std::endl; - find_psi = true; - break; - }; - } - ASSERT_TRUE(find_psi); - - while (ifs_new >> line_new) { - if (regex_search(line_new, match_new, this->pattern)) { //check if the line contains a complex number tuple - double real = stod(match_new[1]); //extract the real part of the complex number - double imag = stod(match_new[2]); //extract the imaginary part of the complex number - reals_read[this->n_match] = real; - imags_read[this->n_match] = imag; - n_match++; - if (this->n_match == this->n_match_max) break; - } - } - this->ifs_new.close(); - for (int i = 0; i < reals_ref.size(); i++) { - EXPECT_NEAR(reals_read[i], reals_ref[i], 1e-5); - EXPECT_NEAR(imags_read[i], imags_ref[i], 1e-5); - } - delete[] reals_read; - delete[] imags_read; -} - -#include -int main(int argc, char **argv) -{ - char cwd[1024]; - std::string return_val = getcwd(cwd, sizeof(cwd)); - std::cout << "present directory: " << cwd << std::endl; - testing::InitGoogleTest(&argc, argv); - int result = RUN_ALL_TESTS(); - return result; -} \ No newline at end of file diff --git a/source/module_psi/test/psi_initializer_unit_test.cpp b/source/module_psi/test/psi_initializer_unit_test.cpp new file mode 100644 index 0000000000..2a2000648b --- /dev/null +++ b/source/module_psi/test/psi_initializer_unit_test.cpp @@ -0,0 +1,633 @@ +#include +#include "../psi_initializer.h" +#include "../psi_initializer_random.h" +#include "../psi_initializer_atomic.h" +#include "../psi_initializer_nao.h" +#include "../psi_initializer_atomic_random.h" +#include "../psi_initializer_nao_random.h" + +/* +========================= +psi initializer unit test +========================= +- Tested functions: + - psi_initializer_random::psi_initializer_random + - constructor of psi_initializer_random + - psi_initializer_atomic::psi_initializer_atomic + - constructor of psi_initializer_atomic + - psi_initializer_atomic_random::psi_initializer_atomic_random + - constructor of psi_initializer_atomic_random + - psi_initializer_nao::psi_initializer_nao + - constructor of psi_initializer_nao + - psi_initializer_nao_random::psi_initializer_nao_random + - constructor of psi_initializer_nao_random + - psi_initializer_random::allocate + - allocate wavefunctions with random-specific method + - psi_initializer_atomic::allocate + - allocate wavefunctions with atomic-specific method + - psi_initializer_atomic_random::allocate + - allocate wavefunctions with atomic-specific method + - psi_initializer_nao::allocate + - allocate wavefunctions with nao-specific method + - psi_initializer_nao_random::allocate + - allocate wavefunctions with nao-specific method + - psi_initializer_random::cal_psig + - calculate wavefunction initial guess (before diagonalization) by randomly generating numbers + - psi_initializer_atomic::cal_psig + - calculate wavefunction initial guess (before diagonalization) with atomic pseudo wavefunctions + - nspin = 4 case + - nspin = 4 with has_so case + - psi_initializer_atomic_random::cal_psig + - calculate wavefunction initial guess (before diagonalization) with atomic pseudo wavefunctions and random numbers + - psi_initializer_nao::cal_psig + - calculate wavefunction initial guess (before diagonalization) with numerical atomic orbital wavefunctions + - psi_initializer_nao_random::cal_psig + - calculate wavefunction initial guess (before diagonalization) with numerical atomic orbital wavefunctions and random numbers +*/ + +// here are many empty functions because we have included many headers in file to be tested +// but it does not mean all functions can only be defined for only once in those corresponding files +// so here we define them again to avoid undefined reference error +Atom_pseudo::Atom_pseudo() {} +Atom_pseudo::~Atom_pseudo() {} +#ifdef __MPI +void Atom_pseudo::bcast_atom_pseudo() {} +#endif +pseudo::pseudo() {} +pseudo::~pseudo() {} + +pseudopot_cell_vnl::pseudopot_cell_vnl() {} +pseudopot_cell_vnl::~pseudopot_cell_vnl() {} +pseudopot_cell_vl::pseudopot_cell_vl() {} +pseudopot_cell_vl::~pseudopot_cell_vl() {} +Magnetism::Magnetism() {} +Magnetism::~Magnetism() {} +void output::printM3(std::ofstream &ofs, const std::string &description, const ModuleBase::Matrix3 &m) {} +#ifdef __LCAO +ORB_gaunt_table::ORB_gaunt_table() {} +ORB_gaunt_table::~ORB_gaunt_table() {} +InfoNonlocal::InfoNonlocal() {} +InfoNonlocal::~InfoNonlocal() {} +#endif +Structure_Factor::Structure_Factor() {} +Structure_Factor::~Structure_Factor() {} +void Structure_Factor::setup_structure_factor(UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis) {} +std::complex* Structure_Factor::get_sk(int ik, int it, int ia, ModulePW::PW_Basis_K const*wfc_basis) const +{ + int npw = wfc_basis->npwk[ik]; + std::complex *sk = new std::complex[npw]; + for(int ipw = 0; ipw < npw; ++ipw) sk[ipw] = std::complex(0.0, 0.0); + return sk; +} + +class PsiIntializerUnitTest : public ::testing::Test { + public: + Structure_Factor* p_sf = nullptr; + ModulePW::PW_Basis_K* p_pw_wfc = nullptr; + UnitCell* p_ucell = nullptr; + pseudopot_cell_vnl* p_pspot_vnl = nullptr; + #ifdef __MPI + Parallel_Kpoints* p_parakpts = nullptr; + #endif + int random_seed = 1; + + psi_initializer* psi_init; + private: + protected: + void SetUp() override + { + // allocate + this->p_sf = new Structure_Factor(); + this->p_pw_wfc = new ModulePW::PW_Basis_K(); + this->p_ucell = new UnitCell(); + this->p_pspot_vnl = new pseudopot_cell_vnl(); + #ifdef __MPI + this->p_parakpts = new Parallel_Kpoints(); + #endif + // mock + GlobalV::NBANDS = 1; + GlobalV::NSPIN = 1; + GlobalV::global_orbital_dir = "./support/"; + GlobalV::global_pseudo_dir = "./support/"; + GlobalV::NPOL = 1; + GlobalV::CALCULATION = "scf"; + GlobalV::init_wfc = "random"; + GlobalV::KS_SOLVER = "cg"; + GlobalV::DOMAG = false; + GlobalV::DOMAG_Z = false; + // lattice + this->p_ucell->a1 = {10.0, 0.0, 0.0}; + this->p_ucell->a2 = {0.0, 10.0, 0.0}; + this->p_ucell->a3 = {0.0, 0.0, 10.0}; + this->p_ucell->lat0 = 1.0; + this->p_ucell->latvec.e11 = 10.0; this->p_ucell->latvec.e12 = 0.0; this->p_ucell->latvec.e13 = 0.0; + this->p_ucell->latvec.e21 = 0.0; this->p_ucell->latvec.e22 = 10.0; this->p_ucell->latvec.e23 = 0.0; + this->p_ucell->latvec.e31 = 0.0; this->p_ucell->latvec.e32 = 0.0; this->p_ucell->latvec.e33 = 10.0; + this->p_ucell->GT = this->p_ucell->latvec.Inverse(); + this->p_ucell->G = this->p_ucell->GT.Transpose(); + this->p_ucell->GGT = this->p_ucell->G * this->p_ucell->GT; + this->p_ucell->tpiba = 2.0 * M_PI / this->p_ucell->lat0; + this->p_ucell->tpiba2 = this->p_ucell->tpiba * this->p_ucell->tpiba; + // atom + if(this->p_ucell->atom_label != nullptr) delete[] this->p_ucell->atom_label; + this->p_ucell->atom_label = new std::string[1]; + this->p_ucell->atom_label[0] = "Si"; + // atom properties + this->p_ucell->nat = 1; + this->p_ucell->ntype = 1; + this->p_ucell->atoms = new Atom[1]; + this->p_ucell->atoms[0].label = "Si"; + this->p_ucell->atoms[0].mass = 28.0855; + this->p_ucell->atoms[0].na = 1; + this->p_ucell->atoms[0].angle1 = new double[1]; + this->p_ucell->atoms[0].angle1[0] = 0.0; + this->p_ucell->atoms[0].angle2 = new double[1]; + this->p_ucell->atoms[0].angle2[0] = 0.0; + // atom position + this->p_ucell->atoms[0].tau[0] = {0.0, 0.0, 0.0}; + this->p_ucell->atoms[0].taud[0] = {0.25, 0.25, 0.25}; + this->p_ucell->atoms[0].mbl[0] = {0, 0, 0}; + // atom pseudopotential + if(this->p_ucell->pseudo_fn != nullptr) delete[] this->p_ucell->pseudo_fn; + this->p_ucell->pseudo_fn = new std::string[1]; + this->p_ucell->pseudo_fn[0] = "Si_NCSR_ONCVPSP_v0.5_dojo.upf"; + this->p_ucell->natomwfc = 4; + this->p_ucell->atoms[0].ncpp.nchi = 2; + this->p_ucell->atoms[0].ncpp.mesh = 10; + this->p_ucell->atoms[0].ncpp.msh = 10; + this->p_ucell->atoms[0].ncpp.lmax = 2; + //if(this->p_ucell->atoms[0].ncpp.rab != nullptr) delete[] this->p_ucell->atoms[0].ncpp.rab; + this->p_ucell->atoms[0].ncpp.rab = new double[10]; + for(int i = 0; i < 10; ++i) this->p_ucell->atoms[0].ncpp.rab[i] = 0.01; + //if(this->p_ucell->atoms[0].ncpp.r != nullptr) delete[] this->p_ucell->atoms[0].ncpp.r; + this->p_ucell->atoms[0].ncpp.r = new double[10]; + for(int i = 0; i < 10; ++i) this->p_ucell->atoms[0].ncpp.r[i] = 0.01*i; + this->p_ucell->atoms[0].ncpp.chi.create(2, 10); + for(int i = 0; i < 2; ++i) for(int j = 0; j < 10; ++j) this->p_ucell->atoms[0].ncpp.chi(i, j) = 0.01; + //if(this->p_ucell->atoms[0].ncpp.lchi != nullptr) delete[] this->p_ucell->atoms[0].ncpp.lchi; + this->p_ucell->atoms[0].ncpp.lchi = new int[2]; + this->p_ucell->atoms[0].ncpp.lchi[0] = 0; + this->p_ucell->atoms[0].ncpp.lchi[1] = 1; + this->p_ucell->lmax_ppwf = 1; + this->p_ucell->atoms[0].ncpp.oc = new double[2]; + this->p_ucell->atoms[0].ncpp.oc[0] = 1.0; + this->p_ucell->atoms[0].ncpp.oc[1] = 1.0; + + this->p_ucell->atoms[0].ncpp.has_so = false; + this->p_ucell->atoms[0].ncpp.jchi = new double[2]; + this->p_ucell->atoms[0].ncpp.jchi[0] = 0.5; + this->p_ucell->atoms[0].ncpp.jchi[1] = 1.5; + // atom numerical orbital + this->p_ucell->lmax = 2; + if(this->p_ucell->orbital_fn != nullptr) delete[] this->p_ucell->orbital_fn; + this->p_ucell->orbital_fn = new std::string[1]; + this->p_ucell->orbital_fn[0] = "Si_gga_8au_60Ry_2s2p1d.orb"; + this->p_ucell->atoms[0].nwl = 2; + delete[] this->p_ucell->atoms[0].l_nchi; + this->p_ucell->atoms[0].l_nchi = new int[3]; + this->p_ucell->atoms[0].l_nchi[0] = 2; + this->p_ucell->atoms[0].l_nchi[1] = 2; + this->p_ucell->atoms[0].l_nchi[2] = 1; + + + // can support function PW_Basis::getfftixy2is + this->p_pw_wfc->nks = 1; + this->p_pw_wfc->npwk_max = 1; + if(this->p_pw_wfc->npwk != nullptr) delete[] this->p_pw_wfc->npwk; + this->p_pw_wfc->npwk = new int[1]; + this->p_pw_wfc->npwk[0] = 1; + this->p_pw_wfc->fftnxy = 1; + this->p_pw_wfc->fftnz = 1; + this->p_pw_wfc->nst = 1; + if(this->p_pw_wfc->is2fftixy != nullptr) delete[] this->p_pw_wfc->is2fftixy; + this->p_pw_wfc->is2fftixy = new int[1]; + this->p_pw_wfc->is2fftixy[0] = 0; + if(this->p_pw_wfc->fftixy2ip != nullptr) delete[] this->p_pw_wfc->fftixy2ip; + this->p_pw_wfc->fftixy2ip = new int[1]; + this->p_pw_wfc->fftixy2ip[0] = 0; + if(this->p_pw_wfc->igl2isz_k != nullptr) delete[] this->p_pw_wfc->igl2isz_k; + this->p_pw_wfc->igl2isz_k = new int[1]; + this->p_pw_wfc->igl2isz_k[0] = 0; + if(this->p_pw_wfc->gcar != nullptr) delete[] this->p_pw_wfc->gcar; + this->p_pw_wfc->gcar = new ModuleBase::Vector3[1]; + this->p_pw_wfc->gcar[0] = {0.0, 0.0, 0.0}; + if(this->p_pw_wfc->igl2isz_k != nullptr) delete[] this->p_pw_wfc->igl2isz_k; + this->p_pw_wfc->igl2isz_k = new int[1]; + this->p_pw_wfc->igl2isz_k[0] = 0; + if(this->p_pw_wfc->gk2 != nullptr) delete[] this->p_pw_wfc->gk2; + this->p_pw_wfc->gk2 = new double[1]; + this->p_pw_wfc->gk2[0] = 0.0; + this->p_pw_wfc->latvec.e11 = this->p_ucell->latvec.e11; this->p_pw_wfc->latvec.e12 = this->p_ucell->latvec.e12; this->p_pw_wfc->latvec.e13 = this->p_ucell->latvec.e13; + this->p_pw_wfc->latvec.e21 = this->p_ucell->latvec.e21; this->p_pw_wfc->latvec.e22 = this->p_ucell->latvec.e22; this->p_pw_wfc->latvec.e23 = this->p_ucell->latvec.e23; + this->p_pw_wfc->latvec.e31 = this->p_ucell->latvec.e31; this->p_pw_wfc->latvec.e32 = this->p_ucell->latvec.e32; this->p_pw_wfc->latvec.e33 = this->p_ucell->latvec.e33; + this->p_pw_wfc->G = this->p_ucell->G; + this->p_pw_wfc->GT = this->p_ucell->GT; + this->p_pw_wfc->GGT = this->p_ucell->GGT; + this->p_pw_wfc->lat0 = this->p_ucell->lat0; + this->p_pw_wfc->tpiba = 2.0 * M_PI / this->p_ucell->lat0; + this->p_pw_wfc->tpiba2 = this->p_pw_wfc->tpiba * this->p_pw_wfc->tpiba; + if(this->p_pw_wfc->kvec_c != nullptr) delete[] this->p_pw_wfc->kvec_c; + this->p_pw_wfc->kvec_c = new ModuleBase::Vector3[1]; + this->p_pw_wfc->kvec_c[0] = {0.0, 0.0, 0.0}; + if(this->p_pw_wfc->kvec_d != nullptr) delete[] this->p_pw_wfc->kvec_d; + this->p_pw_wfc->kvec_d = new ModuleBase::Vector3[1]; + this->p_pw_wfc->kvec_d[0] = {0.0, 0.0, 0.0}; + + #ifdef __MPI + if(this->p_parakpts->startk_pool != nullptr) delete[] this->p_parakpts->startk_pool; + this->p_parakpts->startk_pool = new int[1]; + this->p_parakpts->startk_pool[0] = 0; + #endif + + } + void TearDown() override + { + delete this->p_sf; + delete this->p_pw_wfc; + delete this->p_ucell; + delete this->p_pspot_vnl; + #ifdef __MPI + delete this->p_parakpts; + #endif + } +}; + +TEST_F(PsiIntializerUnitTest, ConstructorRandom) { + #ifdef __MPI + this->psi_init = new psi_initializer_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + EXPECT_EQ("random", this->psi_init->get_method()); + EXPECT_EQ(this->p_sf, this->psi_init->get_interface_sf()); + EXPECT_EQ(this->p_pw_wfc, this->psi_init->get_interface_pw_wfc()); + EXPECT_EQ(this->p_ucell, this->psi_init->get_interface_ucell()); + EXPECT_EQ(this->random_seed, this->psi_init->get_random_seed()); + #ifdef __MPI + EXPECT_EQ(this->p_parakpts, this->psi_init->get_interface_parakpts()); + #endif +} + +TEST_F(PsiIntializerUnitTest, ConstructorAtomic) { + #ifdef __MPI + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + EXPECT_EQ("atomic", this->psi_init->get_method()); + EXPECT_EQ(this->p_sf, this->psi_init->get_interface_sf()); + EXPECT_EQ(this->p_pw_wfc, this->psi_init->get_interface_pw_wfc()); + EXPECT_EQ(this->p_ucell, this->psi_init->get_interface_ucell()); + EXPECT_EQ(this->random_seed, this->psi_init->get_random_seed()); + #ifdef __MPI + EXPECT_EQ(this->p_parakpts, this->psi_init->get_interface_parakpts()); + #endif +} + +TEST_F(PsiIntializerUnitTest, ConstructorAtomicRandom) { + #ifdef __MPI + this->psi_init = new psi_initializer_atomic_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_atomic_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + EXPECT_EQ("atomic+random", this->psi_init->get_method()); + EXPECT_EQ(this->p_sf, this->psi_init->get_interface_sf()); + EXPECT_EQ(this->p_pw_wfc, this->psi_init->get_interface_pw_wfc()); + EXPECT_EQ(this->p_ucell, this->psi_init->get_interface_ucell()); + EXPECT_EQ(this->random_seed, this->psi_init->get_random_seed()); + #ifdef __MPI + EXPECT_EQ(this->p_parakpts, this->psi_init->get_interface_parakpts()); + #endif +} + +TEST_F(PsiIntializerUnitTest, ConstructorNao) { + #ifdef __MPI + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + EXPECT_EQ("nao", this->psi_init->get_method()); + EXPECT_EQ(this->p_sf, this->psi_init->get_interface_sf()); + EXPECT_EQ(this->p_pw_wfc, this->psi_init->get_interface_pw_wfc()); + EXPECT_EQ(this->p_ucell, this->psi_init->get_interface_ucell()); + EXPECT_EQ(this->random_seed, this->psi_init->get_random_seed()); + #ifdef __MPI + EXPECT_EQ(this->p_parakpts, this->psi_init->get_interface_parakpts()); + #endif +} + +TEST_F(PsiIntializerUnitTest, ConstructorNaoRandom) { + #ifdef __MPI + this->psi_init = new psi_initializer_nao_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_nao_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + EXPECT_EQ("nao+random", this->psi_init->get_method()); + EXPECT_EQ(this->p_sf, this->psi_init->get_interface_sf()); + EXPECT_EQ(this->p_pw_wfc, this->psi_init->get_interface_pw_wfc()); + EXPECT_EQ(this->p_ucell, this->psi_init->get_interface_ucell()); + EXPECT_EQ(this->random_seed, this->psi_init->get_random_seed()); + #ifdef __MPI + EXPECT_EQ(this->p_parakpts, this->psi_init->get_interface_parakpts()); + #endif +} + +TEST_F(PsiIntializerUnitTest, AllocateRandom) { + GlobalV::init_wfc = "random"; + #ifdef __MPI + this->psi_init = new psi_initializer_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->initialize_only_once(); + psi::Psi>* psi = this->psi_init->allocate(); + EXPECT_EQ(0, this->psi_init->get_nbands_complem()); + EXPECT_EQ(1, psi->get_nk()); + EXPECT_EQ(1, psi->get_nbands()); + EXPECT_EQ(1, psi->get_nbasis()); + EXPECT_EQ(1, this->psi_init->psig->get_nk()); + EXPECT_EQ(1, this->psi_init->psig->get_nbands()); + EXPECT_EQ(1, this->psi_init->psig->get_nbasis()); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, AllocateAtomic) { + GlobalV::init_wfc = "atomic"; + #ifdef __MPI + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->initialize_only_once(this->p_pspot_vnl); + psi::Psi>* psi = this->psi_init->allocate(); + EXPECT_EQ(0, this->psi_init->get_nbands_complem()); + EXPECT_EQ(1, psi->get_nk()); + EXPECT_EQ(1, psi->get_nbands()); + EXPECT_EQ(1, psi->get_nbasis()); + EXPECT_EQ(1, this->psi_init->psig->get_nk()); + EXPECT_EQ(4, this->psi_init->psig->get_nbands()); + EXPECT_EQ(1, this->psi_init->psig->get_nbasis()); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, AllocateAtomicRandom) { + GlobalV::init_wfc = "atomic+random"; + #ifdef __MPI + this->psi_init = new psi_initializer_atomic_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_atomic_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->initialize_only_once(this->p_pspot_vnl); + psi::Psi>* psi = this->psi_init->allocate(); + EXPECT_EQ(0, this->psi_init->get_nbands_complem()); + EXPECT_EQ(1, psi->get_nk()); + EXPECT_EQ(1, psi->get_nbands()); + EXPECT_EQ(1, psi->get_nbasis()); + EXPECT_EQ(1, this->psi_init->psig->get_nk()); + EXPECT_EQ(4, this->psi_init->psig->get_nbands()); + EXPECT_EQ(1, this->psi_init->psig->get_nbasis()); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, AllocateNao) { + GlobalV::init_wfc = "nao"; + #ifdef __MPI + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->set_orbital_files(this->p_ucell->orbital_fn); + this->psi_init->initialize_only_once(); + psi::Psi>* psi = this->psi_init->allocate(); + EXPECT_EQ(0, this->psi_init->get_nbands_complem()); + EXPECT_EQ(1, psi->get_nk()); + EXPECT_EQ(1, psi->get_nbands()); + EXPECT_EQ(1, psi->get_nbasis()); + EXPECT_EQ(1, this->psi_init->psig->get_nk()); + EXPECT_EQ(13, this->psi_init->psig->get_nbands()); + EXPECT_EQ(1, this->psi_init->psig->get_nbasis()); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, AllocateNaoRandom) { + GlobalV::init_wfc = "nao+random"; + #ifdef __MPI + this->psi_init = new psi_initializer_nao_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_nao_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->set_orbital_files(this->p_ucell->orbital_fn); + this->psi_init->initialize_only_once(); + psi::Psi>* psi = this->psi_init->allocate(); + EXPECT_EQ(0, this->psi_init->get_nbands_complem()); + EXPECT_EQ(1, psi->get_nk()); + EXPECT_EQ(1, psi->get_nbands()); + EXPECT_EQ(1, psi->get_nbasis()); + EXPECT_EQ(1, this->psi_init->psig->get_nk()); + EXPECT_EQ(13, this->psi_init->psig->get_nbands()); + EXPECT_EQ(1, this->psi_init->psig->get_nbasis()); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, CalPsigRandom) { + GlobalV::init_wfc = "random"; + #ifdef __MPI + this->psi_init = new psi_initializer_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, CalPsigAtomic) { + GlobalV::init_wfc = "atomic"; + #ifdef __MPI + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->initialize_only_once(this->p_pspot_vnl); + this->psi_init->cal_ovlp_pswfcjlq(); + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, CalPsigAtomicSoc) { + GlobalV::init_wfc = "atomic"; + GlobalV::NSPIN = 4; + GlobalV::NPOL = 2; + this->p_ucell->atoms[0].ncpp.has_so = false; + this->p_ucell->natomwfc *= 2; + #ifdef __MPI + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->initialize_only_once(this->p_pspot_vnl); + this->psi_init->cal_ovlp_pswfcjlq(); + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + GlobalV::NSPIN = 1; + GlobalV::NPOL = 1; + this->p_ucell->atoms[0].ncpp.has_so = false; + this->p_ucell->natomwfc /= 2; + delete psi; +} + +TEST_F(PsiIntializerUnitTest, CalPsigAtomicSocHasSo) { + GlobalV::init_wfc = "atomic"; + GlobalV::NSPIN = 4; + GlobalV::NPOL = 2; + this->p_ucell->atoms[0].ncpp.has_so = true; + this->p_ucell->natomwfc *= 2; + #ifdef __MPI + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_atomic(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->initialize_only_once(this->p_pspot_vnl); + this->psi_init->cal_ovlp_pswfcjlq(); + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + GlobalV::NSPIN = 1; + GlobalV::NPOL = 1; + this->p_ucell->atoms[0].ncpp.has_so = false; + this->p_ucell->natomwfc /= 2; + delete psi; +} + +TEST_F(PsiIntializerUnitTest, CalPsigAtomicRandom) { + GlobalV::init_wfc = "atomic+random"; + #ifdef __MPI + this->psi_init = new psi_initializer_atomic_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_atomic_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->initialize_only_once(this->p_pspot_vnl); + this->psi_init->cal_ovlp_pswfcjlq(); + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, CalPsigNao) { + GlobalV::init_wfc = "nao"; + #ifdef __MPI + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->set_orbital_files(this->p_ucell->orbital_fn); + this->psi_init->initialize_only_once(); + this->psi_init->cal_ovlp_flzjlq(); + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, CalPsigNaoRandom) { + GlobalV::init_wfc = "nao+random"; + #ifdef __MPI + this->psi_init = new psi_initializer_nao_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_nao_random(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->set_orbital_files(this->p_ucell->orbital_fn); + this->psi_init->initialize_only_once(); + this->psi_init->cal_ovlp_flzjlq(); + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + delete psi; +} +/* +TEST_F(PsiIntializerUnitTest, CalPsigNaoSoc) { + GlobalV::init_wfc = "nao"; + GlobalV::NSPIN = 4; + GlobalV::NPOL = 2; + this->p_ucell->atoms[0].ncpp.has_so = false; + GlobalV::DOMAG = false; + GlobalV::DOMAG_Z = false; + #ifdef __MPI + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->set_orbital_files(this->p_ucell->orbital_fn); + this->psi_init->initialize_only_once(); + this->psi_init->cal_ovlp_flzjlq(); + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, CalPsigNaoSocHasSo) { + GlobalV::init_wfc = "nao"; + GlobalV::NSPIN = 4; + GlobalV::NPOL = 2; + this->p_ucell->atoms[0].ncpp.has_so = true; + GlobalV::DOMAG = false; + GlobalV::DOMAG_Z = false; + #ifdef __MPI + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->set_orbital_files(this->p_ucell->orbital_fn); + this->psi_init->initialize_only_once(); + this->psi_init->cal_ovlp_flzjlq(); + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + delete psi; +} + +TEST_F(PsiIntializerUnitTest, CalPsigNaoSocHasSoDOMAG) { + GlobalV::init_wfc = "nao"; + GlobalV::NSPIN = 4; + GlobalV::NPOL = 2; + this->p_ucell->atoms[0].ncpp.has_so = true; + GlobalV::DOMAG = true; + GlobalV::DOMAG_Z = false; + #ifdef __MPI + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->p_parakpts, this->random_seed); + #else + this->psi_init = new psi_initializer_nao(this->p_sf, this->p_pw_wfc, this->p_ucell, this->random_seed); + #endif + this->psi_init->set_orbital_files(this->p_ucell->orbital_fn); + this->psi_init->initialize_only_once(); + this->psi_init->cal_ovlp_flzjlq(); + psi::Psi>* psi = this->psi_init->allocate(); + psi::Psi>* psig = this->psi_init->cal_psig(0); + EXPECT_NEAR(0, psi->operator()(0,0,0).real(), 1e-12); + delete psi; +} +*/ +int main(int argc, char** argv) +{ + +#ifdef __MPI + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &GlobalV::NPROC); + MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK); +#endif + + testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + +#ifdef __MPI + MPI_Finalize(); +#endif + + return result; +} \ No newline at end of file diff --git a/tests/integrate/102_PW_PINT_RKS/INPUT b/tests/integrate/102_PW_PINT_RKS/INPUT new file mode 100644 index 0000000000..285ec0a5fa --- /dev/null +++ b/tests/integrate/102_PW_PINT_RKS/INPUT @@ -0,0 +1,22 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../../tests/PP_ORB +orbital_dir ../../../tests/PP_ORB +nbands 16 +#Parameters (Accuracy) +ecutwfc 50 +scf_nmax 20 + +basis_type pw +gamma_only 1 + +psi_initializer 1 +init_wfc nao + +smearing_method gauss +smearing_sigma 0.01 + +calculation relax +relax_nmax 50 +force_thr_ev 1.0e-3 diff --git a/tests/integrate/102_PW_PINT_RKS/KPT b/tests/integrate/102_PW_PINT_RKS/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/integrate/102_PW_PINT_RKS/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/integrate/102_PW_PINT_RKS/STRU b/tests/integrate/102_PW_PINT_RKS/STRU new file mode 100644 index 0000000000..c08fb59d58 --- /dev/null +++ b/tests/integrate/102_PW_PINT_RKS/STRU @@ -0,0 +1,24 @@ +#This is the atom file containing all the information +#about the lattice structure. + +ATOMIC_SPECIES +Si 1.000 Si.pz-vbc.UPF #Element, Mass, Pseudopotential + +NUMERICAL_ORBITAL +Si_gga_8au_100Ry_2s2p1d.orb + +LATTICE_CONSTANT +10.2 #Lattice constant + +LATTICE_VECTORS +0.5 0.5 0.0 #Lattice vector 1 +0.5 0.0 0.5 #Lattice vector 2 +0.0 0.5 0.5 #Lattice vector 3 + +ATOMIC_POSITIONS +Cartesian #Cartesian(Unit is LATTICE_CONSTANT) +Si #Name of element +0.0 #Magnetic for this element. +2 #Number of atoms +0.00 0.00 0.00 0 0 0 #x,y,z, move_x, move_y, move_z +0.25 0.25 0.25 1 1 1 diff --git a/tests/integrate/102_PW_PINT_RKS/jd b/tests/integrate/102_PW_PINT_RKS/jd new file mode 100644 index 0000000000..5630f3c6fe --- /dev/null +++ b/tests/integrate/102_PW_PINT_RKS/jd @@ -0,0 +1 @@ +examples/scf/pw_Si2 with psi_initializer init_wfc nao (nspin = 1, relax) \ No newline at end of file diff --git a/tests/integrate/102_PW_PINT_RKS/result.ref b/tests/integrate/102_PW_PINT_RKS/result.ref new file mode 100644 index 0000000000..889ddac358 --- /dev/null +++ b/tests/integrate/102_PW_PINT_RKS/result.ref @@ -0,0 +1,3 @@ +etotref -198.3545030347299303 +etotperatomref -99.1772515174 +totaltimeref 1.71 diff --git a/tests/integrate/102_PW_PINT_UKS/INPUT b/tests/integrate/102_PW_PINT_UKS/INPUT new file mode 100644 index 0000000000..34cf9a8dc4 --- /dev/null +++ b/tests/integrate/102_PW_PINT_UKS/INPUT @@ -0,0 +1,25 @@ +INPUT_PARAMETERS +suffix autotest +ntype 2 + +calculation scf +ecutwfc 20 +scf_thr 1.0e-8 +scf_nmax 50 +out_chg 0 + +smearing_method gaussian +smearing_sigma 0.02 + +mixing_type pulay +mixing_beta 0.4 + +ks_solver cg +basis_type pw +gamma_only 1 +nspin 2 +pseudo_dir ../../../tests/PP_ORB +orbital_dir ../../../tests/PP_ORB + +psi_initializer 1 +init_wfc nao diff --git a/tests/integrate/102_PW_PINT_UKS/KPT b/tests/integrate/102_PW_PINT_UKS/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/integrate/102_PW_PINT_UKS/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/integrate/102_PW_PINT_UKS/STRU b/tests/integrate/102_PW_PINT_UKS/STRU new file mode 100644 index 0000000000..1c346f66f1 --- /dev/null +++ b/tests/integrate/102_PW_PINT_UKS/STRU @@ -0,0 +1,29 @@ +ATOMIC_SPECIES +Fe1 1.000 Fe_ONCV_PBE-1.0.upf +Fe2 1.000 Fe_ONCV_PBE-1.0.upf + +NUMERICAL_ORBITAL +Fe_gga_9au_100Ry_4s2p2d1f.orb +Fe_gga_9au_100Ry_4s2p2d1f.orb + +LATTICE_CONSTANT +15 + +LATTICE_VECTORS + 1.00 0.50 0.50 + 0.50 1.00 0.50 + 0.50 0.50 1.00 +ATOMIC_POSITIONS +Direct + +Fe1 +1.0 +1 +0.00 0.00 0.00 1 1 1 + +Fe2 +-1.0 +1 +0.50 0.50 0.50 1 1 1 + + diff --git a/tests/integrate/102_PW_PINT_UKS/jd b/tests/integrate/102_PW_PINT_UKS/jd new file mode 100644 index 0000000000..0314ed9249 --- /dev/null +++ b/tests/integrate/102_PW_PINT_UKS/jd @@ -0,0 +1 @@ +examples/smearing/lcao_fe with psi_initializer init_wfc nao (nspin = 2, scf) diff --git a/tests/integrate/102_PW_PINT_UKS/result.ref b/tests/integrate/102_PW_PINT_UKS/result.ref new file mode 100644 index 0000000000..5ac2ef8604 --- /dev/null +++ b/tests/integrate/102_PW_PINT_UKS/result.ref @@ -0,0 +1,3 @@ +etotref -6147.553111352082 +etotperatomref -3073.7765556760 +totaltimeref 16.68 diff --git a/tests/integrate/108_PW_RE_PINT_RKS/INPUT b/tests/integrate/108_PW_RE_PINT_RKS/INPUT new file mode 100644 index 0000000000..46a5ca1e42 --- /dev/null +++ b/tests/integrate/108_PW_RE_PINT_RKS/INPUT @@ -0,0 +1,27 @@ +INPUT_PARAMETERS +# Created by Atomic Simulation Enviroment +calculation cell-relax +suffix autotest +relax_nmax 100 +force_thr_ev 0.01 +stress_thr 1 +out_stru 1 +out_level ie +pseudo_rcut 10.0 +pseudo_mesh 1 +ecutwfc 30 +basis_type pw +ks_solver cg +smearing_method gaussian +smearing_sigma 0.01 +mixing_type broyden +mixing_beta 0.7 +scf_thr 1e-08 +cal_force 1 +cal_stress 1 +out_stru 1 +chg_extrap second-order #atomic; first-order; second-order; dm:coefficients of SIA +pseudo_dir ../../../tests/PP_ORB +orbital_dir ../../../tests/PP_ORB +psi_initializer 1 +init_wfc nao diff --git a/tests/integrate/108_PW_RE_PINT_RKS/KPT b/tests/integrate/108_PW_RE_PINT_RKS/KPT new file mode 100644 index 0000000000..f5f7f4ec34 --- /dev/null +++ b/tests/integrate/108_PW_RE_PINT_RKS/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/tests/integrate/108_PW_RE_PINT_RKS/STRU b/tests/integrate/108_PW_RE_PINT_RKS/STRU new file mode 100644 index 0000000000..56932f990f --- /dev/null +++ b/tests/integrate/108_PW_RE_PINT_RKS/STRU @@ -0,0 +1,24 @@ +ATOMIC_SPECIES +Al 26.982 Al.PD04.PBE.UPF upf201 + +NUMERICAL_ORBITAL +Al_gga_10au_100Ry_3s3p2d.orb + +LATTICE_CONSTANT +1.88972612546 + +LATTICE_VECTORS +4.83198403231 -0.389937941056 0.421696660303 #latvec1 +-0.224122113303 2.94867832278 -0.221983075925 #latvec2 +0.446072377462 -0.396375785848 4.82801496974 #latvec3 + +ATOMIC_POSITIONS +Direct + +Al #label +0 #magnetism +4 #number of atoms +0.00825359188008 0.99670712136 0.00678934462388 m 1 1 1 +0.506607213055 0.503316792705 0.00321624481765 m 1 1 1 +0.506746323934 0.00329313381946 0.503210864168 m 1 1 1 +0.00839287113151 0.496682952117 0.50678354639 m 1 1 1 diff --git a/tests/integrate/108_PW_RE_PINT_RKS/jd b/tests/integrate/108_PW_RE_PINT_RKS/jd new file mode 100644 index 0000000000..d5f1e84d62 --- /dev/null +++ b/tests/integrate/108_PW_RE_PINT_RKS/jd @@ -0,0 +1 @@ +integrated test 108_PW_RE_NEW with psi_initializer init_wfc nao (nspin = 1, cell-relax) diff --git a/tests/integrate/108_PW_RE_PINT_RKS/result.ref b/tests/integrate/108_PW_RE_PINT_RKS/result.ref new file mode 100644 index 0000000000..02f28ba88d --- /dev/null +++ b/tests/integrate/108_PW_RE_PINT_RKS/result.ref @@ -0,0 +1,5 @@ +etotref -253.8609235706759648 +etotperatomref -63.4652308927 +totalforceref 0.073112 +totalstressref 1.843207 +totaltimeref 9.15 diff --git a/tests/integrate/Autotest.sh b/tests/integrate/Autotest.sh index 833724ee3f..bfae98f45a 100755 --- a/tests/integrate/Autotest.sh +++ b/tests/integrate/Autotest.sh @@ -9,7 +9,8 @@ threshold=0.0000001 # check accuracy ca=8 # regex of case name -case="^[^#].*_.*$" +case="^[^#].*PINT.*$" +#case="^[^#].*_.*$" # enable AddressSanitizer sanitize=false diff --git a/tests/integrate/CASES b/tests/integrate/CASES index 4a87bd11ed..29267976a5 100644 --- a/tests/integrate/CASES +++ b/tests/integrate/CASES @@ -14,6 +14,8 @@ 101_PW_VW_pseudopots 102_PW_DA_davidson 102_PW_BPCG +102_PW_PINT_RKS +102_PW_PINT_UKS 103_PW_15_CF_CS_S1_smallg 103_PW_15_CF_CS_S2_smallg 103_PW_15_CS_CF @@ -47,6 +49,7 @@ 108_PW_RE_MB 108_PW_RE_MG 108_PW_RE_NEW +108_PW_RE_PINT_RKS 109_PW_CR 109_PW_CR_fix_a 109_PW_CR_fix_ab