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

# david.gunn@stfc.ac.uk

import re
import os.path
import math
import chemsh
from numpy       import full, zeros
from ..qm        import _QM
from ...base     import errors
from ...dl_py2f  import utils
from ...data     import atomicdata, unitconvert
from ...parallel import parallel

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

    _kwargs = {
                'auxbasis'              :'',
                'basisspec'             : False,
                # YL 09/03/2021: (suggested by HMS) with ORCA it's always encouraged to use "TightSCF" although "NormalSCF" is the default
                # YL NB 09/03/2021: but we may change the option's keyword in the future to make it uniform for all interfaces (currently other interfaces use "threshold=1.0e-6")
#                'convergence'           :"TightSCF", # YL 16/02/2024: refactored to meet the Py-ChemShell standard
                'cosmo'                 : None,
                'd3'                    : False,
                'enerflag'              : 1,
                'energy_correction'     : 0.0,
                'eroots'                : 1,
                'estate'                : 1,
#                'excited'               : False,     # YL 16/02/2024: refactored to meet the Py-ChemShell standard
                'explicit'              : False,
                'exp_list'              : [],
                'filename_bq'           :'_orca.bq',
                'finalgridsize'         : None,
                'gridsize'              : None,
                'guess'                 :'hueckel',
                'hamout'                :'HF',
                'image'                 : None,
                'input'                 :"_orca.inp",
                'jobname'               :"_orca",
                'list_option'           : "medium",
                'listing_index'         : 1,
#                'moinp'                 : "save.gbw", # YL 16/02/2024: refactored to meet the Py-ChemShell standard
                'newgto'                : None,
                'nprocs'                : 0, # allow for user override, see note under parallel execution
                'num_explicit'          : 0,
#                'old_moinp'             :"_orca.gbw", # YL 16/02/2024: refactored to meet the Py-ChemShell standard
                'optstr'                :'',
                'output'                :'_orca.out',
                'output_bqs'            :"_orca.pcgrad",
                'output_eandg'          :"_orca.engrad",
                'path'                  :'orca',
                'simple_input'          : True,
                'stddft_ethresh'        : 10.0,       # energy threshold up to which CSFs are included (in eV)
                'stddft_pthresh'        : 1.0e-04,    # selection criterion to include CSF beyond EThresh (in Eh)
                'stddft_ptlimit'        : 30,         # energy threshold up to which CSFs beyond EThresh may be selected (in eV)
                'tda'                   : True,
                'threshold'             : 1.0e-08,    # YL 16/02/2024: defaulted to TightSCF
                'use_finalgrid'         : True,
#                'use_ri'                : False,     # YL 16/02/2024: refactored to meet the Py-ChemShell standard
                 # v4.0 changed some basis set input conventions
                'version'               : 4.0,
                'vecs_inp'              :'_orca_restarted.gbw',
              }

    _internals = {
                   '_basis_element' : [],
                   '_basis_list'    : '_orca.engrad',
                   '_basis_sym'     : [],
                   '_ecp_list'      : '',
                   '_ecp_ncore'     : [],
                   '_ecp_sym'       : [],
                   '_extra_outputs' : [''],
                   '_exec_run_shell': True,
                   '_syscmd'        : 'orca',
                 }

    _synons = {
                'blyp'   : 'becke88 lyp',
                'lda'    : 'slater vwn_5',
                's-vwn'  : 'slater vwn_5',
                'hf'     : 'hf',
                'dft'    : 'dft',
                'mp2'    : 'mp2',
                'rimp2'  : 'ri-mp2',
                'ri-mp2' : 'ri-mp2',
                'stddft' : 'stddft',
                'std-dft': 'stddft',
                'tddft'  : 'tddft',
                'td-dft' : 'tddft',
              }

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

        if self.frag.bqs:
            self.frag.nbqs = len(self.frag.bqs.coords)
            ncent = self.frag.natoms + self.frag.nbqs

        method = self._synons.get(self.method.lower())
        self.hamout = method

        strbuff = '#ORCA input file generated by ChemShell\n'
        if self.d3:
            if self.ri:
                riout = "RI"
            else:
                riout = "NoRI"
            if self.hamout == 'hf':
                strbuff += "! {} HF {} D3\n".format(self._precision, riout)
            if self.hamout == 'dft':
                strbuff += "! %s %s %s D3\n" % (self._precision, self.functional, riout)
        else:
             strbuff += "! %s\n" % (self._precision)

        if self.cosmo:
            strbuff += "! COSMO(%s)\n" % (self.cosmo)

        #Parallel execution
        # TWK 15/5/2024 Re-enable user control of parallelisation of ORCA
        # to allow for the use case of serial ChemShell + parallel ORCA
        if self.nprocs > 0:
            input_nprocs = self.nprocs
        else:
            input_nprocs = parallel.nprocs_per_wg

        if input_nprocs > 1:
            strbuff +=  '# parallel\n'
            strbuff +=  '%pal\n'
            strbuff += f'    nprocs {input_nprocs}\n'
            strbuff +=  'end\n'

        #Specify the point charge file
        if self.frag.nbqs > 0:
            strbuff += "%pointcharges \"pointcharges.xyz\" \n"

        #Method
        strbuff += "%method\n"
        strbuff += "  Method %s\n" % (self.hamout)
        if self.hamout == 'dft':
            strbuff += "  Functional %s\n" % (self.functional)
            if self.gridsize:
                strbuff += "  Grid %s\n" % (self.gridsize)
            if self.use_finalgrid:
                if self.finalgridsize:
                    strbuff += "  FinalGrid %s\n" % (self.finalgridsize)
            else:
              strbuff += "  UseFinalGrid false\n"

        if self.ri:
            strbuff += "  RI on\n"
        else:
            strbuff += "  RI off\n"

        if self.enerflag:
            if self._gradients:
                strbuff += "  RunTyp Gradient\n"
            else:
                strbuff += "  RunTyp Energy\n"
        strbuff += "end\n"

        #Basis set
