diff --git a/monomial/ndrat_monomial.c b/monomial/ndrat_monomial.c
index a79e9b738..77c7cb02d 100644
--- a/monomial/ndrat_monomial.c
+++ b/monomial/ndrat_monomial.c
@@ -71,153 +71,6 @@ void nd_set_global_parameter(monomial * const mnl) {
   return;
 }
 
-
-
-void print_spinor_updow_shift( double** up, double** down,   int odd, int num_shifts) {
-  tm_stopwatch_push(&g_timers, __func__, "");
-
-  const int change_sign[4] = {-1, 1, 1, -1};
-  const int change_spin[4] = {3, 2, 1, 0};
-  const int Vh = VOLUME/2;
-
-  // now copy and reorder from tempSpinor to spinor
-  for(int shift = 0; shift < num_shifts; shift++){
-  for(int ud = 0; ud < 2; ud++){
-    double **sp;
-    if (ud==0) sp=up;
-    if (ud==1) sp=down;
-    for( int x0=0; x0<T; x0++ )
-      for( int x1=0; x1<LX; x1++ )
-        for( int x2=0; x2<LY; x2++ )
-          for( int x3=0; x3<LZ; x3++ ) {
-            const int q_eo_idx = (x1 + LX*x2 + LY*LX*x3 + LZ*LY*LX*x0)/2;
-            const int tm_eo_idx = (x3 + LZ*x2 + LY*LZ*x1 + LX*LY*LZ*x0)/2;
-
-            const int oddBit = (x0+x1+x2+x3) & 1;
-            if( oddBit == odd ){
-              for(int tm_spin = 0; tm_spin < 4; tm_spin++){
-                for(int col = 0; col < 3; col++){
-                  if (g_proc_id==0) printf("MARCOfrom TMLQCD (%d %d %d %d),  %d %d, s=%d  ud=%d    %g  %g\n",x0,x1,x2,x3, tm_spin, col ,
-                    shift, ud,
-                    sp[shift][24*tm_eo_idx + 6*tm_spin + 2*col + 0],
-                    sp[shift][24*tm_eo_idx + 6*tm_spin + 2*col + 1]);
-                }
-              }
-            }
-          }
-  }
-  }
-  tm_stopwatch_pop(&g_timers, 0, 0, "TM_QUDA");
-}      
-
-void print_spinor_updow( double* up, double* down,   int odd , int shift) {
-  tm_stopwatch_push(&g_timers, __func__, "");
-
-  const int change_sign[4] = {-1, 1, 1, -1};
-  const int change_spin[4] = {3, 2, 1, 0};
-  const int Vh = VOLUME/2;
-
-  // now copy and reorder from tempSpinor to spinor
-  
-  for(int ud = 0; ud < 2; ud++){
-    double *sp;
-    if (ud==0) sp=up;
-    if (ud==1) sp=down;
-    for( int x0=0; x0<T; x0++ )
-      for( int x1=0; x1<LX; x1++ )
-        for( int x2=0; x2<LY; x2++ )
-          for( int x3=0; x3<LZ; x3++ ) {
-            const int q_eo_idx = (x1 + LX*x2 + LY*LX*x3 + LZ*LY*LX*x0)/2;
-            const int tm_eo_idx = (x3 + LZ*x2 + LY*LZ*x1 + LX*LY*LZ*x0)/2;
-
-            const int oddBit = (x0+x1+x2+x3) & 1;
-            if( oddBit == odd ){
-              for(int tm_spin = 0; tm_spin < 4; tm_spin++){
-                for(int col = 0; col < 3; col++){
-                  if (g_proc_id==0) printf("MARCOfrom TMLQCD (%d %d %d %d),  %d %d, s=%d  ud=%d    %g  %g\n",x0,x1,x2,x3, tm_spin, col ,
-                    shift, ud,
-                    sp[24*tm_eo_idx + 6*tm_spin + 2*col + 0],
-                    sp[24*tm_eo_idx + 6*tm_spin + 2*col + 1]);
-                }
-              }
-            }
-  
-  }
-  }
-  tm_stopwatch_pop(&g_timers, 0, 0, "TM_QUDA");
-}      
-
-
-void set_spinor_updow_shift( double** up, double** down,   int odd, int num_shifts) {
-  tm_stopwatch_push(&g_timers, __func__, "");
-
-  const int change_sign[4] = {-1, 1, 1, -1};
-  const int change_spin[4] = {3, 2, 1, 0};
-  const int Vh = VOLUME/2;
-
-  // now copy and reorder from tempSpinor to spinor
-  for(int shift = 0; shift < num_shifts; shift++){
-  for(int ud = 1; ud < 2; ud++){
-    double **sp;
-    if (ud==0) sp=up;
-    if (ud==1) sp=down;
-    for( int x0=0; x0<T; x0++ )
-      for( int x1=0; x1<LX; x1++ )
-        for( int x2=0; x2<LY; x2++ )
-          for( int x3=0; x3<LZ; x3++ ) {
-            const int q_eo_idx = (x1 + LX*x2 + LY*LX*x3 + LZ*LY*LX*x0)/2;
-            const int tm_eo_idx = (x3 + LZ*x2 + LY*LZ*x1 + LX*LY*LZ*x0)/2;
-
-            const int oddBit = (x0+x1+x2+x3) & 1;
-            if( oddBit == odd ){
-              for(int q_spin = 0; q_spin < 4; q_spin++){
-                const int tm_spin = change_spin[q_spin];
-                for(int col = 0; col < 3; col++){
-                    sp[shift][24*tm_eo_idx + 6*tm_spin + 2*col + 0]=0;
-                    sp[shift][24*tm_eo_idx + 6*tm_spin + 2*col + 1]=0;
-                }
-              }
-            }
-          }
-  }
-  }
-  // for(int shift = 0; shift < num_shifts; shift++){
-  // for(int ud = 0; ud < 1; ud++){
-  //   double **sp;
-  //   if (ud==0) sp=up;
-  //   if (ud==1) sp=down;
-  //   for( int x0=0; x0<T; x0++ )
-  //     for( int x1=0; x1<LX; x1++ )
-  //       for( int x2=0; x2<LY; x2++ )
-  //         for( int x3=0; x3<LZ; x3++ ) {
-  //           const int q_eo_idx = (x1 + LX*x2 + LY*LX*x3 + LZ*LY*LX*x0)/2;
-  //           const int tm_eo_idx = (x3 + LZ*x2 + LY*LZ*x1 + LX*LY*LZ*x0)/2;
-
-  //           const int oddBit = (x0+x1+x2+x3) & 1;
-  //           if( oddBit == odd ){
-  //             for(int q_spin = 0; q_spin < 4; q_spin++){
-  //               const int tm_spin = change_spin[q_spin];
-  //               for(int col = 0; col < 3; col++){
-  //                   sp[shift][24*tm_eo_idx + 6*3 + 2*0 + 0]=1;
-  //                   sp[shift][24*tm_eo_idx + 6*tm_spin + 2*col + 1]=0;
-  //               }
-  //             }
-  //           }
-  //         }
-  // }
-  // }
-  // up[0][24*0 + 6*0 + 2*0 + 0]=1;
-  // up[0][24*0 + 6*0 + 2*0 + 1]=2;
-  // up[0][24*0 + 6*0 + 2*1 + 0]=3;
-  // up[0][24*1 + 6*0 + 2*0 + 0]=1;
-  // up[0][24*2 + 6*0 + 2*0 + 0]=1;
-  // up[0][24*3 + 6*0 + 2*0 + 0]=1;
-  // up[0][24*4 + 6*0 + 2*0 + 0]=1;
-  // up[0][24*5 + 6*0 + 2*0 + 0]=1;
-  // down[0][24*0 + 6*0 + 2*0 + 0]=2;
-  tm_stopwatch_pop(&g_timers, 0, 0, "TM_QUDA");
-}      
-
 /********************************************
  *
  * Here \delta S_b is computed
@@ -267,8 +120,6 @@ void ndrat_derivative(const int id, hamiltonian_field_t * const hf) {
   mnl->iter1 += solve_mms_nd(g_chi_up_spinor_field, g_chi_dn_spinor_field,
                              mnl->pf, mnl->pf2, &(mnl->solver_params) );
 
-  // set_spinor_updow_shift(  g_chi_up_spinor_field,  g_chi_dn_spinor_field,   1, mnl->rat.np);
-
   if ( mnl->external_library == QUDA_LIB){
     if(!mnl->even_odd_flag) {
       fatal_error("QUDA support only even_odd_flag",__func__);
@@ -304,9 +155,6 @@ void ndrat_derivative(const int id, hamiltonian_field_t * const hf) {
     if (g_debug_level > 3){
       su3adj **given = hf->derivative;
       hf->derivative = debug_derivative;
-      // int store_solver = mnl->solver;
-      // mnl->solver_params.external_inverter = NO_EXT_INV;
-      // mnl->solver = CGMMSND;
       mnl->external_library = NO_EXT_LIB;
       tm_debug_printf( 0, 3, "Recomputing the derivative from tmLQCD\n");
       ndrat_derivative(id, hf);
@@ -320,7 +168,6 @@ void ndrat_derivative(const int id, hamiltonian_field_t * const hf) {
   #endif // no other option, TM_USE_QUDA already checked by solver
   }
   else{
-    // print_spinor_updow_shift( g_chi_up_spinor_field, g_chi_dn_spinor_field,   1, mnl->rat.np);
     // TODO: check that we did not breack anything rearranging the loop
     // for(int j = (mnl->rat.np-1); j > -1; j--) {
     for(int j = 0; j < mnl->rat.np; j++) {
@@ -331,9 +178,6 @@ void ndrat_derivative(const int id, hamiltonian_field_t * const hf) {
         Qsw_tau1_sub_const_ndpsi(mnl->w_fields[0], mnl->w_fields[1],
               g_chi_up_spinor_field[j], g_chi_dn_spinor_field[j], 
               -I*mnl->rat.mu[j], 1., mnl->EVMaxInv);
-        printf("MARCO:printpar mnl->rat.mu[%d] =%g\n",j,mnl->rat.mu[j]);
-        printf("MARCO:printpar mnl->rat.rmu[%d] =%g\n",j,mnl->rat.rmu[j]);
-        printf("MARCO:printpar mnl->EVMaxInv =%g\n",mnl->EVMaxInv);
         tm_stopwatch_pop(&g_timers, 0, 1, "");
         
         /* Get the even parts X_j,e */
