From 4fc41efc6c843501b036f60239ea14b6b4aad4c7 Mon Sep 17 00:00:00 2001 From: "Christian B. Mendl" Date: Tue, 3 Sep 2024 23:22:16 +0200 Subject: [PATCH] molecular Hamiltonian MPO construction for a spin orbital basis available via Python interface --- examples/chemistry/water_molecule.ipynb | 160 ++++++++------------ examples/python/chemtensor_python.ipynb | 4 +- pymodule/pymodule.c | 186 +++++++++++++++++++----- 3 files changed, 211 insertions(+), 139 deletions(-) diff --git a/examples/chemistry/water_molecule.ipynb b/examples/chemistry/water_molecule.ipynb index 9f9bf7c..df205ff 100644 --- a/examples/chemistry/water_molecule.ipynb +++ b/examples/chemistry/water_molecule.ipynb @@ -28,7 +28,7 @@ "sys.path.append(\"../../build/\")\n", "import chemtensor\n", "\n", - "# pyscf (https://pyscf.org/) can define a molecular basis, compute overlap integrals and run other computational methods for comparison.\n", + "# pyscf (https://pyscf.org/) defines a molecular orbital basis, computes overlap integrals and runs other computational methods for comparison.\n", "import pyscf" ] }, @@ -62,7 +62,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "converged SCF energy = -74.9307084821\n" + "converged SCF energy = -74.9307084820999\n" ] } ], @@ -80,13 +80,13 @@ "name": "stdout", "output_type": "stream", "text": [ - "E(CCSD) = -74.97016403895398 E_corr = -0.03945555685402184\n" + "E(CCSD) = -74.97016403895393 E_corr = -0.03945555685402184\n" ] }, { "data": { "text/plain": [ - "-74.97016403895398" + "-74.97016403895393" ] }, "execution_count": 4, @@ -95,7 +95,6 @@ } ], "source": [ - "\n", "# run coupled-cluster with single and double excitations (CCSD), for comparison\n", "ccsd = pyscf.cc.CCSD(hf).run()\n", "ccsd.e_tot" @@ -135,40 +134,6 @@ "cell_type": "code", "execution_count": 6, "metadata": {}, - "outputs": [], - "source": [ - "def spatial_to_spin_overlap_integrals(h1, h2):\n", - " \"\"\"\n", - " Enlarge the single- and two-particle electron overlap integral tensors\n", - " from an orbital basis without spin to a spin-orbital basis.\n", - " \"\"\"\n", - " h1 = np.asarray(h1)\n", - " h2 = np.asarray(h2)\n", - "\n", - " n = h1.shape[0]\n", - " assert h1.shape == (n, n)\n", - " assert h2.shape == (n, n, n, n)\n", - "\n", - " # single-particle integrals\n", - " h1_so = np.kron(np.eye(2), h1)\n", - "\n", - " # two-particle integrals\n", - " tmp = np.zeros((2*n, 2*n, n, n))\n", - " for i in range(n):\n", - " for j in range(n):\n", - " tmp[:, :, i, j] = np.kron(np.eye(2), h2[:,:, i, j])\n", - " h2_so = np.zeros((2*n, 2*n, 2*n, 2*n))\n", - " for i in range(2*n):\n", - " for j in range(2*n):\n", - " h2_so[i, j, :, :] = np.kron(np.eye(2), tmp[i, j, :, :])\n", - "\n", - " return h1_so, h2_so" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, "outputs": [ { "name": "stdout", @@ -186,7 +151,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -208,34 +173,13 @@ }, { "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "(14, 14)\n", - "(14, 14, 14, 14)\n" - ] - } - ], - "source": [ - "# extend to spin-orbitals\n", - "h1_so, eri_so = spatial_to_spin_overlap_integrals(h1_mo, eri_mo)\n", - "print(h1_so.shape)\n", - "print(eri_so.shape)" - ] - }, - { - "cell_type": "code", - "execution_count": 10, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# convert to physicists' convention\n", - "tkin = h1_so\n", - "vint = np.transpose(eri_so, (0, 2, 1, 3))" + "tkin = h1_mo\n", + "vint = np.transpose(eri_mo, (0, 2, 1, 3))" ] }, { @@ -247,44 +191,67 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[1, 4, 16, 31, 38, 43, 52, 51, 52, 43, 38, 31, 16, 4, 1]" + "[1, 36, 58, 96, 96, 58, 36, 1]" ] }, - "execution_count": 11, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "hamiltonian = chemtensor.construct_molecular_hamiltonian_mpo(tkin, vint, optimize=True)\n", + "hamiltonian = chemtensor.construct_spin_molecular_hamiltonian_mpo(tkin, vint)\n", "# virtual bond dimensions\n", "hamiltonian.bond_dims" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[(0, 0), (1, -1), (1, 1), (2, 0)]" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# local physical quantum numbers (number of electrons and spin)\n", + "[chemtensor.decode_quantum_number_pair(qnum) for qnum in hamiltonian.qsite]" + ] + }, + { + "cell_type": "code", + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ - "# number of electrons (determines quantum number sector)\n", - "pnum = 10" + "# overall quantum number sector of quantum state (number of electrons and spin)\n", + "q_pnum = 10\n", + "q_spin = 0\n", + "qnum_sector = chemtensor.encode_quantum_number_pair(q_pnum, q_spin)" ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# run two-site DMRG\n", - "psi, en_sweeps, entropy = chemtensor.dmrg(hamiltonian, num_sweeps=6, maxiter_lanczos=25, tol_split=1e-9, qnum_sector=pnum)" + "psi, en_sweeps, entropy = chemtensor.dmrg(hamiltonian, num_sweeps=6, maxiter_lanczos=25, tol_split=1e-9, qnum_sector=qnum_sector)" ] }, { @@ -296,16 +263,16 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[1, 2, 4, 7, 10, 14, 19, 19, 20, 16, 11, 7, 4, 2, 1]" + "[1, 4, 13, 27, 19, 16, 4, 1]" ] }, - "execution_count": 14, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } @@ -317,21 +284,21 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[-84.88903443027432,\n", - " -84.8890344313337,\n", - " -84.88903443132047,\n", - " -84.88903443131844,\n", - " -84.88903443131815,\n", - " -84.88903443131794]" + "[-84.88903443537285,\n", + " -84.88903443537973,\n", + " -84.88903443537977,\n", + " -84.88903443537987,\n", + " -84.88903443537993,\n", + " -84.8890344353799]" ] }, - "execution_count": 15, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -343,16 +310,16 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "-74.97027263097696" + "-74.97027263503892" ] }, - "execution_count": 16, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -365,16 +332,16 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "0.00010859202298263426" + "0.0001085960849991352" ] }, - "execution_count": 17, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -386,16 +353,16 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "5.8616635101316206e-08" + "5.455467544379644e-08" ] }, - "execution_count": 18, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -404,13 +371,6 @@ "# difference to FCI energy\n", "en_dmrg - en_fci" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/examples/python/chemtensor_python.ipynb b/examples/python/chemtensor_python.ipynb index 0340b13..9f49b84 100644 --- a/examples/python/chemtensor_python.ipynb +++ b/examples/python/chemtensor_python.ipynb @@ -114,7 +114,7 @@ ], "source": [ "# local physical quantum numbers (particle number and spin)\n", - "[chemtensor.fermi_hubbard_decode_quantum_numbers(qnum) for qnum in hamiltonian.qsite]" + "[chemtensor.decode_quantum_number_pair(qnum) for qnum in hamiltonian.qsite]" ] }, { @@ -154,7 +154,7 @@ "# overall quantum number sector of quantum state (particle number and spin)\n", "q_pnum = 7\n", "q_spin = 1\n", - "qnum_sector = chemtensor.fermi_hubbard_encode_quantum_numbers(q_pnum, q_spin)" + "qnum_sector = chemtensor.encode_quantum_number_pair(q_pnum, q_spin)" ] }, { diff --git a/pymodule/pymodule.c b/pymodule/pymodule.c index 7452576..dae9d09 100644 --- a/pymodule/pymodule.c +++ b/pymodule/pymodule.c @@ -671,6 +671,41 @@ static PyTypeObject PyMPOType = { // +static PyObject* Py_encode_quantum_number_pair(PyObject* Py_UNUSED(self), PyObject* args) +{ + // individual quantum numbers + qnumber qa, qb; + + // parse input arguments + if (!PyArg_ParseTuple(args, "ii", &qa, &qb)) { + PyErr_SetString(PyExc_SyntaxError, "error parsing input; syntax: encode_quantum_number_pair(qa: int, qb: int)"); + return NULL; + } + + qnumber qnum = encode_quantum_number_pair(qa, qb); + + return PyLong_FromLong(qnum); +} + + +static PyObject* Py_decode_quantum_number_pair(PyObject* Py_UNUSED(self), PyObject* args) +{ + qnumber qnum; + + // parse input arguments + if (!PyArg_ParseTuple(args, "i", &qnum)) { + PyErr_SetString(PyExc_SyntaxError, "error parsing input; syntax: decode_quantum_number_pair(qnum: int)"); + return NULL; + } + + // decode quantum numbers + qnumber qa, qb; + decode_quantum_number_pair(qnum, &qa, &qb); + + return PyTuple_Pack(2, PyLong_FromLong(qa), PyLong_FromLong(qb)); +} + + static PyObject* Py_construct_random_mps(PyObject* Py_UNUSED(self), PyObject* args, PyObject* kwargs) { // data type @@ -922,45 +957,116 @@ static PyObject* Py_construct_fermi_hubbard_1d_mpo(PyObject* Py_UNUSED(self), Py } -static PyObject* Py_fermi_hubbard_encode_quantum_numbers(PyObject* Py_UNUSED(self), PyObject* args) +static PyObject* Py_construct_molecular_hamiltonian_mpo(PyObject* Py_UNUSED(self), PyObject* args, PyObject* kwargs) { - // particle number quantum number - qnumber q_pnum; - // spin quantum number - qnumber q_spin; + PyObject* py_obj_tkin; + PyObject* py_obj_vint; + int optimize = 0; // parse input arguments - if (!PyArg_ParseTuple(args, "ii", &q_pnum, &q_spin)) { - PyErr_SetString(PyExc_SyntaxError, "error parsing input; syntax: fermi_hubbard_encode_quantum_numbers(q_pnum: int, q_spin: int)"); + char* kwlist[] = { "", "", "optimize", NULL }; + if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|p", kwlist, &py_obj_tkin, &py_obj_vint, &optimize)) { + PyErr_SetString(PyExc_SyntaxError, "error parsing input; syntax: construct_molecular_hamiltonian_mpo(tkin, vint, optimize=False)"); return NULL; } - // encode particle and spin quantum numbers - qnumber q = encode_quantum_number_pair(q_pnum, q_spin); - - return PyLong_FromLong(q); -} + PyArrayObject* py_tkin = (PyArrayObject*)PyArray_ContiguousFromObject(py_obj_tkin, NPY_DOUBLE, 2, 2); + if (py_tkin == NULL) { + PyErr_SetString(PyExc_ValueError, "converting input argument 'tkin' to a NumPy array with degree 2 failed"); + return NULL; + } + PyArrayObject* py_vint = (PyArrayObject*)PyArray_ContiguousFromObject(py_obj_vint, NPY_DOUBLE, 4, 4); + if (py_vint == NULL) { + PyErr_SetString(PyExc_ValueError, "converting input argument 'vint' to a NumPy array with degree 4 failed"); + return NULL; + } + // argument checks + if (PyArray_NDIM(py_tkin) != 2) { + char msg[1024]; + sprintf(msg, "'tkin' must have degree 2, received %i", PyArray_NDIM(py_tkin)); + PyErr_SetString(PyExc_ValueError, msg); + return NULL; + } + if (PyArray_DIM(py_tkin, 0) != PyArray_DIM(py_tkin, 1)) { + char msg[1024]; + sprintf(msg, "'tkin' must be a square matrix, received a %li x %li matrix", PyArray_DIM(py_tkin, 0), PyArray_DIM(py_tkin, 1)); + PyErr_SetString(PyExc_ValueError, msg); + return NULL; + } + if (PyArray_DIM(py_tkin, 0) == 0) { + PyErr_SetString(PyExc_ValueError, "'tkin' cannot be an empty matrix"); + return NULL; + } + if (PyArray_TYPE(py_tkin) != NPY_DOUBLE) { + PyErr_SetString(PyExc_TypeError, "'tkin' must have 'double' format entries"); + return NULL; + } + if (!(PyArray_FLAGS(py_tkin) & NPY_ARRAY_C_CONTIGUOUS)) { + PyErr_SetString(PyExc_SyntaxError, "'tkin' does not have contiguous C storage format"); + return NULL; + } + if (PyArray_NDIM(py_vint) != 4) { + char msg[1024]; + sprintf(msg, "'vint' must have degree 4, received %i", PyArray_NDIM(py_vint)); + PyErr_SetString(PyExc_ValueError, msg); + return NULL; + } + if ((PyArray_DIM(py_vint, 0) != PyArray_DIM(py_vint, 1)) || + (PyArray_DIM(py_vint, 0) != PyArray_DIM(py_vint, 2)) || + (PyArray_DIM(py_vint, 0) != PyArray_DIM(py_vint, 3))) { + char msg[1024]; + sprintf(msg, "all dimensions of 'vint' must agree, received a %li x %li x %li x %li tensor", PyArray_DIM(py_vint, 0), PyArray_DIM(py_vint, 1), PyArray_DIM(py_vint, 2), PyArray_DIM(py_vint, 3)); + PyErr_SetString(PyExc_ValueError, msg); + return NULL; + } + if (PyArray_TYPE(py_vint) != NPY_DOUBLE) { + PyErr_SetString(PyExc_TypeError, "'vint' must have 'double' format entries"); + return NULL; + } + if (!(PyArray_FLAGS(py_vint) & NPY_ARRAY_C_CONTIGUOUS)) { + PyErr_SetString(PyExc_SyntaxError, "'vint' does not have contiguous C storage format"); + return NULL; + } + if (PyArray_DIM(py_tkin, 0) != PyArray_DIM(py_vint, 0)) { + PyErr_SetString(PyExc_SyntaxError, "'tkin' and 'vint' must have the same axis dimensions"); + return NULL; + } -static PyObject* Py_fermi_hubbard_decode_quantum_numbers(PyObject* Py_UNUSED(self), PyObject* args) -{ - qnumber qnum; + const long nsites = (long)PyArray_DIM(py_tkin, 0); + long dim_tkin[2] = { nsites, nsites }; + struct dense_tensor tkin = { + .data = PyArray_DATA(py_tkin), + .dim = dim_tkin, + .dtype = CT_DOUBLE_REAL, + .ndim = 2, + }; + long dim_vint[4] = { nsites, nsites, nsites, nsites }; + struct dense_tensor vint = { + .data = PyArray_DATA(py_vint), + .dim = dim_vint, + .dtype = CT_DOUBLE_REAL, + .ndim = 4, + }; - // parse input arguments - if (!PyArg_ParseTuple(args, "i", &qnum)) { - PyErr_SetString(PyExc_SyntaxError, "error parsing input; syntax: fermi_hubbard_decode_quantum_numbers(qnum: int)"); + PyMPOObject* py_mpo = (PyMPOObject*)PyMPO_new(&PyMPOType, NULL, NULL); + if (py_mpo == NULL) { + PyErr_SetString(PyExc_RuntimeError, "error creating PyMPO object"); return NULL; } - // decode quantum numbers - qnumber q_pnum, q_spin; - decode_quantum_number_pair(qnum, &q_pnum, &q_spin); + // actually construct the assembly and MPO + construct_molecular_hamiltonian_mpo_assembly(&tkin, &vint, (bool)optimize, &py_mpo->assembly); + mpo_from_assembly(&py_mpo->assembly, &py_mpo->mpo); - return PyTuple_Pack(2, PyLong_FromLong(q_pnum), PyLong_FromLong(q_spin)); + Py_DECREF(py_tkin); + Py_DECREF(py_vint); + + return (PyObject*)py_mpo; } -static PyObject* Py_construct_molecular_hamiltonian_mpo(PyObject* Py_UNUSED(self), PyObject* args, PyObject* kwargs) +static PyObject* Py_construct_spin_molecular_hamiltonian_mpo(PyObject* Py_UNUSED(self), PyObject* args, PyObject* kwargs) { PyObject* py_obj_tkin; PyObject* py_obj_vint; @@ -969,7 +1075,7 @@ static PyObject* Py_construct_molecular_hamiltonian_mpo(PyObject* Py_UNUSED(self // parse input arguments char* kwlist[] = { "", "", "optimize", NULL }; if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|p", kwlist, &py_obj_tkin, &py_obj_vint, &optimize)) { - PyErr_SetString(PyExc_SyntaxError, "error parsing input; syntax: construct_molecular_hamiltonian_mpo(tkin, vint, optimize=False)"); + PyErr_SetString(PyExc_SyntaxError, "error parsing input; syntax: construct_spin_molecular_hamiltonian_mpo(tkin, vint, optimize=False)"); return NULL; } @@ -1059,7 +1165,7 @@ static PyObject* Py_construct_molecular_hamiltonian_mpo(PyObject* Py_UNUSED(self } // actually construct the assembly and MPO - construct_molecular_hamiltonian_mpo_assembly(&tkin, &vint, (bool)optimize, &py_mpo->assembly); + construct_spin_molecular_hamiltonian_mpo_assembly(&tkin, &vint, (bool)optimize, &py_mpo->assembly); mpo_from_assembly(&py_mpo->assembly, &py_mpo->mpo); Py_DECREF(py_tkin); @@ -1284,6 +1390,18 @@ static PyObject* Py_operator_average_coefficient_gradient(PyObject* Py_UNUSED(se static PyMethodDef methods[] = { + { + .ml_name = "encode_quantum_number_pair", + .ml_meth = Py_encode_quantum_number_pair, + .ml_flags = METH_VARARGS, + .ml_doc = "Encode a pair of quantum numbers into a single quantum number.\nSyntax: encode_quantum_number_pair(qa: int, qb: int)", + }, + { + .ml_name = "decode_quantum_number_pair", + .ml_meth = Py_decode_quantum_number_pair, + .ml_flags = METH_VARARGS, + .ml_doc = "Decode a quantum number into two separate quantum numbers.\nSyntax: decode_quantum_number_pair(qnum: int)", + }, { .ml_name = "construct_random_mps", .ml_meth = (PyCFunction)Py_construct_random_mps, @@ -1314,24 +1432,18 @@ static PyMethodDef methods[] = { .ml_flags = METH_VARARGS, .ml_doc = "Construct an MPO representation of the Fermi-Hubbard Hamiltonian with nearest-neighbor hopping on a one-dimensional lattice.\nSyntax: construct_fermi_hubbard_1d_mpo(nsites: int, t: float, u: float, mu: float)", }, - { - .ml_name = "fermi_hubbard_encode_quantum_numbers", - .ml_meth = Py_fermi_hubbard_encode_quantum_numbers, - .ml_flags = METH_VARARGS, - .ml_doc = "Encode a particle number and spin quantum number for the Fermi-Hubbard Hamiltonian into a single quantum number.\nSyntax: fermi_hubbard_encode_quantum_numbers(q_pnum: int, q_spin: int)", - }, - { - .ml_name = "fermi_hubbard_decode_quantum_numbers", - .ml_meth = Py_fermi_hubbard_decode_quantum_numbers, - .ml_flags = METH_VARARGS, - .ml_doc = "Decode a quantum number of the Fermi-Hubbard Hamiltonian into separate particle number and spin quantum numbers.\nSyntax: fermi_hubbard_decode_quantum_numbers(qnum: int)", - }, { .ml_name = "construct_molecular_hamiltonian_mpo", .ml_meth = (PyCFunction)Py_construct_molecular_hamiltonian_mpo, .ml_flags = METH_VARARGS | METH_KEYWORDS, .ml_doc = "Construct a molecular Hamiltonian as MPO, using physicists' convention for the interaction term:\nH = sum_{i,j} t_{i,j} ad_i a_j + 1/2 sum_{i,j,k,l} v_{i,j,k,l} ad_i ad_j a_l a_k\nSyntax: construct_molecular_hamiltonian_mpo(tkin, vint, optimize=False)", }, + { + .ml_name = "construct_spin_molecular_hamiltonian_mpo", + .ml_meth = (PyCFunction)Py_construct_spin_molecular_hamiltonian_mpo, + .ml_flags = METH_VARARGS | METH_KEYWORDS, + .ml_doc = "Construct a molecular Hamiltonian as MPO assuming a spin orbital basis, using physicists' convention for the interaction term:\nH = sum_{i,j,sigma} t_{i,j} ad_{i,sigma} a_{j,sigma} + 1/2 sum_{i,j,k,l,sigma,tau} v_{i,j,k,l} ad_{i,sigma} ad_{j,tau} a_{l,tau} a_{k,sigma}\nSyntax: construct_spin_molecular_hamiltonian_mpo(tkin, vint, optimize=False)", + }, { .ml_name = "dmrg", .ml_meth = (PyCFunction)Py_dmrg,