#  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/>.

from sys import stdout
import chemsh
from ..qm        import _QM
from ...dl_py2f  import utils
from ...parallel import parallel

class NWChem(_QM):
    '''Interface to NWChem'''

    # FIXME: options of opt, e.g. by passing in the dlf command
    # TODO: Python doesn't allow *args to be after **kwargs, but we need more flexibility for input
    _kwargs = { 
                'accaofunc'     : 25, # this is not an option in the manual, see: http://www.nwchem-sw.org/index.php/Release66:Density_Functional_Theory_for_Molecules
                'bq_grads'      :'_nwchem.bq_grads',
                'dampon'        : 0.0,
                'dampoff'       : 0.0,
                'diison'        : 0.0,
                'diisoff'       : 0.0,
                'formal_charges':{},
                'fragmovecs'    :'',
                'hl_tol'        : 0.1,
                'incore'        : False,
                'input'         :'_nwchem.inp',
                'integral_io'   :'',
                'levlon'        : 0.0,
                'levloff'       : 0.0,
                'lindep_n'      : -1,
                'lindep_tol'    : 1.E-5,
                'lshift'        : 0.5,
#                'maxnbqs'      : 25000,
                'mem_heap'      : 100,
                'mem_stack'     : 100,
                'mem_global'    : 200,
                'memory'        : 400,
                'ncydp'         : 2,
                'ncyds'         : 30,
                'ncysh'         : 30,
                'necps'         : 0,            # required by NWChem Fortran code
                'ndamp'         : 0,
                'nfock'         : 10,
                'output'        :'_nwchem.out',
                'pre_hf'        : False,
                'project'       :'',
                # YL 18/11/2019 NWChem only supports separate ECPs
                # TWK 23/06/2020 Note separate ECPs are only applicable to the ionic embedding model, will be overriden for covalent coupling
                'separate_ecps' : True,
                'tol_energy'    : 1.E-6,
                'tol_density'   : 1.E-5,
                'tol_gradient'  : 5.E-4,
                'tol_acccoul'   : 8,
                'tol_pre_hf'    : 1.E-5,
                'tol_rho'       : 1.E-10,
                'tol_radius'    : 25.0,
                'vecs_in'       :'',
                'vecs_out'      :'_nwchem.movecs'
              }

    _internals = {
                   '_ecp_charge'   :  1.0E-20, # the default 1.0E-8 isn't small enough
                   '_syscmd'       : 'nwchem', # default system command to run
                   '_extra_outputs': ['_nwchem.bq_grads'], # must be a list
                 }

    # TWK: translate common functional names into NWChem format
    # YL: place all synonyms here
    _synons = {
                'blyp'  : 'becke88 lyp',
                'lda'   : 'slater vwn_5',
                's-vwn' : 'slater vwn_5',
                'hf'    : 'scf',
                'scf'   : 'scf',
                'dft'   : 'dft',
                'tddft' : 'tddft',
                'td-dft': 'tddft',
              }

   
    # FIXME: change to @property?
    # must be a procedure rather than a direct definition, otherwise frag.natoms is unknown
    def setOutKeys(self):
        '''Define {task:keyword} for searching in output file'''

        # each keyword should contain an array of integers specifying:
        #     [ row_0, # starting row
        #       row_n, # ending row
        #       row_i, # row increment
        #       col_0, # starting col
        #       col_n, # ending col
        #       col_i, # col increment
        #       skip , # skip the previous (skip-1) match,
        #       coeff, # unit conversion coefficient
        #     ]
        # referencing to the line where the key is matched
        self._outkeys = {
                          'energy'   :{ 'Total SCF energy': [0,0,1,4,4,1,1,1.0],
                                        'Total DFT energy': [0,0,1,4,4,1,1,1.0]},
                          'gradients':{ 'ENERGY GRADIENTS': [4,4+self.frag.natoms-1,1,5,7,1,1,1.0]},
                          'ecp_grads':{ 'ENERGY GRADIENTS': [4+self.frag.natoms,4+self.frag.natoms+self.necps-1,1,5,7,1,1,1.0]},
                        }

        # YL 14/12/2023: multistate calculation
        # KS - 9/2/24 to accommodate open-shell systems
        if self.mult == 1:
            for i in range(self.nstates):
                # typically: Root   1 singlet a              0.283365076 a.u.                7.7108 eV
                # internally we start from 0 but the NWChem labels start from 1
                self._outkeys.update({'excitation {}'.format(i):{f'Root{i+1:4d}': [0,0,1,2,6,1,0,1.0]}})
                self._outkeys.update({'oscillator strength {}'.format(i):{f'Root{i+1:4d}': [8,8,1,3,3,1,0,1.0]}})
        else:
            for i in range(self.nstates):
                # typically: Root   1 singlet a              0.283365076 a.u.                7.7108 eV
                # internally we start from 0 but the NWChem labels start from 1
                self._outkeys.update({'excitation {}'.format(i):{f'Root{i+1:4d}': [0,0,1,2,5,1,0,1.0]}})
                self._outkeys.update({'oscillator strength {}'.format(i):{f'Root{i+1:4d}': [9,9,1,3,3,1,0,1.0]}})
          

    def parseMainOutput(self, librun=False):
        '''Parse main output file self.output'''

        from ...io   import file
        from ...data import unitconvert

        self.setOutKeys()

        # only for system call
        if not librun:
            self._result.energy = file.getArrayFromFile(self.output, 'out', self._result.energy.shape, regex=self._outkeys['energy'])
            # in excited state calculation we use E0 for ground state energy and result.energy for the state of interest
            if self._excited_state:
                self._result.E0 = self._result.energy[0]

        # excited state calculations
        for i in range(self.nstates):
            abuff = file.getArrayFromFile(self.output,
                                         'out',
                                          self._result.energy.shape,
                                          regex=self._outkeys[f'excitation {i}'],
                                          number_only=False)

            # KS - 20/02/24 fix the index for self._result.excitation[i] 
            #               since NWCHEM is inconsistent with printing spin state with roots, adding in the ground spin state
            spinstates = {2:'doublet', 3:'triplet', 4:'quartet', 5:'quintet', 6:'sextet'}
            if len(abuff) == 4:
                # KS TODO - 20/02/24: in case of open shell we are unsure yet what spin of a excited state NWChem outputs; we should take <S2> values in the result object in the future
