#  Copyright (C) 2017 The authors of Py-ChemShell
#
#  This file is part of Py-ChemShell.
#
#  Py-ChemShell is free software: you can redistribute it and/or modify
#  it under the terms of the GNU Lesser General Public License as
#  published by the Free Software Foundation, either version 3 of the
#  License, or (at your option) any later version.
#
#  Py-ChemShell is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with Py-ChemShell.  If not, see
#  <http://www.gnu.org/licenses/>.

import chemsh
from numpy            import zeros, full, float64, int64, dtype, count_nonzero, ma
from sys              import stdout
from .fragbase        import _FragmentBase
from .bqs             import BQs
from .cell            import Cell
from .xyz             import _XYZ
from .shells          import Shells
from ..data           import atomicdata, unitconvert
from ..dl_py2f        import utils
from ..base.constants import *
from ..io.file        import dump

class Fragment(_FragmentBase):
    '''Chemical fragment class'''

    # default size (natoms) is 0, but the shapes (of charges, conn, coords, etc.) need to be (1,...) to hold the default values
    _kwargs = {
                'conn'         : full((1,MAXCONN) , -1, dtype=int64),
                'frag_name'    : '',
# FIXME: bool is only 8-bit that cannot be aligned in Fortran 2D-array
#        so frozen will be passed as int64 for now
# YL: we may not need this flag (can treat this list as for qm_region)
                'frozen'       : zeros(shape=(1,) ,     dtype=int64),
                'grid'         : None,
                'map'          : zeros(shape=(1,),      dtype=int64),
                'masses'       : zeros(shape=(1,) ,     dtype=float64),
                'natoms'       : 0,
# FIXME: nshells and nbqs are set here to avoid seg fault in cluster routines
                'nbqs'         : 0,
                'nshells'      : 0,
                'regions'      : full((1)         , '0', dtype='S8'),
                'segIDs'       : zeros(shape=(1,) ,     dtype=int64),
                # YL 15/04/2018: shellatoms in Fragment is a list of shells indices and coreatoms in Shells is one of core indices
                'shellatoms'   : full((1,)        , -1, dtype=int64),
                'shells'       : Shells(),
                'spacegroup'   : 1,
                'totalcharge'  : 0.0,
                # YL 03/11/2021: added types for _FragmentBase for we need them for assigning forcefield
                'types'        : full((1)         , '', dtype='S8'),
                'znums'        : zeros(shape=(1,),      dtype=int64),
              }

    # fields of <nputils.RecArray _master>: the order matters!
    _fields = 'names', 'types', 'regions', 'znums', 'map', 'coords', 'charges', 'conn', 'fragIDs', 'frozen', 'masses', 'shellatoms'

    # make sure 'unit' is processed before 'coords' and
    #           'shellatoms' before 'coords' (connectivity depends on shell charges!)
    # YL 02/07/2020: attach more items to make sure they are dealth with in a correct order when the Python version is <3.6 (for example on ARCHER)
    # YL 13/09/2020: dependencies are cell -> cell.radius -> frag.coords
    _priorities = 'connect_mode', 'coords', 'cell', 'shellatoms', 'names', 'znums', 'charges'

    _internals = {
                 }

# FIXME: move to _synons
    synons = { 'natoms':'npts',
               'coords':'xyz',
             }

    _synons = {
                'coords':'xyz',
                'natoms':'npts',
                'q'     :'charges'
              }


    def __new__(cls, *args, **kwargs):
        ''''''

        from os.path   import splitext
        from ctypes    import Structure
        from .fragment import Fragment
        from .fragbio  import FragmentBio

        selectcases = { '.pdb': FragmentBio,
                        '.pqr': FragmentBio,
                        '.psf': FragmentBio,
                      }

        extension = ''
        # distinguish file types to return specific Fragment class
        c = kwargs.get('coords', None)
        if type(c) is str:
            filename, extension = splitext(kwargs['coords'])

        # default class: Fragment
        inst = Structure.__new__(selectcases.get(extension.lower(), Fragment), *args, **kwargs)

        inst.__init__(*args, **kwargs)

        return inst


# FRAGMENT PROPERTIES #
# setters decorated by @nputils.resize will allow users to resize through input data


    @property
    def frag_name(self):
        return self.__frag_name


    @frag_name.setter
    def frag_name(self, frag_nam):

        self.__frag_name = frag_nam


    @property
    def types(self):
        '''Forcefield types'''

        return self._master.types


    @types.setter
    @_XYZ._resize1D
    def types(self, val):
        '''Setter of types'''

        self._master.setField(val, field='types')


    @property
    def regions(self):
        ''''''

        return self._master.regions


    @property
    def nregions(self):
        '''Number of regions'''

        from numpy import char, unique

        return unique(char.strip(self._master.regions)).size


    @property
    def frozen(self):
        '''Frozen or not'''

        return self._master.frozen


    @frozen.setter
    def frozen(self, val):
        '''Setter of frozen'''

        self._master.setField(val, field='frozen')


    @property
    def species(self):
        '''Atom species'''

        from numpy                   import unique, __version__
        # TD 18/06/2024 : Numpy v2.0.0 version dependence
        from numpy.lib import NumpyVersion
        if NumpyVersion(__version__) >= '2.0.0b1':
            from numpy._core.defchararray import title
        else:
            from numpy.core.defchararray import title

        return unique(title(self.symbols.astype(str)))


    @property
    def nspecies(self):
        '''Number of species'''

        return self.species.size


    @property
    def map(self):
        '''Mapping from periodic to fragment'''

        return self._master.map

 
    @map.setter
    @_XYZ._resize1D
    def map(self,val):
        ''''setter of map'''

        self._master.setField(val, field='map')


    @property
    def segIDs(self):
        '''Array of nonbonded segment IDs'''

        try:
            return self.__segIDs
        except:
            pass

        from numpy import any, argmax, full, hstack
        from sys   import getrecursionlimit, setrecursionlimit

        self.__segIDs = full((self.natoms,), -1, dtype=int)

        # YL 07/10/2019: Python has a very low limit of max recursion depth (~1000). but changing it can be dangeous
        setrecursionlimit(100000)

        try:
            bonds = self.bonds
        except:
            return self.__segIDs

        def findChain(_i, residue=0):

            # reckon with all I--J and J--I cases
            js = hstack([bonds[bonds[:,0]==_i][:,1], bonds[bonds[:,1]==_i][:,0]])
            # only calculate undetermined J atoms
            for j in js[self.__segIDs[js]==-1]:
                # assign for current J
                self.__segIDs[j] = residue
                # recurse
                findChain(j, residue=residue)

        segID = 0
        # loop until all atoms determined
        while any(self.__segIDs==-1):

             # find the first undetermined atom I
             i = argmax(self.__segIDs==-1)
             self.__segIDs[i] = segID
             # find the whole chain lead by I
             findChain(i, residue=segID)
             segID += 1

        return self.__segIDs


    # YL 25/11/2021: might feebly improve the performance of _connectCovalentIIs() and _connectCovalentIJs()
    #                by <10%
    # YL NB 29/11/2021: this should NOT be inherited from <class FragmentBase> but must be defined
    #                   independently (same in <class FragmentBio>) otherwise self.__covalent_radii
    #                   becomse self._FragmentBase_covalent__radii but not self._Fragment_covalent__radii
    @property
    def _covalent_radii(self):
        '''Return the covalent radii'''

        try:
            return self.__covalent_radii
        except:
            self.__covalent_radii = atomicdata.covalent_radii[self.znums]
            return self.__covalent_radii


