From 634a8af1997ce4d3baaa1fb5fa21a16f44393b1a Mon Sep 17 00:00:00 2001
From: SaltyChiang <SaltyChiang@users.noreply.github.com>
Date: Thu, 12 Dec 2024 15:40:43 +0800
Subject: [PATCH] Fix a bug in GPT interface. Add a package called `pyquda_io`.

---
 .gitignore                                    |  1 +
 pyproject.toml                                |  2 +-
 pyquda_io/__init__.py                         | 46 +++++++++++++
 .../_field_utils.py                           | 68 ++++++++++++-------
 .../io/mpi_file.py => pyquda_io/_mpi_file.py  |  0
 {pyquda_utils/io => pyquda_io}/chroma.py      |  2 +-
 {pyquda_utils/io => pyquda_io}/io_general.py  |  0
 {pyquda_utils/io => pyquda_io}/kyu.py         |  5 +-
 {pyquda_utils/io => pyquda_io}/lime.py        |  0
 {pyquda_utils/io => pyquda_io}/milc.py        |  2 +-
 {pyquda_utils/io => pyquda_io}/nersc.py       |  4 +-
 {pyquda_utils/io => pyquda_io}/npy.py         |  2 +-
 {pyquda_utils/io => pyquda_io}/openqcd.py     |  4 +-
 {pyquda_utils/io => pyquda_io}/xqcd.py        |  5 +-
 pyquda_utils/gpt.py                           |  2 +-
 pyquda_utils/io/__init__.py                   | 64 +++++++----------
 setup.py                                      |  2 +
 17 files changed, 135 insertions(+), 74 deletions(-)
 create mode 100644 pyquda_io/__init__.py
 rename pyquda_utils/io/gauge_utils.py => pyquda_io/_field_utils.py (82%)
 rename pyquda_utils/io/mpi_file.py => pyquda_io/_mpi_file.py (100%)
 rename {pyquda_utils/io => pyquda_io}/chroma.py (98%)
 rename {pyquda_utils/io => pyquda_io}/io_general.py (100%)
 rename {pyquda_utils/io => pyquda_io}/kyu.py (88%)
 rename {pyquda_utils/io => pyquda_io}/lime.py (100%)
 rename {pyquda_utils/io => pyquda_io}/milc.py (98%)
 rename {pyquda_utils/io => pyquda_io}/nersc.py (97%)
 rename {pyquda_utils/io => pyquda_io}/npy.py (97%)
 rename {pyquda_utils/io => pyquda_io}/openqcd.py (94%)
 rename {pyquda_utils/io => pyquda_io}/xqcd.py (89%)