#                abuff.insert(0, spinstates[self.mult])
                abuff.insert(0, ' ')

            # oscillator strength (f) will be determined next
            if not librun:
                self._result.excitations[i] = self._result.E0+float(abuff[2]),    \
                                              0.0,                                \
                                              unitconvert.eV2nm(float(abuff[4])), \
                                              abuff[1].lower(),                   \
                                              abuff[0].lower()
            # do not overwrite the excitation energies passed from rtdb
            else:
                self._result.excitations[i] = self._result.excitations['energies'][i], \
                                              0.0,                                     \
                                              unitconvert.eV2nm(float(abuff[4])),      \
                                              abuff[1].lower(),                        \
                                              abuff[0].lower()
            f = file.getArrayFromFile(self.output,
                                     'out',
                                      self._result.energy.shape,
                                      regex=self._outkeys[f'oscillator strength {i}'],
                                      number_only=True)
            self._result.excitations['f'][i] = f[0]

            # state 0: ground state; target_state counts from the ground state!
            # KS FIXME - 20/02/24: is this a duplicate of same lines in runLib()?
            if self.target_state:
                self._result.energy = self._result.excitations['energies'][self.target_state-1]
            else:
                self._result.energy = self._result.energy[0]
                
        try:
            if self._gradients:
                self._result.gradients = file.getArrayFromFile(self.output, 'out', self._result.gradients.shape, regex=self._outkeys['gradients'])
                if self.necps and self.separate_ecps:
                    self.ecps.gradients = file.getArrayFromFile(self.output, 'out', self.ecps.gradients.shape, regex=self._outkeys['ecp_grads'])
        except:
            print("\n >>> ChemShell NWChem: failed parsing gradients from output:\n")
            raise


    def getInpBuff(self):
        '''String buff containing the input information'''

        from numpy       import arange, array, char, __version__
        # TD 18/06/2024 : Numpy v2.0.0 version dependence
        from numpy.lib import NumpyVersion
        if NumpyVersion(__version__) >= '2.0.0b1':
            from numpy._core  import defchararray, full
        else:
            from numpy.core  import defchararray, full
        from .fragmovecs import getReorderedAtoms, getFragments, getSCFs, getPreSCF
        from ...data     import basis
        from ...objects  import fragment, fragbio

        # method (default: scf)
        self.method = method = self._synons.get(self.method.lower().strip(), 'scf')

        # headers
        strbuff = 'title \"ChemShell NWCHEM interface\"' +'\n'

        # memory management
        if self.memory != self.mem_stack + self.mem_heap + self.mem_global:
            self.mem_global = self.memory/2
            self.mem_stack  = self.memory/4
            self.mem_heap   = self.memory/4
            strbuff += 'memory stack %d mb heap %d mb global %d mb verify\n'%(self.mem_stack, self.mem_heap, self.mem_global)
        else:
            strbuff += 'memory ' + str(self.memory) + ' mb verify\n'

        # print level
        strbuff += 'print %s\n'%self.print_level

        # new or a restart
        if self.restart:
            # gradient to switch off fragmovecs in a restart job
            self.fragmovecs = ''
            strbuff += 'restart _nwchem\n'
        else:
            strbuff += 'start _nwchem\n'

        # coords header
        strbuff += 'geometry {} print units bohr nocenter noautoz\n'.format(self.frag.formula)

        # symmetry 
        if not self.symmetry:
            strbuff += '    symmetry c1 tol 0.0\n'
        else:
            strbuff += f'    symmetry {self.symmetry}\n'

        # coords
