feat(nttda): add analytic gradients for ROKS and AOCSCF

README.md

@@ -21,7 +21,7 @@ It is built on top of the PySCF electronic structure package.
         │                          │                          │
         ▼                          ▼                          ▼
 Noncollinear TDDFT         Noncollinear SF-TDDFT      Noncollinear Tensor TDA
-                                   ├── Gradients            ├── Gradients (in progress)
+                                   ├── Gradients            ├── Gradients
                                    ├── Oscillator strengths └── SOC (in progress)
                                    ├── NADC
                                    └── SOC
@@ -50,9 +50,18 @@ Noncollinear TDDFT         Noncollinear SF-TDDFT      Noncollinear Tensor TDA
 
   NT-TDA is a spin-consistent extension of noncollinear TDDFT that provides a unified treatment of spin-conserving and spin-flip excitations. For an open-shell reference state with total spin S, it can describe target states with total spins S−1 (except for S = 1/2), S, and S+1. All resulting states are free from spin contamination.
 
-  - Analytic gradients *(in progress)*
+  - Analytic gradients for ΔS = −1 and 0 on ROKS and average-occupation ROKS (AOCSCF) references
+  - Finite-difference gradients for all three spin channels
   - SOC *(in progress)*
 
+  `td.Gradients().kernel(state=n)` returns the gradient of the selected total
+  energy (`state=1` is the lowest NTTDA root; `state=0` is the reference).
+  `td.Gradients().as_scanner(state=n)` provides the energy and gradient for geometry optimization.
+  Analytic DFT gradients omit numerical grid response; finite-difference checks
+  use fixed grids by default. Density fitting, X2C, solvent response, and NLC
+  are not supported by the NTTDA gradient driver.
+  See the [AOCSCF gradient example](examples/nttda/02_nttda_aocscf_grad.py).
+
 
 ## Authors
 

examples/grad/02_aocscf_grad.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python
+# Copyright 2026 The NEST Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+'''
+Analytic nuclear gradient of the AOCSCF (average-occupation) reference.
+
+AOCSCF drives one set of orbitals with average occupations 2/1/0 for a
+high-spin open-shell reference, and takes the high-spin ROKS energy evaluated
+on those orbitals as the reference energy.  ``nuc_grad_method()`` returns the
+analytic gradient of that reference energy, which is the zero state used by
+the NTTDA excited-state gradients (see examples/nttda/02_nttda_aocscf_grad.py).
+
+The reference is non-stationary on the average-occupation orbitals, so the
+driver solves a Z-vector equation for the orbital response; a diffuse enough
+integration grid is required for the force sum to vanish.
+'''
+
+from pyscf import gto
+from nest import aocscf  # necessary import
+
+atom = '''
+N   0.000000  -0.040000   0.000000
+H   0.000000   0.780000   0.590000
+H   0.000000  -0.860000   0.520000
+'''
+mol = gto.M(atom=atom, charge=0, spin=1, basis='6-31g', verbose=3)
+fun = 'PBE'  # try also 'SVWN', 'B3LYP', 'M06-2X', etc.
+mf = mol.ROKS(xc=fun).average_occ()
+mf.conv_tol = 1e-12
+mf.conv_tol_grad = 1e-9
+mf.max_cycle = 120
+mf.grids.level = 5      # dense grid: the force sum is grid-sensitive
+mf.grids.prune = None
+mf.small_rho_cutoff = 0.0
+mf.kernel()
+
+print('AOCSCF reference energy: %.12f' % mf.e_tot)
+
+grad = mf.nuc_grad_method().kernel()
+print('Analytic reference gradient (Eh/Bohr):\n', grad)
+print('Force sum (should be ~0):\n', grad.sum(axis=0))
+
+# Gradients can also be restricted to selected atoms:
+grad_n = mf.nuc_grad_method().kernel(atmlst=[0])
+print('Gradient on the nitrogen atom only:\n', grad_n)

examples/nttda/02_nttda_aocscf_grad.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+# Copyright 2026 The NEST Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+'''
+NTTDA excited-state gradients on an AOCSCF (average-occupation) reference.
+
+An AOCSCF reference gives a common average-occupation orbital set for every
+spin channel.  NTTDA built on that reference can target the same-spin channel
+(``deltaS=0``) and the spin-lowering channel (``deltaS=-1``); the total state
+energy is the high-spin reference energy plus the NTTDA excitation energy, so
+its gradient is the sum of the reference gradient and the excitation-energy
+gradient.
+
+``td.Gradients().kernel(state=n)`` returns the analytic gradient of
+``E_reference + omega_n`` for state ``n`` (1 for the lowest root); ``state=0``
+returns the reference gradient.  Analytic gradients are available for
+``deltaS = -1`` and ``0``; ``deltaS = +1`` uses ``method='finite_diff'``.
+'''
+
+from pyscf import gto
+from nest import aocscf, nttda  # necessary imports
+
+atom = '''
+C    0.020000  -0.030000   0.010000
+H   -0.020000   0.800000   0.620000
+H    0.030000  -0.910000   0.500000
+'''
+mol = gto.M(atom=atom, charge=0, spin=2, basis='sto-3g', verbose=3)
+fun = 'B3LYP'
+mf = mol.ROKS(xc=fun).average_occ()
+mf.conv_tol = 1e-12
+mf.conv_tol_grad = 1e-9
+mf.max_cycle = 150
+mf.grids.level = 5      # dense grid: the force sum is grid-sensitive
+mf.grids.prune = None
+mf.small_rho_cutoff = 0.0
+mf.kernel()
+
+for delta_s in (-1, 0):
+    td = mf.NTTDA().set(
+        deltaS=delta_s,     # Sf = Si + deltaS
+        nstates=3,
+        conv_tol=1e-9,
+        max_cycle=200,
+        verbose=0,
+    ).run()
+    print('deltaS = %+d' % delta_s)
+    print('  NTTDA excitation energies:', td.e)
+    print('  total energies (E_ref + omega):', td.e_tot)
+
+    grad = td.Gradients().kernel(state=1)
+    print('  state-1 analytic gradient (Eh/Bohr):\n', grad)
+    print('  force sum (should be ~0):\n', grad.sum(axis=0))
+
+    ref_grad = td.Gradients().kernel(state=0)
+    print('  reference (state-0) gradient:\n', ref_grad)
+    print('  excitation-only contribution (state 1 - state 0):\n',
+          grad - ref_grad)