#  Copyright (C) 2019 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/>.

__author__ = 'You Lu <you.lu@stfc.ac.uk>'

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

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

    # default size (natoms) is 0, but the shapes (of charges, conn, coords, etc.) need to be (1,...) to hold the default values
    _kwargs = {
                'bfactors'   : zeros(shape=(1,) ,     dtype=float64),
                'bqs'        : BQs(),
                'chainIDs'   : full((1),         '' , dtype='S8'),
                'frozen'     : zeros(shape=(1,) ,     dtype=int64),
                'grid'       : None,
                'natoms'     : 0,
                'nbqs'       : 0,
                'occ'        : ones(shape=(1,)  ,     dtype=float64),
                'pH'         :-1.0,
                'radii'      : zeros(shape=(1,) ,     dtype=float64),
                'residues'   : full((1)         , '', dtype='S8'),
                'resIDs'     : zeros(shape=(1,) ,     dtype=int64),
                'segments'   : full((1)         , '', dtype='S8'),
                'segIDs'     : zeros(shape=(1,) ,     dtype=int64),
                'shellatoms' : full((1,)        , -1, dtype=int64),
                'shells'     : Shells(),
                'totalcharge': 0.0,
                'types'      : full((1)         , '', dtype='S8'),
                'zmatrix'    : zeros(shape=(1,),      dtype=dtype([('j',int64),('bond',float64),('k',int64),('angle',float64),('l',int64),('dihedral',float64)])),
                'znums'      : zeros(shape=(1,),      dtype=int64),
              }

    # fields of <nputils.RecArray _master>: the order matters!
    # YL 31/10/2019: _FragmentBase._fields = ('names', 'znums', 'coords', 'charges', 'conn', 'fragIDs', 'frozen', 'masses', 'shellatoms')
    _fields = tuple(list(_FragmentBase._fields)+['types', 'residues', 'resIDs', 'segments', 'chainIDs', 'occ', 'bfactors', 'radii', 'zmatrix'])

    # make sure 'unit' is processed before 'coords' and
    #           'shellatoms' before 'coords' (connectivity depends on shell charges!)
    # YL 19/02/2021: we need to set `charges` prior to coords because when coords are defined by
    #                PDB parsing charges are the only missing information and due to residues connected
    #                to the main segment the atoms sequence (FragmentBio._master) may be reordered
    # YL 03/03/2021: added names prior to `coords` for the same reason
    _priorities = 'connect_mode', 'charges', 'names', 'coords', 'shellatoms'

    _internals = {
                   '_calcZMatrix'     : True,
                   '_renamed_residues': {},
                 }

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

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


    @property
    def zmatrix(self):
        '''Internal coordinates (Z-matrix)'''

        from numpy import isfinite, ma, stack

        if self._calcZMatrix and self.nfragments == 1:
            # cartesian to z-matrix
            self.c2z()

        b = self._master.zmatrix.bond
        a = self._master.zmatrix.angle
        d = self._master.zmatrix.dihedral
        j = self._master.zmatrix.j
        k = self._master.zmatrix.k
        l = self._master.zmatrix.l

        # array of masks
        mask = stack([j==-1, ~(isfinite(b)),
                      k==-1, ~(isfinite(a)),
                      l==-1, ~(isfinite(d))],
                     axis=-1)


        # construct the masked array
        marray = self._master.zmatrix.view(ma.MaskedArray)

        # apply the masking rules
        for i in range(marray.mask.size):
            marray.mask[i] = tuple(mask[i])

        return marray


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

        self._zmatrix = val


    @property
    def chainIDs(self):
        '''Chain IDs'''

        return self._master.chainIDs


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

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


    @property
    def nchains(self):
        '''Number of unique chains'''

        from numpy import unique

        return unique(self.chains).size


    # segment names (normally defined by PDB)
    @property
    def segments(self):
        '''Segment names'''

        return self._master.segments


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

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


    @property
    def all_segments(self):
        '''Segment names without repeat'''

        # YL 12/09/2021: we don't recalculate this every time from now on but once
        #                segments are reassigned self._all_segments MUST be recalculated!
        try:
            return self._all_segments
        except:
            from numpy import unique
            self._all_segments, self._all_segments_inverse_indices = unique(self.segments, return_inverse=True)
            return self._all_segments


    @property
    def all_segments_inverse_indices(self):
        '''Indices to reconstruct `segments` from `all_segments`'''

        # we don't recalculate this every time from now on but once
        # segments are reassigned self._all_segments MUST be recalculated!
        try:
            return self._all_segments_inverse_indices
        except:
            from numpy import unique
            self._all_segments, self._all_segments_inverse_indices = unique(self.segments, return_inverse=True)
            return self._all_segments_inverse_indices


    @property
    def nsegments(self):
        '''Number of unique segments'''

        from numpy import unique

        return unique(self.segments).size


    @property
    def residues(self):
        '''Residue names'''

        return self._master.residues


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

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


    @property
    def nresidues(self):
        '''Number of unique residues'''

        from numpy import unique

        return unique(self.residues).size


    @property
    def nresIDs(self):
        '''Number of unique residues'''

        from numpy import unique

        return unique(self.resIDs).size


    @property
    def ncoordinations(self):
        '''Coordination number'''

        from numpy import sum

        return sum(self.conn > -1, axis=1)


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

        return self._master.types


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

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


    @property
    def resIDs(self):
        '''Residue IDs'''

        return self._master.resIDs


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

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


    @property
    def occ(self):
        '''Occupancies'''

        return self._master.occ


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

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


    @property
    def bfactors(self):
        '''Temperature factors (b-factors)'''

        return self._master.bfactors


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

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


    @property
    def radii(self):
        '''Radii'''

        return self._master.radii


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

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


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

        return self._pH


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

        self._pH = val


    @property
    def resIDs2residues(self):
        '''Mapping between resIDs and residues'''

        # there may be same resIDs in different segments so that we can't use unique()
        return self._master[['resIDs','residues']]


    def resID2residue(self, resID, segment='', mute=False, decode=True):
        '''Return the matching residue name of a given resID'''

        from numpy import unique

        result = {}

        # if segment is not specified: loop over all segments
        if not segment:
            for segment in self.all_segments:
                selected = self.select(segments=segment)
                slices = dict(unique(self.resIDs2residues[selected]))
                val = slices.get(resID, None)
                if val:
                    result.update({segment:val})

        # if segment is specified
        else:
            selected = self.select(segments=segment)
            slices = dict(unique(self.resIDs2residues[selected]))
            val = slices.get(resID, None)
            if val:
                result.update({segment:val})

        # no matching
        if len(result) == 0:
            if not mute:
                print("\n WARNING: resID2residue() has not found a matching residue with the given resID {}\n Returning None...\n".format(resID))
            return None

        # unique matching
        if len(result) == 1:
            if decode:
                # YL TODO 15/03/2024: should apply .strip(), too, but I'm afraid of breaking things
                return list(result.values())[0].decode()
            else:
                return list(result.values())[0]

        # multiple matching
        if len(result) > 1:
            if not mute:
                print("\n WARNING: resID2residue() has found more than one matching residues with the given resID {}\n Returning a dictionary...\n".format(resID))
            if not mute and segment:
                print("\n WARNING: resID2residue() has found more than one matching residues with the given resID {}\n There seems to be something wrong with segment {}\n Returning a dictionary...\n".format(resID, segment))
            return result


    class _Terminal():
        '''Container of terminal'''

        @property
        def residue(self):

            return self._residue


        @residue.setter
        def residue(self, val):

            try:
                self._residue = str(val.decode())
            except:
                self._residue = str(val)


        @property
        def resID(self):

            return self._resID


        @resID.setter
        def resID(self, val):

            self._resID = int(val)


        @property
        def segment(self):

            return self._segment


        @segment.setter
        def segment(self, val):

            self._segment = val


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

        try:
            return self._terminal
        except:
            self._terminal = self._Terminal()
            return self._terminal


    @terminal.setter
    def terminal(self, val):
        ''''''

        self._terminal = val


    def deleteIons(self, ff='charmm', mute=False, reconnect=True):
        '''Remove all ionic atoms from the current fragment and return the deleted part'''

        return self.delete(self.selectByIons(ff=ff), mute=mute, reconnect=reconnect)


    def delete(self, *args, mute=False, reconnect=True, **kwargs):
        '''An advanced version of _XYZ.delete(): also returns the deleted part'''

        from numpy   import delete

        # YL 25/09/2020: optional positioned arguments as indices
        if args == ():
            selected = self.select(**kwargs, return_masks=False)
        else:
            indices = utils.iterutils.getFlattened(list(args))
            selected = self.select(indices=indices, **kwargs, return_masks=False)

        # if this entry is also available (from PDB)
        try:
            self._pdb_entries = delete(self._pdb_entries, selected)
        except:
            pass

        self._psf.delete(selected)

        deletedFrag = _XYZ.delete(self, selected, mute=mute)
        # YL 30/01/2022: all_segments_inverse_indices needs to be recalculated
        try:
            del self._all_segments_inverse_indices
        except:
            pass

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

        return deletedFrag


    def reorder(self, order):
        '''Wrapper function of _FragmentBase.reorder() with rearrangement of Z-matrix'''

        from numpy import arange, argsort, array, count_nonzero, searchsorted

        _FragmentBase.reorder(self, order)

        # YL 28/12/2021: rearrange Z-matrix
        if count_nonzero(self._master.zmatrix):
            new = arange(self.natoms)
            old = array(order)
            mapping = dict(zip(old,new))
            m = self._master.zmatrix
            # update the indices
            for item in [ 'j', 'k', 'l' ]:
                sorter = argsort(old)
                indices = searchsorted(old, m[item], sorter=sorter)
                m[item] = new[sorter][indices]


    def merge(self, *args, noself=False, reset_bfactors=False, reconnect=True):
        '''Biological version of _XYZ.merge()'''

        from numpy import array

        inst = _XYZ.merge(self, *args, noself=noself)

        if noself:
            frags = list(args)
        else:
            frags = [self] + list(args)

        # merge Z-matrix
        for i, (k, v) in enumerate(inst.stacks.items()):
            if 'zmatrix' not in frags[i]._master.dtype.fields:
                if chemsh.CHEMSH_DEBUG > 2:
                    print("\n >>> fragmentBio.merge(): fragment object {} does not have Z-matrix. Skipping...".format(frags[i]))
                continue
            try:
                offset = list(v)[0]
            # list(v) could be []
            except:
                offset = 0
            inst._master['zmatrix'][v] = frags[i]._master['zmatrix'][:]
            for item in [ 'j', 'k', 'l' ]:
                valid = array(v)[inst._master['zmatrix'][item][v] != -1]
                if valid.size:
                    inst._master['zmatrix'][item][valid] += offset

        inst.bfactors = 0.0
        if reconnect:
            inst.connect()

        # inherit all the titles
        if noself:
            inst.title = '+'.join([ '_'.join(arg.title.strip().split()) for arg in args if arg.title.strip() != 'title' ])
        else:
            inst.title = self.title.strip() + '+'.join([ '_'.join(arg.title.strip().split()) for arg in args if arg.title.strip() != 'title' ])

        # inherit the terminal
        inst.terminal = self.terminal

        # YL 26/08/2021: fragIDs must be recalculated
        # YL 15/03/2024: I never found this wrong but now inst.fragIDs must be firstly reset to -1
        inst.fragIDs = -1
        inst._recalc_fragIDs = True

        return inst


    # YL 02/12/2023: rewritten based on the new method of backbone selection using species check rather than FF names
    def selectByTruncation(self, truncate_by='alpha-beta', side_chain=True, return_masks=False):
        '''Return indices of selected atoms according to the truncating rules given'''

        from numpy import isin, ones

        # return all if truncation method or forcefield isn't supported
        masks = ones((self.natoms,), dtype=bool)

        # truncate by, e.g., alpha-beta bond
        #         HA (N)
        #    \    |   |
        # CB--\---CA--C==O
        #      \  |
        #         N--HN

        # YL 02/12/2023: this implicitly calls self.selectByBackbone(return_masks=True)
        indices_backbone = self.selectByBackbone(protein_only=True)

        resIDs_A = self.resIDs[self.bonds[:,0]]
        resIDs_B = self.resIDs[self.bonds[:,1]]
        seg_A = self.segments[self.bonds[:,0]]
        seg_B = self.segments[self.bonds[:,1]]
        # shape: (self.bonds.shape[0],)
        masks_intra = (resIDs_A == resIDs_B)*(seg_A == seg_B)
        # all intra-residue H atoms connected to the backbone
        masks_backbone_in_bonds = isin(self.bonds[:,0], indices_backbone) + isin(self.bonds[:,1], indices_backbone)
        indices_all_bonded_to_backbone = unique(self.bonds[masks_backbone_in_bonds])
        indices_H_bonded_to_backbone = indices_all_bonded_to_backbone[self.znums[indices_all_bonded_to_backbone]==1]

        # all residues involving the backbone atoms
        masks = self.selectBySameResidues(indices_backbone, return_masks=True)

        if truncate_by.strip().lower() in [ 'alpha-beta', 'alpha_beta', 'a-b', 'a_b' ]:
#            if ff.strip().lower() in [ 'charmm' ]:
#                masks = self.select(atoms=['O','N','C','CA','HN','HA'], strict=True, invert=side_chain, return_masks=True)
            if side_chain:
                # exclude atoms beyond C_alpha (except for H, for example in glycine)
                masks *= isin(self._indices_all, indices_backbone.tolist()+indices_H_bonded_to_backbone.tolist(), invert=True)
            else:
                masks *= isin(self._indices_all, indices_backbone.tolist()+indices_H_bonded_to_backbone.tolist())

        if return_masks:
            return masks
        else:
            # numerical indices
            selected = self._indices_all[masks]
            return selected


    # YL 02/12/2023
    def selectBySameResidues(self, selected, return_masks=False, invert=False):
         '''Return indices/masks of all atoms in the residues involved in the given indices of selection'''

         from numpy import ndarray, isin, unique, zeros

         # in case given a masks array
         if type(selected) is ndarray:
             if selected.dtype == bool:
                 selected = self._indices_all[selected]

         segs_selected = unique(self.segments[selected])
         masks_selected = isin(self._indices_all, selected)

         masks = zeros((self.natoms,), dtype=bool)
         # go through the selected segments
         for seg in segs_selected:
              masks_seg = self.segments == seg
              resIDs_selected_seg = unique(self.resIDs[masks_selected*masks_seg])
              masks += self.select(segments=seg, resIDs=resIDs_selected_seg, return_masks=True)

         if invert:
             masks = ~masks

         if return_masks:
             return masks
         else:
             return self._indices_all[masks]


    # YL 02/12/2023: rewritten based on the new method of backbone selection using species check rather than FF names
    def selectByAminoAcids(self,
                           truncate_by='alpha-beta',
                           side_chain=True,
                           return_masks=False,
                           **kwargs):
        '''Return indices of all amino acid residues: It is to be noticed that this is different
from selectByProtein() which selects the whole segments involving peptides'''

        # get the indices of backbone atoms
        indices_backbone = self.selectByBackbone(protein_only=True, return_masks=False)

        # all residues containing these backbone atoms are considered amino acids
        masks = self.selectBySameResidues(indices_backbone, return_masks=True)

        # apply truncation method
        # (in selectByTruncation() there is also a call to selectBySameResidues() but this is not redudant thinking of future tracation of nucleotides)
        masks *= self.selectByTruncation(truncate_by=truncate_by, side_chain=side_chain, return_masks=True)

        if return_masks:
            return masks
        # numerical indices
        else:
            return self._indices_all[masks]


    # YL 05/08/2021: changed the default behaviour of strict to True
    def selectBio(self,
                  chainIDs      = [],
                  resIDs        = [],
                  segIDs        = [],
                  residues      = [],
                  segments      = [],
#                  hybridisation = '',
                  truncate_by   ='',       # YL 15/12/2021: do not truncate by default
                  side_chain    = True,     # YL 15/12/2021: defaulted to choose the side chain (won't take effect unless truncate_by is defined)
                  ff            ='charmm',
                  return_masks  = False,
                  invert        = False,
                  union         = False,
                  strict        = True,
                  visualise     ='',
                  **kwargs):
        '''Return indices of selection'''

        # check keywords
        from inspect import getfullargspec, getcallargs
        printkeys = False
        for key in kwargs:
            # YL NB 29/06/2020: this method isn't always safe because the function being inspected cannot be wrapped by a decorator
#           and key not in getfullargspec(self.selectByHybridisation).args \
            if key not in getfullargspec(self.selectBio).args   \
           and key not in getfullargspec(self.selectByXYZ).args \
           and key not in getfullargspec(self.selectByRadius).args \
           and key != 'kwargs':
                print(' >>> WARNING: \"{}\" is not a valid keyword of function `select`. Nothing has been selected.\n'.format(key))
                printkeys = True
