Merge branch 'develop' into fft_float2

Author: A-006

docs/advanced/input_files/input-main.md

@@ -153,6 +153,7 @@
     - [out\_mat\_hs2](#out_mat_hs2)
     - [out\_mat\_t](#out_mat_t)
     - [out\_mat\_dh](#out_mat_dh)
+    - [out\_mat\_ds](#out_mat_ds)
     - [out\_mat\_xc](#out_mat_xc)
     - [out\_mat\_xc2](#out_mat_xc2)
     - [out\_mat\_l](#out_mat_l)
@@ -1801,6 +1802,13 @@ These variables are used to control the output of properties.
 - **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix (in Ry/Bohr). The format will be the same as the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ as mentioned in [out_mat_hs2](#out_mat_hs2). The name of the files will be `data-dHRx-sparse_SPIN0.csr` and so on. Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
 - **Default**: False
 
+### out_mat_ds
+
+- **Type**: Boolean
+- **Availability**: Numerical atomic orbital basis (not gamma-only algorithm)
+- **Description**: Whether to print files containing the derivatives of the Overlap matrix (in Ry/Bohr). The format will be the same as the Overlap matrix $dH(R)$ as mentioned in [out_mat_dh](#out_mat_dh). The name of the files will be `data-dSRx-sparse_SPIN0.csr` and so on. Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag). This feature can be used with `calculation get_S`.
+- **Default**: False
+
 ### out_mat_xc
 
 - **Type**: Boolean

source/module_cell/klist.cpp

@@ -1413,12 +1413,6 @@ void K_Vectors::set_after_vc(const int& nspin_in,
         ModuleBase::Matrix3 RT = latvec.Transpose();
         for (int i = 0; i < nks; i++)
         {
-            //			std::cout << " ik=" << i
-            //				<< " kvec.x=" << kvec_c[i].x
-            //				<< " kvec.y=" << kvec_c[i].y
-            //				<< " kvec.z=" << kvec_c[i].z << std::endl;
-            // wrong!            kvec_d[i] = RT * kvec_c[i];
-            // mohan fixed bug 2011-03-07
             kvec_d[i] = kvec_c[i] * RT;
         }
         kd_done = true;

source/module_esolver/esolver_gets.cpp

@@ -137,6 +137,20 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
         r_matrix.out_rR(ucell, gd, istep);
     }
 
+    if (PARAM.inp.out_mat_ds)
+    {
+        LCAO_HS_Arrays HS_Arrays; // store sparse arrays
+        //! Print out sparse matrix
+        ModuleIO::output_dSR(istep,
+            ucell,
+            this->pv,
+            HS_Arrays,
+            gd, // mohan add 2024-04-06
+            two_center_bundle_,
+            orb_,
+            kv);
+    }
+
     ModuleBase::timer::tick("ESolver_GetS", "runner");
 }
 

source/module_esolver/lcao_after_scf.cpp

@@ -305,6 +305,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
         //! Print out sparse matrix
         ModuleIO::output_mat_sparse(PARAM.inp.out_mat_hs2,
                                     PARAM.inp.out_mat_dh,
+                                    PARAM.inp.out_mat_ds,
                                     PARAM.inp.out_mat_t,
                                     PARAM.inp.out_mat_r,
                                     istep,

source/module_esolver/lcao_others.cpp

@@ -345,10 +345,10 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     else if (cal_type == "get_wf")
     {
         std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting wave function");
-        IState_Envelope wavefunc(this->pelec);
+        Get_wf_lcao get_wf(this->pelec);
         if (PARAM.globalv.gamma_only_local)
         {
-            wavefunc.begin(ucell,
+            get_wf.begin(ucell,
                       this->psi,
                       this->pw_rhod,
                       this->pw_wfc,
@@ -370,7 +370,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
         }
         else
         {
-            wavefunc.begin(ucell,
+            get_wf.begin(ucell,
                       this->psi,
                       this->pw_rhod,
                       this->pw_wfc,

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp

@@ -68,16 +68,6 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(const Grid_Driv
 
     this->adjs_all.clear();
     this->adjs_all.reserve(this->ucell->nat);
-    bool pre_cal_nlm = false;
-    if (ucell->nat < 100) // less than 100 atom , cost memory for high performance
-    {                     // pre calculate nlm in initialization
-        this->nlm_tot.resize(ucell->nat);
-        pre_cal_nlm = true;
-    }
-    else
-    { // calculate nlm on the fly
-        this->nlm_tot.resize(1);
-    }
 
     for (int iat0 = 0; iat0 < ucell->nat; iat0++)
     {
@@ -135,10 +125,6 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(const Grid_Driv
                 // }
             }
         }
-        if (pre_cal_nlm)
-        {
-            this->pre_calculate_nlm(iat0, nlm_tot[iat0]);
-        }
     }
     // allocate the memory of BaseMatrix in HR, and set the new values to zero
     // if (std::is_same<TK, double>::value)
@@ -234,70 +220,15 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
 
 #ifdef __DEEPKS
 
-template <typename TK, typename TR>
-void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::pre_calculate_nlm(
-    const int iat0,
-    std::vector<std::unordered_map<int, std::vector<double>>>& nlm_in)
-{
-    const Parallel_Orbitals* paraV = this->hR->get_paraV();
-    const int npol = this->ucell->get_npol();
-    auto tau0 = ucell->get_tau(iat0);
-    int T0 = 0;
-    int I0 = 0;
-    ucell->iat2iait(iat0, &I0, &T0);
-    AdjacentAtomInfo& adjs = this->adjs_all[iat0];
-    nlm_in.resize(adjs.adj_num + 1);
-
-    for (int ad = 0; ad < adjs.adj_num + 1; ++ad)
-    {
-        const int T1 = adjs.ntype[ad];
-        const int I1 = adjs.natom[ad];
-        const int iat1 = ucell->itia2iat(T1, I1);
-        const ModuleBase::Vector3<double>& tau1 = adjs.adjacent_tau[ad];
-        const Atom* atom1 = &ucell->atoms[T1];
-
-        auto all_indexes = paraV->get_indexes_row(iat1);
-        auto col_indexes = paraV->get_indexes_col(iat1);
-        // insert col_indexes into all_indexes to get universal set with no repeat elements
-        all_indexes.insert(all_indexes.end(), col_indexes.begin(), col_indexes.end());
-        std::sort(all_indexes.begin(), all_indexes.end());
-        all_indexes.erase(std::unique(all_indexes.begin(), all_indexes.end()), all_indexes.end());
-        for (int iw1l = 0; iw1l < all_indexes.size(); iw1l += npol)
-        {
-            const int iw1 = all_indexes[iw1l] / npol;
-            std::vector<std::vector<double>> nlm;
-            // nlm is a vector of vectors, but size of outer vector is only 1 here
-            // If we are calculating force, we need also to store the gradient
-            // and size of outer vector is then 4
-            // inner loop : all projectors (L0,M0)
-
-            int L1 = atom1->iw2l[iw1];
-            int N1 = atom1->iw2n[iw1];
-            int m1 = atom1->iw2m[iw1];
-
-            // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l)
-            int M1 = (m1 % 2 == 0) ? -m1 / 2 : (m1 + 1) / 2;
-
-            ModuleBase::Vector3<double> dtau = tau0 - tau1;
-            intor_orb_alpha_->snap(T1, L1, N1, M1, 0, dtau * ucell->lat0, false /*calc_deri*/, nlm);
-            nlm_in[ad].insert({all_indexes[iw1l], nlm[0]});
-            if (npol == 2)
-            {
-                nlm_in[ad].insert({all_indexes[iw1l + 1], nlm[0]});
-            }
-        }
-    }
-}
-
 template <typename TK, typename TR>
 void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
 {
     ModuleBase::TITLE("DeePKS", "calculate_HR");
+    ModuleBase::timer::tick("DeePKS", "calculate_HR");
     if (this->V_delta_R->size_atom_pairs() == 0)
     {
         return;
     }
-    ModuleBase::timer::tick("DeePKS", "calculate_HR");
 
     const Parallel_Orbitals* paraV = this->V_delta_R->get_paraV();
     const int npol = this->ucell->get_npol();
@@ -361,18 +292,6 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
         const int trace_alpha_size = trace_alpha_row.size();
         //--------------------------------------------------
 
-        // if nlm_tot is not calculated already, calculate it on the fly now
-        std::vector<std::unordered_map<int, std::vector<double>>> nlm_on_the_fly;
-        const bool is_on_the_fly = (nlm_tot.size() != this->ucell->nat);
-
-        if (is_on_the_fly)
-        {
-            this->pre_calculate_nlm(iat0, nlm_on_the_fly);
-        }
-
-        std::vector<std::unordered_map<int, std::vector<double>>>& nlm_iat
-            = is_on_the_fly ? nlm_on_the_fly : nlm_tot[iat0];
-
         // 2. calculate <phi_I|beta>D<beta|phi_{J,R}> for each pair of <IJR> atoms
         for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1)
         {
@@ -386,11 +305,17 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
             {
                 continue;
             }
+            ModuleBase::Vector3<int> dR1(adjs.box[ad1].x, adjs.box[ad1].y, adjs.box[ad1].z);
+            if (this->ld->phialpha[0]->find_matrix(iat0, iat1, dR1.x, dR1.y, dR1.z) == nullptr)
+            {
+                continue;
+            }
 
             std::vector<double> s_1t(trace_alpha_size * row_size);
             for (int irow = 0; irow < row_size; irow++)
             {
-                const double* row_ptr = nlm_iat[ad1][row_indexes[irow]].data();
+                const hamilt::BaseMatrix<double>* overlap_1 = this->ld->phialpha[0]->find_matrix(iat0, iat1, dR1);
+                const double* row_ptr = overlap_1->get_pointer() + row_indexes[irow] * overlap_1->get_col_size();
                 double* ps1t = &s_1t[irow * trace_alpha_size];
                 for (int i = 0; i < trace_alpha_size; i++)
                 {
@@ -415,11 +340,23 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
                 }
                 auto col_indexes = paraV->get_indexes_col(iat2);
                 const int col_size = col_indexes.size();
+
+                if (col_size == 0)
+                {
+                    continue;
+                }
+                ModuleBase::Vector3<int> dR2(adjs.box[ad2].x, adjs.box[ad2].y, adjs.box[ad2].z);
+                if (this->ld->phialpha[0]->find_matrix(iat0, iat2, dR2.x, dR2.y, dR2.z) == nullptr)
+                {
+                    continue;
+                }
+
                 std::vector<double> hr_current(row_size * col_size, 0);
                 std::vector<double> s_2t(trace_alpha_size * col_size);
                 for (int icol = 0; icol < col_size; icol++)
                 {
-                    const double* col_ptr = nlm_iat[ad2][col_indexes[icol]].data();
+                    const hamilt::BaseMatrix<double>* overlap_2 = this->ld->phialpha[0]->find_matrix(iat0, iat2, dR2);
+                    const double* col_ptr = overlap_2->get_pointer() + col_indexes[icol] * overlap_2->get_col_size();
                     double* ps2t = &s_2t[icol * trace_alpha_size];
                     for (int i = 0; i < trace_alpha_size; i++)
                     {

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h

@@ -101,8 +101,6 @@ class DeePKS<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
      */
     void cal_HR_IJR(const double* hr_in, const int& row_size, const int& col_size, TR* data_pointer);
 
-    void pre_calculate_nlm(const int iat0, std::vector<std::unordered_map<int, std::vector<double>>>& nlm_in);
-    std::vector<std::vector<std::unordered_map<int, std::vector<double>>>> nlm_tot;
     /**
      * @brief initialize V_delta_R, search the nearest neighbor atoms
      * used for calculate the DeePKS real space Hamiltonian correction with specific <I,J,R> atom-pairs

source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp

@@ -4,6 +4,52 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
 #include <vector>
 
+void sparse_format::cal_dS(const UnitCell& ucell,
+    const Parallel_Orbitals& pv,
+    LCAO_HS_Arrays& HS_Arrays,
+    const Grid_Driver& grid,
+    const TwoCenterBundle& two_center_bundle,
+    const LCAO_Orbitals& orb,
+    const double& sparse_thr)
+{
+ModuleBase::TITLE("sparse_format", "cal_dS");
+
+sparse_format::set_R_range(HS_Arrays.all_R_coor, grid);
+const int nnr = pv.nnr;
+
+ForceStressArrays fsr_dh;
+fsr_dh.DHloc_fixedR_x = new double[nnr];
+fsr_dh.DHloc_fixedR_y = new double[nnr];
+fsr_dh.DHloc_fixedR_z = new double[nnr];
+ModuleBase::GlobalFunc::ZEROS(fsr_dh.DHloc_fixedR_x, nnr);
+ModuleBase::GlobalFunc::ZEROS(fsr_dh.DHloc_fixedR_y, nnr);
+ModuleBase::GlobalFunc::ZEROS(fsr_dh.DHloc_fixedR_z, nnr);
+// the pointers of dS is different from dH, use the dh pointers to reuse the print functions
+fsr_dh.DSloc_Rx = fsr_dh.DHloc_fixedR_x;
+fsr_dh.DSloc_Ry = fsr_dh.DHloc_fixedR_y;
+fsr_dh.DSloc_Rz = fsr_dh.DHloc_fixedR_z;
+// cal dS=<phi|dphi> in LCAO
+const bool cal_deri = true;
+const bool cal_stress = false;
+LCAO_domain::build_ST_new(fsr_dh,
+       'S',
+       cal_deri,
+       cal_stress,
+       ucell,
+       orb,
+       pv,
+       two_center_bundle,
+       &grid,
+       nullptr,
+       false); // delete unused parameter lm.Hloc_fixedR
+
+sparse_format::cal_dSTN_R(ucell,pv, HS_Arrays, fsr_dh, grid, orb.cutoffs(), 0, sparse_thr);
+delete[] fsr_dh.DHloc_fixedR_x;
+delete[] fsr_dh.DHloc_fixedR_y;
+delete[] fsr_dh.DHloc_fixedR_z;
+return;
+}
+
 void sparse_format::cal_dH(const UnitCell& ucell,
                            const Parallel_Orbitals& pv,
                            LCAO_HS_Arrays& HS_Arrays,

source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.h

@@ -21,6 +21,15 @@ void cal_dH(const UnitCell& ucell,
             const double& sparse_thr,
             Gint_k& gint_k);
 
+// calculated the derivative of the overlap matrix: <phi|dphi>
+void cal_dS(const UnitCell& ucell,
+            const Parallel_Orbitals& pv,
+            LCAO_HS_Arrays& HS_Arrays,
+            const Grid_Driver& grid,
+            const TwoCenterBundle& two_center_bundle,
+            const LCAO_Orbitals& orb,
+            const double& sparse_thr);
+
 // be called by 'cal_dH_sparse'
 void set_R_range(std::set<Abfs::Vector3_Order<int>>& all_R_coor, const Grid_Driver& grid);
 

source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp

@@ -64,6 +64,8 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
     {
         // this part is for integrated test of deepks
         // so it is printed no matter even if deepks_out_labels is not used
+        DeePKS_domain::update_dmr(kvec_d, dm->get_DMK_vector(), ucell, orb, *ParaV, GridD, dmr);
+        
         DeePKS_domain::cal_pdm<
             TK>(init_pdm, inlmax, lmaxd, inl2l, inl_index, kvec_d, dmr, phialpha, ucell, orb, GridD, *ParaV, pdm);