Feature: term-separated LCAO Hamiltonian atomic derivatives (dH)

CMakeLists.txt

@@ -32,6 +32,7 @@ option(ENABLE_MLALGO "Enable the machine learning algorithms" OFF)
 option(ENABLE_LCAO "Enable LCAO algorithm" ON)
 option(USE_ELPA "Enable ELPA for LCAO" ON)
 option(ENABLE_LIBRI "Enable LibRI for hybrid functional" OFF)
+option(EXX_DEV "Enable LibRI developing features" OFF)
 option(ENABLE_LIBCOMM "Enable LibComm" OFF)
 option(ENABLE_PEXSI "Enable PEXSI for LCAO" OFF)
 option(ENABLE_DFTD4 "Enable DFT-D4 dispersion correction" OFF)
@@ -759,6 +760,9 @@ if(ENABLE_LIBRI)
   target_link_libraries(${ABACUS_BIN_NAME} ri module_exx_symmetry)
   add_compile_definitions(__EXX EXX_DM=3 EXX_H_COMM=2 TEST_EXX_LCAO=0
                           TEST_EXX_RADIAL=1)
+  if(EXX_DEV)
+   add_compile_definitions(__EXX_DEV)
+  endif()
 endif()
 
 if(ENABLE_LIBRI OR DEFINED LIBCOMM_DIR)

docs/advanced/input_files/input-main.md