#        if printkeys:
#            print(' Valid keywords are:\n {}\n'.format(getfullargspec(self.selectBio).args))

        # get the keyword arguments (NB: including unwanted 'return_masks', 'kwargs', and so)
        kwargs = locals()

        # do the importing after referencing locals()
        from numpy    import arange, array, asarray, bytes_, char, integer, isin, ones, str_, where, zeros
        from ..io     import format

        # select nothing in the beginning if union
        if union:
            masks = zeros((self.natoms,), dtype=bool)
        # select all in the beginning if not union
        else:
            masks = ones((self.natoms,), dtype=bool)

        nothing_selected = True

        # loop over the keyword arguments 
        for key, val in kwargs.items():

            if not hasattr(self, key):
                continue

            # `[val]` can't be directly used because val may be a tuple
            if type(val) is tuple:
                val = list(val)

            # if there is no iterable argument value (see below masks += or *= iter_masks)
            iter_masks = None

            # it should be always "union" within an iterable argument value (e.g., residues=['A','B'])
            # YL 07/03/2021: this can make when [val]==[] iter_masks is None
            for j in utils.iterutils.getFlattened([val]):
                iter_masks = zeros((self.natoms,), dtype=bool)
                break

            # loop over the given values
            # has NOT to be `list(val)` because list() converts bytes to numbers
            for j in utils.iterutils.getFlattened([val]):

                nothing_selected = False
                # unspecified parameter will appears as numpy.ndarray([])
                try:
                    if not j.size:
                        continue
                except:
                    pass

                # YL 07/03/2021: the following code refactored
                dtt = getattr(self, key).dtype.type

                # the attribute array is of bytes type
                if issubclass(dtt, bytes_):
                    # always treat j an array
                    flattened = array(utils.iterutils.getFlattened([j]))
                    # non-strict: match the beginning
                    if not strict:
                        for k in flattened:
                            # k must be bytes, too
                            kbytes = utils.strutils.str2bytes(k)
                            iter_masks += char.startswith(char.lower(char.strip(getattr(self, key))), kbytes.strip().lower())
                    # strict matching
                    else:
                        for k in flattened:
                            # k must be bytes, too
                            kbytes = utils.strutils.str2bytes(k)
                            iter_masks += char.equal(char.lower(char.strip(getattr(self, key))), kbytes.strip().lower())

                # attribute array is str
                elif issubclass(dtt, str_):
                    # always treat j an array
                    flattened = array(utils.iterutils.getFlattened([j]))
                    # non-strict: match the beginning
                    if not strict:
                        for k in flattened:
                            # k must be str, too
                            kstr = utils.strutils.bytes2str(k)
                            iter_masks += char.startswith(char.lower(char.strip(getattr(self, key))), kstr.strip().lower())
                    # strict matching
                    else:
                        for k in flattened:
                            # k must be str, too
                            kstr = utils.strutils.bytes2str(k)
                            iter_masks += char.equal(char.lower(char.strip(getattr(self, key))), kstr.strip().lower())

                # the attribute array is of other types
                # we currently only support integer and exact matching of float numbers
                else:
                    # YL 07/03/2021: removed the warnings filter
                    flattened = array(utils.iterutils.getFlattened([j]))
                    # [j] is array of integers
                    if issubclass(flattened.dtype.type, integer):
                        iter_masks += isin(getattr(self, key), flattened)
                    # [j] is array of float
                    # YL TODO 07/03/2021: comparison with a tolerence
                    elif issubclass(flattened.dtype.type, float):
                        iter_masks += isin(getattr(self, key), flattened)
                    # YL 07/03/2021: do nothing because such comparison doesn't make sense or it
                    #                means: TODO more types should be covered
#                    else:
#                        iter_masks += getattr(self, key) == j

            # YL 15/12/2021: apply truncation method
            if truncate_by:
                masks *= self.selectByTruncation(truncate_by=truncate_by, side_chain=side_chain, return_masks=True)

            # if there is iter_masks
            try:
                if union:
                    masks += iter_masks
                else:
                    masks *= iter_masks
            # do nothing if there is no iter_masks
            except:
                pass

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

#        # select by hybridisation
#        hybrid_kwargs = { k:v for k, v in kwargs.items() if k not in [ 'self', 'hybridisation', 'return_masks', 'invert' ] and k in getfullargspec(self.selectByHybridisation).args }
#        if hybridisation:
#            if union:
#                masks += self.selectByHybridisation(hybridisation, return_masks=True, invert=invert, **hybrid_kwargs)
#            else:
#                masks *= self.selectByHybridisation(hybridisation, return_masks=True, invert=invert, **hybrid_kwargs)

        # YL 08/07/2020: return an all-False array if nothing actually selected
        if nothing_selected:
            masks = zeros((self.natoms,), dtype=bool)

        # invert selection
        if invert:
            masks = ~masks

        # return masks
        if return_masks:
            return masks
        # return array of indices
        else:
            return self._indices_all[masks]


    def selectByCommonNames(self, names, ff='charmm', side_chain=True, truncate_by='alpha-beta', return_masks=False, invert=False):
        '''Select via common names such as water, amino acid, etc'''

        from numpy   import arange, str_, zeros
        from ..data  import forcefield

        # flatten
        names = utils.iterutils.getFlattened([names])
   
        # get the forcefield object
        ffObj = forcefield.getFFObj(ff)

        masks_selected = zeros((self.natoms,), dtype=bool)
        for name in names:
            if type(name) not in [ str, str_ ]:
                print(f'\n WARNING: selectByCommonNames() cannot process option value {name}')
                continue
            # solvent: currently water
            if name.strip().lower().replace('_','').replace('-','').startswith('solv') or \
               name.strip().lower().replace('_','').replace('-','').startswith('water'):
                masks_selected += self.select(residues=ffObj.waters, return_masks=True)
            # non-solvent: currently water+ions
            elif name.strip().lower().replace('_','').replace('-','').startswith('nonsolv'):
                residues_not_to_wrap = ffObj.waters+ffObj.ions
                masks_selected += self.select(residues=residues_not_to_wrap, invert=True, return_masks=True)
            # amino acids
            elif name.strip().lower().replace('_','').replace('-','').startswith('amino'):
                masks_selected += self.selectByAminoAcids(side_chain=side_chain, truncate_by=truncate_by, return_masks=True)
            # protein
            elif name.strip().lower().replace('_','').replace('-','').startswith('prot'):
                masks_selected += self.selectByProtein(return_masks=True)
            # ions
            elif name.strip().lower().replace('_','').replace('-','').startswith('ion'):
                masks_selected += self.selectByIons(ff=ff, return_masks=True)
            # try explicit residue/segment names
            else:
                masks_selected += self.select(residues=name, segments=name, union=True, return_masks=True)

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


    def selectByShell(self, around=[], convex=True, padding=5.0, unit='a.u.', boundary='inclusive', return_masks=False):
        '''A biological version of function FragmentBase.selectByShell with residue-based boundary treatment'''

        from numpy import arange

        masks = _FragmentBase.selectByShell(self, around=around, convex=convex, padding=padding, boundary=boundary, unit=unit, return_masks=True)

        masks_selected = self.selectByBoundary(masks, boundary=boundary, unit=unit, return_masks=True, mute=True)

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


    def selectByBoundary(self, masks, boundary='inclusive', unit='a.u.', return_masks=False, mute=True):
        '''Deal with boundary residues for selected atoms'''

        from numpy   import all, arange, argsort, argwhere, copy, isin, unique
        from ..data  import unitconvert

        masks0 = copy(masks)
        segments_unique_selected = unique(self.segments[masks])
        segments_unique, segment_inverse = self.all_segments, self.all_segments_inverse_indices

        # boundary residues included in selection
        # NB: there could be repeated resIDs in a biological fragment
        if str(boundary)[:4].lower() == 'incl':
            # YL 12/09/2021: improved algorithm 400 faster than before
            for i, segment in enumerate(segments_unique):
                if segment not in segments_unique_selected:
                    continue
                # select the current segment: segment_inverse==i is faster than self.segments==segment
                masks_segment = segment_inverse==i
                # selected resIDs in the current segment (total time doesn't change if in1d(self.resIDs[masks_segment], self.resIDs[masks_segment*masks0]) withouth applying unique)
                resIDs_unique_segment = unique(self.resIDs[masks_segment*masks0])
                # any other atoms in the current segment belonging to the selected resIDs
                masks[masks_segment] += isin(self.resIDs[masks_segment], resIDs_unique_segment)
            if not mute:
                print("\n {} atoms have been selected by radius {:.3f} a.u. with all boundary residues included".format(masks.sum(), radius*unitconvert.unit2au.get(unit, 1.0)))
            
        # boundary residues excluded in selection
        elif str(boundary)[:4].lower() == 'excl':
            # YL 09/09/2021: improved algorithm 200 faster than before
            for i, segment in enumerate(segments_unique):
                if segment not in segments_unique_selected:
                    continue
                # select the current segment: segment_inverse==i is faster than self.segments==segment
                masks_segment = segment_inverse==i
                # resIDs of unselected atoms in the current segment
                resIDs_unique_segment = unique(self.resIDs[masks_segment*~masks0])
                # remove resIDs found in the unselected atoms
                masks[masks_segment] *= isin(self.resIDs[masks_segment], resIDs_unique_segment, invert=True)
            if not mute:
                print("\n {} atoms have been selected by radius {:.3f} a.u. with all boundary residues excluded".format(masks.sum(), radius*unitconvert.unit2au.get(unit, 1.0)))
        else:
            print(" >>> FragmentBio.selectByRadius(): invalid `residues` argument value: {}".format(boundary))

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


    # NB: the built-in Python function range() is overridden!
    def selectByRadius(self, radius, range=[], centre=None, boundary='inclusive', unit='a.u.', return_masks=False, mute=True):
        '''Biological version of FragmentBase.selectByRadius()'''

        from numpy import all, arange, argsort, argwhere, copy, isin, unique

        # boundary=None to skip including/excluding fragments
        masks = _FragmentBase.selectByRadius(self, radius, range=range, centre=centre, unit=unit, boundary=None, return_masks=True)

        # YL 17/01/2022: this function has been made an independent one
        masks = self.selectByBoundary(masks, boundary=boundary, unit=unit, return_masks=True, mute=mute)

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


    # YL: do NOT use! this method is WRONG! there is no easy way to describe lone pairs or charge
    def selectByHybridisation(self, hybridisation, connected_hydrogens=True, return_masks=False, invert=False, **kwargs):
        '''Select atoms by types of hybridisation'''

        from numpy import arange, isin, sum, where, unique, zeros

        # C
        masks_carbon   = self.selectByAtoms( 6, return_masks=True)
        # N
        masks_nitrogen = self.selectByAtoms( 7, return_masks=True)
        # O
        masks_oxygen   = self.selectByAtoms( 8, return_masks=True)
        # S
        masks_sulphur  = self.selectByAtoms(16, return_masks=True)

        masks_selected = zeros((self.natoms,), dtype=bool)

        if hybridisation.strip().lower() in [ 'sp3' ]:
            masks_selected[masks_carbon]   += sum(self.conn[masks_carbon]   > -1, axis=1) == 4
            masks_selected[masks_nitrogen] += sum(self.conn[masks_nitrogen] > -1, axis=1) == 3
            masks_selected[masks_oxygen]   += sum(self.conn[masks_oxygen]   > -1, axis=1) == 2
            masks_selected[masks_sulphur]  += sum(self.conn[masks_sulphur]  > -1, axis=1) == 2

        elif hybridisation.strip().lower() in [ 'sp2' ]:
            masks_selected[masks_carbon]   += sum(self.conn[masks_carbon]   > -1, axis=1) == 3
            masks_selected[masks_nitrogen] += sum(self.conn[masks_nitrogen] > -1, axis=1) == 3
            masks_selected[masks_oxygen]   += sum(self.conn[masks_oxygen]   > -1, axis=1) == 2
            masks_selected[masks_sulphur]  += sum(self.conn[masks_sulphur]  > -1, axis=1) == 2

        elif hybridisation.strip().lower() in [ 'sp' ]:
            masks_selected[masks_carbon]   += sum(self.conn[masks_carbon]   > -1, axis=1) == 2
            masks_selected[masks_nitrogen] += sum(self.conn[masks_nitrogen] > -1, axis=1) == 1

        else:
            if chemsh.CHEMSH_DEBUG > 4:
                print("\n WARNING: no atoms selected for hybridisation type \"{}\"\n".format(hybridisation))

        # add connected H atoms
        if connected_hydrogens:
            indices_selected_connected_atoms = unique(self.conn[masks_selected][self.conn[masks_selected]>-1])
            indices_selected_connected_hydrogens = indices_selected_connected_atoms[self.znums[indices_selected_connected_atoms]==1]
            masks_selected += isin(arange(self.natoms), indices_selected_connected_hydrogens)

        # invert selection
        if invert:
            masks_selected = ~masks_selected

        if return_masks:
            return masks_selected

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


    @miscutils.time
    def detectNucleotides(self, aliases={}, delete_terminal_hydrogen=False, delete_terminal_phosphate=False):
        '''Tell if each residue is a deoxyribonucleic (DNA) or ribonucleic (RNA) acid'''

        from numpy import arange, char, unique, zeros

        indices_all = arange(self.natoms)
        names       = char.decode(char.strip(self.names))
        znums       = self.znums
        ids_segment = unique(self.all_segments_inverse_indices)

        masks_terminal_h         = zeros((self.natoms,), dtype=bool)
        masks_terminal_phosphate = zeros((self.natoms,), dtype=bool)

        has_nucleotide = False

        # loop over segments
        for id_segment in ids_segment:

            segname = self.all_segments[id_segment]
            self.optimiseResIDs(segname)

            masks_segment = self.all_segments_inverse_indices == id_segment
            ids_residue = unique(self.resIDs[masks_segment])
            id_residue_max = ids_residue.max()
            # YL 07/02/2022: determine to non-nucleotide if no chain "A"
            try:
                id_residue_max_A = self.resIDs[self.chainIDs==b'A'].max()
            except:
                return False

            for id_residue in ids_residue:
                masks_residue = self.resIDs == id_residue
                masks = masks_segment*masks_residue
                formula = self.getSelected(indices_all[masks], flatten=False).formula

                # 3'- or 5'-terminal
#                masks_c3p = masks*(names=="C3'")
#                index_c3p = indices_all[masks_c3p][0]
                masks_o3p = masks*(names=="O3'")
                conn_o3p = self.conn[masks_o3p]
#                masks_phosphodiester = (names=="P")+(names=="O1P")+(names=="O2P")+(names=="O5'")
#                masks_phosphodiester *= masks
#                # include the neighbour O3' except for terminal bases
#                if id_residue not in [1, id_residue_max_A, id_residue_max_A+1, id_residue_max_A*2]:
#                    masks_phosphodiester += (names=="O3'")*(self.resIDs==id_residue-1)
#                conn_phosphodiester = self.conn[masks_phosphodiester]
#                znums_connected_phosphodiester = znums[conn_phosphodiester[conn_phosphodiester>-1]]
#                ncarbons_connected_phosphodiester = znums_connected_phosphodiester.tolist().count(6)
                suffix = ''
#                if ncarbons_connected_phosphodiester != 2:
                # NB: assuming resIDs are correct
                if id_residue in [1, id_residue_max_A, id_residue_max_A+1, id_residue_max_A*2]:
                    # 5'
                    znums_connected_o3p = znums[conn_o3p[conn_o3p>-1]]
                    if 15 in znums_connected_o3p or set(znums_connected_o3p) == set():
                        suffix = '5'
                        masks_terminal_phosphate += masks*(names=="O3'")
                    # 3'
                    else:
                        suffix = '3'
                        masks_terminal_h += masks*(names=="H")

                # DNA or RNA
                masks_C2p = names == "C2'"
                conn_C2p = self.conn[masks_C2p*masks]
                znums_connected_C2p = znums[conn_C2p[conn_C2p>-1]]
                if 8 in znums_connected_C2p:
                    prefix = ''
                else:
                    prefix = 'd'

                resname = self.residues[masks][0]
                # A or G
                if 'C10' in formula and 'N5' in formula:
                    masks_C2 = names == 'C2'
                    conn_C2 = self.conn[masks_C2*masks]
                    znums_connected_C2 = znums[conn_C2[conn_C2>-1]]
                    masks_C6 = names == 'C6'
                    conn_C6 = self.conn[masks_C6*masks]
                    znums_connected_C6 = znums[conn_C6[conn_C6>-1]]
                    # G
                    if set(znums_connected_C2) == {7} and set(znums_connected_C6) == {6,7,8}:
                        resname = f'{prefix}G{suffix}'
                        has_nucleotide = True
                    elif set(znums_connected_C2) == {1,7} and set(znums_connected_C6) == {6,7}:
                        resname = f'{prefix}A{suffix}'
                        has_nucleotide = True
                # T
                elif 'C10' in formula and 'N2' in formula:
                    masks_C4 = names == 'C4'
                    conn_C4 = self.conn[masks_C4*masks]
                    znums_connected_C4 = znums[conn_C4[conn_C4>-1]]
                    masks_C5 = names == 'C5'
                    conn_C5 = self.conn[masks_C5*masks]
                    znums_connected_C5 = znums[conn_C5[conn_C5>-1]]
                    if set(znums_connected_C4) == {6,7,8} and set(znums_connected_C5) == {6}:
                        resname = f'{prefix}T{suffix}'
                        has_nucleotide = True
                # C
                elif 'C9' in formula and 'N3' in formula:
                    masks_C4 = names == 'C4'
                    conn_C4 = self.conn[masks_C4*masks]
                    znums_connected_C4 = znums[conn_C4[conn_C4>-1]]
                    masks_C5 = names == 'C5'
                    conn_C5 = self.conn[masks_C5*masks]
                    znums_connected_C5 = znums[conn_C5[conn_C5>-1]]
                    if set(znums_connected_C4) == {6,7} and set(znums_connected_C5) == {1,6}:
                        resname = f'{prefix}C{suffix}'
                        has_nucleotide = True
                # U
                elif 'C9' in formula and 'N2' in formula:
                    masks_C4 = names == 'C4'
                    conn_C4 = self.conn[masks_C4*masks]
                    znums_connected_C4 = znums[conn_C4[conn_C4>-1]]
                    masks_C5 = names == 'C5'
                    conn_C5 = self.conn[masks_C5*masks]
                    znums_connected_C5 = znums[conn_C5[conn_C5>-1]]
                    if set(znums_connected_C4) == {6,7,8} and set(znums_connected_C5) == {1,6}:
                        resname = f'{prefix}T{suffix}'
                        has_nucleotide = True

                self.residues[masks] = aliases.get(resname, resname)

            # delete terminal atoms as requested
            indices_to_delete = []
            if delete_terminal_hydrogen:
                indices_terminal_h = indices_all[masks_terminal_h].tolist()
                indices_to_delete += indices_terminal_h
            if delete_terminal_phosphate:
                indices_terminal_phosphate = indices_all[masks_terminal_phosphate].tolist()
                indices_to_delete += indices_terminal_phosphate
            self.delete(indices_to_delete)

        return has_nucleotide


    def insertDummy(self, recalculate=True):
        '''Insert three dummy atoms to calculate the Z-matrix'''

        from numpy import nan

        # insert dummy atoms and give them coordinates
        self.insert(['X','X','X'],0)
        self.coords[0] = 0.000, 0.000, 0.000
        self.coords[1] = 0.000, 0.000, 1.000
        self.coords[2] = 0.000, 1.000,-0.000
        # define the dummy atoms' internal coordinates
        self._master.zmatrix[0] = -1, nan               , -1,  nan, -1,  nan
        self._master.zmatrix[1] =  0, 1.88972666351     , -1,  nan, -1,  nan
        self._master.zmatrix[2] =  1, 2.6724770767139003,  0, 45.0, -1,  nan
        # define the relative position of the first original atom: 2.83459 Bohr = 1.5 Angstrom which is an arbitrary distance
        self._master.zmatrix[3] =  0, 2.834589995265    ,  2, 90.0,  1,-90.0
        # if there is a second original atom
        if self.natoms > 4:
            self._master.zmatrix[4] = 3, self.dist(4,3), 0, 180.0, 1, 0.0
        # if there is a third original atom
        if self.natoms > 5:
            self._master.zmatrix[5] = 4, self.dist(5,4), 3, self.angle(5,4,3), 0, self.dihedral(5,4,3,0)

        # by now we have already got the Z-matrix we want and the following calculation will be redundant
        # but to keep the structure of the code unchanged we recalculate the Catesian cooridnates and connectivity
        if recalculate:
            self.z2c()
            self.connect(force_reconnect=True, mute=True)


    # YL 24/04/2019