# END OF FRAGMENT PROPERTIES #
# GRID FOR CHARGE FITTING


    @property
    def grid(self):
        '''Grid for charges fitting'''

        return self._grid


    @grid.setter
    def grid(self, val):
        '''Setter of grid'''

        self._grid = val


# END OF GRID FOR CHARGE FITTING
# BKGQ PROPERTIES #


    @property
    def ngridpoints(self):
        '''Number of grid points'''

        if self._grid:
            return self._grid._master.shape[0]
        else:
            return 0


    @ngridpoints.setter
    def ngridpoints(self, val):
        '''Setter of ngridpoints'''

        self._grid.npoints = val


    @property
    def totalcharge(self):
        '''Total charge'''

        return self._totalcharge


    @totalcharge.setter
    def totalcharge(self, val):
        '''Setter of totalcharge'''

        self._totalcharge = float(val)


    # YL 28/08/2019: rewritten
    def addCharges(self, dictObj):
        """
        Add charge values to the atoms in the Fragment object
            Parameters:
                dictObj (Dictionary)
                    Dictionary of type { label: charge } to that values can be appropriately added
        """

        for key, val in dictObj.items():
            selected = self.names==key.encode()
            self.charges[selected] = val

# END OF BKGQ PROPERTIES #
# SHELL PROPERTIES #

    # depends on self.shells
    @property
    def nshells(self):
        '''Number of shells'''

        from numpy import arange

        # get a "private" handle
        self._shells._Shells__corecoords = self.coords

        # AJL, Aug2017: Added this to fix the scenario where you load a fragment but shellatoms
        # hasn't been updated yet with the information from the shells module.
        # e.g. mgo = Fragment(coords='mgo_regions.pun'); print(mgo.nshells)
        # will return 0 even though the shells are in loaded Fragment!
        # However, do you even need self.shellatoms? Couldn't you just return the length of self_shells.coreatoms 
        # from this routine, and return the array itself if someone requests self.shellatoms?
        # This would save on this additional need for book-keeping.
        # YL 17/04/2018: refactered so that Shell.coreatoms is consistent with Fragment.shellatoms
        if self._shells:
            coreatoms = self._shells.coreatoms
            if count_nonzero(coreatoms >= 0):
                self.shellatoms[coreatoms] = arange(self._shells.natoms)

#        # YL 29/03/2018: coreatoms entry may be all -1 in a punch file
        return self._shells.natoms
        # YL FIXME: len(self.shellatoms) may not be up-to-date
#        return len(self.shellatoms)


    @property
    def shells(self):
        '''Wrapped shells data'''

        from numpy import add, arange, count_nonzero, subtract

        # get a "private" handle
        self._shells._Shells__corecoords = self.coords

        # YL 17/04/2018: refactered so that Shell.coreatoms is consistent with Fragment.shellatoms
        coreatoms  = self._shells.coreatoms
        shellatoms = self.shellatoms
        if count_nonzero(coreatoms >= 0):
            self.shellatoms[coreatoms] = arange(self._shells.natoms)

        # check size
        # should NOT use expand(), for example: self._shells._master.expand(self.nshells)
        if self.nshells:
            self._shells.natoms = self.nshells
        else:
            return self._shells

        # see: Shells.coords
        if count_nonzero(coreatoms >= 0) and count_nonzero(self.coords):
            # if displacements are not just default and frag.coords are not just zeros: calculate new positions
            if count_nonzero(self._shells.displace-self._shells._Shells__displace):
                self._shells._master.setField(add(self._shells.displace, self.coords.take(coreatoms, axis=0)), field='xyz')
            # if displacements are just default: keep absolute positions and calculate new displacements
            elif count_nonzero(self._shells.coords):
                self._shells._master.setField(subtract(self._shells.coords, self.coords.take(coreatoms, axis=0)), field='displace')
            # if frag.shells.coords == displacements (in case shells are not read from a file but defined with a list of coreatoms)
            else:
                self._shells._master.setField(self.coords.take(coreatoms, axis=0), field='xyz')

            # names
            self._shells.names = self.names.take(coreatoms)

        return self._shells


    @shells.setter
    def shells(self, val):
        '''Setter of shells'''

        self._shells = val


    @property
    def effectcharges(self):
        '''Total core and shell charge'''

        from numpy import copy as npcopy

        self._effectcharges = npcopy(self.charges)

        # YL 21/10/2017: added a `nshells` check to prevent error
        if self.nshells:
            self._effectcharges[self._shells.coreatoms] += self.shells.charges

        return self._effectcharges


    @property
    def shellatoms(self):
        '''List of shell atom indices'''

        from numpy import all

        # YL 17/08/2019: if not yet initialised
        if all(self._master.shellatoms==-1) and not all(self._shells.coreatoms==-1):
            self._master.shellatoms[self._shells.coreatoms] = range(self._shells.coreatoms.size)

        return self._master.shellatoms


    @shellatoms.setter
    def shellatoms(self, val):
        '''Setter of shellatoms'''

        from numpy import argwhere, count_nonzero

        # YL 16/04/2018: now Fragment.shellatoms is a list of indices of the associated shell centres
        self._master.setField(val, field='shellatoms')
        self._shells.natoms = count_nonzero(self._master.shellatoms+1)
        self._shells.coreatoms[:] = argwhere(self._master.shellatoms+1).ravel()


# END OF SHELL PROPERTIES #
# CLASSMETHODS #





# END OF CLASSMETHODS #
# CLUSTER CUTTING METHODS #

# AJL, 13/06/18. Duplicate of adjustCharge
#    def cluster_adjust_charge(self, *args, **kwargs):
#        ''''''
#
#        from .. import cluster
#
#        return cluster.adjust_charge(self, *args, **kwargs)

    def adjustCharge(self, periodicFrag, *args, **kwargs):
        '''Adjust the charge of the cluster'''

        from .. import cluster

        return cluster.adjustCharge(periodicFrag, self, *args, **kwargs)

    def calcElectrostatics(self, periodicFrag, *args, **kwargs):
        '''Calculate the periodic electrostatics'''

        from .. import cluster

        return cluster.calcElectrostatics(periodicFrag, self, *args, **kwargs)

# AJL, 13/06/2018. Edited to duplicate behaviour of adjustCharge
# For consistency, should this then be renamed just fitBQs?
    def cluster_fit_bqs(self, fragObj, *args, **kwargs):
        '''Fit charges on BQs'''

        from .. import cluster

# AJL, 13/06/18. Call to non-existent function
#        return cluster.fit_bqs(self, *args, **kwargs)

        return cluster.fitBQs(fragObj, self, *args, **kwargs)


    def construct_cluster(self, *args, **kwargs):
        '''Construct the cluster'''

        from .. import cluster

        return cluster.construct(self, *args, **kwargs)


    def cut_cluster(self, *args, **kwargs):
        '''Cut the cluster'''

        from .. import cluster

        clusterObj = Fragment(connect_mode=self.connect_mode)
        cluster.cut(self, clusterObj, *args, **kwargs)
        return clusterObj


    def partition(self, *args, **kwargs):
        '''Divide a cluster according to the 5 region model'''

        from ..cluster import partition

        return partition.partition(self, *args, **kwargs)


    def makeFF(self, *args, **kwargs):
        '''Generate force field data (with DL_FIELD)'''

        from .. import dl_field 

        return dl_field.run(self, *args, **kwargs)


