diff --git a/src/projections.F90 b/src/projections.F90
index 4432cfd..2f32632 100644
--- a/src/projections.F90
+++ b/src/projections.F90
@@ -476,8 +476,26 @@ subroutine point_volume(x0, tu, tv, tw, ku, kv, kw, coef, nctlu, nctlv, nctlw, n
         ! Calculate the Hessian
         do j = 1, 3
             do i = 1, 3
-                hessian(i, j) = dot_product(deriv(:, i), deriv(:, j)) + &
-                                dot_product(R, deriv2(:, i, j))
+                ! TODO clean up this routine
+                ! 3 ways of doing the hessian. Second term causes problems, so we either:
+                ! omit it, only add the diagonal, keep adding it. Needs further investigation
+                ! just not adding them at all seems to be the most robust...
+
+                hessian(i, j) = dot_product(deriv(:, i), deriv(:, j))
+
+                ! if (i .eq. j) then
+                !     hessian(i, j) = hessian(i, j) + dot_product(R, deriv2(:, i, j))
+                ! end if
+
+                ! hessian(i, j) = dot_product(deriv(:, i), deriv(:, j)) + &
+                !                 dot_product(R, deriv2(:, i, j))
+
+                ! The hessian gets terribly ill conditioned when the second derivatives are included.
+                ! In one case, this shows up as a 2x2 sub block getting almost zero determinant.
+                ! This is solving an optimization problem, so check if there are any features that
+                ! would cause the optimization problem to be non-convex!
+
+                ! write(*,*) i, j, dot_product(deriv(:, i), deriv(:, j)), dot_product(R, deriv2(:, i, j))
             end do
         end do