Merge branch 'develop' into ucell3

Author: A-006

CMakeLists.txt

@@ -725,7 +725,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/elec_properties/hs_matrix.md

@@ -61,6 +61,11 @@ We also offer the option of only calculating the overlap matrix without running
 
 A file named `SR.csr` will be generated in the working directory, which contains the overlap matrix.
 
+> When `nspin` is set to 1 or 2, the dimension of the overlap matrix is nlocal $\times$ nlocal, where nlocal is the total number of numerical atomic orbitals. 
+These numerical atomic orbitals are ordered from outer to inner loop as atom, angular quantum number $l$, zeta (multiple radial orbitals corresponding to each $l$), and magnetic quantum number $m$. 
+When `nspin` is set to 4, the dimension of the overlap matrix is (2 $\times$ nlocal) $\times$ (2 $\times$ nlocal). In this case, the numerical atomic orbitals are ordered from outer to inner loop as atom, angular quantum number $l$, zeta (multiple radial orbitals corresponding to each $l$), magnetic quantum number $m$, and npol (index of spin, ranges from 0 to 1).
+
+
 ## examples
 We provide [examples](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/matrix_hs) of outputting the matrices. There are four examples:
 

docs/advanced/elec_properties/position_matrix.md

@@ -18,4 +18,7 @@ Each file or each section of the appended file starts with "STEP: " followed by
 
 Each block here contains the matrix for the corresponding cell. There are three columns in each block, giving the matrix elements in x, y, z directions, respectively. There are altogether nbasis * nbasis lines in each block, which emulates the matrix elements.
 
