diff --git a/Z4co/src/constraint.cxx b/Z4co/src/constraint.cxx
index e3d3204c..d85662a8 100644
--- a/Z4co/src/constraint.cxx
+++ b/Z4co/src/constraint.cxx
@@ -9,7 +9,6 @@
 #endif
 
 #include "derivs.hxx"
-#include "physics.hxx"
 
 #include <loop_device.hxx>
 #include <simd.hxx>
diff --git a/Z4co/src/enforce.cxx b/Z4co/src/enforce.cxx
index 7542e015..428cd9bc 100644
--- a/Z4co/src/enforce.cxx
+++ b/Z4co/src/enforce.cxx
@@ -1,5 +1,3 @@
-#include "physics.hxx"
-
 #include <loop_device.hxx>
 #include <mat.hxx>
 #include <simd.hxx>
diff --git a/Z4co/src/initial1.cxx b/Z4co/src/initial1.cxx
index 7fd3d6af..c65262e0 100644
--- a/Z4co/src/initial1.cxx
+++ b/Z4co/src/initial1.cxx
@@ -1,5 +1,3 @@
-#include "physics.hxx"
-
 #include <loop_device.hxx>
 #include <mat.hxx>
 #include <simd.hxx>
diff --git a/Z4co/src/initial2.cxx b/Z4co/src/initial2.cxx
index fb87e2de..74a2fa12 100644
--- a/Z4co/src/initial2.cxx
+++ b/Z4co/src/initial2.cxx
@@ -1,5 +1,4 @@
 #include "derivs.hxx"
-#include "physics.hxx"
 
 #include <loop_device.hxx>
 #include <simd.hxx>
@@ -16,6 +15,29 @@ using namespace Arith;
 using namespace Loop;
 using namespace std;
 
