#  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 absolute, add, all, amax, any, arange, array, count_nonzero, dtype, float64, full
from numpy            import int32, int64, linalg, ma, multiply, nonzero, sign, subtract, where, zeros
from sys              import stdout
from time             import time
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 ..dl_py2f.utils  import miscutils


class _FragmentBase(_XYZ):
    '''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 = {
                'bqs'          : BQs(),
                'cell'         : Cell(),
                'charges'      : zeros(shape=(1,) ,     dtype=float64),
                'conn'         : full((1,MAXCONN) , -1, dtype=int64),
                'connect_mode' :'covalent',
                'connect_scale': 1.0,
                'connect_toler': 0.5,
                'coords'       : zeros(shape=(1,3),     dtype=float64),
# FIXME: bool is only 8-bit that cannot be aligned in Fortran 2D-array
#        so frozen will be passed as int64 for now
                'fragIDs'      : full((1,)        , -1, dtype=int64),
# YL: we may not need this flag (can treat this list as for qm_region)
                'frozen'       : zeros(shape=(1,) ,     dtype=int64),
                'grid'         : None,
                'masses'       : zeros(shape=(1,) ,     dtype=float64),
                'names'        : full((1),         '',  dtype='S8'),
                'natoms'       : 0,
                'nbqs'         : 0,
                'nshells'      : 0,
                # YL 15/04/2018: shellatoms in Fragment is a list of shells indices and coreatoms in Shells is one of core indices
                'trajectory'   : zeros(shape=(1,3),     dtype=float64),  # YL 24/06/2021: added this option to override the coordinates with an MD trajectory/restart file
                'totalcharge'  : 0.0,
                'znums'        : zeros(shape=(1,),      dtype=int64),
              }

    # fields of <utils.nputils.RecArray _master>: the order matters!
    # YL 24/06/2021: `trajectory` should override any previously defined coordinates
    _fields = 'names', 'znums', 'coords', 'charges', 'conn', 'fragIDs', 'frozen', 'masses', 'shellatoms', 'trajectory'

    # make sure 'unit' is processed before 'coords' and
    #           'shellatoms' before 'coords' (connectivity depends on shell charges!)
    # YL 13/09/2020: dependencies are cell -> cell.radius -> frag.coords
    _priorities = 'connect_mode', 'coords', 'cell', 'shellatoms'

    _internals = {
                   '_recalc_fragIDs': True,
                   '_base_cutoff'   : 10.0,                         # a base constant for cutoff of cell lists (see FragmentBase.cell_lists_cutoff)
                   '_nconns_aux'    : MAXCONN,                      # maxinum number of connections (2nd dimension of _conn_aux)
                   '_conn_slots'    : zeros(shape=(1,), dtype=int), # available slots in `_conn_aux`
                   # TWK 24/06/2020: same charge sign scaling value
                   #                 originally from the UCL cluster construct program
                   '_scale_covalent': 0.5
                 }

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

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


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

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

        # YL 27/03/2017: brought it prior to the self._cell.natoms because the reference should always be there
        # cell.cartesian is a reference to frag.coords
        self._cell._cartesian = self.coords
        self._cell._Cell__radius = self.radius

        # selectively return cell
        if self._cell.natoms:

            # find equivalent coordinates
            # YL 02/05/2018: removed use of findCloserSites() here because Cell.frac2cart() will do so instead
            #                and if coordinates are read in from a file they shouldn't be changed
#            self._cell.findCloserSites()

            # cell.connect is a reference to method frag.connect
            self._cell.connect = self.connect

        return self._cell


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

        self._cell = val
        self._cell._cartesian = self.coords


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

        # NOTE: do not sync cell any more
        # if containing cell, sync the coords
        # it should be safe to use self._cell (rather than self.cell) since cell.fractional.setter now assigns cell.cartesian as well
#        if self._cell:
#            self._master.setField(self._cell.cartesian, field='xyz')

        return self._master.xyz


    @coords.setter
    @_XYZ._resize
    def coords(self, val):
        '''Setter of cartesian coordinates'''

        # since July 2016 there is no more 'unit' in the class
        self._master.setField(val, field='xyz')

        # flag of recalculating connectivity by new coords
        if self.connect_mode:
            self._reconnect = True


    @property
    def charges(self):
        '''Charges: (natoms,1)'''

        return self._master.charges


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

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


    @property
    def conn(self):
        '''Connectivity'''

        # check if reconnecting needed
        if self.connect_mode and self._reconnect:
            self.connect()
#            self._reconnect = False

        return self._master.conn


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

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


    @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 symbols(self):
        '''Atom symbols (types)'''

        from numpy   import array
        from ..data  import atomicdata

        return array(atomicdata.symbols).take(self.znums)


    @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 zspecies(self):
        '''Atom species in terms of znum'''

        from numpy             import asarray
        from ..data.atomicdata import species

        # YL 18/06/2023: return asarray([ int(species[s]) for s in self.species ]) causes such warning
        #                DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated      NumPy 1.25.)
        return asarray([ int(species[s][0]) for s in self.species ])


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

        from numpy import unique

        return self.species.size


    @property
    def znums(self):
        '''Atomic numbers: (natoms,1)'''

        from ..base import forcefield
        from ..data import atomicdata

        # if znums no initialised, try to convert from self.names
        # YL 26/03/2019: a bug found here that a name of b'' leads to znum of 89 and out of bound in masses data array
        if not count_nonzero(self._master.znums) and count_nonzero(self._master.names):
            self._master.setField(forcefield.names2znums(self.names), field='znums')
            # YL 19/03/2019: we forgot to set masses!
            self._master.setField(atomicdata.abundant_masses[self._master.znums], field='masses')

        return self._master.znums


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

        from numpy  import array, char, unique
        from ..base import forcefield
        from ..data import atomicdata

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

        # also set masses
        self._master.setField(atomicdata.abundant_masses[self._master.znums], field='masses')

        # also set names
        # AL, 2019: This is a bit unpleasant if we don't want names changed - can we make this
        #           optional somehow? The only way at present is to reset names afterwards
        # YL 13/09/2020: only overwrite names if not yet initialised
        if not any(char.str_len(self.names)):
            self._master.setField(forcefield.znums2names(self._master.znums), field='names')


    @property
    def masses(self):
        '''Abundant masses'''

        return self._master.masses


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

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


    @property
    def formula(self):
        '''Chemical formula'''

        # elements in atomicdata.order (if count is 1, do not print count)
        l1 = [ element+str(self.symbols.tolist().count(element)) if self.symbols.tolist().count(element) > 1 else
               element
               for element in atomicdata.order if element in self.symbols ]
        # elements not in atomicdata.order (if count is 1, do not print count)
        l2 = [ element+str(self.symbols.tolist().count(element)) if self.symbols.tolist().count(element) > 1 else
               element
               for element in set(self.symbols) - set(atomicdata.order) ]

        return ''.join(l2)+''.join(l1)


    @property
    def _conn_aux(self):
        '''Auxiliary connectivity table for containing more than MAXCONN connections, for example computed using a unitary bond distance'''

        try:
            return self.__conn_aux
        except:
            # default size to save memory
            self.__conn_aux = full((self.natoms,1), -1, dtype=int32)
            return self.__conn_aux


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

        self.__conn_aux = val


    @property
    def _conn_slots(self):
        '''Indices of available slots in `_conn_aux` (namely, the lowest position the number remains -1)'''

        try:
            return self.__conn_slots
        except:
            # default size to save memory
            self.__conn_slots = full(self.natoms, 0, dtype=int32)
            return self.__conn_slots


# 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 bqs(self):
        '''Point charges'''

        return self._bqs


    @bqs.setter
    def bqs(self, val):
        '''Setter of point charges'''

        self._bqs = val


    @property
    def nbqs(self):
        '''Number of point charges'''

        return self._bqs._master.shape[0]


    # do nothing
    @nbqs.setter
    def nbqs(self, val):
        '''Setter of nbqs'''


    @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)
    

# 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)

    # do nothing
    @nshells.setter
    def nshells(self, val):
        ''''''


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

        from numpy import add, arange, argwhere, 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
        # YL 29/01/2019: the following `unassociated` buffer are used to protect centres from unassoicated core-shell pairs (where coreatoms is -1)
        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):
                unassociated = self._shells._master.xyz[argwhere(coreatoms==-1)]
                self._shells._master.setField(add(self._shells.displace, self.coords.take(coreatoms, axis=0)), field='xyz')
                self._shells._master.xyz[argwhere(coreatoms==-1)] = unassociated[:]
            # if displacements are just default: keep absolute positions and calculate new displacements
            elif count_nonzero(self._shells.coords):
                unassociated = self._shells._master.displace[argwhere(coreatoms==-1)]
                self._shells._master.setField(subtract(self._shells.coords, self.coords.take(coreatoms, axis=0)), field='displace')
                self._shells._master.displace[argwhere(coreatoms==-1)] = unassociated[:]
            # if frag.shells.coords == displacements (in case shells are not read from a file but defined with a list of coreatoms)
            else:
                unassociated = self._shells._master.xyz[argwhere(coreatoms==-1)]
                self._shells._master.setField(self.coords.take(coreatoms, axis=0), field='xyz')
                self._shells._master.xyz[argwhere(coreatoms==-1)] = unassociated[:]

            # 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 title(self):
#        '''Title of the fragment'''
#
#        return self._title
#
#
#    @title.setter
#    def title(self, val):
#        '''Setter of title'''
#
#        if type(val) is str:
#            self._title = val
#        else:
#            self._title = ' '.join(val)

    @property
    def net_charge(self):
        '''Net total charge'''

        return self.charges.sum()


    @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'''

        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 #


    def recentre(self):
        '''Recentre the fragment to its centroid'''

        self.coords -= self.centroid


    def reorder(self, order):
        '''Wrapper function of _XYZ.reorder() with reconnect'''

        from .xyz import _XYZ

        _XYZ.reorder(self, order)

        # YL 29/11/2021: __covalent_radii must be reset!
        try:
            delattr(self, '_'+self.__class__.__name__+'__covalent_radii')
        except:
            pass

        # reconnect!
        self.connect(force_reconnect=True)


    def changeConnectOptions(self, connect_scale=None, connect_toler=None, connect_mode=""):
        '''Set new connectivity test parameters and rerun connectivity calc.'''

        # Set new parameters as requested
        if connect_scale is not None:
            self.connect_scale = connect_scale
        if connect_toler is not None:
            self.connect_toler = connect_toler
        # NB None is a valid option for connect_mode so test for empty string instead
        if connect_mode != '':
            self.connect_mode = connect_mode

        # Force compute of connecitivity
        self._reconnect = True
        self.connect()

   
    def connect(self, force_reconnect=False, mute=True):
        '''Compute the connectivity'''

        from ..base  import connectivity

        starttime = time()

        if force_reconnect:
            self._reconnect = True

        multiproc = False

        if self.natoms > 10000 and not mute:
            print("\n Calculating connectivity table...")
            multiproc = True
            stdout.flush()

        if self.connect_mode and self._reconnect:

            # YL 29/11/2021: __covalent_radii must be reset before connect!
            try:
                delattr(self, '_'+self.__class__.__name__+'__covalent_radii')
            except:
                pass

            self._master.conn = -1
            connectivity.connect(self, multiproc=multiproc)

            # all dependant properties must be wiped out
            try:
                del self.__dihedrals
                del self.__angles
                del self.__bonds
                del self.__vdwtypes
            except AttributeError:
                pass

        if self.connect_mode and self._reconnect:
            if not mute or self.natoms > 10000:
                print("\n >>> Time used for calculating connectivity: %.3f s"%(time()-starttime))
                stdout.flush()

        # normally, no more need to calculate connectivity once it is done 
        self._reconnect = False


    @property
    def _max_coords(self):
        '''Greatest values of coords in the format of [[x_min,x_max],[y_min,y_max],[z_min,z_max]]'''

        from numpy import array

        _max_coords = []

        for i in range(3):
            _max_coords.append([self.coords[:,i].min(),self.coords[:,i].max()])

        return array(_max_coords)


    @property
    def _max_abs_coords(self):
        '''Greatest absolute values of coords'''

        from numpy import array, fabs

        _max_abs_coords = []

        for i in range(3):
            _max_abs_coords.append(fabs(self.coords)[:,i].max())

        return array(_max_abs_coords)


    @property
    def _xyz_spanning(self):
        '''Spanning distance along X, Y, and Z, respectively'''

        from numpy import array, fabs

        _xyz_spanning = []

        for i in range(3):
            _xyz_spanning.append(fabs(self.coords[:,i].max()-self.coords[:,i].min()))

        return array(_xyz_spanning)


    # YL 12/03/2021: a vectorised version of connectIJ(), faster than the previous implementation of connectIJ()
    def connectIJs(self, iatoms, jatoms):
        '''Connect iatoms to jatoms: the both must have a same shape to perform one-on-one mapping'''

        from numpy import array, unique

        # find the lowest index of avalable slot (where it's still -1) for each I
        slots = self._master.conn[iatoms].argmin(axis=1)

        # masks array of where atom J already in I's conn list
        masks_j_not_in_i = all(self._master.conn[iatoms]-jatoms[:,None]!=0, axis=1)

        # must apply increments to the slot positions because there may be duplicates in iatoms!
        increments = utils.nputils.getIncrementsFromCounts(iatoms[masks_j_not_in_i])

        # must use tuples for fancy indexing otherwise: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated
        # YL 25/08/2021: we can't connect to more than MAXCONN (14) centres
        positions = tuple(iatoms[masks_j_not_in_i]), tuple(slots[masks_j_not_in_i]+increments)

        # assign filtered values to the slots
        # NB: we can't use array[positions][masks|indices] = values because the left bound doesn't return a reference to the original array!
        try:
            self._master.conn[positions] = jatoms[masks_j_not_in_i]
        # YL 26/05/2024: added debugging print
        except IndexError:
            print(' >>> ERROR: FragmentBase.connectIJs() exceeds max number of conns (14) when')
            masks_exceeding = array(positions[1]) > 13
            iatoms_exceeding = unique(iatoms[masks_exceeding])
            for iatom_exceeding in iatoms_exceeding:
                print(f'            connecting {iatom_exceeding} ({self.names[iatom_exceeding].decode()}) to', flush=True)
                jatoms_connected_to_i = jatoms[array(positions[0])==iatom_exceeding]
                print('            ', jatoms_connected_to_i, flush=True)
            raise


    # YL 26/08/2021: when using a unitary distance which may break MAXCONN we
    #                must know the maximum number of connections beforehand
    # YL 23/07/2023: the previously implemented algorithm wasn't valid because self._conn_aux didn't 
    #                exit and there was no way to exclude already-connect jatoms using masks_j_not_in_i;
    #                however, the current modified method can guarantee self._nconns_aux is greater than
    #                needed
    def _updateNConnsUnitaryDist(self, iatoms, jatoms):
        '''Prepare for connecting iatoms to jatoms: the both must have a same shape to perform one-on-one mapping'''

        # find the lowest index of avalable slot (where it's still -1) for each I
