diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py
index d02d6978..eecf1130 100644
--- a/pyerrors/__init__.py
+++ b/pyerrors/__init__.py
@@ -1,5 +1,6 @@
 from .pyerrors import *
 from . import correlators
+from . import dirac
 from . import fits
 from . import linalg
 from . import misc
diff --git a/pyerrors/dirac.py b/pyerrors/dirac.py
new file mode 100644
index 00000000..f5350b94
--- /dev/null
+++ b/pyerrors/dirac.py
@@ -0,0 +1,61 @@
+import numpy as np
+
+
+gammaX = np.array(
+    [[0, 0, 0, 1j], [0, 0, 1j, 0], [0, -1j, 0, 0], [-1j, 0, 0, 0]],
+    dtype=complex)
+gammaY = np.array(
+    [[0, 0, 0, -1], [0, 0, 1, 0], [0, 1, 0, 0], [-1, 0, 0, 0]],
+    dtype=complex)
+gammaZ = np.array(
+    [[0, 0, 1j, 0], [0, 0, 0, -1j], [-1j, 0, 0, 0], [0, 1j, 0, 0]],
+    dtype=complex)
+gammaT = np.array(
+    [[0, 0, 1, 0], [0, 0, 0, 1], [1, 0, 0, 0], [0, 1, 0, 0]],
+    dtype=complex)
+gamma = np.array([gammaX, gammaY, gammaZ, gammaT])
+gamma5 = np.array(
+    [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, -1]],
+    dtype=complex)
+identity = np.array(
+    [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]],
+    dtype=complex)
+
+
+def Grid_gamma(gamma_tag):
+    """Returns gamma matrix in Grid labeling."""
+    if gamma_tag == 'Identity':
+        g = identity
+    elif gamma_tag == 'Gamma5':
+        g = gamma5
+    elif gamma_tag == 'GammaX':
+        g = gamma[0]
+    elif gamma_tag == 'GammaY':
+        g = gamma[1]
+    elif gamma_tag == 'GammaZ':
+        g = gamma[2]
+    elif gamma_tag == 'GammaT':
+        g = gamma[3]
+    elif gamma_tag == 'GammaXGamma5':
+        g = gamma[0] @ gamma5
+    elif gamma_tag == 'GammaYGamma5':
+        g = gamma[1] @ gamma5
+    elif gamma_tag == 'GammaZGamma5':
+        g = gamma[2] @ gamma5
+    elif gamma_tag == 'GammaTGamma5':
+        g = gamma[3] @ gamma5
+    elif gamma_tag == 'SigmaXT':
+        g = 0.5 * (gamma[0] @ gamma[3] - gamma[3] @ gamma[0])
+    elif gamma_tag == 'SigmaXY':
+        g = 0.5 * (gamma[0] @ gamma[1] - gamma[1] @ gamma[0])
+    elif gamma_tag == 'SigmaXZ':
+        g = 0.5 * (gamma[0] @ gamma[2] - gamma[2] @ gamma[0])
+    elif gamma_tag == 'SigmaYT':
+        g = 0.5 * (gamma[1] @ gamma[3] - gamma[3] @ gamma[1])
+    elif gamma_tag == 'SigmaYZ':
+        g = 0.5 * (gamma[1] @ gamma[2] - gamma[2] @ gamma[1])
+    elif gamma_tag == 'SigmaZT':
+        g = 0.5 * (gamma[2] @ gamma[3] - gamma[3] @ gamma[2])
+    else:
+        raise Exception('Unkown gamma structure', gamma_tag)
+    return g
diff --git a/pyerrors/npr.py b/pyerrors/npr.py
index 98c99481..ad7d7de8 100644
--- a/pyerrors/npr.py
+++ b/pyerrors/npr.py
@@ -2,70 +2,12 @@ import warnings
 import numpy as np
 import autograd.numpy as anp
 from .linalg import mat_mat_op
+from .dirac import gamma, gamma5
 
 
 L = None
 T = None
 
