Skip to content

DOCI 3RDM implementation #83

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
msricher merged 15 commits into from
Dec 9, 2024
Merged

DOCI 3RDM implementation #83

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
94ff745
Adding and computing d3 and d4 matrices for the 3rdm
DanielCalero1 Dec 3, 2024
bd91d4b
Updatind python binding in rmd.cpp
DanielCalero1 Dec 3, 2024
5fe3111
Updating pyci.h
DanielCalero1 Dec 3, 2024
af06c92
Updating spinze_rdms to compute DOCI 3rdm
DanielCalero1 Dec 3, 2024
453e6b3
Adding spin_free_rdms function
DanielCalero1 Dec 4, 2024
5403be9
Updating __init__.py to include spin_free_rdms function
DanielCalero1 Dec 4, 2024
a719ae5
Separating the previous function compute_rdms from the new compute_rd…
DanielCalero1 Dec 6, 2024
4e5be15
Updating bindings
DanielCalero1 Dec 6, 2024
17e31bc
Correcting typo in rdm.cpp
DanielCalero1 Dec 6, 2024
34ac07c
updating pyci.h to include compute_rdms_34 function
DanielCalero1 Dec 6, 2024
6bdb6a6
Including function compute_rdms_34 in __init__.py
DanielCalero1 Dec 6, 2024
998e0bd
Separating functions spinize_rdms and spinize_rdms_34
DanielCalero1 Dec 6, 2024
0cee298
Creating new tests routines for compute_rdms_34 and spinize_rdms_34
DanielCalero1 Dec 6, 2024
6267ab9
Improving code style
DanielCalero1 Dec 9, 2024
ca1f3af
Changing function name to compute_rdms_1234
DanielCalero1 Dec 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions pyci/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,10 @@
from pyci._pyci import secondquant_op, wavefunction, one_spin_wfn, two_spin_wfn
from pyci._pyci import doci_wfn, fullci_wfn, genci_wfn, sparse_op
from pyci._pyci import get_num_threads, set_num_threads, popcnt, ctz
from pyci._pyci import compute_overlap, compute_rdms, compute_transition_rdms
from pyci._pyci import compute_overlap, compute_rdms, compute_transition_rdms,compute_rdms_1234
from pyci._pyci import add_hci, compute_enpt2

from pyci.utility import make_senzero_integrals, reduce_senzero_integrals, spinize_rdms
from pyci.utility import make_senzero_integrals, reduce_senzero_integrals, spinize_rdms,spinize_rdms_1234,spin_free_rdms
from pyci.utility import odometer_one_spin, odometer_two_spin

from pyci.excitation_ci import add_excitations
Expand Down
4 changes: 4 additions & 0 deletions pyci/include/pyci.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -221,6 +221,8 @@ void clearbit_det(const long, ulong *);

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

void compute_rdms_1234(const DOCIWfn &, const double *, double *, double *, double *, double *);

void compute_rdms(const FullCIWfn &, const double *, double *, double *);

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

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

pybind11::tuple py_compute_rdms_1234_doci(const DOCIWfn &, const Array<double>);

pybind11::tuple py_compute_rdms_fullci(const FullCIWfn &, const Array<double>);