diff --git a/.gitignore b/.gitignore
index 825fc2f..c49ff08 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,6 +22,7 @@ pyquda_core/pyquda/enum_quda.py
 pyquda_core/pyquda/src/*.c
 pyquda_core/pyquda/src/*.cpp
 pyquda_utils/_version.py
+pyquda_io/_version.py
 
 # symbolic link
 /pyquda
diff --git a/pyproject.toml b/pyproject.toml
index 074daf0..f2e2296 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -2,7 +2,7 @@
 requires = ["setuptools", "wheel"]
 
 [tool.setuptools.packages.find]
-include = ["pyquda_utils*"]
+include = ["pyquda_utils*", "pyquda_io*"]
 
 [tool.setuptools.dynamic]
 version = { attr = "pyquda_utils._version.__version__" }
diff --git a/pyquda_io/__init__.py b/pyquda_io/__init__.py
new file mode 100644
index 0000000..b501efa
--- /dev/null
+++ b/pyquda_io/__init__.py
@@ -0,0 +1,46 @@
+# flake8: noqa
+
+from .io_general import (
+    IOGeneral,
+    read as readIOGeneral,
+    write as writeIOGeneral,
+)
+from .lime import (
+    Lime,
+)
+from .npy import (
+    readGauge as readNPYGauge,
+    writeGauge as writeNPYGauge,
+    readPropagator as readNPYPropagator,
+    writePropagator as writeNPYPropagator,
+)
+from .chroma import (
+    readQIOGauge as readChromaQIOGauge,
+    readQIOPropagator as readChromaQIOPropagator,
+    readILDGBinGauge,
+)
+from .milc import (
+    readGauge as readMILCGauge,
+    writeGauge as writeMILCGauge,
+    readQIOPropagator as readMILCQIOPropagator,
+)
+from .kyu import (
+    readGauge as readKYUGauge,
+    writeGauge as writeKYUGauge,
+    readPropagator as readKYUPropagator,
+    writePropagator as writeKYUPropagator,
+)
+from .xqcd import (
+    readPropagator as readXQCDPropagator,
+    writePropagator as writeXQCDPropagator,
+    readPropagatorFast as readXQCDPropagatorFast,
+    writePropagatorFast as writeXQCDPropagatorFast,
+)
+from .nersc import (
+    readGauge as readNERSCGauge,
+    writeGauge as writeNERSCGauge,
+)
+from .openqcd import (
+    readGauge as readOpenQCDGauge,
+    writeGauge as writeOpenQCDGauge,
+)
diff --git a/pyquda_utils/io/gauge_utils.py b/pyquda_io/_field_utils.py
similarity index 82%
rename from pyquda_utils/io/gauge_utils.py
rename to pyquda_io/_field_utils.py
index 3b91587..4ff5a58 100644
--- a/pyquda_utils/io/gauge_utils.py
+++ b/pyquda_io/_field_utils.py
@@ -3,7 +3,7 @@
 from mpi4py import MPI
 import numpy
 
-from .mpi_file import getSublatticeSize, getNeighbourRank
+from ._mpi_file import getSublatticeSize, getNeighbourRank
 
 Nd, Nc = 4, 3
 
@@ -72,30 +72,12 @@ def gaugeLexicoPlaquette(latt_size: List[int], grid_size: List[int], gauge: nump
         extended[:, -1, :-1, :-1, :-1] = buf
 
     plaq = numpy.zeros((6))
-    plaq[0] = numpy.vdot(
-        numpy.linalg.matmul(gauge[0], extended[1, :-1, :-1, :-1, 1:]),
-        numpy.linalg.matmul(gauge[1], extended[0, :-1, :-1, 1:, :-1]),
-    ).real
-    plaq[1] = numpy.vdot(
-        numpy.linalg.matmul(gauge[0], extended[2, :-1, :-1, :-1, 1:]),
-        numpy.linalg.matmul(gauge[2], extended[0, :-1, 1:, :-1, :-1]),
-    ).real
-    plaq[2] = numpy.vdot(
-        numpy.linalg.matmul(gauge[1], extended[2, :-1, :-1, 1:, :-1]),
-        numpy.linalg.matmul(gauge[2], extended[1, :-1, 1:, :-1, :-1]),
-    ).real
-    plaq[3] = numpy.vdot(
-        numpy.linalg.matmul(gauge[0], extended[3, :-1, :-1, :-1, 1:]),
-        numpy.linalg.matmul(gauge[3], extended[0, 1:, :-1, :-1, :-1]),
-    ).real
-    plaq[4] = numpy.vdot(
-        numpy.linalg.matmul(gauge[1], extended[3, :-1, :-1, 1:, :-1]),
-        numpy.linalg.matmul(gauge[3], extended[1, 1:, :-1, :-1, :-1]),
-    ).real
-    plaq[5] = numpy.vdot(
-        numpy.linalg.matmul(gauge[2], extended[3, :-1, 1:, :-1, :-1]),
-        numpy.linalg.matmul(gauge[3], extended[2, 1:, :-1, :-1, :-1]),
-    ).real
+    plaq[0] = numpy.vdot(gauge[0] @ extended[1, :-1, :-1, :-1, 1:], gauge[1] @ extended[0, :-1, :-1, 1:, :-1]).real
+    plaq[1] = numpy.vdot(gauge[0] @ extended[2, :-1, :-1, :-1, 1:], gauge[2] @ extended[0, :-1, 1:, :-1, :-1]).real
+    plaq[2] = numpy.vdot(gauge[1] @ extended[2, :-1, :-1, 1:, :-1], gauge[2] @ extended[1, :-1, 1:, :-1, :-1]).real
+    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)
@@ -186,3 +168,39 @@ def gaugeEvenShiftBackward(latt_size: List[int], grid_size: List[int], gauge: nu
         MPI.COMM_WORLD.Sendrecv_replace(buf, dest=neighbour_rank[3], source=neighbour_rank[7])
         gauge_shift[3, 1, 0, :, :, :] = buf
     return gauge_shift
+
+
+# matrices to convert gamma basis bewteen DeGrand-Rossi and Dirac-Pauli
+# DP for Dirac-Pauli, DR for DeGrand-Rossi
+# \psi(DP) = _DR_TO_DP \psi(DR)
+# \psi(DR) = _DP_TO_DR \psi(DP)
+_DP_TO_DR = numpy.array(
+    [
+        [0, 1, 0, -1],
+        [-1, 0, 1, 0],
+        [0, 1, 0, 1],
+        [-1, 0, -1, 0],
+    ]
+)
+_DR_TO_DP = numpy.array(
+    [
+        [0, -1, 0, -1],
+        [1, 0, 1, 0],
+        [0, 1, 0, -1],
+        [-1, 0, 1, 0],
+    ]
+)
+
+
+def propagatorDeGrandRossiToDiracPauli(propagator: numpy.ndarray):
+    P = _DR_TO_DP
+    Pinv = _DP_TO_DR / 2
+
+    return numpy.ascontiguousarray(numpy.einsum("ij,tzyxjkab,kl->tzyxilab", P, propagator.data, Pinv, optimize=True))
+
+
+def propagatorDiracPauliToDeGrandRossi(propagator: numpy.ndarray):
+    P = _DP_TO_DR
+    Pinv = _DR_TO_DP / 2
+
+    return numpy.ascontiguousarray(numpy.einsum("ij,tzyxjkab,kl->tzyxilab", P, propagator.data, Pinv, optimize=True))
diff --git a/pyquda_utils/io/mpi_file.py b/pyquda_io/_mpi_file.py
similarity index 100%
rename from pyquda_utils/io/mpi_file.py
rename to pyquda_io/_mpi_file.py
diff --git a/pyquda_utils/io/chroma.py b/pyquda_io/chroma.py
similarity index 98%
rename from pyquda_utils/io/chroma.py
rename to pyquda_io/chroma.py
index 4ea1872..49eca91 100644
--- a/pyquda_utils/io/chroma.py
+++ b/pyquda_io/chroma.py
@@ -4,7 +4,7 @@
 
 from mpi4py import MPI
 
-from .mpi_file import getSublatticeSize, getGridCoord, readMPIFile
+from ._mpi_file import getSublatticeSize, getGridCoord, readMPIFile
 
 Nd, Ns, Nc = 4, 4, 3
 _precision_map = {"D": 8, "F": 4, "S": 4}
diff --git a/pyquda_utils/io/io_general.py b/pyquda_io/io_general.py
similarity index 100%
rename from pyquda_utils/io/io_general.py
rename to pyquda_io/io_general.py
diff --git a/pyquda_utils/io/kyu.py b/pyquda_io/kyu.py
similarity index 88%
rename from pyquda_utils/io/kyu.py
rename to pyquda_io/kyu.py
index 19baede..4412e4c 100644
--- a/pyquda_utils/io/kyu.py
+++ b/pyquda_io/kyu.py
@@ -3,7 +3,8 @@
 
 import numpy
 
-from .mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
+from ._mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
+from ._field_utils import propagatorDiracPauliToDeGrandRossi, propagatorDeGrandRossiToDiracPauli
 
 Nd, Ns, Nc = 4, 4, 3
 
@@ -48,6 +49,7 @@ def readPropagator(filename: str, latt_size: List[int], grid_size: List[int]):
         .reshape(Lt, Lz, Ly, Lx, Ns, Ns, Nc, Nc * 2)
         .view("<c16")
     )
+    propagator = propagatorDiracPauliToDeGrandRossi(propagator)
     return propagator
 
 
@@ -56,6 +58,7 @@ def writePropagator(filename: str, latt_size: List[int], grid_size: List[int], p
     Lx, Ly, Lz, Lt = getSublatticeSize(latt_size, grid_size)
     dtype, offset = ">f8", 0
 
+    propagator = propagatorDeGrandRossiToDiracPauli(propagator)
     propagator = (
         propagator.view("<f8")
         .reshape(Lt, Lz, Ly, Lx, Ns, Ns, Nc, Nc, 2)
diff --git a/pyquda_utils/io/lime.py b/pyquda_io/lime.py
similarity index 100%
rename from pyquda_utils/io/lime.py
rename to pyquda_io/lime.py
diff --git a/pyquda_utils/io/milc.py b/pyquda_io/milc.py
similarity index 98%
rename from pyquda_utils/io/milc.py
rename to pyquda_io/milc.py
index edb249b..9d668fb 100644
--- a/pyquda_utils/io/milc.py
+++ b/pyquda_io/milc.py
@@ -7,7 +7,7 @@
 import numpy
 from mpi4py import MPI
 
-from .mpi_file import getGridCoord, getSublatticeSize, readMPIFile, writeMPIFile
+from ._mpi_file import getGridCoord, getSublatticeSize, readMPIFile, writeMPIFile
 
 Nd, Ns, Nc = 4, 4, 3
 _precision_map = {"D": 8, "F": 4, "S": 4}
diff --git a/pyquda_utils/io/nersc.py b/pyquda_io/nersc.py
similarity index 97%
rename from pyquda_utils/io/nersc.py
rename to pyquda_io/nersc.py
index b619e8e..8cccb60 100644
--- a/pyquda_utils/io/nersc.py
+++ b/pyquda_io/nersc.py
@@ -5,8 +5,8 @@
 import numpy
 from mpi4py import MPI
 
-from .mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
-from .gauge_utils import gaugeLexicoPlaquette
+from ._mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
+from ._field_utils import gaugeLexicoPlaquette
 
 Nd, Ns, Nc = 4, 4, 3
 
diff --git a/pyquda_utils/io/npy.py b/pyquda_io/npy.py
similarity index 97%
rename from pyquda_utils/io/npy.py
rename to pyquda_io/npy.py
index 667d128..75cdbd9 100644
--- a/pyquda_utils/io/npy.py
+++ b/pyquda_io/npy.py
@@ -5,7 +5,7 @@
 from numpy.lib.format import dtype_to_descr, read_magic, read_array_header_1_0, write_array_header_1_0
 from mpi4py import MPI
 
-from .mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
+from ._mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
 
 Nd, Ns, Nc = 4, 4, 3
 
diff --git a/pyquda_utils/io/openqcd.py b/pyquda_io/openqcd.py
similarity index 94%
rename from pyquda_utils/io/openqcd.py
rename to pyquda_io/openqcd.py
index c49195d..7c9151e 100644
--- a/pyquda_utils/io/openqcd.py
+++ b/pyquda_io/openqcd.py
@@ -5,8 +5,8 @@
 import numpy
 from mpi4py import MPI
 
-from .mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
-from .gauge_utils import gaugeEvenOdd, gaugeLexico, gaugeLexicoPlaquette, gaugeOddShiftForward, gaugeEvenShiftBackward
+from ._mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
+from ._field_utils import gaugeEvenOdd, gaugeLexico, gaugeLexicoPlaquette, gaugeOddShiftForward, gaugeEvenShiftBackward
 
 Nd, Ns, Nc = 4, 4, 3
 
diff --git a/pyquda_utils/io/xqcd.py b/pyquda_io/xqcd.py
similarity index 89%
rename from pyquda_utils/io/xqcd.py
rename to pyquda_io/xqcd.py
index 25e75bd..6e6c18d 100644
--- a/pyquda_utils/io/xqcd.py
+++ b/pyquda_io/xqcd.py
@@ -3,7 +3,8 @@
 
 import numpy
 
-from .mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
+from ._mpi_file import getSublatticeSize, readMPIFile, writeMPIFile
+from ._field_utils import propagatorDiracPauliToDeGrandRossi, propagatorDeGrandRossiToDiracPauli
 
 Ns, Nc = 4, 3
 
@@ -16,6 +17,7 @@ def readPropagator(filename: str, latt_size: List[int], grid_size: List[int], st
     if not staggered:
         propagator = readMPIFile(filename, dtype, offset, (Ns, Nc, Lt, Lz, Ly, Lx, Ns, Nc), (5, 4, 3, 2), grid_size)
         propagator = propagator.transpose(2, 3, 4, 5, 6, 0, 7, 1).astype("<c16")
+        propagator = propagatorDiracPauliToDeGrandRossi(propagator)
     else:
         # QDP_ALIGN16 makes the last Nc to be aligned with 16 Bytes.
         propagator_align16 = readMPIFile(filename, dtype, offset, (Nc, Lt, Lz, Ly, Lx, 4), (4, 3, 2, 1), grid_size)
@@ -32,6 +34,7 @@ def writePropagator(
     dtype, offset = "<c8", 0
 
     if not staggered:
+        propagator = propagatorDeGrandRossiToDiracPauli(propagator)
         propagator = propagator.astype(dtype).transpose(5, 7, 0, 1, 2, 3, 4, 6).copy()
         writeMPIFile(filename, dtype, offset, (Ns, Nc, Lt, Lz, Ly, Lx, Ns, Nc), (5, 4, 3, 2), grid_size, propagator)
     else:
diff --git a/pyquda_utils/gpt.py b/pyquda_utils/gpt.py
index 951e765..d356bba 100644
--- a/pyquda_utils/gpt.py
+++ b/pyquda_utils/gpt.py
@@ -49,7 +49,7 @@ def LatticeGaugeGPT(lattice: List[g.lattice], gen_simd_width: int, gauge: Lattic
     else:
         assert latt_info.size == gauge.latt_info.size
         for index in range(latt_info.Nd):
-            gpt_shape = [i for sl in zip(gpt_simd, gpt_latt) for i in sl]
+            gpt_shape = [i for sl in zip(gpt_simd[::-1], gpt_latt[::-1]) for i in sl]
             lattice[index].mview()[0][:] = (
                 gauge[index]
                 .lexico()
diff --git a/pyquda_utils/io/__init__.py b/pyquda_utils/io/__init__.py
index 61cbb13..1ff98ed 100644
--- a/pyquda_utils/io/__init__.py
+++ b/pyquda_utils/io/__init__.py
@@ -71,7 +71,7 @@ def rotateToDeGrandRossi(propagator: LatticePropagator):
 
 def readChromaQIOGauge(filename: str, checksum: bool = True):
     from pyquda import getGridSize
-    from .chroma import readQIOGauge as read
+    from pyquda_io.chroma import readQIOGauge as read
 
     latt_size, gauge_raw = read(filename, getGridSize(), checksum)
     return LatticeGauge(LatticeInfo(latt_size), evenodd(gauge_raw, [1, 2, 3, 4]))
@@ -79,7 +79,7 @@ def readChromaQIOGauge(filename: str, checksum: bool = True):
 
 def readILDGBinGauge(filename: str, dtype: str, latt_size: List[int]):
     from pyquda import getGridSize
-    from .chroma import readILDGBinGauge as read
+    from pyquda_io.chroma import readILDGBinGauge as read
 
     gauge_raw = read(filename, dtype, latt_size, getGridSize())
     return LatticeGauge(LatticeInfo(latt_size), evenodd(gauge_raw, [1, 2, 3, 4]))
@@ -87,7 +87,7 @@ def readILDGBinGauge(filename: str, dtype: str, latt_size: List[int]):
 
 def readChromaQIOPropagator(filename: str, checksum: bool = True):
     from pyquda import getGridSize
-    from .chroma import readQIOPropagator as read
+    from pyquda_io.chroma import readQIOPropagator as read
 
     latt_size, staggered, propagator_raw = read(filename, getGridSize(), checksum)
     if not staggered:
@@ -98,21 +98,21 @@ def readChromaQIOPropagator(filename: str, checksum: bool = True):
 
 def readMILCGauge(filename: str, checksum: bool = True):
     from pyquda import getGridSize
-    from .milc import readGauge as read
+    from pyquda_io.milc import readGauge as read
 
     latt_size, gauge_raw = read(filename, getGridSize(), checksum)
     return LatticeGauge(LatticeInfo(latt_size), evenodd(gauge_raw, [1, 2, 3, 4]))
 
 
 def writeMILCGauge(filename: str, gauge: LatticeGauge):
-    from .milc import writeGauge as write
+    from pyquda_io.milc import writeGauge as write
 
     write(filename, gauge.latt_info.global_size, gauge.latt_info.grid_size, gauge.lexico())
 
 
 def readMILCQIOPropagator(filename: str):
     from pyquda import getGridSize
-    from .milc import readQIOPropagator as read
+    from pyquda_io.milc import readQIOPropagator as read
 
     latt_size, staggered, propagator_raw = read(filename, getGridSize())
     if not staggered:
@@ -123,62 +123,50 @@ def readMILCQIOPropagator(filename: str):
 
 def readKYUGauge(filename: str, latt_size: List[int]):
     from pyquda import getGridSize
-    from .kyu import readGauge as read
+    from pyquda_io.kyu import readGauge as read
 
     gauge_raw = read(filename, latt_size, getGridSize())
     return LatticeGauge(LatticeInfo(latt_size), evenodd(gauge_raw, [1, 2, 3, 4]))
 
 
 def writeKYUGauge(filename: str, gauge: LatticeGauge):
-    from .kyu import writeGauge as write
+    from pyquda_io.kyu import writeGauge as write
 
     write(filename, gauge.latt_info.global_size, gauge.latt_info.grid_size, gauge.lexico())
 
 
 def readKYUPropagator(filename: str, latt_size: List[int]):
     from pyquda import getGridSize
-    from .kyu import readPropagator as read
+    from pyquda_io.kyu import readPropagator as read
 
     propagator_raw = read(filename, latt_size, getGridSize())
-    propagator = LatticePropagator(LatticeInfo(latt_size), evenodd(propagator_raw, [0, 1, 2, 3]))
-    propagator = rotateToDeGrandRossi(propagator)
-    return propagator
+    return LatticePropagator(LatticeInfo(latt_size), evenodd(propagator_raw, [0, 1, 2, 3]))
 
 
 def writeKYUPropagator(filename: str, propagator: LatticePropagator):
-    from .kyu import writePropagator as write
+    from pyquda_io.kyu import writePropagator as write
 
-    propagator = rotateToDiracPauli(propagator)
-    propagator_raw = propagator.lexico()
-    write(filename, propagator.latt_info.global_size, propagator.latt_info.grid_size, propagator_raw)
+    write(filename, propagator.latt_info.global_size, propagator.latt_info.grid_size, propagator.lexico())
 
 
 def readXQCDPropagator(filename: str, latt_size: List[int], staggered: bool):
     from pyquda import getGridSize
-    from .xqcd import readPropagator as read
+    from pyquda_io.xqcd import readPropagator as read
 
     propagator_raw = read(filename, latt_size, getGridSize(), staggered)
-    if not staggered:
-        return rotateToDeGrandRossi(LatticePropagator(LatticeInfo(latt_size), evenodd(propagator_raw, [0, 1, 2, 3])))
-    else:
-        return LatticeStaggeredPropagator(LatticeInfo(latt_size), evenodd(propagator_raw, [0, 1, 2, 3]))
+    return LatticeStaggeredPropagator(LatticeInfo(latt_size), evenodd(propagator_raw, [0, 1, 2, 3]))
 
 
 def writeXQCDPropagator(filename: str, propagator: Union[LatticePropagator, LatticeStaggeredPropagator]):
-    from .xqcd import writePropagator as write
+    from pyquda_io.xqcd import writePropagator as write
 
-    latt_size = propagator.latt_info.global_size
-    grid_size = propagator.latt_info.grid_size
     staggered = isinstance(propagator, LatticeStaggeredPropagator)
-    if not staggered:
-        write(filename, latt_size, grid_size, rotateToDiracPauli(propagator).lexico(), staggered)
-    else:
-        write(filename, latt_size, grid_size, propagator.lexico(), staggered)
+    write(filename, propagator.latt_info.global_size, propagator.latt_info.grid_size, propagator.lexico(), staggered)
 
 
 def readXQCDPropagatorFast(filename: str, latt_size: List[int]):
     from pyquda import getGridSize
-    from .xqcd import readPropagatorFast as read
+    from pyquda_io.xqcd import readPropagatorFast as read
 
     latt_info = LatticeInfo(latt_size)
     Lx, Ly, Lz, Lt = latt_info.size
@@ -192,7 +180,7 @@ def readXQCDPropagatorFast(filename: str, latt_size: List[int]):
 
 
 def writeXQCDPropagatorFast(filename: str, propagator: LatticePropagator):
-    from .xqcd import writePropagatorFast as write
+    from pyquda_io.xqcd import writePropagatorFast as write
 
     latt_info = propagator.latt_info
     Lx, Ly, Lz, Lt = latt_info.size
@@ -206,7 +194,7 @@ def writeXQCDPropagatorFast(filename: str, propagator: LatticePropagator):
 
 def readNPYGauge(filename: str):
     from pyquda import getGridSize
-    from .npy import readGauge as read
+    from pyquda_io.npy import readGauge as read
 
     filename = filename if filename.endswith(".npy") else filename + ".npy"
     latt_size, gauge_raw = read(filename, getGridSize())
@@ -214,7 +202,7 @@ def readNPYGauge(filename: str):
 
 
 def writeNPYGauge(filename: str, gauge: LatticeGauge):
-    from .npy import writeGauge as write
+    from pyquda_io.npy import writeGauge as write
 
     filename = filename if filename.endswith(".npy") else filename + ".npy"
     write(filename, gauge.latt_info.global_size, gauge.latt_info.grid_size, gauge.lexico())
@@ -222,42 +210,42 @@ def writeNPYGauge(filename: str, gauge: LatticeGauge):
 
 def readNPYPropagator(filename: str):
     from pyquda import getGridSize
-    from .npy import readPropagator as read
+    from pyquda_io.npy import readPropagator as read
 
     latt_size, propagator_raw = read(filename, getGridSize())
     return LatticePropagator(LatticeInfo(latt_size), evenodd(propagator_raw, [0, 1, 2, 3]))
 
 
 def writeNPYPropagator(filename: str, propagator: LatticePropagator):
-    from .npy import writePropagator as write
+    from pyquda_io.npy import writePropagator as write
 
     write(filename, propagator.latt_info.global_size, propagator.latt_info.grid_size, propagator.lexico())
 
 
 def readOpenQCDGauge(filename: str, plaquette: bool = True):
     from pyquda import getGridSize
-    from .openqcd import readGauge as read
+    from pyquda_io.openqcd import readGauge as read
 
     latt_size, gauge = read(filename, getGridSize(), plaquette, False)
     return LatticeGauge(LatticeInfo(latt_size), gauge)
 
 
 def writeOpenQCDGauge(filename: str, gauge: LatticeGauge, plaquette: float = None):
-    from .openqcd import writeGauge as write
+    from pyquda_io.openqcd import writeGauge as write
 
     write(filename, gauge.latt_info.global_size, gauge.latt_info.grid_size, gauge.getHost(), plaquette, False)
 
 
 def readNERSCGauge(filename: str, plaquette: bool = True, link_trace: bool = True, checksum: bool = True):
     from pyquda import getGridSize
-    from .nersc import readGauge as read
+    from pyquda_io.nersc import readGauge as read
 
     latt_size, gauge_raw = read(filename, getGridSize(), plaquette, link_trace, checksum)
     return LatticeGauge(LatticeInfo(latt_size), evenodd(gauge_raw, [1, 2, 3, 4]))
 
 
 def writeNERSCGauge(filename: str, gauge: LatticeGauge, plaquette: float = None, use_fp32: bool = False):
-    from .nersc import writeGauge as write
+    from pyquda_io.nersc import writeGauge as write
 
     write(filename, gauge.latt_info.global_size, gauge.latt_info.grid_size, gauge.lexico(), plaquette, use_fp32)
 
diff --git a/setup.py b/setup.py
index 0d6257a..967f0e4 100644
--- a/setup.py
+++ b/setup.py
@@ -11,5 +11,7 @@
             tag, post = describe, 0
         with open(os.path.join(os.path.dirname(__file__), "pyquda_utils", "_version.py"), "w") as f:
             f.write(f'__version__ = "{tag[1:]}.dev{post}"\n')
+        with open(os.path.join(os.path.dirname(__file__), "pyquda_io", "_version.py"), "w") as f:
+            f.write(f'__version__ = "{tag[1:]}.dev{post}"\n')
 
 setup()