#  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 dot, full, count_nonzero
from numpy.linalg import norm
from .            import Objects
from ..dl_py2f    import utils

class _XYZ(Objects):
    '''Base class of all cartesian-based objects (fragments, gradients, etc.)'''


    _kwargs = {
                # YL: do NOT include `_master` in initialisation items nor
                #     `def _master(self):` as @property, or it will bring
                #      many troubles, for example `npts` (i.e., natoms, nbqs,
                #      etc.) as a prerequisite
#               '_master': utils.nputils.RecArray(shape=(0,))

                'names': full((1,), '', dtype='S8'),
                'title': 'title'+' '*123,            # the initial lenghth determines the max size `title` string holds
                'unit' :'angstrom',
              }

    _inherit = '_fields'

    _init    = '_fields'

    _fields  = 'names'


# YL 05/09/2023: brought above the fragment properties so to use these methods
# WRAPPER METHODS #


    def _resize(setter):
        '''Resizer as a decorator'''

        # wrapper of setter
        def wrapper(self, buffer):

            # convert to ndarray
            abuff = utils.objutils.checkType(buffer, list)

            # this handles the case that buffer is a 1-D array
            if abuff.ndim > 1:
                self.natoms = abuff.shape[0]

            setter(self, abuff)

        return wrapper


    def _resize1D(setter):
        '''Resizer as a decorator'''

        # wrapper of setter
        def wrapper(self, buffer):

            from numpy       import ndarray, ravel
            # YL 20/10/2020: using or importing the ABCs from 'collections' instead
            #                of from 'collections.abc' is deprecated since Python 3.3,
            #                and in 3.9 it will stop working
            from collections import abc

            # array-like
            if isinstance(buffer, (ndarray, abc.Sequence)):
                # convert to flattened ndarray
                abuff = ravel(utils.objutils.checkType(buffer, list))
                self.natoms = abuff.size
                setter(self, abuff)
            # scalar
            else:
                setter(self, buffer)

        return wrapper


# END OF WRAPPER METHODS #

# PROPERTIES

    @property
    def title(self):
        '''Title of the object'''

        return self._title


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

        if type(val) is str:
            self._title = val
        else:
            try:
                self._title = ' '.join([ word for word in val.astype(str) ])
            except:
                try:
                    self._title = ' '.join([ word for word in val ])
                except:
                    # hack to fix 0-d array exception for one word punch file titles - IMPROVE!
                    self._title = str(val)


    @property
    def names(self):
        '''Point names'''

        return self._master.names


    # YL 05/09/2023: added a resize method similarly to FragBase.charges so that names can also be defined prior to `coords` (priority: charges > names > coords)
    @names.setter
    @_resize1D
    def names(self, val):
        '''Setter of names'''

        from numpy  import count_nonzero
        from ..base import forcefield

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

        # also set znums
        # YL 30/05/2018: znums has been removed from Shell
        try:
            if not count_nonzero(self.znums):
                self._master.setField(forcefield.names2znums(self.names), field='znums')
        except:
            pass


    @property
    def unit(self):
        '''Unit'''

        return self._unit


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

        self._unit = val.lower()


    @property
    def coords(self):
        '''Cartesian coordinates'''

        return self._master.xyz


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

        self._master.setField(val, field='xyz')
        self._reconnect = True


    @property
    def x(self):
        '''X coordinates'''

        return self._master.xyz[:,0]


#    @x.setter
#    def x(self, val):
#        '''Setter of X coordinates'''
        

    @property
    def y(self):
        '''Y coordinates'''

        return self._master.xyz[:,1]


    @property
    def z(self):
        '''Z coordinates'''

        return self._master.xyz[:,2]


    def dist(self, *args):
        '''Distance between two points'''

        from numpy import linalg

        return linalg.norm(self.diff(*args))


    def unitvec(self, *args):
        '''Unit vector formed by two points'''

        return self.diff(*args)/self.dist(*args)


    def vec(self, *args):
        '''Vector formed by two points'''

        return self.diff(*args)


    # TODO: need jatom - iatom?
    def diff(self, *args):
        '''ipoint - jpoint'''
#        if len(args) == 1 and type(args[0]) is tuple:
        # tuple
        try:
            return self._master.xyz[args[0][0]] - self._master.xyz[args[0][1]]
