00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00031 class em_field: public mdp_complex_field {
00032 public:
00033 int ndim, nc, nem;
00034 em_field() {
00035 ndim=nc=nem=0;
00036 reset_field();
00037 }
00038 em_field(mdp_lattice &a, int nc_) {
00039 reset_field();
00040 ndim=a.ndim;
00041 nc=nc_;
00042 nem=(ndim*(ndim-1))/2;
00043 allocate_field(a, nem*nc*nc);
00044 }
00045 em_field(em_field &em) {
00046 reset_field();
00047 ndim=em.ndim;
00048 nc=em.nc;
00049 nem=(ndim*(ndim-1))/2;
00050 allocate_field(em.lattice(), nem*nc*nc);
00051 }
00052 void allocate_em_field(mdp_lattice &a, int nc_) {
00053 deallocate_field();
00054 ndim=a.ndim;
00055 nc=nc_;
00056 nem=(ndim*(ndim-1))/2;
00057 allocate_field(a, nem*nc*nc);
00058 }
00059
00060 inline int ordered_index(int mu, int nu) const {
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071 if(mu<nu) return nu+(mu*(2*ndim-mu-3))/2-1;
00072 if(mu>nu) return mu+(nu*(2*ndim-nu-3))/2-1+nem;
00073
00074 return -1;
00075 };
00076
00077 inline mdp_matrix operator()(site x, int mu, int nu) {
00078 #ifdef CHECK_ALL
00079 if(mu>=nu) error("em(x,mu,nu) for mu>=nu is not defined");
00080 #endif
00081 register int k=ordered_index(mu,nu);
00082 return mdp_matrix(address(x,k*nc*nc),nc,nc);
00083 }
00084 inline mdp_complex& operator()(site x, int mu, int nu, int i, int j) {
00085 #ifdef CHECK_ALL
00086 if(mu>=nu) error("em(x,mu,nu) for mu>=nu is not defined");
00087 #endif
00088 register int k=ordered_index(mu,nu);
00089 return *(address(x,(k*nc+i)*nc+j));
00090 }
00091 inline const mdp_complex& operator()(site x, int mu, int nu,
00092 int i, int j) const {
00093 #ifdef CHECK_ALL
00094 if(mu>=nu) error("em(x,mu,nu) for mu>=nu is not defined");
00095 #endif
00096 register int k=ordered_index(mu,nu);
00097 return *(address(x,(k*nc+i)*nc+j));
00098 }
00099
00100 };
00101
00121 class gauge_field: public mdp_complex_field {
00122 public:
00123
00124 em_field em;
00125 mdp_nmatrix_field long_links;
00126 mdp_field<long> i_jump;
00127 mdp_matrix_field swirls;
00128
00129 int ndim,nc;
00130 gauge_field() {
00131 reset_field();
00132 }
00133 gauge_field(const gauge_field &U) : mdp_complex_field() {
00134 ndim=U.ndim;
00135 nc=U.nc;
00136 mdp_complex_field::operator=(U);
00137 }
00138 void operator=(const gauge_field &U) {
00139 ndim=U.ndim;
00140 nc=U.nc;
00141 mdp_complex_field::operator=(U);
00142 }
00143 gauge_field(mdp_lattice &a, int nc_) {
00144 reset_field();
00145 ndim=a.ndim;
00146 nc=nc_;
00147 allocate_field(a, a.ndim*nc*nc);
00148 }
00149 void allocate_gauge_field(mdp_lattice &a, int nc_) {
00150 deallocate_field();
00151 ndim=a.ndim;
00152 nc=nc_;
00153 allocate_field(a, a.ndim*nc*nc);
00154 }
00155 inline mdp_matrix operator() (site x, int mu) {
00156 #ifndef TWIST_BOUNDARY
00157 return mdp_matrix(address(x,mu*nc*nc),nc,nc);
00158 #else
00159 mdp_matrix tmp(address(x,mu*nc*nc),nc,nc);
00160 if(!in_block(x)) {
00161 mdp_matrix a;
00162 a=tmp;
00163 #ifdef ADVANCED_TWIST
00164 twist_eat_links(a,x,mu);
00165 #else
00166 twist_boundary(a,x);
00167 #endif
00168 return(a);
00169 }
00170 return tmp;
00171 #endif
00172 }
00173
00174 inline const mdp_matrix operator() (site x, int mu) const {
00175 #ifndef TWISTED_BOUNDARY
00176 return mdp_matrix(address(x,mu*nc*nc),nc,nc);
00177 #else
00178 mdp_matrix tmp(address(x,mu*nc*nc),nc,nc);
00179 if(!in_block(x)) {
00180 mdp_matrix a;
00181 a=tmp;
00182 #ifdef ADVANCED_TWIST
00183 twist_eat_links(a,x,mu);
00184 #else
00185 twist_boundary(a,x);
00186 #endif
00187 return(a);
00188 }
00189 return tmp;
00190 #endif
00191 }
00192
00193 inline mdp_complex& operator() (site x, int mu, int i, int j) {
00194 #ifdef TWISTED_BOUNDARY
00195 if(!in_block(x))
00196 error("call to &operator=(site x) per x out of block");
00197 #endif
00198 return *(address(x,(mu*nc+i)*nc+j));
00199 }
00200
00201 inline const mdp_complex& operator() (site x, int mu, int i, int j) const {
00202 #ifdef TWISTED_BOUNDARY
00203 if(!in_block(x))
00204 error("call to &operator=(site x) per x out of block");
00205 #endif
00206 return *(address(x,(mu*nc+i)*nc+j));
00207 }
00208
00209 inline mdp_matrix operator() (site x, int sign, int mu) {
00210 mdp_matrix tmp;
00211 if(sign==+1) tmp=(*this)(x,mu);
00212 if(sign==-1) tmp=hermitian((*this)(x-mu,mu));
00213 return tmp;
00214 }
00215 inline const mdp_matrix operator() (site x, int sign, int mu) const {
00216 mdp_matrix tmp;
00217 if(sign==+1) tmp=(*this)(x,mu);
00218 if(sign==-1) tmp=hermitian((*this)(x-mu,mu));
00219 return tmp;
00220 }
00221
00222
00223 #ifndef TWISTED_BOUNDARY
00224 inline const mdp_complex operator() (site x, int sign, int mu,
00225 int i, int j) const {
00226 if(sign==+1) return *(address(x,(mu*nc+i)*nc+j));
00227 if(sign==-1) return conj(*(address(x-mu,(mu*nc+j)*nc+i)));
00228 error("call to U(x,0,mu,i,j)");
00229 return mdp_complex(0,0);
00230 };
00231 #endif
00232
00233
00234
00235
00236
00237
00238
00239
00240 #ifdef TWISTED_BOUNDARY
00241 inline friend void twist_boundary(mdp_matrix &M, site &x) {
00242 static int mu, block;
00243 static mdp_complex z=exp(mdp_complex(0,2.0*Pi/3.0));
00244 static mdp_complex a,b,c,d,e,f,g,h,i;
00245
00246 for(mu=1; mu<x.lattice().ndim; mu++) {
00247 block=x.block[mu];
00248 if(block!=0) {
00249 a=M(0,0); b=M(0,1); c=M(0,2);
00250 d=M(1,0); e=M(1,1); f=M(1,2);
00251 g=M(2,0); h=M(2,1); i=M(2,2);
00252 if(mu==1 && block==1) {
00253 M(0,0)=e; M(0,1)=f; M(0,2)=d;
00254 M(1,0)=h; M(1,1)=i; M(1,2)=g;
00255 M(2,0)=b; M(2,1)=c; M(2,2)=a;
00256 }
00257 if(mu==1 && block==-1) {
00258 M(0,0)=i; M(0,1)=g; M(0,2)=h;
00259 M(1,0)=c; M(1,1)=a; M(1,2)=b;
00260 M(2,0)=f; M(2,1)=d; M(2,2)=e;
00261 }
00262 if(mu==2 && block==1) {
00263 M(0,0)=a; M(0,1)=b/z; M(0,2)=c/z/z;
00264 M(1,0)=d*z; M(1,1)=e; M(1,2)=f/z;
00265 M(2,0)=g*z*z; M(2,1)=h*z; M(2,2)=i;
00266 }
00267 if(mu==2 && block==-1) {
00268 M(0,0)=a; M(0,1)=b*z; M(0,2)=c*z*z;
00269 M(1,0)=d/z; M(1,1)=e; M(1,2)=f*z;
00270 M(2,0)=g/z/z; M(2,1)=h/z; M(2,2)=i;
00271 }
00272 if(mu==3 && block==1) {
00273 M(0,0)=i; M(0,1)=g/z; M(0,2)=h/z/z;
00274 M(1,0)=c*z; M(1,1)=a; M(1,2)=b/z;
00275 M(2,0)=f*z*z; M(2,1)=d*z; M(2,2)=e;
00276 }
00277 if(mu==3 && block==-1) {
00278 M(0,0)=e; M(0,1)=f*z; M(0,2)=d*z*z;
00279 M(1,0)=h/z; M(1,1)=i; M(1,2)=g*z;
00280 M(2,0)=b/z/z; M(2,1)=c/z; M(2,2)=a;
00281 }
00282 if(block*block>1) error("two blocks out\n");
00283 }
00284 }
00285 }
00286
00287 friend void define_twist_matrices() {
00288 begin_function("gauge_field__define_twist_matrices");
00289 mdp_complex z=exp(mdp_complex(0, 2.0*Pi/3.0));
00290 OmegaTwist[0]=identity(3);
00291 OmegaTwist[1]=zero(3);
00292 OmegaTwist[1](0,1)=1;
00293 OmegaTwist[1](1,2)=1;
00294 OmegaTwist[1](2,0)=1;
00295 OmegaTwist[2]=zero(3);
00296 OmegaTwist[2](0,0)=1.0/z;
00297 OmegaTwist[2](1,1)=1;
00298 OmegaTwist[2](2,2)=z;
00299 OmegaTwist[3]=zero(3);
00300 OmegaTwist[3](0,2)=1.0/z;
00301 OmegaTwist[3](1,0)=1;
00302 OmegaTwist[3](2,1)=z;
00303 end_function("gauge_field__define_twist_matrices");
00304 }
00305
00306 inline friend void twist_eat_fields(mdp_matrix &M, site &x,
00307 gauge_field &omega) {
00308 begin_function("gauge_field__twist_eat_matrices");
00309 static int mu, block;
00310 static mdp_complex z=exp(mdp_complex(0,2.0*Pi/3.0));
00311 static mdp_complex a,b,c,d,e,f,g,h,i;
00312
00313 for(mu=1; mu<x.lattice().ndim; mu++) {
00314 block=x.block[mu];
00315 if(block!=0) {
00316 if(mu==1 && block==1)
00317 M=omega(x,mu)*M*hermitian(omega(x,mu));
00318 if(mu==1 && block==-1)
00319 M=hermitian(omega(x,mu))*M*omega(x,mu);
00320 if(mu==2 && block==1)
00321 M=omega(x,mu)*M*hermitian(omega(x,mu));
00322 if(mu==2 && block==-1)
00323 M=hermitian(omega(x,mu))*M*omega(x,mu);
00324 if(mu==3 && block==1)
00325 M=omega(x,mu)*M*hermitian(omega(x,mu));
00326 if(mu==3 && block==-1)
00327 M=hermitian(omega(x,mu))*M*omega(x,mu);
00328 if(block*block>1) error("two blocks out\n");
00329 }
00330 }
00331 end_function("gauge_field__twist_eat_matrices");
00332 }
00333 #endif
00334 };
00335