#        slots = self._conn_aux[iatoms].argmin(axis=1)

        # masks array of where atom J already in I's conn list
#        masks_j_not_in_i = all(self._conn_aux[iatoms]-jatoms[:,None]!=0, axis=1)

        # must apply increments to the slot positions because there may be duplicates in iatoms!
        increments = utils.nputils.getIncrementsFromCounts(iatoms)

        # must use tuples for fancy indexing otherwise: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated
        # YL 25/08/2021: we can't connect to more than MAXCONN (14) centres
#        positions = tuple(iatoms[masks_j_not_in_i]), tuple(slots[masks_j_not_in_i]+increments)
        positions = tuple(iatoms), tuple(self._conn_slots[iatoms]+increments)
        # YL 22/07/2023: we must use the put() function otherwise the values don't change
        self._conn_slots.put(iatoms, self._conn_slots[iatoms]+increments)

        # we don't connect but only get the max number of connections
        self._nconns_aux = max(amax(positions[1])+1, self._nconns_aux)


    # YL 26/08/2021: 
    def connectIJsUnitaryDist(self, iatoms, jatoms):
        '''Connect iatoms to jatoms: the both must have a same shape to perform one-on-one mapping'''

        # find the lowest index of avalable slot (where it's still -1) for each I
        slots = self._conn_aux[iatoms].argmin(axis=1)

        # masks array of where atom J already in I's conn list
        masks_j_not_in_i = all(self._conn_aux[iatoms]-jatoms[:,None]!=0, axis=1)

        # must apply increments to the slot positions because there may be duplicates in iatoms!
        increments = utils.nputils.getIncrementsFromCounts(iatoms[masks_j_not_in_i])

        # must use tuples for fancy indexing otherwise: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated
        # YL 25/08/2021: we can't connect to more than MAXCONN (14) centres
        positions = tuple(iatoms[masks_j_not_in_i]), tuple(slots[masks_j_not_in_i]+increments)

        # assign filtered values to the slots
        # NB: we can't use array[positions][masks|indices] = values because the left bound doesn't return a reference to the original array!
        self._conn_aux[positions] = jatoms[masks_j_not_in_i]


    def connectIJ(self, iatom, jatom):
        '''Connect jatom to iatom'''
   
        def _connectI2J(_iatom, _jatom):
            # find position of next slot
            _idx = self._master.conn[_iatom].argmin()
            # set connections
            if _jatom not in self._master.conn[_iatom]:
                self._master.conn[_iatom][_idx] = _jatom
   
        _connectI2J(iatom, jatom)
        _connectI2J(jatom, iatom)
   
   
    def connPairs(self, pairs):
        '''Connect atom pairs provided in an array'''

        from numpy   import array

        # YL 03/09/2020: reshape anyway otherwise a fragment's conn cannot be read in that contains only two atoms
        abuff = array(pairs)
        abuff = abuff.reshape((int(abuff.size/2),2))
        for pair in abuff:
            self.connectIJ(int(pair[0]), int(pair[1]))


    def sortConnection(self):
        # sort self.conn in place: e.g., to [1,3,6,-1,-1,-1]
        for conn in self.conn:
            conn[:conn.argmin()].sort()


    def save(self, filename='', fmt='', format='', mute=False, **kwargs):
        '''Write out fragment data to disk using specified file format.'''

        from ..io.file import dump

        if not filename:
            filename = self.title.rstrip() + '.pun'

        if format:
            fmt = format

        dump(filename, fmt, self, mute=mute, **kwargs)


    def _updateNConnsUnitaryDistIIs(self, iatoms, unitary_dist=0.0):
        '''Calculate the connectivity within subset iatoms by a unitary distance'''

        # no need to proceed
        if not len(iatoms) or not unitary_dist:
            return

        # dimensions: (I,I,XYZ)
        # we construct arrays of all data using NumPy instead of using Python for loops which are slower (see the previous implementation in base/connectivity.py)
        _vectors = self.coords[iatoms][:,None,:].__rsub__(self.coords[iatoms])

        # positions of qualified atoms in array _cutoffs: exclude I-I
        pos_i, pos_j = where((linalg.norm(_vectors, axis=2) > 0.0)*(linalg.norm(_vectors, axis=2) <= unitary_dist))

        # no need to proceed
        if not len(pos_i)*len(pos_j):
            return

        # get the indices
        indices_i, indices_j = iatoms[pos_i], iatoms[pos_j]

        # connect
        self._updateNConnsUnitaryDist(indices_i, indices_j)
        # also do it vice sersa
        self._updateNConnsUnitaryDist(indices_j, indices_i)


    def _updateNConnsUnitaryDistIJs(self, iatoms, jatoms, unitary_dist=0.0):
        '''Calculate the connectivity between subset iatoms and jatoms by a unitary distance'''

        # no need to proceed
        if not len(iatoms)*len(jatoms) or not unitary_dist:
            return

        # dimensions: (I,J,XYZ)
        # we construct arrays of all data using NumPy instead of using Python for loops which are slower (see the previous implementation in base/connectivity.py)
        _vectors   = self.coords[iatoms][:,None,:].__rsub__(self.coords[jatoms])

        # positions of qualified atoms in array _cutoffs
        pos_i, pos_j = where(linalg.norm(_vectors, axis=2) <= unitary_dist)

        # no need to proceed
        if not len(pos_i)*len(pos_j):
            return

        # get the indices
        indices_i, indices_j = iatoms[pos_i], jatoms[pos_j]

        # connect
        self._updateNConnsUnitaryDist(indices_i, indices_j)
        # also do it vice sersa
        self._updateNConnsUnitaryDist(indices_j, indices_i)


    def _connectUnitaryDistIIs(self, iatoms, unitary_dist=0.0):
        '''Calculate the connectivity within subset iatoms by a unitary distance'''

        # no need to proceed
        if not len(iatoms) or not unitary_dist:
            return

        # dimensions: (I,I,XYZ)
        # we construct arrays of all data using NumPy instead of using Python for loops which are slower (see the previous implementation in base/connectivity.py)
        _vectors = self.coords[iatoms][:,None,:].__rsub__(self.coords[iatoms])

        # positions of qualified atoms in array _cutoffs: exclude I-I
        pos_i, pos_j = where((linalg.norm(_vectors, axis=2) > 0.0)*(linalg.norm(_vectors, axis=2) <= unitary_dist))

        # no need to proceed
        if not len(pos_i)*len(pos_j):
            return

        # get the indices
        indices_i, indices_j = iatoms[pos_i], iatoms[pos_j]

        # connect
        self.connectIJsUnitaryDist(indices_i, indices_j)
        # also do it vice sersa
        self.connectIJsUnitaryDist(indices_j, indices_i)


    def _connectUnitaryDistIJs(self, iatoms, jatoms, unitary_dist=0.0):
        '''Calculate the connectivity between subset iatoms and jatoms by a unitary distance'''

        # no need to proceed
        if not len(iatoms)*len(jatoms) or not unitary_dist:
            return

        # dimensions: (I,J,XYZ)
        # we construct arrays of all data using NumPy instead of using Python for loops which are slower (see the previous implementation in base/connectivity.py)
        _vectors   = self.coords[iatoms][:,None,:].__rsub__(self.coords[jatoms])

        # positions of qualified atoms in array _cutoffs
        pos_i, pos_j = where(linalg.norm(_vectors, axis=2) <= unitary_dist)

        # no need to proceed
        if not len(pos_i)*len(pos_j):
            return

        # get the indices
        indices_i, indices_j = iatoms[pos_i], jatoms[pos_j]

        # connect
        self.connectIJsUnitaryDist(indices_i, indices_j)
        # also do it vice sersa
        self.connectIJsUnitaryDist(indices_j, indices_i)


    def _connectCovalentIIs(self, iatoms, **kwargs):
        '''Calculate the covalent connectivity within subset iatoms'''

        # no need to proceed
        if not len(iatoms):
            return