#        elif len(args) == 2 and utils.objutils.getType(args[0], list_modules=chemsh.CHEMSH_MODS) is int and utils.objutils.getType(args[1]) is int:
        # pair of int indices
        except:
            return self._master.xyz[args[0]] - self._master.xyz[args[1]]
        else:
            print("### xyz.diff: ERROR! args =", args)
            return None


    # ensure natoms is always up-to-date
    @property
    def natoms(self):
        '''Number of atoms'''

        return self._master.shape[0]


    @natoms.setter
    def natoms(self, val):
        '''Number of atoms'''

        # if _master exists
        try:
            self._master.expand(val)
        # create one if not
        except:
            self._master = utils.nputils.RecArray(shape=(val,), dtype=self.__dtype)

        # resizing fragment should require reconnecting
        self._reconnect = True


    @property
    def _indices_all(self):
        '''Indices of all centres'''

        from numpy import arange

        return arange(self.natoms)


    @property
    def __dtype(self):
        '''Datatype of _master'''

        fields = [ self._synons[field] if field in self._synons else
                   field for field in utils.iterutils.getFlattened((self._fields,)) ]

        # construct dtypes according to _kwargs and _fields
        dtypes = [ [field, self._kwargs[field].dtype] if self._kwargs[field].ndim == 1 else
                   [field, self._kwargs[field].dtype, self._kwargs[field].shape[1]] if self._kwargs[field].ndim == 2 else
                   [field, self._kwargs[field].dtype, 0]
                   for field in utils.iterutils.getFlattened((self._fields,)) ]


# YL: version B
#        # construct dtypes according to _kwargs and _fields
#        dtypes = [ [field, self._kwargs[field].dtype] if self._kwargs[field].ndim is 1 else
#                   [field, self._kwargs[field].dtype, self._kwargs[field].shape[1]] if self._kwargs[field].ndim is 2 else
#                   [field, self._kwargs[field].dtype, 0]
#                   # skip embedded structured array for now, e.g., zmatrix
#                   for field in utils.iterutils.getFlattened((self._fields,)) if not self._kwargs[field].dtype.names ]
#
#        # decompose embedded structured array, e.g., zmatrix
#        for field in utils.iterutils.getFlattened((self._fields,)):
#            if self._kwargs[field].dtype.names:
#                dtypes += [ [name, self._kwargs[field].dtype[name]] for name in self._kwargs[field].dtype.names ]

        # check synonyms
        dtypes = [ tuple([self._synons[dtype[0]]]+dtype[1:]) if dtype[0] in self._synons else
                   tuple(dtype) for dtype in dtypes ]

        # initialisers
        values = [ self._kwargs[field] for field in utils.iterutils.getFlattened((self._fields,)) ]

        # final dtypes
        return list(zip(dtypes, values))


# YL: since `_fields` is already (more conveniently) defined at the beginning of the class,
#     do NOT define it as an on-the-fly @property
#    def _fields(self):
#        '''All fields items will enter self._master array'''
#        from ..dl_py2f import utils
##        return utils.iterutils.getFlattened((self.__fields,))
#        return self.__fields
#    @_fields.setter
#    def _fields(self, val):
#        '''Setter of _fields'''
#        print("### Setter of _fields", self, val)
#        # to ensure _fields is still a tuple rather than a string given class-level, for example, _fields = 'bla',
#        # we have to apply the method in _fields.setter rather than in @property definition, because
#        # self._fields in defined prior to the property
#        self._fields = utils.iterutils.getFlattened((val,))


#    @property
#    def _master(self):
#        '''Master array'''
#
#        # if gets AttributeError, for example, 'Cell' object has no attribute '__XYZ__master', check initialisation in utils.objutils and @_master.setter
#        return self.__master
#
#    # this makes sure `_master` is initialised by creation
#    @_master.setter
#    def _master(self, val):
#        '''Setter of _master'''
##        self.__master = utils.nputils.RecArray(shape=(0,), dtype=self.__dtype)
#
#        self.__master = val

        # do nothing: _master can only be resized but never redefined
