fermiqcd: /Users/mdipierro/fermiqcd/development/Libraries/fermiqcd_staggered

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       // only for goldstone pion if(x!=y || phases(x,i)!=1) error("what!");
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;  // 2^0
00144 const int wall_source =2;  // 2^1
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, // input gauge field
00216                          gauge_field &V, // input fat and link links
00217                          coefficients &coeff,    // input quark mass
00218                          float precision=1e-6)        // color source index
00219 {
00220   // // Local variable /////////////  
00221   int i,j;                        //
00222   unsigned int t_source=0;        // location of the source (timeslice)
00223   unsigned int t, nt              // number of timeslices
00224     =U.lattice().size(0);         //
00225   mdp_matrix tmp(nt,U.nc);            // auxiliary var
00226   mdp_matrix prop(2, nt);             // output vector
00227   site x(U.lattice());            // 
00228   // ///////////////////////////////
00229 
00230   // // Local fields ////////////////////////////////
00231   staggered_field quark_source(U.lattice(), U.nc); //
00232   staggered_field quark_prop(U.lattice(), U.nc);   //
00233   // ////////////////////////////////////////////////
00234 
00235   // // Initialization of parameters ////////////////
00236   prop=0;                                          //
00237   // ////////////////////////////////////////////////
00238 
00239   printf("Starting to make 5x5 pion...\n");
00240   for(i=0; i<U.nc; i++) {
00241     // // Make wall source /////////////////////////////
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     // // Make propagator - antipropagator - sink //////
00252     mul_invQ(quark_prop,quark_source,V,coeff,1e-5);   //
00253     // /////////////////////////////////////////////////
00254     
00255     // // trace with local sink ////////////////////////
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     // // trace with wall sink /////////////////////////
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   // // Normalize output and retrun //////////////////
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 /* EXAMPLE:
00296 
00297 // #define BLOCKSITE 4
00298 // #define TWISTED_BOUNDARY
00299 
00300 #include "/home/mdp/mdpqcd2/MDP_Lib2.h"
00301 #include "/home/mdp/mdpqcd2/MDP_MPI.h"
00302 #include "/home/mdp/mdpqcd2/MDP_Prompt.h"
00303 #include "/home/mdp/mdpqcd2/MDP_Gauge.h"
00304 #include "/home/mdp/mdpqcd2/MDP_iGauge.h"
00305 #include "/home/mdp/mdpqcd2/MDP_GaugeFix.h"
00306 #include "/home/mdp/mdpqcd2/MDP_Fermi.h"
00307 #include "/home/mdp/mdpqcd2/MDP_Staggered.h"
00308 #include "/home/mdp/mdpqcd2/MDP_StaggeredMesons.h"
00309 
00310 int main(int argc, char **argv) {
00311   mpi.open_wormholes(argc, argv);
00312   define_base_matrices("FERMILAB");
00313 
00314  
00315   //  define_twist_matrices();
00316 
00317   printf("COMPUTING C2(t) FOR STAGGERED MESON: %s\n", argv[1]);
00318   
00319   float seed=0;
00320   int t, Nt=24, Ns=8;
00321   int Nc=3;
00322   int box[]={Nt, Ns, Ns, Ns};
00323 
00324   float pU, u0=0.8629;
00325   float *c;
00326   char s[256];
00327 
00328   mdp_lattice space_time(4,box,
00329                              default_partitioning<0>,
00330                              torus_topology,
00331                              seed,3);
00332 
00333   int conf, nconf=(int) val(argv[2]);
00334   
00335   gauge_field     U(space_time,Nc);
00336   gauge_field     V(space_time,Nc);
00337   eta_field       eta(space_time);
00338   site x(space_time);
00339   mdp_matrix prop;
00340   float mass_a, mass_b;
00341 
00342   mass_a=0.03;
00343   mass_b=0.03;
00344 
00345   default_staggered_action=staggered_mul_Q_TWO;
00346   default_staggered_inversion_method=staggered_BiCG_UML;
00347 
00348   u0=0.8629;
00349 
00350   JackBoot c2(nconf, Nt);
00351 
00352   for(conf=0; conf<nconf; conf++) {
00353     c2.conf=conf;
00354     sprintf(s, "/home/mdp/data/gauge_improved_b7.4_u0.8629/gauge32x08_b7.4_u0.8629_n%.6i", conf+1);
00355  
00356     U.load(s);
00357     sprintf(s, "/home/mdp/data/gauge_improved_b7.4_u0.8629/gauge24x08_b7.4_u0.8629_n%.6i", conf+1);
00358     U.save(s);
00359 
00360     pU=pow(u0,4);
00361     c=lepage_coeffiecients(pU, "Full");
00362     set_antiperiodic_phases(U,0,true);
00363     lepage_improved_links(V,U,c);
00364 
00365     if(strcmp(argv[1],"5x5")==0) 
00366       prop= make_meson(U,V, Gamma5, Gamma5,mass_a, mass_b, 
00367                        wall_source, wall_source, 1e-5);
00368     if(strcmp(argv[1],"05x05")==0) 
00369       prop= make_meson(U,V, Gamma[0]*Gamma5, Gamma[0]*Gamma5,mass_a, mass_b, 
00370                        wall_source, wall_source, 1e-5);
00371     if(strcmp(argv[1],"5x35")==0) 
00372       prop= make_meson(U,V, Gamma5, Gamma[3]*Gamma5,mass_a, mass_b, 
00373                        wall_source, wall_source, 1e-5);
00374     if(strcmp(argv[1],"05x12")==0) 
00375       prop= make_meson(U,V, Gamma[0]*Gamma5, Gamma[1]*Gamma[2],mass_a, mass_b, 
00376                        wall_source, wall_source, 1e-5);
00377     if(strcmp(argv[1],"5x03")==0) 
00378       prop= make_meson(U,V, Gamma5, Gamma[0]*Gamma[3],mass_a, mass_b, wall_source, 
00379                        wall_source, 1e-5);
00380     if(strcmp(argv[1],"05x3")==0) 
00381       prop= make_meson(U,V, Gamma[0]*Gamma5, Gamma[3],mass_a, mass_b, 
00382                        wall_source, wall_source, 1e-5);
00383     if(strcmp(argv[1],"5x0")==0) 
00384       prop= make_meson(U,V, Gamma5, Gamma[0],mass_a, mass_b, 
00385                        wall_source, wall_source, 1e-5);
00386     if(strcmp(argv[1],"05xI")==0) 
00387       prop= make_meson(U,V, Gamma[0]*Gamma5, Gamma1,mass_a, mass_b, 
00388                        wall_source, wall_source, 1e-5);
00389 
00390     for(t=0; t<Nt; t++) c2(t)=real(prop(0,t));
00391 
00392     float x, dx;
00393 
00394     if(ME==0) {
00395       for(t=0; t<Nt; t=t+2) {
00396         c2.plain(t);
00397         x=c2.mean();
00398         dx=c2.j_err();
00399         printf("%i,\t%e,\t%e\n", t, log(x), dx/x);
00400       }
00401       printf("\n");
00402     }
00403   }
00404   mpi.close_wormholes();
00405   return 0;
00406 }
00407 
00408 
00409 
00410 
00411 NOTES:
00412 make meson dow not work is you enable twist. The rest does work!
00413 
00414 
00415 */
00416 
00417
/Users/mdipierro/fermiqcd/development/Libraries/fermiqcd_staggered_mesons.h
