diff --git a/abipy/core/wyckoff.py b/abipy/core/site_symmetries.py
similarity index 57%
rename from abipy/core/wyckoff.py
rename to abipy/core/site_symmetries.py
index 429e2296..3da1155f 100644
--- a/abipy/core/wyckoff.py
+++ b/abipy/core/site_symmetries.py
@@ -8,214 +8,11 @@ import numpy as np
 import sympy as sp
 
 #from monty.functools import lazy_property
-from collections import Sequence, OrderedDict
+from collections import OrderedDict
 from monty.termcolor import cprint
 from abipy.core.mixins import Has_Structure
 
 
-def _validate_params(params, wpos):
-    """
-    Check keys in params. Raises ValueError.
-    """
-    check = params.copy()
-    for k in params:
-        if k not in wpos.names:
-            raise ValueError("Cannot find key: %s in symbol_names: %s" % (k, wpos.names))
-        check.pop(k)
-    if check:
-        raise ValueError("Found unknown symbol names in %s" % check)
-
-
-class WyckoffPositions(Sequence):
-    """
-    A numpy array with sympy expression.
-    """
-    @classmethod
-    def from_string(cls, s):
-        """
-        Initialize the object from string `s`.
-        """
-        # Build a matrix of strings with expressions.
-        s = s.replace("\n", "")
-        s = s.replace("(", "").replace(")", "")
-        s = s.replace("[", "").replace("]", "")
-
-        try:
-            matrix_str = np.reshape(s.split(","), (-1, 3))
-        except ValueError:
-            raise ValueError("Wrong string in input, perhaps missing comma. Input string:\n %s" % str(s))
-
-        return cls(matrix_str)
-
-    def __init__(self, matrix_str, letter=None, site_symmetry=None):
-        """
-
-        Args:
-            matrix_str:
-            letter:
-            site_symmetry:
-        """
-        self.letter, self.site_symmetry = None, None
-
-        # Build numpy matrix of sympy expressions from the matrix of strings.
-        from sympy import Symbol
-        from sympy.parsing.sympy_parser import parse_expr
-        exprs = []
-        for vec in np.reshape(matrix_str, (-1, 3)):
-            exprs.extend([parse_expr(e) for e in vec])
-
-        self._expr_mat = np.reshape(exprs, (-1, 3))
-
-        # Get set of sympy symbols.
-        symbols = []
-        for expr in self.ravel():
-            for arg in expr.args:
-                if isinstance(arg, Symbol): symbols.append(arg)
-
-        # SymPy symbols and names.
-        self.symbols = list(set(symbols))
-        self.names = [s.name for s in self.symbols]
-
-    def __len__(self):
-        return len(self._expr_mat)
-
-    def __getitem__(self, key):
-        return self._expr_mat.__getitem__(key)
-
-    def __str__(self):
-        """String representation."""
-        lines = []
-        for vec in self:
-            lines.append("(" + ",".join(str(e) for e in vec) + ")")
-        return ", ".join(lines)
-
-    def ravel(self):
-        """Flat list, similar to np.ravel"""
-        return self._expr_mat.ravel()
-
-    @property
-    def mult(self):
-        """Multiplicity"""
-        return len(self)
-
-    @property
-    def num_params(self):
-        """The number of parameters."""
-        return len(self.symbols)
-
-    def params_from_frac_coords(self, frac_coords):
-        """
-        Compute the Wyckoff parameters from the fractional coordinates.
-        """
-        frac_coords = np.reshape(frac_coords, (-1, 3))
-        assert len(self) == len(frac_coords)
-        equations = []
-        for i in range(len(self)):
-            eqs = self[i, :] - frac_coords[i, :]
-            equations.extend(eqs)
-
-        # Solve the problem.
-        from sympy.solvers import solve
-        d = solve(equations, self.symbols)
-        #print("solution: ", d)
-
-        # Return results in a dict with symbol names.
-        if not d:
-            return None
-
-        return {symbol.name: d[symbol] for symbol in self.symbols}
-
-    def error_of_params(self, params, frac_coords):
-        """
-        Return a numpy array with the difference between the
-        wyckoff positions computed from the sympy expressions and the input frac_coords
-        """
-        _validate_params(params, self)
-        frac_coords = np.reshape(frac_coords, (-1, 3))
-        assert len(self) == len(frac_coords)
-
-        frac_errors = np.empty((self.mult, 3))
-        for i in range(len(self)):
-            for j in range(3):
-                frac_errors[i,j] = self[i,j].subs(params) - frac_coords[i,j]
-
-        # Modulo lattice vector (with tolerance?)
-        return frac_errors % 1
-
-
-class Wyckoff(object):
-    """
-    """
-    #@classmethod
-    #def from_site2wpos_ordict(cls, site2wpos):
-    #    return cls(site2wpos)
-
-    def __init__(self, site2wpos):
-        # TODO: Ordered dict.
-        d = {}
-        for k, s in site2wpos.items():
-            d[k] = WyckoffPositions.from_string(s)
-
-        self.site2wpos = d
-
-    def __str__(self):
-        """String representation."""
-        lines = []
-        app = lines.append
-        app("Wyckoff positions:")
-        for site, wpos in self.site2wpos.items():
-            app("%s site:\n\t%s" % (site, wpos))
-
-        return "\n".join(lines)
-
-    def generate_structure(self, lattice, site2params, to_unit_cell=False, struct_cls=None):
-        """
-        Generate structure object from `lattice` and dictionary `site2params`
-        """
-        if not isinstance(site2params, OrderedDict):
-            raise ValueError("Please use a OrderedDict for site2params so that\n" +
-                             "the algorithm used to build the structure is deterministic.")
-
-        # Build species and coords arrays.
-        species, coords = [], []
-        for elsym, wpos in self.site2wpos.items():
-            params = site2params[elsym]
-            _validate_params(params, wpos)
-            for vec_expr in wpos:
-                species.append(elsym)
-                # Evaluate sympy expression with params.
-                coords.append([float(expr.subs(params)) for expr in vec_expr])
-
-        #for sp, c in zip(species, coords):
-        #    print(sp, c, type(c), type(c[0]))
-
-        from abipy.core.structure import Structure
-        cls = struct_cls if struct_cls is not None else Structure
-        return cls(lattice, species, coords, validate_proximity=True,
-                   to_unit_cell=to_unit_cell, coords_are_cartesian=False, site_properties=None)
-
-    def find_params(self, structure):
-        """
-        Compute the value of the parameter given a Structure object.
-        """
-        frac_coords = structure.get_symbol2coords()
-
-        # Solve the problem.
-        params = {}
-        for elsym, wpos in self.site2wpos.items():
-            params[elsym] = wpos.params_from_frac_coords(frac_coords[elsym])
-
-        # Return results in a dict with symbol names.
-        return params
-
-    def error_of_params(self, params, structure):
-        frac_coords = structure.get_symbol2coords()
-
-        for elsym, wpos in self.site2wpos.items():
-            frac_errors = wpos.error_of_params(params[elsym], frac_coords[elsym])
-            print(frac_errors)
-
-
 class SiteSymmetries(Has_Structure):
 
     def __init__(self, structure):