# YL FIXME: unsafe! don't use util we have to rearrange mass
        # ghost atoms
        original_names = self.frag.names[self.ghost].copy()
        for g in self.ghost:
            self.frag.names[g] = self.ghost_prefix.encode() + self.frag.names[g]
        if self.fragmovecs:
            strbuff += getReorderedAtoms(self)
        else:
# YL set a default
            if isinstance(self.frag, fragbio.FragmentBio):
                strbuff += self.frag.coords2str(atoms='symbols', charges=None, coords=14)
            else:
                strbuff += self.frag.coords2str(atoms='names', charges=None, coords=14)
        # ghost atoms
            self.frag.names[self.ghost] = original_names
# END YL FIXME


        # separate boundary ECPs (if qmmm.separate_ecps)
        if self.ecps:
            # prefix with 'bq_' and suffix with '_e'
            ecpnames = defchararray.add(full((self.necps,), 'bq_', dtype='S3'), self.ecps.names)
            ecpnames = defchararray.add(ecpnames, full((self.ecps.natoms,), '_e', dtype='S2'))
            # charges have to be scientific notation!!! otherwise NWChem does not work!!!
            strbuff += self.ecps.coords2str(coords=14,
                                            charges=1,
                                            atoms=None,
                                            indent=4*' ',
                                            delimiter=' '*1,
                                            notations={'q':'e'},
                                            inserts=[(0,ecpnames),(4,'charge')])

        strbuff += 'end\n'

        # YL 08/04/2019: superposition of fragment molecular orbitals
        #                see: http://www.nwchem-sw.org/index.php/Release66:Hartree-Fock_Theory_for_Molecules#Superposition_of_fragment_molecular_orbitals
        if self.fragmovecs:
            strbuff += getFragments(self)

        # default max num of bqs is only 25000
        if self.frag.nbqs > 25000:
            nbqs = self.frag.nbqs
            strbuff += 'set bq:max_nbq %d\n'%((int(str(nbqs)[0])+1)*10**(len(str(nbqs))-1))

        # background point charges
        if self.frag.nbqs:

            # region 2 BQ charges subtracted with self._ecp_charge
            if self.ecps:
                # YL 19/11/2019: ecps should exclude interface region shells but it doesn't really matter since the charges are neglectable
                # YL 13/09/2020: removed the general method select() from <class BQs>
                ecprange = self.frag.bqs.selectByRegions(self.ecp_region)
                self.frag.bqs.charges[ecprange] -= self._ecp_charge

            #Old chemshell had this line below
            #strbuff += 'end \n'
            strbuff += 'bq units au\n'
            if self._gradients:
                strbuff +=f'    # ChemShell: please note the signs in {self.bq_grads} because it is \"gradient\" but not \"force\" dumped\n'
                strbuff += '    force %s\n'%self.bq_grads
            #AJL, April 2018. Introduced higher accuracy for charges as found otherwise didn't get
            #correct agreement with NWChem embedded cluster calculations in TCL. And this accuracy is free!
            strbuff += self.frag.bqs.coords2str(coords=14, charges=8, atoms=None, indent=' ')

            strbuff += 'end\n'

            # restore the region 2 BQ charges
            if self.ecps:
                self.frag.bqs.charges[ecprange] += self._ecp_charge
  
        # basis set
        # write basis using library keyword if new job and do nothing if restart
        if not self.restart:

            # YL: use just a keyword `basis` to replace all `basis`, `basisfile`, `basisdeck`, `basisspec` in Tcl-ChemShell
            # AJL, Aug 2017: This doesn't work if a user wants to control the basis definition though? Perhaps a suitable workaround is
            # to do a basis file check, and if doesn't exist then to define it all internally?

            if self.harmonic:
                spherical = 'spherical '
            else:
                spherical = ''

            if method.endswith('dft'):
                # AJL, Aug2017. Introduced a local variable so we can edit the basis headline
                #strbuff += 'basis print\n'
                basis_headline = 'basis {}print'.format(spherical)
            else:
                #strbuff += 'basis \"ao basis\" print\n'
                basis_headline = 'basis \"ao basis\" {}print'.format(spherical)

            # basis set name
            if not '\n' in self.basis:
                strbuff += basis_headline + '\n'
                strbuff += '    * library ' + self.basis + '\n'
                strbuff += 'end\n'
            # content of basis set file or basis set deck
            else:
                # AJL, Aug 2017
                # Made this consistent with Gamess-UK, so we can define the basis settings (e.g. spherical)
                # YL 28.02.2019: improved 
                # leading "basis" headline
