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

fermiqcd_gauge_fixing.h

Go to the documentation of this file.
00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00014 class gaugefixing_stats {
00015  public:
00016   uint max_steps;
00017   mdp_real target_precision;
00018   uint steps;
00019   mdp_real precision;
00020   mdp_real action;
00021 };
00022 
00033 class GaugeFixing {
00034  public:
00035   static const int Coulomb = 0;
00036   static const int Landau = 10;
00037     
00038   static void hit(gauge_field &U,
00039                   int mu, 
00040                   int parity, 
00041                   int i, int j,
00042                   mdp_real overrelaxation_boost = 1) {
00043 
00044   // This also works for twisted boundary even if I use generic_field<> */
00045   
00046     static mdp_real a0,a1,a2,a3,b,c,d;
00047     static mdp_real a0_sq,ai_sq;
00048     static mdp_complex x0,x1;
00049     int k,nu,nc=U.nc;
00050     int opposite_parity=EVENODD;
00051     generic_field<mdp_complex> W(U.lattice(),4);
00052     mdp_matrix U_up(nc,nc), U_dw(nc,nc);
00053     mdp_matrix A;
00054     site x(U.lattice());
00055     site y(U.lattice());
00056     switch(parity) {
00057     case EVEN: opposite_parity=ODD; break;
00058     case ODD:  opposite_parity=EVEN; break;
00059     }
00060     forallsitesofparity(x,parity) {
00061       a0=a1=a2=a3=0;
00062       for(nu=0; nu<U.ndim; nu++) if(nu!=mu) {
00063         U_up=U(x,nu);
00064         U_dw=U(x-nu,nu);
00065         a0+=
00066           +real(U_dw(i,i))+real(U_up(i,i))
00067           +real(U_dw(j,j))+real(U_up(j,j));
00068         a1+=
00069           +imag(U_dw(j,i))-imag(U_up(j,i))
00070           +imag(U_dw(i,j))-imag(U_up(i,j));
00071         a2+=
00072           -real(U_dw(j,i))+real(U_up(j,i))
00073           +real(U_dw(i,j))-real(U_up(i,j));
00074         a3+=
00075           +imag(U_dw(i,i))-imag(U_up(i,i))
00076           -imag(U_dw(j,j))+imag(U_up(j,j));
00077       }
00078       ai_sq=a1*a1+a2*a2+a3*a3;
00079       a0_sq=a0*a0;
00080       b=(overrelaxation_boost*a0_sq+ai_sq)/(a0_sq+ai_sq);
00081       c=sqrt(a0_sq+b*b*ai_sq);
00082       d=b/c;
00083       a0/=c;
00084       a1*=d;
00085       a2*=d;
00086       a3*=d;
00087       W(x,0)=mdp_complex(a0,a3);
00088       W(x,1)=mdp_complex(a2,a1);
00089       W(x,2)=mdp_complex(-a2,a1);
00090       W(x,3)=mdp_complex(a0,-a3);
00091     }
00092     W.update();
00093     forallsitesofparity(x,parity)
00094       for(nu=0; nu<U.ndim; nu++) 
00095         for(k=0; k<U.nc; k++) {
00096           x0=U(x,nu,i,k);
00097           x1=U(x,nu,j,k);
00098           U(x,nu,i,k)=W(x,0)*x0+W(x,1)*x1;
00099           U(x,nu,j,k)=W(x,2)*x0+W(x,3)*x1;
00100         }
00101     forallsitesofparity(x,opposite_parity)
00102       for(nu=0; nu<U.ndim; nu++) {
00103         y=x+nu;
00104 #ifndef TWISTED_BOUNDARY
00105         for(k=0; k<U.nc; k++) {
00106           x0=U(x,nu,k,i);
00107           x1=U(x,nu,k,j);
00108           U(x,nu,k,i)=conj(W(y,0))*x0+conj(W(y,1))*x1;
00109           U(x,nu,k,j)=conj(W(y,2))*x0+conj(W(y,3))*x1;
00110         }
00111 #else
00112         if(in_block(y)) 
00113           for(k=0; k<U.nc; k++) {
00114             x0=U(x,nu,k,i);
00115             x1=U(x,nu,k,j);
00116             U(x,nu,k,i)=conj(W(y,0))*x0+conj(W(y,1))*x1;
00117             U(x,nu,k,j)=conj(W(y,2))*x0+conj(W(y,3))*x1;
00118           } else {
00119             A=identity(nc);
00120             A(i,i)=W(y,0);
00121             A(i,j)=W(y,1);
00122             A(j,i)=W(y,2);
00123             A(j,j)=W(y,3);
00124             twist_boundary(A,y);
00125             U(x,nu)=U(x,nu)*hermitian(A);
00126           }
00127 #endif
00128       }
00129     U.update();  
00130   }
00131 
00132   static void z3_fix(gauge_field &U, int mu) {
00133     int     i=0,t;
00134     site    x(U.lattice());
00135     mdp_matrix  A;
00136     mdp_complex phase[3]= {mdp_complex(1,0), 
00137                        exp(2.0*Pi*I/3.0), 
00138                        exp(4.0*Pi*I/3.0) };
00139     mdp_real alpha[3];
00140     for(t=0; t<U.lattice().size(mu)-1; t++) {
00141       A=mdp_zero(U.nc);
00142       forallsites(x) if(x(mu)==t) A+=U(x,mu);
00143       mpi.add(A);
00144       alpha[0]=real(trace(A*conj(phase[0])));
00145       alpha[1]=real(trace(A*conj(phase[1])));
00146       alpha[2]=real(trace(A*conj(phase[2])));
00147       if(alpha[0]>=alpha[1] && alpha[0]>=alpha[2]) i=0;
00148       if(alpha[1]> alpha[0] && alpha[1]>=alpha[2]) i=1;
00149       if(alpha[2]> alpha[0] && alpha[2]>=alpha[1]) i=2;
00150       forallsites(x) 
00151         if(x(mu)==t) U(x,mu)*=conj(phase[i]);
00152         else if(x(mu)==t+1) U(x,mu)*=phase[i];
00153       U.update();
00154     }
00155   }
00164   static gaugefixing_stats fix(gauge_field &U, 
00165                   int mu=0, 
00166                   int max_steps=1, 
00167                   mdp_real target_precision=1e-5, 
00168                   mdp_real overrelaxation_boost = 1,
00169                   bool z3=false) {
00170     
00171     gaugefixing_stats stats;
00172     int step,nu,parity,i,j;
00173     site x(U.lattice());
00174     double action=0;
00175     double precision=0;
00176     mdp_matrix M(U.nc,U.nc);
00177 
00178     stats.max_steps=max_steps;
00179     stats.target_precision=target_precision;
00180 
00181     mdp << "step\taction\tprecision\n";
00182     
00183     for(step=0; step<max_steps; step++) {
00184       
00185       for(parity=EVEN; parity<=ODD; parity++) {
00186         for(i=0;     i<U.nc-1; i++)
00187           for(j=i+1; j<U.nc;   j++) {
00188             hit(U,mu,parity,i,j,overrelaxation_boost);
00189           }
00190       }   
00191       // ******** What is the gauge action ? *********
00192       action=0;
00193       precision=0;
00194       forallsites(x) {
00195         M=0;
00196         for(nu=0; nu<U.ndim; nu++) if(nu!=mu) {
00197           M+=U(x,nu)-U(x-nu,nu);
00198           action+=real(trace(U(x,nu)+U(x-nu,nu)));
00199         }
00200         M=(M-trace(M)/U.nc);
00201         M=M-hermitian(M);
00202         for(i=0; i<U.nc; i++)
00203         for(j=0; j<U.nc; j++)
00204           precision+=(double) pow(abs(M(i,j)),2);
00205       }
00206       mpi.add(precision);
00207       mpi.add(action);
00208 
00209       precision=sqrt(precision/(U.nc*U.nc*U.lattice().global_volume()));
00210       action=action/(2.0*U.nc*U.lattice().global_volume());
00211 
00212       mdp << step << "\t" << action << "\t" << precision << '\n';
00213       
00214       if(step !=0 && precision<target_precision) break;
00215       // ***********************************************
00216       
00217     }
00218     stats.steps=step;
00219     stats.precision=precision;
00220     stats.action=action;
00221 
00222     if(z3) z3_fix(U,mu);
00223 
00224     mdp << "steps=" << step
00225         << ", action=" << action
00226         << ", precison+" << precision << '\n';
00227 
00228     return stats;
00229   }
00230 };
00231 
00232 
00233 
00234 
00235 
00236     
00237 
00238 
00239 
00240 
00241 

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