Merge upstream develop into hsolver refactor

Author: zxy.monado

.github/workflows/build_test_cmake.yml

@@ -24,7 +24,7 @@ jobs:
 
           - tag: gnu
             external_toolchain_args: ""
-            build_args: "-DENABLE_LIBXC=ON -DENABLE_MLALGO=ON -DENABLE_LIBRI=ON -DENABLE_DFTD4=ON"
+            build_args: "-DENABLE_LIBXC=ON -DENABLE_MLALGO=ON -DENABLE_LIBRI=ON -DENABLE_DFTD4=ON -DENABLE_PEXSI=ON"
             name: "Build extra components with GNU toolchain"
           - tag: intel
             external_toolchain_args: "--with-intel"

.github/workflows/test.yml

@@ -67,6 +67,7 @@ jobs:
             -DENABLE_RAPIDJSON=ON \
             -DENABLE_FLOAT_FFTW=ON \
             -DENABLE_DFTD4=ON \
+            -DENABLE_PEXSI=ON \
             -Werror=dev
 
 # Temporarily removed because no one maintains this now.

AGENTS.md

@@ -34,6 +34,10 @@ rules. Read the complete governance document before making or reviewing changes:
   is required.
 - Report the exact verification performed. Do not claim completion without
   fresh test or check output.
+- For multi-step refactors (e.g., splitting a large `.cpp` into several
+  files), build and commit after each step rather than batching all changes
+  before verification. This keeps the blast radius small when a step
+  surfaces a missing include or instantiation error.
 - Prefer `std::vector` over raw `new`/`delete` for dynamic arrays; before
   converting class members, confirm no external code consumes them as raw
   pointers (e.g., `std::vector<bool>` has no `.data()`), and use
@@ -99,13 +103,18 @@ rules. Read the complete governance document before making or reviewing changes:
 - Member -> free function: inventory `this->` reads; pass as params (const
   for config, ref for mutable state); move only when body is `this`-free;
   keep thin wrapper; compile each step.
+- Extract a base-class nested-vector member in three steps (hold + forward,
+  switch writers, delete legacy) so no commit mixes old-storage writes with
+  new-storage reads.
 
 ## Local Commands
 
 ```bash
 python3 tools/03_code_analysis/agent_governance_check.py --staged
 python3 tools/03_code_analysis/agent_governance_check.py --base upstream/develop --head HEAD --format text
 pre-commit run abacus-agent-governance --all-files
+# Score changed C++ files for quality debt (pass line is 60):
+python3 tools/03_code_analysis/code_quality_score.py $(git diff --name-only upstream/develop...HEAD | grep -E '\.(cpp|h)$')
 ```
 
 The repository text files have been normalized to LF once. Day-to-day line

docs/advanced/acceleration/cuda.md

