import re
import warnings
from copy import deepcopy
from os.path import dirname, isfile
from pathlib import Path

import numpy as np

from ase.atoms import Atoms
from ase.calculators.singlepoint import SinglePointDFTCalculator
from ase.cell import Cell
from ase.units import Bohr

no_positions_error = (
    'no positions can be read from this onetep output '
    'if you wish to use ASE to read onetep outputs '
    'please use uppercase block positions in your calculations'
)

unable_to_read = 'unable to read this onetep output file, ending'

# taken from onetep source code,
# does not seem to be from any known NIST data
units = {'Hartree': 27.2116529, 'Bohr': 1 / 1.889726134583548707935}

# Want to add a functionality? add a global constant below
ONETEP_START = re.compile(
    r'(?i)^\s*\|\s*Linear-Scaling\s*Ab\s*'
    r'Initio\s*Total\s*Energy\s*Program\s*\|\s*$'
)
ONETEP_STOP = re.compile(r'(?i)^\s*-+\s*TIMING\s*INFORMATION\s*-+\s*$')
ONETEP_TOTAL_ENERGY = re.compile(
    r'(?i)^\s*\|\s*\*{3}\s*NGWF\s*' r'optimisation\s*converged\s*\*{3}\s*\|\s*$'
)
ONETEP_FORCE = re.compile(r'(?i)^\s*\*+\s*Forces\s*\*+\s*$')
ONETEP_MULLIKEN = re.compile(r'(?i)^\s*Mulliken\s*Atomic\s*Populations\s*$')
ONETEP_SPIN = re.compile(r'(?i)^\s*Down\s*spin\s*density')
ONETEP_POSITION = re.compile(r'(?i)^\s*Cell\s*Contents\s*$')
ONETEP_FIRST_POSITION = re.compile(
    r'^\s*%BLOCK\s*POSITIONS\s*_?\s*(ABS|FRAC)\s*:?\s*([*#!].*)?$'
)
ONETEP_WRONG_FIRST_POSITION = re.compile(
    r'^\s*%block\s*positions\s*_?\s*(abs|frac)\s*:?\s*([*#!].*)?$'
)
ONETEP_RESUMING_GEOM = re.compile(
    r'(?i)^\s*<{16}\s*Resuming\s*previous'
    r'\s*ONETEP\s*Geometry\s*Optimisation\s*>{16}\s*$'
)

ONETEP_ATOM_COUNT = re.compile(r'(?i)^\s*Totals\s*:\s*(\d+\s*)*$')
ONETEP_IBFGS_ITER = re.compile(r'(?i)^\s*BFGS\s*:\s*starting\s*iteration')
ONETEP_IBFGS_IMPROVE = re.compile(r'(?i)^\s*BFGS\s*:\s*improving\s*iteration')
ONETEP_START_GEOM = re.compile(
    r'(?i)^<+\s*Starting\s*ONETEP\s*Geometry\s*Optimisation\s*>+$'
)
ONETEP_END_GEOM = re.compile(r'(?i)^\s*BFGS\s*:\s*Final\s*Configuration:\s*$')

ONETEP_SPECIES = re.compile(r'(?i)^\s*%BLOCK\s*SPECIES\s*:?\s*([*#!].*)?$')

ONETEP_FIRST_CELL = re.compile(
    r'(?i)^\s*%BLOCK\s*LATTICE\s*_?\s*CART\s*:?\s*([*#!].*)?$'
)
ONETEP_STRESS_CELL = re.compile(
    r'(?i)^\s*stress_calculation:\s*cell\s*geometry\s*$'
)


def get_onetep_keywords(path):
    if isinstance(path, str):
        with open(path) as fd:
            results = read_onetep_in(fd, only_keywords=True)
    else:
        results = read_onetep_in(path, only_keywords=True)

    # If there is an include file, the entire
    # file keyword's will be included in the dict
    # and the include_file keyword will be deleted
    if 'include_file' in results['keywords']:
        warnings.warn('include_file will be deleted from the dict')
        del results['keywords']['include_file']
    return results['keywords']


