Spn smearing

.gitignore

@@ -28,3 +28,4 @@ Hirep-git.code-workspace
 **/clangd
 .ycm_extra_conf.py*
 compile_commands.json
+Spectrum/measure_spectrum

Include/global.h

@@ -143,6 +143,15 @@ GLB_VAR(suNg_field_flt,*u_gauge_flt,=NULL);
 GLB_VAR(suNf_field,*u_gauge_f,=NULL);
 GLB_VAR(suNg_field,*u_gauge_s,=NULL);
 GLB_VAR(suNf_field_flt,*u_gauge_f_flt,=NULL);
+
+/*APE smearing*/
+GLB_VAR(suNg_field,*u_gauge_APE,=NULL);
+GLB_VAR(suNg_scalar_field,*u_scalar_APE,=NULL);
+GLB_VAR(suNg_field_flt,*u_gauge_APE_flt,=NULL);
+GLB_VAR(suNf_field,*u_gauge_APE_f,=NULL);
+GLB_VAR(suNg_field,*u_gauge_APE_s,=NULL);
+GLB_VAR(suNf_field_flt,*u_gauge_APE_f_flt,=NULL);
+
 #if defined(GAUGE_SPN) && defined(REPR_FUNDAMENTAL)
 GLB_VAR(suNffull_field,*cl_term,=NULL);
 GLB_VAR(suNffull_field,*cl_force,=NULL);
@@ -161,6 +170,14 @@ GLB_VAR(int,gauge_field_active,=0); // whether gauge field interactions is activ
 #define pu_gauge_f(ix,mu) ((u_gauge_f->ptr)+coord_to_index(ix,mu))
 #define pu_gauge_f_flt(ix,mu) ((u_gauge_f_flt->ptr)+coord_to_index(ix,mu))
 
+/*APE smearing*/
+#define pu_gauge_APE(ix,mu) ((u_gauge_APE->ptr)+coord_to_index(ix,mu))
+#define pu_scalar_APE(ix) ((u_scalar_APE->ptr)+ix)
+#define pu_gauge_APE_flt(ix,mu) ((u_gauge_APE_flt->ptr)+coord_to_index(ix,mu))
+#define pu_gauge_APE_f(ix,mu) ((u_gauge_APE_f->ptr)+coord_to_index(ix,mu))
+#define pu_gauge_APE_f_flt(ix,mu) ((u_gauge_APE_f_flt->ptr)+coord_to_index(ix,mu))
+#define pu_gauge_tmp(ix,mu) ((u_gauge_tmp->ptr)+coord_to_index(ix,mu))
+
 /* input parameters */
 #include "input_par.h"
 GLB_VAR(input_glb,glb_var,=init_input_glb(glb_var));

Include/observables.h

@@ -296,4 +296,32 @@ void init_triplet_discon_correlators();
 
 void ff_observables();
 
+//APE smearing
+double plaq_APE(int ix,int mu,int nu);
+void cplaq_APE(complex *ret,int ix,int mu,int nu);
+double avr_spatial_plaquette_APE();
+void full_plaquette_APE();
+void polyakov_APE();
+
+/* source smearing*/
+void smearing_function(spinor_field *source, int tau, double epsilon);
+void smearing_function_with_APE(spinor_field *source, int tau, double epsilon);
+void smearing_function_volume(spinor_field *source, double epsilon);
+void smearing_function_volume_with_APE(spinor_field *source, double epsilon);
+void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear);
+void create_smeared_source_with_APE(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear);
+void create_noise_source_equal_eo_smeared_with_APE(spinor_field *source, double epsilon, int Nsmear);
+void create_noise_source_equal_eo_smeared(spinor_field *source, double epsilon, int Nsmear);
+
+/* Sink smearing*/
+void smeared_propagator(spinor_field* psi, int nm , double epsilon);
+void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon);
+void smeared_propagator_volume(spinor_field* psi, int nm , double epsilon);
+void smeared_propagator_volume_with_APE(spinor_field* psi, int nm , double epsilon);
+
+/* Measurements*/
+void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source_max, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff);
+void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff);
+
+
 #endif

Include/representation.h

@@ -21,4 +21,6 @@ void _group_represent2_flt(suNf_flt* v, suNg_flt *u);
 void represent_gauge_field();
 void represent_gauge_field_measure();
 
+/* APE smearing*/
+void represent_gauge_field_APE();
 #endif

Include/spectrum.h

@@ -19,6 +19,8 @@
 
 void measure_spectrum_semwall(int nm, double* m, int nhits,int conf_num, double precision,storage_switch swc, data_storage_array **ret);
 void measure_spectrum_discon_semwall(int nm, double* m, int nhits,int conf_num, double precision,storage_switch swc, data_storage_array **ret);
