#  Copyright (C) 2017 The authors of Py-ChemShell
#
#  This file is part of Py-ChemShell.
#
#  Py-ChemShell is free software: you can redistribute it and/or modify
#  it under the terms of the GNU Lesser General Public License as
#  published by the Free Software Foundation, either version 3 of the
#  License, or (at your option) any later version.
#
#  Py-ChemShell is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with Py-ChemShell.  If not, see
#  <http://www.gnu.org/licenses/>.

import chemsh
from .               import _Interfaces
from ..base          import run
from ..io            import keyword
#from ..objects import Fragment
from ..dl_py2f       import utils
from ..dl_py2f.utils import miscutils

# inherit <ctypes.Structure> type for convenience of making callback
class _QM(_Interfaces):
    '''A container of QM methods: all QM interfaces should be of this class'''

    from ctypes import RTLD_GLOBAL
    from ..cluster import partition
    from numpy import zeros, float64

    # default mode for ctypes.CDLL() function, see modutils
    _mode = RTLD_GLOBAL

    _kwargs = {
                'analysis'         :'',       # wavefunction analysis
                'atom_names'       :'ff',
                'basis'            :'3-21g',  # also use `basis` to control `basisspec`
                'charge'           :(0, 0.0, ''),
                'damping'          : True,
                'diis'             : True,
                'direct'           : True,
                'dispersion'       : '',
                'ecp'              : '',      # AJL, Aug2017. Used by more than one calculator so included here rather than in child classes
                'ecps'             : None,    # needed by NWChem Fortran code, otherwise could've been kept private
                 # YL 18/11/2019: changed to a user option at QM interface level again and removed from QMMM
                'ecp_region'       : partition._kwargs['label_interface'],
                'freetext'         :'',      # a mult-line string containing a block of user-defined input
                'freetext_only'    : False,  # content of freetext should replace ChemShell-generated part or just be used as an insertion
                'functional'       :'blyp',
                'ghost'            :[],
                'ghost_prefix'     :'X_',
                'grads_qmme'       : zeros(shape=(1,3), dtype=float64, order='C'),
                'grid'             :'medium',
                'guess'            :'',
                'harmonic'         : False,
                'input'            : None,    # dummy definition, to be overwritten by calculator (Needed for FHI-aims interface)
                'level_shift'      : False,
                'level_e1'         : 1.0,
                'level_e2'         : 0.3,
                'level_cycles'     : 999,
                'logger'           :'',       # error logger for stderr
                'maxiter'          : 100,
                'method'           :'dft',
                'mult'             : 0,
                'necps'            : 0,
                'nstates'          : 0,
                'path'             :'',
                'pop_analysis'     :'',       # HNT, Jul2020. Added population analysis keyword for mulliken.
                'properties'       :'',
                'pseudopotential'  :'',
                'restart'          : False,
                'ri'               :'',       # YL 16/02/2024: a str val of RI basis
                'rpa'              : True,
                'target_state'     : 0,
                'threshold'        : 1.0e-6,
                #Rajany: keeping as it is for now.
                'scale_nonperiodic': 3.0,
                # TD 16/05/2024: Must be left as '' to enable the default to depend on mult (mult=1 ==> scftype='rhf', mult>1 ==> scftype='uhf')
                'scftype'          :'',
                'separate_ecps'    : False,   # YL 18/11/2019: changed to a user option at QM interface level again and removed from QMMM
                'soscf'            : True,    # YL 21/02/2024: approximate second order SCF
                'spin'             :'singlet',
                'sysmpi'           :'',
              }

    _internals = {
                   '_ecp_charge'   : 1.0E-8,
                   '_syscmd'       : '',
                   '_extra_outputs': [],
                   '_recalc_cell'  : True,
                   '_qmme'         : 0,
                   '_basis_str'    : '',
                 }

    # set up at the QM/MM level
    @property
    def linkatoms(self):
        '''List of link atom indices'''

        return self._linkatoms

    @property
    def grads_qmme(self):
        return self._grads_qmme

    @grads_qmme.setter
    def grads_qmme(self, val):
        self._grads_qmme = val


    # called by run.runExec()
    def genInput(self, **kwargs):
        '''Generate input'''

        # single input file
        if type(self.input) in [ str ]:
            utils.fileutils.write(self.input, self.getInpBuff())

        # multiple input files
        # JCRT: Added tuple/list of tuples/lists
        elif type(self.input) in [ list, tuple ]:
            for inpt in self.input:
                if type(inpt) in [list, tuple]:
                    for inpt2 in inpt:
                        if inpt2 != None:
                            utils.fileutils.write(inpt2, self.getInpBuff(input=inpt2))
                elif inpt != None:
                    utils.fileutils.write(inpt, self.getInpBuff(input=inpt))



    def parseMainOutput(self):
        '''Parse main output file self.output'''

        # do nothing in QM template


    def parseAdditionalOutput(self):
        '''Parse additional output data, e.g., nwchem.bq_grads. This procedure could be overloaded, otherwise ignored.'''


    def setOutKeys(self):
        '''Set up {task:keyword} for output parsing. Shall be overloaded if needed.'''


    def parseOutput(self):
        '''Parse output'''

        self.parseMainOutput()

        self.parseAdditionalOutput()


    @miscutils.time
    def run(self, **kwargs):
        '''Single-point procedure'''


        from sys     import stderr, stdout
        from ..io    import banner, format
        from ..data  import references, unitconvert

        # JCRT: We have to tell the qm method we are using gradients here, that way if multiple executables are required we can add them now via @_gradients.setter
        try:
            self._gradients = kwargs["_gradients"]
        except:
            pass

        if not self.frag.natoms:
            print("\n\n >>> ChemShell ERROR: no valid fragment has been assigned to theory {}!\n\n".format(self.__class__.__name__))
            exit(999)

        if type(self.input) == tuple:
            for inpt in self.input:
                if inpt != None:
                    utils.fileutils.archiveFile(inpt)
        else:
            utils.fileutils.archiveFile(self.input)
        utils.fileutils.archiveFile(self.output)

        banner.printModuleBanner(f'ChemShell/{self.__class__.__name__} Interface', level=2)
