Merge branch 'develop' into fft5

Author: A-006

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)
@@ -1389,7 +1390,7 @@ These variables are used to control the geometry relaxation.
 
 - **Type**: Integer
 - **Description**: The maximal number of ionic iteration steps, the minimum value is 1.
-- **Default**: 1
+- **Default**: 1 for SCF, 50 for relax and cell-relax calcualtions
 
 ### relax_cg_thr
 
@@ -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

@@ -228,7 +228,6 @@ OBJS_ELECSTAT=elecstate.o\
     pot_xc.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
-      elecstate_lcao_tddft.o\
       elecstate_lcao_cal_tau.o\
       density_matrix.o\
       density_matrix_io.o\
@@ -246,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\
-      lcao_gets.o\
+      esolver_gets.o\
       lcao_others.o\
-      lcao_init_after_vc.o\
 
 OBJS_GINT=gint.o\
       gint_gamma_env.o\
@@ -318,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\
@@ -670,7 +664,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\
@@ -681,7 +675,9 @@ OBJS_VDW=vdw.o\
     vdwd3_parameters.o\
     vdwd2.o\
     vdwd3.o\
-    vdwd3_parameters_tab.o
+    vdwd3_parameters_tab.o\
+    vdwd3_autoset_xcname.o\
+    vdwd3_autoset_xcparam.o
 
 OBJS_DFTU=dftu.o\
       dftu_force.o\
@@ -693,14 +689,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\
 
@@ -727,7 +720,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_run.cpp

@@ -67,13 +67,16 @@ void Driver::driver_run() {
         Relax_Driver rl_driver;
         rl_driver.relax_driver(p_esolver);
     }
+    else if (cal_type == "get_S")
+    {
+        p_esolver->runner(0, GlobalC::ucell);
+    }
     else
     {
         //! supported "other" functions:
         //! get_pchg(LCAO),
         //! test_memory(PW,LCAO),
         //! test_neighbour(LCAO),
-        //! get_S(LCAO),
         //! gen_bessel(PW), et al.
         const int istep = 0;
         p_esolver->others(istep);

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)

source/module_base/vector3.h

@@ -36,6 +36,9 @@ template <class T> class Vector3
     Vector3(const Vector3<T> &v) : x(v.x), y(v.y), z(v.z){}; // Peize Lin add 2018-07-16   
     explicit Vector3(const std::array<T,3> &v) :x(v[0]), y(v[1]), z(v[2]){}
 
+    template <typename U>
+    explicit Vector3(const Vector3<U>& other) : x(static_cast<T>(other.x)), y(static_cast<T>(other.y)), z(static_cast<T>(other.z)) {}
+
     Vector3(Vector3<T> &&v) noexcept : x(v.x), y(v.y), z(v.z) {}
 
     /**

source/module_basis/module_nao/two_center_bundle.cpp

@@ -63,7 +63,7 @@ void TwoCenterBundle::tabulate()
 {
     ModuleBase::SphericalBesselTransformer sbt(true);
     orb_->set_transformer(sbt);
-    beta_->set_transformer(sbt);
+    if (beta_) { beta_->set_transformer(sbt); }
     if (alpha_) {
         alpha_->set_transformer(sbt);
 }
@@ -75,22 +75,17 @@ void TwoCenterBundle::tabulate()
     //              build two-center integration tables
     //================================================================
     // set up a universal radial grid
-    double rmax = std::max(orb_->rcut_max(), beta_->rcut_max());
-    if (alpha_) {
-        rmax = std::max(rmax, alpha_->rcut_max());
-}
+    double rmax = orb_->rcut_max();
+    if (beta_) { rmax = std::max(rmax, beta_->rcut_max()); }
+    if (alpha_) { rmax = std::max(rmax, alpha_->rcut_max()); }
     double dr = 0.01;
     double cutoff = 2.0 * rmax;
     int nr = static_cast<int>(rmax / dr) + 1;
 
     orb_->set_uniform_grid(true, nr, cutoff, 'i', true);
-    beta_->set_uniform_grid(true, nr, cutoff, 'i', true);
-    if (alpha_) {
-        alpha_->set_uniform_grid(true, nr, cutoff, 'i', true);
-}
-    if (orb_onsite_) {
-        orb_onsite_->set_uniform_grid(true, nr, cutoff, 'i', true);
-}
+    if (beta_) { beta_->set_uniform_grid(true, nr, cutoff, 'i', true); }
+    if (alpha_) { alpha_->set_uniform_grid(true, nr, cutoff, 'i', true);}
+    if (orb_onsite_) {  orb_onsite_->set_uniform_grid(true, nr, cutoff, 'i', true);}
 
     // build TwoCenterIntegrator objects
     kinetic_orb = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
@@ -101,9 +96,12 @@ void TwoCenterBundle::tabulate()
     overlap_orb->tabulate(*orb_, *orb_, 'S', nr, cutoff);
     ModuleBase::Memory::record("TwoCenterTable: Overlap", overlap_orb->table_memory());
 
-    overlap_orb_beta = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
-    overlap_orb_beta->tabulate(*orb_, *beta_, 'S', nr, cutoff);
-    ModuleBase::Memory::record("TwoCenterTable: Nonlocal", overlap_orb_beta->table_memory());
+    if (beta_)
+    {
+        overlap_orb_beta = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
+        overlap_orb_beta->tabulate(*orb_, *beta_, 'S', nr, cutoff);
+        ModuleBase::Memory::record("TwoCenterTable: Nonlocal", overlap_orb_beta->table_memory());
+    }
 
     if (alpha_)
     {