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

# this is an interface to C. Yong's DL_FIELD
# see: http://www.ccp5.ac.uk/DL_FIELD
#      DOI: 10.1021/acs.jcim.6b00323

# you.lu@stfc.ac.uk
# March 2020
# YL revised Feb 2021

# ChemShell data structure of forcefield

# { 'scheme name': { 'residues': [...],
#                    'residue name 1': { 'atom': ndarray(...),
#                                        ....
#                                      },
#                    'residue name 2': {...},
#                    'bonds'    : ndarray(...),
#                    'angles'   : ndarray(...),
#                    'dihedrals': ndarray(...),
#                    ....
#                   }
# }
from sys import stdout
import chemsh
from ..driver         import _Driver
from ...dl_py2f.utils import dictutils, fileutils, iterutils, miscutils, modutils, nputils, strutils

class DL_FIELD(_Driver):

    # YL NB: keep this order!
    _kwargs = {
                'dl_poly'         : True,
                'ff'              :'CHARMM',
                'unit'            :'kcal/mol',
                'conversion'      :'normal',
                'bond'            :'harmonic',
                'angle'           :'harmonic',
                'udff'            :'',
                'verbosity'       : True,
                'input'           :'_dl_field.inp',
                'output'          :'_dl_field.out.pdb',
                'optimise'        : True,
                'override_builtin': False,
                'solution'        : False,
                'solution_den'    : 1.05,
                'solution_unit'   :'g/cm^3',
                'solution_cutoff' : 1.9,
                'atom_display'    :'standard',
                'vdw_format'      :'12-6',
                'eps_org'         :'arithmatic',
                'sig_org'         :'arithmatic',
                'eps_inorg'       :'geometric',
                'sig_inorg'       :'arithmatic',
                'eps_diff'        :'arithmatic',
                'sig_diff'        :'arithmatic',
                'extra_info'      : False,
                'freeze'          : False,
                'tether'          : False,
                'constrain'       : False,
                'rigid_body'      : False,
                'pbc'             : 0,           # YL 24/11/2023: this was fixed on the branch for the inorganic Solvation workflow but let's make a hotfix here
                'scale_coul'      : 0,
                'scale_vdw'       : 0,
                'velocity'        : False,
                'temperature'     : 300.0,
                'recentre'        : False,
                'solvent'         : None,
                'distance'        : 1.6,
                'counterions'     : False,
                'min_distance'    : 10.0,
                'mm_energy'       : False,
                'cutoff_coulomb'  : 100.0,
                'cutoff_vdw'      : 10.0,
                'path_lib'        :'',
                'path_sol'        :'',
                'path_out'        :'',
                'filename'        :'_dl_field.pdb',
                'exclude'         :[],
                'qmatoms'         :[],
                'solvents'        :[],
                'remove_unused_definitions': True,
                'original_tip3p'  : False,
                'directory'       :'_dl_field',
                # YL 03/11/2021: we default to run DL_FIELD only at the first cycle of optimisation/MD while keep this optional
                'field_by_chemsh' : True,
                'run_only_once'   : True,
                'patch_target'    :{},
              }

    _internals = {
#                   '_udff_data': dict(),
                   '_reskeys'  : dict(), # foreign residue keys
                   '_dlfkeys'  : dict(), # built-in DL_FIELD residue keys
                   '_dlftypes' : dict(),
                   '_rerun'    : True,
                 }

    _priorities = 'ff', 'patch_target'


    @miscutils.memusage
    def runLib(self, **kwargs):
    
        from os         import getcwd, path
        from sys        import stdout
        from ctypes     import byref, c_int
        from numpy      import all, array, array2string, char, isin, ndarray, unique
        from shutil     import copyfile
        from ...objects import fragbio
        from ...dl_py2f import py2f
        from ...io      import file, format, banner
        from ...data    import forcefield, references
        from .          import dl_f_path

        banner.printModuleBanner('ChemShell/DL_FIELD Interface')