@@ -265,13 +62,27 @@ class SiteSymmetries(Has_Structure):
             if len(rotations) != 0:
                 herm_symbol, ptg_num, trans_mat = spglib.get_pointgroup(rotations)
 
-            self.sitesym_labels.append("%s (#%d) ord: %d" % (herm_symbol.strip(), ptg_num, len(rotations)))
+            self.sitesym_labels.append("%s (#%d,nsym:%d)" % (herm_symbol.strip(), ptg_num, len(rotations)))
+
+        #for irred_isite in self.irred_isites:
+        #    for isite_eq, rm1, tau, l0 in self.eq_sites[irred_isite]:
+        #        # isite_eq = rm1(irred_site - tau) + l0
 
     @property
     def structure(self):
         """|Structure| object."""
         return self._structure
 
+    def __str__(self):
+        return self.to_string()
+
+    def to_string(self, verbose=0):
+        """String representation with verbosity level verbose."""
+        lines = []; app = lines.append
+        app(self.structure.to_string(verbose=verbose, title="Structure"))
+
+        return "\n".join(lines)
+
     def get_wyckoff_dataframe(self, view="all", select_symbols=None, decimals=5, verbose=0):
         """
         Find Wyckoff positions.
@@ -304,7 +115,13 @@ class SiteSymmetries(Has_Structure):
                 l0 = indsym[iatom, isym, :3]
                 l0 = sp.Matrix(l0)
                 m = rm1 * (vector - tau) - (l0 + vector)
-                if verbose: sp.pprint(m)
+                if verbose:
+                    print(92 * "=")
+                    print("System of linear equations for iatom %d, %s" % (iatom, repr(site)))
+                    sp.pprint(m)
+                    print("Rotation:")
+                    sp.pprint(rm1)
+
                 system.append(m)
 
             # Solve system of linear equations.
@@ -322,17 +139,17 @@ class SiteSymmetries(Has_Structure):
                 d["wyckoff"] = wlabel
                 d["site_symmetry"] = sitesym_labels[iatom]
                 for s in frac_symbols:
-                    d[s] = s
+                    d[str(s)] = str(s)
                 if len(solutions) > 1:
                     cprint("Found multiple solutions for iatom %d" % iatom, "red")
-                d.update(**solutions[0])
+                d.update({str(k): str(v) for k, v in solutions[0].items()})
                 rows.append(d)
 
         import pandas as pd
         df = pd.DataFrame(rows, index=None, columns=list(rows[0].keys()) if rows else None)
         return df
 
-    def get_tensor_rank2_dataframe(self, view="all", select_symbols=None, decimal=5, verbose=0):
+    def get_tensor_rank2_dataframe(self, view="all", select_symbols=None, decimals=5, verbose=0):
         """
 
         Args:
@@ -349,9 +166,11 @@ class SiteSymmetries(Has_Structure):
 
         sitesym_labels = self.structure.spget_site_symmetries()
 
