Skip to content

Commit 263450f

Browse files
alexvakimovalexvakimov
committed
Added the option to use the active space in SD, CSF, and CI overlap calculations, as well as in the MOPAC workflow, added the unittest examples to show that
1 parent 9cad9ee commit 263450f

4 files changed

Lines changed: 294 additions & 3 deletions

File tree

‎src/libra_py/citools/ci.py‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,9 @@ def overlap(st_mo, data1, data2, params):
6565
Total number of electrons
6666
nstates : int
6767
Total number of electronic states, including the ground state
68+
active_space : iterable of int, optional
69+
Spatial orbital indices (1-based) defining the active space for
70+
spin adaptation. If None, all available orbitals are used.
6871
6972
Returns
7073
-------
@@ -94,6 +97,7 @@ def overlap(st_mo, data1, data2, params):
9497
nvirt = params["nvirt"]
9598
nelec = params["nelec"]
9699
nstates = params["nstates"]
100+
active_space = params["active_space"]
97101

98102
if nstates <= 1:
99103
raise ValueError("nstates must include at least one excited state")
@@ -131,6 +135,7 @@ def overlap(st_mo, data1, data2, params):
131135
nelec,
132136
homo_indx,
133137
common_sd_basis,
138+
active_space
134139
)
135140

136141
# ------------------------------------------------------------------

‎src/libra_py/packages/mopac/methods.py‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -534,7 +534,8 @@ def mopac_compute_adi(q, params, full_id):
534534
coords = q.col(itraj)
535535