# YL FIXME: don't use
                if len(self.ghost) > 0:
                    basisbuff = self.getBasisWithGhost()
                else:
                    basisbuff = self.basis
# END YL FIXME

                # YL 21/02/2021: refactored with numpy.char
                abuff = array(basisbuff.splitlines())
                # comment out all lines starting with 'basis' and 'end'
                masks_basis = char.startswith(char.lower(char.strip(abuff)), 'basis')
                masks_end   = char.startswith(char.lower(char.strip(abuff)), 'end')
                abuff[masks_basis+masks_end] = char.add('#', abuff[masks_basis+masks_end])
                # attach a headline
                strbuff += f'\n{basis_headline}\n'
                strbuff += '\n'.join(abuff)
                # attach an "end"
                strbuff += '\nend\n'

#                if not basisbuff.strip().lower().startswith('basis'):
#                    strbuff += basis_headline + '\n'
#                strbuff += basisbuff + '\n\n'
#                # trailing "end"
#                if not basisbuff.lower().strip().split('\n')[-1].endswith('end'):
#                    strbuff += 'end\n\n'

            # ecp
            strbuff += utils.fileutils.file2text(self.ecp, prefix=parallel.prefix) + '\n'

        # guess file
        strbuff += utils.fileutils.file2text(self.guess, prefix=parallel.prefix) + '\n'

        # YL 08/04/2019: superposition of fragment molecular orbitals
        #                see: http://www.nwchem-sw.org/index.php/Release66:Hartree-Fock_Theory_for_Molecules#Superposition_of_fragment_molecular_orbitals
        if self.fragmovecs:
            strbuff += getSCFs(self)

        strbuff += '\ncharge {}\n'.format(self.charge)
        strbuff += 'set geometry {}\n'.format(self.frag.formula)

        # YL 14/12/2023: TD-DFT has its own grad keyword
        # KS 20/02/24: brought up and will be used below
        selectcases = {  True   :'gradient',
                        'shells':'gradient' }
        kw_gradients = selectcases.get(self._gradients, '')

        # SCF
        if method == 'scf':
            strbuff += self.__getSCFBlock()

        # DFT
        elif method == 'dft':
            strbuff += self.__getDFTBlock()

        # YL 14/12/2023: TD-DFT
        elif method == 'tddft':
            strbuff += self.__getDFTBlock() # added on suggestion by TK 4/1/24
            # KS - 20/02/24: we keep this feature for a future use case for example that during an MD traject on the ground state
            #                with vertical excitation energies calculated (for simulating absorption spectra)
            if self._gradients and not self.target_state:
                # task: always calculate energy
                strbuff += f'task dft energy {kw_gradients}\n'
            strbuff += self.__getTDDFTBlock()

        # other methods
        else:
            print(' >>> FATAL ERROR: unknown NWChem method!')
            quit(1)

        # set the number of linearly dependent eigenvectors
        if self.lindep_n >= 0:
            strbuff += 'set lindep:n_dep     {}\n'.format(self.lindep_n)

        # properties
        if self.properties:
            strbuff += 'property\n'
            for i in utils.iterutils.getFlattened([self.properties]):
                strbuff += '    {}\n'.format(i)
            strbuff += 'end\n'

        # KS - 20/02/24: in case we want the ground-state gradients but TD-DFT energies at the same time
        #                the gradient task for the TD-DFT block must be turned off
        if self.method == 'tddft' and not self.target_state:
            kw_gradients = ''
        # task: always calculate energy
        strbuff += f'task {method} energy {kw_gradients}'
        if self.properties:
            strbuff += ' property'
        strbuff += '\n'

        return strbuff


    def parseAdditionalOutput(self):
        '''Overload: Parse file for BQ gradients'''

        from ...io import file

        if self.frag.nbqs and self._gradients:
            self.frag.bqs.gradients = file.getArrayFromFile(self.bq_grads, 'nwchem_bq', self.frag.bqs.gradients.shape)


    def runLib(self, **kwargs):
        '''Run a library call of NWChem (called by run.runLib())'''

        from ctypes     import c_bool, c_int, byref
        from sys        import stdout
        from ...dl_py2f import py2f

        # overload self.parseOutput()
        self.parseOutput = lambda : 0

        print("\n >>> running NWChem as a compiled library...")

        # load NWChem library
        libnwchem = utils.modutils.getLinkedLib(self.__class__.__name__, chemsh.CHEMSH_LIB_DIR)

        theoryRef = py2f(self)
        resultRef = py2f(self._result)

        if libnwchem.interface_nwchem(byref(theoryRef),
                                      byref(resultRef),
                                      byref(c_bool(self._gradients)),
                                      byref(c_int(chemsh.CHEMSH_DEBUG))):

            print(" >>> NWChem: libnwchem.interface_nwchem failed.")
            return 1

        else:
            if self.method == 'tddft':
                self._result.E0 = self._result.energy[0]
                # KS - modify prining traget energy state when calculatin only VE
                if self.target_state:
                    self._result.energy = self._result.excitations['energies'][self.target_state-1]
                else:
                    self._result.energy = self._result.energy[0]

                # YL 16/12/2023: because some information is missing in rtdb, we have to parse the output file to complement the result
                self.parseMainOutput(librun=True)
            return 0


    def __getDFTBlock(self):
        ''''''

        dftstrbuff = ''

        if self.pre_hf:
            dftstrbuff += self.__getSCFBlock()

        # pre-SCF: since in NWChem DFT doesn't support superposition of fragment MOs
        if self.fragmovecs:
            dftstrbuff += self.__getSCFBlock()

        dftstrbuff += '\ndft\n'

        # HNT, Jul 2020. Adds population analysis keyword. 
        if self.pop_analysis:
             dftstrbuff += self.pop_analysis + '\n'
        else: 
             dftstrbuff += ''
        
        # AJL, June 2018. This limits the user to only pre-defined XC functionals
        # Changing to allow use of XC functionals not defined in synons
        # The default for functional is BLYP anyway, so this doesn't change default functionality.
        #dftstrbuff += '    xc ' + self._synons.get(self.functional, 'becke88 lyp') + '\n'
        dftstrbuff += '    xc ' + self._synons.get(self.functional, self.functional) + '\n'

