adjust order with custom fermionic ops

block-hczhai · Dec 4, 2023 · 88bf4ec · 88bf4ec
1 parent c729757
commit 88bf4ec
Show file tree

Hide file tree

Showing 4 changed files with 174 additions and 108 deletions.
diff --git 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"
- ]
- }
  ]
  }
- ]
-}
+ ],
+ "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
@@ -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
@@ -313,6 +313,59 @@ template <typename FL> struct GeneralFCIDUMP {
  return r;
  }
  };
+ // abelian symmetry case
+ shared_ptr<GeneralFCIDUMP> adjust_order(const string &fermionic_ops,
+ bool merge = true,
+ FP cutoff = (FP)0.0) const {
+ unordered_map<string, uint64_t> r_exprs;
+ unordered_map<char, int8_t> is_op_f;
+ for (size_t it = 0; it < fermionic_ops.length(); it++)
+ is_op_f[fermionic_ops[it]] = 1;
+ vector<vector<uint16_t>> r_indices;
+ vector<vector<FL>> r_data;
+ for (size_t ix = 0; ix < exprs.size(); ix++) {
+ int nn = (int)exprs[ix].length();
+ vector<uint16_t> 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<uint16_t>());
+ r_data.push_back(vector<FL>());
+ }
+ 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<GeneralFCIDUMP> r = make_shared<GeneralFCIDUMP>();
+ 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 <typename T, typename FLX>
  shared_ptr<GeneralFCIDUMP>
  adjust_order_impl(const vector<shared_ptr<T>> &schemes, bool merge = true,

diff --git a/src/pybind/pybind_core.hpp b/src/pybind/pybind_core.hpp
@@ -3665,6 +3665,12 @@ template <typename FL> void bind_general_fcidump(py::module &m) {
  &GeneralFCIDUMP<FL>::initialize_from_qc, py::arg("fcidump"),
  py::arg("elem_type"),
  py::arg("cutoff") = (typename GeneralFCIDUMP<FL>::FP)0.0)
+ .def("adjust_order",
+ (shared_ptr<GeneralFCIDUMP<FL>>(GeneralFCIDUMP<FL>::*)(
+ const string &, bool, typename GeneralFCIDUMP<FL>::FP) const) &
+ GeneralFCIDUMP<FL>::adjust_order,
+ py::arg("fermionic_ops"), py::arg("merge") = true,
+ py::arg("cutoff") = (typename GeneralFCIDUMP<FL>::FP)0.0)
  .def("adjust_order",
  (shared_ptr<GeneralFCIDUMP<FL>>(GeneralFCIDUMP<FL>::*)(
  bool, bool, typename GeneralFCIDUMP<FL>::FP) const) &