#        print(' Please cite:')
#        format.printwrap(references.references['dl_field'], prepend_blanc_line=False)

        # archive existing dir
        fileutils.archiveDir(self.directory)

        # YL TODO: .PDB files are also used for materials science
        ### PDB ###

        if isinstance(self.frag, fragbio.FragmentBio) and self._format == 'pdb':

            # back up
            fileutils.archiveFile(self.filename, dir=self.directory)

            # YL 24/11/2023: patch residues
            self.patchFragment()

            # YL 24/03/2021: choose between TIP3P versions:
            if self.original_tip3p:
                print()
                format.printwrap('ChemShell/DL_FIELD: Using the original TIP3P model by renaming water residue names to \"TP3O\"')
                indices_all_water = self.frag.select(residues=self._ffObj.waters)
                self.frag.residues[indices_all_water] = 'TP3O'

            # assign forcefield schemes: typically received from ruturned dictionary of charmm2udff()
            # TODO: in the future we could support such as {'protein':'CHARMM22_pot', 'HEMC':'user_defined.rtf'}
            if type(self.ff) is dict:
                self.writeFFsToUDFF([self.ff])

            # str: UDFF or scheme name
            elif issubclass(type(self.ff), str):
                # this implicitly updates self._reskeys
                self.getDataFromFF(self.ff)

            # multiple FFs (not necessarily multiple schemes)
            elif type(self.ff) is not str and hasattr(self.ff, '__iter__'):
                self.writeFFsToUDFF(self.ff)

            # this judge is based on self._reskeys
            all_residues_found = self.all_residues_found

            # if any residue missing
            # METHOD ONE
            # YL 13/02/2021: we screen the built-in CHARMM library to find missing residues
            #                and write out a UDFF to be used together with DL_FIELD's built-in library
            if not all_residues_found:

                # note: relies on self._reskeys
                residues_not_in_reskeys = self.all_resnames[~self.masks_residues_found]
                print("\n Residues (i.e. DL_FIELD molecule types) not found:\n {}".format(', '.join(residues_not_in_reskeys)))

                # find RTF/PRM files containing the missing residues from the ChemShell library
                rtfs, prms, nfounds = self._associatedFFLib.getFilesOfResidues(residues_not_in_reskeys)

                # form a list of the FFs if all residues found
                if len(residues_not_in_reskeys) >= nfounds:
                    all_residues_found = True

                    # NB: dict type also has __iter__
                    # dict of forcefield data (typically by charmm2udff)
                    if type(self.ff) is dict:
                        # YL 20/02/2021: we don't yet deal with multiple schemes
                        # get the data out of the RTF/PRM files (and write out an intermediate UDFF file)
                        dbuff = self.charmm2udff(rtf=rtfs, prm=prms, output='_charmm.udff', version=self._scheme_name, charmm=self._associatedFFLib)

                        # NB: if override_builtin: user-specified data dict self.ff overrides ChemShell library [dbuff]. see: self.writeFFsToUDFF()
                        if self.override_builtin:
                            list_ffs = [dbuff] + [self.ff]
                        else:
                            list_ffs = [self.ff] + [dbuff]

                    # list
                    elif type(self.ff) is not str and hasattr(self.ff, '__iter__'):
                        # YL 20/02/2021: we don't yet deal with multiple schemes
                        # get the data out of the RTF/PRM files (and write out an intermediate UDFF file)
                        dbuff = self.charmm2udff(rtf=rtfs, prm=prms, output='_charmm.udff', version=self._scheme_name, charmm=self._associatedFFLib)

                        # NB: if override_builtin: ChemShell library [dbuff] overrides DL_FIELD built-in library self.ff. see: self.writeFFsToUDFF()
                        if self.override_builtin:
                            list_ffs = list(self.ff) + [dbuff]
                        else:
                            list_ffs = [dbuff] + list(self.ff)

                    # str: UDFF or scheme name
                    elif issubclass(type(self.ff), str):
                        # get the data out of the RTF/PRM files (and write out an intermediate UDFF file)
                        # override_default=True is default
                        dbuff = self.charmm2udff(rtf=rtfs, prm=prms, output='_charmm.udff', version=self.ff, charmm=self._associatedFFLib, override_default=self.override_builtin)
                        list_ffs = [self.ff, dbuff]


                    # merge with DL_FIELD scheme and write out a final UDFF
                    self.writeFFsToUDFF(list_ffs)

            # write out a PDB as the DL_FIELD input
            if all_residues_found:
                # uppercase=False is default: we don't write capital elements to be consistent with the udff.write()
                self.frag.save(self.filename, self._format, uppercase=False)

            # METHOD TWO
            # YL 13/02/2021: if still no matching we have to split the original fragment to
            #                use the multiple-scheme method of DL_FIELD by writing out a new
            #                PDB with each FF scheme name followed by segment
            # YL 15/01/2021: the split() function sorts the segment! we hence 
            else:

                # YL TODO: this method is a plan for the future. we raise an error for now
                residues_not_in_reskeys = self.all_resnames[~self.masks_residues_found]
                print("\n ChemShell/DL_FIELD ERROR: there are still undefined molecule keys:\n ", residues_not_in_reskeys)
                exit(999)
                self.splitFragment()


        ### END PDB ###
        ### non-template-based forcefields (e.g., .XYZ) ###

        # write out a new structure file to be read in by DL_FIELD
        else:
            self.frag.save(self.filename)

        # QM/MM
        if kwargs.get('qmatoms', []):
            self.qmatoms = kwargs.get('qmatoms', [])
        # see: interfaces.QMMM.qmmm_pairs_by_junctions
        if kwargs.get('exclude', []):
            self.exclude = kwargs.get('exclude', [])

        # YL TODO: self.path_lib/path_sol/path_out should overload dl_f_path
        if len(self.qmatoms):
            print("\n ChemShell/DL_FIELD QM/MM treatments")
            print("\n QM atoms:\n {}".format(array2string(array(self.qmatoms), \
                                                          formatter={'int':lambda x:'{:>7d}'.format(x)}, \
                                                          separator=', ', \
                                                          max_line_width=72)[1:-1]))
            print("\n QM-MM interactions to exclude about QM/MM junction (QM-JQ-MM):")
            for pair in self.exclude:
                print(" {:>7} - {:>7}".format(pair[0], pair[1]))

        # transfer data with PY2F and get a handle
        py2f_handle = py2f(self, byref=True)

        # write out a 'dl_f_path' file needed by DL_FIELD
        dl_f_path.write_path_file(dl_field_path=self.dl_field_path, input=self.input, directory=self.directory)

        # if field_by_chemsh (default: True) we deal with the QM/MM special treatments with ChemShell instead of DL_FIELD
        if self.field_by_chemsh:
            self.qmatoms = []
            self.exclude = []

        # run
        ierror = self.lib.dl_field_main(py2f_handle,
                                        byref(c_int(self.frag.natoms)),
                                        byref(c_int(len(self.qmatoms))),
                                        byref(c_int(array(self.exclude).shape[0])))
        stdout.flush()

        del py2f_handle

#        print("\n Elapsed time of procedure chemsh.interfaces.dl_field.runLib: {} sec\n".format(time()-t0))
        stdout.flush()

        for output in [ 'CONFIG', 'FIELD' ]:
            # back up
            fileutils.archiveFile(output, dir=self.directory)