#        radii  = atomicdata.covalent_radii

        # dimensions: (I,I,XYZ)
        # we construct arrays of all data using NumPy instead of using Python for loops which are slower (see the previous implementation in base/connectivity.py)
        coords = self.coords[iatoms]
        _vectors = coords[:,None,:].__rsub__(coords)

        # get scale and tolerance for determining connectivity according to the formula Rij < scale * (radius_i + radius_j) + toler
        radii = self._covalent_radii[iatoms]
        _cutoffs = radii[:,None].__radd__(radii)*self.connect_scale + self.connect_toler

        # positions of qualified atoms in array _cutoffs: exclude I-I
        pos_i, pos_j = where((linalg.norm(_vectors, axis=2) > 0.0)*(linalg.norm(_vectors, axis=2) <= _cutoffs))

        # no need to proceed
        if not len(pos_i)*len(pos_j):
            return

        # get the indices
        indices_i, indices_j = iatoms[pos_i], iatoms[pos_j]

        # connect
        self.connectIJs(indices_i, indices_j)
        # also do it vice sersa
        self.connectIJs(indices_j, indices_i)


    #  YL 12/03/2021: called by base.connectivity.connectNonperiodicCellLists()
    #                 iatoms and jatoms don't have to have the same size
    def _connectCovalentIJs(self, iatoms, jatoms, **kwargs):
        '''Calculate the covalent connectivity between subset iatoms and jatoms'''

        # no need to proceed
        if not len(iatoms)*len(jatoms):
            return

#        radii = atomicdata.covalent_radii

        # dimensions: (I,J,XYZ)
        # we construct arrays of all data using NumPy instead of using Python for loops which are slower (see the previous implementation in base/connectivity.py)
        _vectors = self.coords[iatoms][:,None,:].__rsub__(self.coords[jatoms])

        # get scale and tolerance for determining connectivity according to the formula Rij < scale * (radius_i + radius_j) + toler
        _cutoffs = self._covalent_radii[iatoms][:,None].__radd__(self._covalent_radii[jatoms])*self.connect_scale + self.connect_toler
        # positions of qualified atoms in array _cutoffs
        pos_i, pos_j = where(linalg.norm(_vectors, axis=2) <= _cutoffs)

        # no need to proceed
        if not len(pos_i)*len(pos_j):
            return

        # get the indices
        indices_i, indices_j = iatoms[pos_i], jatoms[pos_j]

        # connect
        self.connectIJs(indices_i, indices_j)
        # also do it vice sersa
        self.connectIJs(indices_j, indices_i)


    def _connectIonicIIs(self, iatoms, **kwargs):
        '''Calculate the ionic connectivity within subset iatoms'''

        # no need to proceed
        if not len(iatoms):
            return

        echarges = self.effectcharges
        radii    = atomicdata.ionic_radii
        natoms   = len(iatoms)**2

        _vectors  = self.coords[iatoms][:,None,:].__rsub__(self.coords[iatoms])
# YL NB 12/03/2021: since we currently only use the current method for performing the cell lists algorithm, there is
#                   no longer need to do getMinImages(); but if the developer use the method for Verlet list or other
#                   algorithms the getMinImages line must be uncommented!
        # unfold the rank 3 array of vectors to be used by connectivity.getMinImages()
        # PBC on mockcell will be '3D' although we currently only run cell lists algorithm on non-periodic systems (see base.connectivity.connect())
#        _vectors  = connectivity.getMinImages(_vectors.reshape(natoms,3), self.mockcell).reshape(len(iatoms),len(iatoms),3)

        _sumradii = radii[self.znums[iatoms]][:,None].__radd__(radii[self.znums[iatoms]])
        # get scale and tolerance for determining connectivity according to the formula Rij < scale * (radius_i + radius_j) + toler
        _cutoffs  = _sumradii*self.connect_scale + self.connect_toler

        icharges  = echarges[iatoms]
        # TWK 24/06/2020: For ions of the same sign the distance cut-off is scaled down
        #                 to exclude "unintentional" coordination. Only covalent pairing
        #                 is thus allowed as in O22- etc.
        masks_same_sign = (sign(icharges[:,None])==-1)==(sign(icharges)==-1)
        _cutoffs[masks_same_sign] -= _sumradii[masks_same_sign]*self._scale_covalent

        # positions of qualified atoms in array _cutoffs: exclude I-I
        pos_i, pos_j = where((linalg.norm(_vectors, axis=2) > 0.0)*(linalg.norm(_vectors, axis=2) <= _cutoffs))

        # no need to proceed
        if not len(pos_i)*len(pos_j):
            return

        # get the indices
        indices_i, indices_j = iatoms[pos_i], iatoms[pos_j]

        # connect
        self.connectIJs(indices_i, indices_j)
        # also do it vice sersa
        self.connectIJs(indices_j, indices_i)


    #  YL 12/03/2021: called by base.connectivity.connectNonperiodicCellLists()
    #                 iatoms and jatoms don't have to have the same size
    def _connectIonicIJs(self, iatoms, jatoms, **kwargs):
        '''Calculate the ionic connectivity within subset iatoms'''

        # no need to proceed
        if not len(iatoms)*len(jatoms):
            return

        echarges = self.effectcharges
        radii    = atomicdata.ionic_radii
        natoms   = len(iatoms)*len(jatoms)

        _vectors = self.coords[iatoms][:,None,:].__rsub__(self.coords[jatoms])
