feat: add deltaSCF methods

examples/deltascf/01_sgm.py

@@ -0,0 +1,90 @@
+#!/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.
+
+"""CH4 core excitation with ordinary SGM and a canonical-orbital guess.
+
+This is a restricted triplet calculation.
+Only use the excitation energy when the excited-state calculation converges.
+"""
+# Q-Chem reference input for this CH4 core excitation benchmark:
+#
+# $rem
+# METHOD BHHLYP
+# BASIS aug-cc-pvtz
+# SYMMETRY false
+# SYM_IGNORE true
+# NO_REORIENT True
+# XC_GRID 000075000302
+# UNRESTRICTED false
+# delta_scf true
+# $end
+#
+# $delta_scf
+#   triplet restricted
+#   triplet_SCF_algorithm SGM
+#   somo_1 1
+#   somo_2 6
+# $end
+#
+# Q-Chem output:
+#   Restricted open-shell triplet state = -29.947480 Ha
+#   SOMO(1): best initial orbital 5 overlap 0.999002
+#   SOMO(2): best initial orbital 6 overlap 0.879378
+#   Excitation energy = 10.554376 Ha / 287.199215 eV
+
+from pyscf import gto
+from pyscf.data import nist
+from nest.deltascf import sgm
+
+# 1. Calculate the closed-shell ground state.
+mol = gto.M(
+    atom="""
+    C  0.0000  0.0000  0.0000
+    H  0.6276  0.6276  0.6276
+    H  0.6276 -0.6276 -0.6276
+    H -0.6276  0.6276 -0.6276
+    H -0.6276 -0.6276  0.6276
+    """,
+    basis='aug-cc-pvtz',
+    spin=0,
+    symmetry=False,
+)
+ground = mol.RKS(xc='BHandHLYP')
+ground.grids.atom_grid = (75, 302)
+ground.kernel()
+
+# 2. Make a triplet guess by changing two occupations from 2/0 to 1/1.
+# Indices start at zero: orbital 0 is C 1s, orbital 5 is the ground-state LUMO.
+hole_index = 0
+particle_index = 5
+initial_orbitals = ground.mo_coeff.copy()
+initial_occupations = ground.mo_occ.copy()
+initial_occupations[hole_index] = 1
+initial_occupations[particle_index] = 1
+
+# 3. Optimize the excited-state orbitals using the original SGM preconditioner.
+triplet_mol = mol.copy()
+triplet_mol.spin = 2  # PySCF spin = N_alpha - N_beta = 2S.
+triplet_scf = triplet_mol.ROKS(xc=ground.xc)
+triplet_scf.grids.atom_grid = (75, 302)
+excited = sgm.SGM(triplet_scf)
+excited.verbose = 4
+excited.kernel(initial_orbitals, initial_occupations)
+
+print('SGM converged:', excited.converged)
+print('Last energy (Ha):', excited.e_tot)
+if excited.converged:
+    excitation_energy = (excited.e_tot - ground.e_tot) * nist.HARTREE2EV
+    print('Triplet excitation energy (eV):', excitation_energy)

examples/deltascf/02_fr.py

@@ -0,0 +1,66 @@
+#!/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.
+
+"""CH4 core excitation with ordinary FR and a canonical-orbital guess.
+
+This is a restricted triplet calculation.
+Only use the excitation energy when the excited-state calculation converges.
+"""
+
+from pyscf import gto
+from pyscf.data import nist
+from nest.deltascf import fr
+
+# 1. Calculate the closed-shell ground state.
+mol = gto.M(
+    atom="""
+    C  0.0000  0.0000  0.0000
+    H  0.6276  0.6276  0.6276
+    H  0.6276 -0.6276 -0.6276
+    H -0.6276  0.6276 -0.6276
+    H -0.6276 -0.6276  0.6276
+    """,
+    basis='aug-cc-pvtz',
+    spin=0,
+    symmetry=False,
+)
+ground = mol.RKS(xc='BHandHLYP')
+ground.grids.atom_grid = (75, 302)
+ground.kernel()
+
+# 2. Make a triplet guess by changing two occupations from 2/0 to 1/1.
+# Indices start at zero: orbital 0 is C 1s, orbital 5 is the ground-state LUMO.
+hole_index = 0
+particle_index = 5
+initial_orbitals = ground.mo_coeff.copy()
+initial_occupations = ground.mo_occ.copy()
+initial_occupations[hole_index] = 1
+initial_occupations[particle_index] = 1
+
+# 3. Freeze both SOMOs, relax the other orbitals, then release with IMOM.
+triplet_mol = mol.copy()
+triplet_mol.spin = 2  # PySCF spin = N_alpha - N_beta = 2S.
+triplet_scf = triplet_mol.ROKS(xc=ground.xc)
+triplet_scf.grids.atom_grid = (75, 302)
+excited = fr.FR(triplet_scf)
+excited.verbose = 4
+excited.kernel(initial_orbitals, initial_occupations)
+
+print('Freeze converged:', excited.freeze_converged, 'cycles:', excited.freeze_cycles)
+print('Release converged:', excited.converged, 'cycles:', excited.cycles)
+print('Last energy (Ha):', excited.e_tot)
+if excited.converged:
+    excitation_energy = (excited.e_tot - ground.e_tot) * nist.HARTREE2EV
+    print('Triplet excitation energy (eV):', excitation_energy)