@@ -166,11 +166,23 @@
     - [out\_stru](#out_stru)
     - [out\_level](#out_level)
     - [out\_mat\_hs](#out_mat_hs)
+    - [out\_mat\_h\_t](#out_mat_h_t)
+    - [out\_mat\_h\_vl](#out_mat_h_vl)
+    - [out\_mat\_h\_vnl](#out_mat_h_vnl)
+    - [out\_mat\_h\_vh](#out_mat_h_vh)
+    - [out\_mat\_h\_vxc](#out_mat_h_vxc)
+    - [out\_mat\_h\_exx](#out_mat_h_exx)
     - [out\_mat\_hs2](#out_mat_hs2)
     - [out\_mat\_tk](#out_mat_tk)
     - [out\_mat\_r](#out_mat_r)
     - [out\_mat\_t](#out_mat_t)
     - [out\_mat\_dh](#out_mat_dh)
+    - [out\_mat\_dh\_t](#out_mat_dh_t)
+    - [out\_mat\_dh\_vl](#out_mat_dh_vl)
+    - [out\_mat\_dh\_vnl](#out_mat_dh_vnl)
+    - [out\_mat\_dh\_vh](#out_mat_dh_vh)
+    - [out\_mat\_dh\_vxc](#out_mat_dh_vxc)
+    - [out\_mat\_dh\_exx](#out_mat_dh_exx)
     - [out\_mat\_ds](#out_mat_ds)
     - [out\_mat\_xc](#out_mat_xc)
     - [out\_mat\_xc2](#out_mat_xc2)
@@ -1865,7 +1877,7 @@
 
   The corresponding sequence of the orbitals can be seen in Basis Set.
 
-  Also controled by out_freq_ion and out_app_flag.
+  Also controlled by out_freq_ion and out_app_flag.
 
   > Note: In the 3.10-LTS version, the file names are WFC_NAO_GAMMA1_ION1.txt and WFC_NAO_K1_ION1.txt, etc.
 - **Default**: 0
@@ -1926,7 +1938,7 @@
 
 - **Type**: Boolean \[Integer\](optional)
 - **Availability**: *Numerical atomic orbital basis*
-- **Description**: Whether to print the upper triangular part of the Hamiltonian matrices and overlap matrices for each k-point into files in the directory OUT.${suffix}. The second number controls precision. For more information, please refer to hs_matrix.md. Also controled by out_freq_ion and out_app_flag.
+- **Description**: Whether to print the upper triangular part of the Hamiltonian matrices and overlap matrices for each k-point into files in the directory OUT.${suffix}. The second number controls precision. For more information, please refer to hs_matrix.md. Also controlled by out_freq_ion and out_app_flag.
   - For gamma only case:
    - nspin = 1: hks1_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix;
    - nspin = 2: hks1_nao.txt and hks2_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix. Note that the code will not output sks2_nao.txt because it is the same as sks1_nao.txt;
@@ -1940,6 +1952,54 @@
 - **Default**: False 8
 - **Unit**: Ry
 
+### out_mat_h_t
+
+- **Type**: Boolean \[Integer\](optional)
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print the kinetic energy matrix $T_{\mu\nu}(k) = \langle\phi_\mu|\hat{T}|\phi_\nu\rangle(k)$ for each k-point. The output format and file naming (e.g. `tks1_nao.txt`, `tks1k1_nao.txt`) follow [`out_mat_hs`](#out_mat_hs).
+- **Default**: False 8
+- **Unit**: Ry
+
+### out_mat_h_vl
+
+- **Type**: Boolean \[Integer\](optional)
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print the local pseudopotential matrix $V^L_{\mu\nu}(k) = \langle\phi_\mu|\hat{V}^L|\phi_\nu\rangle(k)$ for each k-point. The output format and file naming (e.g. `vlks1_nao.txt`, `vlks1k1_nao.txt`) follow [`out_mat_hs`](#out_mat_hs).
+- **Default**: False 8
+- **Unit**: Ry
+
+### out_mat_h_vnl
+
+- **Type**: Boolean \[Integer\](optional)
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print the nonlocal pseudopotential (Kleinman–Bylander) matrix $V^{NL}_{\mu\nu}(k) = \langle\phi_\mu|\hat{V}^{NL}|\phi_\nu\rangle(k)$ for each k-point. The output format and file naming (e.g. `vnlks1_nao.txt`, `vnlks1k1_nao.txt`) follow [`out_mat_hs`](#out_mat_hs).
+- **Default**: False 8
+- **Unit**: Ry
+
+### out_mat_h_vh
+
+- **Type**: Boolean \[Integer\](optional)
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print the Hartree matrix $V^H_{\mu\nu}(k) = \langle\phi_\mu|\hat{V}^H|\phi_\nu\rangle(k)$ for each k-point. The output format and file naming (e.g. `vhks1_nao.txt`, `vhks1k1_nao.txt`) follow [`out_mat_hs`](#out_mat_hs).
+- **Default**: False 8
+- **Unit**: Ry
+
+### out_mat_h_vxc
+
+- **Type**: Boolean \[Integer\](optional)
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print the exchange-correlation matrix $V^{XC}_{\mu\nu}(k) = \langle\phi_\mu|\hat{V}^{XC}|\phi_\nu\rangle(k)$ for each k-point. The output format and file naming (e.g. `vxcks1_nao.txt`, `vxcks1k1_nao.txt`) follow [`out_mat_hs`](#out_mat_hs).
+- **Default**: False 8
+- **Unit**: Ry
+
+### out_mat_h_exx
+
+- **Type**: Boolean \[Integer\](optional)
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4, hybrid functional only)*
+- **Description**: Whether to print the exact-exchange matrix $V^{EXX}_{\mu\nu}(k) = \langle\phi_\mu|\hat{V}^{EXX}|\phi_\nu\rangle(k)$ for each k-point. The output format and file naming (e.g. `vexxks1_nao.txt`, `vexxks1k1_nao.txt`) follow [`out_mat_hs`](#out_mat_hs). Requires a hybrid functional (`cal_exx = true`).
+- **Default**: False 8
+- **Unit**: Ry
+
 ### out_mat_hs2
 
 - **Type**: Boolean \[Integer\](optional)
@@ -1974,7 +2034,7 @@
 
 - **Type**: Boolean \[Integer\](optional)
 - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)*
-- **Description**: Generate files containing the kinetic energy matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag.
+- **Description**: Generate files containing the kinetic energy matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be trs1_nao.csr and so on. Also controlled by out_freq_ion and out_app_flag.
 
   > Note: In the 3.10-LTS version, the file name is data-TR-sparse_SPIN0.csr.
 - **Default**: False 8
@@ -1983,18 +2043,81 @@
 ### out_mat_dh
 
 - **Type**: Integer
-- **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)*
-- **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag.
-
+- **Availability**: *Numerical atomic orbital basis*
+- **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix $dH(k)/d\tau_I=d\braket{\phi|\hat{H}|\phi}(k)/d\tau_I$ where $\tau_I$ is the Ith atom position with the dense format as `out_mat_dh`. The names are dhk[x/y/z]_iat[I][_ik]_nao.txt.
+  - See also the term-separated output parameters: [`out_mat_dh_t`](#out_mat_dh_t), [`out_mat_dh_vl`](#out_mat_dh_vl), [`out_mat_dh_vnl`](#out_mat_dh_vnl), [`out_mat_dh_vh`](#out_mat_dh_vh), [`out_mat_dh_vxc`](#out_mat_dh_vxc) and [`out_mat_dh_exx`](#out_mat_dh_exx).
+  - If not gamma-only, also $\braket{\nabla\phi|\hat{H}\phi}(R)$ of sparse format as `out_mat_hs2` will also be output. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controlled by out_freq_ion and out_app_flag.
   > Note: In the 3.10-LTS version, the file name is data-dHRx-sparse_SPIN0.csr and so on.
+  - **Format**: `<enable> [precision] [iat1 iat2 ...]`
+  - The first value (0/1) enables or disables output.
+  - The second optional value sets the output precision (number of significant digits, default: 8).
+  - Starting from the third value, **1-based atom indices** can be listed to restrict the output to derivatives with respect to those specific atoms only. If no atom indices are given, derivatives are written for all atoms.
+
+  For example, `out_mat_dh 1 8 1 3` writes dH/dR for atoms 1 and 3 only (1-based indexing).
+
+- **Default**: 0 8
+- **Unit**: Ry/Bohr
+
+### out_mat_dh_t
+
+- **Type**: Integer
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print files containing the kinetic energy contribution to the Hamiltonian derivative, see [`out_mat_dh`](#out_mat_dh) for the same format. Output files: dhk[x/y/z]_iat[I][_ik]_nao.txt.
+
+- **Default**: 0 8
+- **Unit**: Ry/Bohr
+
+### out_mat_dh_vl
+
+- **Type**: Integer
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print files containing the local pseudopotential contribution to the Hamiltonian derivative, see [`out_mat_dh`](#out_mat_dh) for the same format. Output files: dvlk[x/y/z]_iat[I][_ik]_nao.txt.
+
+- **Default**: 0 8
+- **Unit**: Ry/Bohr
+
+### out_mat_dh_vnl
+
+- **Type**: Integer
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print files containing the nonlocal pseudopotential contribution to the Hamiltonian derivative,  see [`out_mat_dh`](#out_mat_dh) for the same format. Output files: dvnlk[x/y/z]_iat[I][_ik]_nao.txt.
+
+- **Default**: 0 8
+- **Unit**: Ry/Bohr
+
+### out_mat_dh_vh
+
+- **Type**: Integer
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print files containing the Hartree contribution to the Hamiltonian derivative,  see [`out_mat_dh`](#out_mat_dh) for the same format. Output files: dvhk[x/y/z]_iat[I][_ik]_nao.txt.
+
+- **Default**: 0 8
+- **Unit**: Ry/Bohr
+
+### out_mat_dh_vxc
+
+- **Type**: Integer
+- **Availability**: *Numerical atomic orbital basis (nspin ≠ 4)*
+- **Description**: Whether to print files containing the exchange-correlation contribution to the Hamiltonian derivative,  see [`out_mat_dh`](#out_mat_dh) for the same format. Output files: dvxck[x/y/z]_iat[I][_ik]_nao.txt.
+
+- **Default**: 0 8
+- **Unit**: Ry/Bohr
+
+### out_mat_dh_exx
+
+- **Type**: Integer
+- **Availability**: *Numerical atomic orbital basis, hybrid functional only (nspin ≠ 4)*
+    - Currently only availablewhen compiled with the personal developing branch of LibRI and -DEXX_DEV flag, waiting for the new release of LibRI to remove the flag.
+- **Description**: Whether to print files containing the exact-exchange contribution to the Hamiltonian derivative,  see [`out_mat_dh`](#out_mat_dh) for the same format. Output files: dvexxk[x/y/z]_iat[I][_ik]_nao.txt.
+
 - **Default**: 0 8
 - **Unit**: Ry/Bohr
 
 ### out_mat_ds
 
 - **Type**: Boolean \[Integer\](optional)
 - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)*
-- **Description**: Whether to print files containing the derivatives of the overlap matrix. The format will be the same as the overlap matrix as mentioned in out_mat_dh. The name of the files will be dsxrs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. This feature can be used with calculation get_s.
+- **Description**: Whether to print files containing the derivatives of the overlap matrix. The format will be the same as the overlap matrix as mentioned in out_mat_dh. The name of the files will be dsxrs1_nao.csr and so on. Also controlled by out_freq_ion and out_app_flag. This feature can be used with calculation get_s.
 
   > Note: In the 3.10-LTS version, the file name is data-dSRx-sparse_SPIN0.csr and so on.
 - **Default**: False 8

source/source_esolver/esolver_ks_lcao.cpp

@@ -544,6 +544,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
             PARAM.inp, this->kv, this->pelec, this->dmat.dm, this->pv,
             this->gd, this->psi, hamilt_lcao, this->dftu, this->two_center_bundle_,
             this->orb_, this->pw_wfc, this->pw_rho, this->pw_big, this->sf,
+            this->pw_rhod, this->locpp.vloc, this->solvent,
             this->rdmft_solver, this->deepks, this->exx_nao,
             this->conv_esolver, this->scf_nmax_flag, istep);
 

source/source_estate/init_scf.cpp

@@ -17,7 +17,10 @@ void init_scf(const UnitCell& ucell,
     pelec->charge->set_rho_core(ucell, strucfac, numeric);
 
     //! renormalize the charge density
-    pelec->charge->renormalize_rho();
+    if(PARAM.inp.init_chg != "dm")
+    {
+        pelec->charge->renormalize_rho();
+    }
 
     //! initialize the potential
     pelec->pot->init_pot(pelec->charge);

source/source_estate/module_charge/charge.cpp

@@ -54,7 +54,7 @@ void Charge::set_rhopw(ModulePW::PW_Basis* rhopw_in)
 }
 
 // mohan add 2025-12-02
-bool Charge::kin_density()
+bool Charge::kin_density() const
 {
 	if (XC_Functional::get_ked_flag() || PARAM.inp.out_elf[0] > 0)
 	{

source/source_estate/module_charge/charge.h

@@ -84,7 +84,7 @@ class Charge
                   const void* wfcpw = nullptr);
 
     // mohan add 2025-12-02
-    bool kin_density();
+    bool kin_density() const;
 
     void allocate(const int &nspin_in, const bool kin_den);
 

source/source_estate/module_pot/potential_new.h

@@ -175,6 +175,11 @@ class Potential : public PotBase
     {
         return this->rho_basis_;
     }
+    // get the local pseudopotential vloc(it, G) table; used by the dH module (out_mat_dh_vl)
+    const ModuleBase::matrix* get_vloc() const
+    {
+        return this->vloc_;
+    }
     // What about adding a function to get the wfc?
     // This is useful for the calculation of the exx energy
 

source/source_io/CMakeLists.txt

@@ -82,10 +82,13 @@ if(ENABLE_LCAO)
       module_mulliken/output_mulliken.cpp
       module_ml/io_npz.cpp
       module_hs/cal_pLpR.cpp
+      module_dhs/write_dH.cpp
+      module_dhs/write_dH_terms.cpp
   )
   list(APPEND objects_advanced	   
       module_unk/unk_overlap_lcao.cpp
       module_hs/write_HS_R.cpp
+      module_hs/write_H_terms.cpp
       module_hs/write_HS_sparse.cpp
       module_hs/single_R_io.cpp
       module_hs/cal_r_overlap_R.cpp