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

_authors__ = ['You Lu <you.lu@stfc.ac.uk>']
# uplifted Nov 2020

import chemsh
from numpy      import full
from ...dl_py2f import py2f, utils
from ...base    import errors
from ...io      import format
from ..qm       import _QM
from .filenames import basis_filenames, potential_filenames

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

    # use mode=1 to prevent OSError: build/gnu/lib/libcp2k.so: undefined symbol: libderiv1_storage_required
    _mode = 1

    # Will try and keep keywords in alphabetical order. 
    _kwargs = { 
                'basis'            :'DZVP',
                'basis_file'       :'',
                'diag_method'      :'STANDARD',
                'functional'       :'PADE',
                'grid_cutoff'      : 280.0,
                'grid_real_cutoff' : 40.0,
                'guess'            :'ATOMIC',
                'input'            :'_cp2k.inp',
                'kpoints_method'   :'MONKHORST-PACK',
                #Rajany: keeping as is now
                'lattice_scale'    : 1.0000,
                #'lattice_scale'    : 3.0000,
                'maxiter_outer'    : 50,
                'mixing_alpha'     : 0.3,
                'mixing_method'    :'DIRECT_P_MIXING',
                'mixing_gmix_p'    : False,
                'ngrids'           : 4,
                'nmos_added'       : 0,
                'outer_scf'        : False,
                'output'           :'_cp2k.out',
                'periodic'         :'XYZ',
                'poisson_solver'   :'PERIODIC',
                'potential_file'   :'',
                'precision_all'    : 1.0E-10,
                'precision_gvg'    : 1.0E-8,
                'precision_rho'    : 1.0E-10,
                'smear'            : False,
                'smear_method'     :'ENERGY_WINDOW',
                'smear_temperature': 300.0,
                'smear_window'     : 0.02,
                'symmetry'         :'auto',
                'wavefunction'     :'_cp2k.wfn',
                'xc_derivative'    :'SPLINE2_SMOOTH',
                'xc_finer_grid'    : False,
                'xc_smooth_rho'    :'SPLINE2',
              }

    _internals = {
#                   '_ecpsection':'bq',
                   '_syscmd'        :'cp2k.sopt',      # default system command to run CP2K
                   '_cell_redundant': False,      # CP2K does NOT allow symmetry-redundant centres in cell
                 }


    # 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 (>=1)
        #       col_0, # starting col
        #       col_n, # ending col
        #       col_i, # col increment (>=1)
        #       skip , # skip the previous (skip-1) match,
        #       coeff, # unit conversion coefficient
        #     ]
        # referencing to the line where the key is matched
        #Rajany: This is the final energy extrapolated for TS --> 0
        energy_string = {r'Total FORCE_EVAL \( QS \) energy' : [0,0,1,8,8,1,1,1.0]}
        #energy_string = {'Total energy:':[0,0,1,2,2,1,1,1.0]}
        self._outkeys = {
                          'energy'    : energy_string,
                          'gradients' :{'ATOMIC FORCES in':[3,3+self.frag.natoms-1,1,3,5,1,1,-1.0]},
                          'density'   :{'Electronic density on regular grids:':[0,0,1,6,6,1,1,1.0]},
                        }


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

        from numpy       import unique, __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
        else:
            from numpy.core  import defchararray
        from ...base     import lattice, math
        from ...data     import unitconvert
        from ...io       import format
        from ...parallel import parallel

        utils.fileutils.archiveFile(self.input)

        # prevent reading result from old results
        frag   = self.frag
        cell   = frag.cell
