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, AssertionError) 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_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]