-        # Symmetric tensor.
-        axx, ayy, azz, axy, axz, ayz = symbols = sp.symbols('axx, ayy, azz, axy, axz, ayz')
-        tensor = sp.Matrix(([axx, axy, axz], [axy, ayy, ayz], [axz, ayz, azz]))
+        # Symmetric tensor in reduced coordinates (direct lattice)
+        # Operations in reduced coords are given by integer matrices and this facilitates the solution
+        # of the system of equations with sympy.
+        Txx, Tyy, Tzz, Txy, Txz, Tyz = symbols = sp.symbols('Txx, Tyy, Tzz, Txy, Txz, Tyz')
+        tensor = sp.Matrix(([Txx, Txy, Txz], [Txy, Tyy, Tyz], [Txz, Tyz, Tzz]))
 
         indsym = self.structure.indsym
         rows = []
@@ -361,7 +180,14 @@ class SiteSymmetries(Has_Structure):
             for isym, rotf in enumerate(self.sp_symrel):
                 if indsym[iatom, isym, 3] != iatom: continue
                 m = rotf * tensor * rotf.T - tensor
-                if verbose: sp.pprint(m)
+                if verbose:
+                    print(92 * "=")
+                    print("System of linear equations for iatom %d, %s" % (iatom, repr(site)))
+                    sp.pprint(m)
+                    print("Rotation:")
+                    sp.pprint(rotf)
+                    print(92 * "=")
+
                 system.append(m)
 
             # Solve system of linear equations.
@@ -379,10 +205,10 @@ class SiteSymmetries(Has_Structure):
                 d["wyckoff"] = wlabel
                 d["site_symmetry"] = sitesym_labels[iatom]
                 for s in symbols:
-                    d[s] = s
+                    d[str(s)] = str(s)
                 if len(solutions) > 1:
                     cprint("Found multiple solutions for iatom %d" % iatom, "red")
-                d.update(**solutions[0])
+                d.update({str(k): str(v) for k, v in solutions[0].items()})
                 rows.append(d)
 
         import pandas as pd
diff --git a/abipy/core/structure.py b/abipy/core/structure.py
index 6308213e..1c75dcc7 100644
--- a/abipy/core/structure.py
+++ b/abipy/core/structure.py
@@ -604,8 +604,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
 
         return("\n".join(lines))
 
-    def get_conventional_standard_structure(self, international_monoclinic=True,
-                                            symprec=1e-3, angle_tolerance=5):
+    def get_conventional_standard_structure(self, international_monoclinic=True, symprec=1e-3, angle_tolerance=5):
         """
         Gives a structure with a conventional cell according to certain
 	standards. The standards are defined in :cite:`Setyawan2010`
@@ -638,6 +637,18 @@ class Structure(pymatgen.Structure, NotebookWriter):
 
         return self.__class__.as_structure(standardized_structure)
 
+    def refine(self, symprec=1e-3, angle_tolerance=5):
+        """
+        Get the refined structure based on detected symmetry. The refined
+        structure is a *conventional* cell setting with atoms moved to the
+        expected symmetry positions.
+
+        Returns: Refined structure.
+        """
+        sym_finder = SpacegroupAnalyzer(structure=self, symprec=symprec, angle_tolerance=angle_tolerance)
+        new = sym_finder.get_refined_structure()
+        return self.__class__.as_structure(new)
+
     def abi_sanitize(self, symprec=1e-3, angle_tolerance=5, primitive=True, primitive_standard=False):
         """
         Returns a new structure in which:
@@ -658,8 +669,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
 
         # Refine structure
         if symprec is not None and angle_tolerance is not None:
-            sym_finder = SpacegroupAnalyzer(structure=structure, symprec=symprec, angle_tolerance=angle_tolerance)
-            structure = sym_finder.get_refined_structure()
+            structure = structure.refine(symprec=symprec, angle_tolerance=angle_tolerance)
 
         # Convert to primitive structure.
         if primitive:
@@ -932,13 +942,14 @@ class Structure(pymatgen.Structure, NotebookWriter):
                               "Use `overwrite=True` to allow modification."))
 
         msg = ("Structure object does not have symmetry operations computed from Abinit.\n"
-               "Will call spglib to get symmetry operations.")
+               "Calling spglib to get symmetry operations.")
         cprint(msg, "magenta")
 
         spglib_data = SpacegroupAnalyzer(self).get_symmetry_dataset()
         spgid = spglib_data["number"]
         symrel, tnons = spglib_data["rotations"], spglib_data["translations"]
-        symafm = [1] * len(symrel)  # TODO: Anti-ferromagnetic symmetries are not supported by spglib
+        # TODO: Anti-ferromagnetic symmetries are not supported by spglib
+        symafm = [1] * len(symrel)
 
         abispg = AbinitSpaceGroup(spgid, symrel, tnons, symafm, has_timerev, inord="C")
         self.set_abi_spacegroup(abispg)
@@ -972,7 +983,7 @@ class Structure(pymatgen.Structure, NotebookWriter):
 
     @lazy_property
     def site_symmetries(self):
-        from abipy.core.wyckoff import SiteSymmetries
+        from abipy.core.site_symmetries import SiteSymmetries
         return SiteSymmetries(self)
 
     # TODO: site_symmetry or spget_site_symmetries?