pybind11::tuple py_compute_rdms_genci(const GenCIWfn &, const Array<double>);
Expand Down
40 changes: 40 additions & 0 deletions pyci/src/binding.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1246,6 +1246,46 @@ For Generalized CI wave functions, ``rdm1`` and ``rdm2`` are the full 1-RDM and
)""",
py::arg("wfn"), py::arg("coeffs"));

m.def("compute_rdms_1234", &py_compute_rdms_1234_doci, R"""(
Compute the one-, two-, three-, and four-particle(not implemented yet) reduced density matrices (RDMs) of a wave function.

Parameters
----------
wfn : pyci.wavefunction
Wave function.
coeffs : numpy.ndarray
Coefficient vector.

Returns
-------
For DOCI wave functions, this method returns two nbasis-by-nbasis matrices and two nbasis-by-nbasis-by-nbasis,
which include the unique seniority-zero and seniority-two terms from the full 2-RDMs:

.. math::

D_0 = \left<pp|qq\right>

.. math::

D_2 = \left<pq|pq\right>

.. math::

D_3 = \left<pqr|pqr\right>

.. math::

D_4 = \left<pqq|prr\right>


Notes
-----
This method calculates all RDMs from 1-RDM to 4-RDM (so far implemented only to 3RDM) for a given DOCI wave function.

)""",
py::arg("wfn"), py::arg("coeffs"));


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

m.def("compute_rdms", &py_compute_rdms_genci, py::arg("wfn"), py::arg("coeffs"));
Expand Down
102 changes: 102 additions & 0 deletions pyci/src/rdm.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,95 @@ void compute_rdms(const DOCIWfn &wfn, const double *coeffs, double *d0, double *
}
}


void compute_rdms_1234(const DOCIWfn &wfn, const double *coeffs, double *d0, double *d2, double *d3, double *d4) {
// prepare working vectors
AlignedVector<ulong> v_det(wfn.nword);
AlignedVector<long> v_occs(wfn.nocc_up);
AlignedVector<long> v_virs(wfn.nvir_up);
ulong *det = &v_det[0];
long *occs = &v_occs[0], *virs = &v_virs[0];
// fill rdms with zeros
long i = wfn.nbasis * wfn.nbasis, j = 0;
while (j < i) {
d0[j] = 0;
d2[j++] = 0;
}
i = wfn.nbasis * wfn.nbasis * wfn.nbasis;
j = 0;
while (j < i) {
d3[j]=0;
d4[j++]=0;
}
// iterate over determinants
for (long idet = 0, jdet, mdet, k, l, m, n; idet < wfn.ndet; ++idet) {
double val1, val2;
// fill working vectors
wfn.copy_det(idet, det);
fill_occs(wfn.nword, det, occs);
fill_virs(wfn.nword, wfn.nbasis, det, virs);
// diagonal elements
val1 = coeffs[idet] * coeffs[idet];
for (i = 0; i < wfn.nocc_up; ++i) {
k = occs[i];
d0[k * (wfn.nbasis + 1)] += val1;
for (j = i + 1; j < wfn.nocc_up; ++j) {
l = occs[j];
d2[wfn.nbasis * k + l] += val1;
d2[wfn.nbasis * l + k] += val1;
for (m= j + 1; m < wfn.nocc_up; ++m){
n = occs[m];
d3[(wfn.nbasis * wfn.nbasis) * k + wfn.nbasis * l + n] += val1;
d3[(wfn.nbasis * wfn.nbasis) * k + wfn.nbasis * n + l] += val1;
d3[(wfn.nbasis * wfn.nbasis) * l + wfn.nbasis * k + n] += val1;
d3[(wfn.nbasis * wfn.nbasis) * l + wfn.nbasis * n + k] += val1;
d3[(wfn.nbasis * wfn.nbasis) * n + wfn.nbasis * k + l] += val1;
d3[(wfn.nbasis * wfn.nbasis) * n + wfn.nbasis * l + k] += val1;
}
// pair excitation elements 3rdm j>i
for (m = 0; m < wfn.nvir_up; ++m) {
n = virs[m];
excite_det(l, n, det);
mdet = wfn.index_det(det);
excite_det(n, l, det);
// check if excited determinant is in wfn
if (mdet > idet) {
val2 = coeffs[mdet] * coeffs[idet];
d4[(wfn.nbasis * wfn.nbasis) * k + wfn.nbasis * l + n] += val2;
d4[(wfn.nbasis * wfn.nbasis) * k + wfn.nbasis * n + l] += val2;
}
} //pair excitation elements 3rdm i > j
for (m = 0; m < wfn.nvir_up; ++m) {
n = virs[m];
excite_det(k, n, det);
mdet = wfn.index_det(det);
excite_det(n, k, det);
// check if excited determinant is in wfn
if (mdet > idet) {
val2 = coeffs[mdet] * coeffs[idet];
d4[(wfn.nbasis * wfn.nbasis) * l + wfn.nbasis * k + n] += val2;
d4[(wfn.nbasis * wfn.nbasis) * l + wfn.nbasis * n + k] += val2;
}
}

}
// pair excitation elements
for (j = 0; j < wfn.nvir_up; ++j) {
l = virs[j];
excite_det(k, l, det);
jdet = wfn.index_det(det);
excite_det(l, k, det);
// check if excited determinant is in wfn
if (jdet > idet) {
val2 = coeffs[idet] * coeffs[jdet];
d0[wfn.nbasis * k + l] += val2;
d0[wfn.nbasis * l + k] += val2;
}
}
}
}
}

void compute_rdms(const FullCIWfn &wfn, const double *coeffs, double *rdm1, double *rdm2) {
long n1 = wfn.nbasis;
long n2 = wfn.nbasis * wfn.nbasis;
Expand Down Expand Up @@ -815,6 +904,19 @@ pybind11::tuple py_compute_rdms_doci(const DOCIWfn &wfn, const Array<double> coe
return pybind11::make_tuple(d0, d2);
}

pybind11::tuple py_compute_rdms_1234_doci(const DOCIWfn &wfn, const Array<double> coeffs) {
Array<double> d0({wfn.nbasis, wfn.nbasis});
Array<double> d2({wfn.nbasis, wfn.nbasis});
Array<double> d3({wfn.nbasis, wfn.nbasis, wfn.nbasis});
Array<double> d4({wfn.nbasis, wfn.nbasis, wfn.nbasis});
compute_rdms_1234(wfn, reinterpret_cast<const double *>(coeffs.request().ptr),
reinterpret_cast<double *>(d0.request().ptr),
reinterpret_cast<double *>(d2.request().ptr),
reinterpret_cast<double *>(d3.request().ptr),
reinterpret_cast<double *>(d4.request().ptr));
return pybind11::make_tuple(d0, d2, d3, d4);
}

pybind11::tuple py_compute_rdms_fullci(const FullCIWfn &wfn, const Array<double> coeffs) {
Array<double> rdm1({static_cast<long>(2), wfn.nbasis, wfn.nbasis});
Array<double> rdm2({static_cast<long>(3), wfn.nbasis, wfn.nbasis, wfn.nbasis, wfn.nbasis});
Expand Down
194 changes: 194 additions & 0 deletions pyci/test/data/BH_sto-3g_eq.fcidump
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,194 @@
&FCI NORB= 6,NELEC= 6,MS2=0,
ORBSYM=1,1,1,1,1,1,
ISYM=1,
&END
2.890694185551886 1 1 1 1
-0.2285478644008276 1 1 2 1
0.5953083364634656 1 1 2 2
-0.2162016548582802 1 1 3 1
0.09104871057910127 1 1 3 2
0.6900652953132707 1 1 3 3
0.7419889025330183 1 1 4 4
0.7419889025330181 1 1 5 5
0.1413257075902383 1 1 6 1
-0.1861902551212619 1 1 6 2
0.09008046930667073 1 1 6 3
0.6699597694264849 1 1 6 6
-0.2285478644008276 2 1 1 1
0.02992401570960048 2 1 2 1
-0.001146336251477235 2 1 2 2
0.02347037645652453 2 1 3 1
-0.008237570230799412 2 1 3 2
-0.01537296998854835 2 1 3 3
-0.008814449236465278 2 1 4 4
-0.008814449236465275 2 1 5 5
-0.02314657520492044 2 1 6 1
0.003050017645768165 2 1 6 2
-0.00419211917379873 2 1 6 3
-0.004379253657977653 2 1 6 6
0.5953083364634657 2 2 1 1
-0.001146336251477227 2 2 2 1
0.5804858409021905 2 2 2 2
-0.01900821896163675 2 2 3 1
-0.07241056645830383 2 2 3 2
0.4170856180274074 2 2 3 3
0.4593408467391383 2 2 4 4
0.4593408467391382 2 2 5 5
-0.01238748160972113 2 2 6 1
0.0002084163025978247 2 2 6 2
-0.05254431697407079 2 2 6 3
0.5754116421934867 2 2 6 6
-0.2162016548582803 3 1 1 1
0.02347037645652455 3 1 2 1
-0.01900821896163674 3 1 2 2
0.03397697977407958 3 1 3 1
0.003163633505588129 3 1 3 2
0.00404893468739087 3 1 3 3
-0.007412212106220699 3 1 4 4
-0.007412212106220697 3 1 5 5
-0.005247409917599474 3 1 6 1
0.01158307545322849 3 1 6 2
0.006990498122004285 3 1 6 3
-0.01879243966504147 3 1 6 6
0.09104871057910127 3 2 1 1
-0.008237570230799426 3 2 2 1
-0.07241056645830381 3 2 2 2
0.003163633505588093 3 2 3 1
0.0676156185213167 3 2 3 2
0.05019628841726888 3 2 3 3
0.04621563081453394 3 2 4 4
0.04621563081453393 3 2 5 5
0.01250097777451356 3 2 6 1
-0.05223141555214424 3 2 6 2
0.05747267715242108 3 2 6 3
-0.06327529964411363 3 2 6 6
0.6900652953132709 3 3 1 1
-0.01537296998854832 3 3 2 1
0.4170856180274075 3 3 2 2
0.004048934687390883 3 3 3 1
0.05019628841726886 3 3 3 2
0.5635811972005504 3 3 3 3
0.5041862576414032 3 3 4 4
0.5041862576414031 3 3 5 5
0.0231276130620684 3 3 6 1
-0.03192699042073559 3 3 6 2
0.07541684244416474 3 3 6 3
0.4565326809201556 3 3 6 6
0.02169507828733871 4 1 4 1
0.01921155457995589 4 1 4 2
0.01851865219375422 4 1 4 3
-0.01161252859365546 4 1 6 4
0.01921155457995589 4 2 4 1
0.06039821130266281 4 2 4 2
0.04342107961102947 4 2 4 3
-0.03633847993724627 4 2 6 4
0.01851865219375421 4 3 4 1
0.04342107961102945 4 3 4 2
0.06291076116021634 4 3 4 3
-0.01322878606419047 4 3 6 4
0.7419889025330183 4 4 1 1
-0.008814449236465242 4 4 2 1
0.4593408467391384 4 4 2 2
-0.00741221210622064 4 4 3 1
0.04621563081453399 4 4 3 2
0.5041862576414033 4 4 3 3
0.5867727439789322 4 4 4 4
0.5235134814676911 4 4 5 5
0.004168026575771488 4 4 6 1
-0.09238076313811244 4 4 6 2
0.04177944228679972 4 4 6 3
0.4740773398260829 4 4 6 6
0.0216950782873387 5 1 5 1
0.01921155457995588 5 1 5 2
0.01851865219375421 5 1 5 3
-0.01161252859365545 5 1 6 5
0.01921155457995588 5 2 5 1
0.06039821130266281 5 2 5 2
0.04342107961102946 5 2 5 3
-0.03633847993724625 5 2 6 5
0.01851865219375421 5 3 5 1
0.04342107961102945 5 3 5 2
0.06291076116021634 5 3 5 3
-0.01322878606419046 5 3 6 5
0.03162963125562031 5 4 5 4
0.7419889025330181 5 5 1 1
-0.008814449236465207 5 5 2 1
0.4593408467391384 5 5 2 2
-0.007412212106220606 5 5 3 1
0.04621563081453402 5 5 3 2
0.5041862576414033 5 5 3 3
0.5235134814676911 5 5 4 4
0.586772743978932 5 5 5 5
0.004168026575771471 5 5 6 1
-0.09238076313811253 5 5 6 2
0.04177944228679967 5 5 6 3
0.4740773398260828 5 5 6 6
0.1413257075902383 6 1 1 1
-0.02314657520492046 6 1 2 1
-0.01238748160972115 6 1 2 2
-0.005247409917599475 6 1 3 1
0.01250097777451353 6 1 3 2
0.02312761306206843 6 1 3 3
0.004168026575771572 6 1 4 4
0.00416802657577157 6 1 5 5
0.02880391396937533 6 1 6 1
0.005524338027076478 6 1 6 2
0.01082364008038884 6 1 6 3
-0.008812192154253081 6 1 6 6
-0.186190255121262 6 2 1 1
0.003050017645768176 6 2 2 1
0.0002084163025976651 6 2 2 2
0.01158307545322851 6 2 3 1
-0.05223141555214433 6 2 3 2
-0.03192699042073567 6 2 3 3
-0.0923807631381124 6 2 4 4
-0.09238076313811237 6 2 5 5
0.005524338027076496 6 2 6 1
0.1212730361895562 6 2 6 2
-0.05403162334477736 6 2 6 3
0.01293203355693659 6 2 6 6
0.09008046930667071 6 3 1 1
-0.004192119173798695 6 3 2 1
-0.05254431697407061 6 3 2 2
0.006990498122004348 6 3 3 1
0.05747267715242119 6 3 3 2
0.07541684244416491 6 3 3 3
0.0417794422867997 6 3 4 4
0.04177944228679968 6 3 5 5
0.01082364008038878 6 3 6 1
-0.05403162334477741 6 3 6 2
0.07634911021780025 6 3 6 3
-0.04207173170262291 6 3 6 6
-0.01161252859365545 6 4 4 1
-0.03633847993724626 6 4 4 2
-0.01322878606419046 6 4 4 3
0.03249943702929813 6 4 6 4
-0.01161252859365545 6 5 5 1
-0.03633847993724625 6 5 5 2
-0.01322878606419045 6 5 5 3
0.03249943702929813 6 5 6 5
0.6699597694264852 6 6 1 1
-0.004379253657977538 6 6 2 1
0.5754116421934872 6 6 2 2
-0.01879243966504134 6 6 3 1
-0.06327529964411353 6 6 3 2
0.4565326809201556 6 6 3 3
0.4740773398260831 6 6 4 4
0.474077339826083 6 6 5 5
-0.008812192154253214 6 6 6 1
0.01293203355693684 6 6 6 2
-0.04207173170262306 6 6 6 3
0.6101411445973578 6 6 6 6
-12.73694583468452 1 1 0 0
0.2636037738346674 2 1 0 0
-3.08122009866082 2 2 0 0
0.2419315877376549 3 1 0 0
-0.1364127650686597 3 2 0 0
-2.922828142214635 3 3 0 0
-2.996085190675585 4 4 0 0
-2.996085190675584 5 5 0 0
-0.1527654545391036 6 1 0 0
0.4703521766204644 6 2 0 0
-0.2079679737989178 6 3 0 0
-2.475899716280425 6 6 0 0
2.147634784577923 0 0 0 0
Loading
Loading