#  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 sys              import stdout
from ...io            import format
from ...io.files      import rtf
from ...dl_py2f.utils import dictutils, fileutils, iterutils, miscutils, strutils

# see NAMD User's Guild Chapter 4
@miscutils.time
def psfgen(bioFrag,
           allow_undefined=False,
           atom_aliases={},
           dryrun=False,
           ff='charmm',
           ff_dir='',
           guess_coord=True,
           haem='',                # the user can choose a different segment name such as 'HEME' or even a sequence like [ 'HEMB', 'HEMC' ]
           input='_psfgen.inp',
           ions_in_solvent=False,
           merge_solvents=False,
           no_auto=[],
           output='_psfgen.out',
           output_pdb='',
           output_psf='',
           output_dir='_psfgen',
           override_standard=False,
           patch=True,
           patch_disulphides=True,
           patch_titratables=True,
           path='',
           residue_aliases={},
           solvents=[],
           topologies=[],
           use_newer_version=True):
    '''Run psfgen and generate PSF from PDB'''

    from numpy     import array, char, hstack, isin, ndarray, setdiff1d, str_, unique, zeros
    from copy      import copy
    from os.path   import abspath, isdir, isfile, join, relpath
    from itertools import groupby
    from operator  import itemgetter
    from .         import haems
    from ...io     import banner, format
    from ...data   import forcefield

    # fetch the FF data module
    ffObj = forcefield.getFFObj(ff)

    # YL 26/11/2021: strangely we have to recalculate the connectivity now
    bioFrag.connect(force_reconnect=True)

    # YL 20/09/2021: according to PSFGEN (https://www.ks.uiuc.edu/Research/vmd/plugins/psfgen/ug.pdf)
    #                Many PDB files in the PDB databank contain multiple chains, corresponding to protein subunits, water,
    #                and other miscellaneous groups. Protein subunits are often identified by their chain ID in the PDB file. In
    #                psfgen, each of these groups must be assigned to their own segment. This applies most strictly in the case
    #                of protein chains, each of which must be assigned to its own segment so that N-terminal and C-terminal
    #                patches can be applied. You are free to group water molecules into whatever segments you choose.
    # select all non-amino-acid atoms
    masks_non_protein = bioFrag.selectByProtein(return_masks=True, invert=True)
    segments_created = []
    residue_patches  = {}
    # loop over segments
    for segment in bioFrag.all_segments:
        # select the current segment
        masks_segment = bioFrag.segments==segment
        # proceed iff the current segment contains amino acid
        if (masks_segment*(~masks_non_protein)).sum():
            masks_non_aa_segment = masks_segment*masks_non_protein
            residues_non_protein_unique = unique(bioFrag.residues[masks_non_aa_segment])
            segments_created += char.decode(residues_non_protein_unique).tolist()
            # rename
            # YL TODO 09/02/2022: DNA/RNA?
            for residue_non_protein in residues_non_protein_unique:
                # YL 09/02/2022: added but not in use; in theory ions should be optionally excluded in the solute segment
                if ions_in_solvent and residue_non_protein.decode().strip() in ffObj.ions: 
                    continue
                masks_seleted = masks_non_aa_segment*(bioFrag.residues==residue_non_protein)
                bioFrag.segments[masks_seleted] = residue_non_protein
    # recalculate unique segments
    bioFrag._all_segments, bioFrag._all_segments_inverse_indices = unique(bioFrag.segments, return_inverse=True)
    # resIDs must be rearranged
    for segment in bioFrag.all_segments:
        bioFrag.optimiseResIDs(segment)
    if len(segments_created):
        print('\n NB: PSFGEN requires proteins to be assigned to their own segment')
        print('     Segments were created for the following non-amino-acid residues:\n     {}'.format(' '.join(char.decode(residues_non_protein_unique))))
        stdout.flush()

    # result to return
    additional_prm = []

    # welcome
    banner.printModuleBanner('ChemShell/NAMD Interface: PSFGEN')

    # YL FIXME 20/09/2021: this is hardcoded
    if not solvents:
        if bioFrag.select(segments='CW').size:
            solvents += ['CW']
        # segments starting with 'SOL' (solvated by Py-ChemShell)
        try:
            solvents += char.decode(unique(bioFrag.segments[char.startswith(char.lower(bioFrag.segments), b'sol')])).tolist()
        except TypeError:
            solvents += unique(bioFrag.segments[char.startswith(char.lower(bioFrag.segments), 'sol')]).tolist()
        print("\n Solvent segments detected:")
        format.printwrap(' '.join(solvents), prepend_blanc_line=False)

    # group the solvents, according to topology files, to match the returned dict by bioFrag.split()
    group_dict = dict(zip([ solv.strip() for solv in solvents ], [ True for solv in solvents ]))

    solv_segnames = dict()
    if merge_solvents:
        solv_groups = solvents
        print("\n Solvent segments are grouped:\n", solv_groups)
        print("\n WARNING: merging solvent segments may cause resIDs exceeding the 9999 limit and deletion by psfgen")
        # update topologies dict
        solv_dicts = [ {'_'.join(s.strip() for s in list(map(itemgetter(0),v))):k}
                                for k,v in groupby(sorted(group_dict.items(),key=itemgetter(1)),itemgetter(1)) ]
        # however segname can't exceed 4 characters
        for i,d in enumerate(solv_dicts):