#        print("\n ChemShell >>> WARNING: _master array of", self.__class__.__name__, "instance must not be redefined (but can be initialised/resized via natoms, etc.)! Quitting...\n")
#        quit(1)


    # YL: Python cannot keep the order of kwargs yet, see: http://legacy.python.org/dev/peps/pep-0468
    def coords2str(self,
                   coords    =  9,        # <int>   precision of float number printing, i.e. 9 means '%.9f', None means no printing
                                          # YL 10/07/2019: it is accepted now, e.g., coords = '%8.3f'
                   charges   =  6,        # <int>   precision of float number printing, i.e. 9 means '%.9f', None means no printing
                                          # YL 10/07/2019: it is accepted now, e.g., charges = '%8.3f'
                   inserts   = [],        # <list>  list of tuples, [(xyz, content/, fmt/)], e.g.:inserts=[(1,'core','%8s'), (6, 3.1415926)]
                                          #         a properly sized array can also be inserted, e.g.: inserts=[(2,frag.names)]
                                          #         in practice this can be used to generate any free-formatted blocks
                   select    = [],        # <list>  a list of select atoms/points to print
                   indent    = '',        # <str>   indentation
                   delimiter = ' '*2,     # <str>   separator
                   atoms     = 'symbols', # <str>   'symbols' | 'znums' | 'names' | 'types'; None means no printing
                   order     = 'acq',     # <str>   'a': atoms; 'c': coords; 'q': charges; 'z': znums
                   unit      = 'a.u.',    # <str>   target unit for printing
                   notations = {},        # <dict>  notations: 'f' for float, 'e' for scientific
                   align     = 'left',
                  ):
        '''Convert coords to strings (possibly with column-wise inserted contents)'''

        from numpy  import array, dtype, ndarray, savetxt
        from re     import sub         as re_sub
        from io     import BytesIO
        from ..data import unitconvert

        # `natoms` may have a different definition in some derived classes (for example in <class Grid>)
        npoints = self._master.shape[0]

        synonyms = { 'atoms'  :'a',
                     'coords' :'c',
                     'charges':'q',
                     'x'      :'x',
                     'y'      :'y',
                     'z'      :'z'
                   }

        def _getInsert(_field, fmt=False):
            # try to return formats provided in 'inserts'
            if fmt:
                try:
                    return inserts[int(_field.split('insert_')[1])][2]
                except:
                    None
            else:
                return inserts[int(_field.split('insert_')[1])][1]

        # ordered fields (with inserts)
        fields = [  atoms    if i.lower() == 'a' else
                   'coords'  if i.lower() == 'c' else
                   'charges' if i.lower() == 'q' else
                   'znums'   if i.lower() == 'z' else
                    None for i in order ]

        # replace 'coords' with 'x', 'y', and 'z'
        if 'coords' in fields:
            pos = fields.index('coords')
            fields.remove('coords')
            fields.insert(pos, 'z')
            fields.insert(pos, 'y')
            fields.insert(pos, 'x')

        # exclude unwanted
        for i in fields:
            # for example, charges are not printed if charges=0 or charges=None
            if not i or ( i in locals() and not locals()[i] ):
                fields.remove(i)

        # insert in the fields
        for i in range(len(inserts)):
            fields.insert(inserts[i][0], 'insert_'+str(i))

        # original dtype
        dtp = dtype(sum([ self.getDtype(f) if not f.startswith('insert_') else
                          utils.objutils.type2dtype(_getInsert(f), descr=True)
                          for f in fields ],
                        [ ]))

        # name the dtype's fields
        dtp.names = fields

        descr = dtp.descr
        for f in fields:
            # replace all bytes to (unicode) strings: cannot replace a dtype item but can modify the descriptor
            if dtp[f].char == 'S':
                # keep the bytesize
                descr[fields.index(f)] = (descr[fields.index(f)][0], 'U'+str(dtp[f].itemsize))

        # construct a new dtype according to descriptor
        dtp = dtype(descr)
       

        # returns 0 if x==0; returns y otherwise
        getlen = lambda x, y: y if x !=0 else 0

        # default notations: float
        default_notations = { 'c':'f',
                              'q':'f' }

        # customised notations
        default_notations.update(notations)

        # prepare formats
        # fixed-width float number %w.mf: w = m + m/2
        # int(m or 0) can avoid error when m is None, while str(None) yielding 'None' does not matter
        # YL 12/03/2024: to right align float in scientific notation at leat 7 extra positions are needed, for example %21.14e
        try:
            c_fmt = '%'+str(int(coords  or 0)+getlen(int(coords  or 0),7))+'.'+str(coords) +default_notations['c']
        except:
            c_fmt = coords
        try:
            q_fmt = '%'+str(int(charges or 0)+getlen(int(charges or 0),7))+'.'+str(charges)+default_notations['q']
        except:
            q_fmt = charges

        if align.lower().strip().startswith('l'):
            symbol_align = '-'
        else:
            symbol_align = ''

        selectcases = {
                        'charges': q_fmt,
                        'coords' :(c_fmt)*3,
                        'x'      : c_fmt,
                        'y'      : c_fmt,
                        'z'      : c_fmt,
                        'symbols': f'%{symbol_align}4s',
                        'znums'  : f'%{symbol_align}4d',
                        'names'  : f'%{symbol_align}6s',          # six characters for names
                        'types'  : f'%{symbol_align}6s',          # six characters for types
                         str     : lambda x: f'%{symbol_align}s',
                         int     : lambda x: '%d',
                         float   : lambda x: '%f',
                         ndarray : lambda x: '%'+utils.nputils.dtype2fmt(x.dtype)
                      }

        fmt = indent + delimiter.join([ selectcases[f]          if f in selectcases        else 
                                        _getInsert(f, fmt=True) if _getInsert(f, fmt=True) else
                                        selectcases[type(_getInsert(f))](_getInsert(f))
                                        for f in fields ])

        if chemsh.CHEMSH_DEBUG > 5:
            print("\n >>> CHEMSHELL _XYZ.coords2str: format =\n {}\n".format(fmt))

        # customised array
        abuff = ndarray((npoints,), dtype=dtp)

        # take selected atoms
        if len(select):
            abuff = abuff.take(select, axis=0)
        # take all
        else:
            select = range(npoints)

        for f in fields:
            # first, convert bytes to string
            try:
                abuff[f] = getattr(self, f).take(select).decode()
            except AttributeError:
                # second, try to assign other attributes
                try:
                    abuff[f] = getattr(self, f).take(select)
                # last, insert non-existing attributes
                except AttributeError:
                    abuff[f] = _getInsert(f)

            # unit conversion: from a.u. to any requested
            if f in [ 'x', 'y', 'z', 'coords' ]:
                abuff[f] *= unitconvert.au2unit[unit.lower()]


            # rounding