# YL NB 12/03/2021: since we currently only use the current method for performing the cell lists algorithm, there is
#                   no longer need to do getMinImages(); but if the developer use the method for Verlet list or other
#                   algorithms the getMinImages line must be uncommented!
        # unfold the rank 3 array to be used by connectivity.getMinImages()
        # PBC on mockcell will be '3D' although we currently only run cell lists algorithm on non-periodic systems (see base.connectivity.connect())
#        _vectors = connectivity.getMinImages(_vectors.reshape(natoms,3), self.mockcell).reshape(len(iatoms),len(jatoms),3)

        _sumradii = radii[self.znums[iatoms]][:,None].__radd__(radii[self.znums[jatoms]])
        # get scale and tolerance for determining connectivity according to the formula Rij < scale * (radius_i + radius_j) + toler
        _cutoffs  = _sumradii*self.connect_scale + self.connect_toler

        icharges  = echarges[iatoms]
        jcharges  = echarges[jatoms]
        # TWK 24/06/2020: For ions of the same sign the distance cut-off is scaled down
        #                 to exclude "unintentional" coordination. Only covalent pairing
        #                 is thus allowed as in O22- etc.
        masks_same_sign = (sign(icharges[:,None])==-1)==(sign(jcharges)==-1)
        _cutoffs[masks_same_sign] -= _sumradii[masks_same_sign]*self._scale_covalent

        # positions of qualified atoms in array _cutoffs: exclude I-I
        pos_i, pos_j = where((linalg.norm(_vectors, axis=2) > 0.0)*(linalg.norm(_vectors, axis=2) <= _cutoffs))

        # no need to proceed
        if not len(pos_i)*len(pos_j):
            return

        # get the indices
        indices_i, indices_j = iatoms[pos_i], jatoms[pos_j]

        # connect
        self.connectIJs(indices_i, indices_j)
        # also do it vice sersa
        self.connectIJs(indices_j, indices_i)


# END OF CLASSMETHODS #
# CLUSTER CUTTING METHODS #


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

        from .. import cluster

        return cluster.adjustCharge(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)


    # TODO: to move to the proposed <class Cluster>
    def setRegion(self, list_of_atoms, region=1):
        '''Set atom names as XX? (region ?)'''

        from string import digits

        for i in list_of_atoms:

            # remove any trailing numbers
            self.names[i] = self.names[i].decode().rstrip(digits).encode()

            self.names[i] += str(region).encode()


    def unsetRegion(self, list_of_atoms):
        '''Unset atom names from XX? to XX (region ?)'''

        from string import digits

        for i in list_of_atoms:

            # remove any trailing numbers
            self.names[i] = self.names[i].decode().rstrip(digits).encode()


# END OF CLUSTER CUTTING METHODS #
# INQUIRIES #


    @property
    def centroid(self):
        '''Centroid of the fragment'''

        from numpy import array, mean

        # YL 13/09/2020: `natoms` as a property will be processed by utils.objutils.init() with 0 resulting in
        #                zero-sized fragment (it is sometimes 1 if _resize1D() is triggered by another property)
        if self.natoms:
            return mean(self.coords, axis=0)
        else:
            return array([0.0, 0.0, 0.0])


    @property
    def radius(self):
        '''Radius of the smallest sphere which, when centred at the centroid of the fragment, encloses all the atoms'''

        from numpy import linalg, subtract

        # YL 13/09/2020: `natoms` as a property will be processed by utils.objutils.init() with 0 resulting in
        #                zero-sized fragment (it is sometimes 1 if _resize1D() is triggered by another property)
        if self.natoms:
            return max(linalg.norm(subtract(self.coords, self.centroid), axis=1))
        else:
            return 0.0


    # YL 01/02/2021: moved from base/connectivity.py
    # see linkcells2() in tcl-chemsh/src/fragment/connect.c
    @property
    def cell_lists_cutoff(self):
        '''Cutoff for connectivity'''

        from ..data import atomicdata

        # avoid recalculating
        try:
            return self._cell_lists_cutoff

        except:
            selectcases = { 'covalent': atomicdata.covalent_radii,
                            'ionic'   : atomicdata.ionic_radii
                          }
       
            scale = self.connect_scale
            toler = self.connect_toler
    
            # r_cutoff = max(10.0, 2.0*scale*r_max + toler)
            try:
                max_atomic_radius = selectcases[self.connect_mode.strip().lower()][self.zspecies].max()
            except KeyError:
                raise errors.ChemShValError('connect_mode', selectcases.keys())
       
            self._cell_lists_cutoff = max(self._base_cutoff, 2.0*scale*max_atomic_radius + toler)

            return self._cell_lists_cutoff


    @property
    def mockcell(self):
        '''Get a mock cell'''

        # avoid recalculating
        try:
            return self._mockcell

        except:
            # make a mock cell
            self._mockcell = utils.objutils.getDeepCopy(self.cell)

            # apply periodic boundary condition
            self._mockcell.cell_lists_cutoff = self.cell_lists_cutoff

            self._mockcell.applyConnectPBC(version='Tcl-ChemShell')

            return self._mockcell


    def getClosest(self, i, exclude=[], select=[], return_dist=False):
        '''Return the closest atom'''

        from numpy   import arange, argmin, argsort, delete, isin, integer, linalg, subtract

        allatoms = arange(self.natoms)

        # select all by default
        if select == []:
            select = allatoms

        # excluding atom I
        if isinstance(i, int) or isinstance(i, integer):
            exclude = utils.iterutils.getFlattenedIter([i, exclude])

        # apply the selection and exclusion rules
        selected = isin(arange(self.natoms), select) & isin(arange(self.natoms), exclude, invert=True)

        # get sorted indices in ascending order of distance
        if isinstance(i, int) or isinstance(i, integer):
            dists = linalg.norm(subtract(self.coords, self.coords[i]), axis=1)
        # YL 05/09/2020: allow i to be coords [...,...,...]
        else:
            dists = linalg.norm(subtract(self.coords, i), axis=1)
        rankings = argsort(dists)

        # the closest atom is the 1st in the rankings with all unselected and excluded removed
        closest = rankings[isin(rankings, allatoms[selected])][0]

        if chemsh.CHEMSH_DEBUG > 4:
            nselect  = len(select)
            nexclude = len(exclude)
            if nexclude > 1:
                print("\n >>> getClosest(): found the closest atom {} at {} a.u., selecting\n{}{} and excluding {} atoms\n".format(closest, dists[closest], ' '*19, nselect, nexclude))
            else:
                print("\n >>> getClosest(): found the closest atom {} at {} a.u., selecting\n{}{} atoms\n".format(closest, dists[closest], ' '*19, nselect))

        if return_dist:
            return closest, dists[closest]

        return closest


    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