#            topologies.update(d)
            solv_segnames.update(dict.fromkeys(d, 'SOL{:1}'.format(str(i))))
    else:
        solv_groups = []

    # check forcefield directory
    ff_dir = forcefield.getFFDir(ff, dir=ff_dir)

    # update residue alias rules
    built_in_rules = ffObj.residue_alias
    rules_to_add   = {}
    for res, alias in built_in_rules.items():
        selected = bioFrag.select(residues=res)
        if len(selected):
            rules_to_add.update({res:alias})
    residue_aliases.update(rules_to_add)
    if residue_aliases:
        print("\n Residue aliases of forcefield scheme {}:\n".format(ff.upper()))
    for res, alias in residue_aliases.items():
        print(" {:5} -> {}".format(res, alias))

    # segments without repeat
    segments = unique(bioFrag.segments)

    # default: no automatic angles/dihedrals on solvents, unless specified otherwise
    if no_auto == []:
        no_auto = solvents

    # split into segments
    # TODO: undefined
    undefined = []
    split_frags = bioFrag.split(exclude=undefined, groups=solv_groups, return_dict=False, sort=False)

    proteins = bioFrag.selectByProtein()
    if proteins.size:
        segnames_proteins = char.lower(char.decode(unique(bioFrag.segments[proteins])))
    else:
        segnames_proteins = []

    # YL TODO 18/12/2023: the function bioFrag.selectByNucleotides() is pending reimplementation
    nucleotides = bioFrag.selectByNucleotides(ff=ff)
    if nucleotides.size:
        segnames_nucleotides = char.lower(char.decode(unique(bioFrag.segments[nucleotides])))
    else:
        segnames_nucleotides = []

    # a list of segnames: place the solvent segments in the end; protein segment(s) must be dealt with first
    segnames = [ i for i in split_frags.dtype.names if i not in solvents and i.lower().strip() in segnames_proteins ] \
             + [ i for i in split_frags.dtype.names if i not in solvents and i.lower().strip() in segnames_nucleotides ] \
             + [ i for i in split_frags.dtype.names if i not in solvents and \
                (i.lower().strip() not in segnames_proteins and   \
                 i.lower().strip() not in segnames_nucleotides) ] \
             + [ i for i in split_frags.dtype.names if i in solvents ]

    # get a list of absolute paths to topology files
    topol_abspaths = getAbsPathsList(topologies, ff_dir)

    # get a dict of dicts containing residue-atoms of all topology files
    rtf_res_dicts = rtf.getResidueDicts(topol_abspaths)

    # detect special segments (such as haem) and add missing RTF/PRM
    haem_results = haems._detectHaem(bioFrag, ff, split_frags, segnames, haem, \
                                     rtf_res_dicts, topol_abspaths, residue_aliases, additional_prm)

    # special rules for haems
    if haem_results is not None:
        res_haem = haem_results[0]
        segnames_haems = char.strip(char.lower(char.decode(unique(bioFrag.segments[bioFrag.select(residues=res_haem)]))))
        # YL 16/12/2021: this was a substantial bug!--we must place haems after other segments (excl. solvents)
        #                especially proteins and ligands otherwise patching in PSFGEN cannot work
        segnames = [ i for i in split_frags.dtype.names if i not in solvents and i.lower().strip() in segnames_proteins ] \
                 + [ i for i in split_frags.dtype.names if i not in solvents and i.lower().strip() in segnames_nucleotides ] \
                 + [ i for i in split_frags.dtype.names if i not in solvents and \
                    (i.lower().strip() not in segnames_proteins and   \
                     i.lower().strip() not in segnames_haems    and   \
                     i.lower().strip() not in segnames_nucleotides) ] \
                 + [ i for i in split_frags.dtype.names if i.lower().strip() in segnames_haems ] \
                 + [ i for i in split_frags.dtype.names if i in solvents ]
        # switch off disulphides patching
        if 'HEMC' in residue_aliases.values():
            patch_disulphides = False
            print("\n Switching off automatic disulphides patching to protect the CYS residues")

    # we must know patched residue names beforehand to scan through the RTF files
    if patch:
        for iseg, seg in enumerate(segnames):
            segFrag = split_frags[seg][0]
            # YL 27/05/2024: label all segments with chainID
            chainID = unique(segFrag.chainIDs)[0].decode()
            segname = seg.rstrip(f'_{chainID}')
            __patch(bioFrag, segname, segFrag, ffObj, patch_disulphides, haem, haem_results, residue_patches, mute=True)

    # get a dict of mapping from segments to toplogy files
    seg_top_dict = __getSegTopDict(ff, split_frags, topol_abspaths, rtf_res_dicts, additional_prm, residue_patches, _residue_aliases=residue_aliases, _override_standard=override_standard, _use_newer_version=use_newer_version)

    print("\n\n Topology files associated with segements:")
    max_strlen_seg = max(char.str_len(list(seg_top_dict.keys())))
    for seg in segnames:
        if type(seg_top_dict[seg]) is str:
            print((' {:'+str(max_strlen_seg)+'}: {}').format(seg, seg_top_dict[seg]))
        elif type(seg_top_dict[seg]) is list:
            print((' {:'+str(max_strlen_seg)+'}: {}').format(seg, seg_top_dict[seg][0]))
            for i in seg_top_dict[seg][1:]:
                print((' {:'+str(max_strlen_seg)+'}  {}').format(' ', i))
    stdout.flush()

    ### PSFGEN INPUT ###

    # string buffer
    strbuff = '# psfgen input generated by Py-ChemShell\n'
    strbuff += 'resetpsf\n'

    # read in all topology files
    strbuff += '\n# read in topology files\n'
    strbuff += '\n'.join(['topology {}'.format(abspath(t)) for t in topol_abspaths]) + '\n'

    segments_to_archive = []

    # loop over segments
    for iseg, seg in enumerate(segnames):

        segFrag = split_frags[seg][0]
        if not segFrag:
            print(" >>> WARNING: No segment matched the requested \"{}\". Skipping...\n".format(seg))
            continue

        print('\n '+'-'*70)
        print(" Processing segment \"{}\"...".format(seg))
        print(' '+'-'*70)
        stdout.flush()

        # YL 18/09/2021: allow multiple topology files for a segment (rtf_res_dict is only needed for aliases)
        try:
            topology_file = seg_top_dict[seg]
            # single file
            if type(topology_file) is str:
                rtf_res_dict = rtf_res_dicts[topology_file]
                print(" RTF file: {}".format(abspath(topology_file)))
                stdout.flush()
            # multiple files (saved in a list)
            elif type(topology_file) is list:
                # construct a new rtf_res_dict for later use
                rtf_res_dict = dict()
                # first print all files
                print(" RTF file: {}\n".format(abspath(topology_file[0])), end='')
                stdout.flush()
                for top in topology_file[1:]:
                    print('           {}'.format(abspath(top)))
                # take the first entry's data anyway
                rtf_res_dict.update(rtf_res_dicts[topology_file[0]])
                # loop over files
                for top in topology_file[1:]:
                    # if the current file is a user-defined one
                    if top not in ffObj.default_rtfs and top not in ffObj.builtin_rtfs:
                        user_specified = True
                    else:
                        user_specified = False
                    # for printing later
                    residues_duplicate = []
                    entries_duplicate = {}
                    # loop over source data dict items to check: k is keyword for entry name and v is data
                    for k, v in rtf_res_dicts[top].items():
                        # array entries should be handled carefully
                        if type(v) is ndarray:
                            # select data before merging if keyword already exists in dict (e.g., 'mass')
                            if k in rtf_res_dict:
                                # first entry such as 'residues', 'mass', 'decl', 'defa'
                                field = rtf_res_dict[k].dtype.names[0]
                                # user-specified file has priority
                                if user_specified:
                                    # keep in mind: rtf_res_dict[k] is destination and v is source
                                    masks_existing_in_user = isin(rtf_res_dict[k][field], v[field])
                                    if rtf_res_dict[k][field].shape == ():
                                        # for printing only
                                        if rtf_res_dict[k][field].dtype.type is str_:
                                            entries_duplicate.update({k:copy(rtf_res_dict[k][field][None][~masks_existing_in_user])})
                                        # remove duplicate existing entries from rtf_res_dict and merge
                                        rtf_res_dict.update({k:hstack((rtf_res_dict[k][None][~masks_existing_in_user],v))})
                                    else:
                                        if rtf_res_dict[k][field].dtype.type is str_:
                                            entries_duplicate.update({k:copy(rtf_res_dict[k][field][~masks_existing_in_user])})
                                        # remove duplicate existing entries from rtf_res_dict and merge
                                        rtf_res_dict.update({k:hstack((rtf_res_dict[k][~masks_existing_in_user],v))})
                                # other files have lower priority
                                else:
                                    masks_current_in_existing = isin(v[field], rtf_res_dict[k][field])
                                    if v[field].shape == ():
