introduced a cloverdet monomial which is running.

It is also giving correct results for C_sw = 0

Makefile.in

@@ -63,7 +63,7 @@ MODULES = read_input gamma hybrid_update measure_gauge_action start \
 	xchange_lexicfield xchange_2fields online_measurement \
 	monomial det_monomial detratio_monomial update_momenta \
 	integrator gauge_monomial ndpoly_monomial phmc \
-	clover_trlog_monomial \
+	clover_trlog_monomial cloverdet_monomial \
 	little_D block Dov_psi operator poly_monomial measurements pion_norm Dov_proj \
 	xchange_field_tslice temporalgauge spinor_fft X_psi P_M_eta \
 	xchange_jacobi jacobi init_jacobi_field

clover.c

@@ -71,7 +71,7 @@ void Qsw_psi(spinor * const l, spinor * const k) {
   clover_inv(EO, g_spinor_field[DUM_MATRIX+1]);
   Hopping_Matrix(OE, g_spinor_field[DUM_MATRIX], g_spinor_field[DUM_MATRIX+1]);
   /* temporarily set to 0 -> needs to be fixed */
-  mul_one_pm_imu_sub_mul_gamma5(l, k, g_spinor_field[DUM_MATRIX], +1.);
+  /*   mul_one_pm_imu_sub_mul_gamma5(l, k, g_spinor_field[DUM_MATRIX], +1.); */
   clover_gamma5(OE, l, k, g_spinor_field[DUM_MATRIX], 0.);
 }
 
@@ -96,6 +96,13 @@ void Msw_psi(spinor * const l, spinor * const k) {
   clover(OE, l, k, g_spinor_field[DUM_MATRIX], 0.);
 }
 