examples/deltascf/03_tda_guess.py

@@ -0,0 +1,109 @@
+#!/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.
+
+"""C2H3F C1 core excitation: a core-TDA guess followed by SGM or FR.
+
+The chosen TDA root defines an initial guess; it does not guarantee the
+identity of the final SCF solution. Check convergence before using energies.
+"""
+
+import numpy as np
+from scipy.linalg import null_space
+from pyscf import gto, tdscf
+from pyscf.data import nist
+from nest.deltascf import fr, sgm
+
+# 1. Calculate the closed-shell ground state.
+mol = gto.M(
+    atom="""
+    C1  0.000000  -0.246412  -1.271068
+    C2  0.000000   0.457081  -0.154735
+    F1  0.000000  -0.119195   1.052878
+    H1  0.000000   0.272328  -2.210194
+    H2  0.000000  -1.319906  -1.249847
+    H3  0.000000   1.530323  -0.095954
+    """,
+    basis={'C1': 'aug-cc-pvtz', 'C2': 'aug-cc-pvtz',
+           'F1': 'cc-pvtz', 'H': 'cc-pvdz'},
+    spin=0,
+    symmetry=False,
+)
+ground = mol.RKS(xc='BHandHLYP')
+ground.kernel()
+
+# 2. Obtain a particle orbital from single-hole, triplet core-TDA.
+# These zero-based indices apply to THIS molecule and basis.
+hole_index = 2      # C1 1s orbital.
+particle_index = 12 # First ground-state virtual; also the slot for our particle orbital.
+tda_root = 1       # Second TDA root: the even-parity root with a large LUMO component.
+
+# Allow excitations only from the chosen core orbital, to all virtual orbitals.
+# This freezing applies to TDA only, not to the subsequent SCF calculation.
+frozen_occupied = []
+for orbital_index in range(particle_index):
+    if orbital_index != hole_index:
+        frozen_occupied.append(orbital_index)
+
+tda = tdscf.TDA(ground)
+tda.singlet = False
+tda.frozen = frozen_occupied
+tda.kernel(nstates=3)
+
+for root_index, (energy, amplitudes) in enumerate(zip(tda.e, tda.xy)):
+    excitation_amplitudes, _ = amplitudes
+    particle_vector = excitation_amplitudes[0].copy()
+    particle_vector /= np.linalg.norm(particle_vector)
+    lumo_weight = particle_vector[0]**2
+    print(f'TDA root {root_index}: {energy * nist.HARTREE2EV:.6f} eV, '
+          f'ground-state LUMO weight = {lumo_weight:.4f}')
+
+# Only one occupied orbital participates, so X has a single row. Its entries
+# are the coefficients of the particle orbital in the OLD virtual basis.
+# Normalize explicitly because PySCF's restricted TDA uses 2 * ||X||^2 = 1.
+excitation_amplitudes, _ = tda.xy[tda_root]
+particle_vector = excitation_amplitudes[0].copy()
+particle_vector /= np.linalg.norm(particle_vector)
+
+# Put the particle orbital in the original LUMO slot. Complete the remaining
+# virtual columns with orthonormal combinations perpendicular to that vector.
+# These other columns are a basis completion, not higher TDA excited states.
+old_virtual_orbitals = ground.mo_coeff[:, particle_index:]
+remaining_virtual_vectors = null_space(particle_vector.reshape(1, -1))
+initial_orbitals = ground.mo_coeff.copy()
+initial_orbitals[:, particle_index] = old_virtual_orbitals @ particle_vector
+initial_orbitals[:, particle_index + 1:] = old_virtual_orbitals @ remaining_virtual_vectors
+
+initial_occupations = ground.mo_occ.copy()
+initial_occupations[hole_index] = 1
+initial_occupations[particle_index] = 1
+
+# 3. Optimize the restricted triplet. Choose FR or SGM on the next line.
+triplet_mol = mol.copy()
+triplet_mol.spin = 2
+triplet_scf = triplet_mol.ROKS(xc=ground.xc)
+excited = fr.FR(triplet_scf)  # Or: sgm.SGM(triplet_scf)
+excited.verbose = 4
+excited.conv_tol_grad = 1e-5
+# FR first freezes both singly occupied orbitals, then releases them with IMOM.
+# verbose=4 prints both stage boundaries and the final SOMO-subspace overlaps.
+excited.conv_tol = 1e-9
+excited.max_cycle = 200
+excited.kernel(initial_orbitals, initial_occupations)
+
+print('Optimizer converged:', excited.converged, 'cycles:', excited.cycles)
+print('Last energy (Ha):', excited.e_tot)
+if excited.converged:
+    excitation_energy = (excited.e_tot - ground.e_tot) * nist.HARTREE2EV
+    print('Triplet excitation energy (eV):', excitation_energy)

