Feature: Calculate and output electron localization function (ELF) wi…

…th KSDFT and OFDFT (#5139)

* Feature: Calculate and output electron localization function (ELF) with KSDFT in both PW and LCAO basis.

* Feature: Calculate and output ELF with OFDFT

* Test: Add three integrate tests for out_elf
119_PW_out_elf, 219_NO_out_elf, 919_OF_out_elf

* Test: Update tools/catch_properties.sh

* Feature: Support ELF in SPIN2 case.

* Doc: Update the document of input parameters

* Fix: Update Makefile.Objects

* [pre-commit.ci lite] apply automatic fixes

* Fix: Fix the bug caused by check_value

* Fix: Update the name of several parameters, fix a memory leak

* Feature: support user defined precision of elf

* Doc: Update document of out_elf

* Doc/Test: Update the doc and test

* Test: Update unit test

* Refactor: Remove GlobalV::NPSIN from source/module_hamilt_pw/hamilt_ofdft/

* [pre-commit.ci lite] apply automatic fixes

* Fix: Add a small number (1e-5) onto the numerator of chi to suppress the numerical instability caused by LCAO basis, as suggested by 卢天, 陈飞武. 电子定域化函数的含义与函数形式[J]. 物理化学学报, 2011, 27(12): 2786-2792.

* Test: Update 119_PW_out_elf and 919_OF_out_elf

* Test: Correct the reference of 119_PW_out_elf

* Test: Update 219_LCAO_out_elf

* Refactor: Move cal_tau to new files.

* [pre-commit.ci lite] apply automatic fixes

* Test: Update abacus-develop/source/module_elecstate/test/CMakeLists.txt

---------

Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>
Co-authored-by: Mohan Chen <mohanchen@pku.edu.cn>

docs/advanced/input_files/input-main.md

@@ -161,6 +161,7 @@
     - [nbands\_istate](#nbands_istate)
     - [bands\_to\_print](#bands_to_print)
     - [if\_separate\_k](#if_separate_k)
+    - [out\_elf](#out_elf)
   - [Density of states](#density-of-states)
     - [dos\_edelta\_ev](#dos_edelta_ev)
     - [dos\_sigma](#dos_sigma)
@@ -1521,7 +1522,7 @@ These variables are used to control the output of properties.
 - **Type**: Integer \[Integer\](optional)
 - **Description**: 
   The first integer controls whether to output the charge density on real space grids:
-  - 1. Output the charge density (in Bohr^-3) on real space grids into the density files in the folder `OUT.${suffix}`. The files are named as:
+  - 1: Output the charge density (in Bohr^-3) on real space grids into the density files in the folder `OUT.${suffix}`. The files are named as:
     - nspin = 1: SPIN1_CHG.cube;
     - nspin = 2: SPIN1_CHG.cube, and SPIN2_CHG.cube;
     - nspin = 4: SPIN1_CHG.cube, SPIN2_CHG.cube, SPIN3_CHG.cube, and SPIN4_CHG.cube.
@@ -1801,6 +1802,23 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t
 - **Description**: Specifies whether to write the partial charge densities for all k-points to individual files or merge them. **Warning**: Enabling symmetry may produce incorrect results due to incorrect k-point weights. Therefore, when calculating partial charge densities, it is strongly recommended to set `symmetry = -1`.
 - **Default**: false
 
+### out_elf
+
+- **Type**: Integer \[Integer\](optional)
+- **Availability**: Only for Kohn-Sham DFT and Orbital Free DFT.
+- **Description**: Whether to output the electron localization function (ELF) in the folder `OUT.${suffix}`. The files are named as 
+    - nspin = 1:
+      - ELF.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$;
+    - nspin = 2:
+      - ELF_SPIN1.cube, ELF_SPIN2.cube: ${\rm{ELF}}_\sigma = \frac{1}{1+\chi_\sigma^2}$, $\chi_\sigma = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i,\sigma}|^2} - \frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}$;
+      - ELF.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i,\sigma}{f_i |\nabla\psi_{i,\sigma}|^2} - \sum_{\sigma}{\frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}}{\sum_{\sigma}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}}$;
+
+  The second integer controls the precision of the kinetic energy density output, if not given, will use `3` as default. For purpose restarting from this file and other high-precision involved calculation, recommend to use `10`.
+
+  ---
+  In molecular dynamics calculations, the output frequency is controlled by [out_interval](#out_interval).
+- **Default**: 0 3
+
 [back to top](#full-list-of-input-keywords)
 
 ## Density of states

source/Makefile.Objects

@@ -213,6 +213,7 @@ OBJS_ELECSTAT=elecstate.o\
     elecstate_print.o\
     elecstate_pw.o\
     elecstate_pw_sdft.o\
+    elecstate_pw_cal_tau.o\
     elecstate_op.o\
     efield.o\
     gatefield.o\
@@ -226,6 +227,7 @@ OBJS_ELECSTAT=elecstate.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
       elecstate_lcao_tddft.o\
+      elecstate_lcao_cal_tau.o\
       density_matrix.o\
       cal_dm_psi.o\
 
@@ -496,6 +498,7 @@ OBJS_IO=input_conv.o\
     winput.o\
     write_cube.o\
     write_elecstat_pot.o\
+    write_elf.o\
     write_dipole.o\
     td_current_io.o\
     write_wfc_r.o\

source/module_elecstate/CMakeLists.txt

@@ -7,6 +7,7 @@ list(APPEND objects
     elecstate_print.cpp
     elecstate_pw.cpp
     elecstate_pw_sdft.cpp
+    elecstate_pw_cal_tau.cpp
     potentials/gatefield.cpp
     potentials/efield.cpp
     potentials/H_Hartree_pw.cpp
@@ -31,6 +32,7 @@ if(ENABLE_LCAO)
   list(APPEND objects
       elecstate_lcao.cpp
       elecstate_lcao_tddft.cpp
+      elecstate_lcao_cal_tau.cpp
       potentials/H_TDDFT_pw.cpp
       module_dm/density_matrix.cpp
       module_dm/cal_dm_psi.cpp

source/module_elecstate/elecstate.h

@@ -53,6 +53,14 @@ class ElecState
     }
     // virtual void updateRhoK(const psi::Psi<std::complex<double>> &psi) = 0;
     // virtual void updateRhoK(const psi::Psi<double> &psi)=0
+    virtual void cal_tau(const psi::Psi<std::complex<double>>& psi)
+    {
+        return;
+    }
+    virtual void cal_tau(const psi::Psi<double>& psi)
+    {
+        return;
+    }
 
     // update charge density for next scf step
     // in this function, 1. input rho for construct Hamilt and 2. calculated rho from Psi will mix to 3. new charge

source/module_elecstate/elecstate_lcao.cpp

@@ -67,12 +67,7 @@ void ElecStateLCAO<std::complex<double>>::psiToRho(const psi::Psi<std::complex<d
 
     if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
-        for (int is = 0; is < PARAM.inp.nspin; is++)
-        {
-            ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
-        }
-        Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau);
-        this->gint_k->cal_gint(&inout1);
+        this->cal_tau(psi);
     }
 
     this->charge->renormalize_rho();
@@ -124,12 +119,7 @@ void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi)
 
     if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
-        for (int is = 0; is < PARAM.inp.nspin; is++)
-        {
-            ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
-        }
-        Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau);
-        this->gint_gamma->cal_gint(&inout1);
+        this->cal_tau(psi);
     }
 
     this->charge->renormalize_rho();

source/module_elecstate/elecstate_lcao.h

@@ -46,6 +46,7 @@ class ElecStateLCAO : public ElecState
     // virtual void psiToRho(const psi::Psi<double>& psi) override;
     //  return current electronic density rho, as a input for constructing Hamiltonian
     //  const double* getRho(int spin) const override;
+    virtual void cal_tau(const psi::Psi<TK>& psi) override;
 
     // update charge density for next scf step
     // void getNewRho() override;

source/module_elecstate/elecstate_lcao_cal_tau.cpp

@@ -0,0 +1,41 @@
+#include "elecstate_lcao.h"
+
+#include "module_base/timer.h"
+
+namespace elecstate
+{
+
+// calculate the kinetic energy density tau, multi-k case
+template <>
+void ElecStateLCAO<std::complex<double>>::cal_tau(const psi::Psi<std::complex<double>>& psi)
+{
+    ModuleBase::timer::tick("ElecStateLCAO", "cal_tau");
+
+    for (int is = 0; is < PARAM.inp.nspin; is++)
+    {
+        ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
+    }
+    Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau);
+    this->gint_k->cal_gint(&inout1);
+
+    ModuleBase::timer::tick("ElecStateLCAO", "cal_tau");
+    return;
+}
+
+// calculate the kinetic energy density tau, gamma-only case
+template <>
+void ElecStateLCAO<double>::cal_tau(const psi::Psi<double>& psi)
+{
+    ModuleBase::timer::tick("ElecStateLCAO", "cal_tau");
+
+    for (int is = 0; is < PARAM.inp.nspin; is++)
+    {
+        ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
+    }
+    Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau);
+    this->gint_gamma->cal_gint(&inout1);
+
+    ModuleBase::timer::tick("ElecStateLCAO", "cal_tau");
+    return;
+}
+}

