refactor Exx_Opt_Orb

docs/advanced/input_files/input-main.md

@@ -512,6 +512,7 @@ These variables are used to control general system parameters.
   - get_wf: obtain real space wave functions (for LCAO basis only). See `out_wfc_norm` and `out_wfc_re_im` for more information
   - get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is `SR.csr` with file format being the same as that generated by [out_mat_hs2](#out_mat_hs2). Note: in the 3.10-LTS version, the command was named `get_S`
   - gen_bessel: generates projectors, i.e., a series of Bessel functions, for the DeePKS method (for LCAO basis only); see also keywords `bessel_descriptor_lmax`, `bessel_descriptor_rcut` and `bessel_descriptor_tolerence`. A file named `jle.orb` will be generated which contains the projectors. An example is provided in examples/H2O-deepks-pw
+  - gen_opt_abfs: generate opt-ABFs as discussed in this [article](https://pubs.acs.org/doi/abs/10.1021/acs.jpclett.0c00481).
   - test_memory: obtain a rough estimation of memory consuption for the calculation
   - test_neighbour: obtain information of neighboring atoms (for LCAO basis only), please specify a positive [search_radius](#search_radius) manually
 - **Default**: scf
@@ -540,7 +541,7 @@ These variables are used to control general system parameters.
 - **Default**:
   - 0:
     - if [calculation](#calculation)==md/nscf/get_pchg/get_wf/get_s or [gamma_only](#gamma_only)==True;
-    - If ([dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb or [rpa](#rpa)==True). 
+    - If ([dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0 or [rpa](#rpa)==True). 
     - If [efield_flag](#efield_flag)==1
   - 1: else
 
@@ -763,7 +764,7 @@ These variables are used to control parameters related to input files.
 - **Availability**: Used only when numerical atomic orbitals are employed as basis set.
 - **Description**: If [restart_save](#restart_save) is set to true and an electronic iteration is finished, calculations can be restarted from the charge density file, which are saved in the former calculation. Please ensure [read_file_dir](#read_file_dir) is correct, and  the charge density file exist.
 
-  If EXX(exact exchange) is calculated (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*), the Hexx(R) files in the same folder for each processor will also be read.
+  If EXX(exact exchange) is calculated (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0* or *[rpa](#rpa)==True*), the Hexx(R) files in the same folder for each processor will also be read.
 - **Default**: False
 
 ### wannier_card
@@ -1057,7 +1058,6 @@ calculations.
 
   Furthermore, the old INPUT parameter exx_hybrid_type for hybrid functionals has been absorbed into dft_functional. Options are `hf` (pure Hartree-Fock), `pbe0`(PBE0), `hse` (Note: in order to use HSE functional, LIBXC is required). Note also that HSE has been tested while PBE0 has NOT been fully tested yet, and the maximum CPU cores for running exx in parallel is $N(N+1)/2$, with N being the number of atoms.
 
-  If set to `opt_orb`, the program will not perform hybrid functional calculation. Instead, it is going to generate opt-ABFs as discussed in this [article](https://pubs.acs.org/doi/abs/10.1021/acs.jpclett.0c00481).
 - **Default**: Used the same as DFT functional as specified in the pseudopotential files.
 
 ### xc_temperature
@@ -1645,7 +1645,7 @@ These variables are used to control the output of properties.
   ---
   The circle order of the charge density on real space grids is: x is the outer loop, then y and finally z (z is moving fastest).
 
-  In EXX(exact exchange) calculations, (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*), the Hexx(R) files will be output in the folder `OUT.${suffix}` too, which can be read in NSCF calculation.
+  In EXX(exact exchange) calculations, (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0* or *[rpa](#rpa)==True*), the Hexx(R) files will be output in the folder `OUT.${suffix}` too, which can be read in NSCF calculation.
 
   In molecular dynamics simulations, the output frequency is controlled by [out_interval](#out_interval).
 - **Default**: 0 3
@@ -1965,7 +1965,7 @@ These variables are used to control the output of properties.
   - auto: These files are saved in folder `OUT.${suffix}/restart/`;
   - other: These files are saved in folder `${read_file_dir}/restart/`.
 
-  If EXX(exact exchange) is calculated (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*), the Hexx(R) files for each processor will also be saved in the above folder, which can be read in EXX calculation with *[restart_load](#restart_load)==True*.
+  If EXX(exact exchange) is calculated (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0* or *[rpa](#rpa)==True*), the Hexx(R) files for each processor will also be saved in the above folder, which can be read in EXX calculation with *[restart_load](#restart_load)==True*.
 - **Default**: False
 
 ### rpa (Under Development Feature)
@@ -2927,22 +2927,22 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b
 ### exx_opt_orb_lmax
 
 - **Type**: Integer
-- **Availability**: *[dft_functional](#dft_functional)==opt_orb*
+- **Availability**: *[calculation](#calculation)==gen_opt_abfs*
 - **Description**: The maximum l of the spherical Bessel functions, when the radial part of opt-ABFs are generated as linear combinations of spherical Bessel functions. A reasonable choice is 2.
 - **Default**: 0
 
 ### exx_opt_orb_ecut
 
 - **Type**: Real
-- **Availability**: *[dft_functional](#dft_functional)==opt_orb*
+- **Availability**: *[calculation](#calculation)==gen_opt_abfs*
 - **Description**: The cut-off of plane wave expansion, when the plane wave basis is used to optimize the radial ABFs. A reasonable choice is 60.
 - **Default**: 0
 - **Unit**: Ry
 
 ### exx_opt_orb_tolerence
 
 - **Type**: Real
-- **Availability**: *[dft_functional](#dft_functional)==opt_orb*
+- **Availability**: *[calculation](#calculation)==gen_opt_abfs*
 - **Description**: The threshold when solving for the zeros of spherical Bessel functions. A reasonable choice is 1e-12.
 - **Default**: 0
 
@@ -2975,7 +2975,7 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b
 ### exx_symmetry_realspace
 
 - **Type**: Boolean
-- **Availability**: *[symmetry](#symmetry)==1* and exx calculation  (*[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*)
+- **Availability**: *[symmetry](#symmetry)==1* and exx calculation  (*[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0* or *[rpa](#rpa)==True*)
 - **Description**: 
   - False: only rotate k-space density matrix D(k) from irreducible k-points to accelerate diagonalization
   - True: rotate both D(k) and Hexx(R) to accelerate both diagonalization and EXX calculation

source/source_esolver/esolver_ks_lcao.cpp

@@ -58,6 +58,10 @@
 #include "source_lcao/hamilt_lcao.h"
 #include "source_hsolver/hsolver_lcao.h"
 
+#ifdef __EXX
+#include "../source_lcao/module_ri/exx_opt_orb.h"
+#endif
+
 // test RDMFT
 #include "source_lcao/module_rdmft/rdmft.h"
 
@@ -132,6 +136,17 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
                                  two_center_bundle_,
                                  orb_);
 
+    if (PARAM.inp.calculation == "gen_opt_abfs")
+    {
+      #ifdef __EXX
+        Exx_Opt_Orb exx_opt_orb;
+        exx_opt_orb.generate_matrix(GlobalC::exx_info.info_opt_abfs, this->kv, ucell, this->orb_);
+      #else
+        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_all_runners", "calculation=gen_opt_abfs must compile __EXX");
+      #endif
+        return;
+    }
+
     // 4) initialize electronic wave function psi
     if (this->psi == nullptr)
     {

source/source_esolver/lcao_others.cpp

@@ -74,6 +74,10 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
         std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "testing neighbour");
         return;
     }
+    else if (cal_type == "gen_opt_abfs")
+    {
+        return;
+    }
 
     // 1. prepare HS matrices, prepare grid integral
     // (1) Find adjacent atoms for each atom.

source/source_hamilt/module_xc/exx_info.h

@@ -77,6 +77,15 @@ struct Exx_Info
     };
     Exx_Info_RI info_ri;
 
+    struct Exx_Info_Opt_ABFs
+    {
+        int abfs_Lmax = 0; // tmp
+        double ecut_exx = 60;
+        double tolerence = 1E-2;
+        double kmesh_times = 4;
+    };
+    Exx_Info_Opt_ABFs info_opt_abfs;
+
     Exx_Info() : info_lip(this->info_global), info_ri(this->info_global)
     {
     }

source/source_io/input_conv.cpp

@@ -271,7 +271,6 @@ void Input_Conv::Convert()
         ModuleBase::GlobalFunc::MAKE_DIR(GlobalC::restart.folder);
         if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0"
             || dft_functional_lower == "hse"
-            || dft_functional_lower == "opt_orb"
             || dft_functional_lower == "scan0") {
             GlobalC::restart.info_save.save_charge = true;
             GlobalC::restart.info_save.save_H = true;
@@ -297,7 +296,6 @@ void Input_Conv::Convert()
         GlobalC::restart.folder = PARAM.globalv.global_readin_dir + "restart/";
         if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0"
             || dft_functional_lower == "hse"
-            || dft_functional_lower == "opt_orb"
             || dft_functional_lower == "scan0"
             || dft_functional_lower == "lc_pbe"
             || dft_functional_lower == "lc_wpbe" 
@@ -422,13 +420,6 @@ void Input_Conv::Convert()
             }
         }
     }
-#ifdef __EXX
-    else if (dft_functional_lower == "opt_orb")
-    {
-        GlobalC::exx_info.info_global.cal_exx = false;
-        Exx_Abfs::Jle::generate_matrix = true;
-    }
-#endif
     else
     {
         GlobalC::exx_info.info_global.cal_exx = false;
@@ -456,7 +447,6 @@ void Input_Conv::Convert()
 
     if (GlobalC::exx_info.info_global.cal_exx
 #ifdef __EXX
-        || Exx_Abfs::Jle::generate_matrix
         || PARAM.inp.rpa
 #endif
         )
@@ -483,11 +473,9 @@ void Input_Conv::Convert()
         GlobalC::exx_info.info_ri.ccp_rmesh_times = std::stod(PARAM.inp.exx_ccp_rmesh_times);
         GlobalC::exx_info.info_ri.exx_symmetry_realspace = PARAM.inp.exx_symmetry_realspace;
 
-#ifdef __EXX
-        Exx_Abfs::Jle::Lmax = PARAM.inp.exx_opt_orb_lmax;
-        Exx_Abfs::Jle::Ecut_exx = PARAM.inp.exx_opt_orb_ecut;
-        Exx_Abfs::Jle::tolerence = PARAM.inp.exx_opt_orb_tolerence;
-#endif
+        GlobalC::exx_info.info_opt_abfs.abfs_Lmax  = PARAM.inp.exx_opt_orb_lmax;
+        GlobalC::exx_info.info_opt_abfs.ecut_exx = PARAM.inp.exx_opt_orb_ecut;
+        GlobalC::exx_info.info_opt_abfs.tolerence = PARAM.inp.exx_opt_orb_tolerence;
 
         // EXX does not support symmetry for nspin==4
         if (PARAM.inp.calculation != "nscf" && PARAM.inp.symmetry == "1" && PARAM.inp.nspin == 4 && PARAM.inp.basis_type == "lcao")

source/source_io/read_input_item_system.cpp

@@ -62,7 +62,7 @@ void ReadInput::item_system()
     }
     {
         Input_Item item("calculation");
-        item.annotation = "test; scf; relax; nscf; get_wf; get_pchg";
+        item.annotation = "scf; relax; md; cell-relax; nscf; get_s; get_wf; get_pchg; gen_bessel; gen_opt_abfs; test_memory; test_neighbour";
         item.read_value = [](const Input_Item& item, Parameter& para) {
             para.input.calculation = strvalue;
             std::string& calculation = para.input.calculation;
@@ -73,13 +73,14 @@ void ReadInput::item_system()
                                                 "relax",
                                                 "md",
                                                 "cell-relax",
-                                                "test_memory",
-                                                "test_neighbour",
                                                 "nscf",
                                                 "get_s",
                                                 "get_wf",
                                                 "get_pchg",
-                                                "gen_bessel"};
+                                                "gen_bessel",
+                                                "gen_opt_abfs",
+                                                "test_memory",
+                                                "test_neighbour"};
             if (std::find(callist.begin(), callist.end(), calculation) == callist.end())
             {
                 const std::string warningstr = nofound_str(callist, "calculation");

source/source_lcao/module_ri/Exx_LRI_interface.hpp

@@ -4,7 +4,6 @@
 
 #include "Exx_LRI_interface.h"
 #include "source_lcao/module_ri/exx_abfs-jle.h"
-#include "source_lcao/module_ri/exx_opt_orb.h"
 #include "source_lcao/hamilt_lcao.h"
 #include "source_lcao/module_operator_lcao/op_exx_lcao.h"
 #include "source_base/parallel_common.h"
@@ -183,15 +182,6 @@ void Exx_LRI_Interface<T, Tdata>::exx_beforescf(const int istep,
         this->cal_exx_ions(ucell,PARAM.inp.out_ri_cv);
     }
 
-    if (Exx_Abfs::Jle::generate_matrix)
-    {
-        //program should be stopped after this judgement
-        Exx_Opt_Orb exx_opt_orb;
-        exx_opt_orb.generate_matrix(kv, ucell,orb);
-        ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf");
-        return;
-    }
-
     // set initial parameter for mix_DMk_2D
     if(GlobalC::exx_info.info_global.cal_exx)
     {

source/source_lcao/module_ri/Matrix_Orbs22.cpp

@@ -55,7 +55,6 @@ void Matrix_Orbs22::init(const int mode,
     Lmax = 2 * Lmax + 1;
 
     ModuleBase::timer::tick("Matrix_Orbs22", "init");
-    std::cout << "Matrix_Orbs22::init()::done" << std::endl;
 }
 
 void Matrix_Orbs22::init_radial(const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_A1,

source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp

@@ -45,16 +45,7 @@ std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> Exx_Abfs::Construct_
 		}
 	}
 
-	for (int T = 0;  T < orbs.size() ; T++)
-	{
-		for (int L=orbs[T].size()-1; L >= 0  ; L--)
-		{
-			if (orbs[T][L].size()>0)
-				break;
-			else
-				orbs[T].resize(L);
-		}
-	}
+	Exx_Abfs::Construct_Orbs::filter_empty_orbs(orbs);
 	return orbs;
 }
 
@@ -548,4 +539,18 @@ std::vector<double> Exx_Abfs::Construct_Orbs::get_Rcut(const std::vector<std::ve
     }
 
     return Rcut;
-}
+}
+
+void Exx_Abfs::Construct_Orbs::filter_empty_orbs(std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orbs)
+{
+	for (int T=0;  T<orbs.size(); ++T)
+	{
+		for (int L=orbs[T].size()-1; L>=0 ; --L)
+		{
+			if (orbs[T][L].size()>0)
+				{ break; }
+			else
+				{ orbs[T].resize(L); }
+		}
+	}
+}

source/source_lcao/module_ri/exx_abfs-construct_orbs.h

@@ -57,6 +57,9 @@ class Exx_Abfs::Construct_Orbs
         return get_Rmax(rcut);
     }
 
+	static void filter_empty_orbs(
+		std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orbs);
+
   private:
 	static std::vector<std::vector<std::vector<std::vector<double>>>> psi_mult_psi( 
 		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &lcaos );

source/source_lcao/module_ri/exx_abfs-jle.cpp

@@ -1,45 +1,42 @@
 #include "exx_abfs-jle.h"
 
 #include "source_io/module_parameter/parameter.h"
-#include "../../source_pw/module_pwdft/global.h"
 #include "../../source_basis/module_ao/ORB_read.h"
-#include "../../source_base/global_function.h"
+#include "../../source_cell/unitcell.h"
 #include "../../source_base/mathzone.h"
 #include "../../source_base/math_sphbes.h" // mohan add 2021-05-06
 
-bool Exx_Abfs::Jle::generate_matrix = false;
-int Exx_Abfs::Jle::Lmax = 2;
-double Exx_Abfs::Jle::Ecut_exx = 60;
-double Exx_Abfs::Jle::tolerence = 1.0e-12;	
-
-void Exx_Abfs::Jle::init_jle(const double kmesh_times, 
-							 const UnitCell& ucell,
-							 const LCAO_Orbitals& orb)
+std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>
+Exx_Abfs::Jle::init_jle(
+	const Exx_Info::Exx_Info_Opt_ABFs &info,
+	const double kmesh_times, 
+	const UnitCell& ucell,
+	const LCAO_Orbitals& orb)
 {
-	jle.resize( ucell.ntype );
+	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> jle( ucell.ntype );
 
-	for (int T = 0;  T < ucell.ntype ; T++)
+	for(int T=0; T<ucell.ntype; ++T)
 	{
-		jle[T].resize( Lmax+1 );
-		for (int L=0; L <= Lmax ; ++L)
+		jle[T].resize( info.abfs_Lmax+1 );
+		for(int L=0; L<=info.abfs_Lmax; ++L)
 		{
 			const size_t ecut_number 
-				= static_cast<size_t>( sqrt( Ecut_exx ) * orb.Phi[T].getRcut() / ModuleBase::PI ); // Rydberg Unit.
+				= static_cast<size_t>( std::sqrt( info.ecut_exx ) * orb.Phi[T].getRcut() / ModuleBase::PI ); // Rydberg Unit.
 
 			jle[T][L].resize( ecut_number );
 
 			std::vector<double> en(ecut_number, 0.0);
-			ModuleBase::Sphbes::Spherical_Bessel_Roots(ecut_number, L, tolerence, ModuleBase::GlobalFunc::VECTOR_TO_PTR(en), orb.Phi[T].getRcut());
+			ModuleBase::Sphbes::Spherical_Bessel_Roots(ecut_number, L, info.tolerence, en.data(), orb.Phi[T].getRcut());
 
-			for(size_t E=0; E!=ecut_number; ++E)
+			for(size_t E=0; E<ecut_number; ++E)
 			{
 				std::vector<double> jle_r( orb.Phi[T].PhiLN(0,0).getNr() );
 				ModuleBase::Sphbes::Spherical_Bessel(
 					orb.Phi[T].PhiLN(0,0).getNr(), 
 					orb.Phi[T].PhiLN(0,0).getRadial(), 
 					en[E], 
 					L, 
-					ModuleBase::GlobalFunc::VECTOR_TO_PTR(jle_r));
+					jle_r.data());
 				jle[T][L][E].set_orbital_info(
 					orb.Phi[T].PhiLN(0,0).getLabel(),
 					orb.Phi[T].PhiLN(0,0).getType(),
@@ -49,13 +46,16 @@ void Exx_Abfs::Jle::init_jle(const double kmesh_times,
 					orb.Phi[T].PhiLN(0,0).getRab(),
 					orb.Phi[T].PhiLN(0,0).getRadial(),
 					Numerical_Orbital_Lm::Psi_Type::Psi,
-					ModuleBase::GlobalFunc::VECTOR_TO_PTR(jle_r),
+					jle_r.data(),
 					static_cast<int>(orb.Phi[T].PhiLN(0,0).getNk() * kmesh_times) | 1,
 					orb.Phi[T].PhiLN(0,0).getDk(),
 					orb.Phi[T].PhiLN(0,0).getDruniform(),
 					false,
-					true, PARAM.inp.cal_force);
+					true,
+					PARAM.inp.cal_force);
 			}
 		}
 	}
+
+	return jle;
 }