-In molecular dynamics (MD) calculations, if [out_app_flag](../input_files/input-main.md#out_app_flag) is set to true, then `data-rR-tr` is written in an append manner. Otherwise, output files will be put in a separate directory, `matrix`, and named as `$x`_data-rR-tr, where `$x` is the number of MD step. In addition, The output frequency is controlled by [out_interval](../input_files/input-main.md#out_interval). For example, if we are running a 10-step MD with out_interval = 3, then `$x` will be 0, 3, 6, and 9.
+In molecular dynamics (MD) calculations, if [out_app_flag](../input_files/input-main.md#out_app_flag) is set to true, then `data-rR-tr` is written in an append manner. Otherwise, output files will be put in a separate directory, `matrix`, and named as `$x`_data-rR-tr, where `$x` is the number of MD step. In addition, the output frequency is controlled by [out_interval](../input_files/input-main.md#out_interval). For example, if we are running a 10-step MD with out_interval = 3, then `$x` will be 0, 3, 6, and 9.
+
+## get_S
+We also offer the option of only calculating the position matrix without running SCF. For that purpose, in `INPUT` file we need to set the keyword [calculation](../input_files/input-main.md#calculation) to `get_S`, and [out_mat_r](../input_files/input-main.md#out_mat_r) to `true`.

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
@@ -1717,7 +1720,7 @@ These variables are used to control the output of properties.
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital basis (not gamma-only algorithm)
-- **Description**: Whether to print the matrix representation of the position matrix (in Bohr) into a file named `data-rR-tr` in the directory `OUT.${suffix}`. For more information, please refer to [position_matrix.md](../elec_properties/position_matrix.md#extracting-position-matrices).
+- **Description**: Whether to print the matrix representation of the position matrix (in Bohr) into a file named `data-rR-tr` in the directory `OUT.${suffix}`. If [calculation](#calculation) is set to `get_S`, the position matrix can be obtained without scf iterations. For more information, please refer to [position_matrix.md](../elec_properties/position_matrix.md#extracting-position-matrices).
 - **Default**: False
 
 ### out_mat_hs2
@@ -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\
@@ -436,7 +438,6 @@ OBJS_SURCHEM=surchem.o\
     cal_totn.o\
     cal_vcav.o\
     cal_vel.o\
-    corrected_energy.o\
     minimize_cg.o\
     sol_force.o\
 
@@ -726,3 +727,8 @@ OBJS_TENSOR=tensor.o\
     lr_spectrum.o\
     hamilt_casida.o\
     esolver_lrtd_lcao.o\
+
+OBJS_RDMFT=rdmft.o\
+        rdmft_tools.o\
+        rdmft_pot.o\
+        update_state_rdmft.o\

source/driver.cpp

@@ -19,9 +19,6 @@ Driver::Driver()
 
 Driver::~Driver()
 {
-    // Release the device memory within singleton object GlobalC::ppcell
-    // before the main function exits.
-    GlobalC::ppcell.release_memory();
 }
 
 void Driver::init()

source/module_cell/cell_index.h

@@ -12,7 +12,7 @@
  *    the maximum L of a given atom, the number of chi of a given atom and a given L, the atom label of a given atom,
  * etc. The class provides the interface to get the information of the atoms and orbitals indices in the unit cell. It
  * also provides the interface to write the orbital information into a file. It is now used only in the Output_Mulliken
- * class. However, it is supposed to be able to replace GlobalC::ucell in LCAO codes where only cell indices are needed.
+ * class. However, it is supposed to be able to replace ucell in LCAO codes where only cell indices are needed.
  * Take care that the three key vectors atomCounts, lnchiCounts, and atomLabels should be set from the UnitCell class.
  *    It depends on nspin because the functions get_nw() and get_iwt() are related to nspin, and can
  *    be used in the LCAO Hamiltonian construction.

source/module_cell/klist.cpp

@@ -59,7 +59,8 @@ int K_Vectors::get_ik_global(const int& ik, const int& nkstot)
     }
 }
 
-void K_Vectors::set(const ModuleSymmetry::Symmetry& symm,
+void K_Vectors::set(const UnitCell& ucell,
+                    const ModuleSymmetry::Symmetry& symm,
                     const std::string& k_file_name,
                     const int& nspin_in,
                     const ModuleBase::Matrix3& reciprocal_vec,
@@ -93,7 +94,7 @@ void K_Vectors::set(const ModuleSymmetry::Symmetry& symm,
     this->nspin = (this->nspin == 4) ? 1 : this->nspin;
 
     // read KPT file and generate K-point grid
-    bool read_succesfully = this->read_kpoints(k_file_name);
+    bool read_succesfully = this->read_kpoints(ucell,k_file_name);
 #ifdef __MPI
     Parallel_Common::bcast_bool(read_succesfully);
 #endif
@@ -113,7 +114,7 @@ void K_Vectors::set(const ModuleSymmetry::Symmetry& symm,
     {
         bool match = true;
         // calculate kpoints in IBZ and reduce kpoints according to symmetry
-        this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, GlobalC::ucell, match);
+        this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, ucell, match);
 #ifdef __MPI
         Parallel_Common::bcast_bool(match);
 #endif
@@ -128,7 +129,7 @@ void K_Vectors::set(const ModuleSymmetry::Symmetry& symm,
                 std::cout << "Automatically set symmetry to 0 and continue ..." << std::endl;
                 ModuleSymmetry::Symmetry::symm_flag = 0;
                 match = true;
-                this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, GlobalC::ucell, match);
+                this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, ucell, match);
             } else {
                 ModuleBase::WARNING_QUIT("K_Vectors::ibz_kpoint",
                                          "Possible solutions: \n \
@@ -209,7 +210,8 @@ void K_Vectors::renew(const int& kpoint_number)
 
 // Read the KPT file, which contains K-point coordinates, weights, and grid size information
 // Generate K-point grid according to different parameters of the KPT file
-bool K_Vectors::read_kpoints(const std::string& fn)
+bool K_Vectors::read_kpoints(const UnitCell& ucell,
+                             const std::string& fn)
 {
     ModuleBase::TITLE("K_Vectors", "read_kpoints");
     if (GlobalV::MY_RANK != 0)
@@ -236,16 +238,16 @@ bool K_Vectors::read_kpoints(const std::string& fn)
             ModuleBase::WARNING_QUIT("K_Vectors", "kspacing should > 0");
         };
         // number of K points = max(1,int(|bi|/KSPACING+1))
-        ModuleBase::Matrix3 btmp = GlobalC::ucell.G;
+        ModuleBase::Matrix3 btmp = ucell.G;
         double b1 = sqrt(btmp.e11 * btmp.e11 + btmp.e12 * btmp.e12 + btmp.e13 * btmp.e13);
         double b2 = sqrt(btmp.e21 * btmp.e21 + btmp.e22 * btmp.e22 + btmp.e23 * btmp.e23);
         double b3 = sqrt(btmp.e31 * btmp.e31 + btmp.e32 * btmp.e32 + btmp.e33 * btmp.e33);
         int nk1
-            = std::max(1, static_cast<int>(b1 * ModuleBase::TWO_PI / PARAM.inp.kspacing[0] / GlobalC::ucell.lat0 + 1));
+            = std::max(1, static_cast<int>(b1 * ModuleBase::TWO_PI / PARAM.inp.kspacing[0] / ucell.lat0 + 1));
         int nk2
-            = std::max(1, static_cast<int>(b2 * ModuleBase::TWO_PI / PARAM.inp.kspacing[1] / GlobalC::ucell.lat0 + 1));
+            = std::max(1, static_cast<int>(b2 * ModuleBase::TWO_PI / PARAM.inp.kspacing[1] / ucell.lat0 + 1));
         int nk3
-            = std::max(1, static_cast<int>(b3 * ModuleBase::TWO_PI / PARAM.inp.kspacing[2] / GlobalC::ucell.lat0 + 1));
+            = std::max(1, static_cast<int>(b3 * ModuleBase::TWO_PI / PARAM.inp.kspacing[2] / ucell.lat0 + 1));
 
         GlobalV::ofs_warning << " Generate k-points file according to KSPACING: " << fn << std::endl;
         std::ofstream ofs(fn.c_str());

source/module_cell/klist.h

@@ -50,7 +50,8 @@ class K_Vectors
      *       it will output a warning and suggest possible solutions.
      * @note Only available for nspin = 1 or 2 or 4.
      */
-    void set(const ModuleSymmetry::Symmetry& symm,
+    void set(const UnitCell& ucell,
+        const ModuleSymmetry::Symmetry& symm,
         const std::string& k_file_name,
         const int& nspin,
         const ModuleBase::Matrix3& reciprocal_vec,
@@ -204,7 +205,8 @@ class K_Vectors
      * @note If the k-points type is Line mode and the symmetry flag is 1, it will quit with a warning.
      * @note If the number of k-points is greater than 100000, it will quit with a warning.
      */
-    bool read_kpoints(const std::string& fn); // return 0: something wrong.
+    bool read_kpoints(const UnitCell& ucell,
+                      const std::string& fn); // return 0: something wrong.
 
     /**
      * @brief Adds k-points linearly between special points.