{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "owC5h_4Zs9le" }, "source": [ "# Vibrational 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/vibrational-hamiltonians.ipynb)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "id": "7rG2UFfWl74_" }, "outputs": [], "source": [ "!pip install block2==0.5.3rc16 -qq --progress-bar off --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/" ] }, { "cell_type": "markdown", "metadata": { "id": "Mx_FmbJamLjq" }, "source": [ "In this tutorial, we explain how to perform DMRG for vibrational Hamiltonians to find vibrational energy levels and vibrational mode assignments. We use the Watson Hamiltonian and expand the Potential Energy Surface (PES) as a sextic force field in terms of normal mode coordinates about the equilibrium geometry:\n", "\n", "$$\n", "\\hat{H} = -\\frac{1}{2}\\sum_{i=1} \\frac{\\partial^2}{\\partial Q_i^2}\n", " + \\frac{1}{2!} \\sum_{ij} F_{ij} Q_i Q_j\n", " + \\frac{1}{3!} \\sum_{ijk} F_{ijk} Q_i Q_j Q_k\n", " + \\frac{1}{4!} \\sum_{ijkl} F_{ijkl} Q_i Q_j Q_k Q_l \\\\\n", " + \\frac{1}{5!} \\sum_{ijklm} F_{ijklm} Q_i Q_j Q_k Q_l Q_m\n", " + \\frac{1}{6!} \\sum_{ijklmn} F_{ijklmn} Q_i Q_j Q_k Q_l Q_m Q_n\n", "$$\n", "\n", "For details of the definition of the Hamiltonian and the wavefunction ansatz, see [M. Sibaev, D. L. Crittenden. PyVCI: A flexible open-source code for calculating accurate molecular infrared spectra. *Computer Physics Communications*, **203**, 290-297 (2016)](https://doi.org/10.1016/j.cpc.2016.02.026)." ] }, { "cell_type": "markdown", "metadata": { "id": "4i0vUf6zodRR" }, "source": [ "## Force Field Parameters\n", "\n", "First we need to load the force field parameters in the Hamiltonian. We consider $\\mathrm{N_2H_2}$ as an example, where the PES can be found in the PyPES package (https://github.com/dlc62/pypes-lib) and the reference energy levels computed using Vibrational Configuration Interaction (VCI) are given in [M. Sibaev, D. L. Crittenden. The PyPES library of high quality semi‐global potential energy surfaces.\" *Journal of Computational Chemistry* **29**, 2200-2207 (2015)](https://doi.org/10.1002/jcc.24192).\n", "\n", "We first download the force field file ``force_field.tar`` (without decompressing it). It will be automaticaly decompressed in the following Python script." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "jezxBzi5oc7X", "outputId": "099d2081-fd1e-4a9e-9d53-bc53694099ca" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--2024-08-01 02:09:02-- https://github.com/dlc62/pypes-lib/raw/master/force_field.tar\n", "Resolving github.com (github.com)... 140.82.113.4\n", "Connecting to github.com (github.com)|140.82.113.4|:443... connected.\n", "HTTP request sent, awaiting response... 302 Found\n", "Location: https://raw.githubusercontent.com/dlc62/pypes-lib/master/force_field.tar [following]\n", "--2024-08-01 02:09:03-- https://raw.githubusercontent.com/dlc62/pypes-lib/master/force_field.tar\n", "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...\n", "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.\n", "HTTP request sent, awaiting response... 200 OK\n", "Length: 21156352 (20M) [application/octet-stream]\n", "Saving to: ‘force_field.tar’\n", "\n", "force_field.tar 100%[===================>] 20.18M --.-KB/s in 0.1s \n", "\n", "2024-08-01 02:09:03 (173 MB/s) - ‘force_field.tar’ saved [21156352/21156352]\n", "\n" ] } ], "source": [ "!wget https://github.com/dlc62/pypes-lib/raw/master/force_field.tar" ] }, { "cell_type": "markdown", "metadata": { "id": "xt-zRdWyryDD" }, "source": [ "## Perform Vibrational DMRG\n", "\n", "The following script will solve the 10 low-energy roots for the vibrational Hamiltonian of the $\\mathrm{N_2H_2}$ molecule (rotational effects are not included)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "4wN92_7eoYy5", "outputId": "e0541a6c-a76b-44e1-9dc8-033f4c5b20be" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PES: N2H2\n", " trans-N2H2 molecule. Saved in normal mode coordinates from PESN 2.\n", " Reference:\n", " Martin, J. M. L.; Taylor, P. R. Spectrochim. Acta Part A: Mol. Biomol. Spectrosc. 1997, 53(8), 1039-1050.\n", "\n", "Harmonic Frequencies (cm-1): 1328.4271 1350.2811 1558.4338 1621.8428 3269.8469 3301.5852\n", "Harmonic ZPE (cm-1): 6215.2084\n", "\n", "Build MPO | Nsites = 6 | Nterms = 348 | Algorithm = FastBIP | Cutoff = 1.00e-75\n", " Site = 0 / 6 .. Mmpo = 4 DW = 0.00e+00 NNZ = 4 SPT = 0.0000 Tmvc = 0.000 T = 0.014\n", " Site = 1 / 6 .. Mmpo = 11 DW = 0.00e+00 NNZ = 16 SPT = 0.6364 Tmvc = 0.000 T = 0.027\n", " Site = 2 / 6 .. Mmpo = 23 DW = 0.00e+00 NNZ = 66 SPT = 0.7391 Tmvc = 0.000 T = 0.034\n", " Site = 3 / 6 .. Mmpo = 19 DW = 0.00e+00 NNZ = 104 SPT = 0.7620 Tmvc = 0.000 T = 0.013\n", " Site = 4 / 6 .. Mmpo = 7 DW = 0.00e+00 NNZ = 36 SPT = 0.7293 Tmvc = 0.000 T = 0.016\n", " Site = 5 / 6 .. Mmpo = 1 DW = 0.00e+00 NNZ = 7 SPT = 0.0000 Tmvc = 0.000 T = 0.007\n", "Ttotal = 0.112 Tmvc-total = 0.001 MPO bond dimension = 23 MaxDW = 0.00e+00\n", "NNZ = 233 SIZE = 878 SPT = 0.7346\n", "\n", "Rank = 0 Ttotal = 0.196 MPO method = FastBipartite bond dimension = 23 NNZ = 233 SIZE = 878 SPT = 0.7346\n", "\n", "==== ROOT 0 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.162 | E = 0.0278175746 | DW = 1.50629e-20\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.230 | E = 0.0278175746 | DE = -4.86e-17 | DW = 2.94238e-20\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.303 | E = 0.0278175746 | DE = -4.54e-15 | DW = 9.66411e-21\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.406 | E = 0.0278175746 | DE = 4.55e-15 | DW = 1.94996e-20\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.591 | E = 0.0278175746 | DE = -4.53e-15 | DW = 1.11346e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.681 | E = 0.0278175746 | DE = 4.52e-15 | DW = 5.49483e-21\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.789 | E = 0.0278175746 | DE = -4.55e-15 | DW = 1.43430e-20\n", "\n", "\n", "==== ROOT 1 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.510 | E = 0.0337897813 | DW = 8.91046e-16\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.609 | E = 0.0337897813 | DE = -5.29e-12 | DW = 7.20577e-18\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.717 | E = 0.0337897813 | DE = 0.00e+00 | DW = 1.43698e-20\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.819 | E = 0.0337897813 | DE = -2.08e-17 | DW = 3.82610e-18\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.923 | E = 0.0337897813 | DE = 1.39e-17 | DW = 1.50304e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.030 | E = 0.0337897813 | DE = 9.71e-17 | DW = 2.46714e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.124 | E = 0.0337897813 | DE = 1.39e-17 | DW = 1.36100e-20\n", "\n", "\n", "==== ROOT 2 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.372 | E = 0.0336294440 | DW = 6.32649e-12\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.482 | E = 0.0336294440 | DE = -3.26e-12 | DW = 1.32102e-17\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.612 | E = 0.0336294440 | DE = 0.00e+00 | DW = 8.63491e-21\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.718 | E = 0.0336294440 | DE = -6.94e-18 | DW = 1.09199e-17\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.840 | E = 0.0336294440 | DE = -6.94e-18 | DW = 1.24511e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.955 | E = 0.0336294440 | DE = 6.94e-18 | DW = 1.33919e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.074 | E = 0.0336294440 | DE = 1.39e-17 | DW = 9.75384e-21\n", "\n", "\n", "==== ROOT 3 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.336 | E = 0.0347691843 | DW = 2.71296e-15\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.453 | E = 0.0347691843 | DE = -9.19e-13 | DW = 4.24397e-20\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.546 | E = 0.0347691843 | DE = -6.94e-18 | DW = 1.57276e-20\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.633 | E = 0.0347691843 | DE = -3.47e-17 | DW = 2.62031e-20\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.740 | E = 0.0347691843 | DE = 0.00e+00 | DW = 1.64701e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.824 | E = 0.0347691843 | DE = 3.47e-17 | DW = 1.24025e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.913 | E = 0.0347691843 | DE = -1.39e-17 | DW = 1.95745e-20\n", "\n", "\n", "==== ROOT 4 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.221 | E = 0.0350118953 | DW = 2.73902e-17\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.320 | E = 0.0350118953 | DE = -1.27e-12 | DW = 2.21093e-19\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.426 | E = 0.0350118953 | DE = 0.00e+00 | DW = 8.90919e-21\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.520 | E = 0.0350118953 | DE = -6.94e-18 | DW = 2.29959e-18\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.624 | E = 0.0350118953 | DE = 0.00e+00 | DW = 1.89815e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.738 | E = 0.0350118953 | DE = 0.00e+00 | DW = 1.07261e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.846 | E = 0.0350118953 | DE = 6.94e-18 | DW = 9.05114e-21\n", "\n", "\n", "==== ROOT 5 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.313 | E = 0.0393996116 | DW = 2.56369e-15\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.424 | E = 0.0393996116 | DE = -5.93e-12 | DW = 4.63254e-18\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.538 | E = 0.0393996116 | DE = 0.00e+00 | DW = 1.62961e-20\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.647 | E = 0.0393996116 | DE = 4.16e-17 | DW = 1.40207e-17\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.773 | E = 0.0393996116 | DE = -6.94e-18 | DW = 1.72877e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.890 | E = 0.0393996116 | DE = -2.78e-17 | DW = 1.85333e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.998 | E = 0.0393996116 | DE = 0.00e+00 | DW = 1.03031e-20\n", "\n", "\n", "==== ROOT 6 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.511 | E = 0.0397070198 | DW = 1.26937e-13\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.628 | E = 0.0397070198 | DE = -9.73e-12 | DW = 1.34647e-17\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.770 | E = 0.0397070198 | DE = 1.39e-17 | DW = 1.20041e-20\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.884 | E = 0.0397070198 | DE = -6.94e-18 | DW = 1.12444e-17\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.003 | E = 0.0397070198 | DE = 6.94e-18 | DW = 1.36156e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.124 | E = 0.0397070198 | DE = 0.00e+00 | DW = 1.99115e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.240 | E = 0.0397070198 | DE = -6.94e-18 | DW = 1.60524e-20\n", "\n", "\n", "==== ROOT 7 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.430 | E = 0.0398511574 | DW = 3.70941e-15\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.561 | E = 0.0398511574 | DE = -4.99e-12 | DW = 2.45878e-20\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.691 | E = 0.0398511574 | DE = -1.39e-17 | DW = 2.44316e-20\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.823 | E = 0.0398511574 | DE = -1.39e-17 | DW = 3.23054e-19\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.957 | E = 0.0398511574 | DE = -6.94e-18 | DW = 2.15729e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.083 | E = 0.0398511574 | DE = 3.47e-17 | DW = 1.74024e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.218 | E = 0.0398511574 | DE = 0.00e+00 | DW = 1.91998e-20\n", "\n", "\n", "==== ROOT 8 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.476 | E = 0.0406183665 | DW = 2.45735e-15\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.607 | E = 0.0406183665 | DE = -5.63e-12 | DW = 2.86690e-17\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.754 | E = 0.0406183665 | DE = 6.94e-18 | DW = 1.81005e-20\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.888 | E = 0.0406183665 | DE = -6.94e-18 | DW = 2.33180e-17\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.028 | E = 0.0406183665 | DE = 2.08e-17 | DW = 2.37587e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.194 | E = 0.0406183665 | DE = -9.02e-17 | DW = 2.19444e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.340 | E = 0.0406183665 | DE = 1.39e-17 | DW = 2.59923e-20\n", "\n", "\n", "==== ROOT 9 ====\n", "\n", "Sweep = 0 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.588 | E = 0.0395978303 | DW = 5.53041e-11\n", "\n", "Sweep = 1 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.748 | E = 0.0395978303 | DE = -8.83e-12 | DW = 4.70079e-20\n", "\n", "Sweep = 2 | Direction = forward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 0.903 | E = 0.0395978303 | DE = 6.94e-18 | DW = 1.44256e-20\n", "\n", "Sweep = 3 | Direction = backward | Bond dimension = 50 | Noise = 1.00e-04 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.048 | E = 0.0395978303 | DE = 6.94e-18 | DW = 3.29833e-20\n", "\n", "Sweep = 4 | Direction = forward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.203 | E = 0.0395978303 | DE = -1.39e-17 | DW = 1.81447e-20\n", "\n", "Sweep = 5 | Direction = backward | Bond dimension = 100 | Noise = 1.00e-05 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.353 | E = 0.0395978303 | DE = 0.00e+00 | DW = 1.83984e-20\n", "\n", "Sweep = 6 | Direction = forward | Bond dimension = 100 | Noise = 0.00e+00 | Dav threshold = 1.00e-10\n", "Time elapsed = 1.505 | E = 0.0395978303 | DE = 0.00e+00 | DW = 1.61162e-20\n", "\n", "ZPE = 6105.2523 cm-1 Max weight = 0.9705 Assignment = [0 0 0 0 0 0]\n", "VIB = 1310.7479 cm-1 Max weight = 0.9699 Assignment = [0 1 0 0 0 0]\n", "VIB = 1275.5580 cm-1 Max weight = 0.9732 Assignment = [1 0 0 0 0 0]\n", "VIB = 1525.7021 cm-1 Max weight = 0.9334 Assignment = [0 0 1 0 0 0]\n", "VIB = 1578.9710 cm-1 Max weight = 0.9524 Assignment = [0 0 0 1 0 0]\n", "VIB = 2541.9634 cm-1 Max weight = 0.9322 Assignment = [2 0 0 0 0 0]\n", "VIB = 2609.4317 cm-1 Max weight = 0.9226 Assignment = [0 2 0 0 0 0]\n", "VIB = 2641.0663 cm-1 Max weight = 0.8964 Assignment = [0 3 0 0 0 0]\n", "VIB = 2809.4492 cm-1 Max weight = 0.7351 Assignment = [0 1 1 0 0 0]\n", "VIB = 2585.4674 cm-1 Max weight = 0.9687 Assignment = [1 1 0 0 0 0]\n" ] } ], "source": [ "import numpy as np\n", "import math\n", "import tarfile, bz2\n", "\n", "pes_name = 'N2H2'\n", "pes_n = None\n", "nroots = 10\n", "nexcit = 4 # max excitation number per mode\n", "\n", "# Read force field parameters from PyPES\n", "with tarfile.open(name='force_field.tar') as tf:\n", " fn = [fn for fn in tf.getnames() if fn.endswith('/molecule_%s.py.bz2' % pes_name)][0]\n", " pess, ipes, ff = {}, -1, None\n", " for l in bz2.decompress(tf.extractfile(fn).read()).decode('utf-8').split('\\n'):\n", " if 'if state.pesn ==' in l:\n", " if ff is not None:\n", " pess[ipes]['ff'] = ff\n", " ipes = int(l.split('\\\"')[-2])\n", " pess[ipes], ff = {'nm': False, 'cmt': ''}, None\n", " elif 'state.normal_mode_ff' in l:\n", " pess[ipes]['nm'] = True\n", " elif 'state.ff_comment.append' in l:\n", " pess[ipes]['cmt'] += l.split('\\'')[1] + '\\n'\n", " elif 'ff.ff =' in l:\n", " ff = []\n", " elif 'else:' in l:\n", " if ff is not None:\n", " pess[ipes]['ff'] = ff\n", " ff = None\n", " elif ff is not None:\n", " xs = l.split('[')[1].split(']')[0].split(',')\n", " ff.append([int(x) for x in xs[:-1]] + [float(xs[-1])])\n", " pes = pess[pes_n if pes_n is not None else ipes]\n", " assert pes['nm']\n", "\n", "derivs = [{tuple(l[:-1]): l[-1] for l in pes['ff'] if len(l) == i + 1} for i in range(max([len(l) for l in pes['ff']]))]\n", "n_sites = max(derivs[2])[0] + 1\n", "omega = np.array([derivs[2][(i, i)] for i in range(n_sites)]) ** 0.5\n", "au2cm = 219474.6435\n", "\n", "print('PES:', pes_name)\n", "print(pes['cmt'])\n", "print('Harmonic Frequencies (cm-1): %s' % ''.join(['%12.4f' % x for x in (omega * au2cm)]))\n", "print('Harmonic ZPE (cm-1): %12.4f' % np.sum(omega * au2cm / 2))\n", "\n", "from pyblock2.driver.core import DMRGDriver, SymmetryTypes\n", "driver = DMRGDriver(scratch=\"./tmp\", symm_type=SymmetryTypes.SAny, n_threads=4)\n", "driver.set_symmetry_groups()\n", "q = driver.bw.SX()\n", "\n", "# Table 1 in Computer Physics Communications 203, 290-297 (2016).\n", "site_basis, site_ops = [], []\n", "for k in range(n_sites):\n", " basis = [(q, nexcit)]\n", " pmat, qmat = np.zeros((nexcit, nexcit)), np.zeros((nexcit, nexcit))\n", " for i, j, n, dn in [(i, j, max(i, j), abs(i - j)) for i in range(nexcit) for j in range(nexcit)]:\n", " if dn == 0:\n", " pmat[i, j] = -omega[k] * (i + 0.5)\n", " elif dn == 2:\n", " pmat[i, j] = 0.5 * omega[k] * (n * (n - 1)) ** 0.5\n", " elif dn == 1:\n", " qmat[i, j] = (n / (2 * omega[k])) ** 0.5\n", " ops = {\"\": np.identity(nexcit), \"P\": pmat, \"Q\": qmat}\n", " site_basis.append(basis)\n", " site_ops.append(ops)\n", "\n", "# Set up vibrational Hamiltonian\n", "driver.initialize_system(n_sites=n_sites, vacuum=q, target=q, hamil_init=False)\n", "driver.ghamil = driver.get_custom_hamiltonian(site_basis, site_ops)\n", "b = driver.expr_builder()\n", "\n", "b.add_term(\"P\", list(range(n_sites)), -0.5)\n", "for p, idx, f in [(p, idx, f) for p, d in enumerate(derivs) for idx, f in d.items()]:\n", " fac = np.prod([math.factorial(g) for g in np.unique(idx, return_counts=True)[1]])\n", " b.add_term(\"Q\" * p, idx, f / fac)\n", "\n", "mpo = driver.get_mpo(b.finalize(adjust_order=True, fermionic_ops=''), cutoff=1E-75, iprint=1)\n", "\n", "# Run DMRG\n", "kets = [None for _ in range(nroots)]\n", "energies = []\n", "for ir in range(nroots):\n", " print('\\n==== ROOT %4d ====' % ir)\n", " kets[ir] = driver.get_random_mps(tag=\"KET-%d\" % ir, bond_dim=50, nroots=1)\n", " energies.append(driver.dmrg(mpo, kets[ir], n_sweeps=10, bond_dims=[50] * 4 + [100] * 4,\n", " noises=[1e-4] * 4 + [1e-5] * 2 + [0], thrds=[1e-10] * 8, dav_max_iter=300,\n", " proj_weights=[1.0] * ir, proj_mpss=kets[:ir], iprint=1, tol=1E-11))\n", "\n", "# Energy levels and assignment\n", "for ir in range(nroots):\n", " if ir == 0:\n", " print('ZPE = %12.4f cm-1' % (energies[ir] * au2cm), end=' ', flush=True)\n", " else:\n", " print('VIB = %12.4f cm-1' % ((energies[ir] - energies[0]) * au2cm), end=' ', flush=True)\n", " kets[ir] = driver.adjust_mps(kets[ir], dot=1)[0]\n", " csfs, coeffs = driver.get_csf_coefficients(kets[ir], cutoff=0.1, iprint=0)\n", " csf, weight = csfs[np.argmax(coeffs ** 2)], coeffs[np.argmax(coeffs ** 2)] ** 2\n", " print(' Max weight = %10.4f Assignment = %s' % (weight, csf))" ] } ], "metadata": { "colab": { "provenance": [] }, "kernelspec": { "display_name": "Python 3", "name": "python3" }, "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 0 }