Skip to content

4RDM implementation #105

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
merged 32 commits into from
Apr 10, 2025
Merged

4RDM implementation #105

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
32 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
b137aae
First attempt to include 4RDM elements in compute_rdms_1234
DanielCalero1 Dec 10, 2024
b97f92a
Chaing python biding to receive 4RDM vectors
DanielCalero1 Feb 27, 2025
4707357
Adding H8 electron integrals to test 4RDM
DanielCalero1 Mar 25, 2025
acb3b97
Adding new test file(delete later)
DanielCalero1 Mar 25, 2025
1ee1e25
Adding computation of general 4RDM with the appropiate symmetries
DanielCalero1 Mar 25, 2025
80df38e
Adding new test_routines for the 4RDM
DanielCalero1 Mar 25, 2025
490ef93
Modifying spin-free RDMs function
DanielCalero1 Mar 25, 2025
0441aaf
Modyfing _init_ function to include spinize_RDM_1234 and spin_free_rdms
DanielCalero1 Mar 25, 2025
655bf73
Updating the test rotine for the 4RDM
DanielCalero1 Mar 25, 2025
e74e342
Deleting auxiliary test file
DanielCalero1 Mar 25, 2025
560e785
Fixing typo in d7 definition
DanielCalero1 Apr 2, 2025
ca23e4d
Fixing mixed-spin block definitions
DanielCalero1 Apr 2, 2025
32b92a7
Adding option to select if compute 4rdm or not
DanielCalero1 Apr 2, 2025
e4caf7a
Fixing typo in test routine
DanielCalero1 Apr 2, 2025
4960628
Merge branch 'master' into High-order-RDM
DanielCalero1 Apr 7, 2025
24a1fb0
setting up bigmem mark
DanielCalero1 Apr 8, 2025
87a3ba9
Updating Makefile to ignore 34 RDM test unless specified
DanielCalero1 Apr 10, 2025
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: 4 additions & 0 deletions Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,10 @@ all: pyci/_pyci.so.$(PYCI_VERSION) pyci/_pyci.so.$(VERSION_MAJOR) pyci/_pyci.so
test:
@set -e; $(PYTHON) -m pytest -sv ./pyci

.PHONY: test-bigmem
test-bigmem:
@set -e; $(PYTHON) -m pytest -sv ./pyci --bigmem