#            if abuff[f].dtype.kind == 'f':
#                abuff[f] = abuff[f].round(utils.miscutils.getPrecision(selectcases[f]))
#                abuff[f] = [ round(num, utils.miscutils.getPrecision(selectcases[f])) for num in abuff[f] ]

        # convert to bytes: BytesIO, not StringIO! 
        bbuff = BytesIO()

        savetxt(
                 bbuff,
                 abuff,
#                 delimiter='',
#                 newline='\n ',
                 fmt=fmt
                )

        # convert bytes to string
        outstr = bbuff.getvalue().decode().rstrip() + '\n'

        return outstr


    def getSelected(self, selected):
        '''Return a new instance with selected centres'''

        selected = utils.iterutils.getFlattened(selected)

        # create a new instance of the current class
        inst = super().__self_class__()

        # determine the size
        inst.natoms = len(selected)

        # copy data
        # numpy.take() is faster than fancy indexing!
        # YL 24/10/2019: changed to fancy indexing to allow masks
        inst._master.setField(self._master[selected])

        return inst


    def _getDistsPoints2Plane(self, normal_ABC, d, range=...):
        '''Return the distances from (selected) points to the plane (a*x+b*y+c*z+d=0) defined by the normal, which must be a (3,)-shaped array, and d'''

        dists = (dot(self.coords[range], normal_ABC)+d)/norm(normal_ABC)

        return dists


# YL FIXME: to move and merge into `select()`
    def getRegion(self, *args):
        '''Return atom indices of a labelled region according to `names` (e.g., 4 for all 'O4' and 'Mg4' atoms)'''

        from numpy      import array, concatenate, where, __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

        # allows frag.getRegion(4), frag.getRegion(1,2,3,...), and frag.getRegion([1,2,3])
        int2region = lambda x: where(defchararray.count(defchararray.decode(self.names), str(x)))[0]

        try:
            return concatenate(tuple([ int2region(arg) for arg in utils.iterutils.getFlattenedIter(args) ]))
        except:
            return []