def read_onetep_in(fd, **kwargs):
    """
    Read a single ONETEP input.

    This function can be used to visually check ONETEP inputs,
    using the ase gui. It can also be used to get access to
    the input parameters attached to the ONETEP calculator
    returned

    The function should work on inputs which contain
    'include_file' command(s), (possibly recursively
    but untested)

    The function should work on input which contains
    exotic element(s) name(s) if the specie block is
    present to map them back to real element(s)

    Parameters
    ----------
    fd : io-object
        File to read.

    Return
    ------
    structure: Atoms
        Atoms object with cell and a Onetep calculator
        attached which contains the keywords dictionary
    """

    fdi_lines = fd.readlines()

    try:
        fd_path = Path(fd.name).resolve()
        fd_parent = fd_path.parent
        include_files = [fd_path]
    except AttributeError:
        # We are in a StringIO or something similar
        fd_path = Path().cwd()
        fd_parent = fd_path
        include_files = [Path().cwd()]

    def clean_lines(lines):
        """
        Remove indesirable line from the input
        """
        new_lines = []
        for line in lines:
            sep = re.split(r'[!#]', line.strip())[0]
            if sep:
                new_lines.append(sep)
        return new_lines

    # Skip comments and empty lines
    fdi_lines = clean_lines(fdi_lines)

    # Are we in a block?
    block_start = 0

    keywords = {}
    atoms = Atoms()
    cell = np.zeros((3, 3))
    fractional = False
    positions = False
    symbols = False

    # Main loop reading the input
    for n, line in enumerate(fdi_lines):
        line_lower = line.lower()
        if re.search(r'^\s*%block', line_lower):
            block_start = n + 1
            if re.search(r'lattice_cart$', line_lower):
                if re.search(r'^\s*ang\s*$', fdi_lines[block_start]):
                    cell = np.loadtxt(fdi_lines[n + 2 : n + 5])
                else:
                    cell = np.loadtxt(fdi_lines[n + 1 : n + 4])
                    cell *= Bohr

        if not block_start:
            if 'devel_code' in line_lower:
                warnings.warn('devel_code is not supported')
                continue
            # Splits line on any valid onetep separator
            sep = re.split(r'[:=\s]+', line)
            keywords[sep[0]] = ' '.join(sep[1:])
            # If include_file is used, we open the included file
            # and insert it in the current fdi_lines...
            # ONETEP does not work with cascade
            # and this SHOULD NOT work with cascade
            if re.search(r'^\s*include_file$', sep[0]):
                name = sep[1].replace("'", '')
                name = name.replace('"', '')
                new_path = fd_parent / name
                for path in include_files:
                    if new_path.samefile(path):
                        raise ValueError('invalid/recursive include_file')
                new_fd = open(new_path)
                new_lines = new_fd.readlines()
                new_lines = clean_lines(new_lines)
                for include_line in new_lines:
                    sep = re.split(r'[:=\s]+', include_line)
                    if re.search(r'^\s*include_file$', sep[0]):
                        raise ValueError('nested include_file')
                fdi_lines[:] = (
                    fdi_lines[: n + 1] + new_lines + fdi_lines[n + 1 :]
                )
                include_files.append(new_path)
                continue

        if re.search(r'^\s*%endblock', line_lower):
            if re.search(r'\s*positions_', line_lower):
                head = re.search(r'(?i)^\s*(\S*)\s*$', fdi_lines[block_start])
                head = head.group(1).lower() if head else ''
                conv = 1 if head == 'ang' else units['Bohr']
                # Skip one line if head is True
                to_read = fdi_lines[block_start + int(bool(head)) : n]
                positions = np.loadtxt(to_read, usecols=(1, 2, 3))
                positions *= conv
                symbols = np.loadtxt(to_read, usecols=(0), dtype='str')
                if re.search(r'.*frac$', line_lower):
                    fractional = True
            elif re.search(r'^\s*%endblock\s*species$', line_lower):
                els = fdi_lines[block_start:n]
                species = {}
                for el in els:
                    sep = el.split()
                    species[sep[0]] = sep[1]
                to_read = [i.strip() for i in fdi_lines[block_start:n]]
                keywords['species'] = to_read
            elif re.search(r'lattice_cart$', line_lower):
                pass
            else:
                to_read = [i.strip() for i in fdi_lines[block_start:n]]
                block_title = line_lower.replace('%endblock', '').strip()
                keywords[block_title] = to_read
            block_start = 0

    # We don't need a fully valid onetep
    # input to read the keywords, just
    # the keywords
    if kwargs.get('only_keywords', False):
        return {'keywords': keywords}
    # Necessary if we have only one atom
    # Check if the cell is valid (3D)
    if not cell.any(axis=1).all():
        raise ValueError('invalid cell specified')

    if positions is False:
        raise ValueError('invalid position specified')

    if symbols is False:
        raise ValueError('no symbols found')

    positions = positions.reshape(-1, 3)
    symbols = symbols.reshape(-1)
    tags = []
    info = {'onetep_species': []}
    for symbol in symbols:
        label = symbol.replace(species[symbol], '')
        if label.isdigit():
            tags.append(int(label))
        else:
            tags.append(0)
        info['onetep_species'].append(symbol)
    atoms = Atoms(
        [species[i] for i in symbols], cell=cell, pbc=True, tags=tags, info=info
    )
    if fractional:
        atoms.set_scaled_positions(positions / units['Bohr'])
    else:
        atoms.set_positions(positions)
    results = {'atoms': atoms, 'keywords': keywords}
    return results