## YL FIXME: to delete--this algorithm is terribly memory-consuming!
#    @property
#    def dihedrals2(self):
#        '''Generate a list containing all dihedrals'''
#
#        from numpy     import any, all, asarray, array, extract, ceil, copy, floor, fix, meshgrid, inf, intersect1d, fromiter, nditer, ndarray, unique, where
#        from numpy     import count_nonzero, argmin, flatnonzero, remainder, modf, trunc, int16, hstack, vstack, in1d, isin, set_printoptions, equal, not_equal, vectorize
#        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(1):
#            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


    @property
    def bonds(self):
        '''Generate a list containing all bonds (excluding duplicates such as I-J and J-I)'''

        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 int(count_nonzero(conn+1)/2)
        return self.bonds.shape[0]


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

        from numpy  import array, copy, count_nonzero, logical_and, unique
        from ..base import constants

        try:
            return self.__angles
        except:
            lbuff = []

        # get a handle
        conn = self.conn

        # 1st atom: i
        for i in range(self.natoms):
            # 2nd atom: j
            for j in conn[i][conn[i]>=0]:
                # 3nd atom: k, should not be i
                for k in conn[j][logical_and(conn[j]>=0, conn[j]!=i)]:
                    lbuff.append([i,j,k])

        # make a copy and sort in place
        ijk = array(lbuff)
        sorted = copy(ijk)
        sorted.sort(axis=1)

        # get the indices for reconstructing ijk
        # 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.__angles = unique(ijk[non_repeateds], axis=0)

        return self.__angles


    @property
    def nangles(self):
        '''Number of angles'''

        return self.angles.shape[0]


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

        from numpy  import array, copy, count_nonzero, logical_and, unique
        from ..base import constants

        try:
            return self.__dihedrals
        except:
            lbuff = []

        # get a handle
        conn = self.conn

        # 1st atom: i
        for i in range(self.natoms):
            # 2nd atom: j
            for j in conn[i][conn[i]>=0]:
                # 3nd atom: k, should not be i
                for k in conn[j][logical_and(conn[j]>=0, conn[j]!=i)]:
                    # 4th atom: l, should not be i or j
                    for l in conn[k][logical_and(logical_and(conn[k]>=0, conn[k]!=i), conn[k]!=j)]:
                        lbuff.append([i,j,k,l])

        # make a copy and sort in place
        ijkl = array(lbuff)
        sorted = copy(ijkl)
        sorted.sort(axis=1)

        # get the indices for reconstructing ijkl
        # 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.__dihedrals = unique(ijkl[non_repeateds], axis=0)

        return self.__dihedrals


    @property
    def ndihedrals(self):
        '''Number of dihedral bonds'''

        return self.dihedrals.shape[0]


    @property
    def dihedral_fftypes(self):
        '''Types of dihedral force field'''

        from numpy import unique

        try:
            # YL NB: axis not available utill NumPy v1.13!
            return unique(self.types[self.__dihedrals], axis=0)
        except:
            return unique(self.types[self.dihedrals], axis=0)


    @property
    def dihedral_ztypes(self):
        '''Types of dihedral force field'''

        from numpy import unique

        try:
            # YL NB: axis not available utill NumPy v1.13!
            return unique(self.znums[self.__dihedrals], axis=0)
        except:
            return unique(self.znums[self.dihedrals], axis=0)


    @property
    def nvdwtypes(self):
        '''Number of types of van der Waals interaction pair'''

        return int((1+self.ntypes)*self.ntypes/2)


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

        from numpy     import array, copy, unique
        from itertools import product

        try:
            return self.__vdwtypes
        except:
            pass

        # what is the dtype if we use fromiter()
        ijtypes = array(list(product(unique(self.types.astype(str, copy=False)),
                                     unique(self.types.astype(str, copy=False)))))

        # make a copy and sort in place
        sorted = copy(ijtypes)
        sorted.sort(axis=1)

        # get the indices for reconstructing ijkl
        # 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.__vdwtypes = unique(ijtypes[non_repeateds], axis=0).astype('U8', copy=False)

        return self.__vdwtypes


    # YL 14/12/2021: improved to make it 5+x faster
    # chain indices
    @property
    def fragIDs(self):
        '''Array of fragment IDs those are not connected to each other via bonds (calculated property rather than assigned from structure file)'''

        from numpy import all, any, argmax, argmin, array, asfortranarray, full, hstack, isin
        from sys   import getrecursionlimit, setrecursionlimit, stdout

        t0 = time()

        # YL 13/12/2021: revised by using asfortranarray() to accelerate the current function by ~25%
        fragIDs = asfortranarray(self._master['fragIDs'])

        # no recalc
        if not self._recalc_fragIDs and all(fragIDs!=-1):
            return self._master['fragIDs']

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

        # YL 30/01/2022: avoid using `.` in a loop body
        conn = self.conn
        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]])
            # YL 14/12/2021: using self.conn is 10x faster than the original implementation using bonds
            js = conn[_i][:argmin(conn[_i])]
            # only calculate undetermined J atoms
            for j in js[fragIDs[js]==-1]:
                # assign for current J
                fragIDs[j] = residue
                # recurse
                findChain(j, residue=residue)

            return

        fragID = 0
        # loop until all atoms determined
        while -1 in fragIDs:

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

        # copy values
        self._master['fragIDs'][:] = fragIDs[:]

        # toggle recalc
        self._recalc_fragIDs = False

        if self.natoms >= 10000:
            print(f'\n >>> Time used for calculating non-bond-connected fragments: {time()-t0:.3f} s')
            stdout.flush()

        return self._master['fragIDs']


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

        self._master.setField(val, field='fragIDs')
        self._recalc_fragIDs = False


    def getIndicesOfFragID(self, fragID):
        '''Return the indices of a given fragID'''

        return self._indices_all[self.fragIDs == fragID]


    def getSelectedFragIDs(self, fragIDs):
        '''Return a new frag of the given fragID'''

        from numpy import zeros

        # fragIDs can be an iterable series
        try:
            masks = zeros((self.natoms,), dtype=bool)
            for fragID in fragIDs:
                masks += self.fragIDs == fragID
            indices = self._indices_all[masks]
            return self.getSelected(indices)

        # or a singular
        except:
            return self.getSelected(self.getIndicesOfFragID(fragIDs))


    @property
    def nfragments(self):
        '''Number of unique fragIDs'''

        from numpy import unique

        return unique(self.fragIDs).size


# END OF INQUIRIES #
# SELECTION METHODS #


    # YL 15/10/2019: new from branch issue45
    def select(self, radius=0.0, atoms='', indices=[], regions=[], return_masks=False, invert=False, union=False, **kwargs):
        '''Select atoms by radius, names/znums, indices, or regions'''

        from numpy     import arange, array, ones, concatenate, ndarray, unique, intersect1d, prod, zeros
        from functools import reduce
        from inspect   import getcallargs, getfullargspec


        # force to return mask arrays because we only use booleans before returning the final result
        if radius == 0.0:
            masks_radius = zeros(shape=(self.natoms,), dtype=bool)
        else:
            masks_radius = self.selectByRadius(radius, return_masks=True, **kwargs)
        masks_atoms   = self.selectByAtoms  (atoms   , return_masks=True, **kwargs)
        masks_indices = self.selectByIndices(indices , return_masks=True)
        masks_regions = self.selectByRegions(regions , return_masks=True)

        # locals() in a generator means that within the latter's domain only
        dlocals = locals()
 
        # check if XYZ mode
        if any({'min_x','min_y','min_z','max_x','max_y','max_z'}.intersection(set(kwargs.keys()))):
            dlocals.update({'xyz':True})
        # YL 04/12/2020: we don't do this unless necessary because selectByXYZ could be relatively slow
        if dlocals.get('xyz', False):
            masks_xyz = self.selectByXYZ(**kwargs, return_masks=True)
        else:
            masks_xyz = zeros(shape=(self.natoms,), dtype=bool)

        # remove exclusive keywords (we need selectBio to issue a warning if a keyword is invalid)
        # save strict which is also needed by selectBio() below
        strict = kwargs.get('strict', True)
        for kw in getfullargspec(self.selectByAtoms).kwonlyargs&kwargs.keys():
            kwargs.pop(kw, None)

        # selectBio is a class method of <class FragBio>
        # NB: invert has to be switched off because the current method will check `invert` anyway just before returning
        # TODO: check if union works for selectBio()
        try:
            masks_bio = self.selectBio(**kwargs, strict=strict, return_masks=True, invert=False, union=union)
            bio_off = False
        # if the current class is not <class FragmentBio>
        except AttributeError:
            # defer because it depends on `union` which will be determined later
            bio_off = True

        #     selected_bio = self.selectBio(**kwargs, return_masks=return_masks)
        #     print(" >>> ERROR: Option \"segments\" is not a valid method of {}".format(self.__class__))
 
        # effective keyword arguments of FragmentBio.selectBio()
        try:
            biokwargs  = getcallargs(self.selectBio, **kwargs)
            thiskwargs = getfullargspec(self.select).args + ['kwargs']
            # remove those belonging to this _FragmentBase.select()
            for arg in thiskwargs:
                biokwargs.pop(arg, None)
            if any(set(biokwargs).intersection(set(kwargs.keys()))):
                dlocals.update({'bio':True})
        # if the current class is not <class FragmentBio>
        except AttributeError:
            pass

        # how many methods of selecting
        nways = sum([radius!=0.0,
                     atoms != '',
                     dlocals.get('xyz', False),
                     dlocals.get('bio', False)])
        # YL 16/02/2021: we can't use indices != [] otherwise:
        #                DeprecationWarning: elementwise != comparison failed; this
        #                will raise an error in the future.
        if not (type(indices) is list and not len(indices)):
            nways += 1
        if not (type(regions) is list and not len(regions)):
            nways += 1

        # enforce `union` if only one method
        if type(nways) is int and nways == 1:
            union = True

        # deferred masks_bio: has to be set after nways is calculated
        # YL 23/02/2021: I have fixed this before but it was messed up by Git auto-merging
        if bio_off:
            if union:
                masks_bio = zeros(shape=(self.natoms,), dtype=bool)
            else:
                masks_bio = ones(shape=(self.natoms,), dtype=bool)

        selectcases = { 'radius' : masks_radius,
                        'atoms'  : masks_atoms,
                        'indices': masks_indices,
                        'regions': masks_regions,
                        'xyz'    : masks_xyz,
                        'bio'    : masks_bio,
                      }

        if union:
            masks = masks_radius + masks_atoms + masks_indices + masks_regions + masks_xyz + masks_bio
        else:
            # take only those chosen as parameters
            if chemsh.CHEMSH_DEBUG > 5:
                print(" >>> {}.select(): all_arrays = {}".format(str(__class__), [k for k in selectcases.keys() if dlocals.get(k, None)]))
            # YL 04/09/2020: can't use a nested generator due to complicated types
            all_arrays = []
            for k in selectcases.keys():
                argument = dlocals.get(k, None)
                if type(argument) is ndarray:
                    if argument.size:
                        all_arrays.append(selectcases[k])
                else:
                    if argument:
                        all_arrays.append(selectcases[k])
            # a kwarg must appear in selectcases.keys() and dlocals.keys() at the same time