#        consts = cell.consts

        # QM/Me periodic embedding
        try:
            surface = self._surface
        except:
            surface = 'ab'

        # periodicity for special cases (e.g., QM/Me)
        periodic  = lattice.xyz[surface].strip().upper()
        # non-periodic dimension
        non_periodic_dim = lattice.dimensions[surface]

        booleans = { True:'.TRUE.', False:'.FALSE.' }

        selectcases = {-1: 'SILENT', # almost no output
#                        0: 'LOW'   , # little output
                        0: 'MEDIUM', # little output
                        1: 'MEDIUM', # quite some output
                        2: 'MEDIUM', # quite some output
                        3: 'MEDIUM', # quite some output
                        4: 'HIGH'  , # lots of output
                        5: 'MEDIUM' , # everything is written out, useful for debugging purposes only
                      }

        strbuff =  '&GLOBAL\n'
        strbuff += '    PROJECT     {}\n'.format(self.title)
        if self._gradients:
            strbuff += '    RUN_TYPE    ENERGY_FORCE\n'
        else:
            strbuff += '    RUN_TYPE    ENERGY\n'
        strbuff += '    PRINT_LEVEL {}\n'.format(selectcases.get(chemsh.CHEMSH_DEBUG, 'MEDIUM'))
        strbuff += '    &DBCSR\n'
        # enable memory pools on the CPU
        strbuff += '        USE_MEMPOOLS_CPU  {}\n'.format(booleans[parallel.nprocs_per_wg>1])
        # use MPI allocator to allocate buffers used in MPI communications
        strbuff += '        USE_MPI_ALLOCATOR {}\n'.format(booleans[parallel.nprocs_per_wg>1])
        strbuff += '    &END DBCSR\n'
        strbuff += '&END GLOBAL\n'

        ### FORCE_EVAL ###
        strbuff += '&FORCE_EVAL\n'
        strbuff += '    METHOD        QUICKSTEP\n'
        if self._qmmm_scheme != 'qm/me':
            strbuff += '    STRESS_TENSOR ANALYTICAL\n'

        ###### DFT ######
        strbuff += '    &DFT\n'

        # restart
        if self.restart:
            strbuff += '        WFN_RESTART_FILE_NAME {}\n'.format(self.wavefunction)
        # this keyword can be repeated
        if type(self.basis_file) is str:
            strbuff += '        BASIS_SET_FILE_NAME {}\n'.format(self.basis_file)
        elif type(self.basis_file) is list:
            for basis_file in self.basis_file:
                strbuff += '        BASIS_SET_FILE_NAME {}\n'.format(basis_file)
        # this keyword cannot be repeated
        strbuff += '        POTENTIAL_FILE_NAME {}\n'.format(self.potential_file)

        ######### QS #########
        strbuff += '        &QS\n'
        strbuff += '            EPS_DEFAULT {:e}\n'.format(self.precision_all)
        strbuff += '            EPS_GVG     {:e}\n'.format(self.precision_gvg)
        strbuff += '            EPS_RHO     {:e}\n'.format(self.precision_rho)
        strbuff += '        &END QS\n'

        ######### MGRID #########
        strbuff += '        &MGRID\n'
        strbuff += '            COMMENSURATE .FALSE.\n'
        strbuff += '            NGRIDS        {:d}\n'.format(self.ngrids)
        strbuff += '            CUTOFF        {:f}\n'.format(self.grid_cutoff)
        strbuff += '            REL_CUTOFF    {:f}\n'.format(self.grid_real_cutoff)
        strbuff += '        &END MGRID\n'

        ######### SCF #########
        strbuff += '        &SCF\n'

        if self.level_shift:
            try:
                strbuff += '            LEVEL_SHIFT {:f}\n'.format(self.level_shift)
            except:
                print("\n ChemShell/CP2K ERROR: value of keyword `level_shift` must be a real number!\n")
                exit(112)
        # restart
        if self.restart:
            strbuff += '            SCF_GUESS   RESTART\n'
        # default
        else:
            strbuff += '            SCF_GUESS   {}\n'.format(self.guess.upper().strip())
        # target accuracy for the SCF convergence
        strbuff += '            EPS_SCF     {}\n'.format(self.threshold)
        strbuff += '            MAX_SCF     {}\n'.format(self.maxiter)
        # mixing
        strbuff += '            &MIXING .TRUE.\n'
        strbuff += '                ALPHA   {}\n'.format(self.mixing_alpha)

        # 2nd QM of QM/Me scheme 
        # for 2D systems without adsorbate Broyden mixing causes incomplete convergence
        #Rajany: did not remove the default mixing when no smearing is used
        if (self._qmme == 2) and (not self.smear):
            mixing_method = 'DIRECT_P_MIXING'
        else:
            mixing_method = self.mixing_method
        #mixing_method = self.mixing_method
        mixing_methods = [ 'BROYDEN_MIXING', 'BROYDEN_MIXING_NEW', 'DIRECT_P_MIXING', 'KERKER_MIXING', 'MULTISECANT_MIXING', 'NONE', 'PULAY_MIXING' ]
        if mixing_method.upper().strip() not in mixing_methods:
            print("\n ChemShell/CP2K ERROR: DFT SCF mixing method \"{}\" not valid. Available methods are:\n {}".format(mixing_method.upper().strip(), mixing_methods))
            exit(112)

        strbuff += '                METHOD  {}\n'.format(mixing_method.upper().strip())

        strbuff += '                GMIX_P  {}\n'.format(booleans[self.mixing_gmix_p])
        strbuff += '            &END MIXING\n'

        # 
        smear = self.smear
        # Rajany:adhoc using higher basis set causes slow convergence.
        #if self._qmme ==2:
        #    smear = False
        strbuff += '            &SMEAR {}\n'.format(booleans[smear])
        smear_methods = [ 'ENERGY_WINDOW', 'FERMI_DIRAC', 'LIST' ]
        # Rajany: Print these messages only if smearing is used. Messages are confusing for other tests which don't use it.
        if smear:

            if self.smear_method.upper().strip() not in smear_methods:
                print("\n ChemShell/CP2K ERROR: DFT SCF smearing method \"{}\" not valid. Available methods are:\n {}".format(self.smear_method.upper().strip(), smear_methods), flush=True)
                exit(112)
            if self._qmmm_scheme == 'qm/me' or len(self.kpoints):
                if self.smear_method.strip().upper() != 'FERMI_DIRAC':
                    print('\n'+format.wrap("ChemShell/CP2K NB: Resetting MO smearing method to \"FERMI_DIRAC\"..."))
                    self.smear_method = 'FERMI_DIRAC'

        strbuff += '                METHOD                 {}\n'.format(self.smear_method.upper().strip())
        strbuff += '                ELECTRONIC_TEMPERATURE {}\n'.format(self.smear_temperature)
        strbuff += '                WINDOW_SIZE            {}\n'.format(self.smear_window)
        strbuff += '            &END SMEAR\n'