#        if self.__class__.__name__.lower() in references.references:
#            print(' Please cite:')
#            format.printwrap(references.references[self.__class__.__name__.lower()], prepend_blanc_line=False)

        if chemsh.CHEMSH_DEBUG > 4:
            print(" {} debug mode:".format(__class__))
            print(" path         = {}".format(self.path))
            print(" getLinkedLib = {}".format(utils.modutils.getLinkedLib(self.__class__.__name__, chemsh.CHEMSH_LIB_DIR, mode=self._mode)))
            print(" syscmd       = {}".format(self._syscmd))

        # YL: we allow users to override linked-in libs with `path`: check path first then isLinked then _syscmd
        cmdpath = self.path or utils.modutils.getLinkedLib(self.__class__.__name__, chemsh.CHEMSH_LIB_DIR, mode=self._mode) or self._syscmd

        # don't use higher level task routines here (need to keep lower level separate)
        # all keyword arguments at the run() step will be temporary
        # YL 11/2020: run twice for the QM/Me periodic embedding scheme
        if self._qmmm_scheme == "qm/me":
            ierror = self.__runQMMe(cmdpath, **kwargs)

        # standard methods
        else:
            ierror = run.run(self, cmdpath, stderr=self.logger, shell=self._exec_run_shell, **kwargs)
            stdout.flush()
            stderr.flush()

        # combine the BQ gradients with contribution from separate ECPs
        if self.necps and self.separate_ecps and self._gradients and kwargs.get('_gradients', False) != 'shells':
            # YL 25/07/2021: planned to change to `self.frag.bqs.gradients[self._ecps._range] += self._ecps.gradients`
            #                this should be safe because self._ecps depends on self.bq.coords only which is not changed during QM.run()
            self.frag.bqs.gradients[self.ecps._range] += self.ecps.gradients

        # excitation energies
        if self._result.excitations.size:

            print('\n Summary of QM multi-state calculation:')
            print(' ' + '-'*70)
            sbuff = f'(E{self.target_state})'
            sbuff2 = ' '*(len(sbuff)-6)
            print(f' Ground state       (E0){sbuff2} = {self._result.E0:25.12f} a.u.')
            if self.target_state:
                print(f' State of interest* {sbuff} = {self._result.energy[0]:25.12f} a.u.')
            print(' ' + '-'*70)

            print('\n Excited states calculated:')
            print(' ' + '-'*70)
            print(' State Symm. Spin  deltaE/eV          E/a.u.      lambda/nm   Osc. str.')
            print(' ' + '-'*70)
            # 4 decimals may not be sufficient but let's keep it for now as it's the best print by Gaussian
            for i, e in enumerate(self._result.excitations):
                if i == self.target_state-1:
                    s_index = str(i+1)+'*'
                else:
                    s_index = str(i+1)+' '
                dE = (e['energies']-self._result.E0)*unitconvert.Ha2eV
                print(' {:>4}   {:<5} {:>1} {:>12.4f} {:>19.12f} {:>10.2f} {:>11.6f}'.format(s_index,
                                                                                      e['symmetries'].decode().upper(),
                                                                                      e['spins'].decode().upper()[:1],
                                                                                      dE,
                                                                                      e['energies'],
                                                                                      e['wavelengths'], e['f']))
            print(' ' + '-'*70)
            sbuff = f'{format.integer2english(self.target_state)} excited state'
            if self.target_state:
                print(f' * State of interest: {sbuff} (i.e., target_state={self.target_state})\n')
            else:
                print()

        return ierror


    @property
    def charge(self):
        '''QM charge (must be an integer)'''

        if self._charge == 'automatic':
            return round(self.frag.charges.sum())

        return self._charge


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

        # YL 02/06/2022: deffered QM charge, determined by the MM mechanism later (see: dl_poly/__init__.py)
        if type(val) is str:
            if val.lower().strip().startswith('auto') or val.lower().strip() in [ 'dl_field', 'dl-field', 'dlfield' ]:
                self._charge = 'automatic'

        # NB: default val must be assigned explicitly because in the `_attr` dict it is a tuple of multiple legal types
        elif type(val) is tuple:
            self._charge = val[0]

        # int or float
        else:
            self._charge = round(float(val))
            if float(val) != self._charge:
                print(f'\n >>> Given QM charge value {val} rounded to {self._charge}')


    @property
    def logger(self):
        '''Error logger filename'''

        if not self._logger:
            return '_{}.err'.format(self.__class__.__name__.lower())
        else:
            return self._logger


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

        self._logger = val


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

        try:
            return self._runargs
        except:
            # default is just self.input
            return self.input


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

        self._runargs = val


    @property
    def restart(self):
        '''Restart from a previous calculation'''

        return self._restart


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

        self._restart = val


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

        return '_' + self.__class__.__name__.lower()


    @property
    def _extrafiles(self):
        '''Additional output files'''

        return []


    @property
    def path(self):
        '''Path of command to execute the code'''

        return self._path


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

        self._path = val


    # YL 22/08/2019: we hope to deal with ghost atoms' basis set automatically
    def getBasisWithGhost(self):
        '''Automatically adapt basis set for ghost atoms'''

        from numpy import char, unique, setxor1d

        strbuff = self.basis

        if len(self.ghost) > 0:
            lstrbuff = strbuff.splitlines()
            ghostbuff = '\n# ghost basis set adapted by Py-ChemShell\n'
            ghost_species = char.decode(unique(self.frag.names[self.ghost]))

            lastline = len(lstrbuff)
            for i, l in enumerate(lstrbuff):
                if l.strip().lower().startswith('end'):
                    lastline = i
                    break
            defined_ghosts = []
            # check if already defined
            for g in ghost_species:
                for i, l in enumerate(lstrbuff):
                    if l.strip().lower().startswith(self.ghost_prefix.lower()+g.lower()):
                        print(" \n Basis set for ghost atom \"{}\" detected.\n".format(self.ghost_prefix+g))
                        defined_ghosts += [g]
                        break

            # loop over ghost atoms
            for g in setxor1d(ghost_species, defined_ghosts):
                print(" Adapting basis set for ghost atom \"{}\"...".format(self.ghost_prefix+g))
                ghostbuff += '\n# {}\n'.format(self.ghost_prefix+g)
                # find starting line
                other_species = setxor1d(char.decode(unique(self.frag.names)), [g])
                for i, l in enumerate(lstrbuff):
                    if l.strip().lower().startswith(g.lower()):
                        startline = i
                        break
                # find ending line by matching the first other species
                endline = len(lstrbuff)
                for o in other_species:
                    for i, l in enumerate(lstrbuff[startline+1:]):
                        if l.strip().lower().startswith(o.lower()) and endline>startline+i+1:
                            endline = startline+i+1
                        if l.strip().lower().startswith('end') and endline>startline+i+1:
                            endline = startline+i+1
                            break
                for i, l in enumerate(lstrbuff[startline:endline]):
                    if l.strip().lower().startswith(g.lower()):
                        ghostbuff += self.ghost_prefix + l.strip() + '\n'
                    # get rid of comments
                    elif l.strip().startswith("#"):
                        continue
                    elif l.strip() == '':
                        continue
                    else:
                        ghostbuff += l+'\n'
            lstrbuff.insert(lastline, ghostbuff)
            strbuff = '\n'.join(lstrbuff)

        return strbuff


    @property
    def basis(self):
        '''Basis set: returns a standard name of basis set or the content of basisdeck or whats given'''

        from ..data  import basis

        # YL 21/07/2022: using basis set engine such as MolSSI-BSE
        if hasattr(self._basis, 'run'):
            
            if not self._basis_str:
                self._basis_str = self._basis.run(frag    = self.frag,
                                                  fmt     = self.__class__.__name__.lower(),
                                                  default = self._kwargs['basis'])
            return self._basis_str

        else:
            # search in the standard basis set names, if not found return the content of string or file
            return basis.basissets.get(self._basis.strip().lower().replace('-',''),
                                       utils.fileutils.file2text(self._basis, preceding_newline=True, error_return=self._basis))