#                                        # for printing
#                                        if v[field].dtype.type is str_:
#                                            entries_duplicate.update({k:copy(v[field][None][~masks_current_in_existing])})
                                        # remove duplicate entries from rtf_res_dicts[top] and merge
                                        rtf_res_dict.update({k:hstack((v[None][~masks_current_in_existing],rtf_res_dict[k]))})
                                    else:
#                                        if v[field].dtype.type is str_:
#                                            entries_duplicate.update({k:copy(v[field][~masks_current_in_existing])})
#                                        # remove duplicate entries from rtf_res_dicts[top] and merge
                                        rtf_res_dict.update({k:hstack((v[~masks_current_in_existing],rtf_res_dict[k]))})
                            # merge directly
                            else:
                                rtf_res_dict.update({k:v})
                        # mainly dict type
                        else:
                            if k in rtf_res_dict:
                                # YL 18/09/2021: residues defined in user-specified RTF should always be respected
                                if user_specified:
                                    rtf_res_dict.update({k:v})
                                    residues_duplicate += [k]
                            else:
                                rtf_res_dict.update({k:v})
                    if len(residues_duplicate):
                        print(' NB: Duplicate residues')
                        format.printwrap('{}'.format(' '.join(residues_duplicate)),
                                         indent=' '*5, prepend_blanc_line=False)
                        print(f'     kept from file {top}')
                    if entries_duplicate:
                        print(' NB: Duplicate data overwritten by file')
                        print(f'     {top}')
                        print('     Entry      Data')
                        for k, v in entries_duplicate.items():
                            print('     {:<8}'.format(k), end='   ')
                            format.printwrap(' '.join(v), indent='', initial_indent='', subsequent_indent=' '*16, prepend_blanc_line=False)
        # invalid RTF
        except:
            if allow_undefined:
                print(" No topology file specified for segment \"{}\". Skipping...".format(seg))
                continue
            else:
                print(" No topology file specified for segment \"{}\". Exiting...".format(seg))
                stdout.flush()
                quit(999)
        stdout.flush()

        # apply residue aliases
        aliased_res = iterutils.findMapped(char.strip(unique(segFrag.residues)), residue_aliases, case_sensitive=False)
        if len(aliased_res):
            strbuff += '\n# residue aliases specified via `residue_aliases`\n'
        for r in aliased_res:
            strbuff += 'pdbalias residue {} {}\n'.format(r, residue_aliases[r])

        # build segment
        if merge_solvents and solv_segnames.get(seg, seg) != seg:
            print(" Segment {} renamed as {}".format(seg, solv_segnames.get(seg, seg)))
            stdout.flush()
        strbuff += '\n# build segment \"{}\"\n'.format(seg)
        strbuff += 'segment {} {{\n'.format(solv_segnames.get(seg, seg))

        if strutils.inOrEq(seg, no_auto):
            print(" Automatic generating of angles and dihedrals for segment \"{}\" switched off".format(seg))
            stdout.flush()
            strbuff += '    auto none\n'

        # build segment (continued)
        strbuff += '    pdb {}\n'.format('_segment_{}.pdb'.format(seg))
        strbuff += '}\n'

        segments_to_archive.append('_segment_{}.pdb'.format(seg))

        # patch residues
        if patch:
            # YL 27/05/2024: label all segments with chainID
            chainID = unique(segFrag.chainIDs)[0].decode()
            segname = seg.rstrip(f'_{chainID}')
            strbuff += __patch(bioFrag, segname, segFrag, ffObj, patch_disulphides, haem, haem_results, residue_patches)

        # apply automatic atom aliases
        if not atom_aliases:
            aliases, patches = autoAlias(ff, segFrag, rtf_res_dict, residue_patches, residue_aliases=residue_aliases, patch_titratables=patch_titratables)
            for patch_res, patch_resIDs in patches.items():
                for patch_resID in patch_resIDs:
                    strbuff += 'patch {} {}:{}\n'.format(patch_res, seg, patch_resID)
            for k, v in aliases.items():
                strbuff += '\n# automatic alias generated by ChemShell\n'
                try:
                    strbuff += 'pdbalias atom {} {} {}\n'.format(k.decode(), v[0], v[1])
                except:
                    strbuff += 'pdbalias atom {} {} {}\n'.format(k, v[0], v[1])

        # read coordinates from PDB
        strbuff += '\n# read {} coordinates from PDB file\n'.format(seg)
        strbuff += 'coordpdb {} {}\n'.format('_segment_{}.pdb'.format(seg), solv_segnames.get(seg, seg))

        # end loop over segments
        strbuff += '\n'
        print(' '+'-'*70)
        stdout.flush()


    if guess_coord:
        strbuff += '\n# guess missing coordinate\n'
        strbuff += 'guesscoord\n'

    # write PSF and PDB
    if not output_psf:
        output_psf = '_'.join(bioFrag.title.strip().split())+'_psfgen.psf'
    if not output_pdb:
        output_pdb = '_'.join(bioFrag.title.strip().split())+'_psfgen.pdb'
    strbuff += '\n# write out PSF and PDB files\n'
    strbuff += 'writepsf {}\n'.format(join(output_dir, output_psf))
    strbuff += 'writepdb {}\n'.format(join(output_dir, output_pdb))

    # archive old output dir
    fileutils.archiveDir(output_dir)
    with open(input, 'w') as fp:
        fp.write(strbuff)

    if not dryrun:
        ierror = __run(path, input, output, output_dir, output_psf, output_pdb)

    # archive input/output
    for segment_to_archive in segments_to_archive:
        fileutils.archiveFile(segment_to_archive, dir=output_dir, suffix='')
    fileutils.archiveFile(input , dir=output_dir, suffix='')
    fileutils.archiveFile(output, dir=output_dir, suffix='')
    fileutils.archiveFile(f'{output_dir}.err', dir=output_dir, suffix='')


    stdout.flush()
    # NAMD needs locations of generated PSF and PDB
    return join(output_dir,output_psf), join(output_dir,output_pdb), additional_prm, solvents

    ### END PSFGEN INPUT ###