#Rajany: moved here
        nmos_added = self.nmos_added
        
        if self.smear is True and nmos_added ==0:
            print("\n WARNING! : ADDED MOS cannot be zero when smearing is used! ")
            #exit(100)
        #Rajany : use all available orbitals
        #if self._qmme == 2:
        #    nmos_added = -1

        strbuff += '            ADDED_MOS   {}\n'.format(nmos_added)


        # OT not available with k-points
        if len(self.kpoints):
            strbuff += '            &OT .FALSE.\n'
        else:
            strbuff += '            &OT .TRUE.\n'
        strbuff += '                MINIMIZER      DIIS\n'
        strbuff += '                PRECONDITIONER FULL_ALL\n'
        strbuff += '                ENERGY_GAP         {}\n'.format(0.001)
        strbuff += '            &END OT\n'
        strbuff += '            &OUTER_SCF {}\n'.format(booleans[self.outer_scf])
        strbuff += '                OPTIMIZER DIIS\n'
        strbuff += '                MAX_SCF   {}\n'.format(self.maxiter_outer)
        strbuff += '            &END OUTER_SCF\n'
        strbuff += '            &DIAGONALIZATION {}\n'.format('.TRUE.')
        diag_methods = [ 'DAVIDSON', 'FILTER_MATRIX', 'LANCZOS', 'OT', 'STANDARD' ]
        if self.diag_method.upper().strip() not in diag_methods:
            print("\n ChemShell/CP2K ERROR: DFT SCF diagonalisation method \"{}\" not valid. Available methods are:\n {}".format(self.diag_method.upper().strip(), diag_methods))
            exit(112)
        strbuff += '                ALGORITHM {}\n'.format(self.diag_method.upper().strip())
        strbuff += '            &END DIAGONALIZATION\n'
        strbuff += '            &PRINT\n'
        strbuff += '                &RESTART MEDIUM\n'
        strbuff += '                    BACKUP_COPIES {}\n'.format(0)
        strbuff += '                    FILENAME ={}\n'.format(self.wavefunction)
        strbuff += '                &END RESTART\n'
        strbuff += '                &DIIS_INFO MEDIUM\n'
        strbuff += '                &END DIIS_INFO\n'
        strbuff += '            &END PRINT\n'
        strbuff += '        &END SCF\n'

        ######### XC #########

        # exchange-correlation functional supported by CP2k 7.1
        xc = ['B3LYP','BEEFVDW','BLYP','BP','HCTH120','LDA','NONE','NO_SHORTCUT','OLYP','PADE','PBE','PBE0','TPSS']
        if self.functional.upper().strip() not in xc:
            print("\n >>> ERROR: \"{}\" is not a valid CP2K exchange-correlation functional".format(self.functional.upper()))
            print("\n Available functionals are:\n {}\n".format(xc))
            exit(112)
        strbuff += '        &XC\n'
        strbuff += '            &XC_FUNCTIONAL {}\n'.format(self.functional.upper())
        strbuff += '            &END XC_FUNCTIONAL\n'
        strbuff += '            &XC_GRID\n'
        strbuff += '                USE_FINER_GRID {}\n'.format(booleans[self.xc_finer_grid])
        strbuff += '                XC_DERIV        {}\n'.format(self.xc_derivative) 
        strbuff += '                XC_SMOOTH_RHO   {}\n'.format(self.xc_smooth_rho)
        strbuff += '            &END XC_GRID\n'
        strbuff += '        &END XC\n'
        strbuff += '        &KPOINTS\n'
        strbuff += '            SCHEME    {} {}\n'.format(self.kpoints_method.strip().upper(), ' '.join([ str(i) for i in self.kpoints ]))
        strbuff += '            VERBOSE  .TRUE.\n'
        # Raj: kpoint symmetry leads to reduced grid, hence error-abort when spglib is used
        # Relevant error:k-points need FULL_GRID currently. Reduced grids not yet working correctly,kpoint_methods.F:1501
        strbuff += '            SYMMETRY .FALSE.\n'
        strbuff += '        &END KPOINTS\n'
        strbuff += '        &POISSON\n'
        # this in general should be identical with &SUBSYS&CELL&PERIODIC
        if self._qmmm_scheme == 'qm/me':
            # forced 2D for QM/Me
            # the "periodic" Poisson solver is only valid for 3D
            solver = self.poisson_solver.strip().upper()
            if solver == 'WAVELET':
                if periodic not in [ 'XZ' ]:
                    print('\n'+format.wrap("ChemShell/CP2K NB: Poisson solver \"WAVELET\" only valid for XZ periodicity. Resetting to \"ANALYTIC\"..."))
                    self.poisson_solver = 'ANALYTIC'
            elif solver != 'ANALYTIC':
                if periodic not in [ 'XZ' ]:
                    print('\n'+format.wrap("ChemShell/CP2K NB: Resetting Poisson solver to \"ANALYTIC\" for periodicity {}...\n".format(periodic)))
                    self.poisson_solver = 'ANALYTIC'
                else:
                    print('\n'+format.wrap("ChemShell/CP2K NB: Resetting Poisson solver to \"WAVELET\" for periodicity {}...\n".format(periodic)))
                    self.poisson_solver = 'WAVELET'
            strbuff += '            PERIODIC       {}\n'.format(periodic)
        else:
            strbuff += '            PERIODIC       {}\n'.format(self.periodic.strip().upper())
        strbuff += '            POISSON_SOLVER {}\n'.format(self.poisson_solver.strip().upper())
        strbuff += '        &END POISSON\n'

        ######### PRINT #########
        strbuff += '        &PRINT\n'
        strbuff += '            &MO\n'
        # they are all zeros in CP2K (bug?)
        if chemsh.CHEMSH_DEBUG > 5:
            strbuff += '                EIGENVALUES\n'
            strbuff += '                OCCUPATION_NUMBERS\n'
        strbuff += '            &END MO\n'
        strbuff += '        &END PRINT\n'

        ###### END DFT ######
        strbuff += '    &END DFT\n'
        ###### SUBSYS ######
        strbuff += '    &SUBSYS\n'

        strbuff += '        &TOPOLOGY\n'
        strbuff += '            &DUMP_PDB\n'
        strbuff += '                FILENAME =_cp2k.pdb\n'
        strbuff += '            &END DUMP_PDB\n'
        strbuff += '            &DUMP_PSF\n'
        strbuff += '                FILENAME =_cp2k.psf\n'
        strbuff += '            &END DUMP_PSF\n'
        strbuff += '        &END TOPOLOGY\n'
        strbuff += '        &CELL\n'
        if self._qmmm_scheme == 'qm/me':
            #Rajany: keeping. Examine later
            lattice_scale = self.lattice_scale
            # NB: unit conversion factor in ChemShell and CP2K is both 1.88972666351
            #     however for some reason the ChemShell lattice constant may exceed the lattice size calculated by CP2K
            # scale the non-periodic vector
            factors = [ unitconvert.Bohr2Angstrom if i != non_periodic_dim else unitconvert.Bohr2Angstrom*lattice_scale for i in range(3) ]
        else:
            factors = [ unitconvert.Bohr2Angstrom, unitconvert.Bohr2Angstrom, unitconvert.Bohr2Angstrom ]

        # have to round up to avoid CP2K failure
        a, b, c = math.round_up([self.frag._cell.a*factors[0], self.frag._cell.b*factors[1], self.frag._cell.c*factors[2]], ndecimals=9)

        # even rounding up is not safe!
