diff --git a/.gitignore b/.gitignore index 82a86fead..9489a23c2 100644 --- a/.gitignore +++ b/.gitignore @@ -39,5 +39,6 @@ doc/*.ps doc/*.txt doc/*.gpl doc/wrapper* +doc/omeas_heavy_mesons.html history_hmc_tm .vscode/* \ No newline at end of file diff --git a/doc/acta-ecologica-sinica.csl b/doc/acta-ecologica-sinica.csl new file mode 100644 index 000000000..1d78dcb63 --- /dev/null +++ b/doc/acta-ecologica-sinica.csl @@ -0,0 +1,14 @@ + + diff --git a/doc/bibliography.bib b/doc/bibliography.bib index 616d45e69..1e4ee7f9e 100644 --- a/doc/bibliography.bib +++ b/doc/bibliography.bib @@ -27,6 +27,18 @@ @article{Boyle:2017xcy SLACcitation = "%%CITATION = ARXIV:1711.04883;%%" } +@article{baron2011computing, + title={{Computing K and D meson masses with $N_f= 2 + 1 + 1$ twisted mass lattice QCD}}, + author={Baron, Remi and Boucaud, Philippe and Carbonell, Jaume and Drach, Vincent and Farchioni, Federico and Herdoiza, Gregorio and Jansen, Karl and Michael, Chris and Montvay, Istvan and Pallante, Elisabetta and others}, + journal={Computer Physics Communications}, + volume={182}, + number={2}, + pages={299--316}, + year={2011}, + publisher={Elsevier}, + doi = {https://doi.org/10.1016/j.cpc.2010.10.004} +} + @article{Babich:2011np, archiveprefix = "arXiv", author = "Babich, R. and Clark, M.A. and Joo, B. and Shi, G. and Brower, R.C. and others", @@ -139,6 +151,18 @@ @article{Abdel-Rehim:2004gx year = "2005" } +@article{foley2005practical, + title={Practical all-to-all propagators for lattice QCD}, + author={Foley, Justin and Juge, K Jimmy and Cais, Alan {\'O} and Peardon, Mike and Ryan, Sin{\'e}ad M and Skullerud, Jon-Ivar and TrinLat Collaboration and others}, + journal={Computer physics communications}, + volume={172}, + number={3}, + pages={145--162}, + year={2005}, + publisher={Elsevier}, + doi = {10.1016/j.cpc.2005.06.008} +} + @article{Abdel-Rehim:2005gz, author = "Abdel-Rehim, Abdou M. and Lewis, Randy and Woloshyn, R. M.", eprint = "hep-lat/0503007", @@ -620,6 +644,22 @@ @article{Becher:1999he year = "1999" } +@article{PhysRevD.59.074503, + title = {Quark mass dependence of hadron masses from lattice QCD}, + author = {Foster, M. and Michael, C.}, + collaboration = {UKQCD Collaboration}, + journal = {Phys. Rev. D}, + volume = {59}, + issue = {7}, + pages = {074503}, + numpages = {10}, + year = {1999}, + month = {Feb}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevD.59.074503}, + url = {https://link.aps.org/doi/10.1103/PhysRevD.59.074503} +} + @article{Bietenholz:2004sa, author = "Bietenholz, W. and others", collaboration = "\xlf", diff --git a/doc/mathjax-config.html b/doc/mathjax-config.html new file mode 100644 index 000000000..8fc201d95 --- /dev/null +++ b/doc/mathjax-config.html @@ -0,0 +1,31 @@ + + \ No newline at end of file diff --git a/doc/omeas_heavy_mesons.qmd b/doc/omeas_heavy_mesons.qmd new file mode 100644 index 000000000..b1ef91910 --- /dev/null +++ b/doc/omeas_heavy_mesons.qmd @@ -0,0 +1,381 @@ +--- +title: Online measurements of heavy meson correlators +author: Simone Romiti +date: last-modified +jupyter: python3 +format-links: false # hide the existence of other formats +format: + html: + toc: true + toc-location: body + embed-resources: true # standalone html + html-math-method: mathjax # use the common latex commands + include-in-header: mathjax-config.html # equation numbering with sections + fontsize: 14pt + pdf: + pdf-engine: pdflatex + fontsize: 12pt +bibliography: bibliography.bib +csl: acta-ecologica-sinica.csl +fontsize: 16pt +--- + +### Preamble and notation +
+ Show Content + +`tmLQCD`can compute the correlators for the heavy mesons $K$ and $D$. +The twisted mass formulation of the heavy doublet $(s, c)$ is such that $K$ and $D$ mix in the spectral decomposition @baron2011computing. +In fact, the Dirac operator is not diagonal in flavor, +and the extraction of their masses requires the evaluation of a $(2 \times 2) \otimes (2 \times 2)$ matrix of correlators. +In the following we use the twisted basis $\chi$ for fermion fields: + +\begin{align} +\label{eq:LightFlavorSpinorFields} +\chi_{\ell} &= + \begin{bmatrix} + {\chi}_{u} \\ + 0 + \end{bmatrix} + + + \begin{bmatrix} + 0 \\ + {\chi}_{d} + \end{bmatrix} +\, , \\[1em] +\label{eq:HeavyFlavorSpinorFields} +\chi_{h} &= + \begin{bmatrix} + {\chi}_{c} \\ + 0 + \end{bmatrix} + + + \begin{bmatrix} + 0 \\ + {\chi}_{s} + \end{bmatrix} +\, . +\end{align} + + +
+ +## Correlators for the $K$ and $D$ + +
+ Show Content + +The required correlators are given by the following +expectation values on the interacting vacuum +${\bra{\Omega} \cdot \ket{\Omega} = \braket{\cdot}}$: + +\begin{equation} +\begin{split} +{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t) &= +- \eta_1 \frac{1}{V} +\sum_{\vec{x}} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) +\\ +&= +\eta_{1} \eta_{2} +\frac{1}{V} +\sum_{\vec{x}} +\braket{ +[\bar{\chi}_{d} \Gamma_1 \chi_{h_i}](x) +[\bar{\chi}_{d} \Gamma_2 \chi_{h_j}]^{\dagger}(0) +} +\, , +\end{split} +\end{equation} + +where $h_i=c,s$, $\Gamma_i = 1, \gamma_5$. +The minus sign will vanish when doing the Wick contractions, +while +$\eta_{\mathbb{1}} = 1$ and $\eta_{\gamma_5} = -1$ +are used just to compensate the sign change when daggering the spinor bilinears at the source: + +\begin{equation} +\eta_\Gamma (\bar{\chi}_a \Gamma \chi_b)^\dagger = +\eta_\Gamma \bar{\chi}_b \gamma_0 \Gamma^\dagger \gamma_0 \chi_a = +\bar{\chi}_b \Gamma \chi_a +\, . +\end{equation} + + +We now write the explicit Wick contractions for the object: + +\begin{equation} +\begin{split} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) &= +- \eta_{2} +\braket{ +[\bar{\chi}_{d} \Gamma_1 \chi_{h_i}](x) +[\bar{\chi}_{d} \Gamma_2 \chi_{h_j}]^{\dagger}(0) +} +\\ +&= - \braket{ +[\bar{\chi}_{d} \Gamma_1 \chi_{h_i}](x) +[\bar{\chi}_{h_j} \Gamma_2 \chi_{d}](0) +} +\end{split} +\end{equation} + + +::: {.remark} + +1. The spinor fields $\chi_{u,d}$ and $\chi_{c,s}$ of eq. +\eqref{eq:LightFlavorSpinorFields} and +\eqref{eq:HeavyFlavorSpinorFields} +can be obtained with the projector: + +\begin{equation} +P_{i} = \frac{1}{2} [1 + (-1)^i \sigma_3] \, , \, i=0,1 +\, , +\end{equation} + +where ${\sigma_3 = + \begin{bmatrix} + 1 & 0 \\ + 0 & -1 + \end{bmatrix} + }$. + As a projector, it satisfies $P_i^2 = P_i$. + Moreover $\sigma_3 = \sigma_3^\dagger$, + hence $P_i^\dagger = P_i$. + +2. Obviously, $P_i$ commutes with the $\gamma$ matrices as they act on different spaces. + +3. The action of $\sigma_1$ swaps the up and down flavor components of a spinor. +::: + +
+ +## Wick contractions + +
+ Show Content + +Using the above remarks, +we can write: + +\begin{equation} +\begin{split} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) &= +- \braket{ +[\bar{\chi}_{d} \Gamma_1 \chi_{h_i}](x) +[\bar{\chi}_{h_j} \Gamma_2 \chi_{d}](0) +} +\\ +&= +- \braket{ +[\bar{\chi}_{\ell} \Gamma_1 (\sigma_1)^{1-i} P_i \chi_{h}](x) +[\bar{\chi}_{h} P_j (\sigma_1)^{1-j} \Gamma_2 \chi_{\ell}](0) +} +\, . +\end{split} +\end{equation} + +An implicit summation on flavor indices is understood. + +We now call $S = D^{-1}$ the inverse of the Dirac operator +and use Wick's theorem: + +\begin{equation} +\braket{\Psi_{a}(x) \bar{\Psi}_{b}(0)} = S_{ab}(x|0) \, , +\end{equation} + +where $a$ and $b$ are a shortcut for the other indices of the spinor (spin, flavor, color, etc.). +The correlator is: + +\begin{equation} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) += +\operatorname{Tr} \left[ +S_\ell(0|x) +(\sigma_1)^{1-i} P_i \Gamma_1 +S_h(x|0) +P_j (\sigma_1)^{1-j} \Gamma_2 +\right] +\end{equation} + + +We now use $\gamma_5$ hermiticity: +$(S_\ell) = \sigma_1 \gamma_5 S_\ell^\dagger \gamma_5 \sigma_1$. +Please note that the $\dagger$ is not on flavor space indices, +but here it makes no difference because $S_\ell$ is diagonal in flavor. +Our correlator becomes: + +\begin{equation} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) += +\operatorname{Tr} +\left[ + S_h(x|0) + \Gamma_2 \gamma_5 + P_j (\sigma_1)^{j} + (S_\ell)^\dagger (x|0) + (\sigma_1)^i P_i + \gamma_5 \Gamma_1 +\right] +\, . +\end{equation} + +Upon a careful calculation for all values $i,j = 0,1$ we find, equivalently (using spacetime translational symmetry): + +\begin{equation} \label{eq:C.hihj.Gamma1Gamma2} +\begin{split} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) +&= +\operatorname{Tr} +\left[ + (S_h)_{f_i f_j} (x|0) + \Gamma_2 \gamma_5 + (S_u)^\dagger (x|0) + \gamma_5 \Gamma_1 +\right] +\end{split} +\end{equation} + +This is a generalized case of eq. (A9) of @PhysRevD.59.074503. + +
+ + +## Stochastic approximation of the correlators + +### Index dilution + +
+ Show Content + + +We now make the following remark. If we want to invert numerically the system $D_{ij} \psi_{j} = \eta_{i}$ (where the $i,j$ indices include all internal indices), we have: + +\begin{equation} +D_{ij} \psi_{j} = \eta_{i} \, \implies +\psi_{i} = S_{ij} \eta_{j} \, . +\end{equation} + +If $\eta_i \eta_j^{*} = \delta_{ij}$ we have: + +\begin{equation} +S_{ij} = \psi_{i} \eta_{j}^{*} \, . +\end{equation} + +--- + +We can also use **index dilution** in order to select the components we want. In fact, if we define: $\eta_i^{(a)} = \eta \delta_{i}^{a}$, with $\eta^{*} \eta =1$ we have: + +\begin{equation} +S_{ab} = S_{ij} \delta_{i}^{a} \delta_{j}^{b} = +S_{ij} (\eta^{*} \delta_{i}^{a}) (\eta \delta_{j}^{b}) = +\psi_{i}^{(b)} (\eta_{i}^{a})^{*} += \eta^{(a)} \cdot \psi^{(b)} \, . +\end{equation} + +
+ +### Stochastic expression of the correlators + +
+ Show Content + +We now approximate the propagator using stochastic sources. +Additionally, we use: + +- The one-end-trick @Boucaud:2008xu: non vanishing source only at $x=0$. +- Same source for light and heavy doublet +- Spin dilution (cf. @foley2005practical): one source for each Dirac index $\beta$, with unit norm in color space. The components with index different from the dilution index are set to zero: + + \begin{equation} + \eta^{(\alpha)}_{\beta, c} = + \eta_c \, \delta^\alpha_\beta \,\, , \, \, + \eta_c \otimes \eta^\dagger_c = \mathbb{1}_{N_c \times N_c} \, . + \end{equation} + +- Flavor dilution: sources have an additional flavor index $\phi$, such that their flavor component different from the index vanish. The value of the flavor component however is the same, it changes only its position in the doublet: + +\begin{equation} +\eta^{(\phi=0)} = + \begin{bmatrix} + \eta \\ + 0 + \end{bmatrix} +\, , \, +\eta^{(\phi=1)} = + \begin{bmatrix} + 0 \\ + \eta + \end{bmatrix} +\, . +\end{equation} + + +Therefore, we can approximate the correlators of eq. \eqref{eq:C.hihj.Gamma1Gamma2} above with: + + +\begin{equation} +\begin{split} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) +&= +\left( + \eta^{(f_i, \alpha_1)}(0) + \cdot + \psi_h^{(f_j, \alpha_2)}(x) +\right) +(\Gamma_2 \gamma_5)_{\alpha_2 \alpha_3} +\left( + \eta^{(0, \alpha_3)}(0) + \cdot + \psi_u^{(0, \alpha_4)}(x) +\right)^{*} +(\gamma_5 \Gamma_1)_{\alpha_4 \alpha_1} +\\ +&= +(\psi_h)_{f_i, \alpha_1}^{(f_j, \alpha_2)}(x) +(\Gamma_2 \gamma_5)_{\alpha_2 \alpha_3} +(\psi_u)_{\alpha_3}^{(0, \alpha_4)}(x)^{*} +(\gamma_5 \Gamma_1)_{\alpha_4 \alpha_1} +\end{split} +\end{equation} + + +In the end we have: + +\begin{equation} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) += +[\Gamma_2 \gamma_5 (\psi_u)^{(0, \alpha_2)}(x)^{*}]_{\alpha_1} +[\gamma_5 \Gamma_1 (\psi_h)_{f_i}^{(f_j, \alpha_1)}(x) ]_{\alpha_2} +\end{equation} + +In our case $\Gamma_{1,2} = 1,\gamma_5$. Since we use the chiral basis ${\Gamma_{1,2}^{*} = (\Gamma_{1,2}^{T})^\dagger = \Gamma_{1,2}}$. +Therefore we can equivalently write the correlator as: + +\begin{equation} +\mathcal{C}^{h_i, h_j}_{\Gamma_1, \Gamma_2}(t, \vec{x}) += +[\Gamma_2 \gamma_5 (\psi_u)^{(0, \alpha_2)}(x)]^{*}_{\alpha_1} +\cdot +[\gamma_5 \Gamma_1 (\psi_h)_{f_i}^{(f_j, \alpha_1)}(x) ]_{\alpha_2} \, +\end{equation} + +where $\cdot$ is the dot product in color space (NOTE: it complex-conjugates the 1st vector). + + + +
+ + + + + + +## Implementation + +See comments in the code: `tmLQCD/meas/correlators.c`. + +**NOTE**: test if `invert_doublet_eo_quda` works (do one inversion and check the residual) + + + diff --git a/meas/correlators.c b/meas/correlators.c index de7334edf..8aa929f3c 100644 --- a/meas/correlators.c +++ b/meas/correlators.c @@ -40,6 +40,8 @@ #include "gettime.h" +#define TM_OMEAS_FILENAME_LENGTH 100 + /****************************************************** * * This routine computes the correlators @@ -52,15 +54,12 @@ * * ******************************************************/ - -#define TM_OMEAS_FILENAME_LENGTH 100 - -void correlators_measurement(const int traj, const int id, const int ieo) { +void light_correlators_measurement(const int traj, const int id, const int ieo) { tm_stopwatch_push(&g_timers, __func__, ""); int i, j, t, tt, t0; double *Cpp = NULL, *Cpa = NULL, *Cp4 = NULL; double res = 0., respa = 0., resp4 = 0.; - double atime, etime; + // double atime, etime; float tmp; operator * optr; #ifdef TM_USE_MPI @@ -204,7 +203,7 @@ void correlators_measurement(const int traj, const int id, const int ieo) { /* and write everything into a file */ if(g_mpi_time_rank == 0 && g_proc_coords[0] == 0) { ofs = fopen(filename, "w"); - fprintf( ofs, "1 1 0 %e %e\n", Cpp[t0], 0.); + fprintf( ofs, "1 1 0 %e %e\n", Cpp[t0], 0.0); for(t = 1; t < g_nproc_t*T/2; t++) { tt = (t0+t)%(g_nproc_t*T); fprintf( ofs, "1 1 %d %e ", t, Cpp[tt]); @@ -212,9 +211,9 @@ void correlators_measurement(const int traj, const int id, const int ieo) { fprintf( ofs, "%e\n", Cpp[tt]); } tt = (t0+g_nproc_t*T/2)%(g_nproc_t*T); - fprintf( ofs, "1 1 %d %e %e\n", t, Cpp[tt], 0.); + fprintf( ofs, "1 1 %d %e %e\n", t, Cpp[tt], 0.0); - fprintf( ofs, "2 1 0 %e %e\n", Cpa[t0], 0.); + fprintf( ofs, "2 1 0 %e %e\n", Cpa[t0], 0.0); for(t = 1; t < g_nproc_t*T/2; t++) { tt = (t0+t)%(g_nproc_t*T); fprintf( ofs, "2 1 %d %e ", t, Cpa[tt]); @@ -222,9 +221,9 @@ void correlators_measurement(const int traj, const int id, const int ieo) { fprintf( ofs, "%e\n", Cpa[tt]); } tt = (t0+g_nproc_t*T/2)%(g_nproc_t*T); - fprintf( ofs, "2 1 %d %e %e\n", t, Cpa[tt], 0.); + fprintf( ofs, "2 1 %d %e %e\n", t, Cpa[tt], 0.0); - fprintf( ofs, "6 1 0 %e %e\n", Cp4[t0], 0.); + fprintf( ofs, "6 1 0 %e %e\n", Cp4[t0], 0.0); for(t = 1; t < g_nproc_t*T/2; t++) { tt = (t0+t)%(g_nproc_t*T); fprintf( ofs, "6 1 %d %e ", t, Cp4[tt]); @@ -232,7 +231,7 @@ void correlators_measurement(const int traj, const int id, const int ieo) { fprintf( ofs, "%e\n", Cp4[tt]); } tt = (t0+g_nproc_t*T/2)%(g_nproc_t*T); - fprintf( ofs, "6 1 %d %e %e\n", t, Cp4[tt], 0.); + fprintf( ofs, "6 1 %d %e %e\n", t, Cp4[tt], 0.0); fclose(ofs); } #ifdef TM_USE_MPI @@ -252,3 +251,707 @@ void correlators_measurement(const int traj, const int id, const int ieo) { tm_stopwatch_pop(&g_timers, 0, 1, ""); return; } + + +void spinor_dirac_array(su3_vector* phi, spinor psi){ + phi[0] = psi.s0; + phi[1] = psi.s1; + phi[2] = psi.s2; + phi[3] = psi.s3; +} + + +/****************************************************** + * + * This routine computes the correlator matrix (), + * where G1 and G2 are 2 gamma matrices from the set {1, gamma5} --> {S,P} (scalar and pseudoscalar + *current), , , for all the choices of h1 and h2 (see eq. 18 and 20 of + *https://arxiv.org/pdf/1005.2042.pdf): + * + * C = Tr[S_l * Gamma_1 * S_h * Gamma_2] + * where S_l, S_h are the propagators of the light and heavy quark of the meson (e.g. K --> u, ) + * + * i1, i2 = operator indices from the list of the operators in the input file. + * + * + ******************************************************/ +void heavy_correlators_measurement(const int traj, const int id, const int ieo, const int i1, + const int i2) { + tm_stopwatch_push(&g_timers, __func__, ""); // timer for profiling + printf("Called: %s\n", __func__); + + spinor phi; // dummy spinor + // dummy spinors for the heavy doublet inversion + spinor* phi1 = calloc(VOLUME/2, sizeof(spinor)); + spinor* phi2 = calloc(VOLUME/2, sizeof(spinor)); + int i, j, t, tt, t0; // dummy indices + float tmp; // dummy variable + + /* light-light correlators: dummy variables */ + + // array of the values of the correlator C(t) + double *Cpp = NULL, *Cpa = NULL, *Cp4 = NULL; + // result of the accumulation of MPI partial sums + double res = 0., respa = 0., resp4 = 0.; +#ifdef TM_USE_MPI + double mpi_res = 0., mpi_respa = 0., mpi_resp4 = 0.; + // send buffer for MPI_Gather + double *sCpp = NULL, *sCpa = NULL, *sCp4 = NULL; +#endif + + /* heavy-light correlators variables */ + // sign change bilinear^\dagger, with Gamma = 1,gamma_5 + double eta_Gamma[2] = {1.0, -1.0}; + + // even-odd spinor fields for the light and heavy doublet correlators + // INDICES: source+propagator, doublet, spin dilution index, flavor projection, even-odd, + // flavor, position + // (+ Dirac, color) psi = (psi[(s,p)][db][beta][F][eo][f][x])[alpha][c] + // last 2 indices come from spinor struct + // Note: propagator in the sense that it is D^{-1}*source after the inversion + spinor* arr_eo_spinor[2][2][4][2][2][2]; + + // (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] + spinor* arr_spinor[2][2][4][2][2]; + + // allocating memory + for (size_t i_sp = 0; i_sp < 2; i_sp++) { // source or propagator + for (size_t db = 0; db < 2; db++) { // doublet: light or heavy + for (size_t beta = 0; beta < 4; beta++) { // spin dilution index + for (size_t F = 0; F < 2; F++) { // flavor dilution index + for (size_t i_eo = 0; i_eo < 2; i_eo++) { // even-odd index + for (size_t i_f = 0; i_f < 2; i_f++){// flavor index of the doublet + arr_eo_spinor[i_sp][db][beta][F][i_eo][i_f] = (spinor*) calloc(VOLUME/2, sizeof(spinor)); + zero_spinor_field(arr_eo_spinor[i_sp][db][beta][F][i_eo][i_f], VOLUME/2); + + if (i_eo == 0){ // (trick) doing it only once -> no eo index + arr_spinor[i_sp][db][beta][F][i_f] = (spinor*) calloc(VOLUME, sizeof(spinor)); + zero_spinor_field(arr_spinor[i_sp][db][beta][F][i_f], VOLUME); + } + } + } + } + } + } + } + + if (g_proc_id == 0){ + printf("Checkpoint 1\n"); + } + + if (g_proc_id == 0){ + printf("Checkpoint 2\n"); + } + + double res_hihj_g1g2[2][2][2][2]; + double mpi_res_hihj_g1g2[2][2][2][2]; + + if (g_proc_id == 0){ + printf("Checkpoint 3\n"); + } + + FILE *ofs; // output file stream + char filename[TM_OMEAS_FILENAME_LENGTH]; // file path + + /* building the stochastic propagator spinors needed for the correlators */ + + // MPI_Barrier(MPI_COMM_WORLD); + init_operators(); // initialize operators (if not already done) + + if (g_proc_id == 0){ + printf("Checkpoint 4\n"); + } + + if (no_operators < 1) { + if (g_proc_id == 0) { + // we don't want the simulation to stop, we can do the measurements offline afterwards + fprintf(stderr, + "Warning! no operators defined in input file, cannot perform online correlator " + "mesurements!\n"); + } + tm_stopwatch_pop(&g_timers, 0, 0, ""); // timer for profiling + return; + } + + // selecting the operators i1 and i2 from the list of operators initialized before + operator* optr1 = & operator_list[i1]; // light doublet + operator* optr2 = & operator_list[i2]; // heavy c,s doublet + printf("Initialized the operators\n"); + + bool b1 = (optr1->type != TMWILSON && optr1->type != WILSON && optr1->type != CLOVER); + bool b2 = (optr2->type != DBTMWILSON && optr2->type != DBCLOVER); + if (b1 || b2) { + if (g_proc_id == 0) { + if (b1) { + fprintf(stderr, + "Warning! optr1 should be one of the following: TMWILSON, WILSON and CLOVER\n"); + fprintf(stderr, "Cannot perform correlator online measurement!\n"); + } + if (b2) { + fprintf(stderr, "Warning! optr2 should be one of the following: DBTMWILSON, DBCLOVER\n"); + } + fprintf(stderr, "Cannot perform correlator online measurement!\n"); + } + tm_stopwatch_pop(&g_timers, 0, 0, ""); // timer for profiling + return; + } + + if (ranlxs_init == 0) { + rlxs_init(1, 123456); // initializing random number generator RANLUX + } + + // there are three modes of operation + // 1) one single time-slice source (default) + // 2) no_samples time-slice sources on random time-slices + // 3) one sample on all time-slices + int max_samples = measurement_list[id].all_time_slices ? 1 : measurement_list[id].no_samples; + int max_time_slices = measurement_list[id].all_time_slices ? measurement_list[id].max_source_slice : 1; + for (int sample = 0; sample < max_samples; sample++) { + for (int ts = 0; ts < max_time_slices; ts++) { + // setting output filename + if (max_samples == 1 && max_time_slices == 1) { + snprintf(filename, TM_OMEAS_FILENAME_LENGTH, "%s%06d", "heavy_mesons.", traj); + } else if (max_samples == 1 && max_time_slices > 1) { + snprintf(filename, TM_OMEAS_FILENAME_LENGTH, "%s.t%03d.%06d", "heavy_mesons", ts, traj); + } else { + snprintf(filename, TM_OMEAS_FILENAME_LENGTH, "%s.s%03d.%06d", "heavy_mesons", sample, traj); + } + + /* generate random timeslice */ + t0 = ts; + if (!measurement_list[id].all_time_slices) { + ranlxs(&tmp, 1); + t0 = (int)(measurement_list[id].max_source_slice * tmp); + } +#ifdef TM_USE_MPI + MPI_Bcast(&t0, 1, MPI_INT, 0, MPI_COMM_WORLD); // broadcast t0 to all the ranks +#endif + if (g_debug_level > 1 && g_proc_id == 0) { + printf("# timeslice set to %d (T=%d) for online measurement\n", t0, g_nproc_t * T); + printf("# online measurements parameters: kappa = %.12f, mubar = %.12f, epsbar = %.12f\n", + optr1->kappa, optr1->mubar, optr1->epsbar); + } + +#ifdef TM_USE_MPI + sCpp = (double *)calloc(T, sizeof(double)); + sCpa = (double *)calloc(T, sizeof(double)); + sCp4 = (double *)calloc(T, sizeof(double)); + if (g_mpi_time_rank == 0) { + Cpp = (double *)calloc(g_nproc_t * T, sizeof(double)); + Cpa = (double *)calloc(g_nproc_t * T, sizeof(double)); + Cp4 = (double *)calloc(g_nproc_t * T, sizeof(double)); + } +#else + Cpp = (double *)calloc(T, sizeof(double)); + Cpa = (double *)calloc(T, sizeof(double)); + Cp4 = (double *)calloc(T, sizeof(double)); +#endif + + // the number of independent correlators is 16 = (2*2)_{h_1 h_2} * (2*2)_{Gamma_1 Gamma_2} + // hi: c,s + // Gamma = 1, gamma_5 + double* C_hihj_g1g2[2][2][2][2]; + double* sC_hihj_g1g2[2][2][2][2]; + + // allocating memory + for (int h1=0; h1<2; h1++){ + for (int h2=0; h2<2; h2++){ + for (int g1=0; g1<2; g1++){ + for (int g2=0; g2<2; g2++){ +#ifdef TM_USE_MPI + mpi_res_hihj_g1g2[h1][h2][g1][g2] = 0.0; + res_hihj_g1g2[h1][h2][g1][g2] = 0.0; + sC_hihj_g1g2[h1][h2][g1][g2] = (double*) calloc(T, sizeof(double)); + if (g_mpi_time_rank == 0) { + C_hihj_g1g2[h1][h2][g1][g2] = (double*) calloc(g_nproc_t*T, sizeof(double)); + } +#else + res_hihj_g1g2[h1][h2][g1][g2] = 0.0; + C_hihj_g1g2[h1][h2][g1][g2] = (double*) calloc(T, sizeof(double)); +#endif + } + } + } + } + + /* init random sources: one for each Dirac index beta --> spin dilution */ + const unsigned int seed_i = measurement_list[id].seed; // has to be same seed + for (size_t db = 0; db < 2; db++) { // doublet: light or heavy + for (size_t beta = 0; beta < 4; beta++) { // spin dilution index + for (size_t F = 0; F < 2; F++) { // flavor dilution index + for (size_t i_f = 0; i_f < 2; i_f++) { // flavor index of the doublet + eo_source_spinor_field_spin_diluted_oet_ts( + arr_eo_spinor[0][db][beta][F][0][i_f], + arr_eo_spinor[0][db][beta][F][1][i_f], t0, + beta, sample, traj, seed_i); + } + } + } + } + zero_spinor_field(phi1, VOLUME/2); + zero_spinor_field(phi2, VOLUME/2); + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 5\n"); + } + + // heavy doublet: for each even-odd block, I multiply by (1 + i*tau_2)/sqrt(2) + // the basis for the inversion should be the same as for the light + // --> I will multiply later by the inverse, namely at the end of the inversion + for (size_t beta = 0; beta < 4; beta++) { // spin dilution index + for (size_t F = 0; F < 2; F++) { // flavor projection index + for (size_t i_eo = 0; i_eo < 2; i_eo++) { // even-odd + // (c,s) --> [(1-i*tau_2)/sqrt(2)] * (c,s) + // stored temporarely in the propagator spinors (used as dummy) + if (F==0){ + mul_one_pm_itau2_and_div_by_sqrt2( + arr_eo_spinor[1][1][beta][F][i_eo][0], arr_eo_spinor[1][1][beta][F][i_eo][1], + arr_eo_spinor[0][1][beta][F][i_eo][0], phi1, +1.0, + VOLUME / 2); + } + if (F==1){ + mul_one_pm_itau2_and_div_by_sqrt2( + arr_eo_spinor[1][1][beta][F][i_eo][0], arr_eo_spinor[1][1][beta][F][i_eo][1], + phi2, arr_eo_spinor[0][1][beta][F][i_eo][1], +1.0, + VOLUME / 2); + } + + for (size_t i_f = 0; i_f < 2; i_f++) { + // assigning the result to the first components (the sources). + // The propagators will be overwritten with the inversion + assign(arr_eo_spinor[0][1][beta][F][i_eo][i_f], + arr_eo_spinor[1][1][beta][F][i_eo][i_f], VOLUME / 2); + zero_spinor_field(arr_eo_spinor[1][1][beta][F][i_eo][i_f], VOLUME / 2); + } + } + } + } + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 6\n"); + } + + /* + assign the sources and propagator pointes for the operators + these need to be know when calling the inverter() method */ + for (size_t beta = 0; beta < 4; beta++) { // spin dilution index + /* light doublet */ + + // up + + // sources + optr1->sr0 = arr_eo_spinor[0][0][beta][0][0][0]; + optr1->sr1 = arr_eo_spinor[0][0][beta][0][1][0]; + + MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 7\n"); + } + + // propagators + optr1->prop0 = arr_eo_spinor[1][0][beta][0][0][0]; + optr1->prop1 = arr_eo_spinor[1][0][beta][0][1][0]; + + printf("inverting the light doublet\n"); + optr1->inverter(i1, 0, 0); // inversion for the up flavor + + // PLEASE KEEP THESE LINES COMMENTED, MAY BE USEFUL IN THE FUTURE + // // inversion of the light doublet only inverts the up block (the operator is + // // diagonal in flavor) down components in flavor will be empty + // optr1->DownProp = 1; + + // // sources + // optr1->sr0 = arr_eo_spinor[0][0][beta][0][1]; + // optr1->sr1 = arr_eo_spinor[0][0][beta][1][1]; + // // propagator + // optr1->prop0 = arr_eo_spinor[1][0][beta][0][1]; + // optr1->prop1 = arr_eo_spinor[1][0][beta][1][1]; + + // optr1->inverter(i1, 0, 0); // inversion for the down + + // optr1->DownProp = 0; // restoring to default + + /* heavy doublet */ + printf("Inverting the heavy doublet\n"); + for (size_t F = 0; F < 2; F++) { // flavor projection index + printf("F=%d\n", F); + + optr2->sr0 = arr_eo_spinor[0][1][beta][F][0][0]; + optr2->sr1 = arr_eo_spinor[0][1][beta][F][1][0]; + optr2->sr2 = arr_eo_spinor[0][1][beta][F][0][1]; + optr2->sr3 = arr_eo_spinor[0][1][beta][F][1][1]; + + optr2->prop0 = arr_eo_spinor[1][1][beta][F][0][0]; + optr2->prop1 = arr_eo_spinor[1][1][beta][F][1][0]; + optr2->prop2 = arr_eo_spinor[1][1][beta][F][0][1]; + optr2->prop3 = arr_eo_spinor[1][1][beta][F][1][1]; + + SourceInfo.no_flavours = 1; // only for the heavy + optr2->inverter(i2, 0, 0); // inversion for both flavor components + SourceInfo.no_flavours = 0; // reset to previous value + } + } + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 8\n"); + } + + // conclude the change of basis for the heavy doublet + for (size_t beta = 0; beta < 4; beta++) { // spin dilution index + for (size_t F = 0; F < 2; F++) { // flavor projection index + for (size_t i_eo = 0; i_eo < 2; i_eo++) { // even-odd + // (c,s) --> [(1-i*tau_2)/sqrt(2)] * (c,s) + // stored temporarely in phi1, phi2 + mul_one_pm_itau2_and_div_by_sqrt2(phi1, phi2, + arr_eo_spinor[1][1][beta][F][i_eo][0], arr_eo_spinor[1][1][beta][F][i_eo][1], -1.0, + VOLUME / 2); + + assign(arr_eo_spinor[1][1][beta][F][i_eo][0], phi1, VOLUME / 2); + assign(arr_eo_spinor[1][1][beta][F][i_eo][1], phi2, VOLUME / 2); + } + } + } + // reset to zero the dummy spinors + zero_spinor_field(phi1, VOLUME/2); + zero_spinor_field(phi2, VOLUME/2); + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 9\n"); + } + + // now we switch from even-odd representation to standard + for (size_t i_sp = 0; i_sp < 2; i_sp++) { // source or sink + for (size_t db = 0; db < 2; db++) { // doublet: light of heavy + for (size_t beta = 0; beta < 4; beta++) { // spin dilution index + for (size_t F = 0; F < 2; F++) { // flavor projection index + for (size_t i_f = 0; i_f < 2; i_f++) { // flavor index + convert_eo_to_lexic(arr_spinor[i_sp][db][beta][F][i_f], + arr_eo_spinor[i_sp][db][beta][F][0][i_f], + arr_eo_spinor[i_sp][db][beta][F][1][i_f]); + } + } + } + } + } + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 10\n"); + } + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 11.0\n"); + } + + /* + Now that I have all the propagators (all in the basis of + https://arxiv.org/pdf/1005.2042.pdf) I can build the correlators of eq. (20) + */ + const int f0 = 0; // flavor index of the up + + /* now we sum only over local space for every t */ + const int j_ts = t0-g_proc_coords[0]*T; // checkerboard index of the time at the source + + for (t = 0; t < T; t++) { + j = g_ipt[t][0][0][0]; // source and sink separated by "t" + + // dummy variables + res = 0.0; + respa = 0.0; + resp4 = 0.0; + 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++) { + res_hihj_g1g2[hi][hj][g1][g2] = 0.0; + } + } + } + } + for (i = j; i < j + LX * LY * LZ; i++) { + + // light correlators + for (size_t beta = 0; beta < 4; beta++) { // spin dilution + spinor psi_u = arr_spinor[1][0][beta][f0][f0][i]; + + res += _spinor_prod_re(psi_u, psi_u); + _gamma0(phi, psi_u); + respa += _spinor_prod_re(psi_u, phi); + _gamma5(phi, phi); + resp4 += _spinor_prod_im(psi_u, phi); + } + + // heavy correlators + 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 complex dum_tot = 0.0; + for (int alpha_1 = 0; alpha_1 < 4; alpha_1++){ + spinor psi_h = arr_spinor[1][1][alpha_1][hj][hi][i]; + if (g1 == 0){ // Gamma_1 = Id + _gamma5(psi_h, psi_h); + } + su3_vector psi_h_su3[4]; + spinor_dirac_array(&psi_h_su3[0], psi_h); + for (int alpha_2 = 0; alpha_2 < 4; alpha_2++){ + spinor psi_u_star = arr_spinor[1][0][alpha_2][f0][f0][i]; + if (g2 == 0){ // Gamma_2 = Id. NOTE: works because Gamma_2=Gamma_2* for Gamma_2=1,gamma_5 + _gamma5(psi_u_star, psi_u_star); + } + su3_vector psi_u_star_su3[4]; + spinor_dirac_array(&psi_u_star_su3[0], psi_u_star); + complex double dum_12 = 0.0; + _colorvec_scalar_prod(dum_12, psi_u_star_su3[alpha_1], psi_h_su3[alpha_2]); + dum_tot += dum_12; + } + } + // if (g_proc_id == 0){ + // printf("dum_tot = %d %d %d %d %d %e %e \n", t, hi, hj, g1, g2, creal(dum_tot), cimag(dum_tot)); + // } + if (hi != hj){ + dum_tot *= -1; + } + if (g1==g2){ + res_hihj_g1g2[hi][hj][g1][g2] += creal(dum_tot); // correlators is real + } + else{ + if (g2==1){ + dum_tot *= -1; + } + res_hihj_g1g2[hi][hj][g1][g2] += cimag(dum_tot); // correlator is imaginary + } + } + } + } + + } + } + + const double vol_fact = (g_nproc_x * LX) * (g_nproc_y * LY) * (g_nproc_z * LZ); +#if defined TM_USE_MPI + MPI_Reduce(&res, &mpi_res, 1, MPI_DOUBLE, MPI_SUM, 0, g_mpi_time_slices); + res = mpi_res; + MPI_Reduce(&respa, &mpi_respa, 1, MPI_DOUBLE, MPI_SUM, 0, g_mpi_time_slices); + respa = mpi_respa; + MPI_Reduce(&resp4, &mpi_resp4, 1, MPI_DOUBLE, MPI_SUM, 0, g_mpi_time_slices); + resp4 = mpi_resp4; + + sCpp[t] = +res / vol_fact / 2.0 / optr1->kappa / optr1->kappa; + sCpa[t] = -respa / vol_fact / 2.0 / optr1->kappa / optr1->kappa; + sCp4[t] = +resp4 / vol_fact / 2.0 / optr1->kappa / optr1->kappa; +#else + Cpp[t] = +res / vol_fact / 2.0 / optr1->kappa / optr1->kappa; + Cpa[t] = -respa / vol_fact / 2.0 / optr1->kappa / optr1->kappa; + Cp4[t] = +resp4 / vol_fact / 2.0 / optr1->kappa / optr1->kappa; +#endif + + + 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++) { +#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]; + + sC_hihj_g1g2[hi][hj][g1][g2][t] = - eta_Gamma[g1] * res_hihj_g1g2[hi][hj][g1][g2] / vol_fact; +#else + C_hihj_g1g2[hi][hj][g1][g2][t] = - eta_Gamma[g1] * res_hihj_g1g2[hi][hj][g1][g2] / vol_fact; +#endif + } + } + } + } + + } // end 1st loop over "t" + + + if (g_proc_id == 0){ + printf("Checkpoint 11.6\n"); + } + +#ifdef TM_USE_MPI + /* some gymnastics needed in case of parallelisation */ + if (g_mpi_time_rank == 0) { + // light correlators + MPI_Gather(sCpp, T, MPI_DOUBLE, Cpp, T, MPI_DOUBLE, 0, g_mpi_SV_slices); + MPI_Gather(sCpa, T, MPI_DOUBLE, Cpa, T, MPI_DOUBLE, 0, g_mpi_SV_slices); + MPI_Gather(sCp4, T, MPI_DOUBLE, Cp4, T, MPI_DOUBLE, 0, g_mpi_SV_slices); + + // heavy mesons + 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++) { + MPI_Gather(sC_hihj_g1g2[hi][hj][g1][g2], T, MPI_DOUBLE, C_hihj_g1g2[hi][hj][g1][g2], T, MPI_DOUBLE, 0, g_mpi_SV_slices); + } + } + } + } + + } +#endif + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 12\n"); + } + + printf("check: %d, %d\n", g_mpi_time_rank, g_proc_coords[0]); + + /* and write everything into a file */ + if (g_mpi_time_rank == 0 && g_proc_coords[0] == 0) { + printf("Checkpoint 12.1 : -%s-\n", filename); + ofs = fopen(filename, "w"); + printf("Checkpoint 12.2\n"); + + + fprintf(ofs, "1 1 0 %e %e\n", Cpp[t0], 0.0); + for (t = 1; t < g_nproc_t * T / 2; t++) { + tt = (t0 + t) % (g_nproc_t * T); + fprintf(ofs, "1 1 %d %e ", t, Cpp[tt]); + tt = (t0 + g_nproc_t * T - t) % (g_nproc_t * T); + fprintf(ofs, "%e\n", Cpp[tt]); + } + tt = (t0 + g_nproc_t * T / 2) % (g_nproc_t * T); + fprintf(ofs, "1 1 %d %e %e\n", t, Cpp[tt], 0.0); + + fprintf(ofs, "2 1 0 %e %e\n", Cpa[t0], 0.0); + for (t = 1; t < g_nproc_t * T / 2; t++) { + tt = (t0 + t) % (g_nproc_t * T); + fprintf(ofs, "2 1 %d %e ", t, Cpa[tt]); + tt = (t0 + g_nproc_t * T - t) % (g_nproc_t * T); + fprintf(ofs, "%e\n", Cpa[tt]); + } + tt = (t0 + g_nproc_t * T / 2) % (g_nproc_t * T); + fprintf(ofs, "2 1 %d %e %e\n", t, Cpa[tt], 0.0); + + fprintf(ofs, "6 1 0 %e %e\n", Cp4[t0], 0.0); + for (t = 1; t < g_nproc_t * T / 2; t++) { + tt = (t0 + t) % (g_nproc_t * T); + fprintf(ofs, "6 1 %d %e ", t, Cp4[tt]); + tt = (t0 + g_nproc_t * T - t) % (g_nproc_t * T); + fprintf(ofs, "%e\n", Cp4[tt]); + } + tt = (t0 + g_nproc_t * T / 2) % (g_nproc_t * T); + fprintf(ofs, "6 1 %d %e %e\n", t, Cp4[tt], 0.0); + + printf("Checkpoint 12.3\n"); + // heavy mesons correlators + 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++) { + printf("check matrix: %e\n", C_hihj_g1g2[hi][hj][g1][g2][3]); + fprintf(ofs, "%u %u %u %u 0 %e %e\n", hi, hj, g1, g2, C_hihj_g1g2[hi][hj][g1][g2][t0], 0.0); + for (t = 1; t < g_nproc_t * T / 2; t++) { + tt = (t0 + t) % (g_nproc_t * T); + fprintf(ofs, "%u %u %u %u %d %e ", hi, hj, g1, g2, t, C_hihj_g1g2[hi][hj][g1][g2][tt]); + tt = (t0 + g_nproc_t * T - t) % (g_nproc_t * T); + fprintf(ofs, "%e\n", C_hihj_g1g2[hi][hj][g1][g2][tt]); + } + tt = (t0 + g_nproc_t * T / 2) % (g_nproc_t * T); + fprintf(ofs, "%u %u %u %u %d %e %e\n", hi, hj, g1, g2, t, C_hihj_g1g2[hi][hj][g1][g2][tt], 0.0); + } + } + } + } + fclose(ofs); + } + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 13\n"); + } + + // freeing memory: light correlators +#ifdef TM_USE_MPI + if (g_mpi_time_rank == 0) { + free(Cpp); + free(Cpa); + free(Cp4); + } + free(sCpp); + free(sCpa); + free(sCp4); +#else + free(Cpp); + free(Cpa); + free(Cp4); +#endif + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 14\n"); + } + + + // freeing memory: heavy-light correlators + for (int h1=0; h1<2; h1++){ + for (int h2=0; h2<2; h2++){ + for (int g1=0; g1<2; g1++){ + for (int g2=0; g2<2; g2++){ +#ifdef TM_USE_MPI + free(sC_hihj_g1g2[h1][h2][g1][g2]); + if (g_mpi_time_rank == 0) { + free(C_hihj_g1g2[h1][h2][g1][g2]); + } +#else + free(C_hihj_g1g2[h1][h2][g1][g2]); +#endif + } + } + } + } + + //MPI_Barrier(MPI_COMM_WORLD); + if (g_proc_id == 0){ + printf("Checkpoint 15\n"); + } + + } // for(max_time_slices) + } // for(max_samples) + + // freeing up memory + for (size_t i_sp = 0; i_sp < 2; i_sp++) { // source or propagator + for (size_t db = 0; db < 2; db++) { // doublet: light or heavy + for (size_t beta = 0; beta < 4; beta++) { // spin dilution index + for (size_t F = 0; F < 2; F++) { // flavor dilution index + for (size_t i_eo = 0; i_eo < 2; i_eo++) { // even-odd index + for (size_t i_f = 0; i_f < 2; i_f++){// flavor index of the doublet + free(arr_eo_spinor[i_sp][db][beta][F][i_eo][i_f]); + if(i_eo == 0){ // only once: no "eo" index for arr_spinor + free(arr_spinor[i_sp][db][beta][F][i_f]); + } + } + } + } + } + } + } + + free(phi1); + free(phi2); + + tm_stopwatch_pop(&g_timers, 0, 1, ""); + + return; +} + +void correlators_measurement(const int traj, const int id, const int ieo) { + + // ??? maybe add a double check? does i1 --> light and i2 --> heavy? + if (measurement_list[id].measure_heavy_mesons == 1) { + const int i1 = 0, i2 = 1; + heavy_correlators_measurement(traj, id, ieo, i1, i2); + } + + light_correlators_measurement(traj, id, ieo); +} diff --git a/meas/correlators.h b/meas/correlators.h index c9a1c4ac0..3f817734f 100644 --- a/meas/correlators.h +++ b/meas/correlators.h @@ -21,6 +21,24 @@ #ifndef _ONLINE_MEASUREMENT_H #define _ONLINE_MEASUREMENT_H -void correlators_measurement(const int traj, const int t0, const int ieo); +#include + +#include "init/init_spinor_field.h" +#include "linalg/assign.h" +#include "operator/tm_operators.h" + + +/* measurement of the correlators involving the light doublet (see tmLQCD documentation)*/ +void light_correlators_measurement(const int traj, const int id, const int ieo); + +/* measurement of the 4x4 matrix for the heavy doublet: eq. 20 of https://arxiv.org/pdf/1005.2042.pdf */ +void heavy_correlators_measurement(const int traj, const int id, const int ieo, const int i1, const int i2); + +/* + Function that is called when the correlator measurement is specified in the input file. + Internally, it calls the functions for the measure of the light and (optionally) the heavy mesons correlators + */ +void correlators_measurement(const int traj, const int id, const int ieo); + #endif diff --git a/meas/measurements.c b/meas/measurements.c index b18d263f9..1b9eca94b 100644 --- a/meas/measurements.c +++ b/meas/measurements.c @@ -43,7 +43,7 @@ int no_measurements = 0; int add_measurement(const enum MEAS_TYPE meas_type) { if(no_measurements == max_no_measurements) { - fprintf(stderr, "maximal number of measurementss %d exceeded!\n", max_no_measurements); + fprintf(stderr, "maximal number of measurements %d exceeded!\n", max_no_measurements); exit(-1); } measurement_list[no_measurements].measurefunc = &dummy_meas; diff --git a/meas/measurements.h b/meas/measurements.h index 0f5d6d6ab..f795e7e04 100644 --- a/meas/measurements.h +++ b/meas/measurements.h @@ -68,8 +68,12 @@ typedef struct { /* functions for the measurement */ void (*measurefunc) (const int traj, const int id, const int ieo); - + ExternalLibrary external_library; + + // measuring the 4x4 matrix of correlators in the 1+1 sector: https://arxiv.org/abs/1005.2042 + int measure_heavy_mesons; // 1 or 0: heavy mesons measured or not + } measurement; diff --git a/offline_measurement.c b/offline_measurement.c index f9d4b1215..4911ecd14 100644 --- a/offline_measurement.c +++ b/offline_measurement.c @@ -75,7 +75,7 @@ static void set_default_filenames(char ** input_filename, char ** filename); int main(int argc, char *argv[]) { FILE *parameterfile = NULL; - int j, i, ix = 0, isample = 0, op_id = 0; + int j, i; //, ix = 0, isample = 0, op_id = 0; char datafilename[206]; char parameterfilename[206]; char conf_filename[CONF_FILENAME_LENGTH]; @@ -83,6 +83,7 @@ int main(int argc, char *argv[]) char * filename = NULL; double plaquette_energy; + init_critical_globals(TM_PROGRAM_OFFLINE_MEASUREMENT); #ifdef _KOJAK_INST diff --git a/operator.c b/operator.c index 8ceb8d4dd..154e6eaea 100644 --- a/operator.c +++ b/operator.c @@ -477,9 +477,9 @@ void op_invert(const int op_id, const int index_start, const int write_prop) { /* this requires multiplication of source with */ /* (1+itau_2)/sqrt(2) and the result with (1-itau_2)/sqrt(2) */ - mul_one_pm_itau2(optr->prop0, optr->prop2, g_spinor_field[DUM_DERI], + mul_one_pm_itau2_and_div_by_sqrt2(optr->prop0, optr->prop2, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+2], -1., VOLUME/2); - mul_one_pm_itau2(optr->prop1, optr->prop3, g_spinor_field[DUM_DERI+1], + mul_one_pm_itau2_and_div_by_sqrt2(optr->prop1, optr->prop3, g_spinor_field[DUM_DERI+1], g_spinor_field[DUM_DERI+3], -1., VOLUME/2); /* write propagator */ if(write_prop) optr->write_prop(op_id, index_start, i); @@ -495,11 +495,11 @@ void op_invert(const int op_id, const int index_start, const int write_prop) { fprintf(stdout, "# Inversion done in %d iterations, squared residue = %e!\n", optr->iterations, optr->reached_prec); } - mul_one_pm_itau2(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+2], optr->sr0, optr->sr2, -1., VOLUME/2); - mul_one_pm_itau2(g_spinor_field[DUM_DERI+1], g_spinor_field[DUM_DERI+3], optr->sr1, optr->sr3, -1., VOLUME/2); + mul_one_pm_itau2_and_div_by_sqrt2(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+2], optr->sr0, optr->sr2, -1., VOLUME/2); + mul_one_pm_itau2_and_div_by_sqrt2(g_spinor_field[DUM_DERI+1], g_spinor_field[DUM_DERI+3], optr->sr1, optr->sr3, -1., VOLUME/2); - mul_one_pm_itau2(optr->sr0, optr->sr2, g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI], +1., VOLUME/2); - mul_one_pm_itau2(optr->sr1, optr->sr3, g_spinor_field[DUM_DERI+3], g_spinor_field[DUM_DERI+1], +1., VOLUME/2); + mul_one_pm_itau2_and_div_by_sqrt2(optr->sr0, optr->sr2, g_spinor_field[DUM_DERI+2], g_spinor_field[DUM_DERI], +1., VOLUME/2); + mul_one_pm_itau2_and_div_by_sqrt2(optr->sr1, optr->sr3, g_spinor_field[DUM_DERI+3], g_spinor_field[DUM_DERI+1], +1., VOLUME/2); } /* volume sources need only one inversion */ else if(SourceInfo.type == SRC_TYPE_VOL) i++; diff --git a/operator/tm_operators_nd.c b/operator/tm_operators_nd.c index 5cdc4f385..41729dc36 100644 --- a/operator/tm_operators_nd.c +++ b/operator/tm_operators_nd.c @@ -722,7 +722,7 @@ void Q_test_epsilon(spinor * const l_strange, spinor * const l_charm, } -void mul_one_pm_itau2(spinor * const p, spinor * const q, +void mul_one_pm_itau2_and_div_by_sqrt2(spinor * const p, spinor * const q, spinor * const r, spinor * const s, const double sign, const int N) { double fac = 1./sqrt(2.); diff --git a/operator/tm_operators_nd.h b/operator/tm_operators_nd.h index 138e9b93b..35ff75d67 100644 --- a/operator/tm_operators_nd.h +++ b/operator/tm_operators_nd.h @@ -22,7 +22,8 @@ #ifndef _TM_OPERATTORS_ND_H #define _TM_OPERATTORS_ND_H -void mul_one_pm_itau2(spinor * const p, spinor * const q, +// \frac{1}{\sqrt(2)} (1 \pm \tau_2) +void mul_one_pm_itau2_and_div_by_sqrt2(spinor * const p, spinor * const q, spinor * const r, spinor * const s, const double sign, const int N); diff --git a/prepare_source.c b/prepare_source.c index b4d726ba3..9974df558 100644 --- a/prepare_source.c +++ b/prepare_source.c @@ -311,8 +311,8 @@ void prepare_source(const int nstore, const int isample, const int ix, const int gaussian_volume_source(g_spinor_field[2], g_spinor_field[3], isample, nstore, 2); } - mul_one_pm_itau2(g_spinor_field[4], g_spinor_field[6], g_spinor_field[0], g_spinor_field[2], +1., VOLUME/2); - mul_one_pm_itau2(g_spinor_field[5], g_spinor_field[7], g_spinor_field[1], g_spinor_field[3], +1., VOLUME/2); + mul_one_pm_itau2_and_div_by_sqrt2(g_spinor_field[4], g_spinor_field[6], g_spinor_field[0], g_spinor_field[2], +1., VOLUME/2); + mul_one_pm_itau2_and_div_by_sqrt2(g_spinor_field[5], g_spinor_field[7], g_spinor_field[1], g_spinor_field[3], +1., VOLUME/2); assign(g_spinor_field[0], g_spinor_field[4], VOLUME/2); assign(g_spinor_field[1], g_spinor_field[5], VOLUME/2); assign(g_spinor_field[2], g_spinor_field[6], VOLUME/2); diff --git a/read_input.l b/read_input.l index 3b39524dc..a4972beff 100644 --- a/read_input.l +++ b/read_input.l @@ -3106,6 +3106,14 @@ static inline double fltlist_next_token(int * const list_end){ meas->external_library = QUDA_LIB; if(myverbose) printf(" UseExternalLibrary set to quda in line %d, monomial %d\n", line_of_file, current_monomial); } + {SPC}*HeavyMesons{EQL}yes { + meas->measure_heavy_mesons = 1; + if(myverbose) printf(" HeavyMesons set to yes in line %d, monomial %d\n", line_of_file, current_monomial); + } + {SPC}*HeavyMesons{EQL}no { + meas->measure_heavy_mesons = 0; + if(myverbose) printf(" HeavyMesons set to no in line %d, monomial %d\n", line_of_file, current_monomial); + } } { diff --git a/solver/dirac_operator_eigenvectors.h b/solver/dirac_operator_eigenvectors.h index b4ef212e1..6b7454533 100644 --- a/solver/dirac_operator_eigenvectors.h +++ b/solver/dirac_operator_eigenvectors.h @@ -265,6 +265,14 @@ int cyclicDiff(int a,int b, int period); conj((r).s3.c2) * (s).s3.c2; + +#define _colorvec_scalar_prod(proj,r,s)\ + (proj) = conj((r).c0) * (s).c0 + \ + conj((r).c1) * (s).c1; + \ + conj((r).c2) * (s).c2; + + + #define PROJECTSPLIT(p_plus,up_plus,col_proj,phi_o,phi_plus,col_phi)\ p_plus = 0; \ p_plus += conj(up_plus->s0.col_proj) * (phi_o->s0.col_phi); \ diff --git a/su3.h b/su3.h index b66652fe5..d8f238161 100644 --- a/su3.h +++ b/su3.h @@ -47,6 +47,7 @@ typedef struct _Complex float c00, c01, c02, c10, c11, c12, c20, c21, c22; } su3_32; +// vector in color space typedef struct { _Complex double c0,c1,c2; @@ -57,6 +58,7 @@ typedef struct _Complex float c0,c1,c2; } su3_vector32; +// 4 Dirac components of the (colored) spinor (vector in color space) typedef struct { su3_vector s0,s1,s2,s3;