diff --git a/.gitignore b/.gitignore index 8ee19dc..c9ace05 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,72 @@ -*.pyc -__pycache__/ .DS_Store -.ipynb_checkpoints/ *lastfailed *coverage + +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +*.egg-info/ +.installed.cfg +*.egg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +.hypothesis/ + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# pyenv +.python-version + +# Environments +.env +.venv +env/ +venv/ +ENV/ + +# mkdocs documentation +/site diff --git a/README.rst b/README.rst index 43d011b..607667d 100644 --- a/README.rst +++ b/README.rst @@ -14,13 +14,13 @@ FermiLib - An open source software for analyzing quantum simulation algorithms FermiLib is an open source effort for analyzing quantum simulation algorithms. -The first version (v0.1a0) is an alpha release which features data structures and tools for obtaining and manipulating representations of fermionic Hamiltonians. FermiLib is designed as a library on top of ProjectQ and leverages ProjectQ to compile, emulate and simulate quantum circuits. +The first version (v0.1a0) is an alpha release which features data structures and tools for obtaining and manipulating representations of fermionic Hamiltonians. FermiLib is designed as a library on top of `ProjectQ `__ and leverages ProjectQ to compile, emulate and simulate quantum circuits. There are also `plugins `__ available for FermiLib. Getting started --------------- To start using FermiLib, simply follow the installation instructions in the `intro `__. There, you will also find `code examples `__. Also, make sure to check out the `ProjectQ -website `__ and the detailed `code documentation `__. +website `__ and the detailed `code documentation `__. Moreover, take a look at the available `plugins `__ for FermiLib. How to contribute ----------------- diff --git a/examples/fermilib_demo.ipynb b/examples/fermilib_demo.ipynb index a3867e8..cfc8313 100644 --- a/examples/fermilib_demo.ipynb +++ b/examples/fermilib_demo.ipynb @@ -2,7 +2,10 @@ "cells": [ { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "# Introduction to FermiLib alpha version 0.1a0\n", "### Ryan Babbush, Jarrod McClean, Damian Steiger, Ian Kivlichan, Thomas Haner, Vojtech Havlicek, Matthew Neeley and Wei Sun\n", @@ -11,42 +14,60 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "## Initializing the FermionOperator data structure\n", "\n", "Fermionic systems are often treated in second quantization where arbitrary operators can be expressed using the fermionic creation and annihilation operators, $a^\\dagger_k$ and $a_k$. The fermionic ladder operators play a similar role to their qubit ladder operator counterparts, $\\sigma^+_k$ and $\\sigma^-_k$ but are distinguished by the cannonical fermionic anticommutation relations, $\\{a^\\dagger_i, a^\\dagger_j\\} = \\{a_i, a_j\\} = 0$ and $\\{a_i, a_j^\\dagger\\} = \\delta_{ij}$. Any weighted sums of products of these operators are represented with the FermionOperator data structure in FermiLib. The following are examples of valid FermionOperators:\n", + "\n", + "$$\n", "\\begin{align}\n", - "& a_1 \\\\\n", - "& 1.7 a^\\dagger_3\\\\\n", - "&-1.7 \\, a^\\dagger_3 a_1\\\\\n", - "&(1 + 2i) \\, a^\\dagger_4 a^\\dagger_3 a_9 a_1\\\\\n", - "&(1 + 2i) \\, a^\\dagger_4 a^\\dagger_3 a_9 a_1 - 1.7 \\, a^\\dagger_3 a_1\n", + "& a_1 \\nonumber \\\\\n", + "& 1.7 a^\\dagger_3 \\nonumber \\\\\n", + "&-1.7 \\, a^\\dagger_3 a_1 \\nonumber \\\\\n", + "&(1 + 2i) \\, a^\\dagger_4 a^\\dagger_3 a_9 a_1 \\nonumber \\\\\n", + "&(1 + 2i) \\, a^\\dagger_4 a^\\dagger_3 a_9 a_1 - 1.7 \\, a^\\dagger_3 a_1 \\nonumber\n", "\\end{align}\n", + "$$\n", "\n", "The FermionOperator class is contained in $\\textrm{ops/_fermion_operators.py}$. In order to support fast addition of FermionOperator instances, the class is implemented as hash table (python dictionary). The keys of the dictionary encode the strings of ladder operators and values of the dictionary store the coefficients. The strings of ladder operators are encoded as a tuple of 2-tuples which we refer to as the \"terms tuple\". Each ladder operator is represented by a 2-tuple. The first element of the 2-tuple is an int indictating the tensor factor on which the ladder operator acts. The second element of the 2-tuple is Boole: 1 represents raising and 0 represents lowering. For instance, $a^\\dagger_8$ is represented in a 2-tuple as $(8, 1)$. Note that indices start at 0 and the identity operator is an empty list. Below we give some examples of operators and their terms tuple:\n", + "\n", + "$$\n", "\\begin{align}\n", - "I & \\mapsto ()\\\\\n", - "a_1 & \\mapsto ((1, 0),)\\\\\n", - "a^\\dagger_3 & \\mapsto ((3, 1),)\\\\\n", - "a^\\dagger_3 a_1 & \\mapsto ((3, 1), (1, 0))\\\\\n", - "a^\\dagger_4 a^\\dagger_3 a_9 a_1 & \\mapsto ((4, 1), (3, 1), (9, 0), (1, 0))\n", + "I & \\mapsto () \\nonumber \\\\\n", + "a_1 & \\mapsto ((1, 0),) \\nonumber \\\\\n", + "a^\\dagger_3 & \\mapsto ((3, 1),) \\nonumber \\\\\n", + "a^\\dagger_3 a_1 & \\mapsto ((3, 1), (1, 0)) \\nonumber \\\\\n", + "a^\\dagger_4 a^\\dagger_3 a_9 a_1 & \\mapsto ((4, 1), (3, 1), (9, 0), (1, 0)) \\nonumber\n", "\\end{align}\n", + "$$\n", + "\n", "Note that when initializing a single ladder operator one should be careful to add the comma after the inner pair. This is because in python ((1, 2)) = (1, 2) whereas ((1, 2),) = ((1, 2),). The \"terms tuple\" is usually convenient when one wishes to initialize a term as part of a coded routine. However, the terms tuple is not particularly intuitive. Accordingly, FermiLib also supports another user-friendly, string notation below. This representation is rendered when calling \"print\" on a FermionOperator.\n", + "\n", + "$$\n", "\\begin{align}\n", - "I & \\mapsto \"\"\\\\\n", - "a_1 & \\mapsto \"1\"\\\\\n", - "a^\\dagger_3 & \\mapsto \"3\\textrm{^}\"\\\\\n", - "a^\\dagger_3 a_1 & \\mapsto \"3\\textrm{^} \\,\\, 1\"\\\\\n", - "a^\\dagger_4 a^\\dagger_3 a_9 a_1 & \\mapsto \"4\\textrm{^}\\,\\, 3\\textrm{^} \\,\\, 9 \\,\\,1\"\n", + "I & \\mapsto \\textrm{\"\"} \\nonumber \\\\\n", + "a_1 & \\mapsto \\textrm{\"1\"} \\nonumber \\\\\n", + "a^\\dagger_3 & \\mapsto \\textrm{\"3^\"} \\nonumber \\\\\n", + "a^\\dagger_3 a_1 & \\mapsto \\textrm{\"3^}\\;\\textrm{1\"} \\nonumber \\\\\n", + "a^\\dagger_4 a^\\dagger_3 a_9 a_1 & \\mapsto \\textrm{\"4^}\\;\\textrm{3^}\\;\\textrm{9}\\;\\textrm{1\"} \\nonumber\n", "\\end{align}\n", + "$$\n", + "\n", "Let's initialize our first term! We do it two different ways below." ] }, { "cell_type": "code", "execution_count": 1, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -69,7 +90,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "The preferred way to specify the coefficient in FermiLib is to provide an optional coefficient argument. If not provided, the coefficient defaults to 1. In the code below, the first method is preferred. The multiplication in the second method actually creates a copy of the term, which introduces some additional cost. All inplace operands (such as +=) modify classes whereas binary operands such as + create copies. Important caveats are that the empty tuple FermionOperator(()) and the empty string FermionOperator('') initializes identity. The empty initializer FermionOperator() initializes the zero operator." ] @@ -77,7 +101,11 @@ { "cell_type": "code", "execution_count": 2, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -106,7 +134,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "Note that FermionOperator has only one attribute: .terms. This attribute is the dictionary which stores the term tuples." ] @@ -114,7 +145,11 @@ { "cell_type": "code", "execution_count": 3, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -133,7 +168,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "## Manipulating the FermionOperator data structure\n", "So far we have explained how to initialize a single FermionOperator such as $-1.7 \\, a^\\dagger_3 a_1$. However, in general we will want to represent sums of these operators such as $(1 + 2i) \\, a^\\dagger_4 a^\\dagger_3 a_9 a_1 - 1.7 \\, a^\\dagger_3 a_1$. To do this, just add together two FermionOperators! We demonstrate below." @@ -142,7 +180,11 @@ { "cell_type": "code", "execution_count": 4, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -173,7 +215,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "The print function prints each term in the operator on a different line. Note that the line my_operator = term_1 + term_2 creates a new object, which involves a copy of term_1 and term_2. The second block of code uses the inplace method +=, which is more efficient. This is especially important when trying to construct a very large FermionOperator. FermionOperators also support a wide range of builtins including, str(), repr(), *=, *, /, /=, +, +=, -, -=, - and **. Note that instead of supporting != and ==, we have the method .isclose(), since FermionOperators involve floats. We demonstrate some of these methods below." ] @@ -181,7 +226,11 @@ { "cell_type": "code", "execution_count": 5, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -227,7 +276,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "Additionally, there are a variety of methods that act on the FermionOperator data structure. We demonstrate a small subset of those methods here." ] @@ -235,7 +287,11 @@ { "cell_type": "code", "execution_count": 6, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -279,16 +335,23 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "## The QubitOperator data structure\n", - "The QubitOperator data structure is another essential part of FermiLib. While the QubitOperator was originally developed for FermiLib, it is now part of the core ProjectQ library so that it can be interpreted by the ProjectQ compiler using the TimeEvolution gate. As the name suggests, QubitOperator is used to store qubit operators in almost exactly the same way that FermionOperator is used to store fermion operators. For instance $X_0 Z_3 Y_4$ is a QubitOperator. The internal representation of this as a terms tuple would be $((0, \"X\"), (3, \"Z\"), (4, \"Y\"))$. Note that one important difference between QubitOperator and FermionOperator is that the terms in QubitOperator are always sorted in order of tensor factor. In some cases, this enables faster manipulation. We initialize some QubitOperators below." + "The QubitOperator data structure is another essential part of FermiLib. While the QubitOperator was originally developed for FermiLib, it is now part of the core ProjectQ library so that it can be interpreted by the ProjectQ compiler using the TimeEvolution gate. As the name suggests, QubitOperator is used to store qubit operators in almost exactly the same way that FermionOperator is used to store fermion operators. For instance $X_0 Z_3 Y_4$ is a QubitOperator. The internal representation of this as a terms tuple would be $((0, \\textrm{\"X\"}), (3, \\textrm{\"Z\"}), (4, \\textrm{\"Y\"}))$. Note that one important difference between QubitOperator and FermionOperator is that the terms in QubitOperator are always sorted in order of tensor factor. In some cases, this enables faster manipulation. We initialize some QubitOperators below." ] }, { "cell_type": "code", "execution_count": 7, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -317,7 +380,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "## Jordan-Wigner and Bravyi-Kitaev\n", "FermiLib provides functions for mapping FermionOperators to QubitOperators." @@ -326,7 +392,11 @@ { "cell_type": "code", "execution_count": 8, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -372,7 +442,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "We see that despite the different representation, these operators are iso-spectral. We can also apply the Jordan-Wigner transform in reverse to map arbitrary QubitOperators to FermionOperators. Note that we also demonstrate the .compress() method (a method on both FermionOperators and QubitOperators) which removes zero entries." ] @@ -380,7 +453,11 @@ { "cell_type": "code", "execution_count": 9, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -429,7 +506,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "## Sparse matrices and the Hubbard model\n", "Often, one would like to obtain a sparse matrix representation of an operator which can be analyzed numerically. There is code in both fermilib.transforms and fermilib.utils which facilitates this. The function get_sparse_operator converts either a FermionOperator, a QubitOperator or other more advanced classes such as InteractionOperator to a scipy.sparse.csc matrix. There are numerous functions in fermilib.utils which one can call on the sparse operators such as \"get_gap\", \"get_hartree_fock_state\", \"get_ground_state\", ect. We show this off by computing the ground state energy of the Hubbard model. To do that, we use code from the fermilib.utils module which constructs lattice models of fermions such as Hubbard models." @@ -438,7 +518,11 @@ { "cell_type": "code", "execution_count": 10, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -566,7 +650,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "## Hamiltonians in the plane wave basis\n", "A user can write plugins to FermiLib which allow for the use of, e.g., third-party electronic structure package to compute molecular orbitals, Hamiltonians, energies, reduced density matrices, coupled cluster amplitudes, etc using Gaussian basis sets. We may provide scripts which interface between such packages and FermiLib in future but do not discuss them in this tutorial.\n", @@ -577,7 +664,11 @@ { "cell_type": "code", "execution_count": 11, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -611,23 +702,21 @@ } ], "source": [ - "from fermilib.utils import eigenspectrum, fourier_transform, jellium_model\n", + "from fermilib.utils import eigenspectrum, fourier_transform, jellium_model, Grid\n", "from fermilib.transforms import jordan_wigner\n", "\n", "# Let's look at a very small model of jellium in 1D.\n", - "n_dimensions = 1\n", - "grid_length = 3\n", - "length_scale = 1.\n", + "grid = Grid(dimensions=1, length=3, scale=1.0)\n", "spinless = True\n", "\n", "# Get the momentum Hamiltonian.\n", - "momentum_hamiltonian = jellium_model(n_dimensions, grid_length, length_scale, spinless)\n", + "momentum_hamiltonian = jellium_model(grid, spinless)\n", "momentum_qubit_operator = jordan_wigner(momentum_hamiltonian)\n", "momentum_qubit_operator.compress()\n", "print(momentum_qubit_operator)\n", "\n", "# Fourier transform the Hamiltonian to the position basis.\n", - "position_hamiltonian = fourier_transform(momentum_hamiltonian, n_dimensions, grid_length, length_scale, spinless)\n", + "position_hamiltonian = fourier_transform(momentum_hamiltonian, grid, spinless)\n", "position_qubit_operator = jordan_wigner(position_hamiltonian)\n", "position_qubit_operator.compress()\n", "print('')\n", @@ -641,7 +730,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "## Basics of MolecularData class\n", "\n", @@ -657,7 +749,11 @@ { "cell_type": "code", "execution_count": 12, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -700,7 +796,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "If we had previously computed this molecule using an electronic structure package, we can call molecule.load() to populate all sorts of interesting fields in the data structure. Though we make no assumptions about what electronic structure packages users might install, we assume that the calculations are saved in Fermilib's MolecularData objects. There may be plugins available in future. For the purposes of this example, we will load data that ships with FermiLib to make a plot of the energy surface of hydrogen. Note that helper functions to initialize some interesting chemical benchmarks are found in fermilib.utils." ] @@ -708,7 +807,11 @@ { "cell_type": "code", "execution_count": 13, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -985,7 +1088,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "## InteractionOperator and InteractionRDM for efficient numerical representations\n", "\n", @@ -993,7 +1099,7 @@ "\n", "Because the RDMs and molecular Hamiltonians are both compactly represented and manipulated as 2- and 4- index tensors, we can represent them in a particularly efficient form using similar data structures. The InteractionOperator data structure can be initialized for a Hamiltonian by passing the constant $h_0$ (or 0), as well as numpy arrays representing $h_{pq}$ (or $\\rho_{pq}$) and $h_{pqrs}$ (or $\\rho_{pqrs}$). Importantly, InteractionOperators can also be obtained by calling MolecularData.get_molecular_hamiltonian() or by calling the function get_interaction_operator() (found in fermilib.utils) on a FermionOperator. The InteractionRDM data structure is similar but represents RDMs. For instance, one can get a molecular RDM by calling MolecularData.get_molecular_rdm(). When generating Hamiltonians from the MolecularData class, one can choose to restrict the system to an active space.\n", "\n", - "These classes inherit from the same base class, InteractionTensor. This data structure overloads the slice operator [] so that one can get or set the key attributes of the InteractionOperator: .constant, .one_body_coefficients and .two_body_coefficients. For instance, InteractionOperator[p,q,r,s] would return $h_{pqrs}$ and InteractionRDM would return $\\rho_{pqrs}$. Importantly, the class supports fast basis transformations using the method InteractionTensor.rotate_basis(rotation_matrix).\n", + "These classes inherit from the same base class, InteractionTensor. This data structure overloads the slice operator [] so that one can get or set the key attributes of the InteractionOperator: $\\textrm{.constant}$, $\\textrm{.one_body_coefficients}$ and $\\textrm{.two_body_coefficients}$ . For instance, InteractionOperator[p,q,r,s] would return $h_{pqrs}$ and InteractionRDM would return $\\rho_{pqrs}$. Importantly, the class supports fast basis transformations using the method InteractionTensor.rotate_basis(rotation_matrix).\n", "But perhaps most importantly, one can map the InteractionOperator to any of the other data structures we've described here.\n", "\n", "Below, we load MolecularData from a saved calculation of LiH. We then obtain an InteractionOperator representation of this system in an active space. We then map that operator to qubits. We then demonstrate that one can rotate the orbital basis of the InteractionOperator using random angles to obtain a totally different operator that is still iso-spectral." @@ -1002,7 +1108,11 @@ { "cell_type": "code", "execution_count": 14, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -1138,7 +1248,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "## Simulating a variational quantum eigensolver using ProjectQ\n", "We now demonstrate how one can use both FermiLib and ProjectQ to run a simple VQE example using a Unitary Coupled Cluster ansatz. It demonstrates a simple way to evaluate the energy, optimize the energy with respect to the ansatz and build the corresponding compiled quantum circuit. It utilizes ProjectQ to build and simulate the circuit." @@ -1148,7 +1261,9 @@ "cell_type": "code", "execution_count": 6, "metadata": { - "collapsed": true + "collapsed": true, + "deletable": true, + "editable": true }, "outputs": [], "source": [ @@ -1166,15 +1281,22 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ - "Here we load H$_2$ from a precomputed molecule file found in the test data directory, and initialize the ProjectQ circuit compiler to a standard setting that uses a first-order Trotter decomposition to break up the exponentials of non-commuting operators. " + "Here we load $\\textrm{H}_2$ from a precomputed molecule file found in the test data directory, and initialize the ProjectQ circuit compiler to a standard setting that uses a first-order Trotter decomposition to break up the exponentials of non-commuting operators. " ] }, { "cell_type": "code", "execution_count": 7, - "metadata": {}, + "metadata": { + "collapsed": true, + "deletable": true, + "editable": true + }, "outputs": [], "source": [ "# Load the molecule.\n", @@ -1190,7 +1312,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "The Variational Quantum Eigensolver (or VQE), works by parameterizing a wavefunction $| \\Psi(\\theta) \\rangle$ through some quantum circuit, and minimzing the energy with respect to that angle, which is defined by\n", "\n", @@ -1205,7 +1330,9 @@ "cell_type": "code", "execution_count": 8, "metadata": { - "collapsed": true + "collapsed": true, + "deletable": true, + "editable": true }, "outputs": [], "source": [ @@ -1239,7 +1366,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "While we could plug this objective function into any optimizer, SciPy offers a convenient framework within the Python ecosystem. We'll choose as starting amplitudes the classical CCSD values that can be loaded from the molecule if desired. The optimal energy is found and compared to the exact values to verify that our simulation was successful." ] @@ -1247,7 +1377,11 @@ { "cell_type": "code", "execution_count": 9, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -1286,7 +1420,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "As we can see, the optimization terminates extremely quickly because the classical coupled cluster amplitudes were (for this molecule) already optimal. We can now use ProjectQ to compile this simulation circuit to a set of two-body quanutm gates." ] @@ -1294,7 +1431,11 @@ { "cell_type": "code", "execution_count": 10, - "metadata": {}, + "metadata": { + "collapsed": false, + "deletable": true, + "editable": true + }, "outputs": [ { "name": "stdout", @@ -1481,24 +1622,25 @@ } ], "metadata": { + "anaconda-cloud": {}, "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.12" + "pygments_lexer": "ipython3", + "version": "3.5.3" } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 0 } diff --git a/requirements.txt b/requirements.txt index c00aec8..c74ad79 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,4 +3,4 @@ numpy>=1.11.0 pytest>=3.0 h5py future -projectq>=0.3.0 +projectq>=0.3.1 diff --git a/src/fermilib/_version.py b/src/fermilib/_version.py index c647195..e03592d 100644 --- a/src/fermilib/_version.py +++ b/src/fermilib/_version.py @@ -11,4 +11,4 @@ # limitations under the License. """Define version number here and read it from setup.py automatically""" -__version__ = "0.1a0" +__version__ = "0.1a1" diff --git a/src/fermilib/data/H1-Li1_sto-3g_singlet.hdf5 b/src/fermilib/data/H1-Li1_sto-3g_singlet.hdf5 deleted file mode 100644 index 0f0ebc7..0000000 Binary files a/src/fermilib/data/H1-Li1_sto-3g_singlet.hdf5 and /dev/null differ diff --git a/src/fermilib/data/H1-Li1_sto-3g_singlet_1.45.hdf5 b/src/fermilib/data/H1-Li1_sto-3g_singlet_1.45.hdf5 new file mode 100644 index 0000000..003836c Binary files /dev/null and b/src/fermilib/data/H1-Li1_sto-3g_singlet_1.45.hdf5 differ diff --git a/src/fermilib/data/H1-Li1_sto-3g_singlet_cisd_rdm.npy b/src/fermilib/data/H1-Li1_sto-3g_singlet_cisd_rdm.npy deleted file mode 100644 index 88d92d3..0000000 Binary files a/src/fermilib/data/H1-Li1_sto-3g_singlet_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H1-Li1_sto-3g_singlet_eri.npy b/src/fermilib/data/H1-Li1_sto-3g_singlet_eri.npy deleted file mode 100644 index bafd6ee..0000000 Binary files a/src/fermilib/data/H1-Li1_sto-3g_singlet_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H1-Li1_sto-3g_singlet_fci_rdm.npy b/src/fermilib/data/H1-Li1_sto-3g_singlet_fci_rdm.npy deleted file mode 100644 index e8543a9..0000000 Binary files a/src/fermilib/data/H1-Li1_sto-3g_singlet_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet.hdf5 b/src/fermilib/data/H2_sto-3g_singlet.hdf5 deleted file mode 100644 index 75c3470..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet.hdf5 and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.1.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.1.hdf5 index d9b24ad..29c8d27 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.1.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.1.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.1_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.1_cisd_rdm.npy deleted file mode 100644 index 6bcfeb7..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.1_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.1_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_0.1_eri.npy deleted file mode 100644 index f99d960..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.1_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.1_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.1_fci_rdm.npy deleted file mode 100644 index 701a31d..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.1_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.2.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.2.hdf5 index 2be85f0..40d6cdf 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.2.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.2.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.2_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.2_cisd_rdm.npy deleted file mode 100644 index 200aab7..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.2_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.2_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_0.2_eri.npy deleted file mode 100644 index cbf9515..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.2_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.2_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.2_fci_rdm.npy deleted file mode 100644 index 1137407..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.2_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.3.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.3.hdf5 index b166533..4290a22 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.3.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.3.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.3_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.3_cisd_rdm.npy deleted file mode 100644 index 130c8fa..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.3_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.3_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_0.3_eri.npy deleted file mode 100644 index 891cdea..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.3_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.3_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.3_fci_rdm.npy deleted file mode 100644 index 3d5fb1a..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.3_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.4.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.4.hdf5 index 3ab05d5..a8e3af6 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.4.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.4.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.4_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.4_cisd_rdm.npy deleted file mode 100644 index cafc25a..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.4_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.4_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_0.4_eri.npy deleted file mode 100644 index be76c91..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.4_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.4_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.4_fci_rdm.npy deleted file mode 100644 index 0d6e8f9..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.4_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.5.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.5.hdf5 index 29d35f5..4bba283 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.5.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.5.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.5_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.5_cisd_rdm.npy deleted file mode 100644 index e23bf7f..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.5_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.5_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_0.5_eri.npy deleted file mode 100644 index c126d9c..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.5_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.5_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.5_fci_rdm.npy deleted file mode 100644 index 78dd3c4..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.5_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.6.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.6.hdf5 index 2886a70..8199091 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.6.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.6.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.6_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.6_cisd_rdm.npy deleted file mode 100644 index a8603e6..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.6_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.6_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_0.6_eri.npy deleted file mode 100644 index bec88be..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.6_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.6_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.6_fci_rdm.npy deleted file mode 100644 index e9b7493..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.6_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.7.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.7.hdf5 index c9948b1..0e68f61 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.7.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.7.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.7414.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.7414.hdf5 index 6967042..9db7414 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.7414.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.7414.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.7_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.7_cisd_rdm.npy deleted file mode 100644 index a4152df..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.7_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.7_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_0.7_eri.npy deleted file mode 100644 index 0c9714c..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.7_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.7_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.7_fci_rdm.npy deleted file mode 100644 index a2275f8..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.7_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.8.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.8.hdf5 index ae54c75..74dead8 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.8.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.8.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.8_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.8_cisd_rdm.npy deleted file mode 100644 index f0fdb4b..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.8_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.8_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_0.8_eri.npy deleted file mode 100644 index 541c6e6..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.8_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.8_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.8_fci_rdm.npy deleted file mode 100644 index 6ddc4f1..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.8_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.9.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_0.9.hdf5 index 81a7dbb..858e8ee 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.9.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_0.9.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.9_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.9_cisd_rdm.npy deleted file mode 100644 index 745338f..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.9_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.9_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_0.9_eri.npy deleted file mode 100644 index 358eb80..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.9_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_0.9_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_0.9_fci_rdm.npy deleted file mode 100644 index 417c1a5..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_0.9_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.0.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.0.hdf5 index da3ebfd..a22c99b 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.0.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.0.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.0_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.0_cisd_rdm.npy deleted file mode 100644 index 1301cb1..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.0_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.0_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.0_eri.npy deleted file mode 100644 index 4616a49..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.0_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.0_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.0_fci_rdm.npy deleted file mode 100644 index 71c09f9..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.0_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.1.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.1.hdf5 index a6a9939..8a40567 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.1.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.1.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.1_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.1_cisd_rdm.npy deleted file mode 100644 index fe9861b..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.1_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.1_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.1_eri.npy deleted file mode 100644 index 0dd231f..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.1_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.1_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.1_fci_rdm.npy deleted file mode 100644 index f783bc4..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.1_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.2.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.2.hdf5 index 5898b5c..a476580 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.2.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.2.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.2_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.2_cisd_rdm.npy deleted file mode 100644 index ff6ad20..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.2_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.2_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.2_eri.npy deleted file mode 100644 index 298f3af..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.2_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.2_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.2_fci_rdm.npy deleted file mode 100644 index 1a5d6d7..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.2_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.3.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.3.hdf5 index 86cc133..1902aa8 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.3.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.3.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.3_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.3_cisd_rdm.npy deleted file mode 100644 index e21c209..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.3_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.3_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.3_eri.npy deleted file mode 100644 index 81453e7..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.3_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.3_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.3_fci_rdm.npy deleted file mode 100644 index 491445f..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.3_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.4.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.4.hdf5 index 413c6e0..5a4837e 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.4.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.4.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.4_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.4_cisd_rdm.npy deleted file mode 100644 index 1cab1b6..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.4_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.4_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.4_eri.npy deleted file mode 100644 index 56de894..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.4_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.4_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.4_fci_rdm.npy deleted file mode 100644 index 277ccd3..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.4_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.5.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.5.hdf5 index 39b6d9d..b3ba604 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.5.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.5.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.5_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.5_cisd_rdm.npy deleted file mode 100644 index 1f2f015..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.5_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.5_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.5_eri.npy deleted file mode 100644 index c48bcf9..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.5_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.5_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.5_fci_rdm.npy deleted file mode 100644 index 688d201..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.5_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.6.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.6.hdf5 index 746ce65..b5f1250 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.6.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.6.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.6_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.6_cisd_rdm.npy deleted file mode 100644 index cb447ff..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.6_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.6_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.6_eri.npy deleted file mode 100644 index bb05a24..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.6_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.6_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.6_fci_rdm.npy deleted file mode 100644 index 6fbaf18..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.6_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.7.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.7.hdf5 index 48424d7..e3ddf26 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.7.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.7.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.7_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.7_cisd_rdm.npy deleted file mode 100644 index 2d47c2d..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.7_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.7_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.7_eri.npy deleted file mode 100644 index be69908..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.7_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.7_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.7_fci_rdm.npy deleted file mode 100644 index 8612d61..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.7_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.8.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.8.hdf5 index 44e5064..c94e748 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.8.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.8.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.8_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.8_cisd_rdm.npy deleted file mode 100644 index 2fb0105..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.8_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.8_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.8_eri.npy deleted file mode 100644 index 45b9c6a..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.8_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.8_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.8_fci_rdm.npy deleted file mode 100644 index a0b7496..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.8_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.9.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_1.9.hdf5 index 1786ad7..af88ba1 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.9.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_1.9.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.9_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.9_cisd_rdm.npy deleted file mode 100644 index e5bce86..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.9_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.9_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_1.9_eri.npy deleted file mode 100644 index 2268cc6..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.9_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_1.9_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_1.9_fci_rdm.npy deleted file mode 100644 index 3ee92ac..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_1.9_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.0.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.0.hdf5 index 1ee08b4..7c261d4 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.0.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.0.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.0_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.0_cisd_rdm.npy deleted file mode 100644 index c1db9a1..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.0_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.0_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.0_eri.npy deleted file mode 100644 index 34219ad..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.0_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.0_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.0_fci_rdm.npy deleted file mode 100644 index ea742e7..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.0_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.1.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.1.hdf5 index ed0fd17..ef6ebc3 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.1.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.1.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.1_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.1_cisd_rdm.npy deleted file mode 100644 index 3b8aa9f..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.1_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.1_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.1_eri.npy deleted file mode 100644 index d9f7a46..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.1_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.1_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.1_fci_rdm.npy deleted file mode 100644 index d9c04fd..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.1_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.2.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.2.hdf5 index 1519d28..90cbe61 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.2.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.2.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.2_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.2_cisd_rdm.npy deleted file mode 100644 index f37089c..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.2_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.2_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.2_eri.npy deleted file mode 100644 index 39ad057..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.2_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.2_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.2_fci_rdm.npy deleted file mode 100644 index fcc64e9..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.2_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.3.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.3.hdf5 index cb1a2e0..aefa672 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.3.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.3.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.3_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.3_cisd_rdm.npy deleted file mode 100644 index d905926..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.3_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.3_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.3_eri.npy deleted file mode 100644 index 4950da7..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.3_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.3_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.3_fci_rdm.npy deleted file mode 100644 index 3ba8a4a..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.3_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.4.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.4.hdf5 index 4acffcb..939c358 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.4.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.4.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.4_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.4_cisd_rdm.npy deleted file mode 100644 index cb3db25..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.4_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.4_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.4_eri.npy deleted file mode 100644 index 3a6fef8..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.4_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.4_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.4_fci_rdm.npy deleted file mode 100644 index c3ae695..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.4_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.5.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.5.hdf5 index dcb6feb..7b5e972 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.5.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.5.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.5_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.5_cisd_rdm.npy deleted file mode 100644 index d8cbd3d..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.5_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.5_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.5_eri.npy deleted file mode 100644 index 4368738..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.5_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.5_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.5_fci_rdm.npy deleted file mode 100644 index d8b9db5..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.5_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.6.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.6.hdf5 index c4ee440..99b209b 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.6.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.6.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.6_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.6_cisd_rdm.npy deleted file mode 100644 index 71087cb..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.6_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.6_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.6_eri.npy deleted file mode 100644 index b65a445..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.6_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.6_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.6_fci_rdm.npy deleted file mode 100644 index 0dd5f91..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.6_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.7.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.7.hdf5 index 6ed8101..7e08745 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.7.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.7.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.7_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.7_cisd_rdm.npy deleted file mode 100644 index a77f894..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.7_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.7_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.7_eri.npy deleted file mode 100644 index 5a9fc86..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.7_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.7_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.7_fci_rdm.npy deleted file mode 100644 index 7961b85..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.7_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.8.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.8.hdf5 index f706a45..085611b 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.8.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.8.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.8_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.8_cisd_rdm.npy deleted file mode 100644 index a2154dc..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.8_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.8_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.8_eri.npy deleted file mode 100644 index 395d589..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.8_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.8_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.8_fci_rdm.npy deleted file mode 100644 index 7d1fdb2..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.8_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.9.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_2.9.hdf5 index 33a0c9d..f9a71e9 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.9.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_2.9.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.9_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.9_cisd_rdm.npy deleted file mode 100644 index 1103b0b..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.9_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.9_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_2.9_eri.npy deleted file mode 100644 index a4c8ba1..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.9_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_2.9_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_2.9_fci_rdm.npy deleted file mode 100644 index 052d3dc..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_2.9_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_3.0.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_3.0.hdf5 index 196c7b2..1209525 100644 Binary files a/src/fermilib/data/H2_sto-3g_singlet_3.0.hdf5 and b/src/fermilib/data/H2_sto-3g_singlet_3.0.hdf5 differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_3.0_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_3.0_cisd_rdm.npy deleted file mode 100644 index ca74ddf..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_3.0_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_3.0_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_3.0_eri.npy deleted file mode 100644 index 2daa7c5..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_3.0_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_3.0_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_3.0_fci_rdm.npy deleted file mode 100644 index ca74ddf..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_3.0_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_cisd_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_cisd_rdm.npy deleted file mode 100644 index a476ec3..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_cisd_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_eq.hdf5 b/src/fermilib/data/H2_sto-3g_singlet_eq.hdf5 deleted file mode 100644 index 460fa2d..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_eq.hdf5 and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_eri.npy b/src/fermilib/data/H2_sto-3g_singlet_eri.npy deleted file mode 100644 index 6012aba..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_eri.npy and /dev/null differ diff --git a/src/fermilib/data/H2_sto-3g_singlet_fci_rdm.npy b/src/fermilib/data/H2_sto-3g_singlet_fci_rdm.npy deleted file mode 100644 index 40f3a7f..0000000 Binary files a/src/fermilib/data/H2_sto-3g_singlet_fci_rdm.npy and /dev/null differ diff --git a/src/fermilib/data/dummy_molecule.hdf5 b/src/fermilib/data/dummy_molecule.hdf5 deleted file mode 100644 index c3cc85b..0000000 Binary files a/src/fermilib/data/dummy_molecule.hdf5 and /dev/null differ diff --git a/src/fermilib/ops/_fermion_operator.py b/src/fermilib/ops/_fermion_operator.py index e828397..4e1dca7 100644 --- a/src/fermilib/ops/_fermion_operator.py +++ b/src/fermilib/ops/_fermion_operator.py @@ -39,6 +39,8 @@ def number_operator(n_orbitals, orbital=None, coefficient=1.): orbital (int, optional): The orbital on which to return the number operator. If None, return total number operator on all sites. coefficient (float): The coefficient of the term. + Returns: + operator (FermionOperator) """ if orbital is None: operator = FermionOperator() @@ -134,6 +136,30 @@ def normal_ordered(fermion_operator): return ordered_operator +def _parse_ladder_operator(ladder_operator_text): + """ + Args: + ladder_operator_text (str): + A ladder operator term like '4' or '5^', or an invalid string. + Returns: + tuple[int, int]: The mode then raise-vs-lower. + Raises: + FermionOperatorError: Given invalid text that doesn't match /\d+^?/ . + """ + inverted = 1 if ladder_operator_text.endswith('^') else 0 + mode_text = ladder_operator_text[:-1] if inverted else ladder_operator_text + + try: + mode = int(mode_text) + if mode < 0: + raise ValueError() # Merge with not-an-int failure case. + except ValueError: + raise FermionOperatorError( + "Invalid ladder operator term '{}'.".format(ladder_operator_text)) + + return mode, inverted + + class FermionOperator(object): """FermionOperator stores a sum of products of fermionic ladder operators. @@ -212,16 +238,8 @@ def __init__(self, term=None, coefficient=1.): # String input. if isinstance(term, str): - ladder_operators = [] - for ladder_operator in term.split(): - if ladder_operator[-1] == '^': - ladder_operators.append((int(ladder_operator[:-1]), 1)) - else: - try: - ladder_operators.append((int(ladder_operator), 0)) - except ValueError: - raise ValueError( - 'Invalid action provided to FermionTerm.') + ladder_operators = tuple(_parse_ladder_operator(e) + for e in term.split()) self.terms[tuple(ladder_operators)] = coefficient # Tuple input. @@ -245,13 +263,33 @@ def __init__(self, term=None, coefficient=1.): 'Invalid action in FermionOperator: ' 'Must be 0 (lowering) or 1 (raising).') + @staticmethod + def zero(): + """ + Returns: + additive_identity (FermionOperator): + A fermion operator o with the property that o+x = x+o = x for + all fermion operators x. + """ + return FermionOperator(term=None) + + @staticmethod + def identity(): + """ + Returns: + multiplicative_identity (FermionOperator): + A fermion operator u with the property that u*x = x*u = x for + all fermion operators x. + """ + return FermionOperator(term=()) + def compress(self, abs_tol=EQ_TOLERANCE): """ Eliminates all terms with coefficients close to zero and removes imaginary parts of coefficients that are close to zero. Args: - abs_tol(float): Absolute tolerance, must be at least 0.0 + abs_tol (float): Absolute tolerance, must be at least 0.0 """ new_terms = {} for term in self.terms: @@ -357,6 +395,8 @@ def __imul__(self, multiplier): Args: multiplier(complex float, or FermionOperator): multiplier + Returns: + product (FermionOperator): Mutated self. """ # Handle scalars. if isinstance(multiplier, (int, float, complex)): @@ -388,7 +428,7 @@ def __mul__(self, multiplier): multiplier: A scalar, or a FermionOperator. Returns: - product: A FermionOperator. + product (FermionOperator) Raises: TypeError: Invalid type cannot be multiply with FermionOperator. @@ -405,14 +445,14 @@ def __rmul__(self, multiplier): """Return multiplier * self for a scalar. We only define __rmul__ for scalars because the left multiply - exist for FermionOperator and left multiply + exist for FermionOperator and left multiply is also queried as the default behavior. Args: multiplier: A scalar to multiply by. Returns: - product: A new instance of FermionOperator. + product (FermionOperator) Raises: TypeError: Object of invalid type cannot multiply FermionOperator. @@ -428,10 +468,10 @@ def __truediv__(self, divisor): Note that this is always floating point division. Args: - divisor: A scalar to divide by. + divisor (int|float|complex): A scalar to divide by. Returns: - A new instance of FermionOperator. + quotient (FermionOperator) Raises: TypeError: Cannot divide local operator by non-scalar type. @@ -441,17 +481,33 @@ def __truediv__(self, divisor): return self * (1.0 / divisor) def __div__(self, divisor): - """For compatibility with Python 2. """ + """For compatibility with Python 2. + Args: + divisor (int|float|complex): A scalar to divide by. + Returns: + quotient (FermionOperator) + """ return self.__truediv__(divisor) def __itruediv__(self, divisor): + """ + Args: + divisor (int|float|complex): A scalar to divide by. + Returns: + quotient (FermionOperator): Mutated self. + """ if not isinstance(divisor, (int, float, complex)): raise TypeError('Cannot divide QubitOperator by non-scalar type.') self *= (1.0 / divisor) return self def __idiv__(self, divisor): - """For compatibility with Python 2. """ + """For compatibility with Python 2. + Args: + divisor (int|float|complex): A scalar to divide by. + Returns: + quotient (FermionOperator): Mutated self. + """ return self.__itruediv__(divisor) def __iadd__(self, addend): @@ -460,6 +516,9 @@ def __iadd__(self, addend): Args: addend: A FermionOperator. + Returns: + sum (FermionOperator): Mutated self. + Raises: TypeError: Cannot add invalid type. """ @@ -478,28 +537,43 @@ def __iadd__(self, addend): return self def __add__(self, addend): - """Return self + addend for a FermionOperator. """ + """ + Args: + addend (FermionOperator): The operator to add. + + Returns: + sum (FermionOperator) + """ summand = copy.deepcopy(self) summand += addend return summand def __sub__(self, subtrahend): - """Return self - subtrahend for a FermionOperator.""" + """ + Args: + subtrahend (FermionOperator): The operator to subtract. + Returns: + difference (FermionOperator) + """ if not isinstance(subtrahend, FermionOperator): raise TypeError('Cannot subtract invalid type to FermionOperator.') return self + (-1. * subtrahend) def __neg__(self): + """ + Returns: + negation (FermionOperator) + """ return -1 * self def __pow__(self, exponent): """Exponentiate the FermionOperator. Args: - exponent: An int, the exponent with which to raise the operator. + exponent (int): The exponent with which to raise the operator. Returns: - exponentiated: The exponentiated operator. + exponentiated (FermionOperator) Raises: ValueError: Can only raise FermionOperator to non-negative @@ -508,7 +582,8 @@ def __pow__(self, exponent): # Handle invalid exponents. if not isinstance(exponent, int) or exponent < 0: raise ValueError( - 'Can only raise FermionOperator to positive integer powers.') + 'exponent must be a non-negative int, but was {} {}'.format( + type(exponent), repr(exponent))) # Initialized identity. exponentiated = FermionOperator(()) diff --git a/src/fermilib/ops/_fermion_operator_test.py b/src/fermilib/ops/_fermion_operator_test.py index 4216039..0276185 100644 --- a/src/fermilib/ops/_fermion_operator_test.py +++ b/src/fermilib/ops/_fermion_operator_test.py @@ -56,27 +56,79 @@ def test_init_tuple_npcomplex128_coefficient(self): self.assertEqual(len(fermion_op.terms), 1) self.assertEqual(fermion_op.terms[loc_op], coefficient) + def test_identity_is_multiplicative_identity(self): + u = FermionOperator.identity() + f = FermionOperator(((0, 1), (5, 0), (6, 1)), 0.6j) + g = FermionOperator(((0, 0), (5, 0), (6, 1)), 0.3j) + h = f + g + self.assertTrue(f.isclose(u * f)) + self.assertTrue(f.isclose(f * u)) + self.assertTrue(g.isclose(u * g)) + self.assertTrue(g.isclose(g * u)) + self.assertTrue(h.isclose(u * h)) + self.assertTrue(h.isclose(h * u)) + + u *= h + self.assertTrue(h.isclose(u)) + self.assertFalse(f.isclose(u)) + + # Method always returns new instances. + self.assertFalse(FermionOperator.identity().isclose(u)) + + def test_zero_is_additive_identity(self): + o = FermionOperator.zero() + f = FermionOperator(((0, 1), (5, 0), (6, 1)), 0.6j) + g = FermionOperator(((0, 0), (5, 0), (6, 1)), 0.3j) + h = f + g + self.assertTrue(f.isclose(o + f)) + self.assertTrue(f.isclose(f + o)) + self.assertTrue(g.isclose(o + g)) + self.assertTrue(g.isclose(g + o)) + self.assertTrue(h.isclose(o + h)) + self.assertTrue(h.isclose(h + o)) + + o += h + self.assertTrue(h.isclose(o)) + self.assertFalse(f.isclose(o)) + + # Method always returns new instances. + self.assertFalse(FermionOperator.zero().isclose(o)) + + def test_zero_is_multiplicative_nil(self): + o = FermionOperator.zero() + u = FermionOperator.identity() + f = FermionOperator(((0, 1), (5, 0), (6, 1)), 0.6j) + g = FermionOperator(((0, 0), (5, 0), (6, 1)), 0.3j) + self.assertTrue(o.isclose(o * u)) + self.assertTrue(o.isclose(o * f)) + self.assertTrue(o.isclose(o * g)) + self.assertTrue(o.isclose(o * (f + g))) + def test_init_str(self): fermion_op = FermionOperator('0^ 5 12^', -1.) correct = ((0, 1), (5, 0), (12, 1)) self.assertIn(correct, fermion_op.terms) self.assertEqual(fermion_op.terms[correct], -1.0) + def test_merges_multiple_whitespace(self): + fermion_op = FermionOperator(' \n ') + self.assertEqual(fermion_op.terms, {(): 1}) + def test_init_str_identity(self): fermion_op = FermionOperator('') self.assertIn((), fermion_op.terms) def test_init_bad_term(self): with self.assertRaises(ValueError): - fermion_op = FermionOperator(list()) + _ = FermionOperator(list()) def test_init_bad_coefficient(self): with self.assertRaises(ValueError): - fermion_op = FermionOperator('0^', "0.5") + _ = FermionOperator('0^', "0.5") def test_init_bad_action_str(self): - with self.assertRaises(ValueError): - fermion_op = FermionOperator('0-') + with self.assertRaises(FermionOperatorError): + _ = FermionOperator('0-') def test_init_bad_action_tuple(self): with self.assertRaises(ValueError): @@ -84,15 +136,15 @@ def test_init_bad_action_tuple(self): def test_init_bad_tuple(self): with self.assertRaises(ValueError): - fermion_op = FermionOperator(((0, 1, 1),)) + _ = FermionOperator(((0, 1, 1),)) def test_init_bad_str(self): - with self.assertRaises(ValueError): - fermion_op = FermionOperator('^') + with self.assertRaises(FermionOperatorError): + _ = FermionOperator('^') def test_init_bad_mode_num(self): with self.assertRaises(FermionOperatorError): - fermion_op = FermionOperator('-1^') + _ = FermionOperator('-1^') def test_FermionOperator(self): op = FermionOperator((), 3.) @@ -423,7 +475,7 @@ def test_add(self): def test_add_bad_addend(self): op = FermionOperator((), 1.0) with self.assertRaises(TypeError): - op = op + "0.5" + _ = op + "0.5" def test_sub(self): term_a = ((1, 1), (3, 1), (8, 1)) @@ -442,7 +494,7 @@ def test_sub(self): def test_sub_bad_subtrahend(self): op = FermionOperator((), 1.0) with self.assertRaises(TypeError): - op = op - "0.5" + _ = op - "0.5" def test_isub_different_term(self): term_a = ((1, 1), (3, 1), (8, 0)) @@ -464,7 +516,7 @@ def test_isub_bad_addend(self): def test_neg(self): op = FermionOperator(((1, 1), (3, 1), (8, 1)), 0.5) - -op + _ = -op # out of place self.assertTrue(op.isclose(FermionOperator(((1, 1), (3, 1), (8, 1)), 0.5))) diff --git a/src/fermilib/ops/_interaction_operator.py b/src/fermilib/ops/_interaction_operator.py index c8ccd13..9fe5a4b 100644 --- a/src/fermilib/ops/_interaction_operator.py +++ b/src/fermilib/ops/_interaction_operator.py @@ -75,8 +75,10 @@ def unique_iter(self, complex_valued=False): 2. pqrs = rqps = psrq = srqp = qpsr = rspq = spqr = qrsp. Args: - complex_valued: Bool, whether the operator has complex - coefficients. + complex_valued (bool): + Whether the operator has complex coefficients. + Yields: + tuple[int] """ # Constant. if self.constant: @@ -86,23 +88,24 @@ def unique_iter(self, complex_valued=False): for p in range(self.n_qubits): for q in range(p + 1): if self.one_body_tensor[p, q]: - yield [p, q] + yield p, q # Two-body terms. - record_map = {} - for p in range(self.n_qubits): - for q in range(self.n_qubits): - for r in range(self.n_qubits): - for s in range(self.n_qubits): - if self.two_body_tensor[p, q, r, s] and \ - (p, q, r, s) not in record_map: - yield [p, q, r, s] - record_map[(p, q, r, s)] = [] - record_map[(s, r, q, p)] = [] - record_map[(q, p, s, r)] = [] - record_map[(r, s, p, q)] = [] - if not complex_valued: - record_map[(p, s, r, q)] = [] - record_map[(s, p, q, r)] = [] - record_map[(q, r, s, p)] = [] - record_map[(r, q, p, s)] = [] + seen = set() + for quad in itertools.product(range(self.n_qubits), repeat=4): + if self[quad] and quad not in seen: + seen |= set(_symmetric_two_body_terms(quad, complex_valued)) + yield quad + + +def _symmetric_two_body_terms(quad, complex_valued): + p, q, r, s = quad + yield p, q, r, s + yield q, p, s, r + yield s, r, q, p + yield r, s, p, q + if not complex_valued: + yield p, s, r, q + yield q, r, s, p + yield s, p, q, r + yield r, q, p, s diff --git a/src/fermilib/ops/_interaction_rdm.py b/src/fermilib/ops/_interaction_rdm.py index 818588e..3d4e3c3 100644 --- a/src/fermilib/ops/_interaction_rdm.py +++ b/src/fermilib/ops/_interaction_rdm.py @@ -35,7 +35,6 @@ class InteractionRDM(InteractionTensor): one_body_tensor: The expectation values . two_body_tensor: The expectation values . - n_qubits: An int giving the number of qubits. """ def __init__(self, one_body_tensor, two_body_tensor): @@ -70,9 +69,11 @@ def expectation(self, operator): InteractionRDMError: Invalid operator provided. """ if isinstance(operator, QubitOperator): - expectation_value = 0. - for qubit_term in operator: - expectation += qubit_term_expectation(self, qubit_term) + expectation_op = self.get_qubit_expectations(operator) + expectation = 0.0 + for qubit_term in operator.terms: + expectation += (operator.terms[qubit_term] * + expectation_op.terms[qubit_term]) elif isinstance(operator, InteractionOperator): expectation = operator.constant expectation += numpy.sum(self.one_body_tensor * @@ -99,13 +100,12 @@ def get_qubit_expectations(self, qubit_operator): """ from fermilib.transforms import reverse_jordan_wigner qubit_operator_expectations = copy.deepcopy(qubit_operator) - del qubit_operator_expectations.terms[()] for qubit_term in qubit_operator_expectations.terms: expectation = 0. # Map qubits back to fermions. reversed_fermion_operators = reverse_jordan_wigner( - QubitOperator(qubit_term), self.n_qubits) + QubitOperator(qubit_term)) reversed_fermion_operators = normal_ordered( reversed_fermion_operators) diff --git a/src/fermilib/ops/_interaction_rdm_test.py b/src/fermilib/ops/_interaction_rdm_test.py index c9a0d5d..0fe3c6d 100644 --- a/src/fermilib/ops/_interaction_rdm_test.py +++ b/src/fermilib/ops/_interaction_rdm_test.py @@ -26,7 +26,8 @@ def setUp(self): geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))] basis = 'sto-3g' multiplicity = 1 - filename = os.path.join(THIS_DIRECTORY, 'data', 'H2_sto-3g_singlet') + filename = os.path.join(THIS_DIRECTORY, 'data', + 'H2_sto-3g_singlet_0.7414') self.molecule = MolecularData( geometry, basis, multiplicity, filename=filename) self.molecule.load() @@ -38,7 +39,7 @@ def test_get_qubit_expectations(self): qubit_operator = jordan_wigner(self.hamiltonian) qubit_expectations = self.rdm.get_qubit_expectations(qubit_operator) - test_energy = qubit_operator.terms[()] + test_energy = 0.0 for qubit_term in qubit_expectations.terms: term_coefficient = qubit_operator.terms[qubit_term] test_energy += (term_coefficient * diff --git a/src/fermilib/tests/_hydrogen_integration_test.py b/src/fermilib/tests/_hydrogen_integration_test.py index 3db1ced..5397074 100644 --- a/src/fermilib/tests/_hydrogen_integration_test.py +++ b/src/fermilib/tests/_hydrogen_integration_test.py @@ -31,7 +31,8 @@ def setUp(self): geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))] basis = 'sto-3g' multiplicity = 1 - filename = os.path.join(THIS_DIRECTORY, 'data', 'H2_sto-3g_singlet') + filename = os.path.join(THIS_DIRECTORY, 'data', + 'H2_sto-3g_singlet_0.7414') self.molecule = MolecularData( geometry, basis, multiplicity, filename=filename) self.molecule.load() @@ -192,7 +193,7 @@ def test_rdm_numerically(self): # Confirm expectation on qubit Hamiltonian using reverse JW matches. qubit_rdm = self.fci_rdm.get_qubit_expectations(self.qubit_hamiltonian) - qubit_energy = self.qubit_hamiltonian.terms[()] + qubit_energy = 0.0 for term, expectation in qubit_rdm.terms.items(): qubit_energy += expectation * self.qubit_hamiltonian.terms[term] self.assertAlmostEqual(qubit_energy, self.molecule.fci_energy) @@ -272,9 +273,6 @@ def test_ucc(self): self.hamiltonian_matrix.dot(ccsd_state_r))[0, 0] self.assertAlmostEqual(expected_ccsd_energy, self.molecule.fci_energy) - def test_version(self): - self.assertEqual(__version__, '0.1a0') - if __name__ == '__main__': unittest.main() diff --git a/src/fermilib/tests/_lih_integration_test_.py b/src/fermilib/tests/_lih_integration_test.py similarity index 89% rename from src/fermilib/tests/_lih_integration_test_.py rename to src/fermilib/tests/_lih_integration_test.py index c8c171a..1a29c6e 100644 --- a/src/fermilib/tests/_lih_integration_test_.py +++ b/src/fermilib/tests/_lih_integration_test.py @@ -27,20 +27,24 @@ class LiHIntegrationTest(unittest.TestCase): def setUp(self): - # Set up molecule. geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.45))] basis = 'sto-3g' multiplicity = 1 - filename = os.path.join(THIS_DIRECTORY, 'data', 'H1-Li1_sto-3g_singlet') + filename = os.path.join(THIS_DIRECTORY, 'data', + 'H1-Li1_sto-3g_singlet_1.45') self.molecule = MolecularData( geometry, basis, multiplicity, filename=filename) self.molecule.load() # Get molecular Hamiltonian. self.molecular_hamiltonian = self.molecule.get_molecular_hamiltonian() - self.molecular_hamiltonian_no_core = self.molecule.\ - get_molecular_hamiltonian(active_space_start=1) + self.molecular_hamiltonian_no_core = ( + self.molecule. + get_molecular_hamiltonian(occupied_indices=[0], + active_indices=range(1, + self.molecule. + n_orbitals))) # Get FCI RDM. self.fci_rdm = self.molecule.get_molecular_rdm(use_fci=1) @@ -90,7 +94,7 @@ def test_all(self): # Confirm expectation on qubit Hamiltonian using reverse JW matches. qubit_rdm = self.fci_rdm.get_qubit_expectations(self.qubit_hamiltonian) - qubit_energy = self.qubit_hamiltonian.terms[()] + qubit_energy = 0.0 for term, coefficient in qubit_rdm.terms.items(): qubit_energy += coefficient * self.qubit_hamiltonian.terms[term] self.assertAlmostEqual(qubit_energy, self.molecule.fci_energy) @@ -121,8 +125,8 @@ def test_all(self): # Check that frozen core result matches frozen core FCI from psi4. # Recore frozen core result from external calculation. self.frozen_core_fci_energy = -7.8807607374168 - no_core_fci_energy = scipy. \ - linalg.eigh(self.hamiltonian_matrix_no_core.todense())[0][0] + no_core_fci_energy = scipy.linalg.eigh( + self.hamiltonian_matrix_no_core.todense())[0][0] self.assertAlmostEqual(no_core_fci_energy, self.frozen_core_fci_energy) diff --git a/src/fermilib/transforms/_bravyi_kitaev.py b/src/fermilib/transforms/_bravyi_kitaev.py index 0a6d02b..fb787c1 100644 --- a/src/fermilib/transforms/_bravyi_kitaev.py +++ b/src/fermilib/transforms/_bravyi_kitaev.py @@ -22,7 +22,11 @@ def bravyi_kitaev(operator, n_qubits=None): """Apply the Bravyi-Kitaev transform and return qubit operator. Args: - operator: A FermionOperator to transform. + operator (fermilib.ops.FermionOperator): + A FermionOperator to transform. + n_qubits (int|None): + Can force the number of qubits in the resulting operator above the + number that appear in the input operator. Returns: transformed_operator: An instance of the QubitOperator class. @@ -37,53 +41,102 @@ def bravyi_kitaev(operator, n_qubits=None): if n_qubits < count_qubits(operator): raise ValueError('Invalid number of qubits specified.') + # Build the Fenwick tree. + fenwick_tree = FenwickTree(n_qubits) + # Compute transformed operator. - transformed_operator = QubitOperator() - for term in operator.terms: - - # Initialize identity matrix. - coefficient = operator.terms[term] - transformed_term = QubitOperator((), coefficient) - - # Build the Fenwick Tree - fenwick_tree = FenwickTree(n_qubits) - - # Build the Bravyi-Kitaev transformed operators. - for ladder_operator in term: - index = ladder_operator[0] - - # Parity set. Set of nodes to apply Z to. - parity_set = [node.index for node in - fenwick_tree.get_parity_set(index)] - - # Update set. Set of ancestors to apply X to. - ancestors = [node.index for node in - fenwick_tree.get_update_set(index)] - - # The C(j) set. - ancestor_children = [node.index for node in - fenwick_tree.get_remainder_set(index)] - - # Switch between lowering/raising operators. - d_coefficient = .5j - if ladder_operator[1]: - d_coefficient *= -1. - - # The fermion lowering operator is given by - # a = (c+id)/2 where c, d are the majoranas. - d_majorana_component = QubitOperator( - (((ladder_operator[0], 'Y'),) + - tuple((index, 'Z') for index in ancestor_children) + - tuple((index, 'X') for index in ancestors)), - d_coefficient) - - c_majorana_component = QubitOperator( - (((ladder_operator[0], 'X'),) + - tuple((index, 'Z') for index in parity_set) + - tuple((index, 'X') for index in ancestors)), - 0.5) - - # Update term. - transformed_term *= c_majorana_component + d_majorana_component - transformed_operator += transformed_term - return transformed_operator + transformed_terms = ( + _transform_operator_term(term=term, + coefficient=operator.terms[term], + fenwick_tree=fenwick_tree) + for term in operator.terms + ) + return inline_sum(seed=QubitOperator(), summands=transformed_terms) + + +def _transform_operator_term(term, coefficient, fenwick_tree): + """ + Args: + term (list[tuple[int, int]]): + A list of (mode, raising-vs-lowering) ladder operator terms. + coefficient (float): + fenwick_tree (FenwickTree): + Returns: + QubitOperator: + """ + + # Build the Bravyi-Kitaev transformed operators. + transformed_ladder_ops = ( + _transform_ladder_operator(ladder_operator, fenwick_tree) + for ladder_operator in term + ) + return inline_product(seed=QubitOperator((), coefficient), + factors=transformed_ladder_ops) + + +def _transform_ladder_operator(ladder_operator, fenwick_tree): + """ + Args: + ladder_operator (tuple[int, int]): + fenwick_tree (FenwickTree): + Returns: + QubitOperator: + """ + index = ladder_operator[0] + + # Parity set. Set of nodes to apply Z to. + parity_set = [node.index for node in + fenwick_tree.get_parity_set(index)] + + # Update set. Set of ancestors to apply X to. + ancestors = [node.index for node in + fenwick_tree.get_update_set(index)] + + # The C(j) set. + ancestor_children = [node.index for node in + fenwick_tree.get_remainder_set(index)] + + # Switch between lowering/raising operators. + d_coefficient = -.5j if ladder_operator[1] else .5j + + # The fermion lowering operator is given by + # a = (c+id)/2 where c, d are the majoranas. + d_majorana_component = QubitOperator( + (((ladder_operator[0], 'Y'),) + + tuple((index, 'Z') for index in ancestor_children) + + tuple((index, 'X') for index in ancestors)), + d_coefficient) + + c_majorana_component = QubitOperator( + (((ladder_operator[0], 'X'),) + + tuple((index, 'Z') for index in parity_set) + + tuple((index, 'X') for index in ancestors)), + 0.5) + + return c_majorana_component + d_majorana_component + + +def inline_sum(seed, summands): + """Computes a sum, using the __iadd__ operator. + Args: + seed (T): The starting total. The zero value. + summands (iterable[T]): Values to add (with +=) into the total. + Returns: + T: The result of adding all the factors into the zero value. + """ + for r in summands: + seed += r + return seed + + +def inline_product(seed, factors): + """Computes a product, using the __imul__ operator. + Args: + seed (T): The starting total. The unit value. + factors (iterable[T]): Values to multiply (with *=) into the total. + Returns: + T: The result of multiplying all the factors into the unit value. + """ + for r in factors: + seed *= r + return seed diff --git a/src/fermilib/transforms/_bravyi_kitaev_test.py b/src/fermilib/transforms/_bravyi_kitaev_test.py index 3cf2470..99fd2b7 100644 --- a/src/fermilib/transforms/_bravyi_kitaev_test.py +++ b/src/fermilib/transforms/_bravyi_kitaev_test.py @@ -190,8 +190,7 @@ def test_bk_jw_integration_original(self): jw_spectrum = eigenspectrum(jw_qubit_operator) bk_spectrum = eigenspectrum(bk_qubit_operator) self.assertAlmostEqual(0., numpy.amax(numpy.absolute(jw_spectrum - - bk_spectrum)), - places=5) + bk_spectrum)), places=5) if __name__ == '__main__': diff --git a/src/fermilib/transforms/_conversion.py b/src/fermilib/transforms/_conversion.py index 664044f..8b33aa6 100644 --- a/src/fermilib/transforms/_conversion.py +++ b/src/fermilib/transforms/_conversion.py @@ -12,15 +12,15 @@ """Transformations acting on operators and RDMs.""" from __future__ import absolute_import -from future.utils import iteritems -import copy import itertools + import numpy +from future.utils import iteritems +from projectq.ops import QubitOperator from fermilib.ops import (FermionOperator, normal_ordered, - number_operator, InteractionOperator, InteractionRDM) from fermilib.ops._interaction_operator import InteractionOperatorError @@ -28,15 +28,13 @@ jordan_wigner_sparse, qubit_operator_sparse) -from projectq.ops import QubitOperator - def get_sparse_operator(operator, n_qubits=None): """Map a Fermion, Qubit, or InteractionOperator to a SparseOperator.""" if isinstance(operator, InteractionOperator): return get_sparse_interaction_operator(operator) elif isinstance(operator, FermionOperator): - return jordan_wigner_sparse(operator) + return jordan_wigner_sparse(operator, n_qubits) elif isinstance(operator, QubitOperator): if n_qubits is None: n_qubits = count_qubits(operator) @@ -161,17 +159,15 @@ def get_fermion_operator(interaction_operator): n_qubits = count_qubits(interaction_operator) # Add one-body terms. - for p in range(n_qubits): - for q in range(n_qubits): - coefficient = interaction_operator[p, q] - fermion_operator += FermionOperator(((p, 1), (q, 0)), coefficient) - - # Add two-body terms. - for r in range(n_qubits): - for s in range(n_qubits): - coefficient = interaction_operator[p, q, r, s] - fermion_operator += FermionOperator(((p, 1), (q, 1), - (r, 0), (s, 0)), - coefficient) + for p, q in itertools.product(range(n_qubits), repeat=2): + fermion_operator += FermionOperator( + ((p, 1), (q, 0)), coefficient=interaction_operator[p, q]) + + # Add two-body terms. + for p, q, r, s in itertools.product(range(n_qubits), repeat=4): + coefficient = interaction_operator[p, q, r, s] + fermion_operator += FermionOperator(((p, 1), (q, 1), + (r, 0), (s, 0)), + coefficient) return fermion_operator diff --git a/src/fermilib/transforms/_conversion_test.py b/src/fermilib/transforms/_conversion_test.py index 89e0785..fc09780 100644 --- a/src/fermilib/transforms/_conversion_test.py +++ b/src/fermilib/transforms/_conversion_test.py @@ -115,5 +115,14 @@ def test_sparse_matrix_linearity(self): [0, 3, 5, 6, 9, 10, 12, 15]) +class GetSparseOperatorFermionTest(unittest.TestCase): + + def test_sparse_matrix_zero_n_qubit(self): + sparse_operator = get_sparse_operator(FermionOperator.zero(), 4) + sparse_operator.eliminate_zeros() + self.assertEqual(len(list(sparse_operator.data)), 0) + self.assertEqual(sparse_operator.shape, (16, 16)) + + if __name__ == '__main__': unittest.main() diff --git a/src/fermilib/transforms/_jordan_wigner.py b/src/fermilib/transforms/_jordan_wigner.py index 8683c9d..e8edb69 100644 --- a/src/fermilib/transforms/_jordan_wigner.py +++ b/src/fermilib/transforms/_jordan_wigner.py @@ -151,8 +151,8 @@ def jordan_wigner_two_body(p, q, r, s): elif len(set([p, q, r, s])) == 4: # Loop through different operators which act on each tensor factor. - for operator_p, operator_q, operator_r in \ - itertools.product(['X', 'Y'], repeat=3): + for operator_p, operator_q, operator_r in itertools.product( + ['X', 'Y'], repeat=3): if [operator_p, operator_q, operator_r].count('X') % 2: operator_s = 'X' else: diff --git a/src/fermilib/utils/__init__.py b/src/fermilib/utils/__init__.py index 4000d28..e9ebcc0 100644 --- a/src/fermilib/utils/__init__.py +++ b/src/fermilib/utils/__init__.py @@ -13,6 +13,9 @@ from ._chemical_series import (make_atomic_ring, make_atomic_lattice, make_atom) +from ._graph import Graph, Node + +from ._grid import Grid from ._hubbard import fermi_hubbard @@ -30,7 +33,8 @@ from ._plane_wave_hamiltonian import (inverse_fourier_transform, fourier_transform, - plane_wave_hamiltonian) + plane_wave_hamiltonian, + jordan_wigner_dual_basis_hamiltonian) from ._sparse_tools import (expectation, get_density_matrix, diff --git a/src/fermilib/utils/_graph.py b/src/fermilib/utils/_graph.py new file mode 100644 index 0000000..2a83456 --- /dev/null +++ b/src/fermilib/utils/_graph.py @@ -0,0 +1,228 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from __future__ import absolute_import + +import queue + + +class Node: + """Graph node + + These graph nodes may be used to store data or other attributes within + the Graph structure. + """ + def __init__(self, value=None): + """Build a graph node initialized with generic value""" + self.value = value + + +class Graph: + """An undirected graph of nodes of arbitrary connectivity + + A generic graph class for undirected graphs that holds Nodes and + edges that connect them. + + Attributes: + nodes(list): A list of Node objects containing the nodes of the graph + node_uids(list): A list of unique IDs (uids) for each node + uid_to_index(dict): A dictionary that maps UIDs currently present + in the graph to their node index + neigbors(list of sets): A list of sets that enumerate the neighbors + of each node. For example the neighbors of node i are the set + neighbors[i] + next_uid(int): The next unique ID to be assigned to a node on addition + """ + + def __init__(self): + """Set up an empty graph with no nodes""" + self.nodes = [] + self.node_uids = [] + self.uid_to_index = {} + self.neighbors = [] + self.next_uid = 0 + + def add_node(self, node=Node()): + """Add a Node to the graph. + + Args: + node(Node): A Node object to be added to the graph + + Returns(int): The unique identified that was given to the node + """ + self.nodes += [node] + self.node_uids += [self.next_uid] + self.neighbors += [set()] + self.uid_to_index[self.next_uid] = self.node_count() - 1 + self.next_uid += 1 + + def remove_node(self, node_id): + """Remove a graph node + + This removes a node from the graph and leverages the unique ID + system used internally to avoid having to modify the edges for all + nodes in the graph. + + Args: + node_id(int): Index of the node to be removed. + """ + if node_id >= self.node_count(): + raise IndexError("Node ID out of range for graph") + node_uid = self.node_uids[node_id] + for i in range(self.node_count()): + if node_uid in self.neighbors[i]: + self.neighbors[i].remove(node_uid) + for i in range(node_id + 1, self.node_count()): + self.uid_to_index[self.node_uids[i]] -= 1 + + del self.uid_to_index[node_uid] + del self.nodes[node_id] + del self.neighbors[node_id] + del self.node_uids[node_id] + + def node_count(self): + """Number of nodes in the graph + + Returns(int): Number of nodes currently in the graph + """ + return len(self.nodes) + + def add_edge(self, node_id1, node_id2): + """Add an edge between node 1 and node2 + + Args: + node_id1(int): Index of first node on edge + node_id2(int): Index of second node on edge + """ + if node_id1 >= self.node_count() or node_id2 >= self.node_count(): + raise IndexError("Node ID out of range for graph.") + if node_id1 == node_id2: + raise IndexError("Error, self-loops not supported.") + + self.neighbors[node_id1].add(self.node_uids[node_id2]) + self.neighbors[node_id2].add(self.node_uids[node_id1]) + + def remove_edge(self, node_id1, node_id2): + """Remove an edge between node1 and node2 + + Args: + node_id1(int): Index of first node on edge + node_id2(int): Index of second node on edge + """ + if node_id1 >= self.node_count() or node_id2 >= self.node_count(): + raise IndexError("Node ID out of range for graph") + + self.neighbors[node_id1].remove(self.node_uids[node_id2]) + self.neighbors[node_id2].remove(self.node_uids[node_id1]) + + def find_index(self, value, starting_node=0): + """Find the index of the first node that matches value in a BFS + + Performs a breadth-first search of the graph starting at node index + starting_node. Returns the index or None if no match is found + + Args: + value(Node Value) - Value to match against in the graph + starting_node(int) - Node index to start search from + """ + if starting_node > self.node_count(): + raise IndexError("Node ID out of range.") + + node_queue = queue.LifoQueue() + node_queue.put(starting_node) + visited = [starting_node] + while not node_queue.empty(): + next_id = node_queue.get() + if (self.nodes[next_id].value == value): + return next_id # Success + for uid in self.neighbors[next_id]: + if (self.uid_to_index[uid] not in visited): + node_queue.put(self.uid_to_index[uid]) + visited += [self.uid_to_index[uid]] + return None + + def is_adjacent(self, node_id1, node_id2): + """Test for adjacency between node1 and node2 + + Args: + node_id1(int): Index of first node + node_id2(int): Index of second node + """ + if node_id1 >= self.node_count() or node_id2 >= self.node_count(): + raise IndexError("Node ID out of range for graph") + + return (self.node_uids[node_id1] in self.neighbors[node_id2]) + + def get_neighbors(self, node_id): + """Return list of neighbors of the specified node + + Args: + node_id: Index of node to examine the neighbors of + + Returns(list): List of current node IDs that are neighbors of node_id. + """ + if node_id >= self.node_count(): + raise IndexError("Node ID out of range for graph") + return [self.uid_to_index[i] for i in self.neighbors[node_id]] + + def shortest_path(self, node_id1, node_id2): + """Find the shortest path between node1 and node2 on the graph + + Args: + node_id1(int): Index of first node + node_id2(int): Index of second node + + Returns(list): List of nodes from node_id1 to node_id2 that constitute + the shortest possible path in the graph between those two nodes. + """ + if node_id1 >= self.node_count() or node_id2 >= self.node_count(): + raise IndexError("Node ID out of range for graph") + + # Treat special case of equal inputs + if node_id1 == node_id2: + return [node_id1] + + # Initialize two arrays for backtracking the shortest path found + previous = [None] * self.node_count() # Tracks the moves + distances = [None] * self.node_count() # Records distance to nodes + distances[node_id1] = 0 + + node_queue = queue.LifoQueue() + node_queue.put(node_id1) + + while not node_queue.empty(): + next_id = node_queue.get() + if next_id == node_id2: + break # Success + new_distance = distances[next_id] + 1 + # On each iteration, check if going to a neighbor node was + # shortest path to that node, if so, add that node to the queue + # and record distance and path that is being taken. + for uid in self.neighbors[next_id]: + if ((distances[self.uid_to_index[uid]] is None) or + (distances[self.uid_to_index[uid]] > new_distance)): + distances[self.uid_to_index[uid]] = new_distance + previous[self.uid_to_index[uid]] = next_id + node_queue.put(self.uid_to_index[uid]) + + # Backtrack to get path + if previous[node_id2] is None: + raise IndexError("Unable to find target node, " + "possibly disconnected graph.") + + path = [] + next = node_id2 + while next is not None: + path.append(next) + next = previous[next] + # Return reversed backtrack path to get in order from node 1 to 2 + return path[::-1] diff --git a/src/fermilib/utils/_graph_test.py b/src/fermilib/utils/_graph_test.py new file mode 100644 index 0000000..c4acfa4 --- /dev/null +++ b/src/fermilib/utils/_graph_test.py @@ -0,0 +1,104 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from __future__ import absolute_import + +import unittest + +from fermilib.utils import Graph, Node + + +class GraphTest(unittest.TestCase): + + def test_two_node_graph(self): + """Build a graph with two nodes and one edge""" + two_node = Graph() + two_node.add_node() + two_node.add_node() + two_node.add_edge(0, 1) + + # Check to see if it built normally + self.assertEqual(two_node.neighbors, + [{1}, {0}]) + self.assertEqual(two_node.node_uids, + [0, 1]) + + self.assertEqual(two_node.uid_to_index[0], + 0) + self.assertEqual(two_node.uid_to_index[1], + 1) + self.assertEqual(two_node.node_count(), 2) + + # Now remove an edge, then a node, see if it + two_node.remove_edge(0, 1) + self.assertEqual(two_node.neighbors, + [set([]), set([])]) + two_node.remove_node(0) + self.assertEqual(two_node.node_uids, [1]) + + # Add a node and edge back, verify indexing + two_node.add_node() + two_node.add_edge(0, 1) + self.assertEqual(two_node.neighbors, + [{2}, {1}]) + self.assertEqual(two_node.uid_to_index[1], 0) + self.assertEqual(two_node.uid_to_index[2], 1) + self.assertEqual(two_node.node_count(), 2) + + self.assertEqual(two_node.get_neighbors(0), [1]) + self.assertEqual(two_node.get_neighbors(1), [0]) + + self.assertTrue(two_node.is_adjacent(0, 1)) + + # Check path between nodes + self.assertEqual(two_node.shortest_path(0, 1), [0, 1]) + self.assertEqual(two_node.shortest_path(1, 0), [1, 0]) + self.assertEqual(two_node.shortest_path(0, 0), [0]) + + # Remove just the node and check neighbors + two_node.remove_node(1) + self.assertEqual(two_node.neighbors, + [set([])]) + + def test_ring(self): + """Build a ring of 8 nodes and test path finding""" + eight_node = Graph() + for i in range(8): + eight_node.add_node(Node(value=i)) + for i in range(8): + eight_node.add_edge(i, (i + 1) % 8) + self.assertEqual(eight_node.neighbors, + [{1, 7}, {0, 2}, {1, 3}, {2, 4}, + {3, 5}, {4, 6}, {5, 7}, {0, 6}]) + + self.assertEqual(eight_node.shortest_path(0, 4), + [0, 7, 6, 5, 4]) + self.assertEqual(eight_node.shortest_path(0, 6), + [0, 7, 6]) + self.assertTrue(eight_node.is_adjacent(0, 7)) + + # Look for node with value 6 and value not present + found_index = eight_node.find_index(6) + self.assertEqual(found_index, 6) + found_index = eight_node.find_index(10) + self.assertIsNone(found_index) + + # Make a hole in the ring and check new distance + eight_node.remove_node(7) + self.assertTrue(eight_node.is_adjacent(0, 1)) + self.assertFalse(eight_node.is_adjacent(0, 6)) + self.assertEqual(eight_node.shortest_path(0, 6), + [0, 1, 2, 3, 4, 5, 6]) + +# Run test. +if __name__ == '__main__': + unittest.main() diff --git a/src/fermilib/utils/_grid.py b/src/fermilib/utils/_grid.py new file mode 100644 index 0000000..872242d --- /dev/null +++ b/src/fermilib/utils/_grid.py @@ -0,0 +1,67 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from __future__ import absolute_import + +import itertools + + +class Grid: + """ + A multi-dimensional grid of points. + """ + + def __init__(self, dimensions, length, scale): + """ + Args: + dimensions (int): The number of dimensions the grid lives in. + length (int): The number of points along each grid axis. + scale (float): The total length of each grid dimension. + """ + if not isinstance(dimensions, int) or dimensions < 0: + raise ValueError( + 'dimensions must be a non-negative int but was {} {}'.format( + type(dimensions), repr(dimensions))) + if not isinstance(length, int) or length < 0: + raise ValueError( + 'length must be a non-negative int but was {} {}'.format( + type(length), repr(length))) + if not isinstance(scale, float) or not scale > 0: + raise ValueError( + 'scale must be a positive float but was {} {}'.format( + type(scale), repr(scale))) + + self.dimensions = dimensions + self.length = length + self.scale = scale + + def volume_scale(self): + """ + Returns: + float: The volume of a length-scale hypercube within the grid. + """ + return self.scale ** float(self.dimensions) + + def num_points(self): + """ + Returns: + int: The number of points in the grid. + """ + return self.length ** self.dimensions + + def all_points_indices(self): + """ + Returns: + iterable[tuple[int]]: + The index-coordinate tuple of each point in the grid. + """ + return itertools.product(range(self.length), repeat=self.dimensions) diff --git a/src/fermilib/utils/_grid_test.py b/src/fermilib/utils/_grid_test.py new file mode 100644 index 0000000..1222290 --- /dev/null +++ b/src/fermilib/utils/_grid_test.py @@ -0,0 +1,64 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from __future__ import absolute_import + +import unittest + +from fermilib.utils import Grid + + +class GridTest(unittest.TestCase): + + def test_preconditions(self): + nan = float('nan') + + # No exception + _ = Grid(dimensions=0, length=0, scale=1.0) + _ = Grid(dimensions=1, length=1, scale=1.0) + _ = Grid(dimensions=2, length=3, scale=0.01) + _ = Grid(dimensions=234, length=345, scale=456.0) + + with self.assertRaises(ValueError): + _ = Grid(dimensions=1, length=1, scale=1) + with self.assertRaises(ValueError): + _ = Grid(dimensions=1, length=1, scale=0.0) + with self.assertRaises(ValueError): + _ = Grid(dimensions=1, length=1, scale=-1.0) + with self.assertRaises(ValueError): + _ = Grid(dimensions=1, length=1, scale=nan) + + with self.assertRaises(ValueError): + _ = Grid(dimensions=1, length=-1, scale=1.0) + with self.assertRaises(ValueError): + _ = Grid(dimensions=-1, length=1, scale=1.0) + + def test_properties(self): + g = Grid(dimensions=2, length=3, scale=5.0) + self.assertEqual(g.num_points(), 9) + self.assertEqual(g.volume_scale(), 25) + self.assertEqual(list(g.all_points_indices()), [ + (0, 0), + (0, 1), + (0, 2), + (1, 0), + (1, 1), + (1, 2), + (2, 0), + (2, 1), + (2, 2), + ]) + + +# Run test. +if __name__ == '__main__': + unittest.main() diff --git a/src/fermilib/utils/_jellium.py b/src/fermilib/utils/_jellium.py index 442886e..0782dad 100644 --- a/src/fermilib/utils/_jellium.py +++ b/src/fermilib/utils/_jellium.py @@ -13,24 +13,22 @@ """This module constructs Hamiltonians for the uniform electron gas.""" from __future__ import absolute_import -import itertools import numpy +from projectq.ops import QubitOperator from fermilib.ops import FermionOperator -from projectq.ops import QubitOperator - # Exceptions. class OrbitalSpecificationError(Exception): pass -def orbital_id(grid_length, grid_coordinates, spin=None): +def orbital_id(grid, grid_coordinates, spin=None): """Return the tensor factor of a orbital with given coordinates and spin. Args: - grid_length: Int, the number of points in one dimension of the grid. + grid (Grid): The discretization to use. grid_coordinates: List or tuple of ints giving coordinates of grid element. Acceptable to provide an int (instead of tuple or list) for 1D case. @@ -38,10 +36,8 @@ def orbital_id(grid_length, grid_coordinates, spin=None): If None, assume spinless model. Returns: - tensor_factor: tensor factor associated with provided orbital label. - - Raises: - OrbitalSpecificiationError: Invalid orbital coordinates provided. + tensor_factor (int): + tensor factor associated with provided orbital label. """ # Initialize. if isinstance(grid_coordinates, int): @@ -52,8 +48,8 @@ def orbital_id(grid_length, grid_coordinates, spin=None): for dimension, grid_coordinate in enumerate(grid_coordinates): # Make sure coordinate is an integer in the correct bounds. - if isinstance(grid_coordinate, int) and grid_coordinate < grid_length: - tensor_factor += grid_coordinate * (grid_length ** dimension) + if isinstance(grid_coordinate, int) and grid_coordinate < grid.length: + tensor_factor += grid_coordinate * (grid.length ** dimension) else: # Raise for invalid model. @@ -69,17 +65,17 @@ def orbital_id(grid_length, grid_coordinates, spin=None): return tensor_factor -def grid_indices(qubit_id, n_dimensions, grid_length, spinless): +def grid_indices(qubit_id, grid, spinless): """This function is the inverse of orbital_id. Args: - qubit_id: The tensor factor to map to grid indices. - n_dimensions: An int giving the number of dimensions for the model. - grid_length (int): The number of points in one dimension of the grid. + qubit_id (int): The tensor factor to map to grid indices. + grid (Grid): The discretization to use. spinless (bool): Whether to use the spinless model or not. Returns: - grid_indices: The location of the qubit on the grid. + grid_indices (numpy.ndarray[int]): + The location of the qubit on the grid. """ # Remove spin degree of freedom. orbital_id = qubit_id @@ -90,94 +86,72 @@ def grid_indices(qubit_id, n_dimensions, grid_length, spinless): # Get grid indices. grid_indices = [] - for dimension in range(n_dimensions): - remainder = orbital_id % (grid_length ** (dimension + 1)) - grid_index = remainder // (grid_length ** dimension) + for dimension in range(grid.dimensions): + remainder = orbital_id % (grid.length ** (dimension + 1)) + grid_index = remainder // (grid.length ** dimension) grid_indices += [grid_index] return grid_indices -def position_vector(position_indices, grid_length, length_scale): +def position_vector(position_indices, grid): """Given grid point coordinate, return position vector with dimensions. Args: - position_indices: List or tuple of integers giving grid point - coordinate. Allowed values are ints in [0, grid_length). - grid_length (int): The number of points in one dimension of the grid. - length_scale (float): The real space length of a box dimension. + position_indices (int|iterable[int]): + List or tuple of integers giving grid point coordinate. + Allowed values are ints in [0, grid_length). + grid (Grid): The discretization to use. Returns: - position_vector: A numpy array giving the position vector with - dimensions. - - Raises: - orbitalSpecificationError: Position indices must be integers - in [0, grid_length). + position_vector (numpy.ndarray[float]) """ # Raise exceptions. if isinstance(position_indices, int): position_indices = [position_indices] - if (not isinstance(grid_length, int) or - max(position_indices) >= grid_length or - min(position_indices) < 0.): - raise orbitalSpecificationError( + if not all(0 <= e < grid.length for e in position_indices): + raise OrbitalSpecificationError( 'Position indices must be integers in [0, grid_length).') # Compute position vector. - shift = float(grid_length - 1) / 2. - adjusted_vector = numpy.array(position_indices, float) - shift - position_vector = length_scale * adjusted_vector / float(grid_length) - return position_vector + adjusted_vector = numpy.array(position_indices, float) - grid.length // 2 + return grid.scale * adjusted_vector / float(grid.length) -def momentum_vector(momentum_indices, grid_length, length_scale): +def momentum_vector(momentum_indices, grid): """Given grid point coordinate, return momentum vector with dimensions. Args: momentum_indices: List or tuple of integers giving momentum indices. Allowed values are ints in [0, grid_length). - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. + grid (Grid): The discretization to use. Returns: momentum_vector: A numpy array giving the momentum vector with dimensions. - - Raises: - OrbitalSpecificationError: Momentum indices must be integers - in [0, grid_length). """ # Raise exceptions. if isinstance(momentum_indices, int): momentum_indices = [momentum_indices] - if (not isinstance(grid_length, int) or - max(momentum_indices) >= grid_length or - min(momentum_indices) < 0.): + if not all(0 <= e < grid.length for e in momentum_indices): raise OrbitalSpecificationError( 'Momentum indices must be integers in [0, grid_length).') # Compute momentum vector. - shift = float(grid_length - 1) / 2. - adjusted_vector = numpy.array(momentum_indices, float) - shift - momentum_vector = 2. * numpy.pi * adjusted_vector / length_scale - return momentum_vector + adjusted_vector = numpy.array(momentum_indices, float) - grid.length // 2 + return 2. * numpy.pi * adjusted_vector / grid.scale -def momentum_kinetic_operator(n_dimensions, grid_length, - length_scale, spinless=False): +def momentum_kinetic_operator(grid, spinless=False): """Return the kinetic energy operator in momentum second quantization. Args: - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - spinless: Bool, whether to use the spinless model or not. + grid (fermilib.utils.Grid): The discretization to use. + spinless (bool): Whether to use the spinless model or not. Returns: - operator: An instance of the FermionOperator class. + FermionOperator: The kinetic momentum operator. """ # Initialize. - n_points = grid_length ** n_dimensions operator = FermionOperator() if spinless: spins = [None] @@ -185,14 +159,13 @@ def momentum_kinetic_operator(n_dimensions, grid_length, spins = [0, 1] # Loop once through all plane waves. - for grid_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector(grid_indices, grid_length, length_scale) + for momenta_indices in grid.all_points_indices(): + momenta = momentum_vector(momenta_indices, grid) coefficient = momenta.dot(momenta) / 2. # Loop over spins. for spin in spins: - orbital = orbital_id(grid_length, grid_indices, spin) + orbital = orbital_id(grid, momenta_indices, spin) # Add interaction term. operators = ((orbital, 1), (orbital, 0)) @@ -201,81 +174,58 @@ def momentum_kinetic_operator(n_dimensions, grid_length, return operator -def momentum_potential_operator(n_dimensions, grid_length, - length_scale, spinless=False): +def momentum_potential_operator(grid, spinless=False): """Return the potential operator in momentum second quantization. Args: - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - spinless: Boole, whether to use the spinless model or not. + grid (Grid): The discretization to use. + spinless (bool): Whether to use the spinless model or not. Returns: - operator: An instance of the FermionOperator class. - - Raises: - OrbitalSpecificationError: 'Must use an odd number of momentum modes.' + operator (FermionOperator) """ - # Make sure number of orbitals is odd. - if not (grid_length % 2): - raise OrbitalSpecificationError( - 'Must use an odd number of momentum modes.') - # Initialize. - n_points = grid_length ** n_dimensions - volume = length_scale ** float(n_dimensions) + volume = grid.volume_scale() prefactor = 2. * numpy.pi / volume operator = FermionOperator((), 0.0) - if spinless: - spins = [None] - else: - spins = [0, 1] + spins = [None] if spinless else [0, 1] # Loop once through all plane waves. - for omega_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - shifted_omega_indices = [index - grid_length // 2 for + for omega_indices in grid.all_points_indices(): + shifted_omega_indices = [index - grid.length // 2 for index in omega_indices] # Get the momenta vectors. - omega_momenta = momentum_vector( - omega_indices, grid_length, length_scale) + omega_momenta = momentum_vector(omega_indices, grid) # Skip if omega momentum is zero. if not omega_momenta.any(): continue # Compute coefficient. - coefficient = prefactor / \ - omega_momenta.dot(omega_momenta) + coefficient = prefactor / omega_momenta.dot(omega_momenta) - for grid_indices_a in itertools.product(range(grid_length), - repeat=n_dimensions): + for grid_indices_a in grid.all_points_indices(): shifted_indices_d = [ - (grid_indices_a[i] - shifted_omega_indices[i]) % - grid_length for i in range(n_dimensions)] - for grid_indices_b in itertools.product(range(grid_length), - repeat=n_dimensions): + (grid_indices_a[i] - shifted_omega_indices[i]) % grid.length + for i in range(grid.dimensions)] + for grid_indices_b in grid.all_points_indices(): shifted_indices_c = [ (grid_indices_b[i] + shifted_omega_indices[i]) % - grid_length for i in range(n_dimensions)] + grid.length + for i in range(grid.dimensions)] # Loop over spins. for spin_a in spins: - orbital_a = orbital_id( - grid_length, grid_indices_a, spin_a) - orbital_d = orbital_id( - grid_length, shifted_indices_d, spin_a) + orbital_a = orbital_id(grid, grid_indices_a, spin_a) + orbital_d = orbital_id(grid, shifted_indices_d, spin_a) for spin_b in spins: - orbital_b = orbital_id( - grid_length, grid_indices_b, spin_b) - orbital_c = orbital_id( - grid_length, shifted_indices_c, spin_b) + orbital_b = orbital_id(grid, grid_indices_b, spin_b) + orbital_c = orbital_id(grid, shifted_indices_c, spin_b) # Add interaction term. - if (orbital_a != orbital_b) and \ - (orbital_c != orbital_d): + if ((orbital_a != orbital_b) and + (orbital_c != orbital_d)): operators = ((orbital_a, 1), (orbital_b, 1), (orbital_c, 0), (orbital_d, 0)) operator += FermionOperator(operators, coefficient) @@ -284,44 +234,32 @@ def momentum_potential_operator(n_dimensions, grid_length, return operator -def position_kinetic_operator(n_dimensions, grid_length, - length_scale, spinless=False): +def position_kinetic_operator(grid, spinless=False): """Return the kinetic operator in position space second quantization. Args: - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - spinless: Bool, whether to use the spinless model or not. + grid (Grid): The discretization to use. + spinless (bool): Whether to use the spinless model or not. Returns: - operator: An instance of the FermionOperator class. + operator (FermionOperator) """ # Initialize. - n_points = grid_length ** n_dimensions + n_points = grid.num_points() operator = FermionOperator() - if spinless: - spins = [None] - else: - spins = [0, 1] + spins = [None] if spinless else [0, 1] # Loop once through all lattice sites. - for grid_indices_a in itertools.product(range(grid_length), - repeat=n_dimensions): - coordinates_a = position_vector( - grid_indices_a, grid_length, length_scale) - for grid_indices_b in itertools.product(range(grid_length), - repeat=n_dimensions): - coordinates_b = position_vector( - grid_indices_b, grid_length, length_scale) + for grid_indices_a in grid.all_points_indices(): + coordinates_a = position_vector(grid_indices_a, grid) + for grid_indices_b in grid.all_points_indices(): + coordinates_b = position_vector(grid_indices_b, grid) differences = coordinates_b - coordinates_a # Compute coefficient. coefficient = 0. - for momenta_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector( - momenta_indices, grid_length, length_scale) + for momenta_indices in grid.all_points_indices(): + momenta = momentum_vector(momenta_indices, grid) if momenta.any(): coefficient += ( numpy.cos(momenta.dot(differences)) * @@ -329,8 +267,8 @@ def position_kinetic_operator(n_dimensions, grid_length, # Loop over spins and identify interacting orbitals. for spin in spins: - orbital_a = orbital_id(grid_length, grid_indices_a, spin) - orbital_b = orbital_id(grid_length, grid_indices_b, spin) + orbital_a = orbital_id(grid, grid_indices_a, spin) + orbital_b = orbital_id(grid, grid_indices_b, spin) # Add interaction term. operators = ((orbital_a, 1), (orbital_b, 0)) @@ -340,46 +278,33 @@ def position_kinetic_operator(n_dimensions, grid_length, return operator -def position_potential_operator(n_dimensions, grid_length, - length_scale, spinless=False): +def position_potential_operator(grid, spinless=False): """Return the potential operator in position space second quantization. Args: - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - spinless: Boole, whether to use the spinless model or not. + grid (Grid): The discretization to use. + spinless (bool): Whether to use the spinless model or not. Returns: - operator: An instance of the FermionOperator class. + operator (FermionOperator) """ # Initialize. - n_points = grid_length ** n_dimensions - volume = length_scale ** float(n_dimensions) + volume = grid.volume_scale() prefactor = 2. * numpy.pi / volume operator = FermionOperator() - if spinless: - spins = [None] - else: - spins = [0, 1] + spins = [None] if spinless else [0, 1] # Loop once through all lattice sites. - for grid_indices_a in itertools.product(range(grid_length), - repeat=n_dimensions): - coordinates_a = position_vector( - grid_indices_a, grid_length, length_scale) - for grid_indices_b in itertools.product(range(grid_length), - repeat=n_dimensions): - coordinates_b = position_vector( - grid_indices_b, grid_length, length_scale) + for grid_indices_a in grid.all_points_indices(): + coordinates_a = position_vector(grid_indices_a, grid) + for grid_indices_b in grid.all_points_indices(): + coordinates_b = position_vector(grid_indices_b, grid) differences = coordinates_b - coordinates_a # Compute coefficient. coefficient = 0. - for momenta_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector( - momenta_indices, grid_length, length_scale) + for momenta_indices in grid.all_points_indices(): + momenta = momentum_vector(momenta_indices, grid) if momenta.any(): coefficient += ( prefactor * numpy.cos(momenta.dot(differences)) / @@ -387,9 +312,9 @@ def position_potential_operator(n_dimensions, grid_length, # Loop over spins and identify interacting orbitals. for spin_a in spins: - orbital_a = orbital_id(grid_length, grid_indices_a, spin_a) + orbital_a = orbital_id(grid, grid_indices_a, spin_a) for spin_b in spins: - orbital_b = orbital_id(grid_length, grid_indices_b, spin_b) + orbital_b = orbital_id(grid, grid_indices_b, spin_b) # Add interaction term. if orbital_a != orbital_b: @@ -400,71 +325,54 @@ def position_potential_operator(n_dimensions, grid_length, return operator -def jellium_model(n_dimensions, grid_length, length_scale, - spinless=False, momentum_space=True): +def jellium_model(grid, spinless=False, momentum_space=True): """Return jellium Hamiltonian as FermionOperator class. Args: - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - spinless: Bool, whether to use the spinless model or not. - momentum_space: Boole, whether to return in momentum space (True) + grid (fermilib.utils.Grid): The discretization to use. + spinless (bool): Whether to use the spinless model or not. + momentum_space (bool): Whether to return in momentum space (True) or position space (False). Returns: - hamiltonian: An instance of the FermionOperator class. + FermionOperator: The Hamiltonian of the model. """ + if grid.length & 1 == 0 and grid.length & (grid.length - 1): + raise OrbitalSpecificationError( + 'Must use an odd number or a power of 2 for momentum modes.') + if momentum_space: - hamiltonian = momentum_kinetic_operator(n_dimensions, - grid_length, - length_scale, - spinless) - hamiltonian += momentum_potential_operator(n_dimensions, - grid_length, - length_scale, - spinless) + hamiltonian = momentum_kinetic_operator(grid, spinless) + hamiltonian += momentum_potential_operator(grid, spinless) else: - hamiltonian = position_kinetic_operator(n_dimensions, - grid_length, - length_scale, - spinless) - hamiltonian += position_potential_operator(n_dimensions, - grid_length, - length_scale, - spinless) + hamiltonian = position_kinetic_operator(grid, spinless) + hamiltonian += position_potential_operator(grid, spinless) return hamiltonian -def jordan_wigner_position_jellium(n_dimensions, grid_length, - length_scale, spinless=False): +def jordan_wigner_position_jellium(grid, spinless=False): """Return the position space jellium Hamiltonian as QubitOperator. Args: - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - spinless: Bool, whether to use the spinless model or not. + grid (Grid): The discretization to use. + spinless (bool): Whether to use the spinless model or not. Returns: - hamiltonian: An instance of the QubitOperator class. + hamiltonian (QubitOperator) """ # Initialize. - n_orbitals = grid_length ** n_dimensions - volume = length_scale ** float(n_dimensions) + n_orbitals = grid.num_points() + volume = grid.volume_scale() if spinless: - spins = [None] n_qubits = n_orbitals else: - spins = [0, 1] n_qubits = 2 * n_orbitals hamiltonian = QubitOperator() # Compute the identity coefficient. identity_coefficient = 0. - for k_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector(k_indices, grid_length, length_scale) + for k_indices in grid.all_points_indices(): + momenta = momentum_vector(k_indices, grid) if momenta.any(): identity_coefficient += momenta.dot(momenta) / 2. identity_coefficient -= (numpy.pi * float(n_orbitals) / @@ -478,9 +386,8 @@ def jordan_wigner_position_jellium(n_dimensions, grid_length, # Compute coefficient of local Z terms. z_coefficient = 0. - for k_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector(k_indices, grid_length, length_scale) + for k_indices in grid.all_points_indices(): + momenta = momentum_vector(k_indices, grid) if momenta.any(): z_coefficient += numpy.pi / (momenta.dot(momenta) * volume) z_coefficient -= momenta.dot(momenta) / (4. * float(n_orbitals)) @@ -493,19 +400,18 @@ def jordan_wigner_position_jellium(n_dimensions, grid_length, # Add ZZ terms. prefactor = numpy.pi / volume for p in range(n_qubits): - index_p = grid_indices(p, n_dimensions, grid_length, spinless) - position_p = position_vector(index_p, grid_length, length_scale) + index_p = grid_indices(p, grid, spinless) + position_p = position_vector(index_p, grid) for q in range(p + 1, n_qubits): - index_q = grid_indices(q, n_dimensions, grid_length, spinless) - position_q = position_vector(index_q, grid_length, length_scale) + index_q = grid_indices(q, grid, spinless) + position_q = position_vector(index_q, grid) differences = position_p - position_q # Loop through momenta. zpzq_coefficient = 0. - for k_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector(k_indices, grid_length, length_scale) + for k_indices in grid.all_points_indices(): + momenta = momentum_vector(k_indices, grid) if momenta.any(): zpzq_coefficient += prefactor * numpy.cos( momenta.dot(differences)) / momenta.dot(momenta) @@ -517,25 +423,25 @@ def jordan_wigner_position_jellium(n_dimensions, grid_length, # Add XZX + YZY terms. prefactor = .25 / float(n_orbitals) for p in range(n_qubits): - index_p = grid_indices(p, n_dimensions, grid_length, spinless) - position_p = position_vector(index_p, grid_length, length_scale) + index_p = grid_indices(p, grid, spinless) + position_p = position_vector(index_p, grid) for q in range(p + 1, n_qubits): if not spinless and (p + q) % 2: continue - index_q = grid_indices(q, n_dimensions, grid_length, spinless) - position_q = position_vector(index_q, grid_length, length_scale) + index_q = grid_indices(q, grid, spinless) + position_q = position_vector(index_q, grid) differences = position_p - position_q # Loop through momenta. term_coefficient = 0. - for k_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector(k_indices, grid_length, length_scale) + for k_indices in grid.all_points_indices(): + momenta = momentum_vector(k_indices, grid) if momenta.any(): - term_coefficient += prefactor * momenta.dot(momenta) * \ - numpy.cos(momenta.dot(differences)) + term_coefficient += (prefactor * + momenta.dot(momenta) * + numpy.cos(momenta.dot(differences))) # Add term. z_string = tuple((i, 'Z') for i in range(p + 1, q)) diff --git a/src/fermilib/utils/_jellium_test.py b/src/fermilib/utils/_jellium_test.py index baffc73..632311b 100644 --- a/src/fermilib/utils/_jellium_test.py +++ b/src/fermilib/utils/_jellium_test.py @@ -12,13 +12,24 @@ from __future__ import absolute_import -import itertools -import numpy import unittest +import numpy + from fermilib.transforms import jordan_wigner -from fermilib.utils._jellium import * -from fermilib.utils import count_qubits, eigenspectrum +from fermilib.utils import count_qubits, eigenspectrum, Grid +from fermilib.utils._jellium import ( + jellium_model, + jordan_wigner_position_jellium, + momentum_kinetic_operator, + momentum_potential_operator, + momentum_vector, + orbital_id, + OrbitalSpecificationError, + position_kinetic_operator, + position_potential_operator, + position_vector, +) class JelliumTest(unittest.TestCase): @@ -26,89 +37,76 @@ class JelliumTest(unittest.TestCase): def test_orbital_id(self): # Test in 1D with spin. - grid_length = 5 + grid = Grid(dimensions=1, length=5, scale=1.0) input_coords = [0, 1, 2, 3, 4] tensor_factors_up = [1, 3, 5, 7, 9] tensor_factors_down = [0, 2, 4, 6, 8] - test_output_up = [orbital_id( - grid_length, i, 1) for i in input_coords] - test_output_down = [orbital_id( - grid_length, i, 0) for i in input_coords] + test_output_up = [orbital_id(grid, i, 1) for i in input_coords] + test_output_down = [orbital_id(grid, i, 0) for i in input_coords] self.assertEqual(test_output_up, tensor_factors_up) self.assertEqual(test_output_down, tensor_factors_down) with self.assertRaises(OrbitalSpecificationError): - orbital_id(5, 6, 1) + orbital_id(grid, 6, 1) # Test in 2D without spin. - grid_length = 3 + grid = Grid(dimensions=2, length=3, scale=1.0) input_coords = [(0, 0), (0, 1), (1, 2)] tensor_factors = [0, 3, 7] - test_output = [orbital_id( - grid_length, i) for i in input_coords] + test_output = [orbital_id(grid, i) for i in input_coords] self.assertEqual(test_output, tensor_factors) def test_position_vector(self): # Test in 1D. - grid_length = 4 - length_scale = 4. - test_output = [position_vector(i, grid_length, length_scale) - for i in range(grid_length)] - correct_output = [-1.5, -.5, .5, 1.5] + grid = Grid(dimensions=1, length=4, scale=4.) + test_output = [position_vector(i, grid) + for i in range(grid.length)] + correct_output = [-2, -1, 0, 1] self.assertEqual(correct_output, test_output) - grid_length = 11 - length_scale = 2. * numpy.pi - for i in range(grid_length): + grid = Grid(dimensions=1, length=11, scale=2. * numpy.pi) + for i in range(grid.length): self.assertAlmostEqual( - -position_vector(i, grid_length, length_scale), - position_vector( - grid_length - i - 1, grid_length, length_scale)) + -position_vector(i, grid), + position_vector(grid.length - i - 1, grid)) # Test in 2D. - grid_length = 3 - length_scale = 3. + grid = Grid(dimensions=2, length=3, scale=3.) test_input = [] test_output = [] for i in range(3): for j in range(3): test_input += [(i, j)] - test_output += [position_vector( - (i, j), grid_length, length_scale)] + test_output += [position_vector((i, j), grid)] correct_output = numpy.array([[-1., -1.], [-1., 0.], [-1., 1.], [0., -1.], [0., 0.], [0., 1.], [1., -1.], [1., 0.], [1., 1.]]) self.assertAlmostEqual(0., numpy.amax(test_output - correct_output)) def test_momentum_vector(self): - grid_length = 3 - length_scale = 2. * numpy.pi - test_output = [momentum_vector(i, grid_length, length_scale) - for i in range(grid_length)] + grid = Grid(dimensions=1, length=3, scale=2. * numpy.pi) + test_output = [momentum_vector(i, grid) + for i in range(grid.length)] correct_output = [-1., 0, 1.] self.assertEqual(correct_output, test_output) - grid_length = 11 - length_scale = 2. * numpy.pi - for i in range(grid_length): + grid = Grid(dimensions=1, length=11, scale=2. * numpy.pi) + for i in range(grid.length): self.assertAlmostEqual( - -momentum_vector(i, grid_length, length_scale), - momentum_vector( - grid_length - i - 1, grid_length, length_scale)) + -momentum_vector(i, grid), + momentum_vector(grid.length - i - 1, grid)) # Test in 2D. - grid_length = 3 - length_scale = 2. * numpy.pi + grid = Grid(dimensions=2, length=3, scale=2. * numpy.pi) test_input = [] test_output = [] for i in range(3): for j in range(3): test_input += [(i, j)] - test_output += [momentum_vector( - (i, j), grid_length, length_scale)] + test_output += [momentum_vector((i, j), grid)] correct_output = numpy.array([[-1, -1], [-1, 0], [-1, 1], [0, -1], [0, 0], [0, 1], [1, -1], [1, 0], [1, 1]]) @@ -117,14 +115,10 @@ def test_momentum_vector(self): def test_kinetic_integration(self): # Compute kinetic energy operator in both momentum and position space. - n_dimensions = 2 - grid_length = 2 - length_scale = 3. - spinless = 0 - momentum_kinetic = momentum_kinetic_operator( - n_dimensions, grid_length, length_scale, spinless) - position_kinetic = position_kinetic_operator( - n_dimensions, grid_length, length_scale, spinless) + grid = Grid(dimensions=2, length=2, scale=3.) + spinless = False + momentum_kinetic = momentum_kinetic_operator(grid, spinless) + position_kinetic = position_kinetic_operator(grid, spinless) # Diagonalize and confirm the same energy. jw_momentum = jordan_wigner(momentum_kinetic) @@ -139,16 +133,11 @@ def test_kinetic_integration(self): def test_potential_integration(self): - # Compute potential energy operator in both momentum and position - # space. - n_dimensions = 2 - grid_length = 3 - length_scale = 2. + # Compute potential energy operator in momentum and position space. + grid = Grid(dimensions=2, length=3, scale=2.) spinless = 1 - momentum_potential = momentum_potential_operator( - n_dimensions, grid_length, length_scale, spinless) - position_potential = position_potential_operator( - n_dimensions, grid_length, length_scale, spinless) + momentum_potential = momentum_potential_operator(grid, spinless) + position_potential = position_potential_operator(grid, spinless) # Diagonalize and confirm the same energy. jw_momentum = jordan_wigner(momentum_potential) @@ -164,14 +153,10 @@ def test_potential_integration(self): def test_model_integration(self): # Compute Hamiltonian in both momentum and position space. - n_dimensions = 2 - grid_length = 3 - length_scale = 1. - spinless = 1 - momentum_hamiltonian = jellium_model( - n_dimensions, grid_length, length_scale, spinless, 1) - position_hamiltonian = jellium_model( - n_dimensions, grid_length, length_scale, spinless, 0) + grid = Grid(dimensions=2, length=3, scale=1.0) + spinless = True + momentum_hamiltonian = jellium_model(grid, spinless, True) + position_hamiltonian = jellium_model(grid, spinless, False) # Diagonalize and confirm the same energy. jw_momentum = jordan_wigner(momentum_hamiltonian) @@ -187,22 +172,18 @@ def test_model_integration(self): def test_coefficients(self): # Test that the coefficients post-JW transform are as claimed in paper. - n_dimensions = 2 - grid_length = 3 - length_scale = 2. + grid = Grid(dimensions=2, length=3, scale=2.) spinless = 1 - n_orbitals = grid_length ** n_dimensions + n_orbitals = grid.num_points() n_qubits = (2 ** (1 - spinless)) * n_orbitals - volume = length_scale ** n_dimensions + volume = grid.volume_scale() # Kinetic operator. - kinetic = position_kinetic_operator( - n_dimensions, grid_length, length_scale, spinless) + kinetic = position_kinetic_operator(grid, spinless) qubit_kinetic = jordan_wigner(kinetic) # Potential operator. - potential = position_potential_operator( - n_dimensions, grid_length, length_scale, spinless) + potential = position_potential_operator(grid, spinless) qubit_potential = jordan_wigner(potential) # Check identity. @@ -212,10 +193,8 @@ def test_coefficients(self): paper_kinetic_coefficient = 0. paper_potential_coefficient = 0. - for indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector( - indices, grid_length, length_scale) + for indices in grid.all_points_indices(): + momenta = momentum_vector(indices, grid) paper_kinetic_coefficient += float( n_qubits) * momenta.dot(momenta) / float(4. * n_orbitals) @@ -237,10 +216,8 @@ def test_coefficients(self): paper_kinetic_coefficient = 0. paper_potential_coefficient = 0. - for indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector( - indices, grid_length, length_scale) + for indices in grid.all_points_indices(): + momenta = momentum_vector(indices, grid) paper_kinetic_coefficient -= momenta.dot( momenta) / float(4. * n_orbitals) @@ -260,27 +237,21 @@ def test_coefficients(self): else: spins = [0, 1] - for indices_a in itertools.product(range(grid_length), - repeat=n_dimensions): - for indices_b in itertools.product(range(grid_length), - repeat=n_dimensions): + for indices_a in grid.all_points_indices(): + for indices_b in grid.all_points_indices(): paper_kinetic_coefficient = 0. paper_potential_coefficient = 0. - position_a = position_vector( - indices_a, grid_length, length_scale) - position_b = position_vector( - indices_b, grid_length, length_scale) + position_a = position_vector(indices_a, grid) + position_b = position_vector(indices_b, grid) differences = position_b - position_a for spin_a in spins: for spin_b in spins: - p = orbital_id( - grid_length, indices_a, spin_a) - q = orbital_id( - grid_length, indices_b, spin_b) + p = orbital_id(grid, indices_a, spin_a) + q = orbital_id(grid, indices_b, spin_b) if p == q: continue @@ -291,11 +262,8 @@ def test_coefficients(self): else: potential_coefficient = 0. - for indices_c in \ - itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector( - indices_c, grid_length, length_scale) + for indices_c in grid.all_points_indices(): + momenta = momentum_vector(indices_c, grid) if momenta.any(): potential_contribution = numpy.pi * numpy.cos( @@ -309,19 +277,15 @@ def test_coefficients(self): def test_jordan_wigner_position_jellium(self): # Parameters. - n_dimensions = 2 - grid_length = 3 - length_scale = 1. - spinless = 1 + grid = Grid(dimensions=2, length=3, scale=1.) + spinless = True # Compute fermionic Hamiltonian. - fermion_hamiltonian = jellium_model( - n_dimensions, grid_length, length_scale, spinless, 0) + fermion_hamiltonian = jellium_model(grid, spinless, False) qubit_hamiltonian = jordan_wigner(fermion_hamiltonian) # Compute Jordan-Wigner Hamiltonian. - test_hamiltonian = jordan_wigner_position_jellium( - n_dimensions, grid_length, length_scale, spinless) + test_hamiltonian = jordan_wigner_position_jellium(grid, spinless) # Make sure Hamiltonians are the same. self.assertTrue(test_hamiltonian.isclose(qubit_hamiltonian)) diff --git a/src/fermilib/utils/_molecular_data.py b/src/fermilib/utils/_molecular_data.py index af236be..c336fd3 100644 --- a/src/fermilib/utils/_molecular_data.py +++ b/src/fermilib/utils/_molecular_data.py @@ -222,7 +222,7 @@ class MolecularData(object): ccsd_amplitudes: Molecular operator holding coupled cluster amplitude. """ def __init__(self, geometry=None, basis=None, multiplicity=None, - charge=0, description="", filename=""): + charge=0, description="", filename="", data_directory=None): """Initialize molecular metadata which defines class. Args: @@ -241,6 +241,8 @@ def __init__(self, geometry=None, basis=None, multiplicity=None, (e.g. 0.7414). filename: An optional string giving name of file. If filename is not provided, one is generated automatically. + data_directory: Optional data directory to change from default + data directory specified in config file. """ # Check appropriate data as been provided and autoload if requested. if ((geometry is None) or @@ -277,7 +279,10 @@ def __init__(self, geometry=None, basis=None, multiplicity=None, filename = filename[:(len(filename) - 5)] self.filename = filename else: - self.filename = DATA_DIRECTORY + '/' + self.name + if data_directory is None: + self.filename = DATA_DIRECTORY + '/' + self.name + else: + self.filename = data_directory + '/' + self.name # Attributes generated automatically by class. self.n_atoms = len(geometry) @@ -299,6 +304,7 @@ def __init__(self, geometry=None, basis=None, multiplicity=None, # Attributes generated from integrals. self.orbital_overlaps = None self.one_body_integrals = None + self.two_body_integrals = None # Attributes generated from MP2 calculation. self.mp2_energy = None @@ -306,17 +312,34 @@ def __init__(self, geometry=None, basis=None, multiplicity=None, # Attributes generated from CISD calculation. self.cisd_energy = None self.cisd_one_rdm = None + self.cisd_two_rdm = None # Attributes generated from exact diagonalization. self.fci_energy = None self.fci_one_rdm = None + self.fci_two_rdm = None # Attributes generated from CCSD calculation. self.ccsd_energy = None self.ccsd_amplitudes = None def save(self): - "Method to save the class under a systematic name." + """Method to save the class under a systematic name.""" + + # Load two body integrals/rdms from file before re-saving, since + # they aren't loaded by default + if (os.path.isfile("{}.hdf5".format(self.filename))): + if (self.one_body_integrals is not None and + self.two_body_integrals is None): + self.one_body_integrals, self.two_body_integrals = ( + self.get_integrals()) + if self.cisd_one_rdm is not None and self.cisd_two_rdm is None: + rdm = self.get_molecular_rdm() + self.cisd_two_rdm = rdm.two_body_tensor + if self.fci_one_rdm is not None and self.fci_two_rdm is None: + rdm = self.get_molecular_rdm(use_fci=True) + self.fci_two_rdm = rdm.two_body_tensor + with h5py.File("{}.hdf5".format(self.filename), "w") as f: # Save geometry (atoms and positions need to be separate): d_geom = f.create_group("geometry") @@ -333,7 +356,7 @@ def save(self): # Save description: f["description"] = numpy.string_(self.description) # Save name: - f["name"] = self.name + f["name"] = numpy.string_(self.name) # Save n_atoms: f["n_atoms"] = self.n_atoms # Save atoms: @@ -366,19 +389,26 @@ def save(self): f["one_body_integrals"] = (self.one_body_integrals if self.one_body_integrals is not None else False) + f["two_body_integrals"] = (self.two_body_integrals if + self.two_body_integrals is + not None else False) # Save attributes generated from MP2 calculation. f["mp2_energy"] = (self.mp2_energy if self.mp2_energy is not None else False) # Save attributes generated from CISD calculation. f["cisd_energy"] = (self.cisd_energy if self.cisd_energy is not None else False) - f["cisd_one_rmd"] = (self.cisd_one_rdm if + f["cisd_one_rdm"] = (self.cisd_one_rdm if self.cisd_one_rdm is not None else False) + f["cisd_two_rdm"] = (self.cisd_two_rdm if + self.cisd_two_rdm is not None else False) # Save attributes generated from exact diagonalization. f["fci_energy"] = (self.fci_energy if self.fci_energy is not None else False) f["fci_one_rdm"] = (self.fci_one_rdm if self.fci_one_rdm is not None else False) + f["fci_two_rdm"] = (self.fci_two_rdm if + self.fci_two_rdm is not None else False) # Save attributes generated from CCSD calculation. f["ccsd_energy"] = (self.ccsd_energy if self.ccsd_energy is not None else False) @@ -398,18 +428,18 @@ def load(self): # Load geometry: for atom, pos in zip(f["geometry/atoms"][...], f["geometry/positions"][...]): - geometry.append((str(atom), list(pos))) + geometry.append((atom.tobytes().decode('utf-8'), list(pos))) self.geometry = geometry # Load basis: - self.basis = str(f["basis"][...]) + self.basis = f["basis"][...].tobytes().decode('utf-8') # Load multiplicity: self.multiplicity = int(f["multiplicity"][...]) # Load charge: self.charge = int(f["charge"][...]) # Load description: - self.description = str(f["description"][...]) + self.description = f["description"][...].tobytes().decode('utf-8') # Load name: - self.name = str(f["name"][...]) + self.name = f["name"][...].tobytes().decode('utf-8') # Load n_atoms: self.n_atoms = int(f["n_atoms"][...]) # Load atoms: @@ -443,7 +473,7 @@ def load(self): # Load attributes generated from CISD calculation. d_9 = f["cisd_energy"][...] self.cisd_energy = d_9 if d_9.dtype.num != 0 else None - d_10 = f["cisd_one_rmd"][...] + d_10 = f["cisd_one_rdm"][...] self.cisd_one_rdm = d_10 if d_10.dtype.num != 0 else None # Load attributes generated from exact diagonalization. d_11 = f["fci_energy"][...] @@ -484,25 +514,30 @@ def get_integrals(self): # Make sure integrals have been computed. if self.hf_energy is None: raise MissingCalculationError( - 'Missing file {}. Run SCF before loading integrals.'.format( - self.filename + '_eri.npy')) + 'Missing SCF in {}, run before loading integrals.'.format( + self.filename)) # Get integrals and return. - two_body_integrals = numpy.load(self.filename + '_eri.npy') - return self.one_body_integrals, two_body_integrals + with h5py.File("{}.hdf5".format(self.filename), "r") as f: + tmp = f["two_body_integrals"][...] + self.two_body_integrals = tmp if tmp.dtype.num != 0 else None + return self.one_body_integrals, self.two_body_integrals - def get_active_space_integrals(self, active_space_start, - active_space_stop=None): - """Restricts a molecule at a spatial orbital level to the active space - defined by active_space=[start,stop]. Note that one_body_integrals and - two_body_integrals must be defined in an orthonormal basis set, which - is typically the case when defining an active space. + def get_active_space_integrals(self, + occupied_indices=None, + active_indices=None): + """Restricts a molecule at a spatial orbital level to an active space + + This active space may be defined by a list of active indices and + doubly occupied indices. Note that one_body_integrals and + two_body_integrals must be defined + n an orthonormal basis set. Args: - active_space_start(int): spatial orbital index defining active - space start. - active_space_stop(int): spatial orbital index defining active - space stop. + occupied_indices(list): A list of spatial orbital indices + indicating which orbitals should be considered doubly occupied. + active_indices(list): A list of spatial orbital indices indicating + which orbitals should be considered active. Returns: tuple: Tuple with the following entries: @@ -516,62 +551,63 @@ def get_active_space_integrals(self, active_space_start, **two_body_integrals_new**: two-electron integrals over active space. """ + # Fix data type for a few edge cases + occupied_indices = [] if occupied_indices is None else occupied_indices + if (len(active_indices) < 1): + raise ValueError('Some active indices required for reduction.') + # Get integrals. one_body_integrals, two_body_integrals = self.get_integrals() - n_orbitals = one_body_integrals.shape[0] - if active_space_stop is None: - active_space_stop = n_orbitals # Determine core constant core_constant = 0.0 - for i in range(active_space_start): + for i in occupied_indices: core_constant += 2 * one_body_integrals[i, i] - for j in range(active_space_start): + for j in occupied_indices: core_constant += (2 * two_body_integrals[i, j, j, i] - two_body_integrals[i, j, i, j]) # Modified one electron integrals one_body_integrals_new = numpy.copy(one_body_integrals) - for u in range(active_space_start, active_space_stop): - for v in range(active_space_start, active_space_stop): - for i in range(active_space_start): + for u in active_indices: + for v in active_indices: + for i in occupied_indices: one_body_integrals_new[u, v] += ( 2 * two_body_integrals[i, u, v, i] - two_body_integrals[i, u, i, v]) # Restrict integral ranges and change M appropriately return (core_constant, - one_body_integrals_new[active_space_start: active_space_stop, - active_space_start: active_space_stop], - two_body_integrals[active_space_start: active_space_stop, - active_space_start: active_space_stop, - active_space_start: active_space_stop, - active_space_start: active_space_stop]) + one_body_integrals_new[numpy.ix_(active_indices, + active_indices)], + two_body_integrals[numpy.ix_(active_indices, + active_indices, + active_indices, + active_indices)]) def get_molecular_hamiltonian(self, - active_space_start=None, - active_space_stop=None): + occupied_indices=None, + active_indices=None): """Output arrays of the second quantized Hamiltonian coefficients. Args: rotation_matrix: A square numpy array or matrix having dimensions of n_orbitals by n_orbitals. Assumed real and invertible. - active_space_start: An optional int giving the first orbital - in the active space. - active_space stop: An optional int giving the last orbital - in the active space. + occupied_indices(list): A list of spatial orbital indices + indicating which orbitals should be considered doubly occupied. + active_indices(list): A list of spatial orbital indices indicating + which orbitals should be considered active. Returns: molecular_hamiltonian: An instance of the MolecularOperator class. """ # Get active space integrals. - if active_space_start is None: + if occupied_indices is None and active_indices is None: one_body_integrals, two_body_integrals = self.get_integrals() constant = self.nuclear_repulsion else: core_adjustment, one_body_integrals, two_body_integrals = self.\ - get_active_space_integrals( - active_space_start, active_space_stop) + get_active_space_integrals(occupied_indices, active_indices) constant = self.nuclear_repulsion + core_adjustment n_qubits = 2 * one_body_integrals.shape[0] @@ -640,22 +676,23 @@ def get_molecular_rdm(self, use_fci=False): if use_fci: if self.fci_energy is None: raise MissingCalculationError( - 'Missing {}. '.format( - self.filename + '_fci_rdm.npy') + + 'Missing FCI RDM in {}'.format(self.filename) + 'Run FCI calculation before loading FCI RDMs.') else: - rdm_name = self.filename + '_fci_rdm.npy' + rdm_name = "fci_two_rdm" one_rdm = self.fci_one_rdm else: if self.cisd_energy is None: raise MissingCalculationError( - 'Missing {}. '.format( - self.filename + '_cisd_rdm.npy') + + 'Missing CISD RDM in {}'.format(self.filename) + 'Run CISD calculation before loading CISD RDMs.') else: - rdm_name = self.filename + '_cisd_rdm.npy' + rdm_name = "cisd_two_rdm" one_rdm = self.cisd_one_rdm - two_rdm = numpy.load(rdm_name) + + with h5py.File("{}.hdf5".format(self.filename), "r") as f: + tmp = f[rdm_name][...] + two_rdm = self.two_rdm = tmp if tmp.dtype.num != 0 else None # Truncate. one_rdm[numpy.absolute(one_rdm) < EQ_TOLERANCE] = 0. diff --git a/src/fermilib/utils/_molecular_data_test.py b/src/fermilib/utils/_molecular_data_test.py index 6226d7f..787cd31 100644 --- a/src/fermilib/utils/_molecular_data_test.py +++ b/src/fermilib/utils/_molecular_data_test.py @@ -28,19 +28,20 @@ def setUp(self): self.geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))] self.basis = 'sto-3g' self.multiplicity = 1 - filename = os.path.join(THIS_DIRECTORY, 'data', 'H2_sto-3g_singlet') + filename = os.path.join(THIS_DIRECTORY, 'data', + 'H2_sto-3g_singlet_0.7414') self.molecule = MolecularData( self.geometry, self.basis, self.multiplicity, filename=filename) self.molecule.load() def test_name_molecule(self): charge = 0 - correct_name = 'H2_sto-3g_singlet' + correct_name = str('H2_sto-3g_singlet_0.7414') computed_name = name_molecule(self.geometry, self.basis, self.multiplicity, charge, - description=None) + description="0.7414") self.assertEqual(correct_name, computed_name) self.assertEqual(correct_name, self.molecule.name) @@ -85,26 +86,49 @@ def test_dummy_save(self): molecule.orbital_energies = [5, 6, 7, 8] molecule.orbital_overlaps = [1, 2, 3, 4] molecule.one_body_integrals = [5, 6, 7, 8] + molecule.two_body_integrals = [5, 6, 7, 8] molecule.mp2_energy = -12. molecule.cisd_energy = 32. molecule.cisd_one_rdm = numpy.arange(10) + molecule.cisd_two_rdm = numpy.arange(10) molecule.fci_energy = 232. molecule.fci_one_rdm = numpy.arange(11) + molecule.fci_two_rdm = numpy.arange(11) molecule.ccsd_energy = 88. # Save molecule. molecule.save() - # Change attributes and load. - molecule.ccsd_energy = -2.232 + try: + # Change attributes and load. + molecule.ccsd_energy = -2.232 - # Load molecule. - new_molecule = MolecularData(filename=filename) - molecule.load() + # Load molecule. + new_molecule = MolecularData(filename=filename) + molecule.load() - # Check CCSD energy. - self.assertAlmostEqual(new_molecule.ccsd_energy, molecule.ccsd_energy) - self.assertAlmostEqual(molecule.ccsd_energy, 88.) + # Tests re-load functionality + molecule.save() + + # Check CCSD energy. + self.assertAlmostEqual(new_molecule.ccsd_energy, + molecule.ccsd_energy) + self.assertAlmostEqual(molecule.ccsd_energy, 88.) + finally: + os.remove(filename + '.hdf5') + + def test_active_space(self): + """Test simple active space truncation features""" + + # Start w/ no truncation + core_const, one_body_integrals, two_body_integrals = ( + self.molecule.get_active_space_integrals(active_indices=[0, 1])) + + self.assertAlmostEqual(core_const, 0.0) + self.assertAlmostEqual(scipy.linalg.norm(one_body_integrals - + self.molecule.one_body_integrals), 0.0) + self.assertAlmostEqual(scipy.linalg.norm(two_body_integrals - + self.molecule.two_body_integrals), 0.0) def test_energies(self): self.assertAlmostEqual(self.molecule.hf_energy, -1.1167, places=4) diff --git a/src/fermilib/utils/_plane_wave_hamiltonian.py b/src/fermilib/utils/_plane_wave_hamiltonian.py index 2dfb533..4da0da3 100644 --- a/src/fermilib/utils/_plane_wave_hamiltonian.py +++ b/src/fermilib/utils/_plane_wave_hamiltonian.py @@ -13,62 +13,55 @@ """Construct Hamiltonians in plan wave basis and its dual in 3D.""" from __future__ import absolute_import -import itertools import numpy from fermilib.config import * from fermilib.ops import FermionOperator +from fermilib.utils._grid import Grid from fermilib.utils._jellium import (orbital_id, grid_indices, position_vector, - momentum_vector, jellium_model) + momentum_vector, jellium_model, + jordan_wigner_position_jellium) +from fermilib.utils._molecular_data import periodic_hash_table from projectq.ops import QubitOperator -def dual_basis_u_operator(n_dimensions, grid_length, length_scale, - nuclear_charges, spinless): +def dual_basis_u_operator(grid, geometry, spinless): """Return the external potential operator in plane wave dual basis. Args: - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - nuclear_charges: 3D int array, the nuclear charges. - spinless: Bool, whether to use the spinless model or not. + grid (Grid): The discretization to use. + geometry: A list of tuples giving the coordinates of each atom. + example is [('H', (0, 0, 0)), ('H', (0, 0, 0.7414))]. + Distances in atomic units. Use atomic symbols to specify atoms. + spinless (bool): Whether to use the spinless model or not. Returns: - operator: An instance of the FermionOperator class. + FermionOperator: The dual basis operator. """ - n_points = grid_length ** n_dimensions - volume = length_scale ** float(n_dimensions) - prefactor = -4.0 * numpy.pi / volume + prefactor = -4.0 * numpy.pi / grid.volume_scale() operator = None if spinless: spins = [None] else: spins = [0, 1] - for grid_indices_p in itertools.product(range(grid_length), - repeat=n_dimensions): - coordinate_p = position_vector(grid_indices_p, grid_length, - length_scale) - for grid_indices_j in itertools.product(range(grid_length), - repeat=n_dimensions): - coordinate_j = position_vector(grid_indices_j, grid_length, - length_scale) - for momenta_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momenta = momentum_vector(momenta_indices, grid_length, - length_scale) - momenta_squred = momenta.dot(momenta) - if momenta_squred < EQ_TOLERANCE: + for pos_indices in grid.all_points_indices(): + coordinate_p = position_vector(pos_indices, grid) + for nuclear_term in geometry: + coordinate_j = numpy.array(nuclear_term[1], float) + for momenta_indices in grid.all_points_indices(): + momenta = momentum_vector(momenta_indices, grid) + momenta_squared = momenta.dot(momenta) + if momenta_squared < EQ_TOLERANCE: continue exp_index = 1.0j * momenta.dot(coordinate_j - coordinate_p) - coefficient = prefactor / momenta_squred * \ - nuclear_charges[grid_indices_j] * numpy.exp(exp_index) + coefficient = (prefactor / momenta_squared * + periodic_hash_table[nuclear_term[0]] * + numpy.exp(exp_index)) for spin_p in spins: - orbital_p = orbital_id( - grid_length, grid_indices_p, spin_p) + orbital_p = orbital_id(grid, pos_indices, spin_p) operators = ((orbital_p, 1), (orbital_p, 0)) if operator is None: operator = FermionOperator(operators, coefficient) @@ -78,56 +71,47 @@ def dual_basis_u_operator(n_dimensions, grid_length, length_scale, return operator -def plane_wave_u_operator(n_dimensions, grid_length, length_scale, - nuclear_charges, spinless): +def plane_wave_u_operator(grid, geometry, spinless): """Return the external potential operator in plane wave basis. Args: - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - nuclear_charges: 3D int array, the nuclear charges. + grid (Grid): The discretization to use. + geometry: A list of tuples giving the coordinates of each atom. + example is [('H', (0, 0, 0)), ('H', (0, 0, 0.7414))]. + Distances in atomic units. Use atomic symbols to specify atoms. spinless: Bool, whether to use the spinless model or not. Returns: - operator: An instance of the FermionOperator class. + FermionOperator: The plane wave operator. """ - n_points = grid_length ** n_dimensions - volume = length_scale ** float(n_dimensions) - prefactor = -4.0 * numpy.pi / volume + prefactor = -4.0 * numpy.pi / grid.volume_scale() operator = None if spinless: spins = [None] else: spins = [0, 1] - for grid_indices_p in itertools.product(range(grid_length), - repeat=n_dimensions): - for grid_indices_q in itertools.product(range(grid_length), - repeat=n_dimensions): - shift = grid_length // 2 + for indices_p in grid.all_points_indices(): + for indices_q in grid.all_points_indices(): + shift = grid.length // 2 grid_indices_p_q = [ - (grid_indices_p[i] - grid_indices_q[i] + shift) % grid_length - for i in range(n_dimensions)] - momenta_p_q = momentum_vector(grid_indices_p_q, grid_length, - length_scale) + (indices_p[i] - indices_q[i] + shift) % grid.length + for i in range(grid.dimensions)] + momenta_p_q = momentum_vector(grid_indices_p_q, grid) momenta_p_q_squared = momenta_p_q.dot(momenta_p_q) if momenta_p_q_squared < EQ_TOLERANCE: continue - for grid_indices_j in itertools.product(range(grid_length), - repeat=n_dimensions): - coordinate_j = position_vector(grid_indices_j, grid_length, - length_scale) + for nuclear_term in geometry: + coordinate_j = numpy.array(nuclear_term[1]) exp_index = 1.0j * momenta_p_q.dot(coordinate_j) - coefficient = prefactor / momenta_p_q_squared * \ - nuclear_charges[grid_indices_j] * numpy.exp(exp_index) + coefficient = (prefactor / momenta_p_q_squared * + periodic_hash_table[nuclear_term[0]] * + numpy.exp(exp_index)) for spin in spins: - orbital_p = orbital_id( - grid_length, grid_indices_p, spin) - orbital_q = orbital_id( - grid_length, grid_indices_q, spin) + orbital_p = orbital_id(grid, indices_p, spin) + orbital_q = orbital_id(grid, indices_q, spin) operators = ((orbital_p, 1), (orbital_q, 0)) if operator is None: operator = FermionOperator(operators, coefficient) @@ -137,41 +121,43 @@ def plane_wave_u_operator(n_dimensions, grid_length, length_scale, return operator -def plane_wave_hamiltonian(n_dimensions, grid_length, length_scale, - nuclear_charges, spinless=False, - momentum_space=True): +def plane_wave_hamiltonian(grid, geometry=None, + spinless=False, momentum_space=True): """Returns Hamiltonian as FermionOperator class. Args: - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - nuclear_charges: 3D int array, the nuclear charges. - spinless: Bool, whether to use the spinless model or not. - momentum_space: Boole, whether to return in plane wave basis (True) + grid (Grid): The discretization to use. + geometry: A list of tuples giving the coordinates of each atom. + example is [('H', (0, 0, 0)), ('H', (0, 0, 0.7414))]. + Distances in atomic units. Use atomic symbols to specify atoms. + spinless (bool): Whether to use the spinless model or not. + momentum_space (bool): Whether to return in plane wave basis (True) or plane wave dual basis (False). Returns: - hamiltonian: An instance of the FermionOperator class. + FermionOperator: The hamiltonian. """ - if len(nuclear_charges.shape) != n_dimensions: - raise ValueError('Invalid nuclear charges array shape.') + jellium_op = jellium_model(grid, spinless, momentum_space) + + if geometry is None: + return jellium_op + + for item in geometry: + if len(item[1]) != grid.dimensions: + raise ValueError("Invalid geometry coordinate.") + if item[0] not in periodic_hash_table: + raise ValueError("Invalid nuclear element.") if momentum_space: - return jellium_model(n_dimensions, grid_length, length_scale, spinless, - True) + \ - plane_wave_u_operator(n_dimensions, grid_length, length_scale, - nuclear_charges, spinless) + external_potential = plane_wave_u_operator(grid, geometry, spinless) else: - return jellium_model(n_dimensions, grid_length, length_scale, spinless, - False) + \ - dual_basis_u_operator(n_dimensions, grid_length, length_scale, - nuclear_charges, spinless) + external_potential = dual_basis_u_operator(grid, geometry, spinless) + + return jellium_op + external_potential -def fourier_transform(hamiltonian, n_dimensions, grid_length, length_scale, - spinless): - """Apply Fourier tranform to change hamiltonian in plane wave basis. +def fourier_transform(hamiltonian, grid, spinless): + """Apply Fourier transform to change hamiltonian in plane wave basis. .. math:: @@ -179,129 +165,137 @@ def fourier_transform(hamiltonian, n_dimensions, grid_length, length_scale, c_v = \sqrt{1/N} \sum_m {a_m \exp(i k_v r_m)} Args: - hamiltonian: The hamiltonian in plane wave basis. - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - spinless: Bool, whether to use the spinless model or not. + hamiltonian (FermionOperator): The hamiltonian in plane wave basis. + grid (Grid): The discretization to use. + spinless (bool): Whether to use the spinless model or not. + + Returns: + FermionOperator: The fourier-transformed hamiltonian. + """ + return _fourier_transform_helper(hamiltonian=hamiltonian, + grid=grid, + spinless=spinless, + phase_factor=+1, + vec_func_1=momentum_vector, + vec_func_2=position_vector) + + +def inverse_fourier_transform(hamiltonian, grid, spinless): + """Apply inverse Fourier transform to change hamiltonian in plane wave dual + basis. + + .. math:: + + a^\dagger_v = \sqrt{1/N} \sum_m {c^\dagger_m \exp(i k_v r_m)} + a_v = \sqrt{1/N} \sum_m {c_m \exp(-i k_v r_m)} + + Args: + hamiltonian (FermionOperator): + The hamiltonian in plane wave dual basis. + grid (Grid): The discretization to use. + spinless (bool): Whether to use the spinless model or not. Returns: - hamiltonian_t: An instance of the FermionOperator class. + FermionOperator: The inverse-fourier-transformed hamiltonian. """ - hamiltonian_t = None + + return _fourier_transform_helper(hamiltonian=hamiltonian, + grid=grid, + spinless=spinless, + phase_factor=-1, + vec_func_1=position_vector, + vec_func_2=momentum_vector) + + +def _fourier_transform_helper(hamiltonian, + grid, + spinless, + phase_factor, + vec_func_1, + vec_func_2): + hamiltonian_t = FermionOperator.zero() + normalize_factor = numpy.sqrt(1.0 / float(grid.num_points())) for term in hamiltonian.terms: - transformed_term = None - for ladder_operator in term: - momentum_indices = grid_indices(ladder_operator[0], n_dimensions, - grid_length, spinless) - momentum_vec = momentum_vector(momentum_indices, grid_length, - length_scale) - new_basis = None - for position_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - position_vec = position_vector(position_indices, grid_length, - length_scale) - if spinless: - spin = None - else: - spin = ladder_operator[0] % 2 - orbital = orbital_id(grid_length, position_indices, spin) - exp_index = 1.0j * numpy.dot(momentum_vec, position_vec) - if ladder_operator[1] == 1: + transformed_term = FermionOperator.identity() + for ladder_op_mode, ladder_op_type in term: + indices_1 = grid_indices(ladder_op_mode, grid, spinless) + vec1 = vec_func_1(indices_1, grid) + new_basis = FermionOperator.zero() + for indices_2 in grid.all_points_indices(): + vec2 = vec_func_2(indices_2, grid) + spin = None if spinless else ladder_op_mode % 2 + orbital = orbital_id(grid, indices_2, spin) + exp_index = phase_factor * 1.0j * numpy.dot(vec1, vec2) + if ladder_op_type == 1: exp_index *= -1.0 - element = FermionOperator(((orbital, ladder_operator[1]),), + element = FermionOperator(((orbital, ladder_op_type),), numpy.exp(exp_index)) - if new_basis is None: - new_basis = element - else: - new_basis += element - - new_basis *= numpy.sqrt(1.0/float(grid_length**n_dimensions)) - - if transformed_term is None: - transformed_term = new_basis - else: - transformed_term *= new_basis - if transformed_term is None: - continue + new_basis += element + + new_basis *= normalize_factor + transformed_term *= new_basis # Coefficient. transformed_term *= hamiltonian.terms[term] - if hamiltonian_t is None: - hamiltonian_t = transformed_term - else: - hamiltonian_t += transformed_term + hamiltonian_t += transformed_term return hamiltonian_t -def inverse_fourier_transform(hamiltonian, n_dimensions, grid_length, - length_scale, spinless): - """Apply Fourier tranform to change hamiltonian in plane wave dual basis. - - .. math:: - - a^\dagger_v = \sqrt{1/N} \sum_m {c^\dagger_m \exp(i k_v r_m)} - a_v = \sqrt{1/N} \sum_m {c_m \exp(-i k_v r_m)} +def jordan_wigner_dual_basis_hamiltonian(grid, geometry=None, spinless=False): + """Return the dual basis Hamiltonian as QubitOperator. Args: - hamiltonian: The hamiltonian in plane wave dual basis. - n_dimensions: An int giving the number of dimensions for the model. - grid_length: Int, the number of points in one dimension of the grid. - length_scale: Float, the real space length of a box dimension. - spinless: Bool, whether to use the spinless model or not. + grid (Grid): The discretization to use. + geometry: A list of tuples giving the coordinates of each atom. + example is [('H', (0, 0, 0)), ('H', (0, 0, 0.7414))]. + Distances in atomic units. Use atomic symbols to specify atoms. + spinless (bool): Whether to use the spinless model or not. Returns: - hamiltonian_t: An instance of the FermionOperator class. + hamiltonian (QubitOperator) """ - hamiltonian_t = None + jellium_op = jordan_wigner_position_jellium(grid, spinless) - for term in hamiltonian.terms: - transformed_term = None - for ladder_operator in term: - position_indices = grid_indices(ladder_operator[0], n_dimensions, - grid_length, spinless) - position_vec = position_vector(position_indices, grid_length, - length_scale) - new_basis = None - for momentum_indices in itertools.product(range(grid_length), - repeat=n_dimensions): - momentum_vec = momentum_vector(momentum_indices, grid_length, - length_scale) - if spinless: - spin = None - else: - spin = ladder_operator[0] % 2 - orbital = orbital_id(grid_length, momentum_indices, spin) - exp_index = -1.0j * numpy.dot(position_vec, momentum_vec) - if ladder_operator[1] == 1: - exp_index *= -1.0 + if geometry is None: + return jellium_op - element = FermionOperator(((orbital, ladder_operator[1]),), - numpy.exp(exp_index)) - if new_basis is None: - new_basis = element - else: - new_basis += element - - new_basis *= numpy.sqrt(1.0/float(grid_length**n_dimensions)) - - if transformed_term is None: - transformed_term = new_basis - else: - transformed_term *= new_basis - if transformed_term is None: + for item in geometry: + if len(item[1]) != grid.dimensions: + raise ValueError("Invalid geometry coordinate.") + if item[0] not in periodic_hash_table: + raise ValueError("Invalid nuclear element.") + + n_orbitals = grid.num_points() + volume = grid.volume_scale() + if spinless: + n_qubits = n_orbitals + else: + n_qubits = 2 * n_orbitals + prefactor = -2 * numpy.pi / volume + external_potential = QubitOperator() + + for k_indices in grid.all_points_indices(): + momenta = momentum_vector(k_indices, grid) + momenta_squared = momenta.dot(momenta) + if momenta_squared < EQ_TOLERANCE: continue - # Coefficient. - transformed_term *= hamiltonian.terms[term] + for p in range(n_qubits): + index_p = grid_indices(p, grid, spinless) + coordinate_p = position_vector(index_p, grid) - if hamiltonian_t is None: - hamiltonian_t = transformed_term - else: - hamiltonian_t += transformed_term + for nuclear_term in geometry: + coordinate_j = numpy.array(nuclear_term[1], float) - return hamiltonian_t + exp_index = 1.0j * momenta.dot(coordinate_j - coordinate_p) + coefficient = (prefactor / momenta_squared * + periodic_hash_table[nuclear_term[0]] * + numpy.exp(exp_index)) + external_potential += (QubitOperator((), coefficient) - + QubitOperator(((p, 'Z'),), coefficient)) + + return jellium_op + external_potential diff --git a/src/fermilib/utils/_plane_wave_hamiltonian_test.py b/src/fermilib/utils/_plane_wave_hamiltonian_test.py index e502dc9..4c1bc35 100644 --- a/src/fermilib/utils/_plane_wave_hamiltonian_test.py +++ b/src/fermilib/utils/_plane_wave_hamiltonian_test.py @@ -13,109 +13,96 @@ """Tests for plane_wave_hamiltonian.py""" from __future__ import absolute_import -import numpy import unittest -from fermilib.ops import * +import numpy + +from fermilib.ops import normal_ordered from fermilib.transforms import jordan_wigner -from fermilib.utils import eigenspectrum -from fermilib.utils._plane_wave_hamiltonian import * +from fermilib.utils import eigenspectrum, Grid +from fermilib.utils._plane_wave_hamiltonian import ( + plane_wave_hamiltonian, + fourier_transform, + inverse_fourier_transform, + plane_wave_u_operator, + dual_basis_u_operator, + jordan_wigner_dual_basis_hamiltonian, +) class PlaneWaveHamiltonianTest(unittest.TestCase): def test_fourier_transform(self): - n_dimensions = 1 - length_scale = 1.5 - grid_length = 3 + grid = Grid(dimensions=1, scale=1.5, length=3) spinless_set = [True, False] - nuclear_charges = numpy.empty((3)) - nuclear_charges[0] = 1 - nuclear_charges[1] = -3 - nuclear_charges[2] = 2 + geometry = [('H', (0,)), ('H', (0.5,))] for spinless in spinless_set: h_plane_wave = plane_wave_hamiltonian( - n_dimensions, grid_length, length_scale, nuclear_charges, - spinless, True) + grid, geometry, spinless, True) h_dual_basis = plane_wave_hamiltonian( - n_dimensions, grid_length, length_scale, nuclear_charges, - spinless, False) - h_plane_wave_t = fourier_transform( - h_plane_wave, n_dimensions, grid_length, length_scale, - spinless) + grid, geometry, spinless, False) + h_plane_wave_t = fourier_transform(h_plane_wave, grid, spinless) self.assertTrue(normal_ordered(h_plane_wave_t).isclose( normal_ordered(h_dual_basis))) def test_inverse_fourier_transform_1d(self): - n_dimensions = 1 - length_scale = 1.5 - grid_length = 3 + grid = Grid(dimensions=1, scale=1.5, length=4) spinless_set = [True, False] - nuclear_charges = numpy.empty((3)) - nuclear_charges[0] = 1 - nuclear_charges[1] = -3 - nuclear_charges[2] = 2 + geometry = [('H', (0,)), ('H', (0.5,))] for spinless in spinless_set: h_plane_wave = plane_wave_hamiltonian( - n_dimensions, grid_length, length_scale, nuclear_charges, - spinless, True) + grid, geometry, spinless, True) h_dual_basis = plane_wave_hamiltonian( - n_dimensions, grid_length, length_scale, nuclear_charges, - spinless, False) + grid, geometry, spinless, False) h_dual_basis_t = inverse_fourier_transform( - h_dual_basis, n_dimensions, grid_length, length_scale, - spinless) + h_dual_basis, grid, spinless) self.assertTrue(normal_ordered(h_dual_basis_t).isclose( normal_ordered(h_plane_wave))) def test_inverse_fourier_transform_2d(self): - n_dimensions = 2 - length_scale = 1.5 - grid_length = 3 + grid = Grid(dimensions=2, scale=1.5, length=3) spinless = True - nuclear_charges = numpy.empty((3, 3)) - nuclear_charges[0][1] = 1 - nuclear_charges[0][2] = 1 - nuclear_charges[1][2] = 2 - nuclear_charges[1][0] = -3 - nuclear_charges[2][2] = 2 - nuclear_charges[2][1] = -3 - h_plane_wave = plane_wave_hamiltonian( - n_dimensions, grid_length, length_scale, nuclear_charges, - spinless, True) - h_dual_basis = plane_wave_hamiltonian( - n_dimensions, grid_length, length_scale, nuclear_charges, - spinless, False) + geometry = [('H', (0, 0)), ('H', (0.5, 0.8))] + h_plane_wave = plane_wave_hamiltonian(grid, geometry, spinless, True) + h_dual_basis = plane_wave_hamiltonian(grid, geometry, spinless, False) h_dual_basis_t = inverse_fourier_transform( - h_dual_basis, n_dimensions, grid_length, length_scale, - spinless) + h_dual_basis, grid, spinless) self.assertTrue(normal_ordered(h_dual_basis_t).isclose( normal_ordered(h_plane_wave))) - def test_u_operator_integration(self): - n_dimensions = 1 - length_scale = 1 - grid_length = 3 + def test_plane_wave_hamiltonian_integration(self): + length_set = [3, 4] spinless_set = [True, False] - nuclear_charges = numpy.empty((3)) - nuclear_charges[0] = 1 - nuclear_charges[1] = -3 - nuclear_charges[2] = 2 - for spinless in spinless_set: - u_plane_wave = plane_wave_u_operator( - n_dimensions, grid_length, length_scale, nuclear_charges, - spinless) - u_dual_basis = dual_basis_u_operator( - n_dimensions, grid_length, length_scale, nuclear_charges, - spinless) - jw_u_plane_wave = jordan_wigner(u_plane_wave) - jw_u_dual_basis = jordan_wigner(u_dual_basis) - u_plane_wave_spectrum = eigenspectrum(jw_u_plane_wave) - u_dual_basis_spectrum = eigenspectrum(jw_u_dual_basis) - - diff = numpy.amax(numpy.absolute( - u_plane_wave_spectrum - u_dual_basis_spectrum)) - self.assertAlmostEqual(diff, 0) + geometry = [('H', (0,)), ('H', (0.8,))] + for l in length_set: + for spinless in spinless_set: + grid = Grid(dimensions=1, scale=1.0, length=l) + h_plane_wave = plane_wave_hamiltonian(grid, geometry, spinless, + True) + h_dual_basis = plane_wave_hamiltonian(grid, geometry, spinless, + False) + jw_h_plane_wave = jordan_wigner(h_plane_wave) + jw_h_dual_basis = jordan_wigner(h_dual_basis) + h_plane_wave_spectrum = eigenspectrum(jw_h_plane_wave) + h_dual_basis_spectrum = eigenspectrum(jw_h_dual_basis) + + diff = numpy.amax(numpy.absolute( + h_plane_wave_spectrum - h_dual_basis_spectrum)) + self.assertAlmostEqual(diff, 0) + + def test_jordan_wigner_dual_basis_hamiltonian(self): + grid = Grid(dimensions=2, length=3, scale=1.) + spinless = True + geometry = [('H', (0, 0)), ('H', (0.5, 0.8))] + + fermion_hamiltonian = plane_wave_hamiltonian(grid, geometry, spinless, + False) + qubit_hamiltonian = jordan_wigner(fermion_hamiltonian) + + test_hamiltonian = jordan_wigner_dual_basis_hamiltonian(grid, geometry, + spinless) + self.assertTrue(test_hamiltonian.isclose(qubit_hamiltonian)) + # Run test. if __name__ == '__main__': diff --git a/src/fermilib/utils/_sparse_tools.py b/src/fermilib/utils/_sparse_tools.py index 164f632..98fdb84 100644 --- a/src/fermilib/utils/_sparse_tools.py +++ b/src/fermilib/utils/_sparse_tools.py @@ -246,7 +246,7 @@ def sparse_eigenspectrum(sparse_operator): eigenspectrum = numpy.linalg.eigvalsh(dense_operator) else: eigenspectrum = numpy.linalg.eigvals(dense_operator) - return eigenspectrum + return numpy.sort(eigenspectrum) def expectation(sparse_operator, state): diff --git a/src/fermilib/utils/_trotter_error.py b/src/fermilib/utils/_trotter_error.py index c689a43..bab6ee3 100644 --- a/src/fermilib/utils/_trotter_error.py +++ b/src/fermilib/utils/_trotter_error.py @@ -32,20 +32,18 @@ def trivially_commutes(term_a, term_b): position_b = 0 commutes = True - len_a = len(list(term_a.terms.keys())[0]) - len_b = len(list(term_b.terms.keys())[0]) + term_op_a, = term_a.terms.keys() + term_op_b, = term_b.terms.keys() - while position_a < len_a and position_b < len_b: - qubit_a = list(term_a.terms.keys())[0][position_a][0] - qubit_b = list(term_b.terms.keys())[0][position_b][0] + while position_a < len(term_op_a) and position_b < len(term_op_b): + qubit_a, action_a = term_op_a[position_a] + qubit_b, action_b = term_op_b[position_b] if qubit_a > qubit_b: position_b += 1 elif qubit_a < qubit_b: position_a += 1 else: - action_a = list(term_a.terms.keys())[0][position_a][1] - action_b = list(term_b.terms.keys())[0][position_b][1] if action_a != action_b: commutes = not commutes position_a += 1 @@ -66,12 +64,13 @@ def trivially_double_commutes(term_a, term_b, term_c): qubits) is empty, then the double commutator is trivially zero. """ # determine the set of qubits each term acts on - qubits_a = set([list(term_a.terms.keys())[0][i][0] - for i in range(len(list(term_a.terms.keys())[0]))]) - qubits_b = set([list(term_b.terms.keys())[0][i][0] - for i in range(len(list(term_b.terms.keys())[0]))]) - qubits_c = set([list(term_c.terms.keys())[0][i][0] - for i in range(len(list(term_c.terms.keys())[0]))]) + term_op_a, = term_a.terms.keys() + term_op_b, = term_b.terms.keys() + term_op_c, = term_c.terms.keys() + + qubits_a = set([index for index, _ in term_op_a]) + qubits_b = set([index for index, _ in term_op_b]) + qubits_c = set([index for index, _ in term_op_c]) return (trivially_commutes(term_b, term_c) or not qubits_a.intersection(set(qubits_b.union(qubits_c)))) @@ -153,13 +152,13 @@ def error_bound(terms, tight=False): elif not tight: for alpha in range(len(terms)): term_a = terms[alpha] - coefficient_a = list(term_a.terms.values())[0] + coefficient_a, = term_a.terms.values() if coefficient_a: error_a = 0. for beta in range(alpha + 1, len(terms)): term_b = terms[beta] - coefficient_b = list(term_b.terms.values())[0] + coefficient_b, = term_b.terms.values() if not (trivially_commutes(term_a, term_b) or commutator(term_a, term_b).isclose(zero)): error_a += abs(coefficient_b) diff --git a/src/fermilib/utils/_unitary_cc.py b/src/fermilib/utils/_unitary_cc.py index 7cf8000..d843b39 100644 --- a/src/fermilib/utils/_unitary_cc.py +++ b/src/fermilib/utils/_unitary_cc.py @@ -11,7 +11,6 @@ # limitations under the License. """Module to create and manipulate unitary coupled cluster operators.""" -from __future__ import absolute_import import itertools import numpy @@ -26,24 +25,30 @@ import projectq.ops import projectq.setups import projectq.setups.decompositions +import projectq.types def uccsd_operator(single_amplitudes, double_amplitudes, anti_hermitian=True): """Create a fermionic operator that is the generator of uccsd. This a the most straight-forward method to generate UCCSD operators, - however it is slightly inefficient. In particular, it parameterizes - all possible exictations, so it represents a generalized unitary coupled + however it is slightly inefficient. In particular, it parameterizes + all possible excitations, so it represents a generalized unitary coupled cluster ansatz, but also does not explicitly enforce the uniqueness - in paramterization, so it is redundant. For example there will be a linear + in parametrization, so it is redundant. For example there will be a linear dependency in the ansatz of single_amplitudes[i,j] and single_amplitudes[j,i]. Args: - single_amplitudes(ndarray): [NxN] array storing single excitation - amplitudes corresponding to t[i,j] * (a_i^\dagger a_j + H.C.) - double_amplitudes(ndarray): [NxNxNxN] array storing double excitation + single_amplitudes(list or ndarray): list of lists with each sublist + storing a list of indices followed by single excitation amplitudes + i.e. [[[i,j],t_ij], ...] OR [NxN] array storing single excitation amplitudes corresponding to + t[i,j] * (a_i^\dagger a_j + H.C.) + double_amplitudes(list or ndarray): list of lists with each sublist + storing a list of indices followed by double excitation amplitudes + i.e. [[[i,j,k,l],t_ijkl], ...] OR [NxNxNxN] array storing double + excitation amplitudes corresponding to t[i,j,k,l] * (a_i^\dagger a_j a_k^\dagger a_l + H.C.) anti_hermitian(Bool): Flag to generate only normal CCSD operator rather than unitary variant, primarily for testing @@ -52,33 +57,61 @@ def uccsd_operator(single_amplitudes, double_amplitudes, anti_hermitian=True): uccsd_generator(FermionOperator): Anti-hermitian fermion operator that is the generator for the uccsd wavefunction. """ - n_orbitals = single_amplitudes.shape[0] - assert(n_orbitals == double_amplitudes.shape[0]) + uccsd_generator = FermionOperator() + # Re-format inputs (ndarrays to lists) if necessary + if (isinstance(single_amplitudes, numpy.ndarray) or + isinstance(double_amplitudes, numpy.ndarray)): + single_amplitudes, double_amplitudes = convert_amplitude_format( + single_amplitudes, + double_amplitudes) + # Add single excitations - for i, j in itertools.product(range(n_orbitals), repeat=2): - if single_amplitudes[i, j] == 0.: - continue + for (i, j), t_ij in single_amplitudes: + i, j = int(i), int(j) + uccsd_generator += FermionOperator(((i, 1), (j, 0)), t_ij) + if anti_hermitian: + uccsd_generator += FermionOperator(((j, 1), (i, 0)), -t_ij) + + # Add double excitations + for (i, j, k, l), t_ijkl in double_amplitudes: + i, j, k, l = int(i), int(j), int(k), int(l) uccsd_generator += FermionOperator( - ((i, 1), (j, 0)), single_amplitudes[i, j]) - if (anti_hermitian): + ((i, 1), (j, 0), (k, 1), (l, 0)), t_ijkl) + if anti_hermitian: uccsd_generator += FermionOperator( - ((j, 1), (i, 0)), -single_amplitudes[i, j]) + ((l, 1), (k, 0), (j, 1), (i, 0)), -t_ijkl) + return uccsd_generator - # Add double excitations - for i, j, k, l in itertools.product(range(n_orbitals), repeat=4): - if double_amplitudes[i, j, k, l] == 0.: - continue - uccsd_generator += FermionOperator( - ((i, 1), (j, 0), (k, 1), (l, 0)), - double_amplitudes[i, j, k, l]) - if (anti_hermitian): - uccsd_generator += FermionOperator( - ((l, 1), (k, 0), (j, 1), (i, 0)), - -double_amplitudes[i, j, k, l]) - return uccsd_generator +def convert_amplitude_format(single_amplitudes, double_amplitudes): + """Re-format single_amplitudes and double_amplitudes from ndarrays to lists. + + Args: + single_amplitudes(ndarray): [NxN] array storing single excitation + amplitudes corresponding to t[i,j] * (a_i^\dagger a_j + H.C.) + double_amplitudes(ndarray): [NxNxNxN] array storing double excitation + amplitudes corresponding to + t[i,j,k,l] * (a_i^\dagger a_j a_k^\dagger a_l + H.C.) + + Returns: + single_amplitudes_list(list): list of lists with each sublist storing + a list of indices followed by single excitation amplitudes + i.e. [[[i,j],t_ij], ...] + double_amplitudes_list(list): list of lists with each sublist storing + a list of indices followed by double excitation amplitudes + i.e. [[[i,j,k,l],t_ijkl], ...] + """ + single_amplitudes_list, double_amplitudes_list = [], [] + + for i, j in zip(*single_amplitudes.nonzero()): + single_amplitudes_list.append([[i, j], single_amplitudes[i, j]]) + + for i, j, k, l in zip(*double_amplitudes.nonzero()): + double_amplitudes_list.append([[i, j, k, l], + double_amplitudes[i, j, k, l]]) + return single_amplitudes_list, double_amplitudes_list def uccsd_singlet_paramsize(n_qubits, n_electrons): @@ -126,7 +159,6 @@ def uccsd_singlet_operator(packed_amplitudes, n_occupied = int(numpy.ceil(n_electrons / 2.)) n_virtual = int(n_qubits / 2 - n_occupied) # Virtual Spatial Orbitals n_t1 = int(n_occupied * n_virtual) - n_t2 = int(n_t1 ** 2) t1 = packed_amplitudes[:n_t1] t2 = packed_amplitudes[n_t1:] @@ -135,48 +167,76 @@ def t1_ind(i, j): return i * n_occupied + j def t2_ind(i, j, k, l): - return i * n_occupied * n_virtual * n_occupied \ - + j * n_virtual * n_occupied \ - + k * n_occupied \ - + l + return (i * n_occupied * n_virtual * n_occupied + + j * n_virtual * n_occupied + + k * n_occupied + + l) uccsd_generator = FermionOperator() - for i in range(n_virtual): - for j in range(n_occupied): - for s1 in range(2): - uccsd_generator += FermionOperator(( - (2 * (i + n_occupied) + s1, 1), - (2 * j + s1, 0)), - t1[t1_ind(i, j)]) - - uccsd_generator += FermionOperator(( - (2 * j + s1, 1), - (2 * (i + n_occupied) + s1, 0)), - -t1[t1_ind(i, j)]) - - for i in range(n_virtual): - for j in range(n_occupied): - for s1 in range(2): - for k in range(n_virtual): - for l in range(n_occupied): - for s2 in range(2): - uccsd_generator += FermionOperator(( - (2 * (i + n_occupied) + s1, 1), - (2 * j + s1, 0), - (2 * (k + n_occupied) + s2, 1), - (2 * l + s2, 0)), - t2[t2_ind(i, j, k, l)]) - - uccsd_generator += FermionOperator(( - (2 * l + s2, 1), - (2 * (k + n_occupied) + s2, 0), - (2 * j + s1, 1), - (2 * (i + n_occupied) + s1, 0)), - -t2[t2_ind(i, j, k, l)]) + spaces = range(n_virtual), range(n_occupied), range(2) + + for i, j, s in itertools.product(*spaces): + uccsd_generator += FermionOperator( + ( + (2 * (i + n_occupied) + s, 1), + (2 * j + s, 0), + ), + coefficient=t1[t1_ind(i, j)]) + + uccsd_generator += FermionOperator( + ( + (2 * j + s, 1), + (2 * (i + n_occupied) + s, 0), + ), + coefficient=-t1[t1_ind(i, j)]) + + for i, j, s, i2, j2, s2 in itertools.product(*spaces, repeat=2): + uccsd_generator += FermionOperator(( + (2 * (i + n_occupied) + s, 1), + (2 * j + s, 0), + (2 * (i2 + n_occupied) + s2, 1), + (2 * j2 + s2, 0)), + t2[t2_ind(i, j, i2, j2)]) + + uccsd_generator += FermionOperator(( + (2 * j2 + s2, 1), + (2 * (i2 + n_occupied) + s2, 0), + (2 * j + s, 1), + (2 * (i + n_occupied) + s, 0)), + -t2[t2_ind(i, j, i2, j2)]) + return uccsd_generator +def uccsd_evolution(fermion_generator, fermion_transform=jordan_wigner): + """Create a ProjectQ evolution operator for a UCCSD circuit + + Args: + fermion_generator(FermionOperator): UCCSD generator to evolve. + fermion_transform(fermilib.transform): The transformation that + defines the mapping from Fermions to QubitOperator. + + Returns: + evoution_operator(projectq.ops.TimeEvolution): The unitary operator + that constructs the UCCSD state. + """ + + # Transform generator to qubits + qubit_generator = fermion_transform(fermion_generator) + + # Cast to real part only for compatibility with current ProjectQ routine + for key in qubit_generator.terms: + qubit_generator.terms[key] = float(qubit_generator.terms[key].imag) + qubit_generator.compress() + + # Allocate wavefunction and act evolution on gate according to compilation + evolution_operator = ( + projectq.ops.TimeEvolution(time=1., hamiltonian=qubit_generator)) + + return evolution_operator + + def uccsd_singlet_evolution(packed_amplitudes, n_qubits, n_electrons, fermion_transform=jordan_wigner): """Create a ProjectQ evolution operator for a UCCSD singlet circuit @@ -200,18 +260,9 @@ def uccsd_singlet_evolution(packed_amplitudes, n_qubits, n_electrons, fermion_generator = uccsd_singlet_operator(packed_amplitudes, n_qubits, n_electrons) - # Transform generator to qubits - qubit_generator = fermion_transform(fermion_generator) - # Cast to real part only for compatibility with current ProjectQ routine - for key in qubit_generator.terms: - qubit_generator.terms[key] = float(qubit_generator.terms[key].imag) - qubit_generator.compress() - - # Allocate wavefunction and act evolution on gate according to compilation - evolution_operator = projectq.ops.\ - TimeEvolution(time=1., - hamiltonian=qubit_generator) + evolution_operator = uccsd_evolution(fermion_generator, + fermion_transform) return evolution_operator @@ -237,8 +288,9 @@ def _identify_non_commuting(cmd): for term in hamiltonian.terms: test_op = projectq.ops.QubitOperator(term, hamiltonian.terms[term]) for other in hamiltonian.terms: - other_op = projectq.\ - ops.QubitOperator(other, hamiltonian.terms[other]) + other_op = ( + projectq.ops.QubitOperator(other, + hamiltonian.terms[other])) commutator = test_op * other_op - other_op * test_op if not commutator.isclose(id_op, rel_tol=1e-9, @@ -247,6 +299,100 @@ def _identify_non_commuting(cmd): return False +def _non_adjacent_filter(self, cmd, qubit_graph, flip=False): + """A ProjectQ filter to identify when swaps are needed on a graph + + This flags any gates that act on two non-adjacent qubits with respect to + the qubit_graph that has been given + + Args: + self(Dummy): Dummy parameter to meet function specification. + cmd(projectq.command): Command to be checked for decomposition into + additional swap gates. + qubit_graph(Graph): Graph object specifying connectivity of + qubits. The values of the nodes of this graph are unique qubit ids. + flip(Bool): Flip for switching if identifying a gate is in this class + by true or false. Designed to meet the specification of ProjectQ + InstructionFilter and DecompositionRule with one function. + + Returns: + bool: When flip is False, this returns True when a 2 qubit command + acts on non-adjacent qubits or when it acts only on a single qubit. + This is reversed when flip is used. + + """ + if qubit_graph is None: + return True ^ flip + + total_qubits = (cmd.control_qubits + + [item for qureg in cmd.qubits for item in qureg]) + + # Check for non-connected gate on 2 qubits + if ((len(total_qubits) == 1) or + (len(total_qubits) == 2 and + qubit_graph.is_adjacent( + qubit_graph.find_index(total_qubits[0].id), + qubit_graph.find_index(total_qubits[1].id)))): + return True ^ flip + return False ^ flip + + +def _direct_graph_swap(cmd, qubit_graph): + """Define a naive direct swap sequence to respect qubit_graph connectivity + + Uses the connectivity of qubit_graph to find the shortest path between + two non-adjacent qubits, and swaps/unswaps qubits appropriately. Baseline + for more sophisticated algorithms + + Args: + cmd(projectq.command): A command from ProjectQ that needs to be + broken down due to non-adjacent terms + qubit_graph(Graph): Graph object specifying connectivity of qubits. + The values of the nodes of this graph are unique qubit ids + """ + total_qubits = (cmd.control_qubits + + [item for qureg in cmd.qubits for item in qureg]) + + gate = cmd.gate + engine = cmd.engine + graph_path = qubit_graph.shortest_path( + qubit_graph.find_index(total_qubits[0].id), + qubit_graph.find_index(total_qubits[1].id)) + swap_path = [(graph_path[i], graph_path[i + 1]) + for i in range(len(graph_path) - 2)] + + # SWAP qubit 1 into position adjacent to qubit 2 + for pair in swap_path: + projectq.ops.Swap | (projectq.types. + WeakQubitRef(engine, + qubit_graph.nodes[pair[0]].value), + projectq.types. + WeakQubitRef(engine, + qubit_graph.nodes[pair[1]].value)) + + # Perform original gate + if len(cmd.control_qubits) > 0: + projectq.ops.C(gate) | (projectq.types. + WeakQubitRef(engine, + qubit_graph. + nodes[graph_path[-2]].value), + total_qubits[1]) + else: + gate | (projectq.types. + WeakQubitRef(engine, + qubit_graph.nodes[graph_path[-2]].value), + total_qubits[1]) + + # Reverse the swaps to put qubits back in place + for pair in reversed(swap_path): + projectq.ops.Swap | (projectq.types. + WeakQubitRef(engine, + qubit_graph.nodes[pair[0]].value), + projectq.types. + WeakQubitRef(engine, + qubit_graph.nodes[pair[1]].value)) + + def _first_order_trotter(cmd): """Define a Trotter splitting for non-commuting Pauli in ProjectQ @@ -255,7 +401,7 @@ def _first_order_trotter(cmd): Args: cmd(projectq.command): A command from ProjectQ that needs to be - factorized due to non-commuting terms + factorized due to non-commuting time evolution terms """ qureg = cmd.qubits eng = cmd.engine @@ -264,8 +410,9 @@ def _first_order_trotter(cmd): with projectq.meta.Control(eng, cmd.control_qubits): # First order Trotter splitting for term in hamiltonian.terms: - ind_operator = projectq.\ - ops.QubitOperator(term, hamiltonian.terms[term]) + ind_operator = (projectq. + ops. + QubitOperator(term, hamiltonian.terms[term])) projectq.ops.TimeEvolution(time, ind_operator) | qureg @@ -289,7 +436,8 @@ def _two_gate_filter(self, cmd): return False -def uccsd_trotter_engine(compiler_backend=projectq.backends.Simulator()): +def uccsd_trotter_engine(compiler_backend=projectq.backends.Simulator(), + qubit_graph=None): """Define a ProjectQ compiler engine that is common for use with UCCSD This defines a ProjectQ compiler engine that decomposes time evolution @@ -301,26 +449,49 @@ def uccsd_trotter_engine(compiler_backend=projectq.backends.Simulator()): circuit compiler, so that it may either simulate gates numerically or alternatively print a gate sequence, e.g. using projectq.backends.CommandPrinter() + qubit_graph(Graph): Graph object specifying connectivity of qubits. + The values of the nodes of this unique qubit ids. If None, + all-to-all connectivity is assumed. Returns: projectq.cengine that is the compiler engine set up with these rules and decompostions. """ - rule_set = \ - projectq.cengines. \ - DecompositionRuleSet(modules=[projectq.setups.decompositions]) - trotter_rule_set = projectq.cengines. \ - DecompositionRule(gate_class=projectq.ops.TimeEvolution, - gate_decomposer=_first_order_trotter, - gate_recognizer=_identify_non_commuting) + rule_set = ( + projectq.cengines. + DecompositionRuleSet(modules=[projectq.setups.decompositions])) + + # Set rules for splitting non-commuting operators + trotter_rule_set = (projectq.cengines.DecompositionRule( + gate_class=projectq.ops.TimeEvolution, + gate_decomposer=_first_order_trotter, + gate_recognizer=_identify_non_commuting)) rule_set.add_decomposition_rule(trotter_rule_set) + + # Set rules for 2 qubit gates that act on non-adjacent qubits + if qubit_graph is not None: + connectivity_rule_set = ( + projectq.cengines.DecompositionRule( + gate_class=projectq.ops.NOT.__class__, + gate_decomposer=(lambda x: _direct_graph_swap(x, qubit_graph)), + gate_recognizer=(lambda x: _non_adjacent_filter(None, x, + qubit_graph, + True)))) + rule_set.add_decomposition_rule(connectivity_rule_set) + + # Build the full set of engines that will be applied to qubits replacer = projectq.cengines.AutoReplacer(rule_set) + compiler_engine_list = [replacer, + projectq. + cengines. + InstructionFilter( + lambda x, y: + (_non_adjacent_filter(x, y, qubit_graph) and + _two_gate_filter(x, y))), + projectq.cengines.LocalOptimizer(5)] # Start the compiler engine with these rules - compiler_engine = projectq.\ - MainEngine(backend=compiler_backend, - engine_list=[replacer, - projectq. - cengines.InstructionFilter(_two_gate_filter), - projectq.cengines.LocalOptimizer(3)]) + compiler_engine = ( + projectq.MainEngine(backend=compiler_backend, + engine_list=compiler_engine_list)) return compiler_engine diff --git a/src/fermilib/utils/_unitary_cc_test.py b/src/fermilib/utils/_unitary_cc_test.py index bb19e91..6ab3408 100644 --- a/src/fermilib/utils/_unitary_cc_test.py +++ b/src/fermilib/utils/_unitary_cc_test.py @@ -16,14 +16,15 @@ import unittest from numpy.random import randn +import itertools import fermilib import fermilib.ops from fermilib.ops import FermionOperator import fermilib.utils +from fermilib.utils._graph import (Graph, Node) from fermilib.utils._unitary_cc import * -from projectq import MainEngine from projectq.ops import (All, Measure, TimeEvolution, QubitOperator, X) @@ -85,31 +86,6 @@ def test_uccsd_singlet_build(self): (-0.0565340614) * FermionOperator("0^ 2 1^ 3")) self.assertTrue(test_generator.isclose(generator)) - # Skip this test for now. - @unittest.skip - def test_projectq_filters(self): - """Verify ProjectQ filters work as intended""" - eng = MainEngine() - op = QubitOperator("X0 X1 X2", -0.5) - wavefunction = eng.allocate_qureg(4) - - command = TimeEvolution(time=1., hamiltonian=op) | wavefunction - self.assertFalse(_identify_non_commuting(command)) - - op = QubitOperator("X0 X1 X2", -0.5) + QubitOperator("Z1", 1.0) - command = TimeEvolution(time=1., hamiltonian=op) | wavefunction - - self.assertTrue(_identify_non_commuting(command)) - - op = QubitOperator("X0 X1 X2", -0.5) - command = TimeEvolution(time=1., hamiltonian=op) | wavefunction - - self.assertFalse(_two_gate_filter(None, command)) - - command = X | wavefunction - - self.assertTrue(_two_gate_filter(None, command)) - def test_simulation_energy(self): """Test UCCSD Singlet Energy for H2""" @@ -146,3 +122,99 @@ def test_simulation_energy(self): wavefunction) All(Measure) | wavefunction self.assertAlmostEqual(energy, -1.13727017463) + + def test_simulation_with_graph(self): + """Test UCCSD Singlet Energy for H2 using a restricted qubit_graph""" + + # Define H2 Hamiltonian inline + hamiltonian = ((-0.0453222020986) * QubitOperator("X0 X1 Y2 Y3") + + (0.165867023964) * QubitOperator("Z0 Z3") + + (0.174348441706) * QubitOperator("Z2 Z3") + + (0.120544821866) * QubitOperator("Z0 Z2") + + (3.46944695195e-18) * QubitOperator("X0 Y1 X2 Y3") + + (0.165867023964) * QubitOperator("Z1 Z2") + + (0.171197748533) * QubitOperator("Z0") + + (-0.222785928901) * QubitOperator("Z3") + + (3.46944695195e-18) * QubitOperator("X0 X1 X2 X3") + + (0.168622191433) * QubitOperator("Z0 Z1") + + (0.120544821866) * QubitOperator("Z1 Z3") + + (3.46944695195e-18) * QubitOperator("Y0 Y1 Y2 Y3") + + (-0.0988639735178) * QubitOperator("") + + (0.171197748533) * QubitOperator("Z1") + + (0.0453222020986) * QubitOperator("Y0 X1 X2 Y3") + + (3.46944695195e-18) * QubitOperator("Y0 X1 Y2 X3") + + (-0.0453222020986) * QubitOperator("Y0 Y1 X2 X3") + + (-0.222785928901) * QubitOperator("Z2") + + (0.0453222020986) * QubitOperator("X0 Y1 Y2 X3")) + hamiltonian.compress() + + # Create a star graph of 4 qubits, all connected through qubit 0 + qubit_graph = Graph() + compiler_engine = uccsd_trotter_engine(qubit_graph=qubit_graph) + wavefunction = compiler_engine.allocate_qureg(4) + for i in range(4): + qubit_graph.add_node(Node(value=wavefunction[i].id)) + for i in range(1, 4): + qubit_graph.add_edge(0, i) + + test_amplitudes = [-1.14941450e-08, 5.65340614e-02] + for i in range(2): + X | wavefunction[i] + evolution_operator = uccsd_singlet_evolution(test_amplitudes, 4, 2) + evolution_operator | wavefunction + compiler_engine.flush() + + energy = compiler_engine.backend.get_expectation_value(hamiltonian, + wavefunction) + All(Measure) | wavefunction + self.assertAlmostEqual(energy, -1.13727017463) + + def test_sparse_uccsd_operator_numpy_inputs(self): + """Test numpy ndarray inputs to uccsd_operator that are sparse""" + test_orbitals = 30 + sparse_single_amplitudes = numpy.zeros((test_orbitals, test_orbitals)) + sparse_double_amplitudes = numpy.zeros((test_orbitals, test_orbitals, + test_orbitals, test_orbitals)) + + sparse_single_amplitudes[3, 5] = 0.12345 + sparse_single_amplitudes[12, 4] = 0.44313 + + sparse_double_amplitudes[0, 12, 6, 2] = 0.3434 + sparse_double_amplitudes[1, 4, 6, 13] = -0.23423 + + generator = uccsd_operator(sparse_single_amplitudes, + sparse_double_amplitudes) + + test_generator = (0.12345 * FermionOperator("3^ 5") + + (-0.12345) * FermionOperator("5^ 3") + + 0.44313 * FermionOperator("12^ 4") + + (-0.44313) * FermionOperator("4^ 12") + + 0.3434 * FermionOperator("0^ 12 6^ 2") + + (-0.3434) * FermionOperator("2^ 6 12^ 0") + + (-0.23423) * FermionOperator("1^ 4 6^ 13") + + 0.23423 * FermionOperator("13^ 6 4^ 1")) + self.assertTrue(test_generator.isclose(generator)) + + def test_sparse_uccsd_operator_list_inputs(self): + """Test list inputs to uccsd_operator that are sparse""" + sparse_single_amplitudes = [[[3, 5], 0.12345], + [[12, 4], 0.44313]] + sparse_double_amplitudes = [[[0, 12, 6, 2], 0.3434], + [[1, 4, 6, 13], -0.23423]] + + generator = uccsd_operator(sparse_single_amplitudes, + sparse_double_amplitudes) + + test_generator = (0.12345 * FermionOperator("3^ 5") + + (-0.12345) * FermionOperator("5^ 3") + + 0.44313 * FermionOperator("12^ 4") + + (-0.44313) * FermionOperator("4^ 12") + + 0.3434 * FermionOperator("0^ 12 6^ 2") + + (-0.3434) * FermionOperator("2^ 6 12^ 0") + + (-0.23423) * FermionOperator("1^ 4 6^ 13") + + 0.23423 * FermionOperator("13^ 6 4^ 1")) + self.assertTrue(test_generator.isclose(generator)) + +# Run test. +if __name__ == '__main__': + unittest.main()