def __patch(_bioFrag, _seg, _segFrag, _ffObj, _patch_disulphides, _haem, _haem_results, _residue_patches, mute=False):
    '''Patch residues'''

    from numpy import char, isin, unique
    from .     import haems, nucleotides

    _strbuff = ''

    # YL 27/05/2024: label all segments with chainID
    _chainID = unique(_segFrag.chainIDs)[0].decode()
    # patch special residues
    residues_seg = char.decode(char.strip(unique(_segFrag.residues)))
    # prescreen the residues because looping over a segment may take a lot of time
    if isin(residues_seg, _ffObj.patchable_residues).sum():
        resIDs_seg = unique(_segFrag.resIDs)
        for resID_seg in resIDs_seg:
            # str result is returned
            residue_seg = _bioFrag.resID2residue(resID_seg, segment=_seg).strip()
            if residue_seg in _ffObj.patchable_residues:
                if not mute:
                    print(f"\n PSFGEN: Found patchable residue:\n segment: {_seg}  residue: {residue_seg}  resID: {resID_seg}")
                # YL 04/12/2021: take patched residue names as aliases for RTF scan
                _residue_patches.update({residue_seg:_ffObj.patchable_residues[residue_seg]})
                connected_resIDs, connected_segments = _bioFrag.getConnectedResIDsFromResID(resIDs_seg, segment=_seg)
                for i, connected_resID in enumerate(connected_resIDs):
                    _strbuff += 'patch {} {}:{} {}:{}\n'.format(_ffObj.patchable_residues[residue_seg],
                                                                 connected_segments[i]+f'_{_chainID}',
                                                                 connected_resID,
                                                                _seg+f'_{_chainID}',
                                                                 resID_seg)

    # patch disulphides
    # TODO support user specified
    if _patch_disulphides == True:
        disulphides = _segFrag.getDisulphides(returnResIDs=True)
        if disulphides.size:
            _strbuff += '\n# automatic residue patch generated by ChemShell\n'
            # YL 04/12/2021: take patched residue names as aliases for RTF scan
            _residue_patches.update({'CYS':'DISU'})
            for disulphide in disulphides:
                _strbuff += 'patch DISU {}:{} {}:{}\n'.format(_seg+f'_{_chainID}', disulphide[0], _seg+f'_{_chainID}', disulphide[1])

    # patch haem
    _strbuff += haems._patchHaem(_haem, _bioFrag, _ffObj, _haem_results, _seg, _chainID, residues_seg, _residue_patches, patch_disulphides=_patch_disulphides, mute=mute)

    # patch nucleotides
    _strbuff += nucleotides._patchNucleotides(_bioFrag, _ffObj, _seg, _segFrag, _residue_patches, mute=mute)

    return _strbuff


def __run(path, input, output, output_dir, output_psf, output_pdb):
    '''Execute psfgen'''

    from os.path     import abspath, isdir, isfile, join, relpath
    from subprocess  import run 
    from ...io       import format
    from ...parallel import parallel

    sysmpi = ''

    if not path:
        try:
            from . import namd_path
        except:
            format.printwrap('>>> ChemShell/PSFGEN ERROR: NAMD not built in the current ChemShell installtion. Please recompile with the `./setup --namd ...` option, or specify a full path to the PSFGEN binary via `psfgen(path=...)`')
            exit(112)
        path = join(namd_path.namd_path, 'psfgen')
        sysmpi = namd_path.sysmpi

    if input and not isfile(input):
        format.printwrap(">>> ChemShell/PSFGEN ERROR:", filename, "does not exist. Please check your input.\n")
        exit(999)

    # YL 23/11/2021: on some systems (Cray's ARCHER2!!) calling to PSFGEN in parallel results in error:
    #                ...:_pmi_mmap_tmp: Warning bootstrap barrier failed: num_syncd=X, pes_this_node=YY, timeout=180 secs
    #                then we have to force an `srun -n 1` before the PSFGEN command
#    if parallel.nprocs_per_workgroup > 1 and sysmpi:
#        command = ' '.join([ f'{sysmpi} -n 1 -N 1 --oversubscribe', path, input ])
#    else:
    # YL 20/12/2021: on the ARCHER2 full system the problem is gone
    command = ' '.join([ path, input ])

    miscutils.runSysCmd(command, stdout=output, stderr=f'{output_dir}.err', shell=False, check=True)
# alternatively
#    command = path
#    miscutils.runSysCmd(command, stdinp=input, stdout=output, stderr=f'{output_dir}.err', shell=False, check=True, comm=True)

    # YL 16/12/2021: check errors in output
    with open(output) as fp:
        lines = fp.readlines()
        header = ''
        lines_failed = []
        for line in lines:
            if 'failed' in line:
                header = '\n WARNING: There are at least one error in the PSFGEN run:\n'
                lines_failed += [' '+line]
        if header:
            print(''.join([header]+lines_failed))
        stdout.flush()

    return abspath(join(output_dir,output_psf)), abspath(join(output_dir,output_pdb))


def getUndefinedResidues(bioFrag, segment, topfile, residue_aliases):
    ''''''

    from numpy      import any, append, array, char, isin, unique
    from ...io.file import parse

    # get RTF defined residues
    rtfresdict = dict()
    parse(topfile, fmt='rtf', receiver=rtfresdict)

    # get residues of the current segment
    residues = unique(bioFrag.residues[bioFrag.select(segments=segment)])

    undefined = array([])

    for res in residues:
        # override exsisting definitions if residue aliases are defined
        try:
            oldres = iterutils.findMapped(res, residue_aliases, case_sensitive=False)[0]
            newres = residue_aliases[oldres]
            if not dictutils.isKeyOf(rtfresdict, newres, case_sensitive=False):
                append(undefined, [newres])
        except:
            if not dictutils.isKeyOf(rtfresdict, res, case_sensitive=False):
                undefined = append(undefined, [res])

    return undefined


