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

* initial version of writing H and dH terms (except exx)

* correct 1-electron terms: T, Vl, Vnl

* fix pot register

* Gint_drho

* add EXX H&dH, dH (H-F term) of Veff Hartree, and refactor

* add dm_container_to_Ds and exx nscf-from-dm workflow; minor refactor OperatorEXX

* fix parallel segfaults and kinetic sign

* XC H-F term (FD)

* enable value-adding HContainers, total dH and atom-specified output; enable nspin=2

* dH test cases

* exclude cal_exx_dHs with __EXX_DEV flag, waiting for LibRI PR merged

* small fixes

* change scf_in_dmr ref due to not renormalizing charge at init_chg==dm

* fix conflict #7462 #7475 #7478

* Revert reference change by adding option dm_no_renormalize to init_chg

This reverts commit ed14d9f.

* fix pyabacus CI

Author: maki49

CMakeLists.txt

@@ -35,6 +35,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)
@@ -651,6 +652,9 @@ if(ENABLE_LIBRI)
   endif()
   abacus_add_feature_definitions(__EXX EXX_DM=3 EXX_H_COMM=2 TEST_EXX_LCAO=0
                                  TEST_EXX_RADIAL=1)
+  if(EXX_DEV)
+    abacus_add_feature_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)
@@ -1869,7 +1881,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
@@ -1930,7 +1942,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;
@@ -1944,6 +1956,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)
@@ -1978,27 +2038,90 @@
 
 - **Type**: Boolean \[Integer\](optional)
 - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)*
-- **Description**: Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. 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 optional second parameter controls text output precision. 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
 - **Unit**: Ry
 
 ### out_mat_dh
 
-- **Type**: Boolean \[Integer\](optional)
-- **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)*
-- **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix. The optional second parameter controls text output precision. 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.
-
+- **Type**: Integer
+- **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 optional second parameter controls text output precision. 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 optional second parameter controls text output precision. 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

python/pyabacus/pyproject.toml

@@ -30,7 +30,7 @@ test = ["contourpy","cycler","exceptiongroup","fonttools","importlib-metadata","
 
 [tool.scikit-build]
 wheel.expand-macos-universal-tags = true
-cmake.verbose = true
+build.verbose = true
 logging.level = "INFO"
 
 [tool.scikit-build.cmake.define]

source/source_esolver/esolver_ks_lcao.cpp

@@ -180,7 +180,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     if(istep == 0)//if the first scf step, readin DMR from file,
     {
         //calculate or readin the density matrix DMR
-        if(PARAM.inp.init_chg == "dm")
+        if(PARAM.inp.init_chg == "dm" || PARAM.inp.init_chg == "dm_no_renormalize")
         {
             //! 13.1.1) init charge density from density matrix file
             LCAO_domain::init_chg_dm<TK>(PARAM.globalv.global_readin_dir, PARAM.inp.nspin,
@@ -546,6 +546,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_no_renormalize")
+    {
+        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

@@ -179,6 +179,11 @@ class Potential : public PotBase
     {
         return this->ucell_;
     }
+    // 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

@@ -81,10 +81,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/rr_sparse_writer.cpp