diff --git a/src/quaternion.hh b/src/quaternion.hh
index 048b6e063cfc3e166348f81b9abd846be29bd1e6..e68760e188e7735226e45adfa1148afa7fd4b7d5 100644
--- a/src/quaternion.hh
+++ b/src/quaternion.hh
@@ -260,66 +260,6 @@ public:
 
     }
 
-    void getFirstDerivativesOfDirectors(Dune::array<Dune::array<Dune::FieldVector<double,3>, 3>, 3>& dd_dvj) const
-    {
-        const Quaternion<T>& q = (*this);
-
-        // Contains \partial q / \partial v^i_j  at v = 0
-        Dune::array<Quaternion<double>,3> dq_dvj;
-        for (int j=0; j<3; j++) {
-            for (int m=0; m<4; m++)
-                dq_dvj[j][m]    = (j==m) * 0.5;
-        }
-        
-        // Contains \parder d \parder v_j
-        
-        for (int j=0; j<3; j++) {
-            
-            // d1
-            dd_dvj[0][j][0] = q[0]*(q.mult(dq_dvj[j]))[0] - q[1]*(q.mult(dq_dvj[j]))[1] 
-                - q[2]*(q.mult(dq_dvj[j]))[2] + q[3]*(q.mult(dq_dvj[j]))[3];
-            
-            dd_dvj[0][j][1] = (q.mult(dq_dvj[j]))[0]*q[1] + q[0]*(q.mult(dq_dvj[j]))[1]
-                + (q.mult(dq_dvj[j]))[2]*q[3] + q[2]*(q.mult(dq_dvj[j]))[3];
-            
-            dd_dvj[0][j][2] = (q.mult(dq_dvj[j]))[0]*q[2] + q[0]*(q.mult(dq_dvj[j]))[2]
-                - (q.mult(dq_dvj[j]))[1]*q[3] - q[1]*(q.mult(dq_dvj[j]))[3];
-            
-            // d2
-            dd_dvj[1][j][0] = (q.mult(dq_dvj[j]))[0]*q[1] + q[0]*(q.mult(dq_dvj[j]))[1]
-                - (q.mult(dq_dvj[j]))[2]*q[3] - q[2]*(q.mult(dq_dvj[j]))[3];
-            
-            dd_dvj[1][j][1] = - q[0]*(q.mult(dq_dvj[j]))[0] + q[1]*(q.mult(dq_dvj[j]))[1] 
-                - q[2]*(q.mult(dq_dvj[j]))[2] + q[3]*(q.mult(dq_dvj[j]))[3];
-            
-            dd_dvj[1][j][2] = (q.mult(dq_dvj[j]))[1]*q[2] + q[1]*(q.mult(dq_dvj[j]))[2]
-                + (q.mult(dq_dvj[j]))[0]*q[3] + q[0]*(q.mult(dq_dvj[j]))[3];
-            
-            // d3
-            dd_dvj[2][j][0] = (q.mult(dq_dvj[j]))[0]*q[2] + q[0]*(q.mult(dq_dvj[j]))[2]
-                + (q.mult(dq_dvj[j]))[1]*q[3] + q[1]*(q.mult(dq_dvj[j]))[3];
-            
-            dd_dvj[2][j][1] = (q.mult(dq_dvj[j]))[1]*q[2] + q[1]*(q.mult(dq_dvj[j]))[2]
-                - (q.mult(dq_dvj[j]))[0]*q[3] - q[0]*(q.mult(dq_dvj[j]))[3];
-            
-            dd_dvj[2][j][2] = - q[0]*(q.mult(dq_dvj[j]))[0] - q[1]*(q.mult(dq_dvj[j]))[1] 
-                + q[2]*(q.mult(dq_dvj[j]))[2] + q[3]*(q.mult(dq_dvj[j]))[3];
-            
-            
-            dd_dvj[0][j] *= 2;
-            dd_dvj[1][j] *= 2;
-            dd_dvj[2][j] *= 2;
-            
-        }
-    
-    // Check: The derivatives of the directors must be orthogonal to the directors
-    for (int i=0; i<3; i++)
-        for (int j=0; j<3; j++)
-            assert (std::abs(q.director(i) * dd_dvj[i][j]) < 1e-7);
-
-}
-
-
     /** \brief Turn quaternion into a unit quaternion by dividing by its Euclidean norm */
     void normalize() {
         (*this) /= this->two_norm();