#         Determine if self.basis contains one or more basis, or if file type with newGTO info
#         Individual newGTO Basis info written after the atom in the 'coords' section
#hms      The auxbasis is not handled correctly below. There is another "end" keyword inserted into the %basis block, which is nonsense.
#hms      Instead, we insert "end" only at the very end, after both the basis and auxbasis have been processed.
        words = self.basis.split()
#         case 1: one item in self.basis with global basis info 
        if (len(words) == 1):
          strbuff += "%basis\n"
          if self.version < 4.0:
            strbuff += "  Basis %s\n" % (self.basis)
          else:
            strbuff += "  Basis \"%s\" \n" % (self.basis)
#hms         strbuff += "end\n"
#         case 2: more than one item in self.basis but not read in from file
        if (len(words) > 1) and "\"" not in self.basis.lower():
          strbuff += "%basis\n"
          strbuff += "  %s\n" % (self.basis)
#hms          strbuff += "end\n"
#         case 3: basis set in file
        if (len(words) > 1) and "\"" in self.basis.lower():
          self.basisspec = True
          lines = self.basis.splitlines()
          num_basis = 0
          self._basis_sym = [None] * (self.frag.natoms+1) 
          self._basis_element = [None] * (self.frag.natoms+1)
          self._basis_list = [None] * (self.frag.natoms+1)
          for line in lines:
            if "\"" in line:
              num_basis += 1
              word = line.split("\"")
              self._basis_sym[num_basis] = word[1]
              self._basis_element[num_basis] = ''.join([k for k in word[1] if not k.isdigit()])
              self._basis_list[num_basis] = "newGTO \n"
            else:
              if (num_basis > 0):
                if ("newgto" not in line.lower()):
                    self._basis_list[num_basis] +="%s" % (line.rstrip())
                if ("end" not in line.lower()):
                    self._basis_list[num_basis] +="\n"
                
        #ECPs
        if self.ecp:
            num_ecp = 0
            j = ''
            self._ecp_ncore = [None] * (self.frag.natoms+1)
            self._ecp_sym = [None] * (self.frag.natoms+1)
            self._ecp_list = [None] * (self.frag.natoms+1)
            if os.path.isfile(os.path.abspath(self.ecp)):
