qpms/symmetries.py

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