#            all_arrays = [ selectcases[k] for k in selectcases.keys() if dlocals.get(k, None) ]
            # exclude any empty option when intersections (union=False) of conditions
            masks = prod(all_arrays, axis=0, dtype=bool)
            # YL 08/07/2020: if no valid keyword given
            if not len(all_arrays):
                masks = zeros((self.natoms,), dtype=bool)

        if invert:
            masks = ~masks

        if return_masks:
            return masks
        else:
            return arange(self.natoms)[masks]

#        return selected


    # YL 05/08/2021: changed the default behaviour of strict to True
    # Yl 26/05/2024: added names_only
    def selectByAtoms(self, *args, names_only=False, return_masks=False, strict=True, **kwargs):
        ''''''

        from numpy  import arange, char, isin, issubdtype, ndarray, str_, zeros
        from ..base import forcefield

        # we use a bool array as data buffer
        selected = zeros(shape=(self.natoms,), dtype=bool)

        def _parseArg(_arg):
            '''Parse an arg and return array of bool'''

            # 'X1' needs a strict match but 'X' matches all 'X1', 'X2',...
            if type(_arg) is str or type(_arg) is bytes or issubdtype(type(_arg), str_):
                # '' or b'': nothing selected
                if not _arg:
                    return zeros(shape=self._master.shape, dtype=bool)
                if not strict:
                    try:
                        return char.startswith(char.lower(self.names), _arg.lower().encode())
                    except:
                        return char.startswith(char.lower(self.names), _arg.lower())
                else:
                    # YL 06/11/2023: bug fixed: char.equal(), which element-wise compares two arrays of the same shape, was a wrong function to use!
                    #                we can use in1d(a, x) here which works only when both a and x are str (dtype='U') or both are bytes (dtype='S')
                    try:
                        return isin(char.lower(self.names), _arg.lower().strip().encode())
                    except:
                        return isin(char.lower(self.names), _arg.lower().strip())
            elif miscutils.isInteger(_arg):
                if names_only:
                    return zeros(shape=self._master.shape, dtype=bool)
                else:
                    return self.znums == _arg
            else:
                return zeros(shape=self._master.shape, dtype=bool)

        for arg in args:
            if type(arg) is list or type(arg) is tuple or type(arg) is ndarray:
                for a in arg:
                    selected += _parseArg(a)
            else:
                selected += _parseArg(arg)

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


    def selectByRegions(self, *args, return_masks=False):
        '''Return atom indices of a labelled region according to `names` (e.g., 4 for all 'O4' and 'Mg4' atoms)'''

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


        int2region = lambda x: where(defchararray.count(defchararray.decode(self.names), str(x)))[0]

        try:
            # allows frag.selectByRegions(4), frag.selectByRegions(1,2,3,...), and frag.selectByRegions([1,2,3])
            selected = concatenate(tuple([ int2region(arg) for arg in utils.iterutils.getFlattenedIter(args) ]))
            if return_masks:
                return isin(range(self.natoms), selected)
            else:
                return selected
        # except 
        except:
            if return_masks:
                return zeros(shape=self._master.shape, dtype=bool)
            else:
                return array([], dtype=int)


    def selectByIndices(self, *args, return_masks=False):
        ''''''

        from numpy   import arange, array, isin, zeros, integer

        selected = zeros(shape=(self.natoms,), dtype=bool)
        indices_all = arange(self.natoms)

        for arg in args:
            flattened = array(utils.iterutils.getFlattened([arg]))
            if issubclass(flattened.dtype.type, integer):
                selected += isin(indices_all, flattened)

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


    def selectByRadius(self, radius, range=[], centre=None, unit='a.u.', boundary=None, return_masks=False, **kwargs):
        '''Return an indices list containing all atoms within given radius from centre (default centre is centroid)'''

        from numpy  import arange, array, asarray, isin, linalg, subtract, where, unique, zeros
        from ..data import unitconvert

        # default: all atoms
        if not len(range):
            range = arange(self.natoms)
        else:
            range = asarray(utils.iterutils.getFlattened([range]))

        if radius == 0.0:
            if return_masks:
                return zeros(shape=(self.natoms,), dtype=bool)
            else:
                return array([], dtype=int)

        # default centre is the centroid
        if isinstance(centre, type(None)):
            centre = self.centroid

        # centre could be the position of a given atom index (integer)
        if isinstance(centre, int):
            centre = self.coords[centre]

        # table of distances
        dists = linalg.norm(subtract(self.coords[range], centre), axis=1)

        # radius (distance from centre)
        radius *= unitconvert.unit2au.get(unit.lower(), 'a.u.')

        masks_in_range = dists < radius

        # map the range back to the original sequence
        masks = isin(arange(self.natoms), range[masks_in_range])

        selected = where(masks)[0]

        # boundary molecules/fragments included in selection
        if str(boundary)[:4].lower() == 'incl':
            masks = isin(self.fragIDs, unique(self.fragIDs[masks]))
            selected = arange(self.natoms)[masks]
        # boundary molecules/fragments excluded in selection
        elif str(boundary)[:4].lower() == 'excl':
            masks = isin(self.fragIDs, unique(self.fragIDs[~masks]), invert=True)
            selected = arange(self.natoms)[masks]

        if return_masks:
            return masks
        else:
            return selected


    def selectByXYZ(self,
                    min_x=0.0, min_y=0.0, min_z=0.0,
                    max_x=0.0, max_y=0.0, max_z=0.0,
                    origin=[0.0,0.0,0.0], # TO IMPLEMENT
                    unit='a.u.',
                    boundary=None,
                    return_masks=False,
                    **kwargs):
        ''''''

        from numpy import all, arange, isin, linalg, subtract, unique

        masks = all([max_x,max_y,max_z]>self.coords, axis=1)*all(self.coords>[min_x,min_y,min_z], axis=1)

        # boundary molecules/fragments included in selection
        if str(boundary)[:4].lower() == 'incl':
            masks = isin(self.fragIDs, unique(self.fragIDs[masks]))
        # boundary molecules/fragments excluded in selection
        elif str(boundary)[:4].lower() == 'excl':
            masks = isin(self.fragIDs, unique(self.fragIDs[~masks]), invert=True)

        if return_masks:
            return masks
        else:
            selected = arange(self.natoms)[masks]
            return selected


    def selectUnbonded(self):
        '''Return a indices list of all unbonded atoms'''

        from numpy import all, where

        return where(all(self.conn==-1, axis=1))[0]


    def selectByConvex(self, range=[], padding=0.0, unit='a.u.', invert=False, return_inside=False, return_masks=False):
        '''Select the convex hull atoms'''

        from numpy   import append, arange, asarray, isin, ones, zeros
        from time    import time
        from ..base  import geometry
        from ..data  import unitconvert

        # default: all atoms
        if not len(range):
            range = arange(self.natoms)
        else:
            range = asarray(utils.iterutils.getFlattened([range]))

        if return_inside:
            # not including the selected range of atoms
            masks_range_inverted = isin(arange(self.natoms), range, invert=True)

        # YL 13/01/2022: also get the oriented normals (see base.geometry.getConvexHull())
        coords, vertices, edges, faces, normals, ds = geometry.getConvexHull(self.coords[range], return_vertices=True, return_edges_and_faces=True)

        if return_inside:
            # atoms inside the convex set have positive distances to all the oriented faces
            if invert:
                # return all atoms at a distance from any face further than padding along the negative (exterior) direction
                masks_dists_positive = zeros(shape=(masks_range_inverted.sum(),), dtype=bool)
                for (normal, d) in zip(normals, ds):
                    masks_outside = self._getDistsPoints2Plane(normal, d, range=masks_range_inverted) < -padding*unitconvert.unit2au.get(unit, 1.0)
                    masks_dists_positive += masks_outside
            else:
                # return atoms within a padding distance from all faces along the positive (innerior) direction
                masks_dists_positive = ones(shape=(masks_range_inverted.sum(),), dtype=bool)
                for (normal, d) in zip(normals, ds):
                    masks_dists_positive *= self._getDistsPoints2Plane(normal, d, range=masks_range_inverted) >= padding*unitconvert.unit2au.get(unit, 1.0)

            if return_masks:
                masks_inside = isin(arange(self.natoms), arange(self.natoms)[masks_range_inverted][masks_dists_positive])
                return isin(arange(self.natoms), range[vertices]), masks_inside
            else:
                indices_inside = arange(self.natoms)[masks_range_inverted][masks_dists_positive]
                return range[vertices], indices_inside

        else:
            if return_masks:
                return isin(arange(self.natoms), range[vertices])
            else:
                return range[vertices]


    # YL 18/01/2022: revised to have a fast and convenient convex-based method
    @miscutils.time
    def selectByShell(self, around=[], convex=True, padding=5.0, unit='a.u.', boundary=None, return_masks=False):
        '''Select atoms as a shell around chosen atoms with a given distance'''

        from numpy   import arange, array, asarray, isin, ones, zeros
        from time    import time
        from ..base  import connectivity, geometry
        from ..io    import format

        # return all atoms if no references are given 
        if not len(around):
            format.printwrap('NB: All atoms are selected because no valid reference atoms are defined via the keyword option around=',
                              initial_indent=' ', subsequent_indent=' '*5, indent='')
            if return_masks:
                return ones((self.natoms,), dtype=bool)
            else:
                return arange(self.natoms)
        else:
            around = asarray(utils.iterutils.getFlattened([around]))

        # YL 18/01/2022: corrected the convex algorithm and kept it default as it runs way faster than the radius-based selection
        if convex:
            t0 = time()
            convexes, outside = self.selectByConvex(range=around, padding=padding, unit=unit, invert=True, return_inside=True, return_masks=True)
            masks_selected = ~outside
        else:
            looper = around
            # exclude the atoms around those we're selecting
            masks_selected = isin(arange(self.natoms), around)
            remainings = arange(self.natoms)[~masks_selected]

            # loop over the atoms around those we're selecting
            for i in looper:
                # exclude the atoms already selected to speed up
                masks_selected += self.selectByRadius(padding, range=remainings, unit=unit, centre=self.coords[i], boundary=boundary, return_masks=True)
                remainings = arange(self.natoms)[~masks_selected]

        if return_masks:
            return masks_selected
        else:
            return arange(self.natoms)[masks_selected]


