kluophysics
diff --git a/‎examples/hse/pw_Si2/INPUT
+32 b/‎examples/hse/pw_Si2/INPUT
+32
diff --git a/‎examples/hse/pw_Si2/KPT
+4 b/‎examples/hse/pw_Si2/KPT
+4
diff --git a/‎examples/hse/pw_Si2/STRU
+24 b/‎examples/hse/pw_Si2/STRU
+24
diff --git a/‎source/Makefile.Objects
+4 b/‎source/Makefile.Objects
+4
diff --git a/‎source/module_elecstate/elecstate.h
-4 b/‎source/module_elecstate/elecstate.h
-4
diff --git a/‎source/module_elecstate/elecstate_exx.cpp
-4 b/‎source/module_elecstate/elecstate_exx.cpp
-4
diff --git a/‎source/module_elecstate/potentials/potential_new.h
+6 b/‎source/module_elecstate/potentials/potential_new.h
+6
diff --git a/‎source/module_esolver/CMakeLists.txt
+2 b/‎source/module_esolver/CMakeLists.txt
+2
diff --git a/‎source/module_esolver/esolver_ks.cpp
+3 b/‎source/module_esolver/esolver_ks.cpp
+3
diff --git a/‎source/module_esolver/esolver_ks_pw.cpp
+77 b/‎source/module_esolver/esolver_ks_pw.cpp
+77
diff --git a/‎source/module_esolver/esolver_ks_pw.h
+4 b/‎source/module_esolver/esolver_ks_pw.h
+4
diff --git a/‎source/module_hamilt_general/module_xc/xc_functional.cpp
+15-6 b/‎source/module_hamilt_general/module_xc/xc_functional.cpp
+15-6
diff --git a/‎source/module_hamilt_general/operator.h
+6-1 b/‎source/module_hamilt_general/operator.h
+6-1
diff --git a/‎source/module_hamilt_lcao/module_tddft/propagator.h
+2 b/‎source/module_hamilt_lcao/module_tddft/propagator.h
+2
@@ -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          hse
+
+# 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
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+3 3 3 0 0 0
@@ -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
@@ -50,6 +50,7 @@ VPATH=./src_global:\
 ./module_hamilt_pw/hamilt_stodft:\
 ./module_hamilt_pw/hamilt_pwdft/operator_pw:\
 ./module_hamilt_pw/hamilt_pwdft/kernels:\
+./module_hamilt_pw/hamilt_pwdft/module_exx_helper:\
 ./module_hamilt_pw/hamilt_stodft/kernels:\
 ./module_hamilt_lcao/module_hcontainer:\
 ./module_hamilt_lcao/hamilt_lcaodft:\
@@ -295,6 +296,7 @@ OBJS_HAMILT=hamilt_pw.o\
     hamilt_sdft_pw.o\
     operator.o\
     operator_pw.o\
+    op_exx_pw.o\
     ekinetic_pw.o\
     ekinetic_op.o\
     hpsi_norm_op.o\
@@ -306,6 +308,7 @@ OBJS_HAMILT=hamilt_pw.o\
     meta_op.o\
     velocity_pw.o\
     radial_proj.o\
+    exx_helper.o\
 
 OBJS_HAMILT_OF=kedf_tf.o\
     kedf_vw.o\
@@ -671,6 +674,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     sto_stress_pw.o\
     stress_func_cc.o\
     stress_func_ewa.o\
+    stress_func_exx.o\
     stress_func_gga.o\
     stress_func_mgga.o\
     stress_func_har.o\
 
@@ -131,12 +131,8 @@ class ElecState
     bool vnew_exist = false;
     void cal_converged();
     void cal_energies(const int type);
-#ifdef __EXX
-#ifdef __LCAO
     void set_exx(const double& Eexx);
     void set_exx(const std::complex<double>& Eexx);
-#endif //__LCAO
-#endif //__EXX
 
     double get_hartree_energy();
     double get_etot_efield();
 
@@ -3,8 +3,6 @@
 namespace elecstate
 {
 
-#ifdef __EXX
-#ifdef __LCAO
 /// @brief calculation if converged
 /// @date Peize Lin add 2016-12-03
 void ElecState::set_exx(const double& Eexx)
@@ -17,7 +15,5 @@ void ElecState::set_exx(const double& Eexx)
     }
     return;
 }
-#endif //__LCAO
-#endif //__EXX
 
 }
@@ -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;
 