#               ECP is in a file
                ecp_list = utils.fileutils.file2text(self.ecp, prefix=parallel.prefix)
                lines = ecp_list.splitlines()
                for line in lines:
                  if "\"" in line:
                    num_ecp += 1
                    word = line.split("\"")
                    self._ecp_sym[num_ecp] = word[1]
                    j.join([k for k in word[1] if not k.isdigit()])
                    self._ecp_list[num_ecp] = "NewECP \n"
                  else:
                    if ("newecp" not in line.lower()):
                      self._ecp_list[num_ecp] +="%s\n" % (line)
                  if "core" in line.lower():
                    word = line.split()
                    self._ecp_ncore[num_ecp] = word[1]
            else:
#               ECP was defined in the input script
                lines = self.ecp.splitlines()
                for line in lines:
                    strbuff += line + '\n'
                    if "\"" in line:
                      num_ecp += 1
                      word = line.split("\"")
                      self._ecp_sym[num_ecp] = word[1]
                      j.join([k for k in word[1] if not k.isdigit()])
                      self._ecp_list[num_ecp] = "NewECP \n"
                    else:
                      if ("newecp" not in line.lower()):
                        self._ecp_list[num_ecp] +="%s\n" % (line)
                    if "core" in line.lower():
                      word = line.split()
                      self._ecp_ncore[num_ecp] = word[1]

        #Auxiliary basis for RI method
        if self.auxbasis:
            if self.version < 4.0:
                strbuff += "  Aux %s\n" % (self.auxbasis)
            else:
                strbuff += "  AuxJ \"%s\" \n" % (self.auxbasis)
        # YL/HMS 04/04/2021: this "end" is redundant
        #hms 11/04/2022: To the contrary, this is the only one needed.
        strbuff += "end\n"

        #SCF control
        strbuff += "%scf\n"
        strbuff += "  HFTyp %s\n" % (self.scftype)
        if self.restart:
            utils.fileutils.rename(self.old_moinp,self.moinp)
            strbuff += "  Guess MORead\n"
            strbuff += "  MOInp \"%s\" \n" % (self.moinp)
        else:
#           Patom and Pmodel seem to fail when defining own Basis or ECP so defaulting to hueckel
            if self.guess:
                strbuff += f'  Guess {self.guess}\n'
            strbuff += "  AutoStart false\n"
    
        if self.direct:
            strbuff += "  SCFMode Direct\n"
        else:
            strbuff += "  SCFMode Conventional\n"
    
        if self.maxiter:
            strbuff += "  MaxIter %s\n" % (self.maxiter)
    
        strbuff += "end\n"

        #Excited states
        if self._excited_state:
            if self.hamout == 'hf':
                strbuff += "%cis\n"
            elif self.hamout == 'dft':
                strbuff += "%tddft\n"
        
            # Must calculate as least as many roots as the state of interest
            if self.eroots > self.estate:
                strbuff += "  NRoots %s\n" % (self.eroots)
            else:
                strbuff += "  NRoots %s\n" % (self.estate)
            strbuff += "  IRoot %s\n" % (self.estate)
        
            if (self.hamout == 'dft') and not tda:
                strbuff += "  TDA false\n"
            strbuff += "end\n"

        #Geometry
        strbuff += "%coords\n"
        strbuff += "  CTyp xyz\n"
        strbuff += "  Mult %s\n" % (self.mult)
        strbuff += "  Units bohrs\n"
        strbuff += "  coords\n"

        #Point charge file
        if self.frag.nbqs:

            fp = open('pointcharges.xyz', 'w+')
            self.num_explicit = 0
            self.exp_list.clear() 
            num_file_bqs = 0
            ecprange = self.frag.bqs.getRegion(2)
            for i in range(self.frag.bqs.nbqs): 
#               This initial loop is to determine the number of bqs to be written in the pointcharges file, and removes the pointcharge(s)
#               that have to be treated differently (i.e. when requiring an 'explicit' interface)
                name = str(self.frag.bqs.names[i])
#               There's probably a neater way to strip the extraneous bits from the name but it'll do in a pinch
                name = name.replace(' ','')
                name = name.replace('\'','')
                name = name.replace('b','')
                coords = list(self.frag.bqs.coords[i])
                if (name not in self._ecp_sym):
                  num_file_bqs += 1
            fp.write(str(num_file_bqs))
            fp.write('\n')
            for i in range(self.frag.bqs.nbqs):
