qpms/symmetries.py

   1 import sympy
   2 from sympy.combinatorics import Permutation, PermutationGroup
   3 sympy.init_printing(perm_cyclic = True)
   4 import cmath
   5 from cmath import exp, pi
   6 from math import sqrt
   7 import numpy as np
   8 np.set_printoptions(linewidth=200)
   9 import numbers
  10 import re
  11 ň = None
  12 from .tmatrices import zflip_tyty, xflip_tyty, yflip_tyty, zrotN_tyty, WignerD_yy_fromvector, identity_tyty, apply_ndmatrix_left
  13 from .cyquaternions import IRot3
  14 from .cycommon import get_mn_y
  15
  16 s3long = np.sqrt(np.longdouble(3.))
  17
  18 def grouprep_try(tdict, src, im, srcgens, imgens, immultop = None, imcmp = None):
  19     tdict[src] = im
  20     for i in range(len(srcgens)):
  21         new_src = src * srcgens[i]
  22         new_im = (im * imgens[i]) if (immultop is None) else immultop(im, imgens[i])
  23         if new_src not in tdict.keys():
  24             grouprep_try(tdict, new_src, new_im, srcgens, imgens, immultop, imcmp)
  25         elif ((new_im != tdict[new_src]) if (imcmp is None) else (not imcmp(new_im, tdict[new_src]))): # check consistency
  26             print(src, ' * ', srcgens[i], ' --> ', new_src)
  27             print(im)
  28             print(' * ')
  29             print(imgens[i])
  30             print(' --> ')
  31             print(new_im)
  32             print(' != ')
  33             print(tdict[new_src])
  34             raise ValueError("Homomorphism inconsistency detected")
  35     return
  36
  37 def group_dps_try(elemlist, elem, gens):
  38     '''Deterministic group depth-first search'''
  39     elemlist.append(elem)
  40     for i in range(len(gens)):
  41         newelem = elem * gens[i]
  42         if newelem not in elemlist:
  43             group_dps_try(elemlist, newelem, gens)
  44     return
  45
  46 class SVWFPointGroupInfo: # only for point groups, coz in svwf_rep() I use I_tyty, not I_ptypty or something alike
  47     def __init__(self,
  48                  name,
  49                  permgroupgens, # permutation group generators
  50                  irrepgens_dict, # dictionary with irrep generators,
  51                  svwf_rep_gen_func, # function that generates a tuple with svwf representation generators
  52                  rep3d_gens = None, # 3d (quaternion) representation generators of a point group: sequence of qpms.irep3 instances
  53                 ):
  54         self.name = name
  55         self.permgroupgens = permgroupgens
  56         self.permgroup = PermutationGroup(*permgroupgens)
  57         self.irrepgens_dict = irrepgens_dict
  58         self.svwf_rep_gen_func = svwf_rep_gen_func
  59         self.irreps = dict()
  60         for irrepname, irrepgens in irrepgens_dict.items():
  61             is1d = isinstance(irrepgens[0], (int,float,complex))
  62             irrepdim = 1 if is1d else irrepgens[0].shape[0]
  63             self.irreps[irrepname] = generate_grouprep(self.permgroup,
  64                                                        1 if is1d else np.eye(irrepdim),
  65                                                        permgroupgens, irrepgens,
  66                                                        immultop = None if is1d else np.dot,
  67                                                        imcmp = None if is1d else np.allclose
  68                                                       )
  69         self.rep3d_gens = rep3d_gens
  70         self.rep3d = None if rep3d_gens is None else generate_grouprep(
  71                 self.permgroup,
  72                 IRot3(),
  73                 permgroupgens, rep3d_gens,
  74                 immultop = None, imcmp = (lambda x, y: x.isclose(y))
  75                 )
  76
  77     def deterministic_elemlist(self):
  78         thelist = list()
  79         group_dps_try(thelist, self.permgroup.identity, self.permgroupgens)
  80         return thelist
  81
  82     def svwf_rep(self, lMax, *rep_gen_func_args, **rep_gen_func_kwargs):
  83         '''
  84         This method generates full SVWF (reducible) representation of the group.
  85         '''
  86         svwfgens = self.svwf_rep_gen_func(lMax, *rep_gen_func_args, **rep_gen_func_kwargs)
  87         my, ny = get_mn_y(lMax)
  88         nelem = len(my)
  89         I_tyty = np.moveaxis(np.eye(2)[:,:,ň,ň] * np.eye(nelem), 2,1)
  90         return generate_grouprep(self.permgroup, I_tyty, self.permgroupgens, svwfgens, immultop = mmult_tyty, imcmp = np.allclose)
  91
  92     def svwf_irrep_projectors(self, lMax, *rep_gen_func_args, **rep_gen_func_kwargs):
  93         return gen_point_group_svwfrep_projectors(self.permgroup, self.irreps, self.svwf_rep(lMax, *rep_gen_func_args, **rep_gen_func_kwargs))
  94
  95     # alternative, for comparison and testing; should give the same results
  96     def svwf_irrep_projectors2(self, lMax, *rep_gen_func_args, **rep_gen_func_kwargs):
  97         return gen_point_group_svwfrep_projectors2(self.permgroup, self.irreps, self.svwf_rep(lMax, *rep_gen_func_args, **rep_gen_func_kwargs))
  98
  99     def svwf_irrep_projectors2_w_bases(self, lMax, *rep_gen_func_args, **rep_gen_func_kwargs):
 100         return gen_point_group_svwfrep_projectors2_w_bases(self.permgroup, self.irreps, self.svwf_rep(lMax, *rep_gen_func_args, **rep_gen_func_kwargs))
 101
 102     def generate_c_source(self):
 103         '''
 104         Generates a string with a chunk of C code with a definition of a qpms_finite_group_t instance.
 105         See also groups.h.
 106         '''
 107         permlist = self.deterministic_elemlist()
 108         order = len(permlist)
 109         permindices = {perm: i for i, perm in enumerate(permlist)} # 'invert' permlist
 110         identity = self.permgroup.identity
 111         s = "{\n"
 112         # char *name
 113         s += '  "%s", // name\n' % self.name
 114         # size_t order;
 115         s += '  %d, // order\n' % order
 116         # qpms_gmi_t idi
 117         s += '  %d, // idi\n' % permindices[identity]
 118         # qpms_gmi_t *mt
 119         s += '  (qpms_gmi_t[]) { // mt\n'
 120         for i in range(order):
 121             ss = ', '.join([str(permindices[permlist[i]*permlist[j]]) for j in range(order)])
 122             s += '    ' + ss + ',\n'
 123         s += '  },\n'
 124         # qpms_gmi_t *invi
 125         s += '  (qpms_gmi_t[]) { // invi\n'
 126         s += '    ' + ', '.join([str(permindices[permlist[j]**-1]) for j in range(order)])
 127         s += '\n  },\n'
 128         # qpms_gmi_t *gens
 129         s += '  (qpms_gmi_t[]) {' + ', '.join([str(permindices[g]) for g in self.permgroupgens]) + '}, // gens\n'
 130         # int ngens
 131         s += '  %d, // ngens\n' % len(self.permgroupgens)
 132         # qpms_permutation_t permrep[]
 133         s += '  (qpms_permutation_t[]){ // permrep\n'
 134         for i in range(order):
 135             s += '    "%s",\n' % str(permlist[i])
 136         s += '  },\n'
 137         # char **elemlabels
 138         s += '  NULL, // elemlabels\n'
 139         # int permrep_nelem
 140         s += '  %d, // permrep_nelem\n' % self.permgroup.degree
 141         # qpms_irot3_t rep3d[]
 142         if self.rep3d is None:
 143             s += '  NULL, // rep3d TODO!!!\n'
 144         else:
 145             s += '  (qpms_irot3_t[]) { // rep3d\n'
 146             for i in range(order):
 147                 s += '   ' + self.rep3d[permlist[i]].crepr() + ',\n'
 148             s += '  },\n'
 149         # int nirreps
 150         s += '  %d, // nirreps\n' % len(self.irreps)
 151         # struct qpms_finite_grep_irrep_t irreps[]
 152         s += '  (struct qpms_finite_group_irrep_t[]) { // irreps\n'
 153         for irname  in sorted(self.irreps.keys()):
 154             irrep = self.irreps[irname]
 155             s += '    {\n'
 156             is1d = isinstance(irrep[identity], (int, float, complex))
 157             dim = 1 if is1d else irrep[identity].shape[0]
 158             # int dim
 159             s += '      %d, // dim\n' % dim
 160             # char name[]
 161             s += '      "%s", //name\n' % re.escape(irname)
 162
 163             # complex double *m
 164             if (is1d):
 165                 s += '      (complex double []) {' + ', '.join([str(irrep[permlist[i]]) for i in range(order)]) + '} // m\n'
 166             else:
 167                 s += '      (complex double []) {\n'
 168                 for i in range(order):
 169                     s += '        // %s\n' % str(permlist[i])
 170                     for row in range(dim):
 171                         s += '        '
 172                         for col in range(dim):
 173                             s += '%s, ' % re.sub('j', '*I', str(irrep[permlist[i]][row,col]))
 174                         s += '\n'
 175                     mat = irrep[permlist[i]]
 176                 s += '      }\n'
 177
 178             #s += '       %d, // dim\n' %
 179             s += '    },\n'
 180         s += '  } // end of irreps\n'
 181         s += '}'
 182         return s
 183
 184 # srcgroup is expected to be PermutationGroup and srcgens of the TODO
 185 # imcmp returns True if two elements of the image group are 'equal', otherwise False
 186 def generate_grouprep(srcgroup, im_identity, srcgens, imgens, immultop = None, imcmp = None):
 187     sz = srcgens[0].size
 188     for g in srcgens:
 189         if g.size != sz:
 190             raise ValueError('All the generators must have the same "size"')
 191     tdict = dict()
 192     grouprep_try(tdict, Permutation(sz-1), im_identity, srcgens, imgens, immultop = immultop, imcmp = imcmp)
 193     if(srcgroup.order() != len(tdict.keys())): # basic check
 194         raise ValueError('The supplied "generators" failed to generate the preimage group: ',
 195                          srcgroup.order(), " != ", len(tdict.keys()))
 196     return tdict
 197
 198 # matrices appearing in 2d representations of common groups as used in Bradley, Cracknell p. 61 (with arabic names instead of greek, because lambda is a keyword)
 199 epsilon = np.eye(2)
 200 alif = np.array(((-1/2,-s3long/2),(s3long/2,-1/2)))
 201 bih = np.array(((-1/2,s3long/2),(-s3long/2,-1/2)))
 202 kaf = np.array(((0,1),(1,0)))
 203 lam = np.array(((1,0),(0,-1)))
 204 ra = np.array(((0,-1),(1,0)))
 205 mim =  np.array(((-1/2,-s3long/2),(-s3long/2,1/2)))
 206 nun =  np.array(((-1/2,s3long/2),(s3long/2,1/2)))
 207
 208
 209
 210
 211 def mmult_tyty(a, b):
 212         return(apply_ndmatrix_left(a, b, (-4,-3)))
 213 def mmult_ptypty(a, b):
 214     return(apply_ndmatrix_left(a, b, (-6,-5,-4)))
 215
 216 def gen_point_group_svwfrep_irreps(permgroup, matrix_irreps_dict, sphrep_full):
 217     '''
 218     Gives the projection operators $P_kl('\Gamma')$ from Dresselhaus (4.28)
 219     for all irreps $\Gamma$ of D3h.;
 220     as an array with indices [k,l,t,y,t,y]
 221
 222     Example of creating last argument:
 223     sphrep_full = generate_grouprep(D3h_permgroup, I_tyty, D3h_srcgens, [C3_tyty, vfl_tyty, zfl_tyty],
 224                            immultop = mmult_tyty, imcmp = np.allclose)
 225     '''
 226     order = permgroup.order()
 227     sphreps = dict()
 228     nelem = sphrep_full[permgroup[0]].shape[-1] # quite ugly hack
 229     for repkey, matrixrep in matrix_irreps_dict.items():
 230         arepmatrix = matrixrep[permgroup[0]] # just one of the matrices to get the shape etc
 231         if isinstance(arepmatrix, numbers.Number):
 232             dim = 1 # repre dimension
 233             preprocess = lambda x: np.array([[x]])
 234         elif isinstance(arepmatrix, np.ndarray):
 235             if(len(arepmatrix.shape)) != 2 or arepmatrix.shape[0] != arepmatrix.shape[1]:
 236                 raise ValueError("Arrays representing irrep matrices must be of square shape")
 237             dim = arepmatrix.shape[0]
 238             preprocess = lambda x: x
 239         else:
 240             raise ValueError("Irrep is not a square array or number")
 241         sphrep = np.zeros((dim,dim,2,nelem,2,nelem), dtype=complex)
 242         for i in permgroup.elements:
 243             sphrep += preprocess(matrixrep[i]).conj().transpose()[:,:,ň,ň,ň,ň] * sphrep_full[i]
 244         sphrep *= dim / order
 245         # clean the nonexact values here
 246         for x in [0, 0.5, -0.5, 0.5j, -0.5j]:
 247             sphrep[np.isclose(sphrep,x)]=x
 248         sphreps[repkey] = sphrep
 249     return sphreps
 250
 251
 252 def gen_point_group_svwfrep_projectors(permgroup, matrix_irreps_dict, sphrep_full):
 253     '''
 254     The same as gen_point_group_svwfrep_irreps, but summed over the kl diagonal, so
 255     one gets single projector onto each irrep space and the arrays have indices
 256     [t, y, t, y]
 257     '''
 258     summedprojs = dict()
 259     for repi, W in gen_point_group_svwfrep_irreps(permgroup, matrix_irreps_dict, sphrep_full).items():
 260         irrepd = W.shape[0]
 261         if irrepd == 1:
 262             mat = np.reshape(W, W.shape[-4:])
 263         else:
 264             mat = np.zeros(W.shape[-4:], dtype=complex) # TODO the result should be real — check!
 265             for d in range(irrepd):
 266                 mat += W[d,d]
 267         if not np.allclose(mat.imag, 0):
 268             raise ValueError("The imaginary part of the resulting projector should be zero, damn!")
 269         else:
 270             summedprojs[repi] = mat.real
 271     return summedprojs
 272
 273
 274 def gen_point_group_svwfrep_projectors2_w_bases(permgroup, matrix_irreps_dict, sphrep_full):
 275     return gen_point_group_svwfrep_projectors2(permgroup, matrix_irreps_dict, sphrep_full, do_bases = True)
 276
 277 def gen_point_group_svwfrep_projectors2(permgroup, matrix_irreps_dict, sphrep_full, do_bases = False):
 278     '''
 279     an approach as in gen_hexlattice_Kpoint_svwf_rep_projectors; for comparison and testing
 280     '''
 281     if (do_bases):
 282         bases = dict()
 283     projectors = dict()
 284     for repi, W in gen_point_group_svwfrep_irreps(permgroup, matrix_irreps_dict, sphrep_full).items():
 285         nelem = W.shape[-1] # however, this should change between iterations
 286         totalvecs = 0
 287         tmplist = list()
 288         for t in (0,1):
 289          for y in range(nelem):
 290            for ai in range(W.shape[0]):
 291             for bi in range(W.shape[1]):
 292                 v = np.zeros((2, nelem))
 293                 v[t,y] = 1
 294                 v1 = np.tensordot(W[ai,bi], v, axes = ([-2,-1],[0,1]))
 295
 296                 if not np.allclose(v1,0):
 297                     v1 = normalize(v1)
 298                     for v2 in tmplist:
 299                         dot = np.tensordot(v1.conjugate(),v2, axes=([-2,-1],[0,1]))
 300                         if not (np.allclose(dot,0)):
 301                             if not np.allclose(np.abs(dot),1):
 302                                 raise ValueError('You have to fix this piece of code.')
 303                             break
 304                     else:
 305                         totalvecs += 1
 306                         tmplist.append(v1)
 307         theprojector = np.zeros((2,nelem, 2, nelem), dtype = float)
 308         if do_bases:
 309             thebasis = np.zeros((len(tmplist), 2, nelem), dtype=complex)
 310             for i, v in enumerate(tmplist):
 311                 thebasis[i] = v
 312             bases[repi] = thebasis
 313         for v in tmplist:
 314             theprojector += (v[:,:,ň,ň] * v.conjugate()[ň,ň,:,:]).real
 315         for x in [0, 1, -1, sqrt(.5), -sqrt(.5), .5, -.5]:
 316             theprojector[np.isclose(theprojector,x)] = x
 317         projectors[repi] = theprojector
 318     if do_bases:
 319         return projectors, bases
 320     else:
 321         return projectors
 322
 323
 324 # Group D3h; mostly legacy code (kept because of the the honeycomb lattice K-point code, whose generalised version not yet implemented)
 325 # Note that the size argument of permutations is necessary, otherwise e.g. c*c and  b*b would not be evaluated equal
 326 # N.B. the weird elements as Permutation(N) – it means identity permutation of size N+1.
 327 rot3_perm = Permutation(0,1,2, size=5) # C3 rotation
 328 xflip_perm = Permutation(0,2, size=5) # vertical mirror
 329 zflip_perm = Permutation(3,4, size=5) # horizontal mirror
 330 D3h_srcgens = [rot3_perm,xflip_perm,zflip_perm]
 331 D3h_permgroup = PermutationGroup(*D3h_srcgens) # D3h
 332
 333 D3h_irreps = {
 334     # Bradley, Cracknell p. 61
 335     "E'" : generate_grouprep(D3h_permgroup, epsilon, D3h_srcgens, [alif, lam, epsilon], immultop = np.dot, imcmp = np.allclose),
 336     "E''" : generate_grouprep(D3h_permgroup, epsilon, D3h_srcgens, [alif, lam, -epsilon], immultop = np.dot, imcmp = np.allclose),
 337     # Bradley, Cracknell p. 59, or Dresselhaus, Table A.14 (p. 482)
 338     "A1'" : generate_grouprep(D3h_permgroup, 1, D3h_srcgens, [1,1,1]),
 339     "A2'" : generate_grouprep(D3h_permgroup, 1, D3h_srcgens, [1,-1,1]),
 340     "A1''" : generate_grouprep(D3h_permgroup, 1, D3h_srcgens, [1,-1,-1]),
 341     "A2''" : generate_grouprep(D3h_permgroup, 1, D3h_srcgens, [1,1,-1]),
 342 }
 343
 344 #TODO lepší název fce; legacy, use group_info['D3h'].generate_grouprep() instead
 345 def gen_point_D3h_svwf_rep(lMax, vflip = 'x'):
 346     '''
 347     Gives the projection operators $P_kl('\Gamma')$ from Dresselhaus (4.28)
 348     for all irreps $\Gamma$ of D3h.;
 349     as an array with indices [k,l,t,y,t,y]
 350     '''
 351
 352     my, ny = get_mn_y(lMax)
 353     nelem = len(my)
 354     C3_yy = WignerD_yy_fromvector(lMax, np.array([0,0,2*pi/3]))
 355     C3_tyty = np.moveaxis(np.eye(2)[:,:,ň,ň] * C3_yy, 2,1)
 356     zfl_tyty = zflip_tyty(lMax)
 357     #yfl_tyty = yflip_tyty(lMax)
 358     #xfl_tyty = xflip_tyty(lMax)
 359     vfl_tyty = yflip_tyty(lMax) if vflip == 'y' else xflip_tyty(lMax)
 360     I_tyty = np.moveaxis(np.eye(2)[:,:,ň,ň] * np.eye(nelem), 2,1)
 361     order = D3h_permgroup.order()
 362     sphrep_full = generate_grouprep(D3h_permgroup, I_tyty, D3h_srcgens, [C3_tyty, vfl_tyty, zfl_tyty],
 363                            immultop = mmult_tyty, imcmp = np.allclose)
 364     sphreps = dict()
 365     for repkey, matrixrep in D3h_irreps.items():
 366         arepmatrix = matrixrep[rot3_perm] # just one of the matrices to get the shape etc
 367         if isinstance(arepmatrix, numbers.Number):
 368             dim = 1 # repre dimension
 369             preprocess = lambda x: np.array([[x]])
 370         elif isinstance(arepmatrix, np.ndarray):
 371             if(len(arepmatrix.shape)) != 2 or arepmatrix.shape[0] != arepmatrix.shape[1]:
 372                 raise ValueError("Arrays representing irrep matrices must be of square shape")
 373             dim = arepmatrix.shape[0]
 374             preprocess = lambda x: x
 375         else:
 376             raise ValueError("Irrep is not a square array or number")
 377         sphrep = np.zeros((dim,dim,2,nelem,2,nelem), dtype=complex)
 378         for i in D3h_permgroup.elements:
 379             sphrep += preprocess(matrixrep[i]).conj().transpose()[:,:,ň,ň,ň,ň] * sphrep_full[i]
 380         sphrep *= dim / order
 381         # clean the nonexact values here
 382         for x in [0, 0.5, -0.5, 0.5j, -0.5j]:
 383             sphrep[np.isclose(sphrep,x)]=x
 384         sphreps[repkey] = sphrep
 385     return sphreps
 386
 387 def gen_hexlattice_Kpoint_svwf_rep(lMax, psi, vflip = 'x'):
 388     my, ny = get_mn_y(lMax)
 389     nelem = len(my)
 390     C3_yy = WignerD_yy_fromvector(lMax, np.array([0,0,2*pi/3]))
 391     C3_tyty = np.moveaxis(np.eye(2)[:,:,ň,ň] * C3_yy, 2,1)
 392     zfl_tyty = zflip_tyty(lMax)
 393     #yfl_tyty = yflip_tyty(lMax)
 394     #xfl_tyty = xflip_tyty(lMax)
 395     vfl_tyty = yflip_tyty(lMax) if vflip == 'y' else xflip_tyty(lMax)
 396     I_tyty = np.moveaxis(np.eye(2)[:,:,ň,ň] * np.eye(nelem), 2,1)
 397     hex_C3_K_ptypty = np.diag([exp(-psi*1j*2*pi/3),exp(+psi*1j*2*pi/3)])[:,ň,ň,:,ň,ň] * C3_tyty[ň,:,:,ň,:,:]
 398     hex_zfl_ptypty = np.eye(2)[:,ň,ň,:,ň,ň] * zfl_tyty[ň,:,:,ň,:,:]
 399     #hex_xfl_ptypty = np.array([[0,1],[1,0]])[:,ň,ň,:,ň,ň] * xfl_tyty[ň,:,:,ň,:,:]
 400     hex_vfl_ptypty = np.array([[0,1],[1,0]])[:,ň,ň,:,ň,ň] * vfl_tyty[ň,:,:,ň,:,:]
 401     hex_I_ptypty = np.eye((2*2*nelem)).reshape((2,2,nelem,2,2,nelem))
 402     order = D3h_permgroup.order()
 403     hex_K_sphrep_full = generate_grouprep(D3h_permgroup, hex_I_ptypty, D3h_srcgens, [hex_C3_K_ptypty, hex_vfl_ptypty, hex_zfl_ptypty],
 404                            immultop = mmult_ptypty, imcmp = np.allclose)
 405     hex_K_sphreps = dict()
 406     for repkey, matrixrep in D3h_irreps.items():
 407         arepmatrix = matrixrep[rot3_perm] # just one of the matrices to get the shape etc
 408         if isinstance(arepmatrix, numbers.Number):
 409             dim = 1 # repre dimension
 410             preprocess = lambda x: np.array([[x]])
 411         elif isinstance(arepmatrix, np.ndarray):
 412             if(len(arepmatrix.shape)) != 2 or arepmatrix.shape[0] != arepmatrix.shape[1]:
 413                 raise ValueError("Arrays representing irrep matrices must be of square shape")
 414             dim = arepmatrix.shape[0]
 415             preprocess = lambda x: x
 416         else:
 417             raise ValueError("Irrep is not a square array or number")
 418         sphrep = np.zeros((dim,dim,2,2,nelem,2,2,nelem), dtype=complex)
 419         for i in D3h_permgroup.elements:
 420             sphrep += preprocess(matrixrep[i]).conj().transpose()[:,:,ň,ň,ň,ň,ň,ň] * hex_K_sphrep_full[i]
 421         sphrep *= dim / order
 422         # clean the nonexact values here
 423         for x in [0, 0.5, -0.5, 0.5j, -0.5j]:
 424             sphrep[np.isclose(sphrep,x)]=x
 425         hex_K_sphreps[repkey] = sphrep
 426     return hex_K_sphreps
 427
 428 def normalize(v):
 429     norm = np.linalg.norm(v.reshape((np.prod(v.shape),)), ord=2)
 430     if norm == 0:
 431        return v*np.nan
 432     return v / norm
 433
 434 def gen_hexlattice_Kpoint_svwf_rep_projectors(lMax, psi, vflip='x', do_bases=False):
 435     nelem = lMax * (lMax+2)
 436     projectors = dict()
 437     if do_bases:
 438         bases = dict()
 439     for repi, W in gen_hexlattice_Kpoint_svwf_rep(lMax,psi,vflip=vflip).items():
 440         totalvecs = 0
 441         tmplist = list()
 442         for p in (0,1):
 443          for t in (0,1):
 444           for y in range(nelem):
 445            for ai in range(W.shape[0]):
 446             for bi in range(W.shape[1]):
 447                 v = np.zeros((2,2,nelem))
 448                 v[p,t,y] = 1
 449                 #v = np.ones((2,2,nelem))
 450                 v1 = np.tensordot(W[ai,bi],v, axes = ([-3,-2,-1],[0,1,2]))
 451
 452
 453                 if not np.allclose(v1,0):
 454                     v1 = normalize(v1)
 455                     for v2 in tmplist:
 456                         dot = np.tensordot(v1.conjugate(),v2,axes = ([-3,-2,-1],[0,1,2]))
 457                         if not np.allclose(dot,0):
 458                             if not np.allclose(np.abs(dot),1):
 459                                 raise ValueError('You have to fix this piece of code.')# TODO maybe I should make sure that the absolute value is around 1
 460                             break
 461                     else:
 462                         totalvecs += 1
 463                         tmplist.append(v1)
 464                         #for index, x in np.ndenumerate(v1):
 465                         #    if x!=0:
 466                         #        print(index, x)
 467                         #print('----------')
 468         theprojector = np.zeros((2,2,nelem,2,2,nelem), dtype = float)
 469         if do_bases:
 470             thebasis = np.zeros((len(tmplist), 2,2,nelem), dtype=complex)
 471             for i, v in enumerate(tmplist):
 472                 thebasis[i] = v
 473             bases[repi] = thebasis
 474         for v in tmplist:
 475             theprojector += (v[:,:,:,ň,ň,ň] * v.conjugate()[ň,ň,ň,:,:,:]).real # TODO check is it possible to have imaginary elements?
 476         for x in [0, 1, -1,sqrt(0.5),-sqrt(0.5),0.5,-0.5]:
 477             theprojector[np.isclose(theprojector,x)]=x
 478         projectors[repi] = theprojector
 479     if do_bases:
 480         return projectors, bases
 481     else:
 482         return projectors
 483
 484
 485
 486 point_group_info = { # representation info of some useful point groups
 487     # TODO real trivial without generators
 488     'trivial_g' : SVWFPointGroupInfo('trivial_g',
 489                  # permutation group generators
 490                                (   # I put here the at least the identity for now (it is reduntant, but some functions are not robust enough to have an empty set of generators
 491                                    Permutation(),
 492                                ),
 493                  # dictionary with irrep generators
 494                                {
 495                                     "A" : (1,),
 496                                },
 497                  # function that generates a tuple with svwf representation generators
 498                                lambda lMax : (identity_tyty(lMax),),
 499                  # quaternion rep generators
 500                                rep3d_gens = (
 501                                    IRot3.identity(),
 502                                 )
 503     ),
 504     'C2' : SVWFPointGroupInfo('C2',
 505                  # permutation group generators
 506                                (Permutation(0,1), # 180 deg rotation around z axis
 507                                ),
 508                  # dictionary with irrep generators
 509                                {
 510                                     # Bradley, Cracknell p. 57;
 511                                     'A': (1,),
 512                                     'B': (-1,),
 513                                },
 514                  # function that generates a tuple with svwf representation generators
 515                                lambda lMax : (zrotN_tyty(2, lMax),),
 516                  # quaternion rep generators
 517                                 rep3d_gens = (
 518                                     IRot3.zrotN(2),
 519                                 )
 520
 521     ),
 522     'C2v' : SVWFPointGroupInfo('C2v',
 523                  # permutation group generators
 524                                (Permutation(0,1, size=4)(2,3), # x -> - x mirror operation (i.e. yz mirror plane)
 525                                 Permutation(0,3, size=4)(1,2), # y -> - y mirror operation (i.e. xz mirror plane)
 526                                ),
 527                  # dictionary with irrep generators
 528                                {
 529                                     # Bradley, Cracknell p. 58; not sure about the labels / axes here
 530                                     'A1': (1,1),
 531                                     'B2': (-1,1),
 532                                     'A2': (-1,-1),
 533                                     'B1': (1,-1),
 534                                },
 535                  # function that generates a tuple with svwf representation generators
 536                                lambda lMax : (xflip_tyty(lMax), yflip_tyty(lMax)),
 537                  # quaternion rep generators
 538                                 rep3d_gens = (
 539                                     IRot3.xflip(),
 540                                     IRot3.yflip(),
 541                                 )
 542
 543     ),
 544     'D2h' : SVWFPointGroupInfo('D2h',
 545                  # permutation group generators
 546                                (Permutation(0,1, size=6)(2,3), # x -> - x mirror operation (i.e. yz mirror plane)
 547                                 Permutation(0,3, size=6)(1,2), # y -> - y mirror operation (i.e. xz mirror plane)
 548                                 # ^^^ btw, I guess that Permutation(0,1, size=6) and Permutation(2,3, size=6) would
 549                                 # do exactly the same job (they should; CHECK)
 550                                 Permutation(4,5, size=6)       # z -> - z mirror operation (i.e. xy mirror plane)
 551                                ),
 552                  # dictionary with irrep generators
 553                                {
 554                                     # Product of C2v and zflip; not sure about the labels / axes here
 555                                     "A1'": (1,1,1),
 556                                     "B2'": (-1,1,1),
 557                                     "A2'": (-1,-1,1),
 558                                     "B1'": (1,-1,1),
 559                                     "A1''": (-1,-1,-1),
 560                                     "B2''": (1,-1,-1),
 561                                     "A2''": (1,1,-1),
 562                                     "B1''": (-1,1,-1),
 563                                },
 564                  # function that generates a tuple with svwf representation generators
 565                                lambda lMax : (xflip_tyty(lMax), yflip_tyty(lMax), zflip_tyty(lMax)),
 566                  # quaternion rep generators
 567                                rep3d_gens = (
 568                                    IRot3.xflip(),
 569                                    IRot3.yflip(),
 570                                    IRot3.zflip(),
 571                                 )
 572     ),
 573     'C4' : SVWFPointGroupInfo('C4',
 574                  # permutation group generators
 575                                (Permutation(0,1,2,3, size=4),), #C4 rotation
 576                  # dictionary with irrep generators
 577                                {
 578                                     # Bradley, Cracknell p. 58
 579                                     'A': (1,),
 580                                     'B': (-1,),
 581                                     '1E': (-1j,),
 582                                     '2E': (1j,),
 583                                },
 584                  # function that generates a tuple with svwf representation generators
 585                                lambda lMax : (zrotN_tyty(4, lMax), ),
 586                  # quaternion rep generators
 587                                rep3d_gens = (
 588                                    IRot3.zrotN(4),
 589                                 )
 590     ),
 591     'C4v' : SVWFPointGroupInfo('C4v',
 592                  # permutation group generators
 593                                (Permutation(0,1,2,3, size=4), #C4 rotation
 594                                 Permutation(0,1, size=4)(2,3)), # x -> - x mirror operation (i.e. yz mirror plane)
 595                  # dictionary with irrep generators
 596                                {
 597                                     # Bradley, Cracknell p. 62
 598                                     'E': (ra, -lam),
 599                                     # Bradley, Cracknell p. 59, or Dresselhaus, Table A.18
 600                                     'A1': (1,1),
 601                                     'A2': (1,-1),
 602                                     'B1': (-1,1),
 603                                     'B2': (-1,-1),
 604                                },
 605                  # function that generates a tuple with svwf representation generators
 606                                lambda lMax : (zrotN_tyty(4, lMax), xflip_tyty(lMax)),
 607                  # quaternion rep generators
 608                                rep3d_gens = (
 609                                    IRot3.zrotN(4),
 610                                    IRot3.xflip(),
 611                                 )
 612     ),
 613     'D4h' : SVWFPointGroupInfo('D4h',
 614                  # permutation group generators
 615                                (Permutation(0,1,2,3, size=6),  # C4 rotation
 616                                 Permutation(0,1, size=6)(2,3), # x -> - x mirror operation (i.e. yz mirror plane)
 617                                 Permutation(4,5, size=6), # horizontal mirror operation z -> -z (i.e. xy mirror plane)
 618                                ),
 619                  # dictionary with irrep generators
 620                                {    # product of C4v and zflip
 621                                     "E'": (ra, -lam, epsilon),
 622                                     "E''":(ra, -lam, -epsilon),
 623                                     "A1'": (1,1,1),
 624                                     "A2'": (1,-1,1),
 625                                     "A1''": (1,-1,-1),
 626                                     "A2''": (1,1,-1),
 627                                     "B1'": (-1,1,1),
 628                                     "B2'": (-1,-1,1),
 629                                     "B1''": (-1,-1,-1),
 630                                     "B2''": (-1,1,-1),
 631                                },
 632                  # function that generates a tuple with svwf representation generators
 633                                lambda lMax : (zrotN_tyty(4, lMax), xflip_tyty(lMax), zflip_tyty(lMax)),
 634                  # quaternion rep generators
 635                                rep3d_gens = (
 636                                    IRot3.zrotN(4),
 637                                    IRot3.xflip(),
 638                                    IRot3.zflip(),
 639                                 )
 640     ),
 641     'D3h' : SVWFPointGroupInfo('D3h',
 642                  # permutation group generators
 643                                ( Permutation(0,1,2, size=5),  # C3 rotation
 644                                  Permutation(0,2, size=5), # vertical mirror
 645                                  Permutation(3,4, size=5), # horizontal mirror z -> -z (i.e. xy mirror plane)
 646                                ),
 647                  # dictionary with irrep generators
 648                                {    # Bradley, Cracknell p. 61
 649                                     "E'" : (alif, lam, epsilon),
 650                                     "E''" : (alif, lam, -epsilon),
 651                                     # Bradley, Cracknell p. 59, or Dresselhaus, Table A.14 (p. 482)
 652                                     "A1'" : (1,1,1),
 653                                     "A2'" : (1,-1,1),
 654                                     "A1''" : (1,-1,-1),
 655                                     "A2''" : (1,1,-1),
 656                                },
 657                  # function that generates a tuple with svwf representation generators
 658                                lambda lMax, vflip: (zrotN_tyty(3, lMax), yflip_tyty(lMax) if vflip == 'y' else xflip_tyty(lMax), zflip_tyty(lMax)),
 659                  # quaternion rep generators
 660                                rep3d_gens = (
 661                                    IRot3.zrotN(3),
 662                                    IRot3.xflip(), # if vflip == 'y' else IRot3.xflip(), # FIXME enable to choose
 663                                    IRot3.zflip(),
 664                                 )
 665     ),
 666     'x_and_z_flip': SVWFPointGroupInfo(
 667         'x_and_z_flip',
 668         (
 669             Permutation(0,1, size=4), # x -> -x mirror op
 670             Permutation(2,3, size=4), # z -> -z mirror op
 671         ),
 672         {
 673             "P'": (1, 1),
 674             "R'": (-1, 1),
 675             "P''": (-1,-1),
 676             "R''": (1, -1),
 677         },
 678         lambda lMax : (xflip_tyty(lMax), zflip_tyty(lMax)),
 679         rep3d_gens = (
 680             IRot3.xflip(),
 681             IRot3.zflip(),
 682         )
 683
 684     ),
 685     'y_and_z_flip': SVWFPointGroupInfo(
 686         'y_and_z_flip',
 687         (
 688             Permutation(0,1, size=4), # y -> -y mirror op
 689             Permutation(2,3, size=4), # z -> -z mirror op
 690         ),
 691         {
 692             "P'": (1, 1),
 693             "R'": (-1, 1),
 694             "P''": (-1,-1),
 695             "R''": (1, -1),
 696         },
 697         lambda lMax : (yflip_tyty(lMax), zflip_tyty(lMax)),
 698         rep3d_gens = (
 699             IRot3.yflip(),
 700             IRot3.zflip(),
 701         )
 702     ),
 703 }