Fix a bug in openQCD I/O.

Author: SaltyChiang

pyquda_io/_field_utils.py

@@ -78,10 +78,9 @@ def gaugePlaquette(latt_size: List[int], grid_size: List[int], gauge: numpy.ndar
     plaq[3] = numpy.vdot(gauge[0] @ extended[3, :-1, :-1, :-1, 1:], gauge[3] @ extended[0, 1:, :-1, :-1, :-1]).real
     plaq[4] = numpy.vdot(gauge[1] @ extended[3, :-1, :-1, 1:, :-1], gauge[3] @ extended[1, 1:, :-1, :-1, :-1]).real
     plaq[5] = numpy.vdot(gauge[2] @ extended[3, :-1, 1:, :-1, :-1], gauge[3] @ extended[2, 1:, :-1, :-1, :-1]).real
-
     plaq /= int(numpy.prod(latt_size)) * Nc
     plaq = MPI.COMM_WORLD.allreduce(plaq, MPI.SUM)
-    return [plaq.mean().item(), plaq[:3].mean().item(), plaq[3:].mean().item()]
+    return numpy.array([plaq.mean(), plaq[:3].mean(), plaq[3:].mean()])
 
 
 def gaugeOddShiftForward(latt_size: List[int], grid_size: List[int], gauge: numpy.ndarray):
@@ -229,8 +228,8 @@ def gaugeReconstruct12(gauge: numpy.ndarray):
     gauge_ = gauge.transpose(5, 6, 0, 1, 2, 3, 4)
     gauge = numpy.empty((Nc, *gauge_.shape[1:]), "<c16")
     gauge[:2] = gauge_
-    gauge[2] = numpy.cross(gauge[..., 0, :], gauge[..., 1, :], axis=-1).conjugate()
-    return gauge
+    gauge[2] = numpy.cross(gauge[0], gauge[1], axis=0).conjugate()
+    return gauge.transpose(2, 3, 4, 5, 6, 0, 1)
 
 
 # matrices to convert gamma basis bewteen DeGrand-Rossi and Dirac-Pauli

pyquda_io/openqcd.py

@@ -15,7 +15,7 @@ def readGauge(filename: str, grid_size: List[int], plaquette: bool = True, lexic
     filename = path.expanduser(path.expandvars(filename))
     with open(filename, "rb") as f:
         latt_size = struct.unpack("<iiii", f.read(16))[::-1]
-        plaquette = struct.unpack("<d", f.read(8))[0] / Nc
+        plaquette_ = struct.unpack("<d", f.read(8))[0] / Nc
         offset = f.tell()
     Lx, Ly, Lz, Lt = getSublatticeSize(latt_size, grid_size)
     dtype = "<c16"
@@ -35,10 +35,9 @@ def readGauge(filename: str, grid_size: List[int], plaquette: bool = True, lexic
     if lexico:
         gauge = gaugeLexico([Lx, Ly, Lz, Lt], gauge)
         if plaquette:
-            assert numpy.isclose(gaugePlaquette(latt_size, grid_size, gauge)[0], plaquette)
+            assert numpy.isclose(gaugePlaquette(latt_size, grid_size, gauge)[0], plaquette_)
     elif plaquette:
-        gauge_lexico = gaugeLexico([Lx, Ly, Lz, Lt], gauge)
-        assert numpy.isclose(gaugePlaquette(latt_size, grid_size, gauge_lexico)[0], plaquette)
+        assert numpy.isclose(gaugePlaquette(latt_size, grid_size, gaugeLexico(latt_size, gauge))[0], plaquette_)
     gauge = gauge.astype("<c16")
 
     return latt_size, gauge
@@ -62,8 +61,7 @@ def writeGauge(
             plaquette = gaugePlaquette(latt_size, grid_size, gauge)[0]
         gauge = gaugeEvenOdd([Lx, Ly, Lz, Lt], gauge)
     elif plaquette is None:
-        gauge_lexico = gaugeLexico([Lx, Ly, Lz, Lt], gauge)
-        plaquette = gaugePlaquette(latt_size, grid_size, gauge_lexico)[0]
+        plaquette = gaugePlaquette(latt_size, grid_size, gaugeLexico([Lx, Ly, Lz, Lt], gauge))[0]
     gauge = gaugeEvenShiftBackward(latt_size, grid_size, gauge)
     gauge_reorder = numpy.zeros((Lt, Lx, Ly, Lz // 2, Nd, 2, Nc, Nc), dtype)
     for t in range(Lt):