#            copyfile(path.join(self.dl_field_path, 'output', 'dl_poly.'+output), output)
            # YL 07/07/2022: directly generated to the working directory
            copyfile(path.join(getcwd(), '_dl_field', 'dl_poly.'+output), output)

        # archive old log file and rename the current log file
        fileutils.rename('dl_field.output', '_dl_field.out')
        fileutils.archiveFile([ '_dl_field.inp',  \
                                '_dl_field.out',  \
                                '_dl_field.pdb',  \
                                '_dl_field.udff', \
                                'dl_f_path' ],
                              dir=self.directory, suffix='')

        # save content to fftext
        with open('FIELD', 'r') as fp:
            self.fftext = fp.read()
    
        ierror = self.lib.dl_field_finalise()

        print()


    @miscutils.memusage
    def getInpBuff(self, **kwargs):
        ''''''

        from numpy      import any, array, char, isin
        from .          import schemes
        from ...io      import format, file
        from ...objects import fragbio
        from sys        import stdout

        # check if DL_FIELD is available before proceeding
        if not self.path and not modutils.getLinkedLib(self.__class__.__name__, chemsh.CHEMSH_LIB_DIR) and not self._syscmd:
            format.printwrap('>>> ERROR: DL_FIELD not compiled with ChemShell nor a system path to DL_FIELD executable has been provided. Exiting...',
                             indent='', initial_indent=' ', subsequent_indent=' '*12)
            stdout.flush()
            exit(999)

        # switch off PBC for biological fragments: getInpBuff() is called before runLib()
        if isinstance(self.frag, fragbio.FragmentBio) and self._format == 'pdb':
            self.pbc = None

        # get correct file format if not overridden by user
        if self.filename == self.__class__._kwargs['filename'] and issubclass(type(self.ff), str):
            self.filename = '_dl_field.' + schemes.schemes.get(self.ff.lower(), 'pdb')

        fragObj = self.frag
    
        # YL NB: some lines are only read in with 80 characters
    
        DWIDTH = SWIDTH = 12
    
        ISWITCH = {True:1,False:0}

        # title line
        strbuff = '{:80}\n'.format('Py-ChemShell interface to DL_FIELD v4.8')
    
        # line 1
        strbuff += '{}'.format({True:1}.get(self.dl_poly, 0))
        strbuff += ' * Construct DL_POLY output files\n'
    
        # line 2
        strbuff += '1'
        strbuff += ' * For ChemShell\n'

        strbuff += '{:{sw}}'.format(self._scheme_name, sw=SWIDTH)
        strbuff += ' * Type of force field require\n'
    
        # line 2
        strbuff += '{:{sw}}'.format(self.unit, sw=SWIDTH)
        strbuff += ' * Energy unit: kcal/mol, kJ/mol, eV or default\n'
    
        # line 3
        strbuff += '{:{sw}}'.format(self.conversion, sw=SWIDTH)
        strbuff += ' * Conversion criteria (strict, normal, loose)\n'
    
        strbuff += '{:<{dw}d}'.format({'default':0,'harmonic':1,'morse':2}[self.bond.lower()], dw=DWIDTH)
        strbuff += ' * Bond type (0=default, 1=harmonic , 2=Morse)\n'
    
        strbuff += '{:<{dw}d}'.format({'default':0,'harmonic':1,'cos':2}[self.angle.lower()], dw=DWIDTH)
        strbuff += ' * Angle type (0=default, 1=harmonic, 2=harmonic cos)\n'
    
        # UDFF
        # YL 14/02/2021: changed self.udff to only str type
        strbuff += '{:{sw}}'.format({'':'none'}.get(self.udff, self.udff), sw=SWIDTH)
        strbuff += ' * Include user-defined information. Put \"none\" or a .udff filename\n'
    
        strbuff += '{:<{dw}d}'.format(ISWITCH[self.verbosity], dw=DWIDTH)
        strbuff += ' * Verbosity mode: 1 = on, 0 = off\n'
    
        strbuff += '{:{sw}}'.format(self.filename, sw=SWIDTH)
        strbuff += ' * Filename for the user\'s atomic configuration.\n'

        strbuff += '{:{sw}}'.format({None:'none'}.get(self.output, self.output), sw=SWIDTH)
        strbuff += ' * Output file in PDB. Put \"none\" if not needed.\n'

        strbuff += '{:<2d} {:<6.3} {:8} {:<6.3}'.format(ISWITCH[self.solution],
                                                        self.solution_den,
                                                        self.solution_unit,
                                                        self.solution_cutoff)
        strbuff += ' * Solution Maker: on/off, density, unit, cutoff)\n'

        strbuff += '{:<{dw}d}'.format(ISWITCH[self.optimise], dw=DWIDTH)
        strbuff += ' * Optimise FIELD output size, if possible? 1=Yes 2=No\n'

        strbuff += '{:<{dw}d}'.format({'dl_field':1,'standard':2}[self.atom_display.lower()], dw=DWIDTH)
        strbuff += ' * Atom display: 1 = DL_FIELD format. 2 = Standard format\n'

        strbuff += '{:<{dw}d}'.format({'12-6':1,'lj':2,'l-j':2}[self.vdw_format.lower()], dw=DWIDTH)
        strbuff += ' * Vdw display format: 1 = 12-6 format 2 = LJ format\n'

        strbuff += '{:{sw}}'.format({'geometric':1,'arithmatic':'default'}[self.eps_org.lower()], sw=SWIDTH)
        strbuff += ' * Epsilon mixing rule (organic FF only) : 1 = geometric 2 = arithmatic or default\n'

        strbuff += '{:{sw}}'.format({'geometric':1,'arithmatic':'default'}[self.sig_org.lower()], sw=SWIDTH)
        strbuff += ' * Sigma mixing rule (organic FF only) : 1 = geometric 2 = arithmatic or default\n'

        strbuff += '{:<{dw}d}'.format({'geometric':1,'arithmatic':2}[self.eps_inorg.lower()], dw=DWIDTH)
        strbuff += ' * Epsilon mixing rule (inorganic FF only) : 1 = geometric 2 = arithmatic\n'

        strbuff += '{:<{dw}d}'.format({'geometric':1,'arithmatic':2}[self.sig_inorg.lower()], dw=DWIDTH)
        strbuff += ' * Sigma mixing rule (inorganic FF only) : 1 = geometric 2 = arithmatic\n'

        strbuff += '{:<{sw}}'.format({'geometric':1,'arithmatic':2}[self.eps_diff.lower()], sw=SWIDTH)
        strbuff += ' * Epsilon mixing rule (between different FF) : 1 = geometric 2 = arithmatic or default\n'

        strbuff += '{:<{sw}}'.format({'geometric':1,'arithmatic':2}[self.sig_diff.lower()], sw=SWIDTH)
        strbuff += ' * Sigma mixing rule (between different FF) : 1 = geometric 2 = arithmatic or default\n'

        strbuff += '{:<{dw}d}'.format(ISWITCH[self.extra_info], dw=DWIDTH)
        strbuff += ' * Display additional info. for protein 1=Yes 0=No\n'

        strbuff += '{:<{dw}d}'.format(ISWITCH[self.freeze], dw=DWIDTH)
        strbuff += ' * Freeze atoms? 1=Yes 0=No\n'

        strbuff += '{:<{dw}d}'.format(ISWITCH[self.tether], dw=DWIDTH)
        strbuff += ' * Tether atoms? 1=Yes 0=No\n'

        if self.original_tip3p:
            print('\n ChemShell/DL_FIELD: rigid water is to be used for the original TIP3P model\n')
            self.constrain = True
        strbuff += '{:<{dw}d}'.format(ISWITCH[self.constrain], dw=DWIDTH)
        strbuff += ' * Constrain bonds? 1=Yes 0=No\n'

        strbuff += '{:<{dw}d}'.format(ISWITCH[self.rigid_body], dw=DWIDTH)
        strbuff += ' * Apply rigid body? 1=Yes 0=No\n'


        strbuff += '{:<{dw}}'.format({'no'                :0,
                                       0                  :0,
                                       None               :0,
                                       False              :0,
                                      'cubic'             :1,
                                      'orthorhombic'      :2,
                                      'parallelepiped'    :3,
                                      'x-y_parallelepiped':6,
                                      'slab'              :6,
                                      '2D'                :6}.get(self.pbc, self.pbc), dw=DWIDTH)
        strbuff += ' * Periodic condition ? 0=No, other number = type of box\n'
    
        strbuff += '{0[0]:<14.9} {0[1]:<14.9} {0[2]:<14.9}'.format((fragObj.cell.vectors[0][0], fragObj.cell.vectors[0][1], fragObj.cell.vectors[0][2]))
        strbuff += ' * Cell vector a (x, y, z)\n'
    
        strbuff += '{0[0]:<14.9} {0[1]:<14.9} {0[2]:<14.9}'.format((fragObj.cell.vectors[1][0], fragObj.cell.vectors[1][1], fragObj.cell.vectors[1][2]))
        strbuff += ' * Cell vector b (x, y, z)\n'
    
        strbuff += '{0[0]:<14.9} {0[1]:<14.9} {0[2]:<14.9}'.format((fragObj.cell.vectors[2][0], fragObj.cell.vectors[2][1], fragObj.cell.vectors[2][2]))
        strbuff += ' * Cell vector c (x, y, z)\n'
    
        strbuff += '{:<{sw}}'.format({0:'default',1:1,2:2,3:3,4:4}[self.scale_coul], sw=SWIDTH)
        strbuff += ' * 1-4 scaling for coulombic (put default or x for scaling=x)\n'
    
        strbuff += '{:<{sw}}'.format({0:'default',1:1,2:2,3:3,4:4}[self.scale_vdw], sw=SWIDTH)
        strbuff += ' * 1-4 scaling for vdw (put default or x for scaling=x)\n'
    
        strbuff += '{:<2d} '.format(ISWITCH[self.velocity], dw=DWIDTH)
        strbuff += '{:<9.2f}'.format(self.temperature)
        strbuff += ' * Include velocity? 1=Yes, 0=No and scaling temperature.\n'
    
        strbuff += '{:<{dw}d}'.format(ISWITCH[self.recentre], dw=DWIDTH)
        strbuff += ' * Position solute at origin 1=Yes, 0=No\n'
    
        strbuff += '{:5} '.format({None:'none'}.get(self.solvent, self.solvent), sw=SWIDTH)
        strbuff += '{:<6.3}'.format(self.distance)
        strbuff += ' * Solvate model? \"none\" or specify solvent and distance criteria.\n'
    
        strbuff += '{:<2d} '.format(ISWITCH[self.counterions], dw=DWIDTH)
        strbuff += '{:<9.3}'.format(self.min_distance)
        strbuff += ' * Add counter ions? 1=Yes, 0=No, minimum distance from solute\n'
    
        strbuff += '{:<{dw}d}'.format(ISWITCH[self.mm_energy], dw=DWIDTH)
        strbuff += ' * MM energy calculation. 1=Yes, 0=No\n'
    
        strbuff += '{:<12.3}'.format(self.cutoff_coulomb)
        strbuff += ' * Cut off for electrostatic energy calculation (angstrom)\n'
    
        strbuff += '{:<12.3}'.format(self.cutoff_vdw)
        strbuff += ' * Cut off for vdw energy calculation (angstrom)\n'

        strbuff += '\n\n########################################################\n'
        # YL 24/03/2021: rigid_water must be used for the original TIP3P model
        if self.constrain and self.original_tip3p:
            for segment in self.frag.all_segments:
                residues = char.decode(self.frag.residues[self.frag.select(segments=segment)])
                if any(isin(residues, self._ffObj.waters)):
                    strbuff += f'CONSTRAIN {segment.decode()} rigid_water\n'
        strbuff += '\n\n########################################################\n'
    
        return strbuff


    def parseOutput(self):
        '''Dummy parseOutput'''



    def getFF(self):
        ''''''


    @property
    def qmatoms(self):
        '''List of QM atoms'''

        return self._qmatoms


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

        self._qmatoms = val


    @property
    def exclude(self):
        '''List of pairs to exclude'''

        return self._exclude


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

        self._exclude = val


    # YL: "strip" not in use
    def getDataFromLib(self, ff, receiver=dict(), strip=''):
        '''Return the data defined in a SF or UDFF file'''

        from numpy import array, char, argmax, genfromtxt
        from os    import listdir, path
        from ...io import file

        # typically UDFF
        if path.isfile(ff):
            file.parse(ff, 'dlp_sf', receiver=receiver)
 
        path_lib = path.join(self.dl_field_path, 'lib')

        # all files/dirs
        files = listdir(path_lib)

        # loop over built-in files in DL_FIELD lib
        for f in files:
            fullpath = path.join(path_lib, f)
            # YL 15/01/2021: receiver dict with FF scheme (potential) name as top-level key
            if path.isfile(fullpath):
                # YL 22/11/2023: since DL_FIELD 4.10 the filenames are without a "DLPOLY" prefix
                if ('DLPOLY_'+ff+'.sf' ).lower() == f.lower() or (ff+'.sf' ).lower() == f.lower():
                    file.parse(fullpath, 'dlp_sf' , receiver=receiver)
                if ('DLPOLY_'+ff+'.par').lower() == f.lower() or (ff+'.par').lower() == f.lower():
                    file.parse(fullpath, 'dlp_par', receiver=receiver)

        return receiver


    def getDataFromFF(self, ff):
        '''Return a dict holding the forcefield data according to the given single forcefield name'''

        # load FF file from the DL_FIELD built-in library
        _udff_data = dict()

        # from a UDFF file
        # note: lower=True is default
        # TODO: UDFF parser isn't yet complete
        if fileutils.getExtension(ff, lower=True) == 'udff':
            file.parse(ff, 'udff', receiver=_udff_data)
            # loop over schemes to update self._reskeys
            for scheme, res_dict in _udff_data.items():
                self.updateReskeys(scheme, dbuff)

        # from DL_FIELD library
        else:
            self.getDataFromLib(ff, receiver=_udff_data)
            dbuff = _udff_data[ff]
            # merge and update self._reskeys and self._dlfkeys (built-in keys)
            self.updateReskeys(ff, dbuff)
            self.updateReskeys(ff, dbuff, builtin=True)

            return _udff_data


    def getTypesFromSF(self, ff, strip='', receiver=dict()):
        '''Return the molecular types defined in a SF or UDFF file'''

        from numpy import array, char, argmax, genfromtxt, unique
        from os    import listdir, path
        from ...io import file

        data = self.getDataFromLib(ff, receiver=receiver, strip=strip)

        try:
            return char.strip(data['residues']['key'], chars=strip)
        except:
            return []


    def updateReskeys(self, scheme, res_dict, builtin=False):
        '''Merge and update self._reskeys (or self._dlfkeys/self._dlftypes) for given scheme'''

        from numpy import array, char, unique

        # merge with exsiting keys
        if builtin:
            residue_keys = self._dlfkeys.get(scheme, []) + char.strip(res_dict['residues']['key'], chars='-').tolist()
            residue_types = self._dlftypes.get(scheme, []) + res_dict['residues']['type'].tolist()
        else:
            residue_keys = self._reskeys.get(scheme, []) + char.strip(res_dict['residues']['key'], chars='-').tolist()
        # remove duplicates
        residue_keys = unique(residue_keys).tolist()
        # update
        if builtin:
            self._dlfkeys.update({scheme:residue_keys})
            self._dlftypes.update({scheme:residue_types})
        else:
            self._reskeys.update({scheme:residue_keys})


    @staticmethod
    def _getSchemeNameFromUDFF(udff):
        '''Return the POTENTIAL entry value from a UDFF file'''

        from os    import path 
        from numpy import array
        from ...io import file

        if not path.isfile(udff):
            return udff

        dbuff = {}
        file.parse(udff, 'udff', receiver=dbuff)
        udff_potentials = dbuff.get('potential', array([]))

        if udff_potentials.size == 0:
            print(" >>> ERROR: UDFF file does not contain a POTENTIAL entry")
            exit(112)

        elif udff_potentials.size == 1:
            return udff_potentials[0]

        elif udff_potentials.size > 1:
            #return 'multiple'
            return udff_potentials[0]


    @staticmethod
    def _getSchemeNameFromDict(_ff):
        '''Return the forcefield scheme name from a given data dictionary'''

        if len(_ff.keys()) == 0:
             _scheme = self._kwargs['ff']
             print("\n Input forcefield is a dictionary with no scheme named. Using {} by default...".format(_scheme))
        elif len(_ff.keys()) == 1:
             _scheme = list(_ff.keys())[0]
        elif len(_ff.keys()) > 1:
             _scheme = 'multiple'

        return _scheme


    @staticmethod
    def _getSchemeNameFromStr(_ff):
        '''Return the forcefield scheme name from a given string'''

        # UDFF
        if fileutils.getExtension(_ff) == 'udff':
            return DL_FIELD._getSchemeNameFromUDFF(_ff)
        # plain scheme name
        else:
            return _ff


    @staticmethod
    def _getSchemeNameFromFF(_ff):
        '''Get the forcefield scheme name from given ff which can be of any type'''

        # dict ff: typically returned data dict from charmm2udff()
        if type(_ff) is dict:
            scheme = DL_FIELD._getSchemeNameFromDict(_ff)

        # multiple FFs in list/array/tuple
        elif type(_ff) is not str and hasattr(_ff, '__iter__'):
            # defaults
            scheme = 'multiple'
            # we have to merge the FF schemes into a single UDFF if any in the list is UDFF
            for _f in _ff:
                schemes = []
                # dict
                if type(_f) is dict:
                    schemes += [DL_FIELD._getSchemeNameFromDict(_f)]
                # str (either udff filename or DL_FIELD built-in scheme name)
                elif issubclass(type(_f), str):
                    schemes += [DL_FIELD._getSchemeNameFromStr(_f)]
                schemes = list(set(schemes))
                if len(schemes) == 1:
                    scheme = schemes[0]
                elif len(schemes) == 0:
                    scheme = self._kwargs['_f']
                    print("\n Input forcefield \"{}\" found with no scheme named. Using {} by default...".format(_ff, scheme))

        # str (either udff filename or DL_FIELD built-in scheme name)
        else:
            scheme = DL_FIELD._getSchemeNameFromStr(_ff)

        return scheme


    def writeFFsToUDFF(self, ffs, reset_reskeys=True, reset_builtin=True):
        '''Write out a UDFF file of multiple forcefields (any combination of other UDFF data dict or DL_FIELD scheme name)'''

        from ...io   import file, format
        from ...base import forcefield
        from numpy   import arange, array, char, delete, isin, ndarray, unique

        # start with a new self._reskeys
        if reset_reskeys:
            self._reskeys = dict()
        if reset_builtin:
            self._dlfkeys = dict()
            self._dlftypes = dict()

        _udff_data = dict()
        # dict buffer for non-built-in data
        dbuff_foreign = dict()
        # dict buffer for built-in DL_FIELD data
        dbuff_builtin = dict()

        # loop over FFs (UDFF or DL_FIELD schemes)
        for ff in ffs:

            # if there is a merged UDFF file: keep returned data
            if type(ff) is dict:

                for scheme, dbuff in ff.items():

                    # remove repeated:
                    # if self.override_builtin: built-in dict is in the end
                    # if not self.override_builtin: built-in dict is in the beginning
                    # so that we always delete entries from the dict added lately (a1 in forcefield.removeRepeats(a1, a2))
                    for section, abuff in dbuff.items():
                        if scheme in dbuff_foreign and section in dbuff_foreign[scheme] and type(abuff) is ndarray:
                            ff[scheme][section] = forcefield.removeRepeats(abuff, dbuff_foreign[scheme][section])

                    # merge and update list of residue keys
                    self.updateReskeys(scheme, dbuff)

                # merge dict into dbuff_foreign
                dictutils.mergeDicts(dbuff_foreign, ff)

            # NB: if there is any generated UDFF
            elif issubclass(type(ff), str):

                dbuff = {}

                # this implicitly updates self._reskeys
                dbuff = self.getDataFromFF(ff)

                # merge dict into dbuff_builtin
                dictutils.mergeDicts(dbuff_builtin, dbuff)

        # remove repeated entries
        # respect the built-in DL_FIELD data
        if not self.override_builtin:
            for ff, dbuff in dbuff_foreign.items():
                for section, abuff in dbuff.items():
                    # YL 19/02/2021: it is possible that a UDFF file/dictionary is specified without
                    #                accompanying DL_FIELD built-in scheme name
                    if not ff in dbuff_builtin:
                        dbuff = {}
                        dbuff = self.getDataFromFF(ff)
                        dictutils.mergeDicts(dbuff_builtin, dbuff)
                    if section in dbuff_foreign[ff] and section in dbuff_builtin[ff]:
                        dbuff_foreign[ff][section] = forcefield.removeRepeats(abuff, dbuff_builtin[ff][section])
        # delete the built-in DL_FIELD data which are redefined by supplemented data
        else:
            # YL 23/03/2021: likewise, it is possible that a UDFF file/dictionary is specified without
            #                accompanying DL_FIELD built-in scheme name
            for ff, dbuff in dbuff_foreign.items():
                for section, abuff in dbuff.items():
                    if not ff in dbuff_builtin:
                        dbuff = {}
                        dbuff = self.getDataFromFF(ff)
                        dictutils.mergeDicts(dbuff_builtin, dbuff)
            # now screen the built-in data
            for ff, data in dbuff_builtin.items():
                for section, abuff in dbuff.items():
                    if section in dbuff_foreign[ff] and section in dbuff_builtin[ff]:
                        dbuff_builtin[ff][section] = forcefield.removeRepeats(abuff, dbuff_foreign[ff][section])

        # merge dicts: now _udff_data holds everything
        dictutils.mergeDicts(_udff_data, dbuff_foreign)
        dictutils.mergeDicts(_udff_data, dbuff_builtin)

        # YL 15/02/2021: remove duplicated keys
        for ff, dlftypes in self._dlftypes.items():
            nresidues = _udff_data[ff]['residues'].size
            # keep built-in DL_FIELD entries
            if not self.override_builtin:
                masks_types_not_in_dlf = isin(_udff_data[ff]['residues']['type'], dlftypes, invert=True)
                masks_keys_in_dlf      = isin(_udff_data[ff]['residues']['key'] , self._dlfkeys[ff])
                # remove non-built-in entries that have keys occuring in the built-in DL_FIELD library
                indices_to_remove = arange(nresidues)[masks_types_not_in_dlf*masks_keys_in_dlf]
                nindices_to_remove = indices_to_remove.sum()
                if nindices_to_remove:
                    if chemsh.CHEMSH_DEBUG:
                        print('\n'+format.wrap("Removing supplemented {} with type/key already defined in the built-in DL_FIELD library:".format({1:'entry'}.get(nindices_to_remove, 'entries'))))
                        lbuff = []
                        for i in indices_to_remove:
                            lbuff += ['{}/{}'.format(_udff_data[ff]['residues'][i]['type'], _udff_data[ff]['residues'][i]['key'])]
                        print(format.wrap(', '.join(lbuff)))
                    _udff_data[ff]['residues'] = delete(_udff_data[ff]['residues'], indices_to_remove)
            # override built-in DL_FIELD entries
            # YL TODO 24/03/2021: STILL A DRAFT VERSION! a solution to completely overriding the built-in library hasn't been worked out!
            else:
                masks_types_in_dlf = isin(_udff_data[ff]['residues']['type'], dlftypes)
                # remove the '-' for comparing
                _reskeys_udff = char.strip(_udff_data[ff]['residues']['key'], '-')
                # YL TODO: not yet in use
                masks_keys_in_dlf  = isin(_reskeys_udff, self._dlfkeys[ff])
                # keys that occur more than once (assuming self._dlfkeys[ff] contains no repeat)
                masks_repeated_dlf_keys = nputils.nparray.getCountsA1InA2(_reskeys_udff, self._dlfkeys[ff]) > 1
                # remove built-in DL_FIELD entries that have keys also occuring in foreign source
                indices_to_remove = arange(nresidues)[masks_types_in_dlf*masks_repeated_dlf_keys]
                nindices_to_remove = indices_to_remove.sum()

                if nindices_to_remove:
                    if chemsh.CHEMSH_DEBUG:
                        format.printwrap("Removing built-in DL_FIELD {} with type/key also defined in the supplemented forcefield:".format({1:'entry'}.get(nindices_to_remove, 'entries')))
                        lbuff = []
                        for i in indices_to_remove:
                            lbuff += ['{}/{}'.format(_udff_data[ff]['residues'][i]['type'], _udff_data[ff]['residues'][i]['key'])]
                        print(format.wrap(', '.join(lbuff)))
                    _udff_data[ff]['residues'] = delete(_udff_data[ff]['residues'], indices_to_remove)

        # YL 15/02/2021: remove unused residue entries
        if self.remove_unused_definitions:
            for ff, dlftypes in self._dlftypes.items():
                nresidues = _udff_data[ff]['residues'].size
                # keep the -XXX- ones
                masks_residues_not_in_fragment = isin(char.strip(_udff_data[ff]['residues']['key'], '-'), self.all_resnames, invert=True)