# END OF CLUSTER CUTTING METHODS #
# INQUIRIES #


    def getClosest(self, centre, nclosests=1, species=None):
        '''Index/indices of the closest atom or nclosests atoms'''

        from numpy  import ones, array, argsort, char, linalg, ndarray, subtract
        from ..base import forcefield

        if array(centre).shape == ():
            vects = subtract(self.coords, self.coords[centre])
        elif array(centre).shape == (3,):
            vects = subtract(self.coords, centre)

        # sorted array of closest sites
        closests = argsort(argsort(linalg.norm(vects, axis=1)))

        # filter
        if species:
            masks = self.select(types=species, return_masks=True)[closests]
        else:
            masks = ones((self.natoms,), dtype=bool)

        if closests.size == 0:
            return 0
        elif closests[masks][1:][:nclosests].size == 1:
            return closests[masks][1:][:nclosests][0]
        elif closests[masks][1:][:nclosests].size > 1:
            return closests[masks][1:][:nclosests]


    def angle(self, iatom, jatom, katom):
        '''Angle I-J-K in degrees'''

        from numpy import arccos, clip, dot, rad2deg

        return rad2deg(arccos(dot(self.diff(iatom,jatom),
                                  self.diff(katom,jatom))
                                 /self.dist(iatom,jatom)
                                 /self.dist(jatom,katom)))


    def dihedral(self, iatom, jatom, katom, latom):
        '''Torsion angle I-J-K-L in degrees'''

        from numpy import arctan2, clip, cross, dot, rad2deg

        # normalised vectors
        vec_ij = self.unitvec(iatom,jatom)
        vec_jk = self.unitvec(jatom,katom)
        vec_kl = self.unitvec(katom,latom)

        plane_ijk = cross(vec_ij, vec_jk)
        plane_jkl = cross(vec_jk, vec_kl)
        orthonorm = cross(plane_ijk, vec_jk)

        x = dot(plane_ijk, plane_jkl)
        y = dot(orthonorm, plane_jkl)

        # use arctan2 to get right signs
        return rad2deg(arctan2(y, x))


    # AJL, Aug2017: Check for charges as requested in io/punch.py
    def hasCharges(self):
        '''Check if the cluster has charge on the atoms'''

        from numpy import count_nonzero

        if count_nonzero(self.charges):
            return True
        else:
            return False


    @property
    def bonds(self):
        '''Generate a list containing all bonds'''

        from numpy  import flatnonzero, trunc, vstack, copy, unique
        from ..base import constants

        try:
            return self.__bonds
        except:
            pass

        # get a handle
        conn = self.conn
        MAXCONN = constants.MAXCONN

        iatoms = trunc(flatnonzero(conn+1)/MAXCONN).astype(int, copy=False)
        jatoms = conn[conn>=0]

        # make a copy and sort in place
        ij = vstack((iatoms, jatoms)).T
        sorted = copy(ij)
        sorted.sort(axis=1)

        # get the indices for reconstructing ij
        # YL NB: axis not available utill NumPy v1.13!
        non_repeateds = unique(sorted, axis=0, return_index=True)[1]

        # reconstruct
        # YL NB: axis not available utill NumPy v1.13!
        self.__bonds = unique(ij[non_repeateds], axis=0)

        return self.__bonds


    @property
    def nbonds(self):
        '''Number of bonds'''

        return self.bonds.shape[0]


    @property
    def alldihedrals(self):
        '''Generate a list containing all dihedrals'''

        from numpy     import any, all, asarray, array, extract, ceil, copy, fix, meshgrid, inf, intersect1d, fromiter, nditer, ndarray, where
        from numpy     import count_nonzero, argmin, flatnonzero, remainder, modf, trunc, int16, hstack, vstack, isin, set_printoptions, equal, __version__
        # TD 18/06/2024 : Numpy v2.0.0 version dependence
        from numpy.lib import NumpyVersion
        if NumpyVersion(__version__) >= '2.0.0b1':
            from numpy._core import defchararray
        else:
            from numpy.core import defchararray
        from itertools import permutations
        from ..base    import constants

        from time import time

        ss = time()
        try:
            return self.__dihedrals
        except:
            pass

        maxconn = constants.MAXCONN

        allatoms = range(self.natoms)

        conn = self.conn

        # an array of all permutations of 4 atoms
        klmnatoms = fromiter(permutations(allatoms, 4), dtype=('i,i,i,i'))
        klmnatoms = klmnatoms.view('i').reshape((klmnatoms.size,4)).astype(str, copy=False)
#        klmnatoms = nditer(permutations(allatoms, 4), ['refs_ok'])

        # all bonds (converted to str joined by '-')
        iatoms = trunc(flatnonzero(conn+1)/maxconn).astype(int, copy=False).astype(str, copy=False)
        jatoms = conn[conn>=0].astype(str, copy=False)
        allbonds  = defchararray.add(defchararray.add(iatoms, '-'), jatoms)

        # check if the bonds in klmnatoms contain `allbonds`
        for i in range(3):
            alldiheds = defchararray.add(klmnatoms[:,i], '-')
            alldiheds = defchararray.add(alldiheds, klmnatoms[:,i+1])
            bonds, diheds = meshgrid(allbonds, alldiheds, indexing='ij')[:]
            bondsindiheds = defchararray.find(diheds, bonds)
            klmnatoms = klmnatoms[any(bondsindiheds+1, axis=0)]

        # remove repeated dihedrals
        for i in reversed(range(klmnatoms.shape[0])):
            if any(all(klmnatoms[i][::-1] == klmnatoms[:], axis=1)):
                klmnatoms[i] = '-1'

        dihedrals = klmnatoms.astype(int, copy=False)
        self.__dihedrals = dihedrals[all(dihedrals+1, axis=1)]

        return self.__dihedrals


