corr_nonlocal_fn *= -1;
void wallNuclFormFac(XMLWriter& xml,		     const multi1d<LatticeColorMatrix>& u, 		     const LatticePropagator& forw_u_prop,		     const LatticePropagator& back_u_prop, 		     const LatticePropagator& forw_d_prop,		     const LatticePropagator& back_d_prop, 		     const SftMom& phases,		     int t0, int t_sink){  START_CODE("wallNuclFormFac");  // Length of lattice in j_decay direction and 3pt correlations fcns  int length = phases.numSubsets();  multi1d<Complex> local_cur3ptfn(length);  multi1d<Complex> nonlocal_cur3ptfn(length);  int G5 = Ns*Ns-1;    // Project propagator onto zero momentum: Do a slice-wise sum.  Propagator u_x2 = sum(forw_u_prop, phases.getSet()[t_sink]);  Propagator d_x2 = sum(forw_d_prop, phases.getSet()[t_sink]);  LatticePropagator anti_u_prop = adj(Gamma(G5)*back_u_prop*Gamma(G5));  LatticePropagator anti_d_prop = adj(Gamma(G5)*back_d_prop*Gamma(G5));  // Loop over appropriate form-factor contractions for this system  XMLArrayWriter xml_seq_src(xml, 4);  push(xml_seq_src, "FormFac");  for (int seq_src = 0; seq_src < 4; ++seq_src)   {    push(xml_seq_src);    write(xml_seq_src, "seq_src", seq_src);    // Loop over gamma matrices of the insertion current of insertion current    XMLArrayWriter xml_array(xml_seq_src, Nd);    push(xml_array, "Insertions");    for(int mu = 0; mu < Nd; ++mu)    {      int gamma_value = 1 << mu;      push(xml_array);      write(xml_array, "mu", mu);      write(xml_array, "gamma_value", gamma_value);      LatticeComplex corr_local_fn;      LatticeComplex corr_nonlocal_fn;      switch (seq_src)      {      case 0:      case 2:      {	// "\bar u O u" insertion in proton, ie. "(u C gamma_5 d) u"	// The local non-conserved current contraction	LatticeSpinMatrix local_contract = 	  wallNuclUContract(LatticePropagator(anti_u_prop*Gamma(gamma_value)*forw_u_prop), 			    u_x2, d_x2);	// Construct the non-local (possibly conserved) current contraction	LatticeSpinMatrix nonlocal_contract = 	  wallNuclUContract(nonlocalCurrentProp(u, mu, forw_u_prop, anti_u_prop), 			    u_x2, d_x2);	if (seq_src == 0)	{	  QDPIO::cerr << "do seq_src 0" << endl;	  /* "\bar u O u" insertion in proton, ie. "(u C gamma_5 d) u" */	  /* T = (1 + gamma_4) / 2 = (1 + Gamma(8)) / 2 */	  /* C gamma_5 = Gamma(5) = - (C gamma_5)^T */		  // The local non-conserved vector-current matrix element 	  corr_local_fn = 0.5 * traceSpin(local_contract + Gamma(8)*local_contract);	  // The nonlocal (possibly conserved) current matrix element 	  corr_nonlocal_fn = 0.25 * traceSpin(nonlocal_contract + Gamma(8)*nonlocal_contract);	}	else	{	  QDPIO::cerr << "do seq_src 2" << endl;	  /* "\bar u O u" insertion in proton, ie. "(u C gamma_5 d) u" */	  /* T = \Sigma_3 (1 + gamma_4) / 2 = -i (Gamma(3) + Gamma(11)) / 2 */	  /* C gamma_5 = Gamma(5) = - (C gamma_5)^T */	  // The local non-conserved vector-current matrix element 	  corr_local_fn = 0.5 * timesMinusI(traceSpin(Gamma(3)*local_contract + Gamma(11)*local_contract));	  // The nonlocal (possibly conserved) current matrix element 	  corr_nonlocal_fn = 0.25 * timesMinusI(traceSpin(Gamma(3)*nonlocal_contract + Gamma(11)*nonlocal_contract));	}      }      break;      case 1:      case 3:      {	// "\bar d O d" insertion in proton, ie. "(u C gamma_5 d) u"	// The local non-conserved current contraction	LatticeSpinMatrix local_contract = 	  wallNuclUContract(LatticePropagator(anti_d_prop*Gamma(gamma_value)*forw_d_prop), 			    u_x2, d_x2);	// Construct the non-local (possibly conserved) current contraction	LatticeSpinMatrix nonlocal_contract = 	  wallNuclUContract(nonlocalCurrentProp(u, mu, forw_d_prop, anti_d_prop), 			    u_x2, d_x2);	if (seq_src == 1)	{	  QDPIO::cerr << "do seq_src 1" << endl;	  /* "\bar d O d" insertion in proton, ie. "(u C gamma_5 d) u" */	  /* T = (1 + gamma_4) / 2 = (1 + Gamma(8)) / 2 */	  // The local non-conserved vector-current matrix element 	  corr_local_fn = 0.5 * traceSpin(local_contract + Gamma(8)*local_contract);	  // The nonlocal (possibly conserved) current matrix element 	  corr_nonlocal_fn = 0.25 * traceSpin(nonlocal_contract + Gamma(8)*nonlocal_contract);	}	else	{	  QDPIO::cerr << "do seq_src 3" << endl;	  /* "\bar d O d" insertion in proton, ie. "(u C gamma_5 d) u" */	  /* T = \Sigma_3 (1 + gamma_4) / 2 = -i (Gamma(3) + Gamma(11)) / 2 */	  /* C gamma_5 = Gamma(5) = - (C gamma_5)^T */	  // The local non-conserved vector-current matrix element 	  corr_local_fn = 0.5 * timesMinusI(traceSpin(Gamma(3)*local_contract + Gamma(11)*local_contract));	  // The nonlocal (possibly conserved) current matrix element 	  corr_nonlocal_fn = 0.25 * timesMinusI(traceSpin(Gamma(3)*nonlocal_contract + Gamma(11)*nonlocal_contract));	}      }      break;      default:	QDP_error_exit("Unknown sequential source type", seq_src);      }      multi2d<DComplex> hsum_local = phases.sft(corr_local_fn);      multi2d<DComplex> hsum_nonlocal = phases.sft(corr_nonlocal_fn);        XMLArrayWriter xml_inser_mom(xml_array, phases.numMom());      push(xml_inser_mom, "Momenta");      // Loop over insertion momenta and print out results      for(int inser_mom_num=0; inser_mom_num<phases.numMom(); ++inser_mom_num)       {	push(xml_inser_mom);	write(xml_inser_mom, "inser_mom_num", inser_mom_num);	write(xml_inser_mom, "inser_mom", phases.numToMom(inser_mom_num)) ;	for (int t=0; t < length; ++t) 	{	  int t_eff = (t - t0 + length) % length;	  local_cur3ptfn[t_eff] = Complex(hsum_local[inser_mom_num][t]);	  nonlocal_cur3ptfn[t_eff] = Complex(hsum_nonlocal[inser_mom_num][t]);	} // end for(t)	// Print out the results	write(xml_inser_mom, "local_cur3ptfn", local_cur3ptfn);	write(xml_inser_mom, "nonlocal_cur3ptfn", nonlocal_cur3ptfn);	pop(xml_inser_mom);  // elem      } // end for(inser_mom_num)      pop(xml_inser_mom);    // Momenta      pop(xml_array);        // elem    } // end for(gamma_value)                                pop(xml_array);          // FormFac    pop(xml_seq_src);        // elem  } // end for(seq_src)                              pop(xml_seq_src);          // WallNuclFormFac  END_CODE("wallNuclFormFac");}