From 2f4a0c5de4fc64ca5006e62f809e8efc1dbf4121 Mon Sep 17 00:00:00 2001 From: Azusa Yamaguchi Date: Wed, 10 Jun 2015 11:24:14 +0100 Subject: [PATCH] Successful generation of general SU N generators This class will move into a utils class and gain exponentiation and ta projection to support HMC and heatbath --- tests/Test_lie_generators.cc | 229 +++++++++++++++++++++++++++++++++++ 1 file changed, 229 insertions(+) create mode 100644 tests/Test_lie_generators.cc diff --git a/tests/Test_lie_generators.cc b/tests/Test_lie_generators.cc new file mode 100644 index 00000000..b29682a0 --- /dev/null +++ b/tests/Test_lie_generators.cc @@ -0,0 +1,229 @@ +#include + +#include +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +class suN { +public: + + static int generators(int ncolour) { return ncolour*ncolour-1; } + static int su2subgroups(int ncolour) { return (ncolour*(ncolour-1))/2; } + + template using suNmatrix = iScalar > > ; + + //////////////////////////////////////////////////////////////////////// + // There are N^2-1 generators for SU(N). + // + // We take a traceless hermitian generator basis as follows + // + // * Normalisation: trace ta tb = 1/2 delta_ab + // + // * Off diagonal + // - pairs of rows i1,i2 behaving like pauli matrices signma_x, sigma_y + // + // - there are (Nc-1-i1) slots for i2 on each row [ x 0 x ] + // direct count off each row + // + // - Sum of all pairs is Nc(Nc-1)/2: proof arithmetic series + // + // (Nc-1) + (Nc-2)+... 1 ==> Nc*(Nc-1)/2 + // 1+ 2+ + + Nc-1 + // + // - There are 2 x Nc (Nc-1)/ 2 of these = Nc^2 - Nc + // + // - We enumerate the row-col pairs. + // - for each row col pair there is a (sigma_x) and a (sigma_y) like generator + // + // + // t^a_ij = { in 0.. Nc(Nc-1)/2 -1} => delta_{i,i1} delta_{j,i2} + delta_{i,i1} delta_{j,i2} + // t^a_ij = { in Nc(Nc-1)/2 ... Nc^(Nc-1) -1} => i delta_{i,i1} delta_{j,i2} - i delta_{i,i1} delta_{j,i2} + // + // * Diagonal; must be traceless and normalised + // - Sequence is + // N (1,-1,0,0...) + // N (1, 1,-2,0...) + // N (1, 1, 1,-3,0...) + // N (1, 1, 1, 1,-4,0...) + // + // where 1/2 = N^2 (1+.. m^2)etc.... for the m-th diagonal generator + // NB this gives the famous SU3 result for su2 index 8 + // + // N= sqrt(1/2 . 1/6 ) = 1/2 . 1/sqrt(3) + // + // ( 1 ) + // ( 1 ) / sqrt(3) /2 = 1/2 lambda_8 + // ( -2) + //////////////////////////////////////////////////////////////////////// + template + static void suNgenerator(int lieIndex,suNmatrix &ta){ + // map lie index to which type of generator + int diagIndex; + int su2Index; + int sigxy; + int NNm1 = Ncolour*(Ncolour-1); + if ( lieIndex>= NNm1 ) { + diagIndex = lieIndex -NNm1; + suNgeneratorDiagonal(diagIndex,ta); + return; + } + sigxy = lieIndex&0x1; + su2Index= lieIndex>>1; + if ( sigxy ) suNgeneratorSigmaY(su2Index,ta); + else suNgeneratorSigmaX(su2Index,ta); + } + template + static void suNgeneratorSigmaX(int su2Index,suNmatrix &ta){ + ta=zero; + int i1,i2; + su2SubGroupIndex(i1,i2,su2Index); + ta()()(i1,i2)=1.0; + ta()()(i2,i1)=1.0; + ta= ta *0.5; + } + template + static void suNgeneratorSigmaY(int su2Index,suNmatrix &ta){ + ta=zero; + Complex i(0.0,1.0); + int i1,i2; + su2SubGroupIndex(i1,i2,su2Index); + ta()()(i1,i2)=-i; + ta()()(i2,i1)= i; + ta= ta *0.5; + } + template + static void suNgeneratorDiagonal(int diagIndex,suNmatrix &ta){ + ta=zero; + int trsq=0; + int last=diagIndex+1; + for(int i=0;i<=diagIndex;i++){ + ta()()(i,i) = 1.0; + trsq++; + } + ta()()(last,last) = -last; + trsq+=last*last; + RealD nrm = 1.0/std::sqrt(2.0*trsq); + ta = ta *nrm; + } + //////////////////////////////////////////////////////////////////////// + // Map a + // + //////////////////////////////////////////////////////////////////////// + template + static void su2SubGroupIndex(int &i1,int &i2,int su2_index){ + + assert( (su2_index>=0) && (su2_index< (Ncolour*(Ncolour-1))/2) ); + + int spare=su2_index; + for(i1=0;spare >= (Ncolour-1-i1);i1++ ){ + spare = spare - (Ncolour-1-i1); // remove the Nc-1-i1 terms + } + i2=i1+1+spare; + } + template + static void su2Extract(std::vector &r,const Lattice > &source, int su2_index) + { + assert(r.size() == 4); // store in 4 real parts + + for(int i=0;i<4;i++){ + conformable(r[i],source); + } + + int i1,i2; + su2SubGroupIndex(i1,i2,su2_index); + + /* Compute the b(k) of A_SU(2) = b0 + i sum_k bk sigma_k */ + r[0] = real(source()()(i1,i1) + source()()(i2,i2)); + r[1] = imag(source()()(i1,i2) + source()()(i2,i1)); + r[2] = real(source()()(i1,i2) - source()()(i2,i1)); + r[3] = imag(source()()(i1,i1) - source()()(i2,i2)); + } + + + template static void printGenerators(void) + { + for(int gen=0;gen ta; + suN::suNgenerator(gen,ta); + std::cout<< "Nc = "< static void testGenerators(void){ + suNmatrix ta; + suNmatrix tb; + std::cout<<"Checking trace ta tb is 0.5 delta_ab"< simd_layout = GridDefaultSimd(4,vComplexF::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + std::vector latt_size ({4,4,4,4}); + + GridCartesian Fine(latt_size,simd_layout,mpi_layout); + + LatticeGaugeField Umu(&Fine); + + std::cout<<"*********************************************"<(); + suN::testGenerators<2>(); + std::cout<<"*********************************************"<(); + suN::testGenerators<3>(); + std::cout<<"*********************************************"<(); + suN::testGenerators<4>(); + std::cout<<"*********************************************"<(); + suN::testGenerators<5>(); + + Grid_finalize(); +} + +