         Mat33 PED =
             Potential_Energy_Derivative_Deviatoric(marker_Fe_hat[p], marker_Fp[p], device_settings.mu,
                                                    device_settings.lambda, device_settings.hardening_coefficient);
+
+        // Print(PED, "PED");
+
         Mat33 vPEDFepT = marker_volume[p] * MultTranspose(PED, marker_Fe[p]);
         real JE = Determinant(marker_Fe[p]);  //
         real JP = Determinant(marker_Fp[p]);

        Mat33 Fe_hat_inv_transpose = InverseTranspose(marker_Fe_hat[p]);

        real dJ = J * InnerProduct(Fe_hat_inv_transpose, delta_F);
        Mat33 dF_inverse_transposed = -Fe_hat_inv_transpose * Transpose(delta_F) * Fe_hat_inv_transpose;
        Mat33 dJF_inverse_transposed = dJ * Fe_hat_inv_transpose + J * dF_inverse_transposed;
        Mat33 RD = Rotational_Derivative(marker_Fe_hat[p], delta_F);

            real3 vnew(v_array[current_node * 3 + 0], v_array[current_node * 3 + 1], v_array[current_node * 3 + 2]);
            real3 vold(old_vel_node_mpm[current_node * 3 + 0], old_vel_node_mpm[current_node * 3 + 1],
                       old_vel_node_mpm[current_node * 3 + 2]);
            real3 v0 = vold + vnew;  //

            real3 v1 = dN(xi - current_node_location, inv_bin_edge);  //

            delta_F += OuterProduct(v0, v1);  //
            )

        real plastic_determinant = Determinant(marker_Fp[p]);
