From 159bf49ddd746f1960b2ad290c3206092b1a6d4e Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Tue, 18 Jul 2023 11:50:38 +0200 Subject: [PATCH] C(t) volume average over spinor scalar products --- meas/correlators.c | 54 +++++++++++++++++++++------------------------- 1 file changed, 25 insertions(+), 29 deletions(-) diff --git a/meas/correlators.c b/meas/correlators.c index 69205e1ef..6f2f56c01 100644 --- a/meas/correlators.c +++ b/meas/correlators.c @@ -546,9 +546,9 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, } } - // (no even-odd): source+propagator, doublet, spin dilution index, even-odd, flavor, position (+ Dirac, color) - // (psi[i_sp][phi][beta][f][x])[alpha][c] // last 3 indices are proper of the spinor struct - spinor *****arr_spinor = (spinor *****) callocMultiDimensional({2,2,4,2,VOLUME}, 5, sizeof(spinor)); + // (no even-odd): source+propagator, doublet, spin dilution index, flavor proj, flavor index, position (+ Dirac, color) + // (psi[i_sp][d][beta][F][f][x])[alpha][c] // last 3 indices are proper of the spinor struct + spinor ******arr_spinor = (spinor ******) callocMultiDimensional({2,2,4,2,2,VOLUME}, 6, sizeof(spinor)); // now we switch from even-odd representation to standard for (size_t i_sp = 0; i_sp < 2; i_sp++) { // source or sink @@ -581,40 +581,36 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, for (i = j; i < j + LX * LY * LZ; i++) { // light correlators for (size_t beta = 0; beta < 4; beta++) { // spin dilution - double sign = 1.0; // the light source sould be multiplied by gamma_5 - if(beta >= 2){ + double sign = 1.0; // the light source sould be multiplied by gamma_5 + if (beta >= 2) { sign = -1.0; } - res += sign*_spinor_prod_re(arr_spinor[1][0][beta][0][i], arr_spinor[1][0][beta][0][i]); + res += + sign * _spinor_prod_re(arr_spinor[1][0][beta][0][i], arr_spinor[1][0][beta][0][i]); _gamma0(phi, arr_spinor[1][0][beta][0][i]); - respa += sign*_spinor_prod_re(arr_spinor[1][0][beta][0][i], phi); + respa += sign * _spinor_prod_re(arr_spinor[1][0][beta][0][i], phi); _gamma5(phi, phi); - resp4 += sign*_spinor_prod_im(arr_spinor[1][0][beta][0][i], phi); + resp4 += sign * _spinor_prod_im(arr_spinor[1][0][beta][0][i], phi); } - for (size_t i_f = 0; i_f < 2; i_f++) { - for (size_t j_f = 0; j_f < 2; j_f++) { - for (size_t g1 = 0; g1 < 2; g1++) { - for (size_t g2 = 0; g2 < 2; g2++) { - for (size_t beta1 = 0; beta1 < 4; beta1++) { - for (size_t beta2 = 0; beta2 < 4; beta2++) { - double dum1, dum2 = 0.0, 0.0; // dummy variables - - // 1st contribution - phi = arr_spinor[1][1][beta_2][i_f][i]; - if (Gamma_1 == 1) { // Gamma_1 = 1 + const int f0 = 0; // flavor index of the up + for (size_t F = 0; F < 4; F++) { // flavor projection + for (size_t beta = 0; beta < 4; beta++) { // spin dilution + spinor psi_u = arr_spinor[1][0][beta][f0][f0][i]; + for (size_t hi = 0; hi < 2; hi++) { + for (size_t hj = 0; hj < 2; hj++) { + for (size_t g1 = 0; g1 < 2; g1++) { + for (size_t g2 = 0; g2 < 2; g2++) { + double dum = 0.0; // dummy variable + + // heavy doublet spinor propagator + phi = arr_spinor[1][1][beta][hj][hi][i]; + if (g1 == 1) { // Gamma_1 = 1 _gamma5(phi, phi); } - _spinor_scalar_prod(dum1, arr_spinor[0][0][beta_1][1 - i_f][i], phi); + _spinor_scalar_prod(dum, psi_u, phi); - // 2nd contribution - phi = arr_spinor[1][0][beta_1][1 - j_f][i]; - if (Gamma_2 == 1) { // Gamma_1 = gamma_5 - _gamma5(phi, phi); - } - _spinor_scalar_prod(dum2, arr_spinor[0][1][beta_2][j_f][i], phi); - - res_hihj_g1g2[i_f][j_f][beta_1][beta_2] += dum1 * dum2; + res_hihj_g1g2[hi][hj][g1][g2] += dum; } } } @@ -652,7 +648,7 @@ void heavy_correlators_measurement(const int traj, const int id, const int ieo, #if defined TM_USE_MPI MPI_Reduce(&res_hihj_g1g2[hi][hj][g1][g2], &mpi_res_hihj_g1g2[hi][hj][g1][g2], 1, MPI_DOUBLE, MPI_SUM, 0, g_mpi_time_slices); res_hihj_g1g2[hi][hj][g1][g2] = mpi_res_hihj_g1g2[hi][hj][g1][g2]; - C_hihj_g1g2[hi][hj][g1][g2][t] = -eta_Gamma[g1] * res / (g_nproc_x * LX) / (g_nproc_y * LY) / (g_nproc_z * LZ); + sC_hihj_g1g2[hi][hj][g1][g2][t] = -eta_Gamma[g1] * res / (g_nproc_x * LX) / (g_nproc_y * LY) / (g_nproc_z * LZ); #else C_hihj_g1g2[hi][hj][g1][g2][t] = -eta_Gamma[g1] * res / (g_nproc_x * LX) / (g_nproc_y * LY) / (g_nproc_z * LZ); #endif