@@ -28,6 +28,8 @@ add_library(
     esolver
     OBJECT
     ${objects}
+        ../module_hamilt_pw/hamilt_pwdft/module_exx_helper/exx_helper.cpp
+        ../module_hamilt_pw/hamilt_pwdft/module_exx_helper/exx_helper.h
 )
 
 if(ENABLE_COVERAGE)
 
@@ -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"
@@ -24,6 +25,8 @@
 #include "module_cell/module_paw/paw_cell.h"
 #endif
 
+#include "esolver_ks_pw.h"
+
 namespace ModuleESolver
 {
 
 
@@ -225,6 +225,26 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
                                  this->pelec->wg,
                                  this->pelec->skip_weights);
     }
+
+
+    // 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.set_firstiter();
+        }
+
+        if (GlobalC::exx_info.info_global.cal_exx)
+        {
+            exx_helper.set_wg(&this->pelec->wg);
+        }
+    }
+
 }
 
 template <typename T, typename Device>
@@ -258,6 +278,19 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     // allocate HamiltPW
     this->allocate_hamilt(ucell);
 
+    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);
+        }
+
+    }
+
     //----------------------------------------------------------
     //! calculate the total local pseudopotential in real space
     //----------------------------------------------------------
@@ -333,6 +366,18 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
         this->already_initpsi = true;
     }
 
+    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")
+        {
+            exx_helper.set_psi(kspw_psi);
+        }
+
+    }
+
     ModuleBase::timer::tick("ESolver_KS_PW", "before_scf");
 }
 
@@ -504,6 +549,11 @@ void ESolver_KS_PW<T, Device>::update_pot(UnitCell& ucell, const int istep, cons
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver)
 {
+    if (GlobalC::exx_info.info_global.cal_exx && !exx_helper.op_exx->first_iter)
+    {
+        this->pelec->set_exx(exx_helper.cal_exx_energy(kspw_psi));
+    }
+
     // deband is calculated from "output" charge density calculated
     // in sum_band
     // need 'rho(out)' and 'vr (v_h(in) and v_xc(in))'
@@ -522,6 +572,33 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
         this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell);
     }
 
+    if (GlobalC::exx_info.info_global.cal_exx)
+    {
+        if (GlobalC::exx_info.info_global.separate_loop)
+        {
+            if (conv_esolver)
+            {
+                exx_helper.set_psi(this->kspw_psi);
+
+                conv_esolver = exx_helper.exx_after_converge(iter);
+
+                if (!conv_esolver)
+                {
+                    std::cout << " Setting Psi for EXX PW Inner Loop" << std::endl;
+                    exx_helper.op_exx->first_iter = false;
+                    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+                    update_pot(ucell, istep, iter, conv_esolver);
+                }
+            }
+        }
+        else
+        {
+//            std::cout << "setting psi for each iter" << std::endl;
+            exx_helper.set_psi(this->kspw_psi);
+        }
+
+    }
+
     // 3) Print out electronic wavefunctions in pw basis
     if (PARAM.inp.out_wfc_pw == 1 || PARAM.inp.out_wfc_pw == 2)
     {
 
@@ -3,6 +3,8 @@
 #include "./esolver_ks.h"
 #include "module_hamilt_pw/hamilt_pwdft/operator_pw/velocity_pw.h"
 #include "module_psi/psi_init.h"
+#include "module_hamilt_pw/hamilt_pwdft/module_exx_helper/exx_helper.h"
+#include "module_hamilt_pw/hamilt_pwdft/global.h"
 
 #include <memory>
 #include <module_base/macros.h>
@@ -31,6 +33,8 @@ class ESolver_KS_PW : public ESolver_KS<T, Device>
 
     void after_all_runners(UnitCell& ucell) override;
 
+    Exx_Helper<T, Device> exx_helper;
+
   protected:
     virtual void before_scf(UnitCell& ucell, const int istep) override;
 
 
@@ -40,6 +40,9 @@ method. */
     else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") {
         XC_Functional::set_xc_type("scan");
     }
+    else if (ucell.atoms[0].ncpp.xc_func == "B3LYP") {
+        XC_Functional::set_xc_type("blyp");
+    }
 }
 
 // The setting values of functional id according to the index in LIBXC
@@ -244,6 +247,12 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
         func_type = 2;
         use_libxc = true;
     }