@@ -984,13 +995,13 @@ class Structure(pymatgen.Structure, NotebookWriter):
         sitesym_labels = []
         for iatom, site in enumerate(self):
             rotations = [symrel[isym] for isym in range(nsym) if
-                         indsym[iatom, isym, 3] == iatom and symafm[isym] == 1]
+                         indsym[iatom, isym, 3] == iatom and symafm[isym] == +1]
             # Passing a 0-length rotations list to spglib can segfault.
             herm_symbol, ptg_num = "1", 1
             if len(rotations) != 0:
                 herm_symbol, ptg_num, trans_mat = spglib.get_pointgroup(rotations)
 
-            sitesym_labels.append("%s (#%d) nsym:%d" % (herm_symbol.strip(), ptg_num, len(rotations)))
+            sitesym_labels.append("%s (#%d,nsym:%d)" % (herm_symbol.strip(), ptg_num, len(rotations)))
 
         return sitesym_labels
 
@@ -1857,32 +1868,6 @@ class Structure(pymatgen.Structure, NotebookWriter):
 
         return "\n".join(lines)
 
-    #def ksampling_from_jhudb(self, **kwargs):
-    #    from pymatgen.ext.jhu import get_kpoints
-    #    __doc__ = get_kpoints.__doc__
-    #    kpoints = get_kpoints(self, **kwargs)
-    #    print(kpoints)
-    #    print(kpoints.style)
-    #    print("num_kpts", kpoints.num_kpts)
-
-    #    d = {"kptopt": 0,
-    #         "kpt": kpoints.kpts,
-    #         "nkpt": kpoints.num_kpts,
-    #         #"kptnrm": kptnrm,
-    #         "wtk": kpoints.kpts_weights,
-    #         "shiftk": kpoints.kpts_shift,
-    #         "chksymbreak": 0,
-    #    }
-    #    print(d)
-
-    #    #from pymatgen.io.abinit.abiobjects import KSampling
-    #    #return KSampling(mode=KSamplingModes.automatic,
-    #    #         num_kpts= 0,
-    #    #         kpts=((1, 1, 1),),
-    #    #         kpt_shifts=(0.5, 0.5, 0.5),
-    #    #         kpts_weights=None, use_symmetries=True, use_time_reversal=True, chksymbreak=None,
-    #    #         comment=None)
-
     def calc_ksampling(self, nksmall, symprec=0.01, angle_tolerance=5):
         """
         Return the k-point sampling from the number of divisions ``nksmall`` to be used for
diff --git a/abipy/core/tests/test_site_symmetries.py b/abipy/core/tests/test_site_symmetries.py
new file mode 100644
index 00000000..05d6743c
--- /dev/null
+++ b/abipy/core/tests/test_site_symmetries.py
@@ -0,0 +1,56 @@
+"""Tests for core.site_symmetries module"""
+from __future__ import print_function, division, unicode_literals, absolute_import
+
+import unittest
+import os
+import numpy as np
+
+from abipy.core.testing import AbipyTest
+from abipy.core.structure import Structure
+from abipy.core.site_symmetries import SiteSymmetries
+
+import abipy.data as abidata
+
+class TestSiteSymmetries(AbipyTest):
+    """Unit tests for SiteSymmetries."""
+
+    def test_si(self):
+        """Testing wyckoff positions for Si2"""
+        si = Structure.from_file(abidata.cif_file("si.cif"))
+        ss = si.site_symmetries
+        repr(ss); str(ss)
+        assert ss.to_string(verbose=2)
+        df = ss.get_wyckoff_dataframe(verbose=2)
+        self.assert_equal(np.array(df["xfrac"].values, dtype=float), [0, 0.25])
+        self.assert_equal(np.array(df["yfrac"].values, dtype=float), [0, 0.25])
+        self.assert_equal(np.array(df["zfrac"].values, dtype=float), [0, 0.25])
+        #0  -43m (#31) nsym:24                  0                  0                  0
+        #1  -43m (#31) nsym:24  0.250000000000000  0.250000000000000  0.250000000000000
+
+        df = ss.get_tensor_rank2_dataframe(verbose=2)
+        ref = ["Tzz", "Tzz"]
+        self.assert_equal(df["Txx"].values, ref)
+        self.assert_equal(df["Tyy"].values, ref)
+        ref = ["-Tzz/3", "-Tzz/3"]
+        self.assert_equal(df["Txy"].values, ref)
+        self.assert_equal(df["Txz"].values, ref)
+        self.assert_equal(df["Tyz"].values, ref)
+
+    def test_alpha_sio2(self):
+        """Testing wyckoff positions for alpha-SiO2"""
+        asi02 = Structure.from_file(os.path.join(abidata.dirpath, "refs", "mp-7000_DDB.bz2"))
+        ss = asi02.site_symmetries
+        df = ss.get_wyckoff_dataframe(verbose=2)
+        df[df["element"] == "Si"]
+        self.assert_equal(df["xfrac"].values, ["xfrac", "yfrac", "0.0"])
+        self.assert_equal(df["yfrac"].values, ["0.0", "yfrac", "yfrac"])
+        self.assert_equal(np.array(df["zfrac"].values, dtype=float), [0.833335, 0.5, 0.166665])
+
+        #df = ss.get_tensor_rank2_dataframe()
+        #ref = ["Tzz", "Tzz"]
+        #self.assert_equal(df["Txx"].values, ref)
+        #self.assert_equal(df["Tyy"].values, ref)
+        #ref = ["-Tzz/3", "-Tzz/3"]
+        #self.assert_equal(df["Txy"].values, ref)
+        #self.assert_equal(df["Txz"].values, ref)
+        #self.assert_equal(df["Tyz"].values, ref)
diff --git a/abipy/core/tests/test_wyckoff.py b/abipy/core/tests/test_wyckoff.py
deleted file mode 100644
index aa21e09a..00000000
--- a/abipy/core/tests/test_wyckoff.py
+++ /dev/null
@@ -1,61 +0,0 @@
-"""Tests for core.density module"""
-from __future__ import print_function, division, unicode_literals, absolute_import
-
-import unittest
-
-from collections import OrderedDict
-from abipy.core.testing import AbipyTest
-from abipy.core.structure import Lattice
-from abipy.core.wyckoff import Wyckoff
-
-
-class TestWyckoff(AbipyTest):
-    """Unit tests for Wyckoff."""
-
-    def test_sio2(self):
-        """Testing Wyckoff positions with SiO2 (requires sympy)"""
-        try:
-            import sympy
-        except ImportError:
-            raise unittest.SkipTest("sympy is not installed")
-
-        # http://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-wp-list?gnum=152
-        site2wpos = {}
-        #3	b	.2.
-        #site2wpos["Si"] = "(x,0,5/6)  (0,x,1/6)   (-x,-x,1/2)"
-        #3	a	.2.
-        site2wpos["Si"] = "[x,0,1/3], (0,x,2/3), (-x,-x,0)"
-        #6	c	1
-        site2wpos["O"] = """(x,y,z), (-y,x-y,z+1/3), (-x+y,-x,z+2/3), (y,x,-z),
-                            (x-y,-y,-z+2/3), (-x,-x+y,-z+1/3)"""
-
-        wyckoff = Wyckoff(site2wpos)
-        repr(wyckoff); str(wyckoff)
-
-        site2params = OrderedDict([
-            ("Si", dict(x=0.4763)),
-            ("O", dict(x=0.1588, y=0.7439, z=0.4612)),
-            #("O", dict(x=0.1588, y=0.7439, z=0.4612, yyy=3)),
-        ])
-
-        lattice = Lattice.hexagonal(a=4.971, c=5.473)
-        for to_unit_cell in [False]:
-        #for to_unit_cell in [False, True]:
-            structure = wyckoff.generate_structure(lattice, site2params, to_unit_cell=to_unit_cell)
-            print("to_unit_cell %s\n" % to_unit_cell, structure)
-
-            #from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
-            #spga = SpacegroupAnalyzer(structure)
-            #structure = spga.get_refined_structure()
-            #print("Refined:", structure)
-
-            wyck_params = wyckoff.find_params(structure)
-            #structure.perturb(distance=0.01)
-            #wyck_params = wyckoff.find_params(structure)
-
-            print("Error of params")
-            wyckoff.error_of_params(wyck_params, structure)
-
-            for site, params in site2params.items():
-                print(wyck_params[site])
-                assert wyck_params[site] == params
diff --git a/abipy/dfpt/msqdos.py b/abipy/dfpt/msqdos.py
index e95a31d0..43ad4454 100644
--- a/abipy/dfpt/msqdos.py
+++ b/abipy/dfpt/msqdos.py
@@ -46,11 +46,11 @@ class _Component(object):
 
 class MsqDos(Has_Structure):
     """