#    @miscutils.time
    def c2z(self, atom0=-1, atom1=-1, atom2=-1):
        '''Cartesian to Z-matrix conversion'''

        from numpy import all, any, arange, argmax, concatenate, delete, copy, flip, isin, isin, argwhere, nan, nonzero, ones, set_printoptions

        #set_printoptions(threshold=1000000)

        # YL 28/12/2021: we need at least 4 atoms to represent the Z-matrix
        if self.natoms < 4:
            self.insertDummy()

        # natoms >= 4
        else:
            # all bonds, angles, and dihedrals including inversed redundancies
            def setBADs(self):
                self.__bonds__     = concatenate((flip(self.bonds    ,1), self.bonds)    , axis=0)
                self.__angles__    = concatenate((flip(self.angles   ,1), self.angles)   , axis=0)
                self.__dihedrals__ = concatenate((flip(self.dihedrals,1), self.dihedrals), axis=0)
    
            def getJ(self, _i, readonly=False, skip_H=False, favoured_J=-1):
    
                # skip atoms of coordination number being 1, such as hydrogens
                if skip_H:
                    #nonhydrogens = self.znums[self.__bonds__[:,1]] > 1
                    nonhydrogens = self.ncoordinations[self.__bonds__[:,1]] > 1
                else:
                    nonhydrogens = ones(shape=(self.__bonds__.shape[0],), dtype=int)
    
                # look for a specific J
                if favoured_J > -1:
                    favoured = self.__bonds__[:,1] == favoured_J
                else:
                    favoured = ones(shape=(self.__bonds__.shape[0],), dtype=int)
    
                try:
                    # get the first match of current atom as the leading centre
                    matches = (self.__bonds__[:,0]==_i)&(nonhydrogens)&(favoured)
                    idxl = not any(matches) or argmax(matches)
                    j = self.__bonds__[idxl][1]
                    # remove the used bonds (and its inverse to prevent redundant dependency)
                    # benchmarked: it won't slow down the loops
                    if not readonly:
                        self.__bonds__ = delete(self.__bonds__, idxl, axis=0)
                        # there has got to be a match
                        idxr = argmax((self.__bonds__[:,0]==j)
                                     &(self.__bonds__[:,1]==_i))
                        self.__bonds__ = delete(self.__bonds__, idxr, axis=0)
                # get the closest atom as J if no bond
                except:
                    j = self.getClosest(_i)
    
                return j
    
            def getJK(self, _i, readonly=False, skip_H=True, favoured_J=-1, favoured_K=-1):
    
                # skip Js of coordination number being 1, such as hydrogens
                if skip_H:
                    #nonhydrogens = self.znums[self.__angles__[:,1]] > 1
                    nonhydrogens = self.ncoordinations[self.__angles__[:,1]] > 1
                else:
                    nonhydrogens = ones(shape=(self.__angles__.shape[0],), dtype=int)
    
                # look for a specific J or K
                if favoured_J > -1 and favoured_K == -1:
                    favoured = self.__angles__[:,1] == favoured_J
                elif favoured_J == -1 and favoured_K > -1:
                    favoured = self.__angles__[:,1] == favoured_K
                elif favoured_J > -1 and favoured_K > -1:
                    favoured = (self.__angles__[:,1] == favoured_J) & (self.__angles__[:,2] == favoured_K)
                else:
                    favoured = ones(shape=(self.__angles__.shape[0],), dtype=int)
    
                # get the first bonded angle from the remaining list of angles
                try:
                    # get the first match of current atom as the leading centre
                    matches = (self.__angles__[:,0]==_i)&(nonhydrogens)&(favoured)
                    # raise error when idxl is True
                    idxl = not any(matches) or argmax(matches)
                    j, k = self.__angles__[idxl][1:]
                    # remove the used angles (and its inverse to prevent redundant dependency)
                    # benchmarked: it won't slow down the loops
                    if not readonly:
                        self.__angles__ = delete(self.__angles__, idxl, axis=0)
                        # there has got to be a match
                        idxr = argmax((self.__angles__[:,0]==k)
                                     &(self.__angles__[:,1]==j)
                                     &(self.__angles__[:,2]==_i))
                        self.__angles__ = delete(self.__angles__, idxr, axis=0)
                # get the closest K if no bonded angle
                except:
                    j = getJ(self, _i)
                    k = self.getClosest(j, exclude=_i)
    
                return j, k
    
            def getL(self, _i, _j):
                '''Get the 4th centre L'''
    
                # we need to find a 4th centre L by searching I-J-L, if there is no bonded dihedral via I-J-K
                #     (M)
                #   K--J--I
                #     /
                #    L
                m, l = getJK(self, _i, readonly=True)
                while m != _j:
                    m, l = getJK(self, _i, readonly=True)
    
                return l
    
    
            def getJKL(self, _i, readonly=False, skip_H=True):
                '''Get the latter thre centres of a dihedral'''
    
                # skip Js and Ks of coordination number being 1, such as hydrogens
                if skip_H:
                    #nonhydrogens = (self.znums[self.__dihedrals__[:,1]] > 1) & (self.znums[self.__dihedrals__[:,2]] > 1)
                    nonhydrogens = (self.ncoordinations[self.__dihedrals__[:,1]] > 1) & (self.ncoordinations[self.__dihedrals__[:,2]] > 1)
                else:
                    nonhydrogens = ones(shape=(self.__dihedrals__.shape[0],), dtype=int)
    
                # bonded
                try:
                    # get the first match of current atom as the leading centre
                    # NB: returns: the index (including 0) if matches; True if no match
                    # YL 12/02/2020: there is a risk in this method that the starting end may not take part in a dihedral
                    #                but only form impropers. but since self.__dihedrals__ does not contain improper,
                    #                the fourth atom (index 3) has to be connected to index>3 atoms to form a dihedral.
                    #                for example, in:
                    #                  O0     O3
                    #                    \   /
                    #                O2 -- P1
                    #                      |
                    #                      O4--C5--...
                    #                Z-matrix of O3 will be (1 x.xx 4 xx.xx 5 xx.xx) and in this case z2c() cannot be calculated
                    matches = (self.__dihedrals__[:,0]==_i)&(nonhydrogens)
                    idxl = not any(matches) or argmax(matches)
                    j, k, l = self.__dihedrals__[idxl][1:]
                    # remove the used dihedrals (and its inverse to prevent redundant dependency)
                    # benchmarked: it won't slow down the loops
                    if not readonly:
                        self.__dihedrals__ = delete(self.__dihedrals__, idxl, axis=0)
                        # there has got to be a match
                        idxr = argmax((self.__dihedrals__[:,0]==l)
                                     &(self.__dihedrals__[:,1]==k)
                                     &(self.__dihedrals__[:,2]==j)
                                     &(self.__dihedrals__[:,3]==_i))
                        self.__dihedrals__ = delete(self.__dihedrals__, idxr, axis=0)
                # get improper if possible when no bonded dihedral
                except:
                    j, k = getJK(self, _i)
                    l = getL(self, _i, j)
    
                return j, k, l
    
            # get pool of bonds, angles, and dihedrals
            setBADs(self)
    
            calculated = zeros(shape=(self.natoms,), dtype=int)
    
            # 1st atom: skip atoms of coordination number being 1 such as hydrogens
            if not atom0 > -1:
                atom0 = argmax((self.ncoordinations>1)&(calculated==0))
            self._master.zmatrix[atom0] = -1, nan, -1, nan, -1, nan
            calculated[atom0] = 1
    
            # 2nd atom: skip atoms of coordination number being 1 such as hydrogens
            if not atom1 > -1:
                atom1 = argmax((self.ncoordinations>1)&(calculated==0)&any(self.conn==atom0,axis=1))
            j = getJ(self, atom1, favoured_J=atom0)
            self._master.zmatrix[atom1] =  j, self.dist(atom1,j), -1, nan, -1, nan
            calculated[atom1] = 1
    
            # 3rd atom: skip atoms of coordination number being 1 such as hydrogens
            if not atom2 > -1:
                atom2 = argmax((self.ncoordinations>1)&(calculated==0)&any(self.conn==atom1,axis=1))
            j, k = getJK(self, atom2, favoured_J=atom1, favoured_K=atom0)
            self._master.zmatrix[atom2] =  j, self.dist(atom2,j), k, self.angle(atom2,j,k), -1, nan
            calculated[atom2] = 1
    
            # further atoms
            indices_atoms = arange(self.natoms)
            for i in indices_atoms[calculated==0]:
                j, k, l = getJKL(self, i)
                self._master.zmatrix[i] =  j, self.dist(i,j),  k, self.angle(i,j,k), l, self.dihedral(i,j,k,l)
    
            # non-bonded fragments
            setBADs(self)
            for fragment in arange(1, self.nfragments):
                i = argmax(self.fragIDs==fragment)
                # get the closest atom from Segment 0
                j = self.getClosest(i, select=self._indices_all[self.fragIDs==0])
                k, l = getJK(self, j)
                if chemsh.CHEMSH_DEBUG > 3:
                    print("\n >>> c2z(): resetting the leading atom {} of fragment {}\n".format(i, fragment))
                self._master.zmatrix[i] =  j, self.dist(i,j),  k, self.angle(i,j,k), l, self.dihedral(i,j,k,l)



            # debug print
            if chemsh.CHEMSH_DEBUG > 3:
                print("\n >>> Converting Cartesians to Z-matrix:\n")
                for i in arange(self.natoms):
                    j = self._master.zmatrix[i].j
                    k = self._master.zmatrix[i].k
                    l = self._master.zmatrix[i].l
                    b = self._master.zmatrix[i].bond
                    a = self._master.zmatrix[i].angle
                    d = self._master.zmatrix[i].dihedral
                    if i == atom0:
                        print(" {:5d} {:8s}  {:>5} {:<8} {:>6}  {:>5} {:<8} {:>6}  {:>5} {:<8} {:>6}".format(atom0, self.names[i].decode(),
                                                                                                              '--', '--', '--',
                                                                                                              '--', '--', '--',
                                                                                                              '--', '--', '--'))
                    elif i == atom1:
                        print(" {:5d} {:8s}  {:5d} {:8s} {:6.3f}  {:>5} {:<8} {:>6}  {:>5} {:<8} {:>6}".format(atom1,  self.names[i].decode(),
                                                                                                                   j,  self.names[j].decode(), b,
                                                                                                                '--', '--', '--',
                                                                                                               '--', '--', '--'))
                    elif i == atom2:
                        print(" {:5d} {:8s}  {:5d} {:8s} {:6.3f}  {:5d} {:8s} {:6.2f}  {:>5} {:<8} {:>6}".format(atom2,  self.names[i].decode(),
                                                                                                                     j,  self.names[j].decode(), b,
                                                                                                                    k,  self.names[k].decode(), a,
                                                                                                                  '--', '--', '--'))
                    else:
                        print(" {:5d} {:8s}  {:5d} {:8s} {:6.3f}  {:5d} {:8s} {:6.2f}  {:5d} {:8s} {:6.2f}".format( i, self.names[i].decode(),
                                                                                                                    j, self.names[j].decode(), b,
                                                                                                                    k, self.names[k].decode(), a,
                                                                                                                    l, self.names[l].decode(), d))

            # YL 03/05/2024: we do a drill of z2c here to make sure back transformation
            #                can happen especially for structures with rings; this procedure
            #                will update Z-matrix if it finds it necessary but doesn't not
            #                amend the Cartesian coordinates
            self.z2c(coords_atom0=self.coords[atom0],
                     coords_atom1=self.coords[atom1],
                     coords_atom2=self.coords[atom2],
                     update_coords=False)

        # don't need to recalculate in the future
        self._calcZMatrix = False
    # TODO: dummy atoms are needed if there are ring sturctures (bonds are left over)


    def z2c(self, inverse_z=False, coords_atom0=None, coords_atom1=None, coords_atom2=None, update_coords=True):
        '''Z-matrix to cartesian conversion (NeRF algorithm)
           see: Parsons, et al., J. Comput. Chem. 2005, 26, 1063-1068
        '''

        from numpy import all, arange, argmin, argmax, array, array_equal, cos, cross, dot
        from numpy import linalg, nonzero, pi, sin, stack, unique, zeros
        from math  import radians, cos, sin
        from warnings import filterwarnings, resetwarnings, simplefilter

        # raise error if coordinates are wrong
        simplefilter('error')

        calculated = zeros(shape=(self.natoms,), dtype=int)

        def _NeRF(self, _i, _j, _k, _l):
            '''An implementation of the Self-Normalizing Natural Extension Reference Frame (SN-NeRF) method'''

            # normalised vectors
            kl = self.coords[_k]-self.coords[_l]
            kl = kl/linalg.norm(kl)
            jk = self.coords[_j]-self.coords[_k]
            jk = jk/linalg.norm(jk)
       
            # z and y vectors
            # see: eq. (1)
            z = cross(kl, jk)
            z = z/linalg.norm(z)
            y = cross(z, jk)
       
            # rotation matrix
            # see: eq. (5)
            M = stack([jk, y, z], axis=-1)
       
            # bond, angle (in rad), and dihedral (in rad)
            b = self._master.zmatrix[_i].bond
            a = self._master.zmatrix[_i].angle
            d = self._master.zmatrix[_i].dihedral
            a *= pi/180.0
            d *= pi/180.0
       
            # apply the rotation
            return self.coords[_j] + dot(M, [ -b*cos(a), b*sin(a)*cos(d), b*sin(a)*sin(d)])


        # recursively determine the atoms from 4th on
        def _determineSequence(self, _calculated):
            ''''''

            Z = self._master.zmatrix

            # YL 02/05/2024
            def _findNewNeighbour(_self, _i, __calculated):
                '''Recursively find a new neighbour if _i's current neighbour is bad'''

                _Z = _self._master.zmatrix
                neighbours = _self.conn[_i][:argmin(_self.conn[_i])]
                neighbours_calculated = neighbours[__calculated[neighbours]==1]
                for _j in neighbours_calculated:
                    # redefine the chain
                    _k, _l = _Z[_j].j, _Z[_j].k
                    if all(__calculated[[_j,_k,_l]]):
                        # update the Z-matrix
                        if not update_coords:
                            _Z[_i].j, _Z[_i].k, _Z[_i].l = _j, _k, _l
                            _Z[_i].bond     = self.dist    (_i,_j)
                            _Z[_i].angle    = self.angle   (_i,_j,_k)
                            _Z[_i].dihedral = self.dihedral(_i,_j,_k,_l)
                            __calculated[_i] = 1
                        break
                    # _j is bad and it needs a better neeighbour for itself
                    else:
                        # we shouldn't touch the three first determined atoms
                        if _j in [atom0,atom1,atom2]:
                            continue
                        _findNewNeighbour(_self, _j, __calculated)

            # loop over the atoms of those the Cartesian positions haven't been determined
            for i in self._indices_all[_calculated==0]:
                j, k, l = Z[i].j, Z[i].k, Z[i].l
                # we can use NeRF to calculate the position of I if J, K, and L are already determined
                if all(_calculated[[j,k,l]]):
                    if update_coords:
                        self.coords[i] = _NeRF(self,i,j,k,l)
                    _calculated[i] = 1
                else:
                    # YL 02/05/2024: we must change I's beta- and gamma-neighbours using the determined J (update_coords==False only!)
                    if _calculated[j]:
                        k, l = Z[j].j, Z[j].k
                        # update the Z-matrix
                        if all(_calculated[[j,k,l]]) and not update_coords:
                            Z[i].k, Z[i].l = k, l
                            Z[i].angle    = self.angle   (i,j,k)
                            Z[i].dihedral = self.dihedral(i,j,k,l)
                            _calculated[i] = 1
                    # YL 03/05/2024: if j is undetermined
                    else:
                        # we shouldn't touch the three first determined atoms
                        if i in [atom0,atom1,atom2]:
                            continue
                        # recursively find a new neighbour for i
                        _findNewNeighbour(self, i, _calculated)
            # finished
            if all(_calculated):
                return 0
            # start over if there are remaining undetermined positions
            else:
                _determineSequence(self, _calculated)

        # if the first three atoms aren't defined try to determine the positions, otherwise use the existing Cartesians

        # determin three atoms' positions first that are not necessarily the first three ones
        atom0 = argmax((self._master.zmatrix[:].j==-1)&(self._master.zmatrix[:].k==-1)&(self._master.zmatrix[:].l==-1))
        atom1 = argmax((self._master.zmatrix[:].j!=-1)&(self._master.zmatrix[:].k==-1)&(self._master.zmatrix[:].l==-1))
        atom2 = argmax((self._master.zmatrix[:].j!=-1)&(self._master.zmatrix[:].k!=-1)&(self._master.zmatrix[:].l==-1))

        # atom 0
        if coords_atom0 is not None:
            self.coords[atom0] = coords_atom0[:]
        else:
            self.coords[atom0] = 0.0, 0.0, 0.0

        # atom 1
        if self._master.zmatrix[atom1].j == atom0:
            if inverse_z:
                if coords_atom1 is not None:
                    self.coords[atom1] = coords_atom1[:]
                else:
                    self.coords[atom1] = 0.0, 0.0, -self._master.zmatrix[atom1].bond
            else:
                if coords_atom1 is not None:
                    self.coords[atom1] = coords_atom1[:]
                else:
                    self.coords[atom1] = 0.0, 0.0, self._master.zmatrix[atom1].bond
        else:
            pass

        # atom 2
        j = self._master.zmatrix[atom2].j
        k = self._master.zmatrix[atom2].k
        b = self._master.zmatrix[atom2].bond
        a = radians(self._master.zmatrix[atom2].angle) # degrees to radians
        # I--J--K(origin)
        if j == atom1 and k == atom0:
            if inverse_z:
                if coords_atom2 is not None:
                    self.coords[atom2] = coords_atom2[:]
                else:
                    self.coords[atom2] = 0.0, b*sin(a), b*cos(a)+self.coords[atom1][2]
            else:
                if coords_atom2 is not None:
                    self.coords[atom2] = coords_atom2[:]
                else:
                    self.coords[atom2] = 0.0, b*sin(a), -b*cos(a)+self.coords[atom1][2]
        # I--J(origin)--K
        elif j == atom0 and k == atom1:
            if coords_atom2 is not None:
                self.coords[atom2] = coords_atom2[:]
            else:
                self.coords[atom2] = 0.0, b*sin(a), -b*cos(a)

        calculated[atom0] = calculated[atom1] = calculated[atom2] = 1

        _determineSequence(self, calculated)

        # reset warning setting
        simplefilter('default')
        resetwarnings()

        return self.coords


    def reorientDihedral(self, iatom, jatom, katom, latom):
        '''Re-orient Z-matrix of an atom'''

        self._master['zmatrix']['j']       [iatom] = jatom
        self._master['zmatrix']['k']       [iatom] = katom
        self._master['zmatrix']['l']       [iatom] = latom
        self._master['zmatrix']['bond']    [iatom] = self.dist    (iatom, jatom)
        self._master['zmatrix']['angle']   [iatom] = self.angle   (iatom, jatom, katom)
        self._master['zmatrix']['dihedral'][iatom] = self.dihedral(iatom, jatom, katom, latom)


    def printZMatrix(self, start=0):
        '''Print formatted Z-matrix'''

        print("\n Z-matrix (in a.u. and degrees):\n")
        fmt = '{:>4s} {:<6s} {:>4s} {:<6s} {:>9s} {:>4s} {:<6s} {:>9s} {:>4s} {:<6s} {:>9s}\n'
        print(fmt.format('I', 'name', 'J', 'name', 'bond', 'K', 'name', 'angle', 'L', 'name', 'dihedral'))
        fmt = '{:>4d} {:<6s} {:>4} {:<6s} {:9.3f} {:>4} {:<6s} {:9.3f} {:>4} {:<6s} {:9.3f}'
        for i in range(self.natoms):
            j = self._master['zmatrix']['j'][i]
            k = self._master['zmatrix']['k'][i]
            l = self._master['zmatrix']['l'][i]
            iname = self.names[i].decode().strip()
            if j == -1:
                jname = jatom = '-'
            else:
                jname = self.names[j].decode().strip()
                jatom = j+start
            if k == -1:
                kname = katom = '-'
            else:
                kname = self.names[k].decode().strip()
                katom = k+start
            if l == -1:
                lname = latom = '-'
            else:
                lname = self.names[l].decode().strip()
                latom = l+start
            b = self._master['zmatrix']['bond'][i]
            a = self._master['zmatrix']['angle'][i]
            d = self._master['zmatrix']['dihedral'][i]
            print(fmt.format(i+start, iname, jatom, jname, b, katom, kname, a, latom, lname, d))
            stdout.flush()


    def split(self, save_to_pdb=True, exclude=[], groups=[], return_dict=True, sort=False):
        '''Split into segments, save to separate files, and return a list of the filenames'''

        from collections import OrderedDict
        from numpy       import arange, array, char, recarray, unique, zeros
        from numpy       import sort as np_sort
        from ..io        import format

        segments, unique_indices = unique(self.segments, return_index=True)
        # an unsorted list
        segments = self.segments[np_sort(unique_indices)].tolist()

        def _inSubgroup(_seg, _subgroup):
            ''''''
            try:
                _subgroup = char.strip(array(_subgroup))
            except:
                pass
            return _seg.decode() in _subgroup or _seg.decode().strip() in _subgroup

        # NB: list is not hashable while set is unordered
        segments = [ tuple(groups)
                         if seg.decode() in groups or seg.decode().strip() in groups else
                     tuple(*[subgroup for subgroup in groups if _inSubgroup(seg,subgroup)])
                         if any([True for subgroup in groups if _inSubgroup(seg,subgroup)]) else
                     seg for seg in segments ]

        # remove repeated elements but maintain the order (we can't use set() which is unsorted or unique() which doesn't work here)
        segments = list(OrderedDict.fromkeys(segments))

        # a sorted list (list.sort() won't work)
        if sort:
            lbuffs = [ i for i in segments if type(i) in [ str, bytes ] ]
            lbuffs.sort()
            lbufft = [ i for i in segments if type(i) is tuple ]
            segments = lbuffs + lbufft

        if len(segments) > 1:
            print("\n Splitting into segments:")
            format.printwrap(' '.join(char.strip([seg.decode() if hasattr(seg,'decode') else
                             '_'.join(seg) for seg in segments])),
                             prepend_blanc_line=False)

        dict_segFrags = dict()
        dtype = []

        # loop over segments
        for seg in segments:
            try:
                segstr = seg.decode().strip()
            except:
                segstr = '_'.join([s.strip() for s in seg])
            # select by segment name
            masks_segment = self.segments == seg
            masks_protein = self.selectByProtein(return_masks=True)
            # invert select by residue names to exclude
            excluded = self.select(residues=exclude, invert=True, return_masks=True)
            # get a new fragment from the selected segment
            # YL 25/02/2021: we can't use masks array in FragmentBase.getSelected()
            segFrag = self.getSelected(self._indices_all[masks_segment*excluded])
            # YL 26/05/2024: all segments should be split per chain (not just for PSFGEN but also as a generic splitting function)
            if masks_segment.sum():
                masks_selected = zeros((segFrag.natoms), dtype=bool)
                chainIDs = unique(segFrag.chainIDs)
                for chainID in chainIDs:
                    masks_chainID = segFrag.chainIDs == chainID
                    masks_selected += masks_chainID
                    protFrag = segFrag.getSelected(segFrag._indices_all[masks_chainID])
                    sbuff = segstr + '_' + chainID.decode().strip()
                    protFrag.title = sbuff
                    if save_to_pdb:
                        seg_filename = f'_segment_{sbuff}.pdb'
                        protFrag.save(seg_filename)
                    dict_segFrags.update({sbuff:protFrag})
                    dtype += [ (sbuff,object) ]
                # Yl 04/05/2024: I don't have an example yet but I'm afraid some misc
                #                residues in the protein segment could be left over
                if (~masks_selected).sum():
                    protFrag = segFrag.getSelected(segFrag._indices_all[~masks_selected])
                    sbuff = segstr + '_MISC'
                    protFrag.title = sbuff
                    if save_to_pdb:
                        seg_filename = f'_segment_{sbuff}.pdb'
                        protFrag.save(seg_filename)
                    dict_segFrags.update({sbuff:protFrag})
                    dtype += [ (sbuff,object) ]
            else:
                segFrag.title = segstr
                # save to a new PDB file
                if save_to_pdb:
                    seg_filename = f'_segment_{segstr}.pdb'
                    segFrag.save(seg_filename)
                dict_segFrags.update({segstr:segFrag})
                dtype += [ (segstr,object) ]

        if return_dict:
            return dict_segFrags
        else:
            arr_segFrags = recarray((1,), dtype=dtype)
            for k, v in dict_segFrags.items():
                arr_segFrags[k][0] = v
            return arr_segFrags


    def getSelected(self, selected, flatten=True):
        '''Overridden _FragmentBase.getSelected for biosimulation purpose'''


        inst = _FragmentBase.getSelected(self, selected, flatten=flatten)

        # YL 26/10/2021: we must perform a reordering because by selecting the new fragment's segments are easily made discontiguous
