diff --git a/Base/Spin/SpinMatrix.cpp b/Base/Spin/SpinMatrix.cpp
index 8553d1045e8ea3d15e0cfc10926ab0b8d4113024..2714fa233f53cd29591f4180e9795478b603fa90 100644
--- a/Base/Spin/SpinMatrix.cpp
+++ b/Base/Spin/SpinMatrix.cpp
@@ -38,6 +38,11 @@ SpinMatrix SpinMatrix::One()
     return Diag(1, 1);
 }
 
+SpinMatrix SpinMatrix::Null()
+{
+    return Diag(0, 0);
+}
+
 SpinMatrix SpinMatrix::FromBlochVector(const R3& v)
 {
     return {(1.0 + v.z()) / 2.0, complex_t(v.x(), -v.y()) / 2.0, complex_t(v.x(), v.y()) / 2.0,
diff --git a/Base/Spin/SpinMatrix.h b/Base/Spin/SpinMatrix.h
index caa89b5b8495302ce9f5dffab2ad937c608b9c7e..da1741b07ac83f1855a7a1d5495203bdb18b1edd 100644
--- a/Base/Spin/SpinMatrix.h
+++ b/Base/Spin/SpinMatrix.h
@@ -29,6 +29,7 @@ public:
 
     static SpinMatrix Diag(complex_t a_, complex_t d_);
     static SpinMatrix One();
+    static SpinMatrix Null();
 
     //! Constructs matrix (I+v*s)/2, where s is the Pauli vector.
     static SpinMatrix FromBlochVector(const R3& v);
diff --git a/Resample/Specular/ComputeFluxMagnetic.cpp b/Resample/Specular/ComputeFluxMagnetic.cpp
index e270c781dec2eff8fd3c68fde0e5a62f94e72b04..591d764006ec03f0aa35ee4df080bee1573b67e5 100644
--- a/Resample/Specular/ComputeFluxMagnetic.cpp
+++ b/Resample/Specular/ComputeFluxMagnetic.cpp
@@ -81,7 +81,7 @@ std::vector<MatrixFlux> computeTR(const SliceStack& slices, const std::vector<co
 
     if (N == 1) {
         TR[0].m_T = SpinMatrix::One();
-        TR[0].m_R = SpinMatrix();
+        TR[0].m_R = SpinMatrix::Null();
         return TR;
     }
 
@@ -89,8 +89,8 @@ std::vector<MatrixFlux> computeTR(const SliceStack& slices, const std::vector<co
         TR[0].m_T = SpinMatrix::One();
         TR[0].m_R = SpinMatrix::Diag(-1, -1);
         for (size_t i = 1; i < N; ++i) {
-            TR[i].m_T = SpinMatrix();
-            TR[i].m_R = SpinMatrix();
+            TR[i].m_T = SpinMatrix::Null();
+            TR[i].m_R = SpinMatrix::Null();
         }
         return TR;
     }
@@ -111,7 +111,7 @@ std::vector<MatrixFlux> computeTR(const SliceStack& slices, const std::vector<co
 
     // bottom boundary condition
     TR.back().m_T = SpinMatrix::One();
-    TR.back().m_R = SpinMatrix();
+    TR.back().m_R = SpinMatrix::Null();
 
     for (size_t i = N - 1; i > 0; --i) {
         double sigma = 0.;