Feature: Hybrid functionals for planewave basis

Dockerfile.gnu

@@ -10,7 +10,7 @@
 FROM ubuntu:22.04
 RUN apt update && apt install -y --no-install-recommends \
     libopenblas-openmp-dev liblapack-dev libscalapack-mpi-dev libelpa-dev libfftw3-dev libcereal-dev \
-    libxc-dev libgtest-dev libgmock-dev libbenchmark-dev python3-numpy \
+    libxc-dev libgtest-dev libgmock-dev libbenchmark-dev python3-numpy libparmetis-dev libsuperlu-dist-dev \
     bc cmake git g++ make bc time sudo unzip vim wget gfortran
     # If you wish to use the LLVM compiler, replace 'g++' above with 'clang libomp-dev'.
 
@@ -26,6 +26,15 @@ RUN wget https://download.pytorch.org/libtorch/cpu/libtorch-cxx11-abi-shared-wit
         --no-check-certificate --quiet -O libtorch.zip && \
     unzip -q libtorch.zip -d /opt && rm libtorch.zip
 
+RUN cp -r /usr/include/superlu-dist/* /usr/include && \
+    wget https://bitbucket.org/berkeleylab/pexsi/downloads/pexsi_v2.0.0.tar.gz --no-check-certificate --quiet && \
+    tar -xzf pexsi_v2.0.0.tar.gz && \
+    cd pexsi_v2.0.0 && \
+    cmake -B build -DPEXSI_ENABLE_FORTRAN=OFF -DPEXSI_ENABLE_OPENMP=ON && \
+    cmake --build build -j`nproc` && \
+    cmake --install build && \
+    cd .. && rm -rf pexsi_v2.0.0
+
 ENV CMAKE_PREFIX_PATH=/opt/libtorch/share/cmake
 
 ADD https://api.github.com/repos/deepmodeling/abacus-develop/git/refs/heads/develop /dev/null
@@ -34,7 +43,7 @@ ADD https://api.github.com/repos/deepmodeling/abacus-develop/git/refs/heads/deve
 
 RUN git clone https://github.com/deepmodeling/abacus-develop.git --depth 1 && \
     cd abacus-develop && \
-    cmake -B build -DENABLE_DEEPKS=ON -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON -DENABLE_RAPIDJSON=ON && \
+    cmake -B build -DENABLE_DEEPKS=ON -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON -DENABLE_RAPIDJSON=ON -DENABLE_PEXSI=ON && \
     cmake --build build -j`nproc` && \
     cmake --install build && \
     rm -rf build

examples/hse/pw_Si2/INPUT

@@ -0,0 +1,32 @@
+INPUT_PARAMETERS
+pseudo_dir              ../../../tests/PP_ORB
+orbital_dir             ../../../tests/PP_ORB
+nbands                  4
+nspin                   1
+calculation             scf
+basis_type              pw
+ks_solver               dav
+ecutwfc                 50
+scf_thr                 1e-9
+scf_nmax                100
+gamma_only              0
+symmetry                -1
+smearing_method         fixed
+mixing_type             broyden
+mixing_beta             0.7
+
+dft_functional          pbe0
+
+# init_wfc file
+# init_chg file
+
+pseudo_mesh            1
+pseudo_rcut            10
+
+# out_wfc_pw 1
+# out_chg 1
+
+exx_hybrid_alpha        0.25
+
+cal_stress 1
+# cal_force 1

examples/hse/pw_Si2/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+3 3 3 0 0 0

examples/hse/pw_Si2/STRU

@@ -0,0 +1,24 @@
+ATOMIC_SPECIES
+Si	1	Si_ONCV_PBE_FR-1.1.upf					#Pseudopotentials are downloaded from http://www.quantum-simulation.org/potentials/sg15_oncv/upf/
+
+NUMERICAL_ORBITAL
+orb_Si.dat
+
+LATTICE_CONSTANT
+1.889766
+
+LATTICE_VECTORS
+0.0	2.708337	2.708337
+2.708337	0.0	2.708337
+2.708337	2.708337	0.0
+
+
+
+ATOMIC_POSITIONS
+Direct
+
+Si
+0
+2
+0.125	0.125	0.125	0 0 0
+0.875	0.875	0.875	0 0 0

source/module_elecstate/CMakeLists.txt

@@ -8,6 +8,7 @@ list(APPEND objects
     elecstate_pw_sdft.cpp
     elecstate_pw_cal_tau.cpp
     potentials/gatefield.cpp
+    # potentials/e_exx.cpp
     potentials/efield.cpp
     potentials/H_Hartree_pw.cpp
     potentials/pot_xc.cpp

source/module_elecstate/module_charge/charge_init.cpp

@@ -53,6 +53,42 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             {
                 rhopw->recip2real(rhog[is], rho[is]);
             }
+
+//            // also read from abacus type file, and compare with the binary file
+//            for (int is = 0; is < PARAM.inp.nspin; ++is)
+//            {
+//                std::stringstream ssc;
+//                ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_CHG.cube";
+//                double& ef_tmp = eferm_iout.get_ef(is);
+//                if (ModuleIO::read_rho(
+//#ifdef __MPI
+//                        &(GlobalC::Pgrid),
+//#endif
+//                        GlobalV::MY_RANK,
+//                        PARAM.inp.esolver_type,
+//                        GlobalV::RANK_IN_STOGROUP,
+//                        is,
+//                        GlobalV::ofs_running,
+//                        PARAM.inp.nspin,
+//                        ssc.str(),
+//                        this->rho_save[is],
+//                        this->rhopw->nx,
+//                        this->rhopw->ny,
+//                        this->rhopw->nz,
+//                        ef_tmp,
+//                        &(GlobalC::ucell),
+//                        this->prenspin))
+//                {
+//                    GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl;
+//                }
+//            Todo: Fix Comparison
+
+//            }
+
+//            // compare the two charge density
+//            Charge_Mixing chm;
+//            chm.set_rhopw(rhopw, rhopw);
+//            std::cout << chm.get_drho(this, GlobalV::nelec) << std::endl;
         }
         else
         {
@@ -219,7 +255,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             }
             catch (const std::exception& e)
             {
-                // try to load from the output of `out_chg` 
+                // try to load from the output of `out_chg`
                 std::stringstream ssc;
                 ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_CHG.cube";
                 if (ModuleIO::read_vdata_palgrid(pgrid,
@@ -256,7 +292,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
 // computes the core charge on the real space 3D mesh.
 //==========================================================
 void Charge::set_rho_core(const UnitCell& ucell,
-                          const ModuleBase::ComplexMatrix& structure_factor, 
+                          const ModuleBase::ComplexMatrix& structure_factor,
                           const bool* numeric)
 {
     ModuleBase::TITLE("Charge","set_rho_core");

source/module_elecstate/potentials/e_exx_pw.cpp

@@ -0,0 +1,26 @@
+#include "e_exx_pw.h"
+
+namespace elecstate
+{
+
+double E_exx::exx_energy = 0.0;
+
+E_exx::E_exx()
+{
+}
+
+E_exx::~E_exx()
+{
+}
+
+double E_exx::cal_exx(const ModulePW::PW_Basis_K* wfcpw,
+                      const ModulePW::PW_Basis* rhopw,
+                      const int* isk,
+                      const double* const* const rho,
+                      const double& Eexx)
+{
+    double exx = 0.0;
+    return exx;
+}
+
+} // namespace elecstate

source/module_elecstate/potentials/e_exx_pw.h

@@ -0,0 +1,31 @@
+#ifndef E_EXX_PW_H
+#define E_EXX_PW_H
+
+#include "module_basis/module_pw/pw_basis.h"
+#include "module_basis/module_pw/pw_basis_k.h"
+
+namespace elecstate
+{
+
+class E_exx
+{
+  public:
+    E_exx();
+    ~E_exx();
+
+    static double cal_exx(const ModulePW::PW_Basis_K* wfcpw,
+                        const ModulePW::PW_Basis* rhopw,
+                        const int* isk,
+                        const double* const* const rho,
+                        const double& Eexx);
+
+    static double get_exx_energy() { return exx_energy; }
+
+  private:
+    static double exx_energy;
+
+};
+
+} // namespace elecstate
+
+#endif // E_EXX_PW_H

source/module_elecstate/potentials/potential_new.h

@@ -170,6 +170,12 @@ class Potential : public PotBase
     {
         return this->v_effective_fixed.data();
     }
+    const ModulePW::PW_Basis *get_rho_basis() const
+    {
+        return this->rho_basis_;
+    }
+    // What about adding a function to get the wfc? 
+    // This is useful for the calculation of the exx energy
 
 
     /// @brief get the value of vloc at G=0;

source/module_esolver/CMakeLists.txt

@@ -26,6 +26,7 @@ add_library(
     esolver
     OBJECT
     ${objects}
+        module_exx_helper/exx_helper.cpp
 )
 
 if(ENABLE_COVERAGE)

source/module_esolver/esolver_ks.cpp

@@ -2,6 +2,7 @@
 
 #include "module_base/timer.h"
 #include "module_cell/cal_atoms_info.h"
+#include "module_hamilt_general/module_xc/xc_functional.h"
 #include "module_io/cube_io.h"
 #include "module_io/json_output/init_info.h"
 #include "module_io/json_output/output_info.h"
@@ -22,6 +23,8 @@
 #include "module_cell/module_paw/paw_cell.h"
 #endif
 
+#include "esolver_ks_pw.h"
+
 namespace ModuleESolver
 {
 
@@ -649,6 +652,22 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
     this->pelec->cal_energies(1);
     this->pelec->cal_energies(2);
 
+    // for separate loop in hybrid functionals in pw
+    // this gives correct energy for the first "pure" functional loop
+    // and update the functional for the next loop
+    if (PARAM.inp.basis_type == "pw")
+    {
+        auto p_esolver_ks_pw = dynamic_cast<ESolver_KS_PW<T, Device>*>(this);
+
+        if (this->conv_esolver && GlobalC::exx_info.info_global.cal_exx && GlobalC::exx_info.info_global.separate_loop && p_esolver_ks_pw->exx_helper.first_iter)
+        {
+            this->conv_esolver = false;
+            XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+            this->update_pot(ucell, istep, iter);
+            this->conv_esolver = true;
+        }
+    }
+
     if (iter == 1)
     {
         this->pelec->f_en.etot_old = this->pelec->f_en.etot;

source/module_esolver/esolver_ks_pw.cpp

@@ -230,6 +230,30 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
     {
         this->pelec->fixed_weights(PARAM.inp.ocp_kb, PARAM.globalv.nbands_l, PARAM.inp.nelec);
     }
+
+
+    // EXX Todo: Fix the control flow
+#ifdef __EXX
+    // 10) initialize exx pw
+    if (PARAM.inp.calculation == "scf"
+        || PARAM.inp.calculation == "relax"
+        || PARAM.inp.calculation == "cell-relax"
+        || PARAM.inp.calculation == "md")
+    {
+        if (GlobalC::exx_info.info_global.cal_exx && GlobalC::exx_info.info_global.separate_loop == true)
+        {
+            XC_Functional::set_xc_first_loop(ucell);
+            exx_helper.first_iter = true;
+        }
+
+        if (GlobalC::exx_info.info_global.cal_exx)
+        {
+            exx_helper.wf_wg = &(this->pelec->wg);
+            exx_helper.set_psi(this->kspw_psi[0]);
+        }
+    }
+#endif
+
 }
 
 template <typename T, typename Device>
@@ -273,6 +297,22 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     // allocate HamiltPW
     this->allocate_hamilt(ucell);
 
+#ifdef __EXX
+    if (PARAM.inp.calculation == "scf"
+        || PARAM.inp.calculation == "relax"
+        || PARAM.inp.calculation == "cell-relax"
+        || PARAM.inp.calculation == "md")
+    {
+        if (GlobalC::exx_info.info_global.cal_exx && PARAM.inp.basis_type == "pw")
+        {
+            auto hamilt_pw = reinterpret_cast<hamilt::HamiltPW<T, Device>*>(this->p_hamilt);
+            hamilt_pw->set_exx_helper(&exx_helper);
+        }
+
+    }
+#endif
+
+
     //----------------------------------------------------------
     // about vdw, jiyy add vdwd3 and linpz add vdwd2
     //----------------------------------------------------------
@@ -545,6 +585,13 @@ void ESolver_KS_PW<T, Device>::hamilt2density_single(UnitCell& ucell,
         srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm);
     }
 
+#ifdef __EXX
+    if (GlobalC::exx_info.info_global.cal_exx && !exx_helper.first_iter)
+    {
+        this->pelec->set_exx(exx_helper.cal_exx_energy(this->kspw_psi[0], this));
+    }
+#endif
+
     // deband is calculated from "output" charge density calculated
     // in sum_band
     // need 'rho(out)' and 'vr (v_h(in) and v_xc(in))'
@@ -589,6 +636,31 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
         this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell);
     }
 
+#ifdef __EXX
+    if (GlobalC::exx_info.info_global.cal_exx)
+    {
+        if (GlobalC::exx_info.info_global.separate_loop)
+        {
+            if (this->conv_esolver)
+            {
+                std::cout << " setting psi for exx inner loop" << std::endl;
+                exx_helper.set_psi(this->kspw_psi[0]);
+
+                this->conv_esolver = exx_helper.exx_after_converge(iter);
+                exx_helper.first_iter = false;
+
+                XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+            }
+        }
+        else
+        {
+//            std::cout << "setting psi for each iter" << std::endl;
+            exx_helper.set_psi(this->kspw_psi[0]);
+        }
+
+    }
+#endif // __EXX
+
     // 4) Print out electronic wavefunctions
     if (PARAM.inp.out_wfc_pw == 1 || PARAM.inp.out_wfc_pw == 2)
     {
@@ -721,6 +793,9 @@ void ESolver_KS_PW<T, Device>::after_scf(UnitCell& ucell, const int istep)
                                   this->pelec->wg);
     }
 
+    // 10) reset exx_helper
+    this->exx_helper.reset_div();
+
     ModuleBase::timer::tick("ESolver_KS_PW", "after_scf");
 }