#        for segment in inst.all_segments:
#            inst.keepSegmentContiguous(segment)
#            inst.optimiseResIDs(segment)

        # YL 26/11/2021: inherit renaming history
        inst._renamed_residues.update(self._renamed_residues)

        return inst


    def getDisulphides(self, returnResIDs=False, unique=True):
        '''Find disulphides (e.g., CYS--CYS)'''

        from numpy import any, array

        # list of all sulphurs
        sulphurs = self.select(atoms=16)

        # return if no sulphur
        if not sulphurs.size:
            return array([])

        # znums of the connected atoms (including conn=-1)
        znums = self.znums[self.conn[sulphurs]]

        # znums those are 16
        znums_sulphurs = znums==16

        # exclude conn=-1 items
        znums_sulphurs *= self.conn[sulphurs]>-1

        # any connected sulphurs?
        isulphurs = sulphurs[any(znums_sulphurs, axis=1)]
        jsulphurs = self.conn[sulphurs][znums_sulphurs]

        disulphides = array([isulphurs,jsulphurs]).T
        # return if no sulphur
        if not disulphides.size:
            return array([])

        # remove inverted duplicates
        if unique: 
            disulphides = { tuple(sorted(i)) for i in disulphides }
            disulphides = array(list(disulphides))

        # if requested to return resIDs containing disulphides
        if returnResIDs:
            disulphides = array([self.resIDs[disulphides.T[0]], self.resIDs[disulphides.T[1]]]).T

        return disulphides


    def residue2resIDs(self, residue):
        '''Find all resIDs given a residue name'''

        from numpy import unique

        resIDs = unique(self.resIDs[self.select(residues=residue)])

        return resIDs


    def getConnectedResIDsFromResID(self, resID, segment=''):
        '''Return residue(s) connected to the given resID'''

        from numpy import char, isin, unique

        # atom indices of the current residue
        indices_residue = self.select(resIDs=resID, segments=segment)

        # indices of all atoms connected to the current residue
        conn = self.conn[indices_residue].ravel()

        # exclude atoms of self
        indices_connected = unique(conn[conn>-1])
        indices_connected = indices_connected[~isin(indices_connected, indices_residue)]

        # must return the segments at the same time to ensure the uniqueness
        return self.resIDs[indices_connected], char.strip(char.decode(self.segments[indices_connected]))


    # we only return atom indices for convenience because residues may share same names and resIDs
    def getInterconnectedResidues(self, return_indices=False, invert=False):
        '''Return atom indices of residues with interconnected bonds (residues connected to other residues)'''

        from ..base import constants
        from numpy  import arange, equal, isin, repeat, searchsorted, unique, vstack

        # expand the global resIDs array by repeating MAXCONN times
        global_resIDs = repeat(self.global_resIDs , constants.MAXCONN)
        indices       = repeat(arange(self.natoms), constants.MAXCONN)
        conn          = self.conn.ravel()

        # we don't need a transpose because NumPy is row-major
        abuff = vstack((global_resIDs, indices, conn))
        # remove unconnected rows
        abuff = abuff[:,abuff[2]>-1]

        # use the connected atoms' indices to index the global resIDs
        connected_resIDs = self.global_resIDs[abuff[2]]

        # only use the unmatched pairs
        masks = ~equal(abuff[0], connected_resIDs)

        # indices of interconnected atoms
        indices_interconnected = unique(abuff[1][masks])

        interconnected_residues = unique(self.residues[indices_interconnected])
        # return intraconnected residues
        if invert:
            if return_indices:
                return self.select(residues=interconnected_residues, invert=True)
            else:
                return unique(self.residues[self.select(residues=interconnected_residues, invert=True)])
        # return interconnected residues
        else:
            if return_indices:
                return self.select(residues=interconnected_residues)
            else:
                return interconnected_residues
                

    @property
    def global_resIDs(self):
        '''Global residue IDs of all segments'''

        from numpy import arange, array, searchsorted, unique, vstack

        # avoid recomputing
        try:
            return self._global_resIDs

        except:

            segments = unique(self.segments)
            resIDs   = unique(self.resIDs)
            classified_segments = searchsorted(segments, self.segments)
            classified_resIDs   = searchsorted(resIDs  , self.resIDs)

            # use a list of strings because searchsorted cannot handle 2D arrays (nor even array/list of tuples)
            lbuff = list(map(str, vstack((classified_segments, classified_resIDs)).T))

            unique_abuff = unique(lbuff)
            unique_lbuff = list(map(str, unique_abuff))

            self._global_resIDs = searchsorted(unique_lbuff, lbuff)

            return self._global_resIDs


    @property
    def _pdb_entries(self):
        '''PDB entries ATOM or HETATOM'''

        return self.__pdb_entreis


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

        self.__pdb_entreis = val


    @property
    def _old_order(self):
        '''Keep a reference to the original indices in case the sequence was reorderd'''

        from numpy import arange

        try:
            return self.__old_order
        except:
            self.__old_order = arange(self.natoms)
            return self.__old_order


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

        self.__old_order = val


    @property
    def _old2new_order(self):
        '''Mapping between the original and new indices in case the sequence was reorderd'''

        from numpy import arange

        try:
            return self.__old2new_order
        except:
            all_atoms = arange(self.natoms)
            return dict(zip(all_atoms,all_atoms))


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

        self.__old2new_order = val


    # YL 25/11/2021: this is a biomolecular version of _FragmentBase._covalent_radii
    #                which uses the PQR radii (saved in bfactors) whenever possible
    @property
    def _covalent_radii(self):
        '''Return the covalent radii'''

        try:
            return self.__covalent_radii
        except:
            # pre-defined PQR radii
            try:
                assert all(self.radii)
                self.__covalent_radii = self.radii
            except AssertionError:
                self.__covalent_radii = atomicdata.covalent_radii[self.znums]
            finally:
                return self.__covalent_radii


    def selectTarget(self, target, mute=False, strict=True):
        ''''''

        from numpy import array, append, char, unique

        kwargs = {}
        selected = array([], dtype=int)
        for k in [ 'atoms', 'indices', 'chainIDs', 'resIDs', 'segIDs', 'residues', 'segments' ]:
            kwargs.update({k:target})
            # YL 26/05/2024: added names_only (default: False) in FragmentBase.selectByAtoms() to prevent selecting indices as znums
            s = self.select(**{k:target}, strict=strict, names_only=True)
            selected = append(selected, s)
            if s.size == 1 and not mute:
                print(f'\n Target atom has been selected by keyword {k}={target}:\n {s[0]} ({char.decode(self.names[selected])[0]})')
                print( ' which belongs to')
            if s.size > 1 and not mute:
                print(f'\n Target atoms have been selected by keyword {k}={target}:\n {s} ({char.decode(self.names[selected])})')
                print(' which belong to')

        self.printSelected(selected)
        return selected


    # YL 02/05/2024
    def printSelected(self, selected):
        '''Print a summary of selected atoms's segments:residues/resIDs'''

        from numpy import array, char, isin, unique

        masks = isin(self._indices_all, selected)
        if not masks.sum():
            return
        # we do this for flattening the selection (indexing with indices is faster than with the full masks)
        indices = self._indices_all[masks]
        segments = char.decode(self.segments[indices])
        residues = char.strip(char.decode(self.residues[indices]))
        resIDs   = self.resIDs[indices]
        abuff = char.add(residues, array(resIDs).astype(str))
        abuff = char.add([':']*indices.size, abuff)
        abuff = unique(char.add(segments, abuff))
        for residue in abuff:
            print(f' {residue}')


    # YL 17/05/2024: replaced the previous angle and dihedral with two bond and two
    #                dihedral angles ikj, jkl, ijkl, and jklm to precisely control the binding position
    def bind(self, ligand, target=[], donor='', distance=3.0, ijk=180.0, jkl=90.0, ijkl=180.0, jklm=180.0, flip=True, unit='a.u.'):
        '''Bind a ligand/substrate to the target position
           Ligand      |      Acceptor
               I       |       L -- M
                \\     |     /
                   J --|-- K
              (donor)  | (target)
'''
    
        from numpy       import all, any, arange, argmax, array, append, char, isin, intersect1d, nan, searchsorted, setdiff1d, unique, vectorize
        from collections import OrderedDict
        from sys         import stdout
        from ..data      import ligands, unitconvert
        from ..io        import format
        from .           import xyz
    
        def _reorderFrag(fragObj, chosen):
            '''Return a list of new indices with chosen atom(s) leading followed by connected atoms and the remainder'''
    
            # atom(s) connected to the chosen atom
            conn = fragObj.conn[chosen]
            connected = conn[conn>-1]
            # make a new indexing list of [target, connected, remaining]
            leading = append(chosen, connected)
            remaining = setdiff1d(arange(fragObj.natoms), leading)
    
            order = append(leading, remaining)
            # remove duplicates whilst preserving the order
            order = list(OrderedDict.fromkeys(order))
            fragObj.reorder(order)
    
            return order
    
        print("\n Performing ChemShell ligand/substrate binding...")
        stdout.flush()
    
        # first loop up in the built-in ligands libraries
        try:
            ligandFrag = getattr(ligands, ligand.lower())
        # YL 03/05/2024: enabled user-customised ligand/substrate
        except:
            if issubclass(type(ligand), xyz._XYZ):
                ligandFrag = ligand
                if hasattr(ligand, 'optimiseResIDs'):
                    segments = unique(ligand.segments)
                    for segment in segments:
                        ligand.optimiseResIDs(segment)
            else:
                print(' >>> ERROR: Unsupported ligand', ligand, flush=True)
                exit(999)
    
        # remove all dummy atoms
        dummy = ligandFrag.select(atoms='X')
        ligandFrag.delete(indices=dummy, reconnect=True, mute=True)
        # reset Z-matrix (don't do this in self.delete()!)
        ligandFrag._master.zmatrix = 0
    
        # get the binding atom's index in the ligandFrag
        # YL NB 23/12/2021: bound is scalar!
        if type(donor) is str:
            bound = ligandFrag.select(atoms=donor, strict=False)
        elif type(donor) is int:
            bound = array([donor])
        else:
            print('\n ERROR: Illegal donor. Must be str (atom name) or int (atom index)')
            stdout.flush()
            exit(112)
    
        if bound.size == 1:
            bound = bound[0]
        # TODO multiple donor and target atoms
        elif bound.size > 1:
            format.printwrap("WARNING: More than one donor atoms chosen while the 1st is to be used as we do not yet support multiple donors/targets.",
                             indent='', initial_indent=' ', subsequent_indent=' '*10)
            bound = bound[0]
        elif bound.size == 0:
            format.printwrap("ERROR: No valid donor atom has been defined. Exiting...",
                             indent='', initial_indent=' ', subsequent_indent=' '*10)
            stdout.flush()
            exit(112)
    
        # reorder the Fragment to place the binding atom the first
        _reorderFrag(ligandFrag, bound)
    
        # we must do it with the help of dummy atoms
        if ligandFrag.natoms < 4:
            # this contains implict adding of dummy atoms
            ligandFrag.c2z()
        else:
            # get the Z-matrix first then insert dummy atoms
            ligandFrag.c2z(atom0=bound)
            # must place the atoms with redundant freedoms on the top to replace the dummy atoms later
            indices_masked = arange(ligandFrag.natoms)[ligandFrag.zmatrix['dihedral'].mask]
            _reorderFrag(ligandFrag, indices_masked)
            # insert dummy atoms if the first 3 atoms are not dummies already
            if not all(ligandFrag.znums[:3]==0):
                ligandFrag.insertDummy(recalculate=False)
                M = ligandFrag._master.zmatrix
                for item in [ 'j', 'k', 'l' ]:
                    # shift indices due to added dummy: excluding the 3 dummy atoms and 3 atoms containing redundant freedoms
                    M[item][6:] += 3

        # YL 28/01/2021: changed to target residue from target segment
        selected = self.selectTarget(target)
        selected = unique(selected)
        # YL 02/05/2024: we must confine the selection to a unique local residue
        indices_target = self.select(segments=self.segments[selected],
                                     residues=self.residues[selected],
                                     resIDs  =self.resIDs  [selected],
                                     chainIDs=self.chainIDs[selected])
        targetFrag = self.getSelected(indices_target)
    
        # target atom(s)
        targeted = searchsorted(indices_target, selected)
    
        # reorder the target fragment to make the targeted atom on the top
        order = _reorderFrag(targetFrag, targeted)
    
        # keep a reference (for printing only)
        old2new_order = dict(zip(order,targetFrag._indices_all))
        new2old_order = dict(zip(targetFrag._indices_all,order))
    
        # YL 18/05/2024: find a proper starting atom for c2z() to avoid using target atom as atom0 or atom1
        atom0s = targetFrag._indices_all[targetFrag.ncoordinations>1]
        for atom0 in atom0s:
            atom1 = -1
            # the first atom is the target which cannot be used as atom0 which doesn't possess j or k
            if not atom0:
                continue
            atom1s = targetFrag._indices_all[(targetFrag.ncoordinations>1)&any(targetFrag.conn==atom0,axis=1)]
            for atom1 in atom1s:
                # the first atom is the target which cannot be used as atom1 which doesn't possess j or k
                if not atom1:
                    continue
                else:
                    break
            if atom1 > 0:
                break
        # calculate the Z-matrix only when atom0 and atom1 can avoid atom 0
        if atom1 > 0:
            targetFrag.c2z(atom0=atom0, atom1=atom1)
        else:
            print(' >>> ERROR: Failed to construct internal coordinates using the the chosen target atom', flush=True)
            exit(112)

        # merge with the ligand fragment
        mergedFrag = targetFrag.merge(ligandFrag)
    
        # get indices of dummy and ligand atoms
        dummy = [ targetFrag.natoms, targetFrag.natoms+1, targetFrag.natoms+2 ]
        ndummies = 3
    
        # there should always be 3 dummy atoms
        ligand_atoms = setdiff1d(mergedFrag.stacks['FragmentBio1'], dummy)
    
        # index of donor atom in the merged fragment excluding dummies
        bound = targetFrag.natoms

        Z = mergedFrag._master['zmatrix']

        # connect donor and target
        for item in [ 'j', 'k', 'l' ]:
            for i in ligand_atoms:
                entry = Z[item][i]
                # skip dummy atoms
                if entry in dummy:
                    Z[item][i] -= targetFrag.natoms
                # ligand atoms
                if entry in ligand_atoms:
                    Z[item][i] -= ndummies
    
        # delete dummy atoms
        mergedFrag.delete(indices=dummy, mute=True)

        Z = mergedFrag._master['zmatrix']
    
        # adjust distance (J-K)
        Z['bond'][bound] = distance*unitconvert.unit2au[unit.strip().lower()]

        # Yl 18/05/2024: determine the critical indices according to the Z-matrix dependencies -- K is the target atom
        J = bound
        K = Z['j'][J]
        L = Z['j'][K]
        M = Z['k'][K]
        Z['k'][J] = L
        Z['l'][J] = M

        I = mergedFrag._indices_all[Z['j']==J][0]
        Z['k'][I] = Z['j'][J]
        Z['l'][I] = Z['k'][J]
        
        # YL 03/05/2024: update the angle (J-K-L) and dihedral (J-K-L-M) inherited from dummies (which were 90.0 and -90.0 degrees)
        # adjust angle (I-J-K)
        Z['angle']   [I] = ijk
        Z['angle']   [J] = jkl
        # adjust dihedral (I-J-K-L): NB, we always use the 2nd atom but there is no certain definition which dihedral this will be
        Z['dihedral'][I] = ijkl
        Z['dihedral'][J] = jklm

        # fix the first three atoms and calculate the whole coordinates using Z-matrix
        i = argmax((Z[:].j==-1)&(Z[:].k==-1)&(Z[:].l==-1))
        j = argmax((Z[:].j!=-1)&(Z[:].k==-1)&(Z[:].l==-1))
        k = argmax((Z[:].j!=-1)&(Z[:].k!=-1)&(Z[:].l==-1))
        
        mergedFrag.z2c(coords_atom0=mergedFrag.coords[i],
                       coords_atom1=mergedFrag.coords[j],
                       coords_atom2=mergedFrag.coords[k])
    
        # test the closest atom in the original fragment
        closests = []
        for index in ligand_atoms-ndummies:
            closest = self.getClosest(mergedFrag.coords[index], return_dist=True)[1]
            closests += [closest]
        closest = min(closests)
    
        ligandFrag = mergedFrag.getSelected(ligand_atoms-ndummies)
    
        # for later printing only
        # restore indexing information from self.select() and targetFrag.reorder()
        # because we don't calculate the Z-Matrix of boundFrag and we won't have a
        # chance to get K/L/M directly from Z-Matrix later
        dI = I - J # relative offset
        K = indices_target[new2old_order[Z['j'][J]]]
        L = indices_target[new2old_order[Z['k'][J]]]
        M = indices_target[new2old_order[Z['l'][J]]]
    
        # flip
        if flip:
            Z['angle'][J] *= -1.0
            i = argmax((Z[:].j==-1)&(Z[:].k==-1)&(Z[:].l==-1))
            j = argmax((Z[:].j!=-1)&(Z[:].k==-1)&(Z[:].l==-1))
            k = argmax((Z[:].j!=-1)&(Z[:].k!=-1)&(Z[:].l==-1))
            mergedFrag.z2c(coords_atom0=mergedFrag.coords[i],
                           coords_atom1=mergedFrag.coords[j],
                           coords_atom2=mergedFrag.coords[k])
        
            # test the closest atom again in the original fragment
            flipped_closests = []
            for index in ligand_atoms-ndummies:
                flipped_closest = self.getClosest(mergedFrag.coords[index], return_dist=True)[1]
                flipped_closests += [flipped_closest]
            flipped_closest = min(flipped_closests)
        
            # TODO: rotate about the dihedral to avoid close atoms
        
            # choose the coordinates having no close atoms
            if flipped_closest > closest:
                print(' Flipping to avoid clashing...')
                ligandFrag = mergedFrag.getSelected(ligand_atoms-ndummies)

        boundFrag = self.merge(ligandFrag, reconnect=True)
    
        # YL 28/01/2021: residues connected via bonds must belong to one DL_FIELD/DL_POLY segment ("group")
        boundFrag.connect(force_reconnect=True, mute=True)
        natoms_ligand = ligandFrag.natoms
        fragID_ligand = boundFrag.fragIDs[-natoms_ligand:]
        connected_segment = unique(boundFrag.segments[isin(boundFrag.fragIDs, fragID_ligand)])
        target_segment    = unique(targetFrag.segments)[0]
        # if connected to the target residue via bond(s)
        if target_segment in connected_segment:
            boundFrag.segments[-natoms_ligand:] = target_segment
            # shift the resIDs to avoid conflict
            boundFrag.resIDs[-natoms_ligand:] += boundFrag.resIDs[boundFrag.select(segments=target_segment)].max()
    
            # also make sure the segment's indices are contiguous
            boundFrag.keepSegmentContiguous(target_segment)
    
        if type(ligand) is str:
            name_ligand = ligand
        else:
            # PDB
            try:
                name_ligand = ''.join(unique(char.decode(ligandFrag.segments)))
            # XYZ
            except:
                name_ligand = ligandFrag.formula
    
        # assign a title
        boundFrag.title = self.title.strip()+'+'+str(name_ligand)
    
        # print a summary of the geometrical positon of binding
        offset = len(boundFrag.stacks['FragmentBio0'])
        # for printing only: get the new indices after keepSegmentContiguous() of the ligand's first two atoms
        J = boundFrag._old2new_order[offset]
        I = boundFrag._old2new_order[offset+dI]
        K = boundFrag._old2new_order[K]
        L = boundFrag._old2new_order[L]
        M = boundFrag._old2new_order[M]
        names = char.decode(boundFrag.names)
        name_I, name_J, name_K, name_L, name_M = names[I], names[J], names[K], names[L], names[M]
        jk   = boundFrag.dist(J,K)
        ijk  = boundFrag.angle(I,J,K)
        jkl  = boundFrag.angle(J,K,L)
        ijkl = boundFrag.dihedral(I,J,K,L)
        jklm = boundFrag.dihedral(J,K,L,M)
        print(f'\n {name_ligand} bound to the current fragment in chosen position:')
        print(' Distance {}({})-{}({}): {:.3f} a.u.'.format(J, name_J, K, name_K, jk))
        print(' Angle {}({})-{}({})-{}({}): {:.2f} degrees'.format(I, name_I, J, name_J, K, name_K, ijk))
        print(' Angle {}({})-{}({})-{}({}): {:.2f} degrees'.format(J, name_J, K, name_J, L, name_K, jkl))
        print(' Dihedral {}({})-{}({})-{}({})-{}({}): {:.2f} degrees'.format(I, name_I, J, name_J, K, name_K, L, name_L, ijkl))
        print(' Dihedral {}({})-{}({})-{}({})-{}({}): {:.2f} degrees'.format(J, name_J, K, name_K, L, name_L, M, name_M, jklm))
        stdout.flush()
    
        # rename: after keepSegmentContiguous() boundFrag.stacks can no longer be used!
        indices_new = [ boundFrag._old2new_order[i] for i in boundFrag.stacks['FragmentBio1'] ]
        # YL 03/05/2024: PDB ligand (substrate)
        try:
            resname = ''.join(unique(char.decode(ligandFrag.residues)))[:5]
            segname = ''.join(unique(char.decode(ligandFrag.segments)))[:4]
            resIDs  = ligandFrag.resIDs
        # XYZ ligand
        except:
            resname = name_ligand[:5]
            segname = 'LIG '
            resIDs  = 1
        boundFrag.residues[indices_new] = resname
        boundFrag.segments[indices_new] = segname
        boundFrag.resIDs  [indices_new] = resIDs

        # save to a new file
        boundFrag.save('_bound_{}.pdb'.format(name_ligand))

        return boundFrag


    def retainChainIDs(self, refFrag):
        '''Take the chainIDs from the reference fragment which may contain a different number of atoms'''

        from numpy import char, unique

        for segment in refFrag.all_segments:
            masks_segment = char.lower(char.strip(self.segments)) == segment.strip().lower()
            masks_ref = refFrag.segments == segment
            all_resIDs = unique(refFrag.resIDs[masks_ref])
            for resID in all_resIDs:
                masks_resID = self.resIDs == resID
                # remove blank chainID
                chainID = [ c for c in unique(refFrag.chainIDs[(refFrag.resIDs==resID)*masks_ref]) if c.strip() ]
                if len(chainID) <= 1:
                    chainID = chainID[0]
                    self.chainIDs[masks_resID*masks_segment] = chainID
                else:
                    print(f'\n WARNING: Multiple chainIDs found in segment {segment} resID {resID}')


    def keepSegmentContiguous(self, segname):
        '''Make a given segment contiguous in the atoms sequence if it is not'''

        from numpy import arange, diff, isin, hstack
        from sys   import stdout
        from ..io  import format

        try:
            segname = segname.decode()
        except:
            pass

        selected_segment = self.select(segments=segname, strict=True)
        if not all(diff(selected_segment)==1):
            print('\n\n ' + '*'*70)
            format.printwrap('NB: The atom sequence has been reordered to keep segment "{}" contiguous in the sequence'.format(segname))
            print('\n '+'*'*70+'\n')
            stdout.flush()

            all_atoms  = arange(self.natoms)
            rest_atoms = all_atoms[isin(all_atoms, selected_segment, invert=True)]
            new_order  = hstack((selected_segment, rest_atoms))
            old_order  = all_atoms[new_order]
            # embed a reference for later processing
            self._old_order = old_order
            self._old2new_order = dict(zip(old_order,all_atoms))

            # similar to FragmentBio.delete(), some entries have to be deleted or rearranged
            self._master[:] = self._master[new_order]
            # `.pdb_entries` only exist in fragment read in from .PDB files
            try:
                self._pdb_entries = self._pdb_entries[new_order]
            except:
                pass
            # reconnect
            self.connect(force_reconnect=True, mute=True)
            # if there is a `.psf` object
            try:
                del self._psf
            except:
                pass

            if chemsh.CHEMSH_DEBUG > 4:
                print('\n new index  old index  name  res.  resID  x (a.u.)  y (a.u.)  z (a.u.)  charge\n')
                strbuff = self.coords2str(inserts=[(0,''),
                                                   (1,all_atoms,'%7d'),
                                                   (2,''),
                                                   (3,old_order,'%7d'),
                                                   (5,''),
                                                   (6,self.residues),
                                                   (7,self.resIDs,'%4d')],
                                          coords=3,
                                          charges=2)
                print(strbuff)

            self.save('_rearranged.pdb')


    def optimiseResIDs(self, segname):
        '''Optimise the resIDs of the specified segment in place by reording them from 1 and filling gaps'''

        from numpy import all, any, append, arange, argwhere, char, concatenate, diff, isin, mod, put, searchsorted, split, unique
        from sys   import stdout
        from ..io  import format

        # YL 18/02/2023: corrected this comparison with isin() (NumPy >= 1.13)
