Skip to content

Commit 6da6b57

Browse files
arshadmehmood118Arshad Mehmoodsunqm
authored
Add FOMO-CASCI and CASCI with DFT embedding extensions (#167)
* Add FOMO-CASCI and CASCI with DFT embedding extensions New Features: - FOMO-SCF: pyscf/mcscf/addons.py - DFT-CASCI: pyscf/mcscf/casci_dft.py - Gradients: pyscf/grad/casci_dft.py - Tests: pyscf/mcscf/test/, pyscf/grad/test/ - Examples: examples/mcscf/, examples/grad/ References: - P. Slavíček, T. J. Martínez, J. Chem. Phys. 132, 234102 (2010) - S. Pijeau, E. G. Hohenstein, J. Chem. Theory Comput. 13, 1130 (2017) Author: Arshad Mehmood <arshad.mehmood@stonybrook.edu> * Address reviewer feedback: rename modules and reorganize FOMO-SCF - Move FOMO-SCF from mcscf/addons.py to scf/fomoscf.py (reviewer comment #1) - Rename casci_dft -> dft_corrected_casci throughout codebase (reviewer comment #2) - Update all imports and __init__.py files accordingly - Rename test files to match new module names - Update example files with new naming convention - Use 'DFT-corrected CASCI' terminology consistently in code, comments, and outputs This addresses both reviewer comments from the initial PR review. * Delete pyscf/grad/__init__.py to avoid overwriting pyscf core. * Delete pyscf/mcscf/__init__.py to avoid overwriting pyscf core. * Delete pyscf/scf/__init__.py to avoid overwriting pyscf core. * Adjust pytest configurations * Lint * Rename test files for pytest --------- Co-authored-by: Arshad Mehmood <amehmood@xeonmax.cm.cluster> Co-authored-by: Qiming Sun <osirpt.sun@gmail.com>
1 parent 23e79fa commit 6da6b57

11 files changed

Lines changed: 2879 additions & 0 deletions

‎examples/grad/20_dft_corrected_casci_grad.py‎

Lines changed: 166 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,166 @@
1+
#!/usr/bin/env python
2+
"""
3+
Example 20: Nuclear Gradients for DFT-Corrected CASCI
4+
=====================================================
5+
6+
This example demonstrates nuclear gradient calculations for CASCI with
7+
DFT-corrected core energy using BOTH analytical and numerical methods.
8+
9+
Gradient Methods:
10+
1. Analytical: Fast (~1-2e-2 Ha/Bohr accuracy)
11+
2. Numerical: Slow but accurate (~1e-6 Ha/Bohr accuracy)
12+
3. Auto: Tries analytical, falls back to numerical if needed
13+
14+
Author:
15+
Arshad Mehmood, IACS, Stony Brook University
16+
arshad.mehmood@stonybrook.edu
17+
31 December 2025
18+
"""
19+
20+
from pyscf import gto, scf, mcscf
21+
from pyscf.mcscf import dft_corrected_casci
22+
import numpy as np
23+
import time
24+
25+
# Build water molecule
26+
mol = gto.M(
27+
atom='''
28+
O 0.000 0.000 0.000
29+
H 0.000 0.757 0.587
30+
H 0.000 -0.757 0.587
31+
''',
32+
basis='cc-pvdz',
33+
unit='Angstrom',
34+
verbose=3
35+
)
36+
37+
ncas, nelecas = 6, 6
38+
39+
print("="*70)
40+
print("DFT-corrected CASCI Nuclear Gradients: Analytical vs Numerical")
41+
print("="*70)
42+
43+
# =============================================================================
44+
# Part 1: RHF Reference - Compare Methods
45+
# =============================================================================
46+
print("\n" + "="*70)
47+
print("Part 1: RHF Reference - Method Comparison")
48+
print("="*70)
49+
50+
mf_rhf = scf.RHF(mol).run()
51+
mc = dft_corrected_casci.CASCI(mf_rhf, ncas, nelecas, xc='PBE')
52+
mc.kernel()
53+
54+
print(f"\nCASCI Energy: {mc.e_tot:.10f} Ha")
55+
56+
# Analytical gradient
57+
print("\n1a. Analytical Gradient (fast, ~1-2e-2 Ha/Bohr accuracy):")
58+
t0 = time.time()
59+
g_analytical = mc.Gradients(method='analytical').kernel()
60+
t_analytical = time.time() - t0
61+
print(f" Time: {t_analytical:.2f} seconds")
62+
63+
# Numerical gradient
64+
print("\n1b. Numerical Gradient (slow, ~1e-6 Ha/Bohr accuracy):")
65+
t0 = time.time()
66+
g_numerical = mc.Gradients(method='numerical', step_size=1e-4).kernel()
67+
t_numerical = time.time() - t0
68+
print(f" Time: {t_numerical:.2f} seconds")
69+
print(f" Speedup: {t_numerical/t_analytical:.1f}x faster with analytical")
70+
71+
# Auto mode
72+
print("\n1c. Auto Mode (tries analytical, falls back to numerical):")
73+
g_auto = mc.Gradients(method='auto').kernel()
74+
75+
# =============================================================================
76+
# Part 2: Detailed Comparison
77+
# =============================================================================
78+
print("\n" + "="*70)
79+
print("Part 2: Gradient Comparison (Ha/Bohr)")
80+
print("="*70)
81+
82+
atoms = ['O', 'H1', 'H2']
83+
coords = ['x', 'y', 'z']
84+
85+
print(f"\n{'Atom':<6} {'Coord':<6} {'Analytical':<15} {'Numerical':<15} {'|Difference|':<15}")
86+
print("-"*66)
87+
88+
max_diff = 0.0
89+
for i, atom in enumerate(atoms):
90+
for j, coord in enumerate(coords):
91+
diff = abs(g_analytical[i,j] - g_numerical[i,j])
92+
max_diff = max(max_diff, diff)
93+
print(f"{atom:<6} {coord:<6} {g_analytical[i,j]:>14.8f} {g_numerical[i,j]:>14.8f} {diff:>14.8f}")
94+
95+
print("-"*66)
96+
print(f"Maximum difference: {max_diff:.2e} Ha/Bohr")
97+
print(f"RMS difference: {np.sqrt(np.mean((g_analytical - g_numerical)**2)):.2e} Ha/Bohr")
98+
99+
# =============================================================================
100+
# Part 3: UHF Reference (Numerical Only)
101+
# =============================================================================
102+
print("\n" + "="*70)
103+
print("Part 3: UHF Reference (Numerical Gradients)")
104+
print("="*70)
105+
106+
mf_uhf = scf.UHF(mol).run()
107+
mc_uhf = dft_corrected_casci.UCASCI(mf_uhf, ncas, nelecas, xc='PBE')
108+
mc_uhf.kernel()
109+
110+
print(f"\nUCASCI Energy: {mc_uhf.e_tot:.10f} Ha")
111+
print("\nNote: UCASCI only supports numerical gradients")
112+
113+
t0 = time.time()
114+
g_uhf = mc_uhf.Gradients(method='numerical').kernel()
115+
t_uhf = time.time() - t0
116+
print(f"Time: {t_uhf:.2f} seconds")
117+
118+
print(f"\n{'Atom':<6} {'Coord':<6} {'Gradient':<15}")
119+
print("-"*27)
120+
for i, atom in enumerate(atoms):
121+
for j, coord in enumerate(coords):
122+
print(f"{atom:<6} {coord:<6} {g_uhf[i,j]:>14.8f}")
123+
124+
# =============================================================================
125+
# Part 4: Different XC Functionals
126+
# =============================================================================
127+
print("\n" + "="*70)
128+
print("Part 4: Gradient Comparison Across XC Functionals")
129+
print("="*70)
130+
131+
functionals = ['LDA', 'PBE', 'B3LYP']
132+
results = []
133+
134+
for xc in functionals:
135+
mc_xc = dft_corrected_casci.CASCI(mf_rhf, ncas, nelecas, xc=xc)
136+
mc_xc.kernel()
137+
g_xc = mc_xc.Gradients(method='analytical').kernel()
138+
results.append((xc, mc_xc.e_tot, g_xc))
139+
140+
print(f"\n{'Functional':<12} {'Energy (Ha)':<18} {'Max Gradient':<15}")
141+
print("-"*45)
142+
for xc, energy, grad in results:
143+
max_g = np.abs(grad).max()
144+
print(f"{xc:<12} {energy:<18.10f} {max_g:<15.8f}")
145+
146+
# =============================================================================
147+
# Part 5: Performance Recommendations
148+
# =============================================================================
149+
print("\n" + "="*70)
150+
print("Performance and Accuracy Summary")
151+
print("="*70)
152+
153+
print("\nMethod Characteristics:")
154+
print(f"{'Method':<15} {'Speed':<12} {'Accuracy':<20} {'Use Case'}")
155+
print("-"*70)
156+
print(f"{'Analytical':<15} {'Fast':<12} {'~1-2e-2 Ha/Bohr':<20} {'MD, optimization (loose)'}")
157+
print(f"{'Numerical':<15} {'Slow':<12} {'~1e-6 Ha/Bohr':<20} {'Validation, tight opt'}")
158+
print(f"{'Auto':<15} {'Fast*':<12} {'Variable':<20} {'General purpose'}")
159+
160+
print("\nRecommendations:")
161+
print(" • Use analytical for: Molecular dynamics, exploratory optimization")
162+
print(" • Use numerical for: Final geometries, benchmark calculations")
163+
print(" • Use auto for: General purpose (safe default)")
164+
print(f" • Current speedup: Analytical is {t_numerical/t_analytical:.1f}x faster")
165+
166+
print("\n" + "="*70)

‎examples/grad/21_fomo_casci_grad.py‎

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
#!/usr/bin/env python
2+
"""
3+
Example 21: Nuclear Gradients for FOMO-CASCI (Both Methods)
4+
===========================================================
5+
6+
This example demonstrates gradient calculations for FOMO-CASCI using
7+
both analytical and numerical methods.
8+
9+
Author:
10+
Arshad Mehmood, IACS, Stony Brook University
11+
31 December 2025
12+
"""
13+
14+
from pyscf import gto, scf, mcscf
15+
from pyscf.scf import fomoscf
16+
from pyscf.mcscf import dft_corrected_casci
17+
import numpy as np
18+
import time
19+
20+
mol = gto.M(
21+
atom='O 0 0 0; H 0 0.757 0.587; H 0 -0.757 0.587',
22+
basis='cc-pvdz', unit='Angstrom', verbose=3
23+
)
24+
25+
ncas, nelecas = 6, 6
26+
ncore = (mol.nelectron - nelecas) // 2
27+
28+
print("="*70)
29+
print("FOMO-CASCI Nuclear Gradients: Analytical vs Numerical")
30+
print("="*70)
31+
32+
# Standard RHF
33+
mf_rhf = scf.RHF(mol).run()
34+
35+
# FOMO-SCF
36+
mf_fomo = fomoscf.fomo_scf(
37+
mf_rhf,
38+
temperature=0.25,
39+
method='gaussian',
40+
restricted=(ncore, ncas)
41+
)
42+
mf_fomo.kernel()
43+
print(f"\nFOMO occupations: {mf_fomo.mo_occ}")
44+
45+
# Standard CASCI for comparison
46+
print("\n" + "-"*70)
47+
print("1. Standard CASCI (RHF orbitals)")
48+
print("-"*70)
49+
mc_std = mcscf.CASCI(mf_rhf, ncas, nelecas)
50+
mc_std.kernel()
51+
print(f"Energy: {mc_std.e_tot:.10f} Ha")
52+
53+
# Standard CASCI gradient (only analytical available)
54+
g_std = mc_std.Gradients().kernel()
55+
56+
# FOMO-CASCI
57+
print("\n" + "-"*70)
58+
print("2. FOMO-CASCI (FOMO orbitals, LDA core)")
59+
print("-"*70)
60+
mc_fomo = dft_corrected_casci.CASCI(mf_fomo, ncas, nelecas, xc='LDA')
61+
mc_fomo.kernel()
62+
print(f"Energy: {mc_fomo.e_tot:.10f} Ha")
63+
64+
# Analytical gradient
65+
print("\n2a. Analytical gradient:")
66+
t0 = time.time()
67+
g_fomo_analytical = mc_fomo.Gradients(method='analytical').kernel()
68+
t_analytical = time.time() - t0
69+
print(f" Time: {t_analytical:.2f} seconds")
70+
71+
# Numerical gradient
72+
print("\n2b. Numerical gradient:")
73+
t0 = time.time()
74+
g_fomo_numerical = mc_fomo.Gradients(method='numerical').kernel()
75+
t_numerical = time.time() - t0
76+
print(f" Time: {t_numerical:.2f} seconds")
77+
78+
# Comparison
79+
print("\n" + "="*70)
80+
print("Gradient Comparison")
81+
print("="*70)
82+
83+
atoms = ['O', 'H1', 'H2']
84+
coords = ['x', 'y', 'z']
85+
86+
print(f"\n{'Atom':<6} {'Coord':<6} {'Standard':<15} {'FOMO-Ana':<15} {'FOMO-Num':<15}")
87+
print("-"*57)
88+
89+
for i, atom in enumerate(atoms):
90+
for j, coord in enumerate(coords):
91+
print(f"{atom:<6} {coord:<6} {g_std[i,j]:>14.8f} {g_fomo_analytical[i,j]:>14.8f} {g_fomo_numerical[i,j]:>14.8f}")
92+
93+
# Analysis
94+
diff_std_fomo = np.abs(g_std - g_fomo_analytical).max()
95+
diff_methods = np.abs(g_fomo_analytical - g_fomo_numerical).max()
96+
97+
print("\nAnalysis:")
98+
print(f" Max difference (Standard vs FOMO analytical): {diff_std_fomo:.6f} Ha/Bohr")
99+
print(f" Max difference (Analytical vs Numerical): {diff_methods:.6f} Ha/Bohr")
100+
print(f" Speedup (Analytical over Numerical): {t_numerical/t_analytical:.1f}x")

0 commit comments

Comments
 (0)