Merge branch 'develop' of https://gitee.com/deepmodeling/abacus-develop into rdmft

Author: JGHan7

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)
@@ -2928,7 +2929,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)
 
@@ -3946,6 +3947,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

@@ -247,7 +247,6 @@ 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\
 
@@ -258,7 +257,6 @@ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       lcao_before_scf.o\
       esolver_gets.o\
       lcao_others.o\
-      lcao_init_after_vc.o\
 
 OBJS_GINT=gint.o\
       gint_gamma_env.o\
@@ -319,7 +317,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\
@@ -698,14 +696,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\
 
@@ -732,7 +727,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/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_cell/atom_spec.cpp

@@ -13,55 +13,20 @@ Atom::Atom()
     type = 0;
     stapos_wf = 0;
     mass = 0.0;
-    tau = new ModuleBase::Vector3<double>[1];
-    dis = new ModuleBase::Vector3<double>[1];
-    taud = new ModuleBase::Vector3<double>[1];
-    vel = new ModuleBase::Vector3<double>[1];
-    mag = new double[1];
-    angle1 = new double[1];
-    angle2 = new double[1];
-    m_loc_ = new ModuleBase::Vector3<double>[1];
-    l_nchi = new int[1];
-    iw2l = new int[1];
-    iw2n = new int[1];
-    iw2m = new int[1];
-    iw2_ylm = new int[1];
-    iw2_new = new bool[1];
-    mbl = new ModuleBase::Vector3<int>[1];
 }
 
 Atom::~Atom()
 {
-    delete[] tau;
-    delete[] dis;
-    delete[] taud;
-    delete[] vel;
-    delete[] mag;
-    delete[] angle1;
-    delete[] angle2;
-    delete[] m_loc_;
-    delete[] l_nchi;
-    delete[] iw2l;
-    delete[] iw2n;
-    delete[] iw2m;
-    delete[] iw2_ylm;
-    delete[] iw2_new;
-    delete[] mbl;
 }
 
 void Atom::set_index(void)
 {
     assert(nw != 0);
-    delete[] iw2l;
-    delete[] iw2n;
-    delete[] iw2m;
-    delete[] iw2_ylm;
-    delete[] iw2_new;
-    iw2l = new int[nw];
-    iw2n = new int[nw];
-    iw2m = new int[nw];
-    iw2_ylm = new int[nw];
-    iw2_new = new bool[nw];
+    this->iw2l.resize(nw, 0);
+    this->iw2n.resize(nw, 0);
+    this->iw2m.resize(nw, 0);
+    this->iw2_ylm.resize(nw, 0);
+    this->iw2_new.resize(nw, false); // bool array to check if the local orbital is new
 
     int index = 0;
     for (int L = 0; L <= nwl; L++)
@@ -103,7 +68,7 @@ void Atom::print_Atom(std::ofstream& ofs)
     ModuleBase::GlobalFunc::OUT(ofs, "mass", mass);
     ofs << std::endl;
 
-    output::printv31_d(ofs, "atom_position(cartesian)", tau, na);
+    output::printv31_d(ofs, "atom_position(cartesian)", tau.data(), na);
     /*
     for (int i = 0;i < na;i++)
     {
@@ -131,34 +96,26 @@ void Atom::bcast_atom(void)
     Parallel_Common::bcast_bool(coulomb_potential);
     if (GlobalV::MY_RANK != 0)
     {
-        delete[] l_nchi;
-        l_nchi = new int[nwl + 1];
+        this->l_nchi.resize(nwl + 1, 0);
     }
-    Parallel_Common::bcast_int(l_nchi, nwl + 1);
+    Parallel_Common::bcast_int(l_nchi.data(), nwl + 1);
     Parallel_Common::bcast_bool(this->flag_empty_element);
     Parallel_Common::bcast_double(mass);
 
     if (GlobalV::MY_RANK != 0)
     {
         assert(na != 0);
-        delete[] tau;
-        delete[] dis;
-        delete[] taud;
-        delete[] vel;
-        delete[] mag;
-        delete[] angle1;
-        delete[] angle2;
-        delete[] m_loc_;
-        delete[] mbl;
-        tau = new ModuleBase::Vector3<double>[na];
-        dis = new ModuleBase::Vector3<double>[na];
-        taud = new ModuleBase::Vector3<double>[na];
-        vel = new ModuleBase::Vector3<double>[na];
-        mag = new double[na];
-        angle1 = new double[na];
-        angle2 = new double[na];
-        m_loc_ = new ModuleBase::Vector3<double>[na];
-        mbl = new ModuleBase::Vector3<int>[na];
+        this->tau.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->dis.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->taud.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->vel.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->mag.resize(na, 0);
+        this->angle1.resize(na, 0);
+        this->angle2.resize(na, 0);
+        this->m_loc_.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->mbl.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
+        this->lambda.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
+        this->constrain.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
     }
 
     for (int i = 0; i < na; i++)
@@ -184,6 +141,12 @@ void Atom::bcast_atom(void)
         Parallel_Common::bcast_int(mbl[i].x);
         Parallel_Common::bcast_int(mbl[i].y);
         Parallel_Common::bcast_int(mbl[i].z);
+        Parallel_Common::bcast_double(lambda[i].x);
+        Parallel_Common::bcast_double(lambda[i].y);
+        Parallel_Common::bcast_double(lambda[i].z);
+        Parallel_Common::bcast_int(constrain[i].x);
+        Parallel_Common::bcast_int(constrain[i].y);
+        Parallel_Common::bcast_int(constrain[i].z);
     }
 
     return;

source/module_cell/atom_spec.h

@@ -13,14 +13,14 @@ class Atom
 
     Atom_pseudo ncpp;
     double mass; // the mass of atom
-    ModuleBase::Vector3<int> *mbl; // whether the atoms can move or not
+    std::vector<ModuleBase::Vector3<int>> mbl; // whether the atoms can move or not
 	bool flag_empty_element = false; // whether is the empty element for bsse.	Peize Lin add 2021.04.07
 
-    int *iw2m; // use iw to find m
-    int *iw2n; // use iw to find n
-    int *iw2l; // use iw to find L
-	int *iw2_ylm;
-	bool *iw2_new;
+    std::vector<int> iw2m; // use iw to find m
+    std::vector<int> iw2n; // use iw to find n
+    std::vector<int> iw2l; // use iw to find L
+	std::vector<int> iw2_ylm;
+	std::vector<bool> iw2_new;
     int nw; // number of local orbitals (l,n,m) of this type
 
     void set_index(void);
@@ -30,21 +30,23 @@ class Atom
 
     int nwl; // max L(Angular momentum) (for local basis)
     double Rcut; //pengfei Li 16-2-29
-    int *l_nchi; // number of chi for each L
+    std::vector<int> l_nchi; // number of chi for each L
     int stapos_wf; // start position of wave functions
 
     std::string label; // atomic symbol
-    ModuleBase::Vector3<double> *tau;// Cartesian coordinates of each atom in this type.
-    ModuleBase::Vector3<double> *dis;// direct displacements of each atom in this type in current step  liuyu modift 2023-03-22
-	ModuleBase::Vector3<double> *taud;// Direct coordinates of each atom in this type.
-    ModuleBase::Vector3<double> *vel;// velocities of each atom in this type.
-    ModuleBase::Vector3<double> *force; // force acting on each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> tau;// Cartesian coordinates of each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> dis;// direct displacements of each atom in this type in current step  liuyu modift 2023-03-22
+	std::vector<ModuleBase::Vector3<double>> taud;// Direct coordinates of each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> vel;// velocities of each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> force; // force acting on each atom in this type.
+    std::vector<ModuleBase::Vector3<double>> lambda; // Lagrange multiplier for each atom in this type. used in deltaspin
+    std::vector<ModuleBase::Vector3<int>> constrain; // constrain for each atom in this type. used in deltaspin
     std::string label_orb; // atomic Element symbol in the orbital file of lcao
 
-	double* mag;
-	double* angle1;//spin angle, added by zhengdy-soc
-	double* angle2;
-    ModuleBase::Vector3<double> *m_loc_;
+	std::vector<double> mag;
+	std::vector<double> angle1;//spin angle, added by zhengdy-soc
+	std::vector<double> angle2;
+    std::vector<ModuleBase::Vector3<double>> m_loc_;
     // Coulomb potential v(r) = z/r
     // It is a local potentail, and has no non-local potential parts.
     bool coulomb_potential = false;