# YL 14/09/2021: undefined_res is an array of residue names of the current segment yet to be found in RTF 
def _isAliasedResInRTF(__segname, __bioFrag, __rtf_res_dict, __dict_top_undefined_res, residue_aliases={}, mute=False):
    ''''''

    from numpy    import array, char, isin, unique

    # results
    found = False

    # byte
    residues = unique(char.strip(__bioFrag.residues))

    for res in residues:

        # override existing definitions if residue aliases are defined
        try:
            # bytes
            oldres = iterutils.findMapped(res, residue_aliases, case_sensitive=False)[0]
            newres = residue_aliases[oldres]
            oldres, res = res, newres.encode()
        except:
            pass

        #
        try:
            rtfatoms = dictutils.getValue(__rtf_res_dict, res, case_sensitive=False)
            found = True
            print(" Residue alias {}->{} found in RTF".format(oldre.decode(), res.decode()))
            # YL 14/09/2021: update __dict_top_undefined_res
            unmatched_res = __dict_top_undefined_res.get(__segname, array([], dtype='<U5'))
            unmatched_res = unmatched_res[~isin(unmatched_res, oldres.decode())]
            if len(unmatched_res):
                __dict_top_undefined_res.update({__segname:unmatched_res})
            # remove entry if no more residues unfound
            else:
                __dict_top_undefined_res.pop(__segname)
        except:
            if not mute:
                print(" Residue \"{}\" not defined in RTF. Skipping...".format(res.decode()))
            continue

    return found


def autoAlias(ff, bioFrag, rtf_res_atm_dict, residue_patches, residue_aliases={}, patch_titratables=True, mute=False):
    ''''''

    from difflib    import SequenceMatcher
    from collections import namedtuple
    from numpy      import any, append, array, char, isin, unique
    from ...io.file import parse
    from ...io      import format

    # results
    aliases = dict()
    patches = dict()
    Result = namedtuple('Result', ['aliases', 'patches'])

    # byte
    residues = unique(char.strip(bioFrag.residues))

    for res in residues:

        atomnames = char.decode(bioFrag.names[bioFrag.select(residues=res)])
        oldatoms = unique(atomnames)
        # override exsisting definitions if residue aliases are defined
        try:
            oldres = iterutils.findMapped(res, residue_aliases, case_sensitive=False)[0]
            newres = residue_aliases[oldres]
            if not mute:
                print(" Using residue aliases specified via `residue_aliases`: {}->{}".format(res.decode(), newres))
            # YL 26/11/2021: keep a record
            bioFrag._renamed_residues.update({res.decode():newres})
            stdout.flush()
            oldres, res = res, newres.encode()
            # the old resname should be used if there is an atom alias (see below)
            res_atom_alias = oldres
        except:
            res_atom_alias = res
            oldres = res

        #
        try:
            rtfatoms = dictutils.getValue(rtf_res_atm_dict, res, case_sensitive=False)
            # get the structured array out from the dict
            try:
                rtfatoms = rtfatoms['atom']['key']
                if chemsh.CHEMSH_DEBUG > 4:
                    print(" RTF atom keys defined in residue {}:\n {}".format(res, rtfatoms.decode()))
            except:
                pass
        except:
            if not mute:
                print(" Residue \"{}\" not defined in RTF. Skipping...".format(res.decode()))
            continue

        # YL 26/01/2022: also get the atoms of the patch if there is one associated to the current residue
        patch = residue_patches.get(oldres.decode(), '')
        atoms_patch = array([])
        if patch and patch in rtf_res_atm_dict:
            atoms_patch = dictutils.getValue(rtf_res_atm_dict, patch, case_sensitive=False)
            # get the structured array out from the dict
            try:
                atoms_patch = atoms_patch['atom']['key']
                if chemsh.CHEMSH_DEBUG > 4:
                    print(" RTF atom keys defined in patch residue {}:\n {}".format(patch, atoms_patch.decode()))
            except:
                pass

        rtfatoms = append(rtfatoms, atoms_patch)

        notinrtf = isin(oldatoms, rtfatoms, assume_unique=True, invert=True)
        masks_undefined_atoms = isin(atomnames, rtfatoms, invert=True)

        # if there are undefined atoms in the current residue
        if any(notinrtf):

            notfound = oldatoms[notinrtf]
            # patch titratables
            if patch_titratables == True:
                patch_residue = __getTitratables(ff, res, rtf_res_atm_dict, notfound)
                if patch_residue:
                    patch_resIDs = bioFrag.residue2resIDs(res)
                    patches.update({patch_residue:patch_resIDs})
                    continue

        # YL 17/11/2020: added check for terminal residue which should not be aliased
        # if there are undefined atoms in the current residue
#        if any(masks_undefined_atoms):
            # YL 20/03/2021: should search oldres rather than res
            indices_undefined_atoms = bioFrag.select(atoms=atomnames[masks_undefined_atoms], residues=oldres)
#            indices_undefined_atoms = bioFrag.select(atoms=atomnames[masks_undefined_atoms], residues=res)
            segname = unique(bioFrag.segments[indices_undefined_atoms])[0]
            resID   = unique(bioFrag.resIDs[indices_undefined_atoms])[0]
            selected_segment = bioFrag.select(segments=segname)
            min_res = bioFrag.resIDs[selected_segment].min()
            max_res = bioFrag.resIDs[selected_segment].max()
            
            # NB: check if N or O: in theory only N/O should not be aliased but because
            # other atoms are difficult to find matching names we do not alias at all
#            for i in indices_undefined_atoms:
#                if bioFrag.znums[i] in [ 7, 8 ] and resID in [ min_res, max_res ]:
            if resID in [ min_res, max_res ]:
                continue

            if notinrtf.sum() == 1:
                print(" Atom name \"{}\" of residue \"{}\" not found in RTF which contains:".format(notfound[0], res.decode()), end='\n')
            elif notinrtf.sum() > 1:
                print(" Atom names \"{}\" of residue \"{}\" not found in RTF which contains:".format('\", \"'.join(map(str,notfound)), res.decode()), end='\n')
            stdout.flush()
            if rtfatoms.shape == ():
                format.printwrap("{}".format(' '.join(rtfatoms[None])), prepend_blanc_line=False)
            else:
                format.printwrap("{}".format(' '.join(rtfatoms)), prepend_blanc_line=False)
            for oldatom in notfound:
                if rtfatoms.shape == ():
                    similarities = array([SequenceMatcher(None, oldatom, rtfatom).ratio() for rtfatom in rtfatoms[None]])
                    newatom = rtfatoms[None][similarities.argmax()]
                else:
                    similarities = array([SequenceMatcher(None, oldatom, rtfatom).ratio() for rtfatom in rtfatoms])
                    newatom = rtfatoms[similarities.argmax()]
                print(" \"{}\" used for automatic alias: {}->{}".format(newatom, oldatom, newatom))
                # datatype of result: {str:(str,str)}
                # use old residue name if there is alias for it
                aliases.update({res_atom_alias:(oldatom,newatom)})
#                try:
#                    aliases.update({oldres:(oldatom,newatom)})
#                except:
#                    aliases.update({res.decode():(oldatom,newatom)})

    result = Result(aliases=aliases, patches=patches)

    return result
                

