{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "gEis3BAAaJLD" }, "source": [ "# Custom Hamiltonians\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-hamiltonians.ipynb)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "2FXmEVkEaJxX" }, "outputs": [], "source": [ "!pip install block2==0.5.3rc17 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/" ] }, { "cell_type": "markdown", "metadata": { "id": "aAO2KvxWaSjo" }, "source": [ "\n", "In this tutorial, we provide an example python scripts 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 (the operator name can only be a single lower or upper case character).\n", "\n", "Note the following examples are only supposed to work in the Abelian symmetry modes (``SZ``, ``SZ|CPX``, ``SGF``, ``SGFCPX``, ``SAny``, or ``SAny|CPX``) for Abelian symmetries, and non-Abelian symmetry modes (``SAnySU2`` or ``SAnySU2|CPX``) for non-Abelian symmetries, respectively.\n", "\n", "## The Hubbard Model\n", "\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) are set in ``site_basis``.\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "XuXZpeFTb5e2", "outputId": "d92ee9bc-e0e8-40e5-a4ff-a271948703cd" }, "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.167 | E = -6.2256341447 | DW = 2.65116e-16\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.240 | E = -6.2256341447 | DE = -8.88e-15 | DW = 4.93007e-16\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.305 | E = -6.2256341447 | DE = 1.78e-15 | DW = 9.66767e-17\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.367 | E = -6.2256341447 | DE = 1.78e-15 | DW = 1.20251e-16\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.431 | E = -6.2256341447 | DE = -8.88e-16 | DW = 3.71514e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.493 | E = -6.2256341447 | DE = -3.55e-15 | DW = 6.52472e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.560 | E = -6.2256341447 | DE = -1.78e-15 | DW = 4.11603e-20\n", "\n", "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.626 | E = -6.2256341447 | DE = 0.00e+00 | DW = 4.46198e-20\n", "\n", "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", "Time elapsed = 0.674 | E = -6.2256341447 | DE = 7.11e-15 | DW = 4.96657e-20\n", "\n", "DMRG energy = -6.225634144657917\n" ] } ], "source": [ "from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes\n", "import numpy as np\n", "\n", "L = 8\n", "U = 2\n", "N_ELEC = 8\n", "\n", "driver = DMRGDriver(scratch=\"./tmp\", symm_type=SymmetryTypes.SZ, n_threads=4)\n", "driver.initialize_system(n_sites=L, n_elec=N_ELEC, spin=0)\n", "\n", "# [Part A] Set states and matrix representation of operators in local Hilbert space\n", "site_basis, site_ops = [], []\n", "Q = driver.bw.SX # quantum number wrapper (n_elec, 2 * spin, point group irrep)\n", "\n", "for k in range(L):\n", " basis = [(Q(0, 0, 0), 1), (Q(1, 1, 0), 1), (Q(1, -1, 0), 1), (Q(2, 0, 0), 1)] # [0ab2]\n", " ops = {\n", " \"\": np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]), # identity\n", " \"c\": np.array([[0, 0, 0, 0], [1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0]]), # alpha+\n", " \"d\": np.array([[0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 0, 0]]), # alpha\n", " \"C\": np.array([[0, 0, 0, 0], [0, 0, 0, 0], [1, 0, 0, 0], [0, -1, 0, 0]]), # beta+\n", " \"D\": np.array([[0, 0, 1, 0], [0, 0, 0, -1], [0, 0, 0, 0], [0, 0, 0, 0]]), # beta\n", " \"N\": np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1]]), # cdCD\n", " }\n", " site_basis.append(basis)\n", " site_ops.append(ops)\n", "\n", "# [Part B] Set Hamiltonian terms\n", "driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)\n", "b = driver.expr_builder()\n", "\n", "b.add_term(\"cd\", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)\n", "b.add_term(\"CD\", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)\n", "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, 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": "JqphdawZboo9", "outputId": "bf7815ae-6db6-4831-c1ed-fb0cf2ee63d3" }, "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 = 94.146 | E = -6.9568929542 | DW = 3.61640e-09\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 126.358 | E = -6.9568932112 | DE = -2.57e-07 | DW = 3.01342e-19\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 139.180 | E = -6.9568932112 | DE = 4.44e-15 | DW = 1.23451e-19\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 149.110 | E = -6.9568932112 | DE = 1.78e-15 | DW = 7.41629e-20\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 158.521 | E = -6.9568932112 | DE = 3.55e-15 | DW = 7.71387e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 166.802 | E = -6.9568932112 | DE = -5.33e-15 | DW = 6.44981e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 175.980 | E = -6.9568932112 | DE = 8.88e-16 | DW = 8.11423e-20\n", "\n", "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 184.501 | E = -6.9568932112 | DE = -2.66e-15 | DW = 6.41555e-20\n", "\n", "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", "Time elapsed = 192.172 | E = -6.9568932112 | DE = 0.00e+00 | DW = 7.78301e-20\n", "\n", "DMRG energy = -6.956893211180857\n" ] } ], "source": [ "from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes\n", "import numpy as np\n", "\n", "N_SITES_ELEC, N_SITES_PH, N_ELEC = 4, 4, 4\n", "N_PH, U, OMEGA, G = 11, 2, 0.25, 0.5\n", "L = N_SITES_ELEC + N_SITES_PH\n", "\n", "driver = DMRGDriver(scratch=\"./tmp\", symm_type=SymmetryTypes.SZ, n_threads=4)\n", "driver.initialize_system(n_sites=L, n_elec=N_ELEC, spin=0)\n", "\n", "# [Part A] Set states and matrix representation of operators in local Hilbert space\n", "site_basis, site_ops = [], []\n", "Q = driver.bw.SX # quantum number wrapper (n_elec, 2 * spin, point group irrep)\n", "\n", "for k in range(L):\n", " if k < N_SITES_ELEC:\n", " # electron part\n", " basis = [(Q(0, 0, 0), 1), (Q(1, 1, 0), 1), (Q(1, -1, 0), 1), (Q(2, 0, 0), 1)] # [0ab2]\n", " ops = {\n", " \"\": np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]), # identity\n", " \"c\": np.array([[0, 0, 0, 0], [1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0]]), # alpha+\n", " \"d\": np.array([[0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 0, 0]]), # alpha\n", " \"C\": np.array([[0, 0, 0, 0], [0, 0, 0, 0], [1, 0, 0, 0], [0, -1, 0, 0]]), # beta+\n", " \"D\": np.array([[0, 0, 1, 0], [0, 0, 0, -1], [0, 0, 0, 0], [0, 0, 0, 0]]), # beta\n", " }\n", " else:\n", " # phonon part\n", " basis = [(Q(0, 0, 0), N_PH)]\n", " ops = {\n", " \"\": np.identity(N_PH), # identity\n", " \"E\": np.diag(np.sqrt(np.arange(1, N_PH)), k=-1), # ph+\n", " \"F\": np.diag(np.sqrt(np.arange(1, N_PH)), k=1), # ph\n", " }\n", " site_basis.append(basis)\n", " site_ops.append(ops)\n", "\n", "# [Part B] Set Hamiltonian terms in Hubbard-Holstein model\n", "driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)\n", "b = driver.expr_builder()\n", "\n", "# electron part\n", "b.add_term(\"cd\", np.array([[i, i + 1, i + 1, i] for i in range(N_SITES_ELEC - 1)]).ravel(), -1)\n", "b.add_term(\"CD\", np.array([[i, i + 1, i + 1, i] for i in range(N_SITES_ELEC - 1)]).ravel(), -1)\n", "b.add_term(\"cdCD\", np.array([[i, i, i, i] for i in range(N_SITES_ELEC)]).ravel(), U)\n", "\n", "# phonon part\n", "b.add_term(\"EF\", np.array([[i + N_SITES_ELEC, ] * 2 for i in range(N_SITES_PH)]).ravel(), OMEGA)\n", "\n", "# interaction part\n", "b.add_term(\"cdE\", np.array([[i, i, i + N_SITES_ELEC] for i in range(N_SITES_ELEC)]).ravel(), G)\n", "b.add_term(\"cdF\", np.array([[i, i, i + N_SITES_ELEC] for i in range(N_SITES_ELEC)]).ravel(), G)\n", "b.add_term(\"CDE\", np.array([[i, i, i + N_SITES_ELEC] for i in range(N_SITES_ELEC)]).ravel(), G)\n", "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, 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": "irQr9N_HOQNw" }, "source": [ "## Custom Symmetry Groups\n", "\n", "In the following we show how to set the custom symmetry groups. The symmetry mode ``SymmetryTypes.SAny`` (or ``SymmetryTypes.SAny | SymmetryTypes.CPX`` if complex number is required) should be used for this purpose. Currently we support the definition of symmetry group as an arbitrary direct product of up to six Abelian symmetry sub-groups. Possible sub-group names are \"U1\", \"Z1\", \"Z2\", \"Z3\", ..., \"Z2055\", \"U1Fermi\", \"Z1Fermi\", \"Z2Fermi\", \"Z3Fermi\", ..., \"Z2055Fermi\", \"LZ\", and \"AbelianPG\". The names with the suffix \"Fermi\" should be used for Fermion symmetries. The names without the suffix \"Fermi\" should be used for spin or Boson symmetries. The non-Abelian sub-group names \"SU2\" and \"SU2Fermi\" can be used with some restrictions, see the example in later sections. The ``DMRGDriver.set_symmetry_groups(sub_group_name_1: str, sub_group_name_2: str, ...)`` method can be used to set the symmetry sub-groups. The number of arguments in the quantum number wrapper ``Q`` should then match the number of sub-group names given in ``DMRGDriver.set_symmetry_groups``.\n", "\n", "As a first example, we use the custom symmetry group syntax to recompute the Hubbard model. We first use $U(1) \\times U(1)$ symmetry, which should be equivalent to the previous ``SZ`` mode example." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "J2Ed7zljSOLS", "outputId": "f6f530e3-5270-4c1f-fd38-1dfd599ddb38" }, "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.247 | E = -6.2256341447 | DW = 5.28457e-16\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.344 | E = -6.2256341447 | DE = -8.88e-15 | DW = 6.68285e-16\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.446 | E = -6.2256341447 | DE = -4.44e-15 | DW = 1.63432e-16\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.543 | E = -6.2256341447 | DE = 1.78e-15 | DW = 2.16706e-16\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.623 | E = -6.2256341447 | DE = 3.55e-15 | DW = 2.88426e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.692 | E = -6.2256341447 | DE = -5.33e-15 | DW = 2.28873e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.764 | E = -6.2256341447 | DE = 8.88e-16 | DW = 3.56876e-20\n", "\n", "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.831 | E = -6.2256341447 | DE = 1.78e-15 | DW = 3.86969e-20\n", "\n", "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", "Time elapsed = 0.880 | E = -6.2256341447 | DE = -2.66e-15 | DW = 5.43709e-20\n", "\n", "DMRG energy = -6.225634144658391\n" ] } ], "source": [ "from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes\n", "import numpy as np\n", "\n", "L = 8\n", "U = 2\n", "N_ELEC = 8\n", "TWO_SZ = 0\n", "\n", "driver = DMRGDriver(scratch=\"./tmp\", symm_type=SymmetryTypes.SAny, n_threads=4)\n", "\n", "# quantum number wrapper (U1 / n_elec, U1 / 2*Sz)\n", "driver.set_symmetry_groups(\"U1Fermi\", \"U1\")\n", "Q = driver.bw.SX\n", "\n", "# [Part A] Set states and matrix representation of operators in local Hilbert space\n", "site_basis, site_ops = [], []\n", "\n", "for k in range(L):\n", " basis = [(Q(0, 0), 1), (Q(1, 1), 1), (Q(1, -1), 1), (Q(2, 0), 1)] # [0ab2]\n", " ops = {\n", " \"\": np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]), # identity\n", " \"c\": np.array([[0, 0, 0, 0], [1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0]]), # alpha+\n", " \"d\": np.array([[0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1], [0, 0, 0, 0]]), # alpha\n", " \"C\": np.array([[0, 0, 0, 0], [0, 0, 0, 0], [1, 0, 0, 0], [0, -1, 0, 0]]), # beta+\n", " \"D\": np.array([[0, 0, 1, 0], [0, 0, 0, -1], [0, 0, 0, 0], [0, 0, 0, 0]]), # beta\n", " }\n", " site_basis.append(basis)\n", " site_ops.append(ops)\n", "\n", "# [Part B] Set Hamiltonian terms\n", "driver.initialize_system(n_sites=L, vacuum=Q(0, 0), target=Q(N_ELEC, TWO_SZ), hamil_init=False)\n", "driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)\n", "b = driver.expr_builder()\n", "\n", "b.add_term(\"cd\", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)\n", "b.add_term(\"CD\", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)\n", "b.add_term(\"cdCD\", np.array([i for i in range(L) for _ in range(4)]), U)\n", "\n", "# [Part C] Perform DMRG\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": "GpXYx6P-SUu9" }, "source": [ "As a second example, we recompute the Hubbard model using $Z_2 \\times Z_2$ symmetry. This time we cannot easily target the $N_{\\mathrm{elec}} = 8$ symmetry sector. Instead, we compute a few excited states and compute the $\\langle N\\rangle$ to identify the correct state." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "Sy3eG4KuS9qP", "outputId": "fbaccccb-0674-412b-c2ea-fdb27476f894" }, "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 = 36.478 | E[ 10] = -6.9958183038 -6.5705972881 -6.5705972845 -6.5705972809 -6.2773135425 -6.2256341359 -6.1122535707 -6.1122534910 -6.1122534358 -6.0320900884 | DW = 4.44265e-09\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 40.236 | E[ 10] = -6.9958183059 -6.5705972895 -6.5705972895 -6.5705972895 -6.2773135514 -6.2256341447 -6.1122535739 -6.1122535739 -6.1122535739 -6.0320902499 | DE = -1.62e-07 | DW = 1.12165e-09\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 43.608 | E[ 10] = -6.9958183059 -6.5705972895 -6.5705972895 -6.5705972895 -6.2773135514 -6.2256341447 -6.1122535739 -6.1122535739 -6.1122535739 -6.0320902499 | DE = 3.43e-12 | DW = 1.18398e-09\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 46.245 | E[ 10] = -6.9958183059 -6.5705972895 -6.5705972895 -6.5705972895 -6.2773135514 -6.2256341447 -6.1122535739 -6.1122535739 -6.1122535739 -6.0320902499 | DE = -1.75e-12 | DW = 1.12063e-09\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 49.844 | E[ 10] = -6.9958183059 -6.5705972895 -6.5705972895 -6.5705972895 -6.2773135514 -6.2256341447 -6.1122535739 -6.1122535739 -6.1122535739 -6.0320902499 | DE = -1.88e-12 | DW = 7.51709e-17\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 52.386 | E[ 10] = -6.9958183059 -6.5705972895 -6.5705972895 -6.5705972895 -6.2773135514 -6.2256341447 -6.1122535739 -6.1122535739 -6.1122535739 -6.0320902499 | DE = -8.88e-15 | DW = 1.55303e-17\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 55.928 | E[ 10] = -6.9958183059 -6.5705972895 -6.5705972895 -6.5705972895 -6.2773135514 -6.2256341447 -6.1122535739 -6.1122535739 -6.1122535739 -6.0320902499 | DE = 0.00e+00 | DW = 1.59309e-17\n", "\n", "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 58.510 | E[ 10] = -6.9958183059 -6.5705972895 -6.5705972895 -6.5705972895 -6.2773135514 -6.2256341447 -6.1122535739 -6.1122535739 -6.1122535739 -6.0320902499 | DE = 1.24e-14 | DW = 6.42978e-18\n", "\n", "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", "Time elapsed = 60.677 | E[ 10] = -6.9958183059 -6.5705972895 -6.5705972895 -6.5705972895 -6.2773135514 -6.2256341447 -6.1122535739 -6.1122535739 -6.1122535739 -6.0320902499 | DE = 0.00e+00 | DW = 9.28475e-19\n", "\n", "Root = 0 = -6.995818305899173 = 6.000\n", "Root = 1 = -6.570597289543093 = 6.000\n", "Root = 2 = -6.570597289542542 = 6.000\n", "Root = 3 = -6.570597289541823 = 6.000\n", "Root = 4 = -6.277313551397241 = 6.000\n", "Root = 5 = -6.225634144677118 = 8.000\n", "Root = 6 = -6.112253573866418 = 6.000\n", "Root = 7 = -6.112253573864956 = 6.000\n", "Root = 8 = -6.112253573862747 = 6.000\n", "Root = 9 = -6.032090249939520 = 4.000\n" ] } ], "source": [ "from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes\n", "import numpy as np\n", "\n", "L = 8\n", "U = 2\n", "N_ELEC = 8\n", "TWO_SZ = 0\n", "\n", "driver = DMRGDriver(scratch=\"./tmp\", symm_type=SymmetryTypes.SAny, n_threads=4)\n", "\n", "# quantum number wrapper (Z2 / n_elec, Z2 / 2*Sz)\n", "driver.set_symmetry_groups(\"Z2Fermi\", \"Z2\")\n", "Q = driver.bw.SX\n", "\n", "# [Part A] Set states and matrix representation of operators in local Hilbert space\n", "site_basis, site_ops = [], []\n", "\n", "for k in range(L):\n", " basis = [(Q(0, 0), 2), (Q(1, 1), 2)] # [02ab]\n", " ops = {\n", " # note the order of row and column is different from the U1xU1 case\n", " \"\": np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]), # identity\n", " \"c\": np.array([[0, 0, 0, 0], [0, 0, 0, 1], [1, 0, 0, 0], [0, 0, 0, 0]]), # alpha+\n", " \"d\": np.array([[0, 0, 1, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 1, 0, 0]]), # alpha\n", " \"C\": np.array([[0, 0, 0, 0], [0, 0, -1, 0], [0, 0, 0, 0], [1, 0, 0, 0]]), # beta+\n", " \"D\": np.array([[0, 0, 0, 1], [0, 0, 0, 0], [0, -1, 0, 0], [0, 0, 0, 0]]), # beta\n", " }\n", " site_basis.append(basis)\n", " site_ops.append(ops)\n", "\n", "# [Part B] Set Hamiltonian terms\n", "driver.initialize_system(n_sites=L, vacuum=Q(0, 0), target=Q(N_ELEC % 2, TWO_SZ % 2), hamil_init=False)\n", "driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)\n", "b = driver.expr_builder()\n", "\n", "b.add_term(\"cd\", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)\n", "b.add_term(\"CD\", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), -1)\n", "b.add_term(\"cdCD\", np.array([i for i in range(L) for _ in range(4)]), U)\n", "\n", "# [Part C] Perform state-averaged DMRG\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=10)\n", "energies = 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=200, iprint=1)\n", "\n", "# [Part D] Check particle number expectations\n", "b = driver.expr_builder()\n", "b.add_term(\"cd\", np.array([[i, i] for i in range(L)]).ravel(), 1)\n", "b.add_term(\"CD\", np.array([[i, i] for i in range(L)]).ravel(), 1)\n", "partile_n_mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops=\"cdCD\"), algo_type=MPOAlgorithmTypes.FastBipartite)\n", "\n", "kets = [driver.split_mps(mps, ir, tag=\"KET-%d\" % ir) for ir in range(mps.nroots)]\n", "for ir in range(mps.nroots):\n", " n_expt = driver.expectation(kets[ir], partile_n_mpo, kets[ir])\n", " print(\"Root = %d = %20.15f = %10.3f\" % (ir, energies[ir], n_expt))" ] }, { "cell_type": "markdown", "metadata": { "id": "lMtwBV41VQHF" }, "source": [ "## Bose-Hubbard Model\n", "\n", "In this example we will find the ground state of 1D Bose-Hubbard model. Unlike the Holstein model, in this model the boson number is conserved as an Abelian U(1) symmetry. The definition of the Hamiltonian can be found in https://en.wikipedia.org/wiki/Bose%E2%80%93Hubbard_model, which is given by\n", "\n", "$$\n", "\\hat{H} = -t \\sum_{\\langle ij\\rangle} \\Big( b^\\dagger_i\n", " b_j + b^\\dagger_j b_i \\Big)\n", " + \\frac{U}{2} \\sum_i n_i (n_i - 1) - \\mu \\sum_i n_i\n", "$$\n", "\n", "In the following script, we consider a 1D chain including $L= 10$ sites with open boundary condition, $t = 1, U = 0.1$ and $\\mu = 0$. The total number of boson in the many-body state is set to 10, and the maximal number of boson per site is 5." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "XdwaoNOjXINz", "outputId": "c6a8f9af-10a5-48bc-8fa3-862fcefc2aba" }, "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 = 1.026 | E = -18.5867775673 | DW = 4.33791e-12\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 3.090 | E = -18.5873694219 | DE = -5.92e-04 | DW = 3.77544e-16\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 3.436 | E = -18.5873694219 | DE = -7.11e-15 | DW = 3.97569e-17\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 3.813 | E = -18.5873694219 | DE = -8.30e-12 | DW = 1.83165e-17\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 4.397 | E = -18.5873694219 | DE = 3.55e-15 | DW = 1.35857e-19\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 4.703 | E = -18.5873694219 | DE = -3.55e-15 | DW = 1.43599e-19\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 4.971 | E = -18.5873694219 | DE = -1.07e-14 | DW = 8.50577e-20\n", "\n", "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 5.226 | E = -18.5873694219 | DE = 3.55e-15 | DW = 1.44991e-19\n", "\n", "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", "Time elapsed = 5.415 | E = -18.5873694219 | DE = 0.00e+00 | DW = 8.57382e-20\n", "\n", "DMRG energy = -18.587369421896540 (per site = -1.858737)\n" ] } ], "source": [ "from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes\n", "import numpy as np\n", "\n", "L = 10\n", "T = 1\n", "U = 0.1\n", "MU = 0\n", "NB_MAX = 5 # max n_boson per site\n", "N_BOSON = 10\n", "\n", "driver = DMRGDriver(scratch=\"./tmp\", symm_type=SymmetryTypes.SAny, n_threads=4)\n", "\n", "driver.set_symmetry_groups(\"U1\")\n", "Q = driver.bw.SX\n", "\n", "# [Part A] Set states and matrix representation of operators in local Hilbert space\n", "site_basis, site_ops = [], []\n", "\n", "for k in range(L):\n", " basis = [(Q(i), 1) for i in range(NB_MAX + 1)] # [012..NB_MAX]\n", " ops = {\n", " \"\": np.identity(NB_MAX + 1), # identity\n", " \"C\": np.diag(np.sqrt(np.arange(1, NB_MAX + 1)), k=-1), # b+\n", " \"D\": np.diag(np.sqrt(np.arange(1, NB_MAX + 1)), k=1), # b\n", " \"N\": np.diag(np.arange(0, NB_MAX + 1), k=0), # particle number\n", " }\n", " site_basis.append(basis)\n", " site_ops.append(ops)\n", "\n", "# [Part B] Set Hamiltonian terms\n", "driver.initialize_system(n_sites=L, vacuum=Q(0), target=Q(N_BOSON), hamil_init=False)\n", "driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)\n", "b = driver.expr_builder()\n", "\n", "b.add_term(\"CD\", np.array([j for i in range(L - 1) for j in [i, i + 1, i + 1, i]]), -T)\n", "b.add_term(\"N\", np.array([i for i in range(L)]), -(MU + U / 2))\n", "b.add_term(\"NN\", np.array([j for i in range(L) for j in [i, i]]), U / 2)\n", "\n", "# [Part C] Perform DMRG\n", "mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops=\"\"), 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 (per site = %10.6f)\" % (energy, energy / L))" ] }, { "cell_type": "markdown", "metadata": { "id": "XhU8f8s3TFYK" }, "source": [ "## SU(3) Heisenberg Model\n", "\n", "In this example we will find the ground state of the 1D SU(3) Heisenberg model using the custom symemtry group syntax. We will only use Abelian symmetry groups (for quantum numbers $S_z$ and $Q_z$) for this problem. The model used here can be found in Eq. (2) in *Phys. Rev. B* **79**, 012408 (2009)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "QlYXVVljUKVU", "outputId": "7eff7afe-02e6-4e47-e11c-fb6c2f376ef2" }, "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 = 11.778 | E = -0.5120663491 | DW = 1.30017e-05\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 18.660 | E = -0.5145229771 | DE = -2.46e-03 | DW = 4.16418e-08\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 23.883 | E = -0.5145506979 | DE = -2.77e-05 | DW = 3.51827e-07\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 26.629 | E = -0.5145508464 | DE = -1.48e-07 | DW = 4.79795e-07\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 34.588 | E = -0.5145512453 | DE = -3.99e-07 | DW = 4.81635e-09\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 47.894 | E = -0.5145512513 | DE = -6.07e-09 | DW = 4.36734e-09\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 61.178 | E = -0.5145512510 | DE = 2.96e-10 | DW = 1.36107e-09\n", "\n", "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 74.667 | E = -0.5145512506 | DE = 4.30e-10 | DW = 8.13159e-10\n", "\n", "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", "Time elapsed = 84.889 | E = -0.5145512503 | DE = 2.84e-10 | DW = 9.28021e-19\n", "\n", "DMRG energy (per site) = -0.514551250318804\n" ] } ], "source": [ "from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes\n", "import numpy as np\n", "\n", "L = 72\n", "\n", "driver = DMRGDriver(scratch=\"./tmp\", symm_type=SymmetryTypes.SAny, n_threads=4)\n", "\n", "# quantum number wrapper (2Sz / X, 2Qz / Y)\n", "driver.set_symmetry_groups(\"U1\", \"U1\")\n", "Q = driver.bw.SX\n", "\n", "# [Part A] Set states and matrix representation of operators in local Hilbert space\n", "site_basis, site_ops = [], []\n", "\n", "# Gell Mann operators\n", "lambda_ops = {\n", " \"L1\": np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]]), # lambda_1\n", " \"L2\": np.array([[0, -1j, 0], [1j, 0, 0], [0, 0, 0]]), # lambda_2\n", " \"L3\": np.array([[1, 0, 0], [0, -1, 0], [0, 0, 0]]), # lambda_3\n", " \"L4\": np.array([[0, 0, 1], [0, 0, 0], [1, 0, 0]]), # lambda_4\n", " \"L5\": np.array([[0, 0, -1j], [0, 0, 0], [1j, 0, 0]]), # lambda_5\n", " \"L6\": np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0]]), # lambda_6\n", " \"L7\": np.array([[0, 0, 0], [0, 0, -1j], [0, 1j, 0]]), # lambda_7\n", " \"L8\": np.array([[1, 0, 0], [0, 1, 0], [0, 0, -2]]) / 3 ** 0.5, # lambda_8\n", "}\n", "\n", "for k in range(L):\n", " basis = [(Q(1, 1), 1), (Q(-1, 1), 1), (Q(0, -2), 1)]\n", " ops = {\n", " \"\": np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), # identity\n", " \"T\": (lambda_ops[\"L1\"] + 1j * lambda_ops[\"L2\"]).real, # T+\n", " \"t\": (lambda_ops[\"L1\"] - 1j * lambda_ops[\"L2\"]).real, # T-\n", " \"V\": (lambda_ops[\"L4\"] + 1j * lambda_ops[\"L5\"]).real, # V+\n", " \"v\": (lambda_ops[\"L4\"] - 1j * lambda_ops[\"L5\"]).real, # V-\n", " \"U\": (lambda_ops[\"L6\"] + 1j * lambda_ops[\"L7\"]).real, # U+\n", " \"u\": (lambda_ops[\"L6\"] - 1j * lambda_ops[\"L7\"]).real, # U-\n", " \"L\": lambda_ops[\"L3\"], # L3\n", " \"l\": lambda_ops[\"L8\"], # L8\n", " }\n", " site_basis.append(basis)\n", " site_ops.append(ops)\n", "\n", "# [Part B] Set Hamiltonian terms\n", "driver.initialize_system(n_sites=L, vacuum=Q(0, 0), target=Q(0, 0), hamil_init=False)\n", "driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)\n", "b = driver.expr_builder()\n", "\n", "b.add_term(\"Tt\", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), 0.5 * 0.25)\n", "b.add_term(\"Vv\", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), 0.5 * 0.25)\n", "b.add_term(\"Uu\", np.array([[i, i + 1, i + 1, i] for i in range(L - 1)]).ravel(), 0.5 * 0.25)\n", "b.add_term(\"LL\", np.array([[i, i + 1] for i in range(L - 1)]).ravel(), 0.25)\n", "b.add_term(\"ll\", np.array([[i, i + 1] for i in range(L - 1)]).ravel(), 0.25)\n", "b.iscale(1 / L) # compute energy per site instead of total energy\n", "\n", "# [Part C] Perform DMRG\n", "mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops=\"\"), algo_type=MPOAlgorithmTypes.FastBipartite)\n", "mps = driver.get_random_mps(tag=\"KET\", bond_dim=250, nroots=1)\n", "\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 (per site) = %20.15f\" % energy)" ] }, { "cell_type": "markdown", "metadata": { "id": "xomn3CZSB7gq" }, "source": [ "## SU(2) $t-J$ Model\n", "\n", "As an example of non-Abelian custom symmetry groups, in the following we will find the ground state of the 2D SU(2) $t-J$ model using the custom symemtry group syntax. We will use the $\\mathrm{U(1) \\times SU(2)}$ symmetry group (for particle number $N$ and total spin $S$) for this problem.\n", "\n", "The $t-J$ model Hamiltonian is given by\n", "\n", "$$\n", "\\hat{H} = -t \\sum_{\\langle ij\\rangle, \\sigma} \\Big( a^\\dagger_{i\\sigma}\n", " a_{j\\sigma} + a^\\dagger_{j\\sigma} a_{i\\sigma} \\Big)\n", " + J \\sum_{\\langle ij\\rangle} \\Big( \\mathbf{\\hat{S}}_i \\cdot \\mathbf{\\hat{S}}_j - \\frac{1}{4} n_i n_j \\Big)\n", "$$\n", "\n", "We consider $4\\times 4$ square 2D lattice with open boundary condition, $t = 1, J = 0.4$, and hole doping $=1/8$. $\\langle ij\\rangle$ only conatins one of $i=1, j=2$ and $i=2, j=1$, for example. The reference energy per site can be found in *npj Quantum Materials* **5**, 28 (2020).\n", "\n", "Written in SU(2) notation, the Hamiltonian is\n", "\n", "$$\n", "\\hat{H} = -\\sqrt{2} t \\sum_{\\langle ij\\rangle} \\Big(\n", " \\big(a^\\dagger_i\\big)^{[1/2]} \\otimes_{[0]} \\big(a_j\\big)^{[1/2]} +\n", " \\big(a^\\dagger_j\\big)^{[1/2]} \\otimes_{[0]} \\big(a_i\\big)^{[1/2]} \\Big) \\\\\n", " + J \\sum_{\\langle ij\\rangle} \\Big[ -\\frac{\\sqrt{3}}{2}\n", " \\Big( \\big(a^\\dagger_i\\big)^{[1/2]} \\otimes_{[1]} \\big(a_i\\big)^{[1/2]}\n", " \\Big)\\otimes_{[0]} \\Big(\n", " \\big(a^\\dagger_j\\big)^{[1/2]} \\otimes_{[1]} \\big(a_j\\big)^{[1/2]}\n", " \\Big) \\\\\n", " - \\frac{(\\sqrt{2})^2}{4}\n", " \\Big( \\big(a^\\dagger_i\\big)^{[1/2]} \\otimes_{[0]} \\big(a_i\\big)^{[1/2]}\n", " \\Big)\\otimes_{[0]} \\Big(\n", " \\big(a^\\dagger_j\\big)^{[1/2]} \\otimes_{[0]} \\big(a_j\\big)^{[1/2]}\n", " \\Big)\n", " \\Big]\n", "$$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "bH-6fp6dGij1", "outputId": "53f493a3-7493-4f21-e65f-cbc131851f78" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Build MPO | Nsites = 16 | Nterms = 96 | Algorithm = FastBIP | Cutoff = 1.00e-14\n", " Site = 0 / 16 .. Mmpo = 5 DW = 0.00e+00 NNZ = 5 SPT = 0.0000 Tmvc = 0.000 T = 0.032\n", " Site = 1 / 16 .. Mmpo = 10 DW = 0.00e+00 NNZ = 13 SPT = 0.7400 Tmvc = 0.000 T = 0.021\n", " Site = 2 / 16 .. Mmpo = 14 DW = 0.00e+00 NNZ = 18 SPT = 0.8714 Tmvc = 0.000 T = 0.041\n", " Site = 3 / 16 .. Mmpo = 18 DW = 0.00e+00 NNZ = 22 SPT = 0.9127 Tmvc = 0.000 T = 0.012\n", " Site = 4 / 16 .. Mmpo = 18 DW = 0.00e+00 NNZ = 22 SPT = 0.9321 Tmvc = 0.000 T = 0.009\n", " Site = 5 / 16 .. Mmpo = 18 DW = 0.00e+00 NNZ = 26 SPT = 0.9198 Tmvc = 0.000 T = 0.012\n", " Site = 6 / 16 .. Mmpo = 18 DW = 0.00e+00 NNZ = 26 SPT = 0.9198 Tmvc = 0.000 T = 0.011\n", " Site = 7 / 16 .. Mmpo = 18 DW = 0.00e+00 NNZ = 26 SPT = 0.9198 Tmvc = 0.000 T = 0.012\n", " Site = 8 / 16 .. Mmpo = 18 DW = 0.00e+00 NNZ = 22 SPT = 0.9321 Tmvc = 0.000 T = 0.015\n", " Site = 9 / 16 .. Mmpo = 18 DW = 0.00e+00 NNZ = 26 SPT = 0.9198 Tmvc = 0.000 T = 0.013\n", " Site = 10 / 16 .. Mmpo = 18 DW = 0.00e+00 NNZ = 26 SPT = 0.9198 Tmvc = 0.000 T = 0.013\n", " Site = 11 / 16 .. Mmpo = 18 DW = 0.00e+00 NNZ = 26 SPT = 0.9198 Tmvc = 0.000 T = 0.012\n", " Site = 12 / 16 .. Mmpo = 14 DW = 0.00e+00 NNZ = 22 SPT = 0.9127 Tmvc = 0.000 T = 0.021\n", " Site = 13 / 16 .. Mmpo = 10 DW = 0.00e+00 NNZ = 18 SPT = 0.8714 Tmvc = 0.000 T = 0.018\n", " Site = 14 / 16 .. Mmpo = 5 DW = 0.00e+00 NNZ = 13 SPT = 0.7400 Tmvc = 0.000 T = 0.011\n", " Site = 15 / 16 .. Mmpo = 1 DW = 0.00e+00 NNZ = 5 SPT = 0.0000 Tmvc = 0.000 T = 0.008\n", "Ttotal = 0.263 Tmvc-total = 0.001 MPO bond dimension = 18 MaxDW = 0.00e+00\n", "NNZ = 316 SIZE = 3486 SPT = 0.9094\n", "\n", "Rank = 0 Ttotal = 0.480 MPO method = FastBipartite bond dimension = 18 NNZ = 316 SIZE = 3486 SPT = 0.9094\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 3.850 | E = -8.9796008520 | DW = 5.86154e-07\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 7.778 | E = -9.0284250184 | DE = -4.88e-02 | DW = 1.36125e-06\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 9.742 | E = -9.0298619274 | DE = -1.44e-03 | DW = 1.23409e-06\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 10.604 | E = -9.0298670256 | DE = -5.10e-06 | DW = 1.13415e-06\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 11.898 | E = -9.0298686867 | DE = -1.66e-06 | DW = 2.27179e-10\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 12.869 | E = -9.0298686872 | DE = -4.35e-10 | DW = 1.12668e-10\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 13.782 | E = -9.0298686872 | DE = -8.19e-12 | DW = 2.27110e-10\n", "\n", "Sweep = 7 | Direction = backward | Bond dimension = 500 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 14.710 | E = -9.0298686872 | DE = -1.02e-11 | DW = 1.12643e-10\n", "\n", "Sweep = 8 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-09\n", "Time elapsed = 15.473 | E = -9.0298686872 | DE = 2.91e-11 | DW = 1.78156e-19\n", "\n", "DMRG energy = -9.029868687160421 (per site = -0.564367)\n" ] } ], "source": [ "from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes\n", "import numpy as np\n", "\n", "LX, LY = 4, 4\n", "L = LX * LY\n", "J = 0.4\n", "N_ELEC = 14 # 1/8 doping\n", "TWO_S = 0\n", "\n", "driver = DMRGDriver(scratch=\"./tmp\", symm_type=SymmetryTypes.SAnySU2, n_threads=4)\n", "\n", "driver.set_symmetry_groups(\"U1Fermi\", \"SU2\", \"SU2\")\n", "Q = driver.bw.SX\n", "\n", "# [Part A] Set states and matrix representation of operators in local Hilbert space\n", "site_basis, site_ops = [], []\n", "\n", "for k in range(L):\n", " basis = [(Q(0, 0, 0), 1), (Q(1, 1, 1), 1)] # [01]\n", " ops = {\n", " \"\": np.array([[1, 0], [0, 1]]), # identity\n", " \"C\": np.array([[0, 0], [1, 0]]), # a+\n", " \"D\": np.array([[0, 2**0.5], [0, 0]]), # a\n", " }\n", " site_basis.append(basis)\n", " site_ops.append(ops)\n", "\n", "# [Part B] Set Hamiltonian terms\n", "driver.initialize_system(n_sites=L, vacuum=Q(0, 0, 0), target=Q(N_ELEC, TWO_S, TWO_S), hamil_init=False)\n", "driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)\n", "b = driver.expr_builder()\n", "\n", "f = lambda i, j: i * LY + j if i % 2 == 0 else i * LY + LY - 1 - j\n", "\n", "for i in range(0, LX):\n", " for j in range(0, LY):\n", " if i + 1 < LX:\n", " b.add_term(\"(C+D)0\", [f(i, j), f(i + 1, j), f(i + 1, j), f(i, j)], -(2 ** 0.5))\n", " b.add_term(\"((C+D)2+(C+D)2)0\", [f(i, j), f(i, j), f(i + 1, j), f(i + 1, j)], J * -(3 ** 0.5) / 2)\n", " b.add_term(\"((C+D)0+(C+D)0)0\", [f(i, j), f(i, j), f(i + 1, j), f(i + 1, j)], J * -1 / 2)\n", " if j + 1 < LY:\n", " b.add_term(\"(C+D)0\", [f(i, j), f(i, j + 1), f(i, j + 1), f(i, j)], -(2 ** 0.5))\n", " b.add_term(\"((C+D)2+(C+D)2)0\", [f(i, j), f(i, j), f(i, j + 1), f(i, j + 1)], J * -(3 ** 0.5) / 2)\n", " b.add_term(\"((C+D)0+(C+D)0)0\", [f(i, j), f(i, j), f(i, j + 1), f(i, j + 1)], J * -1 / 2)\n", "\n", "# [Part C] Perform DMRG\n", "mpo = driver.get_mpo(b.finalize(adjust_order=True), algo_type=MPOAlgorithmTypes.FastBipartite, iprint=1)\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 (per site = %10.6f)\" % (energy, energy / L))" ] } ], "metadata": { "colab": { "provenance": [] }, "kernelspec": { "display_name": "Python 3", "name": "python3" }, "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 0 }