def write_onetep_in(
    fd,
    atoms,
    edft=False,
    xc='PBE',
    ngwf_count=-1,
    ngwf_radius=9.0,
    keywords={},
    pseudopotentials={},
    pseudo_path='.',
    pseudo_suffix=None,
    **kwargs,
):
    """
    Write a single ONETEP input.

    This function will be used by ASE to perform
    various workflows (Opt, NEB...) or can be used
    manually to quickly create ONETEP input file(s).

    The function will first write keywords in
    alphabetic order in lowercase. Secondly, blocks
    will be written in alphabetic order in uppercase.

    Two ways to work with the function:

        - By providing only (simple) keywords present in
          the parameters. ngwf_count and ngwf_radius
          accept multiple types as described in the Parameters
          section.

        - If the keywords parameters is provided as a dictionary
          these keywords will be used to write the input file and
          will take priority.

    If no pseudopotentials are provided in the parameters and
    the function will try to look for suitable pseudopotential
    in the pseudo_path.

    Parameters
    ----------
    fd : file
        File to write.
    atoms: Atoms
        Atoms including Cell object to write.
    edft: Bool
        Activate EDFT.
    xc: str
        DFT xc to use e.g (PBE, RPBE, ...)
    ngwf_count: int|list|dict
        Behaviour depends on the type:
            int: every species will have this amount
            of ngwfs.
            list: list of int, will be attributed
            alphabetically to species:
            dict: keys are species name(s),
            value are their number:
    ngwf_radius: int|list|dict
        Behaviour depends on the type:
            float: every species will have this radius.
            list: list of float, will be attributed
            alphabetically to species:
            [10.0, 9.0]
            dict: keys are species name(s),
            value are their radius:
            {'Na': 9.0, 'Cl': 10.0}
    keywords: dict
        Dictionary with ONETEP keywords to write,
        keywords with lists as values will be
        treated like blocks, with each element
        of list being a different line.
    pseudopotentials: dict
        Behaviour depends on the type:
            keys are species name(s) their
            value are the pseudopotential file to use:
            {'Na': 'Na.usp', 'Cl': 'Cl.usp'}
    pseudo_path: str
        Where to look for pseudopotential, correspond
        to the pseudo_path keyword of ONETEP.
    pseudo_suffix: str
        Suffix for the pseudopotential filename
        to look for, useful if you have multiple sets of
        pseudopotentials in pseudo_path.
    """

    label = kwargs.get('label', 'onetep')
    try:
        directory = kwargs.get('directory', Path(dirname(fd.name)))
    except AttributeError:
        directory = '.'
    autorestart = kwargs.get('autorestart', False)
    elements = np.array(atoms.symbols)
    tags = np.array(atoms.get_tags())
    species_maybe = atoms.info.get('onetep_species', False)
    #  We look if the atom.info contains onetep species information
    # If it does, we use it, as it might contains character
    #  which are not allowed in ase tags, if not we fall back
    # to tags and use them instead.
    if species_maybe:
        if set(species_maybe) != set(elements):
            species = np.array(species_maybe)
        else:
            species = elements
    else:
        formatted_tags = np.array(['' if i == 0 else str(i) for i in tags])
        species = np.char.add(elements, formatted_tags)
    numbers = np.array(atoms.numbers)
    tmp = np.argsort(species)
    # We sort both Z and name the same
    numbers = np.take_along_axis(numbers, tmp, axis=0)
    # u_elements = np.take_along_axis(elements, tmp, axis=0)
    u_species = np.take_along_axis(species, tmp, axis=0)
    elements = np.take_along_axis(elements, tmp, axis=0)
    # We want to keep unique but without sort: small trick with index
    idx = np.unique(u_species, return_index=True)[1]
    elements = elements[idx]
    # Unique species
    u_species = u_species[idx]
    numbers = numbers[idx]
    n_sp = len(u_species)

    if isinstance(ngwf_count, int):
        ngwf_count = dict(zip(u_species, [ngwf_count] * n_sp))
    elif isinstance(ngwf_count, list):
        ngwf_count = dict(zip(u_species, ngwf_count))
    elif isinstance(ngwf_count, dict):
        pass
    else:
        raise TypeError('ngwf_count can only be int|list|dict')

    if isinstance(ngwf_radius, float):
        ngwf_radius = dict(zip(u_species, [ngwf_radius] * n_sp))
    elif isinstance(ngwf_radius, list):
        ngwf_radius = dict(zip(u_species, ngwf_radius))
    elif isinstance(ngwf_radius, dict):
        pass
    else:
        raise TypeError('ngwf_radius can only be float|list|dict')

    pp_files = re.sub('\'|"', '', keywords.get('pseudo_path', pseudo_path))
    pp_files = Path(pp_files).glob('*')
    pp_files = [i for i in pp_files if i.is_file()]

    if pseudo_suffix:
        common_suffix = [pseudo_suffix]
    else:
        common_suffix = ['.usp', '.recpot', '.upf', '.paw', '.psp', '.pspnc']

    if keywords.get('species_pot', False):
        pp_list = keywords['species_pot']
    elif isinstance(pseudopotentials, dict):
        pp_list = []
        for idx, el in enumerate(u_species):
            if el in pseudopotentials:
                pp_list.append(f'{el} {pseudopotentials[el]}')
            else:
                for i in pp_files:
                    reg_el_candidate = re.split(r'[-_.:= ]+', i.stem)[0]
                    if (
                        elements[idx] == reg_el_candidate.title()
                        and i.suffix.lower() in common_suffix
                    ):
                        pp_list.append(f'{el} {i.name}')
    else:
        raise TypeError('pseudopotentials object can only be dict')

    default_species = []
    for idx, el in enumerate(u_species):
        tmp = ''
        tmp += u_species[idx] + ' ' + elements[idx] + ' '
        tmp += str(numbers[idx]) + ' '
        try:
            tmp += str(ngwf_count[el]) + ' '
        except KeyError:
            tmp += str(ngwf_count[elements[idx]]) + ' '
        try:
            tmp += str(ngwf_radius[el])
        except KeyError:
            tmp += str(ngwf_radius[elements[idx]])
        default_species.append(tmp)

    positions_abs = ['ang']
    for s, p in zip(species, atoms.get_positions()):
        line = '{s:>5} {0:>12.6f} {1:>12.6f} {2:>12.6f}'.format(s=s, *p)
        positions_abs.append(line)

    lattice_cart = ['ang']
    for axis in atoms.get_cell():
        line = '{:>16.8f} {:>16.8f} {:>16.8f}'.format(*axis)
        lattice_cart.append(line)

    # Default keywords if not provided by the user,
    # most of them are ONETEP default, except write_forces
    # which is always turned on.
    default_keywords = {
        'xc_functional': xc,
        'edft': edft,
        'cutoff_energy': 20,
        'paw': False,
        'task': 'singlepoint',
        'output_detail': 'normal',
        'species': default_species,
        'pseudo_path': pseudo_path,
        'species_pot': pp_list,
        'positions_abs': positions_abs,
        'lattice_cart': lattice_cart,
        'write_forces': True,
        'forces_output_detail': 'verbose',
    }

    # Main loop, fill the keyword dictionary
    keywords = {key.lower(): value for key, value in keywords.items()}
    for value in default_keywords:
        if not keywords.get(value, None):
            keywords[value] = default_keywords[value]

    # No pseudopotential provided, we look for them in pseudo_path
    # If autorestart is True, we look for restart files,
    # and turn on relevant keywords...
    if autorestart:
        keywords['read_denskern'] = isfile(directory / (label + '.dkn'))
        keywords['read_tightbox_ngwfs'] = isfile(
            directory / (label + '.tightbox_ngwfs')
        )
        keywords['read_hamiltonian'] = isfile(directory / (label + '.ham'))

    # If not EDFT, hamiltonian is irrelevant.
    # print(keywords.get('edft', False))
    # keywords['read_hamiltonian'] = \
    # keywords.get('read_hamiltonian', False) & keywords.get('edft', False)

    keywords = dict(sorted(keywords.items()))

    lines = []
    block_lines = []

    for key, value in keywords.items():
        if isinstance(value, (list, np.ndarray)):
            if not all(isinstance(_, str) for _ in value):
                raise TypeError('list values for blocks must be strings only')
            block_lines.append(('\n%block ' + key).upper())
            block_lines.extend(value)
            block_lines.append(('%endblock ' + key).upper())
        elif isinstance(value, bool):
            lines.append(str(key) + ' : ' + str(value)[0])
        elif isinstance(value, (str, int, float)):
            lines.append(str(key) + ' : ' + str(value))
        else:
            raise TypeError('keyword values must be list|str|bool')
    input_header = (
        '!'
        + '-' * 78
        + '!\n'
        + '!'
        + '-' * 33
        + ' INPUT FILE '
        + '-' * 33
        + '!\n'
        + '!'
        + '-' * 78
        + '!\n\n'
    )

    input_footer = (
        '\n!'
        + '-' * 78
        + '!\n'
        + '!'
        + '-' * 32
        + ' END OF INPUT '
        + '-' * 32
        + '!\n'
        + '!'
        + '-' * 78
        + '!'
    )

    fd.write(input_header)
    fd.writelines(line + '\n' for line in lines)
    fd.writelines(b_line + '\n' for b_line in block_lines)

    if 'devel_code' in kwargs:
        warnings.warn('writing devel code as it is, at the end of the file')
        fd.writelines('\n' + line for line in kwargs['devel_code'])

    fd.write(input_footer)