+void measure_spectrum_discon_semwall_smeared(int nm, double* m, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double APE_epsilon, int APE_N, int N_diff);
+void measure_spectrum_discon_semwall_smeared_single_inversion(int nm, double* m, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double APE_epsilon, int APE_N, int N_diff);
 void measure_spectrum_discon_gfwall(int nm, double* m, int conf_num, double precision,storage_switch swc, data_storage_array **ret);
 void measure_spectrum_discon_volume(int nm, double* m, int conf_num, double precision, int dil,storage_switch swc, data_storage_array **ret);
 void measure_spectrum_gfwall(int nm, double* m, int conf_num, double precision,storage_switch swc, data_storage_array **ret);

Include/utils.h

@@ -44,6 +44,7 @@ void init_plaq_open_BCs(double *plaq_weight, double *rect_weight, double ct, dou
 
 void free_BCs();
 void apply_BCs_on_represented_gauge_field();
+void apply_BCs_on_represented_gauge_field_APE();
 void apply_BCs_on_fundamental_gauge_field();
 void apply_BCs_on_momentum_field(suNg_av_field *force);
 void apply_BCs_on_spinor_field(spinor_field *sp);
@@ -176,4 +177,15 @@ int reserve_wrk_space_with_pointers(suNg_field **g_wrk_out, int **i_up_wrk_out,
 void release_wrk_space(int id_release);
 void free_wrk_space();
 
+/* APE smearing */
+void APE_smearing(double smear_val, int Nsmear);
+void assign_ud2u_APE_f(void);
+#ifdef GAUGE_SPN
+void cooling_SPN(suNg* g_out, suNg* g_in, suNg* staple, int cooling);
+void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11);
+void vmxsu2(int i1, int i2, suNgfull* A, complex B[4]);
+void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11);
+void vmxsu2_tau(int n1, int n2, suNgfull* A, complex B[4]);
+#endif
+
 #endif

LibHR/Observables/avr_plaquette.c

@@ -416,3 +416,214 @@ double complex avr_plaquette_wrk()
 
   return pa;
 }
