A beta version of the DFPT program

docs/advanced/input_files/input-main.md

@@ -586,6 +586,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
 
@@ -634,6 +642,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
@@ -5170,3 +5179,55 @@
 - **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)
+- **Availability**: *esolver_type = dfpt*
+- **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
+- **Availability**: *esolver_type = dfpt*
+- **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.
+
+### dfpt_compute_q0
+
+- **Type**: Boolean
+- **Availability**: *esolver_type = dfpt*
+- **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
+- **Availability**: *esolver_type = dfpt*
+- **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
+- **Availability**: *esolver_type = dfpt*
+- **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
+- **Availability**: *esolver_type = dfpt*
+- **Description**: Set the maximum number of self-consistent DFPT iterations for each atomic displacement.
+- **Default**: 100
+
+### dfpt_mix_beta
+
+- **Type**: Real
+- **Availability**: *esolver_type = dfpt*
+- **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

@@ -53,6 +53,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: ""
@@ -2958,6 +2959,62 @@ parameters:
     default_value: "0"
     unit: ""
     availability: ""
+  - 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: esolver_type = dfpt
+  - 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: esolver_type = dfpt
+  - 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: esolver_type = dfpt
+  - 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: esolver_type = dfpt
+  - 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: esolver_type = dfpt
+  - 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: esolver_type = dfpt
+  - 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: esolver_type = dfpt
   - name: plot_istate
     category: Linear Response TDDFT
     type: Integer

source/source_cell/CMakeLists.txt

@@ -22,6 +22,7 @@ add_library(
     read_atoms_helper.cpp
     read_orb.cpp
     klist.cpp
+    reciprocal_grid.cpp
     parallel_kpoints.cpp
     cell_index.cpp
     cell_tools.cpp