+    else if (xc_func == "B3LYP")
+    {
+        func_id.push_back(XC_HYB_GGA_XC_B3LYP);
+        func_type = 4;
+        use_libxc = true;
+    }
 #endif
     else
     {
@@ -268,19 +277,19 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
         std::cerr << "\n OPTX untested please test,";
     }
 
-    if((func_type == 4 || func_type == 5) && PARAM.inp.basis_type == "pw")
-    {
-        ModuleBase::WARNING_QUIT("set_xc_type","hybrid functional not realized for planewave yet");
-    }
+    // if((func_type == 4 || func_type == 5) && PARAM.inp.basis_type == "pw")
+    // {
+    //     ModuleBase::WARNING_QUIT("set_xc_type","hybrid functional not realized for planewave yet");
+    // }
     if((func_type == 3 || func_type == 5) && PARAM.inp.nspin==4)
     {
         ModuleBase::WARNING_QUIT("set_xc_type","meta-GGA has not been implemented for nspin = 4 yet");
     }
 
 #ifndef __EXX
-    if(func_type == 4 || func_type == 5)
+    if((func_type == 4 || func_type == 5) && PARAM.inp.basis_type == "lcao")
     {
-        ModuleBase::WARNING_QUIT("set_xc_type","compile with libri to use hybrid functional");
+        ModuleBase::WARNING_QUIT("set_xc_type","compile with libri to use hybrid functional in lcao basis");
     }
 #endif
 
 
@@ -18,6 +18,7 @@ enum class calculation_type
     pw_veff,
     pw_meta,
     pw_onsite,
+    pw_exx,
     lcao_overlap,
     lcao_fixed,
     lcao_gint,
@@ -41,7 +42,6 @@ class Operator
 
     // this is the core function for Operator
     //  do H|psi> from input |psi> ,
-
     /// as default, different operators donate hPsi independently
     /// run this->act function for the first operator and run all act() for other nodes in chain table
     /// if this procedure is not suitable for your operator, just override this function.
@@ -90,6 +90,11 @@ class Operator
         return this->act_type;
     }
 
+    calculation_type get_cal_type() const
+    {
+        return this->cal_type;
+    }
+
   protected:
     int ik = 0;
     int act_type = 1; ///< determine which act() interface would be called in hPsi()
 
@@ -84,6 +84,7 @@ ct::Tensor create_identity_matrix(const int n, ct::DeviceType device = ct::Devic
             data_ptr[i * n + i] = init_value<T>();
         }
     }
+#if ((defined __CUDA))
     else if (device == ct::DeviceType::GpuDevice)
     {
         // For GPU, we need to use a kernel to set the diagonal elements
@@ -94,6 +95,7 @@ ct::Tensor create_identity_matrix(const int n, ct::DeviceType device = ct::Devic
             ct::kernels::set_memory<T, ct::DEVICE_GPU>()(data_ptr + i * n + i, value, 1);
         }
     }
+#endif
 
     return tensor;
 }
Original file line number	Diff line number	Diff line change
`@@ -3,8 +3,6 @@`
`3`	`3`	`namespace elecstate`
`4`	`4`	`{`
`5`	`5`
`6`		`-#ifdef __EXX`
`7`		`-#ifdef __LCAO`
`8`	`6`	`/// @brief calculation if converged`
`9`	`7`	`/// @date Peize Lin add 2016-12-03`
`10`	`8`	`void ElecState::set_exx(const double& Eexx)`
`@@ -17,7 +15,5 @@ void ElecState::set_exx(const double& Eexx)`
`17`	`15`	`}`
`18`	`16`	`return;`
`19`	`17`	`}`
`20`		`-#endif //__LCAO`
`21`		`-#endif //__EXX`
`22`	`18`
`23`	`19`	`}`
Original file line number	Diff line number	Diff line change
`@@ -170,6 +170,12 @@ class Potential : public PotBase`
`170`	`170`	`{`
`171`	`171`	`return this->v_effective_fixed.data();`
`172`	`172`	`}`
	`173`	`+ const ModulePW::PW_Basis *get_rho_basis() const`
	`174`	`+ {`
	`175`	`+ return this->rho_basis_;`
	`176`	`+ }`
	`177`	`+ // What about adding a function to get the wfc?`
	`178`	`+ // This is useful for the calculation of the exx energy`