#        masks_segment = self.segments==segname
        masks_segment = isin(char.upper(char.strip(char.decode(self.segments))), segname.strip().upper())
        selected_segment = arange(self.natoms)[masks_segment]

        resIDs_unique = unique(self.resIDs[selected_segment])
        chainIDs_unique = unique(self.chainIDs[selected_segment])

        # YL 05/08/2021: check if there are duplicate resIDs within a segment's unique chain
        # YL 19/08/2021: the performance can of course be further improved by pre-chopping the atoms
        #                ranges according to resIDs but it's currently fast enough already
        for chainID_unique in chainIDs_unique:
            counter = 0
            strbuff = ''
            masks_segment_chain = masks_segment*(self.chainIDs==chainID_unique)
            resIDs_unique_in_chain = unique(self.resIDs[masks_segment_chain])
            for resID_unique in resIDs_unique_in_chain:
                # YL: to explan, for example, the following case is illegal within chain "A"
                #     atom        0    1    2    3    4  ...   84   85   86   87  ...  152  153  154  155
                #     resisues  ALA  ALA  ALA  ALA  SER  ...  ALA  ALA  ALA  ALA  ...  ALA  ALA  ALA  ALA
                #     resIDs      1    1    1    1    2  ...    1    1    1    1  ...    1    1    1    1
                #     chainIDs    A    A    A    A    A  ...    A    A    A    A  ...    B    B    B    B
                # YL 27/11/2021: we can directly use self.resIDs==resID_unique which is way faster than the select() function