+
+
+double plaq_APE(int ix,int mu,int nu){
+  /**
+   * @brief Computes the value of the plaquette at the given position in the mu-nu plane and return only the real part.
+   *        Modified from `plaq(int ix, int mu, int nu)` by changing `pu_gauge` to `pu_gauge_APE`
+   * 
+   * @param ix The position of the starting site.
+   * @param mu The first direction of the plaquette.
+   * @param nu The second direction of the plaquette.
+   * 
+   * @return The real part of the computed plaquette value, possibly weighted if PLAQ_WEIGHTS is defined.
+   */
+  int iy,iz;
+  double p;
+  suNg *v1,*v2,*v3,*v4,w1,w2,w3;
+
+  iy=iup(ix,mu);
+  iz=iup(ix,nu);
+
+  v1=pu_gauge_APE(ix,mu);
+  v2=pu_gauge_APE(iy,nu);
+  v3=pu_gauge_APE(iz,mu);
+  v4=pu_gauge_APE(ix,nu);
+
+  _suNg_times_suNg(w1,(*v1),(*v2));
+  _suNg_times_suNg(w2,(*v4),(*v3));
+  _suNg_times_suNg_dagger(w3,w1,w2);
+
+  _suNg_trace_re(p,w3);
+
+#ifdef PLAQ_WEIGHTS
+  if(plaq_weight==NULL) return p;
+  return plaq_weight[ix*16+mu*4+nu]*p;
+#else
+  return p;
+#endif
+}
+
+
+void cplaq_APE(complex *ret,int ix,int mu,int nu){
+  /**
+   * @brief Computes the value of the plaquette at the given position in the mu-nu plane.
+   *        Modified from `cplaq(complex *ret,int ix,int mu,int nu)` by changing `pu_gauge` to `pu_gauge_APE`
+   * 
+   * @param ret A pointer to store the resulting complex plaquette value.
+   * @param ix The position of the starting site.
+   * @param mu The first direction of the plaquette.
+   * @param nu The second direction of the plaquette.
+   * 
+   * @note If PLAQ_WEIGHTS is defined, the result may include a weighting factor.
+   */
+  int iy,iz;
+  suNg *v1,*v2,*v3,*v4,w1,w2,w3;
+  double tmpre = 0.;
+  double tmpim = 0.;
+
+  iy=iup(ix,mu);
+  iz=iup(ix,nu);
+
+  v1=pu_gauge_APE(ix,mu);
+  v2=pu_gauge_APE(iy,nu);
+  v3=pu_gauge_APE(iz,mu);
+  v4=pu_gauge_APE(ix,nu);
+
+  _suNg_times_suNg(w1,(*v1),(*v2));
+  _suNg_times_suNg(w2,(*v4),(*v3));
+  _suNg_times_suNg_dagger(w3,w1,w2);
+
+  _suNg_trace_re(tmpre, w3);
+#ifdef GAUGE_SON
+  _suNg_trace_im(tmpim, w3);
+#else
+  *ret = tmpre + I * tmpim;
+#endif
+
+#ifdef PLAQ_WEIGHTS
+  if(plaq_weight!=NULL) {
+    *ret *= plaq_weight[ix * 16 + mu * 4 + nu];;
+  }
+#endif
+
+}
+
+
+double avr_spatial_plaquette_APE(){
+  /**
+   * @brief Computes the average (over the whole lattice) plaquette value over all spatial directions 
+   *        for the APE-smeared gauge links.
+   *        Modified from `avr_spatial_plaquette()` by changing `plaq()` to `plaq_APE()`
+   * 
+   */
+  double pa=0.;
+
+  _PIECE_FOR(&glattice,ixp) {
+    if(ixp==glattice.inner_master_pieces) {
+      _OMP_PRAGMA( master )
+      /* wait for gauge field to be transfered */
+      complete_gf_sendrecv(u_gauge_APE);
+      _OMP_PRAGMA( barrier )
+    }
+    _SITE_FOR_SUM(&glattice,ixp,ix,pa) {
+
+      pa+=plaq_APE(ix,2,1);
+      pa+=plaq_APE(ix,3,1);
+      pa+=plaq_APE(ix,3,2);
+    }
+  }
+
+  global_sum(&pa, 1);
+
+#ifdef BC_T_OPEN
+    pa /= 3.0*NG*GLB_VOLUME*(GLB_T-1)/GLB_T;
+#else
+    pa /= 3.0*NG*GLB_VOLUME;
+#endif
+
+  return pa;
+
+}
+
+void full_plaquette_APE(){
+  /**
+   * @brief Computes and prints the average plaquette value over all directions, 
+   *        including the temporal direction, for the APE-smeared gauge links.
+   *        Modified from `ull_plaquette())` by changing `cplaq(...)` to `cplaq_APE(...)`
+   * 
+   * @note APE smearing is not applied in the temporal direction. 
+   *       As a result, plaquette values involving the temporal direction are zero.
+   */
+
+    static double complex pa[6];
+    static double complex r0;
+    static double complex r1;
+    static double complex r2;
+    static double complex r3;
+    static double complex r4;
+    static double complex r5;
+
+    _OMP_PRAGMA(single)
+    {
+      r0 = 0.;
+      r1 = 0.;
+      r2 = 0.;
+      r3 = 0.;
+      r4 = 0.;
+      r5 = 0.;
+    }
+
+    _PIECE_FOR(&glattice,ixp)
+    {
+        if(ixp == glattice.inner_master_pieces)
+        {
+            _OMP_PRAGMA( master )
+            /* wait for gauge field to be transfered */
+            complete_gf_sendrecv(u_gauge_APE);
+            _OMP_PRAGMA( barrier )
+        }
+
+        _SITE_FOR_SUM(&glattice,ixp,ix,r0re,r0im,r1re,r1im,r2re,r2im,r3re,r3im,r4re,r4im,r5re,r5im)
+        {
+            complex tmp;
+
+            cplaq_APE(&tmp,ix,1,0);
+            r0 += tmp;
+
+            cplaq_APE(&tmp,ix,2,0);
+            r1 += tmp;
+
+            cplaq_APE(&tmp,ix,2,1);
+            r2 += tmp;
+
+            cplaq_APE(&tmp,ix,3,0);
+            r3 += tmp;
+
+            cplaq_APE(&tmp,ix,3,1);
+            r4 += tmp;
+
+            cplaq_APE(&tmp,ix,3,2);
+            r5 += tmp;
+
+        }
+    }
+
+    _OMP_PRAGMA(single)
+    {
+      pa[0] = r0;
+      pa[1] = r1;
+      pa[2] = r2;
+      pa[3] = r3;
+      pa[4] = r4;
+      pa[5] = r5;
+    }
+
+    global_sum((double*)pa,12);
+    for(int k = 0; k < 6; k++)
+    {
+#ifdef BC_T_OPEN
+        pa[k] /= NG*GLB_VOLUME*(GLB_T-1)/GLB_T;
+#else
+        pa[k] /= NG*GLB_VOLUME;
+#endif
+    }
+
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 1, 0, creal(pa[0]), cimag(pa[0]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 2, 0, creal(pa[1]), cimag(pa[1]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 2, 1, creal(pa[2]), cimag(pa[2]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 3, 0, creal(pa[3]), cimag(pa[3]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 3, 1, creal(pa[4]), cimag(pa[4]));
+    lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 3, 2, creal(pa[5]), cimag(pa[5]));
+}