#               This loop writes the bq info into the pointcharges file.
                name = str(self.frag.bqs.names[i])
                name = name.replace(' ','')
                name = name.replace('\'','')
                name = name.replace('b','')
                coords = list(self.frag.bqs.coords[i])
                if (name not in self._ecp_sym):
#               Not an explicit ECP - added to pointcharges file
                  for j in range(3):
                      coords[j] = coords[j] / unitconvert.Angstrom2Bohr
                  fp.write('{0:3.10f} {1:3.10f} {2:3.10f} {3:3.10f}\n'.format(self.frag.bqs.charges[i], coords[0], coords[1], coords[2]))
                else:
#               Explicit ECP - added to main input file
                  element = '' 
                  element = element.join([k for k in name if not k.isdigit()])
                  atom_charge = self.frag.bqs.charges[i]
                  self.explicit = True
                  self.num_explicit += 1
                  self.exp_list.append(i)
                  strbuff += "%s> %s %s %s %s " % (element, atom_charge, coords[0], coords[1], coords[2])
                  for j in range(1,num_ecp+1):
                    if name in self._ecp_sym[j]:
                      strbuff += "\n%s" % (self._ecp_list[j])
            fp.close()

#           Due to the way the explicit ECPs are treated we need to remove the coulombic interaction between the explicit ECPs and the 
#           remaining bqs.
            self.energy_correction = 0.0
            test_e = 0.0
            for i in self.exp_list:
                coords_i = list(self.frag.bqs.coords[i])
                for j in range(self.frag.bqs.nbqs):
                  dist = 0.0
                  chargeprod = 0.0
                  coords_j = list(self.frag.bqs.coords[j])
                  chargeprod=-1.0 * self.frag.bqs.charges[i] * self.frag.bqs.charges[j]
                  for k in range(3):
                    dist += (coords_i[k] - coords_j[k])**2
                  dist = math.sqrt(dist) 
                  if dist and chargeprod:
                    if (j in self.exp_list) and (j > i):
                      test_e += chargeprod/dist
                    if j not in self.exp_list:
                      self.energy_correction += chargeprod / dist
            self.energy_correction += test_e

        for i in range(self.frag.natoms):
            coords = list(self.frag.coords[i])
            symbol = self.frag.names[i].decode('utf-8')#str(atomicdata.symbols[self.frag.znums[i]])
            element = str(atomicdata.symbols[self.frag.znums[i]])
