Merge upstream develop into source_base coverage branch

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

@@ -103,6 +103,9 @@ 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
 

docs/advanced/input_files/input-main.md

@@ -589,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
 
@@ -643,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
@@ -5192,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

@@ -61,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: ""
@@ -3046,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

source/Makefile.Objects

@@ -211,7 +211,7 @@ OBJS_CELL=atom_pseudo.o\
     print_cell.o\
     setup_nonlocal.o\
     klist.o\
-    k_vector_utils.o\
+    reciprocal_grid.o\
     cell_index.o\
     cell_tools.o\
     check_atomic_stru.o\
@@ -281,6 +281,7 @@ OBJS_ELECSTAT=elecstate.o\
     read_orb.o\
     setup_estate_pw.o\
     update_pot.o\
+    occ_matrix.o
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
       init_dm.o\
@@ -373,12 +374,23 @@ OBJS_HAMILT_OF=kedf_tf.o\
     kedf_manager.o\
     evolve_ofdft.o\
 
-OBJS_DFPT=dfpt_metal.o\
+OBJS_DFPT=dfpt_hamilt_shift.o\
+    dfpt_kq_basis.o\
+    dfpt_metal.o\
     dfpt_pert.o\
+    dfpt_pert_nl.o\
+    dfpt_pert_vkb.o\
     dfpt_phon.o\
+    dfpt_phon_elec.o\
+    dfpt_phon_ewald.o\
     dfpt_pw.o\
     dfpt_pw_data.o\
+    dfpt_pw_init.o\
+    dfpt_pw_q0.o\
+    dfpt_pw_run.o\
+    dfpt_pw_solve.o\
     dfpt_q0.o\
+    dfpt_q0_pos.o\
     dfpt_rho.o\
     dfpt_stern.o
 
@@ -399,13 +411,12 @@ OBJS_HAMILT_LCAO=hamilt_lcao.o\
     veff_lcao.o\
     veff_dh.o\
     meta_lcao.o\
-    dftu_lcao_op.o\
+    dftu_nao_op.o\
     deepks_lcao.o\
     op_exx_lcao.o\
     dspin_lcao.o\
     dspin_fs.o\
     setup_dftu_lcao.o\
-    dftu_fs.o\
     operator_fs_utils.o\
 
 OBJS_HCONTAINER=base_matrix.o\
@@ -564,6 +575,7 @@ OBJS_SYMMETRY=symm_other.o\
     symm_magnetic.o\
     symm_pricell.o\
     symm_rho.o\
+    little_group.o\
     symmetry.o\
 
 OBJS_XC=xc_functional.o\
@@ -674,6 +686,7 @@ OBJS_IO=module_parameter/input_conv.o\
     module_parameter/read_inp_model.o\
     module_parameter/read_inp_postproc.o\
     module_parameter/read_inp_exx_dftu.o\
+    module_parameter/read_inp_dfpt.o\
     module_parameter/read_inp_other.o\
     module_parameter/read_inp_out.o\
     module_parameter/read_set_globalv.o\
@@ -808,12 +821,12 @@ OBJS_SRCPW=h_ewald_pw.o\
     setup_pwwfc.o\
     update_cell_pw.o\
     dftu_base.o\
-    dftu_output.o\
-    dftu_tools_pw.o\
-    dftu_cal_occ_pw.o\
+    dftu_base_io.o\
+    dftu_base_occ.o\
+    dftu_base_tools.o\
+    yukawa_screening.o\
     setup_dftu_pw.o\
     deltaspin_pw.o\
-    deltaspin_pw_impl.o\
     force_pw.o\
     force_pw_us.o\
     force_pw_nl.o\
@@ -885,25 +898,29 @@ OBJS_VDW=vdw.o\
     vdwd3_autoset_xcname.o\
     vdwd3_auto_xcpar.o
 