-    This object stores the generalized phonon DOS with the mean square displacement tensor in cartesian coords.
+    This object stores the generalized phonon DOS with the mean square displacement tensor in CARTESIAN coords.
     This DOS-like quantity allows one to calculate Debye Waller factors as a function of T
     by integration with 1/omega and the Bose-Einstein factor.
     See :cite:`Lee1995` for the further details about the internal implementation and
-    :cite:`Trueblood1996` for the different conventions used by crystallographers
+    :cite:`Trueblood1996` for the different conventions used by crystallographers.
     """
     C = _Component
     ALL_COMPS = OrderedDict([
@@ -73,7 +73,7 @@ class MsqDos(Has_Structure):
         Arg:
             structure: |Structure| object.
             wmesh: Frequency mesh
-            values: (natom, 3, 3, nomega) arrays with generalized DOS.
+            values: (natom, 3, 3, nomega) arrays with generalized DOS in cart coords.
             amu_symbol: Dictionary element.symbol -> mass in atomic units.
         """
         self._structure = structure
@@ -120,7 +120,7 @@ class MsqDos(Has_Structure):
         Compute mean square displacement for each atom in `iatom_list` as a function of T (bohr^2).
 
         Args:
-            tmesh: array-like with temperatures in Kelvin degrees
+            tmesh: array-like with temperatures in Kelvin.
             iatom_list: List of atom sites to comput. None for all.
             what_list:
         """
@@ -252,7 +252,7 @@ class MsqDos(Has_Structure):
 
     def write_cif_file(self, filepath, temp=300):
         """
-        Write CIF file with structure info and anisotropic U tensor in CIF format.
+        Write CIF file with structure and anisotropic U tensor in CIF format.
 
         Args:
             filepath: Name of CIF file. If None, a temporary filepath is created.