-gammaX = np.array(
-    [[0, 0, 0, 1j], [0, 0, 1j, 0], [0, -1j, 0, 0], [-1j, 0, 0, 0]],
-    dtype=complex)
-gammaY = np.array(
-    [[0, 0, 0, -1], [0, 0, 1, 0], [0, 1, 0, 0], [-1, 0, 0, 0]],
-    dtype=complex)
-gammaZ = np.array(
-    [[0, 0, 1j, 0], [0, 0, 0, -1j], [-1j, 0, 0, 0], [0, 1j, 0, 0]],
-    dtype=complex)
-gammaT = np.array(
-    [[0, 0, 1, 0], [0, 0, 0, 1], [1, 0, 0, 0], [0, 1, 0, 0]],
-    dtype=complex)
-gamma = np.array([gammaX, gammaY, gammaZ, gammaT])
-gamma5 = np.array(
-    [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, -1]],
-    dtype=complex)
-identity = np.array(
-    [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]],
-    dtype=complex)
-
-
-def Grid_gamma(gamma_tag):
-    """Returns gamma matrix in Grid labeling."""
-    if gamma_tag == 'Identity':
-        g = identity
-    elif gamma_tag == 'Gamma5':
-        g = gamma5
-    elif gamma_tag == 'GammaX':
-        g = gamma[0]
-    elif gamma_tag == 'GammaY':
-        g = gamma[1]
-    elif gamma_tag == 'GammaZ':
-        g = gamma[2]
-    elif gamma_tag == 'GammaT':
-        g = gamma[3]
-    elif gamma_tag == 'GammaXGamma5':
-        g = gamma[0] @ gamma5
-    elif gamma_tag == 'GammaYGamma5':
-        g = gamma[1] @ gamma5
-    elif gamma_tag == 'GammaZGamma5':
-        g = gamma[2] @ gamma5
-    elif gamma_tag == 'GammaTGamma5':
-        g = gamma[3] @ gamma5
-    elif gamma_tag == 'SigmaXT':
-        g = 0.5 * (gamma[0] @ gamma[3] - gamma[3] @ gamma[0])
-    elif gamma_tag == 'SigmaXY':
-        g = 0.5 * (gamma[0] @ gamma[1] - gamma[1] @ gamma[0])
-    elif gamma_tag == 'SigmaXZ':
-        g = 0.5 * (gamma[0] @ gamma[2] - gamma[2] @ gamma[0])
-    elif gamma_tag == 'SigmaYT':
-        g = 0.5 * (gamma[1] @ gamma[3] - gamma[3] @ gamma[1])
-    elif gamma_tag == 'SigmaYZ':
-        g = 0.5 * (gamma[1] @ gamma[2] - gamma[2] @ gamma[1])
-    elif gamma_tag == 'SigmaZT':
-        g = 0.5 * (gamma[2] @ gamma[3] - gamma[3] @ gamma[2])
-    else:
-        raise Exception('Unkown gamma structure', gamma_tag)
-    return g
-
 
 class Npr_matrix(np.ndarray):
 
diff --git a/tests/test_npr.py b/tests/test_dirac.py
similarity index 58%
rename from tests/test_npr.py
rename to tests/test_dirac.py
index fa5cec63..28222da6 100644
--- a/tests/test_npr.py
+++ b/tests/test_dirac.py
@@ -6,7 +6,7 @@ np.random.seed(0)
 
 
 def test_gamma_matrices():
-    for matrix in pe.npr.gamma:
+    for matrix in pe.dirac.gamma:
         assert np.allclose(matrix @ matrix, np.identity(4))
         assert np.allclose(matrix, matrix.T.conj())
-    assert np.allclose(pe.npr.gamma5, pe.npr.gamma[0] @ pe.npr.gamma[1] @ pe.npr.gamma[2] @ pe.npr.gamma[3])
+    assert np.allclose(pe.dirac.gamma5, pe.dirac.gamma[0] @ pe.dirac.gamma[1] @ pe.dirac.gamma[2] @ pe.dirac.gamma[3])