-OBJS_DFTU=dftu_lcao.o\
-      dftu_force.o\
-      dftu_yukawa.o\
-      dftu_folding.o\
-      dftu_lcao_pots.o\
-      dftu_lcao_energy.o\
-      dftu_lcao_op_legacy.o\
-      dftu_lcao_occ.o\
+OBJS_DFTU=dftu_nao.o\
+      dftu_nao_fs_k.o\
+      dftu_nao_for_r.o\
+      dftu_nao_fs_r.o\
+      dftu_nao_str_r.o\
+      dftu_nao_folding.o\
+      dftu_nao_pots.o\
+      dftu_nao_energy.o\
+      dftu_nao_op_legacy.o\
+      dftu_nao_occ.o\
       dftu_hamilt.o
 
 OBJS_DELTASPIN=basic_funcs.o\
       cal_mw_from_lambda.o\
       cal_mw.o\
-      init_sc.o\
+      deltaspin_init.o\
       lambda_loop_helper.o\
       lambda_loop.o\
       spin_constrain.o\
-      cal_mw_helper.o\
       deltaspin_lcao.o\
+      deltaspin_lcao_mi.o\
+      deltaspin_state.o\
+      deltaspin_pw_mi.o\
       mi_tools.o\
       template_helpers.o\
 

source/source_base/module_container/ATen/core/tensor.cpp

@@ -295,7 +295,8 @@ bool Tensor::AllocateFrom(const Tensor& other, const TensorShape& shape) {
 
 void Tensor::sync(const Tensor& rhs) {
     REQUIRES_OK(this->data_type_ == rhs.data_type_ 
-        && this->device_ == rhs.device_)
+        && this->device_ == rhs.device_,
+        "sync: data_type and device must match between tensors")
 
     if (this->shape_ == rhs.shape_) {
         TEMPLATE_ALL_2(data_type_, device_,

source/source_base/module_container/ATen/kernels/memory.h

@@ -81,7 +81,7 @@ struct synchronize_memory_stride {
         const std::vector<int64_t>& out_size,
         const std::vector<int64_t>& in_size)
     {
-        REQUIRES_OK(in_size.size() == out_size.size() && in_size.size() <= 2);
+        REQUIRES_OK(in_size.size() == out_size.size() && in_size.size() <= 2, "rank mismatch: in_size and out_size must have the same rank <= 2");
         if (in_size.size() == 1) {
             synchronize_memory<T, Device_out, Device_in>()(arr_out, arr_in, in_size[0]);
         }

source/source_base/module_container/base/macros/macros.h

@@ -70,7 +70,7 @@
             __func__,                                   \
             __FILE__,                                   \
             static_cast<uint32_t>(__LINE__),            \
-            CHECK_MSG(expr, ##__VA_ARGS__));            \
+            CHECK_MSG(expr, __VA_ARGS__));             \
     }
 
 // The macro TEMPLATE_1() expands to a switch statement conditioned on

source/source_base/parallel_cell.cpp

@@ -2,13 +2,11 @@
 
 namespace ModuleBase
 {
-CommunicationDomain::CommunicationDomain()
-{
-}
-
 #ifdef __MPI
-CommunicationDomain::CommunicationDomain(MPI_Comm communicator) : communicator_(communicator)
+void CommunicationDomain::initialize(MPI_Comm communicator)
 {
+    communicator_ = communicator;
+    rank_ = 0;
     if (communicator_ != MPI_COMM_NULL)
     {
         MPI_Comm_rank(communicator_, &rank_);
@@ -28,10 +26,10 @@ int CommunicationDomain::rank() const
 
 CommunicationDomain world_communication_domain()
 {
+    CommunicationDomain communication_domain;
 #ifdef __MPI
-    return CommunicationDomain(MPI_COMM_WORLD);
-#else
-    return CommunicationDomain();
+    communication_domain.initialize(MPI_COMM_WORLD);
 #endif
+    return communication_domain;
 }
 } // namespace ModuleBase