Skip to main content
`Ctrl`+`K`

#

# Source code for ase.io.gaussian

    # fmt: off

    import logging
    import re
    import warnings
    from collections.abc import Iterable
    from copy import deepcopy

    import numpy as np

    from ase import Atoms
    from ase.calculators.calculator import Calculator, InputError
    from ase.calculators.gaussian import Gaussian
    from ase.calculators.singlepoint import SinglePointCalculator
    from ase.data import atomic_masses_iupac2016, chemical_symbols
    from ase.io import ParseError
    from ase.io.zmatrix import parse_zmatrix
    from ase.units import Bohr, Debye, Hartree

    logger = logging.getLogger(__name__)

    _link0_keys = [
        'mem',
        'chk',
        'oldchk',
        'schk',
        'rwf',
        'oldmatrix',
        'oldrawmatrix',
        'int',
        'd2e',
        'save',
        'nosave',
        'errorsave',
        'cpu',
        'nprocshared',
        'gpucpu',
        'lindaworkers',
        'usessh',
        'ssh',
        'debuglinda',
    ]

    _link0_special = [
        'kjob',
        'subst',
    ]

    # Certain problematic methods do not provide well-defined potential energy
    # surfaces, because these "composite" methods involve geometry optimization
    # and/or vibrational frequency analysis. In addition, the "energy" calculated
    # by these methods are typically ZPVE corrected and/or temperature dependent
    # free energies.
    _problem_methods = [
        'cbs-4m', 'cbs-qb3', 'cbs-apno',
        'g1', 'g2', 'g3', 'g4', 'g2mp2', 'g3mp2', 'g3b3', 'g3mp2b3', 'g4mp4',
        'w1', 'w1u', 'w1bd', 'w1ro',
    ]

    _xc_to_method = dict(
        pbe='pbepbe',
        pbe0='pbe1pbe',
        hse06='hseh1pbe',
        hse03='ohse2pbe',
        lda='svwn',  # gaussian "knows about" LSDA, but maybe not LDA.
        tpss='tpsstpss',
        revtpss='revtpssrevtpss',
    )

    _nuclear_prop_names = ['spin', 'zeff', 'qmom', 'nmagm', 'znuc',
                           'radnuclear', 'iso']

    def _get_molecule_spec(atoms, nuclear_props):
        ''' Generate the molecule specification section to write
        to the Gaussian input file, from the Atoms object and dict
        of nuclear properties'''
        molecule_spec = []
        for i, atom in enumerate(atoms):
            symbol_section = atom.symbol + '('
            # Check whether any nuclear properties of the atom have been set,
            # and if so, add them to the symbol section.
            nuclear_props_set = False
            for keyword, array in nuclear_props.items():
                if array is not None and array[i] is not None:
                    string = keyword + '=' + str(array[i]) + ', '
                    symbol_section += string
                    nuclear_props_set = True

            # Check whether the mass of the atom has been modified,
            # and if so, add it to the symbol section:
            mass_set = False
            symbol = atom.symbol
            expected_mass = atomic_masses_iupac2016[chemical_symbols.index(
                symbol)]
            if expected_mass != atoms[i].mass:
                mass_set = True
                string = 'iso' + '=' + str(atoms[i].mass)
                symbol_section += string

            if nuclear_props_set or mass_set:
                symbol_section = symbol_section.strip(', ')
                symbol_section += ')'
            else:
                symbol_section = symbol_section.strip('(')

            # Then attach the properties appropriately
            # this formatting was chosen for backwards compatibility reasons, but
            # it would probably be better to
            # 1) Ensure proper spacing between entries with explicit spaces
            # 2) Use fewer columns for the element
            # 3) Use 'e' (scientific notation) instead of 'f' for positions
            molecule_spec.append('{:<10s}{:20.10f}{:20.10f}{:20.10f}'.format(
                symbol_section, *atom.position))

        # unit cell vectors, in case of periodic boundary conditions
        for ipbc, tv in zip(atoms.pbc, atoms.cell):
            if ipbc:
                molecule_spec.append('TV {:20.10f}{:20.10f}{:20.10f}'.format(*tv))

        molecule_spec.append('')

        return molecule_spec

    def _format_output_type(output_type):
        ''' Given a letter: output_type, return
        a string formatted for a gaussian input file'''
        # Change output type to P if it has been set to T, because ASE
        # does not support reading info from 'terse' output files.
        if output_type is None or output_type == '' or 't' in output_type.lower():
            output_type = 'P'

        return f'#{output_type}'

    def _check_problem_methods(method):
        ''' Check method string for problem methods and warn appropriately'''
        if method.lower() in _problem_methods:
            warnings.warn(
                'The requested method, {}, is a composite method. Composite '
                'methods do not have well-defined potential energy surfaces, '
                'so the energies, forces, and other properties returned by '
                'ASE may not be meaningful, or they may correspond to a '
                'different geometry than the one provided. '
                'Please use these methods with caution.'.format(method)
            )

    def _pop_link0_params(params):
        '''Takes the params dict and returns a dict with the link0 keywords
        removed, and a list containing the link0 keywords and options
        to be written to the gaussian input file'''
        params = deepcopy(params)
        out = []
        # Remove keywords from params, so they are not set again later in
        # route section
        for key in _link0_keys:
            if key not in params:
                continue
            val = params.pop(key)
            if not val or (isinstance(val, str) and key.lower() == val.lower()):
                out.append(f'%{key}')
            else:
                out.append(f'%{key}={val}')

        # These link0 keywords have a slightly different syntax
        for key in _link0_special:
            if key not in params:
                continue
            val = params.pop(key)
            if not isinstance(val, str) and isinstance(val, Iterable):
                val = ' '.join(val)
            out.append(f'%{key} L{val}')

        return params, out

    def _format_method_basis(output_type, method, basis, fitting_basis):
        output_string = ""
        if basis and method and fitting_basis:
            output_string = '{} {}/{}/{} ! ASE formatted method and basis'.format(
                output_type, method, basis, fitting_basis)
        elif basis and method:
            output_string = '{} {}/{} ! ASE formatted method and basis'.format(
                output_type, method, basis)
        else:
            output_string = f'{output_type}'
            for value in [method, basis]:
                if value is not None:
                    output_string += f' {value}'
        return output_string

    def _format_route_params(params):
        '''Get keywords and values from the params dictionary and return
        as a list of lines to add to the gaussian input file'''
        out = []
        for key, val in params.items():
            # assume bare keyword if val is falsey, i.e. '', None, False, etc.
            # also, for backwards compatibility: assume bare keyword if key and
            # val are the same
            if not val or (isinstance(val, str) and key.lower() == val.lower()):
                out.append(key)
            elif not isinstance(val, str) and isinstance(val, Iterable):
                out.append('{}({})'.format(key, ','.join(val)))
            else:
                out.append(f'{key}({val})')
        return out

    def _format_addsec(addsec):
        '''Format addsec string as a list of lines to be added to the gaussian
         input file.'''
        out = []
        if addsec is not None:
            out.append('')
            if isinstance(addsec, str):
                out.append(addsec)
            elif isinstance(addsec, Iterable):
                out += list(addsec)
        return out

    def _format_basis_set(basis, basisfile, basis_set):
        '''Format either: the basis set filename (basisfile), the basis set file
        contents (from reading basisfile), or the basis_set text as a list of
        strings to be added to the gaussian input file.'''
        out = []
        # if basis='gen', set basisfile. Either give a path to a basisfile, or
        # read in the provided file and paste it verbatim
        if basisfile is not None:
            if basisfile[0] == '@':
                out.append(basisfile)
            else:
                with open(basisfile) as fd:
                    out.append(fd.read())
        elif basis_set is not None:
            out.append(basis_set)
        else:
            if basis is not None and basis.lower() == 'gen':
                raise InputError('Please set basisfile or basis_set')
        return out

    def write_gaussian_in(fd, atoms, properties=['energy'],
                          method=None, basis=None, fitting_basis=None,
                          output_type='P', basisfile=None, basis_set=None,
                          xc=None, charge=None, mult=None, extra=None,
                          ioplist=None, addsec=None, spinlist=None,
                          zefflist=None, qmomlist=None, nmagmlist=None,
                          znuclist=None, radnuclearlist=None,
                          **params):
        '''
        Generates a Gaussian input file

        Parameters
        -----------
        fd: file-like
            where the Gaussian input file will be written
        atoms: Atoms
            Structure to write to the input file
        properties: list
            Properties to calculate
        method: str
            Level of theory to use, e.g. ``hf``, ``ccsd``, ``mp2``, or ``b3lyp``.
            Overrides ``xc`` (see below).
        xc: str
            Level of theory to use. Translates several XC functionals from
            their common name (e.g. ``PBE``) to their internal Gaussian name
            (e.g. ``PBEPBE``).
        basis: str
            The basis set to use. If not provided, no basis set will be requested,
            which usually results in ``STO-3G``. Maybe omitted if basisfile is set
            (see below).
        fitting_basis: str
            The name of the fitting basis set to use.
        output_type: str
            Level of output to record in the Gaussian
            output file - this may be ``N``- normal or ``P`` -
            additional.
        basisfile: str
            The name of the basis file to use. If a value is provided, basis may
            be omitted (it will be automatically set to 'gen')
        basis_set: str
            The basis set definition to use. This is an alternative
            to basisfile, and would be the same as the contents
            of such a file.
        charge: int
            The system charge. If not provided, it will be automatically
            determined from the ``Atoms`` object’s initial_charges.
        mult: int
            The system multiplicity (``spin + 1``). If not provided, it will be
            automatically determined from the ``Atoms`` object’s
            ``initial_magnetic_moments``.
        extra: str
            Extra lines to be included in the route section verbatim.
            It should not be necessary to use this, but it is included for
            backwards compatibility.
        ioplist: list
            A collection of IOPs definitions to be included in the route line.
        addsec: str
            Text to be added after the molecular geometry specification, e.g. for
            defining masses with ``freq=ReadIso``.
        spinlist: list
            A list of nuclear spins to be added into the nuclear
            propeties section of the molecule specification.
        zefflist: list
            A list of effective charges to be added into the nuclear
            propeties section of the molecule specification.
        qmomlist: list
            A list of nuclear quadropole moments to be added into
            the nuclear propeties section of the molecule
            specification.
        nmagmlist: list
            A list of nuclear magnetic moments to be added into
            the nuclear propeties section of the molecule
            specification.
        znuclist: list
            A list of nuclear charges to be added into the nuclear
            propeties section of the molecule specification.
        radnuclearlist: list
            A list of nuclear radii to be added into the nuclear
            propeties section of the molecule specification.
        params: dict
            Contains any extra keywords and values that will be included in either
            the link0 section or route section of the gaussian input file.
            To be included in the link0 section, the keyword must be one of the
            following: ``mem``, ``chk``, ``oldchk``, ``schk``, ``rwf``,
            ``oldmatrix``, ``oldrawmatrix``, ``int``, ``d2e``, ``save``,
            ``nosave``, ``errorsave``, ``cpu``, ``nprocshared``, ``gpucpu``,
            ``lindaworkers``, ``usessh``, ``ssh``, ``debuglinda``.
            Any other keywords will be placed (along with their values) in the
            route section.
        '''

        params = deepcopy(params)

        if properties is None:
            properties = ['energy']

        output_type = _format_output_type(output_type)

        # basis can be omitted if basisfile is provided
        if basis is None:
            if basisfile is not None or basis_set is not None:
                basis = 'gen'

        # determine method from xc if it is provided
        if method is None:
            if xc is not None:
                method = _xc_to_method.get(xc.lower(), xc)

        # If the user requests a problematic method, rather than raising an error
        # or proceeding blindly, give the user a warning that the results parsed
        # by ASE may not be meaningful.
        if method is not None:
            _check_problem_methods(method)

        # determine charge from initial charges if not passed explicitly
        if charge is None:
            charge = atoms.get_initial_charges().sum()

        # determine multiplicity from initial magnetic moments
        # if not passed explicitly
        if mult is None:
            mult = atoms.get_initial_magnetic_moments().sum() + 1

        # set up link0 arguments
        out = []
        params, link0_list = _pop_link0_params(params)
        out.extend(link0_list)

        # begin route line
        # note: unlike in old calculator, each route keyword is put on its own
        # line.
        out.append(_format_method_basis(output_type, method, basis, fitting_basis))

        # If the calculator's parameter dictionary contains an isolist, we ignore
        # this - it is up to the user to attach this info as the atoms' masses
        # if they wish for it to be used:
        params.pop('isolist', None)

        # Any params left will belong in the route section of the file:
        out.extend(_format_route_params(params))

        if ioplist is not None:
            out.append('IOP(' + ', '.join(ioplist) + ')')

        # raw list of explicit keywords for backwards compatibility
        if extra is not None:
            out.append(extra)

        # Add 'force' iff the user requested forces, since Gaussian crashes when
        # 'force' is combined with certain other keywords such as opt and irc.
        if 'forces' in properties and 'force' not in params:
            out.append('force')

        # header, charge, and mult
        out += ['', 'Gaussian input prepared by ASE', '',
                f'{charge:.0f} {mult:.0f}']

        # make dict of nuclear properties:
        nuclear_props = {'spin': spinlist, 'zeff': zefflist, 'qmom': qmomlist,
                         'nmagm': nmagmlist, 'znuc': znuclist,
                         'radnuclear': radnuclearlist}
        nuclear_props = {k: v for k, v in nuclear_props.items() if v is not None}

        # atomic positions and nuclear properties:
        molecule_spec = _get_molecule_spec(atoms, nuclear_props)
        for line in molecule_spec:
            out.append(line)

        out.extend(_format_basis_set(basis, basisfile, basis_set))

        out.extend(_format_addsec(addsec))

        out += ['', '']
        fd.write('\n'.join(out))

    # Regexp for reading an input file:

    _re_link0 = re.compile(r'^\s*%([^\=\)\(!]+)=?([^\=\)\(!]+)?(!.+)?')
    # Link0 lines are in the format:
    # '% keyword = value' or '% keyword'
    # (with or without whitespaces)

    _re_output_type = re.compile(r'^\s*#\s*([NPTnpt]?)\s*')
    # The start of the route section begins with a '#', and then may
    # be followed by the desired level of output in the output file: P, N or T.

    _re_method_basis = re.compile(
        r"\s*([\w-]+)\s*\/([^/=!]+)([\/]([^!]+))?\s*(!.+)?")
    # Matches method, basis and optional fitting basis in the format:
    # method/basis/fitting_basis ! comment
    # They will appear in this format if the Gaussian file has been generated
    # by ASE using a calculator with the basis and method keywords set.

    _re_chgmult = re.compile(r'^(\s*[+-]?\d+(?:,\s*|\s+)[+-]?\d+\s*){1,}$')
    # This is a bit more complex of a regex than we typically want, but it
    # can be difficult to determine whether a line contains the charge and
    # multiplicity, rather than just another route keyword. By making sure
    # that the line contains exactly an even number of *integers*, separated by
    # either a comma (and possibly whitespace) or some amount of whitespace, we
    # can be more confident that we've actually found the charge and multiplicity.
    # Note that in recent versions of Gaussian, the gjf file can contain fragments,
    # where you can give a charge and multiplicity to each fragment. This pattern
    # will allow ASE to read this line in for the charge and multiplicity, however
    # the _get_charge_mult method will only input the first two integers that
    # always gives the overall charge and multiplcity for the full chemical system.
    # The charge and multiplicity of the fragments will be ignored in the
    # _get_charge_mult method.

    _re_nuclear_props = re.compile(r'\(([^\)]+)\)')
    # Matches the nuclear properties, which are contained in parantheses.

    # The following methods are used in GaussianConfiguration's
    # parse_gaussian_input method:

    def _get_link0_param(link0_match):
        '''Gets link0 keyword and option from a re.Match and returns them
        in a dictionary format'''
        value = link0_match.group(2)
        if value is not None:
            value = value.strip()
        else:
            value = ''
        return {link0_match.group(1).lower().strip(): value.lower()}

    def _get_all_link0_params(link0_section):
        ''' Given a string link0_section which contains the link0
        section of a gaussian input file, returns a dictionary of
        keywords and values'''
        parameters = {}
        for line in link0_section:
            link0_match = _re_link0.match(line)
            link0_param = _get_link0_param(link0_match)
            if link0_param is not None:
                parameters.update(link0_param)
        return parameters

    def _convert_to_symbol(string):
        '''Converts an input string into a format
        that can be input to the 'symbol' parameter of an
        ASE Atom object (can be a chemical symbol (str)
        or an atomic number (int)).
        This is achieved by either stripping any
        integers from the string, or converting a string
        containing an atomic number to integer type'''
        symbol = _validate_symbol_string(string)
        if symbol.isnumeric():
            atomic_number = int(symbol)
            symbol = chemical_symbols[atomic_number]
        else:
            match = re.match(r'([A-Za-z]+)', symbol)
            symbol = match.group(1)
        return symbol

    def _validate_symbol_string(string):
        if "-" in string:
            raise ParseError("ERROR: Could not read the Gaussian input file, as"
                             " molecule specifications for molecular mechanics "
                             "calculations are not supported.")
        return string

    def _get_key_value_pairs(line):
        '''Read line of a gaussian input file with keywords and options
        separated according to the rules of the route section.

        Parameters
        ----------
        line (string)
            A line of a gaussian input file.

        Returns
        ---------
        params (dict)
            Contains the keywords and options found in the line.
        '''
        params = {}
        line = line.strip(' #')
        line = line.split('!')[0]  # removes any comments
        # First, get the keywords and options sepatated with
        # parantheses:
        match_iterator = re.finditer(r'\(([^\)]+)\)', line)
        index_ranges = []
        for match in match_iterator:
            index_range = [match.start(0), match.end(0)]
            options = match.group(1)
            # keyword is last word in previous substring:
            keyword_string = line[:match.start(0)]
            keyword_match_iter = [k for k in re.finditer(
                r'[^\,/\s]+', keyword_string) if k.group() != '=']
            keyword = keyword_match_iter[-1].group().strip(' =')
            index_range[0] = keyword_match_iter[-1].start()
            params.update({keyword.lower(): options.lower()})
            index_ranges.append(index_range)

        # remove from the line the keywords and options that we have saved:
        index_ranges.reverse()
        for index_range in index_ranges:
            start = index_range[0]
            stop = index_range[1]
            line = line[0: start:] + line[stop + 1::]

        # Next, get the keywords and options separated with
        # an equals sign, and those without an equals sign
        # must be keywords without options:

        # remove any whitespaces around '=':
        line = re.sub(r'\s*=\s*', '=', line)
        line = [x for x in re.split(r'[\s,\/]', line) if x != '']

        for s in line:
            if '=' in s:
                s = s.split('=')
                keyword = s.pop(0)
                options = s.pop(0)
                params.update({keyword.lower(): options.lower()})
            else:
                if len(s) > 0:
                    params.update({s.lower(): None})

        return params

    def _get_route_params(line):
        '''Reads keywords and values from a line in
        a Gaussian input file's route section,
        and returns them as a dictionary'''
        method_basis_match = _re_method_basis.match(line)
        if method_basis_match:
            params = {}
            ase_gen_comment = '! ASE formatted method and basis'
            if method_basis_match.group(5) == ase_gen_comment:
                params['method'] = method_basis_match.group(1).strip().lower()
                params['basis'] = method_basis_match.group(2).strip().lower()
                if method_basis_match.group(4):
                    params['fitting_basis'] = method_basis_match.group(
                        4).strip().lower()
                return params

        return _get_key_value_pairs(line)

    def _get_all_route_params(route_section):
        ''' Given a string: route_section which contains the route
        section of a gaussian input file, returns a dictionary of
        keywords and values'''
        parameters = {}

        for line in route_section:
            output_type_match = _re_output_type.match(line)
            if not parameters.get('output_type') and output_type_match:
                line = line.strip(output_type_match.group(0))
                parameters.update(
                    {'output_type': output_type_match.group(1).lower()})
            # read route section
            route_params = _get_route_params(line)
            if route_params is not None:
                parameters.update(route_params)
        return parameters

    def _get_charge_mult(chgmult_section):
        '''return a dict with the charge and multiplicity from
        a list chgmult_section that contains the charge and multiplicity
        line, read from a gaussian input file'''
        chgmult_match = _re_chgmult.match(str(chgmult_section))
        try:
            chgmult = chgmult_match.group(0).split()
            return {'charge': int(chgmult[0]), 'mult': int(chgmult[1])}
        except (IndexError, AttributeError):
            raise ParseError("ERROR: Could not read the charge and "
                             "multiplicity from the Gaussian input "
                             "file. There must be an even number of "
                             "integers separated with whitespace or "
                             "a comma, where the first two integers "
                             "must be the overall charge and overall "
                             "multiplicity of the chemical system, "
                             "respectively.")

    def _get_nuclear_props(line):
        ''' Reads any info in parantheses in the line and returns
        a dictionary of the nuclear properties.'''
        nuclear_props_match = _re_nuclear_props.search(line)
        nuclear_props = {}
        if nuclear_props_match:
            nuclear_props = _get_key_value_pairs(nuclear_props_match.group(1))
            updated_nuclear_props = {}
            for key, value in nuclear_props.items():
                if value.isnumeric():
                    value = int(value)
                else:
                    value = float(value)
                if key not in _nuclear_prop_names:
                    if "fragment" in key:
                        warnings.warn("Fragments are not "
                                      "currently supported.")
                    warnings.warn("The following nuclear properties "
                                  "could not be saved: {}".format(
                                      {key: value}))
                else:
                    updated_nuclear_props[key] = value
            nuclear_props = updated_nuclear_props

        for k in _nuclear_prop_names:
            if k not in nuclear_props:
                nuclear_props[k] = None

        return nuclear_props

    def _get_atoms_info(line):
        '''Returns the symbol and position of an atom from a line
        in the molecule specification section'''
        nuclear_props_match = _re_nuclear_props.search(line)
        if nuclear_props_match:
            line = line.replace(nuclear_props_match.group(0), '')
        tokens = line.split()
        symbol = _convert_to_symbol(tokens[0])
        pos = list(tokens[1:])

        return symbol, pos

    def _get_cartesian_atom_coords(symbol, pos):
        '''Returns the coordinates: pos as a list of floats if they
        are cartesian, and not in z-matrix format'''
        if len(pos) < 3 or (pos[0] == '0' and symbol != 'TV'):
            # In this case, we have a z-matrix definition, so
            # no cartesian coords.
            return None
        elif len(pos) > 3:
            raise ParseError("ERROR: Gaussian input file could "
                             "not be read as freeze codes are not"
                             " supported. If using cartesian "
                             "coordinates, these must be "
                             "given as 3 numbers separated "
                             "by whitespace.")
        else:
            try:
                return list(map(float, pos))
            except ValueError:
                raise ParseError(
                    "ERROR: Molecule specification in"
                    "Gaussian input file could not be read")

    def _get_zmatrix_line(line):
        ''' Converts line into the format needed for it to
        be added to the z-matrix contents '''
        line_list = line.split()
        if len(line_list) == 8 and line_list[7] == '1':
            raise ParseError(
                "ERROR: Could not read the Gaussian input file"
                ", as the alternative Z-matrix format using "
                "two bond angles instead of a bond angle and "
                "a dihedral angle is not supported.")
        return line.strip() + '\n'

    def _read_zmatrix(zmatrix_contents, zmatrix_vars=None):
        ''' Reads a z-matrix (zmatrix_contents) using its list of variables
        (zmatrix_vars), and returns atom positions and symbols '''
        try:
            atoms = parse_zmatrix(zmatrix_contents, defs=zmatrix_vars)
        except (ValueError, RuntimeError) as e:
            raise ParseError("Failed to read Z-matrix from "
                             "Gaussian input file: ", e)
        except KeyError as e:
            raise ParseError("Failed to read Z-matrix from "
                             "Gaussian input file, as symbol: {}"
                             "could not be recognised. Please make "
                             "sure you use element symbols, not "
                             "atomic numbers in the element labels.".format(e))
        positions = atoms.positions
        symbols = atoms.get_chemical_symbols()
        return positions, symbols

    def _get_nuclear_props_for_all_atoms(nuclear_props):
        ''' Returns the nuclear properties for all atoms as a dictionary,
        in the format needed for it to be added to the parameters dictionary.'''
        params = {k + 'list': [] for k in _nuclear_prop_names}
        for dictionary in nuclear_props:
            for key, value in dictionary.items():
                params[key + 'list'].append(value)

        for key, array in params.items():
            values_set = False
            for value in array:
                if value is not None:
                    values_set = True
            if not values_set:
                params[key] = None
        return params

    def _get_atoms_from_molspec(molspec_section):
        ''' Takes a string: molspec_section which contains the molecule
        specification section of a gaussian input file, and returns an atoms
        object that represents this.'''
        # These will contain info that will be attached to the Atoms object:
        symbols = []
        positions = []
        pbc = np.zeros(3, dtype=bool)
        cell = np.zeros((3, 3))
        npbc = 0

        # Will contain a dictionary of nuclear properties for each atom,
        # that will later be saved to the parameters dict:
        nuclear_props = []

        # Info relating to the z-matrix definition (if set)
        zmatrix_type = False
        zmatrix_contents = ""
        zmatrix_var_section = False
        zmatrix_vars = ""

        for line in molspec_section:
            # Remove any comments and replace '/' and ',' with whitespace,
            # as these are equivalent:
            line = line.split('!')[0].replace('/', ' ').replace(',', ' ')
            if (line.split()):
                if zmatrix_type:
                    # Save any variables set when defining the z-matrix:
                    if zmatrix_var_section:
                        zmatrix_vars += line.strip() + '\n'
                        continue
                    elif 'variables' in line.lower():
                        zmatrix_var_section = True
                        continue
                    elif 'constants' in line.lower():
                        zmatrix_var_section = True
                        warnings.warn("Constants in the optimisation are "
                                      "not currently supported. Instead "
                                      "setting constants as variables.")
                        continue

                symbol, pos = _get_atoms_info(line)
                current_nuclear_props = _get_nuclear_props(line)

                if not zmatrix_type:
                    pos = _get_cartesian_atom_coords(symbol, pos)
                    if pos is None:
                        zmatrix_type = True

                    if symbol.upper() == 'TV' and pos is not None:
                        pbc[npbc] = True
                        cell[npbc] = pos
                        npbc += 1
                    else:
                        nuclear_props.append(current_nuclear_props)
                        if not zmatrix_type:
                            symbols.append(symbol)
                            positions.append(pos)

                if zmatrix_type:
                    zmatrix_contents += _get_zmatrix_line(line)

        # Now that we are past the molecule spec. section, we can read
        # the entire z-matrix (if set):
        if len(positions) == 0:
            if zmatrix_type:
                if zmatrix_vars == '':
                    zmatrix_vars = None
                positions, symbols = _read_zmatrix(
                    zmatrix_contents, zmatrix_vars)

        try:
            atoms = Atoms(symbols, positions, pbc=pbc, cell=cell)
        except (IndexError, ValueError, KeyError) as e:
            raise ParseError("ERROR: Could not read the Gaussian input file, "
                             "due to a problem with the molecule "
                             "specification: {}".format(e))

        nuclear_props = _get_nuclear_props_for_all_atoms(nuclear_props)

        return atoms, nuclear_props

    def _get_readiso_param(parameters):
        ''' Returns a tuple containing the frequency
        keyword and its options, if the frequency keyword is
        present in parameters and ReadIso is one of its options'''
        freq_options = parameters.get('freq', None)
        if freq_options:
            freq_name = 'freq'
        else:
            freq_options = parameters.get('frequency', None)
            freq_name = 'frequency'
        if freq_options is not None:
            if ('readiso' or 'readisotopes') in freq_options:
                return freq_name, freq_options
        return None, None

    def _get_readiso_info(line, parameters):
        '''Reads the temperature, pressure and scale from the first line
        of a ReadIso section of a Gaussian input file. Returns these in
        a dictionary.'''
        readiso_params = {}
        if _get_readiso_param(parameters)[0] is not None:
            # when count_iso is 0 we are in the line where
            # temperature, pressure, [scale] is saved
            line = line.replace(
                '[', '').replace(']', '')
            tokens = line.strip().split()
            try:
                readiso_params['temperature'] = tokens[0]
                readiso_params['pressure'] = tokens[1]
                readiso_params['scale'] = tokens[2]
            except IndexError:
                pass
        if readiso_params != {}:

            return readiso_params

    def _delete_readiso_param(parameters):
        '''Removes the readiso parameter from the parameters dict'''
        parameters = deepcopy(parameters)
        freq_name, freq_options = _get_readiso_param(parameters)
        if freq_name is not None:
            if 'readisotopes' in freq_options:
                iso_name = 'readisotopes'
            else:
                iso_name = 'readiso'
            freq_options = [v.group() for v in re.finditer(
                r'[^\,/\s]+', freq_options)]
            freq_options.remove(iso_name)
            new_freq_options = ''
            for v in freq_options:
                new_freq_options += v + ' '
            if new_freq_options == '':
                new_freq_options = None
            else:
                new_freq_options = new_freq_options.strip()
            parameters[freq_name] = new_freq_options
        return parameters

    def _update_readiso_params(parameters, symbols):
        ''' Deletes the ReadIso option from the route section as we
        write out the masses in the nuclear properties section
        instead of the ReadIso section.
        Ensures the masses array is the same length as the
        symbols array. This is necessary due to the way the
        ReadIso section is defined:
        The mass of each atom is listed on a separate line, in
        order of appearance in the molecule spec. A blank line
        indicates not to modify the mass for that atom.
        But you do not have to leave blank lines equal to the
        remaining atoms after you finsihed setting masses.
        E.g. if you had 10 masses and only want to set the mass
        for the first atom, you don't have to leave 9 blank lines
        after it.
        '''
        parameters = _delete_readiso_param(parameters)
        if parameters.get('isolist') is not None:
            if len(parameters['isolist']) < len(symbols):
                for _ in range(len(symbols) - len(parameters['isolist'])):
                    parameters['isolist'].append(None)
            if all(m is None for m in parameters['isolist']):
                parameters['isolist'] = None

        return parameters

    def _validate_params(parameters):
        '''Checks whether all of the required parameters exist in the
        parameters dict and whether it contains any unsupported settings
        '''
        # Check for unsupported settings
        unsupported_settings = {
            "z-matrix", "modredun", "modredundant", "addredundant", "addredun",
            "readopt", "rdopt"}

        for s in unsupported_settings:
            for v in parameters.values():
                if v is not None and s in str(v):
                    raise ParseError(
                        "ERROR: Could not read the Gaussian input file"
                        ", as the option: {} is currently unsupported."
                        .format(s))

        for k in list(parameters.keys()):
            if "popt" in k:
                parameters["opt"] = parameters.pop(k)
                warnings.warn("The option {} is currently unsupported. "
                              "This has been replaced with {}."
                              .format("POpt", "opt"))
                return

    def _get_extra_section_params(extra_section, parameters, atoms):
        ''' Takes a list of strings: extra_section, which contains
        the 'extra' lines in a gaussian input file. Also takes the parameters
        that have been read so far, and the atoms that have been read from the
        file.
        Returns an updated version of the parameters dict, containing details from
        this extra section. This may include the basis set definition or filename,
        and/or the readiso section.'''

        new_parameters = deepcopy(parameters)
        # Will store the basis set definition (if set):
        basis_set = ""
        # This will indicate whether we have a readiso section:
        readiso = _get_readiso_param(new_parameters)[0]
        # This will indicate which line of the readiso section we are reading:
        count_iso = 0
        readiso_masses = []

        for line in extra_section:
            if line.split():
                # check that the line isn't just a comment
                if line.split()[0] == '!':
                    continue
                line = line.strip().split('!')[0]

            if len(line) > 0 and line[0] == '@':
                # If the name of a basis file is specified, this line
                # begins with a '@'
                new_parameters['basisfile'] = line
            elif readiso and count_iso < len(atoms.symbols) + 1:
                # The ReadIso section has 1 more line than the number of
                # symbols
                if count_iso == 0 and line != '\n':
                    # The first line in the readiso section contains the
                    # temperature, pressure, scale. Here we save this:
                    readiso_info = _get_readiso_info(line, new_parameters)
                    if readiso_info is not None:
                        new_parameters.update(readiso_info)
                    # If the atom masses were set in the nuclear properties
                    # section, they will be overwritten by the ReadIso
                    # section
                    readiso_masses = []
                    count_iso += 1
                elif count_iso > 0:
                    # The remaining lines in the ReadIso section are
                    # the masses of the atoms
                    try:
                        readiso_masses.append(float(line))
                    except ValueError:
                        readiso_masses.append(None)
                    count_iso += 1
            else:
                # If the rest of the file is not the ReadIso section,
                # then it must be the definition of the basis set.
                if new_parameters.get('basis', '') == 'gen' \
                        or 'gen' in new_parameters:
                    if line.strip() != '':
                        basis_set += line + '\n'

        if readiso:
            new_parameters['isolist'] = readiso_masses
            new_parameters = _update_readiso_params(new_parameters, atoms.symbols)

        # Saves the basis set definition to the parameters array if
        # it has been set:
        if basis_set != '':
            new_parameters['basis_set'] = basis_set
            new_parameters['basis'] = 'gen'
            new_parameters.pop('gen', None)

        return new_parameters

    def _get_gaussian_in_sections(fd):
        ''' Reads a gaussian input file and returns
        a dictionary of the sections of the file - each dictionary
        value is a list of strings - each string is a line in that
        section. '''
        # These variables indicate which section of the
        # input file is currently being read:
        route_section = False
        atoms_section = False
        atoms_saved = False
        # Here we will store the sections of the file in a dictionary,
        # as lists of strings
        gaussian_sections = {'link0': [], 'route': [],
                             'charge_mult': [], 'mol_spec': [], 'extra': []}

        for line in fd:
            line = line.strip(' ')
            link0_match = _re_link0.match(line)
            output_type_match = _re_output_type.match(line)
            chgmult_match = _re_chgmult.match(line)
            if link0_match:
                gaussian_sections['link0'].append(line)
            # The first blank line appears at the end of the route section
            # and a blank line appears at the end of the atoms section:
            elif line == '\n' and (route_section or atoms_section):
                route_section = False
                atoms_section = False
            elif output_type_match or route_section:
                route_section = True
                gaussian_sections['route'].append(line)
            elif chgmult_match:
                gaussian_sections['charge_mult'] = line
                # After the charge and multiplicty have been set, the
                # molecule specification section of the input file begins:
                atoms_section = True
            elif atoms_section:
                gaussian_sections['mol_spec'].append(line)
                atoms_saved = True
            elif atoms_saved:
                # Next we read the other sections below the molecule spec.
                # This may include the ReadIso section and the basis set
                # definition or filename
                gaussian_sections['extra'].append(line)

        return gaussian_sections

    class GaussianConfiguration:

        def __init__(self, atoms, parameters):
            self.atoms = atoms.copy()
            self.parameters = deepcopy(parameters)

        def get_atoms(self):
            return self.atoms

        def get_parameters(self):
            return self.parameters

        def get_calculator(self):
            # Explicitly set parameters that must for security reasons not be
            # taken from file:
            calc = Gaussian(atoms=self.atoms, command=None, restart=None,
                            ignore_bad_restart_file=Calculator._deprecated,
                            label='Gaussian', directory='.', **self.parameters)
            return calc

        @staticmethod
        def parse_gaussian_input(fd):
            '''Reads a gaussian input file into an atoms object and
            parameters dictionary.

            Parameters
            ----------
            fd: file-like
                Contains the contents of a  gaussian input file

            Returns
            ---------
            GaussianConfiguration
                Contains an atoms object created using the structural
                information from the input file.
                Contains a parameters dictionary, which stores any
                keywords and options found in the link-0 and route
                sections of the input file.
            '''
            # The parameters dict to attach to the calculator
            parameters = {}

            file_sections = _get_gaussian_in_sections(fd)

            # Update the parameters dictionary with the keywords and values
            # from each section of the input file:
            parameters.update(_get_all_link0_params(file_sections['link0']))
            parameters.update(_get_all_route_params(file_sections['route']))
            parameters.update(_get_charge_mult(file_sections['charge_mult']))
            atoms, nuclear_props = _get_atoms_from_molspec(
                file_sections['mol_spec'])
            parameters.update(nuclear_props)
            parameters.update(_get_extra_section_params(
                file_sections['extra'], parameters, atoms))

            _validate_params(parameters)

            return GaussianConfiguration(atoms, parameters)

    def read_gaussian_in(fd, attach_calculator=False):
        '''
        Reads a gaussian input file and returns an Atoms object.

        Parameters
        ----------
        fd: file-like
            Contains the contents of a gaussian input file

        attach_calculator: bool
            When set to ``True``, a Gaussian calculator will be
            attached to the Atoms object which is returned.
            This will mean that additional information is read
            from the input file (see below for details).

        Returns
        ----------
        Atoms
            An Atoms object containing the following information that has been
            read from the input file: symbols, positions, cell.

        Notes
        ----------
        Gaussian input files can be read where the atoms' locations are set with
        cartesian coordinates or as a z-matrix. Variables may be used in the
        z-matrix definition, but currently setting constants for constraining
        a geometry optimisation is not supported. Note that the `alternative`
        z-matrix definition where two bond angles are set instead of a bond angle
        and a dihedral angle is not currently supported.

        If the parameter ``attach_calculator`` is set to ``True``, then the Atoms
        object is returned with a Gaussian calculator attached.

        This Gaussian calculator will contain a parameters dictionary which will
        contain the Link 0 commands and the options and keywords set in the route
        section of the Gaussian input file, as well as:

        • The charge and multiplicity

        • The selected level of output

        • The method, basis set and (optional) fitting basis set.

        • Basis file name/definition if set

        • Nuclear properties

        • Masses set in the nuclear properties section or in the ReadIsotopes
          section (if ``freq=ReadIso`` is set). These will be saved in an array
          with keyword ``isolist``, in the parameters dictionary. For these to
          actually be used in calculations and/or written out to input files,
          the Atoms object's masses must be manually updated to these values
          (this is not done automatically)

        If the Gaussian input file has been generated by ASE's
        ``write_gaussian_in`` method, then the basis set, method and fitting
        basis will be saved under the ``basis``, ``method`` and ``fitting_basis``
        keywords in the parameters dictionary. Otherwise they are saved as
        keywords in the parameters dict.

        Currently this does not support reading of any other sections which would
        be found below the molecule specification that have not been mentioned
        here (such as the ModRedundant section).
        '''
        gaussian_input = GaussianConfiguration.parse_gaussian_input(fd)
        atoms = gaussian_input.get_atoms()

        if attach_calculator:
            atoms.calc = gaussian_input.get_calculator()

        return atoms

    # In the interest of using the same RE for both atomic positions and forces,
    # we make one of the columns optional. That's because atomic positions have
    # 6 columns, while forces only has 5 columns. Otherwise they are very similar.
    _re_atom = re.compile(
        r'^\s*\S+\s+(\S+)\s+(?:\S+\s+)?(\S+)\s+(\S+)\s+(\S+)\s*$'
    )
    _re_forceblock = re.compile(r'^\s*Center\s+Atomic\s+Forces\s+\S+\s*$')
    _re_l716 = re.compile(r'^\s*\(Enter .+l716.exe\)$')

    def _compare_merge_configs(configs, new):
        """Append new to configs if it contains a new geometry or new data.

        Gaussian sometimes repeats a geometry, for example at the end of an
        optimization, or when a user requests vibrational frequency
        analysis in the same calculation as a geometry optimization.

        In those cases, rather than repeating the structure in the list of
        returned structures, try to merge results if doing so doesn't change
        any previously calculated values. If that's not possible, then create
        a new "image" with the new results.
        """
        if not configs:
            configs.append(new)
            return

        old = configs[-1]

        if old != new:
            configs.append(new)
            return

        oldres = old.calc.results
        newres = new.calc.results
        common_keys = set(oldres).intersection(newres)

        for key in common_keys:
            if np.any(oldres[key] != newres[key]):
                configs.append(new)
                return
        oldres.update(newres)

    def _read_charges(fd):
        fd.readline()
        qs = []
        ms = []
        for line in fd:
            if not line.strip()[0].isdigit():
                break
            qs.append(float(line.split()[2]))
            if len(line.split()) > 3:
                ms.append(float(line.split()[3]))
        return {'charges': qs, 'magmoms': ms} if ms else {'charges': qs}

    def read_gaussian_out(fd, index=-1):
        """Reads a gaussian output file and returns an Atoms object."""
        configs = []
        atoms = None
        energy = None
        results = {}
        orientation = None  # Orientation of the coordinates stored in atoms
        for line in fd:
            line = line.strip()
            if line.startswith(r'1\1\GINC'):
                # We've reached the "archive" block at the bottom, stop parsing
                break

            if (line == 'Input orientation:'
                    or line == 'Z-Matrix orientation:'
                    or line == "Standard orientation:"):
                if atoms is not None:
                    # Add configuration to the currently-parsed list
                    #  only after an energy or force has been parsed
                    #  (we assume these are in the output file)
                    if energy is None:
                        continue
                    # "atoms" should store the first geometry encountered
                    #  in the input file which is often the input orientation,
                    #  which is the orientation for forces.
                    #  If there are forces and the orientation of atoms is not
                    #  the input coordinate system, warn the user
                    if orientation != "Input" and 'forces' in results:
                        logger.warning('Configuration geometry is not in the input'
                                       f'orientation. It is {orientation}')
                    calc = SinglePointCalculator(atoms, energy=energy, **results)
                    atoms.calc = calc
                    _compare_merge_configs(configs, atoms)
                atoms = None
                orientation = line.split()[0]  # Store the orientation
                energy = None
                results = {}

                numbers = []
                positions = []
                pbc = np.zeros(3, dtype=bool)
                cell = np.zeros((3, 3))
                npbc = 0
                # skip 4 irrelevant lines
                for _ in range(4):
                    fd.readline()
                while True:
                    match = _re_atom.match(fd.readline())
                    if match is None:
                        break
                    number = int(match.group(1))
                    pos = list(map(float, match.group(2, 3, 4)))
                    if number == -2:
                        pbc[npbc] = True
                        cell[npbc] = pos
                        npbc += 1
                    else:
                        numbers.append(max(number, 0))
                        positions.append(pos)
                atoms = Atoms(numbers, positions, pbc=pbc, cell=cell)
            elif (line.startswith('Energy=')
                    or line.startswith('SCF Done:')):
                # Some semi-empirical methods (Huckel, MINDO3, etc.),
                # or SCF methods (HF, DFT, etc.)
                energy = float(line.split('=')[1].split()[0].replace('D', 'e'))
                energy *= Hartree
            elif (line.startswith('E2 =') or line.startswith('E3 =')
                    or line.startswith('E4(') or line.startswith('DEMP5 =')
                    or line.startswith('E2(')):
                # MP{2,3,4,5} energy
                # also some double hybrid calculations, like B2PLYP
                energy = float(line.split('=')[-1].strip().replace('D', 'e'))
                energy *= Hartree
            elif line.startswith('Wavefunction amplitudes converged. E(Corr)'):
                # "correlated method" energy, e.g. CCSD
                energy = float(line.split('=')[-1].strip().replace('D', 'e'))
                energy *= Hartree
            elif line.startswith('CCSD(T)='):
                # CCSD(T) energy
                energy = float(line.split('=')[-1].strip().replace('D', 'e'))
                energy *= Hartree
            elif (
                line.startswith('Mulliken charges')
                or line.startswith('Lowdin Atomic Charges')
                or line.startswith('Hirshfeld charges, spin densities,')
            ):
                # Löwdin is printed after Mulliken and overwrites `charges`.
                # Hirshfeld is printed after Mulliken and overwrites `charges`.
                results.update(_read_charges(fd))
            elif line.startswith('Dipole moment ') and energy is not None:
                # dipole moment in `l601.exe`, printed unless `Pop=None`
                # Skipped if energy is not printed in the same section.
                # This happens in the last geometry record when `opt` or `irc` is
                # specified. In this case, the record is compared with the previous
                # one in `_compare_merge_configs`, and there the dipole moment
                # from `l601` conflicts with the previous record from `l716`.
                line = fd.readline().strip()
                dipole = np.array([float(_) for _ in line.split()[1:6:2]]) * Debye
                results['dipole'] = dipole
            elif re.match(r"^Recovered energy=.*dipole=.*$", line):
                # energy and dipole moment obtained from an external script
                energy = float(line.split()[2].replace('D', 'e'))
                energy *= Hartree
                dipole = np.array([float(_) for _ in line.split()[4:7]]) * Debye
                results['dipole'] = dipole
            elif _re_l716.match(line):
                # Sometimes Gaussian will print "Rotating derivatives to
                # standard orientation" after the matched line (which looks like
                # "(Enter /opt/gaussian/g16/l716.exe)", though the exact path
                # depends on where Gaussian is installed). We *skip* the dipole
                # in this case, because it might be rotated relative to the input
                # orientation (and also it is numerically different even if the
                # standard orientation is the same as the input orientation).
                line = fd.readline().strip()
                if not line.startswith('Dipole'):
                    continue
                dip = line.split('=')[1].replace('D', 'e')
                tokens = dip.split()
                dipole = []
                # dipole elements can run together, depending on what method was
                # used to calculate them. First see if there is a space between
                # values.
                if len(tokens) == 3:
                    dipole = list(map(float, tokens))
                elif len(dip) % 3 == 0:
                    # next, check if the number of tokens is divisible by 3
                    nchars = len(dip) // 3
                    for i in range(3):
                        dipole.append(float(dip[nchars * i:nchars * (i + 1)]))
                else:
                    # otherwise, just give up on trying to parse it.
                    dipole = None
                    continue
                # this dipole moment is printed in atomic units, e-Bohr
                # ASE uses e-Angstrom for dipole moments.
                results['dipole'] = np.array(dipole) * Bohr
            elif _re_forceblock.match(line):
                # skip 2 irrelevant lines
                fd.readline()
                fd.readline()
                forces = []
                while True:
                    match = _re_atom.match(fd.readline())
                    if match is None:
                        break
                    forces.append(list(map(float, match.group(2, 3, 4))))
                results['forces'] = np.array(forces) * Hartree / Bohr
        if atoms is not None:
            atoms.calc = SinglePointCalculator(atoms, energy=energy, **results)
            _compare_merge_configs(configs, atoms)
        return configs[index]