examples/deltascf/04_nttda.py

@@ -0,0 +1,78 @@
+#!/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.
+
+"""Run deltaS=-1 NTTDA with SGM/FR delta-SCF references for CH4.
+
+Optimize the core-excited triplet with SGM or FR, then use it as the
+reference for singlet NTTDA. Use the CH4 geometry and C 1s/LUMO guess
+from 01_sgm.py and 02_fr.py, with a 6-31G* basis and the default grid.
+"""
+
+import numpy as np
+from pyscf import gto
+from pyscf.data import nist
+from nest import nttda  # Registers NTTDA on ROKS.
+from nest.deltascf import fr, sgm
+
+# 1. Calculate the closed-shell ground state.
+mol = gto.M(
+    atom="""
+    C  0.0000  0.0000  0.0000
+    H  0.6276  0.6276  0.6276
+    H  0.6276 -0.6276 -0.6276
+    H -0.6276  0.6276 -0.6276
+    H -0.6276 -0.6276  0.6276
+    """,
+    basis='631G*',
+    spin=0,
+    symmetry=False,
+)
+ground = mol.RKS(xc='BHandHLYP')
+ground.kernel()
+
+# 2. Make a triplet guess by changing two occupations from 2/0 to 1/1.
+# Indices start at zero: orbital 0 is C 1s, orbital 5 is the ground-state LUMO.
+hole_index = 0
+particle_index = 5
+initial_orbitals = ground.mo_coeff.copy()
+initial_occupations = ground.mo_occ.copy()
+initial_occupations[hole_index] = 1
+initial_occupations[particle_index] = 1
+
+# 3. Optimize the delta-SCF triplet reference with SGM or FR.
+triplet_mol = mol.copy()
+triplet_mol.spin = 2
+
+triplet_scf = triplet_mol.ROKS(xc=ground.xc)
+excited = fr.FR(triplet_scf)  # or sgm.SGM(triplet_scf)
+excited.verbose = 4
+excited.kernel(initial_orbitals, initial_occupations)
+excited.analyze(verbose=4)
+
+# 4. Run deltaS=-1 NTTDA on the optimized delta-SCF reference.
+reference = excited  # for SGM, reference = excited.undo_sgm()
+td = reference.NTTDA()
+td.deltaS = -1  # S_reference=1 -> S_final=0 (singlets).
+td.nstates = 10  # Include core-excited singlets above the lower valence roots.
+td.kernel()
+td.analyze(verbose=4)
+
+# Supplementary excitation-energy comparison.
+# E_triplet - E_RKS defines the delta-SCF gap for the settings used here.
+# omega is relative to the triplet, so add that gap for a common RKS zero.
+delta_scf_ev = (reference.e_tot - ground.e_tot) * nist.HARTREE2EV
+nttda_rks_ev = delta_scf_ev + td.e * nist.HARTREE2EV
+nttda_excitation_ev = (td.e - td.e[0]) * nist.HARTREE2EV
+print(f'delta-SCF triplet excitation = {delta_scf_ev:.6f} eV')