@@ -275,7 +275,7 @@ class MsqDos(Has_Structure):
         In the Vesta GUI, select: Properties -> Atoms -> Show as displament ellipsoids.
         """
         filepath = self.write_cif_file(filepath=None, temp=temp)
-        print("Writing structure + Debye-Waller tensor in CIF format for T = %s to file: %s" % (temp, filepath))
+        cprint("Writing structure + Debye-Waller tensor in CIF format for T = %s to file: %s" % (temp, filepath), "green")
         from abipy.iotools import Visualizer
         visu = Visualizer.from_name("vesta")
 
@@ -318,15 +318,16 @@ _atom_site_aniso_U_12""".splitlines()
         return s + "\n".join(aniso_u)
 
     def check_site_symmetries(self, temp=300, verbose=0):
-
+        """
+        Check symmetries of the displacement tensor at temperature `temp`.
+        Return maximum error.
+        """
         msq = self.get_msq_tmesh([float(temp)], what_list="displ")
-        values = getattr(msq, "displ")
+        values_cart = getattr(msq, "displ")
         natom = len(self.structure)
-        values = np.reshape(values, (natom, 3, 3))
+        values_cart = np.reshape(values_cart, (natom, 3, 3))
 
-        from abipy.core.wyckoff import SiteSymmetries
-        ss = SiteSymmetries(self.structure)
-        return ss.check_site_symmetries(values, verbose=verbose)
+        return self.structure.site_symmetries.check_site_symmetries(values_cart, verbose=verbose)
 
     def _get_components(self, components):
         """
diff --git a/abipy/dfpt/phonons.py b/abipy/dfpt/phonons.py
index 3754455e..de20910b 100644
--- a/abipy/dfpt/phonons.py
+++ b/abipy/dfpt/phonons.py
@@ -2833,7 +2833,7 @@ class PhdosFile(AbinitNcFile, Has_Structure, NotebookWriter):
             msqd_dos = None
         if msqd_dos is not None:
             yield msqd_dos.plot(units=units, show=False)
-            yield msqd_dos.plot_bfactors(show=False)
+            yield msqd_dos.plot_tensor(show=False)
 
     def write_notebook(self, nbpath=None):
         """
diff --git a/abipy/dfpt/phtk.py b/abipy/dfpt/phtk.py
index 2bfb18f3..bd40540b 100644
--- a/abipy/dfpt/phtk.py
+++ b/abipy/dfpt/phtk.py
@@ -13,7 +13,6 @@ from abipy.dfpt.phonons import PhononBands, get_dyn_mat_eigenvec, match_eigenvec
 from abipy.abio.inputs import AnaddbInput
 from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt
 from pymatgen.core.units import bohr_to_angstrom, eV_to_Ha
-from pymatgen.symmetry.bandstructure import HighSymmKpath
 
 
 class SoundVelocity(Has_Structure, NotebookWriter):
@@ -108,8 +107,7 @@ class SoundVelocity(Has_Structure, NotebookWriter):
             )
 
             if not directions:
-                hs = HighSymmKpath(ddb.structure, symprec=1e-2)
-
+                hs = ddb.structure.hsym_kpath
                 kpath = hs.kpath
 
                 directions = []