def read_onetep_out(fd, index=-1, improving=False, **kwargs):
    """
    Read ONETEP output(s).

    !!!
    This function will be used by ASE when performing
    various workflows (opt, NEB...)
    !!!

    Parameters
    ----------
    fd : file
        File to read.
    index: slice
        Which atomic configuration to read
    improving: Bool
        If the output is a geometry optimisation,
        improving = True will keep line search
        configuration from BFGS

    Yields
    ------
    structure: Atoms|list of Atoms
    """
    # Put everything in memory
    fdo_lines = fd.readlines()
    n_lines = len(fdo_lines)

    freg = re.compile(r'-?(?:0|[1-9]\d*)(?:\.\d+)?(?:[eE][+\-]?\d+)?')

    # Used to store index of important elements
    output = {
        ONETEP_START: [],
        ONETEP_STOP: [],
        ONETEP_TOTAL_ENERGY: [],
        ONETEP_FORCE: [],
        ONETEP_SPIN: [],
        ONETEP_MULLIKEN: [],
        ONETEP_POSITION: [],
        ONETEP_FIRST_POSITION: [],
        ONETEP_WRONG_FIRST_POSITION: [],
        ONETEP_ATOM_COUNT: [],
        ONETEP_IBFGS_IMPROVE: [],
        ONETEP_IBFGS_ITER: [],
        ONETEP_START_GEOM: [],
        ONETEP_RESUMING_GEOM: [],
        ONETEP_END_GEOM: [],
        ONETEP_SPECIES: [],
        ONETEP_FIRST_CELL: [],
        ONETEP_STRESS_CELL: [],
    }

    # Index will be treated to get rid of duplicate or improving iterations
    output_corr = deepcopy(output)

    # Core properties that will be used in Yield
    properties = [
        ONETEP_TOTAL_ENERGY,
        ONETEP_FORCE,
        ONETEP_MULLIKEN,
        ONETEP_FIRST_CELL,
    ]

    # Find all matches append them to the dictionary
    breg = '|'.join([i.pattern.replace('(?i)', '') for i in output.keys()])
    prematch = {}

    for idx, line in enumerate(fdo_lines):
        matches = re.search(breg, line)
        if matches:
            prematch[idx] = matches.group(0)

    for key, value in prematch.items():
        for reg in output.keys():
            if re.search(reg, value):
                output[reg].append(key)
                break

    output = {key: np.array(value) for key, value in output.items()}
    # Conveniance notation (pointers: no overhead, no additional memory)
    ibfgs_iter = np.hstack((output[ONETEP_IBFGS_ITER], output[ONETEP_END_GEOM]))
    ibfgs_start = output[ONETEP_START_GEOM]
    ibfgs_improve = output[ONETEP_IBFGS_IMPROVE]
    ibfgs_resume = output[ONETEP_RESUMING_GEOM]

    onetep_start = output[ONETEP_START]
    onetep_stop = output[ONETEP_STOP]

    bfgs_keywords = np.hstack((ibfgs_improve, ibfgs_resume, ibfgs_iter))
    bfgs_keywords = np.sort(bfgs_keywords)

    core_keywords = np.hstack(
        (
            ibfgs_iter,
            ibfgs_start,
            ibfgs_improve,
            ibfgs_resume,
            ibfgs_iter,
            onetep_start,
            onetep_stop,
        )
    )
    core_keywords = np.sort(core_keywords)

    i_first_positions = output[ONETEP_FIRST_POSITION]
    is_frac_positions = [i for i in i_first_positions if 'FRAC' in fdo_lines[i]]

    # In onetep species can have arbritary names,
    # We want to map them to real element names
    # Via the species block
    # species = np.concatenate((output[ONETEP_SPECIES],
    #                          output[ONETEP_SPECIESL])).astype(np.int32)
    species = output[ONETEP_SPECIES]

    icells = np.hstack((output[ONETEP_FIRST_CELL], output[ONETEP_STRESS_CELL]))
    icells = icells.astype(np.int32)
    # Using the fact that 0 == False and > 0 == True
    has_bfgs = (
        len(ibfgs_iter)
        + len(output[ONETEP_START_GEOM])
        + len(output[ONETEP_RESUMING_GEOM])
    )

    # When the input block position is written in lowercase
    # ONETEP does not print the initial position but a hash
    # of it, might be needed
    has_hash = len(output[ONETEP_WRONG_FIRST_POSITION])

    def is_in_bfgs(idx):
        """
        Check if a given index is in a BFGS block
        """
        for past, future in zip(
            output[ONETEP_START],
            np.hstack((output[ONETEP_START][1:], [n_lines])),
        ):
            if past < idx < future:
                if np.any(
                    (past < ibfgs_start) & (ibfgs_start < future)
                ) or np.any((past < ibfgs_resume) & (ibfgs_resume < future)):
                    return True
        return False

    def where_in_bfgs(idx):
        for past, future in zip(
            core_keywords, np.hstack((core_keywords[1:], [n_lines]))
        ):
            if past < idx < future:
                if past in onetep_start:
                    if future in ibfgs_start or future in ibfgs_resume:
                        return 'resume'
                    continue
                # Are we in start or resume or improve
                if past in ibfgs_start:
                    return 'start'
                elif past in ibfgs_resume:
                    return 'resume'
                elif past in ibfgs_improve:
                    return 'improve'

        return False

    ipositions = np.hstack((output[ONETEP_POSITION], i_first_positions)).astype(
        np.int32
    )
    ipositions = np.sort(ipositions)

    n_pos = len(ipositions)

    # Some ONETEP files will not have any positions
    # due to how the software is coded. As a last
    # resort we look for a geom file with the same label.

    if n_pos == 0:
        if has_hash:
            raise RuntimeError(no_positions_error)
        raise RuntimeError(unable_to_read)

    to_del = []

    # Important loop which:
    # - Get rid of improving BFGS iteration if improving == False
    # - Append None to properties to make sure each properties will
    # have the same length and each index correspond to the right
    # atomic configuration (hopefully).
    # Past is the index of the current atomic conf, future is the
    # index of the next one.

    for idx, (past, future) in enumerate(
        zip(ipositions, np.hstack((ipositions[1:], [n_lines])))
    ):
        if has_bfgs:
            which_bfgs = where_in_bfgs(past)

            if which_bfgs == 'resume':
                to_del.append(idx)
                continue

            if not improving:
                if which_bfgs == 'improve':
                    to_del.append(idx)
                    continue

        # We append None if no properties in contained for
        # one specific atomic configurations.
        for prop in properties:
            (tmp,) = np.where((past < output[prop]) & (output[prop] <= future))
            if len(tmp) == 0:
                output_corr[prop].append(None)
            else:
                output_corr[prop].extend(output[prop][tmp[:1]])

    if to_del and len(to_del) != n_pos:
        new_indices = np.setdiff1d(np.arange(n_pos), to_del)
        ipositions = ipositions[new_indices]

    # Bunch of methods to grep properties from output.
    def parse_cell(idx):
        a, b, c = np.loadtxt([fdo_lines[idx + 2]]) * units['Bohr']
        al, be, ga = np.loadtxt([fdo_lines[idx + 4]])
        cell = Cell.fromcellpar([a, b, c, al, be, ga])
        return np.array(cell)

    def parse_charge(idx):
        n = 0
        offset = 4
        while idx + n < len(fdo_lines):
            if not fdo_lines[idx + n].strip():
                tmp_charges = np.loadtxt(
                    fdo_lines[idx + offset : idx + n - 1], usecols=3
                )
                return np.reshape(tmp_charges, -1)
            n += 1
        return None

    #  In ONETEP there is no way to differentiate electronic entropy
    #  and entropy due to solvent, therefore there is no way to
    # extrapolate the energy at 0 K. We return the last energy
    #  instead.

    def parse_energy(idx):
        n = 0
        while idx + n < len(fdo_lines):
            if re.search(r'^\s*\|\s*Total\s*:.*\|\s*$', fdo_lines[idx + n]):
                energy_str = re.search(freg, fdo_lines[idx + n]).group(0)
                return float(energy_str) * units['Hartree']
            n += 1
        return None

    def parse_fermi_level(idx):
        n = 0
        fermi_levels = None
        while idx + n < len(fdo_lines):
            if 'Fermi_level' in fdo_lines[idx + n]:
                tmp = '\n'.join(fdo_lines[idx + n : idx + n + 1])
                fermi_level = re.findall(freg, tmp)
                fermi_levels = [
                    float(i) * units['Hartree'] for i in fermi_level
                ]
            if re.search(
                r'^\s*<{5}\s*CALCULATION\s*SUMMARY\s*>{5}\s*$',
                fdo_lines[idx + n],
            ):
                return fermi_levels
            n += 1
        return None

    def parse_first_cell(idx):
        n = 0
        offset = 1
        while idx + n < len(fdo_lines):
            if re.search(
                r'(?i)^\s*"?\s*ang\s*"?\s*([*#!].*)?$', fdo_lines[idx + n]
            ):
                offset += 1
            if re.search(
                r'(?i)^\s*%ENDBLOCK\s*LATTICE'
                r'\s*_?\s*CART\s*:?\s*([*#!].*)?$',
                fdo_lines[idx + n],
            ):
                cell = np.loadtxt(fdo_lines[idx + offset : idx + n])
                return cell if offset == 2 else cell * units['Bohr']
            n += 1
        return None

    def parse_first_positions(idx):
        n = 0
        offset = 1
        while idx + n < len(fdo_lines):
            if re.search(
                r'(?i)^\s*"?\s*ang\s*"?\s*([*#!].*)?$', fdo_lines[idx + n]
            ):
                offset += 1
            if re.search(r'^\s*%ENDBLOCK\s*POSITIONS_', fdo_lines[idx + n]):
                if 'FRAC' in fdo_lines[idx + n]:
                    conv_factor = 1
                else:
                    conv_factor = units['Bohr']
                tmp = np.loadtxt(
                    fdo_lines[idx + offset : idx + n], dtype='str'
                ).reshape(-1, 4)
                els = np.char.array(tmp[:, 0])
                if offset == 2:
                    pos = tmp[:, 1:].astype(np.float64)
                else:
                    pos = tmp[:, 1:].astype(np.float64) * conv_factor
                try:
                    atoms = Atoms(els, pos)
                # ASE doesn't recognize names used in ONETEP
                # as chemical symbol: dig deeper
                except KeyError:
                    tags, real_elements = find_correct_species(
                        els, idx, first=True
                    )
                    atoms = Atoms(real_elements, pos)
                    atoms.set_tags(tags)
                    atoms.info['onetep_species'] = list(els)
                return atoms
            n += 1
        return None

    def parse_force(idx):
        n = 0
        while idx + n < len(fdo_lines):
            if re.search(r'(?i)^\s*\*\s*TOTAL:.*\*\s*$', fdo_lines[idx + n]):
                tmp = np.loadtxt(
                    fdo_lines[idx + 6 : idx + n - 2],
                    dtype=np.float64,
                    usecols=(3, 4, 5),
                )
                return tmp * units['Hartree'] / units['Bohr']
            n += 1
        return None

    def parse_positions(idx):
        n = 0
        offset = 7
        stop = 0
        while idx + n < len(fdo_lines):
            if re.search(r'^\s*x{60,}\s*$', fdo_lines[idx + n]):
                stop += 1
            if stop == 2:
                tmp = np.loadtxt(
                    fdo_lines[idx + offset : idx + n],
                    dtype='str',
                    usecols=(1, 3, 4, 5),
                )
                els = np.char.array(tmp[:, 0])
                pos = tmp[:, 1:].astype(np.float64) * units['Bohr']
                try:
                    atoms = Atoms(els, pos)
                # ASE doesn't recognize names used in ONETEP
                # as chemical symbol: dig deeper
                except KeyError:
                    tags, real_elements = find_correct_species(els, idx)
                    atoms = Atoms(real_elements, pos)
                    atoms.set_tags(tags)
                    atoms.info['onetep_species'] = list(els)
                return atoms
            n += 1
        return None

    def parse_species(idx):
        n = 1
        element_map = {}
        while idx + n < len(fdo_lines):
            sep = fdo_lines[idx + n].split()
            if re.search(
                r'(?i)^\s*%ENDBLOCK\s*SPECIES\s*:?\s*([*#!].*)?$',
                fdo_lines[idx + n],
            ):
                return element_map
            to_skip = re.search(
                r'(?i)^\s*(ang|bohr)\s*([*#!].*)?$', fdo_lines[idx + n]
            )
            if not to_skip:
                element_map[sep[0]] = sep[1]
            n += 1
        return None

    def parse_spin(idx):
        n = 0
        offset = 4
        while idx + n < len(fdo_lines):
            if not fdo_lines[idx + n].strip():
                # If no spin is present we return None
                try:
                    tmp_spins = np.loadtxt(
                        fdo_lines[idx + offset : idx + n - 1], usecols=4
                    )
                    return np.reshape(tmp_spins, -1)
                except ValueError:
                    return None
            n += 1
        return None

    # This is needed if ASE doesn't recognize the element
    def find_correct_species(els, idx, first=False):
        real_elements = []
        tags = []
        # Find nearest species block in case of
        # multi-output file with different species blocks.
        if first:
            closest_species = np.argmin(abs(idx - species))
        else:
            tmp = idx - species
            tmp[tmp < 0] = 9999999999
            closest_species = np.argmin(tmp)
        elements_map = real_species[closest_species]
        for el in els:
            real_elements.append(elements_map[el])
            tag_maybe = el.replace(elements_map[el], '')
            if tag_maybe.isdigit():
                tags.append(int(tag_maybe))
            else:
                tags.append(False)
        return tags, real_elements

    cells = []
    for idx in icells:
        if idx in output[ONETEP_STRESS_CELL]:
            cell = parse_cell(idx) if idx else None
        else:
            cell = parse_first_cell(idx) if idx else None
        cells.append(cell)

    charges = []
    for idx in output_corr[ONETEP_MULLIKEN]:
        charge = parse_charge(idx) if idx else None
        charges.append(charge)

    energies = []
    for idx in output_corr[ONETEP_TOTAL_ENERGY]:
        energy = parse_energy(idx) if idx else None
        energies.append(energy)

    fermi_levels = []
    for idx in output_corr[ONETEP_TOTAL_ENERGY]:
        fermi_level = parse_fermi_level(idx) if idx else None
        fermi_levels.append(fermi_level)

    magmoms = []
    for idx in output_corr[ONETEP_MULLIKEN]:
        magmom = parse_spin(idx) if idx else None
        magmoms.append(magmom)

    real_species = []
    for idx in species:
        real_specie = parse_species(idx)
        real_species.append(real_specie)

    positions, forces = [], []
    for idx in ipositions:
        if idx in i_first_positions:
            position = parse_first_positions(idx)
        else:
            position = parse_positions(idx)
        if position:
            positions.append(position)
        else:
            n_pos -= 1
            break
    for idx in output_corr[ONETEP_FORCE]:
        force = parse_force(idx) if idx else None
        forces.append(force)

    n_pos = len(positions)

    # Numpy trick to get rid of configuration that are essentially the same
    # in a regular geometry optimisation with internal BFGS, the first
    # configuration is printed three time, we get rid of it
    properties = [energies, forces, charges, magmoms]

    if has_bfgs:
        tmp = [i.positions for i in positions]
        to_del = []
        for i in range(len(tmp[:-1])):
            if is_in_bfgs(ipositions[i]):
                if np.array_equal(tmp[i], tmp[i + 1]):
                    # If the deleted configuration has a property
                    # we want to keep it
                    for prop in properties:
                        if prop[i + 1] is not None and prop[i] is None:
                            prop[i] = prop[i + 1]
                    to_del.append(i + 1)
        c = np.full((len(tmp)), True)
        c[to_del] = False
        energies = [energies[i] for i in range(n_pos) if c[i]]
        forces = [forces[i] for i in range(n_pos) if c[i]]
        charges = [charges[i] for i in range(n_pos) if c[i]]
        magmoms = [magmoms[i] for i in range(n_pos) if c[i]]
        positions = [positions[i] for i in range(n_pos) if c[i]]
        ipositions = [ipositions[i] for i in range(n_pos) if c[i]]

    n_pos = len(positions)

    # We take care of properties that only show up at
    # the beginning of onetep calculation
    whole = np.append(output[ONETEP_START], n_lines)

    if n_pos == 0:
        raise RuntimeError(unable_to_read)

    spin = np.full((n_pos), 1)
    for sp in output[ONETEP_SPIN]:
        tmp = zip(whole, whole[1:])
        for past, future in tmp:
            if past < sp < future:
                p = (past < ipositions) & (ipositions < future)
                spin[p] = 2

    cells_all = np.ones((n_pos, 3, 3))
    for idx, cell in enumerate(output[ONETEP_FIRST_CELL]):
        tmp = zip(whole, whole[1:])
        for past, future in tmp:
            if past < cell < future:
                p = (past < ipositions) & (ipositions < future)
                cells_all[p] = cells[idx]

    # Prepare atom objects with all the properties
    if isinstance(index, int):
        indices = [range(n_pos)[index]]
    else:
        indices = range(n_pos)[index]

    for idx in indices:
        positions[idx].set_cell(cells_all[idx])
        if ipositions[idx] in is_frac_positions:
            positions[idx].set_scaled_positions(positions[idx].get_positions())
        positions[idx].set_initial_charges(charges[idx])
        calc = SinglePointDFTCalculator(
            positions[idx],
            energy=energies[idx] if energies else None,
            free_energy=energies[idx] if energies else None,
            forces=forces[idx] if forces else None,
            charges=charges[idx] if charges else None,
            magmoms=magmoms[idx] if magmoms else None,
        )
        # calc.kpts = [(0, 0, 0) for _ in range(spin[idx])]
        positions[idx].calc = calc
        yield positions[idx]