diff --git a/lyncs_quda/clover_field.py b/lyncs_quda/clover_field.py
index ce89a84..ae4838b 100644
--- a/lyncs_quda/clover_field.py
+++ b/lyncs_quda/clover_field.py
@@ -28,17 +28,41 @@ class CloverField(LatticeField):
      This is designed as an intermediary to QUDA CloverField class
      so that it should have 1-to-1 correspondence to an QUDA instance.
     Note:
-     * This class stores the corresponding gauge field in its "field" attribute
-        to make except-clause of copy() work
      * direct & inverse fields are both allocated upon initialization
      * Only rho is mutable.  To change other params, a new instance should be created
      * QUDA convention for clover field := 1+i ( kappa csw )/4 sigma_mu,nu F_mu,nu (<-sigma_mu,nu: spinor tensor)
      *  so that sigma_mu,nu = i[g_mu, g_nu], F_mu,nu = (Q_mu,nu - Q_nu,mu)/8 (1/2 is missing from sigma_mu,nu)
+     * Apparently, an input to QUDA clover object, coeff = kappa*csw
+     *  wihout a normalization factor of 1/4 or 1/32 (suggested in interface_quda.cpp)
     """
 
+    def __new__(cls, fmunu, **kwargs):
+        # TODO: get dofs and local dims from kwargs, instead of getting them
+        # from self.shape assuming that it has the form (dofs, local_dims)
+        if isinstance(fmunu, CloverField):
+            return fmunu
+        if not isinstance(fmunu, (numpy.ndarray, cupy.ndarray)):
+            raise TypeError(
+                f"Supporting only numpy or cupy for field, got {type(fmunu)}"
+            )
+        is_clover = False
+        if not isinstance(fmunu, GaugeField):
+            if kwargs.get("is_clover", False):
+                is_clover = True
+                field = fmunu
+            else:
+                fmunu = GaugeField(fmunu)
+
+        if not is_clover:  # not copying from a clover-field array
+            idof = int((fmunu.ncol * fmunu.ndims) ** 2 / 2)
+            prec = fmunu.dtype
+            field = fmunu.backend.empty((idof,) + fmunu.dims, dtype=prec)
+
+        return super().__new__(cls, field, **kwargs)
+
     def __init__(
         self,
-        fmunu,
+        obj,
         coeff=0.0,
         twisted=False,
         mu2=0,
@@ -46,69 +70,46 @@ def __init__(
         eps2=0,
         rho=0,
         computeTrLog=False,
+        **kwargs,
     ):
-        # ASSUME: coeff = kappa*csw wihout a normalization factor of 1/4 or 1/32 (suggested in interface_quda.cpp)
-        # ? better to store fmunu.quda_field to _fmunu -> import gauge_tensor to be used in some methods
-        # ? better to put clover into self.field -> need walk-around to make copy() work
-        if not isinstance(fmunu, GaugeField):
-            fmunu = GaugeField(fmunu)
-        self._fmunu = fmunu.compute_fmunu()
-        super().__init__(self._fmunu.field, comm=self._fmunu.comm)
-
-        # QUDA clover field inherently works with real's not with complex's (c.f., include/clover_field_order.h)
-
-        idof = int((self._fmunu.ncol * self._fmunu.ndims) ** 2 / 2)
-        prec = self._fmunu.precision
-        self._direct = (
-            False  # Here, it is a flag to indicate whether the field has been computed
-        )
-        self._inverse = (
-            False  # Here, it is a flag to indicate whether the field has been computed
-        )
-
-        new = lambda idof: LatticeField.create(
-            self._fmunu.global_lattice,
-            dofs=(idof,),
-            dtype=prec,
-            device=self._fmunu.device,
-            empty=True,
-        )
-        self._clover = new(idof)
-        self._cloverInv = new(idof)
-        self.coeff = coeff
-        self._twisted = twisted
-        self._twist_flavor = tf
-        self._mu2 = mu2
-        self._eps2 = eps2
-        self._rho = rho
-        self.computeTrLog = computeTrLog
-
-    # shape, dofs, dtype, iscomlex, isreal are overwriten to report their values for the clover field, instead of _fmunu
-
-    @property
-    def shape(self):
-        "Shape of the clover field"
-        return self._clover.shape
-
-    @property
-    def dofs(self):
-        "Shape of the per-site degrees of freedom"
-        return self._clover.dofs
-
-    @property
-    def dtype(self):
-        "Clover field data type"
-        return self._clover.dtype
-
-    @property
-    def iscomplex(self):
-        "Whether the clover field dtype is complex"
-        return self._clover.iscomplex
+        # WARNING: ndarray object is not supposed to be view-casted to CloverField object
+        #           except in __new__, for which __init__ will be called subsequently,
+        #           as the result won't come with  meta info such as 'coeff' or 'mu2'
+
+        super().__init__(obj, getattr(obj, "comm", None))
+        if isinstance(obj, GaugeField):
+            # explicit construction
+            # QUDA clover field inherently works with real's not with complex's (c.f., include/clover_field_order.h)
+            self._cloverInv = LatticeField.create(
+                self.global_lattice,
+                dofs=self.dofs,
+                dtype=self.dtype,
+                device=self.device_id,
+                empty=True,
+            )
+            self._fmunu = obj.compute_fmunu()
+            self._direct = False  # Here, it is a flag to indicate whether the field has been computed
+            self._inverse = False  # Here, it is a flag to indicate whether the field has been computed
+            self.coeff = coeff
+            self._twisted = twisted
+            self._twist_flavor = tf
+            self._mu2 = mu2
+            self._eps2 = eps2
+            self._rho = rho
+            self.computeTrLog = computeTrLog
+        elif isinstance(obj, CloverField):
+            # upcasting to ndarray or new from template
+            self.__dict__.update(obj.__dict__)
+        elif isinstance(obj, self.backend.ndarray):
+            pass
+        else:
+            raise ValueError(
+                "The input is expected to be ndarray or LatticeField object"
+            )
 
-    @property
-    def isreal(self):
-        "Whether the clover field dtype is real"
-        return self._clover.isreal
+    def _prepare(self, field, copy=False, check=False, **kwargs):
+        # When CloverField object prepares its input, the input is assumed to be of CloverField
+        return super()._prepare(field, copy=copy, check=check, is_clover=True, **kwargs)
 
     # naming suggestion: native_view? default_* saved for dofs+lattice?
     def default_view(self):
@@ -116,7 +117,7 @@ def default_view(self):
         shape = (2,)  # even-odd
         shape += (self.dofs[0] // N, -1, N)
 
-        return self.field.view().reshape(shape)
+        return self.float_view().reshape(shape)
 
     @property
     def twisted(self):
@@ -178,7 +179,7 @@ def quda_params(self):
         params = lib.CloverFieldParam()
         lib.copy_struct(params, super().quda_params)
         params.inverse = True
-        params.clover = to_pointer(self._clover.ptr)
+        params.clover = to_pointer(self.ptr)
         params.cloverInv = to_pointer(self._cloverInv.ptr)
         params.coeff = self.coeff
         params.twisted = self.twisted
@@ -205,7 +206,7 @@ def clover_field(self):
         if not self._direct:
             lib.computeClover(self.quda_field, self._fmunu.quda_field, self.coeff)
             self._direct = True
-        return self._clover.field
+        return self.view(self.backend.ndarray)
 
     @property
     def inverse_field(self):
@@ -213,7 +214,7 @@ def inverse_field(self):
             self.clover_field
             lib.cloverInvert(self.quda_field, self.computeTrLog)
             self._inverse = True
-        return self._cloverInv.field
+        return self._cloverInv.view(self.backend.ndarray)
 
     @property
     def trLog(self):
diff --git a/lyncs_quda/dirac.py b/lyncs_quda/dirac.py
index f35b02a..5de9d98 100644
--- a/lyncs_quda/dirac.py
+++ b/lyncs_quda/dirac.py
@@ -242,9 +242,9 @@ def action(self, phi, **params):
         out = 0
         if not self.full and "CLOVER" in self.type:
             if self.even:
-                out -= 2 * self.clover.trLog[1]
-            else:
                 out -= 2 * self.clover.trLog[0]
+            else:
+                out -= 2 * self.clover.trLog[1]
 
         parity = None
         if not self.full:
diff --git a/lyncs_quda/gauge_field.py b/lyncs_quda/gauge_field.py
index 65c83ba..25f629a 100644
--- a/lyncs_quda/gauge_field.py
+++ b/lyncs_quda/gauge_field.py
@@ -74,11 +74,23 @@ def momentum(lattice, **kwargs):
 class GaugeField(LatticeField):
     "Mimics the quda::GaugeField object"
 
-    @LatticeField.field.setter
-    def field(self, field):
-        LatticeField.field.fset(self, field)
-        if self.reconstruct == "INVALID":
-            raise TypeError(f"Unrecognized field dofs {self.dofs}")
+    def _check_field(self, field=None):
+        # NOTE:
+        #  - For now, we store dofs and local lattice info in the shape of the array
+        #  - The canonical form for GaugeField is [dofs, local_lattice]
+        #     where len(dofs) == 1,2 and len(local_lattice) = 4
+        super()._check_field(field)
+        if field is None:
+            field = self
+        dofs = field.shape[: -self.ndims]
+        pdofs = (
+            prod(dofs)
+            if dofs[0] not in self._geometry_values[1] or dofs[0] == 1
+            else prod(dofs[1:])
+        )
+        pdofs *= 2 ** (self.iscomplex)
+        if not (pdofs in (12, 8, 10) or sqrt(pdofs / 2).is_integer()):
+            raise TypeError(f"Unrecognized field dofs {dofs}")
 
     def new(self, reconstruct=None, geometry=None, **kwargs):
         "Returns a new empty field based on the current"
@@ -139,11 +151,6 @@ def copy(self, other=None, out=None, **kwargs):
 
         return super().copy(other, out, **kwargs)
 
-    def __array_finalize__(self, obj):
-        "Support for __array_finalize__ standard"
-        # need to reset QUDA object when meta data of its Python wrapper is changed
-        self._quda = None
-
     def _prepare(self, field, **kwargs):
         field = super()._prepare(field, **kwargs)
         is_momentum = kwargs.get("is_momentum", self.is_momentum)
@@ -254,7 +261,7 @@ def is_momentum(self):
     def is_momentum(self, value):
         self._is_momentum = value
         if self._quda is not None:
-            self._quda.link_type = self.quda_link_type
+            self._quda = None
 
     @property
     def t_boundary(self):
@@ -523,7 +530,7 @@ def project(self, tol=None):
         if tol is None:
             tol = numpy.finfo(self.dtype).eps
 
-        assert self.device == cupy.cuda.runtime.getDevice()
+        assert self.device_id == cupy.cuda.runtime.getDevice()
         fails = cupy.zeros((1,), dtype="int32")
         lib.projectSU3(self.quda_field, tol, to_pointer(fails.data.ptr, "int *"))
         # return fails.get()[0]  # shouldn't we reduce?
diff --git a/lyncs_quda/lattice_field.py b/lyncs_quda/lattice_field.py
index 15c2d92..24476c6 100644
--- a/lyncs_quda/lattice_field.py
+++ b/lyncs_quda/lattice_field.py
@@ -165,39 +165,36 @@ def new(self, empty=True, **kwargs):
             self.global_lattice,
             dofs=kwargs.pop("dofs", self.dofs),
             dtype=kwargs.pop("dtype", self.dtype),
-            device=kwargs.pop("device", self.device),
+            device=kwargs.pop("device", self.device_id),
             comm=kwargs.pop("comm", self.comm),
             empty=empty,
             **kwargs,
         )
-        out.__array_finalize__(self)
         return out
 
-    def __array_finalize__(self, obj):
-        "Support for __array_finalize__ standard"
-        if self.comm is None:
-            self.comm = obj.comm
-
     def copy(self, other=None, out=None, **kwargs):
         "Returns out, a copy+=kwargs, of other if given, else of self"
         # src: other if given; otherwise self
         # dst: out (created anew, if not explicitly given)
         # ASSUME: if out != None, out <=> other/self+=kwargs
         #
-        # IF other and out are both given, this behaves like a classmethod
+        # If other and out are both given, this behaves like a classmethod
         # where out&other are casted into type(self)
+        # TODO:
+        #  * check if this dose not cause any bugs if it overwrites ndarray.copy
+        #     - For now, we cast self to ndarray before performing ndarray methods like flatten
+        #     - the second arg should be "order='C'" to match the signiture?
 
         # check=False => here any output is accepted
         out = self.prepare_out(out, check=False, **kwargs)
-
         if other is None:
             other = self
         # we prepare other without copying because we do copy here!
         other = out.prepare_in(other, copy=False, check=False, **kwargs)
         try:
             out.quda_field.copy(other.quda_field)
-        except:  # NotImplementedError:  #raised if self is LatticeField# at least, serial version calls exit(1) from qudaError, which is not catched by this
-            # As last resort trying to copy elementwise
+        except:  # NotImplementedError:  #raised if self is LatticeField# at least, serial version calls exit(1) from qudaError, which is not caught by this
+            # As a last resort, trying to copy elementwise
             out.default_view()[:] = other.default_view()
 
         return out
@@ -211,8 +208,8 @@ def equivalent(self, other, **kwargs):
         dtype = kwargs.get("dtype", self.dtype)
         if other.dtype != dtype:
             return False
-        device = kwargs.get("device", self.device)
-        if other.device != device:
+        device_id = kwargs.get("device", self.device_id)
+        if other.device_id != device_id:
             return False
         dofs = kwargs.get("dofs", self.dofs)
         if other.dofs != dofs:
@@ -226,12 +223,11 @@ def _prepare(self, field, copy=False, check=False, **kwargs):
             return self.new(**kwargs)
         cls = type(self)
         if not isinstance(field, cls):
-            field = cls(field)
+            field = cls(field, **kwargs)
         if check and not self.equivalent(field, **kwargs):
             if copy:
                 return self.copy(other=field, **kwargs)
             raise ValueError("The given field is not appropriate")
-        field.__array_finalize__(self)
         return field
 
     def prepare(self, fields, **kwargs):
@@ -259,8 +255,45 @@ def prepare_in(self, fields, **kwargs):
         kwargs.setdefault("copy", True)
         return self.prepare(fields, **kwargs)
 
+    _children = {}
+
+    def __new__(cls, field, **kwargs):
+        # TODO: get dofs and local dims from kwargs, instead of getting them
+        # from self.shape assuming that it has the form (dofs, local_dims)
+
+        if isinstance(field, cls):
+            return field
+        if not isinstance(field, (numpy.ndarray, cupy.ndarray)):
+            raise TypeError(
+                f"Supporting only numpy or cupy for field, got {type(field)}"
+            )
+        parent = numpy.ndarray if isinstance(field, numpy.ndarray) else cupy.ndarray
+        child = (cls, parent)
+        if child not in cls._children:
+            cls._children[child] = type(cls.__name__ + "ext", child, {})
+        obj = field.view(type=cls._children[child])
+        # self._dims = kwargs.get("dims", self.shape[-self.ndims :])
+        # self._dofs = kwargs.get("dofs", field.shape[: -self.ndims])
+
+        return obj
+
+    # field check should be performed
+    def __array_finalize__(self, obj):
+        "Support for __array_finalize__ standard"
+        # Note: this is called when creating a temporary, possibly causing
+        # some issues.  Apparently, this temporary is flattened and loose
+        # shape information
+        # self: newly created instance
+        # obj: input instance
+        self._check_field(obj)
+        if obj is None:
+            return  # can be removed; we are not bypassing ndarray.__new__
+        # This will require some care when we use attr for dims and dofs in _check_field
+        # self._dims = kwargs.get("dims", self.shape[-self.ndims :])
+        # self._dofs = kwargs.get("dofs", field.shape[: -self.ndims])
+        self.__init__(obj, comm=getattr(obj, "comm", None))
+
     def __init__(self, field, comm=None, **kwargs):
-        self.field = field
         if comm is False:
             self.comm = comm
         elif comm is None:
@@ -270,52 +303,41 @@ def __init__(self, field, comm=None, **kwargs):
         self._quda = None
         self.activate()
 
+    def _check_field(self, field=None):
+        if field is None:
+            field = self
+        if isinstance(field, cupy.ndarray) and field.device.id != lib.device_id:
+            raise TypeError("Field is stored on a different device than the quda lib")
+        if len(field.shape) < 4:
+            raise ValueError("A lattice field should not have shape smaller than 4")
+
     def activate(self):
         "Activates the current field. To be called before using the object in quda"
         "to make sure the communicator is set for MPI"
         # if self.comm is None, but #ranks>1, this will mess thigns up
         lib.set_comm(self.comm)
 
-    @property
-    def field(self):
-        "The underlaying lattice field"
-        return self._field
-
-    @field.setter
-    def field(self, field):
-        if isinstance(field, LatticeField):
-            field = field.field
-        if not isinstance(field, (numpy.ndarray, cupy.ndarray)):
-            raise TypeError(
-                f"Supporting only numpy or cupy for field, got {type(field)}"
-            )
-        if isinstance(field, cupy.ndarray) and field.device.id != lib.device_id:
-            raise TypeError("Field is stored on a different device than the quda lib")
-        if len(field.shape) < 4:
-            raise ValueError("A lattice field should not have shape smaller than 4")
-        self._field = field
-
-    def get(self):
+    def get(self, *args, **kwargs):
         "Returns the field as numpy array"
-        return self.__array__()
+        return self.__array__(*args, **kwargs)
 
     def __array__(self, *args, **kwargs):
-        out = self.field
         if self.device is not None:
-            out = out.get()
-        return out.__array__(*args, **kwargs)
+            return super().get(*args, **kwargs)
+        return super().__array__(*args, **kwargs)
 
     def complex_view(self):
         "Returns a complex view of the field"
         if self.iscomplex:
-            return self.field
-        return self.field.view(get_complex_dtype(self.dtype))
+            return self.view(type=self.backend.ndarray)
+        return self.view(get_complex_dtype(self.dtype), self.backend.ndarray)
 
     def float_view(self):
         "Returns a complex view of the field"
+        # don't need to upcast if we keep dofs and dims as attributes
         if not self.iscomplex:
-            return self.field
-        return self.field.view(get_float_dtype(self.dtype))
+            return self.view(type=self.backend.ndarray)
+        return self.view(get_float_dtype(self.dtype), self.backend.ndarray)
 
     def default_view(self):
         "Returns the default view of the field including reshaping"
@@ -324,26 +346,25 @@ def default_view(self):
     @property
     def backend(self):
         "The backend of the field: cupy or numpy"
-        if isinstance(self.field, cupy.ndarray):
+        if isinstance(self, cupy.ndarray):
             return cupy
         return numpy
 
     @property
     def device(self):
         "Device id of the field (None if not on GPUs)"
-        if isinstance(self.field, cupy.ndarray):
-            return self.field.device.id
+        if isinstance(self, cupy.ndarray):
+            return super().device
         return None
 
     @property
-    def shape(self):
-        "Shape of the field"
-        return self.field.shape
+    def device_id(self):
+        return getattr(self.device, "id", None)
 
     @property
     def location(self):
         "Memory location of the field (CPU or CUDA)"
-        return "CPU" if isinstance(self.field, numpy.ndarray) else "CUDA"
+        return "CPU" if isinstance(self, numpy.ndarray) else "CUDA"
 
     @property
     def quda_location(self):
@@ -382,20 +403,15 @@ def dofs(self):
         "Shape of the per-site degrees of freedom"
         return self.shape[: -self.ndims]
 
-    @property
-    def dtype(self):
-        "Field data type"
-        return self.field.dtype
-
     @property
     def iscomplex(self):
         "Whether the field dtype is complex"
-        return self.backend.iscomplexobj(self.field)
+        return self.backend.iscomplexobj(self)
 
     @property
     def isreal(self):
         "Whether the field dtype is real"
-        return self.backend.isrealobj(self.field)
+        return self.backend.isrealobj(self)
 
     @property
     def precision(self):
@@ -425,9 +441,9 @@ def pad(self):
     @property
     def ptr(self):
         "Memory pointer"
-        if isinstance(self.field, numpy.ndarray):
-            return self.field.__array_interface__["data"][0]
-        return self.field.data.ptr
+        if isinstance(self, numpy.ndarray):
+            return self.__array_interface__["data"][0]
+        return self.data.ptr
 
     @staticmethod
     @cache
@@ -480,9 +496,8 @@ def gpu_field(self):
 
     def __array_ufunc__(self, ufunc, method, *args, **kwargs):
         out = kwargs.get("out", (self,))[0]
-
         prepare = (
-            lambda arg: out.prepare_in(arg).field
+            lambda arg: out.prepare_in(arg).view(type=self.backend.ndarray)
             if isinstance(arg, (LatticeField, cupy.ndarray, numpy.ndarray))
             else arg
         )
@@ -495,8 +510,10 @@ def __array_ufunc__(self, ufunc, method, *args, **kwargs):
                 kwargs[key] = prepare(val)
 
         fnc = getattr(ufunc, method)
-
-        return self.prepare_out(fnc(*args, **kwargs), check=False)
+        result = fnc(*args, **kwargs)
+        if not isinstance(result, self.backend.ndarray):
+            return result
+        return self.prepare_out(result, check=False)
 
     def __bool__(self):
-        return bool(self.field.all())
+        return bool(self.all())
diff --git a/lyncs_quda/lib.py b/lyncs_quda/lib.py
index bbd47df..2c68901 100644
--- a/lyncs_quda/lib.py
+++ b/lyncs_quda/lib.py
@@ -295,7 +295,9 @@ def __del__(self):
     @fixture(scope="session")
     def fixlib():
         "A fixture to guarantee that in pytest lib is finalized at the end"
-        if not lib.initialized:
+        if QUDA_MPI and MPI.COMM_WORLD.Get_size() > 1:
+            pass
+        elif not lib.initialized:
             lib.init_quda()
         yield lib
         if lib.initialized:
diff --git a/lyncs_quda/testing.py b/lyncs_quda/testing.py
index 99a25b4..df89fec 100644
--- a/lyncs_quda/testing.py
+++ b/lyncs_quda/testing.py
@@ -9,6 +9,7 @@
     "parallel_loop",
     "mtype_loop",
     "mark_mpi",
+    "epsilon_loop",
 ]
 
 from random import random
diff --git a/test/test_clover_field.py b/test/test_clover_field.py
index 967f80c..6019237 100644
--- a/test/test_clover_field.py
+++ b/test/test_clover_field.py
@@ -31,7 +31,6 @@ def test_default(lattice):
 @lattice_loop  # enables lattice
 def test_params(lib, lattice, device, dtype):
     gf = gauge(lattice, dtype=dtype, device=device)
-    # gf.zero()
     clv = CloverField(gf, computeTrLog=True, coeff=0)
     params = clv.quda_params
 
@@ -117,7 +116,9 @@ def test_unit(lib, lattice, device, dtype):
         gf, coeff=1.0, tf="SINGLET", twisted=True, mu2=mu2, computeTrLog=True
     )
 
-    assert (clv.field == 0).all()
+    clv.fill(0)
+    assert clv == 0
+    assert clv + 0 == 0
     idof = int(((clv.ncol * clv.nspin) ** 2 / 2))
     if dtype == "float64":
         tmp = np.zeros((idof,) + lattice, dtype=dtype).reshape((2, 2, 36, -1))
diff --git a/test/test_dirac.py b/test/test_dirac.py
index af931b3..6917b01 100644
--- a/test/test_dirac.py
+++ b/test/test_dirac.py
@@ -1,7 +1,9 @@
 from random import random
 import numpy as np
+from lyncs_utils import isclose
 from lyncs_quda import (
     gauge,
+    momentum,
     spinor,
     gauge_coarse,
     gauge_scalar,
@@ -15,6 +17,7 @@
     dtype_loop,
     mu_loop,
     gamma_loop,
+    epsilon_loop,
 )
 
 
@@ -65,25 +68,25 @@ def test_zero(lib, lattice, device, gamma, mu, dtype=None):
     sf.uniform()
     kappa = random()
     dirac = gf.Dirac(kappa=kappa)
-    assert (dirac.M(sf).field == sf.field).all()
-    assert (dirac.Mdag(sf).field == sf.field).all()
-    assert (dirac.MdagM(sf).field == sf.field).all()
-    assert (dirac.MMdag(sf).field == sf.field).all()
+    assert (dirac.M(sf) == sf).all()
+    assert (dirac.Mdag(sf) == sf).all()
+    assert (dirac.MdagM(sf) == sf).all()
+    assert (dirac.MMdag(sf) == sf).all()
 
     dirac = gf.Dirac(kappa=kappa, mu=mu)
-    sfmu = (2 * kappa * mu) * 1j * sf.gamma5().field
-    assert np.allclose(dirac.M(sf).field, sf.field + sfmu)
-    assert np.allclose(dirac.Mdag(sf).field, sf.field - sfmu)
-    assert np.allclose(dirac.MdagM(sf).field, (1 + (2 * kappa * mu) ** 2) * sf.field)
-    assert np.allclose(dirac.MMdag(sf).field, (1 + (2 * kappa * mu) ** 2) * sf.field)
+    sfmu = (2 * kappa * mu) * 1j * sf.gamma5()
+    assert np.allclose(dirac.M(sf), sf + sfmu)
+    assert np.allclose(dirac.Mdag(sf), sf - sfmu)
+    assert np.allclose(dirac.MdagM(sf), (1 + (2 * kappa * mu) ** 2) * sf)
+    assert np.allclose(dirac.MMdag(sf), (1 + (2 * kappa * mu) ** 2) * sf)
 
     csw = random()
     dirac = gf.Dirac(kappa=kappa, mu=mu, csw=csw)
-    sfmu = (2 * kappa * mu) * 1j * sf.gamma5().field
-    assert np.allclose(dirac.M(sf).field, sf.field + sfmu)
-    assert np.allclose(dirac.Mdag(sf).field, sf.field - sfmu)
-    assert np.allclose(dirac.MdagM(sf).field, (1 + (2 * kappa * mu) ** 2) * sf.field)
-    assert np.allclose(dirac.MMdag(sf).field, (1 + (2 * kappa * mu) ** 2) * sf.field)
+    sfmu = (2 * kappa * mu) * 1j * sf.gamma5()
+    assert np.allclose(dirac.M(sf), sf + sfmu)
+    assert np.allclose(dirac.Mdag(sf), sf - sfmu)
+    assert np.allclose(dirac.MdagM(sf), (1 + (2 * kappa * mu) ** 2) * sf)
+    assert np.allclose(dirac.MMdag(sf), (1 + (2 * kappa * mu) ** 2) * sf)
 
 
 # @dtype_loop  # enables dtype #Double precision multigrid has not been enabled
@@ -99,10 +102,10 @@ def test_coarse_zero(lib, lattice, device, dtype=None):
     sf.uniform()
 
     dirac = gf.Dirac(clover=gf2)
-    assert (dirac.M(sf).field == sf.field).all()
-    assert (dirac.Mdag(sf).field == sf.field).all()
-    assert (dirac.MdagM(sf).field == sf.field).all()
-    assert (dirac.MMdag(sf).field == sf.field).all()
+    assert (dirac.M(sf) == sf).all()
+    assert (dirac.Mdag(sf) == sf).all()
+    assert (dirac.MdagM(sf) == sf).all()
+    assert (dirac.MMdag(sf) == sf).all()
 
 
 # @dtype_loop  # enables dtype
diff --git a/test/test_field.py b/test/test_field.py
index 3e0fd0b..a600455 100644
--- a/test/test_field.py
+++ b/test/test_field.py
@@ -22,6 +22,7 @@ def test_precision():
 
 def test_numpy():
     field = LatticeField(np.zeros(shape))
+
     assert field.location == "CPU"
     assert field.dims == shape[-4:]
     assert field.dofs == shape[:-4]
diff --git a/test/test_gauge_field.py b/test/test_gauge_field.py
index 980f427..63c2fca 100644
--- a/test/test_gauge_field.py
+++ b/test/test_gauge_field.py
@@ -10,7 +10,7 @@
     epsilon_loop,
 )
 from lyncs_cppyy.ll import addressof
-from lyncs_utils import isclose, allclose
+from lyncs_utils import isclose
 
 
 @lattice_loop
@@ -60,7 +60,6 @@ def test_zero(lib, lattice, device, dtype):
 
     assert gf.project() == 4 * np.prod(lattice)
     gf.zero()
-
     gf2 = gf.new()
     gf2.gaussian()
     assert gf.dot(gf2) == 0
@@ -134,7 +133,7 @@ def test_random(lib, lattice, device, dtype):
 
     gf2 = gf.copy()
     assert gf == gf2
-    assert isclose(gf.norm2(), (gf.field**2).sum(), rel_tol=1e-6)
+    assert isclose(gf.norm2(), (gf**2).sum(), rel_tol=1e-6)
 
 
 @dtype_loop  # enables dtype
@@ -147,7 +146,7 @@ def test_exponential(lib, lattice, device, dtype):
 
     gf.unity()
     mom.copy(out=gf)
-    assert np.allclose(gf.field, 0)
+    assert np.allclose(gf, 0)
     # gf.is_momentum = False
     assert gf == 0
 
@@ -161,7 +160,7 @@ def test_exponential(lib, lattice, device, dtype):
 
     mom.gaussian(epsilon=0)
     gf2 = mom.exponentiate()
-    assert np.allclose(gf.field, gf2.field)
+    assert np.allclose(gf, gf2)
     assert gf2 == gf
 
     gf.gaussian()
@@ -253,9 +252,9 @@ def test_paths_wins(lib, lattice, device):
     rparts = gf.plaquette_field(sum_paths=False, insertion=mom, left=False)
     lparts = tuple(zip(*lparts))
     rparts = tuple(zip(*rparts))
-    assert allclose(lparts[0], rparts[0])
-    assert allclose(lparts[1], out)
-    assert allclose(rparts[1], out)
+    assert all(np.allclose(left, right) for left, right in zip(lparts[0], rparts[0]))
+    assert all(np.allclose(left, out) for left in lparts[1])
+    assert all(np.allclose(left, out) for left in rparts[1])
 
     mom = momentum(lattice, dtype=dtype, device=device)
     mom.gaussian()
@@ -264,7 +263,7 @@ def test_paths_wins(lib, lattice, device):
     rparts = gf.plaquette_field(sum_paths=False, insertion=mom, left=False)
     lparts = tuple(zip(*lparts))
     rparts = tuple(zip(*rparts))
-    assert allclose(lparts[0], rparts[0])
+    assert all(np.allclose(left, right) for left, right in zip(lparts[0], rparts[0]))
     ltracs = tuple(loop.reduce() for loop in lparts[1])
     rtracs = tuple(loop.reduce() for loop in rparts[1])
     assert np.allclose(ltracs, rtracs)
diff --git a/test/test_solver.py b/test/test_solver.py
index 2fc4137..2aa01c3 100644
--- a/test/test_solver.py
+++ b/test/test_solver.py
@@ -20,6 +20,6 @@ def test_solve_random(lib, lattice, device, gamma, dtype=None):
     mat = dirac.M
     out = mat.solve(rhs, delta=1e-4)  # this value allowed convergence for all cases
     res = mat(out)
-    res.field -= rhs.field
+    res -= rhs
     res = res.norm() / rhs.norm()
     assert res < 1e-9
diff --git a/test/test_spinor.py b/test/test_spinor.py
index 2c7d3b1..9ad4a21 100644
--- a/test/test_spinor.py
+++ b/test/test_spinor.py
@@ -51,13 +51,13 @@ def test_params(lib, lattice, device, dtype):
 def test_init(lib, lattice, device, dtype):
     sf = spinor(lattice, dtype=dtype, device=device)
     sf.zero()
-    assert (sf.field == 0).all()
+    assert (sf == 0).all()
     assert sf.norm1() == 0
     assert sf.norm2() == 0
     sf.uniform()
-    assert np.isclose(sf.field.mean(), 0.5 + 0.5j, atol=0.1)
+    assert np.isclose(sf.mean(), 0.5 + 0.5j, atol=0.1)
     sf.gaussian()
-    assert np.isclose(sf.field.mean(), 0.0, atol=0.1)
+    assert np.isclose(sf.mean(), 0.0, atol=0.1)
 
 
 # @dtype_loop  # enables dtype
@@ -68,4 +68,4 @@ def test_gamma5(lib, lattice, device, gamma, dtype=None):
     sf = spinor(lattice, dtype=dtype, device=device, gamma_basis=gamma)
     sf.uniform()
     sf2 = sf.gamma5().apply_gamma5()
-    assert (sf.field == sf2.field).all()
+    assert (sf == sf2).all()