CHROMA: meslate.cc Source File

00001 // $Id: meslate.cc,v 3.0 2006/04/03 04:58:58 edwards Exp $
00002 
00003 #error "NOT FULLY CONVERTED"
00004 
00005 
00006 /* MESLATE - measure the lattice energy and the naive topological charge,  */
00007 /*           with the "clover" definition */
00008 /*           of the naive topological charge density. It also computes the */
00009 /*           ratio of the action to the continuum instaton action. */
00010 
00011 /* u -- gauge field (Read) */
00012 /* i_cool  -- cooling sweep number (used for printing out headers) (Read) */
00013 /* t_slice -- time slice number (Read) */
00014 /* qtop    -- topological charge (Write) */
00015 /* S_ratio -- action to continuum instanton action (Write) */
00016 include(types.mh)
00017 
00018 SUBROUTINE(meslate, u, i_cool, j_decay, t_slice, qtop, S_ratio)
00019 
00020 
00021 multi1d<LatticeColorMatrix> u(Nd);
00022 int i_cool;
00023 int j_decay;
00024 int t_slice;
00025 Double qtop;
00026 Double S_ratio;
00027 
00028 { /* Local Variables */
00029   include(COMMON_DECLARATIONS)
00030   /* local stuff for printing lattice action */
00031   int ix;
00032   int iy;
00033   int iz;
00034   int it;
00035   int x;
00036   int y;
00037   int z;
00038   int t;
00039   int lx;
00040   int ly;
00041   int lz;
00042   int lt;
00043   int i;
00044   multi1d<int> coord(Nd);
00045   LatticeReal lract;
00046   LatticeReal lrqtop;
00047 
00048   /* other local variables */
00049   LatticeColorMatrix u_clov1;
00050   LatticeColorMatrix u_clov2;
00051   LatticeColorMatrix tmp_1;
00052   LatticeColorMatrix tmp_2;
00053   LatticeColorMatrix tmp_3;
00054   LatticeReal qtop_tmp;
00055   LatticeReal plaq_tmp;
00056 
00057   multi1d<Double> action_slice(lt);
00058   multi1d<Double> qtop_slice(lt);
00059   Double plaq;
00060   Real rtmp;
00061   Double tmp;
00062   int c;
00063   int cb;
00064   int cb_0;
00065   int mu;
00066   int nu;
00067   int mu1;
00068   int nu1;
00069   int mu2;
00070   int nu2;
00071 
00072   START_CODE();
00073 
00074   if ( Nd != 4 )
00075     QDP_error_exit("code requires 4 dimensions for printing[only], Nd);
00076 
00077   if ( j_decay != Nd-1 && j_decay > 0 )
00078     QDP_error_exit("code requires j_decay == Nd-1 and > 0", j_decay);
00079 
00080   lx = nrow[0];
00081   ly = nrow[1];
00082   lz = nrow[2];
00083   lt = nrow[3];
00084 
00085                   
00086   lract = 0;
00087   lrqtop = 0;
00088   qtop = 0;
00089   plaq = 0;
00090 
00091   /* Lattice version of S_ratio */
00092   rtmp = TO_REAL(2*Nd*(Nd-1)*Nc);
00093   FILL(lract(0), rtmp);
00094   FILL(lract(1), rtmp);
00095 
00096   /* Loop over checkerboards and triplet of perpendicular planes */
00097   mu1 = 0;
00098   for(cb = 0;cb  <= ( 1); ++cb )
00099     for(nu1 = 1;nu1  < ( Nd); ++nu1 )
00100     {
00101       /* First "plus-plus" plaquette */
00102       /* tmp_1(x) = u(x+mu1,nu1) */
00103       tmp_1(rb[cb]) = shift(u[nu1][1-cb], FORWARD, mu1);
00104 
00105       /* tmp_2(x) = u(x+nu1,mu1) */
00106       tmp_2(rb[cb]) = shift(u[mu1][1-cb], FORWARD, nu1);
00107 
00108       /* tmp_3(x) = tmp_1 * tmp_2^dag = u(x+mu1,nu1)*u_dag(x+nu1,mu1) */
00109       tmp_3 = tmp_1 * adj(tmp_2);
00110 
00111       /* tmp_1(x) = tmp_3 * u_dag(x,nu1)= u(x+mu1,nu1)*u_dag(x+nu1,mu1)*u_dag(x,nu1) */
00112       tmp_1 = tmp_3 * adj(u[nu1][cb]);
00113 
00114       /* u_clov1(x) = u(x,mu1) * tmp_1 = u(x,mu1)*u(x+mu1,nu1)* */
00115       /*                                 u_dag(x+nu1,mu1)*u_dag(x,nu1) */
00116       u_clov1 = u[mu1][cb] * tmp_1;
00117 
00118             plaq_tmp = real(trace(u_clov1));
00119       plaq += sum(plaq_tmp);
00120       lract[cb] -= plaq_tmp;
00121       
00122       /* First "plus-minus" plaquette */
00123       /* tmp_1(x) = u(x+mu1,nu1) */
00124       tmp_1(rb[1-cb]) = shift(u[nu1][cb], FORWARD, mu1);
00125 
00126       /* tmp_2(x) = u(x,mu1) * tmp_1 = u(x,mu1)*u(x+mu1,nu1) */
00127       tmp_2 = u[mu1][1-cb] * tmp_1;
00128 
00129       /* tmp_1(x) = tmp_2_dag * u(x,nu1) = u_dag(x+mu1,nu1)*u_dag(x,mu1)*u(x,nu1) */
00130       tmp_1 = adj(tmp_2) * u[nu1][1-cb];
00131 
00132       /* tmp_2(x) = tmp_1(x-nu1) */
00133       tmp_2(rb[cb]) = shift(tmp_1, BACKWARD, nu1);
00134 
00135       /* tmp_1(x) = u(x,mu1) * tmp_2 = u(x,mu1)*u_dag(x-nu1+mu1,nu1)* */
00136       /*                               u_dag(x-nu1,mu1)*u(x-nu1,nu1) */
00137       tmp_1 = u[mu1][cb] * tmp_2;
00138 
00139       u_clov1 -= tmp_1;
00140 
00141             plaq_tmp = real(trace(tmp_1));
00142       plaq += sum(plaq_tmp);
00143       lract[cb] -= plaq_tmp;
00144       
00145       /* First "minus-minus" plaquette */
00146       /* tmp_1(x) = u(x+mu1,nu1) */
00147       tmp_1(rb[cb]) = shift(u[nu1][1-cb], FORWARD, mu1);
00148 
00149       /* tmp_2(x) = u(x,mu1) * tmp_1 = u(x,mu1)*u(x+mu1,nu1) */
00150       tmp_2 = u[mu1][cb] * tmp_1;
00151 
00152       /* tmp_1(x) = u_dag(x,nu1) * tmp_2 = u_dag(x,nu1)*u(x,mu1)*u(x+mu1,nu1) */
00153       tmp_1 = adj(u[nu1][cb]) * tmp_2;
00154 
00155       /* tmp_2(x) = tmp_1(x-nu1) */
00156       tmp_2(rb[1-cb]) = shift(tmp_1, BACKWARD, nu1);
00157 
00158       /* tmp_1(x) = u_dag(x,mu1) * tmp_2 = u_dag(x,mu1)*u_dag(x-nu1,nu1)* */
00159       /*                                   u(x-nu1,mu1)*u(x-nu1+mu1,nu1) */
00160       tmp_1 = adj(u[mu1][1-cb]) * tmp_2;
00161 
00162       /* tmp_2(x) = tmp_1(x-mu1) */
00163       tmp_2(rb[cb]) = shift(tmp_1, BACKWARD, mu1);
00164 
00165       u_clov1 += tmp_2;
00166 
00167             plaq_tmp = real(trace(tmp_2));
00168       plaq += sum(plaq_tmp);
00169       lract[cb] -= plaq_tmp;
00170       
00171       /* First "minus-plus" plaquette */
00172       /* tmp_1(x) = u_dag(x,mu1) * u(x,nu1) */
00173       tmp_1 = adj(u[mu1][1-cb]) * u[nu1][1-cb];
00174 
00175       /* tmp_2(x) = u(x+nu1,mu1) */
00176       tmp_2(rb[1-cb]) = shift(u[mu1][cb], FORWARD, nu1);
00177 
00178       /* tmp_3(x) = tmp_1 * tmp_2 = u_dag(x,mu1)*u(x,nu1)*u(x+nu1,mu1) */
00179       tmp_3 = tmp_1 * tmp_2;
00180 
00181       /* tmp_1(x) = tmp_3(x-mu1) */
00182       tmp_1(rb[cb]) = shift(tmp_3, BACKWARD, mu1);
00183 
00184       /* tmp_2(x) = tmp_1 * u_dag(x,nu1) = u_dag(x-mu1,mu1)*u(x-mu1,nu1)* */
00185       /*                                   u(x-mu1+nu1,mu1)*u_dag(x,nu1) */
00186       tmp_2 = tmp_1 * adj(u[nu1][cb]);
00187 
00188       u_clov1 -= tmp_2;
00189 
00190             plaq_tmp = real(trace(tmp_2));
00191       plaq += sum(plaq_tmp);
00192       lract[cb] -= plaq_tmp;
00193       
00194       mu2 = INTEGER_MOD_FUNCTION(nu1,3) + 1;
00195       nu2 = INTEGER_MOD_FUNCTION(mu2,3) + 1;
00196 
00197 /* Second "plus-plus" plaquette */
00198   /* tmp_1(x) = u(x+mu2,nu2) */
00199       tmp_1(rb[cb]) = shift(u[nu2][1-cb], FORWARD, mu2);
00200 
00201       /* tmp_2(x) = u(x+nu2,mu2) */
00202       tmp_2(rb[cb]) = shift(u[mu2][1-cb], FORWARD, nu2);
00203 
00204       /* tmp_3(x) = tmp_1 * tmp_2^dag = u(x+mu2,nu2)*u_dag(x+nu2,mu2) */
00205       tmp_3 = tmp_1 * adj(tmp_2);
00206 
00207       /* tmp_1(x) = tmp_3 * u_dag(x,nu2)= u(x+mu2,nu2)*u_dag(x+nu2,mu2)*u_dag(x,nu2) */
00208       tmp_1 = tmp_3 * adj(u[nu2][cb]);
00209 
00210       /* u_clov2(x) = u(x,mu2) * tmp_1 = u(x,mu2)*u(x+mu2,nu2)* */
00211       /*                                 u_dag(x+nu2,mu2)*u_dag(x,nu2) */
00212       u_clov2 = u[mu2][cb] * tmp_1;
00213 
00214             plaq_tmp = real(trace(u_clov2));
00215       plaq += sum(plaq_tmp);
00216       lract[cb] -= plaq_tmp;
00217       
00218       /* Second "plus-minus" plaquette */
00219       /* tmp_1(x) = u(x+mu2,nu2) */
00220       tmp_1(rb[1-cb]) = shift(u[nu2][cb], FORWARD, mu2);
00221 
00222       /* tmp_2(x) = u(x,mu2) * tmp_1 = u(x,mu2)*u(x+mu2,nu2) */
00223       tmp_2 = u[mu2][1-cb] * tmp_1;
00224 
00225       /* tmp_1(x) = tmp_2_dag * u(x,nu2) = u_dag(x+mu2,nu2)*u_dag(x,mu2)*u(x,nu2) */
00226       tmp_1 = adj(tmp_2) * u[nu2][1-cb];
00227 
00228       /* tmp_2(x) = tmp_1(x-nu2) */
00229       tmp_2(rb[cb]) = shift(tmp_1, BACKWARD, nu2);
00230 
00231       /* tmp_1(x) = u(x,mu2) * tmp_2 = u(x,mu2)*u_dag(x-nu2+mu2,nu2)* */
00232       /*                               u_dag(x-nu2,mu2)*u(x-nu2,nu2) */
00233       tmp_1 = u[mu2][cb] * tmp_2;
00234 
00235       u_clov2 -= tmp_1;
00236 
00237             plaq_tmp = real(trace(tmp_1));
00238       plaq += sum(plaq_tmp);
00239       lract[cb] -= plaq_tmp;
00240       
00241       /* Second "minus-minus" plaquette */
00242       /* tmp_1(x) = u(x+mu2,nu2) */
00243       tmp_1(rb[cb]) = shift(u[nu2][1-cb], FORWARD, mu2);
00244 
00245       /* tmp_2(x) = u(x,mu2) * tmp_1 = u(x,mu2)*u(x+mu2,nu2) */
00246       tmp_2 = u[mu2][cb] * tmp_1;
00247 
00248       /* tmp_1(x) = u_dag(x,nu2) * tmp_2 = u_dag(x,nu2)*u(x,mu2)*u(x+mu2,nu2) */
00249       tmp_1 = adj(u[nu2][cb]) * tmp_2;
00250 
00251       /* tmp_2(x) = tmp_1(x-nu2) */
00252       tmp_2(rb[1-cb]) = shift(tmp_1, BACKWARD, nu2);
00253 
00254       /* tmp_1(x) = u_dag(x,mu2) * tmp_2 = u_dag(x,mu2)*u_dag(x-nu2,nu2)* */
00255       /*                                   u(x-nu2,mu2)*u(x-nu2+mu2,nu2) */
00256       tmp_1 = adj(u[mu2][1-cb]) * tmp_2;
00257 
00258       /* tmp_2(x) = tmp_1(x-mu2) */
00259       tmp_2(rb[cb]) = shift(tmp_1, BACKWARD, mu2);
00260 
00261       u_clov2 += tmp_2;
00262 
00263             plaq_tmp = real(trace(tmp_2));
00264       plaq += sum(plaq_tmp);
00265       lract[cb] -= plaq_tmp;
00266       
00267       /* Second "minus-plus" plaquette */
00268       /* tmp_1(x) = u_dag(x,mu2) * u(x,nu2) */
00269       tmp_1 = adj(u[mu2][1-cb]) * u[nu2][1-cb];
00270 
00271       /* tmp_2(x) = u(x+nu2,mu2) */
00272       tmp_2(rb[1-cb]) = shift(u[mu2][cb], FORWARD, nu2);
00273 
00274       /* tmp_3(x) = tmp_1 * tmp_2 = u_dag(x,mu2)*u(x,nu2)*u(x+nu2,mu2) */
00275       tmp_3 = tmp_1 * tmp_2;
00276 
00277       /* tmp_1(x) = tmp_3(x-mu2) */
00278       tmp_1(rb[cb]) = shift(tmp_3, BACKWARD, mu2);
00279 
00280       /* tmp_2(x) = tmp_1 * u_dag(x,nu2) = u_dag(x-mu2,mu2)*u(x-mu2,nu2)* */
00281       /*                                   u(x-mu2+nu2,mu2)*u_dag(x,nu2) */
00282       tmp_2 = tmp_1 * adj(u[nu2][cb]);
00283 
00284       u_clov2 -= tmp_2;
00285 
00286             plaq_tmp = real(trace(tmp_2));
00287       plaq += sum(plaq_tmp);
00288       lract[cb] -= plaq_tmp;
00289       
00290       /* Now comes the contribution to the topological charge */
00291       tmp_1 = 1;
00292       tmp_2 = adj(u_clov1) * tmp_1;
00293       u_clov1 -= tmp_2;
00294       tmp_2 = adj(u_clov2) * tmp_1;
00295       u_clov2 -= tmp_2;
00296       tmp_2 = u_clov1 * u_clov2;
00297 
00298             qtop_tmp = real(trace(tmp_2));
00299       lrqtop[cb] -= qtop_tmp;
00300           }
00301 
00302           
00303   plaq = plaq / TO_DOUBLE(2*(vol)*Nd*(Nd-1)*Nc);
00304 
00305   /* Lattice version of S_ratio */
00306   FILL(rtmp,PI);
00307   rtmp = TO_REAL(WORD_VALUE(WORD_rtmp,ONE) / (16*rtmp*rtmp));
00308   lract[0] = lract[0] * rtmp;
00309   lract[1] = lract[1] * rtmp;
00310 
00311   /* Check */
00312   S_ratio = sum(lract[0]);
00313   S_ratio += sum(lract[1]);
00314 
00315   /* Lattice version of qtop */
00316   FILL(rtmp,PI);
00317   rtmp = TO_REAL(WORD_VALUE(WORD_rtmp,ONE) / (256*rtmp*rtmp));
00318   lrqtop[0] = lrqtop[0] * rtmp;
00319   lrqtop[1] = lrqtop[1] * rtmp;
00320 
00321   /* Check */
00322   qtop = sum(lrqtop[0]);
00323   qtop += sum(lrqtop[1]);
00324 
00325   /* sum the action density over each time slice */
00326   action_slice = sumMulti(lract[0], timeslice);
00327   action_slice += sumMulti(lract[1], timeslice);
00328 
00329   /* sum the topological charge density over each time slice */
00330   qtop_slice = sumMulti(lrqtop[0], timeslice);
00331   qtop_slice += sumMulti(lrqtop[1], timeslice);
00332 
00333 
00334   /*### Print out the action and topological charge density in SCIAN stf format */
00335   
00336   fprintf(unit_21," TIME %4d\n RANK 3\n DIMENSIONS %4d%4d%4d\n\
00337  BOUNDS 0 %4d 0 %4d 0 %4d\n SCALAR\n ORDER COLUMN\n INTERLACED\n DATA\n",
00338           i_cool,lx+1,ly+1,lz+1,lx,ly,lz);
00339   fprintf(unit_22," TIME %4d\n RANK 3\n DIMENSIONS %4d%4d%4d\n\
00340  BOUNDS 0 %4d 0 %4d 0 %4d\n SCALAR\n ORDER COLUMN\n INTERLACED\n DATA\n",
00341           i_cool,lx+1,ly+1,lz+1,lx,ly,lz);
00342 
00343   /* Extend the x,y,z axes for padding. Done for Scian and periodic BC. */
00344   for(z = 0;z  <= ( lz); ++z )
00345     for(y = 0;y  <= ( ly); ++y )
00346       for(x = 0;x  <= ( lx); ++x )
00347       {
00348         ix = INTEGER_MOD_FUNCTION(x,lx);
00349         iy = INTEGER_MOD_FUNCTION(y,ly);
00350         iz = INTEGER_MOD_FUNCTION(z,lz);
00351         it = INTEGER_MOD_FUNCTION(t_slice,lt);
00352 
00353         cb = INTEGER_MOD_FUNCTION(ix+iy+iz+it,2);
00354         i  = (int)(ix/2);
00355 
00356         coord[0] = i;
00357         coord[1] = iy;
00358         coord[2] = iz;
00359         coord[3] = it;
00360 
00361         INDEXING(lract(cb),coord,rtmp);
00362         fprintf(unit_21," %14.7g\n",rtmp);
00363 
00364         INDEXING(lrqtop(cb),coord,rtmp);
00365         fprintf(unit_22," %14.7g\n",rtmp);
00366       }
00367 
00368   fprintf(unit_21," END\n");
00369   fprintf(unit_22," END\n");
00370 
00371     /*### END ############### */
00372 
00373 
00374   /*### Print out the action and topological charge density in SCIAN stf format */
00375   
00376   lx = cb_nrow[0];
00377 
00378 
00379   fprintf(unit_23," TIME %4d\n RANK 4\n DIMENSIONS %4d%4d%4d%4d\n\
00380  BOUNDS 0 %4d 0 %4d 0 %4d 0 %4d\n SCALAR\n ORDER COLUMN\n INTERLACED\n DATA\n",
00381           i_cool,lx+1,ly+1,lz+1,lt+1,lx,ly,lz,lt);
00382 
00383 
00384   /* Extend the x,y,z axes for padding. Done for Scian and periodic BC. */
00385   for(t = 0;t  <= ( lt); ++t )
00386     for(z = 0;z  <= ( lz); ++z )
00387       for(y = 0;y  <= ( ly); ++y )
00388         for(x = 0;x  <= ( lx); ++x )
00389         {
00390           ix = INTEGER_MOD_FUNCTION(x,lx);
00391           iy = INTEGER_MOD_FUNCTION(y,ly);
00392           iz = INTEGER_MOD_FUNCTION(z,lz);
00393           it = INTEGER_MOD_FUNCTION(t,lt);
00394 
00395           coord[0] = ix;
00396           coord[1] = iy;
00397           coord[2] = iz;
00398           coord[3] = it;
00399 
00400           INDEXING(lract(0),coord,rtmp);
00401           fprintf(unit_23," %14.7g\n",rtmp);
00402         }
00403 
00404   fprintf(unit_23," END\n");
00405 
00406     /*### END ############### */
00407 
00408   printf(" %d %d %g %g %g\n",i_cool,t_slice,plaq,qtop,S_ratio);
00409   push(xml_out,"MesLatE");
00410 write(xml_out, "i_cool", i_cool);
00411 write(xml_out, "t_slice", t_slice);
00412 write(xml_out, "plaq", plaq);
00413 write(xml_out, "qtop", qtop);
00414 write(xml_out, "S_ratio", S_ratio);
00415 write(xml_out, "qtop_slice", qtop_slice);
00416 write(xml_out, "action_slice", action_slice);
00417 pop(xml_out);
00418 
00419         
00420   END_CODE();
00421 }