source/module_elecstate/elecstate_pw.cpp

@@ -34,7 +34,7 @@ ElecStatePW<T, Device>::~ElecStatePW()
     if (base_device::get_device_type<Device>(this->ctx) == base_device::GpuDevice)
     {
         delmem_var_op()(this->ctx, this->rho_data);
-        if (get_xc_func_type() == 3)
+        if (get_xc_func_type() == 3 || PARAM.inp.out_elf[0] > 0)
         {
             delmem_var_op()(this->ctx, this->kin_r_data);
         }
@@ -53,7 +53,7 @@ void ElecStatePW<T, Device>::init_rho_data()
         for (int ii = 0; ii < this->charge->nspin; ii++) {
             this->rho[ii] = this->rho_data + ii * this->charge->nrxx;
         }
-        if (get_xc_func_type() == 3)
+        if (get_xc_func_type() == 3 || PARAM.inp.out_elf[0] > 0)
         {
             this->kin_r = new Real*[this->charge->nspin];
             resmem_var_op()(this->ctx, this->kin_r_data, this->charge->nspin * this->charge->nrxx);
@@ -64,7 +64,7 @@ void ElecStatePW<T, Device>::init_rho_data()
     }
     else {
         this->rho = reinterpret_cast<Real **>(this->charge->rho);
-        if (get_xc_func_type() == 3)
+        if (get_xc_func_type() == 3 || PARAM.inp.out_elf[0] > 0)
         {
             this->kin_r = reinterpret_cast<Real **>(this->charge->kin_r);
         }

source/module_elecstate/elecstate_pw.h

@@ -33,6 +33,7 @@ class ElecStatePW : public ElecState
     virtual void psiToRho(const psi::Psi<T, Device>& psi);
     // return current electronic density rho, as a input for constructing Hamiltonian
     // const double* getRho(int spin) const override;
+    virtual void cal_tau(const psi::Psi<T, Device>& psi);
 
     // update charge density for next scf step
     // void getNewRho() override;

source/module_elecstate/elecstate_pw_cal_tau.cpp

@@ -0,0 +1,68 @@
+#include "elecstate_pw.h"
+#include "elecstate_getters.h"
+
+namespace elecstate {
+
+template<typename T, typename Device>
+void ElecStatePW<T, Device>::cal_tau(const psi::Psi<T, Device>& psi)
+{
+    ModuleBase::TITLE("ElecStatePW", "cal_tau");
+    for(int is=0; is<PARAM.inp.nspin; is++)
+	{
+        setmem_var_op()(this->ctx, this->kin_r[is], 0,  this->charge->nrxx);
+	}
+
+    for (int ik = 0; ik < psi.get_nk(); ++ik)
+    {
+        psi.fix_k(ik);
+        int npw = psi.get_current_nbas();
+        int current_spin = 0;
+        if (PARAM.inp.nspin == 2)
+        {
+            current_spin = this->klist->isk[ik];
+        }
+        int nbands = psi.get_nbands();
+        for (int ibnd = 0; ibnd < nbands; ibnd++)
+        {
+            this->basis->recip_to_real(this->ctx, &psi(ibnd,0), this->wfcr, ik);
+
+            const auto w1 = static_cast<Real>(this->wg(ik, ibnd) / get_ucell_omega());
+
+            // kinetic energy density
+            for (int j = 0; j < 3; j++)
+            {
+                setmem_complex_op()(this->ctx, this->wfcr, 0,  this->charge->nrxx);
+
+                meta_op()(this->ctx,
+                            ik,
+                            j,
+                            npw,
+                            this->basis->npwk_max,
+                            static_cast<Real>(get_ucell_tpiba()),
+                            this->basis->template get_gcar_data<Real>(),
+                            this->basis->template get_kvec_c_data<Real>(),
+                            &psi(ibnd, 0),
+                            this->wfcr);
+
+                this->basis->recip_to_real(this->ctx, this->wfcr, this->wfcr, ik);
+
+                elecstate_pw_op()(this->ctx, current_spin, this->charge->nrxx, w1, this->kin_r, this->wfcr);
+            }
+        }
+    }
+    if (GlobalV::device_flag == "gpu" || PARAM.inp.precision == "single") {
+        for (int ii = 0; ii < PARAM.inp.nspin; ii++) {
+            castmem_var_d2h_op()(cpu_ctx, this->ctx, this->charge->kin_r[ii], this->kin_r[ii], this->charge->nrxx);
+        }
+    }
+    this->parallelK();
+    ModuleBase::TITLE("ElecStatePW", "cal_tau");
+}
+
+template class ElecStatePW<std::complex<float>, base_device::DEVICE_CPU>;
+template class ElecStatePW<std::complex<double>, base_device::DEVICE_CPU>;
+#if ((defined __CUDA) || (defined __ROCM))
+template class ElecStatePW<std::complex<float>, base_device::DEVICE_GPU>;
+template class ElecStatePW<std::complex<double>, base_device::DEVICE_GPU>;
+#endif 
+} // namespace elecstate

source/module_elecstate/module_charge/charge.cpp

@@ -80,7 +80,7 @@ void Charge::destroy()
         delete[] _space_rhog_save;
         delete[] _space_kin_r;
         delete[] _space_kin_r_save;
-        if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5)
+        if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0)
         {
             delete[] kin_r;
             delete[] kin_r_save;
@@ -121,7 +121,7 @@ void Charge::allocate(const int& nspin_in)
     _space_rho_save = new double[nspin * nrxx];
     _space_rhog = new std::complex<double>[nspin * ngmc];
     _space_rhog_save = new std::complex<double>[nspin * ngmc];
-    if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5)
+    if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0)
     {
         _space_kin_r = new double[nspin * nrxx];
         _space_kin_r_save = new double[nspin * nrxx];
@@ -130,7 +130,7 @@ void Charge::allocate(const int& nspin_in)
     rhog = new std::complex<double>*[nspin];
     rho_save = new double*[nspin];
     rhog_save = new std::complex<double>*[nspin];
-    if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5)
+    if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0)
     {
         kin_r = new double*[nspin];
         kin_r_save = new double*[nspin];
@@ -151,7 +151,7 @@ void Charge::allocate(const int& nspin_in)
         ModuleBase::GlobalFunc::ZEROS(rhog[is], ngmc);
         ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx);
         ModuleBase::GlobalFunc::ZEROS(rhog_save[is], ngmc);