@@ -342,12 +186,7 @@ void ndrat_derivative(const int id, hamiltonian_field_t * const hf) {
         H_eo_sw_ndpsi(mnl->w_fields[2], mnl->w_fields[3], 
           g_chi_up_spinor_field[j], g_chi_dn_spinor_field[j]);
         tm_stopwatch_pop(&g_timers, 0, 1, "");
-        // mul_gamma5(g_chi_up_spinor_field[j], VOLUME/2);
-        // mul_gamma5(g_chi_dn_spinor_field[j], VOLUME/2);
-        // print_spinor_updow( g_chi_up_spinor_field[j], g_chi_dn_spinor_field[j],   1, j); //x[i][parity]
-        // print_spinor_updow( mnl->w_fields[2], mnl->w_fields[3],   0, j); // x[i][other_parity]
-        // print_spinor_updow( mnl->w_fields[0], mnl->w_fields[1],   1, j); // p[i][parity] 
-      } else {
+       } else {
         // multiply with Q_h * tau^1 + i mu_j to get Y_j,o (odd sites)
         // needs phmc_Cpol = 1 to work for ndrat!
         tm_stopwatch_push(&g_timers, "Q_tau1_sub_const_ndpsi", "");
@@ -381,7 +220,6 @@ void ndrat_derivative(const int id, hamiltonian_field_t * const hf) {
         H_eo_tm_ndpsi(mnl->w_fields[4], mnl->w_fields[5], mnl->w_fields[0], mnl->w_fields[1], EO);
         tm_stopwatch_pop(&g_timers, 0, 1, "");
       }
-      print_spinor_updow( mnl->w_fields[4], mnl->w_fields[5],   0, j); // p[i][other_parity] // p_e
 
       /* X_j,o \delta M_oe Y_j,e */
       deriv_Sb(OE, g_chi_up_spinor_field[j], mnl->w_fields[4],