Repository: irreducible-representations/irrep (language: python)

Issue:
Wannier
Tested Wannier90 interface. Main change: Wannier90 admits 3 different forms of setting the value of a parameter in .win file, namely with symbols =, : or " " (blank space). Only the first case was implemented, now all of them are.


Discussion / hints:

Fix patch:
diff --git a/irrep/__version.py b/irrep/__version.py
index 74faa3e39..453af2b0c 100644
--- a/irrep/__version.py
+++ b/irrep/__version.py
@@ -1,1 +1,1 @@
-__version__="1.5.7"
+__version__="1.6.0"
\ No newline at end of file
diff --git a/irrep/bandstructure.py b/irrep/bandstructure.py
index 4d7826ae7..395a29a14 100644
--- a/irrep/bandstructure.py
+++ b/irrep/bandstructure.py
@@ -67,6 +67,12 @@ class BandStructure:
     spin_channel : str, default=None
         Selection of the spin-channel. 'up' for spin-up, 'dw' for spin-down. 
         Only applied in the interface to Quantum Espresso.
+    refUC : array, default=None
+        3x3 array describing the transformation of vectors defining the 
+        unit cell to the standard setting.
+    shiftUC : array
+        Translation taking the origin of the unit cell used in the DFT 
+        calculation to that of the standard setting.
 
     Attributes
     ----------
@@ -112,7 +118,9 @@ def __init__(
         code="vasp",
         EF='0.0',
         onlysym=False,
-        spin_channel=None
+        spin_channel=None,
+        refUC = None,
+        shiftUC = None
     ):
         code = code.lower()
 