#        padding = [ 0.0 if i != non_periodic_dim else 1.0E-10 for i in range(3) ]
#        a, b, c = [ v if i != non_periodic_dim else v+padding[i] for i, v in enumerate([a,b,c]) ]

        strbuff += '            ABC              {:>16.10f} {:>16.10f} {:>16.10f}\n'.format(a, b, c)
        strbuff += '            ALPHA_BETA_GAMMA {:>16.10f} {:>16.10f} {:>16.10f}\n'.format(self.frag._cell.alpha, self.frag._cell.beta, self.frag._cell.gamma)

        # forced 2D for QM/Me
        if self._qmmm_scheme == 'qm/me':
            strbuff += '            PERIODIC {}\n'.format(periodic)
        else:
            strbuff += '            PERIODIC {}\n'.format(self.periodic.strip().upper())
        strbuff += '            SYMMETRY {}\n'.format(self.symmetry.strip().upper())
        strbuff += '        &END CELL\n'
        strbuff += '        &COORD\n'
#        # TODO: fractional with SCALED .TRUE.
#                coords2str() method does not yet work for <class Cell>
#        for i in range(cell.natoms):
#            strbuff += '            %s %15.9f %15.9f %15.9f\n'%(cell.names.astype(str)[i], cell.fractional[i][0], cell.fractional[i][1], cell.fractional[i][2])
        strbuff += frag.coords2str(indent=12*' ', atoms='symbols', coords=10, charges=None, unit='angstrom')
        # CP2K has a bug! it doesn't accept upper-case unit name!