.PHONY: clean
clean:
rm -rf pyci/src/*.o pyci/_pyci.so*
Expand Down
2 changes: 2 additions & 0 deletions pyci/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,8 @@
"make_senzero_integrals",
"reduce_senzero_integrals",
"spinize_rdms",
"spinize_rdms_1234",
"spin_free_rdms",
"odometer_one_spin",
"odometer_two_spin",
"add_excitations",
Expand Down
2 changes: 1 addition & 1 deletion pyci/include/pyci.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -221,7 +221,7 @@ 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_1234(const DOCIWfn &, const double *, double *, double *, double *, double *, double *, double *, double *);

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

Expand Down
133 changes: 126 additions & 7 deletions pyci/src/rdm.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ 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) {
void compute_rdms_1234(const DOCIWfn &wfn, const double *coeffs, double *d0, double *d2, double *d3, double *d4, double *d5, double *d6, double *d7) {
// prepare working vectors
AlignedVector<ulong> v_det(wfn.nword);
AlignedVector<long> v_occs(wfn.nocc_up);
Expand All @@ -84,8 +84,15 @@ void compute_rdms_1234(const DOCIWfn &wfn, const double *coeffs, double *d0, dou
d3[j]=0;
d4[j++]=0;
}
i = wfn.nbasis * wfn.nbasis * wfn.nbasis * wfn.nbasis;
j = 0;
while (j < i) {
d5[j]=0;
d6[j]=0;
d7[j++]=0;
}
// iterate over determinants
for (long idet = 0, jdet, mdet, k, l, m, n; idet < wfn.ndet; ++idet) {
for (long idet = 0, jdet, mdet, pdet, k, l, m, n, p, q; idet < wfn.ndet; ++idet) {
double val1, val2;
// fill working vectors
wfn.copy_det(idet, det);
Expand All @@ -107,9 +114,86 @@ void compute_rdms_1234(const DOCIWfn &wfn, const double *coeffs, double *d0, dou
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;
d3[(wfn.nbasis * wfn.nbasis) * n + wfn.nbasis * l + k] += val1;
for (p= m + 1; p < wfn.nocc_up; ++p){
q= occs[p];
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * n + q] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * q + n] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * l + q] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * q + l] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * q + (wfn.nbasis) * l + n] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * q + (wfn.nbasis) * n + l] += val1;

d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * n + q] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * q + n] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * k + q] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * q + k] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * q + (wfn.nbasis) * k + n] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * q + (wfn.nbasis) * n + k] += val1;

d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * k + q] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * q + k] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * l + q] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * q + l] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * q + (wfn.nbasis) * l + k] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * q + (wfn.nbasis) * k + l] += val1;

d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * q + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * n + k] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * q + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * k + n] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * q + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * l + k] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * q + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * k + l] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * q + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * l + n] += val1;
d5[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * q + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * n + l] += val1;

}
//pair excitation elements 4rdm n >k,l
for (p= 0; p < wfn.nvir_up; ++p){
q = virs[p];
excite_det(n, q, det);
pdet = wfn.index_det(det);
excite_det(q, n, det);
// check if excited determinant is in wfn
if (pdet > idet) {
val2 = coeffs[pdet] * coeffs[idet];
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * n + q] += val2;
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * q + n] += val2;
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * n + q] += val2;
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * q + n] += val2;
}
}
//pair excitation elements 4rdm l >k,n
for (p= 0; p < wfn.nvir_up; ++p){
q = virs[p];
excite_det(l, q, det);
pdet = wfn.index_det(det);
excite_det(q, l, det);
// check if excited determinant is in wfn
if (pdet > idet) {
val2 = coeffs[pdet] * coeffs[idet];
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * l + q] += val2;
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * q + l] += val2;
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * l + q] += val2;
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * q + l] += val2;
}
}
//pair excitation elements 4rdm k >n,l
for (p= 0; p < wfn.nvir_up; ++p){
q = virs[p];
excite_det(k, q, det);
pdet = wfn.index_det(det);
excite_det(q, k, det);
// check if excited determinant is in wfn
if (pdet > idet) {
val2 = coeffs[pdet] * coeffs[idet];
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * k + q] += val2;
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * q + k] += val2;
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * k + q] += val2;
d6[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * q + k] += val2;
}
}

}
// pair excitation elements 3rdm j>i
// pair excitation elements 3rdm l>k
for (m = 0; m < wfn.nvir_up; ++m) {
n = virs[m];
excite_det(l, n, det);
Expand All @@ -121,7 +205,7 @@ void compute_rdms_1234(const DOCIWfn &wfn, const double *coeffs, double *d0, dou
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
} //pair excitation elements 3rdm k > l
for (m = 0; m < wfn.nvir_up; ++m) {
n = virs[m];
excite_det(k, n, det);
Expand All @@ -134,6 +218,35 @@ void compute_rdms_1234(const DOCIWfn &wfn, const double *coeffs, double *d0, dou
d4[(wfn.nbasis * wfn.nbasis) * l + wfn.nbasis * n + k] += val2;
}
}
//double-pair excitation elements 4rdm
for (m = 0; m < wfn.nvir_up; ++m) {
n = virs[m];
for (p= m + 1; p < wfn.nvir_up; ++p){
q = virs[p];

excite_det(l, q, det);
pdet = wfn.index_det(det);
excite_det(q, l, det);

excite_det(k, n, det);
mdet = wfn.index_det(det);
excite_det(n, k, det);
// check if excited determinants are the same and if it is in wfn
if ((mdet==pdet) && (mdet> idet)) {
val2 = coeffs[mdet] * coeffs[idet];
d7[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * n + q] += val2;
d7[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * k + (wfn.nbasis * wfn.nbasis) * l + (wfn.nbasis) * q + n] += val2;
d7[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * q + n] += val2;
d7[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * l + (wfn.nbasis * wfn.nbasis) * k + (wfn.nbasis) * n + q] += val2;
d7[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * q + (wfn.nbasis) * k + l] += val2;
d7[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * n + (wfn.nbasis * wfn.nbasis) * q + (wfn.nbasis) * l + k] += val2;
d7[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * q + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * l + k] += val2;
d7[(wfn.nbasis * wfn.nbasis * wfn.nbasis ) * q + (wfn.nbasis * wfn.nbasis) * n + (wfn.nbasis) * k + l] += val2;


}
}
}

}
// pair excitation elements
Expand Down Expand Up @@ -909,12 +1022,18 @@ pybind11::tuple py_compute_rdms_1234_doci(const DOCIWfn &wfn, const Array<double
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});
Array<double> d5({wfn.nbasis, wfn.nbasis, wfn.nbasis, wfn.nbasis});
Array<double> d6({wfn.nbasis, wfn.nbasis, wfn.nbasis, wfn.nbasis});
Array<double> d7({wfn.nbasis, 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);
reinterpret_cast<double *>(d4.request().ptr),
reinterpret_cast<double *>(d5.request().ptr),
reinterpret_cast<double *>(d6.request().ptr),
reinterpret_cast<double *>(d7.request().ptr));
return pybind11::make_tuple(d0, d2, d3, d4, d5, d6, d7);
}

pybind11::tuple py_compute_rdms_fullci(const FullCIWfn &wfn, const Array<double> coeffs) {
Expand Down
22 changes: 22 additions & 0 deletions pyci/test/conftest.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
import pytest

def pytest_addoption(parser):
parser.addoption(
"--bigmem",
action="store_true",
default=False,
help="run tests that require a large amount of memory"
)

def pytest_configure(config):
config.addinivalue_line(
"markers",
"bigmem: mark test as requiring a large amount of memory"
)

def pytest_collection_modifyitems(config, items):
if not config.getoption("--bigmem"):
skip_bigmem = pytest.mark.skip(reason="need --bigmem option to run")
for item in items:
if "bigmem" in item.keywords:
item.add_marker(skip_bigmem)
Loading