Merge pull request #28 from ITensor/corbett/one-large-graph-left

Corbett/one large graph left

Author: corbett5

Project.toml

@@ -4,18 +4,19 @@ authors = ["Ben Corbett and contributors"]
 version = "0.1.0"
 
 [deps]
+FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 ITensorMPS = "0.1, 0.2, 0.3"
 ITensors = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9"
 LinearAlgebra = "1.9, 1.10"
-Memoize = "0.4"
 SparseArrays = "1.9, 1.10"
 TimerOutputs = "0.5"
 julia = "1.10"

README.md

@@ -1,38 +1,124 @@
 import numpy as np
 from timeit import default_timer as timer
 from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes
+import random
 
-L = 100
-N = L
-TWOSZ = 0
+def fermi_hubbard(N, alg, t=1, U=4, J=0):
+    start = timer()
+    
+    driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
+    driver.initialize_system(n_sites=N, n_elec=N, spin=0)
 
-driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=4)
-driver.initialize_system(n_sites=L, n_elec=N, spin=TWOSZ)
+    b = driver.expr_builder()
 
-t = 1
-U = 4
+    epsilon = lambda k : 2 * np.cos(2 * np.pi * k / N)
 
-b = driver.expr_builder()
+    # Kinetic term
+    for k in range(N):
+        b.add_term("cd", np.array([k, k]), -t * epsilon(k))
+        b.add_term("CD", np.array([k, k]), -t * epsilon(k))
 
-epsilon = lambda k : 2 * np.cos(2 * np.pi * k / L)
+    # Interaction term
+    for p in range(N):
+        for q in range(N):
+            for k in range(N):
+                if J == 0:
+                    b.add_term("cCDd", np.array([(p - k) % N, (q + k) % N, q, p]), U / N)
+                else:
+                    factor = U - t * ((np.exp(J) - 1) * epsilon(p - k) + (np.exp(-J) - 1) * epsilon(p))
+                    b.add_term("cCDd", np.array([(p - k) % N, (q + k) % N, q, p]), factor / N)
+                    b.add_term("CcdD", np.array([(p - k) % N, (q + k) % N, q, p]), factor / N)
 
-# Kinetic term
-for k in range(L):
-  b.add_term("cd", np.array([k, k]), -t * epsilon(k))
-  b.add_term("CD", np.array([k, k]), -t * epsilon(k))
+    if J != 0:
+        prefactor = 2 * t * (np.cosh(J) - 1) / N**2
+        for p in range(N):
+            for q in range(N):
+                for s in range(N):
+                    for k in range(N):
+                        for kp in range(N):
+                            factor = prefactor * epsilon(p - (k - kp))
+                            b.add_term("cCCDDd", np.array([(p - k) % N, (q + kp) % N, (s + k - kp) % N, s, q, p]), factor)
+                            b.add_term("CccddD", np.array([(p - k) % N, (q + kp) % N, (s + k - kp) % N, s, q, p]), factor)
 
-# Interaction term
-sites = np.zeros(4 * L * L * L, dtype=np.int16)
-for p in range(L):
-  for q in range(L):
-    for k in range(L):
-      b.add_term("cCDd", np.array([(p - k) % L, (q + k) % L, q, p]), U)
+    print("Done constructing representation.")
 
-print("Done constructing representation.")
+    b = b.finalize()
+    print("Done finalizing representation.")
 
-alg = MPOAlgorithmTypes.SVD | MPOAlgorithmTypes.Blocked
-# alg = MPOAlgorithmTypes.SVD | MPOAlgorithmTypes.Fast | MPOAlgorithmTypes.Blocked | MPOAlgorithmTypes.Disjoint
-start = timer()
-mpo = driver.get_mpo(b.finalize(), iprint=1, algo_type=alg)
-stop = timer()
-print(f"N = {L}, time = {stop - start}")
+    mpo = driver.get_mpo(b, iprint=1, algo_type=alg)
+    stop = timer()
+    print(f"N = {N}, time = {stop - start}")
+
+
+def get_coefficients(N):
+    rng = np.random.default_rng(0)
+    one_electron = rng.normal(size=(N, N))
+    two_electron = rng.normal(size=(N, 2, N, 2, N, 2, N, 2))
+
+    return one_electron, two_electron
+
+def electronic_structure(N, alg):
+    one_electron_coeffs, two_electron_coeffs = get_coefficients(N)    
+
+    start = timer()
+    driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ, n_threads=8)
+    driver.initialize_system(n_sites=N, n_elec=N, spin=0, orb_sym=None)
+
+    builder = driver.expr_builder()
+
+    c = lambda spin : "c" if (spin == 0) else "C"
+    d = lambda spin : "d" if (spin == 0) else "D"
+
+    for a in range(N):
+        for b in range(a, N):
+            factor = one_electron_coeffs[a, b]
+            for spin in (0, 1):
+                sites = np.array([a, b])
+                builder.add_term(c(spin) + d(spin), sites, factor)
+
+                if a != b:
+                    sites = [b, a]
+                    builder.add_term(c(spin) + d(spin), sites, 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
+
+                            factor = two_electron_coeffs[j, s_j, l, s_l, m, s_m, k, s_k]
+                            sites = np.array([j, l, m, k])
+                            builder.add_term(c(s_j) + c(s_l) + d(s_m) + d(s_k), sites, factor)
+                            
+                            sites = np.flip(sites)
+                            builder.add_term(c(s_k) + c(s_m) + d(s_l) + d(s_j), sites, factor)
+
+    print("Done constructing representation.")
+
+    builder = builder.finalize()
+    print("Done finalizing representation.")
+
+    mpo = driver.get_mpo(builder, iprint=1, algo_type=alg)
+    stop = timer()
+    print(f"N = {N}, time = {stop - start}")
+
+
+for N in [10]:
+    fermi_hubbard(N, MPOAlgorithmTypes.FastBlockedDisjointSVD)
+    print()
+
+for N in [10, 10, 20, 30, 40, 50, 60, 70]:
+    electronic_structure(N, MPOAlgorithmTypes.FastBipartite)
+    print()