+template <typename T, int D>
+ARITH_INLINE ARITH_DEVICE ARITH_HOST constexpr vec<smat<T, D>, D>
+calc_gammal(const smat<vec<T, D>, D> &dg) {
+  // Gammal_abc
+  return vec<smat<T, D>, D>([&](int a) ARITH_INLINE {
+    return smat<T, D>([&](int b, int c) ARITH_INLINE {
+      return (dg(a, b)(c) + dg(a, c)(b) - dg(b, c)(a)) / 2;
+    });
+  });
+}
+
+template <typename T, int D>
+ARITH_INLINE ARITH_DEVICE ARITH_HOST constexpr vec<smat<T, D>, D>
+calc_gamma(const smat<T, D> &gu, const vec<smat<T, D>, D> &Gammal) {
+  // Gamma^a_bc
+  return vec<smat<T, D>, D>([&](int a) ARITH_INLINE {
+    return smat<T, D>([&](int b, int c) ARITH_INLINE {
+      return sum<D>([&](int x)
+                        ARITH_INLINE { return gu(a, x) * Gammal(x)(b, c); });
+    });
+  });
+}
+
 extern "C" void Z4co_Initial2(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTS_Z4co_Initial2;
 
diff --git a/Z4co/src/physics.hxx b/Z4co/src/physics.hxx
deleted file mode 100644
index 9772b5c2..00000000
--- a/Z4co/src/physics.hxx
+++ /dev/null
@@ -1,99 +0,0 @@
-#ifndef PHYSICS_HXX
-#define PHYSICS_HXX
-
-#include <defs.hxx>
-#include <mat.hxx>
-#include <sum.hxx>
-#include <vec.hxx>
-
-namespace Z4co {
-using namespace Arith;
-
-template <typename T, int D>
-ARITH_INLINE ARITH_DEVICE ARITH_HOST constexpr smat<vec<T, D>, D>
-calc_dgu(const smat<T, D> &gu, const smat<vec<T, D>, D> &dg) {
-  // g_xy = g_xa g_yb g^ab
-  // g_xy,c = (g_xa g_yb g^ab),c
-  //        = g_xa,c g_yb g^ab + g_xa g_yb,c g^ab + g_xa g_yb g^ab,c
-  // g_xa g_yb g^ab,c = g_xy,c - g_xa,c g_yb g^ab - g_xa g_yb,c g^ab
-  //                  = g_xy,c - g_xy,c - g_xy,c
-  //                  = - g_xy,c
-  // g^ab,c = - g^ax g^by g_xy,c
-  return smat<vec<T, D>, D>([&](int a, int b) ARITH_INLINE {
-    return vec<T, D>([&](int c) ARITH_INLINE {
-      // return sum2sym([&](int x, int y) ARITH_INLINE {
-      //   return -gu(a, x) * gu(b, y) * dg(x, y)(c);
-      // });
-      return sum<D>([&](int x) ARITH_INLINE {
-        return -gu(a, x) * sum<D>([&](int y) ARITH_INLINE {
-          return gu(b, y) * dg(x, y)(c);
-        });
-      });
-    });
-  });
-}
-
-template <typename T, int D>
-ARITH_INLINE ARITH_DEVICE ARITH_HOST constexpr smat<vec<T, D>, D>
-calc_dAu(const smat<T, D> &gu, const smat<vec<T, D>, D> &dgu,
-         const smat<T, D> &A, const smat<vec<T, D>, D> &dA) {
-  // A^ab,c = (g^ax g^by A_xy),c
-  //        = g^ax,c g^by A_xy + g^ax g^by,c A_xy + g^ax g^by A_xy,c
-  return smat<vec<T, D>, D>([&](int a, int b) ARITH_INLINE {
-    return vec<T, D>([&](int c) ARITH_INLINE {
-      // return sum2sym([&](int x, int y) ARITH_INLINE {
-      //   return dgu(a, x)(c) * gu(b, y) * A(x, y)   //
-      //          + gu(a, x) * dgu(b, y)(c) * A(x, y) //
-      //          + gu(a, x) * gu(b, y) * dA(x, y)(c);
-      // });
-      return sum<D>([&](int x) ARITH_INLINE {
-        return gu(b, x) * sum<D>([&](int y) ARITH_INLINE {
-                 return dgu(a, y)(c) * A(x, y)   //
-                        + dgu(b, y)(c) * A(x, y) //
-                        + gu(b, y) * dA(x, y)(c);
-               });
-      });
-    });
-  });
-}
-
-template <typename T, int D>
-ARITH_INLINE ARITH_DEVICE ARITH_HOST constexpr vec<smat<T, D>, D>
-calc_gammal(const smat<vec<T, D>, D> &dg) {
-  // Gammal_abc
-  return vec<smat<T, D>, D>([&](int a) ARITH_INLINE {
-    return smat<T, D>([&](int b, int c) ARITH_INLINE {
-      return (dg(a, b)(c) + dg(a, c)(b) - dg(b, c)(a)) / 2;
-    });
-  });
-}
-
-template <typename T, int D>
-ARITH_INLINE ARITH_DEVICE ARITH_HOST constexpr vec<smat<T, D>, D>
-calc_gamma(const smat<T, D> &gu, const vec<smat<T, D>, D> &Gammal) {
-  // Gamma^a_bc
-  return vec<smat<T, D>, D>([&](int a) ARITH_INLINE {
-    return smat<T, D>([&](int b, int c) ARITH_INLINE {
-      return sum<D>([&](int x)
-                        ARITH_INLINE { return gu(a, x) * Gammal(x)(b, c); });
-    });
-  });
-}
-
-template <typename T, int D>
-ARITH_INLINE ARITH_DEVICE ARITH_HOST constexpr vec<vec<vec<T, D>, D>, D>
-calc_gammalu(const smat<T, D> &gu, const vec<smat<T, D>, D> &Gammal) {
-  // Gamma_ab^c
-  return vec<vec<vec<T, D>, D>, D>([&](int a) ARITH_INLINE {
-    return vec<vec<T, D>, D>([&](int b) ARITH_INLINE {
-      return vec<T, D>([&](int c) ARITH_INLINE {
-        return sum<D>([&](int x)
-                          ARITH_INLINE { return Gammal(a)(b, x) * gu(x, c); });
-      });
-    });
-  });
-}
-
-} // namespace Z4co
-
-#endif // #ifndef PHYSICS_HXX
diff --git a/Z4co/src/rhs.cxx b/Z4co/src/rhs.cxx
index 4a6b5fbc..df619023 100644
--- a/Z4co/src/rhs.cxx
+++ b/Z4co/src/rhs.cxx
@@ -11,7 +11,6 @@
 // #define Power(x, y) (Arith::pown((x), (y)))
 
 #include "derivs.hxx"
-#include "physics.hxx"
 
 #include <loop_device.hxx>
 #include <mat.hxx>