+void H_eo_sw_inv_psi(spinor * const l, spinor * const k, const int ieo) {
+  Hopping_Matrix(ieo, l, k);
+  clover_inv(ieo, l);
+  return;
+}
+
+
 /**********************************************************
  *
  * clover_inv applies the inverse of the clover term
@@ -305,10 +312,10 @@ void clover(const int ieo,
 
 su3 ** sw1, ** sw_inv1;
 su3 * _sw, *_sw_inv;
-static int sw_init = 0;
 
 void init_sw_fields(const int V) {
   su3 * tmp;
+  static int sw_init = 0;
 
   if(!sw_init) {
     if((void*)(sw = (su3***)calloc(V, sizeof(su3**))) == NULL) {
@@ -336,8 +343,8 @@ void init_sw_fields(const int V) {
       sw_inv[i] = sw_inv[i-1]+3;
     }
 #    if (defined SSE || defined SSE2 || defined SSE3)
-    sw[0][0] = (su3**)(((unsigned long int)(_sw)+ALIGN_BASE)&~ALIGN_BASE);
-    sw_inv[0][0] = (su3**)(((unsigned long int)(_sw_inv)+ALIGN_BASE)&~ALIGN_BASE);
+    sw[0][0] = (su3*)(((unsigned long int)(_sw)+ALIGN_BASE)&~ALIGN_BASE);
+    sw_inv[0][0] = (su3*)(((unsigned long int)(_sw_inv)+ALIGN_BASE)&~ALIGN_BASE);
 #    else
     sw[0][0] = _sw;
     sw_inv[0][0] = _sw_inv;

clover.h

@@ -31,6 +31,7 @@ extern su3 ** swm, ** swp;
 void Qsw_psi(spinor * const l, spinor * const k);
 void Qsw_sq_psi(spinor * const l, spinor * const k);
 void Msw_psi(spinor * const l, spinor * const k);
+void H_eo_sw_inv_psi(spinor * const l, spinor * const k, const int ieo);
 void init_sw_fields(const int V);
 
 #endif

clover_leaf.c

@@ -789,3 +789,40 @@ void sw_all(hamiltonian_field_t * const hf, const double kappa, const double c_s
   }
   return;
 }
+
+su3 * _swp;
+
+int init_swpm(const int V) {
+  int i=0;
+  static int swpm_init=0;
+
+  if(!swpm_init) {
+    if((void*)(swp = (su3**)calloc(V, sizeof(su3*))) == NULL) {
+      printf ("malloc errno : %d\n",errno); 
+      errno = 0;
+      return(1);
+    }
+    if((void*)(swm = (su3**)calloc(V, sizeof(su3*))) == NULL) {
+      printf ("malloc errno : %d\n",errno); 
+      errno = 0;
+      return(1);
+    }
+    if((void*)(_swp = (su3*)calloc(2*4*V+1, sizeof(su3))) == NULL) {
+      printf ("malloc errno : %d\n",errno); 
+      errno = 0;
+      return(2);
+    }
+#if (defined SSE || defined SSE2 || defined SSE3)
+    swp[0] = (su3*)(((unsigned long int)(_swp)+ALIGN_BASE)&~ALIGN_BASE);
+#else
+    swp[0] = _swp;
+#endif
+    swm[0] = swp[0] + V;
+    for(i = 1; i < V; i++){
+      swp[i] = swp[i-1]+4;
+      swm[i] = swm[i-1]+4;
+    }
+    swpm_init = 1;
+  }
+  return(0);
+}

clover_leaf.h

@@ -32,5 +32,6 @@ void sw_invert(const int ieo);
 void sw_deriv(const int ieo);
 void sw_spinor(const int ieo, spinor * const kk, spinor * const ll);
 void sw_all(hamiltonian_field_t * const hf, const double kappa, const double c_sw);
+int init_swpm(const int V);
 
 #endif

clover_trlog_monomial.c

@@ -29,32 +29,45 @@
 #include "su3.h"
 #include "su3adj.h"
 #include "su3spinor.h"
-#include "ranlxd.h"
-#include "start.h"
-#include "read_input.h"
 #include "clover_leaf.h"
-
+#include "clover.h"
 #include "monomial.h"
 #include "clover_trlog_monomial.h"
 
-void clover_trlog_derivative(const int no, hamiltonian_field_t * const hf) {
-
+void clover_trlog_derivative(const int id, hamiltonian_field_t * const hf) {
+  monomial * mnl = &monomial_list[id];
+  /* this term has no derivative */
+  /* so a dummy function         */
+  if(g_proc_id == 0 && g_debug_level > 4) {
+    printf("called clover_trlog_derivative for id %d %d, which is a dummy function\n", id, mnl->even_odd_flag);
+  }
   return;
 }
 
 
 void clover_trlog_heatbath(const int id, hamiltonian_field_t * const hf) {
   monomial * mnl = &monomial_list[id];
+  mnl->energy0 = 0.;
+
+  init_sw_fields(VOLUME);
   sw_term(hf->gaugefield, mnl->kappa, mnl->c_sw); 
-  /*compute the contribution from the clover-term*/
-  mnl->energy0 = 2.*sw_trace(1);
+  /*compute the contribution from the clover trlog term */
+  mnl->energy0 = 2.*sw_trace(0);
+  if(g_proc_id == 0 && g_debug_level > 3) {
+    printf("called clover_trlog_heatbath for id %d %d\n", id, mnl->even_odd_flag);
+  }
   return;
 }
 
 double clover_trlog_acc(const int id, hamiltonian_field_t * const hf) {
   monomial * mnl = &monomial_list[id];
+  mnl->energy1 = 0.;
   sw_term(hf->gaugefield, mnl->kappa, mnl->c_sw); 
-  /*compute the contribution from the clover-term*/
-  mnl->energy1 = 2.*sw_trace(1);  
+  /*compute the contribution from the clover trlog term */
+  mnl->energy1 = 2.*sw_trace(0);  
+  if(g_proc_id == 0 && g_debug_level > 3) {
+    printf("called clover_trlog_acc for id %d %d dH = %1.4e\n", 
+	   id, mnl->even_odd_flag, mnl->energy1 - mnl->energy0);
+  }
   return(mnl->energy1 - mnl->energy0);
 }