def __getTitratables(_ff, residue, RTFResDict, atoms_not_found):
    ''''''

    from numpy    import all, isin, unique
    from ...data  import forcefield

    ffObj = forcefield.getFFObj(_ff)

    patch = dictutils.getValue(ffObj.titratables, residue, case_sensitive=False, raiseError=False)

    patoms = dictutils.getValue(RTFResDict, patch, case_sensitive=False, raiseError=False)

    # confirm if all missing atoms are found in the patch residue
    if all(isin(atoms_not_found, patoms)):
        return patch

    else:
        return ''


def __getSegTopDict(_ff, _split_frags, _topol_abspaths, _rtf_res_dicts, _additional_prm, _residue_patches, _residue_aliases={}, _override_standard=False, _use_newer_version=True):
    '''Assign segments to topology files and return a dict'''

    from numpy   import all, append, argwhere, array, char, isin, vstack, unique, zeros
    from copy    import copy
    from os      import listdir, sep
    from os.path import dirname, join
    from ...data import forcefield
    from ...io   import format

    ffObj = forcefield.getFFObj(_ff)

    # associate segments to topology files
    _seg_top_dict = dict()
    _dict_seg_undefined_res = dict()
    _top_seg_map = zeros((len(_topol_abspaths),len(_split_frags.dtype)), dtype=bool)

    if chemsh.CHEMSH_DEBUG > 2:
        print("\n Analysing residues of segments:")
    if chemsh.CHEMSH_DEBUG > 4:
        print(" __getSegTopDict()")

    # step 1: search in user-specified files
    # loop over segments
    for _iseg, _seg in enumerate(_split_frags.dtype.names):
        _segFrag = _split_frags[_seg][0]
        _residues = char.strip(char.decode(unique(_segFrag.residues)))
        # YL 13/09/2021: all residues have to be found
        masks_residues_found = zeros(_residues.shape, dtype=bool)
        if chemsh.CHEMSH_DEBUG > 2:
            print("\n {}:\n {}".format(_seg, _residues))
        # loop over topology files
        for _itop, _topol_abspath in enumerate(_topol_abspaths):
            _rtf_res_dict = _rtf_res_dicts[_topol_abspath]
#            if chemsh.CHEMSH_DEBUG > 4:
#                print("\n Residues contained by topology file:\n {}\n {}".format(_topol_abspath, list(_rtf_res_dict.keys())))
            # any matching residue
            masks_residues_found_in_rtf = isin(_residues, array(list(_rtf_res_dict.keys())))
            masks_residues_found += masks_residues_found_in_rtf
            _top_seg_map[_itop,_iseg] = masks_residues_found_in_rtf.sum()
            # YL 16/09/2021: update segment-topology dict
            if masks_residues_found_in_rtf.sum():
                _seg_top_dict.update({_seg:_topol_abspath})
        # save unmatched seg/res to the dict for later use
        if not all(masks_residues_found):
            _dict_seg_undefined_res.update({_seg:_residues[~masks_residues_found]})
    if _dict_seg_undefined_res:
        print('\n The following segment:residue(s) are not found in the given RTF:')
        for k, v in _dict_seg_undefined_res.items():
            format.printwrap('{:4}: '.format(k)+' '.join(v),
                             indent='', initial_indent=' ', subsequent_indent=' '*7, prepend_blanc_line=False)

    _list_seg_user_undefined = list(copy(_dict_seg_undefined_res).keys())

    # number of matches
    matched_segs = _top_seg_map.sum(axis=0)
    matched_tops = _top_seg_map.sum(axis=1)
    # step 2: update the mapping using aliases
    for _iseg, _seg in enumerate(_split_frags.dtype.names):
        _segFrag = _split_frags[_seg][0]
        for _itop, _topol_abspath in enumerate(_topol_abspaths):
            _rtf_res_dict = _rtf_res_dicts[_topol_abspath]
            if _top_seg_map[_itop,_iseg] == False and matched_segs[_iseg] == 0:
                # YL 14/09/2021: update missing residues in the current segment
                if _isAliasedResInRTF(_seg, _segFrag, _rtf_res_dict, _dict_seg_undefined_res, residue_aliases=_residue_aliases, mute=True):
                    _top_seg_map[_itop,_iseg] = True
                    _seg_top_dict.update({_seg:_topol_abspaths[_itop]})
                    continue

