diff --git a/src/algorithm/rnea-derivatives.hpp b/src/algorithm/rnea-derivatives.hpp
index 4c7ebb7822fc62c46a25cac5d8183d2e01903192..cea7fbf43a28d9f46d36739b401e8a2c6d7e2fb4 100644
--- a/src/algorithm/rnea-derivatives.hpp
+++ b/src/algorithm/rnea-derivatives.hpp
@@ -54,7 +54,7 @@ namespace se3
   /// \param[out] rnea_partial_dv Partial derivative of the generalized torque vector with respect to the joint velocity.
   /// \param[out] rnea_partial_da Partial derivative of the generalized torque vector with respect to the joint acceleration.
   ///
-  /// \remark rnea_partial_dq, rnea_partial_dv and rnea_partial_da must be first initialized with zeros (gravity_partial_dq.setZero).
+  /// \remark rnea_partial_dq, rnea_partial_dv and rnea_partial_da must be first initialized with zeros (rnea_partial_dq.setZero(),etc).
   ///         As for se3::crba, only the upper triangular part of rnea_partial_da is filled.
   ///
   /// \sa se3::rnea
@@ -78,11 +78,9 @@ namespace se3
   /// \param[in] v The joint velocity vector (dim model.nv).
   /// \param[in] a The joint acceleration vector (dim model.nv).
   ///
-  /// \returns The results are stored in data.dtau_dq, data.tau_dv and data.dtau_da which respectively correspond
+  /// \returns The results are stored in data.dtau_dq, data.M and data.dtau_da which respectively correspond
   ///          to the partial derivatives of the joint torque vector with respect to the joint configuration, velocity and acceleration.
-  ///          data.dtau_da is a reference on data.M. And as for se3::crba, only the upper triangular part of data.M is filled.
-  ///
-  /// \remark rnea_partial_dq and rnea_partial_dv must be first initialized with zeros (gravity_partial_dq.setZero).
+  ///          And as for se3::crba, only the upper triangular part of data.M is filled.
   ///
   /// \sa se3::rnea, se3::crba, se3::cholesky::decompose
   ///
@@ -93,7 +91,7 @@ namespace se3
                          const Eigen::VectorXd & a)
   {
     computeRNEADerivatives(model,data,q,v,a,
-                           data.dtau_dq, data.dtau_dv, data.dtau_da);
+                           data.dtau_dq, data.dtau_dv, data.M);
   }
 
 
diff --git a/src/multibody/model.hpp b/src/multibody/model.hpp
index 535cb478a77cf9c51ef90968614714a69694fe50..f1ce3eb2e1375faa5c91af6f570f07b534d83233 100644
--- a/src/multibody/model.hpp
+++ b/src/multibody/model.hpp
@@ -522,9 +522,6 @@ namespace se3
     /// \brief Variation of the joint torque vector with respect to the joint velocity.
     Eigen::MatrixXd dtau_dv;
     
-    /// \brief Variation of the joint torque vector with respect to the joint acceleration. This is only a reference on data.M.
-    Eigen::MatrixXd & dtau_da;
-    
     /// \brief Vector of joint placements wrt to algorithm end effector.
     container::aligned_vector<SE3> iMf;
 
diff --git a/src/multibody/model.hxx b/src/multibody/model.hxx
index 95b63a6126a794dd2aa48f0e4a19635d91b28ac9..5d1b697f3b23a8955a4cf834f32473d398942def 100644
--- a/src/multibody/model.hxx
+++ b/src/multibody/model.hxx
@@ -278,7 +278,6 @@ namespace se3
     ,dAdv(6,model.nv)
     ,dtau_dq(Eigen::MatrixXd::Zero(model.nv,model.nv))
     ,dtau_dv(Eigen::MatrixXd::Zero(model.nv,model.nv))
-    ,dtau_da(M)
     ,iMf((std::size_t)model.njoints)
     ,com((std::size_t)model.njoints)
     ,vcom((std::size_t)model.njoints)
diff --git a/unittest/rnea-derivatives.cpp b/unittest/rnea-derivatives.cpp
index 074eca2c800d37c8cf27eb4e3d2375906a3976e6..8d59632d223844dca74df5c0c04e26ff49ff24f3 100644
--- a/unittest/rnea-derivatives.cpp
+++ b/unittest/rnea-derivatives.cpp
@@ -256,7 +256,7 @@ BOOST_AUTO_TEST_CASE(test_rnea_derivatives)
 
   BOOST_CHECK(rnea_partial_dq.isApprox(data2.dtau_dq));
   BOOST_CHECK(rnea_partial_dv.isApprox(data2.dtau_dv));
-  BOOST_CHECK(rnea_partial_da.isApprox(data2.dtau_da));
+  BOOST_CHECK(rnea_partial_da.isApprox(data2.M));
   
 }