clover_trlog_monomial.h

@@ -23,6 +23,7 @@
 
 #include "hamiltonian_field.h"
 
+void clover_trlog_derivative(const int id, hamiltonian_field_t * const hf);
 void clover_trlog_heatbath(const int id, hamiltonian_field_t * const hf);
 double clover_trlog_acc(const int id, hamiltonian_field_t * const hf);
 

hmc_tm.c

@@ -201,8 +201,25 @@ int main(int argc,char *argv[]) {
 
   tmlqcd_mpi_init(argc, argv);
 
+  if(even_odd_flag) {
+    j = init_monomials(VOLUMEPLUSRAND/2, even_odd_flag);
+  }
+  else {
+    j = init_monomials(VOLUMEPLUSRAND, even_odd_flag);
+  }
+  if (j != 0) {
+    fprintf(stderr, "Not enough memory for monomial pseudo fermion fields! Aborting...\n");
+    exit(0);
+  }
+
   init_integrator();
 
+  if(g_proc_id == 0) {
+    for(j = 0; j < no_monomials; j++) {
+      printf("# monomial id %d type = %d timescale %d\n", j, monomial_list[j].type, monomial_list[j].timescale);
+    }
+  }
+
   if(nstore == -1) {
     countfile = fopen(nstore_filename, "r");
     if(countfile != NULL) {
@@ -224,11 +241,6 @@ int main(int argc,char *argv[]) {
   g_dbw2rand = 0;
 #endif
 
-  if(g_proc_id == 0) {
-    for(j = 0; j < no_monomials; j++) {
-      printf("# monomial id %d type = %d timescale %d\n", j, monomial_list[j].type, monomial_list[j].timescale);
-    }
-  }
 
   g_mu = g_mu1;
 
@@ -246,16 +258,6 @@ int main(int argc,char *argv[]) {
     fprintf(stderr, "Not enough memory for geometry_indices! Aborting...\n");
     exit(0);
   }
-  if(even_odd_flag) {
-    j = init_monomials(VOLUMEPLUSRAND/2, even_odd_flag);
-  }
-  else {
-    j = init_monomials(VOLUMEPLUSRAND, even_odd_flag);
-  }
-  if (j != 0) {
-    fprintf(stderr, "Not enough memory for monomial pseudo fermion fields! Aborting...\n");
-    exit(0);
-  }
   if(even_odd_flag) {
     j = init_spinor_field(VOLUMEPLUSRAND/2, NO_OF_SPINORFIELDS);
   }

integrator.c

@@ -79,7 +79,7 @@ int init_integrator() {
 
   for(i = 0; i < no_monomials; i++) {
     ts = monomial_list[i].timescale;
-    if(ts < Integrator.no_timescales && ts > -1 && monomial_list[i].type != CLOVERTRLOG) {
+    if(ts < Integrator.no_timescales && ts > -1) {
       Integrator.mnls_per_ts[ ts ][ Integrator.no_mnls_per_ts[ts] ] = monomial_list[i].id;
       Integrator.no_mnls_per_ts[ ts ]++;
     }

monomial.c

@@ -32,6 +32,9 @@
 #include "su3.h"
 #include "su3adj.h"
 #include "su3spinor.h"
+#include "tm_operators.h"
+#include "clover_leaf.h"
+#include "clover.h"
 #include "ranlxd.h"
 #include "sse.h"
 #include "linalg_eo.h"
@@ -143,9 +146,10 @@ int init_monomials(const int V, const int even_odd_flag) {
 
   no = 0;
   for(i = 0; i < no_monomials; i++) {
-    if(monomial_list[i].c_sw > 0) {
-      clover_trlog_monomial = 1;
-    }
+    /* Should be here if implementation is finished */
+    /*     if(monomial_list[i].c_sw > 0) { */
+    /*       clover_trlog_monomial = 1; */
+    /*     } */
     if((monomial_list[i].type != GAUGE) && (monomial_list[i].type != SFGAUGE)) {
 
       monomial_list[i].pf = __pf+no*V;
@@ -157,6 +161,18 @@ int init_monomials(const int V, const int even_odd_flag) {
 	monomial_list[i].accfunction = &det_acc;
 	monomial_list[i].derivativefunction = &det_derivative;
       }
+      else if(monomial_list[i].type == CLOVERDET) {
+	monomial_list[i].hbfunction = &cloverdet_heatbath;
+	monomial_list[i].accfunction = &cloverdet_acc;
+	monomial_list[i].derivativefunction = &cloverdet_derivative;
+	monomial_list[i].mu = 0.;
+	monomial_list[i].even_odd_flag = 0;
+	monomial_list[i].Qsq = &Qsw_sq_psi;
+	monomial_list[i].Qp = &Qsw_psi;
+	monomial_list[i].Qm = &Qsw_psi;
+	init_swpm(VOLUME);
+	clover_trlog_monomial = 1;
+      }
       else if(monomial_list[i].type == DETRATIO) {
 	monomial_list[i].hbfunction = &detratio_heatbath;
 	monomial_list[i].accfunction = &detratio_acc;
@@ -223,6 +239,8 @@ int init_monomials(const int V, const int even_odd_flag) {
     monomial_list[i].even_odd_flag = even_odd_flag;
   }
   if(clover_trlog_monomial && even_odd_flag) {
+    monomial_list[no_monomials].type = CLOVERTRLOG;
+    strcpy( monomial_list[no_monomials].name, "CLTRLOG");
     add_monomial(CLOVERTRLOG);
     monomial_list[no_monomials-1].pf = NULL;
     monomial_list[no_monomials-1].id = no_monomials-1;
@@ -231,8 +249,12 @@ int init_monomials(const int V, const int even_odd_flag) {
     monomial_list[no_monomials-1].mu = g_kappa;
     monomial_list[no_monomials-1].hbfunction = &clover_trlog_heatbath;
     monomial_list[no_monomials-1].accfunction = &clover_trlog_acc;
-    monomial_list[no_monomials-1].derivativefunction = &dummy_derivative;
+    monomial_list[no_monomials-1].derivativefunction = &clover_trlog_derivative;
     monomial_list[no_monomials-1].timescale = 0;
+    monomial_list[no_monomials-1].even_odd_flag = even_odd_flag;
+    if(g_proc_id == 0) {
+      printf("# Initialised clover_trlog_monomial, no_monomials= %d\n", no_monomials);
+    }
   }
   return(0);
 }

monomial.h

@@ -36,6 +36,7 @@
 #define NDDETRATIO 6
 #define POLYDETRATIO 7
 #define CLOVERTRLOG 8
+#define CLOVERDET 9
 
 #define max_no_monomials 20
 
@@ -97,6 +98,10 @@ typedef struct {
   void (*hbfunction) (const int no, hamiltonian_field_t * const hf);
   double (*accfunction) (const int no, hamiltonian_field_t * const hf);
   void (*derivativefunction) (const int no, hamiltonian_field_t * const hf);
+  /* the operator definitions */
+  void (*Qsq) (spinor * const, spinor * const);
+  void (*Qp) (spinor * const, spinor * const);
+  void (*Qm) (spinor * const, spinor * const);
 } monomial;
 
 #include "det_monomial.h"
@@ -107,6 +112,7 @@ typedef struct {
 #include "gauge_monomial.h"
 #include "sf_gauge_monomial.h"
 #include "clover_trlog_monomial.h"
+#include "cloverdet_monomial.h"
 
 /* list of all monomials */
 extern monomial monomial_list[max_no_monomials];