#    for _seg, _residues_undefined in _dict_seg_undefined_res.items():
#        for _res in _residues_undefined:
#            if _res in _residue_aliases:
#                _dict_seg_undefined_res.update({_seg:append(_residues_undefined, _residue_aliases.get(_res, array([], dtype=_residues_undefined.dtype)))})

    # update number of matches
    matched_segs = _top_seg_map.sum(axis=0)
    matched_tops = _top_seg_map.sum(axis=1)

    # if built-in prior to standard
    if _override_standard:
        # step 3: completment missing RTF using the built-in library if any segment has not been matched with a RTF
        if 0 in matched_segs or _dict_seg_undefined_res:
            print("\n Searching in ChemShell built-in RTF library...")
            stdout.flush()
            _seg_top_dict_prescreened = {}
            _top_seg_map = _complementRTF(ffObj, _split_frags, ffObj.builtin_rtfs, _seg_top_dict, _seg_top_dict_prescreened, _topol_abspaths, _rtf_res_dicts, _top_seg_map, _additional_prm, _residue_patches, _residue_aliases, _dict_seg_undefined_res, _use_newer_version=_use_newer_version, _prescreen=True)
            _top_seg_map = _complementRTF(ffObj, _split_frags, ffObj.builtin_rtfs, _seg_top_dict, _seg_top_dict_prescreened, _topol_abspaths, _rtf_res_dicts, _top_seg_map, _additional_prm, _residue_patches, _residue_aliases, _dict_seg_undefined_res, _use_newer_version=_use_newer_version)
            matched_segs = _top_seg_map.sum(axis=0)
            matched_tops = _top_seg_map.sum(axis=1)
    
        # step 4: completment missing RTF using the standard library if any segment has not been matched with a RTF
        if 0 in matched_segs or _dict_seg_undefined_res:
            print("\n Searching in CHARMM standard RTF library...")
            stdout.flush()
            _seg_top_dict_prescreened = {}
            _top_seg_map = _complementRTF(ffObj, _split_frags, ffObj.default_rtfs, _seg_top_dict, _seg_top_dict_prescreened, _topol_abspaths, _rtf_res_dicts, _top_seg_map, _additional_prm, _residue_patches, _residue_aliases, _dict_seg_undefined_res, _use_newer_version=_use_newer_version, _prescreen=True)
            _top_seg_map = _complementRTF(ffObj, _split_frags, ffObj.default_rtfs, _seg_top_dict, _seg_top_dict_prescreened, _topol_abspaths, _rtf_res_dicts, _top_seg_map, _additional_prm, _residue_patches, _residue_aliases, _dict_seg_undefined_res, _use_newer_version=_use_newer_version)
            matched_segs = _top_seg_map.sum(axis=0)
            matched_tops = _top_seg_map.sum(axis=1)

    # default: if standard prior to built-in
    else:
        # step 3: completment missing RTF using the standard library if any segment has not been matched with a RTF
        if 0 in matched_segs or _dict_seg_undefined_res:
            print("\n Searching in CHARMM standard RTF library...")
            stdout.flush()
            _seg_top_dict_prescreened = {}
            _top_seg_map = _complementRTF(ffObj, _split_frags, ffObj.default_rtfs, _seg_top_dict, _seg_top_dict_prescreened, _topol_abspaths, _rtf_res_dicts, _top_seg_map, _additional_prm, _residue_patches, _residue_aliases, _dict_seg_undefined_res, _use_newer_version=_use_newer_version, _prescreen=True)
            _top_seg_map = _complementRTF(ffObj, _split_frags, ffObj.default_rtfs, _seg_top_dict, _seg_top_dict_prescreened, _topol_abspaths, _rtf_res_dicts, _top_seg_map, _additional_prm, _residue_patches, _residue_aliases, _dict_seg_undefined_res, _use_newer_version=_use_newer_version)
            matched_segs = _top_seg_map.sum(axis=0)
            matched_tops = _top_seg_map.sum(axis=1)

        # step 4: completment missing RTF using the built-in library if any segment has not been matched with a RTF
        if 0 in matched_segs or _dict_seg_undefined_res:
            print("\n Searching in ChemShell built-in RTF library...")
            stdout.flush()
            _seg_top_dict_prescreened = {}
            _top_seg_map = _complementRTF(ffObj, _split_frags, ffObj.builtin_rtfs, _seg_top_dict, _seg_top_dict_prescreened, _topol_abspaths, _rtf_res_dicts, _top_seg_map, _additional_prm, _residue_patches, _residue_aliases, _dict_seg_undefined_res, _use_newer_version=_use_newer_version, _prescreen=True)
            _top_seg_map = _complementRTF(ffObj, _split_frags, ffObj.builtin_rtfs, _seg_top_dict, _seg_top_dict_prescreened, _topol_abspaths, _rtf_res_dicts, _top_seg_map, _additional_prm, _residue_patches, _residue_aliases, _dict_seg_undefined_res, _use_newer_version=_use_newer_version)
            matched_segs = _top_seg_map.sum(axis=0)
            matched_tops = _top_seg_map.sum(axis=1)

    # print
    print("\n Automatic mapping between segments and topology files:")
    filenames_topol = char.rpartition(_topol_abspaths, sep)[:,-1]
    max_strlen_topol = max(char.str_len(filenames_topol))
    format_topol = ' '*6 + ''.join(' {:^'+str(max_strlen_topol)+'}')*len(_topol_abspaths)
    print(' '+'-'*(6+(max_strlen_topol+1)*len(_topol_abspaths)))
    print(' '+format_topol.format(*filenames_topol))
    print(' '+'-'*(6+(max_strlen_topol+1)*len(_topol_abspaths)))
    format_seg = '{:<6}' + ''.join(' {:^'+str(max_strlen_topol)+'}')*len(_topol_abspaths)
    for _iseg, _seg in enumerate(_split_frags.dtype.names):
        mappings = [ '+' if matched else '' for matched in _top_seg_map.T[_iseg] ]
        print(' '+format_seg.format(_seg, *mappings))
    print(' '+'-'*(6+(max_strlen_topol+1)*len(_topol_abspaths)))

#    # step 5: there might be more than one files associated with a same segment
#    #         and we need to remove the unnecessary ones
#    print("### _topol_abspaths =", _topol_abspaths)
#    for _iseg, _seg in enumerate(_split_frags.dtype.names):
#        for _itop, _topol_abspath in enumerate(_topol_abspaths):
#            # exclusively matched file
#            if _top_seg_map[_itop,_iseg] == True and (matched_segs[_iseg] == 1 or matched_tops[_itop] == 1):
#                _seg_top_dict.update({_seg:_topol_abspaths[_itop]})
#            # multiple files
#            elif _top_seg_map[_itop,_iseg] == True and (matched_segs[_iseg] > 1):
#                # YL 22/07/2021: if a segment is associated with multiple RTF files the user-defined file should be taken first 
#                if not _topol_abspath in ffObj.default_rtfs and not _topol_abspath in ffObj.builtin_rtfs:
#                    _seg_top_dict.update({_seg:[_topol_abspaths[_itop]]})
##                else:
##                    _seg_top_dict.update({_seg:_topol_abspaths[_itop]]})

    if 0 in matched_segs:
        format.printwrap('>>> ChemShell/PSFGEN ERROR: There is at least one segment without matching topology:')
        format.printwrap(' '.join(array(_split_frags.dtype.names)[matched_segs==0]))
        stdout.flush()
        exit(999)

    return _seg_top_dict


def __getAbsPath(_topology, _dir):
    '''Check if a given file exists'''

    from os.path import abspath, isdir, isfile, join, relpath

    _topol_abspath = ''

    if isfile(abspath(_topology)):
        _topol_abspath = abspath(_topology)
    elif isfile(join(_dir, _topology)):
        _topol_abspath = join(_dir, _topology)
    else:
        print("\n WARNING: Specified file {} cannot be located as\n {}".format(_topology, abspath(_topology)))
        if _dir:
            print(" nor as\n {}".format(join(_dir, _topology)))
        print(" Skipping...")

    return _topol_abspath

    
def getAbsPathsList(_topologies, _dir):
    '''Get a list of absolute paths to files'''

    _topol_abspaths = []

    if type(_topologies) is str:
        _topologies = [_topologies]

    for _topol in _topologies:
        if not _topol:
            continue
        _topol_abspath = __getAbsPath(_topol, _dir)
        if _topol_abspath:
            _topol_abspaths.append(_topol_abspath)
        # not found
        else:
            continue

    return _topol_abspaths

