Feature: add out_mat_ds to print <phi|dphi> matrix

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_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_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_io/output_mat_sparse.cpp

@@ -9,6 +9,7 @@ namespace ModuleIO
 template <>
 void output_mat_sparse(const bool& out_mat_hsR,
                        const bool& out_mat_dh,
+                       const bool& out_mat_ds,
                        const bool& out_mat_t,
                        const bool& out_mat_r,
                        const int& istep,
@@ -27,6 +28,7 @@ void output_mat_sparse(const bool& out_mat_hsR,
 template <>
 void output_mat_sparse(const bool& out_mat_hsR,
                        const bool& out_mat_dh,
+                       const bool& out_mat_ds,
                        const bool& out_mat_t,
                        const bool& out_mat_r,
                        const int& istep,
@@ -68,6 +70,18 @@ void output_mat_sparse(const bool& out_mat_hsR,
                    orb,
                    kv); // LiuXh add 2019-07-15
     }
+    //! generate a file containing the derivatives of the overlap matrix (in Ry/Bohr)
+    if (out_mat_ds)
+    {
+        output_dSR(istep,
+                   ucell,
+                   pv,
+                   HS_Arrays,
+                   grid, // mohan add 2024-04-06
+                   two_center_bundle,
+                   orb,
+                   kv);
+    }
 
     // add by jingan for out r_R matrix 2019.8.14
     if (out_mat_r)

source/module_io/output_mat_sparse.h

@@ -13,6 +13,7 @@ namespace ModuleIO
 template <typename T>
 void output_mat_sparse(const bool& out_mat_hsR,
                        const bool& out_mat_dh,
+                       const bool& out_mat_ds,
                        const bool& out_mat_t,
                        const bool& out_mat_r,
                        const int& istep,

source/module_io/read_input.cpp

@@ -182,7 +182,7 @@ void ReadInput::create_directory(const Parameter& param)
     //----------------------------------------------------------
     bool out_dir = false;
     if (!param.input.out_app_flag
-        && (param.input.out_mat_hs2 || param.input.out_mat_r || param.input.out_mat_t || param.input.out_mat_dh))
+        && (param.input.out_mat_hs2 || param.input.out_mat_r || param.input.out_mat_t || param.input.out_mat_dh || param.input.out_mat_ds))
     {
         out_dir = true;
     }

source/module_io/read_input_item_output.cpp

@@ -322,6 +322,18 @@ void ReadInput::item_output()
         };
         this->add_item(item);
     }
+    {
+        Input_Item item("out_mat_ds");
+        item.annotation = "output of derivative of S(R) matrix";
+        read_sync_bool(input.out_mat_ds);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.out_mat_ds && para.input.nspin == 4)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_ds is not available for nspin = 4");
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("out_mat_xc");
         item.annotation = "output exchange-correlation matrix in KS-orbital representation";

source/module_io/write_HS_R.cpp

@@ -83,6 +83,37 @@ void ModuleIO::output_HSR(const UnitCell& ucell,
     return;
 }
 
