00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 static const int KS_NDIM=4;
00014
00015 #ifndef TWISTED_BOUNDARY
00016
00017 class phase_field : public mdp_field<int> {
00018 public:
00019 phase_field(mdp_lattice &a) {
00020 if(a.ndim!=4) error("fermiqcd_staggered_mesons/phase_field: ndim!=4");
00021 allocate_field(a,16);
00022 };
00023 int component(site x, site y) {
00024 int i=0, mu;
00025 for(mu=0; mu<KS_NDIM; mu++) {
00026 if(x(mu)/2 != y(mu)/2)
00027 return 0;
00028 i = (i << 1) + (x(mu) % 2);
00029 }
00030 return *address(x,i);
00031 }
00032 void compute(mdp_matrix GAMMA, mdp_matrix ZETA) {
00033 int A[KS_NDIM], B[KS_NDIM];
00034 mdp_matrix G1, G2, ZETA_H;
00035 site x(lattice());
00036 int i,mu,nu;
00037 ZETA_H=hermitian(ZETA);
00038 forallsites(x) {
00039 G1=mdp_identity(KS_NDIM);
00040 for(mu=1; mu<=KS_NDIM; mu++) {
00041 nu=mu % KS_NDIM;
00042 A[nu]=x(nu) % 2;
00043 if(A[nu]!=0) G1=-1.0*Gamma[nu]*G1;
00044 }
00045 for(i=0; i<16; i++) {
00046 B[0]=(i >> 3) & 0x1;
00047 B[1]=(i >> 2) & 0x1;
00048 B[2]=(i >> 1) & 0x1;
00049 B[3]=(i >> 0) & 0x1;
00050 G2=mdp_identity(KS_NDIM);
00051 for(mu=1; mu<=KS_NDIM; mu++) {
00052 nu=mu % KS_NDIM;
00053 if(B[nu]!=0)
00054 G2=G2*Gamma[nu];
00055 }
00056 *address(x,i)=(int) (0.25*real(trace(G1*GAMMA*G2*ZETA_H)));
00057 }
00058 }
00059 }
00060 };
00061
00062
00063 void operator_staggered_meson(staggered_field &out,
00064 staggered_field &in,
00065 phase_field &phases,
00066 gauge_field &U) {
00067
00068 site x(U.lattice());
00069 site x1(U.lattice());
00070 site x2(U.lattice());
00071 site y(U.lattice());
00072 int i, c, mu, A[KS_NDIM], B[KS_NDIM], P[4][2];
00073 int d, i0, i1, i2, i3;
00074 mdp_matrix paths, path;
00075 path=mdp_identity(U.nc);
00076 paths.dimension(U.nc,U.nc);
00077 paths=0;
00078
00079 forallsites(x) {
00080 for(c=0; c<U.nc; c++) out(x,c)=0;
00081 for(mu=0; mu<KS_NDIM; mu++)
00082 A[mu]=x(mu) % 2;
00083 for(i=0; i<16; i++) if(phases(x,i)!=0) {
00084 B[0]=(i >> 3) & 0x1;
00085 B[1]=(i >> 2) & 0x1;
00086 B[2]=(i >> 1) & 0x1;
00087 B[3]=(i >> 0) & 0x1;
00088 y.set(x(0)-A[0]+B[0],x(1)-A[1]+B[1], x(2)-A[2]+B[2],x(3)-A[3]+B[3]);
00089
00090 if(phases(x,i)!=0) {
00091 for(d=0, mu=0; mu<KS_NDIM; mu++)
00092 if(B[mu]-A[mu]!=0) {
00093 P[d][0]=B[mu]-A[mu];
00094 P[d][1]=mu;
00095 d++;
00096 }
00097 if(d==0) paths=mdp_identity(U.nc);
00098 if(d==1) paths=U(x,P[0][0],P[0][1]);
00099 if(d==2) {
00100 paths=
00101 U(x,P[0][0],P[0][1])*U(x.hop(P[0][0],P[0][1]),P[1][0],P[1][1])+
00102 U(x,P[1][0],P[1][1])*U(x.hop(P[1][0],P[1][1]),P[0][0],P[0][1]);
00103 paths/=2;
00104 }
00105 if(d==3) {
00106 x1=x.hop(P[0][0],P[0][1]);
00107 paths=U(x,P[0][0],P[0][1])*U(x1,P[1][0],P[1][1])*U(x1.hop(P[1][0],P[1][1]),P[2][0],P[2][1]);
00108 x1=x.hop(P[0][0],P[0][1]);
00109 paths+=U(x,P[0][0],P[0][1])*U(x1,P[2][0],P[2][1])*U(x1.hop(P[2][0],P[2][1]),P[1][0],P[1][1]);
00110 x1=x.hop(P[1][0],P[1][1]);
00111 paths+=U(x,P[1][0],P[1][1])*U(x1,P[0][0],P[0][1])*U(x1.hop(P[0][0],P[0][1]),P[2][0],P[2][1]);
00112 x1=x.hop(P[1][0],P[1][1]);
00113 paths+=U(x,P[1][0],P[1][1])*U(x1,P[2][0],P[2][1])*U(x1.hop(P[2][0],P[2][1]),P[0][0],P[0][1]);
00114 x1=x.hop(P[2][0],P[2][1]);
00115 paths+=U(x,P[2][0],P[2][1])*U(x1,P[0][0],P[0][1])*U(x1.hop(P[0][0],P[0][1]),P[1][0],P[1][1]);
00116 x1==x.hop(P[2][0],P[2][1]);
00117 paths+=U(x,P[2][0],P[2][1])*U(x1,P[1][0],P[1][1])*U(x1.hop(P[1][0],P[1][1]),P[0][0],P[0][1]);
00118 paths/=(3*2);
00119 }
00120 if(d==4) {
00121 paths=0;
00122 for(i0=0; i0<d; i0++)
00123 for(i1=0; i1<d; i1++) if(i1!=i0)
00124 for(i2=0; i2<d; i2++) if(i2!=i0 && i2!=i1)
00125 for(i3=0; i3<d; i3++) if(i3!=i0 && i3!=i1 && i3!=i2) {
00126 x1=x.hop(P[i0][0],P[i0][1]);
00127 x2=x1.hop(P[i1][0],P[i1][1]);
00128 path=
00129 U(x,P[i0][0],P[i0][1])*
00130 U(x1,P[i1][0],P[i1][1])*
00131 U(x2,P[i2][0],P[i2][1])*
00132 U(x2.hop(P[i2][0],P[i2][1]),P[i3][0],P[i3][1]);
00133 paths+=path;
00134 }
00135 paths/=(4*3*2);
00136 }
00137 out(x)+=phases(x,i)*paths*in(y);
00138 }
00139 }
00140 }
00141 }
00142
00143 const int local_source=1;
00144 const int wall_source =2;
00145
00146 mdp_matrix make_meson(gauge_field &U, gauge_field &V,
00147 mdp_matrix GAMMA,
00148 mdp_matrix ZETA,
00149 coefficients &coeff1,
00150 coefficients &coeff2,
00151 int source1_type=wall_source,
00152 int source2_type=wall_source & local_source,
00153 mdp_real precision =1e-7) {
00154
00155 if(source1_type!=wall_source && source1_type!=local_source)
00156 error("fermiqcd_staggered_mesons/make_meson: source option not implemented");
00157 if(source2_type!=local_source)
00158 error("fermiqcd_staggered_mesons/make_meson: source option not implemented");
00159
00160 int t, t_source, nsources=1, sourcestep=1;
00161 int i, j, nt=U.lattice().size(0);
00162 int nc=U.nc;
00163 mdp_complex c;
00164 site x(U.lattice());
00165 staggered_field tmp(U.lattice(),nc);
00166 staggered_field quark_source(U.lattice(),nc);
00167 staggered_field quark_prop(U.lattice(),nc);
00168 staggered_field anti_prop(U.lattice(),nc);
00169
00170 phase_field phases(U.lattice());
00171 phases.compute(GAMMA, ZETA);
00172
00173 mdp_matrix prop(1,nt);
00174 prop=0;
00175
00176 U.update();
00177 V.update();
00178
00179 for(t_source=0; t_source<nsources*sourcestep; t_source+=sourcestep) {
00180 for(i=0; i<U.nc; i++) {
00181 forallsites(x)
00182 for(j=0; j<U.nc; j++)
00183 if(i==j &&
00184 (x(0)==t_source || x(0)==(t_source+1)) &&
00185 source1_type==wall_source)
00186 quark_source(x,j)=mdp_complex(1,0);
00187 else if(i==j &&
00188 (x(0)==t_source || x(0)==(t_source+1)) &&
00189 (x(1)/2==0) && (x(2)/2==0) && (x(3)/2==0) &&
00190 source1_type==local_source)
00191 quark_source(x,j)=mdp_complex(1,0);
00192 else
00193 quark_source(x,j)=mdp_complex(0,0);
00194
00195 quark_source.update();
00196 mul_invQ(quark_prop,quark_source, V,coeff1,precision);
00197 quark_prop.update();
00198 operator_staggered_meson(tmp, quark_source, phases, U);
00199 tmp.update();
00200 mul_invQ(anti_prop, tmp, V,coeff2,precision);
00201 quark_prop.update();
00202 operator_staggered_meson(tmp, quark_prop, phases, U);
00203 forallsites(x) {
00204 for(j=0, c=mdp_complex(0,0); j<U.nc; j++)
00205 c+=conj(anti_prop(x,j))*tmp(x,j);
00206 t=(x(0)-t_source+nt) % nt; if(t%2==1) t--;
00207 prop(0,t)+=c;
00208 }
00209 }
00210 }
00211 return prop;
00212 }
00213
00214
00215 mdp_matrix GoldstonBoson_5x5(gauge_field &U,
00216 gauge_field &V,
00217 coefficients &coeff,
00218 float precision=1e-6)
00219 {
00220
00221 int i,j;
00222 unsigned int t_source=0;
00223 unsigned int t, nt
00224 =U.lattice().size(0);
00225 mdp_matrix tmp(nt,U.nc);
00226 mdp_matrix prop(2, nt);
00227 site x(U.lattice());
00228
00229
00230
00231 staggered_field quark_source(U.lattice(), U.nc);
00232 staggered_field quark_prop(U.lattice(), U.nc);
00233
00234
00235
00236 prop=0;
00237
00238
00239 printf("Starting to make 5x5 pion...\n");
00240 for(i=0; i<U.nc; i++) {
00241
00242 forallsites(x)
00243 for(j=0; j<U.nc; j++)
00244 if(i==j && x(0)-t_source==0)
00245 quark_source(x,j)=mdp_complex(1,0);
00246 else
00247 quark_source(x,j)=mdp_complex(0,0);
00248 quark_source.update();
00250
00251
00252 mul_invQ(quark_prop,quark_source,V,coeff,1e-5);
00253
00254
00255
00256 tmp=0;
00257 forallsites(x) {
00258 t=(x(0)-t_source+nt) % nt;
00259 for(j=0; j<U.nc; j++)
00260 tmp(t,j)+=conj(quark_prop(x,j))*quark_prop(x,j);
00261 }
00262 mpi.add(tmp);
00263 for(t=0; t<nt; t++)
00264 for(j=0; j<U.nc; j++)
00265 prop(0,t)+=tmp(t,j);
00266
00267
00268
00269 tmp=0;
00270 forallsites(x) {
00271 t=(x(0)-t_source+nt) % nt;
00272 for(j=0; j<U.nc; j++)
00273 tmp(t,j)+=quark_prop(x,j);
00274 }
00275 mpi.add(tmp);
00276 for(t=0; t<nt; t++)
00277 for(j=0; j<U.nc; j++)
00278 prop(1,t)+=conj(tmp(t,j))*tmp(t,j);
00279
00280 }
00281
00282 int volume=U.lattice().size(1)*
00283 U.lattice().size(2)*
00284 U.lattice().size(3);
00285 for(t=0; t<nt; t++) {
00286 prop(0,t)/=U.nc*volume;
00287 prop(1,t)/=U.nc*volume*volume;
00288 }
00289
00290 return prop;
00291 }
00292
00293 #endif
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417