#           Here we are checking if each element has an ECP, BASIS, both or neither and adding the appropriate info to strbuff
#           We need the differentiation as if ECP-only (treating as embedded shell) then the formatting is different
#           ORCA doesn't understand numbers after the element, so individual atom identification (for ECP and Basis assignment)
#           happens in the background with the 'symbol' variable.
            if self.ecp:
              shell_only = True 
              basis_only = True
              for j in range(1,num_ecp+1):
                if symbol in self._ecp_sym[j]:
                  basis_only = False
              if basis_only:
                strbuff += "%s %s %s %s " % (element, coords[0], coords[1], coords[2])
              else:
                if self.basisspec:
                  for k in range (1,num_basis+1):
                    if symbol in self._basis_sym[k]:
                      shell_only = False
                  if shell_only: #ECP only found for this atom (Implicit interface)
                      for j in range(1,num_ecp+1):
                        if symbol in self._ecp_sym[j]:
                            atom_charge = self.frag.znums[i] - int(self._ecp_ncore[j])
                            strbuff += "%s> %s %s %s %s " % (element, atom_charge, coords[0], coords[1], coords[2])
                  else: #ECPs used and basis found for this atom
                    strbuff += "%s %s %s %s " % (element, coords[0], coords[1], coords[2])
                else: #ECPs used but global basis used therefore not shell only
                  strbuff += "%s %s %s %s " % (element, coords[0], coords[1], coords[2])
            else: 
              strbuff += "%s %s %s %s " % (element, coords[0], coords[1], coords[2])
            if self.basisspec:
              for j in range(1,num_basis+1):
                if symbol in self._basis_sym[j]:
                  strbuff += "\n%s" % (self._basis_list[j])
            if self.ecp:
              for j in range(1,num_ecp+1):
                if symbol in self._ecp_sym[j]:
                  strbuff += "\n%s" % (self._ecp_list[j])
            strbuff += "\n"

        strbuff += "  end\n"
        strbuff += "  Charge %s\n" % (self.charge)
        strbuff += "end\n"

        #User-defined options
        if self.optstr:
            strbuff += "%s" % (self.optstr)

        return strbuff


    # YL 16/02/2024: refactored using the latest style of ORCA input
    def getInpBuff(self):
        '''String buff containing the input information'''

        from time import time as time_time
        from time import asctime, localtime

        # let's also support the old format deleveoped by DG which however
        # we don't guarantee to run for advanced functionalities
        # TODO Port all functionality from old to new style to avoid having two separate routines
        if not self.simple_input:
            return self.getInpBuffOLD()

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

        strbuff =   '# ORCA input file generated by ChemShell\n'
        strbuff += f'# {asctime(localtime(time_time()))}\n'
        if not self.freetext_only:
            strbuff +=  '\n! '

        # insert freetext of the current section
        if type(self.freetext) is dict:
            for section in [ '!' ]:
                strbuff += self.freetext.get(section, '').strip() + ' '

        # insert str type freetext
        if type(self.freetext) is str and self.freetext:
            strbuff += '\n# user-defined freetext input\n'
            strbuff += utils.fileutils.file2text(self.freetext).strip() + '\n'

        # get a comment-free strbuff
        strbuff_urtext = '\n'.join([ line.lower().strip().split('#')[0] for line in strbuff.splitlines() if not line.strip().startswith('#') ])

        # the following input is the ChemShell default but won't be used
        # if the user wants to replace with freetext (freetext_only=True)
        strbuff0 = ''

        # method
        if self.method == 'dft':
            strbuff0 += self.functional
        elif self.method == 'tddft' or self.method == 'stddft':
            strbuff0 += self.functional
        else:
            strbuff0 += self.method

        # basis set
        strbuff0 += f' {self.basis} {self.ri}'

        # TWK 24/6/24 Support dispersion option (e.g. D3ZERO, D3BJ)
        # TODO note the non-standard d3 option in the original ORCA interface 
        # can be removed after refactoring is completed 
        if self.dispersion:
            strbuff0 += f' {self.dispersion}'

        strbuff0 += '\n'

        # if no RI we must specify "! NORI"
        if not self.ri:
            strbuff0 += '\n! nori\n'

        # toggle gradients
        if self._gradients:
            strbuff0 += '! engrad\n'

        # printing
        if chemsh.CHEMSH_DEBUG:
            strbuff0 += '! largeprint\n'
            strbuff0 += '! printbasis\n'
            strbuff0 += '! printmos\n'
        else:
            strbuff0 += '! normalprint\n'
            strbuff0 += '! printgap\n'

        # TD-DFT
        if self.method in [ 'tddft', 'stddft' ]:
            section = '%tddft'
            strbuff0 +=  f'{section}\n'
            # insert freetext of the current section
            strbuff0 += self._getFreetext(section)
            if self.method == 'stddft':
                self.nstates = 0
                strbuff0 +=  '    mode stddft\n'
                strbuff0 += f'    ethresh {self.stddft_ethresh}\n'
                strbuff0 += f'    pthresh {self.stddft_pthresh}\n'
                strbuff0 += f'    ptlimit {self.stddft_ptlimit}\n'
            if not self.method == 'stddft':
                strbuff0 += f'    nroots {self.nstates}\n'
                strbuff0 += f'    iroot {self.target_state}\n'
            # in ORCA, singlets are calculated anyway
            if self.spin.strip().lower().startswith('t'):
                strbuff0 += '    triplets true\n'
            strbuff0 +=  'end\n'

        strbuff0 += '\n'

        # SCF control: adapted from DG's version
        section = '%scf'
        strbuff0 +=  f'{section}\n'
        # insert freetext of the current section
        strbuff0 += self._getFreetext(section)
        # numerical precision (suppressed if there is already a simple input keyword in freetext)
        if not any([ k in strbuff_urtext for k in self._precisions.values() ]):
            strbuff0 += f'    # the requested threshold is equivalent/close to ORCA simple input keyword {self._precision}\n'
            strbuff0 += f'    tole {self.threshold}\n'
        # 10/02/2024: "convcheckmode 0" (to check all) is default in ORCA for extremescf; we currently don't
        #             provide other options (1: any one of the criteria; 2: energy)
        if self.threshold <= 1.0E-12:
            strbuff0 += f'    convcheckmode 0\n'
        # post-HF methods need fully converged wavefunctions
        if 'mp2' in self.method:
            strbuff0 += '    convforced true\n'
        strbuff0 += f'    hftyp {self.scftype}\n'
        if not self.restart:
            strbuff0 += f'    guess {self.guess}\n'
            # we don't read from previous calculations by default
            strbuff0 +=  '    autostart false\n'
        if self.direct:
            strbuff0 +=  '    scfmode direct\n'
        else:
            strbuff0 +=  '    scfmode conventional\n'
        strbuff0 += f'    maxiter {self.maxiter}\n'
        strbuff0 +=  'end\n'

        # point charge file
        if self.frag.nbqs:
            with open(self.filename_bq, 'w') as fp:
                fp.write(f'{self.frag.nbqs}\n')
                fp.write(self.frag.bqs.coords2str(order='qc',
                                                  coords=14,
                                                  charges=14,
                                                  atoms=None,
                                                  unit='angstrom',
                                                  notations={'c':self._notation,'q':self._notation}))
            strbuff0 += '# background point charges for electrostatic embedding\n'
            strbuff0 += f'%pointcharges "{self.filename_bq}"\n'
            strbuff0 += '\n'

        # the input block generated by ChemShell is only kept if the user
        # doesn't want to use freetext only (but restart and geometry still
        # must be handled by ChemShell)
        if not self.freetext_only:
            strbuff += strbuff0

        # parallel execution