+void ModuleIO::output_dSR(const int& istep,
+    const UnitCell& ucell,
+    const Parallel_Orbitals& pv,
+    LCAO_HS_Arrays& HS_Arrays,
+    const Grid_Driver& grid, // mohan add 2024-04-06
+    const TwoCenterBundle& two_center_bundle,
+    const LCAO_Orbitals& orb,
+    const K_Vectors& kv,
+    const bool& binary,
+    const double& sparse_thr)
+{
+ModuleBase::TITLE("ModuleIO", "output_dSR");
+ModuleBase::timer::tick("ModuleIO", "output_dSR");
+
+
+sparse_format::cal_dS(ucell,
+            pv,
+            HS_Arrays,
+            grid,
+            two_center_bundle,
+            orb,
+            sparse_thr);
+// mohan update 2024-04-01
+ModuleIO::save_dH_sparse(istep, pv, HS_Arrays, sparse_thr, binary, "S");
+
+sparse_format::destroy_dH_R_sparse(HS_Arrays);
+
+ModuleBase::timer::tick("ModuleIO", "output_dSR");
+return;
+}
+
 void ModuleIO::output_dHR(const int& istep,
                           const ModuleBase::matrix& v_eff,
                           Gint_k& gint_k, // mohan add 2024-04-01

source/module_io/write_HS_R.h

@@ -42,6 +42,17 @@ void output_dHR(const int& istep,
                 const bool& binary = false,
                 const double& sparse_threshold = 1e-10);
 
+void output_dSR(const int& istep,
+                const UnitCell& ucell,
+                const Parallel_Orbitals& pv,
+                LCAO_HS_Arrays& HS_Arrays,
+                const Grid_Driver& grid, // mohan add 2024-04-06
+                const TwoCenterBundle& two_center_bundle,
+                const LCAO_Orbitals& orb,
+                const K_Vectors& kv,
+                const bool& binary = false,
+                const double& sparse_thr = 1e-10);
+
 void output_TR(const int istep,
                const UnitCell& ucell,
                const Parallel_Orbitals& pv,

source/module_io/write_HS_sparse.cpp

@@ -328,7 +328,8 @@ void ModuleIO::save_dH_sparse(const int& istep,
                               const Parallel_Orbitals& pv,
                               LCAO_HS_Arrays& HS_Arrays,
                               const double& sparse_thr,
-                              const bool& binary) {
+                              const bool& binary,
+                              const std::string& fileflag) {
     ModuleBase::TITLE("ModuleIO", "save_dH_sparse");
     ModuleBase::timer::tick("ModuleIO", "save_dH_sparse");
 
@@ -429,24 +430,24 @@ void ModuleIO::save_dH_sparse(const int& istep,
     std::stringstream sshz[2];
     if (PARAM.inp.calculation == "md" && !PARAM.inp.out_app_flag) {
         sshx[0] << PARAM.globalv.global_matrix_dir << step << "_"
-                << "data-dHRx-sparse_SPIN0.csr";
+                << "data-d"<<fileflag<<"Rx-sparse_SPIN0.csr";
         sshx[1] << PARAM.globalv.global_matrix_dir << step << "_"
-                << "data-dHRx-sparse_SPIN1.csr";
+                << "data-d"<<fileflag<<"Rx-sparse_SPIN1.csr";
         sshy[0] << PARAM.globalv.global_matrix_dir << step << "_"
-                << "data-dHRy-sparse_SPIN0.csr";
+                << "data-d"<<fileflag<<"Ry-sparse_SPIN0.csr";
         sshy[1] << PARAM.globalv.global_matrix_dir << step << "_"
-                << "data-dHRy-sparse_SPIN1.csr";
+                << "data-d"<<fileflag<<"Ry-sparse_SPIN1.csr";
         sshz[0] << PARAM.globalv.global_matrix_dir << step << "_"
-                << "data-dHRz-sparse_SPIN0.csr";
+                << "data-d"<<fileflag<<"Rz-sparse_SPIN0.csr";
         sshz[1] << PARAM.globalv.global_matrix_dir << step << "_"
-                << "data-dHRz-sparse_SPIN1.csr";
+                << "data-d"<<fileflag<<"Rz-sparse_SPIN1.csr";
     } else {
-        sshx[0] << PARAM.globalv.global_out_dir << "data-dHRx-sparse_SPIN0.csr";
-        sshx[1] << PARAM.globalv.global_out_dir << "data-dHRx-sparse_SPIN1.csr";
-        sshy[0] << PARAM.globalv.global_out_dir << "data-dHRy-sparse_SPIN0.csr";
-        sshy[1] << PARAM.globalv.global_out_dir << "data-dHRy-sparse_SPIN1.csr";
-        sshz[0] << PARAM.globalv.global_out_dir << "data-dHRz-sparse_SPIN0.csr";
-        sshz[1] << PARAM.globalv.global_out_dir << "data-dHRz-sparse_SPIN1.csr";
+        sshx[0] << PARAM.globalv.global_out_dir << "data-d"<<fileflag<<"Rx-sparse_SPIN0.csr";
+        sshx[1] << PARAM.globalv.global_out_dir << "data-d"<<fileflag<<"Rx-sparse_SPIN1.csr";
+        sshy[0] << PARAM.globalv.global_out_dir << "data-d"<<fileflag<<"Ry-sparse_SPIN0.csr";
+        sshy[1] << PARAM.globalv.global_out_dir << "data-d"<<fileflag<<"Ry-sparse_SPIN1.csr";
+        sshz[0] << PARAM.globalv.global_out_dir << "data-d"<<fileflag<<"Rz-sparse_SPIN0.csr";
+        sshz[1] << PARAM.globalv.global_out_dir << "data-d"<<fileflag<<"Rz-sparse_SPIN1.csr";
     }
     std::ofstream g1x[2];
     std::ofstream g1y[2];

source/module_io/write_HS_sparse.h

@@ -25,7 +25,8 @@ void save_dH_sparse(const int& istep,
                     const Parallel_Orbitals& pv,
                     LCAO_HS_Arrays& HS_Arrays,
                     const double& sparse_thr,
-                    const bool& binary);
+                    const bool& binary,
+                    const std::string& fileflag = "H");
 
 template <typename Tdata>
 void save_sparse(const std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, Tdata>>>& smat,

source/module_parameter/input_parameter.h

@@ -380,6 +380,7 @@ struct Input_para
     bool out_mat_hs2 = false;             ///< LiuXh add 2019-07-16, output H(R) matrix and
                                           ///< S(R) matrix in local basis.
     bool out_mat_dh = false;
+    bool out_mat_ds = false;
     bool out_mat_xc = false;      ///< output exchange-correlation matrix in
                                   ///< KS-orbital representation.
     bool out_mat_xc2 = false;     ///< output exchange-correlation matrix Vxc(R) in NAO representation.

tests/integrate/207_NO_KP_OS/INPUT

@@ -8,3 +8,4 @@ ecutwfc 		10
 scf_thr 		1e-6 
 basis_type 		lcao
 calculation 		get_S 
+out_mat_ds  1