# END OF SELECTION METHODS #
# FACILITIES #

    def castTo(self, inst):
        '''Cast to another class and return (not in place)'''

        from .fragment import Fragment
        from .fragbio  import FragmentBio
        from .xyz      import _XYZ

        selectcases = { 'fragment':Fragment,
                        'bio'     :FragmentBio,
                        'fragbio' :FragmentBio,
                      }

        # YL 10/02/2021: subtractive embedding branch merged with the biological branch so that the merge method must be _XYZ.merge()
        # instance
        if isinstance(inst, Fragment) or isinstance(inst, FragmentBio):
            inst = _XYZ.merge(self, cls=inst.__class__, noself=False)
        # str as class name
        elif isinstance(inst, str):
            inst = _XYZ.merge(self, cls=selectcases.get(inst.lower(), None), noself=False)
        # class
        else:
            inst = _XYZ.merge(self, cls=inst, noself=False)

        # also take the stacks
        inst.stacks = {}
        inst.stacks.update(self.stacks)

        return inst


    # YL 13/09/2020: overload _XYZ.getSelected() with `connect_mode`
    def getSelected(self, selected, flatten=True):
        '''Return a new instance with selected centres'''

        from numpy   import integer, issubdtype

        # a single integer
        if isinstance(selected, (int, integer)):
            selected = [selected]

        # YL 20/08/2021: added a switch if the developer knows the input sequence is
        #                already flattened because when the list is large utils.iterutils.getFlattened() may be very slow
        if flatten:
            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
        inst._master.setField(self._master[selected])

        # YL 25/11/2021: however connectivity must be reset!
        inst._master.conn = -1

        inst.connect_mode = self.connect_mode

#        # YL 26/08/2021: unlike merge(), fragIDs won't need to be recalculated
#        inst._recalc_fragIDs = True

        # inherit the title
        inst.title = self.title

        return inst


    @miscutils.memusage
    def getLatticeVectors(self):
        '''NDLP algorithm:
           Henk Bekker, Jur P. van den Berg, Tsjerk A. Wassenaar
           P.M.A. Sloot et al. (Eds.): ICCS 2003, LNCS 2659, pp. 70−79, 2003
        '''

        from itertools import combinations, product
        from numpy     import add, all, append, argsort, argwhere, array, equal, float64, isin, linalg, meshgrid, sign, sort, sum, unique
        from math      import ceil
        from .         import grid
        from ..base    import geometry

        SMALL = 1.0E-10

        # recentre
        self.recentre()

        # make a copy to avoid recalculating
        origin = array(self.centroid, dtype=float64)

        # convex hull
        hull = geometry.getConvexHull(self.coords)

        # get N by Minkowski sum: N ≡ M ⊕ -M
        M = origin - hull
        # we can afford memory for the whole cartesian product
        # N is symmetric
        N = unique(add(M[:,None,:], -M).reshape((M.shape[0]**2,3)), axis=0)

        with open('N.xyz', 'w') as fp:
            fp.write(str(N.shape[0]))
            fp.write('\nN\n')
            for i in N:
                strbuff = 'H {} {} {}\n'.format(*i)
                fp.write(strbuff)


        def _getIntersections(vertices0, ifaces0, vertices1, ifaces1, save=False):
            ''''''

            nonlocal SMALL


            intersections = array([], dtype=float64).reshape((0,3))
            intersect_faces0 = array([], dtype=int).reshape((0,3))
            intersect_faces1 = array([], dtype=int).reshape((0,3))

            # loop over faces on vertices0
            for iface0 in ifaces0:
                plane0 = geometry.getPlane(*vertices0[iface0])
                iA0, iB0, iC0 = iface0[0], iface0[1], iface0[2]

                # loop over faces on vertices1
                for iface1 in ifaces1:
                    plane1 = geometry.getPlane(*vertices1[iface1])
                    # distances from vertices A_a, B_a, and C_a of face (triangle) on vertices1 to plane (face) on vertices0
                    iA1, iB1, iC1 = iface1[0], iface1[1], iface1[2]

                    dist_A1 = geometry.getDistPoint2Plane([*vertices1[iA1]], *plane0)
                    dist_B1 = geometry.getDistPoint2Plane([*vertices1[iB1]], *plane0)
                    dist_C1 = geometry.getDistPoint2Plane([*vertices1[iC1]], *plane0)
                    # they can't be of same sign to intersect
                    if not ((dist_A1>-SMALL and dist_B1>-SMALL and dist_C1>-SMALL) or
                            (dist_A1< SMALL and dist_B1< SMALL and dist_C1< SMALL)):
                        # distances from vertices A, B, and C of face (triangle) on vertices0 to plane (face) on vertices1
                        dist_A0 = geometry.getDistPoint2Plane([*vertices0[iA0]], *plane1)
                        dist_B0 = geometry.getDistPoint2Plane([*vertices0[iB0]], *plane1)
                        dist_C0 = geometry.getDistPoint2Plane([*vertices0[iC0]], *plane1)
                        # they can't be of same sign to intersect
                        if not ((dist_A0>-SMALL and dist_B0>-SMALL and dist_C0>-SMALL) or
                                (dist_A0< SMALL and dist_B0< SMALL and dist_C0< SMALL)):
                            for iI, iJ in combinations([iA0,iB0,iC0], 2):
                                vec_IJ = (vertices0[iI]-vertices0[iJ])/linalg.norm(vertices0[iI]-vertices0[iJ])
                                intersect = geometry.getRayTriangleIntersect(vertices0[iI], vec_IJ, vertices1[iA1], vertices1[iB1], vertices1[iC1])
                                if intersect is not None:
                                    intersections = append(intersections, [intersect], axis=0)
                                    intersect_faces1 = append(intersect_faces1, [[iA1,iB1,iC1]], axis=0)
                                    if save:
                                        print("### intersect_faces1 =", intersect_faces1.shape, vertices0[iI], vec_IJ, vertices1[iA1], vertices1[iB1], vertices1[iC1])
                                        with open('{}-{}-{}.xyz'.format(iA1,iB1,iC1), 'w') as fp:
                                            fp.write(str(semisphere.shape[0]+3))
                                            fp.write('\n{}-{}-{}\n'.format(iA0,iB0,iC0))
                                            for i in semisphere:
                                                strbuff = 'H {} {} {}\n'.format(*i)
                                                fp.write(strbuff)
                                            strbuff = 'P {} {} {}\n'.format(*N_vertices[iA1])
                                            fp.write(strbuff)
                                            strbuff = 'P {} {} {}\n'.format(*N_vertices[iB1])
                                            fp.write(strbuff)
                                            strbuff = 'P {} {} {}\n'.format(*N_vertices[iC1])
                                            fp.write(strbuff)
                            for iI, iJ in combinations([iA1,iB1,iC1], 2):
                                vec_IJ = (vertices1[iI]-vertices1[iJ])/linalg.norm(vertices1[iI]-vertices1[iJ])
                                intersect = geometry.getRayTriangleIntersect(vertices1[iI], vec_IJ, vertices0[iA0], vertices0[iB0], vertices0[iC0])
                                if intersect is not None:
                                    intersections = append(intersections, [intersect], axis=0)
                                    intersect_faces0 = append(intersect_faces0, [[iA0,iB0,iC0]], axis=0)
                                    if save:
                                        print("### intersect_faces0 =", intersect_faces0.shape, vertices1[iI], vec_IJ, vertices0[iA0], vertices0[iB0], vertices0[iC0])
                                        with open('{}-{}-{}.xyz'.format(iA0,iB0,iC0), 'w') as fp:
                                            fp.write(str(semisphere.shape[0]+3))
                                            fp.write('\n{}-{}-{}\n'.format(iA0,iB0,iC0))
                                            for i in semisphere:
                                                strbuff = 'H {} {} {}\n'.format(*i)
                                                fp.write(strbuff)
                                            strbuff = 'P {} {} {}\n'.format(*N_vertices[iA0])
                                            fp.write(strbuff)
                                            strbuff = 'P {} {} {}\n'.format(*N_vertices[iB0])
                                            fp.write(strbuff)
                                            strbuff = 'P {} {} {}\n'.format(*N_vertices[iC0])
                                            fp.write(strbuff)

            return intersections, intersect_faces0, intersect_faces1


        # the convex hull of N is also symmetric
        N_vertices, iN_edges, iN_faces = geometry.getConvexHull(N, return_edges_and_faces=True)
        print("### N_vertices =", N_vertices.shape)
        ind_semi = argsort(N_vertices[:,0])[:ceil(N_vertices.shape[0]/2)]
        semisphere = N_vertices[ind_semi]
        # edges only on the semisphere
        iN_edges_semi = iN_edges[all(isin(iN_edges, ind_semi).reshape((iN_edges.shape[0],2)), axis=1)]
        iN_faces_semi = iN_faces[all(isin(iN_faces, ind_semi).reshape((iN_faces.shape[0],3)), axis=1)]


        for a in semisphere:

            N_a = a + N_vertices