# YL 19/02/2021: keep the titratables
#                masks_not_titratables = in1d(char.strip(_udff_data[ff]['residues']['key'], '-'), self._associatedFFLib.protonateds, invert=True)
#                masks_residues_not_in_fragment *= masks_not_titratables
                indices_to_remove = arange(nresidues)[masks_residues_not_in_fragment]
                nindices_to_remove = indices_to_remove.sum()
                if nindices_to_remove:
                    if chemsh.CHEMSH_DEBUG:
                        print('\n'+format.wrap("Removing {} with type/key not used in the current fragment:".format({1:'entry'}.get(nindices_to_remove, 'entries'))))
                        lbuff = []
                        for i in indices_to_remove:
                            lbuff += ['{}/{}'.format(_udff_data[ff]['residues'][i]['type'], _udff_data[ff]['residues'][i]['key'])]
                        print(format.wrap(', '.join(lbuff)))
                    _udff_data[ff]['residues'] = delete(_udff_data[ff]['residues'], indices_to_remove)

        # rename as -XXX- if residue XXX is connected to other residues
        # this relies on self._reskeys
        self.checkInterconnectedResidues(_udff_data)

        # YL FIXME 15/02/2021: still necessary?
        # remove repeated residues
#        for ff, data in _udff_data.items():
#            unique_residues, unique_indices = unique(data['residues']['type'], return_index=True)
#            _udff_data[ff]['residues'] = data['residues'][unique_indices]

        if self.override_builtin:
            # YL 24/03/2021: to override the built-in library we must rename the scheme's name
            for ff in self._dlftypes.keys():
                _udff_data['USER_'+ff] = _udff_data.pop(ff)

        # back up and write out
        fileutils.archiveFile(self.udff, dir=self.directory)
        print('\n Writing out merged DL_FIELD and supplimented ChemShell data to:\n {}'.format(self.udff))
        file.dump(self.udff, 'udff', _udff_data)


    def checkInterconnectedResidues(self, udff_data):
        '''Rename as -XXX- if residue XXX is connected to other residues'''

        from numpy import append, array, char, copy, isin, unique

        interconnected_residues = self.frag.getInterconnectedResidues()
        interconnected_residues = char.decode(char.strip(char.upper(unique(char.ljust(interconnected_residues, 4)))))

        for ff, ff_reskeys in self._reskeys.items():

            # YL 27/01/2021: we have to check if a residue is connected to others via bonds (interconnected):
            #                if so the DL_FIELD molecule type must be "-RESNAME-"
            interconnected_reskeys = array(ff_reskeys)[isin(ff_reskeys, interconnected_residues)]
            residue_keys  = udff_data[ff]['residues']['key']
            residue_types = udff_data[ff]['residues']['type']
            copy_interconnected_residues = array([], dtype=udff_data[ff]['residues'].dtype)

            for interconnected_reskey in interconnected_reskeys:

                if interconnected_reskey in residue_keys:

                    if not '-'+interconnected_reskey+'-' in residue_keys:
                        masks_interconnected_keys = residue_keys == interconnected_reskey

                        # YL 17/02/2021: skip single-atom residues
                        is_single_atom_residue = False
                        for i, residue_key in enumerate(residue_types[masks_interconnected_keys]):
                            residue_atom_abuff = udff_data[ff].get(residue_types[masks_interconnected_keys][i], {}).get('atom', [])
                            if len(residue_atom_abuff) == 1:
                                is_single_atom_residue = True
                        if is_single_atom_residue:
                            continue

                        residue_types_renamed = char.add(residue_types[masks_interconnected_keys], '_connected')

                        # first update the dictionary then rename
                        for i, residue_type_renamed in enumerate(residue_types_renamed):
                            udff_data[ff].update({residue_type_renamed:udff_data[ff].get(residue_types[masks_interconnected_keys][i], {})})

                        # YL 17/02/2021: copy rather than just renaming because some of the residue type may be interconnected and some remain unconnected
                        copy_interconnected_residues = append(copy_interconnected_residues, copy(udff_data[ff]['residues'][masks_interconnected_keys]))

                        udff_data[ff]['residues']['key'] [masks_interconnected_keys] = '-'+interconnected_reskey+'-'
                        udff_data[ff]['residues']['type'][masks_interconnected_keys] = residue_types_renamed
                        print("\n ChemShell/DL_FIELD: molecule key for interconnecting \"{}\" -> \"-{}-\"".format(interconnected_reskey, interconnected_reskey))
            # add the copies together
            udff_data[ff]['residues'] = append(udff_data[ff]['residues'], copy_interconnected_residues)

    @property
    def dl_field_path(self):
        '''The path to DL_FIELD'''

        from os import path

        try:
            return self._dl_field_path

        except:
            # read in dl_field_path from a compiletime file
            with open(path.join(path.dirname(__file__), 'dl_field_path'), 'r') as fp:
                self._dl_field_path = fp.read().strip()
            return self._dl_field_path


    @property
    def all_resnames(self):
        '''An array of all residue names in the current fragment object'''

        from numpy import char, unique

        try:
            return self._all_resnames

        except:
            self._all_resnames = char.decode(char.strip(char.upper(unique(char.ljust(self.frag.residues, 4)))))
            return self._all_resnames


    # note: relies on self._reskeys
    @property
    def masks_residues_found(self):
        '''Mask array of self.all_resnames found in the self._reskeys dict'''

        from numpy import isin, zeros

        self._masks_residues_found = zeros((self.all_resnames.size,), dtype=bool)

        for ff, ff_reskeys in self._reskeys.items():
            self._masks_residues_found += isin(self.all_resnames, ff_reskeys)

        return self._masks_residues_found


    @property
    def all_residues_found(self):
        '''Return True if all residues are defined in forcefield; False otherwise'''

        from numpy import char, isin 

        # screen if all residues can be found in a single FF scheme
        all_residues_found = False

        # all residue names (including patched ones) found in .PS file
        # NOT IN USE: remove trailing characters for example 'CAM_CHARMM' that is used to distinguish from DL_POLY keys
        if all(self.masks_residues_found) or \
           all(isin(self.all_resnames[~self.masks_residues_found], list(self._associatedFFLib.titratables.keys()))):# or \
        #   all(in1d(self.all_resnames, char.partition(ff_reskeys, '_'+ff)[:,0])):
            all_residues_found = True

        return all_residues_found



    def splitFragment(self):
        '''Split the whole fragment to use the one-scheme-per-segment method'''

        matched = {}

        split_frags = self.frag.split(save_to_pdb=False,
                                      exclude=[],
                                      groups=kwargs.get('solvents', []),
                                      return_dict=True)

        with open(self.filename, 'a') as fp:

            for key, val in split_frags.items():
                print("\n Processing segment \"{}\"...".format(key))

                for ff, ff_reskeys in self._reskeys.items():

                    # if all residue names are defined in SF
                    res_in_reskeys = isin(self.all_resnames, ff_reskeys)

                    # all residue names (including patched ones) found in .PS file
                    if all(res_in_reskeys) or all(isin(self.all_resnames[~res_in_reskeys], list(self._associatedFFLib.titratables.keys()))):
                        print("\n All residues (molecule types) of segment {} have been found".format(key))
                        matched.update({key:ff.upper()})
                        ff = self._getSchemeNameFromUDFF(ff)
                        val.title = 'POTENTIAL {}'.format(ff.upper())
                        val.save(fp, self._format, uppercase=False)

                    else:
                        # no molecular types defined in SF at all
                        if ff_reskeys is None and key not in matched:
                            print(" \n No molecule type of segment {} is defined".format(key))
                            matched.update({key:ff})
                            val.title = 'POTENTIAL {}'.format(ff)
                            val.save(fp, self._format)
                        # if scheme needs no template (such as CVFF, PCFF, OPLS2005, etc.)
                        if ff.lower() in ['cvff','pcff'] and key not in matched:
                            print(" \n No molecule type of segment {} is defined but forcefield scheme {} does not need a template".format(key, ff.upper()))
                            matched.update({key:ff.upper()})
                            val.title = 'POTENTIAL {}'.format(ff.upper())
                            val.save(fp, self._format)

            print("\n Matched segments and forcefield schemes:")
            for key, val in matched.items():
                print(" {:<6}: {}".format(key, val))


    @staticmethod
    def charmm2udff(**kwargs):
        '''A CHARMM-to-DL_FIELD converter, saving to output file'''

        from . import tools

        return tools.charmm2udff(**kwargs)


    @staticmethod
    def patchRTF(output='_patched.rtf', **kwargs):
        '''A wrapper around tools.patchResidues()'''

        from . import tools

        fileutils.archiveFile(output)
        return tools.patchRTF(output=output, **kwargs)


    @property
    def _format(self):
        '''File format of the self.filename'''

        return fileutils.getExtension(self.filename)


    # YL FIXME: 19/02/2021 it's OK since we only support CHARMM now but for the
    #           future developement with more forcefield schemes we need to expand
    #           this method
    @property
    def _associatedFFLib(self):
        '''Associated forcefield library'''

        from ...data import forcefield, charmm

        selectcases = { 'charmm': charmm,
                      }

        try:
            for k, v in selectcases.items():
               if self.ff.lower().strip().startswith(k):
                   builtin_ff = k
                   break
        except:
            builtin_ff = 'charmm'

        return forcefield.getFFObj(builtin_ff)


    @property
    def fftext(self):
        '''Plain text of FF'''

        try:
            return self._fftext
        except AttributeError:
            self.run()
            return self._fftext


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

        self._fftext = val


    @property
    def udff(self):
        '''UDFF filename to write ('' if no UDFF requested)'''

        from ...objects import fragbio

        # if not defined by the user (default is '')
        if not self._udff:
            # PDB for biological system
            if isinstance(self.frag, fragbio.FragmentBio) and self._format == 'pdb':
                self._udff = '_dl_field.udff'

        return self._udff


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

        self._udff = val

        return self._udff


    @property
    def _scheme_name(self):
        '''Return a forcefield scheme name according to self.ff'''

        return self._getSchemeNameFromFF(self.ff)


    @property
    def _standard_ff_name(self):
        '''Standard forcefield name in lower case, for example charmm or amber'''

        for name in [ 'charmm', 'amber' ]:
            if name in self._scheme_name.lower():
                return name

        return 'unknown'


    @property
    def _ffObj(self):
        '''Return the forcefield module object accoriding to self._standard_ff_name'''

        from ...data import forcefield

        return forcefield.getFFObj(self._standard_ff_name)


    # YL 24/11/2023: patch residues of fragment in place
    def patchFragment(self, prefix='_patched'):
        '''Patch residues in place'''

        from numpy import char, isin, unique, zeros

        dict_rename = {}
        dict_deletion = {}
