diff --git a/src/averagedistanceassembler.hh b/src/averagedistanceassembler.hh
index 17f9c3ee2667bf5d6988c12d3fdb016af6b7580c..0f62f9846c5913e56f35b8948a31ee86e2a1c642 100644
--- a/src/averagedistanceassembler.hh
+++ b/src/averagedistanceassembler.hh
@@ -33,9 +33,16 @@ public:
                           typename TargetSpace::EmbeddedTangentVector& gradient) const
     {
         gradient = 0;
+#if 0  // This old code uses the derivative of dist(,), which is (frequently, at least)
+        // not differentiable at 0
         for (size_t i=0; i<coefficients_.size(); i++)
             gradient.axpy(weights_[i]*TargetSpace::distance(coefficients_[i], x),
                           TargetSpace::derivativeOfDistanceWRTSecondArgument(coefficients_[i], x));
+#else
+        for (size_t i=0; i<coefficients_.size(); i++)
+            gradient.axpy(0.5*weights_[i], 
+                          TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
+#endif
     }
 
     void assembleMatrix(const TargetSpace& x,