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<mdp_int> 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(U) {
00134 ndim=U.ndim;
00135 nc=U.nc;
00136 }
00137 void operator=(const gauge_field &U) {
00138 ndim=U.ndim;
00139 nc=U.nc;
00140 mdp_complex_field::operator=(U);
00141 }
00142 gauge_field(mdp_lattice &a, int nc_) {
00143 reset_field();
00144 ndim=a.ndim;
00145 nc=nc_;
00146 allocate_field(a, a.ndim*nc*nc);
00147 }
00148 void allocate_gauge_field(mdp_lattice &a, int nc_) {
00149 deallocate_field();
00150 ndim=a.ndim;
00151 nc=nc_;
00152 allocate_field(a, a.ndim*nc*nc);
00153 }
00154 inline mdp_matrix operator() (site x, int mu) {
00155 #ifndef TWIST_BOUNDARY
00156 return mdp_matrix(address(x,mu*nc*nc),nc,nc);
00157 #else
00158 mdp_matrix tmp(address(x,mu*nc*nc),nc,nc);
00159 if(!in_block(x)) {
00160 mdp_matrix a;
00161 a=tmp;
00162 #ifdef ADVANCED_TWIST
00163 twist_eat_links(a,x,mu);
00164 #else
00165 twist_boundary(a,x);
00166 #endif
00167 return(a);
00168 }
00169 return tmp;
00170 #endif
00171 }
00172
00173 inline const mdp_matrix operator() (site x, int mu) const {
00174 #ifndef TWISTED_BOUNDARY
00175 return mdp_matrix(address(x,mu*nc*nc),nc,nc);
00176 #else
00177 mdp_matrix tmp(address(x,mu*nc*nc),nc,nc);
00178 if(!in_block(x)) {
00179 mdp_matrix a;
00180 a=tmp;
00181 #ifdef ADVANCED_TWIST
00182 twist_eat_links(a,x,mu);
00183 #else
00184 twist_boundary(a,x);
00185 #endif
00186 return(a);
00187 }
00188 return tmp;
00189 #endif
00190 }
00191
00192 inline mdp_complex& operator() (site x, int mu, int i, int j) {
00193 #ifdef TWISTED_BOUNDARY
00194 if(!in_block(x))
00195 error("call to &operator=(site x) per x out of block");
00196 #endif
00197 return *(address(x,(mu*nc+i)*nc+j));
00198 }
00199
00200 inline const mdp_complex& operator() (site x, int mu, int i, int j) const {
00201 #ifdef TWISTED_BOUNDARY
00202 if(!in_block(x))
00203 error("call to &operator=(site x) per x out of block");
00204 #endif
00205 return *(address(x,(mu*nc+i)*nc+j));
00206 }
00207
00208 inline mdp_matrix operator() (site x, int sign, int mu) {
00209 mdp_matrix tmp;
00210 if(sign==+1) tmp=(*this)(x,mu);
00211 if(sign==-1) tmp=hermitian((*this)(x-mu,mu));
00212 return tmp;
00213 }
00214 inline const mdp_matrix operator() (site x, int sign, int mu) const {
00215 mdp_matrix tmp;
00216 if(sign==+1) tmp=(*this)(x,mu);
00217 if(sign==-1) tmp=hermitian((*this)(x-mu,mu));
00218 return tmp;
00219 }
00220
00221
00222 #ifndef TWISTED_BOUNDARY
00223 inline const mdp_complex operator() (site x, int sign, int mu,
00224 int i, int j) const {
00225 if(sign==+1) return *(address(x,(mu*nc+i)*nc+j));
00226 if(sign==-1) return conj(*(address(x-mu,(mu*nc+j)*nc+i)));
00227 error("call to U(x,0,mu,i,j)");
00228 return mdp_complex(0,0);
00229 };
00230 #endif
00231
00232
00233
00234
00235
00236
00237
00238
00239 #ifdef TWISTED_BOUNDARY
00240 inline friend void twist_boundary(mdp_matrix &M, site &x) {
00241 static int mu, block;
00242 static mdp_complex z=exp(mdp_complex(0,2.0*Pi/3.0));
00243 static mdp_complex a,b,c,d,e,f,g,h,i;
00244
00245 for(mu=1; mu<x.lattice().ndim; mu++) {
00246 block=x.block[mu];
00247 if(block!=0) {
00248 a=M(0,0); b=M(0,1); c=M(0,2);
00249 d=M(1,0); e=M(1,1); f=M(1,2);
00250 g=M(2,0); h=M(2,1); i=M(2,2);
00251 if(mu==1 && block==1) {
00252 M(0,0)=e; M(0,1)=f; M(0,2)=d;
00253 M(1,0)=h; M(1,1)=i; M(1,2)=g;
00254 M(2,0)=b; M(2,1)=c; M(2,2)=a;
00255 }
00256 if(mu==1 && block==-1) {
00257 M(0,0)=i; M(0,1)=g; M(0,2)=h;
00258 M(1,0)=c; M(1,1)=a; M(1,2)=b;
00259 M(2,0)=f; M(2,1)=d; M(2,2)=e;
00260 }
00261 if(mu==2 && block==1) {
00262 M(0,0)=a; M(0,1)=b/z; M(0,2)=c/z/z;
00263 M(1,0)=d*z; M(1,1)=e; M(1,2)=f/z;
00264 M(2,0)=g*z*z; M(2,1)=h*z; M(2,2)=i;
00265 }
00266 if(mu==2 && block==-1) {
00267 M(0,0)=a; M(0,1)=b*z; M(0,2)=c*z*z;
00268 M(1,0)=d/z; M(1,1)=e; M(1,2)=f*z;
00269 M(2,0)=g/z/z; M(2,1)=h/z; M(2,2)=i;
00270 }
00271 if(mu==3 && block==1) {
00272 M(0,0)=i; M(0,1)=g/z; M(0,2)=h/z/z;
00273 M(1,0)=c*z; M(1,1)=a; M(1,2)=b/z;
00274 M(2,0)=f*z*z; M(2,1)=d*z; M(2,2)=e;
00275 }
00276 if(mu==3 && block==-1) {
00277 M(0,0)=e; M(0,1)=f*z; M(0,2)=d*z*z;
00278 M(1,0)=h/z; M(1,1)=i; M(1,2)=g*z;
00279 M(2,0)=b/z/z; M(2,1)=c/z; M(2,2)=a;
00280 }
00281 if(block*block>1) error("two blocks out\n");
00282 }
00283 }
00284 }
00285
00286 friend void define_twist_matrices() {
00287 begin_function("gauge_field__define_twist_matrices");
00288 mdp_complex z=exp(mdp_complex(0, 2.0*Pi/3.0));
00289 OmegaTwist[0]=identity(3);
00290 OmegaTwist[1]=zero(3);
00291 OmegaTwist[1](0,1)=1;
00292 OmegaTwist[1](1,2)=1;
00293 OmegaTwist[1](2,0)=1;
00294 OmegaTwist[2]=zero(3);
00295 OmegaTwist[2](0,0)=1.0/z;
00296 OmegaTwist[2](1,1)=1;
00297 OmegaTwist[2](2,2)=z;
00298 OmegaTwist[3]=zero(3);
00299 OmegaTwist[3](0,2)=1.0/z;
00300 OmegaTwist[3](1,0)=1;
00301 OmegaTwist[3](2,1)=z;
00302 end_function("gauge_field__define_twist_matrices");
00303 }
00304
00305 inline friend void twist_eat_fields(mdp_matrix &M, site &x,
00306 gauge_field &omega) {
00307 begin_function("gauge_field__twist_eat_matrices");
00308 static int mu, block;
00309 static mdp_complex z=exp(mdp_complex(0,2.0*Pi/3.0));
00310 static mdp_complex a,b,c,d,e,f,g,h,i;
00311
00312 for(mu=1; mu<x.lattice().ndim; mu++) {
00313 block=x.block[mu];
00314 if(block!=0) {
00315 if(mu==1 && block==1)
00316 M=omega(x,mu)*M*hermitian(omega(x,mu));
00317 if(mu==1 && block==-1)
00318 M=hermitian(omega(x,mu))*M*omega(x,mu);
00319 if(mu==2 && block==1)
00320 M=omega(x,mu)*M*hermitian(omega(x,mu));
00321 if(mu==2 && block==-1)
00322 M=hermitian(omega(x,mu))*M*omega(x,mu);
00323 if(mu==3 && block==1)
00324 M=omega(x,mu)*M*hermitian(omega(x,mu));
00325 if(mu==3 && block==-1)
00326 M=hermitian(omega(x,mu))*M*omega(x,mu);
00327 if(block*block>1) error("two blocks out\n");
00328 }
00329 }
00330 end_function("gauge_field__twist_eat_matrices");
00331 }
00332 #endif
00333 };
00334