# YL 26/11/2023: not in use
#        list_atoms_to_delete = []

        # reconstruct patching information from the dict data, including
        # deletetion and renaming
        if type(self.ff) is list:
            for ff in self.ff:
                if type(ff) is dict:
                    for k_ff, v_ff in ff.items():
                        if type(v_ff) is dict:
                            # k_item will be the new residue name
                            for k_item, v_item in v_ff.items():
                                if type(v_item) is dict:
# YL 26/11/2023: not in use
#                                    if 'patch' in v_item:
#                                        dict_deletion.update({k_item:v_item['patch']})
                                    if v_item.get('renamed_from', ''):
                                        dict_rename.update({v_item['renamed_from']:k_item})

        if self.patch_target:
            print('\n --------------------')
            print(' Patching fragment...')
            print(' --------------------\n')

        # seg must be str/bytes
        for seg, residues in self.patch_target.items():
            print(f'\n Segment: {seg}')
            # this could save a bit time
            masks_seg = self.frag.select(segments=seg, return_masks=True)
            for res in residues:
                # resID
                if miscutils.isInteger(res):
                    masks_res = masks_seg*self.frag.select(resIDs=res, return_masks=True)
                # resname
                if strutils.isStr(res):
                    masks_res = masks_seg*self.frag.select(residues=res, return_masks=True)
                # could be multiple residues
                resIDs_to_rename = unique(self.frag.resIDs[masks_res])
                # go through one by one
                for i, resID_to_rename in enumerate(resIDs_to_rename):
                    # the returned value is str
                    resname_old = self.frag.resID2residue(resID_to_rename, segment=seg).strip()
                    # NB: the trailing space must be stripped
                    if dict_rename.get(resname_old, ''):
                        print(f' Renaming residue {resname_old}{resID_to_rename} to {dict_rename[resname_old]}{resID_to_rename}')
                        self.frag.residues[masks_res] = dict_rename[resname_old]
# YL 26/11/2023: not in use
#                    # it could be possible that the patched residue is not renamed
#                    resname_to_patch = dict_rename.get(resname_old, resname_old)
#                    if resname_to_patch in dict_deletion:
#                        for operation in dict_deletion[resname_to_patch]:
#                            if operation.startswith('# ATOM REMOVE'):
#                                atom_to_delete = operation.split('# ATOM REMOVE')[1]
#                                masks_to_delete = in1d(char.strip(char.decode(self.frag.names[masks_res])), atom_to_delete)
#                                list_atoms_to_delete += self.frag._indices_all[masks_res][masks_to_delete].tolist()
 
                stdout.flush()

        # save resulting fragment to new files
        fileutils.archiveFile(prefix+'.pdb')
        fileutils.archiveFile(prefix+'.pqr')
        self.frag.save(prefix+'.pdb')
        self.frag.save(prefix+'.pqr')

        if self.patch_target:
            print('\n ---------------------------')
            print(' Finished patching fragment.')
            print(' ---------------------------\n')