diff --git a/abipy/scripts/abistruct.py b/abipy/scripts/abistruct.py
index 594d85cc..5c8f636d 100755
--- a/abipy/scripts/abistruct.py
+++ b/abipy/scripts/abistruct.py
@@ -227,9 +227,11 @@ codes), a looser tolerance of 0.1 (the value used in Materials Project) is often
         help="Extract/Compute Abinit space group from file with structure.")
     p_abispg.add_argument("-t", "--tolsym", type=float, default=None, help="""\
 Gives the tolerance on the atomic positions (reduced coordinates), primitive vectors, or magnetization,
-to be considered equivalent, thanks to symmetry operations. This is used in the recognition of the set
+to be considered equivalent, thanks to symmetry operations. This value is used by ABINIT in the recognition of the set
 of symmetries of the system, or the application of the symmetry operations to generate from a reduced set of atoms,
-the full set of atoms. Note that a value larger than 0.01 is considered to be unacceptable.""")
+The internal default is 1e-8. Setting tolsym to a value larger than 1e-8 will make Abinit detect the spacegroup within
+this tolerance and re-symmetrize the input structure. This option is useful if the structure has been taken from a CIF
+file that does not have enough significant digits.""")
     p_abispg.add_argument("-d", "--diff-mode", type=str, default="table", choices=["table", "diff"],
         help="Select diff output format.")
 
@@ -283,12 +285,16 @@ Has to be all integers. Several options are possible:
     p_proto.add_argument("--angle-tol", default=5, type=float, help="angle tolerance.")
 
     # Subparser for wyckoff.
-    p_wyckoff = subparsers.add_parser('wyckoff', parents=[copts_parser, path_selector],
+    p_wyckoff = subparsers.add_parser('wyckoff', parents=[copts_parser, spgopt_parser, path_selector],
             help="Print wyckoff positions. WARNING: still under development!")
+    p_wyckoff.add_argument("--refine", default=False, action="store_true",
+                            help="Use spglib to refine structure before computation")
 
     # Subparser for tensor_site.
-    p_tensor_site = subparsers.add_parser('tensor_site', parents=[copts_parser, path_selector],
+    p_tensor_site = subparsers.add_parser('tensor_site', parents=[copts_parser, spgopt_parser, path_selector],
             help="Print symmetry properties of tensors due to site-symmetries. WARNING: still under development!")
+    p_tensor_site.add_argument("--refine", default=False, action="store_true",
+                                help="Use spglib to refine structure before computation")
 
     # Subparser for neighbors.
     p_neighbors = subparsers.add_parser('neighbors', parents=[copts_parser, path_selector],
@@ -498,20 +504,25 @@ def main():
     elif options.command == "abispg":
         structure = abilab.Structure.from_file(options.filepath)
         check_ordered_structure(structure)
-        spgrp = structure.abi_spacegroup
+        abi_spg = structure.abi_spacegroup
 
-        if spgrp is not None:
+        if abi_spg is not None and options.tolsym is None:
             print(structure.spget_summary(verbose=options.verbose))
         else:
-            # Here we compare Abinit wrt spglib. If spgrp is None, we create a temporary
+            # Here we compare Abinit wrt spglib. If abi_spg is None, we create a temporary
             # task to run the code in dry-run mode.
-            print("FILE does not contain Abinit symmetry operations.")
-            print("Calling Abinit in --dry-run mode with chkprim = 0 to get space group.")
+            if abi_spg is None:
+                print("FILE does not contain Abinit symmetry operations.")
+            cprint("Calling Abinit in --dry-run mode with chkprim = 0 to get space group.")
+            if options.tolsym is not None and options.tolsym > 1e-8:
+                cprint("Crystal structure will be re-symmetrized by Abinit with tolsym: %s" % options.tolsym, "yellow")
+
             from abipy.data.hgh_pseudos import HGH_TABLE
             gsinp = factories.gs_input(structure, HGH_TABLE, spin_mode="unpolarized")
             gsinp["chkprim"] = 0
             abistructure = gsinp.abiget_spacegroup(tolsym=options.tolsym)
             print(abistructure.spget_summary(verbose=options.verbose))
+            print("")
 
             diff_structures([structure, abistructure], mode=options.diff_mode,
                             headers=["Input structure", "After Abinit symmetrization"], fmt="abivars")
@@ -659,17 +670,23 @@ def main():
 
     elif options.command == "wyckoff":
         structure = abilab.Structure.from_file(options.filepath)
-        from abipy.core.wyckoff import SiteSymmetries
-        ss = SiteSymmetries(structure)
-        df = ss.get_wyckoff_dataframe(view="all", select_symbols=None, verbose=options.verbose)
-        abilab.print_dataframe(df)
+        if options.refine:
+            print("Refining structure with symprec: %s, angle_tolerance: %s" % (options.symprec, options.angle_tolerance))
+            structure = structure.refine(symprec=options.symprec, angle_tolerance=options.angle_tolerance)
+        print(structure.spget_summary(verbose=options.verbose))
+        ss = structure.site_symmetries
+        df = ss.get_wyckoff_dataframe(verbose=options.verbose)
+        abilab.print_dataframe(df, title="\nWyckoff positions in reduced coordinates.")
 
     elif options.command == "tensor_site":
         structure = abilab.Structure.from_file(options.filepath)
-        from abipy.core.wyckoff import SiteSymmetries
-        ss = SiteSymmetries(structure)
-        df = ss.get_tensor_rank2_dataframe(view="all", select_symbols=None, verbose=options.verbose)
-        abilab.print_dataframe(df)
+        if options.refine:
+            print("Refining structure with symprec: %s, angle_tolerance: %s" % (options.symprec, options.angle_tolerance))
+            structure = structure.refine(symprec=options.symprec, angle_tolerance=options.angle_tolerance)
+        print(structure.spget_summary(verbose=options.verbose))
+        ss = structure.site_symmetries
+        df = ss.get_tensor_rank2_dataframe(verbose=options.verbose)
+        abilab.print_dataframe(df, title="\nTensor components in reduced coordinates (rank 2, symmetric)")
 
     elif options.command == "neighbors":
         abilab.Structure.from_file(options.filepath).print_neighbors(radius=options.radius)
@@ -757,8 +774,11 @@ def main():
 
     #elif options.command == "kmesh_jhu":
     #    structure = abilab.Structure.from_file(options.filepath)
-    #    ksampling = structure.ksampling_from_jhudb(kppra=1000)
-    #    #print(ksampling)
+    #    from pymatgen.ext.jhu import get_kpoints
+    #    kpoints = get_kpoints(structure, min_distance=0, min_total_kpoints=1,
+    #                           kppra=None, gap_distance=7, remove_symmetry=None,
+    #                           include_gamma="auto", header="simple", incar=None)
+    #    #print(kpoints)
 
     elif options.command == "lgk":
         structure = abilab.Structure.from_file(options.filepath)
diff --git a/docs/api/core_api.rst b/docs/api/core_api.rst
index d55fefac..c64ae4f5 100644
--- a/docs/api/core_api.rst
+++ b/docs/api/core_api.rst
@@ -124,10 +124,10 @@ core Package
    :undoc-members:
    :show-inheritance:
 
-:mod:`wyckoff` Module
----------------------
+:mod:`site_symmetries` Module
+-----------------------------
 
-.. automodule:: abipy.core.wyckoff
+.. automodule:: abipy.core.site_symmetries
    :members:
    :undoc-members:
    :show-inheritance:
diff --git a/docs/installation.rst b/docs/installation.rst
index aa05c980..707f5efe 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -72,9 +72,9 @@ By default, the installer creates the ``anaconda`` directory in your home.
 Anaconda will add one line to your ``.bashrc`` to enable access to the anaconda executables.
 Once the installation is completed, execute::
 
-    source ~/anaconda/bin/activate root
+    source ~/anaconda/bin/activate base
 
-to activate the ``root`` environment.
+to activate the ``base`` environment.
 The output of ``which python`` should show that you are using the python interpreter provided by anaconda.
 
 Use the conda_ command-line interface to install the packages not included in the official distribution.
@@ -86,17 +86,17 @@ Remember that if a package is not available in the official conda repository, yo
 download the package from one of the conda channels or use ``pip install`` if no conda package is available.
 
 Fortunately there are conda channels providing all dependencies needed by AbiPy.
-To install the pymatgen_ package from the matsci_ channel, use::
+Now add ``conda-forge``, ``matsci`` and ``abinit`` to your conda channels with::
 
-    conda install pymatgen --channel matsci
+    conda config --add channels conda-forge
+    conda config --add channels matsci
+    conda config --add channels abinit
 
-then install Abipy from the abinit-channel_ with::
+These are the channels from which we will download pymatgen, abipy and abinit.
+Finally, install AbiPy from the abinit-channel_ with::
 
     conda install abipy --channel abinit
 
-Visit `materials.sh <http://materials.sh>`_ for instructions on how to use the
-matsci channel to install pymatgen and other packages.
-
 Once you have completed the installation of AbiPy and pymatgen, open the ipython_ shell and type::
 
     # make sure spglib library works
diff --git a/docs/postprocessing_howto.rst b/docs/postprocessing_howto.rst
index 3058d247..1255a64c 100644
--- a/docs/postprocessing_howto.rst
+++ b/docs/postprocessing_howto.rst
@@ -64,7 +64,7 @@ Use::
 to read the structure from ``FILE`` and generate a CIF_ file (default behaviour).
 
 The majority of the netcdf_ files produced by Abinit contain structural information
-so this command can be used with netcdf ouput files as well as Abinit input/output 
+so this command can be used with netcdf output files as well as Abinit input/output 
 files and all the other formats supported by pymatgen e.g. POSCAR files.
 Other formats can be specified with the ``-f`` option.
 For example::
@@ -101,7 +101,7 @@ to validate the input file with Abinit (requires ``manager.yml`` and, obviously,
 The script provides other options to invoke Abinit 
 to get space group information, the list of k-points in the IBZ.
 the list of atomic perturbations for phonons or the list of autoparal configurations.
-See ``abinp.py --help`` for futher info.
+See ``abinp.py --help`` for further info.
 
 Print the warnings in the log file
 ----------------------------------
@@ -232,6 +232,24 @@ to generate a template with the input variables defining the k-path
         +0.62500  +0.25000  +0.62500 # U
         +0.50000  +0.00000  +0.50000 # X
 
+
+Re-symmetrize a structure when Abinit reports less symmetries than expected
+---------------------------------------------------------------------------
+
+Crystalline structures saved in text format (e.g. CIF files downloaded from
+the Materials Project websites) may not have enough significant digits
+and Abinit may not find the same spacegroup as the one reported by the source
+as the default tolerance for symmetry detection in Abinit is tight (tolsym = 1e-8).
+
+In this case, one can use the `abispg` option of abistruct.py to compute the spacegroup
+with Abinit and a tolerance larger that the default value::
+
+    abistruct.py abispg problematic.cif --tolsym=1e-6
+
+Hopefully, the code will detect the correct spacegroup, will re-symmetrize 
+the initial lattice vectors and atomic positions and print the new symmetrized structure to terminal.
+
+
 Get neighbors for each atom in the unit cell out to a distance radius
 ---------------------------------------------------------------------
 
@@ -402,7 +420,7 @@ Use::
 
     abiview.py phbands out_PHBST.nc -web
 
-to start a local webserver and open the HTML page inside the default browser 
+to start a local web server and open the HTML page inside the default browser 
 (the browser can be changed with the ``--browser`` option).
 
 It is also possible to visualize the phonon modes starting directly from a DDB_ file with::
@@ -410,7 +428,7 @@ It is also possible to visualize the phonon modes starting directly from a DDB_
     abiview.py ddb -web
 
 In this case, AbiPy will invoke anaddb to produce the ``PHBST.nc`` file on an automatically 
-generated q-path and then start the webserver.
+generated q-path and then start the web server.
 
 Visualize the results of a structural relaxation
 ------------------------------------------------
@@ -461,7 +479,7 @@ the AbiPy scripts are quite handy for a quick analysis of the results.
 Compare multiple files
 ----------------------
 
-The :ref:`abicomp.py` script is explicitely designed for this kind of task.
+The :ref:`abicomp.py` script is explicitly designed for this kind of task.
 It operates on multiple files (usually files with the same extension) and 
 either produces matplotlib_ plots or creates AbiPy robots providing methods
 to analyze the results, perform convergence studies and build pandas DataFrames_.