Compare initial and final slab configuration 
        to determine if slab reconstructs during relaxation

        xy_cutoff: Allowed xy-movement is determined from 
                   the covalent radii as:
                   xy_cutoff * np.mean(cradii)
        z_cutoff:  Allowed z-movement is determined as
                   z_cutoff * cradii_i

def is_reconstructed(self, xy_cutoff=0.3, z_cutoff=0.4):
        """Compare initial and final slab configuration 
        to determine if slab reconstructs during relaxation

        xy_cutoff: Allowed xy-movement is determined from 
                   the covalent radii as:
                   xy_cutoff * np.mean(cradii)
        z_cutoff:  Allowed z-movement is determined as
                   z_cutoff * cradii_i
        """

        assert self.A, \
            'Initial slab geometry needed to classify reconstruction'

        # remove adsorbate
        A = self.A[:-1].copy()
        B = self.B[:-1].copy()

        # Order wrt x-positions
        x_indices = np.argsort(A.positions[:, 0])
        A = A[x_indices]
        B = B[x_indices]
        a = A.positions
        b = B.positions

        allowed_z_movement = z_cutoff * cradii[A.get_atomic_numbers()]
        allowed_xy_movement = \
            xy_cutoff * np.mean(cradii[A.get_atomic_numbers()])

        D, D_len = get_distances(p1=a, p2=b, cell=A.cell, pbc=True)
        d_xy = np.linalg.norm(np.diagonal(D)[:2], axis=0)
        d_z = np.diagonal(D)[2:][0]

        cond1 = np.all(d_xy < allowed_xy_movement)
        cond2 = np.all([d_z[i] < allowed_z_movement[i]
                        for i in range(len(a))])

        if cond1 and cond2:  # not reconstructed
            return False
        else:
            return True