static int co = 0; co++;
void update_fermion_momenta(double step, const int S) {  int i,mu;  double tmp;  su3adj *xm,*deriv;#ifdef _COMP_FORCE  double sum=0., max=0.;  double sum2=0.;#endif  derivative_psf(S, g_eps_sq_force);#ifdef MPI  xchange_deri();#endif  for(i = 0; i < VOLUME; i++){    for(mu=0;mu<4;mu++){      xm=&moment[i][mu];      deriv=&df0[i][mu];#ifdef _COMP_FORCE      sum2 = _su3adj_square_norm(*deriv);       sum+= sum2;      if(sum2 > max) max = sum2;#endif      /* This 2* is coming from what? */      tmp = 2.*step;      _minus_const_times_mom(*xm,tmp,*deriv);     }  }#ifdef _COMP_FORCE#ifdef MPI  MPI_Reduce(&sum, &sum2, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);  sum = sum2;  MPI_Reduce(&max, &sum2, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);  max = sum2;#endif  if(g_proc_id == 0) {    printf("fermionforce%d %e max %e\n", S, sum/((double)(VOLUME*g_nproc))/4., max);    fflush(stdout);  }#endif}