Merge branch 'develop' into stress_ewa

Author: dzzz2001

docs/advanced/input_files/input-main.md

@@ -302,6 +302,9 @@
   - [Exact Exchange (PW)](#exact-exchange-pw)
     - [exxace](#exxace)
     - [exx\_gamma\_extrapolation](#exx_gamma_extrapolation)
+    - [ecutexx](#ecutexx)
+    - [exx_thr_type](#exx_thr_type)
+    - [exx_ene_thr](#exx_ene_thr)
   - [Molecular Dynamics](#molecular-dynamics)
     - [md\_type](#md_type)
     - [md\_nstep](#md_nstep)
@@ -3101,6 +3104,26 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b
 - **Description**: Whether to use the gamma point extrapolation method to calculate the Fock exchange operator. See [https://doi.org/10.1103/PhysRevB.79.205114](https://doi.org/10.1103/PhysRevB.79.205114) for details. Should be set to true most of the time.
 - **Default**: True
 
+### ecutexx
+- **Type**: Real
+- **Description**: The energy cutoff for EXX (Fock) exchange operator in plane wave basis calculations. Reducing `ecutexx` below `ecutrho` may significantly accelerate EXX computations. This speed improvement comes with a reduced numerical accuracy in the exchange energy calculation.
+- **Default**: same as *[ecutrho](#ecutrho)*
+- **Unit**: Ry
+
+### exx_thr_type
+- **Type**: String
+- **Description**: The type of threshold used to judge whether the outer loop has converged in the separate loop EXX calculation.
+  - energy: use the change of exact exchange energy to judge convergence.
+  - density: if the change of charge density difference between two successive outer loop iterations is seen as converged according to *[scf_thr](#scf_thr)*, then the outer loop is seen as converged.
+- **Default**: `density`
+
+### exx_ene_thr
+- **Type**: Real
+- **Availability**: *[exx_thr_type](#exx_thr_type)*==`energy`
+- **Description**: The threshold for the change of exact exchange energy to judge convergence of the outer loop in the separate loop EXX calculation.
+- **Default**: 1e-5
+- **Unit**: Ry
+
 ## Molecular dynamics
 
 These variables are used to control molecular dynamics calculations. For more information, please refer to [md.md](../md.md#molecular-dynamics) in detail.

source/source_hamilt/module_vdw/vdwd3_autoset_xcparam.cpp

@@ -243,6 +243,9 @@ const std::pair<const char*, std::vector<double>> zero_data[] = {
     {"b97_d", {1.0, 0.892, 0.892, 0.909, 1.0, 1.0, 1.0, 14.0, 0.0}},
     {"b973c", {1.0, 1.06, 1.06, 1.5, 1.0, 1.0, 1.0, 14.0, 0.0}},
     {"pbe", {1.0, 1.217, 1.217, 0.722, 1.0, 1.0, 1.0, 14.0, 0.0}},
+    {"pbesol", {1.0, 1.345, 1.345, 0.612, 1.0, 1.0, 1.0, 14.0, 0.0}}, 
+    // issue#6646, d3 zero-damping support for PBEsol, 
+    // parameters retrived from https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3/zero_damping
     {"rpw86pbe", {1.0, 1.224, 1.224, 0.901, 1.0, 1.0, 1.0, 14.0, 0.0}},
     {"b3lyp", {1.0, 1.261, 1.261, 1.703, 1.0, 1.0, 1.0, 14.0, 0.0}},
     {"tpss", {1.0, 1.166, 1.166, 1.105, 1.0, 1.0, 1.0, 14.0, 0.0}},

source/source_io/cal_pLpR.cpp

@@ -150,12 +150,21 @@ ModuleIO::AngularMomentumCalculator::AngularMomentumCalculator(
 
     // for neighbor list search
     double temp = -1.0;
+    if (search_radius < rcut_max)
+    {
+        *ofs_ << "Find the `search_radius` from the input file being smaller than the \n"
+                 "`rcut_max` of the orbitals.\n"
+              << "Reset the `search_radius` (" << search_radius << ") "
+              << "to `rcut_max` ("<< rcut_max << ")." 
+              << std::endl;
+        // we don't really set, but use std::max to mask :)
+    }
     temp = atom_arrange::set_sr_NL(*ofs_,
                                    PARAM.inp.out_level,
-                                   search_radius,
+                                   std::max(search_radius, rcut_max),
                                    ucell.infoNL.get_rcutmax_Beta(),
                                    PARAM.globalv.gamma_only_local);
-    temp = std::max(temp, search_radius);
+    temp = std::max(temp, std::max(search_radius, rcut_max));
     this->neighbor_searcher_ = std::unique_ptr<Grid_Driver>(new Grid_Driver(tdestructor, tgrid));
     atom_arrange::search(searchpbc,
                          *ofs_,
@@ -190,24 +199,29 @@ void ModuleIO::AngularMomentumCalculator::kernel(
     fmtstr += "%" + std::to_string(precision*2) + "." + std::to_string(precision) + "e\n";
     FmtCore fmt(fmtstr);
 
-    ModuleBase::Vector3<double> ri, rj, dr;
+    // placeholders
+    std::complex<double> val = 0;
+    ModuleBase::Vector3<double> taui; // the origin position
+    ModuleBase::Vector3<double> dtau; // the displacement
+    AdjacentAtomInfo adjinfo; // adjacent atom information carrier
     for (int it = 0; it < ucell.ntype; it++)
     {
         const Atom& atyp_i = ucell.atoms[it];
         for (int ia = 0; ia < atyp_i.na; ia++)
         {
-            ri = atyp_i.tau[ia];
-            neighbor_searcher_->Find_atom(ucell, ri, it, ia);
-            for (int ia_adj = 0; ia_adj < neighbor_searcher_->getAdjacentNum(); ia_adj++)
+            taui = ucell.get_tau(ucell.itia2iat(it, ia));
+            neighbor_searcher_->Find_atom(ucell, taui, it, ia, &adjinfo);
+            for (int ia_adj = 0; ia_adj < adjinfo.adj_num + 1; ia_adj++) // "+1" is to include itself
             {
-                rj = neighbor_searcher_->getAdjacentTau(ia_adj);
-                int jt = neighbor_searcher_->getType(ia_adj);
+                int jt = adjinfo.ntype[ia_adj]; // ityp
+                int ja = adjinfo.natom[ia_adj]; // iat with in atomtype
                 const Atom& atyp_j = ucell.atoms[jt];
-                int ja = neighbor_searcher_->getNatom(ia_adj);
-                dr = (ri - rj) * ucell.lat0;
-                const ModuleBase::Vector3<int> iR = neighbor_searcher_->getBox(ia_adj);
-                // the two-center-integral
+                const ModuleBase::Vector3<int> iR = adjinfo.box[ia_adj];
+                dtau = ucell.cal_dtau(ucell.itia2iat(it, ia), 
+                                      ucell.itia2iat(jt, ja), 
+                                      iR) * ucell.lat0; // convert to unit of Bohr
 
+                // nested loop: calculate the two-center-integral
                 for (int li = 0; li < atyp_i.nwl + 1; li++)
                 {
                     for (int iz = 0; iz < atyp_i.l_nchi[li]; iz++)
@@ -220,21 +234,20 @@ void ModuleIO::AngularMomentumCalculator::kernel(
                                 {
                                     for (int mj = -lj; mj <= lj; mj++)
                                     {
-                                        std::complex<double> val = 0;
                                         if (dir == 'x')
                                         {
                                             val = cal_LxijR(calculator_, 
-                                                it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
+                                                it, ia, li, iz, mi, jt, ja, lj, jz, mj, dtau);
                                         }
                                         else if (dir == 'y')
                                         {
                                             val = cal_LyijR(calculator_, 
-                                                it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
+                                                it, ia, li, iz, mi, jt, ja, lj, jz, mj, dtau);
                                         }
                                         else if (dir == 'z')
                                         {
                                             val = cal_LzijR(calculator_, 
-                                                it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
+                                                it, ia, li, iz, mi, jt, ja, lj, jz, mj, dtau);
                                         }
 
                                         *ofs << fmt.format(
@@ -266,7 +279,7 @@ void ModuleIO::AngularMomentumCalculator::calculate(
     }
     std::ofstream ofout;
     const std::string dir = "xyz";
-    const std::string title = "# it ia il iz im iRx iRy iRz jt ja jl jz jm <a|L|b>\n"
+    const std::string title = "# it ia il iz im iRx iRy iRz jt ja jl jz jm Re[<a|L|b>] Im[<a|L|b>]\n"
                               "# it: atomtype index of the first atom\n"
                               "# ia: atomic index of the first atom within the atomtype\n"
                               "# il: angular momentum index of the first atom\n"
@@ -278,7 +291,8 @@ void ModuleIO::AngularMomentumCalculator::calculate(
                               "# jl: angular momentum index of the second atom\n"
                               "# jz: zeta function index of the second atom\n"
                               "# jm: magnetic quantum number of the second atom\n"
-                              "# <a|L|b>: the value of the matrix element\n";
+                              "# Re[<a|L|b>], Im[<a|L|b>]: the real and imaginary parts "
+                              "of the value of the matrix element\n";
 
     for (char d : dir)
     {

source/source_io/module_parameter/input_parameter.h

@@ -656,6 +656,7 @@ struct Input_para
     bool exx_gamma_extrapolation = true; // gamma point extrapolation for exx, https://doi.org/10.1103/PhysRevB.79.205114
     std::string exx_thr_type = "density"; // threshold type for exx outer loop, energy or density
     double exx_ene_thr = 1e-5; // threshold for exx outer loop when exx_thr_type = energy
+    double ecutexx = 0.0; // energy cutoff for exx calculation, Ry
 
     // ====   #Parameters (23.XC external parameterization) ========
     /*

source/source_io/read_input_item_other.cpp

@@ -571,6 +571,18 @@ void ReadInput::item_others()
         };
         this->add_item(item);
     }
+    {
+        Input_Item item("ecutexx");
+        item.annotation = "energy cutoff for exx calculation, Ry";
+        read_sync_double(input.ecutexx);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.ecutexx < 0)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "ecutexx must >= 0");
+            }
+        };
+        this->add_item(item);
+    }
 
 }
 } // namespace ModuleIO

source/source_pw/module_pwdft/operator_pw/op_exx_pw.cpp

@@ -76,15 +76,22 @@ OperatorEXXPW<T, Device>::OperatorEXXPW(const int* isk_in,
     tpiba = ucell->tpiba;
     Real tpiba2 = tpiba * tpiba;
 
+    // initialize rhopw_dev
+    double ecut_exx = PARAM.inp.ecutexx;
+    if (ecut_exx == 0.0)
+    {
+        ecut_exx = PARAM.inp.ecutrho;
+    }
+
     rhopw_dev = new ModulePW::PW_Basis(wfcpw->get_device(), rhopw->get_precision());
     rhopw_dev->fft_bundle.setfft(wfcpw->get_device(), rhopw->get_precision());
 #ifdef __MPI
     rhopw_dev->initmpi(rhopw->poolnproc, rhopw->poolrank, rhopw->pool_world);
 #endif
     // here we can actually use different ecut to init the grids
-    rhopw_dev->initgrids(rhopw->lat0, rhopw->latvec, rhopw->gridecut_lat * rhopw->tpiba2);
+    rhopw_dev->initgrids(rhopw->lat0, rhopw->latvec, ecut_exx);
     rhopw_dev->initgrids(rhopw->lat0, rhopw->latvec, rhopw->nx, rhopw->ny, rhopw->nz);
-    rhopw_dev->initparameters(rhopw->gamma_only, rhopw->ggecut * rhopw->tpiba2, rhopw->distribution_type, rhopw->xprime);
+    rhopw_dev->initparameters(rhopw->gamma_only, ecut_exx, rhopw->distribution_type, rhopw->xprime);
     rhopw_dev->setuptransform();
     rhopw_dev->collect_local_pw();