From dbb28dcbbfe5c6f50a41dd3352e08af2bfde8fa1 Mon Sep 17 00:00:00 2001 From: Qianrui Liu <76200646+Qianruipku@users.noreply.github.com> Date: Sat, 4 Nov 2023 21:38:54 +0800 Subject: [PATCH] Feature: add initsto_ecut (#3149) * Feature: add initsto_ecut Stochastic wave functions are initialized with a large G-ball with cutoff energy of initsto_ecut. It can make SDFT with different wfcecuts comparable. It makes results of different cores equal. --- docs/advanced/input_files/input-main.md | 9 ++ .../module_mixing/broyden_mixing.h | 1 - source/module_esolver/esolver_sdft_pw.cpp | 13 ++- .../module_hamilt_pw/hamilt_stodft/sto_wf.cpp | 105 ++++++++++++++---- .../module_hamilt_pw/hamilt_stodft/sto_wf.h | 49 ++++---- source/module_hsolver/hsolver_pw_sdft.cpp | 2 +- source/module_io/DEFAULT_TYPE.conf | 1 + source/module_io/DEFAULT_VALUE.conf | 1 + source/module_io/input.cpp | 5 + source/module_io/input.h | 1 + source/module_io/parameter_pool.cpp | 4 + source/module_io/test/input_test.cpp | 2 + source/module_io/test/input_test_para.cpp | 1 + source/module_io/test/support/INPUT | 1 + source/module_io/test/support/witestfile | 1 + source/module_io/test/write_input_test.cpp | 3 + source/module_io/write_input.cpp | 1 + 17 files changed, 150 insertions(+), 50 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index aed16d90c1..95ce07dce7 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -88,6 +88,7 @@ - [emin\_sto](#emin_sto) - [emax\_sto](#emax_sto) - [seed\_sto](#seed_sto) + - [initsto\_ecut](#initsto_ecut) - [initsto\_freq](#initsto_freq) - [npart\_sto](#npart_sto) - [Geometry relaxation](#geometry-relaxation) @@ -1119,6 +1120,14 @@ These variables are used to control the parameters of stochastic DFT (SDFT), mi - -1: the seed is decided by time(NULL). - **Default**: 0 +### initsto_ecut + +- **Type**: Real +- **Availability**: [esolver_type](#esolver_type) = `sdft` +- **Description**: Stochastic wave functions are initialized in a large box generated by "4*`initsto_ecut`". `initsto_ecut` should be larger than [ecutwfc](#ecutwfc). In this method, SDFT results are the same when using different cores. Besides, coefficients of the same G are the same when ecutwfc is rising to initsto_ecut. If it is smaller than [ecutwfc](#ecutwfc), it will be turned off. +- **Default**: 0.0 +- **Unit**: Ry + ### initsto_freq - **Type**: Integer diff --git a/source/module_base/module_mixing/broyden_mixing.h b/source/module_base/module_mixing/broyden_mixing.h index cc94bd471a..4f3db44f1d 100755 --- a/source/module_base/module_mixing/broyden_mixing.h +++ b/source/module_base/module_mixing/broyden_mixing.h @@ -277,7 +277,6 @@ class Broyden_Mixing : public Mixing } else { - beta(0, 0) = inner_product(FP_dF, FP_dF); coef[0] = 1.0; } diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index 808dac0220..c5b033dc59 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -78,8 +78,17 @@ void ESolver_SDFT_PW::Init(Input& inp, UnitCell& ucell) this->Init_GlobalC(inp, ucell); // temporary stowf.init(&kv, pw_wfc->npwk_max); - if (INPUT.nbands_sto != 0) - Init_Sto_Orbitals(this->stowf, inp.seed_sto); + if (inp.nbands_sto != 0) + { + if (inp.initsto_ecut < inp.ecutwfc) + { + Init_Sto_Orbitals(this->stowf, inp.seed_sto); + } + else + { + Init_Sto_Orbitals_Ecut(this->stowf, inp.seed_sto, kv, *pw_wfc, inp.initsto_ecut); + } + } else Init_Com_Orbitals(this->stowf); size_t size = stowf.chi0->size(); diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_wf.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_wf.cpp index 38fef87093..9013dde133 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_wf.cpp +++ b/source/module_hamilt_pw/hamilt_stodft/sto_wf.cpp @@ -47,6 +47,12 @@ void Init_Sto_Orbitals(Stochastic_WF& stowf, const int seed_in) srand((unsigned)std::abs(seed_in) + GlobalV::MY_RANK * 10000); } + Allocate_Chi0(stowf); + Update_Sto_Orbitals(stowf, seed_in); +} + +void Allocate_Chi0(Stochastic_WF& stowf) +{ bool firstrankmore = false; int igroup = 0; // I am not sure which is better. @@ -83,30 +89,6 @@ void Init_Sto_Orbitals(Stochastic_WF& stowf, const int seed_in) { stowf.nchip[ik] = tmpnchip; } - stowf.chi0->fix_k(0); - - if (seed_in >= 0) - { - for (int i = 0; i < size; ++i) - { - const double phi = 2 * ModuleBase::PI * rand() / double(RAND_MAX); - stowf.chi0->get_pointer()[i] = std::complex(cos(phi), sin(phi)) / sqrt(double(nchi)); - } - } - else - { - for (int i = 0; i < size; ++i) - { - if (rand() / double(RAND_MAX) < 0.5) - { - stowf.chi0->get_pointer()[i] = -1.0 / sqrt(double(nchi)); - } - else - { - stowf.chi0->get_pointer()[i] = 1.0 / sqrt(double(nchi)); - } - } - } } void Update_Sto_Orbitals(Stochastic_WF& stowf, const int seed_in) @@ -245,3 +227,78 @@ void Init_Com_Orbitals(Stochastic_WF& stowf) } } #endif +void Init_Sto_Orbitals_Ecut(Stochastic_WF& stowf, + const int seed_in, + const K_Vectors& kv, + const ModulePW::PW_Basis_K& wfcpw, + const int max_ecut) +{ + Allocate_Chi0(stowf); + + ModulePW::PW_Basis pwmax; +#ifdef __MPI + pwmax.initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); +#endif + pwmax.initgrids(wfcpw.lat0, wfcpw.latvec, max_ecut); + const int nx = pwmax.nx; + const int ny = pwmax.ny; + const int nz = pwmax.nz; + const int nkstot = kv.nkstot; + const int nks = kv.nks; + const int nchitot = INPUT.nbands_sto; + bool* updown = new bool[nx * ny * nz]; + int* nrecv = new int[GlobalV::NSTOGROUP]; + const int nchiper = stowf.nchip[0]; +#ifdef __MPI + MPI_Allgather(&nchiper, 1, MPI_INT, nrecv, 1, MPI_INT, PARAPW_WORLD); +#endif + int ichi_start = 0; + for (int i = 0; i < GlobalV::MY_STOGROUP; ++i) + { + ichi_start += nrecv[i]; + } + + for (int ik = 0; ik < nks; ++ik) + { + const int iktot = kv.getik_global(ik); + const int npw = wfcpw.npwk[ik]; + int* ig2ixyz = new int[npw]; + + for (int ig = 0; ig < npw; ++ig) + { + ModuleBase::Vector3 gdirect = wfcpw.getgdirect(ik, ig); + int ix = static_cast(gdirect.x); + int iy = static_cast(gdirect.y); + int iz = static_cast(gdirect.z); + ix = (ix + nx) % nx; + iy = (iy + ny) % ny; + iz = (iz + nz) % nz; + ig2ixyz[ig] = ix * ny * nz + iy * nz + iz; + } + + for (int ichi = 0; ichi < nchiper; ++ichi) + { + unsigned int seed = std::abs(seed_in) * (nkstot * nchitot) + iktot * nchitot + ichi_start + ichi; + srand(seed); + for (int i = 0; i < nx * ny * nz; ++i) + { + updown[i] = (rand() / double(RAND_MAX) < 0.5); + } + + for (int ig = 0; ig < npw; ++ig) + { + if (updown[ig2ixyz[ig]]) + { + stowf.chi0->operator()(ik, ichi, ig) = -1.0 / sqrt(double(nchitot)); + } + else + { + stowf.chi0->operator()(ik, ichi, ig) = 1.0 / sqrt(double(nchitot)); + } + } + } + delete[] ig2ixyz; + } + delete[] nrecv; + delete[] updown; +} diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_wf.h b/source/module_hamilt_pw/hamilt_stodft/sto_wf.h index 04203ed512..f198839f2d 100644 --- a/source/module_hamilt_pw/hamilt_stodft/sto_wf.h +++ b/source/module_hamilt_pw/hamilt_stodft/sto_wf.h @@ -1,23 +1,22 @@ #ifndef STOCHASTIC_WF_H #define STOCHASTIC_WF_H -#include "module_cell/klist.h" #include "module_base/complexmatrix.h" +#include "module_cell/klist.h" #include "module_psi/psi.h" +#include "module_basis/module_pw/pw_basis_k.h" //---------------------------------------------- -// Generate stochastic wave functions +// Generate stochastic wave functions //---------------------------------------------- -class Stochastic_WF +class Stochastic_WF { - public: - + public: Stochastic_WF(); ~Stochastic_WF(); - - void init(K_Vectors* p_kv, const int npwx_in); + void init(K_Vectors* p_kv, const int npwx_in); // origin stochastic wavefunctions in real space psi::Psi>* chi0 = nullptr; @@ -25,22 +24,28 @@ class Stochastic_WF psi::Psi>* chiortho = nullptr; // sqrt(f(H))|chi> psi::Psi>* shchi = nullptr; - int nchi = 0; // Total number of stochatic obitals - int *nchip = nullptr; // The number of stochatic orbitals in current process of each k point. - int nchip_max = 0; // Max number of stochastic orbitals among all k points. - int nks = 0; //number of k-points - int *ngk = nullptr; // ngk in klist - int npwx = 0; // max ngk[ik] in all processors - - int nbands_diag = 0; // number of bands obtained from diagonalization - - int nbands_total = 0; // number of bands in total, nbands_total=nchi+nbands_diag; - + int nchi = 0; ///< Total number of stochatic obitals + int* nchip = nullptr; ///< The number of stochatic orbitals in current process of each k point. + int nchip_max = 0; ///< Max number of stochastic orbitals among all k points. + int nks = 0; ///< number of k-points + int* ngk = nullptr; ///< ngk in klist + int npwx = 0; ///< max ngk[ik] in all processors + int nbands_diag = 0; ///< number of bands obtained from diagonalization + int nbands_total = 0; ///< number of bands in total, nbands_total=nchi+nbands_diag; }; -//init stochastic orbitals +// init stochastic orbitals void Init_Sto_Orbitals(Stochastic_WF& stowf, const int seed_in); -//update stochastic orbitals +// init stochastic orbitals from a large Ecut +// It can test the convergence of SDFT with respect to Ecut +void Init_Sto_Orbitals_Ecut(Stochastic_WF& stowf, + const int seed_in, + const K_Vectors& kv, + const ModulePW::PW_Basis_K& wfcpw, + const int max_ecut); +// allocate chi0 +void Allocate_Chi0(Stochastic_WF& stowf); +// update stochastic orbitals void Update_Sto_Orbitals(Stochastic_WF& stowf, const int seed_in); -//init complete orbitals +// init complete orbitals void Init_Com_Orbitals(Stochastic_WF& stowf); -#endif +#endif diff --git a/source/module_hsolver/hsolver_pw_sdft.cpp b/source/module_hsolver/hsolver_pw_sdft.cpp index 821512c6f8..4b1f6e79ee 100644 --- a/source/module_hsolver/hsolver_pw_sdft.cpp +++ b/source/module_hsolver/hsolver_pw_sdft.cpp @@ -51,7 +51,7 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt>* pHamilt, #ifdef __MPI if(nbands > 0) { - MPI_Bcast(&psi(ik,0,0), npwx*nbands*2, MPI_DOUBLE , 0, PARAPW_WORLD); + MPI_Bcast(&psi(ik,0,0), npwx*nbands, MPI_DOUBLE_COMPLEX , 0, PARAPW_WORLD); MPI_Bcast(&(pes->ekb(ik, 0)), nbands, MPI_DOUBLE, 0, PARAPW_WORLD); } #endif diff --git a/source/module_io/DEFAULT_TYPE.conf b/source/module_io/DEFAULT_TYPE.conf index 1e1402909a..886d3a692b 100644 --- a/source/module_io/DEFAULT_TYPE.conf +++ b/source/module_io/DEFAULT_TYPE.conf @@ -35,6 +35,7 @@ nche_sto int nbands_sto int nbndsto_str string seed_sto int +initsto_ecut double emax_sto double emin_sto double bndpar int diff --git a/source/module_io/DEFAULT_VALUE.conf b/source/module_io/DEFAULT_VALUE.conf index 7e1fe5a2ab..7eebdaa45e 100644 --- a/source/module_io/DEFAULT_VALUE.conf +++ b/source/module_io/DEFAULT_VALUE.conf @@ -20,6 +20,7 @@ emax_sto 0.0 nche_sto 100 seed_sto 0 + initsto_ecut 0.0 bndpar 1 kpar 1 initsto_freq 0 diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index b5ee3041b4..a2dd3b4dd7 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -758,6 +758,10 @@ bool Input::Read(const std::string &fn) { read_value(ifs, seed_sto); } + else if (strcmp("initsto_ecut", word) == 0) + { + read_value(ifs, initsto_ecut); + } else if (strcmp("pw_seed", word) == 0) { read_value(ifs, pw_seed); @@ -2974,6 +2978,7 @@ void Input::Bcast() Parallel_Common::bcast_double(min_dist_coef); Parallel_Common::bcast_int(nche_sto); Parallel_Common::bcast_int(seed_sto); + Parallel_Common::bcast_double(initsto_ecut); Parallel_Common::bcast_int(pw_seed); Parallel_Common::bcast_double(emax_sto); Parallel_Common::bcast_double(emin_sto); diff --git a/source/module_io/input.h b/source/module_io/input.h index d3e45211b7..ef9abe46b9 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -84,6 +84,7 @@ class Input int nbands_sto; // number of stochastic bands //qianrui 2021-2-5 std::string nbndsto_str; // string parameter for stochastic bands int seed_sto; // random seed for sDFT + double initsto_ecut = 0.0; // maximum ecut to init stochastic bands double emax_sto; // Emax & Emin to normalize H double emin_sto; int bndpar; //parallel for stochastic/deterministic bands diff --git a/source/module_io/parameter_pool.cpp b/source/module_io/parameter_pool.cpp index 7858ca25ca..fcced938fa 100644 --- a/source/module_io/parameter_pool.cpp +++ b/source/module_io/parameter_pool.cpp @@ -479,6 +479,10 @@ bool input_parameters_set(std::map input_parameters { INPUT.seed_sto = *static_cast(input_parameters["seed_sto"].get()); } + else if (input_parameters.count("initsto_ecut") != 0) + { + INPUT.initsto_ecut = *static_cast(input_parameters["initsto_ecut"].get()); + } else if (input_parameters.count("emax_sto") != 0) { INPUT.emax_sto = *static_cast(input_parameters["emax_sto"].get()); diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index c3d0113cb0..2b7adea125 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -47,6 +47,7 @@ TEST_F(InputTest, Default) EXPECT_EQ(INPUT.emax_sto,0.0); EXPECT_EQ(INPUT.nche_sto,100); EXPECT_EQ(INPUT.seed_sto,0); + EXPECT_EQ(INPUT.initsto_ecut,0.0); EXPECT_EQ(INPUT.bndpar,1); EXPECT_EQ(INPUT.kpar,1); EXPECT_EQ(INPUT.initsto_freq,0); @@ -403,6 +404,7 @@ TEST_F(InputTest, Read) EXPECT_EQ(INPUT.emax_sto,0.0); EXPECT_EQ(INPUT.nche_sto,100); EXPECT_EQ(INPUT.seed_sto,0); + EXPECT_EQ(INPUT.initsto_ecut,0.0); EXPECT_EQ(INPUT.bndpar,1); EXPECT_EQ(INPUT.kpar,1); EXPECT_EQ(INPUT.initsto_freq,0); diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp index 1998d9cc77..95e908c050 100644 --- a/source/module_io/test/input_test_para.cpp +++ b/source/module_io/test/input_test_para.cpp @@ -54,6 +54,7 @@ TEST_F(InputParaTest,Bcast) EXPECT_EQ(INPUT.emax_sto,0.0); EXPECT_EQ(INPUT.nche_sto,100); EXPECT_EQ(INPUT.seed_sto,0); + EXPECT_EQ(INPUT.initsto_ecut,0.0); EXPECT_EQ(INPUT.bndpar,1); EXPECT_EQ(INPUT.kpar,1); EXPECT_EQ(INPUT.initsto_freq,0); diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT index eceb942dd2..5bb4efadba 100644 --- a/source/module_io/test/support/INPUT +++ b/source/module_io/test/support/INPUT @@ -77,6 +77,7 @@ nche_sto 100 #Chebyshev expansion orders emin_sto 0 #trial energy to guess the lower bound of eigen energies of the Hamitonian operator emax_sto 0 #trial energy to guess the upper bound of eigen energies of the Hamitonian operator seed_sto 0 #the random seed to generate stochastic orbitals +initsto_ecut 0 #maximum ecut to init stochastic bands initsto_freq 0 #frequency to generate new stochastic orbitals when running md cal_cond 0 #calculate electronic conductivities cond_che_thr 1e-08 #control the error of Chebyshev expansions for conductivities diff --git a/source/module_io/test/support/witestfile b/source/module_io/test/support/witestfile index 6d8b4bb95d..9bf440281c 100644 --- a/source/module_io/test/support/witestfile +++ b/source/module_io/test/support/witestfile @@ -74,6 +74,7 @@ nche_sto 100 #Chebyshev expansion orders emin_sto 0 #trial energy to guess the lower bound of eigen energies of the Hamitonian operator emax_sto 0 #trial energy to guess the upper bound of eigen energies of the Hamitonian operator seed_sto 0 #the random seed to generate stochastic orbitals +initsto_ecut 0 #maximum ecut to init stochastic bands initsto_freq 0 #frequency to generate new stochastic orbitals when running md cal_cond 0 #calculate electronic conductivities cond_che_thr 1e-08 #control the error of Chebyshev expansions for conductivities diff --git a/source/module_io/test/write_input_test.cpp b/source/module_io/test/write_input_test.cpp index ba1b2ce25a..164a22c4dc 100644 --- a/source/module_io/test/write_input_test.cpp +++ b/source/module_io/test/write_input_test.cpp @@ -202,6 +202,9 @@ TEST_F(write_input, STO3) EXPECT_THAT( output, testing::HasSubstr("seed_sto 0 #the random seed to generate stochastic orbitals")); + EXPECT_THAT( + output, + testing::HasSubstr("initsto_ecut 0 #maximum ecut to init stochastic bands")); EXPECT_THAT(output, testing::HasSubstr( "initsto_freq 0 #frequency to generate new stochastic orbitals when running md")); diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index a680fb8284..880fdc91c4 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -146,6 +146,7 @@ void Input::Print(const std::string &fn) const ModuleBase::GlobalFunc::OUTP(ofs, "emin_sto", emin_sto, "trial energy to guess the lower bound of eigen energies of the Hamitonian operator"); ModuleBase::GlobalFunc::OUTP(ofs, "emax_sto", emax_sto, "trial energy to guess the upper bound of eigen energies of the Hamitonian operator"); ModuleBase::GlobalFunc::OUTP(ofs, "seed_sto", seed_sto, "the random seed to generate stochastic orbitals"); + ModuleBase::GlobalFunc::OUTP(ofs, "initsto_ecut", initsto_ecut, "maximum ecut to init stochastic bands"); ModuleBase::GlobalFunc::OUTP(ofs, "initsto_freq", initsto_freq, "frequency to generate new stochastic orbitals when running md"); ModuleBase::GlobalFunc::OUTP(ofs, "cal_cond", cal_cond, "calculate electronic conductivities"); ModuleBase::GlobalFunc::OUTP(ofs, "cond_che_thr", cond_che_thr, "control the error of Chebyshev expansions for conductivities");