# END OF INQUIRIES #
# FACILITIES #

    # AL, 08/2019: Doesn't need the kwargs - we don't do anything with the information
    def addShells(self, *args, displace=Shells._Shells__displace, factor=1.5, shift=1.0, charges=None): #**kwargs):
        '''
           If charges is defined as a dictionary, this will be used to set charges.
           
           If charges is undefined, factor and shift will be used to calculate shell charges:

           Core       :  1.1 ->  1.1*factor - shift
           Added shell:          1.1*(1.0-factor) + shift
           Core       : -1.1 ->  1.1*factor + shift
           Added shell:         -1.1*(1.0+factor) - shift
        '''

        from numpy import arange

        # selected atoms
        selected = self.selectByAtoms(*args, return_masks=True)
        selected += self.selectByIndices(*args, return_masks=True)

        # exclude exsiting shells
        shells = self.shellatoms > -1
        nshells_exsiting = sum(shells)
        toadd = selected * ~shells
        print(" >>> Adding shells associated to atoms:\n", arange(self.natoms)[toadd])

        # resize self.shells
        nshells_to_add = sum(toadd)
        self.shells.natoms += nshells_to_add

        # assign values
        indices_added = range(nshells_exsiting, self._shells.natoms)
        self.shellatoms[toadd] = indices_added
        # update self.shells.coreatoms
        self._shells.coreatoms[self.shellatoms[self.shellatoms>-1]] = arange(self.natoms)[self.shellatoms>-1]
        self.shells.coords [indices_added] = self.coords [toadd] + displace
        self.shells.names  [indices_added] = self.names  [toadd]

        # assign the charges. If a dictionary is given, use this to assign charges; otherwise, use algebra.
        if charges:
            shellCharges = []
            # For each of the atomic numbers in the list, append the correct charge to the shellCharges list
            for eachAtom in self.shells.names[indices_added]:
                try:
                    shellCharges.append(charges[eachAtom.decode('UTF-8')])
                except KeyError:
                    # Should this throw an error? For now we just set the value to 0.0
                    shellCharges.append(0.0)

            # Set the Fragment's charges equal to the chargesToAdd list
            self.shells.charges[indices_added] = shellCharges
            # Subtract shell charges from core charges, to keep charge balance.
            self.charges[toadd] -= shellCharges
        else:
            positive = (self.charges > 0.0) * toadd
            # shell: Q*(1.0-factor) + shift
            self.shells.charges[self.shellatoms[positive]] = (1.0-factor)*self.charges[positive]
            self.shells.charges[self.shellatoms[positive]] += abs(shift)
            # core (Q): Q*factor - shift
            self.charges[positive] *= factor
            self.charges[positive] -= abs(shift)
            negative = (self.charges < 0.0) * toadd
            # shell: -Q*(1.0+factor) - shift
            self.shells.charges[self.shellatoms[negative]] = (1.0+factor)*self.charges[negative]
            self.shells.charges[self.shellatoms[negative]] -= abs(shift)
            # core (-Q): Q*factor + shift
            self.charges[negative] *= -factor
            self.charges[negative] += abs(shift)


    # YL 21/08/2019: added parameters to control the returned size
    def getSuperCell(self, nlayers=1, symmetric=True, redundant=False, accuracy=1.0E-2):
        '''Translate the unit cell to all 26 (3*3*3-1) directions and form a super cell'''

        from itertools import product
        from numpy     import all, append, dot, linalg, tile, zeros
        from ..base    import constants, lattice

        selectcases = { 3 : 3,
                        2 : 2,
                        1 : 1,
                        0 : 3 }

        supercell = self.__class__(natoms=0, connect_mode=self.connect_mode)


        coreatoms     = []
        shell_coords  = []
        shell_charges = []
        shell_names   = []

        i = 0
        nshells = 0

        # no translation will be made
        if nlayers == 0:
            nlayers = -1

        if symmetric:
            rangel = -nlayers
        else:
            rangel = 0
        ranger = nlayers+1

        for translation in product(range(rangel,ranger), repeat=selectcases[self._cell.ndimensions]):

            displacement = zeros((3,))
            # YL 24/06/2020: Cell.vectors's shape changed to (3,3):always
            displacement[:selectcases[self._cell.ndimensions]] = dot(self._cell.vectors[:selectcases[self._cell.ndimensions],:selectcases[self._cell.ndimensions]], \
                                                             translation)

            # YL 21/10/2020: there could be duplicates when applying periodicity to the unitcell (shells also changed accordingly below)
            coords = self.coords+displacement
            valids = [ i for i, c in enumerate(coords) if all(linalg.norm(c-supercell.coords, axis=1) > accuracy) ]

            unitCell = self.__class__(coords=coords[valids],        \
                                      znums=self.znums[valids],     \
                                      charges=self.charges[valids], \
                                      shellatoms=self.shellatoms[valids], \
                                      connect_mode=self.connect_mode)

            supercell.insert(unitCell, 0)

            shellatoms = unitCell.shellatoms[unitCell.shellatoms>-1]
#            coreatoms = append(coreatoms, self.natoms*i+self._shells.coreatoms)
            coreatoms = append(coreatoms, supercell.natoms-len(valids)+self._shells.coreatoms[shellatoms])
#            shell_coords  = append(shell_coords, self._shells.coords+displacement)
            shell_coords  = append(shell_coords , self._shells.coords[shellatoms]+displacement)
            shell_charges = append(shell_charges, self._shells.charges[shellatoms])
            shell_names   = append(shell_names  , self._shells.names[shellatoms])

            i += 1
            nshells += len(shellatoms)

        # cell
        supercell.cell.natoms = supercell.natoms
        supercell.cell.a = len(range(rangel,ranger))*self._cell.a
        supercell.cell.b = len(range(rangel,ranger))*self._cell.b
        supercell.cell.c = len(range(rangel,ranger))*self._cell.c
        supercell.cell.alpha = self._cell.alpha
        supercell.cell.beta  = self._cell.beta
        supercell.cell.gamma = self._cell.gamma

        # YL 13/07/2020: vectors must be updated after change to consts!
        supercell.cell.vectors

        # shells
        # more attributes need to be set?
        supercell.shells.resize(nshells)
        supercell.shells.coreatoms = coreatoms
        supercell.shells.coords    = shell_coords
        supercell.shells.charges   = shell_charges
        supercell.shells.names     = shell_names

        # if redundant atoms not allowed, the boundary atoms on the three surfaces (100), (010), and (001) should be removed
        if not redundant:
            redundant = supercell._selectRedundant(accuracy=accuracy)
            deleted = supercell.delete(redundant, mute=True)
            if chemsh.CHEMSH_DEBUG > 4:
                deleted.save('_redundant_boundary.xyz')

        return supercell


    def _selectRedundant(self, accuracy=1.0E-2, returnMask=False, bottom=True):
        '''Return the list of atoms on surfaces that may be redundant by periodic repeat'''

        from numpy import abs, arange, copy, cross, dot, linalg

        vectors = self.cell.vectors

        # translate centroid to origin
        centroid = copy(self.centroid)
        self.coords -= centroid

        # unit normal vectors of planes ab, ac, and bc
        plane_ab = cross(vectors[0,:], vectors[1,:])/linalg.norm(cross(vectors[0,:], vectors[1,:]))
        plane_ac = cross(vectors[0,:], vectors[2,:])/linalg.norm(cross(vectors[0,:], vectors[2,:]))
        plane_bc = cross(vectors[1,:], vectors[2,:])/linalg.norm(cross(vectors[1,:], vectors[2,:]))

        # default: (1¬00), (01¬0), and (001¬)
        # because for a surface calculation we normally add the adsorbate on the top
        if bottom:
            half = -0.5
        # (100), (010), and (001)
        else:
            half =  0.5

        # translate along negative half vector A
        self.coords -= vectors[:,0]*half
        # distances to plane BC (yz)
        distances_to_bc = abs(dot(-self.coords, plane_bc))
        selected_bc = distances_to_bc <= accuracy
        self.coords += vectors[:,0]*half

        # translate along negative half vector B
        self.coords -= vectors[:,1]*half
        # distances to plane AC (xz)
        distances_to_ac = abs(dot(-self.coords, plane_ac))
        selected_ac = distances_to_ac <= accuracy
        self.coords += vectors[:,1]*half

        # translate along negative half vector C
        self.coords -= vectors[:,2]*half
        # distances to plane AB (xy)
        distances_to_ab = abs(dot(-self.coords, plane_ab))
        selected_ab = distances_to_ab <= accuracy
        self.coords += vectors[:,2]*half

        # all selected
        selected = selected_bc + selected_ac + selected_ab

        # translate back
        self.coords += centroid

        if returnMask:
            return selected
        else:
            return arange(self.natoms)[selected]



