From 7bf6d1dc7beee49eb552c4d483f32fc082ce4c12 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Mon, 17 Dec 2018 16:46:51 +0200
Subject: [PATCH] Symmetry projection operator generation

Former-commit-id: 58b0d7b3f2a292c26964571edf159b26f0eef0ed
---
 qpms/symmetries.py | 74 +++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 73 insertions(+), 1 deletion(-)

diff --git a/qpms/symmetries.py b/qpms/symmetries.py
index fcd241f..dc7dc7d 100644
--- a/qpms/symmetries.py
+++ b/qpms/symmetries.py
@@ -63,6 +63,10 @@ class SVWFPointGroupInfo: # only for point groups, coz in svwf_rep() I use I_tyt
     
     def svwf_irrep_projectors(self, lMax, *rep_gen_func_args, **rep_gen_func_kwargs):
         return gen_point_group_svwfrep_projectors(self.permgroup, self.irreps, self.svwf_rep(lMax, *rep_gen_func_args, **rep_gen_func_kwargs))
+    
+    # alternative, for comparison and testing; should give the same results
+    def svwf_irrep_projectors2(self, lMax, *rep_gen_func_args, **rep_gen_func_kwargs):
+        return gen_point_group_svwfrep_projectors(self.permgroup, self.irreps, self.svwf_rep(lMax, *rep_gen_func_args, **rep_gen_func_kwargs))
 
 # srcgroup is expected to be PermutationGroup and srcgens of the TODO
 # imcmp returns True if two elements of the image group are 'equal', otherwise False
@@ -96,7 +100,7 @@ def mmult_tyty(a, b):
 def mmult_ptypty(a, b):
     return(qpms.apply_ndmatrix_left(a, b, (-6,-5,-4)))
 
-def gen_point_group_svwfrep_projectors(permgroup, matrix_irreps_dict, sphrep_full):
+def gen_point_group_svwfrep_irreps(permgroup, matrix_irreps_dict, sphrep_full):
     '''
     Gives the projection operators $P_kl('\Gamma')$ from Dresselhaus (4.28)
     for all irreps $\Gamma$ of D3h.;
@@ -131,6 +135,74 @@ def gen_point_group_svwfrep_projectors(permgroup, matrix_irreps_dict, sphrep_ful
         sphreps[repkey] = sphrep
     return sphreps
 
+
+def gen_point_group_svwfrep_projectors(permgroup, matrix_irreps_dict, sphrep_full):
+    '''
+    The same as gen_point_group_svwfrep_irreps, but summed over the kl diagonal, so
+    one gets single projector onto each irrep space and the arrays have indices
+    [t, y, t, y]
+    '''
+    summedprojs = dict()
+    for repi, W in gen_point_group_svwfrep_irreps(permgroup, matrix_irreps_dict, sphrep_full).items():
+        irrepd = W.shape[0]
+        if irrepd == 1:
+            mat = np.reshape(W, W.shape[-4:])
+        else:
+            mat = np.zeros(W.shape[-4:], dtype=complex) # TODO the result should be real — check!
+            for d in range(irrepd):
+                mat += W[d,d]
+        if not np.allclose(mat.imag, 0):
+            raise ValueError("The imaginary part of the resulting projector should be zero, damn!")
+        else:
+            summedprojs[repi] = mat.real
+    return summedprojs
+
+def gen_point_group_svwfrep_projectors2(permgroup, matrix_irreps_dict, sphrep_full, do_bases = False):
+    '''
+    an approach as in gen_hexlattice_Kpoint_svwf_rep_projectors; for comparison and testing
+    '''
+    nelem = sphrep_full.values[0].shape[-1]
+    if (do_bases):
+        bases = dict()
+    projectors = dict()
+    for repi, W in gen_point_group_svwfrep_irreps(permgroup, matrix_irreps_dict, sphrep_full).items():
+        totalvecs = 0
+        tmplist = list()
+        for t in (0,1):
+         for y in range(nelem):
+           for ai in range(W.shape[0]):
+            for bi in range(W.shape[1]):
+                v = np.zeros((2, nelem))
+                v[t,y] = 1
+                v1 = np.tensordot(W[ai,bi], v, axes = ([-2,-1],[0,1]))
+
+                if not np.allclose(v1,0):
+                    v1 = normalize(v1)
+                    for v2 in tmplist:
+                        dot = np.tensordot(v1.conjugate(),v2, axes=([-2,-1],[0,1]))
+                        if not (np.allclose(dot,0)):
+                            if not np.allclose(np.abs(dot),1):
+                                raise ValueError('You have to fix this piece of code.')
+                            break
+                    else:
+                        totalvecs += 1
+                        tmplist.append(v1)
+        theprojectors = np.zeros((2,nelem, 2, nelem), dtype = float)
+        if do_bases:
+            thebasis = np.zeros((len(tmplist), 2, nelem), dtype=complex)
+            for i, v in enumerate(tmplist):
+                thebasis[i] = v
+            beses[repi] = thebasis
+        for v in tmplist:
+            theprojector += (v[:,:,ň,ň] * v.conjugate()[ň,ň,ň,:,:,:]).real 
+        for x in [0, 1, -1, sqrt(.5), -sqrt(.5), .5, -.5]:
+            theprojector[np.isclose(theprojector,x)] = x
+    if do_bases:
+        return projectors, bases
+    else:
+        return projectors
+
+
 # Group D3h; mostly legacy code (kept because of the the honeycomb lattice K-point code, whose generalised version not yet implemented)
 # Note that the size argument of permutations is necessary, otherwise e.g. c*c and  b*b would not be evaluated equal
 # N.B. the weird elements as Permutation(N) – it means identity permutation of size N+1.