def _complementRTF(_ffObj,
                   __split_frags,
                   __rtfs,
                   __seg_top_dict,
                   __seg_top_dict_prescreened, # YL 05/12/2021: added for pre-screening, because there may be more than one RTFs meeting the criteria
                   __topol_abspaths,
                   __rtf_res_dicts,
                   __top_seg_map,
                   __additional_prm,
                   __residue_patches,
                   __residue_aliases,
                   __dict_seg_undefined_res,
                   _use_newer_version=True,
                   _prescreen=False):
    ''''''

    from re      import sub
    from numpy   import all, array, char, isin, unique, vstack, zeros
    from os.path import basename, join
    from sys     import maxsize

    # we only want the residue names as some RTF files are really big and slow-parsing
    extra_rtf_res_dict = rtf.getResidueDicts(__rtfs, residues_only=True)

    # loop over segments
    for _iseg, _seg in enumerate(__split_frags.dtype.names):

        # YL 13/09/2021: only search for the still missing residues
        if _seg in __dict_seg_undefined_res:
            _segFrag = __split_frags[_seg][0]
            _residues = char.strip(char.decode(unique(_segFrag.residues)))
            # YL 21/09/2021: construct a list of residues with aliased names whereever possible
            _residues_aliased = [ __residue_aliases.get(i, i) for i in _residues ]
            # YL 04/12/2021: patches must be dealt with separately
            _residues_patched = [ v for k, v in __residue_patches.items() if k in _residues ]

            # YL 05/12/2021: prescreen all files
            if _prescreen:
                rtfs_candidate = __rtfs
            # YL 05/12/2021: only prescreened files
            else:
                rtfs_candidate = __seg_top_dict_prescreened.get(_seg, [])
                # YL 05/12/2021: if there are multiple candidates
                if len(rtfs_candidate) > 1:
                    versions = []
                    print_warning = False
                    for rtf_candidate in rtfs_candidate:
                        # get a version number
                        try:
                            version = int(sub('[^0-9]','',basename(rtf_candidate)))
                        # we can't decide then
                        except:
                            version = maxsize
                            print_warning = True
                        if _use_newer_version:
                            if all(version>array(versions)):
                                rtfs_candidate = [rtf_candidate]
                                versions += [version]
                        else:
                            if all(version<array(versions)):
                                rtfs_candidate = [rtf_candidate]
                                versions += [version]
                    if print_warning:
                        print('\n WARNING: Found multiple qualified RTFs but using the first')
                        print('', '\n '.join(rtfs_candidate))
                    else:
                        print(f'\n Found multiple qualified RTFs and version {versions[-1]} is used as requested')
                        print('', rtfs_candidate[0])

            # loop over topology files
            for _itop, _topol_abspath in enumerate(rtfs_candidate):
                _rtf_res_dict = extra_rtf_res_dict[_topol_abspath]
                # any matching residue: we assume there is no duplicate residues in the built-in RTF files
                _masks_residues_found_in_rtf = isin(_residues, list(_rtf_res_dict.keys()))
                # YL 04/12/2021: all aliased and patched residues must be found
                _masks_residues_aliased_found_in_rtf = isin(_residues_aliased, list(_rtf_res_dict.keys()))
                _masks_residues_patched_found_in_rtf = isin(_residues_patched, list(_rtf_res_dict.keys()))
                if not all(_masks_residues_aliased_found_in_rtf) or not all(_masks_residues_patched_found_in_rtf):
                    continue
                _masks_residues_found_in_rtf += _masks_residues_aliased_found_in_rtf
                _residues_found_in_rtf = _residues[_masks_residues_found_in_rtf]
                _masks_residues_found_in_wanted = isin(_residues_found_in_rtf, __dict_seg_undefined_res[_seg])
                _nresidues_found_in_wanted = _masks_residues_found_in_wanted.sum()
                if _masks_residues_found_in_rtf.sum() and _nresidues_found_in_wanted:
                    if not _prescreen:
                        if _nresidues_found_in_wanted == 1:
                            format.printwrap("RTF automatically suggested for segment \"{}\" containing residue:\n {}".format(_seg, _residues_found_in_rtf[0]),
                                             indent='', initial_indent=' ', subsequent_indent=' ')
                        elif _nresidues_found_in_wanted > 1:
                            format.printwrap("RTF automatically suggested for segment \"{}\" containing resiudes:\n {}".format(_seg, ' '.join(_residues_found_in_rtf)),
                                             indent='', initial_indent=' ', subsequent_indent=' ')
                        print(' '+_topol_abspath)
                        stdout.flush()
                    else:
                        # always append to a list
                        __seg_top_dict_prescreened.update({_seg:__seg_top_dict_prescreened.get(_seg, [])+[_topol_abspath]})

                    # YL 16/09/2021: remove the found residues from __dict_seg_undefined_res
                    _masks_wanted_found_in_residues = isin(__dict_seg_undefined_res[_seg], _residues_found_in_rtf)
                    _wanted_remaining = __dict_seg_undefined_res[_seg][~_masks_wanted_found_in_residues]
                    if not _prescreen:
                        if len(_wanted_remaining):
                            __dict_seg_undefined_res.update({_seg:_wanted_remaining})
                        else:
                            __dict_seg_undefined_res.pop(_seg)

                    # update __seg_top_dict
                    if not _prescreen:
                        if _seg in __seg_top_dict:
                            if type(__seg_top_dict[_seg]) is str and _topol_abspath != __seg_top_dict[_seg]:
                                __seg_top_dict.update({_seg:[_topol_abspath,__seg_top_dict[_seg]]})
                            elif type(__seg_top_dict[_seg]) is list and _topol_abspath not in __seg_top_dict[_seg]:
                                __seg_top_dict.update({_seg:[_topol_abspath]+__seg_top_dict[_seg]})
                        else:
                            __seg_top_dict.update({_seg:_topol_abspath})

                    if not _topol_abspath in __topol_abspaths and not _prescreen:
                        __topol_abspaths.append(_topol_abspath)
                        # add the associated PRM file
                        filename = basename(_topol_abspath)
                        if filename in _ffObj.rtfs2prms:
                            __additional_prm += [_ffObj.prms2abspaths[_ffObj.rtfs2prms[filename]]]
                        # YL 17/09/2021: attach an extra row to _top_seg_map iff _topol_abspath not previously added
                        _extra_top_seg_map = zeros((1,len(__split_frags.dtype)), dtype=bool)
                        _extra_top_seg_map[0,_iseg] = True
                        __top_seg_map = vstack((__top_seg_map, _extra_top_seg_map))

                    # update the record
                    if not _prescreen:
                        position = __topol_abspaths.index(_topol_abspath)
                        __top_seg_map[position,_iseg] = True

                        # here we need the real data
                        dbuff = rtf.getResidueDicts([_topol_abspath])
                        _rtf_res_dict = dbuff[_topol_abspath]

                        # update the RTF-residues dict
                        __rtf_res_dicts.update({_topol_abspath:_rtf_res_dict})

                        # YL TODO 03/12/2021: match only one file in the non-prescreening mode (see the above version selecting method)
                        break

    # has to be returned as it cannot be changed in place as a parameter being passed
    return __top_seg_map





