Merge branch 'develop' of https://github.com/deepmodeling/abacus-develop into hotfix

Author: Qianruipku

docs/advanced/input_files/input-main.md

@@ -423,6 +423,7 @@
     - [pexsi\_zero\_thr](#pexsi_zero_thr)
   - [Linear Response TDDFT](#linear-response-tddft)
     - [xc\_kernel](#xc_kernel)
+    - [lr\_init\_xc\_kernel](#lr_init_xc_kernel)
     - [lr\_solver](#lr_solver)
     - [lr\_thr](#lr_thr)
     - [nocc](#nocc)
@@ -988,7 +989,7 @@ calculations.
 
 - **Type**: String
 - **Description**: In our package, the XC functional can either be set explicitly using the `dft_functional` keyword in `INPUT` file. If `dft_functional` is not specified, ABACUS will use the xc functional indicated in the pseudopotential file.
-  On the other hand, if dft_functional is specified, it will overwrite the functional from pseudopotentials and performs calculation with whichever functional the user prefers. We further offer two ways of supplying exchange-correlation functional. The first is using 'short-hand' names such as 'LDA', 'PBE', 'SCAN'. A complete list of 'short-hand' expressions can be found in [the source code](../../../source/module_hamilt_general/module_xc/xc_functional.cpp). The other way is only available when ***compiling with LIBXC***, and it allows for supplying exchange-correlation functionals as combinations of LIBXC keywords for functional components, joined by a plus sign, for example, 'dft_functional='LDA_X_1D_EXPONENTIAL+LDA_C_1D_CSC'. The list of LIBXC keywords can be found on its [website](https://www.tddft.org/programs/libxc/functionals/). In this way, **we support all the LDA,GGA and mGGA functionals provided by LIBXC**.
+  On the other hand, if dft_functional is specified, it will overwrite the functional from pseudopotentials and performs calculation with whichever functional the user prefers. We further offer two ways of supplying exchange-correlation functional. The first is using 'short-hand' names such as 'LDA', 'PBE', 'SCAN'. A complete list of 'short-hand' expressions can be found in [the source code](../../../source/module_hamilt_general/module_xc/xc_functional.cpp). The other way is only available when ***compiling with LIBXC***, and it allows for supplying exchange-correlation functionals as combinations of LIBXC keywords for functional components, joined by a plus sign, for example, dft_functional='LDA_X_1D_EXPONENTIAL+LDA_C_1D_CSC'. The list of LIBXC keywords can be found on its [website](https://libxc.gitlab.io/functionals/). In this way, **we support all the LDA,GGA and mGGA functionals provided by LIBXC**.
 
   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. And forces for hybrid functionals are not supported yet.
 
@@ -1759,14 +1760,14 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital basis
-- **Description**: Whether to print Hamiltonian matrices H(R)/density matrics DM(R) in npz format. This feature does not work for gamma-only calculations. Currently only intended for internal usage.
+- **Description**: Whether to print Hamiltonian matrices $H(R)$/density matrics $DM(R)$ in npz format. This feature does not work for gamma-only calculations. Currently only intended for internal usage.
 - **Default**: False
 
 ### dm_to_rho
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital basis
-- **Description**: Reads density matrix DM(R) in npz format and creates electron density on grids. This feature does not work for gamma-only calculations. Only supports serial calculations. Currently only intended for internal usage.
+- **Description**: Reads density matrix $DM(R)$ in npz format and creates electron density on grids. This feature does not work for gamma-only calculations. Only supports serial calculations. Currently only intended for internal usage.
 - **Default**: False
 
 ### out_app_flag
@@ -2925,7 +2926,7 @@ These variables are used to control DFT+U correlated parameters
 
   - where $\gamma$ is a parameter that adjusts the relative weight of the error function to the derivative error function.
 - **Unit**: Bohr
-- **Default**: 5.0
+- **Default**: 3.0
 
 [back to top](#full-list-of-input-keywords)
 
@@ -3943,6 +3944,15 @@ These parameters are used to solve the excited states using. e.g. LR-TDDFT.
 Currently supported: `RPA`, `LDA`, `PBE`, `HSE`, `HF`.
 - **Default**: LDA
 
+### lr_init_xc_kernel
+
+- **Type**: String
+- **Description**: The method to initalize the xc kernel. 
+  - "default": Calculate xc kerenel ($f_\text{xc}$) from the ground-state charge density.
+  - "file": Read the xc kernel $f_\text{xc}$ on grid from the provided files. The following words should be the paths of ".cube" files, where the first 1 (*[nspin](#nspin)==1*) or 3 (*[nspin](#nspin)==2*, namely spin-aa, spin-ab and spin-bb) will be read in. The parameter [xc_kernel](#xc_kernel) will be invalid. Now only LDA-type kernel is supproted as the potential will be calculated by directly multiplying the transition density.
+  - "from_charge_file": Calculate fxc from the charge density read from the provided files. The following words should be the paths of ".cube" files, where the first [nspin]($nspin) files will be read in. 
+- **Default**: "default"
+
 ### lr_solver
 
 - **Type**: String

examples/bsse/water/result.ref

@@ -1,5 +1,5 @@
--13.31798074740440
-E_H2O: -466.1225988776397
-E_O: -427.6271689751553
-E_H1: -12.58872469295076
-E_H2: -12.58872446212924
+-13.49968292248493
+E_H2O: -466.1225988772539
+E_O: -427.5222287307378
+E_H1: -12.55034372743879
+E_H2: -12.55034349659238

examples/scf/lcao_Cu/INPUT

@@ -4,7 +4,6 @@ orbital_dir                 ../../../tests/PP_ORB
 nbands		10
 
 calculation	scf
-ecutwfc		100 
 ecutwfc		100 ###Energy cutoff needs to be tested to ensure your calculation is reliable.[1]
 scf_thr		1.0e-8
 scf_nmax	100

source/Makefile.Objects

@@ -245,18 +245,13 @@ OBJS_ESOLVER=esolver.o\
     esolver_of.o\
     esolver_of_tool.o\
     esolver_of_interface.o\
-    pw_init_after_vc.o\
-    pw_init_globalc.o\
     pw_others.o\
 
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       esolver_ks_lcao_tddft.o\
-      dpks_cal_e_delta_band.o\
-      set_matrix_grid.o\
       lcao_before_scf.o\
       esolver_gets.o\
       lcao_others.o\
-      lcao_init_after_vc.o\
 
 OBJS_GINT=gint.o\
       gint_gamma_env.o\
@@ -317,7 +312,7 @@ OBJS_HAMILT_LCAO=hamilt_lcao.o\
     op_dftu_lcao.o\
     deepks_lcao.o\
     op_exx_lcao.o\
-    sc_lambda_lcao.o\
+    dspin_lcao.o\
     dftu_lcao.o\
 
 OBJS_HCONTAINER=base_matrix.o\
@@ -671,7 +666,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     symmetry_rhog.o\
     wavefunc.o\
     wf_atomic.o\
-    wfinit.o\
+    psiinit.o\
     elecond.o\
     sto_tool.o\
     sto_elecond.o\
@@ -696,14 +691,11 @@ OBJS_DFTU=dftu.o\
       dftu_hamilt.o
 
 OBJS_DELTASPIN=basic_funcs.o\
-      cal_h_lambda.o\
       cal_mw_from_lambda.o\
-      cal_mw_helper.o\
       cal_mw.o\
       init_sc.o\
       lambda_loop_helper.o\
       lambda_loop.o\
-      sc_parse_json.o\
       spin_constrain.o\
       template_helpers.o\
 
@@ -730,7 +722,7 @@ OBJS_TENSOR=tensor.o\
     dmr_complex.o\
     operator_lr_hxc.o\
     operator_lr_exx.o\
-    kernel_xc.o\
+    xc_kernel.o\
     pot_hxc_lrtd.o\
     lr_spectrum.o\
     hamilt_casida.o\

source/driver.cpp

@@ -183,7 +183,7 @@ void Driver::atomic_world()
     //--------------------------------------------------
 
     // where the actual stuff is done
-    this->driver_run();
+    this->driver_run(GlobalC::ucell);
 
     ModuleBase::timer::finish(GlobalV::ofs_running);
     ModuleBase::Memory::print_all(GlobalV::ofs_running);

source/driver.h

@@ -1,6 +1,8 @@
 #ifndef DRIVER_H
 #define DRIVER_H
 
+#include "module_cell/unitcell.h"
+
 class Driver
 {
   public:
@@ -34,7 +36,7 @@ class Driver
     void atomic_world();
 
     // the actual calculations
-    void driver_run();
+    void driver_run(UnitCell& ucell);
 };
 
 #endif

source/driver_run.cpp

@@ -24,7 +24,8 @@
  * the configuration-changing subroutine takes force and stress and updates the
  * configuration
  */
-void Driver::driver_run() {
+void Driver::driver_run(UnitCell& ucell)
+{
     ModuleBase::TITLE("Driver", "driver_line");
     ModuleBase::timer::tick("Driver", "driver_line");
 
@@ -39,37 +40,35 @@ void Driver::driver_run() {
 #endif
 
     // the life of ucell should begin here, mohan 2024-05-12
-    // delete ucell as a GlobalC in near future
-    GlobalC::ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
-    Check_Atomic_Stru::check_atomic_stru(GlobalC::ucell,
-                                         PARAM.inp.min_dist_coef);
+    ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
+    Check_Atomic_Stru::check_atomic_stru(ucell, PARAM.inp.min_dist_coef);
 
     //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`)
-    ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, GlobalC::ucell);
+    ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell);
 
     //! 3: initialize Esolver and fill json-structure
-    p_esolver->before_all_runners(PARAM.inp, GlobalC::ucell);
+    p_esolver->before_all_runners(ucell, PARAM.inp);
 
     // this Json part should be moved to before_all_runners, mohan 2024-05-12
 #ifdef __RAPIDJSON
-    Json::gen_stru_wrapper(&GlobalC::ucell);
+    Json::gen_stru_wrapper(&ucell);
 #endif
 
     const std::string cal_type = PARAM.inp.calculation;
 
     //! 4: different types of calculations
     if (cal_type == "md")
     {
-        Run_MD::md_line(GlobalC::ucell, p_esolver, PARAM);
+        Run_MD::md_line(ucell, p_esolver, PARAM);
     }
     else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf")
     {
         Relax_Driver rl_driver;
-        rl_driver.relax_driver(p_esolver);
+        rl_driver.relax_driver(p_esolver, ucell);
     }
     else if (cal_type == "get_S")
     {
-        p_esolver->runner(0, GlobalC::ucell);
+        p_esolver->runner(ucell, 0);
     }
     else
     {
@@ -79,11 +78,11 @@ void Driver::driver_run() {
         //! test_neighbour(LCAO),
         //! gen_bessel(PW), et al.
         const int istep = 0;
-        p_esolver->others(istep);
+        p_esolver->others(ucell, istep);
     }
 
     //! 5: clean up esolver
-    p_esolver->after_all_runners();
+    p_esolver->after_all_runners(ucell);
 
     ModuleESolver::clean_esolver(p_esolver);
 

source/module_base/module_mixing/mixing_data.cpp

@@ -3,7 +3,7 @@
 namespace Base_Mixing
 {
 
-Mixing_Data::Mixing_Data(const int& ndim, const int& length, const size_t& type_size)
+Mixing_Data::Mixing_Data(const int& ndim, const std::size_t& length, const size_t& type_size)
 {
     this->ndim_tot = ndim;
     this->length = length;
@@ -15,16 +15,18 @@ Mixing_Data::Mixing_Data(const int& ndim, const int& length, const size_t& type_
 
 Mixing_Data::~Mixing_Data()
 {
-    if (this->data != nullptr)
+    if (this->data != nullptr) {
         free(this->data);
 }
+}
 
-void Mixing_Data::resize(const int& ndim, const int& length, const size_t& type_size)
+void Mixing_Data::resize(const int& ndim, const std::size_t& length, const size_t& type_size)
 {
     this->ndim_tot = ndim;
     this->length = length;
-    if (this->data != nullptr)
+    if (this->data != nullptr) {
         free(this->data);
+}
     if (ndim * length > 0)
     {
         this->data = malloc(ndim * length * type_size);

source/module_base/module_mixing/mixing_data.h

@@ -22,7 +22,7 @@ class Mixing_Data
      * @param type_size size of type
      *
      */
-    Mixing_Data(const int& ndim, const int& length, const size_t& type_size);
+    Mixing_Data(const int& ndim, const std::size_t& length, const size_t& type_size);
 
     /**
      * @brief Destroy the Mixing_Data object
@@ -38,7 +38,7 @@ class Mixing_Data
      * @param type_size size of type
      *
      */
-    void resize(const int& ndim, const int& length, const size_t& type_size);
+    void resize(const int& ndim, const std::size_t& length, const size_t& type_size);
 
     /**
      * @brief push data to the tensor
@@ -54,7 +54,7 @@ class Mixing_Data
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 4096/sizeof(FPTYPE))
 #endif
-        for (int i = 0; i < length; ++i)
+        for (std::size_t i = 0; i < length; ++i)
         {
             FP_startdata[i] = data_in[i];
         }
@@ -86,7 +86,7 @@ class Mixing_Data
     // the number of vectors for mixing
     int ndim_tot = 0;
     // the length of each vector
-    int length = 0;
+    std::size_t length = 0;
     // the start index for vector: start = this->index_move(0)
     int start = -1;
     // the number of used vectors for mixing

source/module_base/module_mixing/pulay_mixing.cpp

@@ -28,13 +28,13 @@ void Pulay_Mixing::tem_push_data(Mixing_Data& mdata,
                                  std::function<void(FPTYPE*, const FPTYPE*, const FPTYPE*)> mix,
                                  const bool& need_calcoef)
 {
-    const size_t length = mdata.length;
+    const std::size_t length = mdata.length;
     std::vector<FPTYPE> F_tmp(length);
 
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
 #endif
-    for (int i = 0; i < length; ++i)
+    for (std::size_t i = 0; i < length; ++i)
     {
         F_tmp[i] = data_out[i] - data_in[i];
     }
@@ -69,7 +69,7 @@ void Pulay_Mixing::tem_push_data(Mixing_Data& mdata,
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
 #endif
-        for (int i = 0; i < length; ++i)
+        for (std::size_t i = 0; i < length; ++i)
         {
             FP_F[i] = F_tmp[i];
         }
@@ -81,7 +81,7 @@ void Pulay_Mixing::tem_push_data(Mixing_Data& mdata,
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 4096 / sizeof(FPTYPE))
 #endif
-        for (int i = 0; i < length; ++i)
+        for (std::size_t i = 0; i < length; ++i)
         {
             FP_startF[i] = F_tmp[i];
         }
@@ -103,7 +103,7 @@ void Pulay_Mixing::tem_cal_coef(const Mixing_Data& mdata, std::function<double(F
         ModuleBase::WARNING_QUIT(
             "Pulay_mixing",
             "One Pulay_Mixing object can only bind one Mixing_Data object to calculate coefficients");
-    const int length = mdata.length;
+    const std::size_t length = mdata.length;
     FPTYPE* FP_F = static_cast<FPTYPE*>(F);
 
     if (mdata.ndim_use > 1)