#    def getRange(self, iter):
#        '''Return atom indices for a given (nested) list'''
#
#        from numpy     import where
#        from ..dl_py2f import utils
#
#        # flatten
#        abuff = utils.iterutils.getFlattenedIter(iter)
#
#        # translate string species to indices
#        # `self.atomtypes` already renamed as `self.symbols`
#        try:
#            abuff1 = where([ self.symbols == i.title() for i in abuff if type(i) is str ])[1].tolist()
#        # if none
#        except:
#            abuff1 = []
#        print("abuff1 =", abuff1)
#
#        # find integer (could be numpy.int!) indices
#        abuff2 = [ i for i in abuff if utils.objutils.getType(i) is int ]
#
#        # merge the two lists
#        abuff1.extend(abuff2)
#
#        return sorted(set(abuff1))


    @property
    def stacks(self):
        '''Stacks that are used by merge() to form the current object'''

        # if stacks is missing
        if not hasattr(self, '_stacks'):
            setattr(self, '_stacks', {self.__class__.__name__:range(self.natoms)})

        # return a dictionary
        return self._stacks


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

        self._stacks = val


    @property
    def _width(self):
        '''Total width of array'''

        from numpy import array

        # for example, width of frag is 1(names)+1(znums)+3(coords)+1(charges)+14(conn)+1(frozen)+1(masses) = 22
        # YL 28/01/2020: allow embedded structure such as zmatrix that contains 6 sub-fields
        width_list = [ self._master[field].shape[1] if self._master[field].ndim == 2 else
                       array(self._master.dtype.fields[field][0].names).size for field in self._master.dtype.fields ]

        return sum(width_list)


# MANIPULATION METHODS #


    def reorder(self, order):
        '''Reorder centres in place according to the given sequence'''

        from numpy import array

        if array(order).size != self.natoms:
            print("\n >>> ERROR: indexing array's shape must match the fragment")

        else:
            self._master[:] = self._master[order][:]

   
    def resize(self, newsize):
        '''Resize: same as redefining natoms'''

        self.natoms = newsize


    def delete(self, todelete, mute=False):
        '''Delete an atom/point or a range of atoms/points'''

        # TODO: cell and shells

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

        if self.__class__.__name__ not in [ 'Shells', 'BQs' ] and not self.__class__.__name__.startswith('Fragment'):
            raise errors.ChemShNAError('delete', self)

        indices = arange(self.natoms)

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

        # save the deleted and construct a new instance for returning
        deleted = self.__class__()
        deleted._master = self._master.take(todelete)

        # delete
        self._master = delete(self._master, todelete, axis=0)

        # YL 05/05/2020: if there is `stacks` also shift all the indices
        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))})

        if not mute and len(todelete):
            try:
                strbuff = deleted.coords2str()
                print("\n >>> the following part has been deleted:\n{}".format(strbuff))
            except:
                pass

        return deleted


    # AJL, Aug2017: Picked up a bug here. shells.shellatoms are not updated to reflect these changes, where the coordinating
    # atom for a shell may have changed. TODO: Automate update - on first attempt it seems shellatoms can't be manually altered?
    def insert(self, toinsert, pos=0):
        '''Insert an atom/point or a range of atoms/points in place'''

        from numpy     import insert as np_insert
        from ..data    import atomicdata

        # to insert one or more atoms by providing atomic symbols and position(s)
        try:
            inst = self.__class__()
            # YL 28/12/2021: must use unique=False to enable adding multiple same species of atoms
            inst.znums = [ atomicdata.symbols.index(symbol.title()) for symbol in utils.iterutils.getFlattened([toinsert], unique=False) ]
            toinsert = inst
        # to insert a same-kind _XYZ-object, pos can be either a int or a listi/slice:
        #     xyzObj.insert(toinsert, 0)
        # for the latter case, pos should be of consistent length, for example:
        #     xyzObj.insert(toinsert, [1,3,5,7]), if _XYZ-object `toinsert` contains four atoms/points 
        except:
            pass

        # insert
        self._master = np_insert(self._master, pos, toinsert._master, axis=0)

