CHROMA
meslate.cc
Go to the documentation of this file.
1 
2 #error "NOT FULLY CONVERTED"
3 
4 
5 /* MESLATE - measure the lattice energy and the naive topological charge, */
6 /* with the "clover" definition */
7 /* of the naive topological charge density. It also computes the */
8 /* ratio of the action to the continuum instaton action. */
9 
10 /* u -- gauge field (Read) */
11 /* i_cool -- cooling sweep number (used for printing out headers) (Read) */
12 /* t_slice -- time slice number (Read) */
13 /* qtop -- topological charge (Write) */
14 /* S_ratio -- action to continuum instanton action (Write) */
15 include(types.mh)
16 
17 SUBROUTINE(meslate, u, i_cool, j_decay, t_slice, qtop, S_ratio)
18 
19 
20 multi1d<LatticeColorMatrix> u(Nd);
21 int i_cool;
22 int j_decay;
23 int t_slice;
24 Double qtop;
25 Double S_ratio;
26 
27 { /* Local Variables */
28  include(COMMON_DECLARATIONS)
29  /* local stuff for printing lattice action */
30  int ix;
31  int iy;
32  int iz;
33  int it;
34  int x;
35  int y;
36  int z;
37  int t;
38  int lx;
39  int ly;
40  int lz;
41  int lt;
42  int i;
43  multi1d<int> coord(Nd);
44  LatticeReal lract;
45  LatticeReal lrqtop;
46 
47  /* other local variables */
48  LatticeColorMatrix u_clov1;
49  LatticeColorMatrix u_clov2;
50  LatticeColorMatrix tmp_1;
51  LatticeColorMatrix tmp_2;
52  LatticeColorMatrix tmp_3;
53  LatticeReal qtop_tmp;
54  LatticeReal plaq_tmp;
55 
56  multi1d<Double> action_slice(lt);
57  multi1d<Double> qtop_slice(lt);
58  Double plaq;
59  Real rtmp;
60  Double tmp;
61  int c;
62  int cb;
63  int cb_0;
64  int mu;
65  int nu;
66  int mu1;
67  int nu1;
68  int mu2;
69  int nu2;
70 
71  START_CODE();
72 
73  if ( Nd != 4 )
74  QDP_error_exit("code requires 4 dimensions for printing[only], Nd);
75 
76  if ( j_decay != Nd-1 && j_decay > 0 )
77  QDP_error_exit("code requires j_decay == Nd-1 and > 0", j_decay);
78 
79  lx = nrow[0];
80  ly = nrow[1];
81  lz = nrow[2];
82  lt = nrow[3];
83 
84 
85  lract = 0;
86  lrqtop = 0;
87  qtop = 0;
88  plaq = 0;
89 
90  /* Lattice version of S_ratio */
91  rtmp = TO_REAL(2*Nd*(Nd-1)*Nc);
92  FILL(lract(0), rtmp);
93  FILL(lract(1), rtmp);
94 
95  /* Loop over checkerboards and triplet of perpendicular planes */
96  mu1 = 0;
97  for(cb = 0;cb <= ( 1); ++cb )
98  for(nu1 = 1;nu1 < ( Nd); ++nu1 )
99  {
100  /* First "plus-plus" plaquette */
101  /* tmp_1(x) = u(x+mu1,nu1) */
102  tmp_1(rb[cb]) = shift(u[nu1][1-cb], FORWARD, mu1);
103 
104  /* tmp_2(x) = u(x+nu1,mu1) */
105  tmp_2(rb[cb]) = shift(u[mu1][1-cb], FORWARD, nu1);
106 
107  /* tmp_3(x) = tmp_1 * tmp_2^dag = u(x+mu1,nu1)*u_dag(x+nu1,mu1) */
108  tmp_3 = tmp_1 * adj(tmp_2);
109 
110  /* tmp_1(x) = tmp_3 * u_dag(x,nu1)= u(x+mu1,nu1)*u_dag(x+nu1,mu1)*u_dag(x,nu1) */
111  tmp_1 = tmp_3 * adj(u[nu1][cb]);
112 
113  /* u_clov1(x) = u(x,mu1) * tmp_1 = u(x,mu1)*u(x+mu1,nu1)* */
114  /* u_dag(x+nu1,mu1)*u_dag(x,nu1) */
115  u_clov1 = u[mu1][cb] * tmp_1;
116 
117  plaq_tmp = real(trace(u_clov1));
118  plaq += sum(plaq_tmp);
119  lract[cb] -= plaq_tmp;
120 
121  /* First "plus-minus" plaquette */
122  /* tmp_1(x) = u(x+mu1,nu1) */
123  tmp_1(rb[1-cb]) = shift(u[nu1][cb], FORWARD, mu1);
124 
125  /* tmp_2(x) = u(x,mu1) * tmp_1 = u(x,mu1)*u(x+mu1,nu1) */
126  tmp_2 = u[mu1][1-cb] * tmp_1;
127 
128  /* tmp_1(x) = tmp_2_dag * u(x,nu1) = u_dag(x+mu1,nu1)*u_dag(x,mu1)*u(x,nu1) */
129  tmp_1 = adj(tmp_2) * u[nu1][1-cb];
130 
131  /* tmp_2(x) = tmp_1(x-nu1) */
132  tmp_2(rb[cb]) = shift(tmp_1, BACKWARD, nu1);
133 
134  /* tmp_1(x) = u(x,mu1) * tmp_2 = u(x,mu1)*u_dag(x-nu1+mu1,nu1)* */
135  /* u_dag(x-nu1,mu1)*u(x-nu1,nu1) */
136  tmp_1 = u[mu1][cb] * tmp_2;
137 
138  u_clov1 -= tmp_1;
139 
140  plaq_tmp = real(trace(tmp_1));
141  plaq += sum(plaq_tmp);
142  lract[cb] -= plaq_tmp;
143 
144  /* First "minus-minus" plaquette */
145  /* tmp_1(x) = u(x+mu1,nu1) */
146  tmp_1(rb[cb]) = shift(u[nu1][1-cb], FORWARD, mu1);
147 
148  /* tmp_2(x) = u(x,mu1) * tmp_1 = u(x,mu1)*u(x+mu1,nu1) */
149  tmp_2 = u[mu1][cb] * tmp_1;
150 
151  /* tmp_1(x) = u_dag(x,nu1) * tmp_2 = u_dag(x,nu1)*u(x,mu1)*u(x+mu1,nu1) */
152  tmp_1 = adj(u[nu1][cb]) * tmp_2;
153 
154  /* tmp_2(x) = tmp_1(x-nu1) */
155  tmp_2(rb[1-cb]) = shift(tmp_1, BACKWARD, nu1);
156 
157  /* tmp_1(x) = u_dag(x,mu1) * tmp_2 = u_dag(x,mu1)*u_dag(x-nu1,nu1)* */
158  /* u(x-nu1,mu1)*u(x-nu1+mu1,nu1) */
159  tmp_1 = adj(u[mu1][1-cb]) * tmp_2;
160 
161  /* tmp_2(x) = tmp_1(x-mu1) */
162  tmp_2(rb[cb]) = shift(tmp_1, BACKWARD, mu1);
163 
164  u_clov1 += tmp_2;
165 
166  plaq_tmp = real(trace(tmp_2));
167  plaq += sum(plaq_tmp);
168  lract[cb] -= plaq_tmp;
169 
170  /* First "minus-plus" plaquette */
171  /* tmp_1(x) = u_dag(x,mu1) * u(x,nu1) */
172  tmp_1 = adj(u[mu1][1-cb]) * u[nu1][1-cb];
173 
174  /* tmp_2(x) = u(x+nu1,mu1) */
175  tmp_2(rb[1-cb]) = shift(u[mu1][cb], FORWARD, nu1);
176 
177  /* tmp_3(x) = tmp_1 * tmp_2 = u_dag(x,mu1)*u(x,nu1)*u(x+nu1,mu1) */
178  tmp_3 = tmp_1 * tmp_2;
179 
180  /* tmp_1(x) = tmp_3(x-mu1) */
181  tmp_1(rb[cb]) = shift(tmp_3, BACKWARD, mu1);
182 
183  /* tmp_2(x) = tmp_1 * u_dag(x,nu1) = u_dag(x-mu1,mu1)*u(x-mu1,nu1)* */
184  /* u(x-mu1+nu1,mu1)*u_dag(x,nu1) */
185  tmp_2 = tmp_1 * adj(u[nu1][cb]);
186 
187  u_clov1 -= tmp_2;
188 
189  plaq_tmp = real(trace(tmp_2));
190  plaq += sum(plaq_tmp);
191  lract[cb] -= plaq_tmp;
192 
193  mu2 = INTEGER_MOD_FUNCTION(nu1,3) + 1;
194  nu2 = INTEGER_MOD_FUNCTION(mu2,3) + 1;
195 
196 /* Second "plus-plus" plaquette */
197  /* tmp_1(x) = u(x+mu2,nu2) */
198  tmp_1(rb[cb]) = shift(u[nu2][1-cb], FORWARD, mu2);
199 
200  /* tmp_2(x) = u(x+nu2,mu2) */
201  tmp_2(rb[cb]) = shift(u[mu2][1-cb], FORWARD, nu2);
202 
203  /* tmp_3(x) = tmp_1 * tmp_2^dag = u(x+mu2,nu2)*u_dag(x+nu2,mu2) */
204  tmp_3 = tmp_1 * adj(tmp_2);
205 
206  /* tmp_1(x) = tmp_3 * u_dag(x,nu2)= u(x+mu2,nu2)*u_dag(x+nu2,mu2)*u_dag(x,nu2) */
207  tmp_1 = tmp_3 * adj(u[nu2][cb]);
208 
209  /* u_clov2(x) = u(x,mu2) * tmp_1 = u(x,mu2)*u(x+mu2,nu2)* */
210  /* u_dag(x+nu2,mu2)*u_dag(x,nu2) */
211  u_clov2 = u[mu2][cb] * tmp_1;
212 
213  plaq_tmp = real(trace(u_clov2));
214  plaq += sum(plaq_tmp);
215  lract[cb] -= plaq_tmp;
216 
217  /* Second "plus-minus" plaquette */
218  /* tmp_1(x) = u(x+mu2,nu2) */
219  tmp_1(rb[1-cb]) = shift(u[nu2][cb], FORWARD, mu2);
220 
221  /* tmp_2(x) = u(x,mu2) * tmp_1 = u(x,mu2)*u(x+mu2,nu2) */
222  tmp_2 = u[mu2][1-cb] * tmp_1;
223 
224  /* tmp_1(x) = tmp_2_dag * u(x,nu2) = u_dag(x+mu2,nu2)*u_dag(x,mu2)*u(x,nu2) */
225  tmp_1 = adj(tmp_2) * u[nu2][1-cb];
226 
227  /* tmp_2(x) = tmp_1(x-nu2) */
228  tmp_2(rb[cb]) = shift(tmp_1, BACKWARD, nu2);
229 
230  /* tmp_1(x) = u(x,mu2) * tmp_2 = u(x,mu2)*u_dag(x-nu2+mu2,nu2)* */
231  /* u_dag(x-nu2,mu2)*u(x-nu2,nu2) */
232  tmp_1 = u[mu2][cb] * tmp_2;
233 
234  u_clov2 -= tmp_1;
235 
236  plaq_tmp = real(trace(tmp_1));
237  plaq += sum(plaq_tmp);
238  lract[cb] -= plaq_tmp;
239 
240  /* Second "minus-minus" plaquette */
241  /* tmp_1(x) = u(x+mu2,nu2) */
242  tmp_1(rb[cb]) = shift(u[nu2][1-cb], FORWARD, mu2);
243 
244  /* tmp_2(x) = u(x,mu2) * tmp_1 = u(x,mu2)*u(x+mu2,nu2) */
245  tmp_2 = u[mu2][cb] * tmp_1;
246 
247  /* tmp_1(x) = u_dag(x,nu2) * tmp_2 = u_dag(x,nu2)*u(x,mu2)*u(x+mu2,nu2) */
248  tmp_1 = adj(u[nu2][cb]) * tmp_2;
249 
250  /* tmp_2(x) = tmp_1(x-nu2) */
251  tmp_2(rb[1-cb]) = shift(tmp_1, BACKWARD, nu2);
252 
253  /* tmp_1(x) = u_dag(x,mu2) * tmp_2 = u_dag(x,mu2)*u_dag(x-nu2,nu2)* */
254  /* u(x-nu2,mu2)*u(x-nu2+mu2,nu2) */
255  tmp_1 = adj(u[mu2][1-cb]) * tmp_2;
256 
257  /* tmp_2(x) = tmp_1(x-mu2) */
258  tmp_2(rb[cb]) = shift(tmp_1, BACKWARD, mu2);
259 
260  u_clov2 += tmp_2;
261 
262  plaq_tmp = real(trace(tmp_2));
263  plaq += sum(plaq_tmp);
264  lract[cb] -= plaq_tmp;
265 
266  /* Second "minus-plus" plaquette */
267  /* tmp_1(x) = u_dag(x,mu2) * u(x,nu2) */
268  tmp_1 = adj(u[mu2][1-cb]) * u[nu2][1-cb];
269 
270  /* tmp_2(x) = u(x+nu2,mu2) */
271  tmp_2(rb[1-cb]) = shift(u[mu2][cb], FORWARD, nu2);
272 
273  /* tmp_3(x) = tmp_1 * tmp_2 = u_dag(x,mu2)*u(x,nu2)*u(x+nu2,mu2) */
274  tmp_3 = tmp_1 * tmp_2;
275 
276  /* tmp_1(x) = tmp_3(x-mu2) */
277  tmp_1(rb[cb]) = shift(tmp_3, BACKWARD, mu2);
278 
279  /* tmp_2(x) = tmp_1 * u_dag(x,nu2) = u_dag(x-mu2,mu2)*u(x-mu2,nu2)* */
280  /* u(x-mu2+nu2,mu2)*u_dag(x,nu2) */
281  tmp_2 = tmp_1 * adj(u[nu2][cb]);
282 
283  u_clov2 -= tmp_2;
284 
285  plaq_tmp = real(trace(tmp_2));
286  plaq += sum(plaq_tmp);
287  lract[cb] -= plaq_tmp;
288 
289  /* Now comes the contribution to the topological charge */
290  tmp_1 = 1;
291  tmp_2 = adj(u_clov1) * tmp_1;
292  u_clov1 -= tmp_2;
293  tmp_2 = adj(u_clov2) * tmp_1;
294  u_clov2 -= tmp_2;
295  tmp_2 = u_clov1 * u_clov2;
296 
297  qtop_tmp = real(trace(tmp_2));
298  lrqtop[cb] -= qtop_tmp;
299  }
300 
301 
302  plaq = plaq / TO_DOUBLE(2*(vol)*Nd*(Nd-1)*Nc);
303 
304  /* Lattice version of S_ratio */
305  FILL(rtmp,PI);
306  rtmp = TO_REAL(WORD_VALUE(WORD_rtmp,ONE) / (16*rtmp*rtmp));
307  lract[0] = lract[0] * rtmp;
308  lract[1] = lract[1] * rtmp;
309 
310  /* Check */
311  S_ratio = sum(lract[0]);
312  S_ratio += sum(lract[1]);
313 
314  /* Lattice version of qtop */
315  FILL(rtmp,PI);
316  rtmp = TO_REAL(WORD_VALUE(WORD_rtmp,ONE) / (256*rtmp*rtmp));
317  lrqtop[0] = lrqtop[0] * rtmp;
318  lrqtop[1] = lrqtop[1] * rtmp;
319 
320  /* Check */
321  qtop = sum(lrqtop[0]);
322  qtop += sum(lrqtop[1]);
323 
324  /* sum the action density over each time slice */
325  action_slice = sumMulti(lract[0], timeslice);
326  action_slice += sumMulti(lract[1], timeslice);
327 
328  /* sum the topological charge density over each time slice */
329  qtop_slice = sumMulti(lrqtop[0], timeslice);
330  qtop_slice += sumMulti(lrqtop[1], timeslice);
331 
332 
333  /*### Print out the action and topological charge density in SCIAN stf format */
334 
335  fprintf(unit_21," TIME %4d\n RANK 3\n DIMENSIONS %4d%4d%4d\n\
336  BOUNDS 0 %4d 0 %4d 0 %4d\n SCALAR\n ORDER COLUMN\n INTERLACED\n DATA\n",
337  i_cool,lx+1,ly+1,lz+1,lx,ly,lz);
338  fprintf(unit_22," TIME %4d\n RANK 3\n DIMENSIONS %4d%4d%4d\n\
339  BOUNDS 0 %4d 0 %4d 0 %4d\n SCALAR\n ORDER COLUMN\n INTERLACED\n DATA\n",
340  i_cool,lx+1,ly+1,lz+1,lx,ly,lz);
341 
342  /* Extend the x,y,z axes for padding. Done for Scian and periodic BC. */
343  for(z = 0;z <= ( lz); ++z )
344  for(y = 0;y <= ( ly); ++y )
345  for(x = 0;x <= ( lx); ++x )
346  {
347  ix = INTEGER_MOD_FUNCTION(x,lx);
348  iy = INTEGER_MOD_FUNCTION(y,ly);
349  iz = INTEGER_MOD_FUNCTION(z,lz);
350  it = INTEGER_MOD_FUNCTION(t_slice,lt);
351 
352  cb = INTEGER_MOD_FUNCTION(ix+iy+iz+it,2);
353  i = (int)(ix/2);
354 
355  coord[0] = i;
356  coord[1] = iy;
357  coord[2] = iz;
358  coord[3] = it;
359 
360  INDEXING(lract(cb),coord,rtmp);
361  fprintf(unit_21," %14.7g\n",rtmp);
362 
363  INDEXING(lrqtop(cb),coord,rtmp);
364  fprintf(unit_22," %14.7g\n",rtmp);
365  }
366 
367  fprintf(unit_21," END\n");
368  fprintf(unit_22," END\n");
369 
370  /*### END ############### */
371 
372 
373  /*### Print out the action and topological charge density in SCIAN stf format */
374 
375  lx = cb_nrow[0];
376 
377 
378  fprintf(unit_23," TIME %4d\n RANK 4\n DIMENSIONS %4d%4d%4d%4d\n\
379  BOUNDS 0 %4d 0 %4d 0 %4d 0 %4d\n SCALAR\n ORDER COLUMN\n INTERLACED\n DATA\n",
380  i_cool,lx+1,ly+1,lz+1,lt+1,lx,ly,lz,lt);
381 
382 
383  /* Extend the x,y,z axes for padding. Done for Scian and periodic BC. */
384  for(t = 0;t <= ( lt); ++t )
385  for(z = 0;z <= ( lz); ++z )
386  for(y = 0;y <= ( ly); ++y )
387  for(x = 0;x <= ( lx); ++x )
388  {
389  ix = INTEGER_MOD_FUNCTION(x,lx);
390  iy = INTEGER_MOD_FUNCTION(y,ly);
391  iz = INTEGER_MOD_FUNCTION(z,lz);
392  it = INTEGER_MOD_FUNCTION(t,lt);
393 
394  coord[0] = ix;
395  coord[1] = iy;
396  coord[2] = iz;
397  coord[3] = it;
398 
399  INDEXING(lract(0),coord,rtmp);
400  fprintf(unit_23," %14.7g\n",rtmp);
401  }
402 
403  fprintf(unit_23," END\n");
404 
405  /*### END ############### */
406 
407  printf(" %d %d %g %g %g\n",i_cool,t_slice,plaq,qtop,S_ratio);
408  push(xml_out,"MesLatE");
409 write(xml_out, "i_cool", i_cool);
410 write(xml_out, "t_slice", t_slice);
411 write(xml_out, "plaq", plaq);
412 write(xml_out, "qtop", qtop);
413 write(xml_out, "S_ratio", S_ratio);
414 write(xml_out, "qtop_slice", qtop_slice);
415 write(xml_out, "action_slice", action_slice);
416 pop(xml_out);
417 
418 
419  END_CODE();
420 }
n
unsigned n
Definition: ldumul_w.cc:36
START_CODE
START_CODE()
cb
int cb
Definition: meslate.cc:62
tmp
Double tmp
Definition: meslate.cc:60
y
int y
Definition: meslate.cc:35
nu1
int nu1
Definition: meslate.cc:67
z
int z
Definition: meslate.cc:36
it
int it
Definition: meslate.cc:33
mu
int mu
Definition: meslate.cc:64
action_slice
multi1d< Double > action_slice(lt)
tmp_2
LatticeColorMatrix tmp_2
Definition: meslate.cc:51
coord
multi1d< int > coord(Nd)
j_decay
int j_decay
Definition: meslate.cc:22
ly
int ly
Definition: meslate.cc:39
c
int c
Definition: meslate.cc:61
i_cool
int i_cool
Definition: meslate.cc:21
iy
int iy
Definition: meslate.cc:27
lz
int lz
Definition: meslate.cc:40
lt
int lt
Definition: meslate.cc:41
x
int x
Definition: meslate.cc:34
tmp_3
LatticeColorMatrix tmp_3
Definition: meslate.cc:52
iz
int iz
Definition: meslate.cc:32
S_ratio
Double S_ratio
Definition: meslate.cc:25
cb_0
int cb_0
Definition: meslate.cc:63
qtop_tmp
LatticeReal qtop_tmp
Definition: meslate.cc:53
lract
LatticeReal lract
Definition: meslate.cc:44
plaq
Double plaq
Definition: meslate.cc:58
plaq_tmp
LatticeReal plaq_tmp
Definition: meslate.cc:54
mu1
int mu1
Definition: meslate.cc:66
lx
int lx
Definition: meslate.cc:38
u_clov2
LatticeColorMatrix u_clov2
Definition: meslate.cc:49
include
include(types.mh) multi1d< LatticeColorMatrix > u(Nd)
t
int t
Definition: meslate.cc:37
qtop
Double qtop
Definition: meslate.cc:24
u_clov1
LatticeColorMatrix u_clov1
Definition: meslate.cc:48
i
int i
Definition: meslate.cc:42
qtop_slice
multi1d< Double > qtop_slice(lt)
mu2
int mu2
Definition: meslate.cc:68
tmp_1
LatticeColorMatrix tmp_1
Definition: meslate.cc:50
nu
int nu
Definition: meslate.cc:65
rtmp
Real rtmp
Definition: meslate.cc:59
t_slice
int t_slice
Definition: meslate.cc:23
lrqtop
LatticeReal lrqtop
Definition: meslate.cc:45
Nd
Nd
Definition: meslate.cc:74
nu2
int nu2
Definition: meslate.cc:69
Chroma::QDP_error_exit
QDP_error_exit("too many BiCG iterations", n_count, rsd_sq, cp, c, re_rvr, im_rvr, re_a, im_a, re_b, im_b)
Chroma::u
static multi1d< LatticeColorMatrix > u
Definition: syssolver_linop_qop_mg_w.cc:39
Chroma::d
DComplex d
Definition: invbicg.cc:99
testing::internal::Double
FloatingPoint< double > Double
Definition: gtest.h:7351
Generated by doxygen 1.9.1