536536
critical_params = ["labels", "timestep"]
537-
default_params = {"mopac_exe": "mopac", "is_first_time": True, "orbital_space": None, "nstates":2,
537+
default_params = {"mopac_exe": "mopac", "is_first_time": True, "orbital_space": None, "active_space": None,
538+
"nstates":2,
538539
"mopac_run_params": "INDO C.I.=(6,3) CHARGE=0 RELSCF=0.000001 ALLVEC WRTCONF=0.00 WRTCI=2",
539540
"mopac_working_directory": "mopac_wd",
540541
"mopac_input_prefix": "input_", "mopac_output_prefix": "output_",
@@ -553,6 +554,7 @@ def mopac_compute_adi(q, params, full_id):
553554
mopac_input_prefix = params[itraj]["mopac_input_prefix"]
554555
mopac_output_prefix = params[itraj]["mopac_output_prefix"]
555556
orbital_space = params[itraj]["orbital_space"]
557+
active_space = params[itraj]["active_space"]
556558
nstates = params[itraj]["nstates"]
557559
dt = params[itraj]["dt"]
558560
do_Lowdin = params[itraj]["do_Lowdin"]
@@ -604,8 +606,9 @@ def mopac_compute_adi(q, params, full_id):
604606
ovlp_params = {"homo_indx":info["nocc"],
605607
"nocc":info["nocc"] - 1,
606608
"nvirt":info["nmo"] - info["nocc"],
607-
"nelec":info["nelec"], "nstates":nstates }
608-
609+
"nelec":info["nelec"], "nstates":nstates,
610+
"active_space":active_space
611+
}
609612
st_ci = ci.overlap(st_mo, data_prev, data_curr, ovlp_params)
610613

611614
#=============== Now, populate the allocated matrices ======================

‎unittests/test_ci_overlaps_equivalence.py‎

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
import numpy as np
2+
import pytest
3+
4+
from libra_py.citools.ci import overlap
5+
6+
7+
# ============================================================
8+
# Example runner (can be used outside pytest)
9+
# ============================================================
10+
11+
def run_examples():
12+
results = []
13+
14+
# ================= Case 1 ======================
15+
print("Case 1")
16+
s = np.array([
17+
[1.0, 0.2],
18+
[0.2, 1.2],
19+
])
20+
S = np.kron(np.eye(2), s)
21+
22+
# lowest = homo - nocc
23+
# highest = homo + nvirt
24+
params = dict(
25+
S=S,
26+
nelec=2,
27+
nocc=0,
28+
nvirt=1,
29+
homo_indx=1,
30+
nstates = 2,
31+
data = [ [],
32+
[ [[1, 2]] ],
33+
[ [1.0 ] ]
34+
],
35+
active_space=None,
36+
)
37+
38+
st_ci = overlap(params["S"], params["data"], params["data"], params)
39+
40+
print("CI overlap:\n", st_ci)
41+
results.append(st_ci)
42+
43+
# ================= Case 2 ======================
44+
print("\nCase 2")
45+
s = np.array([
46+
[-1.0, 0.0, 0.0],
47+
[0.0, 1.0, 0.2],
48+
[0.0, 0.2, 1.2],
49+
])
50+
S = np.kron(np.eye(2), s)
51+
52+
params = dict(
53+
S=S,
54+
nelec=4,
55+
nocc=1,
56+
nvirt=1,
57+
homo_indx=2,
58+
nstates = 2,
59+
data = [ [],
60+
[ [[2, 3]] ],
61+
[ [1.0 ] ]
62+
],
63+
active_space=None,
64+
)
65+
66+
st_ci = overlap(params["S"], params["data"], params["data"], params)
67+
68+
print("CI overlap:\n", st_ci)
69+
results.append(st_ci)
70+
71+
# ================= Case 3 ======================
72+
print("\nCase 3")
73+
s = np.array([
74+
[-1.0, 0.0, 0.0],
75+
[0.0, 1.0, 0.2],
76+
[0.0, 0.2, 1.2],
77+
])
78+
S = np.kron(np.eye(2), s)
79+
80+
params = dict(
81+
S=S,
82+
nelec=2,
83+
nocc=1,
84+
nvirt=1,
85+
homo_indx=2,
86+
nstates = 2,
87+
data = [ [],
88+
[ [[2, 3]] ],
89+
[ [1.0 ] ]
90+
],
91+
active_space=[2, 3],
92+
)
93+
94+
st_ci = overlap(params["S"], params["data"], params["data"], params)
95+
96+
print("CI overlap:\n", st_ci)
97+
results.append(st_ci)
98+
99+
return results
100+
101+
102+
# Uncomment to run the example
103+
#run_examples()
104+
105+
106+
# ============================================================
107+
# Pytest tests
108+
# ============================================================
109+
110+
def test_all_cases_produce_identical_overlaps():
111+
"""
112+
All three examples are physically equivalent constructions
113+
and must produce identical CI overlap matrices.
114+
"""
115+
results = run_examples()
116+
117+
ref_ci = results[0]
118+
119+
for i, st_ci in enumerate(results[1:], start=2):
120+
assert np.allclose(
121+
st_ci, ref_ci, atol=1e-12
122+
), f"CI overlap differs in case {i}"
123+
124+
125+

‎unittests/test_sd_and_csf_overlaps_equivalence.py‎

Lines changed: 158 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,158 @@
1+
import numpy as np
2+
import pytest
3+
4+
from libra_py.citools.interfaces import sd_and_csf_overlaps_singlet
5+
6+
7+
# ============================================================
8+
# Example runner (can be used outside pytest)
9+
# ============================================================
10+
11+
def run_examples():
12+
results = []
13+
14+
# ================= Case 1 ======================
15+
print("Case 1")
16+
s = np.array([
17+
[1.0, 0.2],
18+
[0.2, 1.2],
19+
])
20+
S = np.kron(np.eye(2), s)
21+
22+
params = dict(
23+
S=S,
24+
lowest_orbital=1,
25+
highest_orbital=2,
26+
nelec=2,
27+
homo_indx=1,
28+
common_sd_basis=[[1, 2]],
29+
active=None,
30+
)
31+
32+
st_csf, st_sd = sd_and_csf_overlaps_singlet(
33+
params["S"],
34+
params["lowest_orbital"],
35+
params["highest_orbital"],
36+
params["nelec"],
37+
params["homo_indx"],
38+
params["common_sd_basis"],
39+
params["active"],
40+
)
41+
42+
print("CSF overlap:\n", st_csf)
43+
print("SD overlap:\n", st_sd)
44+
results.append((st_csf, st_sd))
45+
46+
# ================= Case 2 ======================
47+
print("\nCase 2")
48+
s = np.array([
49+
[-1.0, 0.0, 0.0],
50+
[0.0, 1.0, 0.2],
51+
[0.0, 0.2, 1.2],
52+
])
53+
S = np.kron(np.eye(2), s)
54+
55+
params = dict(
56+
S=S,
57+
lowest_orbital=1,
58+
highest_orbital=3,
59+
nelec=4,
60+
homo_indx=2,
61+
common_sd_basis=[[2, 3]],
62+
active=None,
63+
)
64+
65+
st_csf, st_sd = sd_and_csf_overlaps_singlet(
66+
params["S"],
67+
params["lowest_orbital"],
68+
params["highest_orbital"],
69+
params["nelec"],
70+
params["homo_indx"],
71+
params["common_sd_basis"],
72+
params["active"],
73+
)
74+
75+
print("CSF overlap:\n", st_csf)
76+
print("SD overlap:\n", st_sd)
77+
results.append((st_csf, st_sd))
78+
79+
# ================= Case 3 ======================
80+
print("\nCase 3")
81+
s = np.array([
82+
[-1.0, 0.0, 0.0],
83+
[0.0, 1.0, 0.2],
84+
[0.0, 0.2, 1.2],
85+
])
86+
S = np.kron(np.eye(2), s)
87+
88+
params = dict(
89+
S=S,
90+
lowest_orbital=1,
91+
highest_orbital=3,
92+
nelec=2,
93+
homo_indx=2,
94+
common_sd_basis=[[2, 3]],
95+
active=[2, 3],
96+
)
97+
98+
st_csf, st_sd = sd_and_csf_overlaps_singlet(
99+
params["S"],
100+
params["lowest_orbital"],
101+
params["highest_orbital"],
102+
params["nelec"],
103+
params["homo_indx"],
104+
params["common_sd_basis"],
105+
params["active"],
106+
)
107+
108+
print("CSF overlap:\n", st_csf)
109+
print("SD overlap:\n", st_sd)
110+
results.append((st_csf, st_sd))
111+
112+
return results
113+
114+
115+
# ============================================================
116+
# Pytest tests
117+
# ============================================================
118+
119+
def test_all_cases_produce_identical_overlaps():
120+
"""
121+
All three examples are physically equivalent constructions
122+
and must produce identical SD and CSF overlap matrices.
123+
"""
124+
results = run_examples()
125+
126+
ref_csf, ref_sd = results[0]
127+
128+
for i, (st_csf, st_sd) in enumerate(results[1:], start=2):
129+
assert np.allclose(
130+
st_sd, ref_sd, atol=1e-12
131+
), f"SD overlap differs in case {i}"
132+
133+
assert np.allclose(
134+
st_csf, ref_csf, atol=1e-12
135+
), f"CSF overlap differs in case {i}"
136+
137+
138+
def test_overlap_matrices_are_symmetric_and_square():
139+
results = run_examples()
140+
141+
for st_csf, st_sd in results:
142+
assert st_sd.ndim == 2
143+
assert st_csf.ndim == 2
144+
145+
assert st_sd.shape[0] == st_sd.shape[1]
146+
assert st_csf.shape[0] == st_csf.shape[1]
147+
148+
assert np.allclose(st_sd, st_sd.T, atol=1e-12)
149+
assert np.allclose(st_csf, st_csf.T, atol=1e-12)
150+
151+
152+
def test_example_runner_produces_output(capsys):
153+
run_examples()
154+
out = capsys.readouterr().out
155+
assert "Case 1" in out
156+
assert "Case 2" in out
157+
assert "Case 3" in out
158+

0 commit comments

Comments
 (0)