@@ -44,7 +44,7 @@ The ABACUS program will automatically determine whether the current ELPA support
 ## Run with the GPU support by editing the INPUT script:
 
 In `INPUT` file we need to set the input parameter [device](../input_files/input-main.md#device) to `gpu`. If this parameter is not set, ABACUS will try to determine if there are available GPUs.
-- Set `ks_solver`: For the PW basis, CG, BPCG and Davidson methods are supported on GPU; set the input parameter [ks_solver](../input_files/input-main.md#ks_solver) to `cg`, `bpcg` or `dav`. For the LCAO basis, `cusolver`, `cusolvermp` and `elpa` is supported on GPU.
+- Set `ks_solver`: For the PW basis, CG, BPCG, Davidson, and Davidson subspace methods are supported on GPU; set the input parameter [ks_solver](../input_files/input-main.md#ks_solver) to `cg`, `bpcg`, `dav`, or `dav_subspace`. For the LCAO basis, `cusolver`, `cusolvermp`, and `elpa` are supported on GPU.
 - **single-card**: ABACUS allows for single-GPU acceleration. You can run ABACUS without any MPI process by command `abacus`, and `ks_solver cusolver` is recommended for the LCAO basis. *note: avoid using `mpirun -n 1 abacus`*.
 - **multi-cards**: ABACUS allows for multi-GPU acceleration. If you have multiple GPU cards, you can run ABACUS with several MPI processes, and each process will utilize one GPU card. For example, the command `mpirun -n 2 abacus` will by default launch two GPUs for computation. If you only have one card, this command will only start one GPU. *note: the number of MPI processes SHOULD be equal to the number of GPU cards, unless you are using MPS in your computer.*
 
@@ -58,4 +58,4 @@ PW basis:
 
 LCAO basis:
 - Unless there is a specific reason, avoid using multiple GPUs, as it can be slower than using a single GPU. This is because the generalized eigenvalue solution of the LCAO basis set will incur additional communication overhead when calculated on multiple cards. When the memory limit of a GPU card makes it insufficient to complete the task, it is recommended to use multiple cards for calculation.
-- When using elpa on GPUs, some ELPA internal logs will be output.
+- When using elpa on GPUs, some ELPA internal logs will be output.

docs/advanced/input_files/input-main.md

@@ -8,6 +8,7 @@
   - [System variables](#system-variables)
     - [suffix](#suffix)
     - [ntype](#ntype)
+    - [cell\_replica](#cell_replica)
     - [calculation](#calculation)
     - [esolver\_type](#esolver_type)
     - [symmetry](#symmetry)
@@ -356,6 +357,8 @@
     - [md\_restart](#md_restart)
     - [md\_restartfreq](#md_restartfreq)
     - [md\_dumpfreq](#md_dumpfreq)
+    - [md\_neighbor\_skin](#md_neighbor_skin)
+    - [md\_out\_force](#md_out_force)
     - [dump\_force](#dump_force)
     - [dump\_vel](#dump_vel)
     - [dump\_virial](#dump_virial)
@@ -586,6 +589,14 @@
   - [Reduced Density Matrix Functional Theory](#reduced-density-matrix-functional-theory)
     - [rdmft](#rdmft)
     - [rdmft\_power\_alpha](#rdmft_power_alpha)
+  - [Density functional perturbation theory](#density-functional-perturbation-theory)
+    - [dfpt\_qmesh](#dfpt_qmesh)
+    - [dfpt\_qfile](#dfpt_qfile)
+    - [dfpt\_compute\_q0](#dfpt_compute_q0)
+    - [dfpt\_loto](#dfpt_loto)
+    - [dfpt\_conv\_thr](#dfpt_conv_thr)
+    - [dfpt\_max\_iter](#dfpt_max_iter)
+    - [dfpt\_mix\_beta](#dfpt_mix_beta)
 
 ## System variables
 
@@ -601,6 +612,12 @@
 - **Description**: Number of different atom species in the calculation.
 - **Default**: 0
 
+### cell_replica
+
+- **Type**: Three Integers
+- **Description**: Replicate the input STRU by Na, Nb, and Nc along its lattice vectors for distributed MDCell workflows. This parameter is only used for classical potentials or machine-learned interatomic potentials. The default is 1 1 1, which preserves the input structure.
+- **Default**: 1 1 1
+
 ### calculation
 
 - **Type**: String
@@ -634,6 +651,7 @@
   - nep: Neuroevolution Potential
   - ks-lr: Kohn-Sham density functional theory + LR-TDDFT (Under Development Feature)
   - lr: LR-TDDFT with given KS orbitals (Under Development Feature)
+  - dfpt: density functional perturbation theory (Under Development Feature)
 - **Default**: ksdft
 
 ### symmetry
@@ -3504,15 +3522,28 @@
 ### md_restartfreq
 
 - **Type**: Integer
-- **Description**: The output frequency of OUT.{suffix}/STRIU/, which are used to restart molecular dynamics calculations, see md_restart in detail.
+- **Description**: The output frequency of OUT.{suffix}/STRU_MD_*, which are used to restart molecular dynamics calculations, see md_restart in detail. Set to 0 to disable MD restart output.
 - **Default**: 5
 
 ### md_dumpfreq
 
 - **Type**: Integer
-- **Description**: The output frequency of OUT.${suffix}/MD_dump in molecular dynamics calculations, which including the information of lattices and atoms.
+- **Description**: The output frequency of OUT.${suffix}/MD_dump in molecular dynamics calculations, which includes lattice and atomic information. Set to 0 to disable MD_dump output.
 - **Default**: 1
 
+### md_neighbor_skin
+
+- **Type**: Real
+- **Description**: The extra neighbor-list radius in Angstrom for MDCell molecular dynamics. This parameter is only used for classical potentials or machine-learned interatomic potentials. A positive value reuses the cutoff-plus-skin candidate list until an atom has moved by half this distance; 0 rebuilds the list every force evaluation.
+- **Default**: 0.0
+- **Unit**: Angstrom
+
+### md_out_force
+
+- **Type**: Boolean
+- **Description**: Whether to output the TOTAL-FORCE table in OUT.${suffix}/running_md.log for MDCell molecular dynamics. This does not affect force calculation or molecular dynamics integration.
+- **Default**: True
+
 ### dump_force
 
 - **Type**: Boolean
@@ -3535,8 +3566,8 @@
 
 - **Type**: Integer
 - **Description**: The random seed to initialize random numbers used in molecular dynamics calculations.
-  - < 0: No srand() function is called.
-  - >= 0: The function srand(md_seed) is called.
+  - < 0: Each MPI rank uses the default seed 1 plus its rank.
+  - >= 0: Each MPI rank uses md_seed plus its rank.
 - **Default**: -1
 
 ### md_tfreq
@@ -5170,3 +5201,49 @@
 - **Default**: 0.656
 
 [back to top](#full-list-of-input-keywords)
+
+## Density functional perturbation theory
+
+### dfpt_qmesh
+
+- **Type**: Vector of Int (1 or 3 values)
+- **Description**: Set the Monkhorst-Pack q mesh (gamma-centered) for DFPT phonon calculations. The q mesh must be commensurate with the ground-state k mesh: k + q must be a point of the k list (modulo a reciprocal lattice vector). For example, a 4x4x4 KPT mesh is commensurate with dfpt_qmesh values of 1, 2, or 4 along each direction. This parameter is ignored when dfpt_qfile is set.
+- **Default**: 1 1 1
+
+### dfpt_qfile
+
+- **Type**: String
+- **Description**: Set the file containing the q points for DFPT, in the same format as the KPT file (Q_POINTS card: Gamma/Monkhorst-Pack mesh, or an explicit Direct/Cartesian list; symmetry reduction is not applied to file q lists). When set, it overrides dfpt_qmesh. Each q point must still be commensurate with the ground-state k mesh.
+- **Default**: ""
+
+### dfpt_compute_q0
+
+- **Type**: Boolean
+- **Description**: Whether to compute the macroscopic dielectric tensor (epsilon_inf) and the Born effective charges at q = 0 within the same DFPT run. Requires a q point at Gamma (the default dfpt_qmesh 1 1 1).
+- **Default**: false
+
+### dfpt_loto
+
+- **Type**: Boolean
+- **Description**: Whether to apply the Lyddane-Sachs-Teller non-analytic correction to the Gamma-point dynamical matrix, which splits the longitudinal and transverse optical modes. Requires dfpt_compute_q0 to be true, since the correction is built from epsilon_inf and the Born effective charges.
+- **Default**: false
+
+### dfpt_conv_thr
+
+- **Type**: Real
+- **Description**: Set the convergence threshold of the self-consistent DFPT cycle: the iteration stops when the relative residual of the first-order density ||drho_out - drho_in|| / ||drho_out|| drops below this value for every displacement.
+- **Default**: 1.0e-8
+
+### dfpt_max_iter
+
+- **Type**: Integer
+- **Description**: Set the maximum number of self-consistent DFPT iterations for each atomic displacement.
+- **Default**: 100
+
+### dfpt_mix_beta
+
+- **Type**: Real
+- **Description**: Set the plain-mixing coefficient of the first-order density in the self-consistent DFPT cycle. The response Jacobian has strongly negative eigenvalues on the smallest-G shells (Coulomb stiffness), so beta must stay below 2 / (1 + |lambda_min|); the default 0.4 keeps margin up to |lambda_min| ~ 3. A larger value accelerates convergence for weakly screened systems but may diverge.
+- **Default**: 0.4
+
+[back to top](#full-list-of-input-keywords)

docs/parameters.yaml

@@ -17,6 +17,14 @@ parameters:
     default_value: "0"
     unit: ""
     availability: ""
+  - name: cell_replica
+    category: System variables
+    type: Three Integers
+    description: |
+      Replicate the input STRU by Na, Nb, and Nc along its lattice vectors for distributed MDCell workflows. This parameter is only used for classical potentials or machine-learned interatomic potentials. The default is 1 1 1, which preserves the input structure.
+    default_value: 1 1 1
+    unit: ""
+    availability: ""
   - name: calculation
     category: System variables
     type: String
@@ -53,6 +61,7 @@ parameters:
       * nep: Neuroevolution Potential
       * ks-lr: Kohn-Sham density functional theory + LR-TDDFT (Under Development Feature)
       * lr: LR-TDDFT with given KS orbitals (Under Development Feature)
+      * dfpt: density functional perturbation theory (Under Development Feature)
     default_value: ksdft
     unit: ""
     availability: ""
@@ -1449,18 +1458,34 @@ parameters:
     category: Molecular dynamics
     type: Integer
     description: |
-      The output frequency of OUT.{suffix}/STRIU/, which are used to restart molecular dynamics calculations, see md_restart in detail.
+      The output frequency of OUT.{suffix}/STRU_MD_*, which are used to restart molecular dynamics calculations, see md_restart in detail. Set to 0 to disable MD restart output.
     default_value: "5"
     unit: ""
     availability: ""
   - name: md_dumpfreq
     category: Molecular dynamics
     type: Integer
     description: |
-      The output frequency of OUT.${suffix}/MD_dump in molecular dynamics calculations, which including the information of lattices and atoms.
+      The output frequency of OUT.${suffix}/MD_dump in molecular dynamics calculations, which includes lattice and atomic information. Set to 0 to disable MD_dump output.
     default_value: "1"
     unit: ""
     availability: ""
+  - name: md_neighbor_skin
+    category: Molecular dynamics
+    type: Real
+    description: |
+      The extra neighbor-list radius in Angstrom for MDCell molecular dynamics. This parameter is only used for classical potentials or machine-learned interatomic potentials. A positive value reuses the cutoff-plus-skin candidate list until an atom has moved by half this distance; 0 rebuilds the list every force evaluation.
+    default_value: "0.0"
+    unit: Angstrom
+    availability: ""
+  - name: md_out_force
+    category: Molecular dynamics
+    type: Boolean
+    description: |
+      Whether to output the TOTAL-FORCE table in OUT.${suffix}/running_md.log for MDCell molecular dynamics. This does not affect force calculation or molecular dynamics integration.
+    default_value: "True"
+    unit: ""
+    availability: ""
   - name: dump_force
     category: Molecular dynamics
     type: Boolean
@@ -1490,8 +1515,8 @@ parameters:
     type: Integer
     description: |
       The random seed to initialize random numbers used in molecular dynamics calculations.
-      * < 0: No srand() function is called.
-      * >= 0: The function srand(md_seed) is called.
+      * < 0: Each MPI rank uses the default seed 1 plus its rank.
+      * >= 0: Each MPI rank uses md_seed plus its rank.
     default_value: "-1"
     unit: ""
     availability: ""
@@ -3022,6 +3047,62 @@ parameters:
     default_value: "-1 2 -1 2"
     unit: primitive cells
     availability: "lr_solver==plot and exciton_plot_format in [slice, both]"
+  - name: dfpt_qmesh
+    category: Density functional perturbation theory
+    type: Vector of Int (1 or 3 values)
+    description: |
+      Set the Monkhorst-Pack q mesh (gamma-centered) for DFPT phonon calculations. The q mesh must be commensurate with the ground-state k mesh: k + q must be a point of the k list (modulo a reciprocal lattice vector). For example, a 4x4x4 KPT mesh is commensurate with dfpt_qmesh values of 1, 2, or 4 along each direction. This parameter is ignored when dfpt_qfile is set.
+    default_value: 1 1 1
+    unit: ""
+    availability: ""
+  - name: dfpt_qfile
+    category: Density functional perturbation theory
+    type: String
+    description: |
+      Set the file containing the q points for DFPT, in the same format as the KPT file (Q_POINTS card: Gamma/Monkhorst-Pack mesh, or an explicit Direct/Cartesian list; symmetry reduction is not applied to file q lists). When set, it overrides dfpt_qmesh. Each q point must still be commensurate with the ground-state k mesh.
+    default_value: "\"\""
+    unit: ""
+    availability: ""
+  - name: dfpt_compute_q0
+    category: Density functional perturbation theory
+    type: Boolean
+    description: |
+      Whether to compute the macroscopic dielectric tensor (epsilon_inf) and the Born effective charges at q = 0 within the same DFPT run. Requires a q point at Gamma (the default dfpt_qmesh 1 1 1).
+    default_value: "false"
+    unit: ""
+    availability: ""
+  - name: dfpt_loto
+    category: Density functional perturbation theory
+    type: Boolean
+    description: |
+      Whether to apply the Lyddane-Sachs-Teller non-analytic correction to the Gamma-point dynamical matrix, which splits the longitudinal and transverse optical modes. Requires dfpt_compute_q0 to be true, since the correction is built from epsilon_inf and the Born effective charges.
+    default_value: "false"
+    unit: ""
+    availability: ""
+  - name: dfpt_conv_thr
+    category: Density functional perturbation theory
+    type: Real
+    description: |
+      Set the convergence threshold of the self-consistent DFPT cycle: the iteration stops when the relative residual of the first-order density ||drho_out - drho_in|| / ||drho_out|| drops below this value for every displacement.
+    default_value: "1.0e-8"
+    unit: ""
+    availability: ""
+  - name: dfpt_max_iter
+    category: Density functional perturbation theory
+    type: Integer
+    description: |
+      Set the maximum number of self-consistent DFPT iterations for each atomic displacement.
+    default_value: "100"
+    unit: ""
+    availability: ""
+  - name: dfpt_mix_beta
+    category: Density functional perturbation theory
+    type: Real
+    description: |
+      Set the plain-mixing coefficient of the first-order density in the self-consistent DFPT cycle. The response Jacobian has strongly negative eigenvalues on the smallest-G shells (Coulomb stiffness), so beta must stay below 2 / (1 + |lambda_min|); the default 0.4 keeps margin up to |lambda_min| ~ 3. A larger value accelerates convergence for weakly screened systems but may diverge.
+    default_value: "0.4"
+    unit: ""
+    availability: ""
   - name: out_freq_ion
     category: Output information
     type: Integer

generate_build_info.sh

@@ -257,7 +257,8 @@ if [ "${USE_CUDA}" == "ON" ]; then
 fi
 # --- Final File Generation ---
 
-INPUT_FILE="source_io/build_info.h.in"
+SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+INPUT_FILE="${SCRIPT_DIR}/source/source_io/build_info.h.in"
 
 # Use sed to replace all placeholders with detected values
 # Note the use of different delimiters (#) for paths to avoid conflicts with /

python/pyabacus/CONTRIBUTING.md

@@ -375,7 +375,7 @@ To contribute to the `pyabacus` project, follow these steps:
 
        @staticmethod
        def sphbes_zeros(l: int, n: int, zeros: NDArray[np.float64]) -> None: 
-           super().sphbes_zeros(l, n, zeros)
+           _Sphbes.sphbes_zeros(l, n, zeros)
    ```
 
 ## Conclusion

python/pyabacus/src/pyabacus/ModuleBase/_module_base.py

@@ -48,7 +48,7 @@ def dsphbesj(self, *args, **kwargs):
 
     @staticmethod
     def sphbes_zeros(l: int, n: int, zeros: NDArray[np.float64]) -> None: 
-        super().sphbes_zeros(l, n, zeros)
+        _Sphbes.sphbes_zeros(l, n, zeros)
 
 class Integral(_Integral):
     def __init__(self) -> None: 
@@ -81,7 +81,7 @@ def Simpson_Integral_0toall(
         rab: NDArray[np.float64], 
         asum: NDArray[np.float64]
     ) -> None: 
-        super().Simpson_Integral_0toall(mesh, func, rab, asum)
+        _Integral.Simpson_Integral_0toall(mesh, func, rab, asum)
 
     @staticmethod
     def Simpson_Integral_alltoinf(
@@ -90,7 +90,7 @@ def Simpson_Integral_alltoinf(
         rab: NDArray[np.float64], 
         asum: NDArray[np.float64]
     ) -> None: 
-        super().Simpson_Integral_alltoinf(mesh, func, rab, asum)
+        _Integral.Simpson_Integral_alltoinf(mesh, func, rab, asum)
 
     @overload
     @staticmethod