# KS WIP 20/02/24 - define some range separated functionals useful for excited states

#CAM-B3LYP:
#xc xcamb88 1.00 lyp 0.81 vwn_5 0.19 hfexch 1.00
#cam 0.33 cam_alpha 0.19 cam_beta 0.46

#LC-BLYP:
#xc xcamb88 1.00 lyp 1.0 hfexch 1.00
#cam 0.33 cam_alpha 0.0 cam_beta 1.0

#LC-PBE:
#xc xcampbe96 1.0 cpbe96 1.0 HFexch 1.0
#cam 0.30 cam_alpha 0.0 cam_beta 1.0

#LC-PBE0 or CAM-PBE0:
#xc xcampbe96 1.0 cpbe96 1.0 HFexch 1.0
#cam 0.30 cam_alpha 0.25 cam_beta 0.75





####################

        # grid options
        selectcases = { 'high'    :'fine',
                        'veryhigh':'xfine' }
        if self.grid:
            dftstrbuff += '    grid ' + selectcases.get(self.grid, self.grid) + ' nodisk\n'
        else:
            dftstrbuff += '    grid nodisk\n'

        # the inverted charge-density and exchange-correlation matrices for a DFT calculation are normally written to disk storage.
        # the user can prevent this by specifying the keyword noio within the input for the DFT directive.
        dftstrbuff += '    noio\n'

#### adding dispersion functions 
# Q. Wu & W. Yang, J. Chem. Phys. 116 515 (2002) and
# U. Zimmerli, M Parrinello and P. Koumoutsakos J. Chem. Phys. 120 2693 (2004)
        if self.dispersion.lower() == 'd0':
           dftstrbuff += '    disp vdw 0\n'
# Grimme's D1: S. Grimme J. Comp. Chem. 25 1463 (2004)
        elif self.dispersion.lower() == 'd1':
             dftstrbuff += '    disp vdw 1\n'
