fix: Brillouin zone meshes (#30)

* fixed Brillouin zone meshes

* fix reshape for kmesh, and some variable typos

* make G vectors along columns

* fix G to follow docs

---------

Co-authored-by: H. L. Nourse <[email protected]>

docs/how-to/quadratic_terms.md

@@ -28,6 +28,7 @@ First construct a Monkhorst-Pack $k$-point mesh in fractional coordinates
 
 ```python
 import numpy as np
+from itertools import product
 
 # The spatial dimensions of the problem
 d = ...
@@ -49,14 +50,24 @@ shift_1 = ...
 shift_d = ...
 shift_list = [shift_1, ..., shift_d]
 
+# Method 1: Using np.meshgrid
+
 # Create linear array in each spatial dimension
-k_linear_list = [np.linspace(0, 1, nk_list[i], endpoint = False) + shift_list[i] / nk_list[i] for i in range(d)]
+k_linear_list = [np.linspace(0, 1, n_k_list[i], endpoint = False) + shift_list[i] / n_k_list[i] for i in range(d)]
 
 # Create a meshgrid
 k_mesh = np.meshgrid(*k_linear_list, indexing='ij')
 
 # Make it a list indexed as idx, k1, k2, ...
-k_mesh = np.reshape(k_mesh, (n_k, d))
+k_mesh = np.reshape(k_mesh, (d,-1)).T
+
+# Method 2: Using an explicit for loop
+
+# Create empty mesh and populate it
+k_mesh = np.empty(shape=(n_k, d))
+for idx, coords in enumerate(product( *[range(n_k_list[i]) for i in range(len(n_k_list))] )):
+    for dim in range(d):
+        k_mesh[idx,dim] = (coords[dim] + shift_list[dim]) / n_k_list[dim]
 ```
 
 If your function for $\hat{H}_0$ is in a different basis then you will 
@@ -74,11 +85,16 @@ a_d = [float_1, ..., float_d]
 A = np.array([a_1, a_2, ...]).T
 
 # The reciprocol lattice vectors as a matrix where each column is a vector
-G = 2 * np.pi * np.linalg.inv(A)
+G = 2 * np.pi * np.linalg.inv(A).T
 
 # Rotate mesh into this basis as k_cartesian = G @ k_fractional
-# For numpy broadcasting: k_cartesian = (G @ k_fractional.T).T
+
+# Method 1: Using numpy broadcasting
 k_mesh = (G @ k_mesh.T).T
+
+# Method 2: Using an explicit for loop
+for k in range(k_mesh.shape[0]):
+    k_mesh[k,:] = G @ k_mesh[k,:]
 ```
 
 Next, construct `h0_k` as

examples/decorated_honeycomb.py

@@ -20,7 +20,7 @@ def get_h0_k(tg=0.5, tk=1.0, nkx=18, spin_names=['up','dn'], basis='trimer'):
 
     # Build shifted 2D mesh
     mesh = np.empty(shape=(n_k, 2))
-    for idx,coords in enumerate(zip(range(nkx), range(nkx))):
+    for idx,coords in enumerate(product(range(nkx), range(nkx))):
         mesh[idx,0] = coords[0]/nkx + 0.5/nkx
         mesh[idx,1] = coords[1]/nkx + 0.5/nkx
 
@@ -30,7 +30,7 @@ def get_h0_k(tg=0.5, tk=1.0, nkx=18, spin_names=['up','dn'], basis='trimer'):
     R = np.array((R1, R2)).T
 
     # Bravais lattice vectors
-    G = 2.0*np.pi*np.linalg.inv(R)
+    G = 2.0*np.pi*np.linalg.inv(R).T
 
     # Vectors to inter-triangle nearest neighbors
     d0 = ( 1.0, 0.0 )
@@ -42,7 +42,7 @@ def get_h0_k(tg=0.5, tk=1.0, nkx=18, spin_names=['up','dn'], basis='trimer'):
 
     # Construct in inequivalent block matrix structure  
     for k,i,j,m,mm in product(range(n_k), range(na), range(na), range(n_orb), range(n_orb)):
-        kay = np.dot(G.T, mesh[k,:])
+        kay = np.dot(G, mesh[k,:])
 
         # Dispersion terms between clusters
         if (i == 0) and (j == 1):

examples/kagome_hubbard.py

@@ -16,7 +16,7 @@ def get_h0_k(t=1, nkx=18, spin_names=['up','dn']):
     # Build shifted 2D mesh
     n_k = nkx**2
     mesh = np.empty(shape=(n_k, 2))
-    for idx,coords in enumerate(zip(range(nkx), range(nkx))):
+    for idx,coords in enumerate(product(range(nkx), range(nkx))):
         mesh[idx,0] = coords[0]/nkx + 0.5/nkx
         mesh[idx,1] = coords[1]/nkx + 0.5/nkx
 
@@ -26,11 +26,11 @@ def get_h0_k(t=1, nkx=18, spin_names=['up','dn']):
     R = np.array((R1, R2)).T
 
     # Bravais lattice vectors
-    G = 2.0*np.pi*np.linalg.inv(R)
+    G = 2.0*np.pi*np.linalg.inv(R).T
 
     h0_k = np.zeros([n_k, n_orb, n_orb], dtype=complex)
     for k in range(n_k):
-        kay = np.dot(G.T, mesh[k,:])
+        kay = np.dot(G, mesh[k,:])
         k1 = kay[0]
         k2 = 0.5 * kay[0] + 0.5 * np.sqrt(3) * kay[1]
         k3 = -0.5 * kay[0] + 0.5 * np.sqrt(3) * kay[1]