add RDMFT(Reduced Density Matrix Functional Theory) code （Useful Information: for those who what to contribute codes to ABACUS, this is a good example)

CMakeLists.txt

@@ -724,7 +724,8 @@ if(ENABLE_LCAO)
     hcontainer
     deltaspin
     numerical_atomic_orbitals
-    lr)
+    lr
+    rdmft)
   if(USE_ELPA)
     target_link_libraries(${ABACUS_BIN_NAME} genelpa)
   endif()

docs/advanced/input_files/input-main.md

@@ -435,6 +435,9 @@
     - [abs\_broadening](#abs_broadening)
     - [ri\_hartree\_benchmark](#ri_hartree_benchmark)
     - [aims\_nbasis](#aims_nbasis)
+  - [Reduced Density Matrix Functional Theory](#Reduced-Density-Matrix-Functional-Theory)
+    - [rdmft](#rdmft)
+    - [rdmft\_power\_alpha](#rdmft_power_alpha)
 
 [back to top](#full-list-of-input-keywords)
 ## System variables
@@ -4031,4 +4034,21 @@ The output files are `OUT.${suffix}/Excitation_Energy.dat` and `OUT.${suffix}/Ex
 - **Description**: Atomic basis set size for each atom type (with the same order as in `STRU`) in FHI-aims.
 - **Default**: {} (empty list, where ABACUS use its own basis set size)
 
+## Reduced Density Matrix Functional Theory
+
+ab-initio methods and the xc-functional parameters used in RDMFT.
+The physical quantities that RDMFT temporarily expects to output are the kinetic energy, total energy, and 1-RDM of the system in the ground state, etc.
+
+### rdmft
+
+- **Type**: Boolean
+- **Description**: Whether to perform rdmft calculation (reduced density matrix funcional theory)
+- **Default**: false
+
+### rdmft_power_alpha
+
+- **Type**: Real
+- **Description**: The alpha parameter of power-functional(or other exx-type/hybrid functionals) which used in RDMFT, g(occ_number) = occ_number^alpha
+- **Default**: 0.656
+
 [back to top](#full-list-of-input-keywords)

source/CMakeLists.txt

@@ -18,6 +18,9 @@ add_subdirectory(module_ri)
 add_subdirectory(module_parameter)
 add_subdirectory(module_lr)
 
+# add by jghan
+add_subdirectory(module_rdmft)
+
 add_library(
     driver
     OBJECT

source/Makefile.Objects

@@ -74,6 +74,7 @@ VPATH=./src_global:\
 ./module_lr/operator_casida:\
 ./module_lr/potentials:\
 ./module_lr/utils:\
+./module_rdmft:\
 ./\
 
 OBJS_ABACUS_PW=${OBJS_MAIN}\
@@ -115,6 +116,7 @@ ${OBJS_DELTASPIN}\
 ${OBJS_TENSOR}\
 ${OBJS_HSOLVER_PEXSI}\
 ${OBJS_LR}\
+${OBJS_RDMFT}
 
 OBJS_MAIN=main.o\
     driver.o\
@@ -727,3 +729,8 @@ OBJS_TENSOR=tensor.o\
     lr_spectrum.o\
     hamilt_casida.o\
     esolver_lrtd_lcao.o\
+
+OBJS_RDMFT=rdmft.o\
+        rdmft_tools.o\
+        rdmft_cal_hamilt.o\
+        update_state_rdmft.o\

source/module_esolver/esolver_ks_lcao.cpp

@@ -62,6 +62,10 @@
 // #include "module_elecstate/cal_dm.h"
 //---------------------------------------------------
 
+// test RDMFT
+#include "module_rdmft/rdmft.h"
+#include <iostream>
+
 namespace ModuleESolver
 {
 
@@ -249,6 +253,13 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
         }
     }
 
+    // 14) initialize rdmft, added by jghan
+    if( PARAM.inp.rdmft == true )
+    {
+        rdmft_solver.init( this->GG, this->GK, this->pv, ucell, this->kv, *(this->pelec),
+                                this->orb_, two_center_bundle_, PARAM.inp.dft_functional, PARAM.inp.rdmft_power_alpha);
+    }
+
     ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners");
     return;
 }
@@ -838,6 +849,7 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(UnitCell& ucell, const int istep, const
     {
         this->pelec->cal_converged();
     }
+
 }
 
 //------------------------------------------------------------------------------
@@ -1041,6 +1053,23 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
     ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels");
 #endif
 
+    /******** test RDMFT *********/
+    if ( PARAM.inp.rdmft == true ) // rdmft, added by jghan, 2024-10-17
+    {
+        ModuleBase::matrix occ_number_ks(this->pelec->wg);
+        for(int ik=0; ik < occ_number_ks.nr; ++ik) { for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik]; }
+        this->rdmft_solver.update_elec(occ_number_ks, *(this->psi));
+
+        //initialize the gradients of Etotal on occupation numbers and wfc, and set all elements to 0. 
+        ModuleBase::matrix dE_dOccNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
+        psi::Psi<TK> dE_dWfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis()); 
+        dE_dWfc.zero_out();
+
+        double Etotal_RDMFT = this->rdmft_solver.run(dE_dOccNum, dE_dWfc);
+    }
+    /******** test RDMFT *********/
+
+
 #ifdef __EXX
     // 11) write rpa information
     if (PARAM.inp.rpa)
@@ -1227,6 +1256,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
     }
 }
 
+
 //------------------------------------------------------------------------------
 //! the 20th,21th,22th functions of ESolver_KS_LCAO
 //! mohan add 2024-05-11

source/module_esolver/esolver_ks_lcao.h

@@ -12,6 +12,9 @@
 #include "module_basis/module_nao/two_center_bundle.h"
 #include "module_io/output_mat_sparse.h"
 
+// added by jghan for rdmft calculation
+#include "module_rdmft/rdmft.h"
+
 #include <memory>
 
 namespace LR
@@ -68,6 +71,8 @@ class ESolver_KS_LCAO : public ESolver_KS<TK> {
 
     TwoCenterBundle two_center_bundle_;
 
+    rdmft::RDMFT<TK, TR> rdmft_solver;  // added by jghan for rdmft calculation, 2024-03-16
+
     // temporary introduced during removing GlobalC::ORB
     LCAO_Orbitals orb_;
 

source/module_esolver/lcao_before_scf.cpp

@@ -369,6 +369,14 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
 
     this->p_hamilt->non_first_scf = istep;
+
+    // update in ion-step
+    if( PARAM.inp.rdmft == true )
+    {
+        // necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp
+        rdmft_solver.update_ion(ucell, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac);   // add by jghan, 2024-03-16/2024-10-08
+    }
+
     return;
 }
 

source/module_hamilt_general/module_xc/xc_functional.cpp

@@ -18,6 +18,7 @@ std::vector<int> XC_Functional::func_id(1);
 int XC_Functional::func_type = 0;
 bool XC_Functional::use_libxc = true;
 double XC_Functional::hybrid_alpha = 0.25;
+std::map<int, double> XC_Functional::scaling_factor_xc = { {1, 1.0} }; // added by jghan, 2024-10-10
 
 void XC_Functional::set_hybrid_alpha(const double alpha_in)
 {
@@ -61,6 +62,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
     //        func_id.push_back(XC_GGA_C_PBE);
 
     func_id.clear();
+    scaling_factor_xc.clear(); // added by jghan, 2024-07-07
     std::string xc_func = xc_func_in;
     std::transform(xc_func.begin(), xc_func.end(), xc_func.begin(), (::toupper));
 	if( xc_func == "LDA" || xc_func == "PZ" || xc_func == "SLAPZNOGXNOGC") //SLA+PZ
@@ -196,13 +198,60 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
     {
         // not doing anything
     }
+    else if( xc_func == "MULLER" || xc_func == "POWER" ) // added by jghan, 2024-07-06
+    {
+        func_type = 4;
+        use_libxc = false;
+    }
 #ifdef USE_LIBXC
     else if( xc_func == "HSE")
     {
         func_id.push_back(XC_HYB_GGA_XC_HSE06);
         func_type = 4;
         use_libxc = true;
     }
+    // added by jghan, 2024-07-06
+    else if( xc_func == "WP22")
+    {
+        func_id.push_back(XC_GGA_X_ITYH);   // short-range of B88_X, id=529
+        func_id.push_back(XC_GGA_C_LYPR);   // short-range of LYP_C, id=624
+        func_type = 4;
+        use_libxc = true;
+    }
+    else if( xc_func == "CWP22")
+    {   
+        // BLYP_XC_lr = -BLYP_XC_sr + BLYP_XC, the realization of it is in v_xc_libxc() function, xc_functional_libxc_vxc.cpp
+        func_id.push_back(XC_GGA_X_ITYH);   // short-range of B88_X, id=529
+        func_id.push_back(XC_GGA_C_LYPR);   // short-range of LYP_C, id=624
+        func_id.push_back(XC_GGA_X_B88);    // complete B88_X, id=106
+        func_id.push_back(XC_GGA_C_LYP);    // complete LYP_C, id=131
+
+        // the scaling factor of CWP22-functionals
+        scaling_factor_xc[XC_GGA_X_ITYH] = -1.0;
+        scaling_factor_xc[XC_GGA_C_LYPR] = -1.0;
+        scaling_factor_xc[XC_GGA_X_B88] = 1.0;
+        scaling_factor_xc[XC_GGA_X_B88] = 1.0;
+
+        func_type = 4;
+        use_libxc = true;
+    }
+    else if( xc_func == "BLYP_LR")
+    {   
+        // BLYP_XC_lr = -BLYP_XC_sr + BLYP_XC, the realization of it is in v_xc_libxc() function, xc_functional_libxc_vxc.cpp
+        func_id.push_back(XC_GGA_X_ITYH);   // short-range of B88_X, id=529
+        func_id.push_back(XC_GGA_C_LYPR);   // short-range of LYP_C, id=624
+        func_id.push_back(XC_GGA_X_B88);    // complete B88_X, id=106
+        func_id.push_back(XC_GGA_C_LYP);    // complete LYP_C, id=131
+
+        // the scaling factor of BLYP_LR-functionals
+        scaling_factor_xc[XC_GGA_X_ITYH] = -1.0;
+        scaling_factor_xc[XC_GGA_C_LYPR] = -1.0;
+        scaling_factor_xc[XC_GGA_X_B88] = 1.0;
+        scaling_factor_xc[XC_GGA_X_B88] = 1.0;
+
+        func_type = 2;
+        use_libxc = true;
+    }
 #endif
     else
     {
@@ -243,7 +292,8 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
 #endif
 
 #ifndef USE_LIBXC
-    if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0")
+    if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0" 
+        || xc_func == "MULLER" || xc_func == "POWER" || xc_func == "WP22" || xc_func == "CWP22")
     {
         ModuleBase::WARNING_QUIT("set_xc_type","to use SCAN, SCAN0, or HSE, LIBXC is required");
     }

source/module_hamilt_general/module_xc/xc_functional.h

@@ -20,6 +20,8 @@
 #include "module_elecstate/module_charge/charge.h"
 #include "module_cell/unitcell.h"
 
+#include <map> // added by jghan, 2024-10-10
+
 class XC_Functional
 {
 	public:
@@ -80,6 +82,10 @@ class XC_Functional
 	//exx_hybrid_alpha for mixing exx in hybrid functional:
 	static double hybrid_alpha;
 
+	// added by jghan, 2024-07-07
+	// as a scaling factor for different xc-functionals
+	static std::map<int, double> scaling_factor_xc;
+
 	public:
 	static std::vector<int> get_func_id() { return func_id; }
 
@@ -162,7 +168,7 @@ class XC_Functional
                          ModulePW::PW_Basis* rhopw,
                          const UnitCell* ucell,
                          std::vector<double>& stress_gga,
-                         const bool is_stress = 0);
+                         const bool is_stress = false);
 	template <typename T, typename Device,
           typename Real = typename GetTypeReal<T>::type>
 	static void grad_wfc(

source/module_hamilt_general/module_xc/xc_functional_libxc.cpp

@@ -135,6 +135,20 @@ std::vector<xc_func_type> XC_Functional_Libxc::init_func(const std::vector<int>
 				GlobalC::exx_info.info_global.hse_omega };
 			xc_func_set_ext_params(&funcs.back(), parameter_hse);
 		}
+        // added by jghan, 2024-07-06
+		else if( id == XC_GGA_X_ITYH ) // short-range of B88_X
+		{
+			add_func( XC_GGA_X_ITYH );
+			double parameter_omega[1] = {PARAM.inp.exx_hse_omega}; // GlobalC::exx_info.info_global.hse_omega
+			xc_func_set_ext_params(&funcs.back(), parameter_omega);	
+		}
+		else if( id == XC_GGA_C_LYPR ) // short-range of LYP_C
+		{
+			add_func( XC_GGA_C_LYPR );
+            // the first six parameters come from libxc, and may need to be modified in some cases
+			double parameter_lypr[7] = {0.04918, 0.132, 0.2533, 0.349, 0.35/2.29, 2.0/2.29, PARAM.inp.exx_hse_omega};
+			xc_func_set_ext_params(&funcs.back(), parameter_lypr);	
+		}
 #endif
 		else
 		{

source/module_hamilt_general/module_xc/xc_functional_libxc.h

@@ -11,6 +11,9 @@
 #include <tuple>
 #include <vector>
 
+#include <map> // added by jghan, 2024-10-10
+#include <utility>
+
 class Charge;
 
 namespace XC_Functional_Libxc
@@ -33,12 +36,13 @@ namespace XC_Functional_Libxc
 //  xc_functional_libxc_vxc.cpp
 //-------------------
 
-    extern std::tuple<double,double,ModuleBase::matrix> v_xc_libxc(
-        const std::vector<int> &func_id,
-        const int &nrxx, // number of real-space grid
-        const double &omega, // volume of cell
-        const double tpiba,
-        const Charge* const chr); // charge density
+	extern std::tuple<double,double,ModuleBase::matrix> v_xc_libxc(
+		const std::vector<int> &func_id,
+		const int &nrxx, // number of real-space grid
+		const double &omega, // volume of cell
+		const double tpiba,
+		const Charge* const chr, // charge density
+		const std::map<int, double>* scaling_factor = nullptr); // added by jghan, 2024-10-10
 
     // for mGGA functional
     extern std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> v_xc_meta(

source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp

@@ -18,7 +18,8 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
         const int &nrxx, // number of real-space grid
         const double &omega, // volume of cell
         const double tpiba,
-        const Charge* const chr)
+        const Charge* const chr,
+        const std::map<int, double>* scaling_factor)
 {
     ModuleBase::TITLE("XC_Functional_Libxc","v_xc_libxc");
     ModuleBase::timer::tick("XC_Functional_Libxc","v_xc_libxc");
@@ -109,14 +110,25 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
                 break;
         }
 
-        etxc += XC_Functional_Libxc::convert_etxc(nspin, nrxx, sgn, rho, exc);
+        // added by jghan, 2024-10-10
+        double factor = 1.0;
+        if( scaling_factor == nullptr ) { ;
+        } else
+        {
+            auto pair_factor = scaling_factor->find(func.info->number);
+            if( pair_factor != scaling_factor->end() ) { factor = pair_factor->second;
+}
+        }
+
+        // time factor is added by jghan, 2024-10-10
+        etxc += XC_Functional_Libxc::convert_etxc(nspin, nrxx, sgn, rho, exc) * factor;
         const std::pair<double,ModuleBase::matrix> vtxc_v = XC_Functional_Libxc::convert_vtxc_v(
             func, nspin, nrxx,
             sgn, rho, gdr,
             vrho, vsigma,
             tpiba, chr);
-        vtxc += std::get<0>(vtxc_v);
-        v += std::get<1>(vtxc_v);
+        vtxc += std::get<0>(vtxc_v) * factor;
+        v += std::get<1>(vtxc_v) * factor;
     } // end for( xc_func_type &func : funcs )
 
     if(4==PARAM.inp.nspin)
@@ -240,10 +252,12 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
 #endif
         for( int ir=0; ir<nrxx; ++ir )
         {
-            if ( rho[ir*2]<rho_th || sqrt(std::abs(sigma[ir*3]))<grho_th || std::abs(kin_r[ir*2])<tau_th)
+            if ( rho[ir*2]<rho_th || sqrt(std::abs(sigma[ir*3]))<grho_th || std::abs(kin_r[ir*2])<tau_th) {
                 sgn[ir*2] = 0.0;
-            if ( rho[ir*2+1]<rho_th || sqrt(std::abs(sigma[ir*3+2]))<grho_th || std::abs(kin_r[ir*2+1])<tau_th)
+}
+            if ( rho[ir*2+1]<rho_th || sqrt(std::abs(sigma[ir*3+2]))<grho_th || std::abs(kin_r[ir*2+1])<tau_th) {
                 sgn[ir*2+1] = 0.0;
+}
         }
     }