# Grimme's D2: S. Grimme J. Comp. Chem. 271787 (2006)
        elif self.dispersion.lower() == 'd2':
             dftstrbuff += '    disp vdw 2\n'
# Grimme's D3: S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 132, 154104 (2010)
        elif self.dispersion.lower() == 'd3':
             dftstrbuff += '    disp vdw 3\n'
 # D3BJ: only BLYP, B3LYP, BHLYP, TPSS, and TPSSh are supported
        elif self.dispersion.lower() == 'd3bj':
             dftstrbuff += '    disp vdw 4\n'

######

        if self.direct:
            dftstrbuff += '    direct\n'
        else:
            dftstrbuff += '    incore\n'


        # iterations
        dftstrbuff += 4*' ' + 'maxiter ' + str(self.maxiter) +'\n'

        # convergence
        # tolerance of energy, density, and orbital gradient
        dftstrbuff +=  4*' ' + 'convergence '
        dftstrbuff +=          'energy   %.1E \\\n'%self.tol_energy
        dftstrbuff += 16*' ' + 'density  %.1E \\\n'%self.tol_density
        dftstrbuff += 16*' ' + 'gradient %.1E \\\n'%self.tol_gradient
        # damping
        if self.damping:
            dftstrbuff += 16*' ' + 'damp ' + str(self.ndamp) + ' \\\n'
        else:
            dftstrbuff += 16*' ' + 'nodamp \\\n'
        # diis
        if self.diis:
            dftstrbuff += 16*' ' + 'diis ' + str(self.nfock) +' \\\n'
        else:
            dftstrbuff += 16*' ' + 'nodiis \\\n'
        # level shifting
        if self.lshift:
            dftstrbuff += 16*' ' + 'lshift ' + str(self.lshift) + ' \\\n'
        else:
            dftstrbuff += 16*' ' + 'nolevelshifting \\\n'
        # more convergence helpers
        dftstrbuff += 16*' ' + 'dampon  ' + str(self.dampon ) + ' \\\n'
        dftstrbuff += 16*' ' + 'dampoff ' + str(self.dampoff) + ' \\\n'
        dftstrbuff += 16*' ' + 'diison  ' + str(self.diison ) + ' \\\n'
        dftstrbuff += 16*' ' + 'diisoff ' + str(self.diisoff) + ' \\\n'
        dftstrbuff += 16*' ' + 'levlon  ' + str(self.levlon ) + ' \\\n'
        dftstrbuff += 16*' ' + 'levloff ' + str(self.levloff) + ' \\\n'
        # iterations per procedure
        dftstrbuff += 16*' ' + 'ncydp ' + str(self.ncydp) + ' \\\n'
        dftstrbuff += 16*' ' + 'ncyds ' + str(self.ncyds) + ' \\\n'
        dftstrbuff += 16*' ' + 'ncysh ' + str(self.ncysh) + ' \\\n'
        # HOMO-LUMO gap
        dftstrbuff += 16*' ' + 'hl_tol ' + str(self.hl_tol) + '\n'

        # tolerances
        # YL 08/01/2018: removed accAoFunc as it is not in the official manual
        dftstrbuff += '    tolerances tol_rho %.1E \\\n               accCoul %d \\\n               radius %f\n' \
                  %(self.tol_rho, self.tol_acccoul, self.tol_radius)

        # projection of different basis sets
        if self.project:
            project = 'project \"{}\" '.format(self.project)
        else:
            project = ''
        # vectors print
        if self.integral_io:
              dftstrbuff += self.integral_io + '\n'

        if self.print_level.lower() in [ 'debug' ]:
            dftstrbuff += '    print "common" \\\n'                   # debug   dump of common blocks
            dftstrbuff += '          "io info"\n'                    # debug   reading from and writing to disk
        elif self.print_level.lower() in [ 'high' ]:
            dftstrbuff += '    print "all vector symmetries" \\\n'    # high    symmetries of all molecular orbitals
            dftstrbuff += '          "alpha partner info" \\\n'       # high    unpaired alpha orbital analysis
            dftstrbuff += '          "coulomb fit" \\\n'              # high    fitting electronic charge density
            dftstrbuff += '          "dft timings" \\\n'              # high    
            dftstrbuff += '          "final vectors" \\\n'            # high    
            dftstrbuff += '          "final vectors analysis" \\\n'   # high    print all orbital energies and orbitals
            dftstrbuff += '          "initial vectors" \\\n'          # high    
            dftstrbuff += '          "intermediate energy info" \\\n' # high    
            dftstrbuff += '          "intermediate evals" \\\n'       # high    intermediate orbital energies
            dftstrbuff += '          "intermediate fock matrix" \\\n' # high    
            dftstrbuff += '          "intermediate overlap" \\\n'     # high    overlaps between the alpha and beta sets
            dftstrbuff += '          "intermediate S2" \\\n'          # high    values of S2
            dftstrbuff += '          "intermediate vectors" \\\n'     # high    intermediate molecular orbitals
            dftstrbuff += '          "interm vector symm" \\\n'       # high    symmetries of intermediate orbitals
            dftstrbuff += '          "kinetic_energy" \\\n'           # high    kinetic energy      
            dftstrbuff += '          "quadrature" \\\n'               # high    numerical quadrature
            dftstrbuff += '          "schwarz" \\\n'                  # high    integral screening info & stats at completion
            dftstrbuff += '          "screening parameters"\n'       # high    integral accuracies
            if not self.pop_analysis == 'mulliken':                  # high    mulliken atomic orbital population
                dftstrbuff += '          "mulliken ao"\n'    
        elif self.print_level.lower() in [ 'medium', 'default' ]:
            dftstrbuff += '    print "convergence" \\\n'              # default convergence of SCF procedure
            dftstrbuff += '          "final vector symmetries" \\\n'  # default symmetries of final molecular orbitals
            dftstrbuff += '          "multipole" \\\n'                # default moments of alpha, beta, and nuclear charge densities
            dftstrbuff += '          "parameters" \\\n'               # default input parameters
            dftstrbuff += '          "semi-direct info"\n'           # default semi direct algorithm   
        elif self.print_level.lower() in [ 'medium', 'default' ]:
            dftstrbuff += '    print "information"\n'                # low     general information
        elif self.print_level.lower() in [ 'no', 'none' ]:
            dftstrbuff += '    noprint\n'                            # noprint 
        
        pamethods = {'mulliken':'"mulliken ao"', 'fukui':'"Fukui information"'}
        if self.pop_analysis:
            dftstrbuff = dftstrbuff.rstrip() + ' \\\n'
            dftstrbuff += '          ' +  pamethods.get(self.pop_analysis, self.pop_analysis) + '\n'

        # .movecs
        vecs_out = 'output %s'%self.vecs_out
        if self.fragmovecs:
            vecs_in = '_{}_prescf.movecs'.format(self.frag.formula)
            # for future restart
            self.vecs_in = self.vecs_out
        else:
            if self.pre_hf:
                vecs_in = 'input _nwchem_pre_hf.movecs'
            else:
                if self.vecs_in:
                    vecs_in = '{}{}'.format(project, self.vecs_in)
                else:
                    vecs_in = ''
        dftstrbuff += 4*' ' + 'vectors %s %s\n'%(vecs_in, vecs_out)

        # SCF type parsing
        if self.scftype.lower() in [ 'uhf', 'uks' ]:
            dftstrbuff += '    odft\n'
        elif self.scftype.lower() in [ 'rohf', 'roks' ]:
            dftstrbuff += '    rodft\n'

        # end with multiplicity
        dftstrbuff += '    mult ' + str(self.mult) + '\n'
        dftstrbuff += 'end\n'
    
        # fix to prevent mirroring messing up DFT when using GA
        dftstrbuff += 'set dft:xcreplicated false\n'
        dftstrbuff += 'set dft:converged    false\n'
        # set ScaLAPACK to 'x' or 'v', the default 'd' does NOT work!
        dftstrbuff += 'set dft:scaladiag    x\n'
