deepmodeling
diff --git a/‎examples/hse/pw_Si2/INPUT‎
Lines changed: 32 additions & 0 deletions b/‎examples/hse/pw_Si2/INPUT‎
Lines changed: 32 additions & 0 deletions
diff --git a/‎examples/hse/pw_Si2/KPT‎
Lines changed: 4 additions & 0 deletions b/‎examples/hse/pw_Si2/KPT‎
Lines changed: 4 additions & 0 deletions
diff --git a/‎examples/hse/pw_Si2/STRU‎
Lines changed: 24 additions & 0 deletions b/‎examples/hse/pw_Si2/STRU‎
Lines changed: 24 additions & 0 deletions
diff --git a/‎source/Makefile.Objects‎
Lines changed: 4 additions & 0 deletions b/‎source/Makefile.Objects‎
Lines changed: 4 additions & 0 deletions
diff --git a/‎source/module_elecstate/elecstate.h‎
Lines changed: 0 additions & 4 deletions b/‎source/module_elecstate/elecstate.h‎
Lines changed: 0 additions & 4 deletions
diff --git a/‎source/module_elecstate/elecstate_exx.cpp‎
Lines changed: 0 additions & 4 deletions b/‎source/module_elecstate/elecstate_exx.cpp‎
Lines changed: 0 additions & 4 deletions
diff --git a/‎source/module_elecstate/potentials/potential_new.h‎
Lines changed: 6 additions & 0 deletions b/‎source/module_elecstate/potentials/potential_new.h‎
Lines changed: 6 additions & 0 deletions
diff --git a/‎source/module_esolver/CMakeLists.txt‎
Lines changed: 2 additions & 0 deletions b/‎source/module_esolver/CMakeLists.txt‎
Lines changed: 2 additions & 0 deletions
diff --git a/‎source/module_esolver/esolver_ks.cpp‎
Lines changed: 3 additions & 0 deletions b/‎source/module_esolver/esolver_ks.cpp‎
Lines changed: 3 additions & 0 deletions
diff --git a/‎source/module_esolver/esolver_ks_pw.cpp‎
Lines changed: 77 additions & 0 deletions b/‎source/module_esolver/esolver_ks_pw.cpp‎
Lines changed: 77 additions & 0 deletions
@@ -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;
 
@@ -27,6 +27,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)
     {
-Original file line number
+Diff line change
@@ @@ -0,0 +1,4 @@ @@
 +K_POINTS
 +0
 +Gamma
 +3 3 3 0 0 0
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
`@@ -27,6 +27,8 @@ add_library(`
`27`	`27`	`esolver`
`28`	`28`	`OBJECT`
`29`	`29`	`${objects}`
	`30`	`+ ../module_hamilt_pw/hamilt_pwdft/module_exx_helper/exx_helper.cpp`
	`31`	`+ ../module_hamilt_pw/hamilt_pwdft/module_exx_helper/exx_helper.h`
`30`	`32`	`)`
`31`	`33`
`32`	`34`	`if(ENABLE_COVERAGE)`