#        strbuff += '%maxcore 3000\n'
        # TWK 15/5/2024 Re-enable user control of parallelisation of ORCA
        # to allow for the use case of serial ChemShell + parallel ORCA
        if self.nprocs > 0:
            input_nprocs = self.nprocs
        else:
            input_nprocs = parallel.nprocs_per_wg

        if input_nprocs > 1:
            strbuff +=  '# parallel\n'
            strbuff +=  '%pal\n'
            strbuff += f'    nprocs {input_nprocs}\n'
            strbuff +=  'end\n'

        # restart
        if self.restart:
            strbuff +=  '\n# restart\n'
            strbuff +=  '! moread\n'
            strbuff += f'%moinp "{self.vecs_inp}"\n'

        # geometry
        strbuff +=  '\n# geometry\n'
        strbuff +=  '! bohrs\n'
        strbuff += f'*xyz {self.charge} {self.mult}\n'
        strbuff += self.frag.coords2str(order='ac', coords=14, charges=None, atoms='symbols', indent=2*' ', unit='au', notations={'c':self._notation})
        strbuff +=  '*\n'

        return strbuff


    def _getFreetext(self, section):
        '''Return a string containing the freetext content for a specified section'''

        if type(self.freetext) is dict:
            # always attach a new line because we don't expect the user to do so
            strbuff = self.freetext.get(section, '').rstrip()
            strbuff += {'':''}.get(strbuff, '\n')
            return strbuff
        else:
            return ''

    @property
    def runargs(self):
        '''Override the default runargs'''

        # YL 09/11/2021: because ORCA uses OpenMPI and it's very likely ChemShell is driven by Intel MPI
        #                we need --oversubscribe to prevent the error
        #                "There are not enough slots available in the system to satisfy..."
        #                and --mca btl_base_warn_component_unused 0 to suppress the warning
        #                "unable to find any relevant network interfaces"
        # TWK 15/05/2024: in the case of serial ChemShell + parllel ORCA it shouldn't be neceesary
        # to do this so not testing for user set parallelism
        if parallel.nprocs_per_wg > 1:
            return f'{self.input} "--oversubscribe --mca btl_base_warn_component_unused 0"'
        else:
            return self.input


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

        self._runargs = val


    def setOutKeys(self):
        '''Define {task:keyword} for searching in output file'''

        if self.frag.bqs:
            self.frag.nbqs = len(self.frag.bqs.coords)
            ncent = self.frag.natoms + self.frag.nbqs

        # 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 Energy":[0,0,1,3,3,1,1,1.0]},
                          # YL 19/02/2024: it's unsafe to use "Total Energy" yet "FINAL SINGLE POINT ENERGY"
                          #                prints the energy for all ground-state methods (at least in v5.0.4)
                          'energy'   : {'FINAL SINGLE POINT ENERGY':[0,0,1,4,4,1,1,1.0]},
                          # YL 20/02/2024: in excited state calculations "FINAL SINGLE POINT ENERGY" is the energy
                          #                of iroot (state of interest) then we need "TOTAL SCF ENERGY" for the
                          #                ground state (at least for TD-DFT)
                          'energy_0' : {'TOTAL SCF ENERGY':[3,3,1,3,3,1,1,1.0]},
                          'gradients': {"current gradient":[2,(1+3*(self.frag.natoms+self.num_explicit)),1,0,0,1,1,1.0]},
                          'bq_grads' : {"":[1,(self.frag.bqs.nbqs-self.num_explicit),1,0,2,1,1,1.0]},
                        }

        # YL 20/02/2024: multistate calculation
        for i in range(self.nstates):
            # typically:
            # ------------------------------------
            # TD-DFT/TDA EXCITED STATES (SINGLETS)
            # ------------------------------------
            # 
            # the weight of the individual excitations are printed if larger than 1.0e-02
            # 
            # STATE  1:  E=   0.280774 au      7.640 eV    61622.8 cm**-1 <S**2> =   0.000000
            #      4a ->   5a  :     0.993605 (c=  0.99679741) 
            # internally we start from 0 but the NWChem labels start from 1
            self._outkeys.update({'excitation {}'.format(i):{f'STATE {i+1:2d}': [0,0,1,3,11,1,0,1.0]}})
            # YL 20/02/2024: for transition electric dipole moments vs. transition velocity dipole moments, see:
            #                https://mattermodeling.stackexchange.com/questions/8717/different-types-of-transition-dipole-moment-in-orca-output
            self._outkeys.update({'oscillator strength {}'.format(i):{'TRANSITION ELECTRIC DIPOLE MOMENTS': [i+5,i+5,1,3,3,1,1,1.0]}})


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

        from ...io import file
        from numpy import all, array, argmax, argsort, char, sort, vstack

        self.setOutKeys()