#        strbuff += '            UNIT bohr\n'
        strbuff += '        &END COORD\n'
        # single basis set/pseudopotential)
        if type(self.basis) is str and type(self.pseudopotential) is str:
            for species in frag.species:
                strbuff += '        &KIND %s\n'%(species)
                strbuff += '            BASIS_SET %s\n'%self.basis
                strbuff += '            POTENTIAL %s\n'%self.pseudopotential
                strbuff += '        &END KIND\n'
        # multiple basis set/pseudopotential
        elif type(self.basis) is dict or type(self.pseudopotential) is dict:
            for species in frag.species:
                strbuff += '        &KIND {}\n'.format(species)
                # dict
                try:
                    strbuff += '            BASIS_SET {}\n'.format(self.basis.get(species.title(), 'DZVP-ALL-PADE'))
                # str
                except:
                    strbuff += '            BASIS_SET {}\n'.format(self.basis)
                # dict
                try:
                    strbuff += '            POTENTIAL {}\n'.format(self.pseudopotential.get(species.title(), 'GTH-'+self.functional.upper().strip()))
                # str
                except:
                    strbuff += '            POTENTIAL {}\n'.format(self.pseudopotential)
                strbuff += '        &END KIND\n'

        ###### END SUBSYS ######
        strbuff += '    &END SUBSYS\n'
        ###### PRINT ######
        strbuff += '    &PRINT HIGH\n'

        strbuff += '        &FORCES\n'
        strbuff += '        &END\n'
        strbuff += '        &GRID_INFORMATION ON\n'
        strbuff += '        &END\n'
        strbuff += '    &END\n'

        ### END FORCE_EVAL ###
        strbuff += '&END FORCE_EVAL\n'

        return strbuff


