deepmodeling · JGHan7 · Dec 20, 2023 · Dec 20, 2023 · Dec 20, 2023 · Dec 20, 2023
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -730,7 +730,8 @@ if(ENABLE_LCAO)
     hcontainer
     deltaspin
     numerical_atomic_orbitals
-    lr)
+    lr
+    rdmft)
   if(USE_ELPA)
     target_link_libraries(${ABACUS_BIN_NAME} genelpa)
   endif()

diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -427,6 +427,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)
+    - [ab\_initio\_type](#ab_initio_type)
+    - [rdmft\_power\_alpha](#rdmft_power_alpha)
 
 [back to top](#full-list-of-input-keywords)
 ## System variables
@@ -3999,4 +4002,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.
+
+### ab_initio_type
+
+- **Type**: String
+- **Description**: Determine the ab-initio method, ksdft or rdmft (reduced density matrix funcional theory)
+- **Default**: ksdft
+
+### 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)
diff --git a/source/CMakeLists.txt b/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

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -72,6 +72,7 @@ VPATH=./src_global:\
 ./module_lr/operator_casida:\
 ./module_lr/potentials:\
 ./module_lr/utils:\
+./module_rdmft:\
 ./\
 
 OBJS_ABACUS_PW=${OBJS_MAIN}\
@@ -113,6 +114,7 @@ ${OBJS_DELTASPIN}\
 ${OBJS_TENSOR}\
 ${OBJS_HSOLVER_PEXSI}\
 ${OBJS_LR}\
+${OBJS_RDMFT}
 
 OBJS_MAIN=main.o\
     driver.o\
@@ -730,3 +732,8 @@ OBJS_TENSOR=tensor.o\
     lr_spectrum.o\
     hamilt_casida.o\
     esolver_lrtd_lcao.o\
+
+OBJS_RDMFT=rdmft.o\
+        rdmft_tools.o\
+        cal_V_rdmft.o\
+        update_state_rdmft.o\
diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h
@@ -9,7 +9,6 @@
 #include "module_hsolver/hsolver.h"
 #include "module_io/cal_test.h"
 #include "module_psi/psi.h"
-
 #include <fstream>
 #include <cstring>
 namespace ModuleESolver
@@ -94,6 +93,9 @@ class ESolver_KS : public ESolver_FP
 
 		wavefunc wf;
 
+        // // TR need be std::complex<double> when NSPIN = 4
+        // rdmft::RDMFT<T, double> rdmft_solver;  // add by jghan for rdmft calculation
+
         // wavefunction coefficients
         psi::Psi<T>* psi = nullptr;
 

diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -55,6 +55,12 @@
 #include "module_io/write_dmr.h"
 #include "module_io/write_wfc_nao.h"
 
+// test RDMFT
+#include "module_rdmft/rdmft.h"
+#include "module_elecstate/module_dm/cal_dm_psi.h"  //temp ,delete
+
+#include <iostream>
+
 namespace ModuleESolver
 {
 
@@ -256,7 +262,15 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(const Input_para& inp, UnitCell
         this->pelec->fixed_weights(PARAM.inp.ocp_kb, PARAM.inp.nbands, PARAM.inp.nelec);
     }
 
-    // 15) if kpar is not divisible by nks, print a warning
+
+    // 15) initialize rdmft, added by jghan
+    if( PARAM.inp.ab_initio_type == "rdmft" )
+    {
+        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);
+    }
+
+    // 16) if kpar is not divisible by nks, print a warning
     if (GlobalV::KPAR_LCAO > 1)
     {
         if (this->kv.get_nks() % GlobalV::KPAR_LCAO != 0)
@@ -902,6 +916,7 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
     {
         this->pelec->cal_converged();
     }
+
 }
 
 //------------------------------------------------------------------------------
@@ -939,7 +954,24 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
         }
     }
 
+    // 2.5) determine whether rdmft needs to get the initial value, added by jghan, 2024-10-25
+    int one_step_exx = false;
+    bool get_init_value_rdmft = false;
+    if( iter == 1 ) get_init_value_rdmft = true; // the case without hybrid functionals
 #ifdef __EXX
+    if( GlobalC::exx_info.info_global.cal_exx )
+    {
+        if( this->conv_esolver ) one_step_exx = true;
+        // the case with hybrid functionals, calculate rdmft just after updateExx, cal once in  one inner loop
+        if( one_step_exx ) get_init_value_rdmft = true;
+        else get_init_value_rdmft = false;
+    }
+#endif
+
+#ifdef __EXX
+
+    if( GlobalC::exx_info.info_global.cal_exx && this->conv_esolver ) one_step_exx = true;
+
     // 3) save exx matrix
     int two_level_step = GlobalC::exx_info.info_ri.real_number ? this->exd->two_level_step : this->exc->two_level_step;
 
@@ -1068,6 +1100,37 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
     {
         GlobalC::dftu.initialed_locale = true;
     }
+
+    // 7) rdmft, added by jghan, 2024-10-25
+    if ( PARAM.inp.ab_initio_type == "rdmft" )
+    {
+        ModuleBase::TITLE("RDMFT", "E & Egradient");
+        ModuleBase::timer::tick("RDMFT", "E & Egradient");
+
+        if( get_init_value_rdmft )
+        {
+            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));
+            std::cout << "\n******\n" << "update elec in rdmft successfully" << "\n******\n" << std::endl;
+
+            //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);
+
+            ModuleBase::timer::tick("RDMFT", "E & Egradient");
+
+            // break;
+        }
+    }
+
 }
 
 //------------------------------------------------------------------------------
@@ -1162,6 +1225,29 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
     ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels");
 #endif
 
+    /******** test RDMFT *********/
+
+    // rdmft, added by jghan, 2024-10-17
+    if ( PARAM.inp.ab_initio_type == "rdmft" )
+    {
+        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)
@@ -1287,6 +1373,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
     }
 }
 
+
 //------------------------------------------------------------------------------
 //! the 20th,21th,22th functions of ESolver_KS_LCAO
 //! mohan add 2024-05-11

diff --git a/source/module_esolver/esolver_ks_lcao.h b/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
@@ -78,6 +81,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_;
 

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

diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/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,68 @@ 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;
+    }
+    // just for test BLYP in libxc, added by jghan, 2024-10-22
+    else if( xc_func == "BLYP_LIBXC")
+    {
+        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
+        func_type = 2;
+        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 +300,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");
     }

diff --git a/source/module_hamilt_general/module_xc/xc_functional.h b/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(