# YL 21/04/2018: moved to cluster/__init__.py for _XYZ is generic
#        # AJL, Aug2017: shift shellatoms that need updating
#        if self.nshells > 0:
#            self.shellatoms = [(i + toinsert.natoms) if i > pos else i for i in self.shellatoms]
#            # Below is a horrible hack because in calling nshells we overwrite self.shellatoms with self._shells.coreatoms
#            # so essentially here I'm just protecting against that information overwriting our new self.shellatoms
#            # It'd be much nicer if the request for nshells didn't need to overwrite shellatoms (I can't follow
#            # the logic why self.shellatoms is needed anyway?) by just accessing the shell object directly? And then this 
#            # could be removed. Easy test of the effect of this code is to edit out this code, do an insert and then look
#            # self.shells.names, where you'll see the indices haven't increased and you'll get e.g. region 2 atoms with shells.
#            self._shells.coreatoms = self.shellatoms

    # AJL, Aug2017: It's not possible to give pos -1 for insert and get the atoms added at the end, so here I have
    # created a short routine that'll do this for the user. It is similar to merge but users the machinary for insert,
    # which is safer for keeping the characteristics of a certain object type at present (e.g. insert xyz into punch)
    def append(self, toappend):
        '''Append an atom/point or a range of atoms/points'''

        self.insert(toappend, self.natoms)


    def merge(self, *args, noself=False, cls=None, **kwargs):
        '''Stack _XYZ-type data according to their common fields of _master'''

        from numpy     import vstack
        from itertools import accumulate
        from inspect   import getmro
        from .xyz      import _XYZ

        # omit non-_XYZ-class and should NOT omit zero-sized objects (required in, for example, QM/MM additive scheme)
        sources = [ arg for arg in args if _XYZ in getmro(arg.__class__) ]

        # do the merge in place
        if noself:
            new = self
        # create a new instance if not in place
        else:
            if cls:
                new = cls()
            else:
                new = self.__class__()
            # YL 09/10/2018: sources here should contain `self`
            sources = [ self ] + sources
        nsources = len(sources)

        # with repeats
        classes = [ sources[i].__class__.__name__ for i in range(nsources) ]
        # label repeats
        classes = [ classes[i]+str(i) if classes.count(classes[i]) > 1 else classes[i]
                    for i in range(nsources) ]

        # resize
#        self.natoms = sum([ arg.natoms for arg in sources ])
        new.natoms = sum([ source.natoms for source in sources ])

        # information of stacks
        new.stacks = dict(zip( classes,
                              [ range(r[0],r[1]) for r in zip([0]+list(accumulate([ source.natoms for source in sources ]))[:-1],
                                                                  list(accumulate([ source.natoms for source in sources ])))]))

        # data are accessed via _master rather than properties
        for field in self._master.dtype.fields:

            # YL 05/09/2020:
            if field not in new._master.dtype.fields:
                continue
            # only deal with common fields in all sources
            # AJL, Aug2017: If I understand correctly, this means we end up with a structure that is common to the lowest
            # possible filetype. When joining an xyz with a punch file, do we just want to keep the name, x, y and z? Shouldn't we aspire to keep all the information
            # we have, and put holders for everything that is unavailable, rather than throw away good data?
            # Definitely in our experience merge does not work with a punch file and xyz, when we want punch file accuracy on the output - in the scenario
            # explained we lose all the shells from e.g. an mgo surface, which is not a good outcome
            if not all([ field in source._master.dtype.fields for source in sources ]):
                continue
           
            try:
                new._master.setField(vstack([ source._master[field] for source in sources ]), field=field)
            except:
                new._master.setField(vstack([ source._master[field].reshape(-1,1) for source in sources ]), field=field)

        return new


# END OF MANIPULATION METHODS #
# SELECTION METHODS #


    # this method is only valid for classes containing `symbols` attribute
    def getRange(self, iter):
        '''Return atom indices for a given (nested) list'''

        from collections import OrderedDict
        from numpy       import where

        # flatten: coreatoms should not be sorted!
        abuff = utils.iterutils.getFlattenedIter(iter, sort=False)

        # translate string species to indices
        # `self.atomtypes` already renamed as `self.symbols`
        try:
            abuff1 = where([ self.symbols == i.title() for i in abuff if type(i) is str ])[1].tolist()
        # if none
        except:
            abuff1 = []

        # find integer (could be numpy.int!) indices
        abuff2 = [ i for i in abuff if utils.objutils.getType(i) is int ]

        # merge the two lists
        abuff1.extend(abuff2)

        # do not use set() which sorts the sequence!