#    # overloaded QM.runargs
#    @property
#    def runargs(self):
#        '''Process arguments to be passed to rungamess'''
#
#        # do nothing more than just returning! define `runargs` only in the setter! otherwise its length will grow in every cycle!
#        return '-i ' + self.input + ' -o ' + self.output


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

        from os       import rename
        from ...io    import file, format

        self.setOutKeys()

        # check electronic density (on r-space and g-space grids) for full convergence
        density = file.getArrayFromFile(self.output, 'out', (1,), regex=self._outkeys['density'])[0]
        if density > 0.1:
            print()
            print(format.wrap("ChemShell/CP2K WARNING: the calculation may not have fully converged with a total charge density on r-space grids being {}!\n".format(density)))

        self._result.energy = file.getArrayFromFile(self.output, 'out', 
                                       self._result.energy.shape, regex=self._outkeys['energy'])[0]

        # gradients for QM/Me periodic embedding of 2-step QM calculations: with and without 
        if self._gradients:
            if self._qmmm_scheme == 'qm/me':
            # 1st: with adsorbate
                if self.natoms > len(self._non_adsorbate):
                    self._result.gradients = file.getArrayFromFile(self.output, 'out', 
                                        self._result.gradients.shape, regex=self._outkeys['gradients'])

            # 2nd: without adsorbate
                elif self.natoms == len(self._non_adsorbate):
                    self.grads_qmme = file.getArrayFromFile( self.output,
                                                   'out',
                                                   (len(self._non_adsorbate),3),
                                                    regex=self._outkeys['gradients'])

                else:
                    print("\n ChemShell/CP2K ERROR: wrong number of atoms in QM/Me scheme\n")
                    exit(999)

            else:
            # YL 23/02/2021: it's strange in CP2K if multiple threads are used there are multiple times of energy/gradients are printed
                self._result.gradients = file.getArrayFromFile(self.output, 'out', self._result.gradients.shape, regex=self._outkeys['gradients'])


   
    def parseAdditionalOutput(self):
        '''Overloaded: Parse punch file'''



    def runLib(self, options=None, **kwargs):
        '''Run CP2K as a shared library (called by run.runLib())'''

        from ctypes   import byref

        utils.fileutils.archiveFile(self.output)

        # load CP2K library
        # use mode=1 to prevent OSError: build/gnu/lib/libcp2k.so: undefined symbol: libderiv1_storage_required
        #libcp2k = utils.modutils.getLinkedLib(self.__class__.__name__, mode=self._mode)
        libcp2k = utils.modutils.getLinkedLib(self.__class__.__name__, chemsh.CHEMSH_LIB_DIR)

        fragRef    = py2f(self.frag   , byref=True)
        optionsRef = py2f(self.options, byref=True)

        if libcp2k.interface_cp2k(fragRef, optionsRef):
            raise errors.ChemShModuleError('CP2K')


    @property
    def basis(self):
        '''Basis set'''

        if not self._basis or self._basis == _QM._kwargs['basis']:
            self._basis = 'DZVP-ALL-PADE'
            print("\n ChemShell/CP2K WARNING: a basis set has not been assigned. \"{}\" is used as default.\n".format(self._basis))

        return self._basis


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

        # str
        try:
            self._basis = val.strip().upper()
        # dict
        except:
            self._basis = {}
            for k, v in val.items():
                self._basis.update({k.strip().title():v.strip().upper()})


    @property
    def pseudopotential(self):
        '''Pseudopotential'''

        if not self._pseudopotential:
            self._pseudopotential = 'GTH-'+self.functional.upper().strip()

        return self._pseudopotential


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

        # str
        try:
            self._pseudopotential = val.strip().upper()
        # dict
        except:
            self._pseudopotential = {}
            for k, v in val.items():
                self._pseudopotential.update({k.strip().title():v.strip().upper()})


    @property
    def basis_file(self):
        '''Filename of basis set'''

        from numpy import unique

        if not self._basis_file:
            if type(self.basis) is str:
                self._basis_file = basis_filenames.get(self.basis, 'BASIS_SET')
                print("\n ChemShell/CP2K: automatically assigned filename \"{}\" for".format(self._basis_file))
                print("                 basis set \"{}\"\n".format(self.basis))
            elif type(self.basis) is dict:
                filenames = []
                for k, v in self.basis.items():
                    filenames += [basis_filenames.get(v, 'BASIS_SET')]
                filenames = unique(filenames).tolist()
                print()
                if len(filenames) == 0:
                    print(format.wrap("ChemShell/CP2K ERROR: no filename specified for basis set and ChemShell cannot find one"))
                    exit(112)
                elif len(filenames) == 1:
                    print(format.wrap("ChemShell/CP2K: automatically assigned filename for basis set\n {}".format(filenames[0])))
                elif len(filenames) > 1:
                    print(format.wrap("ChemShell/CP2K: automatically assigned filenames for basis sets\n {}".format(' '.join(filenames))))
                self._basis_file = filenames

        return self._basis_file


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

        self._basis_file = val.strip()


    @property
    def potential_file(self):
        '''Filename of pseudopotential'''

        from numpy import unique

        if not self._potential_file:
            if type(self.pseudopotential) is str:
                self._potential_file = potential_filenames.get(self.pseudopotential, 'POTENTIAL')
                print("\n ChemShell/CP2K: automatically assigned filename for")
                print("                 pseudopotential \"{}\"\n {}".format(self.pseudopotential, self._potential_file))
            elif type(self.pseudopotential) is dict:
                filenames = []
                for k, v in self.pseudopotential.items():
                    filenames += [potential_filenames.get(v, 'POTENTIAL')]
                filenames = unique(filenames) 
                if len(filenames) > 1:
                    print(format.wrap("ChemShell/CP2K ERROR: pseudopotential filename must be unique, however multiple candidates are found:"))
                    print("", filenames)
                    exit(112)
                elif len(filenames) == 1:
                    print(format.wrap("ChemShell/CP2K: automatically assigned filename for pseudopotential\n {}".format(filenames[0])))
                elif len(filenames) == 0:
                    print(format.wrap("ChemShell/CP2K ERROR: no filename specified for pseudopotential and ChemShell cannot find one"))
                    exit(112)

        return self._potential_file


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

        self._potential_file = val.strip()


