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

@@ -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

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;
         }
 

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();

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<double>(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<double> gdirect = wfcpw.getgdirect(ik, ig);
+            int ix = static_cast<int>(gdirect.x);
+            int iy = static_cast<int>(gdirect.y);
+            int iz = static_cast<int>(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;
+}

source/module_hamilt_pw/hamilt_stodft/sto_wf.h

@@ -1,46 +1,51 @@
 #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<std::complex<double>>* chi0 = nullptr;
     // stochastic wavefunctions after in reciprocal space orthogonalized with KS wavefunctions
     psi::Psi<std::complex<double>>* chiortho = nullptr;
     // sqrt(f(H))|chi>
     psi::Psi<std::complex<double>>* 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

source/module_hsolver/hsolver_pw_sdft.cpp

@@ -51,7 +51,7 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt<std::complex<double>>* 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

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

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

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);

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

source/module_io/parameter_pool.cpp

@@ -479,6 +479,10 @@ bool input_parameters_set(std::map<std::string, InputParameter> input_parameters
     {
         INPUT.seed_sto = *static_cast<int*>(input_parameters["seed_sto"].get());
     }
+    else if (input_parameters.count("initsto_ecut") != 0)
+    {
+        INPUT.initsto_ecut = *static_cast<double*>(input_parameters["initsto_ecut"].get());
+    }
     else if (input_parameters.count("emax_sto") != 0)
     {
         INPUT.emax_sto = *static_cast<double*>(input_parameters["emax_sto"].get());

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);

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);

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

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

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"));

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");