@@ -122,19 +130,19 @@ def __init__(
 
         if code == "vasp":
             self.__init_vasp(
-                fWAV, fPOS, Ecut, IBstart, IBend, kplist, spinor, EF=EF, onlysym=onlysym
+                fWAV, fPOS, Ecut, IBstart, IBend, kplist, spinor, EF=EF, onlysym=onlysym, refUC=refUC, shiftUC=shiftUC
             )
         elif code == "abinit":
             self.__init_abinit(
-                fWFK, Ecut, IBstart, IBend, kplist, EF=EF, onlysym=onlysym
+                fWFK, Ecut, IBstart, IBend, kplist, EF=EF, onlysym=onlysym, refUC=refUC, shiftUC=shiftUC
             )
         elif code == "espresso":
             self.__init_espresso(
-                prefix, Ecut, IBstart, IBend, kplist, EF=EF, onlysym=onlysym, spin_channel=spin_channel
+                prefix, Ecut, IBstart, IBend, kplist, EF=EF, onlysym=onlysym, spin_channel=spin_channel, refUC=refUC, shiftUC=shiftUC
             )
         elif code == "wannier90":
             self.__init_wannier(
-                prefix, Ecut, IBstart, IBend, kplist, EF=EF, onlysym=onlysym
+                prefix, Ecut, IBstart, IBend, kplist, EF=EF, onlysym=onlysym, refUC=refUC, shiftUC=shiftUC
             )
         else:
             raise RuntimeError("Unknown/unsupported code :{}".format(code))
@@ -150,6 +158,8 @@ def __init_vasp(
         spinor=None,
         EF='0.0',
         onlysym=False,
+        refUC=None,
+        shiftUC=None,
     ):
         """
         Initialization for vasp. Read data and save it in attributes.
@@ -174,12 +184,18 @@ def __init_vasp(
             Fermi-energy.
         onlysym : bool, default=False
             Exit after printing info about space-group.
+        refUC : array, default=None
+            3x3 array describing the transformation of vectors defining the 
+            unit cell to the standard setting.
+        shiftUC : array, default=None
+            Translation taking the origin of the unit cell used in the DFT 
+            calculation to that of the standard setting.
         """
         if spinor is None:
             raise RuntimeError(
                 "spinor should be specified in the command line for VASP bandstructure"
             )
-        self.spacegroup = SpaceGroup(inPOSCAR=fPOS, spinor=spinor)
+        self.spacegroup = SpaceGroup(inPOSCAR=fPOS, spinor=spinor, refUC=refUC, shiftUC=shiftUC)
         self.spinor = spinor
         if onlysym:
             return
@@ -262,7 +278,7 @@ def __init_vasp(
                 Ecut,
                 Ecut0,
                 self.RecLattice,
-                SG=self.spacegroup,
+                symmetries_SG=self.spacegroup.symmetries,
                 spinor=self.spinor,
                 WCF=WCF,
             )
@@ -278,6 +294,8 @@ def __init_abinit(
         kplist=None,
         EF='0.0',
         onlysym=False,
+        refUC=None,
+        shiftUC=None,
     ):
         """
         Initialization for abinit. Read data and store it in attributes.
@@ -298,13 +316,19 @@ def __init_abinit(
             Fermi-energy.
         onlysym : bool, default=False
             Exit after printing info about space-group.
+        refUC : array, default=None
+            3x3 array describing the transformation of vectors defining the 
+            unit cell to the standard setting.
+        shiftUC : array, default=None
+            Translation taking the origin of the unit cell used in the DFT 
+            calculation to that of the standard setting.
         """
 
         header = AbinitHeader(WFKname)
         usepaw = header.usepaw
         self.spinor = header.spinor
         self.spacegroup = SpaceGroup(
-            cell=(header.rprimd, header.xred, header.typat), spinor=self.spinor
+            cell=(header.rprimd, header.xred, header.typat), spinor=self.spinor, refUC=refUC, shiftUC=shiftUC
         )
         if onlysym:
             return
@@ -376,7 +400,7 @@ def __init_abinit(
                 Ecut,
                 Ecut0,
                 self.RecLattice,
-                SG=self.spacegroup,
+                symmetries_SG=self.spacegroup.symmetries,
                 spinor=self.spinor,
                 code="abinit",
                 kpt=header.kpt[ik],
@@ -397,6 +421,8 @@ def __init_wannier(
         kplist=None,
         EF='0.0',
         onlysym=False,
+        refUC=None,
+        shiftUC=None,
     ):
         """
         Initialization for wannier90. Read data and store it in attibutes.
@@ -418,6 +444,12 @@ def __init_wannier(
             Fermi-energy.
         onlysym : bool, default=False
             Exit after printing info about space-group.
+        refUC : array, default=None
+            3x3 array describing the transformation of vectors defining the 
+            unit cell to the standard setting.
+        shiftUC : array, default=None
+            Translation taking the origin of the unit cell used in the DFT 
+            calculation to that of the standard setting.
         """
         if Ecut is None:
             raise RuntimeError("Ecut mandatory for Wannier90")
@@ -642,7 +674,7 @@ def check_end(name):
         )
 
         self.spacegroup = SpaceGroup(
-            cell=(self.Lattice, xred, typat), spinor=self.spinor
+            cell=(self.Lattice, xred, typat), spinor=self.spinor, refUC=refUC, shiftUC=shiftUC
         )
         if onlysym:
             return
@@ -684,7 +716,7 @@ def check_end(name):
                 Ecut,
                 None,
                 self.RecLattice,
-                SG=self.spacegroup,
+                symmetries_SG=self.spacegroup.symmetries,
                 spinor=self.spinor,
                 code="wannier",
                 eigenval=eigenval[ik - 1],
@@ -702,7 +734,9 @@ def __init_espresso(
         kplist=None,
         EF='0.0',
         onlysym=False,
-        spin_channel=None
+        spin_channel=None,
+        refUC=None,
+        shiftUC=None,
     ):
         """
         Initialization for Quantum Espresso. Read data and store in attributes.
@@ -726,6 +760,12 @@ def __init_espresso(
             Exit after printing info about space-group.
         spin_channel : str, default=None
             Selection of the spin-channel. 'up' for spin-up, 'dw' for spin-down.
+        refUC : array, default=None
+            3x3 array describing the transformation of vectors defining the 
+            unit cell to the standard setting.
+        shiftUC : array, default=None
+            Translation taking the origin of the unit cell used in the DFT 
+            calculation to that of the standard setting.
         """
         import xml.etree.ElementTree as ET
 
@@ -760,7 +800,7 @@ def __init_espresso(
         xred = (np.array(xcart, dtype=float) * BOHR).dot(np.linalg.inv(self.Lattice))
         #        print ("xred=",xred)
         self.spacegroup = SpaceGroup(
-            cell=(self.Lattice, xred, typat), spinor=self.spinor
+            cell=(self.Lattice, xred, typat), spinor=self.spinor, refUC=refUC, shiftUC=shiftUC
         )
         if onlysym:
             return
@@ -857,7 +897,7 @@ def __init_espresso(
                 Ecut,
                 Ecut0,
                 self.RecLattice,
-                SG=self.spacegroup,
+                symmetries_SG=self.spacegroup.symmetries,
                 spinor=self.spinor,
                 code="espresso",
                 kptxml=kpall[ik],
@@ -884,8 +924,6 @@ def getNK():
     def write_characters(
         self,
         degen_thresh=1e-8,
-        refUC=None,
-        shiftUC=np.zeros(3),
         kpnames=None,
         symmetries=None,
         preline="",
@@ -929,11 +967,9 @@ def write_characters(
         GAP = np.Inf
         Low = -np.Inf
         Up = np.inf
-        if kpnames is not None and refUC is not None:
+        if kpnames is not None:
             for kpname, KP in zip(kpnames, self.kpoints):
-                irreps = self.spacegroup.get_irreps_from_table(
-                    refUC, shiftUC, kpname, KP.K
-                )
+                irreps = self.spacegroup.get_irreps_from_table(kpname, KP.K)
                 ninv, low, up = KP.write_characters(
                     degen_thresh,
                     irreptable=irreps,
@@ -942,6 +978,9 @@ def write_characters(
                     efermi=self.efermi,
                     plotFile=pFile,
                     kpl=kpline,
+                    symmetries_tables=self.spacegroup.symmetries_tables,
+                    refUC=self.spacegroup.refUC,
+                    shiftUC=self.spacegroup.shiftUC
                 )
                 NBANDINV += ninv
                 GAP = min(GAP, up - low)
@@ -1004,8 +1043,6 @@ def getNbands(self):
     def write_trace(
         self,
         degen_thresh=1e-8,
-        refUC=None,
-        shiftUC=np.zeros(3),
         kpnames=None,
         symmetries=None,
     ):
@@ -1018,12 +1055,6 @@ def write_trace(
         degen_thresh : float, default=1e-8
             Threshold energy used to decide whether wave-functions are
             degenerate in energy.
-        refUC : array, default=None
-            3x3 array describing the transformation of vectors defining the 
-            unit cell to the standard setting.
-        shiftUC : array, default=np.zeros(3)
-            Translation taking the origin of the unit cell used in the DFT 
-            calculation to that of the standard setting.
         kpnames : list, default=None
             Labels of maximal k-points at which irreps of bands must be computed. 
             If it is not specified, only traces will be printed, not irreps.
@@ -1038,7 +1069,9 @@ def write_trace(
             ).format(self.getNbands(), 1 if self.spinor else 0)
         )
 
-        f.write(self.spacegroup.write_trace(refUC=refUC, shiftUC=shiftUC))
+        f.write(
+                self.spacegroup.write_trace()
+                )
         # Number of maximal k-vectors in the space group. In the next files
         # introduce the components of the maximal k-vectors))
         f.write("  {0}  \n".format(len(self.kpoints)))
@@ -1046,7 +1079,7 @@ def write_trace(
             f.write(
                 "   ".join(
                     "{0:10.6f}".format(x)
-                    for x in (refUC.dot(KP.K) if refUC is not None else KP.K)
+                    for x in self.spacegroup.refUC.dot(KP.K)
                 )
                 + "\n"
             )
@@ -1262,8 +1295,6 @@ def write_bands(self, locs=None):
     def write_trace_all(
         self,
         degen_thresh=1e-8,
-        refUC=None,
-        shiftUC=np.zeros(3),
         symmetries=None,
         fname="trace_all.dat",
     ):
@@ -1276,12 +1307,6 @@ def write_trace_all(
         degen_thresh : float, default=1e-8
             Threshold energy used to decide whether wave-functions are
             degenerate in energy.
-        refUC : array, default=None
-            3x3 array describing the transformation of vectors defining the 
-            unit cell to the standard setting.
-        shiftUC : array, default=np.zeros(3)
-            Translation taking the origin of the unit cell used in the DFT 
-            calculation to that of the standard setting.
         symmetries : list, default=None
             Index of symmetry operations whose traces will be printed. 
         fname : str, default=trace_all.dat
@@ -1299,9 +1324,7 @@ def write_trace_all(
         f.write(
             "\n".join(
                 ("#" + l)
-                for l in self.spacegroup.write_trace(
-                    refUC=refUC, shiftUC=shiftUC
-                ).split("\n")
+                for l in self.spacegroup.write_trace().split("\n")
             )
             + "\n\n"
         )
diff --git a/irrep/cli.py b/irrep/cli.py
index 31e3b95b6..065cf6045 100644
--- a/irrep/cli.py
+++ b/irrep/cli.py
@@ -136,6 +136,7 @@
 @click.option(
     "-refUC",
     type=str,
+    default=None,
     help="The lattice vectors of the reference unit cell (as given in the crystallographic tables) "
     "expressed in terms of the unit cell vectors used in the calculation. "
     "Nine comma-separated numbers.",
@@ -143,6 +144,7 @@
 @click.option(
     "-shiftUC",
     type=str,
+    default=None,
     help="The vector to shift the calculated unit cell origin (in units of the calculated lattice), "
     "to get the unit cell as defined in crystallographic tables. Three comma-separated numbers.",
 )
@@ -237,12 +239,6 @@ def cli(
     # for k, v in locals().items():
     #     print("{}\t{}".format(k, v))
 
-    # if supplied, convert refUC and shiftUC from comma-separated lists into arrays
-    if refuc:
-        refuc = np.array(refuc.split(","), dtype=float).reshape((3, 3))
-    if shiftuc:
-        shiftuc = np.array(shiftuc.split(","), dtype=float).reshape(3)
-
     # parse input arguments into lists if supplied
     if symmetries:
         symmetries = str2list(symmetries)
@@ -253,9 +249,6 @@ def cli(
     if kpnames:
         kpnames = kpnames.split(",")
 
-    if onlysym:
-        spinor = False
-
     try:
         print(fwfk.split("/")[0].split("-"))
         preline = " ".join(s.split("_")[1] for s in fwfk.split("/")[0].split("-")[:3])
@@ -263,9 +256,6 @@ def cli(
         print(err)
         preline = ""
 
-    if (refuc is not None) and (shiftuc is None):
-        shiftuc = np.zeros(3)
-
     bandstr = BandStructure(
         fWAV=fwav,
         fWFK=fwfk,
@@ -279,19 +269,18 @@ def cli(
         code=code,
         EF=ef,
         onlysym=onlysym,
+        refUC = refuc,
+        shiftUC = shiftuc
     )
-    bandstr.spacegroup.show(refUC=refuc, shiftUC=shiftuc, symmetries=symmetries)
+    bandstr.spacegroup.show(symmetries=symmetries)
 
     if onlysym:
         exit()
 
-    if refuc is None:
-        refuc = np.eye(3)
-    if shiftuc is None:
-        shiftuc = np.zeros(3)
-
     with open("irreptable-template", "w") as f:
-        f.write(bandstr.spacegroup.str(refUC=refuc, shiftUC=shiftuc))
+        f.write(
+                bandstr.spacegroup.str()
+                )
 
     subbands = {(): bandstr}
 
@@ -310,7 +299,7 @@ def cli(
         for k in subbands:
             print("symmetry eigenvalue : {0} \n Traces are : ".format(k))
             subbands[k].write_characters(
-                degen_thresh=0.001, refUC=refuc, symmetries=symmetries
+                degen_thresh=0.001, symmetries=symmetries
             )
             print("symmetry eigenvalue : {0} \n Zak phases are : ".format(k))
             zak = subbands[k].zakphase()
@@ -336,8 +325,6 @@ def cli(
 
     bandstr.write_trace(
         degen_thresh=degenthresh,
-        refUC=refuc,
-        shiftUC=shiftuc,
         symmetries=symmetries,
     )
     for k, sub in subbands.items():
@@ -351,8 +338,6 @@ def cli(
         plotfile=None # being implemented, not finished yet...
         sub.write_characters(
             degen_thresh=degenthresh,
-            refUC=refuc,
-            shiftUC=shiftuc,
             symmetries=symmetries,
             kpnames=kpnames,
             preline=preline,
diff --git a/irrep/kpoint.py b/irrep/kpoint.py
index c382ba282..717d08b5e 100644
--- a/irrep/kpoint.py
+++ b/irrep/kpoint.py
@@ -51,8 +51,10 @@ class Kpoint:
     RecLattice : array, shape=(3,3)
         Each row contains the cartesian coordinates of a basis vector forming 
         the unit-cell in reciprocal space.
-    SG : class, default=None
-        Instance of `SpaceGroup`.
+    symmetries_SG : list, default=None
+        Each element is an instance of class `SymmetryOperation` corresponding 
+        to a symmetry in the point group of the space-group. The value 
+        passed is the attribute `symmetries` of class `SpaceGroup`.
     spinor : bool, default=None
         `True` if wave functions are spinors, `False` if they are scalars.
     code : str, default='vasp'
@@ -109,8 +111,6 @@ class Kpoint:
         to short plane-waves based on energy (ascending order). Fitfth 
         (sixth) row contains the index of the first (last) plane-wave with 
         the same energy as the plane-wave of the current column.
-    SG : class
-        Instance of `SpaceGroup`.
     K : array, shape=(3,)
         Direct coordinates of the k-point.
     Energy : array
@@ -119,7 +119,10 @@ class Kpoint:
         Energy of the first band above the set of bands whose traces should be 
         calculated. It will be set to `numpy.NaN` if the last band matches the 
         last band in the DFT calculation (default `IBend`).
-
+    symmetries_SG : list
+        Each element is an instance of class `SymmetryOperation` corresponding 
+        to a symmetry in the point group of the space-group. The value 
+        passed is the attribute `symmetries` of class `SpaceGroup`.
     """
     
     # creates attribute symmetries, if it was not created before
@@ -141,9 +144,9 @@ def symmetries(self):
         For more about `lazy-property`, check the documentation `here <https://pypi.org/project/lazy-property/>`_ .
         """
         symmetries = {}
-        #        print ("calculating symmetry eigenvalues for E={0}, WF={1} SG={2}".format(self.Energy,self.WF.shape,self.SG) )
-        if not (self.SG is None):
-            for symop in self.SG.symmetries:
+        #        print ("calculating symmetry eigenvalues for E={0}, WF={1} SG={2}".format(self.Energy,self.WF.shape,symmetries_SG) )
+        if not (self.symmetries_SG is None):
+            for symop in self.symmetries_SG:
                 try:
                     symmetries[symop] = symm_eigenvalues(
                         self.K,
@@ -168,7 +171,7 @@ def __init__(
         Ecut,
         Ecut0,
         RecLattice,
-        SG=None,
+        symmetries_SG=None,
         spinor=None,
         code="vasp",
         kpt=None,
@@ -188,6 +191,7 @@ def __init__(
         self.Nband = IBend - IBstart
         #        self.n=np.arange(IBstart,IBend)+1
         self.RecLattice = RecLattice
+        self.symmetries_SG = symmetries_SG  # lazy_property needs it
 
         if code.lower() == "vasp":
             self.WF, self.ig = self.__init_vasp(
@@ -222,10 +226,6 @@ def __init__(
         self.WF /= (
             np.sqrt(np.abs(np.einsum("ij,ij->i", self.WF.conj(), self.WF)))
         ).reshape(self.Nband, 1)
-        self.SG = SG
-
-    #        self.__calc_sym_eigenvalues()
-    #        print("WF=\n",WF)
 
     def copy_sub(self, E, WF):
         """
@@ -944,6 +944,9 @@ def write_characters(
         efermi=0.0,
         plotFile=None,
         kpl="",
+        symmetries_tables=None,
+        refUC=np.eye(3),
+        shiftUC=np.zeros(3)
     ):
         """
         Calculate traces and determine and print irreps in a k-point. Write them 
@@ -959,7 +962,7 @@ def write_characters(
             Returned by method `get_irreps_from_table` of class `SpaceGroup`. 
             Each key is the label of an irrep, each value another `dict`. Keys 
             of every secondary `dict` are indices of symmetries (starting from 
-            1 and following order of operations as returned by `spglib`) and 
+            1 and following order of operations in tables of BCS) and 
             values are traces of symmetries.
         symmetries : list, default=None
             Index of symmetry operations whose traces will be printed. 
@@ -973,6 +976,16 @@ def write_characters(
         kpl : float, default=''
             Length accumulated until the k-point. Can be used to locate irreps 
             in the x-axis of a band structure plot.
+        symmetries_tables : list, default=None
+            Each component is an instance of class `SymopTable` corresponding to a 
+            symmetry operation in the "point-group" of the space-group. Values 
+            passed are attributes `symmetries` of class `IrrepTable`.
+        refUC : array, default=np.eye(3)
+            3x3 array describing the transformation of vectors defining the 
+            unit cell to the standard setting.
+        shiftUC : array, default=np.zeros(3)
+            Translation taking the origin of the unit cell used in the DFT 
+            calculation to that of the standard setting.
 
         Returns
         -------
@@ -989,7 +1002,9 @@ def write_characters(
             sym = {s.ind: s for s in self.symmetries}
         else:
             sym = {s.ind: s for s in self.symmetries if s.ind in symmetries}
-
+        
+        # Generate array char, where each row corresponds to a sym. op
+        # and every column to a wave function
         char = np.vstack([self.symmetries[sym[i]] for i in sorted(sym)])
         borders = np.hstack(
             [
@@ -1004,7 +1019,7 @@ def write_characters(
         Nirrep = int(round(Nirrep))
         char = np.array(
             [char[:, start:end].sum(axis=1) for start, end in zip(borders, borders[1:])]
-        )
+        )  # Every column in char corresponds to a sym. op.
 
         #        print(" char ",char.shape,"\n",char)
         writeimaginary = np.abs(char.imag).max() > 1e-4
@@ -1037,7 +1052,7 @@ def write_characters(
                 print("irreptable:", irreptable)
                 print([sym.ind for sym in self.symmetries])
                 raise ke
-            # generate str describing irrep corresponding to sets of states
+            # Generate str describing irrep corresponding to sets of states
             irreps = [
                 ", ".join(
                     ir
@@ -1048,39 +1063,75 @@ def write_characters(
                         else ""
                     )
                     + ")"
-                    for ir in irr # irreps of little-group
-                    if abs(irr[ir]) > 1e-3 # check multiplicity
+                    for ir in irr  # Irreps of little-group
+                    if abs(irr[ir]) > 1e-3  # Check multiplicity
                 )
-                for irr in irreps # group of degen. states
+                for irr in irreps  # Group of degen. states
             ]
         #            irreps=[ "None" ]*(len(borders)-1)
-        irreplen = max(len(irr) for irr in irreps) # len of largest line
+
+        # Transfer traces in calculational cell to refUC
+        char_refUC = char.copy()
+        if (not np.allclose(refUC, np.eye(3, dtype=float)) or
+            not np.allclose(shiftUC, np.zeros(3, dtype=float))):
+            # Calculational and reference cells are not identical
+            for i,ind in enumerate(sym):
+                dt = (symmetries_tables[ind-1].t 
+                      - sym[ind].translation_refUC(refUC, shiftUC))
+                k = np.round(refUC.dot(self.K), 5)
+                char_refUC[:,i] *= (sym[ind].sign 
+                                     * np.exp(-2j*np.pi*dt.dot(k)))
+        if np.allclose(char, char_refUC, rtol=0.0, atol=1e-4):
+            write_refUC = False  # Tr identical in both unit cells
+        else:
+            write_refUC = True   # Write traces in refUC
+            print(("For each irrep, traces of symmetries in the calculation "
+                   "unit cell will be printed first and traces of symmetries "
+                   "in tables will be printed in the line below")
+                 )
+
+        irreplen = max(len(irr) for irr in irreps)  # len of largest line
         if irreplen % 2 == 1:
             irreplen += 1
         s2 = " " * int(irreplen / 2 - 3)
+
+        # Header of the block
         print(
             "\n\nk-point {0:3d} :{1} \n number of irreps = {2}".format(
                 self.ik0, self.K, Nirrep
             )
         )
         print("   Energy  | multiplicity |{0} irreps {0}| sym. operations  ".format(s2))
+
+        # Symmetry operations
         print(
             "           |              |{0}        {0}| ".format(s2),
             " ".join(s1 + "{0:4d}    ".format(i) + s1 for i in sorted(sym)),
         )
-        print(
-            "\n".join(
-                (" {0:8.4f}  |    {1:5d}     | {2:" + str(irreplen) + "s} |").format(
-                    e - efermi, d, ir
-                )
-                + " ".join(
-                    "{0:8.4f}".format(c.real)
-                    + ("{0:+7.4f}j".format(c.imag) if writeimaginary else "")
-                    for c in ch
-                )
-                for e, d, ir, ch in zip(E, dim, irreps, char)
-            )
-        )
+
+        # Energy-levels, irrep's label and traces
+        for e, d, ir, ch, ch2 in zip(E, dim, irreps, char, char_refUC):
+            # Print characters in calculational unit cell
+            left_str = (" {0:8.4f}  |    {1:5d}     | {2:{3}s} |"
+                        .format(e - efermi, d, ir, irreplen)
+                       )
+            right_str = " ".join(
+                        "{0:8.4f}".format(c.real)
+                        + ("{0:+7.4f}j".format(c.imag) if writeimaginary else "")
+                        for c in ch
+                    )
+            print(left_str + " " + right_str)
+            # Print characters in reference unit cell
+            if write_refUC:
+                left_str = ("           |              | {0:{1}s} |"
+                            .format(len(ir)*" ", irreplen)
+                           )
+                right_str = " ".join(
+                            "{0:8.4f}".format(c.real)
+                            + ("{0:+7.4f}j".format(c.imag) if writeimaginary else "")
+                            for c in ch2
+                        )
+                print(left_str + " " + right_str)
 
         if plotFile is not None:
             plotFile.write(
@@ -1253,7 +1304,7 @@ def write_trace_all(self, degen_thresh=1e-8, symmetries=None, efermi=0.0, kpline
                         if i in char
                         else (" " * 7 + "X" * 3 + " " * 8 + "X" * 3)
                     )
-                    for i in range(1, len(self.SG.symmetries) + 1)
+                    for i in range(1, len(symmetries_SG) + 1)
                 )
                 for e, d, ib, j in zip(E, dim, IB, np.arange(len(dim)))
             )
diff --git a/irrep/spacegroup.py b/irrep/spacegroup.py
index de11d4012..cb4468190 100644
--- a/irrep/spacegroup.py
+++ b/irrep/spacegroup.py
@@ -77,6 +77,9 @@ class SymmetryOperation():
         `True` if wave-functions are spinors, `False` if they are scalars.
     spinor_rotation : array, shape=(2,2)
         Matrix describing how spinors transform under the symmetry.
+    sign : float
+        Factor needed to match the matrix for the rotation of spinors 
+        to that in tables.
     """
 
     def __init__(self, rot, trans, Lattice, ind=-1, spinor=True):
@@ -94,6 +97,7 @@ def __init__(self, rot, trans, Lattice, ind=-1, spinor=True):
         self.spinor = spinor
         self.spinor_rotation = expm(-0.5j * self.angle *
                                     np.einsum('i,ijk->jk', self.axis, pauli_sigma))
+        self.sign = 1  # May be changed later externally
 
     def get_angle_str(self):
         """
@@ -218,20 +222,28 @@ def translation_refUC(self, refUC, shiftUC):
             self.rotation.dot(shiftUC)).dot(
             np.linalg.inv(refUC))
 
-    def show(self, refUC=None, shiftUC=np.zeros(3)):
+    def show(self, refUC=np.eye(3), shiftUC=np.zeros(3)):
         """
         Print description of symmetry operation.
         
         Parameters
         ----------
-        refUC : array, default=None
+        refUC : array, default=np.eye(3)
             3x3 array describing the transformation of vectors defining the 
             unit cell to the standard setting.
         shiftUC : array, default=np.zeros(3)
             Translation taking the origin of the unit cell used in the DFT 
             calculation to that of the standard setting.
         """
-        print(" # ", self.ind)
+        if (not np.allclose(refUC, np.eye(3)) or
+            not np.allclose(shiftUC, np.zeros(3))):
+            write_ref = True  # To avoid writing this huge block again
+        else:
+            write_ref = False
+
+        # Print header
+        print("\n # ", self.ind)
+        # Print rotation part
         rotstr = [s +
                   " ".join("{0:3d}".format(x) for x in row) +
                   t for s, row, t in zip(["rotation : |", " " *
@@ -239,60 +251,75 @@ def show(self, refUC=None, shiftUC=np.zeros(3)):
                                           "|", " " *
                                           11 +
                                           "|"], self.rotation, [" |", " |", " |"])]
-        if (refUC is not None) or ((shiftUC is not None)
-                                   and np.linalg.norm(shiftUC) > 1e-5):
-            if refUC is None:
-                refUC = np.eye(3)
+        if write_ref:
             fstr = ("{0:3d}")
             R = self.rotation_refUC(refUC)
             rotstr1 = [" " *
                        5 +
                        s +
                        " ".join(fstr.format(x) for x in row) +
-                       t for s, row, t in zip(["in refUC : |", " " *
-                                               11 +
-                                               "|", " " *
-                                               11 +
-                                               "|"], R, [" |", " |", " |"])]
+                       t for s, row, t in zip(["rotation : |",
+                                               " (refUC)   |",
+                                               " " * 11 + "|"
+                                               ],
+                                              R,
+                                              [" |", " |", " |"])]
             rotstr = [r + r1 for r, r1 in zip(rotstr, rotstr1)]
         print("\n".join(rotstr))
 
+        # Print spinor transformation matrix
         if self.spinor:
-            print("\n".join(s +
-                            " ".join("{0:6.3f}{1:+6.3f}j".format(x.real, x.imag) for x in row) +
-                            t for s, row, t in zip(["spinor   : |", " " *
-                                                    11 +
-                                                    "|", " " *
-                                                    11 +
-                                                    "|"], self.spinor_rotation, [" |", " |", " |"])))
-        print(
-            " translation : [ " +
-            " ".join(
-                "{0:10.6f}".format(
-                    x %
-                    1) for x in self.translation.round(6)) +
-            " ] ")
-        if refUC is not None:
-            print("     in the reference unit cell :")
-            print(
-                "     translation : [ " +
-                " ".join(
-                    "{0:10.6f}".format(
-                        x %
-                        1) for x in self.translation_refUC(
-                        refUC,
-                        shiftUC).round(6)) +
-                " ] ")
+            spinstr = [s +
+                       " ".join("{0:6.3f}{1:+6.3f}j".format(x.real, x.imag) for x in row) +
+                       t 
+                       for s, row, t in zip(["spinor rot.  : |",
+                                             " " * 15 + "|",
+                                             ], 
+                                             self.spinor_rotation, 
+                                             [" |", " |"]
+                                           )
+                      ]
+            spinstr1 = [s +
+                       " ".join("{0:6.3f}{1:+6.3f}j".format(x.real, x.imag) for x in row) +
+                       t 
+                       for s, row, t in zip(["spinor rot.  : |",
+                                             "  (refUC)      |"
+                                             ], 
+                                             self.spinor_rotation*self.sign, 
+                                             [" |", " |"]
+                                           )
+                       ]
+
+            print("\n".join([row + "    " + row1 for row,row1 in zip(spinstr, spinstr1)]))
+
+        # Print translation part
+        trastr = ("translation :  [ " 
+                  + " ".join("{0:8.4f}"
+                             .format(x%1) for x in self.translation.round(6)
+                             ) 
+                  + " ] "
+                  )
+        if write_ref:
+            trastr += ("    translation : [ "
+                       + " ".join(
+                          "{0:8.4f}".format(x % 1) 
+                          for x in self.translation_refUC(refUC,shiftUC).round(6)) 
+                       + " ] \n"
+                       + " " * 52 + "(refUC)"  # Lower line
+                       )
+        print(trastr)
+
         print("axis: {0} ; angle = {1}, inversion : {2} ".format(
             self.axis.round(6), self.angle_str, self.inversion))
+        print('sign=',self.sign)  # test
 
-    def str(self, refUC=None, shiftUC=np.zeros(3)):
+    def str(self, refUC=np.eye(3), shiftUC=np.zeros(3)):
         """
         Construct description of symmetry operation.
 
         Parameters
         ----------
-        refUC : array, default=None
+        refUC : array, default=np.eye(3)
             3x3 array describing the transformation of vectors defining the 
             unit cell to the standard setting.
         shiftUC : array, default=np.zeros(3)
@@ -304,9 +331,6 @@ def str(self, refUC=None, shiftUC=np.zeros(3)):
         str
             Description to print.
         """
-        if refUC is None:
-            refUC = np.eye(3, dtype=int)
-#       refUC - row-vectors, expressing the reference unit cell vectors in terms of the lattice used in calculation
 #        print ( "symmetry # ",self.ind )
         R = self.rotation_refUC(refUC)
         t = self.translation_refUC(refUC, shiftUC)
@@ -315,7 +339,7 @@ def str(self, refUC=None, shiftUC=np.zeros(3)):
         return ("   ".join(" ".join(str(x) for x in r) for r in R) + "     " + " ".join(str_(x) for x in t) + ("      " + \
                 "    ".join("  ".join(str_(x) for x in X) for X in (np.abs(S.reshape(-1)), np.angle(S.reshape(-1)) / np.pi))))
 
-    def str2(self, refUC=None, shiftUC=np.zeros(3)):
+    def str2(self, refUC=np.eye(3), shiftUC=np.zeros(3)):
         """
         Print matrix of a symmetry operation in the format: 
         {{R|t}}-> R11,R12,...,R23,R33,t1,t2,t3 and, when SOC was included, the 
@@ -325,7 +349,7 @@ def str2(self, refUC=None, shiftUC=np.zeros(3)):
         
         Parameters
         ----------
-        refUC : array, default=None
+        refUC : array, default=np.eye(3)
             3x3 array describing the transformation of vectors defining the 
             unit cell to the standard setting.
         shiftUC : array, default=np.zeros(3)
@@ -372,6 +396,12 @@ class SpaceGroup():
         `Spglib <https://spglib.github.io/spglib/python-spglib.html#get-symmetry>`_.
     spinor : bool, default=True
         `True` if wave-functions are spinors (SOC), `False` if they are scalars.
+    refUC : array, default=None
+        3x3 array describing the transformation of vectors defining the 
+        unit cell to the standard setting.
+    shiftUC : array, default=None
+        Translation taking the origin of the unit cell used in the DFT 
+        calculation to that of the standard setting.
 
     Attributes
     ----------
@@ -380,6 +410,10 @@ class SpaceGroup():
     symmetries : list
         Each element is an instance of class `SymmetryOperation` corresponding 
         to a symmetry in the point group of the space-group.
+    symmetries_tables : list
+        Attribute `symmetries` of class `IrrepTable`. Each component is an 
+        instance of class `SymopTable` corresponding to a symmetry operation
+        in the "point-group" of the space-group.
     name : str 
         Symbol of the space-group in Hermann-Mauguin notation. 
     number : int 
@@ -490,11 +524,18 @@ def _findsym(self, inPOSCAR, cell):
         array
             3x3 array where cartesian coordinates of basis  vectors **a**, **b** 
             and **c** are given in rows. 
+        array
+            3x3 array describing the transformation of vectors defining the 
+            unit cell to the convenctional setting.
+        array
+            Translation taking the origin of the unit cell used in the DFT 
+            calculation to that of the standard setting of spglib. It may not be
+            the shift taking to the convenctional cell of tables; indeed, in 
+            centrosymmetric groups they adopt origin choice 1 of ITA, rather 
+            than choice 2 (BCS).
         """
         if cell is None:
             cell = self.__cell_vasp(inPOSCAR=inPOSCAR)
-#    cell1=tuple( [cell.lattice,cell.positions,cell.numbers] )
-#    print (cell)
         print('')
         print('\n ----------INFORMATION ABOUT THE UNIT CELL----------- \n')
         print('')
@@ -505,43 +546,71 @@ def _findsym(self, inPOSCAR, cell):
             cell[1],
             '\n Atom type indices: \n',
             cell[2])
-        symmetries = spglib.get_symmetry(cell)
-#    print ("symmetriesreturned by spglib : ",symmetries)
+        dataset = spglib.get_symmetry_dataset(cell)
         symmetries = [
             SymmetryOperation(
                 rot,
-                symmetries['translations'][i],
+                dataset['translations'][i],
                 cell[0],
                 ind=i + 1,
                 spinor=self.spinor) for i,
             rot in enumerate(
-                symmetries['rotations'])]
+                dataset['rotations'])]
         nsym = len(symmetries)
-        s = spglib.get_spacegroup(cell).split(" ")
 
-        return symmetries, s[0], int(s[1].strip("()")), cell[0]
+        return (symmetries, 
+                dataset['international'],
+                dataset['number'], 
+                cell[0], 
+                np.linalg.inv(dataset['transformation_matrix']),
+                dataset['origin_shift']
+                )
 
-    def __init__(self, inPOSCAR=None, cell=None, spinor=True):
+    def __init__(self, inPOSCAR=None, cell=None, spinor=True, refUC=None, shiftUC=None):
         self.spinor = spinor
-        self.symmetries, self.name, self.number, self.Lattice = self._findsym(
-            inPOSCAR, cell)
+        (self.symmetries, 
+         self.name, 
+         self.number, 
+         self.Lattice, 
+         refUC_tmp, 
+         shiftUC_tmp) = self._findsym(inPOSCAR, cell)
         self.RecLattice = np.array([np.cross(self.Lattice[(i + 1) %
                                                           3], self.Lattice[(i + 2) %
                                                                            3]) for i in range(3)]) * 2 * np.pi / np.linalg.det(self.Lattice)
         print("\n Reciprocal lattice:\n", self.RecLattice)
 
-    def show(self, refUC=None, shiftUC=np.zeros(3), symmetries=None):
+        # Determine refUC and shiftUC according to entries in CLI
+        self.symmetries_tables = IrrepTable(self.number, self.spinor).symmetries
+        self.refUC, self.shiftUC = self.determine_basis_transf(
+                                            refUC_cli=refUC, 
+                                            shiftUC_cli=shiftUC,
+                                            refUC_lib=refUC_tmp, 
+                                            shiftUC_lib=shiftUC_tmp
+                                            )
+
+        # Check matching of symmetries in refUC
+        ind, dt, signs = self.match_symmetries(signs=self.spinor)
+        print(
+              "refUC=\n{}"
+              .format(self.refUC)
+              +"\nshiftUC=\n {}"
+              .format(self.shiftUC)
+              )
+
+        # Sort symmetries like in tables
+        args = np.argsort(ind)
+        for i,i_ind in enumerate(args):
+            self.symmetries[i_ind].ind = i+1
+            self.symmetries[i_ind].sign = signs[i_ind]
+            self.symmetries.append(self.symmetries[i_ind])
+        self.symmetries = self.symmetries[i+1:]
+
+    def show(self, symmetries=None):
         """
         Print description of space-group and symmetry operations.
         
         Parameters
         ----------
-        refUC : array, default=None
-            3x3 array describing the transformation of vectors defining the 
-            unit cell to the standard setting.
-        shiftUC : array, default=np.zeros(3)
-            Translation taking the origin of the unit cell used in the DFT 
-            calculation to that of the standard setting.
         symmetries : int, default=None
             Index of symmetry operations whose description will be printed. 
             Run `IrRep` with flag `onlysym` to check the index corresponding 
@@ -550,13 +619,14 @@ def show(self, refUC=None, shiftUC=np.zeros(3), symmetries=None):
         print('')
         print("\n ---------- INFORMATION ABOUT THE SPACE GROUP ---------- \n")
         print('')
-        print(
-            "Space group # {0} has {1} symmetry operations  ".format(
-                self.number, len(
-                    self.symmetries)))
+        print("Space group {0} (# {1}) has {2} symmetry operations  ".format(
+            self.name,
+            self.number, 
+            len(self.symmetries))
+            )
         for symop in self.symmetries:
             if symmetries is None or symop.ind in symmetries:
-                symop.show(refUC=refUC, shiftUC=shiftUC)
+                symop.show(refUC=self.refUC, shiftUC=self.shiftUC)
 
 
 #  def show2(self,refUC=None,shiftUC=np.zeros(3)):
@@ -568,7 +638,7 @@ def show(self, refUC=None, shiftUC=np.zeros(3), symmetries=None):
 #       symop.show2(refUC=refUC,shiftUC=shiftUC)
 
 
-    def write_trace(self, refUC=None, shiftUC=np.zeros(3)):
+    def write_trace(self):
         """
         Construct description of matrices of symmetry operations of the 
         space-group in the format: 
@@ -576,15 +646,6 @@ def write_trace(self, refUC=None, shiftUC=np.zeros(3)):
         elements of the matrix describing the transformation of the spinor in 
         the format:
         Re(S11),Im(S11),Re(S12),...,Re(S22),Im(S22).
-        
-        Parameters
-        ----------
-        refUC : array, default=None
-            3x3 array describing the transformation of vectors defining the 
-            unit cell to the standard setting.
-        shiftUC : array, default=np.zeros(3)
-            Translation taking the origin of the unit cell used in the DFT 
-            calculation to that of the standard setting.
 
         Returns
         -------
@@ -595,21 +656,12 @@ def write_trace(self, refUC=None, shiftUC=np.zeros(3)):
                # In the following lines, one symmetry operation for each operation of the point group n"""
                ).format(len(self.symmetries))
         for symop in self.symmetries:
-            res += symop.str2(refUC=refUC, shiftUC=shiftUC)
+            res += symop.str2(refUC=self.refUC, shiftUC=self.shiftUC)
         return(res)
 
-    def str(self, refUC=np.eye(3), shiftUC=np.zeros(3)):
+    def str(self):
         """
         Print description of space-group and its symmetry operations.
-        
-        Parameters
-        ----------
-        refUC : array, default=np.eye(3)
-            3x3 array describing the transformation of vectors defining the 
-            unit cell to the standard setting.
-        shiftUC : array, default=np.zeros(3)
-            Translation taking the origin of the unit cell used in the DFT 
-            calculation to that of the standard setting.
 
         Returns
         -------
@@ -626,8 +678,8 @@ def str(self, refUC=np.eye(3), shiftUC=np.zeros(3)):
             "symmetries=\n" +
             "\n".join(
                 s.str(
-                    refUC,
-                    shiftUC) for s in self.symmetries) +
+                    self.refUC,
+                    self.shiftUC) for s in self.symmetries) +
             "\n\n")
 
     def __match_spinor_rotations(self, S1, S2):
@@ -725,18 +777,12 @@ def __gen_refUC():
         '''used somewhere?'''
         nmax = 3
 
-    def get_irreps_from_table(self, refUC, shiftUC, kpname, K):
+    def get_irreps_from_table(self, kpname, K):
         """
         Read irreps of the little-group of a maximal k-point. 
         
         Parameters
         ----------
-        refUC : array
-            3x3 array describing the transformation of vectors defining the 
-            unit cell to the standard setting.
-        shiftUC : array
-            Translation taking the origin of the unit cell used in the DFT 
-            calculation to that of the standard setting.
         kpname : str
             Label of the maximal k-point.
         K : array, shape=(3,)
@@ -747,7 +793,7 @@ def get_irreps_from_table(self, refUC, shiftUC, kpname, K):
         tab : dict
             Each key is the label of an irrep, each value another `dict`. Keys 
             of every secondary `dict` are indices of symmetries (starting from 
-            1 and following order of operations as returned by `spglib`) and 
+            1 and following order of operations in tables of BCS) and 
             values are traces of symmetries.
 
         Raises
@@ -768,7 +814,206 @@ def get_irreps_from_table(self, refUC, shiftUC, kpname, K):
         """
         #        self.show()
         table = IrrepTable(self.number, self.spinor)
-#        if self.name!=table.name     : raise RuntimeError(  "names of the symmetry groups do not match : {0} and {1}".format(self.name,table.name) )
+        tab = {}
+        for irr in table.irreps:
+            if irr.kpname == kpname:
+                k1 = np.round(np.linalg.inv(self.refUC).dot(irr.k), 5) % 1
+                k2 = np.round(K, 5) % 1
+                if not all(np.isclose(k1, k2)):
+                    raise RuntimeError(
+                        "the kpoint {0} does not correspond to the point {1} ({2} in refUC / {3} in primUC) in the table".format(
+                            K,
+                            kpname,
+                            np.round(
+                                irr.k,
+                                3),
+                            k1))
+#            print (irr.characters)
+                tab[irr.name] = {}
+                for i,(sym1,sym2) in enumerate(zip(self.symmetries,table.symmetries)):
+                    try:
+                        dt = sym2.t - sym1.translation_refUC(self.refUC, self.shiftUC)
+                        tab[irr.name][i + 1] = irr.characters[i + 1] * \
+                            sym1.sign * np.exp(2j * np.pi * dt.dot(irr.k))
+                    except KeyError as err:
+                        pass
+        if len(tab) == 0:
+            raise RuntimeError(
+                "the k-point with name {0} is not found in the spacegroup {1}. found only :\n{2}".format(
+                    kpname, table.number, "\n ".join(
+                        "{0}({1}/{2})".format(
+                            irr.kpname, irr.k, np.linalg.inv(self.refUC).dot(
+                                irr.k) %
+                            1) for irr in table.irreps)))
+#            raise RuntimeError("the k-point with name {0} is not found in the spacegroup {1}. found only {2}".format(kpname,table.number,set([irr.kpname for irr in table.irreps]) ) )
+        return tab
+
+#            irr.characters[i]
+# return( { irr.name: np.array([irr.characters[i]*signs[j] for j,i in
+# enumerate(ind)]) for irr in table.irreps if irr.kpname==kpname})
+
+    def determine_basis_transf(self, refUC_cli, shiftUC_cli, refUC_lib, shiftUC_lib):
+        """ 
+        Determine basis transformation to convenctional cell. Priority 
+        is given to the transformation set by the user in CLI.
+
+        Parameters
+        ----------
+        refUC_cli : array
+            3x3 array describing the transformation of vectors defining the 
+            unit cell to the standard setting. Set in CLI.
+        shiftUC_cli : array
+            Translation taking the origin of the unit cell used in the DFT 
+            calculation to that of the standard setting. Set in CLI.
+        refUC_lib : array
+            Obtained via spglib.
+        shiftUC_lib : array
+            Obtained via splib. It may not be the shift taking the 
+            origin to the position adopted in the tables (BCS). For 
+            example, origin choice 1 of ITA is adopted in spglib for 
+            centrosymmetric groups, while origin choice 2 in BCS.
+
+        Returns
+        -------
+        array 
+            Transformation of vectors defining the unit cell to the 
+            standard setting.
+        array
+            Shift taking the origin of the unit cell used in the DFT 
+            calculation to that of the convenctional setting used in 
+            BCS.
+
+        Raises
+        ------
+        RuntimeError
+            Could not find a pait (refUC,shiftUC) matching symmetries 
+            obtained from spglib to those in tables (mod. lattice 
+            translations of the primitive cell).
+
+        """
+        # Set the tranformation to convenctional cell.
+        flag = False
+        if refUC_cli and shiftUC_cli:  # Both specified in CLI.
+            refUC = np.array(refUC_cli.split(","), dtype=float).reshape((3, 3))
+            shiftUC = np.array(shiftUC_cli.split(","), dtype=float).reshape(3)
+            print('refUC and shiftUC read from CLI')
+            return refUC, shiftUC
+        elif refUC_cli and not shiftUC_cli:  # shiftUC not given in CLI.
+            refUC = np.array(refUC_cli.split(","), dtype=float).reshape((3, 3))
+            shiftUC = np.zeros(3, dtype=float)
+            print(('refUC was specified in CLI, but shiftUC was not. Taking '
+                   'shiftUC=(0,0,0).'))
+            return refUC, shiftUC
+        elif not refUC_cli and shiftUC_cli:  # refUC not given in CLI.
+            refUC = np.eye(3, dtype=float)
+            shiftUC = np.array(shiftUC_cli.split(","), dtype=float).reshape(3)
+            print(('shitfUC was specified in CLI, but refUC was not. Taking '
+                   '3x3 identity matrix as refUC.'))
+            return refUC, shiftUC
+        else:  # Neither specifiend in CLI.
+            refUC = refUC_lib.T  # IrRep's treats vecs as column array
+            print(refUC, shiftUC_lib)
+            found = False
+
+            # Check if the group is centrosymmetric
+            inv = None
+            for sym in self.symmetries:
+                if (np.allclose(sym.rotation, -np.eye(3))):
+                    inv = sym
+
+            if inv is None:  # Not centrosymmetric
+                for r_cent in self.vecs_centering():
+                    shiftUC = shiftUC_lib + r_cent
+                    print(shiftUC)
+                    try:
+                        ind, dt, signs = self.match_symmetries(
+                                            refUC,
+                                            shiftUC,
+                                            )
+                        found = True
+                        break
+                    except RuntimeError:
+                        pass
+
+            else:  # Centrosymmetric. Origin must sit in an inv. center
+                for center in self.vecs_inv_centers():
+                    shiftUC = - 0.5 * (inv.translation
+                                       + center.dot(refUC)
+                                       )
+                    try:
+                        ind, dt, signs = self.match_symmetries(
+                                            refUC,
+                                            shiftUC,
+                                            )
+                        found = True
+                        break
+                    except RuntimeError:
+                        pass
+
+            # Print info about center or error if it was not found
+            if found:
+                print("Placing origin at ", -shiftUC)
+                if inv is not None:
+                    print(("where {} = {} + {} x refUC".format(
+                                -shiftUC,
+                                0.5*inv.translation, 
+                                0.5*center)
+                          )
+                           )
+                return refUC, shiftUC
+            else:
+                raise RuntimeError(("Could not find any shift "
+                                    "placing the origin in an "
+                                    "inversion center that leads "
+                                    "to the expressions for "
+                                    "symmetries found in tables."))
+
+
+    def match_symmetries(self, refUC=None, shiftUC=None, signs=False):
+        """
+        Matches symmetry operations of two lists. Translational parts 
+        are matched mod. lattice translations (important for centered 
+        structures).
+
+        Parameters
+        ----------
+        refUC : array, default=None
+            3x3 array describing the transformation of vectors defining the 
+            unit cell to the standard setting.
+        shiftUC : array, default=None
+            Translation taking the origin of the unit cell used in the DFT 
+            calculation to that of the standard setting.
+        signs : bool, default=False
+            If `True`, match also rotations of spinors corresponding to 
+            each symmetry.
+        
+        Returns
+        -------
+        list
+            The :math:`i^{th}` element corresponds to the :math:`i^{th}` 
+            symmetry found by `spglib` and it is the position that the 
+            same symmetry has in the tables.
+        list
+            The :math:`i^{th}` element corresponds to the :math:`i^{th}` 
+            symmetry found by `spglib` and it is the phase difference 
+            w.r.t. the same symmetry in the tables, which may arise if 
+            their translational parts differ by a lattice vector.
+        array
+            The :math:`i^{th}` element corresponds to the :math:`i^{th}` 
+            symmetry found by `spglib` and it is the sign needed to make
+            the matrix for spinor rotation identical to that in tables.
+
+        Note
+        ----
+        Arguments symmetries1 and symmetries2 always take values of 
+        attributes self.symmetry and self.symmetries_tables, so they can
+        be removed, but this way keeps the function more generic.
+        """
+
+        if refUC is None:
+            refUC = self.refUC
+        if shiftUC is None:
+            shiftUC = self.shiftUC
         ind = []
         dt = []
         errtxt = ""
@@ -776,7 +1021,7 @@ def get_irreps_from_table(self, refUC, shiftUC, kpname, K):
             R, t = sym.rotation_refUC(
                 refUC), sym.translation_refUC(refUC, shiftUC)
             found = False
-            for i, sym2 in enumerate(table.symmetries):
+            for i, sym2 in enumerate(self.symmetries_tables):
                 t1 = np.dot(sym2.t - t, refUC) % 1
                 # t1=(sym2.t-t)%1
                 t1[1 - t1 < 1e-5] = 0
@@ -787,13 +1032,6 @@ def get_irreps_from_table(self, refUC, shiftUC, kpname, K):
                         found = True
                         break
                     else:
-                        print(
-                            'table t=',
-                            sym2.t,
-                            '\nfound t=',
-                            t,
-                            "\nt(table)-t(spglib) (mod. lattice translation):",
-                            t1)
                         raise RuntimeError(
                             "Error matching translational part for symmetry " +
                             "{}. A symmetry with identical rotational part \n"
@@ -818,49 +1056,74 @@ def get_irreps_from_table(self, refUC, shiftUC, kpname, K):
                 "Error in matching symmetries detected by spglib with the \
                  symmetries in the tables. Try to modify the refUC and shiftUC \
                  parameters")
-        if self.spinor:
+        if signs:
             S1 = [sym.spinor_rotation for sym in self.symmetries]
-            S2 = [table.symmetries[i].S for i in ind]
-            signs = self.__match_spinor_rotations(S1, S2)
-#            print ("signs = ",signs)
+            S2 = [self.symmetries_tables[i].S for i in ind]
+            signs_array = self.__match_spinor_rotations(S1, S2)
         else:
-            signs = np.ones(len(ind), dtype=int)
+            signs_array = np.ones(len(ind), dtype=int)
+        return ind, dt, signs_array
 
-        tab = {}
-        for irr in table.irreps:
-            if irr.kpname == kpname:
-                k1 = np.round(np.linalg.inv(refUC).dot(irr.k), 5) % 1
-                k2 = np.round(K, 5) % 1
-                if not all(np.isclose(k1, k2)):
-                    raise RuntimeError(
-                        "the kpoint {0} does not correspond to the point {1} ({2} in refUC / {3} in primUC) in the table".format(
-                            K,
-                            kpname,
-                            np.round(
-                                irr.k,
-                                3),
-                            k1))
-#            print (irr.characters)
-                tab[irr.name] = {}
-                for j, i in enumerate(ind):
-                    try:
-                        #                    print (i,j)
-                        tab[irr.name][j + 1] = irr.characters[i + 1] * \
-                            signs[j] * np.exp(2j * np.pi * dt[j].dot(irr.k))
-                    except KeyError as err:
-                        pass
-#        print (tab)
-        if len(tab) == 0:
-            raise RuntimeError(
-                "the k-point with name {0} is not found in the spacegroup {1}. found only :\n{2}".format(
-                    kpname, table.number, "\n ".join(
-                        "{0}({1}/{2})".format(
-                            irr.kpname, irr.k, np.linalg.inv(refUC).dot(
-                                irr.k) %
-                            1) for irr in table.irreps)))
-#            raise RuntimeError("the k-point with name {0} is not found in the spacegroup {1}. found only {2}".format(kpname,table.number,set([irr.kpname for irr in table.irreps]) ) )
-        return tab
+    def vecs_centering(self):
+        """ 
+        Check the space group and generate vectors of centering.
 
-#            irr.characters[i]
-# return( { irr.name: np.array([irr.characters[i]*signs[j] for j,i in
-# enumerate(ind)]) for irr in table.irreps if irr.kpname==kpname})
+        Returns
+        -------
+        array
+            Each row is a lattice vector describing the centering.
+        """
+        cent = np.array([[0,0,0]])
+        if self.name[0] == 'P':
+            pass  # Just to make it explicit
+        elif self.name[0] == 'C':
+            cent = np.vstack((cent, cent + [1/2,1/2,0]))
+        elif self.name[0] == 'I':
+            cent = np.vstack((cent, cent + [1/2,1/2,1/2]))
+        elif self.name[0] == 'F':
+            cent = np.vstack((cent,
+                              cent + [0,1/2,1/2],
+                              cent + [1/2,0,1/2],
+                              cent + [1/2,1/2,0],
+                              )
+                             )
+        elif self.name[0] == 'A':  # test this
+            cent = np.vstack((cent, cent + [0,1/2,1/2]))
+        else:  # R-centered
+            cent = np.vstack((cent,
+                              cent + [2/3,1/3,1/3],
+                              cent + [1/3,2/3,2/3],
+                              )
+                             )
+        return cent
+
+    def vecs_inv_centers(self):
+        """
+        Get the positions of all inversion centers in the unit cell.
+
+        Returns
+        -------
+        array
+            Each element is a vector pointing to a position center in 
+            the convenctional unit cell.
+
+        Notes
+        -----
+        If the space group is primitive, there are 8 inversion centers, 
+        but there are more if it is a group with a centering.
+
+        """
+        vecs = np.array(
+                        [
+                        [0,0,0],
+                        [1,0,0],
+                        [0,1,0],
+                        [0,0,1],
+                        [1,1,0],
+                        [1,0,1],
+                        [0,1,1],
+                        [1,1,1]
+                        ]
+                        )
+        vecs = np.vstack([vecs + r for r in self.vecs_centering()])
+        return vecs
diff --git a/irreptables/__init__.py b/irreptables/__init__.py
index a590c1f25..c58add592 100644
--- a/irreptables/__init__.py
+++ b/irreptables/__init__.py
@@ -346,7 +346,7 @@ def __init__(self, SGnumber, spinor, name=None):
             elif l[0].lower() == "spinor":
                 assert str2bool(l[1]) == self.spinor
             elif l[0].lower() == "symmetries":
-                print("reading symmetries")
+                print("reading symmetries from tables")
                 self.symmetries = []
                 while len(self.symmetries) < self.nsym:
                     l = lines.pop()


Test patch:
diff --git a/irrep/tests/test_cli.py b/irrep/tests/test_cli.py
index 19bfe37e1..42aad14cb 100644
--- a/irrep/tests/test_cli.py
+++ b/irrep/tests/test_cli.py
@@ -25,6 +25,31 @@ def test_bi_scalar_example():
     return_code = output.returncode
     stdout = output.stdout
 
+    with open("test_out", "w") as f:
+        f.write(stdout)
+
+    assert return_code == 0, output.stderr
+    assert "GM1+(1.0)" in stdout, stdout
+    assert "GM2-(1.0)" in stdout, stdout
+    assert "GM1+(1.0)" in stdout, stdout
+    assert "GM3+(1.0)" in stdout, stdout
+    
+    command = [
+        "irrep",
+        "-Ecut=50",
+        "-code=abinit",
+        "-fWFK=O_DS2_WFK",
+        "-kpoints=1",
+        "-IBstart=11",
+        "-IBend=15",
+        "-kpnames=GM",
+    ]
+
+    output = subprocess.run(command, capture_output=True, text=True)
+
+    return_code = output.returncode
+    stdout = output.stdout
+
     with open("test_out", "w") as f:
         f.write(stdout)
 
@@ -72,6 +97,31 @@ def test_bi_spinor_example():
     assert "-GM8(1.0)" in stdout, stdout
     assert "-GM9(1.0)" in stdout, stdout
 
+    command = [
+        "irrep",
+        "-spinor",
+        "-Ecut=50",
+        "-code=abinit",
+        "-fWFK=O_DS2_WFK",
+        "-kpoints=1",
+        "-IBend=21",
+        "-IBend=32",
+        "-kpnames=GM",
+    ]
+
+    output = subprocess.run(command, capture_output=True, text=True)
+
+    with open("out") as f:
+        reference_out = f.read()
+
+    return_code = output.returncode
+    stdout = output.stdout
+
+    assert return_code == 0, output.stderr
+
+    assert "-GM8(1.0)" in stdout, stdout
+    assert "-GM9(1.0)" in stdout, stdout
+
     for test_output_file in (
         "irreps.dat",
         "irreptable-template",
@@ -210,19 +260,38 @@ def test_bi_hoti():
         f.write(stdout)
 
     assert return_code == 0, output.stderr
+    assert "-2.7306  |        2     | -GM8(1.0)" in stdout, stdout
+    assert "-0.7762  |        2     | -GM8(1.0)" in stdout, stdout
+    assert "-0.4961  |        2     | -GM4(1.0), -GM5(1.0)" in stdout, stdout
+    assert "-4.8263  |        2     | -F3(1.0), -F4(1.0)" in stdout, stdout
+    assert "-3.6784  |        2     | -F5(1.0), -F6(1.0)" in stdout, stdout
+    assert "-2.4303  |        2     | -F5(1.0), -F6(1.0)" in stdout, stdout
+    assert "-1.7054  |        2     | -L5(1.0), -L6(1.0)" in stdout, stdout
+    assert "-1.6885  |        2     | -L3(1.0), -L4(1.0)" in stdout, stdout
+    assert "-0.1312  |        2     | -L5(1.0), -L6(1.0)" in stdout, stdout
+    assert "-1.5597  |        2     | -T9(1.0)" in stdout, stdout
+    assert "-1.2220  |        2     | -T8(1.0)" in stdout, stdout
+    assert "0.1460  |        2     | -T6(1.0), -T7(1.0)" in stdout, stdout
 
-#     known_output = """k-point   2 :[0. 0. 0.] 
-#  number of irreps = 6
-#    Energy  | multiplicity |        irreps        | sym. operations  
-#            |              |                      |     1        2        3        4        5        6        7        8        9       10       11       12    
-#   -2.7306  |        2     | -GM8(1.0)            |  2.0000   2.0000   1.0000   1.0000   1.0000   1.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000
-#   -0.7762  |        2     | -GM8(1.0)            |  2.0000   2.0000   1.0000   1.0000   1.0000   1.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000
-#   -0.4961  |        2     | -GM4(1.0), -GM5(1.0) |  2.0000   2.0000  -2.0000  -2.0000  -2.0000  -2.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000
-# inversion is # 2
-# number of inversions-odd Kramers pairs :  0"""
+    command = [
+        "irrep",
+        "-spinor",
+        "-code=vasp",
+        "-kpnames=T,GM,F,L",
+        "-Ecut=50",
+        "-EF=5.2156",
+        "-IBstart=5",
+        "-IBend=10"
+    ]
+    output = subprocess.run(command, capture_output=True, text=True)
 
-    # assert known_output in stdout, stdout
+    return_code = output.returncode
+    stdout = output.stdout
 
+    with open("test_out", "w") as f:
+        f.write(stdout)
+
+    assert return_code == 0, output.stderr
     assert "-2.7306  |        2     | -GM8(1.0)" in stdout, stdout
     assert "-0.7762  |        2     | -GM8(1.0)" in stdout, stdout
     assert "-0.4961  |        2     | -GM4(1.0), -GM5(1.0)" in stdout, stdout