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

Issue:
(VASP) Unfolding issue with irrep-1.9.1
Dear Admin,

I am using irrep-1.9.1 with banduppy to unfold the band structure from VASP WAVECAR file. However, I am encountering a few issues:

1. The unfolding process is significantly slower compared to irrep-1.8.3.
2. I receive a non-integer warning for each state: `WARNING - non-integer number of states: 411.4891663397381`.
3. Consequently, the unfolding results are incorrect, with all the Bloch weights accumulating at the top bands.

Do you have any suggestions to address these issues?

Best regards,  
Badal

![unfolded_bandstructure_](https://github.com/stepan-tsirkin/irrep/assets/56228126/0dcd26cf-0d65-4edf-9493-f899e65f6d2f)

Reading and writing of prefix.sym file for quantum espresso 
Implements an option to write a prefix.sym file, or to read the symmetry perations from it.

This is intended to be used with pw2wannier90. Not  with banduppy, and not with identification of names of irreps.
Just get the characters of symmetry operations, in a desired order. 

The .sym file is natively written by quantum espresso, but for trigonal systems, before pw2wannier90, one needs to wrun nscf calculation with nosym=True, and therefore the file contains only the identity, which is useless. Therefore I need to write the file with irrep

On the other hand, reading symmetries form the file allows to customize the order or the set of symmetry operations (exclude some of them) 




Discussion / hints:
Hi @bmondal94 ,
really looks like something's broken. Just to be clear, is it true that with irrep-1.8.3 things work as expected?

second, is it only with VASP, or similar problem arises with other code (e.g. QE or abinit)

@MIraola , any ideas?
Dear Admins,
1. Yes.  irrep-1.8.3 works perfectly as expected. All my tested unfolded band structures came out as expected.
2. I have only tested with VASP so far. Unfortunately, at the moment I do not have access to QE or Abinit. So, I can only confirm that the problem is with VASP.
Hi @bmondal94 ,

As I understand, there is no advantage in using irrep-1.9 or newer with banduppy for now, right?

In that case I would suggest that in banduppy for now you you fix in requirements that irrep<=1.8.3, and also when running banduppy raise a warning if the installed version of irrep happens to be newer.  

After we resolve this compatibility and performance problem, we remove this restriction in banduppy. 
Hi Badal,

thank you for your report. Unfortunately, this week I will not be able to
fix the problem because I am very busy with preparing some documents. I
will look at the problem at the beginning of next week. Meanwhile, please
install irrep version <= 1.8.3 as Stepan suggested.

Sorry for the inconvenience.

Best,
Mikel I.

El jue, 11 jul 2024 a las 14:29, Stepan Tsirkin ***@***.***>)
escribió:

> Hi @bmondal94 <https://github.com/bmondal94> ,
>
> As I understand, there is no advantage in using irrep-1.9 or newer with
> banduppy for now, right?
>
> In that case I would suggest that in banduppy for now you you fix in
> requirements that irrep<=1.8.3, and also when running banduppy raise a
> warning if the installed version of irrep happens to be newer.
>
> After we resolve this compatibility and performance problem, we remove
> this restriction in banduppy.
>
> —
> Reply to this email directly, view it on GitHub
> <https://github.com/stepan-tsirkin/irrep/issues/73#issuecomment-2222808213>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AJUEUPHSFU2BFLOEI4RF2JDZLZ3BLAVCNFSM6AAAAABKQIV3OOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMRSHAYDQMRRGM>
> .
> You are receiving this because you were mentioned.Message ID:
> ***@***.***>
>

@bmondal94 , regarding the issue that the code became slower. Could you run the unfolding with 

```
python -m cProfile script.py
```

both with irrep-1.9.1 and 1.8.3, and attach the outputs here, to identify the new bottlenecks
I think I know where the issue hides. 

In this line of `Kpoint.__init__()`

https://github.com/MIraola/irrep/blob/a2d50b9b2daa794044576d9bbeac2b72a9f2c5f2/irrep/kpoint.py#L201-L202

the traces are calculated (which are not actually needed for unfolding). But in addition self.Energy is substituted by Energy_mean

https://github.com/MIraola/irrep/blob/a2d50b9b2daa794044576d9bbeac2b72a9f2c5f2/irrep/kpoint.py#L600-L602

which is returned by calculate_traces

https://github.com/MIraola/irrep/blob/a2d50b9b2daa794044576d9bbeac2b72a9f2c5f2/irrep/kpoint.py#L615

Now, when `Kpoint.unfold()` calls for `self.get_rhi_spin` here 

https://github.com/MIraola/irrep/blob/a2d50b9b2daa794044576d9bbeac2b72a9f2c5f2/irrep/kpoint.py#L310
The energies are still assumed not-averaged 

So, in these borders 

https://github.com/MIraola/irrep/blob/a2d50b9b2daa794044576d9bbeac2b72a9f2c5f2/irrep/kpoint.py#L366-L372

no degeneracies are detected, because the degenerate states are already grouped. 

Now, for simplicity assume that each band is doubly-degenerate, i.e. we read energies as [E1,E1, E2,E2,E3E3,...] and wavefunctions [WF1_1, WF1_2, WF2_1, WF2_2, ... ]

Now after averaging the Energies become [E1,E2,E3], but the wavefunctions do not change. So, __eval_rho_spin considers the pairs (energy,wavefunction) as 

```python 
[ (E1, WF1_1) , (E2, WF1_1) , (E3, WF2_1) , (E4, WF2_2), ...] 
```

So, a higher energy is attributed to each  wavefunction, therefore all unfolding weights tend to be at the top of the figure.

This issue is not VASP-related, it should be generall for codes, as long as there are some degeneracies in the bands.

So, question  - how to solve it?  I think, we should keep both Energy_raw, and Energy_mean. 

On the other hand, what I **do not** like in the fact that traces are calculated in `Kpoint.__init__()` is that they are calculated once for a fixed threshold, and cannot be recalculated . 

So, what is the point of the getter 

https://github.com/MIraola/irrep/blob/a2d50b9b2daa794044576d9bbeac2b72a9f2c5f2/irrep/kpoint.py#L320
if anyway, it may be called only once, with one threshold ?

Also, I recently learned  that those getters are useless, they can be replaced by `@functools.lru_cache` decorator

@bmondal94 ,

I tried to fix the problem on this branch, could you try it:

https://github.com/stepan-tsirkin/irrep/tree/fix-issue-73
Dear Stepan,

First of all, I like the explanation on why the Bloch weights were appearing on the top bands. 