`173`	`179`
`174`	`180`
`175`	`181`	`/// @brief get the value of vloc at G=0;`
Original file line number	Diff line number	Diff line change
`@@ -28,6 +28,8 @@ add_library(`
`28`	`28`	`esolver`
`29`	`29`	`OBJECT`
`30`	`30`	`${objects}`
	`31`	`+ ../module_hamilt_pw/hamilt_pwdft/module_exx_helper/exx_helper.cpp`
	`32`	`+ ../module_hamilt_pw/hamilt_pwdft/module_exx_helper/exx_helper.h`
`31`	`33`	`)`
`32`	`34`
`33`	`35`	`if(ENABLE_COVERAGE)`
Original file line number	Diff line number	Diff line change
`@@ -40,6 +40,9 @@ method. */`
`40`	`40`	`else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") {`
`41`	`41`	`XC_Functional::set_xc_type("scan");`
`42`	`42`	`}`
	`43`	`+ else if (ucell.atoms[0].ncpp.xc_func == "B3LYP") {`
	`44`	`+ XC_Functional::set_xc_type("blyp");`
	`45`	`+ }`
`43`	`46`	`}`
`44`	`47`
`45`	`48`	`// The setting values of functional id according to the index in LIBXC`
`@@ -244,6 +247,12 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)`
`244`	`247`	`func_type = 2;`
`245`	`248`	`use_libxc = true;`
`246`	`249`	`}`
	`250`	`+ else if (xc_func == "B3LYP")`
	`251`	`+ {`
	`252`	`+ func_id.push_back(XC_HYB_GGA_XC_B3LYP);`
	`253`	`+ func_type = 4;`
	`254`	`+ use_libxc = true;`
	`255`	`+ }`
`247`	`256`	`#endif`
`248`	`257`	`else`
`249`	`258`	`{`
`@@ -268,19 +277,19 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)`
`268`	`277`	`std::cerr << "\n OPTX untested please test,";`
`269`	`278`	`}`
`270`	`279`
`271`		`- if((func_type == 4 \|\| func_type == 5) && PARAM.inp.basis_type == "pw")`
`272`		`- {`
`273`		`- ModuleBase::WARNING_QUIT("set_xc_type","hybrid functional not realized for planewave yet");`
`274`		`- }`
	`280`	`+ // if((func_type == 4 \|\| func_type == 5) && PARAM.inp.basis_type == "pw")`
	`281`	`+ // {`
	`282`	`+ // ModuleBase::WARNING_QUIT("set_xc_type","hybrid functional not realized for planewave yet");`
	`283`	`+ // }`
`275`	`284`	`if((func_type == 3 \|\| func_type == 5) && PARAM.inp.nspin==4)`
`276`	`285`	`{`
`277`	`286`	`ModuleBase::WARNING_QUIT("set_xc_type","meta-GGA has not been implemented for nspin = 4 yet");`
`278`	`287`	`}`
`279`	`288`
`280`	`289`	`#ifndef __EXX`
`281`		`- if(func_type == 4 \|\| func_type == 5)`
	`290`	`+ if((func_type == 4 \|\| func_type == 5) && PARAM.inp.basis_type == "lcao")`
`282`	`291`	`{`
`283`		`- ModuleBase::WARNING_QUIT("set_xc_type","compile with libri to use hybrid functional");`
	`292`	`+ ModuleBase::WARNING_QUIT("set_xc_type","compile with libri to use hybrid functional in lcao basis");`
`284`	`293`	`}`
`285`	`294`	`#endif`
`286`	`295`
Original file line number	Diff line number	Diff line change
`@@ -84,6 +84,7 @@ ct::Tensor create_identity_matrix(const int n, ct::DeviceType device = ct::Devic`
`84`	`84`	`data_ptr[i * n + i] = init_value<T>();`
`85`	`85`	`}`
`86`	`86`	`}`
	`87`	`+#if ((defined __CUDA))`
`87`	`88`	`else if (device == ct::DeviceType::GpuDevice)`
`88`	`89`	`{`
`89`	`90`	`// For GPU, we need to use a kernel to set the diagonal elements`
`@@ -94,6 +95,7 @@ ct::Tensor create_identity_matrix(const int n, ct::DeviceType device = ct::Devic`
`94`	`95`	`ct::kernels::set_memory<T, ct::DEVICE_GPU>()(data_ptr + i * n + i, value, 1);`
`95`	`96`	`}`
`96`	`97`	`}`
	`98`	`+#endif`
`97`	`99`
`98`	`100`	`return tensor;`
`99`	`101`	`}`