Updated README. (#8)

README.md

@@ -114,3 +114,19 @@ opCacheVec = [
 
 return MPO_new(os, sites; basisOpCacheVec=opCacheVec)
 ```
+
+## Examples: Electronic Structure Hamiltonian
+
+After looking at the previous example one might assume that the relative speed of ITensorMPOConstruction over Renormalizer might be due to the fact that for the Fermi-Hubbard Hamiltonian ITensorMPOConstruction is able to construct a much more compact MPO. In the case of the electronic structure Hamiltonian all algorithms produce MPOs of similar bond dimensions.
+
+![](./docs/plot-generators/es.png)
+
+However ITensorMPOConstruction is still an order of magnitude faster than Renormalizer. The code for this example can be found in [examples/electronic-structure.jl](https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/electronic-structure.jl). The run time to generate these MPOs are shown below (again on a 2021 MacBook Pro with the M1 Max CPU and 32GB of memory).
+
+| $N$ | ITensors | Renormalizer | ITensorMPOConstruction |
+|-----|----------|--------------|------------------------|
+| 10  | 3.0s     | 2.1s         | 0.37s                  |
+| 20  | 498s     | 61s          | 7.1s                   |
+| 30  | N/A      | 458s         | 46s                    |
+| 40  | N/A      | 2250s        | 183                    |
+| 50  | N/A      | N/A          | 510s                   |

docs/plot-generators/renormalizer-plots.py

@@ -1,9 +1,10 @@
 import numpy as np
 from timeit import default_timer as timer
 from renormalizer import Model, Op, BasisSimpleElectron, Mpo
+import random
 
 
-def fermiHubbardKS(N, t, U):
+def fermi_hubbard_ks(N, t, U):
     toSite = lambda k, s: 2 * k + s
 
     terms = []
@@ -22,12 +23,58 @@ def fermiHubbardKS(N, t, U):
     model = Model([BasisSimpleElectron(i) for i in range(2 * N)], terms)
     return Mpo(model)
 
+
+def electronic_structure(N):
+    toSite = lambda k, s: 2 * k + s
+    coeff = lambda : random.random()
+
+    terms = []
+    for a in range(N):
+        for b in range(a, N):
+            factor = coeff()
+            for spin in (0, 1):
+                sites = [toSite(a, spin), toSite(b, spin)]
+                terms.append(Op(r"a^\dagger a", sites, factor))
+
+                if a != b:
+                    sites = [toSite(b, spin), toSite(a, spin)]
+                    terms.append(Op(r"a^\dagger a", sites, np.conj(factor)))
+
+    for j in range(N):
+        for s_j in (0, 1):
+
+            for k in range(N):
+                s_k = s_j
+                if (s_k, k) <= (s_j, j):
+                    continue
+
+                for l in range(N):
+                    for s_l in (0, 1):
+                        if (s_l, l) <= (s_j, j):
+                            continue
+
+                        for m in range(N):
+                            s_m = s_l
+                            if (s_m, m) <= (s_k, k):
+                                continue
+
+                            value = coeff()
+                            sites = [toSite(j, s_j), toSite(l, s_l), toSite(m, s_m), toSite(k, s_k)]
+                            terms.append(Op(r"a^\dagger a^\dagger a a", sites, factor))
+
+                            sites = list(reversed(sites))
+                            terms.append(Op(r"a^\dagger a^\dagger a a", sites, np.conj(factor)))
+
+    model = Model([BasisSimpleElectron(i) for i in range(2 * N)], terms)
+    return Mpo(model)
+
+
 numSites = []
 bondDims = []
 times = []
-for N in [10, 20, 30, 40, 50]:
+for N in [40]:
   start = timer()
-  mpo = fermiHubbardKS(N, 1, 4)
+  mpo = electronic_structure(N)
   stop = timer()
 
   numSites.append(N)