The [fix-issue-73](https://github.com/stepan-tsirkin/irrep/tree/fix-issue-73) branch is working as expected. In terms of unfolding time, it is also fast similar to irrep-1.8.3. However, I want to point out a few points:

1. I noticed that there have been modifications made to variable names in latest irrep. For example, in the `Kpoints` module's `unfold()` function, `self.K` has been changed to `self.k`. This breaks BandUPpy, as we were previously accessing the unfolded k-points using `.K`. I have updated our code, but I want to highlight that changing the names of especially the `self` parameters can break the irrep-dependent software if they access these `self` parameters for processing. 
With this update, the unfolding is working, but I need to re-check all other irrep `self` parameters we are accessing in BandUPpy to be sure.

2. With the [fix-issue-73](https://github.com/stepan-tsirkin/irrep/tree/fix-issue-73) branch the Bloch weights in the unfolded band structure have slightly changed (in comparison to irrep-1.8.3) in floating-point precision, although mostly in the 4th-8th decimal places, which may not be significant.

Best regards,
Badal
Hi, If I am correct, we access the Kpoint.K attribute only here in banduppy

https://github.com/band-unfolding/banduppy/blob/67728db444fbc2e8dddd6a6bcb95116d277248a6/banduppy/src/unfolding.py#L163

```python 
             if is_round(KP.K - self.kpointsSBZ[key, :3], prec = 1e-6):
```
I suggest that for now in irrep we keep a wrapper property `K` to keep compatibility:

https://github.com/stepan-tsirkin/irrep/blob/4c592e7e90abe39c61f1d6bfc65c169c81206e4e/irrep/kpoint.py#L226-L246

also the `k_close_mod1` function is added to Kpoint, that can de accessed by future version of banduppy directly. This way we both avoid using .K/.k properties, and the import of the auxiliary function `is_round` from irrep.

@bmondal94 , could you list all the properties/functions in irrep that are accessed directly by banduppy. We may put there a warning which should read smth like "ACCESSED DIRECTLY BY BANDUPPY. PLEASE DO NOT CHANGE WITHOUT NOTIFYING ITS DEVELOPERS"

Well, maybe it is better to put it as a comment, rather than in a docstring?

Also, a good idea is to add a short test (to irrep tests) which uses banduppy, so that any compatibility-breaking change will not be unnoticed.  



Hi @stepan-tsirkin,
Just to let you know, I will be quite busy this week. I may not be able to get back to you guys with the requirements you mentioned. I will try my best to come back to you with a positive response ASAP.
Hi,

1- I think it's a good idea to label the attributes that banduppy accesses directly, as Stepan suggested. I think that the section of "Attributes" in the docstring of the class is the correct place to do it, since if we add a comment in a line of the code then we have to check all instances of the attribute before changing the name.

2- Originally, we used to calculate the traces only for a single energy threshold. This is convenient if we want to get rid of wave functions; in order to recalculate the traces, we need to save the wave functions. I will tackle this point once we implement this branch in the code (probably, I will do it by creating a flag `save_wf`, since this possibility of recalculating the traces is more likely to be used when `BandStructure` or `Kpoints` are called in a jupyter-notebook environment). 

>  I think that the section of "Attributes" in the docstring of the class is the correct place to do it,

:+1: 

>  probably, I will do it by creating a flag `save_wf`

:+1:   yes, it is better to control it with an explicit tag like that


Dear @stepan-tsirkin @MIraola ,
Thanks a lot for looking into the issues. Here is the collections of `irrep` functions that `banduppy` uses:

Notes: 
1. Update require in  `Kpoint.k_close_mod1(self, kpt, prec=1e-6)` return. ==> `return is_round(self.k - kpt, prec = prec)`

__FYI:__ `is_round()` is called within `Kpoint.unfold` as well which can be replaced with `k_close_mod1()` functions.

2. There is a `print` statement in `gvectors.calc_gvectors()` which is called by `bandstructure.BandStructure`.
How about a bit more information about what it is printing? If the print statement is kept there we can do like -
```
if N % 10 == 0:
    print(f'cycle-{N:>3d}: No. of plane waves = {len(igall):>10d}')
```
3. Suggestion (not significant): I noticed `self.RecLattice` is calculated when `spacegroup.SpaceGroup` is called. Within `bandstructure.BandStructure`
`self.RecLattice` is again calculated despite `spacegroup.SpaceGroup` is called earlier within 
`bandstructure.BandStructure`. Recalculating `self.RecLattice` can be avoided here.

===============================================
```
class Kpoint:
    @property
    def K(self):
        """Getter for the redfuced coordinates of the k-point
        needed to keep compatibility with banduppy<=0.3.3.
        
        ACCESSED DIRECTLY BY BANDUPPY<=0.3.3. 
        BANDUPPY>=0.3.4 NO LONGER USES IT.
        - banduppy.src.unfolding._BandUnfolding._gather_generated_ab_calculated_kpts()
        DO NOT CHANGE UNLESS NECESSARY.
        PLEASE NOTIFY THE DEVELOPERS IF ANY CHANGES ARE MADE.
        """
        return self.k
    
    def k_close_mod1(self, kpt, prec=1e-6):
        """
        Check if the k-point is close to another k-point modulo 1. (in reduced coordinates)
        
        ACCESSED DIRECTLY BY BANDUPPY>=0.3.4. 
        - banduppy.src.unfolding._BandUnfolding._gather_generated_ab_calculated_kpts()
        DO NOT CHANGE UNLESS NECESSARY.
        PLEASE NOTIFY THE DEVELOPERS IF ANY CHANGES ARE MADE.

        Parameters
        ----------
        kpt : array
            Coordinates of the k-point to compare.
        prec : float, default=1e-6
            Threshold to consider the k-points as equal.
        """
        return is_round(self.k - kpt, prec = prec)
    
class BandStructure:
    def __init__(...):
        """
        Parses files and organizes info about the whole band structure in 
        attributes. Contains methods to calculate and write traces (and irreps), 
        for the separation of the band structure in terms of a symmetry operation 
        and for the calculation of the Zak phase and wannier charge centers.
        
        'self.kpoints' IS ACCESSED DIRECTLY BY BANDUPPY>=0.3.4. 
        - banduppy.src.unfolding._BandUnfolding._unfold_bandstructure()
        DO NOT CHANGE UNLESS NECESSARY.
        PLEASE NOTIFY THE DEVELOPERS IF ANY CHANGES ARE MADE.
        """
        .
        .
        .
        kp = Kpoint(...)
        self.kpoints.append(kp)
        
    def KPOINTSline(self, kpred=None, supercell=None, breakTHRESH=0.1):
            """
            Calculate cumulative length along a k-path in cartesian coordinate.
            
            o k_vec = u.b1 + v.b2 + w.b3
            o e.g. (u,v,w)=(1/2, 0, 0) -> k_vec==(kx, ky, kz)
            o [kx ky kz] = [u v w].[[b11 b12 b13], [b21 b22 b23], [b31 b32 b33]]
            
            ACCESSED DIRECTLY BY BANDUPPY>=0.3.4. 
            - banduppy.src.unfolding._BandUnfolding._unfold_bandstructure()
            DO NOT CHANGE UNLESS NECESSARY.
            PLEASE NOTIFY THE DEVELOPERS IF ANY CHANGES ARE MADE.
    
            Parameters
            ----------
            kpred : list, default=None
                Each element contains the direct coordinates of a k-point in the
                attribute `kpoints`.
            supercell : array, shape=(3,3), default=None
                Describes how the lattice vectors of the (super)cell used in the 
                calculation are expressed in the basis vectors of the primitive 
                cell. 
            breakTHRESH : float, default=0.1
                If the distance between two neighboring k-points in the path is 
                larger than `break_thresh` break continuity in k-path. Set breakTHRESH 
                to a large value if the unfolded kpoints line is continuous.
                The default is 0.1.
    
            Returns
            -------
            K : array
                Each element is the cumulative distance along the path up to a 
                k-point. The first element is 0, so that the number of elements
                matches the number of k-points in the path.
            """
            if kpred is None:
                kpred = [k.k for k in self.kpoints]
            
            if supercell is None:
                reciprocal_lattice = self.RecLattice
            else:
                reciprocal_lattice = supercell.T @ self.RecLattice 
                
            KPcart = np.dot(kpred, reciprocal_lattice) 
            K = np.zeros(KPcart.shape[0])
            k = np.linalg.norm(KPcart[1:, :] - KPcart[:-1, :], axis=1)
            k[k > breakTHRESH] = 0.0
            K[1:] = np.cumsum(k)
            return K
        
```
Thank you for this info @bmondal94 . Definitely, it is very useful. I will implement the required changes immediately. 
@MIraola 
One more addition to earlier post. There is a 'warning' msg printed during unfolding:
```
Neither refUC nor shiftUC were specified in CLI and searchcell was False.
Taking 3x3 identity matrix as refUC and shiftUC=(0,0,0). 
If you want to calculate the transformation to conventional cell, run IrRep with -searchcell ???
```
The msg is originating from `irrep.spacegroup.SpaceGroup.determine_basis_transf()`. This is probably not of any concerns but does it make sense to suppress this msg during unfolding? Otherwise, this may cause confusion for banduppy users.


Fix patch:
diff --git a/irrep/__init__.py b/irrep/__init__.py
index dfd773aca..ae98917ce 100644
--- a/irrep/__init__.py
+++ b/irrep/__init__.py
@@ -1,2 +1,2 @@
-__version__ = "1.9.2"
+__version__ = "1.9.3"
 #from .bandstructure import BandStructure
diff --git a/irrep/bandstructure.py b/irrep/bandstructure.py
index 3ec88f966..7ba4f8b31 100644
--- a/irrep/bandstructure.py
+++ b/irrep/bandstructure.py
@@ -26,6 +26,7 @@
 from .kpoint import Kpoint
 from .spacegroup import SpaceGroup
 from .gvectors import sortIG, calc_gvectors
+from .utility import log_message
 
 
 class BandStructure:
@@ -85,6 +86,12 @@ class BandStructure:
     save_wf : bool
         Whether wave functions should be kept as attribute after calculating 
         traces.
+    v : int, default=0
+        Number controlling the verbosity. 
+        0: minimalistic printing. 
+        1: print info about decisions taken internally by the code, recommended 
+        when the code runs without errors but the result is not the expected.
+        2: print detailed info, recommended when the code stops with an error
 
     Attributes
     ----------
@@ -113,7 +120,9 @@ class BandStructure:
     kpoints : list
         Each element is an instance of `class Kpoint` corresponding to a 
         k-point specified in input parameter `kpoints`. If this input was not 
-        set, all k-points found in DFT files will be considered.
+        set, all k-points found in DFT files will be considered. ACCESSED 
+        DIRECTLY BY BANDUPPY>=0.3.4. DO NOT CHANGE UNLESS NECESSARY. NOTIFY 
+        THE DEVELOPERS IF ANY CHANGES ARE MADE.
     num_bandinvs : int
         Property that returns the number of inversion odd states in the 
         given TRIM.
@@ -149,6 +158,7 @@ def __init__(
         trans_thresh=1e-5,
         degen_thresh=1e-8,
         save_wf=True,
+        v=0,
         alat=None,
         from_sym_file=None,
     ):
@@ -167,7 +177,7 @@ def __init__(
                 )
             self.spinor = spinor
             parser = ParserVasp(fPOS, fWAV, onlysym)
-            self.Lattice, positions, typat = parser.parse_poscar()
+            self.Lattice, positions, typat = parser.parse_poscar(v)
             if not onlysym:
                 NK, NBin, self.Ecut0, lattice = parser.parse_header()
                 if not np.allclose(self.Lattice, lattice):
@@ -184,7 +194,7 @@ def __init__(
              self.spinor,
              typat,
              positions,
-             EF_in) = parser.parse_header()
+             EF_in) = parser.parse_header(v=v)
             NBin = max(nband)
 
         elif code == "espresso":
@@ -234,15 +244,18 @@ def __init__(
                               shiftUC=shiftUC,
                               search_cell=search_cell,
                               trans_thresh=trans_thresh,
-                              alat=alat, from_sym_file=from_sym_file)
+                              v=v,
+                              alat=alat,
+                              from_sym_file=from_sym_file)
         if onlysym:
             return
 
         # Set Fermi energy
         if EF.lower() == "auto":
             if EF_in is None:
-                print("WARNING : fermi-energy not found. Setting it as zero")
                 self.efermi = 0.0
+                msg = "WARNING : fermi-energy not found. Setting it as 0 eV"
+                log_message(msg, v, 1)
             else:
                 self.efermi = EF_in
         else:
@@ -251,7 +264,8 @@ def __init__(
             except:
                 raise RuntimeError("Invalid value for keyword EF. It must be "
                                    "a number or 'auto'")
-        print("Efermi: {:.4f} eV".format(self.efermi))
+
+        log_message(f"Efermi: {self.efermi:.4f} eV", v, 1)
 
         # Fix indices of bands to be considered
         if IBstart is None or IBstart <= 0:
@@ -277,13 +291,14 @@ def __init__(
         self.RecLattice *= (2.0*np.pi/np.linalg.det(self.Lattice))
 
         # To do: create writer of description for this class
-        print(
-            "WAVECAR contains {0} k-points and {1} bands.\n Saving {2} bands starting from {3} in the output".format(
-                NK, NBin, NBout, IBstart + 1
-            )
-        )
-        print("Energy cutoff in WAVECAR : ", self.Ecut0)
-        print("Energy cutoff reduced to : ", self.Ecut)
+        msg = ("WAVECAR contains {} k-points and {} bands.\n"
+               "Saving {} bands starting from {} in the output"
+               .format(NK, NBin, NBout, IBstart + 1))
+        log_message(msg, v, 1)
+        msg = f"Energy cutoff in WAVECAR : {self.Ecut0}"
+        log_message(msg, v, 1)
+        msg = f"Energy cutoff reduced to : {self.Ecut}"
+        log_message(msg, v, 1)
 
         # Create list of indices for k-points
         if kplist is None:
@@ -297,13 +312,16 @@ def __init__(
         for ik in kplist:
 
             if code == 'vasp':
+                msg = f'Parsing wave functions at k-point #{ik:>3d}'
+                log_message(msg, v, 2)
                 WF, Energy, kpt, npw = parser.parse_kpoint(ik, NBin, self.spinor)
                 kg = calc_gvectors(kpt,
                                    self.RecLattice,
                                    self.Ecut0,
                                    npw,
                                    self.Ecut,
-                                   spinor=self.spinor
+                                   spinor=self.spinor,
+                                   v=v
                                    )
                 if not self.spinor:
                     selectG = kg[3]
@@ -314,11 +332,15 @@ def __init__(
             elif code == 'abinit':
                 NBin = parser.nband[ik]
                 kpt = parser.kpt[ik]
+                msg = f'Parsing wave functions at k-point #{ik:>3d}: {kpt}'
+                log_message(msg, v, 2)
                 WF, Energy, kg = parser.parse_kpoint(ik)
                 WF, kg = sortIG(ik, kg, kpt, WF, self.RecLattice, self.Ecut0, self.Ecut, self.spinor)
 
             elif code == 'espresso':
-                WF, Energy, kg, kpt = parser.parse_kpoint(ik, NBin, spin_channel)
+                msg = f'Parsing wave functions at k-point #{ik:>3d}'
+                log_message(msg, v, 2)
+                WF, Energy, kg, kpt = parser.parse_kpoint(ik, NBin, spin_channel, v=v)
                 WF, kg = sortIG(ik+1, kg, kpt, WF, self.RecLattice/2.0, self.Ecut0, self.Ecut, self.spinor)
 
             elif code == 'wannier90':
@@ -329,9 +351,12 @@ def __init__(
                                    self.RecLattice,
                                    self.Ecut,
                                    spinor=self.spinor,
-                                   nplanemax=np.max([ngx, ngy, ngz]) // 2
+                                   nplanemax=np.max([ngx, ngy, ngz]) // 2,
+                                   v=v
                                    )
                 selectG = tuple(kg[0:3])
+                msg = f'Parsing wave functions at k-point #{ik:>3d}: {kpt}'
+                log_message(msg, v, 2)
                 WF = parser.parse_kpoint(ik+1, selectG)
 
             # Pick energy of IBend+1 band to calculate gaps
@@ -360,7 +385,8 @@ def __init__(
                 refUC=self.spacegroup.refUC,
                 shiftUC=self.spacegroup.shiftUC,
                 symmetries_tables=self.spacegroup.symmetries_tables,
-                save_wf=save_wf
+                save_wf=save_wf,
+                v=v
                 )
             self.kpoints.append(kp)
         del WF
@@ -370,7 +396,7 @@ def num_k(self):
         '''Getter for the number of k points'''
         return len(self.kpoints)
 
-    def identify_irreps(self, kpnames):
+    def identify_irreps(self, kpnames, v=0):
         '''
         Identifies the irreducible representations of wave functions based on 
         the traces of symmetries of the little co-group. Each element of 
@@ -381,12 +407,14 @@ def identify_irreps(self, kpnames):
         ----------
         kpnames : list
             List of labels of the maximal k points.
+        v : int, default=0
+            Verbosity level. Default set to minimalistic printing
         '''
 
         for ik, KP in enumerate(self.kpoints):
             
             if kpnames is not None:
-                irreps = self.spacegroup.get_irreps_from_table(kpnames[ik], KP.k)
+                irreps = self.spacegroup.get_irreps_from_table(kpnames[ik], KP.k, v=v)
             else:
                 irreps = None
             KP.identify_irreps(irreptable=irreps)
@@ -600,7 +628,7 @@ def write_trace(self,):
                 KP.write_trace()
             )
 
-    def Separate(self, isymop, degen_thresh=1e-5, groupKramers=True):
+    def Separate(self, isymop, degen_thresh=1e-5, groupKramers=True, v=0):
         """
         Separate band structure according to the eigenvalues of a symmetry 
         operation.
@@ -613,6 +641,8 @@ def Separate(self, isymop, degen_thresh=1e-5, groupKramers=True):
             Energy threshold used to determine degeneracy of energy-levels.
         groupKramers : bool, default=True
             If `True`, states will be coupled by Kramers' pairs.
+        v : int, default=0
+            Verbosity level. Default is set to minimalistic printing
 
         Returns
         -------
@@ -669,7 +699,7 @@ def Separate(self, isymop, degen_thresh=1e-5, groupKramers=True):
                 return dict({allvalues.mean(): self})
         else:
             allvalues = allvalues[np.argsort(np.angle(allvalues))]
-            print("allvalues:", allvalues)
+            log_message(f'allvalues: {allvalues}', v, 1)
             borders = np.where(abs(allvalues - np.roll(allvalues, 1)) > 0.01)[0]
             nv = len(allvalues)
             if len(borders) > 0:
@@ -679,7 +709,7 @@ def Separate(self, isymop, degen_thresh=1e-5, groupKramers=True):
                         for b1, b2 in zip(borders, np.roll(borders, -1))
                     ]
                 )
-                print("distinct values:", allvalues)
+                log_message(f'Distinct values: {allvalues}', v, 1)
                 subspaces = {}
                 for v in allvalues:
                     other = copy.copy(self)
@@ -687,8 +717,6 @@ def Separate(self, isymop, degen_thresh=1e-5, groupKramers=True):
                     for K in kpseparated:
                         vk = list(K.keys())
                         vk0 = vk[np.argmin(np.abs(v - vk))]
-                        #                    print ("v,vk",v,vk)
-                        #                    print ("v,vk",v,vk[np.argmin(np.abs(v-vk))])
                         if abs(vk0 - v) < 0.05:
                             other.kpoints.append(K[vk0])
                         subspaces[v] = other
@@ -801,29 +829,40 @@ def write_plotfile(self, filename='bands-tognuplot.dat'):
                 file.write(s)
         file.close()
 
-    def KPOINTSline(self, kpred=None, breakTHRESH=0.1):
+    def KPOINTSline(self, kpred=None, supercell=None, breakTHRESH=0.1):
         """
-        Calculate cumulative length along a path in reciprocal space.
+        Calculate cumulative length along a k-path in cartesian coordinates.
+        ACCESSED DIRECTLY BY BANDUPPY>=0.3.4. DO NOT CHANGE UNLESS NECESSARY. 
+        NOTIFY THE DEVELOPERS IF ANY CHANGE IS MADE.
 
         Parameters
         ----------
         kpred : list, default=None
             Each element contains the direct coordinates of a k-point in the
             attribute `kpoints`.
+        supercell : array, shape=(3,3), default=None
+                Describes how the lattice vectors of the (super)cell used in the 
+                calculation are expressed in the basis vectors of the primitive 
+                cell. USED IN BANDUPPY. DO NOT CHANGE UNLESS NECESSARY.
         breakTHRESH : float, default=0.1
             If the distance between two neighboring k-points in the path is 
-            larger than `breakTHRESH`, it is taken to be 0.
+            larger than `breakTHRESH`, it is taken to be 0. Set `breakTHRESH` 
+            to a large value if the unforlded kpoints line is continuous.
 
         Returns
         -------
-        array
+        K : array
             Each element is the cumulative distance along the path up to a 
             k-point. The first element is 0, so that the number of elements
             matches the number of k-points in the path.
         """
         if kpred is None:
             kpred = [k.k for k in self.kpoints]
-        KPcart = np.array(kpred).dot(self.RecLattice)
+        if supercell is None:
+            reciprocal_lattice = self.RecLattice
+        else:
+            reciprocal_lattice = supercell.T @ self.RecLattice 
+        KPcart = np.dot(kpred, reciprocal_lattice)
         K = np.zeros(KPcart.shape[0])
         k = np.linalg.norm(KPcart[1:, :] - KPcart[:-1, :], axis=1)
         k[k > breakTHRESH] = 0.0
diff --git a/irrep/cli.py b/irrep/cli.py
index da173e542..e755bfe0a 100644
--- a/irrep/cli.py
+++ b/irrep/cli.py
@@ -25,7 +25,7 @@
 from monty.serialization import dumpfn, loadfn
 
 from .bandstructure import BandStructure
-from .utility import str2list, short
+from .utility import str2list, short, log_message
 from . import __version__ as version
 
 
@@ -271,6 +271,16 @@ def invoke(self, ctx):
           "Default: 1e-5"
           )
 )
+@click.option("-v",
+              count=True,
+              default=1,
+              help=("Verbosity flag. -vv: very detailed info, recommended "
+                    "when you get an error. -v (default for CLI): info about "
+                    "some decissions taken internaly through the code's "
+                    "execution, recommended when the code runs without "
+                    "errors, but the result is not what you expected. If you "
+                    "don't set this tag, you will get the basic info.")
+)
 def cli(
     ecut,
     fwav,
@@ -301,7 +311,8 @@ def cli(
     config,
     searchcell,
     correct_ecut0,
-    trans_thresh
+    trans_thresh,
+    v
 ):
     """
     Defines the "irrep" command-line tool interface.
@@ -317,17 +328,20 @@ def cli(
     # Warning about kpnames
     if kpnames is not None:
         searchcell = True
+
     elif not searchcell:
-        print(("Warning: transformation to the convenctional unit "
+        msg = ("Warning: transformation to the convenctional unit "
                "cell will not be calculated, nor its validity checked "
                "(if given). See the description of flag -searchcell "
                "on:\n"
                "irrep --help"
-               ))
-        print(("Warning: -kpnames not specified. Only traces of "
+               )
+        log_message(msg, v, 1)
+        msg = ("Warning: -kpnames not specified. Only traces of "
                "symmetry operations will be calculated. Remember that "
                "-kpnames must be specified to identify irreps"
-               ))
+               )
+        log_message(msg, v, 1)
 
     # parse input arguments into lists if supplied
     if symmetries:
@@ -363,7 +377,8 @@ def cli(
         search_cell = searchcell,
         degen_thresh=degenthresh,
         save_wf=save_wf,
-        from_sym_file=from_sym_file,
+        v=v,
+        from_sym_file=from_sym_file
     )
 
     bandstr.spacegroup.show()
@@ -380,7 +395,7 @@ def cli(
                 )
 
     # Identify irreps. If kpnames wasn't set, all will be labelled as None
-    bandstr.identify_irreps(kpnames)
+    bandstr.identify_irreps(kpnames, v=v)
 
     # Temporary, until we make it valid for isymsep
     bandstr.write_characters()
@@ -408,7 +423,7 @@ def cli(
         for isym in isymsep:
             print("\n-------- SEPARATING BY SYMMETRY # {} --------".format(isym))
             for s_old, bs in subbands.items():
-                separated = bs.Separate(isym, degen_thresh=degenthresh, groupKramers=groupkramers)
+                separated = bs.Separate(isym, degen_thresh=degenthresh, groupKramers=groupkramers, v=v)
                 for s_new, bs_separated in separated.items():
                     tmp_subbands[tuple(list(s_old) + [s_new])] = bs_separated
             subbands = tmp_subbands
diff --git a/irrep/gvectors.py b/irrep/gvectors.py
index 2c7cc39d9..85bf2b49e 100644
--- a/irrep/gvectors.py
+++ b/irrep/gvectors.py
@@ -19,6 +19,7 @@
 import numpy as np
 import numpy.linalg as la
 from .readfiles import Hartree_eV
+from .utility import log_message
 
 
 class NotSymmetryError(RuntimeError):
@@ -45,6 +46,7 @@ def calc_gvectors(
     thresh=1e-3,
     spinor=True,
     nplanemax=10000,
+    v=0
 ):
     """ 
     Generates G-vectors taking part in the plane-wave expansion of 
@@ -85,6 +87,10 @@ def calc_gvectors(
         the same energy as the plane-wave of the current column.
     
 """
+
+    msg = ('Generating plane waves at k: ({} )'
+           .format(' '.join([f'{x:6.3f}' for x in K])))
+    log_message(msg, v, 2)
     if Ecut1 <= 0:
         Ecut1 = Ecut
     B = RecLattice
@@ -96,7 +102,8 @@ def calc_gvectors(
     for N in range(nplanemax):
         flag = True
         if N % 10 == 0:
-            print(N, len(igall))
+            msg = f'Cycle {N:>3d}: number of plane waves = {len(igall):>10d}'
+            log_message(msg, v, 2)
         if len(igall) >= nplane / 2:    # Only enters if vasp
             if spinor:
                 break
diff --git a/irrep/kpoint.py b/irrep/kpoint.py
index faf618919..0d6a2d7d8 100644
--- a/irrep/kpoint.py
+++ b/irrep/kpoint.py
@@ -21,7 +21,7 @@
 import numpy.linalg as la
 import copy
 from .gvectors import symm_eigenvalues, symm_matrix
-from .utility import compstr, is_round, format_matrix
+from .utility import compstr, is_round, format_matrix, log_message
 
 class Kpoint:
     """
@@ -82,6 +82,8 @@ class Kpoint:
     save_wf : bool
         Whether wave functions should be kept as attribute after calculating 
         traces.
+    v : int
+        Verbosity level. Default set to minimalistic printing
     
     Attributes
     ----------
@@ -105,8 +107,12 @@ 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.
-    K : array, shape=(3,)
+    k : array, shape=(3,)
         Direct coordinates of the k point in the DFT cell setting.
+    K : array, shape=(3,)
+        Property getter for `self.k`. NEEDED TO PRESERVE COMPATIBILITY WITH
+        BANDUPPY<=0.3.3. DO NOT CHANGE UNLESS NECESSARY. NOTIFY THE DEVELOPERS 
+        IF ANY CHANGES ARE MADE.
     k_refUC : array, shape=(3,)
         Direct coordinates of the k point in the reference cell setting.
     Energy_raw : array
@@ -166,7 +172,8 @@ def __init__(
         refUC=np.eye(3),
         shiftUC=np.zeros(3),
         symmetries_tables=None,  # calculate_traces needs it
-        save_wf=True
+        save_wf=True,
+        v=0
     ):
         self.spinor = spinor
         self.ik0 = ik + 1  # the index in the WAVECAR (count start from 1)
@@ -211,7 +218,7 @@ def __init__(
 
         # Calculate traces
         if calculate_traces:
-            self.char, self.char_refUC, self.Energy_mean = self.calculate_traces(refUC, shiftUC, symmetries_tables, degen_thresh)
+            self.char, self.char_refUC, self.Energy_mean = self.calculate_traces(refUC, shiftUC, symmetries_tables, v)
 
             # Determine number of band inversions based on parity
             found = False
@@ -237,14 +244,18 @@ def K(self):
         """Getter for the redfuced coordinates of the k-point
         needed to keep compatibility with banduppy
         
-        ACCESSED BY BANDUPPY, AVOID CHANGING UNLESS NECESSARY
+        ACCESSED BY BANDUPPY.NEEDED TO PRESERVE COMPATIBILITY WITH 
+        BANDUPPY<=0.3.3. AVOID CHANGING UNLESS NECESSARY.
         """
         return self.k
     
     def k_close_mod1(self, kpt, prec=1e-6):
         """
-        Check if the k-point is close to another k-point modulo 1. (in reduced coordinates)
-        ACCESSED BY BANDUPPY, AVOID CHANGING UNLESS NECESSARY
+        Check if the k-point is close to another k-point modulo 1. (in reduced 
+        coordinates)
+
+        ACCESSED BY BANDUPPY. AVOID CHANGING UNLESS NECESSARY. NOTIFY 
+        DEVELOPERS IF ANY CHANGE IS MADE.
 
         Parameters
         ----------
@@ -253,7 +264,7 @@ def k_close_mod1(self, kpt, prec=1e-6):
         prec : float, default=1e-6
             Threshold to consider the k-points as equal.
         """
-        return is_round(self.k - kpt, prec = 1e-6)
+        return is_round(self.k - kpt, prec = prec)
 
     def copy_sub(self, E, WF, inds):
         """
@@ -323,7 +334,7 @@ def unfold(self, supercell, kptPBZ, degen_thresh=1e-4):
             W - weight of the band(s) projected onto the PBZ kpoint.
             Sx, Sy, Sz - Spin components projected onto the PBZ kpoint.
         """
-        if not is_round(kptPBZ.dot(supercell.T) - self.k, prec=1e-5):
+        if not self.k_close_mod1(kptPBZ.dot(supercell.T), prec=1e-5):
             raise RuntimeError(
                 "unable to unfold {} to {}, withsupercell={}".format(
                     self.k, kptPBZ, supercell
@@ -425,7 +436,7 @@ def get_rho_spin(self, degen_thresh=1e-4):
                 result.append((b1, b2, E, (W,)))
         return result
 
-    def Separate(self, symop, degen_thresh, groupKramers=True):
+    def Separate(self, symop, degen_thresh, groupKramers=True, v=0):
         """
         Separate the band structure in a particular k-point according to the 
         eigenvalues of a symmetry operation.
@@ -438,6 +449,8 @@ def Separate(self, symop, degen_thresh, groupKramers=True):
             Energy threshold used to determine degeneracy of energy-levels.
         groupKramers : bool, default=True
             If `True`, states will be coupled by pairs of Kramers.
+        v : int, default=0
+            Verbosity level. Default set to minimalistic printing
 
         Returns
         -------
@@ -452,11 +465,9 @@ def Separate(self, symop, degen_thresh, groupKramers=True):
         norms = self.WF.conj().dot(self.WF.T)
         check = np.max(abs(norms - np.eye(norms.shape[0])))
         if check > 1e-5:
-            print(
-                "orthogonality (largest of diag. <psi_nk|psi_mk>): {0:7.5} > 1e-5   \n".format(
-                    check
-                )
-            )
+            msg = ("orthogonality (largest of diag. <psi_nk|psi_mk>): "
+                   "{0:7.5} > 1e-5   \n".format(check))
+            log_message(msg, v, 1)
 
         S = symm_matrix(
             self.k,
@@ -476,10 +487,12 @@ def Separate(self, symop, degen_thresh, groupKramers=True):
             Sblock[b1:b2, b1:b2] = 0
         check = np.max(abs(Sblock))
         if check > 0.1:
-            print(("WARNING: matrix of symmetry has non-zero elements between "
-                   "states of different energy:  \n", check))
-            print("Printing matrix of symmetry at k={}".format(self.k))
-            print(format_matrix(Sblock))
+            msg = ("WARNING: matrix of symmetry has non-zero elements between "
+                   "states of different energy:  \n", check)
+            log_message(msg, v, 1)
+            msg = (f"Printing matrix of symmetry at k={self.k}")
+            log_message(msg, v, 1)
+            log_message(format_matrix(Sblock), v, 1)
 
         # Calculate eigenvalues and eigenvectors in each block
         eigenvalues = []
@@ -504,7 +517,7 @@ def Separate(self, symop, degen_thresh, groupKramers=True):
 
         # Check unitarity of the symmetry
         if np.abs((np.abs(w) - 1.0)).max() > 1e-4:
-            print("WARNING : some eigenvalues are not unitary :{0} ".format(w))
+            log_message(f"WARNING: some eigenvalues are not unitary: {w}",v,1)
         if np.abs((np.abs(w) - 1.0)).max() > 3e-1:
             raise RuntimeError(" some eigenvalues are not unitary :{0} ".format(w))
         w /= np.abs(w)
@@ -527,7 +540,6 @@ def Separate(self, symop, degen_thresh, groupKramers=True):
             if len(borders) > 0:
                 for b1, b2 in zip(borders, borders[1:]):
                     v1 = v[:, b1:b2]
-                    print(w[b1:b2].mean())
                     subspaces[w[b1:b2].mean()] = self.copy_sub(E=Eloc[b1:b2], WF=v1.T.dot(self.WF), inds=inds_states[b1:b2])
             else:
                 v1 = v
@@ -555,7 +567,7 @@ def Separate(self, symop, degen_thresh, groupKramers=True):
 
         return subspaces
 
-    def calculate_traces(self, refUC, shiftUC, symmetries_tables, degen_thresh=1e-8):
+    def calculate_traces(self, refUC, shiftUC, symmetries_tables, v=0):
         '''
         Calculate traces of symmetry operations
 
@@ -607,7 +619,8 @@ def calculate_traces(self, refUC, shiftUC, symmetries_tables, degen_thresh=1e-8)
         # Check that number of irreps is int
         Nirrep = np.linalg.norm(char.sum(axis=1)) ** 2 / char.shape[0]
         if abs(Nirrep - round(Nirrep)) > 1e-2:
-            print("WARNING - non-integer number of states : {0}".format(Nirrep))
+            msg = f"WARNING - non-integer number of states : {Nirrep}"
+            log_message(msg, v, 2)
         Nirrep = int(round(Nirrep))
 
         # Sum traces of degenerate states. Rows (cols) correspond to states (syms)
diff --git a/irrep/readfiles.py b/irrep/readfiles.py
index 34d899ef3..5682ee64b 100644
--- a/irrep/readfiles.py
+++ b/irrep/readfiles.py
@@ -21,7 +21,7 @@
 from scipy.io import FortranFile as FF
 from sys import stdout
 from .utility import FortranFileR as FFR
-from .utility import str2bool, BOHR, split
+from .utility import str2bool, BOHR, split, log_message
 import xml.etree.ElementTree as ET
 
 Rydberg_eV = 13.605693  # eV
@@ -126,11 +126,16 @@ def __init__(self, filename):
         self.fWFK = FFR(filename)  # temporary
         self.kpt_count = 0  # index of the next k-point to be read
 
-    def parse_header(self):
+    def parse_header(self, v=0):
         '''
         Parse header of WFK file and save as attributes quantities that 
         will be used in the rest of methods
 
+        Parameters
+        ----------
+        v : int, default=0
+            Verbosity level. Default set to minimal printing
+
         Returns
         -------
         nband : array
@@ -165,10 +170,10 @@ def parse_header(self):
         headform, fform = record[0][1]
         defversion = ['8.6.3', '9.6.2', '8.4.4', '8.10.3']
         if codsvn not in defversion:
-            print(("WARNING, the version {0} of abinit is not in {1} "
+            msg = ("WARNING, the version {0} of abinit is not in {1} "
                    "and may not be fully tested"
                    .format(codsvn, defversion))
-                  )
+            log_message(msg, v, 1)
         if headform < 80:
             raise ValueError(
                 "Head form {0}<80 is not supported".format(headform)
@@ -312,7 +317,6 @@ def parse_kpoint(self, ik):
             used in the expansion of wave functions
         '''
 
-        print("Reading k-point", ik)
         nspinor = 2 if self.spinor else 1
 
         # We need to skip lines in fWFK until we reach the lines of ik
@@ -392,10 +396,15 @@ def __init__(self, fPOS, fWAV, onlysym=False):
         if not onlysym:
             self.fWAV = WAVECARFILE(fWAV)
 
-    def parse_poscar(self):
+    def parse_poscar(self, v=0):
         """
         Parses POSCAR.
 
+        Parameters
+        ----------
+        v : int, default=0
+            Verbosity level. Default set to minimal printing
+
         Returns
         ------
         lattice : array
@@ -407,6 +416,8 @@ def parse_poscar(self):
             Each element is a number identifying the atomic species of an ion.
         """
 
+        msg = f'Reading POSCAR: {self.fPOS}'
+        log_message(msg, v, 1)
         fpos = (l.strip() for l in open(self.fPOS))
         title = next(fpos)  # title
         del title
@@ -437,7 +448,7 @@ def parse_poscar(self):
                 positions[i] = np.array(l.split()[:3])
                 i += 1
             except Exception as err:
-                print(err)
+                log_message(err, msg, 1)
                 pass
         if sum(nat) != i:
             raise RuntimeError(
@@ -576,7 +587,8 @@ def parse_header(self):
             NBin_dw=int(self.bandstr.find('nbnd_dw').text)
             NBin_up=int(self.bandstr.find('nbnd_up').text)
             spinpol=True
-            print ("spin-polarised bandstructure composed of {} up and {} dw states".format(NBin_dw,NBin_up))
+            print("spin-polarised bandstructure composed of {} up and {} dw "
+                  "states".format(NBin_up,NBin_dw))
             NBin_dw+NBin_up
         except AttributeError:
             spinpol=False
@@ -646,7 +658,7 @@ def parse_lattice(self):
         return lattice, positions, typat, alat
 
 
-    def parse_kpoint(self, ik, NBin, spin_channel):
+    def parse_kpoint(self, ik, NBin, spin_channel, v=0):
         '''
         Parse block of a particular k-point from `data-file-schema.xml` file
 
@@ -671,6 +683,8 @@ def parse_kpoint(self, ik, NBin, spin_channel):
             reciprocal lattice vector
         kpt : array
             Direct coords of the k-point w.r.t. DFT cell vectors
+        v : int, default=0
+            Verbosity level. Default set to minimalistic printing
         '''
 
         kptxml = self.bandstr.findall("ks_energies")[ik]
@@ -703,7 +717,8 @@ def parse_kpoint(self, ik, NBin, spin_channel):
         # Parse coefficients of wave functions
         npw = int(kptxml.find("npw").text)
         npwtot = npw * (2 if self.spinor else 1)
-        print('npwtot: {}, igwx: {}'.format(npwtot, igwx))
+        msg = 'npwtot: {}, igwx: {}'.format(npwtot, igwx)
+        log_message(msg, v, 2)
         WF = np.zeros((NBin, npwtot), dtype=complex)
         for ib in range(NBin):
             rec = record_abinit(fWFC, "{}f8".format(npwtot * 2))
diff --git a/irrep/spacegroup.py b/irrep/spacegroup.py
index 6a1847ce4..e37419e02 100644
--- a/irrep/spacegroup.py
+++ b/irrep/spacegroup.py
@@ -23,7 +23,7 @@
 import spglib
 from irreptables import IrrepTable
 from scipy.optimize import minimize
-from .utility import str_, BOHR
+from .utility import str_, log_message, BOHR
 
 pauli_sigma = np.array(
     [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]])
@@ -500,6 +500,8 @@ class SpaceGroup():
     from_sym_file : str, default=None
         If provided, the symmetry operations are read from this file.
         (format of pw2wannier90 prefix.sym  file)
+    v : int, default=0
+        Verbosity level. Default set to minimalistic printing
 
     Attributes
     ----------
@@ -519,9 +521,6 @@ class SpaceGroup():
     Lattice : array, shape=(3,3) 
         Each row contains cartesian coordinates of a basis vector forming the 
         unit-cell in real space.
-    RecLattice : array, shape=(3,3)
-        Each row contains the cartesian coordinates of a basis vector forming 
-        the unit-cell in reciprocal space.
     positions : array
         Direct coordinate of sites in the DFT cell setting.
     typat : list
@@ -561,7 +560,8 @@ def __init__(
             search_cell=False,
             trans_thresh=1e-5,
             alat=None,
-            from_sym_file=None
+            from_sym_file=None,
+            v=0
             ):
         self.spinor = spinor
         self.Lattice = cell[0]
@@ -572,19 +572,19 @@ def __init__(
          self.number, 
          refUC_tmp, 
          shiftUC_tmp) = self._findsym(cell, from_sym_file, alat)
-        self.RecLattice = 2*np.pi*np.linalg.inv(self.Lattice).T
         self.order = len(self.symmetries)
         self.alat=alat
 
         # Determine refUC and shiftUC according to entries in CLI
-        self.symmetries_tables = IrrepTable(self.number, self.spinor).symmetries
+        self.symmetries_tables = IrrepTable(self.number, self.spinor, v=v).symmetries
         self.refUC, self.shiftUC = self.determine_basis_transf(
                                             refUC_cli=refUC, 
                                             shiftUC_cli=shiftUC,
                                             refUC_lib=refUC_tmp, 
                                             shiftUC_lib=shiftUC_tmp,
                                             search_cell=search_cell,
-                                            trans_thresh=trans_thresh
+                                            trans_thresh=trans_thresh,
+                                            v=v
                                             )
 
         # Check matching of symmetries in refUC. If user set transf.
@@ -613,10 +613,11 @@ def __init__(
                 # User specified refUC or shiftUC in CLI. He/She may
                 # want the traces in a cell that is not neither the
                 # one in tables nor the DFT one
-                print(("WARNING: refUC and shiftUC don't transform the cell to "
+                msg = ("WARNING: refUC and shiftUC don't transform the cell to "
                        "one where symmetries are identical to those read from "
                        "tables. If you want to achieve the same cell as in "
-                       "tables, try not specifying refUC and shiftUC."))
+                       "tables, try not specifying refUC and shiftUC.")
+                log_message(msg, v, 1)
                 pass
 
     def _findsym(self, cell, from_sym_file, alat):
@@ -661,30 +662,30 @@ def _findsym(self, cell, from_sym_file, alat):
             centrosymmetric groups they adopt origin choice 1 of ITA, rather 
             than choice 2 (BCS).
         """
-        dataset = spglib.get_symmetry_dataset(cell)
+        dataset: spglib.SpglibDataset | None = spglib.get_symmetry_dataset(cell)
         lattice = cell[0]
         if from_sym_file is not None:
             assert alat is not None, "Lattice parameter must be provided to read symmetries from file"
             rot_cart, trans_cart = read_sym_file(from_sym_file)
             rotations, translations = cart_to_crystal(rot_cart, trans_cart, lattice, alat )
         else:
-            rotations = dataset['rotations']
-            translations = dataset['translations']
+            rotations = dataset.rotations
+            translations = dataset.translations
         symmetries = [
             SymmetryOperation(
                 rot,
                 translations[i],
-                lattice,
+                cell[0],
                 ind=i + 1,
                 spinor=self.spinor) for i,
             rot in enumerate(
                 rotations)]
 
         return (symmetries, 
-                dataset['international'],
-                dataset['number'], 
-                dataset['transformation_matrix'],
-                dataset['origin_shift']
+                dataset.international,
+                dataset.number,
+                dataset.transformation_matrix,
+                dataset.origin_shift
                 )
 
     def json(self, symmetries=None):
@@ -936,7 +937,7 @@ def residue(x):
         return np.array([R1.dot(b).dot(R1.T.conj()).dot(np.linalg.inv(
             a)).diagonal().mean().real.round() for a, b in zip(S1, S2)], dtype=int)
 
-    def get_irreps_from_table(self, kpname, K):
+    def get_irreps_from_table(self, kpname, K, v=0):
         """
         Read irreps of the little-group of a maximal k-point. 
         
@@ -946,6 +947,8 @@ def get_irreps_from_table(self, kpname, K):
             Label of the maximal k-point.
         K : array, shape=(3,)
             Direct coordinates of the k-point.
+        v : int, default=0
+            Verbosity level. Default set to minimalistic printing
 
         Returns
         -------
@@ -972,7 +975,7 @@ def get_irreps_from_table(self, kpname, K):
             given in parameter `kpname`.
         """
 
-        table = IrrepTable(self.number, self.spinor)
+        table = IrrepTable(self.number, self.spinor, v=v)
         tab = {}
         for irr in table.irreps:
             if irr.kpname == kpname:
@@ -1012,7 +1015,8 @@ def determine_basis_transf(
             refUC_lib,
             shiftUC_lib,
             search_cell,
-            trans_thresh
+            trans_thresh,
+            v=0
             ):
         """ 
         Determine basis transformation to conventional cell. Priority
@@ -1038,6 +1042,8 @@ def determine_basis_transf(
             It is `True` if kpnames was specified in CLI.
         trans_thresh : float, default=1e-5
             Threshold to compare translational parts of symmetries.
+        v : int, default=0
+            Verbosity level. Default set to minimalistic printing
 
         Returns
         -------
@@ -1063,31 +1069,34 @@ def determine_basis_transf(
         if refUC_cli_bool and shiftUC_cli_bool:  # Both specified in CLI.
             refUC = refUC_cli.T  # User sets refUC as if it was acting on column
             shiftUC = shiftUC_cli
-            print('refUC and shiftUC read from CLI')
+            log_message('refUC and shiftUC read from CLI', v, 1)
             return refUC, shiftUC
         elif refUC_cli_bool and not shiftUC_cli_bool:  # shiftUC not given in CLI.
             refUC = refUC_cli.T  # User sets refUC as if it was acting on column
             shiftUC = np.zeros(3, dtype=float)
-            print(('refUC was specified in CLI, but shiftUC was not. Taking '
-                   'shiftUC=(0,0,0).'))
+            msg = ('refUC was specified in CLI, but shiftUC was not. Taking '
+                   'shiftUC=(0,0,0)')
+            log_message(msg, v, 1)
             return refUC, shiftUC
         elif not refUC_cli_bool and shiftUC_cli_bool:  # refUC not given in CLI.
             refUC = np.eye(3, dtype=float)
             shiftUC = shiftUC_cli
-            print(('shitfUC was specified in CLI, but refUC was not. Taking '
-                   '3x3 identity matrix as refUC.'))
+            msg = ('shitfUC was specified in CLI, but refUC was not. Taking '
+                   '3x3 identity matrix as refUC.')
+            log_message(msg, v, 1)
             return refUC, shiftUC
         elif not search_cell:
             refUC = np.eye(3, dtype=float)
             shiftUC = np.zeros(3, dtype=float)
-            print(('Neither refUC nor shiftUC were specified in CLI '
-                   'and searchcell was False. Taking 3x3 identity '
-                   'matrix as refUC and shiftUC=(0,0,0). If you want '
-                   'to calculate the transformation to conventional '
-                   'cell, run IrRep with -searchcell'
-                   ))
+            msg = ('Taking 3x3 identity matrix as refUC and shiftUC=(0,0,0). '
+                   'If you want to calculate the transformation to '
+                   'conventional cell, run IrRep with -searchcell')
+            log_message(msg, v, 1)
             return refUC, shiftUC
         else:  # Neither specifiend in CLI.
+            msg = ('Determining transformation to conventional setting '
+                   '(refUC and shiftUC)')
+            log_message(msg, v, 1)
             refUC = np.linalg.inv(refUC_lib)  # from DFT to convenctional cell
 
             # Check if the shift given by spglib works
@@ -1117,9 +1126,8 @@ def determine_basis_transf(
                                             shiftUC,
                                             trans_thresh=trans_thresh
                                             )
-                        print(('ShiftUC achieved with the centering: {}'
-                                   .format(r_center))
-                              )
+                        msg = (f'ShiftUC achieved with the centering: {r_center}')
+                        log_message(msg, v, 1)
                         return refUC, shiftUC
                     except RuntimeError:
                         pass
@@ -1136,14 +1144,13 @@ def determine_basis_transf(
                                             shiftUC,
                                             trans_thresh=trans_thresh
                                             )
-                        print(('ShiftUC achieved in 2 steps:\n'
+                        msg = ('ShiftUC achieved in 2 steps:\n'
                                '  (1) Place origin of primitive cell on '
                                'inversion center: {}\n'
                                '  (2) Move origin of convenctional cell to the '
                                'inversion-center: {}'
-                               .format(0.5 * inv.translation, r_center)
-                               )
-                              )
+                               .format(0.5 * inv.translation, r_center))
+                        log_message(msg, v, 1)
                         return refUC, shiftUC
                     except RuntimeError:
                         pass
@@ -1154,7 +1161,6 @@ def determine_basis_transf(
                                     "shiftUC in command line"))
 
 
-
     def match_symmetries(
             self,
             refUC=None,
diff --git a/irrep/utility.py b/irrep/utility.py
index 0f63a7b1e..17efc4425 100644
--- a/irrep/utility.py
+++ b/irrep/utility.py
@@ -249,4 +249,22 @@ def format_matrix(A):
         "   ".join("{0:+5.2f} {1:+5.2f}".format(x.real, x.imag) for x in a)
         + "\n"
         for a in A
-    )
\ No newline at end of file
+    )
+
+
+def log_message(msg, verbosity, level):
+    '''
+    Logger to decide if a message is printed or not
+
+    Parameters
+    ----------
+    msg : str
+        Message to print
+    verbosity : int
+        Verbosity set for the current run of the code
+    level : int
+        If `verbosity >= level`, the message will be printed
+    '''
+
+    if verbosity >= level:
+        print(msg)
\ No newline at end of file
diff --git a/irreptables/irreptables/__init__.py b/irreptables/irreptables/__init__.py
index c90e7c102..61743d8a6 100644
--- a/irreptables/irreptables/__init__.py
+++ b/irreptables/irreptables/__init__.py
@@ -14,12 +14,12 @@
 ##  e-mail: stepan.tsirkin@physik.uzh.ch                         #
 ##################################################################
 
-__version__="1.0.3"
+__version__="1.1.0"
 
 import os
 import logging
 import numpy as np
-from irrep.utility import str2bool, str2list_space, str_
+from irrep.utility import str2bool, str2list_space, str_, log_message
 
 # using a logger to print useful information during debugging,
 # set to logging.INFO to disable debug messages
@@ -211,6 +211,8 @@ class Irrep:
         info about the space-group and irreps.
     k_point : class instance, default=None
         Instance of class `KPoint`.
+    v : int, default=0
+        Verbosity level. Default set to minimalistic printing
 
     Attributes
     ----------
@@ -232,7 +234,7 @@ class Irrep:
         and the corresponding value is the trace of that symmetry in the irrep.
     """
 
-    def __init__(self, line, k_point):
+    def __init__(self, line, k_point, v=0):
         logger.debug("reading irrep line <{0}> for KP=<{1}> ".format(line, k_point.str()))
         self.k = k_point.k
         self.kpname = k_point.name
@@ -249,7 +251,7 @@ def __init__(self, line, k_point):
                 * np.array(line[2 + self.nsym : 2 + 2 * self.nsym], dtype=float)
             )
         self.characters = {k_point.isym[i]: ch[i] for i in range(self.nsym)}
-        logger.debug("the irrep {0}  ch= {1}".format(self.name, self.characters))
+        log_message(f"## Irrep {self.name}\nCharacter:\n{self.characters}", v, 2)
         assert len(self.characters) == self.nsym
 
     def show(self):
@@ -295,6 +297,8 @@ class IrrepTable:
         Name of the file from which info about the space-group and irreps 
         should be read. If `None`, the code will try to open a file already 
         included in it.
+    v : int, default=0
+        Verbosity level. Default set to minimalistic printing
 
     Attributes
     ----------
@@ -316,7 +320,7 @@ class IrrepTable:
         irrep of the little group of a maximal k-point.
     """
 
-    def __init__(self, SGnumber, spinor, name=None):
+    def __init__(self, SGnumber, spinor, name=None, v=0):
         self.number = SGnumber
         self.spinor = spinor
         if name is None:
@@ -325,10 +329,13 @@ def __init__(self, SGnumber, spinor, name=None):
                 spinor="spin" if self.spinor else "scal",
                 root=os.path.dirname(__file__),
             )
-            logger.debug("reading from a standard irrep table <{0}>".format(name))
+            msg = f"Reading standard irrep table <{name}>"
+            log_message(msg, v, 2)
         else:
-            logger.debug("reading from a user-defined irrep table <{0}>".format(name))
+            msg = f"Reading a user-defined irrep table <{name}>"
+            log_message(msg, v, 2)
 
+        log_message("\n---------- DATA FROM THE TABLE ----------\n", v, 2)
         lines = open(name).readlines()[-1::-1]
         while len(lines) > 0:
             l = lines.pop().strip().split("=")
@@ -342,7 +349,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 from tables")
+                log_message("Reading symmetries from tables", v, 2)
                 self.symmetries = []
                 while len(self.symmetries) < self.nsym:
                     l = lines.pop()
@@ -354,21 +361,27 @@ def __init__(self, SGnumber, spinor, name=None):
                         pass
                 break
 
-        logger.debug("symmetries are:\n" + "\n".join(s.str() for s in self.symmetries))
+        msg = "Symmetries are:\n" + "\n".join(s.str() for s in self.symmetries)
+        log_message(msg, v, 2)
 
         self.irreps = []
         while len(lines) > 0:
             l = lines.pop().strip()
             try:
                 kp = KPoint(line=l)
-                logger.debug("kpoint successfully read:", kp.str())
+                msg = f"k-point successfully read:\n{kp.str()}"
+                log_message(msg, v, 2)
             except Exception as err:
-                logger.debug("error while reading k-point <{0}>".format(l), err)
                 try:
-                    self.irreps.append(Irrep(line=l, k_point=kp))
+                    self.irreps.append(Irrep(line=l, k_point=kp, v=v))
                 except Exception as err:
-                    logger.debug("error while reading irrep <{0}>".format(l), err)
-                    pass
+                    if len(l.split()) > 0:
+                        msg = ("WARNING: could not parse k-point nor irrep from the "
+                               "following line <\n{}>"
+                               .format(l))
+                        log_message(msg, v, 2)
+                    else:
+                        pass
 
     def show(self):
         '''


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 364abb0f8..e7bbb18a9 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -31,6 +31,10 @@ jobs:
         python -m pip install --upgrade pip
         pip install pytest
         pip install -e .
+    - name: Install tables from local repo, not from pypi
+      working-directory: ./irreptables
+      run: |
+        pip install -e .
     - name: pytest
       run: |
         pytest --color=yes
diff --git a/irrep/tests/test_espresso_spinor.py b/irrep/tests/test_espresso_spinor.py
index 93fa9ce0a..d0d659e5e 100644
--- a/irrep/tests/test_espresso_spinor.py
+++ b/irrep/tests/test_espresso_spinor.py
@@ -4,9 +4,6 @@
 from monty.serialization import loadfn
 import numpy as np
 
-import irrep
-import irrep.bandstructure
-
 TEST_FILES_PATH = Path(__file__).parents[2] / "examples"
 
 def test_espresso_spinor_example():