diff --git a/docs/source/tutorial/custom-hamiltonian.ipynb b/docs/source/tutorial/custom-hamiltonian.ipynb index 89fa5ed3..9f3c0d50 100644 --- a/docs/source/tutorial/custom-hamiltonian.ipynb +++ b/docs/source/tutorial/custom-hamiltonian.ipynb @@ -1,44 +1,33 @@ { - "nbformat": 4, - "nbformat_minor": 0, - "metadata": { - "colab": { - "provenance": [] - }, - "kernelspec": { - "name": "python3", - "display_name": "Python 3" - }, - "language_info": { - "name": "python" - } - }, "cells": [ { "cell_type": "markdown", + "metadata": { + "id": "gEis3BAAaJLD" + }, "source": [ "# Custom Hamiltonian\n", "\n", "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/block-hczhai/block2-preview/blob/master/docs/source/tutorial/custom-hamiltonian.ipynb)" - ], - "metadata": { - "id": "gEis3BAAaJLD" - } + ] }, { "cell_type": "code", - "source": [ - "!pip install block2==0.5.2rc13 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/\n", - "!pip install pyscf==2.3.0 -qq --progress-bar off" - ], + "execution_count": 1, "metadata": { "id": "2FXmEVkEaJxX" }, - "execution_count": 1, - "outputs": [] + "outputs": [], + "source": [ + "!pip install block2==0.5.3rc4 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/\n", + "!pip install pyscf==2.3.0 -qq --progress-bar off" + ] }, { "cell_type": "markdown", + "metadata": { + "id": "aAO2KvxWaSjo" + }, "source": [ "\n", "In this tutorial, we provide an example python script for performing DMRG using custom Hamiltonians, where the operators and states at local Hilbert space at every site can be redefined. It is also possible to use different local Hilbert space for different sites. New letters can be introduced for representing new operators.\n", @@ -49,13 +38,55 @@ "\n", "In the following example, we implement a custom Hamiltonian for the Hubbard model. In the standard implementation, the on-site term was represented as ``cdCD``. Here we instead introduce a single letter ``N`` for the ``cdCD`` term. For each letter in ``cdCDN`` (representing elementary operators), we define its matrix representation in the local basis in ``site_ops``. The quantum number and number of states in each quantum number at each site (which defines the local Hilbert space) is set in ``site_basis``.\n", "\n" - ], - "metadata": { - "id": "aAO2KvxWaSjo" - } + ] }, { "cell_type": "code", + "execution_count": 2, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "XuXZpeFTb5e2", + "outputId": "878cdd30-3bb5-4b47-d7de-6d1ae8bfdef3" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", + "Time elapsed = 0.342 | E = -6.2256341447 | DW = 2.65e-16\n", + "\n", + "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", + "Time elapsed = 0.456 | E = -6.2256341447 | DE = -1.07e-14 | DW = 4.93e-16\n", + "\n", + "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", + "Time elapsed = 0.585 | E = -6.2256341447 | DE = -1.78e-15 | DW = 9.81e-17\n", + "\n", + "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", + "Time elapsed = 0.694 | E = -6.2256341447 | DE = -3.55e-15 | DW = 1.20e-16\n", + "\n", + "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", + "Time elapsed = 0.873 | E = -6.2256341447 | DE = -8.88e-16 | DW = 4.50e-20\n", + "\n", + "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", + "Time elapsed = 1.028 | E = -6.2256341447 | DE = -2.66e-15 | DW = 4.87e-20\n", + "\n", + "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", + "Time elapsed = 1.178 | E = -6.2256341447 | DE = 5.33e-15 | DW = 3.86e-20\n", + "\n", + "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", + "Time elapsed = 1.320 | E = -6.2256341447 | DE = 1.78e-15 | DW = 4.94e-20\n", + "\n", + "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", + "Time elapsed = 1.472 | E = -6.2256341447 | DE = 0.00e+00 | DW = 2.86e-20\n", + "\n", + "DMRG energy = -6.225634144657919\n" + ] + } + ], "source": [ "from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes\n", "import numpy as np\n", @@ -93,71 +124,71 @@ "b.add_term(\"N\", np.array([i for i in range(L)]), U)\n", "\n", "# [Part C] Perform DMRG\n", - "mpo = driver.get_mpo(b.finalize(adjust_order=True), algo_type=MPOAlgorithmTypes.FastBipartite)\n", + "mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops=\"cdCD\"), algo_type=MPOAlgorithmTypes.FastBipartite)\n", "mps = driver.get_random_mps(tag=\"KET\", bond_dim=250, nroots=1)\n", "energy = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[250] * 4 + [500] * 4,\n", " noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=30, iprint=1)\n", "print(\"DMRG energy = %20.15f\" % energy)" - ], + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "AoKOfFPEbrHw" + }, + "source": [ + "## The Hubbard-Holstein Model\n", + "\n", + "The above script can be easily extended to treat phonons." + ] + }, + { + "cell_type": "code", + "execution_count": 3, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, - "id": "XuXZpeFTb5e2", - "outputId": "878cdd30-3bb5-4b47-d7de-6d1ae8bfdef3" + "id": "JqphdawZboo9", + "outputId": "919844ac-eed3-44f2-af2b-dd6ac5d5ee55" }, - "execution_count": 2, "outputs": [ { - "output_type": "stream", "name": "stdout", + "output_type": "stream", "text": [ "\n", "Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", - "Time elapsed = 0.342 | E = -6.2256341447 | DW = 2.65e-16\n", + "Time elapsed = 108.945 | E = -6.9568929229 | DW = 3.62e-09\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", - "Time elapsed = 0.456 | E = -6.2256341447 | DE = -1.07e-14 | DW = 4.93e-16\n", + "Time elapsed = 142.030 | E = -6.9568932112 | DE = -2.88e-07 | DW = 3.07e-19\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", - "Time elapsed = 0.585 | E = -6.2256341447 | DE = -1.78e-15 | DW = 9.81e-17\n", + "Time elapsed = 154.346 | E = -6.9568932112 | DE = -1.78e-15 | DW = 1.38e-19\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", - "Time elapsed = 0.694 | E = -6.2256341447 | DE = -3.55e-15 | DW = 1.20e-16\n", + "Time elapsed = 165.954 | E = -6.9568932112 | DE = -1.78e-15 | DW = 6.71e-20\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", - "Time elapsed = 0.873 | E = -6.2256341447 | DE = -8.88e-16 | DW = 4.50e-20\n", + "Time elapsed = 176.280 | E = -6.9568932112 | DE = -8.88e-16 | DW = 7.26e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", - "Time elapsed = 1.028 | E = -6.2256341447 | DE = -2.66e-15 | DW = 4.87e-20\n", + "Time elapsed = 186.605 | E = -6.9568932112 | DE = 2.66e-15 | DW = 6.17e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", - "Time elapsed = 1.178 | E = -6.2256341447 | DE = 5.33e-15 | DW = 3.86e-20\n", + "Time elapsed = 196.209 | E = -6.9568932112 | DE = -1.78e-15 | DW = 9.34e-20\n", "\n", "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", - "Time elapsed = 1.320 | E = -6.2256341447 | DE = 1.78e-15 | DW = 4.94e-20\n", + "Time elapsed = 205.734 | E = -6.9568932112 | DE = 2.66e-15 | DW = 6.36e-20\n", "\n", "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", - "Time elapsed = 1.472 | E = -6.2256341447 | DE = 0.00e+00 | DW = 2.86e-20\n", + "Time elapsed = 215.378 | E = -6.9568932112 | DE = -9.77e-15 | DW = 7.87e-20\n", "\n", - "DMRG energy = -6.225634144657919\n" + "DMRG energy = -6.956893211180049\n" ] } - ] - }, - { - "cell_type": "markdown", - "source": [ - "## The Hubbard-Holstein Model\n", - "\n", - "The above script can be easily extended to treat phonons." ], - "metadata": { - "id": "AoKOfFPEbrHw" - } - }, - { - "cell_type": "code", "source": [ "from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes\n", "import numpy as np\n", @@ -214,57 +245,26 @@ "b.add_term(\"CDF\", np.array([[i, i, i + N_SITES_ELEC] for i in range(N_SITES_ELEC)]).ravel(), G)\n", "\n", "# [Part C] Perform DMRG\n", - "mpo = driver.get_mpo(b.finalize(adjust_order=True), algo_type=MPOAlgorithmTypes.FastBipartite)\n", + "mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops=\"cdCD\"), algo_type=MPOAlgorithmTypes.FastBipartite)\n", "mps = driver.get_random_mps(tag=\"KET\", bond_dim=250, nroots=1)\n", "energy = driver.dmrg(mpo, mps, n_sweeps=10, bond_dims=[250] * 4 + [500] * 4,\n", " noises=[1e-4] * 4 + [1e-5] * 4 + [0], thrds=[1e-10] * 8, dav_max_iter=30, iprint=1)\n", "print(\"DMRG energy = %20.15f\" % energy)" - ], - "metadata": { - "colab": { - "base_uri": "https://localhost:8080/" - }, - "id": "JqphdawZboo9", - "outputId": "919844ac-eed3-44f2-af2b-dd6ac5d5ee55" - }, - "execution_count": 3, - "outputs": [ - { - "output_type": "stream", - "name": "stdout", - "text": [ - "\n", - "Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", - "Time elapsed = 108.945 | E = -6.9568929229 | DW = 3.62e-09\n", - "\n", - "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", - "Time elapsed = 142.030 | E = -6.9568932112 | DE = -2.88e-07 | DW = 3.07e-19\n", - "\n", - "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", - "Time elapsed = 154.346 | E = -6.9568932112 | DE = -1.78e-15 | DW = 1.38e-19\n", - "\n", - "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", - "Time elapsed = 165.954 | E = -6.9568932112 | DE = -1.78e-15 | DW = 6.71e-20\n", - "\n", - "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", - "Time elapsed = 176.280 | E = -6.9568932112 | DE = -8.88e-16 | DW = 7.26e-20\n", - "\n", - "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", - "Time elapsed = 186.605 | E = -6.9568932112 | DE = 2.66e-15 | DW = 6.17e-20\n", - "\n", - "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", - "Time elapsed = 196.209 | E = -6.9568932112 | DE = -1.78e-15 | DW = 9.34e-20\n", - "\n", - "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", - "Time elapsed = 205.734 | E = -6.9568932112 | DE = 2.66e-15 | DW = 6.36e-20\n", - "\n", - "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", - "Time elapsed = 215.378 | E = -6.9568932112 | DE = -9.77e-15 | DW = 7.87e-20\n", - "\n", - "DMRG energy = -6.956893211180049\n" - ] - } ] } - ] -} \ No newline at end of file + ], + "metadata": { + "colab": { + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/pyblock2/driver/core.py b/pyblock2/driver/core.py index ca26a73b..24e60523 100644 --- a/pyblock2/driver/core.py +++ b/pyblock2/driver/core.py @@ -5192,9 +5192,16 @@ def iscale(self, d): self.data.data[i] = self.bw.VectorFL(d * np.array(ix)) return self - def finalize(self, adjust_order=True, merge=True, is_drt=False): + def finalize(self, adjust_order=True, merge=True, is_drt=False, fermionic_ops=None): if adjust_order: - self.data = self.data.adjust_order(merge=merge, is_drt=is_drt) + if fermionic_ops is not None: + assert (SymmetryTypes.SU2 not in self.bw.symm_type + and SymmetryTypes.PHSU2 not in self.bw.symm_type + and SymmetryTypes.SO4 not in self.bw.symm_type + and SymmetryTypes.SO3 not in self.bw.symm_type) + self.data = self.data.adjust_order(fermionic_ops, merge=merge) + else: + self.data = self.data.adjust_order(merge=merge, is_drt=is_drt) elif merge: self.data.merge_terms() return self.data diff --git a/src/core/integral_general.hpp b/src/core/integral_general.hpp index a6138159..62322cb6 100644 --- a/src/core/integral_general.hpp +++ b/src/core/integral_general.hpp @@ -313,6 +313,59 @@ template struct GeneralFCIDUMP { return r; } }; + // abelian symmetry case + shared_ptr adjust_order(const string &fermionic_ops, + bool merge = true, + FP cutoff = (FP)0.0) const { + unordered_map r_exprs; + unordered_map is_op_f; + for (size_t it = 0; it < fermionic_ops.length(); it++) + is_op_f[fermionic_ops[it]] = 1; + vector> r_indices; + vector> r_data; + for (size_t ix = 0; ix < exprs.size(); ix++) { + int nn = (int)exprs[ix].length(); + vector idx_idx(nn); + for (size_t i = 0; i < (nn == 0 ? 1 : indices[ix].size()); + i += (nn == 0 ? 1 : nn)) { + for (int j = 0; j < nn; j++) + idx_idx[j] = j; + string xex = exprs[ix]; + int8_t n = 0; + for (int xi = 0; xi < (int)nn - 1; xi++) + for (int xj = xi; xj >= 0; xj--) + if (indices[ix][i + idx_idx[xj]] > + indices[ix][i + idx_idx[xj + 1]]) { + swap(idx_idx[xj], idx_idx[xj + 1]); + swap(xex[xj], xex[xj + 1]); + n ^= (is_op_f.count(xex[xj]) && + is_op_f.count(xex[xj + 1])); + } + if (!r_exprs.count(xex)) { + r_exprs[xex] = r_data.size(); + r_indices.push_back(vector()); + r_data.push_back(vector()); + } + uint64_t ir = r_exprs.at(xex); + for (int j = 0; j < nn; j++) + r_indices[ir].push_back(indices[ix][i + idx_idx[j]]); + r_data[ir].push_back((FL)(FP)(1 - (n << 1)) * + data[ix][nn == 0 ? i : i / nn]); + } + } + shared_ptr r = make_shared(); + r->params = params; + r->const_e = const_e; + r->elem_type = elem_type; + r->exprs.resize(r_exprs.size()); + for (auto &x : r_exprs) + r->exprs[x.second] = x.first; + r->indices = r_indices; + r->data = r_data; + if (merge) + r->merge_terms(cutoff); + return r; + } template shared_ptr adjust_order_impl(const vector> &schemes, bool merge = true, diff --git a/src/pybind/pybind_core.hpp b/src/pybind/pybind_core.hpp index d93edc29..347127af 100644 --- a/src/pybind/pybind_core.hpp +++ b/src/pybind/pybind_core.hpp @@ -3665,6 +3665,12 @@ template void bind_general_fcidump(py::module &m) { &GeneralFCIDUMP::initialize_from_qc, py::arg("fcidump"), py::arg("elem_type"), py::arg("cutoff") = (typename GeneralFCIDUMP::FP)0.0) + .def("adjust_order", + (shared_ptr>(GeneralFCIDUMP::*)( + const string &, bool, typename GeneralFCIDUMP::FP) const) & + GeneralFCIDUMP::adjust_order, + py::arg("fermionic_ops"), py::arg("merge") = true, + py::arg("cutoff") = (typename GeneralFCIDUMP::FP)0.0) .def("adjust_order", (shared_ptr>(GeneralFCIDUMP::*)( bool, bool, typename GeneralFCIDUMP::FP) const) &