#        dftstrbuff += 'set dft:lrt_tddft    true\n'
#        dftstrbuff += 'set qmmm:init        true\n'
#        dftstrbuff += 'set dft:bqdontcare   true\n'

        return dftstrbuff

    def __getTDDFTBlock(self):
        ''''''

        # method: NWChem defaults to RPA (TD-DFT), otherwise CIS (Tamm-Dancoff TD-DFT)
        rpa = {True:'rpa',False:'cis'}.get(self.rpa, 'rpa')
        # gradients: target_state counts from the ground state (state 0)
        if self._gradients:
            if self.target_state:
                grad = f'''    civecs
        grad
            root {self.target_state}
        end'''
            # KS - 20/02/24: if we want the gradients on the ground state we can't do it here in the TD-DFT block
            else:
                grad = f'''civecs
        '''
        else:
            grad = ''
        # spin
        # KS TODO - 20/02/24: in case of open shell we are unsure yet what spin of a excited state NWChem outputs; we should take <S2> values in the result object in the future
        if self.mult == 1:
            if self.spin.strip().lower().startswith('sing'):
                spin = 'singlet\n    notriplet'
            elif self.spin.strip().lower().startswith('trip'):
                spin = 'triplet\n    nosinglet'
        else:
            spin = ''
        # symmetry
        if self.symmetry:
            symmetry = f'symmetry\n    targetsym none\n'
        else:
            symmetry = ''

        sbuff = f'''
tddft
    {rpa}
    nroots {self.nstates}
    maxvecs 1000
    {spin}
    thresh  {self.tol_energy}
    maxiter {self.maxiter}
    {symmetry}
#    ecut] <-cutoff energy>  
#    ewin] <-lower cutoff energy>  <-higher cutoff energy>  
#    alpha] <integer lower orbital>  <integer upper orbital>  
#    beta] <integer lower orbital>  <integer upper orbital>  
    {grad}
#    cdspectrum
#    giao
#    velocity
#    simpleso
    algorithm 0
#    freeze [[core] (atomic || <integer nfzc default 0>)]
#            [virtual <integer nfzv default 0>]]  
    print {self.print_level}  \\
        "timings"             \\
#        "trial vectors"       \\
        "initial guess"       \\ 
        "general information" \\
        "xc information"      \\
        "memory information"  \\
        "convergence"         \\
        "subspace"            \\
        "transform"           \\
        "diagonalization"     \\
        "iteration"           \\
        "contract"            \\
        "ground state"        \\
        "excited state"
end
'''
        return sbuff


    # YL 15/04/2019: modified
    def __getSCFBlock(self):
        '''Generate an SCF block'''

        scfstrbuff = '\nscf\n'
        scfstrbuff += 4*' ' + 'maxiter %d\n'%self.maxiter
        if self.direct:
            scfstrbuff += 4*' ' + 'direct\n'
        scfstrbuff +=  4*' ' + 'print "information"            \\\n' # low     general information
        scfstrbuff += 10*' ' + '"mo guess"               \\\n' # 
        scfstrbuff += 10*' ' + '"final vectors analysis" \\\n' # 
        scfstrbuff += 10*' ' + '"final evals"            \\\n' # 
        scfstrbuff += 10*' ' + '"vectors i/o"            \\\n' # 
        scfstrbuff += 10*' ' + '"parameters"             \\\n' # 
        scfstrbuff += 10*' ' + '"convergence"              \n' # 
        if self.fragmovecs or self.pre_hf:
            scfstrbuff += 4*' ' + 'thresh {}\n'.format(self.tol_pre_hf)  #
        else:
            scfstrbuff += 4*' ' + 'thresh {}\n'.format(self.threshold)  #

        if self.mult == 1:
            scfstrbuff += '    singlet\n'
            scfstrbuff += '    {}\n'.format(self.scftype)
        elif self.mult == 2:
            scfstrbuff += '    doublet\n'
            scfstrbuff += '    {}\n'.format(self.scftype)
        elif self.mult == 3:
            scfstrbuff += '    triplet\n'
            scfstrbuff += '    {}\n'.format(self.scftype)

        # superposition of fragment MOs
        if self.fragmovecs:
            from .fragmovecs import getPreSCF

            scfstrbuff += getPreSCF(self)

        else:
            # .movecs
            if self.vecs_in and self.pre_hf:
                vecs_in  = 'input %s'%self.vecs_in
                vecs_out = 'output _nwchem_pre_hf.movecs'
            elif self.vecs_in and not self.pre_hf:
                vecs_in  = 'input project \"{}\" {}'.format(self.project, self.vecs_in)
                vecs_out = ''
            elif not self.vecs_in and self.pre_hf:
                vecs_in  = ''
                vecs_out = 'output _nwchem_pre_hf.movecs'
            else:
                vecs_in  = ''
                vecs_out = 'output %s'%self.vecs_out
            scfstrbuff += 4*' ' + 'vectors {} {}\n'.format(vecs_in, vecs_out)

        # Open shells
        nopen = self.mult - 1
        if nopen > 0:
            scfstrbuff += '    nopen ' + str(nopen) + '\n'
       
        scfstrbuff += 'end\n'

        # set the number of linearly dependent eigenvectors
        if self.lindep_n >= 0:
            scfstrbuff += 'set lindep:n_dep     {}\n'.format(self.lindep_n)

        if self.pre_hf or self.fragmovecs:
            scfstrbuff += 'task scf energy\n'

        return scfstrbuff