#        e = file.getArrayFromFile(self.output, 'out', self._result.energy.shape   , regex=self._outkeys['energy'])
        self._result.energy    = file.getArrayFromFile(self.output, 'out', self._result.energy.shape   , regex=self._outkeys['energy'])
        # If Explicit interface, need to remove BQ-BQ interaction
        if self.explicit:
          self._result.energy += self.energy_correction 

        if self._gradients and os.path.isfile(os.path.abspath(self.output_eandg)):
            temp_grad = zeros(shape=(self.frag.natoms+self.num_explicit,3))
            temp_grad = file.getArrayFromFile(self.output_eandg, 'out', (self.frag.natoms+self.num_explicit, 3), regex=self._outkeys['gradients'])
            for i in range (self.frag.natoms+self.num_explicit):
              if (i >= self.num_explicit):
                self._result.gradients[i-self.num_explicit] = temp_grad[i] 

        if self._gradients and self.frag.nbqs:
            if os.path.isfile(os.path.abspath(self.output_bqs)):
               if self.explicit:
#              If this involves explicit interface, then the bq-bq forces need to be removed
                 self.frag.bqs.gradients = 0.0
                 temp_grads = zeros(shape=(self.frag.bqs.nbqs-self.num_explicit,3))
                 temp_grads = file.getArrayFromFile(self.output_bqs, 'out', (self.frag.bqs.nbqs-self.num_explicit,3), regex=self._outkeys['bq_grads'])
                 exp_counter = 0
                 counter = 0
                 for i in range(self.frag.bqs.nbqs):
                   if i in self.exp_list:
                     self.frag.bqs.gradients[i] = temp_grad[exp_counter]
                     exp_counter += 1
                   else:
                     self.frag.bqs.gradients[i] = temp_grads[counter]
                     counter += 1
                 for i in self.exp_list:
                   coords_i = list(self.frag.bqs.coords[i])
                   for j in range(self.frag.bqs.nbqs):
                     chargeprod = 0.0
                     dist = 0.0
                     coords_j = list(self.frag.bqs.coords[j])
                     chargeprod=-1.0 * self.frag.bqs.charges[i] * self.frag.bqs.charges[j]
                     for k in range(3):
                       dist += (coords_i[k] - coords_j[k])**2
                     dist = math.sqrt(dist) 
                     if (j != i):
                       correction = chargeprod / (dist**3)
                       if (j in self.exp_list) and (j > i):
                         for k in range(3):
                           diff = coords_j[k] - coords_i[k]
                           self.frag.bqs.gradients[i][k] += correction  * diff
                           self.frag.bqs.gradients[j][k] -= correction  * diff
                       if j not in self.exp_list:
                         for k in range(3):
                           diff = coords_j[k] - coords_i[k]
                           self.frag.bqs.gradients[i][k] += correction * diff
                           self.frag.bqs.gradients[j][k] -= correction * diff

               else:
                 self.frag.bqs.gradients = file.getArrayFromFile(self.output_bqs, 'out', self.frag.bqs.gradients.shape, regex=self._outkeys['bq_grads'])

        # YL 20/02/2024: excited states
        if self._excited_state:
            self._result.E0 = file.getArrayFromFile(self.output, 'out', self._result.energy.shape, regex=self._outkeys['energy_0'])[0]
            # YL TODO we could use utils.fileutils.getChoppedMasks() which however needs to be modified
            with open(self.output, 'r') as fp:
                lines = fp.readlines()
                # form a list of starting line of singlets and triplets, respectively, with a trailing None: [index_s,index_t,None]
                indices = [None]
                spins = [ 'SINGLETS', 'TRIPLETS' ]
                for ispin, spin in enumerate(spins):
                    afind = char.rfind(lines, f'EXCITED STATES ({spin})')
                    # remove the current spin if it is not found (in fact
                    # triplets are calculated in addition to singlets in ORCA)
                    if all(afind==-1):
                        indices.insert(ispin, None)
                        spins.pop(ispin)
                    else:
                        indices.insert(ispin, argmax(afind))
                for ispin, spin in enumerate(spins):
                    # loop over states to collect information
                    for i in range(self.nstates):
                        abuff = file.getArrayFromFile(lines[indices[ispin]:indices[ispin+1]],
                                                     'out',
                                                      self._result.energy.shape,
                                                      regex=self._outkeys[f'excitation {i}'],
                                                      number_only=False)
                        # we use column 1 to mark the spin
                        abuff[1] = spin
                        # YL FIXME 20/02/2024: placeholder for symmetry
                        abuff[3] = 'a'
                        try:
                            abuff_sing_trip = vstack((abuff_sing_trip,abuff))
                        except:
                            abuff_sing_trip = abuff
                e = array(abuff_sing_trip)[:,0]
                sorter = argsort(e)
                abuff_sing_trip = array(abuff_sing_trip)[sorter]

            # loop over states to collect information
            for i in range(self.nstates):
                abuff = abuff_sing_trip[i]
                self._result.excitations[i] = self._result.E0+float(abuff[0]),    \
                                              0.0,                                \
                                              unitconvert.eV2nm(float(abuff[2])), \
                                              abuff[3].lower(),                   \
                                              abuff[1].lower()
                # oscillator strengths
                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!
                self._result.energy = self._result.excitations['energies'][self.target_state-1]


        # YL 16/02/2024: if _orca.gbw is the gbw file we read in then _orca.inp can not be the name of the inputfile
        utils.fileutils.rename(utils.fileutils.getBaseName(self.input)+'.gbw', self.vecs_inp)
 

    @property
    def _precisions(self):
        '''Return a dict of simple input precision keywords'''

        return { 3.0e-05:'sloppyscf',
                 1.0e-05:'loosescf',
                 1.0e-06:'normalscf',
                 3.0e-07:'strongscf',
                 1.0e-08:'tightscf',
                 1.0e-09:'verytightscf',
                 1.0e-14:'extremescf',
               }


    @property
    def _precision(self):
        '''ORCA keyword for numerical precision in SCF'''

        from math import inf


        k0 = inf
        for i, (k, v) in enumerate(self._precisions.items()):
            if self.threshold == k:
                strbuff = v
                return strbuff
            elif k0 > self.threshold > k:
                strbuff = v
                return strbuff
            else:
                k0 = k

        return strbuff