#                masks_residue = self.select(resIDs=resID_unique, strict=True, return_masks=True)
                masks_residue = self.resIDs==resID_unique
                masks_residue *= masks_segment_chain
                selected_residue = arange(self.natoms)[masks_residue]
                if any(diff(selected_residue)>1):
                    nodes = arange(selected_residue.size)[append([False], diff(selected_residue)>1)]
                    residues_split = split(selected_residue, nodes)
                    for residue_split in residues_split[1:]:
                        resID_max = self.resIDs[selected_segment].max()
                        # shift the number
                        strbuff += f' atoms: {residue_split[0]} to {residue_split[-1]} resID: {unique(self.resIDs[residue_split])[0]}'
                        self.resIDs[residue_split] = resID_max + 1
                        strbuff += f' -> {resID_max + 1}\n'
                        counter += 1
            if counter:
                chainID = chainID_unique.decode().strip()
                format.printwrap(f'NB: {counter} duplicate resIDs reassigned for segment \"{segname.strip()}\" chain \"{chainID}\"')
                if chemsh.CHEMSH_DEBUG:
                    print(strbuff)

        # for example, original resIDs [5,5,5,5,3,3,3,6,6,6,6,6,6,9,9,9,9,9,1,1,1,1,1]
        # will become [3,3,3,3,2,2,2,4,4,4,4,4,4,5,5,5,5,5,1,1,1,1,1]
        resIDs_unique = unique(self.resIDs[selected_segment])
        self.resIDs[selected_segment] = searchsorted(resIDs_unique, self.resIDs[selected_segment]) + 1

        # YL 07/08/2021: there can be at most 9999 resIDs in a segment for the time being
        segments_grouped = (self.resIDs[selected_segment]/9999).astype(int)
        masks_9999 = segments_grouped > 0
        if any(masks_9999):
            # renumber resIDs exceeding 9999
            # self.resIDs[selected_segment][masks_9999] = ... doesn't change the values
            put(self.resIDs, selected_segment[masks_9999], mod(self.resIDs[selected_segment][masks_9999], 9999))
            # rename segment names of resIDs exceeding 9999 by appending a number
            segments_grouped_unique = unique(segments_grouped)
            for segment_grouped_unique in segments_grouped_unique[1:]:
                length_extra = len(str(segment_grouped_unique))
                masks_segment_to_rename = segments_grouped==segment_grouped_unique
                format_with_renumber = f'%-s{segment_grouped_unique}'
                # indexing using self.segments[selected_segment][masks_segment_to_rename] cannot set the values
                put(self.segments, selected_segment[masks_segment_to_rename], format_with_renumber%(segname.strip()[:4-length_extra]))


        # YL 26/10/2021: let [3,3,3,3,2,2,2,4,4,4,4,4,4,5,5,5,5,5,1,1,1,1,1] become
        #                [1,1,1,1,2,2,2,3,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5]
        # find nodes: also attach first (0) and last indices
        dbuff = diff(self.resIDs[selected_segment])
        # only proceed when there are discrete numbers
        if not all((dbuff==0)+(dbuff==1)):
            nodes = concatenate(([0], argwhere(dbuff!=0).ravel()+1, [selected_segment.size]))
            for i, node in enumerate(nodes[:-1]):
                put(self.resIDs, selected_segment[nodes[i]:nodes[i+1]], i+1)

# YL TODO 27/02/2021: also shift the resID of terminal
#        self.terminal.resID


    # YL 28/01/2021: residues connected via bonds must belong to one DL_FIELD/DL_POLY segment ("group")
    def mergeConnectedResiduesToSegment(self, segname, ff='charmm', ions_in_solvent=False):
        '''Merge residues connected via bonds to the selected segment'''

        from numpy  import any, argmin, char, isin, unique
        from ..data import forcefield

        ffObj = forcefield.getFFObj(ff)

        # in case of bytes
        try:
            segname = segname.decode()
        except:
            pass

        # realign resIDs to start from 1 and fill in any gap
        self.optimiseResIDs(segname)

        try:
            masks_segment = char.decode(char.lower(char.strip(self.segments))) == segname.strip().lower().decode()
        except:
            masks_segment = char.decode(char.lower(char.strip(self.segments))) == segname.strip().lower()
        indices_segment = self._indices_all[masks_segment]

        # YL 24/02/2021: residues of a same name may have some connected while some other not
        masks_connected = isin(self.fragIDs, self.fragIDs[masks_segment])*(~masks_segment)

        # YL 26/11/2023: we also need to exclude non-standard amino acids
        masks_other_amino_acids = self.selectByNonProteinogenic(ff=ff, return_masks=True, mute=True)
        residues_other_amino_acids = char.decode(char.strip(unique(self.residues[masks_other_amino_acids*masks_segment])))

        # YL 18/02/2021: only perform the operation when there are residues to connect
        if masks_connected.sum():

            # final list of indices to rename segments, used later
            indices_residues_to_reassign = []

            residues_connected_to_segment = char.strip(char.decode(self.residues[masks_connected]))
            # exclude amino acids, water, and ions
            if ions_in_solvent:
                residues_connected_to_segment = set(residues_connected_to_segment) - set(ffObj.amino_acids) - set(residues_other_amino_acids) - set(ffObj.waters) - set(ffObj.ions)
            # exclude amino acids and water
            else:
                residues_connected_to_segment = set(residues_connected_to_segment) - set(ffObj.amino_acids) - set(residues_other_amino_acids) - set(ffObj.waters)

            masks_residues_connected_to_segment = self.select(residues=list(residues_connected_to_segment), return_masks=True, strict=True)
            masks_residues_connected_to_segment *= masks_connected

            # YL 16/03/2024: it's possible that segments are connected by some residues and have same
            #                fragIDs and such assignment is messed up by fragIDs; now we only consider
            #                direct connection between a residue and the current segment
            # a temporary fragment object of the current segment only
            segFrag = self.getSelected(indices_segment)
            segFrag.save('segFrag.pdb', mute=True)
            lbuff = []
            list_residues_excluded = []
            indices_residues_connected_to_segment = self._indices_all[masks_residues_connected_to_segment]
            # loop over all atoms connected
            for i in indices_residues_connected_to_segment:
                resname = f'{self.segments[i].decode().strip()}:{self.residues[i].decode().strip()}{self.resIDs[i]}'
                # skip if the current residue is already checked or not connected
                if resname in lbuff or resname in list_residues_excluded:
                    continue
                indices_current_residue = self._indices_all[(self.segments==self.segments[i])*(self.resIDs==self.resIDs[i])]
                # a temporary fragment object formed by merging the current residue and segment
                tmpFrag = segFrag.merge(self.getSelected(indices_current_residue))
                tmpFrag.save(f'{resname}.pdb', mute=True)
                # if they share a same fragID
                if any(isin(tmpFrag.fragIDs[tmpFrag.stacks['FragmentBio1']], tmpFrag.fragIDs[tmpFrag.stacks['FragmentBio0']])):
                    lbuff += [ resname ]
                    indices_residues_to_reassign += indices_current_residue.tolist()
                # for skipping
                else:
                    list_residues_excluded += [ resname ]

            # YL 15/03/2024: print with resIDs
            if len(lbuff):
                plural = {1:''}.get(len(lbuff), 's')
                print(f' NB: Residue{plural} (seg:resID) connected to "{segname.strip()}" also included in the segment:')
                print( '     {}'.format(' '.join(unique(lbuff)))+'\n')

            # YL 28/01/2021: shift the resIDs to avoid clash
            # YL 24/02/2021: refactored to avoid exceeding 9999
            # YL 16/03/2024: using indices instead of masks is faster
            resIDs_unique   = unique(self.resIDs  [indices_residues_to_reassign])
            residues_unique = unique(self.residues[indices_residues_to_reassign])
            segments_unique = unique(self.segments[indices_residues_to_reassign])
            counter = 0
            # max residue ID
            max_resID_segment = self.resIDs[masks_segment].max()
            for resID in resIDs_unique:
                for residue in residues_unique:
                    for segment in segments_unique:
                        selected_residue = self.select(resIDs=resID, residues=residue, segments=segment, strict=True)
                        if len(selected_residue):
                            if max_resID_segment + counter <= 9999:
                                counter += 1
                                self.resIDs[selected_residue] = max_resID_segment + counter
                            else:
                                print('\n WARNING: we do not support resIDs greater than 9999!\n')

            # we only rename the segment after the resIDs have been shifted
            self.segments[indices_residues_to_reassign] = segname
            # YL 26/11/2023: we must reset self.all_segments after any change to self.segments
            try:
                del self._all_segments
            except AttributeError:
                pass


    def neutralise(self, cation='Na', anion='Cl', ion_segment='', solvent='SOL', ff='charmm'):
        ''''''

        from numpy  import arange, array2string, random, unique
        from sys    import stdout
        from ..data import forcefield

        ffObj = forcefield.getFFObj(ff)

        # fetch some data
        ion_charges = forcefield.ion_charges
        ion_znums = forcefield.ion_znums
        ion_names = ffObj.ion_symbols2names

        ncations = 0
        nanions  = 0

        charge  = round(self.net_charge)
        icharge = int(round(charge))
        print("\n Net total charge: {}\n".format(charge))

        if icharge > 0:
            nanions = -int(icharge/ion_charges.get(anion.title(), -1))
            ion = anion.title()
            default_name, default_znum, default_charge = 'CLA', 17, -1.0
        elif icharge < 0:
            if icharge%ion_charges.get(cation.title(), 1):
                print("\n The default {} is used instead, as {} is not suitable for the system's total charge".format('Na', cation))
            ncations = -int(icharge/ion_charges.get(cation.title(), 1))
            ion = cation.title()
            default_name, default_znum, default_charge = 'SOD', 11, 1.0
        else:
            print("\n ChemShell neutralisation: The system is already neutral\n")
            return self

        nions = ncations or nanions

        # replace the water oxygens
        solvent_o = self.select(atoms='O', residues='TIP', segments=solvent, union=False, strict=False)
        to_replace = random.choice(solvent_o, nions, replace=False)

        # change the residues first to exclude the oxygens
        self.residues[to_replace] = ion_names.get(ion, default_name)

        # we have to loop over the atoms to replace because the connected atoms may have same resIDs/residues/segments
        to_delete = []
        for i in to_replace:
            resID = self.resIDs[i]
            segment = self.segments[i]
            to_delete += self.selectBio(resIDs=resID, residues='TIP', segments=segment, union=False, strict=False).tolist()

        selectcases = {  1:'+',
                         2:'2+',
                        -1:'-',
                      }
        if nions == 1:
            print(" {} {}({}) ion has been added to neutralise the system. Coordinates in a.u.:\n".format(nions,
                                                                                                          ion,
                                                                                                          selectcases[ion_charges[ion]]))
        elif nions > 1:
            print(" {} {}({}) ions have been added to neutralise the system. Coordinates in a.u.:\n".format(nions,
                                                                                                            ion,
                                                                                                            selectcases[ion_charges[ion]]))
        print(' '+array2string(self.coords[to_replace])[1:-1])
        stdout.flush()

        # replace the entries
        # YL TODO 27/02/2021: we now only support single-atom ions
        self.names[to_replace] = ion_names.get(ion, default_name)
        self.types[to_replace] = ion_names.get(ion, default_name)
        self.znums[to_replace] = ion_znums.get(ion, default_znum)
        self.charges [to_replace] = ion_charges.get(ion, default_charge)
        # YL 11/09/2021: we keep the segment name by default
        if ion_segment:
            self.segments[to_replace] = ion_segment
            # YL 26/11/2023: we must reset self.all_segments after any change to self.segments
            try:
                del self._all_segments
            except AttributeError:
                pass

        # delete the PSF entries of the oxygens that are replaced
        # YL 09/09/2020: we don't need this for now because all entries involving
        #                oxygen also involve hydrogen and will be removed anyway below
#        self._psf.delete(to_replace)

        # delete the unwanted at last (the associated PSF entries are also deleted within the delete() function)
        self.delete(indices=to_delete, reconnect=False, mute=True)

        self._psf.titles += [' REMARKS segment {}'.format(ion_segment) + ' { first NONE; last NONE; auto none }\n']


    def load_psf(self, filename):
        '''Load information from the provided PSF filename to the current fragment object'''

        from ..io import file

        file.parse(filename, 'psf', receiver=self)


    @property
    def _psf(self):
        '''An object to store PSF data'''

        from sys import modules

        try:
            return self.__psf
        except:
            self.__psf = getattr(modules[__name__], '__PSF')()
            return self.__psf


    def analyse(self, dcd='_namd.dcd', xst='_namd.xst', nsnapshots=10, **kwargs):
        '''Load the DCD file and do analysis according to the options provided'''

        from ..tools import analysis

        if dcd:
            analysis.analyseMD(self, dcd=dcd, nsnapshots=nsnapshots, **kwargs)


    # YL 01/07/2023: invented this new algorithm
    # YL 31/10/2023: changed the algorithm; we probably should publish this
#    @miscutils.time
    def _getPBCWrapCentre(self, accuracy=1.0E-2, select=[], check_duplicate=False):
        '''Get a sensible centre for PBC wrapping when the molecule to wrap spans two faces of the bounding box in different directions'''

        from itertools import product
        from numpy     import all, any,  arange, argmax, argmin, append, array, dot, isin, linalg
        from numpy     import mean, subtract, tile, unique, vstack, zeros
        from ..base    import constants, lattice

        if len(select):
            indices = self._indices_all[select]
            unselected = self._indices_all[~isin(self._indices_all,select)]
        else:
            indices = ...
            unselected = self._indices_all

        # create a template
        images = self.__class__(natoms=0, connect_mode=None)

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

        ndims = selectcases[self._cell.ndimensions]

        # no translation will be made

        rangel = 0
        ranger = 2

        # *unselect* the chosen atoms (solute) to get the centre of the solvent box!
        centre_box = self._getBoundingBoxCentre(select=unselected)
        vectors = self._cell.vectors
        volume  = self._cell.volume
        A, B, C = self._cell.A, self._cell.B, self._cell.C

        icell = 0
        natoms_0 = len(indices)
        for translation in product(range(rangel,ranger), repeat=ndims):

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

            # YL TODO 01/11/2023: or this might not be true? let me switch it off for now
            # there could be duplicates when applying periodicity to the unitcell
            coords = self.coords[indices]+displacement
            centre_unitcell = centre_box + displacement
            if icell > 0:
                centres_unitcell = vstack((centres_unitcell, centre_unitcell))
            else:
                centres_unitcell = centre_unitcell
            c = coords[:,None,:]
            if icell > 0 and check_duplicate:
                c_diff = subtract(c, images.coords)
                masks_valid = all(linalg.norm(c_diff, axis=2) > accuracy, axis=1)
                indices_valid = indices[masks_valid]
            else:
                indices_valid = indices
