Main Page | Class Hierarchy | Class List | File List | Class Members | File Members

fermiqcd_gauge_actions.h

Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00014 class gauge_stats {};
00015 
00016 
00029 class WilsonGaugeAction {
00030  public:
00031   static void heatbath_SU2(mdp_prng &random,
00032                            mdp_real beta_eff, mdp_complex *a) {
00033     mdp_real e0,e1,e2,e3, dk, p0;
00034     mdp_real r1,r2,r3,r4;
00035     mdp_real a0,a1,a2,a3;
00036     mdp_complex u0,u1,u2,u3;
00037     mdp_complex v0,v1,v2,v3;
00038     mdp_real delta, phi, sin_alpha, sin_theta, cos_theta;
00039     e0=real(a[0])+real(a[3]);
00040     e1=imag(a[1])+imag(a[2]);
00041     e2=real(a[1])-real(a[2]);
00042     e3=imag(a[0])-imag(a[3]);
00043     dk=sqrt(e0*e0+e1*e1+e2*e2+e3*e3);
00044     p0=(dk*beta_eff);
00045     u0=mdp_complex(e0/dk,-e3/dk);
00046     u2=mdp_complex(e2/dk,-e1/dk);
00047     u1=mdp_complex(-e2/dk,-e1/dk);
00048     u3=mdp_complex(e0/dk,e3/dk);
00049     
00050     if(beta_eff<=0.0) 
00051       error("fermiqcd_gauge_algorithms/heatbath_SU2: beta is zero");
00052     
00053     do {
00054       do; while ((r1 = random.plain()) < 0.0001);
00055       r1 = -log(r1)/p0;
00056       do; while ((r2 = random.plain()) < 0.0001);
00057       r2 = -log(r2)/p0;
00058       r3 = cos(2.0*Pi*random.plain());
00059       r3 = r3*r3;
00060       delta = r2+r1*r3;
00061       r4=random.plain();
00062     } while (r4*r4 > (1.0-0.5*delta));
00063     a0=1.0-delta;
00064     cos_theta=2.0*random.plain()-1.0;
00065     sin_theta=sqrt(1.0-cos_theta*cos_theta);
00066     sin_alpha=sqrt(1-a0*a0);
00067     phi=2.0*Pi*random.plain();
00068     a1=sin_alpha*sin_theta*cos(phi);
00069     a2=sin_alpha*sin_theta*sin(phi);
00070     a3=sin_alpha*cos_theta;
00071     v0=mdp_complex(a0,a3);
00072     v1=mdp_complex(a2,a1);
00073     v2=mdp_complex(-a2,a1);
00074     v3=mdp_complex(a0,-a3);
00075     a[0]=v0*u0+v1*u2;
00076     a[1]=v0*u1+v1*u3;
00077     a[2]=v2*u0+v3*u2;
00078     a[3]=v2*u1+v3*u3;
00079   }
00080   
00081   static gauge_stats heatbath(gauge_field &U, 
00082                               coefficients &coeff, 
00083                               int n_iter=1) {
00084     begin_function("WilsonGaugeAction__heatbath");
00085     if(U.nc==1) error("fermiqcd_gauge_algorithms/heatbath(): U(1)? (use metropolis)");
00086     gauge_stats stats;
00087     mdp_real beta, zeta;
00088     if(coeff.has_key("beta")) beta=coeff["beta"]; 
00089     else error("beta undeclared");
00090     if(coeff.has_key("zeta")) zeta=coeff["zeta"]; 
00091     else zeta=1;
00092 
00093     int i,j,k,iter,mu,parity;
00094     mdp_matrix M;
00095     mdp_complex a[4], tmpUik;
00096     site x(U.lattice());
00097     double time=mpi.time();
00098 
00099     mdp << coeff;
00100 
00101     for(iter=0; iter<n_iter; iter++) 
00102       for(parity=0; parity<2; parity++) 
00103         for(mu=0; mu<U.ndim; mu++) {
00104           forallsitesofparity(x,parity) {
00105             for(i=0; i<U.nc-1; i++)
00106               for(j=i+1; j<U.nc; j++) {
00107                 if(zeta==1)    M=U(x,mu)*staple_H(U,x,mu);
00108                 else if(mu==0) M=zeta*U(x,0)*staple_H(U,x,0);
00109                 else           M=((mdp_real) 1.0/zeta)*U(x,mu)*staple_H(U,x,mu);
00110                 a[0]=M(i,i); 
00111                 a[1]=M(i,j);
00112                 a[2]=M(j,i);
00113                 a[3]=M(j,j);
00114                 heatbath_SU2(U.lattice().random(x),beta/U.nc,a);
00115                 for(k=0; k<U.nc; k++) {
00116                   tmpUik=a[0]*U(x,mu,i,k)+a[1]*U(x,mu,j,k);
00117                   U(x,mu,j,k)=a[2]*U(x,mu,i,k)+a[3]*U(x,mu,j,k);
00118                   U(x,mu,i,k)=tmpUik;
00119                 }
00120               }
00121           }
00122           // The next command does all the communications!
00123           U.update(parity, mu, U.nc*U.nc);
00124         }
00125     mdp << "\t<stats>\n\t\t<time>" << mpi.time()-time << "</time>\n\t</stats>\n";
00126     end_function("WilsonGaugeAction__heatbath");
00127     return stats;
00128   }
00129 };
00130 
00159 class ImprovedGaugeAction : public WilsonGaugeAction {
00160  private:
00161   static mdp_matrix rectangles_0i_H(gauge_field &U, site x, int mu) {
00162     mdp_matrix tmp(U.nc,U.nc);
00163     mdp_matrix b1(U.nc, U.nc);
00164     mdp_matrix b2(U.nc, U.nc);
00165     mdp_matrix b3(U.nc, U.nc);
00166     site y0(U.lattice());
00167     site y1(U.lattice());
00168     site y2(U.lattice());
00169     int nu;
00170     tmp=0;
00171     if(mu==0) {
00172       for(nu=1; nu<U.ndim; nu++) {
00173         y0=x+mu;
00174         y1=y0+nu;
00175         tmp+=U(y0,nu)*U(y1,nu)*hermitian(U(x,nu)*U(x+nu,nu)*U((x+nu)+nu,mu));
00176         y0=(x-nu)-nu;
00177         y1=y0+mu;
00178         y2=y1+nu;
00179         tmp+=hermitian(U(y0,mu)*U(y1,nu)*U(y2,nu))*U(y0,nu)*U(x-nu,nu);
00180       }
00181     } else {
00182       nu=0;
00183       y0=(x-mu)+nu;
00184       tmp+=U(x+mu,nu)*hermitian(U(x-mu,nu)*U(y0,mu)*U(x+nu,mu))*U(x-mu,mu);
00185       y0=x-mu;
00186       y1=y0-nu;
00187       y2=y1+mu;
00188       tmp+=hermitian(U(y1,mu)*U(y2,mu)*U(y2+mu,nu))*U(y1,nu)*U(y0,mu);
00189       y0=x+mu;
00190       y1=y0+nu;
00191       y2=y1-mu;
00192       tmp+=U(y0,mu)*U(y0+mu,nu)*hermitian(U(x,nu)*U(y2,mu)*U(y1,mu));
00193       y0=((x+mu)+mu)-nu;
00194       y1=y0-mu;
00195       y2=y1-mu;
00196       tmp+=U(x+mu,mu)*hermitian(U(y2,mu)*U(y1,mu)*U(y0,nu))*U(y2,nu);
00197       
00198     }
00199     prepare(tmp);
00200     return tmp;
00201   }
00202   
00203 // if min_nu==0 then rectangles_ij computes all 6 rectanges
00204   
00205   static mdp_matrix rectangles_ij_H(gauge_field &U, site x, int mu, int min_nu=1) {
00206     mdp_matrix tmp(U.nc,U.nc);
00207     mdp_matrix b1(U.nc, U.nc);
00208     mdp_matrix b2(U.nc, U.nc);
00209     mdp_matrix b3(U.nc, U.nc);
00210     site y0(U.lattice());
00211     site y1(U.lattice());
00212     site y2(U.lattice());
00213     int nu;
00214     for(nu=min_nu; nu<U.ndim; nu++) if(nu!=mu) {    
00215       y0=x+mu;
00216       y1=y0+nu;
00217       tmp+=U(y0,nu)*U(y1,nu)*hermitian(U(x,nu)*U(x+nu,nu)*U((x+nu)+nu,mu));
00218       
00219       y0=(x-nu)-nu;
00220       y1=y0+mu;
00221       y2=y1+nu;
00222       tmp+=hermitian(U(y0,mu)*U(y1,nu)*U(y2,nu))*U(y0,nu)*U(x-nu,nu);
00223       
00224       y0=(x-mu)+nu;
00225       tmp+=U(x+mu,nu)*hermitian(U(x-mu,nu)*U(y0,mu)*U(x+nu,mu))*U(x-mu,mu);
00226       
00227       y0=x-mu;
00228       y1=y0-nu;
00229       y2=y1+mu;
00230       tmp+=hermitian(U(y1,mu)*U(y2,mu)*U(y2+mu,nu))*U(y1,nu)*U(y0,mu);
00231 
00232       y0=x+mu;
00233       y1=y0+nu;
00234       y2=y1-mu;
00235       tmp+=U(y0,mu)*U(y0+mu,nu)*hermitian(U(x,nu)*U(y2,mu)*U(y1,mu));
00236 
00237       y0=((x+mu)+mu)-nu;
00238       y1=y0-mu;
00239       y2=y1-mu;
00240       tmp+=U(x+mu,mu)*hermitian(U(y2,mu)*U(y1,mu)*U(y0,nu))*U(y2,nu);
00241     }
00242     prepare(tmp);
00243     return tmp;
00244   }
00245 
00246   // //////////////////////////////////////////////////////
00247   // this is slow but should make the chair correcly ...
00248   // see: hep-lat/0712010
00249   // //////////////////////////////////////////////////////
00250 
00251   static mdp_matrix chair_H(gauge_field &U, site x, int mu) {
00252     int ndim=U.ndim;
00253     int nu, rho;
00254     mdp_matrix tmp(U.nc, U.nc);
00255     mdp_matrix b1(U.nc, U.nc);
00256     mdp_matrix b2(U.nc, U.nc);
00257     mdp_matrix b3(U.nc, U.nc);
00258     site y1(U.lattice());
00259     site y2(U.lattice());
00260     site y3(U.lattice());
00261     site y4(U.lattice());
00262     site y5(U.lattice());
00263     tmp=0;
00264     for(nu=0; nu<ndim; nu++) if(nu!=mu)
00265       for(rho=0; rho<ndim; rho++) if((rho!=nu) && (rho!=mu)) { 
00266         y1=x+mu;
00267         y2=y1+nu;
00268         y3=y2+rho;
00269         y4=y3-mu;
00270         y5=y4-nu;
00271         tmp+=U(y1,nu)*U(y2,rho)*hermitian(U(x,rho)*U(y5,nu)*U(y4,mu));
00272         y1=x+mu;
00273         y2=y1-nu;
00274         y3=y2+rho;
00275         y4=y3-mu;
00276         y5=y4+nu;
00277         tmp+=hermitian(U(y2,nu))*U(y2,rho)
00278           *hermitian(U(y4,mu))*U(y4,nu)*hermitian(U(x,rho));
00279         
00280         y1=x+mu;
00281         y2=y1+nu;
00282         y3=y2-rho;
00283         y4=y3-mu;
00284         y5=y4-nu;
00285         tmp+=U(y1,nu)*hermitian(U(y5,nu)*U(y4,mu)*U(y3,rho))*U(y5,rho);
00286         
00287         y1=x+mu;
00288         y2=y1-nu;
00289         y3=y2-rho;
00290         y4=y3-mu;
00291         y5=y4+nu;
00292         tmp+=hermitian(U(y4,mu)*U(y3,rho)*U(y2,nu))*U(y4,nu)*U(y5,rho);
00293 
00294       }
00295     return(tmp);
00296   }
00297   
00298   // ////////////////////////////////////////////////////////////////////
00299   // new_heatbath uses an improved gauge action!
00300   // both isotropic (param.zeta=1) and anisotropic
00301   // ////////////////////////////////////////////////////////////////////
00302 
00303   enum { iGauge_min=4 };
00304 
00305   static int strange_mapping(site &x) {
00306     int mu, type=0;
00307     for(mu=0; mu<x.lattice().ndim; mu++) type+=(int) pow((float) iGauge_min,mu)*(x(mu) % iGauge_min);
00308     return type;
00309   }
00310   
00311  public:
00312   static gauge_stats heatbath(gauge_field &U,
00313                               coefficients &coeff,
00314                               int n_iter=1,
00315                               string model="MILC") { 
00316     
00317     begin_function("ImprovedGaugeAction__heatbath");
00318 
00319     gauge_stats stats;
00320     mdp_real beta, zeta, u_s, u_t;
00321 
00322     if(coeff.has_key("beta")) beta=coeff["beta"]; else error("beta undeclared");
00323     if(coeff.has_key("zeta")) zeta=coeff["zeta"]; else zeta=1;
00324     if(coeff.has_key("u_t"))  u_t=coeff["u_t"];   else u_t=1;
00325     if(coeff.has_key("u_s"))  u_s=coeff["u_s"];   else u_s=1;
00326     
00327     // if(Nproc!=1)  error("improved_heatbath() does not work in parallel!");
00328 
00329     if(U.ndim!=4) error("fermiqcd_gauge_algorithms/improved_heatbath(): ndim!=4 (use heatbath instead)");
00330 
00331     /*
00332       if((U.lattice().size(0) % 3!=0) || (U.lattice().size(1) % 3!=0) ||
00333       (U.lattice().size(2) % 3!=0) || (U.lattice().size(3) % 3!=0))
00334       error("lattice is not divisible by 3");
00335     */
00336 
00337     if(U.nc==1)   error("fermiqcd_gauge_algorithms/improved_heatbath(): U(1)? (use metropolis instead)");
00338     int nc=U.nc;
00339     int ndim=U.ndim;
00340     int i,j,k,iter,mu,type;
00341     mdp_matrix M;
00342     site x(U.lattice());
00343     double time=mpi.time();
00344     mdp_complex a[4],tmpUik;
00345     mdp_real alpha_s;
00346     mdp_real c_tp=0, c_tr=0, c_sp=0, c_sr=0, c_p=0, c_r=0, c_c=0;
00347     
00348     if(model=="Morningstar") {
00349       
00350       c_tp=4.0*zeta/(3.0*pow((double)u_s*u_t,(double)2.0));
00351       c_tr=-1.0*zeta/(12.0*pow((double)u_s,(double)4.0)*pow((double)u_t,(double)2.0));
00352       c_sp=5.0/(3.0*zeta*pow((double)u_s,(double)4.0));
00353       c_sr=-1.0/(12.0*zeta*pow((double)u_s,(double)6.0));
00354       
00355       c_p=1.0*pow((double)u_s,(double)-4.0);
00356       c_r=-0.05*pow((double)u_s,(double)-6.0);
00357       c_c=0;
00358       
00359     } else if(model=="MILC") {
00360       
00361       if(zeta!=1) 
00362         error("fermiqcd_gauge_algorithms/improved_heatbath: zeta!=1");
00363       
00364       alpha_s=-4.0*log(u_s)/3.0684;
00365       c_p=1.0;
00366       c_r=-0.05*pow((double)u_s,(double)-2.0)*(1.0+0.4805*alpha_s);;
00367       c_c=-1.00*pow((double)u_s,(double)-2.0)*(0.03325*alpha_s);;
00368       
00369     } else {
00370       mdp << "Using default non-improved action" << '\n';
00371       stats=WilsonGaugeAction::heatbath(U,coeff,n_iter);
00372       end_function("ImprovedGaugeAction__heatbath");
00373       return stats;
00374     }
00375     
00376     mdp << coeff;
00377     
00378     for(iter=0; iter<n_iter; iter++)
00379       for(type=0; type<(int) pow((float) iGauge_min, U.ndim); type++) {
00380         forallsites(x) {
00381           if(strange_mapping(x)==type) {
00382             for(mu=0; mu<ndim; mu++) 
00383               for(i=0; i<nc-1; i++)
00384                 for(j=i+1; j<nc; j++) { 
00385                   if(zeta!=1) {
00386                     // anisotropic San Diego action
00387                     if(mu==0) M=U(x,mu)*(c_tp*staple_0i_H(U,x,0)+
00388                                          c_tr*rectangles_0i_H(U,x,0));
00389                     else      M=U(x,mu)*(c_sp*staple_ij_H(U,x,mu)+
00390                                          c_tp*staple_0i_H(U,x,mu)+
00391                                          c_sr*rectangles_ij_H(U,x,mu)+
00392                                          c_tr*rectangles_0i_H(U,x,mu));
00393                   } else if(c_c==0) {
00394                     // isotropic San Diego action
00395                     M=U(x,mu)*(c_p*staple_H(U,x,mu)+
00396                                c_r*rectangles_ij_H(U,x,mu,0));
00397                   } else {
00398                     // fully O(a^2) improved isotropic MILC action
00399                     M=U(x,mu)*(c_p*staple_H(U,x,mu)+
00400                                c_r*rectangles_ij_H(U,x,mu,0)+
00401                                c_c*chair_H(U,x,mu));
00402                   }
00403                   a[0]=M(i,i); 
00404                   a[1]=M(i,j);
00405                   a[2]=M(j,i);
00406                   a[3]=M(j,j);
00407                   heatbath_SU2(U.lattice().random(x),beta/U.nc,a);
00408                   for(k=0; k<U.nc; k++) {
00409                     tmpUik=a[0]*U(x,mu,i,k)+a[1]*U(x,mu,j,k);
00410                     U(x,mu,j,k)=a[2]*U(x,mu,i,k)+a[3]*U(x,mu,j,k);
00411                     U(x,mu,i,k)=tmpUik;
00412                   }
00413                 }
00414           }
00415         }
00416         U.update();
00417       }
00418     mdp << "\t<stats>\n\t\t<time>" << mpi.time()-time << "</time>\n\t</stats>\n";
00419     end_function("ImprovedGaugeAction__heatbath");
00420     return stats;
00421   }
00422 };
00423 

Generated on Sun Feb 27 15:12:18 2005 by  doxygen 1.4.1