Feature: LR-TDDFT analytical forces (gradients) by Z-vector method and excited-state geometry relaxation

docs/advanced/input_files/input-main.md

@@ -574,6 +574,11 @@
     - [nocc](#nocc)
     - [nvirt](#nvirt)
     - [lr\_nstates](#lr_nstates)
+    - [lr\_target\_state](#lr_target_state)
+    - [lr\_degen\_thr](#lr_degen_thr)
+    - [lr\_degen\_mode](#lr_degen_mode)
+    - [lr\_grad\_solver](#lr_grad_solver)
+    - [lr\_target\_spin](#lr_target_spin)
     - [lr\_unrestricted](#lr_unrestricted)
     - [abs\_wavelen\_range](#abs_wavelen_range)
     - [out\_wfc\_lr](#out_wfc_lr)
@@ -5131,7 +5136,7 @@
 ### xc_kernel
 
 - **Type**: String
-- **Description**: The exchange-correlation kernel used in the calculation. Currently supported: RPA, LDA, PBE, HSE, HF.
+- **Description**: The exchange-correlation kernel used in the calculation. Currently supported: RPA, LDA, PWLDA, PBE, and the hybrids HF, PBE0, HSE, B3LYP, CAM_PBEH, LC_PBE, LC_WPBE, LRC_WPBE, LRC_WPBEH. A hybrid kernel needs the ground state to use the same functional: the exact-exchange operator $[\alpha+\beta\,\mathrm{erfc}(\mu r)]/r$ is built from exx_fock_alpha ($\alpha$), exx_erfc_alpha ($\beta$) and exx_erfc_omega ($\mu$), which are keyed off dft_functional, not off this parameter.
 - **Default**: LDA
 
 ### lr_init_xc_kernel
@@ -5161,9 +5166,10 @@
 ### nocc
 
 - **Type**: Integer
-- **Description**: The number of occupied orbitals (up to HOMO) used in the LR-TDDFT calculation.
-  - Note: If the value is illegal ( > nelec/2 or <= 0), it will be autoset to nelec/2.
-- **Default**: nband
+- **Description**: The number of occupied orbitals (up to HOMO) retained in the majority-spin LR-TDDFT window. A positive value selects a shared core prefix to discard from both spin channels; it does not change the ground-state occupations.
+  - If omitted, non-positive, or larger than the occupied majority-spin channel, all occupied orbitals are used.
+  - The full occupied window is determined by the effective electron number (including nelec_delta once) and the ground-state spin populations. For nspin=2, the minority-spin window has abs(N_up-N_down) fewer occupied orbitals.
+- **Default**: all occupied orbitals
 
 ### nvirt
 
@@ -5177,6 +5183,73 @@
 - **Description**: The number of 2-particle states to be solved.
 - **Default**: 0
 
+### lr_target_state
+
+- **Type**: Integer
+- **Description**: Index of the excited state whose potential energy surface `calculation = relax` follows, counted from 0 within the spin channel selected by `lr_target_spin`.
+
+  Only the gradient of this one state is computed, since solving the Z-vector equation dominates the cost of an excited-state gradient. It also selects the state whose excitation energy is added to the ground-state total energy, which is the quantity the energy-based relaxation algorithms (`cg`, `bfgs`, `lbfgs`) line-search on.
+
+  Ignored outside `calculation = relax`: a single-point run solves and reports the gradients of every state.
+
+  > Note: The state is followed by index, not by character. If it crosses another state during the relaxation, the optimizer will silently continue on the other surface.
+- **Default**: 0
+
+### lr_degen_thr
+
+- **Type**: Real
+- **Description**: Excited states whose excitation energies lie within this threshold of each other are treated as one degenerate multiplet, and the full gradient matrix $G^{(A\alpha)}_{kl}=\langle X_k|\partial A/\partial R_{A\alpha}|X_l\rangle$ is computed for it in addition to the per-state gradients. Zero (the default) disables this and leaves the per-state gradients as the only output.
+
+  At a $d$-fold degeneracy no single state has a gradient vector: the branch slopes along a displacement $u$ are the eigenvalues of $\sum_{A\alpha}u_{A\alpha}G^{(A\alpha)}$, and the eigenvectors that diagonalise it depend on $u$. The per-state gradients are the diagonal of $G$ in whichever basis the eigensolver happened to return, so only their sum (the trace) is basis-independent, while $G$ itself is the complete first-order information -- it is the linear vibronic coupling Hamiltonian of the multiplet. The extra cost is $d(d-1)/2$ further Z-vector solves per multiplet.
+
+  The threshold proposes candidates; it cannot tell a true degeneracy from an accidental near-degeneracy, where the states have genuinely different excitation energies and the construction does not apply. Each multiplet's actual energy spread and the orthonormality of its eigenvectors are reported in the running log so the distinction can be made there.
+
+  > Note: A sensible value is a few times the eigensolver threshold `lr_thr`, so that states split by real physics are not merged.
+- **Default**: 0
+- **Unit**: Ry
+
+### lr_degen_mode
+
+- **Type**: String
+- **Description**: What `calculation = relax` follows when `lr_target_state` sits inside a degenerate multiplet, as identified by `lr_degen_thr`. It has no effect when the target state is non-degenerate.
+
+  - state: follow the gradient of that one state, as returned by the eigensolver. This is the historical behaviour and is what reproduces earlier results, but inside a multiplet it is not a well-defined quantity: the per-state gradients are the diagonal of the subspace gradient matrix in whichever basis the eigensolver happened to return, so they depend on numerical details of the diagonalisation rather than on physics.
+  - average: follow the multiplet average $\bar\Omega=\frac{1}{d}\sum_k\Omega_k$, whose gradient is $\operatorname{Tr}G/d$. Unlike the individual states this is a smooth, basis-independent surface, and by symmetry its gradient is totally symmetric, so following it keeps the geometry on the symmetric configuration. Both the reported energy and the reported gradient switch to the average together, which the energy-based optimisers (`cg`, `bfgs`, `lbfgs`) require -- a gradient of one surface line-searched against the energy of another does not converge. This mode deliberately does NOT find the Jahn-Teller distortion, which is orthogonal to the totally symmetric average gradient.
+  - jt: descend the Jahn-Teller branch. Solves $\min_{\|u\|=1}\lambda_{\min}(\sum_{A\alpha}u_{A\alpha}G^{(A\alpha)})$ -- a joint optimisation over the displacement and the mixing inside the multiplet, since the two are determined together -- and follows the force of the resulting branch. This needs the off-diagonal part of the gradient matrix, so it costs $d(d-1)/2$ further Z-vector solves per step on top of the $d$ diagonal ones. The running log reports the branch's force, its mixing coefficients, and its split into the part common to the multiplet and the part that actually breaks the degeneracy.
+
+  > Note: The usual sequence is `average` first, to reach the symmetric stationary point, then `jt` from there: at a stationary point of the average surface the common part vanishes and the whole force is Jahn-Teller. `jt` is self-limiting -- once a step has split the multiplet there is no group left and the ordinary single-state gradient takes over.
+
+  > Note: `jt` gives the first-order DIRECTION. The distortion amplitude also needs the harmonic term, and the step norm is Cartesian rather than mass-weighted. A linear molecule has no first-order term at all (the effect is second-order Renner-Teller) and the log says so.
+- **Default**: state
+
+### lr_grad_solver
+
+- **Type**: String
+- **Description**: The method to solve the Z-vector (relaxed-density) equation $(A+B)Z=R$ in LR-TDDFT force and relaxation calculations, the linear-equation counterpart of `lr_solver`. Its dimension is $n_k n_{occ} n_{virt}$ summed over spin, where $n_{virt}$ counts every virtual band of the ground state, not only the `nvirt` window of the excitation.
+  - cg: Solve iteratively with the conjugate-gradient method, applying the orbital Hessian $A+B$ to a vector at each step. The matrix is never built.
+  - lapack: Construct the full matrix and solve directly with LAPACK (LU). Every MPI process holds the whole matrix and solves the same system.
+  - scalapack: Construct the matrix distributed over the MPI processes (2D block-cyclic) and solve with ScaLAPACK (LU).
+  - scalapack_chol: Construct the matrix distributed as for scalapack and solve by a ScaLAPACK Cholesky factorization, about half the flops of the LU and in place.
+  - elpa: Construct the matrix distributed as for scalapack and solve by an ELPA Cholesky factorization.
+
+  > Note: The direct solvers build the matrix column by column, at the cost of one application of $A+B$ per column, which usually dominates the cost of the solve itself; scalapack, scalapack_chol and elpa need an MPI build, elpa also an ELPA build.
+
+  > Note: scalapack_chol and elpa require $A+B$ to be positive definite, which holds at a stable ground state. If it is not, scalapack_chol stops with an error, while elpa does so only in a single-process run and hangs in a multi-process one; scalapack (LU) has no such requirement.
+- **Default**: cg
+
+### lr_target_spin
+
+- **Type**: String
+- **Description**: Which spin channel `lr_target_state` indexes.
+
+  - singlet / triplet: the two closed-shell channels solved at `nspin = 2`. At `nspin = 1` only `singlet` exists.
+  - updown: the single spin-conserving channel of an open-shell calculation (`lr_unrestricted`, or a spin-polarised ground state with a non-zero moment).
+
+  An open-shell calculation has only one channel, so any value is accepted there and relaxes that channel; an explicit `triplet` is reported as ignored. A closed-shell calculation rejects `updown`, since singlet and triplet are separate states with separate gradients.
+
+  Ignored outside `calculation = relax`.
+- **Default**: singlet
+
 ### lr_unrestricted
 
 - **Type**: Boolean

docs/parameters.yaml

@@ -2940,7 +2940,7 @@ parameters:
     category: Linear Response TDDFT
     type: String
     description: |
-      The exchange-correlation kernel used in the calculation. Currently supported: RPA, LDA, PBE, HSE, HF.
+      The exchange-correlation kernel used in the calculation. Currently supported: RPA, LDA, PWLDA, PBE, and the hybrids HF, PBE0, HSE, B3LYP, CAM_PBEH, LC_PBE, LC_WPBE, LRC_WPBE, LRC_WPBEH. A hybrid kernel needs the ground state to use the same functional: the exact-exchange operator $[\alpha+\beta\,\mathrm{erfc}(\mu r)]/r$ is built from exx_fock_alpha ($\alpha$), exx_erfc_alpha ($\beta$) and exx_erfc_omega ($\mu$), which are keyed off dft_functional, not off this parameter.
     default_value: LDA
     unit: ""
     availability: ""
@@ -2978,9 +2978,10 @@ parameters:
     category: Linear Response TDDFT
     type: Integer
     description: |
-      The number of occupied orbitals (up to HOMO) used in the LR-TDDFT calculation.
-      * Note: If the value is illegal ( > nelec/2 or <= 0), it will be autoset to nelec/2.
-    default_value: nband
+      The number of occupied orbitals (up to HOMO) retained in the majority-spin LR-TDDFT window. A positive value selects a shared core prefix to discard from both spin channels; it does not change the ground-state occupations.
+      * If omitted, non-positive, or larger than the occupied majority-spin channel, all occupied orbitals are used.
+      * The full occupied window is determined by the effective electron number (including nelec_delta once) and the ground-state spin populations. For nspin=2, the minority-spin window has abs(N_up-N_down) fewer occupied orbitals.
+    default_value: all occupied orbitals
     unit: ""
     availability: ""
   - name: nvirt
@@ -2999,6 +3000,82 @@ parameters:
     default_value: "0"
     unit: ""
     availability: ""
+  - name: lr_target_state
+    category: Linear Response TDDFT
+    type: Integer
+    description: |
+      Index of the excited state whose potential energy surface `calculation = relax` follows, counted from 0 within the spin channel selected by `lr_target_spin`.
+
+      Only the gradient of this one state is computed, since solving the Z-vector equation dominates the cost of an excited-state gradient. It also selects the state whose excitation energy is added to the ground-state total energy, which is the quantity the energy-based relaxation algorithms (`cg`, `bfgs`, `lbfgs`) line-search on.
+
+      Ignored outside `calculation = relax`: a single-point run solves and reports the gradients of every state.
+
+      [NOTE] The state is followed by index, not by character. If it crosses another state during the relaxation, the optimizer will silently continue on the other surface.
+    default_value: "0"
+    unit: ""
+    availability: ""
+  - name: lr_degen_thr
+    category: Linear Response TDDFT
+    type: Real
+    description: |
+      Excited states whose excitation energies lie within this threshold of each other are treated as one degenerate multiplet, and the full gradient matrix $G^{(A\alpha)}_{kl}=\langle X_k|\partial A/\partial R_{A\alpha}|X_l\rangle$ is computed for it in addition to the per-state gradients. Zero (the default) disables this and leaves the per-state gradients as the only output.
+
+      At a $d$-fold degeneracy no single state has a gradient vector: the branch slopes along a displacement $u$ are the eigenvalues of $\sum_{A\alpha}u_{A\alpha}G^{(A\alpha)}$, and the eigenvectors that diagonalise it depend on $u$. The per-state gradients are the diagonal of $G$ in whichever basis the eigensolver happened to return, so only their sum (the trace) is basis-independent, while $G$ itself is the complete first-order information -- it is the linear vibronic coupling Hamiltonian of the multiplet. The extra cost is $d(d-1)/2$ further Z-vector solves per multiplet.
+
+      The threshold proposes candidates; it cannot tell a true degeneracy from an accidental near-degeneracy, where the states have genuinely different excitation energies and the construction does not apply. Each multiplet's actual energy spread and the orthonormality of its eigenvectors are reported in the running log so the distinction can be made there.
+
+      [NOTE] A sensible value is a few times the eigensolver threshold `lr_thr`, so that states split by real physics are not merged.
+    default_value: "0"
+    unit: Ry
+    availability: ""
+  - name: lr_degen_mode
+    category: Linear Response TDDFT
+    type: String
+    description: |
+      What `calculation = relax` follows when `lr_target_state` sits inside a degenerate multiplet, as identified by `lr_degen_thr`. It has no effect when the target state is non-degenerate.
+
+      * state: follow the gradient of that one state, as returned by the eigensolver. This is the historical behaviour and is what reproduces earlier results, but inside a multiplet it is not a well-defined quantity: the per-state gradients are the diagonal of the subspace gradient matrix in whichever basis the eigensolver happened to return, so they depend on numerical details of the diagonalisation rather than on physics.
+      * average: follow the multiplet average $\bar\Omega=\frac{1}{d}\sum_k\Omega_k$, whose gradient is $\operatorname{Tr}G/d$. Unlike the individual states this is a smooth, basis-independent surface, and by symmetry its gradient is totally symmetric, so following it keeps the geometry on the symmetric configuration. Both the reported energy and the reported gradient switch to the average together, which the energy-based optimisers (`cg`, `bfgs`, `lbfgs`) require -- a gradient of one surface line-searched against the energy of another does not converge. This mode deliberately does NOT find the Jahn-Teller distortion, which is orthogonal to the totally symmetric average gradient.
+      * jt: descend the Jahn-Teller branch. Solves $\min_{\|u\|=1}\lambda_{\min}(\sum_{A\alpha}u_{A\alpha}G^{(A\alpha)})$ -- a joint optimisation over the displacement and the mixing inside the multiplet, since the two are determined together -- and follows the force of the resulting branch. This needs the off-diagonal part of the gradient matrix, so it costs $d(d-1)/2$ further Z-vector solves per step on top of the $d$ diagonal ones. The running log reports the branch's force, its mixing coefficients, and its split into the part common to the multiplet and the part that actually breaks the degeneracy.
+
+      [NOTE] The usual sequence is `average` first, to reach the symmetric stationary point, then `jt` from there: at a stationary point of the average surface the common part vanishes and the whole force is Jahn-Teller. `jt` is self-limiting -- once a step has split the multiplet there is no group left and the ordinary single-state gradient takes over.
+
+      [NOTE] `jt` gives the first-order DIRECTION. The distortion amplitude also needs the harmonic term, and the step norm is Cartesian rather than mass-weighted. A linear molecule has no first-order term at all (the effect is second-order Renner-Teller) and the log says so.
+    default_value: state
+    unit: ""
+    availability: ""
+  - name: lr_grad_solver
+    category: Linear Response TDDFT
+    type: String
+    description: |
+      The method to solve the Z-vector (relaxed-density) equation $(A+B)Z=R$ in LR-TDDFT force and relaxation calculations, the linear-equation counterpart of `lr_solver`. Its dimension is $n_k n_{occ} n_{virt}$ summed over spin, where $n_{virt}$ counts every virtual band of the ground state, not only the `nvirt` window of the excitation.
+      * cg: Solve iteratively with the conjugate-gradient method, applying the orbital Hessian $A+B$ to a vector at each step. The matrix is never built.
+      * lapack: Construct the full matrix and solve directly with LAPACK (LU). Every MPI process holds the whole matrix and solves the same system.
+      * scalapack: Construct the matrix distributed over the MPI processes (2D block-cyclic) and solve with ScaLAPACK (LU).
+      * scalapack_chol: Construct the matrix distributed as for scalapack and solve by a ScaLAPACK Cholesky factorization, about half the flops of the LU and in place.
+      * elpa: Construct the matrix distributed as for scalapack and solve by an ELPA Cholesky factorization.
+
+      [NOTE] The direct solvers build the matrix column by column, at the cost of one application of $A+B$ per column, which usually dominates the cost of the solve itself; scalapack, scalapack_chol and elpa need an MPI build, elpa also an ELPA build.
+
+      [NOTE] scalapack_chol and elpa require $A+B$ to be positive definite, which holds at a stable ground state. If it is not, scalapack_chol stops with an error, while elpa does so only in a single-process run and hangs in a multi-process one; scalapack (LU) has no such requirement.
+    default_value: cg
+    unit: ""
+    availability: ""
+  - name: lr_target_spin
+    category: Linear Response TDDFT
+    type: String
+    description: |
+      Which spin channel `lr_target_state` indexes.
+
+      * singlet / triplet: the two closed-shell channels solved at `nspin = 2`. At `nspin = 1` only `singlet` exists.
+      * updown: the single spin-conserving channel of an open-shell calculation (`lr_unrestricted`, or a spin-polarised ground state with a non-zero moment).
+
+      An open-shell calculation has only one channel, so any value is accepted there and relaxes that channel; an explicit `triplet` is reported as ignored. A closed-shell calculation rejects `updown`, since singlet and triplet are separate states with separate gradients.
+
+      Ignored outside `calculation = relax`.
+    default_value: singlet
+    unit: ""
+    availability: ""
   - name: lr_unrestricted
     category: Linear Response TDDFT
     type: Boolean

source/Makefile.Objects

@@ -129,6 +129,7 @@ ${OBJS_DELTASPIN}\
 ${OBJS_TENSOR}\
 ${OBJS_HSOLVER_PEXSI}\
 ${OBJS_LR}\
+${OBJS_LR_GRAD}\
 ${OBJS_RDMFT}
 
 OBJS_MAIN=main.o\
@@ -1032,7 +1033,6 @@ OBJS_TENSOR=tensor.o\
     refcount.o
 
     OBJS_LR=lr_util.o\
-    lr_util_hcontainer.o\
     utils/lr_io.o\
     utils/exciton_plotter.o\
     ao_to_mo_parallel.o\
@@ -1048,6 +1048,28 @@ OBJS_TENSOR=tensor.o\
     hamilt_casida.o\
     esolver_lr_lcao_tddft.o\
 
+ifdef LIBRI_DIR
+# BSE-related code and DMBand: only compiled with LibRI (__EXX), see module_lr/CMakeLists.txt
+OBJS_LR+=utils/lr_io_krlist.o
+OBJS_LR+=dm_band.o
+endif
+
+OBJS_LR_GRAD=lr_force.o\
+    lr_force_test.o\
+    grad_degen.o\
+    grad_jt.o\
+    gradient_output.o\
+    lr_grad_cs.o\
+    lr_grad_os.o\
+    exx_proj.o\
+    cvcx_serial.o\
+    cvcx_par.o\
+    cal_edm.o\
+    zeqlin_solv.o\
+    pot_grad_xc.o\
+    esolver_lr_grad.o\
+    esolver_lr_rlx.o\
+
 OBJS_RDMFT=rdmft.o\
         rdmft_tools.o\
         rdmft_pot.o\