-        if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5)
+        if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0)
         {
             kin_r[is] = _space_kin_r + is * nrxx;
             ModuleBase::GlobalFunc::ZEROS(kin_r[is], nrxx);
@@ -171,7 +171,7 @@ void Charge::allocate(const int& nspin_in)
     ModuleBase::Memory::record("Chg::rho_save", sizeof(double) * nspin * nrxx);
     ModuleBase::Memory::record("Chg::rhog", sizeof(double) * nspin * ngmc);
     ModuleBase::Memory::record("Chg::rhog_save", sizeof(double) * nspin * ngmc);
-    if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5)
+    if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0)
     {
         ModuleBase::Memory::record("Chg::kin_r", sizeof(double) * nspin * ngmc);
         ModuleBase::Memory::record("Chg::kin_r_save", sizeof(double) * nspin * ngmc);

source/module_elecstate/module_charge/charge.h

@@ -138,6 +138,7 @@ class Charge
     int ngmc; // number of g vectors in this processor
     int nspin; // number of spins
     ModulePW::PW_Basis* rhopw = nullptr;// When double_grid is used, rhopw = rhodpw (dense grid)
+    bool cal_elf = false; // whether to calculate electron localization function (ELF)
   private:
 
     void destroy();    // free arrays  liuyu 2023-03-12

source/module_elecstate/module_charge/charge_mpi.cpp

@@ -132,7 +132,7 @@ void Charge::rho_mpi()
     for (int is = 0; is < PARAM.inp.nspin; ++is)
     {
         reduce_diff_pools(this->rho[is]);
-        if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5)
+        if (elecstate::get_xc_func_type() == 3 || elecstate::get_xc_func_type() == 5 || PARAM.inp.out_elf[0] > 0)
         {
             reduce_diff_pools(this->kin_r[is]);
         }

source/module_elecstate/module_charge/symmetry_rho.cpp

@@ -33,7 +33,7 @@ void Symmetry_rho::begin(const int& spin_now,
         psymmg(CHR.rhog[spin_now], rho_basis, symm); // need to modify
         rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]);
 
-        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5 || CHR.cal_elf)
         {
             // Use std::vector to manage kin_g instead of raw pointer
             std::vector<std::complex<double>> kin_g(CHR.ngmc);

source/module_elecstate/test/CMakeLists.txt

@@ -53,6 +53,7 @@ AddTest(
   LIBS parameter  ${math_libs} planewave_serial base device
   SOURCES elecstate_pw_test.cpp
     ../elecstate_pw.cpp
+    ../elecstate_pw_cal_tau.cpp
     ../elecstate.cpp
     ../occupy.cpp
     ../../module_psi/psi.cpp