# END OF FACILITIES #


    # YL 12/09/2020: this should be used as a wrapper version for FragmentMat
    def delete(self, todelete, mute=False, reconnect=True):
        '''Delete an atom/point or a range of atoms/points'''

        # TODO: cell

        from numpy  import arange, delete, searchsorted
        from copy   import copy
        from ..base import errors

        # original natoms
        natoms = copy(self.natoms)

        # get array of atom(s)/point(s) to delete
        todelete = self.getRange(todelete)

        # YL 17/02/2020: shells
        todelete_shells = self.shellatoms[todelete][self.shellatoms[todelete]>-1]

        # save the deleted and construct a new instance for returning
        deleted = self.__class__()
        deleted._master = self._master.take(todelete)
        # YL 17/02/2020: shells
        if todelete_shells.size:
            deleted.shells = self._shells.__class__()
            deleted._shells._master = self._shells._master.take(todelete_shells)

        if not mute:
            print("\n >>> the following part has been deleted:\n")
            if todelete_shells.size:
                print(deleted.coords2str(atoms='names', inserts=[(1,'core ')]))
                print(deleted._shells.coords2str(atoms=None, inserts=[(0,self.names[self.shells.coreatoms[todelete_shells]]),(1,'shell')]))
            else:
                print(deleted.coords2str(atoms='names'))

        if self._cell.natoms:
            print("\n >>> WARNING: the associated lattice atoms of the current fragment has not been deleted! (to be implemented))\n")

        # shift the indices of shells.coreatoms
        self.shells.coreatoms -= searchsorted(todelete, self.shells.coreatoms)

        # delete
        self._master = delete(self._master, todelete, axis=0)
        if todelete_shells.size:
            self._shells._master = delete(self._shells._master, todelete_shells, axis=0)

        # YL 05/05/2020: if there is `stacks` also shift all the indices
        indices = arange(natoms)
        indices = delete(indices, todelete)
        for key, val in self.stacks.items():
            new = searchsorted(indices, val)
            self.stacks.update({key:range(searchsorted(indices,val.start),searchsorted(indices,val.stop))})

        # YL 20/02/2020: reconnect after deletion
        if reconnect:
            self.connect(force_reconnect=True, mute=mute)

        return deleted

    def bind(self, adsorbate, donor='', distance=3.0, angle=180.0, unit='a.u.', surface='ab'):
        '''Bind an adsorbate to the target position'''

        from numpy  import absolute, arccos, copy, cross, dot, intersect1d, linalg, rad2deg, zeros
        from ..data import ligands, unitconvert
        from ..base import constants, lattice

        print("\n >>> Performing ChemShell automatic adsorbate binding on periodic systems...\n")

        surfaces = {  'a': [ 1, 0, 0],
                     '+a': [ 1, 0, 0],
                     'bc': [ 1, 0, 0],
                     'cb': [ 1, 0, 0],
                     '-a': [-1, 0, 0],
                    '-bc': [-1, 0, 0],
                    '-cb': [-1, 0, 0],
                      'b': [ 0, 1, 0],
                     '+b': [ 0, 1, 0],
                     'ac': [ 0, 1, 0],
                     'ca': [ 0, 1, 0],
                     '-b': [ 0,-1, 0],
                    '-ac': [ 0,-1, 0],
                    '-ca': [ 0,-1, 0],
                      'c': [ 0, 0, 1],
                     '+c': [ 0, 0, 1],
                     'ab': [ 0, 0, 1],
                     'ba': [ 0, 0, 1],
                     '-c': [ 0, 0,-1],
                    '-ab': [ 0, 0,-1],
                    '-ba': [ 0, 0,-1],
                   }

        rotations = { 0:1, 1:2, 2:0 }

        surface = surface.strip().lower()
        ranks = surfaces.get(surface, 'ab')

        # first loop up in the built-in adsorbates library
        adsorbateFrag = getattr(ligands, adsorbate.lower())

        # YL 10/02/2021: merged with the biological branch and now we must cast the type of
        #                adsorbateFrag which is initiated as FragmentBio to the type of self
        adsorbateFrag = adsorbateFrag.castTo(self)

        # choose the adsorbate's donor
        indices_donor = adsorbateFrag.select(atoms=donor, strict=False)
        # choose non-dummy atoms if nothing selected
        if not indices_donor.size:
            indices_donor = adsorbateFrag.select(atoms='X', invert=True)

        # YL 27/07/2023: we need a default index_donor
        index_donor = indices_donor[0]
        # choose the one connected a dummy atom if more than one selected by the user
        if indices_donor.size > 1:
            for i in indices_donor:
                if b'X' in adsorbateFrag.names[adsorbateFrag.conn[i][adsorbateFrag.conn[i]>-1]]:
                    index_donor = i
                    index_dummy = i
                    break

        # find a second atom that is connected and non-dummy
        connected = [adsorbateFrag.conn[index_donor][adsorbateFrag.conn[index_donor]>-1]]
        non_dummy = adsorbateFrag.select(atoms='X', invert=True)
        index_2nd = intersect1d(connected, non_dummy)[0]

        # coordinates of the donor atom
        coords_B = adsorbateFrag.coords[index_donor]
        coords_C = adsorbateFrag.coords[index_2nd]
        print(" Donor atom \"{}\" of adsorbate {} selected, coordinates in a.u.:\n {}".format({'':adsorbateFrag.symbols[index_donor]}.get(donor.strip(), donor.upper()),
                                                                                               adsorbate.upper(),
                                                                                               coords_B))


        # move to the centroid first
        coords_A = self.centroid
        adsorbateFrag.coords -= coords_B-coords_A

        selectcases = { 3 : 3,
                        2 : 2,
                        1 : 1,
                        0 : 3 }

        # displacement is half the lattice vector
        displacement = zeros((3,))
        displacement[:selectcases[self._cell.ndimensions]] = dot(self._cell.vectors[:selectcases[self._cell.ndimensions],:selectcases[self._cell.ndimensions]]*0.5, \
                                                         ranks)

        # add padding distance to the displacement vector
        displacement += displacement/linalg.norm(displacement)*(distance*unitconvert.unit2au.get(unit.strip().lower(), 1.0))
        adsorbateFrag.coords += displacement

        # current angle between the adsorbate molecule and the vector of displacement
        vector_BC = coords_B - coords_C
        vector_BA = coords_B - coords_A
        angle_adsorbate_vector = rad2deg(arccos(dot(vector_BC, displacement)/linalg.norm(displacement)/linalg.norm(vector_BC)))
        angle_to_rotate = angle - angle_adsorbate_vector

        # for convenience purpose, translate to origin before rotating
        translation = copy(coords_B)

        # rotate to re-orient
        normal = cross(vector_BC, vector_BA)
        # we rotate the adsorbate a bit if A, B, and C are on the same straight line already
        if linalg.norm(normal) < constants.LARGE:
            adsorbateFrag.coords -= translation
            # rotate about a lattice vector that is different from the current surface axis
            rank = rotations[absolute(ranks).tolist().index(1)]
            axis = self._cell.vectors[:selectcases[self._cell.ndimensions],:selectcases[self._cell.ndimensions]][rank,:]
            # any arbitrary angle
            adsorbateFrag.rotate(theta=15.0, vector=axis)
            adsorbateFrag.coords += translation
            # recalculate vectors
            vector_BC = coords_B - coords_C
            vector_BA = coords_B - coords_A
            # recalculate angle to rotate and the normal
            angle_adsorbate_vector = rad2deg(arccos(dot(vector_BC, displacement)/linalg.norm(displacement)/linalg.norm(vector_BC)))
            angle_to_rotate = angle - angle_adsorbate_vector
            normal = cross(vector_BC, vector_BA)

        # for convenience purpose, translate to origin before rotating
        adsorbateFrag.coords -= translation
        adsorbateFrag.rotate(theta=angle_to_rotate, vector=normal)
        adsorbateFrag.coords += translation

        # remove dummy atoms and merge
        adsorbateFrag = adsorbateFrag.getSelected(adsorbateFrag.select(atoms='X', invert=True))
        mergedFrag = self.merge(adsorbateFrag)

        # inherit the cell properties
        mergedFrag.cell.vectors[:] = self._cell.vectors[:]
        mergedFrag.cell.natoms = mergedFrag.natoms

        # save useful information
        mergedFrag._adsorbate = adsorbateFrag
        mergedFrag._adsorbate_surface = surface

        return mergedFrag


    def selectPeriodic(self, origin=None, dimensions=(1,1,1), accuracy=1.0E-2, surface='ab', visualise=''):
        ''''''

        from numpy  import abs, any, append, arange, argmin, around, array, cross, dot, floor, hstack, isin, integer, issubdtype, linalg, log10, min
        from numpy  import repeat, searchsorted, stack, tile, unique, vstack
        from ..base import lattice
        from ..io   import format
        
        selectcases = { 3 : 3,
                        2 : 2,
                        1 : 1,
                        0 : 3 }

        # saved symbols for choosing a surface <class dict>
        surfaces = lattice.surfaces

        surface = surfaces.get(surface.lower().strip(), 'ab')

        # check surface
        if surface != self._adsorbate_surface:
            print("\n WARNING: selectPeriodic() surface not same as previous adsorbate binding surface \"{}\"\n".format(self._adsorbate_surface))

        # default origin: centroid
        if origin is None or (type(origin) is str and origin.lower().strip() == 'centroid'):
            origin = self.centroid

        # optional origin:
        if type(origin) is int or issubdtype(type(origin), integer):
            try:
                origin = self.coords[origin]
            except:
                print("\n >>> WARNING: function selectPeriodic()'s keyword argument origin given wrong integer value\n")
                # TODO

        # lattice vectors
        vectors = self._cell.vectors[:selectcases[self._cell.ndimensions],:selectcases[self._cell.ndimensions]]

        # unit normal vectors of planes ab, ac, and bc
        plane_ab = cross(vectors[0,:], vectors[1,:])/linalg.norm(cross(vectors[0,:], vectors[1,:]))
        plane_ac = cross(vectors[0,:], vectors[2,:])/linalg.norm(cross(vectors[0,:], vectors[2,:]))
        plane_bc = cross(vectors[1,:], vectors[2,:])/linalg.norm(cross(vectors[1,:], vectors[2,:]))

        # masks array of centres inside the lattice
        masks_in_cell = self.selectOutOfCell(returnMask=True, invert=True)
        masks_out_cell = ~masks_in_cell

        # distances from all centres to plane AB
        distances_to_ab = dot(-self.coords[masks_in_cell], plane_ab)
        distances_adsorbate_to_ab = dot(-self.coords[~masks_in_cell], plane_ab)
        grouped_ab, tagged_ab = self.__getGroupedDists(distances_to_ab, accuracy)

        # distances from all centres to plane AC
        distances_to_ac = dot(-self.coords[masks_in_cell], plane_ac)
        distances_adsorbate_to_ac = dot(-self.coords[~masks_in_cell], plane_ac)
        grouped_ac, tagged_ac = self.__getGroupedDists(distances_to_ac, accuracy)

        # distances from all centres to plane BC
        distances_to_bc = dot(-self.coords[masks_in_cell], plane_bc)
        distances_adsorbate_to_bc = dot(-self.coords[~masks_in_cell], plane_bc)
        grouped_bc, tagged_bc = self.__getGroupedDists(distances_to_bc, accuracy)

        # distances from the origin point to three planes
        distance_origin_to_ab = dot(-origin, plane_ab)
        distance_origin_to_ac = dot(-origin, plane_ac)
        distance_origin_to_bc = dot(-origin, plane_bc)

        # length
        length = dimensions[0]
        # use the other two dimensions as length and width
        selectcases = { 'ab': [distances_to_ac, grouped_ac-distance_origin_to_ac, tagged_ac],
                        'ac': [distances_to_bc, grouped_bc-distance_origin_to_bc, tagged_bc],
                        'bc': [distances_to_ab, grouped_ab-distance_origin_to_ab, tagged_ab],
                      }
        distances, distance_origin_to_slices, tagged = selectcases[surface]
        min_dists = abs(distance_origin_to_slices).argsort()[:length]
        closest_length = isin(tagged, min_dists)
        # get the lattice parameters from the selection's span
        slices = distances[closest_length]
        length_constant = slices.max() - slices.min()

        # width
        width = dimensions[1]
        # use the other two dimensions as length and width
        selectcases = { 'ab': [distances_to_bc, grouped_bc-distance_origin_to_bc, tagged_bc],
                        'ac': [distances_to_ab, grouped_ab-distance_origin_to_ab, tagged_ab],
                        'bc': [distances_to_ac, grouped_ac-distance_origin_to_ac, tagged_ac],
                      }
        distances, distance_origin_to_slices, tagged = selectcases[surface]
        min_dists = abs(distance_origin_to_slices).argsort()[:width]
        closest_width  = isin(tagged, min_dists)
        # get the lattice parameters from the selection's span
        slices = distances[closest_width]
        width_constant = slices.max() - slices.min()

        # depth (nlayers)
        depth = dimensions[2]
        selectcases = { 'ab': [distances_to_ab, grouped_ab-distance_origin_to_ab, tagged_ab],
                        'ac': [distances_to_ac, grouped_ac-distance_origin_to_ac, tagged_ac],
                        'bc': [distances_to_bc, grouped_bc-distance_origin_to_bc, tagged_bc],
                      }
        distances, distance_origin_to_slices, tagged = selectcases[surface]
        min_dists = abs(distance_origin_to_slices).argsort()[:depth]
        closest_depth  = isin(tagged, min_dists)
        # get the lattice parameters from the selection's span
        slices = distances[closest_depth]

        # depth should plus the adsorbate atoms
        selectcases = { 'ab': distances_adsorbate_to_ab,
                        'ac': distances_adsorbate_to_ac,
                        'bc': distances_adsorbate_to_bc,
                      }
        slices = append(slices, selectcases[surface])
        depth_constant = slices.max() - slices.min()

        # save the calculated lattice constants of the selected atoms
        selectcases = { 'ab': [ width_constant , length_constant, depth_constant , self.cell.alpha, self.cell.beta, self.cell.gamma ],
                        'ac': [ length_constant, depth_constant , width_constant , self.cell.alpha, self.cell.beta, self.cell.gamma ],
                        'bc': [ depth_constant , width_constant , length_constant, self.cell.alpha, self.cell.beta, self.cell.gamma ],
                      }
        self._lattice_periodic_selection = selectcases[surface]

        atoms_in_cell = arange(self.natoms)[masks_in_cell]
        selected = atoms_in_cell[closest_length*closest_width*closest_depth]

        # also save the selected atoms for verification
        self._periodic_selection = selected
        self._surface = surface

        # visualisation
        if visualise.lower().strip() in [ 'vmd' ]:
            vmd_message = "To visualise the selection in VMD, copy the following clause and paste to the VMD \"Selected Atoms\" field:\n"
            vmd_select  = 'index ' + ' or index '.join([ str(s) for s in selected ])
            print('\n' + format.wrap(vmd_message))
            print('\n' + format.wrap(vmd_select ))

        return selected


    # this is a method independant of any previous operations
    def getSurface(self, adsorbate=[], accuracy=1.0E-2):
        '''Detect the surface of interest by finding which is the closest to the adsorbate'''

        from numpy   import abs, arange, around, count_nonzero, floor, linalg, log10, tile, vstack
        from ..base  import lattice

        adsorbate = utils.iterutils.getFlattened([adsorbate])

        if not len(adsorbate):
            adsorbate = self.selectOutOfCell()
            # print warning if there is cell but no adsorbate
            if not len(adsorbate) and count_nonzero(self.cell.vectors):
                print("\n WARNING: detectSurface() has no centre out of lattice given\n")
            return ''

        ndecimals = int(abs(floor(log10(accuracy))))

        selectcases = { 3 : 3,
                        2 : 2,
                        1 : 1,
                        0 : 3 }

        # lattice vectors
        # NB: vec[x/y/z, a/b/c]
        vectors = self._cell.vectors[:selectcases[self._cell.ndimensions],:selectcases[self._cell.ndimensions]]

        # volume of the lattice
        vol_ABC = linalg.det(vectors)

        # vectors of centres
        c = self.coords[adsorbate,:,None].T

        natoms_adsorbate = len(adsorbate)
        # vstack acts on the 1st axis
        stacked_AB = vstack((c, tile(self._cell.vectors[:,[0,1]], (natoms_adsorbate,1)).reshape(natoms_adsorbate,3,2).T)).T
        stacked_AC = vstack((c, tile(self._cell.vectors[:,[0,2]], (natoms_adsorbate,1)).reshape(natoms_adsorbate,3,2).T)).T
        stacked_BC = vstack((c, tile(self._cell.vectors[:,[1,2]], (natoms_adsorbate,1)).reshape(natoms_adsorbate,3,2).T)).T

        vol_AB = abs(around(linalg.det(stacked_AB)/vol_ABC, decimals=ndecimals)) > 1.0
        vol_AC = abs(around(linalg.det(stacked_AC)/vol_ABC, decimals=ndecimals)) > 1.0
        vol_BC = abs(around(linalg.det(stacked_BC)/vol_ABC, decimals=ndecimals)) > 1.0

        if vol_AB.sum():
            surface = 'ab'
        elif vol_AC.sum():
            surface = 'ac'
        elif vol_BC.sum():
            surface = 'bc'

        return surface


    def selectOutOfCell(self, accuracy=1.0E-2, returnMask=False, invert=False):
        '''Return indices or masks array of centres out of periodic lattice'''

        from numpy  import abs, arange, around, floor, linalg, log10, ones, tile, vstack, zeros
        from ..base import lattice

        # select nothing if non-periodic
        if self._cell.ndimensions == 0:
            masks_out_cell = zeros((self.natoms,), dtype=bool)

            if returnMask:
                if invert:
                    return ~masks_out_cell
                else:
                    return masks_out_cell
            else:
                if invert:
                    return arange(self.natoms)[~masks_out_cell]
                else:
                    return arange(self.natoms)[masks_out_cell]


        ndecimals = int(abs(floor(log10(accuracy))))

        selectcases = { 3 : 3,
                        2 : 2,
                        1 : 1,
                        0 : 3 }

        # lattice vectors
        # NB: vec[x/y/z, a/b/c]
        vectors = self._cell.vectors[:selectcases[self._cell.ndimensions],:selectcases[self._cell.ndimensions]]

        # volume of the lattice
        vol_ABC = linalg.det(vectors)

        # do nothing if no cell
        if not vol_ABC:
            if returnMask:
                if invert:
                    return ones((self.natoms,), dtype=bool)
                else:
                    return zeros((self.natoms,), dtype=bool)
            else:
                if invert:
                    return arange(self.natoms)
                else:
                    return arange(0)

        # vectors of centres
        c = self.coords[:,:,None].T

        # vstack acts on the 1st axis
        stacked_AB = vstack((c, tile(self._cell.vectors[:,[0,1]], (self.natoms,1)).reshape(self.natoms,3,2).T)).T
        stacked_AC = vstack((c, tile(self._cell.vectors[:,[0,2]], (self.natoms,1)).reshape(self.natoms,3,2).T)).T
        stacked_BC = vstack((c, tile(self._cell.vectors[:,[1,2]], (self.natoms,1)).reshape(self.natoms,3,2).T)).T

        vol_AB = abs(around(linalg.det(stacked_AB)/vol_ABC, decimals=ndecimals)) > 1.0
        vol_AC = abs(around(linalg.det(stacked_AC)/vol_ABC, decimals=ndecimals)) > 1.0
        vol_BC = abs(around(linalg.det(stacked_BC)/vol_ABC, decimals=ndecimals)) > 1.0

        # centres in-/outside the lattice
        masks_out_cell = vol_AB + vol_AC + vol_BC

        if returnMask:
            if invert:
                return ~masks_out_cell
            else:
                return masks_out_cell
        else:
            if invert:
                return arange(self.natoms)[~masks_out_cell]
            else:
                return arange(self.natoms)[masks_out_cell]


    # many are repeated as selectPeriodic() but we have to make the current method independent
    def getLatticeOfSelectedPeriodic(self, selected, accuracy=1.0E-2, surface='ab', recalculate=False, redundant=True):
        '''Return the lattice constants of atoms chosen by periodic selection'''

        from numpy   import abs, arange, around, array, array_equal, cross, diff, dot, floor, isin, linalg, log10, mean
        from numpy   import searchsorted, sort, std, unique, zeros
        from ..base  import geometry, lattice

        selected = utils.iterutils.getFlattened([selected])

        # if Fragment.selectPeriodic() has already been run and lattice of selection was saved
        try:
            # verify the selected atoms
            if array_equal(utils.iterutils.getFlattened([selected]), self._periodic_selection) and not recalculate:
                return self._lattice_periodic_selection
        except:
            pass

        coords = self.coords[selected]