#            # search in the standard basis set names, if not found return the content of string or file
#            strbuff = basis.basissets.get(self._basis.strip().lower().replace('-',''),
#                                          '\n'+fileutils.file2text(self._basis, error_return=self._basis))
#            return strbuff


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

        self._basis = val


    @property
    def ghost(self):
        '''Ghost atoms'''

        return self._ghost
       

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

        from numpy   import unique

        self._ghost = utils.iterutils.getFlattened([val])

        nghosts = len(self._ghost)

        # can't print names before frag may not be ready
        if nghosts == 1:
            print(" >>> {}: ghost atom {} requested".format(self.__class__.__name__, self._ghost))
        elif nghosts > 1:
            print(" >>> {}: ghost atoms {} requested".format(self.__class__.__name__, self._ghost))


    @property
    def functional(self):
        '''DFT functional'''

        from ..data import dft

        return dft.functionals.get(self._functional.lower().replace('-',''),
                                   self._functional.lower())


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

        self._functional = val


    @property
    def grid(self):
        '''Grid accuracy'''

        selectcases = { 'low'      :'low',
                        'medium'   :'medium',
                        'high'     :'high',
                        'veryhigh' :'veryhigh',
                        'ultrahigh':'ultrahigh',
                      }

        return selectcases.get(self._grid.lower(), '')


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

        self._grid = val


    @property
    def pop_analysis(self):
        '''Defines population analysis methods'''

        analysis_method = { 'mulliken' : 'mulliken',
                            'fukui' : 'fukui', 
                          }

        return analysis_method.get(self._pop_analysis.lower(), '')


    @pop_analysis.setter
    def pop_analysis(self, val):
        '''Setter of population analysis method'''
        self._pop_analysis = val



    @property
    def scftype(self):
        '''Type of SCF'''

        # default is RHF for singlets and UHF for higher spin multiplicities
        if not self._scftype:
            if self.mult == 1:
                return 'rhf'
            else: 
                return 'uhf'

        return self._scftype


    @scftype.setter
    def scftype(self, val):
        '''Setter of SCF type'''

        self._scftype = val


    # YL 18/11/2019: removed from _kwargs
    @property
    def necps(self):
        '''Number of ECP point charges'''

        try:
            return self.ecps._master.shape[0]
        except:
            return 0


    # only used by initialisation
    @necps.setter
    def necps(self, val):
        '''Setter of necps'''


    # YL 18/11/2019: refactored
    @property
    def ecps(self):
        '''An instance of <class BQs> for separate ECPs'''

        from numpy import array, isin

        # return if already defined to avoid adding/subtracting _ecp_charge for more than once
        try:
            #  YL 25/04/2020: update coordinates of ECPs (for separate ECPs where shells are retained)
            self._ecps.coords = self.frag.bqs.coords[self._ecps._range]
            return self._ecps

        except:

            # requested
            if self.separate_ecps:
                try:
                    # BQs.select() has problems with methods such as selectByRadius()
                    ecp_range = self.frag.bqs.selectByRegions(self.ecp_region)
                    # exclude shell centres in separate EPCs
                    shl_range = array(self.frag.bqs.stacks['Shells'])
                    ecp_range = ecp_range[isin(ecp_range, shl_range, invert=True)]
                    if len(ecp_range):
                        self._ecps = self.frag.bqs.getSelected(ecp_range)
                        self._ecps._range = ecp_range
                        self._ecps.charges = self._ecp_charge
                    else:
                        self._ecps = None
                except:
                    self._ecps = None

            # not requested 
            else:
                self._ecps = None

            return self._ecps


    # only used by initialisation
    # AJL, June2020: This doesn't contain any actions? Is this safe?
    # YL 25/07/2021: it was safe thanks to the try-except clause but added `self._ecps = val` now
    @ecps.setter
    def ecps(self, val):
        '''Setter of ecps'''

        self._ecps = val


    @property
    def method(self):
        '''A wrapper for some pre-processing'''

        return self._method


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

        self._method = val.lower().strip()


    @property
    def mult(self):
        '''Multiplicity (automatically calculated but can be overridden)'''

        from ..base import forcefield

        if self._mult > 0:
            return self._mult

        else:
            znum = forcefield.names2znums(self.frag.names).sum()
            return (znum - self.charge)%2 + 1


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

        self._mult = val


    @property
    def _time_dependent(self):
        '''Return True if time-dependent excited state calculations are requested'''

        # supported time-dependent methods
        if self.method.lower().strip() in [ 'td-dft', 'td-hf', 'tddft', 'tdhf' ]:
            return True
        else:
            return False


    @property
    def _excited_state(self):
        '''Return True if time-dependent excited state calculations are requested'''

        # supported methods
        if self._time_dependent or self.method.lower().strip() in [ 'mrci', 'mr-ci', 'ccsd', 'ccsdt', 'casscf', 'caspt2' ]:
            return True
        else:
            return False


    @property
    def nstates(self):
        '''Number of states requested'''

        # defaults to 3 if not set by the user
        if self._excited_state and not self._nstates:
            self._nstates = 3

        return self._nstates


    # YL: it won't work if we initialise self._result.excitations by resizing
    #     because `_kwargs` are processed before `_internals` and `self._result`
    #     will be overwritten anyway
    #     so we need to override the _result() method to include multistate results
    @nstates.setter
    def nstates(self, val):
        '''Setter of nstates'''

        self._nstates = val


    def __runQMMe(self, cmdpath, **kwargs):
        '''Run QM/Me method for periodic embedding'''

        from sys    import stdout
        from copy   import copy as pycopy
        # YL 10/02/2021: gcd is not available until NumPy v1.15
        #                temporarily use math.gcd
        from numpy  import arange, array, copy
        from math   import gcd
        from ..base import lattice
        from ..io   import format

        # choose a surface
        try:
            # should be already determined by <class QMMM>
            surface = self._surface
        except:
            print("\n WARNING: no surface has been chosen so \"{}\" will be used as default\n".format('ab'))
            self._surface = 'ab'
            surface = self._surface

        # non-periodic dimension
        dimension = lattice.dimensions[surface]

        # remember filenames
        input  = pycopy(self.input)
        output = pycopy(self.output)

        print(' '+'-'*70)
        print("\n QM/Me scheme: running the 1st QM calculation with adsorbate\n")
        print(' '+'-'*70)
        stdout.flush()

        # cell.vectors were used implicitly, which were created in QMMM.qm
        # recalculate the lattice constants in the beginning
        consts = array(self.frag.getLatticeOfSelectedPeriodic(arange(self.frag.natoms), surface=surface, recalculate=True))

        if self._recalc_cell:
            # scale the non-periodic vector to make sure the volume is large enough
            # NB: this is required for some codes such as CP2K
            self.frag.cell.consts = [ p*self.scale_nonperiodic if i == dimension else p for i, p in enumerate(consts) ]
        else:
            self.frag.cell.consts = self.__saved_consts[:]

        # calculate padding distance (unit vector length) along periodic dimensions if the code doesn't allow redundant centres in cell (e.g., CP2K)
        # NB: the cell fragment must be correctly prepared without redundant centres otherwise the calculation will be wrong!
        if not self._cell_redundant and self._recalc_cell:
            # always reset vectors before calculating constants! otherwise constants are wiped according to vectors!
            self.frag.cell.vectors[:] = 0.0

            # calculate the unit vector lengths: do NOT include adsorbate atoms!
            padded_consts = array(self.frag.getLatticeOfSelectedPeriodic(arange(self.frag.natoms)[self._non_adsorbate],
                                                                         surface     = surface,
                                                                         accuracy    = 0.5,
                                                                         recalculate = True,
                                                                         redundant   = self._cell_redundant))

            # apply the padding to the periodic vectors
            self.frag.cell.consts = [ p if i != dimension else self.frag.cell.consts[i] for i, p in enumerate(padded_consts) ]

        # save the constants to keep them constant for all following runs because
        if self._recalc_cell:
            a, b, c, alpha, beta, gamma = [ self.frag._cell.a, self.frag._cell.b, self.frag._cell.c, self.frag._cell.alpha, self.frag._cell.beta, self.frag._cell.gamma ]
            self.__saved_consts = [ a, b, c, alpha, beta, gamma ]
            # keep a copy of original cell vectors for use in the second run
            self.__saved_vectors = copy(self.frag.cell.vectors)

        a, b, c, alpha, beta, gamma = [ self.frag._cell.a, self.frag._cell.b, self.frag._cell.c, self.frag._cell.alpha, self.frag._cell.beta, self.frag._cell.gamma ]
        print(format.wrap("Lattice constants (a.u.) of the selected periodic QM region automatically assigned:"))
        print(" a     = {:>15.6f}".format(a))
        print(" b     = {:>15.6f}".format(b))
        print(" c     = {:>15.6f}".format(c))
        print(" alpha = {:>15.6f}".format(alpha))
        print(" beta  = {:>15.6f}".format(beta))
        print(" gamma = {:>15.6f}".format(gamma))

        # the k-points sampling along the two periodic directions
        if self._recalc_kpoints:
            two_points = tuple(k for i, k in enumerate(self.kpoints) if i != dimension)
            # check the greatest common diviser
            # YL 11/02/2021: numpy.gcd() not available until NumPy v1.15
            d = gcd(*two_points)
            kpoints = tuple(k for k in self.kpoints)
            # set the k-points sampling 1 along the unperiodic direction
            self.kpoints = tuple(int(k/d) if i != dimension else 1 for i, k in enumerate(kpoints))

        # first run with adsorbate
        self.wavefunction = '_{}.wfn.w_adsorbate'.format(self.__class__.__name__.lower())
        utils.fileutils.archiveFile(self.input +'.w_adsorbate')
        utils.fileutils.archiveFile(self.output+'.w_adsorbate')

        # temporary filenames (will be restored after use)
        self.input  = input +'.w_adsorbate'
        self.output = output+'.w_adsorbate'

        ierror = run.run(self,
                         cmdpath,
                         stderr=self.logger,
                         # special options per QM interface
                         _qmme=1,
                         **kwargs,
                        )

        # geometry optimisation (DL_FIND) archiving because we cannot do with the temporary filenames there
        try: 
            if chemsh.CHEMSH_DEBUG:
                utils.fileutils.archiveFile([self.input, self.output], dir=self._archive_dir, suffix=self._archive_suffix)
        except:
            pass

        e1 = pycopy(self._result.energy[0])
        print("\n QM/Me scheme: 1st QM energy (with adsorbate): {} a.u.\n".format(e1))

        print(' '+'-'*70)
        print("\n QM/Me scheme: running the 2nd QM calculation without adsorbate\n")
        print(' '+'-'*70)

        self.frag = self.frag.getSelected(self._non_adsorbate)

        # restore cell constants/vectors as they should be kept constant
        self.frag.cell.vectors[:] = self.__saved_vectors
        self.frag.cell.consts = self.__saved_consts[:]

        a, b, c, alpha, beta, gamma = [ self.frag._cell.a, self.frag._cell.b, self.frag._cell.c, self.frag._cell.alpha, self.frag._cell.beta, self.frag._cell.gamma ]
        print("\n Lattice constants (a.u.) of the selected periodic QM region automatically assigned:")
        print(" a     = {:>15.6f}".format(a))
        print(" b     = {:>15.6f}".format(b))
        print(" c     = {:>15.6f}".format(c))
        print(" alpha = {:>15.6f}".format(alpha))
        print(" beta  = {:>15.6f}".format(beta))
        print(" gamma = {:>15.6f}".format(gamma))

        # second run without adsorbate
        self.wavefunction = '_{}.wfn.wo_adsorbate'.format(self.__class__.__name__.lower())
        utils.fileutils.archiveFile(self.input +'.wo_adsorbate')
        utils.fileutils.archiveFile(self.output+'.wo_adsorbate')

        # temporary filenames (will be restored after use)
        self.input  = input +'.wo_adsorbate'
        self.output = output+'.wo_adsorbate'

        # TODO: codes that do not have these keywords
        ierror = run.run(self,
                         cmdpath,
                         stderr=self.logger,
                         # special options per QM interface
                         _qmme=2,
                         **kwargs,
                        )

        # geometry optimisation (DL_FIND) archiving because we cannot do with the temporary filenames there
        try: 
            if chemsh.CHEMSH_DEBUG:
                utils.fileutils.archiveFile([self.input, self.output], dir=self._archive_dir, suffix=self._archive_suffix)
        except:
            pass

        # restore filenames
        self.input  = input
        self.output = output

        self._recalc_cell = False

        e2 = pycopy(self._result.energy)
        self._result.energy = e1-e2
        print("\n QM/Me scheme: 2nd QM energy (without adsorbate): {} a.u.\n".format(e2))
        print("\n QM/Me scheme:  QM energy : {} a.u.\n".format(self._result.energy[0]))

        self._result.gradients[self._non_adsorbate] -= self.grads_qmme

        return ierror


    # YL: override `_Interface._result` to enable multistate results
    @property
    def _result(self):
        '''Result containing energy and gradients'''

        from ..objects import result

        # associated to task.result
        try:
            # initialise the multi-state result array (in-place operation)
            if self.__result.excitations.size != self.nstates:
                self.__result.excitations.resize(self.nstates, refcheck=False)
            return self.__result

        # for example: when a QM/MM task does not initialise task.qmmm.qm._result or task.qmmm.mm._result
        except:
            self.__result = result.Result()
            self.__result.excitations.resize(self.nstates, refcheck=False)
            return self.__result


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

        self.__result = val




        