# YL: the above method isn't faster than this explicit generator--
#            indices_valid = array([ i for i, c in enumerate(coords) if all(linalg.norm(c-images.coords, axis=1) > accuracy) ])

            # create a fragment for the unit cell to insert
            unitCell = self.__class__(natoms=len(indices_valid), connect_mode=None)
            unitCell.coords  = coords[indices_valid]
            unitCell.names   = self.names[indices][indices_valid]
            unitCell.znums   = self.znums[indices][indices_valid]
            unitCell.charges = self.charges[indices][indices_valid]
            unitCell.residues = self.residues[indices][indices_valid]
            unitCell.segments = self.segments[indices][indices_valid]

            # insert unit cell
            images.insert(unitCell, 0)
            icell += 1

        # find atoms in each bounding unit cell and label them
        natoms = images.natoms
        # initialise the region label with -1 (we borrow the `frozen` column)
        images.frozen = -1
        for icell, centre_unitcell in enumerate(centres_unitcell):

            # translate for using half of the volume below
            images.translate(vector=-centre_unitcell)
            c = images._coords_transpose
            # vstack acts on the 1st axis
            stacked_AB = vstack((c, tile(vectors[:,[0,1]], (natoms,1)).reshape(natoms,3,2).T)).T
            stacked_CA = vstack((c, tile(vectors[:,[2,0]], (natoms,1)).reshape(natoms,3,2).T)).T
            stacked_BC = vstack((c, tile(vectors[:,[1,2]], (natoms,1)).reshape(natoms,3,2).T)).T
            # find atoms beyond the unitcell's box
            overflow_AB_p = linalg.det(stacked_AB)/volume >  0.5
            overflow_CA_p = linalg.det(stacked_CA)/volume >  0.5
            overflow_BC_p = linalg.det(stacked_BC)/volume >  0.5
            overflow_AB_n = linalg.det(stacked_AB)/volume < -0.5
            overflow_CA_n = linalg.det(stacked_CA)/volume < -0.5
            overflow_BC_n = linalg.det(stacked_BC)/volume < -0.5
            masks_inside = ~(overflow_AB_p+overflow_CA_p+overflow_BC_p+overflow_AB_n+overflow_CA_n+overflow_BC_n)
            # give the inside atoms a region label
            images.frozen[masks_inside] = icell

            # translate back
            images.translate(vector=centre_unitcell)

            icell += 1

        # do not connect until now
        images.connect_mode = self.connect_mode
        images.conn
        region_bonds = images.frozen[images.bonds]
        # find bonds A-B where A and B belong to different bounding unit cell
        masks_boundary = region_bonds[:,0] != region_bonds[:,1]
        # find bonds A-B where either A or B belongs to bounding unit cell 0
        masks_boundary *= (region_bonds[:,0]==0)+(region_bonds[:,1]==0)
        # find bonds A-B where neither A nor B is outside the bounding box of the whole images
        masks_boundary *= (region_bonds[:,0]!=-1)*(region_bonds[:,1]!=-1)

        indices_boundary_bonds = unique(images.bonds[masks_boundary])
        # we use the centre of mass of the selected atoms if there is no bonds crossing the boundary
        if not len(indices_boundary_bonds):
            centroid_mass_weighed = (self.coords[indices]*self.masses[indices][:,None]).sum(axis=0)/self.masses[indices].sum()
        # use the mass-weighed gravity for PBC wrapping
        else:
            fragIDs = unique(images.fragIDs[indices_boundary_bonds])

#            nfrags = len(fragIDs)
#            # if only one fragment sticks out of the box, we use its whole centroid (to avoid the other end sticks out after wrapped)
#            if nfrags == 1:
#                masks_frag = images.fragIDs == fragIDs[0]
#                centroid_mass_weighed = (images.coords[masks_frag]*images.masses[masks_frag][:,None]).sum(axis=0)/images.masses[masks_frag].sum()
#            # otherwise we use the intersection's centroid
#            else:

            def _setBoundingVectors(_frag, _A_ref, _B_ref, _C_ref):
                '''Set the minimal bounding box's lattice vectors according to given reference vectors (_A/B/C_ref)'''

                # project all atoms onto the original cell vector _A_ref
                on_A = dot(_frag.coords,_A_ref[:,None])/dot(_A_ref,_A_ref)*_A_ref
                # YL TODO: check the correctness!
                # calculate new vectors using the far ends of the points projected on the original vector _A_ref
                A_current_frag = on_A[:,0].max()-on_A[:,0].min(), \
                                 on_A[:,1].max()-on_A[:,1].min(), \
                                 on_A[:,2].max()-on_A[:,2].min()
                # project all atoms onto the original cell vector _B_ref
                on_B = dot(_frag.coords,_B_ref[:,None])/dot(_B_ref,_B_ref)*_B_ref
                # calculate new vectors using the far ends of the points projected on the original vector _B_ref
                B_current_frag = on_B[:,0].max()-on_B[:,0].min(), \
                                 on_B[:,1].max()-on_B[:,1].min(), \
                                 on_B[:,2].max()-on_B[:,2].min()
                # project all atoms onto the original cell vector _C_ref
                on_C = dot(_frag.coords,_C_ref[:,None])/dot(_C_ref,_C_ref)*_C_ref
                # calculate new vectors using the far ends of the points projected on the original vector _C_ref
                C_current_frag = on_C[:,0].max()-on_C[:,0].min(), \
                                 on_C[:,1].max()-on_C[:,1].min(), \
                                 on_C[:,2].max()-on_C[:,2].min()
                # assign values
                _frag.cell.vectors = array([A_current_frag,B_current_frag,C_current_frag]).T

            def _InBoundingBox(_frag, _frag_ref, _fragIDs):
                '''Return a mask array of atoms of _frag lying within the minimal bounding box of _frag_ref (obtained from _setBoundingVectors())'''

                _natoms = _frag.natoms

                # natoms of _frag_ref cannot be 1
                if _frag_ref.natoms == 1:
                    # only the current fragment's single atom is included (the list _fragIDs contains only one element)
                    return images.fragIDs == _fragIDs[0]

                # translate for using half of the volume below
                centre_frag_ref = _frag_ref._getBoundingBoxCentre()
                _frag.translate(vector=-centre_frag_ref)

                # shortcuts
                c = _frag._coords_transpose
                V = _frag_ref.cell.vectors
                v = _frag_ref.cell.volume
                # vstack acts on the 1st axis
                stacked_AB = vstack((c, tile(V[:,[0,1]], (_natoms,1)).reshape(_natoms,3,2).T)).T
                stacked_CA = vstack((c, tile(V[:,[2,0]], (_natoms,1)).reshape(_natoms,3,2).T)).T
                stacked_BC = vstack((c, tile(V[:,[1,2]], (_natoms,1)).reshape(_natoms,3,2).T)).T
                # find atoms beyond the unitcell's box
                # YL NB: an additional distance must be added on to eliminate the incorrect exclusion caused by numerical errors;
                # YL TODO: or shall we use natoms_0 +- 1% to check the number of atoms below?
                overflow_AB_p = linalg.det(stacked_AB)/v >  0.5+constants.MEDIUM
                overflow_CA_p = linalg.det(stacked_CA)/v >  0.5+constants.MEDIUM
                overflow_BC_p = linalg.det(stacked_BC)/v >  0.5+constants.MEDIUM
                overflow_AB_n = linalg.det(stacked_AB)/v < -0.5-constants.MEDIUM
                overflow_CA_n = linalg.det(stacked_CA)/v < -0.5-constants.MEDIUM
                overflow_BC_n = linalg.det(stacked_BC)/v < -0.5-constants.MEDIUM

                # translate back
                images.translate(vector=centre_frag_ref)

                # mask array of atoms that lie within this bounding box
                return ~(overflow_AB_p+overflow_CA_p+overflow_BC_p+overflow_AB_n+overflow_CA_n+overflow_BC_n)


            def _getIndicesAllRelated(_frag, _A_ref, _B_ref, _C_ref, _fragIDs):
                '''Recursively return all related atoms lying within any bounding box of fragments, by expanding the list of _fragIDs'''

                # initialise the result mask array
                _masks_related = zeros((_frag.natoms,), dtype=bool) 

                # we need to create a fragment object given a list of _fragIDs
                _currentFrag = _frag.getSelectedFragIDs(_fragIDs)
                # define the minimal bounding box (lattice vectors) of which the vectors share the same directions of the original lattice
                _setBoundingVectors(_currentFrag, _A_ref, _B_ref, _C_ref)
                # get a mask array of all atoms lying within the mininal bounding box of _currentFrag
                _masks_in_bounding_box = _InBoundingBox(_frag, _currentFrag, _fragIDs)
                # expand _masks_related by including all fragIDs
                _masks_related += _masks_in_bounding_box
                # get the frag IDs contained in this bounding box
                _fragIDs_in_bounding_box = unique(_frag.fragIDs[_masks_in_bounding_box])

                # we return the result when the cluster no longer grows (namely, no new fragID can be included)
                # _fragIDs_in_bounding_box could be empty (especially when the fragment is a single atom)
                if set(_fragIDs_in_bounding_box) == set(_fragIDs) or not len(_fragIDs_in_bounding_box):
                    return _masks_related

                # update the list of _fragIDs
                _fragIDs = unique(_fragIDs+_fragIDs_in_bounding_box.tolist()).tolist()
                # recursively run the function
                _masks_related += _getIndicesAllRelated(_frag, _A_ref, _B_ref, _C_ref, _fragIDs)

                return _masks_related
                

            # loop over the boundary fragments (NB: we cannot loop over all fragments)
            for fragID in fragIDs:
                # use the recursive function defined above to get all atoms of the cluster in which the current fragment is
                masks_all_related = _getIndicesAllRelated(images, A, B, C, [fragID])
                # we can return the result as soon as we find a cluster that contains all atoms of the original solute to wrap around
                if masks_all_related.sum() == natoms_0:
                    centroid_mass_weighed = (images.coords[masks_all_related]*images.masses[masks_all_related][:,None]).sum(axis=0)/images.masses[masks_all_related].sum()

        print('\n Mass-weighed centroid as BPC wrapping centre (in a.u.):')
        print(f' {centroid_mass_weighed[0]:.6f}  {centroid_mass_weighed[1]:.6f}  {centroid_mass_weighed[2]:.6f}\n')
        stdout.flush()
        return centroid_mass_weighed


    @property
    def _coords_transpose(self):
        '''The transpose of coordinates which acts as coordinate vectors: shape (1,3,natoms)'''

        # always calculate on the fly because coordinates may be altered anytime
        return self.coords[:,:,None].T


    @miscutils.time
    def pbc_wrap(self, centre=-1, box_parameters=[], select=[]):
        '''Wrap the system into the simulation cell (box)'''

        from numpy  import any, arange, array, dot, isin, linalg, tile, unique, vstack, zeros
        from ..base import lattice

        # PBC wrap
        if len(select) and len(box_parameters):
            # NB: vec[x/y/z, a/b/c]
            print("\n Applying PBC wrapping...")
            print("\n Cell vectors (a.u.):")
            print("           x               y               z")
            print(" a {:15.6f} {:15.6f} {:15.6f}".format(*box_parameters['A']))
            print(" b {:15.6f} {:15.6f} {:15.6f}".format(*box_parameters['B']))
            print(" c {:15.6f} {:15.6f} {:15.6f}".format(*box_parameters['C']))
            print("\n Origin (a.u.):")
            print("   {:15.6f} {:15.6f} {:15.6f}".format(*box_parameters['origin']))
            print()
            stdout.flush()
            self.cell.vectors[:,0] = box_parameters['A']
            self.cell.vectors[:,1] = box_parameters['B']
            self.cell.vectors[:,2] = box_parameters['C']
            # the centre can be an atom specified by the user
            if miscutils.isInteger(centre):
                if centre >= 0:
                    indice_centre, centre  = centre, self.coords[centre]
                    print(f' User-specified wrapping centre (in a.u.) at atom {indice_centre}:')
                    print(f' {centre[0]:.6f}  {centre[1]:.6f}  {centre[2]:.6f}\n')
                    stdout.flush()
                # if the user doesn't specify a centre, let's automatically get a wrapping centre according to the chosen atoms
                else:
                    centre = self._getPBCWrapCentre(select=select)
            # or a specific cooridate
            elif hasattr(centre, '__iter__'):
                if len(centre) == 3:
                    centre = array(centre, dtype=float)
                    print(' User-specified wrapping centre (in a.u.):')
                    print(f' {centre[0]:.6f}  {centre[1]:.6f}  {centre[2]:.6f}\n')
                    stdout.flush()
            # unknown centre
            else:
                centre = self._getPBCWrapCentre(select=select)
        # do nothing
        else:
            print("\n PBC wrapping cannot be applied because no atoms are selected")
            return

        # for convenience to use half of the volume to filter out overflown atoms, move the whole system to the centre
        self.translate(vector=-centre)

        # YL 01/11/2023: reimplemented by replacing the previous explicit loop over all atoms; now it's 10**3 faster
        def __includeResidues(_self, _masks):
            '''Also include atoms of same residues of the selected atoms'''

            if not any(_masks):
                return _masks

            # result mask array
            _masks_result = zeros((_self.natoms), dtype=bool)

            # all segments _masks involves
            _all_segIDs = unique(_self.segments[_masks])
            # loop over these segments
            for _segID in _all_segIDs:
                # get a mask array of the current segment
                _masks_seg = _self.segments == _segID
                # masks of atoms in the current segment only
                _masks_in_seg = _masks_seg*_masks
                # all resIDs confined by _masks_in_seg
                _all_resIds_in_seg = unique(_self.resIDs[_masks_in_seg])
                # all atoms of these resIDs but only in the current segment
                _masks_result += isin(_self.resIDs, _all_resIds_in_seg)*_masks_seg

            return _masks_result

        def _getDisplacement(_self, _ndims, _translation):
            '''Displacement (translation) along six different directions'''

            _displacement = zeros((3,))
            _displacement[:_ndims] = dot(_self._cell.vectors[:_ndims,:_ndims], _translation)

            return _displacement

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

        ndims = selectcases[self._cell.ndimensions]

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

        c = self._coords_transpose

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

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

        overflow_AB_p = linalg.det(stacked_AB)/volume >  0.5
        overflow_CA_p = linalg.det(stacked_CA)/volume >  0.5
        overflow_BC_p = linalg.det(stacked_BC)/volume >  0.5
        overflow_AB_n = linalg.det(stacked_AB)/volume < -0.5
        overflow_CA_n = linalg.det(stacked_CA)/volume < -0.5
        overflow_BC_n = linalg.det(stacked_BC)/volume < -0.5

        # also include atoms of same residues of the selected atoms
        overflow_AB_p = __includeResidues(self, overflow_AB_p)
        overflow_CA_p = __includeResidues(self, overflow_CA_p)
        overflow_BC_p = __includeResidues(self, overflow_BC_p)
        overflow_AB_n = __includeResidues(self, overflow_AB_n)
        overflow_CA_n = __includeResidues(self, overflow_CA_n)
        overflow_BC_n = __includeResidues(self, overflow_BC_n)

        # apply translation to wrap
        self.coords[overflow_AB_p] += _getDisplacement(self, ndims, [ 0, 0,-1])
        self.coords[overflow_CA_p] += _getDisplacement(self, ndims, [ 0,-1, 0])
        self.coords[overflow_BC_p] += _getDisplacement(self, ndims, [-1, 0, 0])
        self.coords[overflow_AB_n] += _getDisplacement(self, ndims, [ 0, 0, 1])
        self.coords[overflow_CA_n] += _getDisplacement(self, ndims, [ 0, 1, 0])
        self.coords[overflow_BC_n] += _getDisplacement(self, ndims, [ 1, 0, 0])

        # move the whole system back
        self.translate(vector=centre)


    def _getBoundingBoxCentre(self, select=[]):
        '''Return the centre of the bounding box defined by cell vectors'''

        from numpy  import arange, array, cross, dot, linalg
        from ..data import unitconvert
        from ..base import lattice

        if not len(select):
            select = ...

        coords = self.coords[select]

        selectcases = { 3 : 3,
                        2 : 2,
                        1 : 1,
                        0 : 3 }
        ndims = selectcases[self._cell.ndimensions]

        # lattice vectors
        vectors = self._cell.vectors[:ndims,:ndims]

        # unit normal vectors of planes ab, ac, and bc
        # note the right-hand rule because the distances are signed!
        plane_ab = cross(vectors[0,:], vectors[1,:])/linalg.norm(cross(vectors[0,:], vectors[1,:]))
        plane_ca = cross(vectors[2,:], vectors[0,:])/linalg.norm(cross(vectors[0,:], vectors[2,:]))
        plane_bc = cross(vectors[1,:], vectors[2,:])/linalg.norm(cross(vectors[1,:], vectors[2,:]))

        # distances from all centres to plane AB (signed!)
        distances_to_ab = dot(-coords, plane_ab)
        min_ab = distances_to_ab.min()
        max_ab = distances_to_ab.max()
        d_ab = (max_ab+min_ab)*0.5

        # distances from all centres to plane AC (signed!)
        distances_to_ac = dot(-coords, plane_ca)
        min_ac = distances_to_ac.min()
        max_ac = distances_to_ac.max()
        d_ac = (max_ac+min_ac)*0.5

        # distances from all centres to plane BC (signed!)
        distances_to_bc = dot(-coords, plane_bc)
        min_bc = distances_to_bc.min()
        max_bc = distances_to_bc.max()
        d_bc = (max_bc+min_bc)*0.5

        # we need to solve the linear system to find the point which has the distances d_ab, d_ac, and d_bc to planes ab, ac, and bc
        D  = linalg.det([-plane_ab, -plane_ca , -plane_bc])
        Dx = linalg.det([[ d_ab       , -plane_ab[1], -plane_ab[2]],
                         [ d_ac       , -plane_ca[1], -plane_ca[2]],
                         [ d_bc       , -plane_bc[1], -plane_bc[2]]])
        Dy = linalg.det([[-plane_ab[0],  d_ab       , -plane_ab[2]],
                         [-plane_ca[0],  d_ac       , -plane_ca[2]], 
                         [-plane_bc[0],  d_bc       , -plane_bc[2]]])
        Dz = linalg.det([[-plane_ab[0], -plane_ab[1],  d_ab       ], 
                         [-plane_ca[0], -plane_ca[1],  d_ac       ], 
                         [-plane_bc[0], -plane_bc[1],  d_bc       ]])

        centre = array([Dx/D, Dy/D, Dz/D])

        return centre


#    @miscutils.time
    # YL 02/12/2023: rewritten based on the new method of backbone selection using species check rather than FF names
    # YL TODO 01/12/2023: changed to protein_only=True for now because self.selectByNucleotides() also needs a reimplementation
    def selectByBackbone(self, protein_only=True, return_masks=False):
        '''Return indices of protein/nucleotide backbone atoms'''

        from numpy  import any, append, isin, unique, zeros
#        from ..data import forcefield

        # get the forcefield object
#        ffObj = forcefield.getFFObj(ff)

