Skip to content

Commit f79ef9d

Browse files
DanielCalero1DanielCalero1
authored
DOCI 3 and 4 -RDM implementation (theochem#83)
Add routine for computing 3,4-RDMs
1 parent 1c85597 commit f79ef9d

File tree

9 files changed

+3884
-3
lines changed

9 files changed

+3884
-3
lines changed

‎pyci/__init__.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,10 +19,10 @@
1919
from pyci._pyci import secondquant_op, wavefunction, one_spin_wfn, two_spin_wfn
2020
from pyci._pyci import doci_wfn, fullci_wfn, genci_wfn, sparse_op
2121
from pyci._pyci import get_num_threads, set_num_threads, popcnt, ctz
22-
from pyci._pyci import compute_overlap, compute_rdms, compute_transition_rdms
22+
from pyci._pyci import compute_overlap, compute_rdms, compute_transition_rdms,compute_rdms_1234
2323
from pyci._pyci import add_hci, compute_enpt2
2424

25-
from pyci.utility import make_senzero_integrals, reduce_senzero_integrals, spinize_rdms
25+
from pyci.utility import make_senzero_integrals, reduce_senzero_integrals, spinize_rdms,spinize_rdms_1234,spin_free_rdms
2626
from pyci.utility import odometer_one_spin, odometer_two_spin
2727

2828
from pyci.excitation_ci import add_excitations

‎pyci/include/pyci.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -221,6 +221,8 @@ void clearbit_det(const long, ulong *);
221221

222222
void compute_rdms(const DOCIWfn &, const double *, double *, double *);
223223

224+
void compute_rdms_1234(const DOCIWfn &, const double *, double *, double *, double *, double *);
225+
224226
void compute_rdms(const FullCIWfn &, const double *, double *, double *);
225227

226228
void compute_rdms(const GenCIWfn &, const double *, double *, double *);
@@ -252,6 +254,8 @@ long py_ctz(const Array<ulong>);
252254

253255
pybind11::tuple py_compute_rdms_doci(const DOCIWfn &, const Array<double>);
254256

257+
pybind11::tuple py_compute_rdms_1234_doci(const DOCIWfn &, const Array<double>);
258+
255259
pybind11::tuple py_compute_rdms_fullci(const FullCIWfn &, const Array<double>);
256260

257261
pybind11::tuple py_compute_rdms_genci(const GenCIWfn &, const Array<double>);

‎pyci/src/binding.cpp‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1246,6 +1246,46 @@ For Generalized CI wave functions, ``rdm1`` and ``rdm2`` are the full 1-RDM and
12461246
)""",
12471247
py::arg("wfn"), py::arg("coeffs"));
12481248

1249+
m.def("compute_rdms_1234", &py_compute_rdms_1234_doci, R"""(
1250+
Compute the one-, two-, three-, and four-particle(not implemented yet) reduced density matrices (RDMs) of a wave function.
1251+
1252+
Parameters
1253+
----------
1254+
wfn : pyci.wavefunction
1255+
Wave function.
1256+
coeffs : numpy.ndarray
1257+
Coefficient vector.
1258+
1259+
Returns
1260+
-------
1261+
For DOCI wave functions, this method returns two nbasis-by-nbasis matrices and two nbasis-by-nbasis-by-nbasis,
1262+
which include the unique seniority-zero and seniority-two terms from the full 2-RDMs:
1263+
1264+
.. math::
1265+
1266+
D_0 = \left<pp|qq\right>
1267+
1268+
.. math::
1269+
1270+
D_2 = \left<pq|pq\right>
1271+
1272+
.. math::
1273+
1274+
D_3 = \left<pqr|pqr\right>
1275+
1276+
.. math::
1277+
1278+
D_4 = \left<pqq|prr\right>
1279+
1280+
1281+
Notes
1282+
-----
1283+
This method calculates all RDMs from 1-RDM to 4-RDM (so far implemented only to 3RDM) for a given DOCI wave function.
1284+
1285+
)""",
1286+
py::arg("wfn"), py::arg("coeffs"));
1287+
1288+
12491289
m.def("compute_rdms", &py_compute_rdms_fullci, py::arg("wfn"), py::arg("coeffs"));
12501290

12511291
m.def("compute_rdms", &py_compute_rdms_genci, py::arg("wfn"), py::arg("coeffs"));

‎pyci/src/rdm.cpp‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,95 @@ void compute_rdms(const DOCIWfn &wfn, const double *coeffs, double *d0, double *
6464
}
6565
}
6666

67+
68+
void compute_rdms_1234(const DOCIWfn &wfn, const double *coeffs, double *d0, double *d2, double *d3, double *d4) {
69+
// prepare working vectors
70+
AlignedVector<ulong> v_det(wfn.nword);
71+
AlignedVector<long> v_occs(wfn.nocc_up);
72+
AlignedVector<long> v_virs(wfn.nvir_up);
73+
ulong *det = &v_det[0];
74+
long *occs = &v_occs[0], *virs = &v_virs[0];
75+
// fill rdms with zeros
76+
long i = wfn.nbasis * wfn.nbasis, j = 0;
77+
while (j < i) {
78+
d0[j] = 0;
79+
d2[j++] = 0;
80+
}
81+
i = wfn.nbasis * wfn.nbasis * wfn.nbasis;
82+
j = 0;
83+
while (j < i) {
84+
d3[j]=0;
85+
d4[j++]=0;
86+
}
87+
// iterate over determinants
88+
for (long idet = 0, jdet, mdet, k, l, m, n; idet < wfn.ndet; ++idet) {
89+
double val1, val2;
90+
// fill working vectors
91+
wfn.copy_det(idet, det);
92+
fill_occs(wfn.nword, det, occs);
93+
fill_virs(wfn.nword, wfn.nbasis, det, virs);
94+
// diagonal elements
95+
val1 = coeffs[idet] * coeffs[idet];
96+
for (i = 0; i < wfn.nocc_up; ++i) {
97+
k = occs[i];
98+
d0[k * (wfn.nbasis + 1)] += val1;
99+
for (j = i + 1; j < wfn.nocc_up; ++j) {
100+
l = occs[j];
101+
d2[wfn.nbasis * k + l] += val1;
102+
d2[wfn.nbasis * l + k] += val1;
103+
for (m= j + 1; m < wfn.nocc_up; ++m){
104+
n = occs[m];
105+
d3[(wfn.nbasis * wfn.nbasis) * k + wfn.nbasis * l + n] += val1;
106+
d3[(wfn.nbasis * wfn.nbasis) * k + wfn.nbasis * n + l] += val1;
107+
d3[(wfn.nbasis * wfn.nbasis) * l + wfn.nbasis * k + n] += val1;
108+
d3[(wfn.nbasis * wfn.nbasis) * l + wfn.nbasis * n + k] += val1;
109+
d3[(wfn.nbasis * wfn.nbasis) * n + wfn.nbasis * k + l] += val1;
110+
d3[(wfn.nbasis * wfn.nbasis) * n + wfn.nbasis * l + k] += val1;
111+
}
112+
// pair excitation elements 3rdm j>i
113+
for (m = 0; m < wfn.nvir_up; ++m) {
114+
n = virs[m];
115+
excite_det(l, n, det);
116+
mdet = wfn.index_det(det);
117+
excite_det(n, l, det);
118+
// check if excited determinant is in wfn
119+
if (mdet > idet) {
120+
val2 = coeffs[mdet] * coeffs[idet];
121+
d4[(wfn.nbasis * wfn.nbasis) * k + wfn.nbasis * l + n] += val2;
122+
d4[(wfn.nbasis * wfn.nbasis) * k + wfn.nbasis * n + l] += val2;
123+
}
124+
} //pair excitation elements 3rdm i > j
125+
for (m = 0; m < wfn.nvir_up; ++m) {
126+
n = virs[m];
127+
excite_det(k, n, det);
128+
mdet = wfn.index_det(det);
129+
excite_det(n, k, det);
130+
// check if excited determinant is in wfn
131+
if (mdet > idet) {
132+
val2 = coeffs[mdet] * coeffs[idet];
133+
d4[(wfn.nbasis * wfn.nbasis) * l + wfn.nbasis * k + n] += val2;
134+
d4[(wfn.nbasis * wfn.nbasis) * l + wfn.nbasis * n + k] += val2;
135+
}
136+
}
137+
138+
}
139+
// pair excitation elements
140+
for (j = 0; j < wfn.nvir_up; ++j) {
141+
l = virs[j];
142+
excite_det(k, l, det);
143+
jdet = wfn.index_det(det);
144+
excite_det(l, k, det);
145+
// check if excited determinant is in wfn
146+
if (jdet > idet) {
147+
val2 = coeffs[idet] * coeffs[jdet];
148+
d0[wfn.nbasis * k + l] += val2;
149+
d0[wfn.nbasis * l + k] += val2;
150+
}
151+
}
152+
}
153+
}
154+
}
155+
67156
void compute_rdms(const FullCIWfn &wfn, const double *coeffs, double *rdm1, double *rdm2) {
68157
long n1 = wfn.nbasis;
69158
long n2 = wfn.nbasis * wfn.nbasis;
@@ -815,6 +904,19 @@ pybind11::tuple py_compute_rdms_doci(const DOCIWfn &wfn, const Array<double> coe
815904
return pybind11::make_tuple(d0, d2);
816905
}
817906

907+
pybind11::tuple py_compute_rdms_1234_doci(const DOCIWfn &wfn, const Array<double> coeffs) {
908+
Array<double> d0({wfn.nbasis, wfn.nbasis});
909+
Array<double> d2({wfn.nbasis, wfn.nbasis});
910+
Array<double> d3({wfn.nbasis, wfn.nbasis, wfn.nbasis});
911+
Array<double> d4({wfn.nbasis, wfn.nbasis, wfn.nbasis});
912+
compute_rdms_1234(wfn, reinterpret_cast<const double *>(coeffs.request().ptr),
913+
reinterpret_cast<double *>(d0.request().ptr),
914+
reinterpret_cast<double *>(d2.request().ptr),
915+
reinterpret_cast<double *>(d3.request().ptr),
916+
reinterpret_cast<double *>(d4.request().ptr));
917+
return pybind11::make_tuple(d0, d2, d3, d4);
918+
}
919+
818920
pybind11::tuple py_compute_rdms_fullci(const FullCIWfn &wfn, const Array<double> coeffs) {
819921
Array<double> rdm1({static_cast<long>(2), wfn.nbasis, wfn.nbasis});
820922
Array<double> rdm2({static_cast<long>(3), wfn.nbasis, wfn.nbasis, wfn.nbasis, wfn.nbasis});

‎pyci/test/data/BH_sto-3g_eq.fcidump‎

Lines changed: 194 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,194 @@
1+
&FCI NORB= 6,NELEC= 6,MS2=0,
2+
ORBSYM=1,1,1,1,1,1,
3+
ISYM=1,
4+
&END
5+
2.890694185551886 1 1 1 1
6+
-0.2285478644008276 1 1 2 1
7+
0.5953083364634656 1 1 2 2
8+
-0.2162016548582802 1 1 3 1
9+
0.09104871057910127 1 1 3 2
10+
0.6900652953132707 1 1 3 3
11+
0.7419889025330183 1 1 4 4
12+
0.7419889025330181 1 1 5 5
13+
0.1413257075902383 1 1 6 1
14+
-0.1861902551212619 1 1 6 2
15+
0.09008046930667073 1 1 6 3
16+
0.6699597694264849 1 1 6 6
17+
-0.2285478644008276 2 1 1 1
18+
0.02992401570960048 2 1 2 1
19+
-0.001146336251477235 2 1 2 2
20+
0.02347037645652453 2 1 3 1
21+
-0.008237570230799412 2 1 3 2
22+
-0.01537296998854835 2 1 3 3
23+
-0.008814449236465278 2 1 4 4
24+
-0.008814449236465275 2 1 5 5
25+
-0.02314657520492044 2 1 6 1
26+
0.003050017645768165 2 1 6 2
27+
-0.00419211917379873 2 1 6 3
28+
-0.004379253657977653 2 1 6 6
29+
0.5953083364634657 2 2 1 1
30+
-0.001146336251477227 2 2 2 1
31+
0.5804858409021905 2 2 2 2
32+
-0.01900821896163675 2 2 3 1
33+
-0.07241056645830383 2 2 3 2
34+
0.4170856180274074 2 2 3 3
35+
0.4593408467391383 2 2 4 4
36+
0.4593408467391382 2 2 5 5
37+
-0.01238748160972113 2 2 6 1
38+
0.0002084163025978247 2 2 6 2
39+
-0.05254431697407079 2 2 6 3
40+
0.5754116421934867 2 2 6 6
41+
-0.2162016548582803 3 1 1 1
42+
0.02347037645652455 3 1 2 1
43+
-0.01900821896163674 3 1 2 2
44+
0.03397697977407958 3 1 3 1
45+
0.003163633505588129 3 1 3 2
46+
0.00404893468739087 3 1 3 3
47+
-0.007412212106220699 3 1 4 4
48+
-0.007412212106220697 3 1 5 5
49+
-0.005247409917599474 3 1 6 1
50+
0.01158307545322849 3 1 6 2
51+
0.006990498122004285 3 1 6 3
52+
-0.01879243966504147 3 1 6 6
53+
0.09104871057910127 3 2 1 1
54+
-0.008237570230799426 3 2 2 1
55+
-0.07241056645830381 3 2 2 2
56+
0.003163633505588093 3 2 3 1
57+
0.0676156185213167 3 2 3 2
58+
0.05019628841726888 3 2 3 3
59+
0.04621563081453394 3 2 4 4
60+
0.04621563081453393 3 2 5 5
61+
0.01250097777451356 3 2 6 1
62+
-0.05223141555214424 3 2 6 2
63+
0.05747267715242108 3 2 6 3
64+
-0.06327529964411363 3 2 6 6
65+
0.6900652953132709 3 3 1 1
66+
-0.01537296998854832 3 3 2 1
67+
0.4170856180274075 3 3 2 2
68+
0.004048934687390883 3 3 3 1
69+
0.05019628841726886 3 3 3 2
70+
0.5635811972005504 3 3 3 3
71+
0.5041862576414032 3 3 4 4
72+
0.5041862576414031 3 3 5 5
73+
0.0231276130620684 3 3 6 1
74+
-0.03192699042073559 3 3 6 2
75+
0.07541684244416474 3 3 6 3
76+
0.4565326809201556 3 3 6 6
77+
0.02169507828733871 4 1 4 1
78+
0.01921155457995589 4 1 4 2
79+
0.01851865219375422 4 1 4 3
80+
-0.01161252859365546 4 1 6 4
81+
0.01921155457995589 4 2 4 1
82+
0.06039821130266281 4 2 4 2
83+
0.04342107961102947 4 2 4 3
84+
-0.03633847993724627 4 2 6 4
85+
0.01851865219375421 4 3 4 1
86+
0.04342107961102945 4 3 4 2
87+
0.06291076116021634 4 3 4 3
88+
-0.01322878606419047 4 3 6 4
89+
0.7419889025330183 4 4 1 1
90+
-0.008814449236465242 4 4 2 1
91+
0.4593408467391384 4 4 2 2
92+
-0.00741221210622064 4 4 3 1
93+
0.04621563081453399 4 4 3 2
94+
0.5041862576414033 4 4 3 3
95+
0.5867727439789322 4 4 4 4
96+
0.5235134814676911 4 4 5 5
97+
0.004168026575771488 4 4 6 1
98+
-0.09238076313811244 4 4 6 2
99+
0.04177944228679972 4 4 6 3
100+
0.4740773398260829 4 4 6 6
101+
0.0216950782873387 5 1 5 1
102+
0.01921155457995588 5 1 5 2
103+
0.01851865219375421 5 1 5 3
104+
-0.01161252859365545 5 1 6 5
105+
0.01921155457995588 5 2 5 1
106+
0.06039821130266281 5 2 5 2
107+
0.04342107961102946 5 2 5 3
108+
-0.03633847993724625 5 2 6 5
109+
0.01851865219375421 5 3 5 1
110+
0.04342107961102945 5 3 5 2
111+
0.06291076116021634 5 3 5 3
112+
-0.01322878606419046 5 3 6 5
113+
0.03162963125562031 5 4 5 4
114+
0.7419889025330181 5 5 1 1
115+
-0.008814449236465207 5 5 2 1
116+
0.4593408467391384 5 5 2 2
117+
-0.007412212106220606 5 5 3 1
118+
0.04621563081453402 5 5 3 2
119+
0.5041862576414033 5 5 3 3
120+
0.5235134814676911 5 5 4 4
121+
0.586772743978932 5 5 5 5
122+
0.004168026575771471 5 5 6 1
123+
-0.09238076313811253 5 5 6 2
124+
0.04177944228679967 5 5 6 3
125+
0.4740773398260828 5 5 6 6
126+
0.1413257075902383 6 1 1 1
127+
-0.02314657520492046 6 1 2 1
128+
-0.01238748160972115 6 1 2 2
129+
-0.005247409917599475 6 1 3 1
130+
0.01250097777451353 6 1 3 2
131+
0.02312761306206843 6 1 3 3
132+
0.004168026575771572 6 1 4 4
133+
0.00416802657577157 6 1 5 5
134+
0.02880391396937533 6 1 6 1
135+
0.005524338027076478 6 1 6 2
136+
0.01082364008038884 6 1 6 3
137+
-0.008812192154253081 6 1 6 6
138+
-0.186190255121262 6 2 1 1
139+
0.003050017645768176 6 2 2 1
140+
0.0002084163025976651 6 2 2 2
141+
0.01158307545322851 6 2 3 1
142+
-0.05223141555214433 6 2 3 2
143+
-0.03192699042073567 6 2 3 3
144+
-0.0923807631381124 6 2 4 4
145+
-0.09238076313811237 6 2 5 5
146+
0.005524338027076496 6 2 6 1
147+
0.1212730361895562 6 2 6 2
148+
-0.05403162334477736 6 2 6 3
149+
0.01293203355693659 6 2 6 6
150+
0.09008046930667071 6 3 1 1
151+
-0.004192119173798695 6 3 2 1
152+
-0.05254431697407061 6 3 2 2
153+
0.006990498122004348 6 3 3 1
154+
0.05747267715242119 6 3 3 2
155+
0.07541684244416491 6 3 3 3
156+
0.0417794422867997 6 3 4 4
157+
0.04177944228679968 6 3 5 5
158+
0.01082364008038878 6 3 6 1
159+
-0.05403162334477741 6 3 6 2
160+
0.07634911021780025 6 3 6 3
161+
-0.04207173170262291 6 3 6 6
162+
-0.01161252859365545 6 4 4 1
163+
-0.03633847993724626 6 4 4 2
164+
-0.01322878606419046 6 4 4 3
165+
0.03249943702929813 6 4 6 4
166+
-0.01161252859365545 6 5 5 1
167+
-0.03633847993724625 6 5 5 2
168+
-0.01322878606419045 6 5 5 3
169+
0.03249943702929813 6 5 6 5
170+
0.6699597694264852 6 6 1 1
171+
-0.004379253657977538 6 6 2 1
172+
0.5754116421934872 6 6 2 2
173+
-0.01879243966504134 6 6 3 1
174+
-0.06327529964411353 6 6 3 2
175+
0.4565326809201556 6 6 3 3
176+
0.4740773398260831 6 6 4 4
177+
0.474077339826083 6 6 5 5
178+
-0.008812192154253214 6 6 6 1
179+
0.01293203355693684 6 6 6 2
180+
-0.04207173170262306 6 6 6 3
181+
0.6101411445973578 6 6 6 6
182+
-12.73694583468452 1 1 0 0
183+
0.2636037738346674 2 1 0 0
184+
-3.08122009866082 2 2 0 0
185+
0.2419315877376549 3 1 0 0
186+
-0.1364127650686597 3 2 0 0
187+
-2.922828142214635 3 3 0 0
188+
-2.996085190675585 4 4 0 0
189+
-2.996085190675584 5 5 0 0
190+
-0.1527654545391036 6 1 0 0
191+
0.4703521766204644 6 2 0 0
192+
-0.2079679737989178 6 3 0 0
193+
-2.475899716280425 6 6 0 0
194+
2.147634784577923 0 0 0 0

0 commit comments

Comments
 (0)