#        return sorted(set(abuff1))
#        return list(set(abuff1))
        return list(OrderedDict.fromkeys(abuff1))


    # YL 15/08/2019: recursively copy all attributes
    def getCopy(self):
        '''Return a copy of the current instance'''

        from copy import deepcopy

        _kwargs = {}

        for key, val in self._kwargs.items():

            attr = getattr(self, key)

            # _XYZ attributes
            if hasattr(attr, 'getCopy'):
                _kwargs.update({key:attr.getCopy()})
            # ndarray attributes
            if hasattr(attr, 'copy'):
                _kwargs.update({key:attr.copy()})
            # other attributes
            try:
                _kwargs.update({key:deepcopy(attr)})
            except:
                continue

        inst = self.__class__()

        try:
            inst.__init__(**_kwargs)
        except:
            pass

        return inst

        


# END OF SELECTION METHODS #


    def rotate(self, alpha=0.0, beta=0.0, gamma=0.0, theta=0.0, vector=[1.0,1.0,1.0]):
        '''Counterclockwise rotate the cartesian coordinates (in degrees):
               alpha - about X-axis
               beta  - about Y-axis
               gamma - about Z-axis
               theta - about vector
        '''

        from numpy import array, cos, dot, float64, linalg, pi, sin, zeros

        if alpha != 0.0:
            alpha = alpha*pi/180.0
            Ax = array([[1.0,0.0,0.0],[0.0,cos(alpha),-sin(alpha)],[0.0,sin(alpha),cos(alpha)]], dtype=float64)
            self._master.xyz = dot(Ax, self._master.xyz.T).T

        if beta  != 0.0:
            beta  = beta *pi/180.0
            Ay = array([[cos(beta ),0.0,sin(beta )],[0.0,1.0,0.0],[-sin(beta ),0.0,cos(beta )]], dtype=float64)
            self._master.xyz = dot(Ay, self._master.xyz.T).T

        if gamma != 0.0:
            gamma = gamma*pi/180.0
            Az = array([[cos(gamma),-sin(gamma),0.0],[sin(gamma),cos(gamma),0.0],[0.0,0.0,1.0]], dtype=float64)
            self._master.xyz = dot(Az, self._master.xyz.T).T
        
        if theta != 0.0:
            theta = theta*pi/180.0

            # normalise the vector
            vector = array(vector, dtype=float64)/linalg.norm(vector)

            M = array([[ cos(theta) + (1-cos(theta))*vector[0]**2,
                        (1-cos(theta))*vector[0]*vector[1] - sin(theta)*vector[2],
                        (1-cos(theta))*vector[0]*vector[2] + sin(theta)*vector[1]],
                       [(1-cos(theta))*vector[1]*vector[0] + sin(theta)*vector[2],
                         cos(theta) + (1-cos(theta))*vector[1]**2,
                        (1-cos(theta))*vector[1]*vector[2] - sin(theta)*vector[0]],
                       [(1-cos(theta))*vector[2]*vector[0] - sin(theta)*vector[1],
                        (1-cos(theta))*vector[2]*vector[1] + sin(theta)*vector[0],
                         cos(theta) + (1-cos(theta))*vector[2]**2]],
                      dtype=float64)
            self._master.xyz = dot(self._master.xyz, M)


    def translate(self, x=0.0, y=0.0, z=0.0, vector=[0.0,0.0,0.0]):
        '''Translate the cartesian coordinates (in a.u.)'''

        from numpy import array

        self._master.xyz[:,0] += x
        self._master.xyz[:,1] += y
        self._master.xyz[:,2] += z

        self._master.xyz += array(vector)


    # aliases
    def coords2txt(self, **kwargs):
        '''Alias to coords2str'''

        return self.coords2str(**kwargs)


# YL 27/07/2022: moved from base.nparray which has been moved to DL_PY2F as utils/nputils.py
# TODO change name
def xyzarray(buffer, npts=0, factor=1.0):
    '''Create an instance of <ndarray> with data contained in buffer and shape it into (natoms,3)'''

    from numpy import array, float64, require, ndarray, full
#    from .errors import ChemShNoneError

#    if buffer == None:
#        raise ValueError
#        raise ChemShNoneError

    # array buffer
    if hasattr(buffer, '__len__'):
        buffer = array(float64(buffer))
        buffer.__imul__(factor)
        if not npts:
            npts = int(array(buffer).size/3)

        return require(ndarray(shape=(npts,3),
                               dtype=float64,
                               buffer=buffer),
                       requirements=['O'])

    # scalar buffer
    else:
        return require(full((npts,3),
                             buffer*factor,
                             dtype=float64),
                       requirements=['O'])
