diff --git a/CHANGELOG.md b/CHANGELOG.md
index fc312449fb..a904c4f3b1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 ## [Unreleased]
 
+* Added 6D loop contacts to account for contacts between parts of the robot in https://github.com/loco-3d/crocoddyl/pull/1309
 * Fixed the inequality constraints' feasibility computation by incorporating bounds into the calculation in https://github.com/loco-3d/crocoddyl/pull/1307
 * Improved the action factory used for unit testing in https://github.com/loco-3d/crocoddyl/pull/1300
 * Ignore ruff issues in ipython notebook files in https://github.com/loco-3d/crocoddyl/pull/1297
diff --git a/bindings/python/crocoddyl/multibody/contacts/contact-6d-loop.cpp b/bindings/python/crocoddyl/multibody/contacts/contact-6d-loop.cpp
new file mode 100644
index 0000000000..13693e4e3a
--- /dev/null
+++ b/bindings/python/crocoddyl/multibody/contacts/contact-6d-loop.cpp
@@ -0,0 +1,242 @@
+///////////////////////////////////////////////////////////////////////////////
+// BSD 3-Clause License
+//
+// Copyright (C) 2019-2024, LAAS-CNRS, University of Edinburgh
+//                          Heriot-Watt University
+// Copyright note valid unless otherwise stated in individual files.
+// All rights reserved.
+///////////////////////////////////////////////////////////////////////////////
+
+#include "crocoddyl/multibody/contacts/contact-6d-loop.hpp"
+
+#include "python/crocoddyl/multibody/multibody.hpp"
+#include "python/crocoddyl/utils/copyable.hpp"
+
+namespace crocoddyl {
+namespace python {
+
+void exposeContact6DLoop() {
+  bp::register_ptr_to_python<boost::shared_ptr<ContactModel6DLoop> >();
+
+  bp::class_<ContactModel6DLoop, bp::bases<ContactModelAbstract> >(
+      "ContactModel6DLoop",
+      "Rigid 6D contact model.\n\n"
+      "It defines a rigid 6D contact models based on acceleration-based "
+      "holonomic constraints.\n"
+      "The calc and calcDiff functions compute the contact Jacobian and drift "
+      "(holonomic constraint) or\n"
+      "the derivatives of the holonomic constraint, respectively.",
+      bp::init<boost::shared_ptr<StateMultibody>, int, pinocchio::SE3, int,
+               pinocchio::SE3, pinocchio::ReferenceFrame, std::size_t,
+               bp::optional<Eigen::Vector2d> >(
+          bp::args("self", "state", "joint1_id", "joint1_placement",
+                   "joint2_id", "joint2_placement", "type", "nu", "gains"),
+          "Initialize the contact model.\n\n"
+          ":param state: state of the multibody system\n"
+          ":param joint1_id: Parent joint id of the first contact frame\n"
+          ":param joint1_placement: Placement of the first contact frame with "
+          "respect to the parent joint\n"
+          ":param joint2_id: Parent joint id of the second contact frames\n"
+          ":param joint2_placement: Placement of the second contact frame with "
+          "respect to the parent joint\n"
+          ":param type: reference frame of contact (must be pinocchio::LOCAL)\n"
+          ":param nu: dimension of control vector\n"
+          ":param gains: gains of the contact model (default "
+          "np.matrix([0.,0.]))"))
+      .def("calc", &ContactModel6DLoop::calc, bp::args("self", "data", "x"),
+           "Compute the 6D loop-contact Jacobian and drift.\n\n"
+           "The rigid contact model throught acceleration-base holonomic "
+           "constraint\n"
+           "of the contact frame placement.\n"
+           ":param data: contact data\n"
+           ":param x: state point (dim. state.nx)")
+      .def("calcDiff", &ContactModel6DLoop::calcDiff,
+           bp::args("self", "data", "x"),
+           "Compute the derivatives of the 6D loop-contact holonomic "
+           "constraint.\n\n"
+           "The rigid contact model throught acceleration-base holonomic "
+           "constraint\n"
+           "of the contact frame placement.\n"
+           "It assumes that calc has been run first.\n"
+           ":param data: cost data\n"
+           ":param x: state point (dim. state.nx)")
+      .def("updateForce", &ContactModel6DLoop::updateForce,
+           bp::args("self", "data", "force"),
+           "Convert the Lagrangian into a stack of spatial forces.\n\n"
+           ":param data: cost data\n"
+           ":param force: force vector (dimension 6)")
+      .def("updateForceDiff", &ContactModel6DLoop::updateForceDiff,
+           bp::args("self", "data", "df_dx", "df_du"),
+           "Update the force derivatives.\n\n"
+           ":param data: cost data\n"
+           ":param df_dx: Jacobian of the force with respect to the state\n"
+           ":param df_du: Jacobian of the force with respect to the control")
+      .def("createData", &ContactModel6DLoop::createData,
+           bp::with_custodian_and_ward_postcall<0, 2>(),
+           bp::args("self", "data"),
+           "Create the 6D loop-contact data.\n\n"
+           "Each contact model has its own data that needs to be allocated. "
+           "This function\n"
+           "returns the allocated data for a predefined cost.\n"
+           ":param data: Pinocchio data\n"
+           ":return contact data.")
+      .add_property(
+          "joint1_id",
+          bp::make_function(&ContactModel6DLoop::get_joint1_id,
+                            bp::return_value_policy<bp::return_by_value>()),
+          "Parent joint id of the first contact frame")
+      .add_property(
+          "joint1_placement",
+          bp::make_function(&ContactModel6DLoop::get_joint1_placement,
+                            bp::return_value_policy<bp::return_by_value>()),
+          "Placement of the first contact frame with respect to the parent "
+          "joint")
+      .add_property(
+          "joint2_id",
+          bp::make_function(&ContactModel6DLoop::get_joint2_id,
+                            bp::return_value_policy<bp::return_by_value>()),
+          "Parent joint id of the second contact frame")
+      .add_property(
+          "joint2_placement",
+          bp::make_function(&ContactModel6DLoop::get_joint2_placement,
+                            bp::return_value_policy<bp::return_by_value>()),
+          "Placement of the second contact frame with respect to the parent "
+          "joint")
+      .add_property(
+          "gains",
+          bp::make_function(&ContactModel6DLoop::get_gains,
+                            bp::return_value_policy<bp::return_by_value>()),
+          "Baumegarte stabilisation gains (Kp, Kd)")
+      .def(CopyableVisitor<ContactModel6DLoop>());
+
+  bp::register_ptr_to_python<boost::shared_ptr<ContactData6DLoop> >();
+
+  bp::class_<ContactData6DLoop, bp::bases<ContactDataAbstract> >(
+      "ContactData6DLoop", "Data for 6DLoop contact.\n\n",
+      bp::init<ContactModel6DLoop*, pinocchio::Data*>(
+          bp::args("self", "model", "data"),
+          "Create 6D loop-contact data.\n\n"
+          ":param model: 6D loop-contact model\n"
+          ":param data: Pinocchio data")[bp::with_custodian_and_ward<
+          1, 2, bp::with_custodian_and_ward<1, 3> >()])
+      .add_property(
+          "v1_partial_dq",
+          bp::make_getter(&ContactData6DLoop::v1_partial_dq,
+                          bp::return_internal_reference<>()),
+          "Jacobian of the spatial velocity of the first contact frame wrt q")
+      .add_property("a1_partial_dq",
+                    bp::make_getter(&ContactData6DLoop::a1_partial_dq,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the spatial acceleration of the first contact "
+                    "frame wrt q")
+      .add_property("a1_partial_dv",
+                    bp::make_getter(&ContactData6DLoop::a1_partial_dv,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the spatial acceleration of the first contact "
+                    "frame wrt v")
+      .add_property(
+          "v2_partial_dq",
+          bp::make_getter(&ContactData6DLoop::v2_partial_dq,
+                          bp::return_internal_reference<>()),
+          "Jacobian of the spatial velocity of the second contact frame wrt q")
+      .add_property("a2_partial_dq",
+                    bp::make_getter(&ContactData6DLoop::a2_partial_dq,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the spatial acceleration of the second "
+                    "contact frame wrt q")
+      .add_property("a2_partial_dv",
+                    bp::make_getter(&ContactData6DLoop::a2_partial_dv,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the spatial acceleration of the second "
+                    "contact frame wrt v")
+      .add_property("da0_dx",
+                    bp::make_getter(&ContactData6DLoop::da0_dx,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the acceleration drift wrt x")
+      .add_property("da0_dq_t1",
+                    bp::make_getter(&ContactData6DLoop::da0_dq_t1,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the acceleration drift wrt q - part 1")
+      .add_property("da0_dq_t2",
+                    bp::make_getter(&ContactData6DLoop::da0_dq_t2,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the acceleration drift wrt q - part 2")
+      .add_property("da0_dq_t3",
+                    bp::make_getter(&ContactData6DLoop::da0_dq_t3,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the acceleration drift wrt q - part 3")
+      .add_property("f1Xj1",
+                    bp::make_getter(&ContactData6DLoop::f1Xj1,
+                                    bp::return_internal_reference<>()),
+                    "Inverse of the placement of the first contact frame in "
+                    "the joint frame - "
+                    "Action Matrix")
+      .add_property("f2Xj2",
+                    bp::make_getter(&ContactData6DLoop::f2Xj2,
+                                    bp::return_internal_reference<>()),
+                    "Inverse of the placement of the second contact frame in "
+                    "the joint frame "
+                    "- Action Matrix")
+      .add_property("f1Mf2",
+                    bp::make_getter(&ContactData6DLoop::f1Mf2,
+                                    bp::return_internal_reference<>()),
+                    "Relative placement of the contact frames")
+      .add_property("f1Xf2",
+                    bp::make_getter(&ContactData6DLoop::f1Xf2,
+                                    bp::return_internal_reference<>()),
+                    "Relative placement of the contact frames - Action Matrix")
+      .add_property("f1Jf1",
+                    bp::make_getter(&ContactData6DLoop::f1Jf1,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the first contact frame")
+      .add_property("f2Jf2",
+                    bp::make_getter(&ContactData6DLoop::f2Jf2,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the second frame in the joint frame")
+      .add_property("j1Jj1",
+                    bp::make_getter(&ContactData6DLoop::j1Jj1,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the first joint in the joint frame")
+      .add_property("j2Jj2",
+                    bp::make_getter(&ContactData6DLoop::j2Jj2,
+                                    bp::return_internal_reference<>()),
+                    "Jacobian of the second contact frame")
+      .add_property("f1vf1",
+                    bp::make_getter(&ContactData6DLoop::f1vf1,
+                                    bp::return_internal_reference<>()),
+                    "Velocity of the first contact frame")
+      .add_property("f2vf2",
+                    bp::make_getter(&ContactData6DLoop::f2vf2,
+                                    bp::return_internal_reference<>()),
+                    "Velocity of the second contact frame")
+      .add_property(
+          "f1vf2",
+          bp::make_getter(&ContactData6DLoop::f1vf2,
+                          bp::return_internal_reference<>()),
+          "Velocity of the second contact frame in the first contact frame")
+      .add_property("f1af1",
+                    bp::make_getter(&ContactData6DLoop::f1af1,
+                                    bp::return_internal_reference<>()),
+                    "Acceleration of the first contact frame")
+      .add_property("f2af2",
+                    bp::make_getter(&ContactData6DLoop::f2af2,
+                                    bp::return_internal_reference<>()),
+                    "Acceleration of the second contact frame")
+      .add_property(
+          "f1af2",
+          bp::make_getter(&ContactData6DLoop::f1af2,
+                          bp::return_internal_reference<>()),
+          "Acceleration of the second contact frame in the first contact frame")
+      .add_property("joint1_f",
+                    bp::make_getter(&ContactData6DLoop::joint1_f,
+                                    bp::return_internal_reference<>()),
+                    "Force at the first joint")
+      .add_property("joint2_f",
+                    bp::make_getter(&ContactData6DLoop::joint2_f,
+                                    bp::return_internal_reference<>()),
+                    "Force at the second joint")
+      .def(CopyableVisitor<ContactData6DLoop>());
+}
+
+}  // namespace python
+}  // namespace crocoddyl
diff --git a/bindings/python/crocoddyl/multibody/multibody.cpp b/bindings/python/crocoddyl/multibody/multibody.cpp
index 57dcfca1c7..96cf782e06 100644
--- a/bindings/python/crocoddyl/multibody/multibody.cpp
+++ b/bindings/python/crocoddyl/multibody/multibody.cpp
@@ -56,6 +56,7 @@ void exposeMultibody() {
   exposeContact2D();
   exposeContact3D();
   exposeContact6D();
+  exposeContact6DLoop();
   exposeImpulse3D();
   exposeImpulse6D();
 }
diff --git a/bindings/python/crocoddyl/multibody/multibody.hpp b/bindings/python/crocoddyl/multibody/multibody.hpp
index 4eb73b0a1b..6e1f150092 100644
--- a/bindings/python/crocoddyl/multibody/multibody.hpp
+++ b/bindings/python/crocoddyl/multibody/multibody.hpp
@@ -61,6 +61,7 @@ void exposeContact1D();
 void exposeContact2D();
 void exposeContact3D();
 void exposeContact6D();
+void exposeContact6DLoop();
 void exposeImpulse3D();
 void exposeImpulse6D();
 void exposeMultibody();
diff --git a/include/crocoddyl/multibody/contact-base.hpp b/include/crocoddyl/multibody/contact-base.hpp
index 1b705f2cca..846548b15a 100644
--- a/include/crocoddyl/multibody/contact-base.hpp
+++ b/include/crocoddyl/multibody/contact-base.hpp
@@ -94,8 +94,9 @@ class ContactModelAbstractTpl {
    * @param[in] data   Contact data
    * @param[in] force  Contact force
    */
-  void updateForceDiff(const boost::shared_ptr<ContactDataAbstract>& data,
-                       const MatrixXs& df_dx, const MatrixXs& df_du) const;
+  virtual void updateForceDiff(
+      const boost::shared_ptr<ContactDataAbstract>& data, const MatrixXs& df_dx,
+      const MatrixXs& df_du) const;
 
   /**
    * @brief Set the stack of spatial forces to zero
diff --git a/include/crocoddyl/multibody/contacts/contact-6d-loop.hpp b/include/crocoddyl/multibody/contacts/contact-6d-loop.hpp
new file mode 100644
index 0000000000..457b66049c
--- /dev/null
+++ b/include/crocoddyl/multibody/contacts/contact-6d-loop.hpp
@@ -0,0 +1,365 @@
+///////////////////////////////////////////////////////////////////////////////
+// BSD 3-Clause License
+//
+// Copyright (C) 2019-2024, LAAS-CNRS, University of Edinburgh,
+//                          Heriot-Watt University
+// Copyright note valid unless otherwise stated in individual files.
+// All rights reserved.
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef CROCODDYL_MULTIBODY_CONTACTS_CONTACT_6D_LOOP_HPP_
+#define CROCODDYL_MULTIBODY_CONTACTS_CONTACT_6D_LOOP_HPP_
+
+#include <pinocchio/multibody/data.hpp>
+#include <pinocchio/spatial/motion.hpp>
+
+#include "crocoddyl/core/utils/deprecate.hpp"
+#include "crocoddyl/multibody/contact-base.hpp"
+#include "crocoddyl/multibody/fwd.hpp"
+
+namespace crocoddyl {
+
+template <typename _Scalar>
+class ContactModel6DLoopTpl : public ContactModelAbstractTpl<_Scalar> {
+ public:
+  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
+  typedef _Scalar Scalar;
+  typedef MathBaseTpl<Scalar> MathBase;
+  typedef ContactModelAbstractTpl<Scalar> Base;
+  typedef ContactData6DLoopTpl<Scalar> Data;
+  typedef StateMultibodyTpl<Scalar> StateMultibody;
+  typedef ContactDataAbstractTpl<Scalar> ContactDataAbstract;
+  typedef pinocchio::SE3Tpl<Scalar> SE3;
+  typedef pinocchio::ForceTpl<Scalar> Force;
+  typedef typename MathBase::Vector2s Vector2s;
+  typedef typename MathBase::Vector3s Vector3s;
+  typedef typename MathBase::VectorXs VectorXs;
+  typedef typename MathBase::Matrix3s Matrix3s;
+  typedef typename MathBase::Matrix6s Matrix6s;
+  typedef typename MathBase::MatrixXs MatrixXs;
+
+  /**
+   * @brief Initialize the 6d loop-contact model from joint and placements
+   *
+   *
+   * @param[in] state             State of the multibody system
+   * @param[in] joint1_id         Parent joint id of the first contact
+   * @param[in] joint1_placement  Placement of the first contact with
+   *                                  respect to the parent joint
+   * @param[in] joint2_id         Parent joint id of the second contact
+   * @param[in] joint2_placement  Placement of the second contact with
+   *                                  respect to the parent joint
+   * @param[in] type              Reference frame of contact
+   * @param[in] nu                Dimension of the control vector
+   * @param[in] gains             Baumgarte stabilization gains
+   */
+  ContactModel6DLoopTpl(
+      boost::shared_ptr<StateMultibody> state, const int joint1_id,
+      const SE3 &joint1_placement, const int joint2_id,
+      const SE3 &joint2_placement,
+      const pinocchio::ReferenceFrame type const std::size_t nu,
+      const Vector2s &gains = Vector2s::Zero());
+
+  /**
+   * @brief Initialize the 6d loop-contact model from joint and placements
+   *
+   *
+   * @param[in] state             State of the multibody system
+   * @param[in] joint1_id         Parent joint id of the first contact
+   * @param[in] joint1_placement  Placement of the first contact with
+   *                                  respect to the parent joint
+   * @param[in] joint2_id         Parent joint id of the second contact
+   * @param[in] joint2_placement  Placement of the second contact with
+   *                                  respect to the parent joint
+   * @param[in] type              Reference frame of contact
+   * @param[in] gains             Baumgarte stabilization gains
+   */
+  ContactModel6DLoopTpl(boost::shared_ptr<StateMultibody> state,
+                        const int joint1_id, const SE3 &joint1_placement,
+                        const int joint2_id, const SE3 &joint2_placement,
+                        const pinocchio::ReferenceFrame type const Vector2s
+                            &gains = Vector2s::Zero());
+
+  virtual ~ContactModel6DLoopTpl();
+
+  /**
+   * @brief Compute the 6d loop-contact Jacobian and drift
+   *
+   * @param[in] data  6d loop-contact data
+   * @param[in] x     State point \f$\mathbf{x}\in\mathbb{R}^{ndx}\f$
+   * @param[in] u     Control input \f$\mathbf{u}\in\mathbb{R}^{nu}\f$
+   */
+  virtual void calc(const boost::shared_ptr<ContactDataAbstract> &data,
+                    const Eigen::Ref<const VectorXs> &x);
+
+  /**
+   * @brief Compute the derivatives of the 6d loop-contact holonomic constraint
+   *
+   * @param[in] data  6d loop-contact data
+   * @param[in] x     State point \f$\mathbf{x}\in\mathbb{R}^{ndx}\f$
+   * @param[in] u     Control input \f$\mathbf{u}\in\mathbb{R}^{nu}\f$
+   */
+  virtual void calcDiff(const boost::shared_ptr<ContactDataAbstract> &data,
+                        const Eigen::Ref<const VectorXs> &x);
+
+  /**
+   * @brief Convert the force into a stack of spatial forces
+   *
+   * @param[in] data   6d loop-contact data
+   * @param[in] force  6d force
+   */
+  virtual void updateForce(const boost::shared_ptr<ContactDataAbstract> &data,
+                           const VectorXs &force);
+
+  /**
+   * @brief Updates the force differential for the given contact data.
+   *
+   * This function updates the force differential matrices with respect to the
+   * state and control variables.
+   *
+   * @param[in] data  Shared pointer to the contact data abstract.
+   * @param[in] df_dx Matrix representing the differential of the force with
+   * respect to the state variables.
+   * @param[in] df_du Matrix representing the differential of the force with
+   * respect to the control variables.
+   */
+  virtual void updateForceDiff(
+      const boost::shared_ptr<ContactDataAbstract> &data, const MatrixXs &df_dx,
+      const MatrixXs &df_du);
+
+  /**
+   * @brief Create the 6d loop-contact data
+   */
+  virtual boost::shared_ptr<ContactDataAbstract> createData(
+      pinocchio::DataTpl<Scalar> *const data);
+
+  /**
+   * @brief Return the first contact frame parent joint
+   */
+  const int get_joint1_id() const;
+
+  /**
+   * @brief Return the first contact frame placement with respect to the parent
+   * joint
+   */
+  const SE3 &get_joint1_placement() const;
+
+  /**
+   * @brief Return the second contact frame parent joint
+   */
+  const int get_joint2_id() const;
+
+  /**
+   * @brief Return the second contact frame placement with respect to the parent
+   * joint
+   */
+  const SE3 &get_joint2_placement() const;
+
+  /**
+   * @brief Return the Baumgarte stabilization gains
+   */
+  const Vector2s &get_gains() const;
+
+  /**
+   * @brief Set the first contact frame parent joint
+   */
+  void set_joint1_id(const int joint1_id);
+
+  /**
+   * @brief Set the first contact frame placement with respect to the parent
+   * joint
+   */
+  void set_joint1_placement(const SE3 &joint1_placement);
+
+  /**
+   * @brief Set the second contact frame parent joint
+   */
+  void set_joint2_id(const int joint2_id);
+
+  /**
+   * @brief Set the second contact frame placement with respect to the parent
+   * joint
+   */
+  void set_joint2_placement(const SE3 &joint2_placement);
+
+  /**
+   * @brief Set the Baumgarte stabilization gains
+   */
+  void set_gains(const Vector2s &gains);
+
+  /**
+   * @brief Print relevant information of the 6D loop-contact model
+   *
+   * @param[out] os  Output stream object
+   */
+  virtual void print(std::ostream &os) const;
+
+ protected:
+  using Base::id_;
+  using Base::nc_;
+  using Base::nu_;
+  using Base::state_;
+  using Base::type_;
+
+ private:
+  int joint1_id_;         //!< Parent joint id of the first contact
+  SE3 joint1_placement_;  //!< Placement of the first contact with respect to
+                          //!< the parent joint
+  int joint2_id_;         //!< Parent joint id of the second contact
+  SE3 joint2_placement_;  //!< Placement of the second contact with respect to
+                          //!< the parent joint
+  Vector2s gains_;        //!< Baumgarte stabilization gains
+};
+
+template <typename _Scalar>
+struct ContactData6DLoopTpl : public ContactDataAbstractTpl<_Scalar> {
+  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
+  typedef _Scalar Scalar;
+  typedef MathBaseTpl<Scalar> MathBase;
+  typedef ContactDataAbstractTpl<Scalar> Base;
+  typedef typename MathBase::Matrix3s Matrix3s;
+  typedef typename MathBase::Matrix6xs Matrix6xs;
+  typedef typename MathBase::Matrix6s Matrix6s;
+  typedef typename MathBase::MatrixXs MatrixXs;
+  typedef typename pinocchio::SE3Tpl<Scalar> SE3;
+  typedef typename pinocchio::SE3Tpl<Scalar>::ActionMatrixType SE3ActionMatrix;
+  typedef typename pinocchio::MotionTpl<Scalar> Motion;
+  typedef typename pinocchio::ForceTpl<Scalar> Force;
+
+  template <template <typename Scalar> class Model>
+  ContactData6DLoopTpl(Model<Scalar> *const model,
+                       pinocchio::DataTpl<Scalar> *const data)
+      : Base(model, data),
+        v1_partial_dq(6, model->get_state()->get_nv()),
+        f1_v1_partial_dq(6, model->get_state()->get_nv()),
+        a1_partial_dq(6, model->get_state()->get_nv()),
+        a1_partial_dv(6, model->get_state()->get_nv()),
+        v2_partial_dq(6, model->get_state()->get_nv()),
+        f2_v2_partial_dq(6, model->get_state()->get_nv()),
+        f1_v2_partial_dq(6, model->get_state()->get_nv()),
+        a2_partial_dq(6, model->get_state()->get_nv()),
+        f2_a2_partial_dq(6, model->get_state()->get_nv()),
+        f2_a2_partial_dv(6, model->get_state()->get_nv()),
+        a2_partial_dv(6, model->get_state()->get_nv()),
+        __partial_da(6, model->get_state()->get_nv()),
+        da0_dq_t1(6, model->get_state()->get_nv()),
+        da0_dq_t2(6, model->get_state()->get_nv()),
+        da0_dq_t2_tmp(6, model->get_state()->get_nv()),
+        da0_dq_t3(6, model->get_state()->get_nv()),
+        da0_dq_t3_tmp(6, model->get_state()->get_nv()),
+        dpos_dq(6, model->get_state()->get_nv()),
+        dvel_dq(6, model->get_state()->get_nv()),
+        dtau_dq_tmp(model->get_state()->get_nv(), model->get_state()->get_nv()),
+        f1Jf1(6, model->get_state()->get_nv()),
+        f2Jf2(6, model->get_state()->get_nv()),
+        f1Jf2(6, model->get_state()->get_nv()),
+        j1Jj1(6, model->get_state()->get_nv()),
+        j2Jj2(6, model->get_state()->get_nv()),
+        f1Xj1(model->get_joint1_placement().inverse()),
+        f2Xj2(model->get_joint2_placement().inverse()),
+        f1Mf2(SE3::Identity()),
+        f1Xf2(SE3ActionMatrix::Identity()),
+        f1vf1(Motion::Zero()),
+        f2vf2(Motion::Zero()),
+        f1vf2(Motion::Zero()),
+        f1af1(Motion::Zero()),
+        f2af2(Motion::Zero()),
+        f1af2(Motion::Zero()),
+        joint1_f(Force::Zero()),
+        joint2_f(Force::Zero()) {
+    v1_partial_dq.setZero();
+    f1_v1_partial_dq.setZero();
+    a1_partial_dq.setZero();
+    a1_partial_dv.setZero();
+    v2_partial_dq.setZero();
+    f2_v2_partial_dq.setZero();
+    f1_v2_partial_dq.setZero();
+    a2_partial_dq.setZero();
+    f2_a2_partial_dq.setZero();
+    f2_a2_partial_dv.setZero();
+    a2_partial_dv.setZero();
+    __partial_da.setZero();
+    da0_dq_t1.setZero();
+    da0_dq_t2.setZero();
+    da0_dq_t2_tmp.setZero();
+    da0_dq_t3.setZero();
+    da0_dq_t3_tmp.setZero();
+    dpos_dq.setZero();
+    dvel_dq.setZero();
+    dtau_dq_tmp.setZero();
+    f1Jf1.setZero();
+    f2Jf2.setZero();
+    f1Jf2.setZero();
+    j1Jj1.setZero();
+    j2Jj2.setZero();
+    j2Jj1.setZero();
+  }
+
+  using Base::a0;
+  using Base::da0_dx;
+  using Base::df_du;
+  using Base::df_dx;
+  using Base::dtau_dq;
+  using Base::f;
+  using Base::Jc;
+  using Base::pinocchio;
+
+  Matrix6xs v1_partial_dq;
+  Matrix6xs f1_v1_partial_dq;
+  Matrix6xs a1_partial_dq;
+  Matrix6xs a1_partial_dv;
+  Matrix6xs v2_partial_dq;
+  Matrix6xs f2_v2_partial_dq;
+  Matrix6xs f1_v2_partial_dq;
+  Matrix6xs a2_partial_dq;
+  Matrix6xs f2_a2_partial_dq;
+  Matrix6xs f2_a2_partial_dv;
+  Matrix6xs a2_partial_dv;
+  Matrix6xs __partial_da;
+
+  Matrix6xs da0_dq_t1;
+  Matrix6xs da0_dq_t2;
+  Matrix6xs da0_dq_t2_tmp;
+  Matrix6xs da0_dq_t3;
+  Matrix6xs da0_dq_t3_tmp;
+
+  Matrix6xs dpos_dq;
+  Matrix6xs dvel_dq;
+  MatrixXs dtau_dq_tmp;
+  // Placement related data
+  SE3 oMf1;   // Placement of the first contact frame in the world frame
+  SE3 oMf2;   // Placement of the second contact frame in the world frame
+  SE3 f1Mf2;  // Relative placement of the contact frames in the first contact
+              // frame
+  SE3ActionMatrix f1Xj1;
+  SE3ActionMatrix f2Xj2;
+  SE3ActionMatrix f1Xf2;
+  // Jacobian related data
+  Matrix6xs f1Jf1;
+  Matrix6xs f2Jf2;
+  Matrix6xs f1Jf2;
+  Matrix6xs j1Jj1;
+  Matrix6xs j2Jj2;
+  Matrix6xs j2Jj1;
+  // Velocity related data
+  Motion f1vf1;
+  Motion f2vf2;
+  Motion f1vf2;
+  // Acceleration related data
+  Motion f1af1;
+  Motion f2af2;
+  Motion f1af2;
+  // Force related data
+  Force joint1_f;
+  Force joint2_f;
+};
+
+}  // namespace crocoddyl
+/* --- Details -------------------------------------------------------------- */
+/* --- Details -------------------------------------------------------------- */
+/* --- Details -------------------------------------------------------------- */
+#include "crocoddyl/multibody/contacts/contact-6d-loop.hxx"
+
+#endif  // CROCODDYL_MULTIBODY_CONTACTS_CONTACT_6D_LOOP_HPP_
diff --git a/include/crocoddyl/multibody/contacts/contact-6d-loop.hxx b/include/crocoddyl/multibody/contacts/contact-6d-loop.hxx
new file mode 100644
index 0000000000..860f09ade4
--- /dev/null
+++ b/include/crocoddyl/multibody/contacts/contact-6d-loop.hxx
@@ -0,0 +1,294 @@
+
+///////////////////////////////////////////////////////////////////////////////
+// BSD 3-Clause License
+//
+// Copyright (C) 2019-2024, LAAS-CNRS, University of Edinburgh,
+//                          Heriot-Watt University
+// Copyright note valid unless otherwise stated in individual files.
+// All rights reserved.
+///////////////////////////////////////////////////////////////////////////////
+
+#include <pinocchio/algorithm/frames.hpp>
+#include <pinocchio/algorithm/kinematics-derivatives.hpp>
+
+#include "crocoddyl/core/utils/exception.hpp"
+#include "crocoddyl/multibody/contacts/contact-6d.hpp"
+
+namespace crocoddyl {
+
+template <typename Scalar>
+ContactModel6DLoopTpl<Scalar>::ContactModel6DLoopTpl(
+    boost::shared_ptr<StateMultibody> state, const int joint1_id,
+    const SE3 &joint1_placement, const int joint2_id,
+    const SE3 &joint2_placement,
+    const pinocchio::ReferenceFrame type const std::size_t nu,
+    const Vector2s &gains)
+    : Base(state, pinocchio::ReferenceFrame::LOCAL, 6, nu),
+      joint1_id_(joint1_id),
+      joint2_id_(joint2_id),
+      joint1_placement_(joint1_placement),
+      joint2_placement_(joint2_placement),
+      gains_(gains) {
+  if (type != pinocchio::ReferenceFrame::LOCAL) {
+    std::cerr << "Warning: Only reference frame pinocchio::LOCAL is supported "
+                 "for 6D loop contacts"
+              << std::endl;
+  }
+  if (joint1_id == 0 || joint2_id == 0) {
+    std::cerr << "Warning: At least one of the parents joints id is zero"
+                 "you should use crocoddyl::ContactModel6D instead"
+              << std::endl;
+  }
+}
+
+template <typename Scalar>
+ContactModel6DLoopTpl<Scalar>::ContactModel6DLoopTpl(
+    boost::shared_ptr<StateMultibody> state, const int joint1_id,
+    const SE3 &joint1_placement, const int joint2_id,
+    const SE3 &joint2_placement,
+    const pinocchio::ReferenceFrame type const Vector2s &gains)
+    : Base(state, pinocchio::ReferenceFrame::LOCAL, 6),
+      joint1_id_(joint1_id),
+      joint2_id_(joint2_id),
+      joint1_placement_(joint1_placement),
+      joint2_placement_(joint2_placement),
+      gains_(gains) {
+  if (type != pinocchio::ReferenceFrame::LOCAL) {
+    std::cerr << "Warning: Only reference frame pinocchio::LOCAL is supported "
+                 "for 6D loop contacts"
+              << std::endl;
+  }
+  if (joint1_id == 0 || joint2_id == 0) {
+    std::cerr << "Warning: At least one of the parents joints id is zero"
+                 "you should use crocoddyl::ContactModel6D instead"
+              << std::endl;
+  }
+}
+
+template <typename Scalar>
+ContactModel6DLoopTpl<Scalar>::~ContactModel6DLoopTpl() {}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::calc(
+    const boost::shared_ptr<ContactDataAbstract> &data,
+    const Eigen::Ref<const VectorXs> &) {
+  Data *d = static_cast<Data *>(data.get());
+  pinocchio::getJointJacobian(*state_->get_pinocchio().get(), *d->pinocchio,
+                              joint1_id_, pinocchio::LOCAL, d->j1Jj1);
+  pinocchio::getJointJacobian(*state_->get_pinocchio().get(), *d->pinocchio,
+                              joint2_id_, pinocchio::LOCAL, d->j2Jj2);
+  d->f1Jf1.noalias() = d->f1Xj1 * d->j1Jj1;
+  d->f2Jf2.noalias() = d->f2Xj2 * d->j2Jj2;
+
+  d->oMf1 = d->pinocchio->oMi[joint1_id_].act(joint1_placement_);
+  d->oMf2 = d->pinocchio->oMi[joint2_id_].act(joint2_placement_);
+  d->f1Mf2 = d->oMf1.actInv(d->oMf2);
+  d->f1Xf2.noalias() = d->f1Mf2.toActionMatrix();
+
+  d->Jc.noalias() = d->f1Jf1 - d->f1Xf2 * d->f2Jf2;
+  // Compute the acceleration drift
+  if (joint1_id_ > 0) {
+    d->f1vf1 = joint1_placement_.actInv(d->pinocchio->v[joint1_id_]);
+    d->f1af1 = joint1_placement_.actInv(d->pinocchio->a[joint1_id_]);
+  }
+  if (joint2_id_ > 0) {
+    d->f2vf2 = joint2_placement_.actInv(d->pinocchio->v[joint2_id_]);
+    d->f2af2 = joint2_placement_.actInv(d->pinocchio->a[joint2_id_]);
+    d->f1vf2 = d->f1Mf2.act(d->f2vf2);
+    d->f1af2 = d->f1Mf2.act(d->f2af2);
+  }
+  d->a0.noalias() = (d->f1af1 - d->f1af2 + d->f1vf1.cross(d->f1vf2)).toVector();
+
+  if (std::abs<Scalar>(gains_[0]) > std::numeric_limits<Scalar>::epsilon()) {
+    d->a0.noalias() -= gains_[0] * pinocchio::log6(d->f1Mf2).toVector();
+  }
+  if (std::abs<Scalar>(gains_[1]) > std::numeric_limits<Scalar>::epsilon()) {
+    d->a0.noalias() += gains_[1] * (d->f1vf1 - d->f1vf2).toVector();
+  }
+}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::calcDiff(
+    const boost::shared_ptr<ContactDataAbstract> &data,
+    const Eigen::Ref<const VectorXs> &) {
+  Data *d = static_cast<Data *>(data.get());
+  const std::size_t nv = state_->get_nv();
+  if (joint1_id_ > 0) {
+    d->f1af1 = joint1_placement_.actInv(d->pinocchio->a[joint1_id_]);
+  }
+  if (joint2_id_ > 0) {
+    d->f2af2 = joint2_placement_.actInv(d->pinocchio->a[joint2_id_]);
+    d->f1af2 = d->f1Mf2.act(d->f2af2);
+  }
+
+  pinocchio::getJointAccelerationDerivatives(
+      *state_->get_pinocchio().get(), *d->pinocchio, joint1_id_,
+      pinocchio::LOCAL, d->v1_partial_dq, d->a1_partial_dq, d->a1_partial_dv,
+      d->__partial_da);
+  pinocchio::getJointAccelerationDerivatives(
+      *state_->get_pinocchio().get(), *d->pinocchio, joint2_id_,
+      pinocchio::LOCAL, d->v2_partial_dq, d->a2_partial_dq, d->a2_partial_dv,
+      d->__partial_da);
+
+  d->da0_dq_t1.noalias() =
+      joint1_placement_.toActionMatrixInverse() * d->a1_partial_dq;
+
+  d->da0_dq_t2.noalias() = d->f1af2.toActionMatrix() * d->Jc;
+  d->f2_a2_partial_dq.noalias() =
+      joint2_placement_.toActionMatrixInverse() * d->a2_partial_dq;
+  d->da0_dq_t2.noalias() += d->f1Xf2 * d->f2_a2_partial_dq;
+
+  d->f1_v1_partial_dq.noalias() =
+      joint1_placement_.toActionMatrixInverse() * d->v1_partial_dq;
+  d->da0_dq_t3.noalias() = -d->f1vf2.toActionMatrix() * d->f1_v1_partial_dq;
+  d->da0_dq_t3_tmp.noalias() = d->f1vf2.toActionMatrix() * d->Jc;
+  d->da0_dq_t3.noalias() += d->f1vf1.toActionMatrix() * d->da0_dq_t3_tmp;
+  d->da0_dq_t3_tmp.noalias() =
+      joint2_placement_.toActionMatrixInverse() * d->v2_partial_dq;
+  d->da0_dq_t3.noalias() +=
+      d->f1vf1.toActionMatrix() * d->f1Xf2 * d->da0_dq_t3_tmp;
+  d->da0_dx.leftCols(nv).noalias() = d->da0_dq_t1 - d->da0_dq_t2 + d->da0_dq_t3;
+
+  d->f2_a2_partial_dv.noalias() =
+      joint2_placement_.toActionMatrixInverse() * d->a2_partial_dv;
+  d->f1Jf2.noalias() = d->f1Xf2 * d->f2Jf2;
+  d->da0_dx.rightCols(nv).noalias() =
+      joint1_placement_.toActionMatrixInverse() * d->a1_partial_dv;
+  d->da0_dx.rightCols(nv).noalias() -= d->f1Xf2 * d->f2_a2_partial_dv;
+  d->da0_dx.rightCols(nv).noalias() -= d->f1vf2.toActionMatrix() * d->f1Jf1;
+  d->da0_dx.rightCols(nv).noalias() += d->f1vf1.toActionMatrix() * d->f1Jf2;
+
+  if (std::abs<Scalar>(gains_[0]) > std::numeric_limits<Scalar>::epsilon()) {
+    Matrix6s f1Mf2_log6;
+    pinocchio::Jlog6(d->f1Mf2, f1Mf2_log6);
+    d->dpos_dq.noalias() =
+        d->oMf2.toActionMatrixInverse() * d->oMf1.toActionMatrix() * d->f1Jf1;
+    d->dpos_dq.noalias() -= d->f2Jf2;
+    d->da0_dx.leftCols(nv).noalias() += gains_[0] * f1Mf2_log6 * d->dpos_dq;
+  }
+  if (std::abs<Scalar>(gains_[1]) > std::numeric_limits<Scalar>::epsilon()) {
+    d->f2_v2_partial_dq.noalias() =
+        joint2_placement_.toActionMatrixInverse() * d->v2_partial_dq;
+    d->f1_v2_partial_dq.noalias() = d->f1Xf2 * d->f2_v2_partial_dq;
+    d->dvel_dq.noalias() = d->f1_v1_partial_dq - d->f1_v2_partial_dq;
+    d->dvel_dq.noalias() -= d->f1vf2.toActionMatrix() * d->Jc;
+    d->da0_dx.leftCols(nv).noalias() += gains_[1] * d->dvel_dq;
+    d->da0_dx.rightCols(nv).noalias() += gains_[1] * d->Jc;
+  }
+}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::updateForce(
+    const boost::shared_ptr<ContactDataAbstract> &data, const VectorXs &force) {
+  if (force.size() != 6) {
+    throw_pretty("Contact force vector has wrong dimension (expected 6 got "
+                 << force.size() << ")");
+  }
+  Data *d = static_cast<Data *>(data.get());
+  d->f = pinocchio::ForceTpl<Scalar>(-force);
+  d->fext = joint1_placement_.act(d->f);
+  d->joint1_f = -joint1_placement_.act(d->f);
+  d->joint2_f = (joint2_placement_ * d->f1Mf2.inverse()).act(d->f);
+
+  Matrix6s f_cross = Matrix6s::Zero(6, 6);
+  f_cross.template topRightCorner<3, 3>() =
+      pinocchio::skew(d->joint2_f.linear());
+  f_cross.template bottomLeftCorner<3, 3>() =
+      pinocchio::skew(d->joint2_f.linear());
+  f_cross.template bottomRightCorner<3, 3>() =
+      pinocchio::skew(d->joint2_f.angular());
+
+  SE3 j2Mj1 =
+      joint2_placement_.act(d->f1Mf2.actInv(joint1_placement_.inverse()));
+  d->j2Jj1.noalias() = j2Mj1.toActionMatrix() * d->j1Jj1;
+  d->dtau_dq_tmp.noalias() = -f_cross * (d->j2Jj2 - d->j2Jj1);
+  d->dtau_dq.noalias() = d->j2Jj2.transpose() * d->dtau_dq_tmp;
+}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::updateForceDiff(
+    const boost::shared_ptr<ContactDataAbstract> &data, const MatrixXs &df_dx,
+    const MatrixXs &df_du) {
+  if (static_cast<std::size_t>(df_dx.rows()) != nc_ ||
+      static_cast<std::size_t>(df_dx.cols()) != state_->get_ndx())
+    throw_pretty("df_dx has wrong dimension");
+
+  if (static_cast<std::size_t>(df_du.rows()) != nc_ ||
+      static_cast<std::size_t>(df_du.cols()) != nu_)
+    throw_pretty("df_du has wrong dimension");
+
+  data->df_dx = -df_dx;
+  data->df_du = -df_du;
+}
+
+template <typename Scalar>
+boost::shared_ptr<ContactDataAbstractTpl<Scalar>>
+ContactModel6DLoopTpl<Scalar>::createData(
+    pinocchio::DataTpl<Scalar> *const data) {
+  return boost::allocate_shared<Data>(Eigen::aligned_allocator<Data>(), this,
+                                      data);
+}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::print(std::ostream &os) const {
+  os << "ContactModel6D {frame=" << state_->get_pinocchio()->frames[id_].name
+     << ", type=" << type_ << "}";
+}
+
+template <typename Scalar>
+const int ContactModel6DLoopTpl<Scalar>::get_joint1_id() const {
+  return joint1_id_;
+}
+
+template <typename Scalar>
+const int ContactModel6DLoopTpl<Scalar>::get_joint2_id() const {
+  return joint2_id_;
+}
+
+template <typename Scalar>
+const typename pinocchio::SE3Tpl<Scalar> &
+ContactModel6DLoopTpl<Scalar>::get_joint1_placement() const {
+  return joint1_placement_;
+}
+
+template <typename Scalar>
+const typename pinocchio::SE3Tpl<Scalar> &
+ContactModel6DLoopTpl<Scalar>::get_joint2_placement() const {
+  return joint2_placement_;
+}
+
+template <typename Scalar>
+const typename MathBaseTpl<Scalar>::Vector2s &
+ContactModel6DLoopTpl<Scalar>::get_gains() const {
+  return gains_;
+}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::set_joint1_id(const int joint1_id) {
+  joint1_id_ = joint1_id;
+}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::set_joint2_id(const int joint2_id) {
+  joint2_id_ = joint2_id;
+}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::set_joint1_placement(
+    const typename pinocchio::SE3Tpl<Scalar> &joint1_placement) {
+  joint1_placement_ = joint1_placement;
+}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::set_joint2_placement(
+    const typename pinocchio::SE3Tpl<Scalar> &joint2_placement) {
+  joint2_placement_ = joint2_placement;
+}
+
+template <typename Scalar>
+void ContactModel6DLoopTpl<Scalar>::set_gains(
+    const typename MathBaseTpl<Scalar>::Vector2s &gains) {
+  gains_ = gains;
+}
+
+}  // namespace crocoddyl
diff --git a/include/crocoddyl/multibody/contacts/multiple-contacts.hxx b/include/crocoddyl/multibody/contacts/multiple-contacts.hxx
index 8130bb4402..90c40f5fc8 100644
--- a/include/crocoddyl/multibody/contacts/multiple-contacts.hxx
+++ b/include/crocoddyl/multibody/contacts/multiple-contacts.hxx
@@ -7,6 +7,7 @@
 // All rights reserved.
 ///////////////////////////////////////////////////////////////////////////////
 
+#include "contact-6d-loop.hpp"
 namespace crocoddyl {
 
 template <typename Scalar>
@@ -250,14 +251,25 @@ void ContactModelMultipleTpl<Scalar>::updateForce(
         const Eigen::VectorBlock<const VectorXs, Eigen::Dynamic> force_i =
             force.segment(nc, nc_i);
         m_i->contact->updateForce(d_i, force_i);
+        const ContactModel6DLoopTpl<Scalar>* c =
+            dynamic_cast<ContactModel6DLoopTpl<Scalar>*>(m_i->contact.get());
+        if (c != nullptr) {
+          const ContactData6DLoopTpl<Scalar>* dc =
+              static_cast<ContactData6DLoopTpl<Scalar>*>(d_i.get());
+          const pinocchio::JointIndex joint1 = c->get_joint1_id();
+          const pinocchio::JointIndex joint2 = c->get_joint2_id();
+          data->fext[joint1] = dc->joint1_f;
+          data->fext[joint2] = dc->joint2_f;
+        } else {
 #if PINOCCHIO_VERSION_AT_LEAST(3, 0, 0)
-        const pinocchio::JointIndex joint =
-            state_->get_pinocchio()->frames[d_i->frame].parentJoint;
+          const pinocchio::JointIndex joint =
+              state_->get_pinocchio()->frames[d_i->frame].parentJoint;
 #else
-        const pinocchio::JointIndex joint =
-            state_->get_pinocchio()->frames[d_i->frame].parent;
+          const pinocchio::JointIndex joint =
+              state_->get_pinocchio()->frames[d_i->frame].parent;
 #endif
-        data->fext[joint] = d_i->fext;
+          data->fext[joint] = d_i->fext;
+        }
       } else {
         m_i->contact->setZeroForce(d_i);
       }
@@ -276,14 +288,25 @@ void ContactModelMultipleTpl<Scalar>::updateForce(
         const Eigen::VectorBlock<const VectorXs, Eigen::Dynamic> force_i =
             force.segment(nc, nc_i);
         m_i->contact->updateForce(d_i, force_i);
+        const ContactModel6DLoopTpl<Scalar>* c =
+            dynamic_cast<ContactModel6DLoopTpl<Scalar>*>(m_i->contact.get());
+        if (c != nullptr) {
+          const ContactData6DLoopTpl<Scalar>* dc =
+              static_cast<ContactData6DLoopTpl<Scalar>*>(d_i.get());
+          const pinocchio::JointIndex joint1 = c->get_joint1_id();
+          const pinocchio::JointIndex joint2 = c->get_joint2_id();
+          data->fext[joint1] = dc->joint1_f;
+          data->fext[joint2] = dc->joint2_f;
+        } else {
 #if PINOCCHIO_VERSION_AT_LEAST(3, 0, 0)
-        const pinocchio::JointIndex joint =
-            state_->get_pinocchio()->frames[d_i->frame].parentJoint;
+          const pinocchio::JointIndex joint =
+              state_->get_pinocchio()->frames[d_i->frame].parentJoint;
 #else
-        const pinocchio::JointIndex joint =
-            state_->get_pinocchio()->frames[d_i->frame].parent;
+          const pinocchio::JointIndex joint =
+              state_->get_pinocchio()->frames[d_i->frame].parent;
 #endif
-        data->fext[joint] = d_i->fext;
+          data->fext[joint] = d_i->fext;
+        }
         nc += nc_i;
       } else {
         m_i->contact->setZeroForce(d_i);
@@ -399,14 +422,7 @@ void ContactModelMultipleTpl<Scalar>::updateRneaDiff(
     assert_pretty(it_m->first == it_d->first,
                   "it doesn't match the contact name between data and model");
     if (m_i->active) {
-      switch (m_i->contact->get_type()) {
-        case pinocchio::ReferenceFrame::LOCAL:
-          break;
-        case pinocchio::ReferenceFrame::WORLD:
-        case pinocchio::ReferenceFrame::LOCAL_WORLD_ALIGNED:
-          pinocchio.dtau_dq += d_i->dtau_dq;
-          break;
-      }
+      pinocchio.dtau_dq += d_i->dtau_dq;
     }
   }
 }
diff --git a/include/crocoddyl/multibody/fwd.hpp b/include/crocoddyl/multibody/fwd.hpp
index c6e44428d8..1081341975 100644
--- a/include/crocoddyl/multibody/fwd.hpp
+++ b/include/crocoddyl/multibody/fwd.hpp
@@ -191,6 +191,11 @@ class ContactModel6DTpl;
 template <typename Scalar>
 struct ContactData6DTpl;
 
+template <typename Scalar>
+class ContactModel6DLoopTpl;
+template <typename Scalar>
+struct ContactData6DLoopTpl;
+
 // friction
 template <typename Scalar>
 class FrictionConeTpl;
@@ -362,6 +367,8 @@ typedef ContactModel3DTpl<double> ContactModel3D;
 typedef ContactData3DTpl<double> ContactData3D;
 typedef ContactModel6DTpl<double> ContactModel6D;
 typedef ContactData6DTpl<double> ContactData6D;
+typedef ContactModel6DLoopTpl<double> ContactModel6DLoop;
+typedef ContactData6DLoopTpl<double> ContactData6DLoop;
 
 typedef StateMultibodyTpl<double> StateMultibody;
 
diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt
index 1c1efc4dff..b141ba2483 100644
--- a/unittest/CMakeLists.txt
+++ b/unittest/CMakeLists.txt
@@ -28,6 +28,7 @@ set(${PROJECT_NAME}_CPP_TESTS
     test_constraints
     test_constraint_manager
     test_contacts
+    test_contacts_loop
     test_controls
     test_impulses
     test_multiple_contacts
diff --git a/unittest/factory/contact_loop.cpp b/unittest/factory/contact_loop.cpp
new file mode 100644
index 0000000000..21c7e93010
--- /dev/null
+++ b/unittest/factory/contact_loop.cpp
@@ -0,0 +1,90 @@
+///////////////////////////////////////////////////////////////////////////////
+// BSD 3-Clause License
+//
+// Copyright (C) 2019-2023, University of Edinburgh, Heriot-Watt University
+// Copyright note valid unless otherwise stated in individual files.
+// All rights reserved.
+///////////////////////////////////////////////////////////////////////////////
+
+#include "contact_loop.hpp"
+
+#include "crocoddyl/core/utils/exception.hpp"
+#include "crocoddyl/multibody/contacts/contact-6d-loop.hpp"
+
+namespace crocoddyl {
+namespace unittest {
+
+const std::vector<ContactLoopModelTypes::Type> ContactLoopModelTypes::all(
+    ContactLoopModelTypes::init_all());
+
+std::ostream& operator<<(std::ostream& os,
+                         const ContactLoopModelTypes::Type& type) {
+  switch (type) {
+    case ContactLoopModelTypes::ContactModel6DLoop_LOCAL:
+      os << "ContactLoopModel6D_LOCAL";
+      break;
+    case ContactLoopModelTypes::NbContactModelTypes:
+      os << "NbContactModelTypes";
+      break;
+    default:
+      os << "Unknown type";
+      break;
+  }
+  return os;
+}
+
+ContactLoopModelFactory::ContactLoopModelFactory() {}
+ContactLoopModelFactory::~ContactLoopModelFactory() {}
+
+boost::shared_ptr<crocoddyl::ContactModelAbstract>
+ContactLoopModelFactory::create(ContactLoopModelTypes::Type contact_type,
+                                PinocchioModelTypes::Type model_type,
+                                const int joint1_id,
+                                const pinocchio::SE3& joint1_placement,
+                                const int joint2_id,
+                                const pinocchio::SE3& joint2_placement,
+                                Eigen::Vector2d gains, std::size_t nu) const {
+  PinocchioModelFactory model_factory(model_type);
+  boost::shared_ptr<crocoddyl::StateMultibody> state =
+      boost::make_shared<crocoddyl::StateMultibody>(model_factory.create());
+  boost::shared_ptr<crocoddyl::ContactModelAbstract> contact;
+  if (nu == std::numeric_limits<std::size_t>::max()) {
+    nu = state->get_nv();
+  }
+  switch (contact_type) {
+    case ContactLoopModelTypes::ContactModel6DLoop_LOCAL:
+      contact = boost::make_shared<crocoddyl::ContactModel6DLoop>(
+          state, joint1_id, joint1_placement, joint2_id, joint2_placement,
+          pinocchio::ReferenceFrame::LOCAL, nu, gains);
+      break;
+    default:
+      throw_pretty(__FILE__ ": Wrong ContactLoopModelTypes::Type given");
+      break;
+  }
+  return contact;
+}
+
+boost::shared_ptr<crocoddyl::ContactModelAbstract>
+create_random_loop_contact() {
+  static bool once = true;
+  if (once) {
+    srand((unsigned)time(NULL));
+    once = false;
+  }
+  boost::shared_ptr<crocoddyl::ContactModelAbstract> contact;
+  ContactLoopModelFactory factory;
+  ;
+  if (rand() % 1 == 0) {
+    contact = factory.create(
+        ContactLoopModelTypes::ContactModel6DLoop_LOCAL,
+        PinocchioModelTypes::RandomHumanoid, 0,
+        pinocchio::SE3(Eigen::Matrix3d::Identity(), Eigen::Vector3d::Random()),
+        1,
+        pinocchio::SE3(Eigen::Matrix3d::Identity(), Eigen::Vector3d::Random()),
+        Eigen::Vector2d::Random());
+  }
+  return contact;
+}
+
+}  // namespace unittest
+}  // namespace crocoddyl
diff --git a/unittest/factory/contact_loop.hpp b/unittest/factory/contact_loop.hpp
new file mode 100644
index 0000000000..f348cc86d8
--- /dev/null
+++ b/unittest/factory/contact_loop.hpp
@@ -0,0 +1,60 @@
+///////////////////////////////////////////////////////////////////////////////
+// BSD 3-Clause License
+//
+// Copyright (C) 2019-2024, University of Edinburgh, Heriot-Watt University,
+//                          LAAS-CNRS
+// Copyright note valid unless otherwise stated in individual files.
+// All rights reserved.
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef CROCODDYL_CONTACTS_LOOP_FACTORY_HPP_
+#define CROCODDYL_CONTACTS_LOOP_FACTORY_HPP_
+
+#include <iostream>
+#include <limits>
+
+#include "crocoddyl/multibody/contact-base.hpp"
+#include "crocoddyl/multibody/contacts/multiple-contacts.hpp"
+#include "crocoddyl/multibody/numdiff/contact.hpp"
+#include "state.hpp"
+
+namespace crocoddyl {
+namespace unittest {
+
+struct ContactLoopModelTypes {
+  enum Type { ContactModel6DLoop_LOCAL, NbContactModelTypes };
+  static std::vector<Type> init_all() {
+    std::vector<Type> v;
+    v.reserve(NbContactModelTypes);
+    for (int i = 0; i < NbContactModelTypes; ++i) {
+      v.push_back((Type)i);
+    }
+    return v;
+  }
+  static const std::vector<Type> all;
+};
+
+std::ostream& operator<<(std::ostream& os,
+                         const ContactLoopModelTypes::Type& type);
+
+class ContactLoopModelFactory {
+ public:
+  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
+  explicit ContactLoopModelFactory();
+  ~ContactLoopModelFactory();
+
+  boost::shared_ptr<crocoddyl::ContactModelAbstract> create(
+      ContactLoopModelTypes::Type contact_type,
+      PinocchioModelTypes::Type model_type, const int joint1_id,
+      const pinocchio::SE3& joint1_placement, const int joint2_id,
+      const pinocchio::SE3& joint2_placement, Eigen::Vector2d gains,
+      std::size_t nu = std::numeric_limits<std::size_t>::max()) const;
+};
+
+boost::shared_ptr<crocoddyl::ContactModelAbstract> create_random_loop_contact();
+
+}  // namespace unittest
+}  // namespace crocoddyl
+
+#endif  // CROCODDYL_CONTACTS_LOOP_FACTORY_HPP_
diff --git a/unittest/test_contacts_loop.cpp b/unittest/test_contacts_loop.cpp
new file mode 100644
index 0000000000..df3f6427ec
--- /dev/null
+++ b/unittest/test_contacts_loop.cpp
@@ -0,0 +1,270 @@
+///////////////////////////////////////////////////////////////////////////////
+// BSD 3-Clause License
+//
+// Copyright (C) 2019-2023, University of Edinburgh, Heriot-Watt University
+// Copyright note valid unless otherwise stated in individual files.
+// All rights reserved.
+///////////////////////////////////////////////////////////////////////////////
+
+#define BOOST_TEST_NO_MAIN
+#define BOOST_TEST_ALTERNATIVE_INIT_API
+
+#include <boost/make_shared.hpp>
+#include <boost/shared_ptr.hpp>
+#include <pinocchio/algorithm/frames.hpp>
+#include <pinocchio/algorithm/kinematics-derivatives.hpp>
+
+#include "crocoddyl/multibody/contacts/contact-6d-loop.hpp"
+#include "factory/contact_loop.hpp"
+#include "unittest_common.hpp"
+
+using namespace crocoddyl::unittest;
+using namespace boost::unit_test;
+
+//----------------------------------------------------------------------------//
+
+void test_construct_data(ContactLoopModelTypes::Type contact_type,
+                         PinocchioModelTypes::Type model_type) {
+  // create the model
+  ContactLoopModelFactory factory;
+  boost::shared_ptr<crocoddyl::ContactModelAbstract> model =
+      factory.create(contact_type, model_type, 0, pinocchio::SE3::Identity(), 0,
+                     pinocchio::SE3::Identity(), Eigen::Vector2d::Random());
+
+  // Run the print function
+  std::ostringstream tmp;
+  tmp << *model;
+
+  // create the corresponding data object
+  const boost::shared_ptr<pinocchio::Model>& pinocchio_model =
+      model->get_state()->get_pinocchio();
+  pinocchio::Data pinocchio_data(*pinocchio_model.get());
+  boost::shared_ptr<crocoddyl::ContactDataAbstract> data =
+      model->createData(&pinocchio_data);
+}
+
+void test_calc_fetch_jacobians(ContactLoopModelTypes::Type contact_type,
+                               PinocchioModelTypes::Type model_type) {
+  // create the model
+  ContactLoopModelFactory factory;
+  boost::shared_ptr<crocoddyl::ContactModelAbstract> model =
+      factory.create(contact_type, model_type, 1, pinocchio::SE3::Random(), 2,
+                     pinocchio::SE3::Random(), Eigen::Vector2d::Random());
+
+  // create the corresponding data object
+  const boost::shared_ptr<pinocchio::Model>& pinocchio_model =
+      model->get_state()->get_pinocchio();
+  pinocchio::Data pinocchio_data(*pinocchio_model.get());
+  boost::shared_ptr<crocoddyl::ContactDataAbstract> data =
+      model->createData(&pinocchio_data);
+
+  // Compute the jacobian and check that the impulse model fetch it.
+  Eigen::VectorXd x = model->get_state()->rand();
+  crocoddyl::unittest::updateAllPinocchio(pinocchio_model.get(),
+                                          &pinocchio_data, x);
+
+  // Getting the jacobian from the model
+  model->calc(data, x);
+
+  // Check that only the Jacobian has been filled
+  BOOST_CHECK(!data->Jc.isZero());
+  if (model_type !=
+      PinocchioModelTypes::Hector) {  // this is due to Hector is a single rigid
+                                      // body system.
+    BOOST_CHECK(!data->a0.isZero());
+  }
+  BOOST_CHECK(data->da0_dx.isZero());
+  BOOST_CHECK(data->f.toVector().isZero());
+  BOOST_CHECK(data->df_dx.isZero());
+  BOOST_CHECK(data->df_du.isZero());
+}
+
+void test_calc_diff_fetch_derivatives(ContactLoopModelTypes::Type contact_type,
+                                      PinocchioModelTypes::Type model_type) {
+  // create the model
+  ContactLoopModelFactory factory;
+  boost::shared_ptr<crocoddyl::ContactModelAbstract> model =
+      factory.create(contact_type, model_type, 1, pinocchio::SE3::Random(), 2,
+                     pinocchio::SE3::Random(), Eigen::Vector2d::Random());
+
+  // create the corresponding data object
+  const boost::shared_ptr<pinocchio::Model>& pinocchio_model =
+      model->get_state()->get_pinocchio();
+  pinocchio::Data pinocchio_data(*pinocchio_model.get());
+  boost::shared_ptr<crocoddyl::ContactDataAbstract> data =
+      model->createData(&pinocchio_data);
+
+  // Compute the jacobian and check that the impulse model fetch it.
+  Eigen::VectorXd x = model->get_state()->rand();
+  crocoddyl::unittest::updateAllPinocchio(pinocchio_model.get(),
+                                          &pinocchio_data, x);
+
+  // Getting the jacobian from the model
+  model->calc(data, x);
+  model->calcDiff(data, x);
+
+  // Check that nothing has been computed and that all value are initialized to
+  // 0
+  BOOST_CHECK(!data->Jc.isZero());
+  if (model_type !=
+      PinocchioModelTypes::Hector) {  // this is due to Hector is a single rigid
+                                      // body system.
+    BOOST_CHECK(!data->a0.isZero());
+    BOOST_CHECK(!data->da0_dx.isZero());
+  }
+  BOOST_CHECK(data->f.toVector().isZero());
+  BOOST_CHECK(data->df_dx.isZero());
+  BOOST_CHECK(data->df_du.isZero());
+}
+
+void test_update_force(ContactLoopModelTypes::Type contact_type,
+                       PinocchioModelTypes::Type model_type) {
+  // create the model
+  ContactLoopModelFactory factory;
+  boost::shared_ptr<crocoddyl::ContactModelAbstract> model =
+      factory.create(contact_type, model_type, 1, pinocchio::SE3::Random(), 2,
+                     pinocchio::SE3::Random(), Eigen::Vector2d::Random());
+
+  // create the corresponding data object
+  const boost::shared_ptr<pinocchio::Model>& pinocchio_model =
+      model->get_state()->get_pinocchio();
+  pinocchio::Data pinocchio_data(*pinocchio_model.get());
+  boost::shared_ptr<crocoddyl::ContactDataAbstract> data =
+      model->createData(&pinocchio_data);
+
+  // Create a random force and update it
+  Eigen::VectorXd f = Eigen::VectorXd::Random(data->Jc.rows());
+  model->updateForce(data, f);
+
+  // Check that nothing has been computed and that all value are initialized to
+  // 0
+  BOOST_CHECK(data->Jc.isZero());
+  BOOST_CHECK(data->a0.isZero());
+  BOOST_CHECK(data->da0_dx.isZero());
+  BOOST_CHECK(!data->f.toVector().isZero());
+  BOOST_CHECK(!data->fext.toVector().isZero());
+  BOOST_CHECK(data->df_dx.isZero());
+  BOOST_CHECK(data->df_du.isZero());
+}
+
+void test_update_force_diff(ContactLoopModelTypes::Type contact_type,
+                            PinocchioModelTypes::Type model_type) {
+  // create the model
+  ContactLoopModelFactory factory;
+  boost::shared_ptr<crocoddyl::ContactModelAbstract> model =
+      factory.create(contact_type, model_type, 1, pinocchio::SE3::Random(), 2,
+                     pinocchio::SE3::Random(), Eigen::Vector2d::Random());
+
+  // create the corresponding data object
+  const boost::shared_ptr<pinocchio::Model>& pinocchio_model =
+      model->get_state()->get_pinocchio();
+  pinocchio::Data pinocchio_data(*pinocchio_model.get());
+  boost::shared_ptr<crocoddyl::ContactDataAbstract> data =
+      model->createData(&pinocchio_data);
+
+  // Create a random force and update it
+  Eigen::MatrixXd df_dx =
+      Eigen::MatrixXd::Random(data->df_dx.rows(), data->df_dx.cols());
+  Eigen::MatrixXd df_du =
+      Eigen::MatrixXd::Random(data->df_du.rows(), data->df_du.cols());
+  model->updateForceDiff(data, df_dx, df_du);
+
+  // Check that nothing has been computed and that all value are initialized to
+  // 0
+  BOOST_CHECK(data->Jc.isZero());
+  BOOST_CHECK(data->a0.isZero());
+  BOOST_CHECK(data->da0_dx.isZero());
+  BOOST_CHECK(data->f.toVector().isZero());
+  BOOST_CHECK(data->fext.toVector().isZero());
+  BOOST_CHECK(!data->df_dx.isZero());
+  BOOST_CHECK(!data->df_du.isZero());
+}
+
+void test_partial_derivatives_against_numdiff(
+    ContactLoopModelTypes::Type contact_type,
+    PinocchioModelTypes::Type model_type) {
+  using namespace boost::placeholders;
+
+  // create the model
+  ContactLoopModelFactory factory;
+  boost::shared_ptr<crocoddyl::ContactModelAbstract> model =
+      factory.create(contact_type, model_type, 0, pinocchio::SE3::Identity(), 0,
+                     pinocchio::SE3::Identity(), Eigen::Vector2d::Random());
+
+  // create the corresponding data object
+  pinocchio::Model& pinocchio_model =
+      *model->get_state()->get_pinocchio().get();
+  pinocchio::Data pinocchio_data(pinocchio_model);
+  boost::shared_ptr<crocoddyl::ContactDataAbstract> data =
+      model->createData(&pinocchio_data);
+
+  // Create the equivalent num diff model and data.
+  crocoddyl::ContactModelNumDiff model_num_diff(model);
+  const boost::shared_ptr<crocoddyl::ContactDataAbstract>& data_num_diff =
+      model_num_diff.createData(&pinocchio_data);
+
+  // Generating random values for the state
+  const Eigen::VectorXd x = model->get_state()->rand();
+
+  // Compute all the pinocchio function needed for the models.
+  crocoddyl::unittest::updateAllPinocchio(&pinocchio_model, &pinocchio_data, x);
+
+  // set the function that needs to be called at every step of the numdiff
+  std::vector<crocoddyl::ContactModelNumDiff::ReevaluationFunction> reevals;
+  reevals.push_back(boost::bind(&crocoddyl::unittest::updateAllPinocchio,
+                                &pinocchio_model, &pinocchio_data, _1, _2));
+  model_num_diff.set_reevals(reevals);
+
+  // Computing the contact derivatives
+  model->calc(data, x);
+  model->calcDiff(data, x);
+  model_num_diff.calc(data_num_diff, x);
+  model_num_diff.calcDiff(data_num_diff, x);
+  // Tolerance defined as in
+  // http://www.it.uom.gr/teaching/linearalgebra/NumericalRecipiesInC/c5-7.pdf
+  double tol = std::pow(model_num_diff.get_disturbance(), 1. / 3.);
+  BOOST_CHECK((data->da0_dx - data_num_diff->da0_dx).isZero(tol));
+}
+
+//----------------------------------------------------------------------------//
+
+void register_contact_model_unit_tests(ContactLoopModelTypes::Type contact_type,
+                                       PinocchioModelTypes::Type model_type) {
+  boost::test_tools::output_test_stream test_name;
+  test_name << "test_" << contact_type << "_" << model_type;
+  std::cout << "Running " << test_name.str() << std::endl;
+  test_suite* ts = BOOST_TEST_SUITE(test_name.str());
+  ts->add(BOOST_TEST_CASE(
+      boost::bind(&test_construct_data, contact_type, model_type)));
+  ts->add(BOOST_TEST_CASE(
+      boost::bind(&test_calc_fetch_jacobians, contact_type, model_type)));
+  ts->add(BOOST_TEST_CASE(boost::bind(&test_calc_diff_fetch_derivatives,
+                                      contact_type, model_type)));
+  ts->add(BOOST_TEST_CASE(
+      boost::bind(&test_update_force, contact_type, model_type)));
+  ts->add(BOOST_TEST_CASE(
+      boost::bind(&test_update_force_diff, contact_type, model_type)));
+  ts->add(BOOST_TEST_CASE(boost::bind(&test_partial_derivatives_against_numdiff,
+                                      contact_type, model_type)));
+  framework::master_test_suite().add(ts);
+}
+
+bool init_function() {
+  for (size_t contact_type = 0;
+       contact_type < ContactLoopModelTypes::all.size(); ++contact_type) {
+    for (size_t model_type = 0; model_type < PinocchioModelTypes::all.size();
+         ++model_type) {
+      if (model_type == PinocchioModelTypes::Hector) {
+        continue;
+      }
+      register_contact_model_unit_tests(
+          ContactLoopModelTypes::all[contact_type],
+          PinocchioModelTypes::all[model_type]);
+    }
+  }
+  return true;
+}
+
+int main(int argc, char** argv) {
+  return ::boost::unit_test::unit_test_main(&init_function, argc, argv);
+}