Add periodic boundary condition resource estimates and factorizations.

[fdmalone]

Adds functionality to:

• Compute sparse, SF, DF, and THC factorizations using symmetry adapted Bloch orbitals under periodic boundary conditions.
• Compute resource estimates for these factorizations.

This code was used to generate the majority of the data in the recent paper: Fault-tolerant quantum simulation of materials using Bloch orbitals

The code follows the structure of the molecular resource estimation code (see #763), and provides utility functions to factorize the integrals if appropriate, compute the lambda value of the hamiltonian and compute fault tolerant resource estimates. A lot of boilerplate can be reduced by using the generate_costing_table functions which perform the necessary steps.

A minimal example is given below. In depth tutorials are provided in notebooks/isdf.ipynb and notebooks/resource_estimates.ipynb.

The first step is to run a periodic Hartree-Fock calculation for a system of interest (here carbon in the diamond structure)

from ase.build import bulk
import numpy as np

from pyscf.pbc import gto, scf
from pyscf.pbc.tools import pyscf_ase

# Build a 2 atom unit cell for carbon in the diamond structure near it's
# equilibrium lattice constant.
ase_atom = bulk("C", "diamond", a=3.5)
cell = gto.Cell()
cell.atom = pyscf_ase.ase_atoms_to_pyscf(ase_atom)
cell.a = ase_atom.cell[:].copy()
# Using a minimal basis set for expediency.
cell.basis = "gth-szv"
cell.pseudo = "gth-hf-rev"
cell.verbose = 0
cell.build()

# We are using a very small k-point mesh for speed purposes too.
kmesh = [1, 1, 3]
kpts = cell.make_kpts(kmesh)
num_kpts = len(kpts)
mf = scf.KRHF(cell, kpts).rs_density_fit() # RSGDF is required!
mf.kernel()
print("SCF energy: ", mf.e_tot)

From this mean field object we can generate a table of resource estimates given a set of threshold parameters. Here we choose the sparse representation so the thresholds correspond to when matrix elements of the Hamiltonian are set to zero.

thresholds = np.logspace(-1, -5, 5)
sparse_costing_table = sparse.generate_costing_table(mf, thresholds=thresholds)
print(sparse_costing_table.to_string(index=False))

The output of which is

system_name  num_spin_orbitals  num_kpts     dE  chi  exact_energy energy_method  cutoff  approx_energy  lambda_total  lambda_one_body  lambda_two_body  num_sym_unique  toffolis_per_step  total_toffolis  logical_qubits
        pbc                 16         3 0.0016   10     -0.126802           MP2 0.10000      -0.062425    556.289220        23.739303       532.549917            1406               1088       594195968             385
        pbc                 16         3 0.0016   10     -0.126802           MP2 0.01000      -0.125852    940.971717        23.739303       917.232414            9129               2214      2045286558            1145
        pbc                 16         3 0.0016   10     -0.126802           MP2 0.00100      -0.126799   1000.895662        23.739303       977.156359           16901               2752      2704192256            1143
        pbc                 16         3 0.0016   10     -0.126802           MP2 0.00010      -0.126802   1003.045646        23.739303       979.306344           19318               2927      2882328126            1146
        pbc                 16         3 0.0016   10     -0.126802           MP2 0.00001      -0.126802   1003.080141        23.739303       979.340838           19760               2940      2895229680            1146

Note the return type of the costing table is a pandas DataFrame. I decided it was more convenient this way rather than writing to a text file (this can be delegated to pandas to_csv(), to_string(), to_latex() methods), but I'm open to suggestions.

More fine grained control can be achieved through the compute_lambda and compute_cost functions.

The code also implements a k-point dependent THC factorization based upon an interpolative separable density fitting (ISDF) initial guess. This procedure is quite slow and involved but can be accelerated a bit on GPUs. I did not try terribly hard to optimize einsum contractions so there are potentially some optimizations to be had.

Much like the molecular resource estimates, this code is not hooked up to the CI and skips the tests if pyscf is not available. We should probably use markers instead and/or enable the tests on a less frequent / incremental basis if possible.

Major contributions were made by Nick Rubin and Alec White.

[fdmalone]

Looks like numpy deprecated np.float. I'll open a separate PR to fix this so as to not pollute the git history with changes to mainline openfermion.

[mpharrigan]

note that np.float was always sortof an accident and was an alias for the built-in float type. If you wanted the numpy/C type it is spelled np.float_

[fdmalone]

I'm going to split this PR into several smaller ones to make it a bit more digestible.