#        ndecimals = int(abs(floor(log10(accuracy))))

        selectcases = { 3 : 3,
                        2 : 2,
                        1 : 1,
                        0 : 3 }

        # saved symbols for choosing a surface <class dict>
        surfaces = lattice.surfaces

        surface = surfaces.get(surface.lower().strip(), 'ab')

        # lattice vectors
        vectors = self._cell.vectors[:selectcases[self._cell.ndimensions],:selectcases[self._cell.ndimensions]]

        # unit normal vectors of planes ab, ac, and bc
        plane_ab = cross(vectors[0,:], vectors[1,:])/linalg.norm(cross(vectors[0,:], vectors[1,:]))
        plane_ac = cross(vectors[0,:], vectors[2,:])/linalg.norm(cross(vectors[0,:], vectors[2,:]))
        plane_bc = cross(vectors[1,:], vectors[2,:])/linalg.norm(cross(vectors[1,:], vectors[2,:]))

        # YL 18/11/2020: the method based on rounding is very unreliable to reflect the resolution request
        # distances from all centres to plane AB (signed!)
        distances_to_ab = dot(-coords, plane_ab)
        grouped_ab, tagged_ab = self.__getGroupedDists(distances_to_ab, accuracy)

        # distances from all centres to plane AC (signed!)
        distances_to_ac = dot(-coords, plane_ac)
        grouped_ac, tagged_ac = self.__getGroupedDists(distances_to_ac, accuracy)

        # distances from all centres to plane BC (signed!)
        distances_to_bc = dot(-coords, plane_bc)
        grouped_bc, tagged_bc = self.__getGroupedDists(distances_to_bc, accuracy)

        def __getPaddedDists(_tags, _distances, _coords, _plane, _accuracy, _redundant):
            '''Get the unit vector length (minimal distance between symmetry-redundant centres along a vector)'''

            # we don't calculate this distance if redundant centres are allowed in cell
            if _redundant:
                return _distances.max() - _distances.min()

            # number of centres in the first layer
            tags_0      = _tags==0
            ntagged_0   = tags_0.sum()
            paddings    = {}
            others      = {}
            dists0      = _distances[tags_0]
            mean_dists0 = mean(dists0)

            # centres of the first layer projected on plane
            projected_0 = geometry.projectOntoPlane(_coords[tags_0], _plane[0], _plane[1], _plane[2])

            # loop over the following layers until we find symmetry-redundant layers
            for i in range(1, _tags.max()+1):
                tags_i    = tags==i
                ntagged_i = tags_i.sum()
                if ntagged_i == ntagged_0:

                    # TODO: defect
                    # project the points of the current layer to the plane
                    matched = zeros(ntagged_0, dtype=bool)
                    projected_i = geometry.projectOntoPlane(_coords[tags_i], _plane[0], _plane[1], _plane[2])
                    # check the distances
                    for p in projected_i:
                        matched += linalg.norm(projected_0-p, axis=1) <= _accuracy
                    # we call the current layer identical pattern only if all points match
                    if matched.sum() < ntagged_0:
                        continue

                    # distance between the current layer and the first
                    displacements = _distances[tags_i] - mean_dists0

                    if std(displacements) <= _accuracy:
                        paddings.update({i:mean(displacements)})
                # other layers (not in use)
                else:
                    displacements = _distances[tags_i] - mean_dists0
                    others.update({i:mean(displacements)})

            # check the evenness (assuming the lattice is even without adsorbate)
            if len(paddings) > 1:
                try:
                    diff_keys = diff(list(paddings.keys()  ))
                    diff_vals = diff(list(paddings.values()))
                    if std(diff_keys) <= _accuracy and std(diff_vals) <= _accuracy:
                        padding = abs(mean(diff_vals))
                    else:
                        print("\n WARNING: Fragment.getLatticeOfSelectedPeriodic() cannot calculate the cell's vector length\n")
                        padding = 0.0
                except:
                    print("\n WARNING: Fragment.getLatticeOfSelectedPeriodic() cannot calculate the cell's vector length\n")
                    padding = 0.0

                # nearest and furthest layers
                large = max(paddings.values())
                small = min(paddings.values())
                # if same sign
                if large*small > 0:
                    _result = max([abs(large),abs(small)]) + padding
                else:
                    _result = large - small + padding
                
            # only the 0th layer (origin) and one more layer
            elif len(paddings) == 1:
                padding = abs(mean(list(paddings.values())))
                _result = padding*2.0
            # no padding
            else:
                _result = _distances.max() - _distances.min()

            return _result

        # length
        # use the other two dimensions as length and width
        dict_distances  = { 'ab': distances_to_ac, 'ac': distances_to_bc, 'bc': distances_to_ab }
        dict_tags       = { 'ab': tagged_ac      , 'ac': tagged_bc      , 'bc': tagged_ab       }
        dict_planes     = { 'ab': plane_ac       , 'ac': plane_bc       , 'bc': plane_ab        }
        distances       = dict_distances[surface]
        tags            = dict_tags[surface]
        plane           = dict_planes[surface]
        length_constant = __getPaddedDists(tags, distances, coords, plane, accuracy, redundant)
        # get the lattice parameters from the selection's span

        # width
        # use the other two dimensions as length and width
        dict_distances  = { 'ab': distances_to_bc, 'ac': distances_to_ab, 'bc': distances_to_ac }
        dict_tags       = { 'ab': tagged_bc      , 'ac': tagged_ab      , 'bc': tagged_ac       }
        dict_planes     = { 'ab': plane_bc       , 'ac': plane_ab       , 'bc': plane_ac        }
        distances       = dict_distances[surface]
        tags            = dict_tags[surface]
        plane           = dict_planes[surface]
        width_constant  = __getPaddedDists(tags, distances, coords, plane, accuracy, redundant)
        # get the lattice parameters from the selection's span

        # depth (nlayers)
        dict_distances  = { 'ab': distances_to_ab, 'ac': distances_to_ac, 'bc': distances_to_bc }
        dict_tags       = { 'ab': tagged_ab      , 'ac': tagged_ac      , 'bc': tagged_bc       }
        dict_planes     = { 'ab': plane_ab       , 'ac': plane_ac       , 'bc': plane_bc        }
        distances       = dict_distances[surface]
        tags            = dict_tags[surface]
        plane           = dict_planes[surface]
        depth_constant  = __getPaddedDists(tags, distances, coords, plane, accuracy, redundant)
        # get the lattice parameters from the selection's span

        # save the calculated lattice constants of the selected atoms
        selectcases = { 'ab': [ width_constant , length_constant, depth_constant , self.cell.alpha, self.cell.beta, self.cell.gamma ],
                        'ac': [ length_constant, depth_constant , width_constant , self.cell.alpha, self.cell.beta, self.cell.gamma ],
                        'bc': [ depth_constant , width_constant , length_constant, self.cell.alpha, self.cell.beta, self.cell.gamma ],
                      }

        # set vectors to zeros so that when constants are accessed they won't be rewritten by
        # triggering vectors2consts and vectors will be freshly calculated by consts2vectors
        self._cell.vectors[:] = 0.0

        return selectcases[surface]


    def __getGroupedDists(self, _dists, _resolution):
        '''Return grouped distances according to given resolution'''

        from numpy import argsort, argwhere, array, diff, full, mean, vstack

        _grouped = []
        _tagged = full(_dists.shape, -1, dtype=int)

        # NB: _dists may be signed!
        _asort = argsort(_dists)
        _adiff = diff(_dists[_asort])
        _anode = argwhere(abs(_adiff)>=_resolution).T[0] + 1
        _aindices = vstack(([0]+_anode.tolist(), _anode.tolist()+[_adiff.size+1])).T

        for i, (s, e) in enumerate(_aindices):
            _grouped += [ mean(_dists[_asort[s:e]]) ]
            _tagged[_asort[s:e]] = i
        undefined = _tagged == -1

        return _grouped, _tagged


# synonyms
Frag = Fragment

