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
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
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