Skip to content

Commit 2b96a94

Browse files
yideng2025yideng2025
authored
Add generalized active space CI solver with OpenMP C kernels (#203)
* fci: add core GAS FCI C kernels * build: add CMake target for GAS FCI kernels * feat(mcscf): add GASCI solver * fix(mcscf): finalize GASCI kernel and Python interfaces Honor FCISolver.spin when normalizing active electron counts and remove GASSCF-only library compatibility paths. Add PySCF source headers and finalize comments in the frozen GAS C kernels. * test(mcscf): add GASCI tests and examples Cover restriction normalization, GAS-as-CAS behavior, restricted spaces, RDMs, spin guards, multiroot GASCI, and state averaging. Add compact examples for explicit spin-supergroups, restriction formats, and cumulative-occ multiroot analysis. * fix(mcscf): propagate solver spin through GASCI Use FCISolver.spin consistently when normalizing GAS restrictions, reporting determinant spaces, tracking CI reuse signatures, and evaluating spin-square values. Add a GASCI regression test for spin-derived active electron counts. * fix(mcscf): address GASCI review feedback Move the GAS FCI target into mcscf/CMakeLists.txt while preserving the isolated kernel build. Make binomial initialization thread-safe, select tied diagonal addresses deterministically, and add a regression test. * fix(mcscf): label GASCI kernel output Use a GASCI-specific fixed-orbital kernel wrapper so runtime messages and timers identify the GAS calculation rather than the reused CASCI orchestration. Add a regression test requiring Start GASCI and excluding Start CASCI. * examples: validate GASCI energies against OpenMolcas * mcscf: add O2 GASCI energy regression test * mcscf: add OpenMolcas GASCI regression coverage * mcscf: add GAS-FCI projection validation * mcscf: unify GAS Hamiltonian contraction plans * mcscf: move spin-free GAS RDM assembly to Python * mcscf: manage GAS contraction intermediates in Python * mcscf: reuse PySCF electron-count unpacking * mcscf: require integer GAS restriction inputs * mcscf: simplify GAS scalar validation * mcscf: inline GAS CAS block construction * mcscf: remove unused GAS helper * mcscf: align GAS diagonal selection with PySCF * mcscf: reuse GAS space for CI initial guesses * mcscf: reuse CASCI one-electron Hamiltonian helper * mcscf: add transition metal RASCI example * fix(mcscf): update Al2FeO4+ RASCI orbitals * refactor(mcscf): make GAS binomial table static * style(mcscf): clean up GAS comments and spacing
1 parent 7a1bdab commit 2b96a94

24 files changed

Lines changed: 12907 additions & 1 deletion

‎examples/mcscf/90_gasci_spin_supergroup.py‎

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
#!/usr/bin/env python
2+
#
3+
# Author: Yi Deng <yideng@uchicago.edu>
4+
#
5+
6+
"""Fixed-orbital GASCI with an explicit spin-supergroup set D."""
7+
8+
import numpy
9+
10+
from pyscf import gto
11+
from pyscf import scf
12+
from pyscf.mcscf import gasci
13+
14+
15+
mol = gto.M(
16+
atom="N 0 0 0; N 0 0 1.10",
17+
basis="sto-3g",
18+
spin=0,
19+
verbose=0,
20+
)
21+
mf = scf.RHF(mol).run()
22+
23+
# Each row is (alpha occupations, beta occupations). Rows are normalized to
24+
# lexicographic order and duplicates are removed before entering the C kernel.
25+
gas_restr = numpy.asarray([
26+
[0, 2, 2, 0],
27+
[1, 1, 1, 1],
28+
[2, 0, 0, 2],
29+
], dtype=numpy.int32)
30+
31+
mc = gasci.GASCI(
32+
mf, 4, (2, 2), ncore=5,
33+
gas_orbs=(2, 2),
34+
gas_restr=gas_restr,
35+
gas_restr_type="spin-supergroup",
36+
)
37+
mc.verbose = 4
38+
e_tot, e_gas, ci, _, _ = mc.kernel()
39+
40+
print("GASCI total energy = %.12f" % e_tot)
41+
print("GASCI active-space energy = %.12f" % e_gas)
42+
print("GAS determinant count = %d" % ci.size)
43+
print("GAS space information = %s" % mc.gas_space_info())

‎examples/mcscf/91_gas_restriction_types.py‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
#!/usr/bin/env python
2+
#
3+
# Author: Yi Deng <yideng@uchicago.edu>
4+
#
5+
6+
"""Normalize the four supported GAS restriction input formats."""
7+
8+
import numpy
9+
10+
from pyscf.mcscf import addons_gas
11+
12+
13+
gas_orbs = (2, 2)
14+
nelec = (2, 2)
15+
16+
spin_supergroups = numpy.asarray([
17+
[2, 0, 0, 2],
18+
[1, 1, 1, 1],
19+
[0, 2, 2, 0],
20+
[1, 1, 1, 1], # Duplicate rows are removed.
21+
])
22+
supergroups = numpy.asarray([[2, 2]])
23+
cumulative_occ = numpy.asarray([[2, 2], [4, 4]])
24+
25+
for restriction_type, restriction in (
26+
("spin-supergroup", spin_supergroups),
27+
("supergroup", supergroups),
28+
("cumulative-occ", cumulative_occ)):
29+
blocks = addons_gas.normalize_gas_restr(
30+
gas_orbs, nelec, restriction, restriction_type)
31+
print("%-17s ->\n%s" % (restriction_type, blocks))
32+
33+
ras_blocks = addons_gas.normalize_gas_restr(
34+
(1, 2, 1), nelec,
35+
{"max_holes": 1, "max_particles": 1},
36+
"ras",
37+
)
38+
print("ras ->\n%s" % ras_blocks)

‎examples/mcscf/92_gasci_multiroot_rdm.py‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
#!/usr/bin/env python
2+
#
3+
# Author: Yi Deng <yideng@uchicago.edu>
4+
#
5+
6+
"""Multiroot GASCI with cumulative restrictions and state-resolved RDMs."""
7+
8+
import numpy
9+
10+
from pyscf import gto
11+
from pyscf import scf
12+
from pyscf.mcscf import gasci
13+
14+
15+
mol = gto.M(
16+
atom="N 0 0 0; N 0 0 1.10", basis="sto-3g", spin=0, verbose=0)
17+
mf = scf.RHF(mol).run()
18+
19+
# Each row gives the lower and upper cumulative electron occupation after
20+
# adding the corresponding GAS subspace.
21+
mc = gasci.GASCI(
22+
mf, 4, (2, 2), ncore=5,
23+
gas_orbs=(2, 2),
24+
gas_restr=[[1, 3], [4, 4]],
25+
gas_restr_type="cumulative-occ",
26+
)
27+
mc.fcisolver.nroots = 3
28+
mc.verbose = 4
29+
energies, _, _, _, _ = mc.kernel()
30+
31+
for state, energy in enumerate(numpy.asarray(energies).reshape(-1)):
32+
dm1, dm2 = mc.make_gasdm12(state=state)
33+
ss, multiplicity = mc.spin_square(state=state)
34+
pair_trace = numpy.einsum("pprr", dm2).real
35+
print("state %d E = %.12f Tr(D1) = %.8f S^2 = %.8f mult = %.4f" %
36+
(state, energy, numpy.trace(dm1), ss, multiplicity))
37+
print("state %d Tr_pair(D2) = %.8f" % (state, pair_trace))
38+
39+
tdm1, tdm2 = mc.trans_gasdm12(bra_state=0, ket_state=2)
40+
print("transition 0 -> 2: ||D1|| = %.8f ||D2|| = %.8f" %
41+
(numpy.linalg.norm(tdm1), numpy.linalg.norm(tdm2)))
42+
43+
_, natural_occ = mc.get_gas_natorb(state=0)
44+
_, pseudo_occ = mc.get_gas_pseudo_natorb(state=0)
45+
print("natural occupations = %s" % natural_occ)
46+
for igas, occupations in enumerate(pseudo_occ):
47+
print("GAS%d pseudo-natural occupations = %s" %
48+
(igas + 1, occupations))

‎examples/mcscf/gasci_openmolcas/01_o2_multiroot.py‎

Lines changed: 252 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,252 @@
1+
#!/usr/bin/env python
2+
#
3+
# Author: Yi Deng <yideng@uchicago.edu>
4+
#
5+
6+
"""Fixed-orbital multiroot GASCI example using OpenMolcas orbitals."""
7+
8+
from pathlib import Path
9+
10+
import numpy
11+
12+
from pyscf import gto
13+
from pyscf import scf
14+
from pyscf.mcscf import gasci
15+
16+
17+
# OpenMolcas reference input:
18+
#
19+
# &GATEWAY
20+
# Title = O2 GASCI 2S+1=3
21+
# Coord = o2.xyz
22+
# Basis = def2-TZVP
23+
# Group = C1
24+
# noCD
25+
#
26+
# &SEWARD
27+
#
28+
# &RASSCF
29+
# Title = O2 GASCI(12e,8o), GAS=2+4+2, triplet
30+
#
31+
# CIONLY
32+
# FILEORB = o2_guessorb.INPORB
33+
#
34+
# Spin = 3
35+
# Symmetry = 1
36+
#
37+
# Nactel = 12 0 0
38+
# Inactive = 2
39+
#
40+
# GASSCF
41+
# 3
42+
# 2
43+
# 2 4
44+
# 4
45+
# 8 10
46+
# 2
47+
# 12 12
48+
#
49+
# CIRoot = 20 20 1
50+
# CIMX = 200
51+
#
52+
# Notes:
53+
# - OpenMolcas Spin and PySCF spin use multiplicity 2S+1 and 2S,
54+
# respectively.
55+
# - CIONLY keeps the molecular orbitals fixed.
56+
# - noCD avoids the Cholesky approximation to the ERIs.
57+
#
58+
# The exact PySCF molecule and fixed MO coefficients were imported from
59+
# the converged OpenMolcas HDF5 file as follows:
60+
#
61+
# import numpy
62+
# from mrh.my_pyscf.tools.molcas2pyscf import get_mo_from_h5
63+
# from mrh.my_pyscf.tools.molcas2pyscf import get_mol_from_h5
64+
#
65+
# h5_file = "o2_m3.rasscf.h5"
66+
# mol = get_mol_from_h5(
67+
# h5_file,
68+
# charge=0,
69+
# spin=2,
70+
# )
71+
# mo_coeff = get_mo_from_h5(mol, h5_file)
72+
# numpy.savez_compressed(
73+
# "data/o2_mo.npz",
74+
# mol=mol.dumps(),
75+
# mo_coeff=mo_coeff[:, :10],
76+
# )
77+
#
78+
# OpenMolcas output:
79+
#
80+
# :: RASSCF root number 1 Total energy: -149.76505363
81+
# :: RASSCF root number 2 Total energy: -149.53703734
82+
# :: RASSCF root number 3 Total energy: -149.53703734
83+
# :: RASSCF root number 4 Total energy: -149.53215375
84+
# :: RASSCF root number 5 Total energy: -149.38576564
85+
# :: RASSCF root number 6 Total energy: -149.28581775
86+
# :: RASSCF root number 7 Total energy: -149.28581722
87+
# :: RASSCF root number 8 Total energy: -149.15394925
88+
# :: RASSCF root number 9 Total energy: -149.11022925
89+
# :: RASSCF root number 10 Total energy: -149.11022908
90+
# :: RASSCF root number 11 Total energy: -149.07315914
91+
# :: RASSCF root number 12 Total energy: -149.07315679
92+
# :: RASSCF root number 13 Total energy: -149.02767903
93+
# :: RASSCF root number 14 Total energy: -149.02767903
94+
# :: RASSCF root number 15 Total energy: -149.01720780
95+
# :: RASSCF root number 16 Total energy: -149.01720383
96+
# :: RASSCF root number 17 Total energy: -149.00345234
97+
# :: RASSCF root number 18 Total energy: -149.00345234
98+
# :: RASSCF root number 19 Total energy: -148.99747740
99+
# :: RASSCF root number 20 Total energy: -148.99747411
100+
101+
102+
openmolcas_energies = numpy.asarray([
103+
-149.76505363,
104+
-149.53703734,
105+
-149.53703734,
106+
-149.53215375,
107+
-149.38576564,
108+
-149.28581775,
109+
-149.28581722,
110+
-149.15394925,
111+
-149.11022925,
112+
-149.11022908,
113+
-149.07315914,
114+
-149.07315679,
115+
-149.02767903,
116+
-149.02767903,
117+
-149.01720780,
118+
-149.01720383,
119+
-149.00345234,
120+
-149.00345234,
121+
-148.99747740,
122+
-148.99747411,
123+
])
124+
125+
here = Path(__file__).resolve().parent
126+
data_file = here / "data" / "o2_mo.npz"
127+
128+
# Load the exact PySCF molecule and fixed OpenMolcas MO coefficients.
129+
with numpy.load(data_file, allow_pickle=False) as data:
130+
mol = gto.loads(data["mol"].item())
131+
mo_coeff = data["mo_coeff"]
132+
133+
mol.max_memory = 8000 # MB
134+
135+
# One-cycle ROHF for mf generation.
136+
mf = scf.ROHF(mol)
137+
mf.max_cycle = 1
138+
mf.kernel()
139+
mf.mo_coeff = mo_coeff
140+
141+
# Active space
142+
ncore = 2
143+
ncas = 8
144+
nelecas = (7, 5)
145+
146+
gas_orbs = (2, 4, 2)
147+
148+
# OpenMolcas cumulative GAS constraints.
149+
gas_restr = (
150+
(2, 4),
151+
(8, 10),
152+
(12, 12),
153+
)
154+
gas_restr_type = "cumulative-occ"
155+
156+
# GASCI
157+
mc = gasci.GASCI(
158+
mf,
159+
ncas,
160+
nelecas,
161+
ncore=ncore,
162+
gas_orbs=gas_orbs,
163+
gas_restr=gas_restr,
164+
gas_restr_type=gas_restr_type,
165+
)
166+
167+
nroots = 30
168+
169+
mc.verbose = 4
170+
mc.canonicalization = False
171+
mc.fcisolver.spin = 2
172+
mc.fcisolver.nroots = nroots
173+
mc.fcisolver.max_cycle = 300
174+
mc.fcisolver.conv_tol = 1e-10
175+
176+
mc.kernel()
177+
178+
# The PySCF M_S = 1 space contains states with different total spin.
179+
# Select the triplet roots using <S^2> = 2.
180+
pyscf_triplet_energies = []
181+
182+
print()
183+
print(f"PySCF triplet roots among the lowest {nroots} states")
184+
print("state E(PySCF) / Eh <S^2>")
185+
186+
for state, energy in enumerate(numpy.asarray(mc.e_tot).reshape(-1)):
187+
ss, _ = mc.spin_square(state=state)
188+
if abs(ss - 2.0) < 1e-6:
189+
pyscf_triplet_energies.append(energy)
190+
print(f"{state + 1:5d} {energy:22.12f} {ss:12.8f}")
191+
192+
ncompare = min(
193+
len(pyscf_triplet_energies),
194+
len(openmolcas_energies),
195+
)
196+
if ncompare == 0:
197+
raise RuntimeError(
198+
f"No triplet root was found among the lowest {nroots} states"
199+
)
200+
201+
pyscf_triplet_energies = numpy.asarray(
202+
pyscf_triplet_energies[:ncompare]
203+
)
204+
reference_energies = openmolcas_energies[:ncompare]
205+
differences = pyscf_triplet_energies - reference_energies
206+
207+
print()
208+
print("OpenMolcas/PySCF GASCI comparison")
209+
print("root E(OpenMolcas) / Eh E(PySCF) / Eh diff / Eh")
210+
211+
for root, (e_ref, e_pyscf, diff) in enumerate(
212+
zip(reference_energies, pyscf_triplet_energies, differences),
213+
start=1,
214+
):
215+
print(f"{root:4d} {e_ref:22.12f} {e_pyscf:22.12f} {diff: .3e}")
216+
217+
max_diff = numpy.max(numpy.abs(differences))
218+
print()
219+
print(f"Compared roots : {ncompare}")
220+
print(f"max |difference| : {max_diff:.3e} Eh")
221+
222+
# Example output (final comparison section only):
223+
#
224+
# Numerical values in the last few digits may depend on the platform.
225+
#
226+
# GASCI converged
227+
#
228+
# OpenMolcas/PySCF GASCI comparison
229+
# root E(OpenMolcas) / Eh E(PySCF) / Eh diff / Eh
230+
# 1 -149.765053630000 -149.765053626447 3.553e-09
231+
# 2 -149.537037340000 -149.537037344145 -4.145e-09
232+
# 3 -149.537037340000 -149.537037343830 -3.830e-09
233+
# 4 -149.532153750000 -149.532153753145 -3.145e-09
234+
# 5 -149.385765640000 -149.385765640052 -5.247e-11
235+
# 6 -149.285817750000 -149.285817749263 7.374e-10
236+
# 7 -149.285817220000 -149.285817220731 -7.307e-10
237+
# 8 -149.153949250000 -149.153949254242 -4.242e-09
238+
# 9 -149.110229250000 -149.110229248524 1.476e-09
239+
# 10 -149.110229080000 -149.110229082511 -2.511e-09
240+
# 11 -149.073159140000 -149.073159141160 -1.160e-09
241+
# 12 -149.073156790000 -149.073156787594 2.406e-09
242+
# 13 -149.027679030000 -149.027679032858 -2.858e-09
243+
# 14 -149.027679030000 -149.027679032856 -2.856e-09
244+
# 15 -149.017207800000 -149.017207802629 -2.629e-09
245+
# 16 -149.017203830000 -149.017203833468 -3.468e-09
246+
# 17 -149.003452340000 -149.003452343080 -3.080e-09
247+
# 18 -149.003452340000 -149.003452342985 -2.985e-09
248+
# 19 -148.997477400000 -148.997477403623 -3.623e-09
249+
# 20 -148.997474110000 -148.997474109395 6.051e-10
250+
#
251+
# Compared roots : 20
252+
# max |difference| : 4.242e-09 Eh

0 commit comments

Comments
 (0)