#    @property
#    def smear(self):
#        '''Smearing of the MO occupation numbers'''
#
#        if len(self.kpoints) and not self._smear:
#            self._smear = ''
#
#        return self._smear
#
#
#    @smear.setter
#    def smear(self, val):
#        '''Setter of smear'''
#
#        self._smear = val.strip().upper()


    @property
    def poisson_solver(self):
        '''Poisson solver'''

#        if self._qmmm_scheme == 'qm/me':
#            print("\n Resetting Poisson to the default method as CP2K does not yet support 2D periodicity...\n")
#            self._poisson_solver = ''

        return self._poisson_solver


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

        self._poisson_solver = val.strip().upper()

        solvers = [ 'ANALYTIC', 'IMPLICIT', 'MT', 'MULTIPOLE', 'PERIODIC', 'WAVELET' ]
        if self._poisson_solver and self._poisson_solver not in solvers:
            print("\n ChemShell/CP2K WARNING: specified Poisson solver \"{}\" may not be supported by CP2K!\n".format(self._poisson_solver))
            print("\n Available solvers are:\n", solvers)


    # the definition and classification are different from those of Cell.pbc so we don't use the latter
    @property
    def symmetry(self):
        '''Lattice symmetry'''

        from numpy import round

        # switched off by the user
        if self._symmetry in [ '', None, False ]:
            self._symmetry = 'NONE'

        # default: always detect the symmetry
        else:

            consts = [ self.frag._cell.a, self.frag._cell.b, self.frag._cell.c, self.frag._cell.alpha, self.frag._cell.beta, self.frag._cell.gamma ]
            a, b, c, alpha, beta, gamma = round(consts, decimals=2)

            if a == b == c and alpha == beta == gamma == 90.0:
                self._symmetry = 'CUBIC'
            elif a == b != c and alpha == beta == 90.0 and gamma == 60.0:
                self._symmetry = 'HEXAGONAL'
            elif a != b != c != a and alpha == gamma == 90.0 and beta != 90.0:
                self._symmetry = 'MONOCLINIC'
            elif a == b != c and alpha == beta == 90.0 and gamma != 90.0:
                self._symmetry = 'MONOCLINIC_GAMMA_AB'
            elif a != b != c and alpha == beta == gamma == 90.0:
                self._symmetry = 'ORTHORHOMBIC'
            elif a == b == c and alpha == beta == gamma != 90.0:
                self._symmetry = 'RHOMBOHEDRAL'
            elif a == b != c and alpha == beta == gamma == 90.0:
                self._symmetry = 'TETRAGONAL_AB'
            elif a == c != b and alpha == beta == gamma == 90.0:
                self._symmetry = 'TETRAGONAL_AC'
            elif a != b == c and alpha == beta == gamma == 90.0:
                self._symmetry = 'TETRAGONAL_BC'
            elif a != b != c != a and alpha != beta != gamma != alpha != 90.0:
                self._symmetry = 'TRICLINIC'
            else:
                self._symmetry = 'NONE'

        return self._symmetry


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

        self._symmetry = val