#            ind = argwhere(all(N_a==semisphere, axis=1))
#            print("### N_a =", N_a[ind])
#            print("### N_a =", add(N_a[:,None,:],-semisphere))

            intersect_N_Na, intersect_N_Na_faces0, intersect_N_Na_faces1 = _getIntersections(N_vertices, iN_faces_semi, N_a, iN_faces)
            print("### intersect_N_Na =", a, intersect_N_Na.shape[0], intersect_N_Na_faces0.shape[0], intersect_N_Na_faces1.shape[0])


            for b in intersect_N_Na:
                N_b = b + N_vertices
                with open('N_a.xyz', 'w') as fp:
                    fp.write(str(semisphere.shape[0]+N_a.shape[0]+N_b.shape[0]+intersect_N_Na.shape[0]))
                    fp.write('\nN_a\n')
                    for i in semisphere:
                        strbuff = 'H {} {} {}\n'.format(*i)
                        fp.write(strbuff)
                    for i in N_a:
                        strbuff = 'O {} {} {}\n'.format(*i)
                        fp.write(strbuff)
                    for i in N_b:
                        strbuff = 'N {} {} {}\n'.format(*i)
                        fp.write(strbuff)
                    for i in intersect_N_Na:
                        strbuff = 'F {} {} {}\n'.format(*i)
                        fp.write(strbuff)
                    for f in intersect_N_Na_faces0:
                        for j in f:
                            strbuff = 'P {} {} {}\n'.format(*N_vertices[j])
                            fp.write(strbuff)
                intersect_N_Nb, intersect_N_Nb_faces0, intersect_N_Nb_faces1 = _getIntersections(N_vertices, intersect_N_Na_faces0, N_b, iN_faces, save=True)
#                intersect_Na_Nb = _getIntersections(N_a, iN_faces, N_b, iN_faces)
                print("### intersect_N_Nb =", b, intersect_N_Nb.shape)
                quit()

        print("### nvertices =", N_vertices.shape)


    @miscutils.time
    def solvate(self,
                boundary=4.53534,               # 4.53534 Bohr =  2.4 Angstrom
                padding=28.3459,                # 28.3459 Bohr = 15.0 Angstrom
                save_to='_solvated.pdb',
                save_to_solvent='_solvent.pdb',
                shape='cube',
                solvent='TIP3P',
                recentre=True,
                reconnect=True,
                unit='a.u.'):
        '''Solvate the fragment'''

        from numpy  import any, isin, linalg, load, ones, subtract, unique, zeros
        from os     import path
        from ..data import unitconvert
        from ..base import connectivity, constants

        print("\n Solvating...")
        stdout.flush()

        unit_factor = unitconvert.unit2au.get(unit.lower(), 1.0)
        boundary *= unit_factor
        padding  *= unit_factor

        print(f' shape   : {shape.lower()}')
        print(f' model   : {solvent.upper()}')
        print(f' padding : {padding} a.u.')
        print(f' boundary: {boundary} a.u.')
        stdout.flush()

        # solvent model and associated data file
        selectcases = { 'tip3p': 'tip3p_cube_150A.npz' }

        # location of built-in solvent data
        datadir = path.abspath(path.join(path.dirname(__file__), '..', 'data', 'solvents'))

        # load data and save to master array
        abuff = load(path.join(datadir, selectcases[solvent.lower()]), fix_imports=False)['arr_0']
        solvFrag = self.__class__(connect_mode=None)
        solvFrag.resize(abuff.shape[0])

        # similar to the merge() method in <class _XYZ>
        for field in solvFrag._master.dtype.fields:
            if field not in abuff.dtype.fields:
                continue
            try:
                solvFrag._master.setField(abuff[field], field=field)
            except:
                solvFrag._master.setField(abuff[field].reshape(-1,1), field=field)

        # recentre
        if recentre:
            self.recentre()
            solvFrag.recentre()
        # translate solvent fragment
        else:
            solvFrag.coords -= self.centroid

        # select solvent atoms: solute's radius plus padding distance 
        cutoff = self.radius + padding

        # shape sphere
        if shape.lower() == 'sphere' or shape.lower() == 'droplet':
            selected = solvFrag.select(radius=cutoff, boundary='excl')
        # shape box
        elif shape.lower() == 'box' or shape.lower() == 'cuboid':
            selected = solvFrag.select(min_x=self.x.min()-padding,
                                       min_y=self.y.min()-padding,
                                       min_z=self.z.min()-padding,
                                       max_x=self.x.max()+padding,
                                       max_y=self.y.max()+padding,
                                       max_z=self.z.max()+padding,
                                       boundary='excl')
        elif shape.lower() == 'cube':
            max_strech = self._max_abs_coords.max()
            selected = solvFrag.select(min_x=-max_strech-padding,
                                       min_y=-max_strech-padding,
                                       min_z=-max_strech-padding,
                                       max_x= max_strech+padding,
                                       max_y= max_strech+padding,
                                       max_z= max_strech+padding,
                                       boundary='excl')
        # unknown
        else:
            print("\n >>> ERROR: unknown solvent shape: {}\n".format(self.shape))
            exit(112)

        # cut out from the original solvent
        # YL 27/11/2021: this will be extremely slow if flatten=True which is unnecessary
        solvFrag = solvFrag.getSelected(selected, flatten=False)
        solvent_coords = solvFrag.coords

        # solvent atoms to delete
        todels = zeros((solvFrag.natoms), dtype=bool)

        # YL 26/08/2021: improved method 4x faster than calculating the full distances
        #                of each centres pair previously - we create a merged fragment
        #                object and calculate the connectivity using `boundary` as a
        #                unitary distance
        tmpFrag = solvFrag.merge(self, reconnect=False)
        # YL 26/08/2021: see changes in connectivity.connect()
        connectivity.connect(tmpFrag, unitary_dist=boundary)
        if chemsh.CHEMSH_DEBUG > 2:
             print(f'\n >>> solvate: memory used for temporary connectivity: {tmpFrag._conn_aux.nbytes/1024/1024} MB')
        range_solute = range(solvFrag.natoms, solvFrag.natoms+self.natoms)
        for i in range(tmpFrag._nconns_aux):
            todels += isin(tmpFrag._conn_aux[:solvFrag.natoms,i], range_solute)
#        # use explicit loop to save memory
#        # this is rather slow but we can't use the substrate's convex points to solvate
#        for icoords in self.coords:
#            dist = linalg.norm(solvent_coords-icoords, axis=1)
#            todels += dist <= boundary

        # exclude solvent molecules broken by selection
        todels = isin(solvFrag.fragIDs, unique(solvFrag.fragIDs[todels]))

        # cut a new solvent fragment out of the original selection
        # YL 27/11/2021: this will be extremely slow if flatten=True which is unnecessary
        solvFrag = solvFrag.getSelected(arange(solvFrag.natoms)[~todels], flatten=False)

        # save the solvent structure
        solvFrag.save(save_to_solvent)

        # merge solvent fragment with solute
        solvatedFrag = self.merge(solvFrag, reconnect=reconnect)

        # title
        solvatedFrag.title = self.title.strip() + '_solvated'

        # YL 26/11/2021: inherit renamed residues
        solvatedFrag._renamed_residues.update(self._renamed_residues)

        # save the solvated structure
        solvatedFrag.save(save_to)

        return solvatedFrag


# END OF FACILITIES #

#    def remove(self, foo):
#        '''remove atom(s)'''
#
#        from numpy import delete
#
#        if type(foo) is int:
#            bar = [foo]
#
#        elif type(foo) is list:
#            bar = foo
#
#        elif type(foo) is str:
#            if foo.lower() == "all":
#                bar = list(range(self.nat))
#
#        for i in bar:
#            self.coords = delete(self.coords, i, 0)
#            self.znums  = delete(self.znums , i, 0)
#            if self.charges > 0:
#                self.charges = delete(self.charges, i, 0)


# YL TODO improve
# currently FHI-aims needs selectByRadius
BQs.selectByRegions = _FragmentBase.selectByRegions
BQs.selectByRadius  = _FragmentBase.selectByRadius
