diff --git a/python/triqs_dft_tools/converters/elk.py b/python/triqs_dft_tools/converters/elk.py
index 50e6ce82..0b62a9b3 100644
--- a/python/triqs_dft_tools/converters/elk.py
+++ b/python/triqs_dft_tools/converters/elk.py
@@ -142,7 +142,7 @@ class ElkConverter(ConverterTools,Elk_tools,read_Elk):
              for isym in range(n_symm):
            #rotate symmetry matrix into basis defined by T
                mat[isym][ish]=numpy.matmul(T[ish],mat[isym][ish])
-               mat[isym][ish]=numpy.matmul(mat[isym][ish],T[ish].transpose())
+               mat[isym][ish]=numpy.matmul(mat[isym][ish],T[ish].conjugate().transpose())
            #put desired subset of transformed symmetry matrix into temp matrix for symmetry isym
                for id in range(len(ind[ish])):
                  i=ind[ish][id]
@@ -348,15 +348,19 @@ class ElkConverter(ConverterTools,Elk_tools,read_Elk):
         use_rotations = 1
         rot_mat = [numpy.identity(corr_shells[icrsh]['dim'], numpy.complex_) for icrsh in range(n_corr_shells)]
         for icrsh in range(n_corr_shells):
-          incrsh = corr_to_inequiv[icrsh]
-          iatom = corr_shells[incrsh]['atom']
+          #incrsh = corr_to_inequiv[icrsh]
+          #iatom = corr_shells[incrsh]['atom']
+          #want to rotate atom to first inequivalent atom in list
+          iatom = 1
           for isym in range(n_symm):
             jatom=perm[isym][corr_shells[icrsh]['atom']-1]
             #determinant determines if crystal symmetry matrix has inversion symmetry (=-1)
             det = numpy.linalg.det(symmat[isym][:,:])
             if((jatom==iatom)&(det>0.0)):
-              rot_mat[icrsh][:,:]=mat[isym][icrsh][:,:]
-            #used first desired symmetry in crystal symmetry list
+              #local rotation which rotates equivalent atom into its local coordinate system
+              #(inverse of the symmetry operator applied to the projectors in Elk)
+              rot_mat[icrsh][:,:]=mat[isym][icrsh][:,:].conjugate().transpose()
+              #used first desired symmetry in crystal symmetry list
               break
         # Elk does not currently use time inversion symmetry
         rot_mat_time_inv = [0 for i in range(n_corr_shells)]
diff --git a/python/triqs_dft_tools/converters/elktools/elk_converter_tools.py b/python/triqs_dft_tools/converters/elktools/elk_converter_tools.py
index c7686e3d..66f91a65 100644
--- a/python/triqs_dft_tools/converters/elktools/elk_converter_tools.py
+++ b/python/triqs_dft_tools/converters/elktools/elk_converter_tools.py
@@ -92,8 +92,9 @@ class ElkConverterTools:
             v[2]=numpy.sqrt(abs(rotp[2,2]+1.0)/2.0)
             v[0]=(rotp[0,2]+rotp[2,0])/(4.0*v[2])
             v[1]=(rotp[1,2]+rotp[2,1])/(4.0*v[2])
-      # return theta and v
-      return v,th
+      # return -theta and v. -theta is returned as TRIQS does not rotate
+      # the observable (such as the density matrix) which is done in Elk
+      return v,-th
 
     def axangsu2(self,v,th):
       """
@@ -198,9 +199,10 @@ class ElkConverterTools:
             #ang=r.as_euler('zyz')
             ang=self.zyz_euler(rot)
             #Elk uses inverse rotations, i.e. the function is being rotated, not the spherical harmonics
-            angi[0]=-ang[2]
-            angi[1]=-ang[1]
-            angi[2]=-ang[0]
+            #TRIQS rotates the spherical harmonics instead
+            angi[0]=ang[0]
+            angi[1]=ang[1]
+            angi[2]=ang[2]
             #calculate the symmetry in the complex spherical harmonic basis.
             d = self.ylmrot(p,angi,l)
             symmat[isym][ish][:,:] = d[:,:]
diff --git a/test/python/elk_convert.ref.h5 b/test/python/elk_convert.ref.h5
index 2ffe186b..8a100766 100644
Binary files a/test/python/elk_convert.ref.h5 and b/test/python/elk_convert.ref.h5 differ