# YL 01/12/2023: this must be reimplemented!
#        masks_nucleotide  = self.selectByNucleotides(ff=ff, return_masks=True)

        # YL 01/12/2023
        # peptide
        #     HN HA  O  :
        #     |  |   || :
        #   --N--CA--C--:
        #        |      :
        #       (CB)    :
        # NB: do NOT use self.all_segments which makes a self._all_segments saved too soon (though this loophole could have been plugged up)!
        indices_backbone = []
        # reset these hidden feature arrays
        self.__indices_peptide_N = []
        self.__indices_peptide_C_alpha = []
        self.__indices_peptide_C = []
        self.__indices_peptide_O = []
        for seg in unique(self.segments):
            masks_seg = self.segments == seg
            indices_seg = self._indices_all[masks_seg]
            masks_A = isin(self.bonds[:,0], indices_seg)
            masks_B = isin(self.bonds[:,1], indices_seg)
            masks_bonds = masks_A*masks_B
            bonds_seg = self.bonds[masks_bonds]
            resIDs_A = self.resIDs[bonds_seg[:,0]]
            resIDs_B = self.resIDs[bonds_seg[:,1]]
            # masks array of inter- and intra-residue bonds
            masks_inter = resIDs_A != resIDs_B
            masks_intra = resIDs_A == resIDs_B
            # do not continue if no inter-residue bond
            if not any(masks_inter):
                continue
            A_inter = bonds_seg[masks_inter][:,0]
            B_inter = bonds_seg[masks_inter][:,1]
            znums_A_inter = self.znums[A_inter]
            znums_B_inter = self.znums[B_inter]
            A_intra = bonds_seg[masks_intra][:,0]
            B_intra = bonds_seg[masks_intra][:,1]
            znums_A_intra = self.znums[A_intra]
            znums_B_intra = self.znums[B_intra]

            # find inter-residue C-N bonds
            masks_inter_C_N = (znums_A_inter==6)*(znums_B_inter==7)
            masks_inter_N_C = (znums_A_inter==7)*(znums_B_inter==6)

            # all these C-end (carbonyl C) must be connected to an intra-residue O (carbonyl O) and an intra-residue C (C_alpha)
            indices_end_C = unique(A_inter[masks_inter_C_N].tolist() + \
                                   B_inter[masks_inter_N_C].tolist())
            indices_end_N = unique(B_inter[masks_inter_C_N].tolist() + \
                                   A_inter[masks_inter_N_C].tolist())

            # find all intra-residue C-C bonds
            masks_intra_C_C = (znums_A_intra==6)*(znums_B_intra==6)
            indices_C_has_C_C = unique(A_intra[masks_intra_C_C].tolist() + \
                                       B_intra[masks_intra_C_C].tolist())

            # find all intra-residue C-O bonds
            masks_intra_C_O = (znums_A_intra==6)*(znums_B_intra==8)
            masks_intra_O_C = (znums_A_intra==8)*(znums_B_intra==6)
            indices_C_has_C_O = unique(A_intra[masks_intra_C_O].tolist() + \
                                       B_intra[masks_intra_O_C].tolist())
            indices_C = unique(list(set(indices_end_C)&set(indices_C_has_C_C)&set(indices_C_has_C_O))).tolist()
            # keep a reference for potential future use
            self.__indices_peptide_C += indices_C

            # carbonyl O
            indices_O_carbonyl = B_intra[masks_intra_C_O][isin(A_intra[masks_intra_C_O],indices_C)].tolist() + \
                                 A_intra[masks_intra_O_C][isin(B_intra[masks_intra_O_C],indices_C)].tolist()
            self.__indices_peptide_O += indices_O_carbonyl

            # C_alpha
            indices_C_alpha = B_intra[masks_intra_C_C][isin(A_intra[masks_intra_C_C],indices_C)].tolist() + \
                              A_intra[masks_intra_C_C][isin(B_intra[masks_intra_C_C],indices_C)].tolist()
            self.__indices_peptide_C_alpha += indices_C_alpha

            # find all intra-residue C-N bonds
            masks_intra_C_N = (znums_A_intra==6)*(znums_B_intra==7)
            masks_intra_N_C = (znums_A_intra==7)*(znums_B_intra==6)

            # these N atoms are connected to C_alpha
            indices_N_connected_to_C_alpha = B_intra[masks_intra_C_N][isin(A_intra[masks_intra_C_N],indices_C_alpha)].tolist() + \
                                             A_intra[masks_intra_N_C][isin(B_intra[masks_intra_N_C],indices_C_alpha)].tolist()

            # and these N atoms must connect to a inter-residue C
            set_indices_N_connected_to_C_alpha = set(indices_N_connected_to_C_alpha)
            set_indices_end_N = set(indices_end_N)
            indices_N = unique(list(set_indices_N_connected_to_C_alpha & set_indices_end_N)).tolist()
            self.__indices_peptide_N += indices_N

            # now let's inspect if there're ending amino acid
            Ns_other = list(set_indices_end_N ^ set_indices_N_connected_to_C_alpha)
            indices_end = []
            # loop over the above omitted N atoms
            for N in Ns_other:
                for X in self.conn[N]:
                    if X == -1:
                        break
                    # there must be an intra-residue C (C_alpha) connected
                    if self.znums[X] == 6 and self.resIDs[N] == self.resIDs[X] and self.segments[N] == self.segments[X]:
                        for Y in self.conn[X]:
                            if Y == -1:
                                break
                            # C_alpha must be connected to an intra-residue C (carbonyl C)
                            if self.znums[Y] == 6 and self.resIDs[X] == self.resIDs[Y] and self.segments[X] == self.segments[Y]:
                                for Z in self.conn[Y]:
                                    if Z == -1:
                                        break
                                    # this carbonyl C should be connected to an O (carbonyl O)
                                    if self.znums[Z] == 8 and self.resIDs[Y] == self.resIDs[Z] and self.segments[Y] == self.segments[Z]:
                                        # these atoms belong to a peptide
                                        indices_end += [ N, X, Y, Z]
                                        # update the reference to these indices
                                        self.__indices_peptide_N += [N]
                                        self.__indices_peptide_C_alpha += [X]
                                        self.__indices_peptide_C += [Y]
                                        self.__indices_peptide_O += [Z]

            indices_backbone += indices_C + indices_O_carbonyl + indices_C_alpha + indices_N + indices_end

        masks_protein = isin(self._indices_all, unique(indices_backbone))
        stdout.flush()
        
        # YL 27/01/2022: added nucleotide
        # YL TODO 01/12/2023: we need to reimplement for nucleotides; this is a dummy
        masks_nucleotide = zeros((self.natoms,), dtype=bool)

        if protein_only:
            masks = masks_protein
        else:
            masks = masks_protein+masks_nucleotide

        if return_masks:
            return masks
        else:
            return self._indices_all[masks]


    @property
    def _indices_peptide_N(self):
        '''Return an array of the indices of peptide N atoms'''

        try:
            return unique(self.__indices_peptide_N)
        # these lists are initialised only by function selectByBackbone()
        except:
            self.selectByBackbone(protein_only=True)
            return unique(self.__indices_peptide_N)


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

        self.__indices_peptide_N = val


    @property
    def _indices_peptide_C_alpha(self):
        '''Return an array of the indices of peptide C_alpha atoms'''

        try:
            return unique(self.__indices_peptide_C_alpha)
        # these lists are initialised only by function selectByBackbone()
        except:
            self.selectByBackbone(protein_only=True)
            return unique(self.__indices_peptide_C_alpha)


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

        self.__indices_peptide_C_alpha = val


    @property
    def _indices_peptide_C(self):
        '''Return an array of the indices of peptide C atoms'''

        try:
            return unique(self.__indices_peptide_C)
        # these lists are initialised only by function selectByBackbone()
        except:
            self.selectByBackbone(protein_only=True)
            return unique(self.__indices_peptide_C)


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

        self.__indices_peptide_C = val


    @property
    def _indices_peptide_O(self):
        '''Return an array of the indices of peptide O atoms'''

        try:
            return unique(self.__indices_peptide_O)
        # these lists are initialised only by function selectByBackbone()
        except:
            self.selectByBackbone(protein_only=True)
            return unique(self.__indices_peptide_O)


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

        self.__indices_peptide_O = val


    def selectByIons(self, ff='charmm', invert=False, return_masks=False):
        '''Return indices of ionic atoms'''

        from ..data import forcefield
   
        # get the forcefield object
        ffObj = forcefield.getFFObj(ff)
   
        return self.select(residues=ffObj.ions, invert=invert, return_masks=return_masks)


    # YL 26/11/2023
    # YL 02/12/2023: this seems no longer useful since I just implemented the new selectByBackbone()
    def selectByNonProteinogenic(self, ff='charmm', invert=False, return_masks=False, mute=False):
        '''Return a list of residue names that are amino acids but with non-proteinogenic names'''

        from numpy       import any, arange, char, isin, unique, zeros
        from collections import defaultdict
        from ..data      import forcefield

        # result masks array
        masks = zeros(shape=(self.natoms,), dtype=bool)

        dict_seg_resnames = defaultdict(list)

        # get the forcefield object
        ffObj = forcefield.getFFObj(ff)

        # get the backbone masks according to atom names
        masks_backbone = self.selectByBackbone(protein_only=True, return_masks=True)

        # nothing to do if there is even no amino acid
        if not any(masks_backbone):
            if invert:
                masks = ~masks
            if return_masks:
                return masks
            else:
                return self._indices_all[masks]

        residues_backbone = char.strip(char.decode(unique(self.residues[masks_backbone])))
        # names of residues on the backbone but not in the standard set of proteinogenic amino acids
        residues_other = residues_backbone[isin(residues_backbone, ffObj.amino_acids, invert=True)]
        # YL 30/11/2023: do nothing if there is no other amino-acid-like resiude found, otherwise
        #                self.select(segments=seg, residues=residues_other) below causes chaos!
        if not len(residues_other):
            if invert:
                masks = ~masks
            if return_masks:
                return masks
            else:
                return self._indices_all[masks]

        # loop over segments
        # NB: do NOT use self.all_segments which makes a self._all_segments saved too soon (though this loophole could have been plugged up)!
        for seg in unique(self.segments):
            masks_seg = self.segments == seg
            # masks of possibly non-proteinogenic in the current segment
            masks_residues_other = self.select(segments=seg, residues=residues_other, return_masks=True)
            resIDs_other = unique(self.resIDs[masks_residues_other])
            # we should always use resID rather than resnames
            for resID_other in resIDs_other:
                is_amino_acid = True
                masks_res_other = (self.resIDs==resID_other)*masks_seg
                masks_res_other_backbone = masks_backbone*masks_res_other
                # loop over all atoms in the backbone of the current residue
                for i, idx in enumerate(self._indices_all[masks_res_other_backbone]):
                    atomname = self.names[idx].strip().decode()
                    # consider the current residue non-amino-acid if the current atom is not one of recognised peptide backbone
                    if atomname not in ffObj.peptide_bonds:
                        is_amino_acid = False
                        break
                    indices_conneced = [ x for x in self.conn[idx] if x > -1 ]
                    # check the connected atoms one by one (this doesn't exclude the case of end amino acids)
                    for index_connected in indices_conneced:
                        atomname_connected = self.names[index_connected].strip().decode()
                        resID_connected = self.resIDs[index_connected]
                        if atomname_connected in ffObj.peptide_bonds[atomname]:
                            # only intra-residue N-C and C-N are acknowledged
                            if (atomname == 'N' and atomname_connected == 'C') or \
                               (atomname == 'C' and atomname_connected == 'N'):
                                if not resID_other != resID_connected:
                                    is_amino_acid = False
                                    break
                        else:
                            is_amino_acid = False
                            break
                    # if non-amino-acid, no need to continue the loop over atoms
                    if not is_amino_acid:
                        break
                # include the current residue
                if is_amino_acid:
                    masks += masks_res_other
                    dict_seg_resnames[seg].append(self.resID2residue(resID_other, seg).strip()+str(resID_other))

        if not mute and dict_seg_resnames:
            print(f'\n Found amino acids with non-standard {ff.upper()} names (non-proteinogenic amino acids):')
            for k, v in dict_seg_resnames.items():
                print(f' Segment {k.strip().decode()}')
                print(' '+' '.join(v))
        stdout.flush()

        if invert:
            masks = ~masks
        if return_masks:
            return masks
        else:
            return self._indices_all[masks]


    # YL 02/12/2023: rewritten based on the new method of backbone selection using species check rather than FF names
    def selectByProtein(self, proteinogenic_only=False, invert=False, return_masks=False):
        '''Return indices of protein atoms'''

        from numpy  import arange, char, isin, unique
        from ..data import forcefield
   
        # get the forcefield object
#        ffObj = forcefield.getFFObj(ff)

        # get the indices of backbone atoms
        indices_backbone = self.selectByBackbone(protein_only=True)

        # then we consider all segments containing the backbone atoms belong to protein
        # NB: this is different from selectByAminoAcids() which select the residues containing backbone
        segments = unique(self.segments[indices_backbone])

        # do not check further if only proteinogenic requested
        masks = self.select(segments=segments, invert=invert, return_masks=True)

        # YL 02/12/2023: this was implemented several days ago but now seems unnecessary but let's keep it
        if not proteinogenic_only:
            masks += self.selectByNonProteinogenic(return_masks=True)

        if invert:
            masks = ~masks
        if return_masks:
            return masks
        else:
            return self._indices_all[masks]


    def selectByNucleotides(self, ff='charmm', invert=False, return_masks=False):
        '''Return indices of nucleotide atoms'''

        from numpy  import arange
        from ..data import forcefield

        # get the forcefield object
        ffObj = forcefield.getFFObj(ff)
   
        masks = self.select(residues=ffObj.nucleotides, invert=invert, return_masks=return_masks)

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


    # TWK 19/08/2022 Use PDB data to define HDLC residues for DL-FIND optimisations.
    # Somewhat similar to the pdb_to_res function in Tcl-ChemShell
    # TODO: handle cases where segments aren't present etc.
    def getHDLCResidues(self):
        '''Return residue specifications for HDLC optimisations with DL-FIND.

        Builds up HDLC residue definitions based on PDB residue data.
        Residue names are formed by combining segment names with residue numbers.

        Returns a dictionary of the form:

        { 'res1' : [atom1, atom2, atom3, ...],
          'res2' : [atom4, atom5, atom6, ...] }

        for use with the residues= argument in Opt().
        '''

        hdlcres = {}

        for iat in range(self.natoms):
            resname = "{segname}{resid}".format(segname=self.segments[iat].decode().strip(), resid=self.resIDs[iat])
            try:
                # Try appending to an existing residue
                hdlcres[resname].append(iat)
            except KeyError:
                # Residue name does not exist yet
                hdlcres[resname] = [iat]

        return hdlcres


class __PSF():
    '''PSF data'''

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

        return self._natoms


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

        self._natoms = val


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

        return self._bonds


    @bonds.setter
    def bonds(self, val):
        ''''''

        self._bonds = val


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

        try:
            return int(self._bonds.size/2)
        except:
            return 0


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

        return self._angles


    @angles.setter
    def angles(self, val):
        ''''''

        self._angles = val


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

        try:
            return int(self._angles.size/3)
        except:
            return 0


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

        return self._dihedrals


    @dihedrals.setter
    def dihedrals(self, val):
        ''''''

        self._dihedrals = val


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

        try:
            return int(self._dihedrals.size/4)
        except:
            return 0


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

        return self._impropers


    @impropers.setter
    def impropers(self, val):
        ''''''

        self._impropers = val


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

        try:
            return int(self._impropers.size/4)
        except:
            return 0


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

        return self._donors


    @donors.setter
    def donors(self, val):
        ''''''

        self._donors = val


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

        return self._nbs


    @nbs.setter
    def nbs(self, val):
        ''''''

        self._nbs = val


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

        return self._groups


    @groups.setter
    def groups(self, val):
        ''''''

        self._groups = val


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

        try:
            return self._groups.size
        except:
            return 0


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

        return self._cterms


    @cterms.setter
    def cterms(self, val):
        ''''''

        self._cterms = val


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

        try:
            return int(self._cterms.size/4)
        except:
            return 0


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

        return self._titles


    @titles.setter
    def titles(self, val):
        ''''''

        self._titles = val


    # PSF
    def delete(self, todelete):
        '''Deleting method for PSF entries'''

        from numpy import any, isin, searchsorted

        if chemsh.CHEMSH_DEBUG > 3 and len(todelete):
            print(" Deleting PSF entries involving the following atoms:\n {}".format(todelete))

        # delete bonds
        if self.nbonds:
            # remember positions to remove
            masks = ~any(isin(self.bonds,todelete).reshape((self.nbonds,2)),axis=1)
            # shift indices
            self.bonds -= searchsorted(todelete, self.bonds)
            # remove
            self.bonds = self.bonds[masks]

        # delete angles
        if self.nangles:
            masks = ~any(isin(self.angles,todelete).reshape((self.nangles,3)),axis=1)
            self.angles -= searchsorted(todelete, self.angles)
            self.angles = self.angles[masks]
            

        # delete diherals
        if self.ndihedrals:
            masks = ~any(isin(self.dihedrals,todelete).reshape((self.ndihedrals,4)),axis=1)
            self.dihedrals -= searchsorted(todelete, self.dihedrals)
            self.dihedrals = self.dihedrals[masks]
            
        # delete impropers
        if self.nimpropers:
            masks = ~any(isin(self.impropers,todelete).reshape((self.nimpropers,4)),axis=1)
            self.impropers -= searchsorted(todelete, self.impropers)
            self.impropers = self.impropers[masks]
            
        # delete cross-terms
        if self.ncterms:
            masks = ~any(isin(self.cterms,todelete).reshape((self.ncterms,4)),axis=1)
            self.cterms -= searchsorted(todelete, self.cterms)
            self.cterms = self.cterms[masks]
            







