# -*- coding: utf-8 -*-
#
# Copyright (c) 2018, the cclib development team
#
# This file is part of cclib (http://cclib.github.io) and is distributed under
# the terms of the BSD 3-Clause License.

"""Parser for GAMESS(US) output files"""


from __future__ import print_function
import re

import numpy


from cclib.parser import logfileparser
from cclib.parser import utils


class GAMESS(logfileparser.Logfile):
    """A GAMESS/Firefly log file."""

    # Used to index self.scftargets[].
    SCFRMS, SCFMAX, SCFENERGY = list(range(3))

    # Used to extact Dunning basis set names.
    dunningbas = {'CCD': 'cc-pVDZ', \
            'CCT': 'cc-pVTZ', \
            'CCQ': 'cc-pVQZ', \
            'CC5': 'cc-pV5Z', \
            'CC6': 'cc-pV6Z', \
            'ACCD': 'aug-cc-pVDZ', \
            'ACCT': 'aug-cc-pVTZ', \
            'ACCQ': 'aug-cc-pVQZ', \
            'ACC5': 'aug-cc-pV5Z', \
            'ACC6': 'aug-cc-pV6Z', \
            'CCDC': 'cc-pCVDZ', \
            'CCTC': 'cc-pCVTZ', \
            'CCQC': 'cc-pCVQZ', \
            'CC5C': 'cc-pCV5Z', \
            'CC6C': 'cc-pCV6Z', \
            'ACCDC': 'aug-cc-pCVDZ', \
            'ACCTC': 'aug-cc-pCVTZ', \
            'ACCQC': 'aug-cc-pCVQZ', \
            'ACC5C': 'aug-cc-pCV5Z', \
            'ACC6C': 'aug-cc-pCV6Z'}

    def __init__(self, *args, **kwargs):

        # Call the __init__ method of the superclass
        super(GAMESS, self).__init__(logname="GAMESS", *args, **kwargs)

    def __str__(self):
        """Return a string representation of the object."""
        return "GAMESS log file %s" % (self.filename)

    def __repr__(self):
        """Return a representation of the object."""
        return 'GAMESS("%s")' % (self.filename)

    def normalisesym(self, label):
        """Normalise the symmetries used by GAMESS.

        To normalise, two rules need to be applied:
        (1) Occurences of U/G in the 2/3 position of the label
            must be lower-cased
        (2) Two single quotation marks must be replaced by a double
        """

        if label[1:] == "''":
            end = '"'
        else:
            end = label[1:].replace("U", "u").replace("G", "g")
        return label[0] + end

    def before_parsing(self):

        self.firststdorient = True  # Used to decide whether to wipe the atomcoords clean
        self.cihamtyp = "none"  # Type of CI Hamiltonian: saps or dets.
        self.scftype = "none"  # Type of SCF calculation: BLYP, RHF, ROHF, etc.

    def extract(self, inputfile, line):
        """Extract information from the file object inputfile."""
        
        # Extract the version number. If the calculation is from
        # Firefly, its version number comes before a line that looks
        # like the normal GAMESS version number...
        if "Firefly version" in line:
            match = re.search(r"Firefly version\s([\d.]*)\D*(\d*)\s*\*", line)
            if match:
                version, build = match.groups()
                package_version = "{}.b{}".format(version, build)
                self.metadata["package_version"] = package_version
        if "GAMESS VERSION" in line:
            # ...so avoid overwriting it if Firefly already set this field.
            if "package_version" not in self.metadata:
                tokens = line.split()
                self.metadata["package_version"] = ' '.join(tokens[4:-1])

        if line[1:12] == "INPUT CARD>":
            return

        # extract the methods 
        if line[1:7] == "SCFTYP":
            method = line.split()[0][7:]
            if len(self.metadata["methods"]) == 0:
                self.metadata["methods"].append(method)

        # extract the basis set name
        if line[5:11] == "GBASIS":
            basnm1 = line.split()[0][7:]
            if basnm1 in self.dunningbas:
                self.metadata["basis_set"] = self.dunningbas[basnm1]
            else:
                if basnm1 == "PM3" or basnm1 == "AM1":
                    self.metadata["methods"].append(basnm1)
                if basnm1 == "STO" :
                    if line.split()[2] == "2":
                        self.metadata["basis_set"] = "STO-2G"
                    elif line.split()[2] == "3":
                        self.metadata["basis_set"] = "STO-3G"
                    elif line.split()[2] == "4":
                        self.metadata["basis_set"] = "STO-4G"
                    elif line.split()[2] == "5":
                        self.metadata["basis_set"] = "STO-5G"
                if basnm1 == "N21" :
                    if line.split()[2] == "3" and line.split()[3] == "POLAR=COMMON":
                        self.metadata["basis_set"] = "3-21G*"
                    if line.split()[2] == "3" and line.split()[3] == "POLAR=NONE":
                        self.metadata["basis_set"] = "3-21G"
                    if line.split()[2] == "4" and line.split()[3] == "POLAR=NONE":
                        self.metadata["basis_set"] = "4-21G"
                    if line.split()[2] == "6" and line.split()[3] == "POLAR=NONE":
                        self.metadata["basis_set"] = "6-21G"
                if basnm1 == "N31" :
                    if line.split()[2] == "6" and (line.split()[3] == "POLAR=POPN31" \
                            or line.split()[3] == "POLAR=POPLE"):
                        self.metadata["basis_set"] = "6-31G*"
                        line = next(inputfile)
                        if line.split()[-1] == "T":
                            self.metadata["basis_set"] = "6-31+G*"
                            line = next(inputfile)
                            if line.split()[1] == "0" and line.split()[3] == "T":
                                self.metadata["basis_set"] = "6-31++G*"
                            if line.split()[1] == "1" and line.split()[3] == "T":
                                self.metadata["basis_set"] = "6-31++G**"
                        else:
                            line = next(inputfile)
                            if line.split()[1] == "1":  #NPFUNC = 1
                                self.metadata["basis_set"] = "6-31G**"
                    if line.split()[2] == "6" and line.split()[3] == "POLAR=NONE":
                        self.metadata["basis_set"] = "6-31G"
                    if line.split()[2] == "4" and line.split()[3] == "POLAR=NONE":
                        self.metadata["basis_set"] = "4-31G"
                    if line.split()[2] == "4" and line.split()[3] == "POLAR=POPN31":
                        self.metadata["basis_set"] = "4-31G*"
                if basnm1 == "N311" :
                    if line.split()[2] == "6" and line.split()[3] == "POLAR=POPN311":
                        self.metadata["basis_set"] = "6-311G*"
                        line = next(inputfile)
                        if line.split()[-1] == "T":
                            self.metadata["basis_set"] = "6-311+G*"
                            line = next(inputfile)
                            if line.split()[1] == "0" and line.split()[3] == "T":
                                self.metadata["basis_set"] = "6-311++G*"
                            if line.split()[1] == "1" and line.split()[3] == "T":
                                self.metadata["basis_set"] = "6-311++G**"
                        else:
                            line = next(inputfile)
                            if line.split()[1] == "1":  #NPFUNC = 1
                                self.metadata["basis_set"] = "6-311G**"
                    if line.split()[2] == "6" and line.split()[3] == "POLAR=NONE":
                        self.metadata["basis_set"] = "6-311G"

        # We are looking for this line:
        #           PARAMETERS CONTROLLING GEOMETRY SEARCH ARE
        #           ...
        #           OPTTOL = 1.000E-04          RMIN   = 1.500E-03
        if line[10:18] == "OPTTOL =":
            if not hasattr(self, "geotargets"):
                opttol = float(line.split()[2])
                self.geotargets = numpy.array([opttol, 3. / opttol], "d")

        # Has to deal with such lines as:
        #  FINAL R-B3LYP ENERGY IS     -382.0507446475 AFTER  10 ITERATIONS
        #  FINAL ENERGY IS     -379.7594673378 AFTER   9 ITERATIONS
        # ...so take the number after the "IS"
        if line.find("FINAL") == 1:
            if not hasattr(self, "scfenergies"):
                self.scfenergies = []
            temp = line.split()
            self.scfenergies.append(utils.convertor(float(temp[temp.index("IS") + 1]), "hartree", "eV"))

        # For total energies after Moller-Plesset corrections, the output looks something like this:
        #
        # RESULTS OF MOLLER-PLESSET 2ND ORDER CORRECTION ARE
        #         E(0)=      -285.7568061536
        #         E(1)=         0.0
        #         E(2)=        -0.9679419329
        #       E(MP2)=      -286.7247480864
        # where E(MP2) = E(0) + E(2)
        #
        # With GAMESS-US 12 Jan 2009 (R3), the preceding text is different:
        #      DIRECT 4-INDEX TRANSFORMATION
        #      SCHWARZ INEQUALITY TEST SKIPPED          0 INTEGRAL BLOCKS
        #                     E(SCF)=       -76.0088477471
        #                       E(2)=        -0.1403745370
        #                     E(MP2)=       -76.1492222841
        #
        # With GAMESS-US 20 APR 2017 (R1), the following block may be present:
        #       SCHWARZ INEQUALITY TEST SKIPPED          0 INTEGRAL BLOCKS
        # ... END OF INTEGRAL TRANSFORMATION ...

        if line.find("RESULTS OF MOLLER-PLESSET") >= 0 or line[6:37] == "SCHWARZ INEQUALITY TEST SKIPPED":

            if not hasattr(self, "mpenergies"):
                self.mpenergies = []

            line = next(inputfile)
            # Each iteration has a new print-out
            if "END OF INTEGRAL TRANSFORMATION" not in line:
                self.mpenergies.append([])

                # GAMESS-US presently supports only second order corrections (MP2)
                # PC GAMESS also has higher levels (3rd and 4th), with different output
                # Only the highest level MP4 energy is gathered (SDQ or SDTQ)
                # Loop breaks when substring "DONE WITH MPn ENERGY" is encountered,
                # where n=2, 3 or 4.
                while "DONE WITH MP" not in line:

                    if len(line.split()) > 0:
                        # Only up to MP2 correction
                        if line.split()[0] == "E(MP2)=":
                            self.metadata["methods"].append("MP2")
                            mp2energy = float(line.split()[1])
                            self.mpenergies[-1].append(utils.convertor(mp2energy, "hartree", "eV"))
                        # MP2 before higher order calculations
                        if line.split()[0] == "E(MP2)":
                            mp2energy = float(line.split()[2])
                            self.mpenergies[-1].append(utils.convertor(mp2energy, "hartree", "eV"))
                        if line.split()[0] == "E(MP3)":
                            self.metadata["methods"].append("MP3")
                            mp3energy = float(line.split()[2])
                            self.mpenergies[-1].append(utils.convertor(mp3energy, "hartree", "eV"))
                        if line.split()[0] in ["E(MP4-SDQ)", "E(MP4-SDTQ)"]:
                            self.metadata["methods"].append("MP4")
                            mp4energy = float(line.split()[2])
                            self.mpenergies[-1].append(utils.convertor(mp4energy, "hartree", "eV"))
                    line = next(inputfile)

        # Total energies after Coupled Cluster calculations
        # Only the highest Coupled Cluster level result is gathered
        if line[12:23] == "CCD ENERGY:":
            self.metadata["methods"].append("CCD")
            if not hasattr(self, "ccenergies"):
                self.ccenergies = []
            ccenergy = float(line.split()[2])
            self.ccenergies.append(utils.convertor(ccenergy, "hartree", "eV"))
        if line.find("CCSD") >= 0 and line.split()[0:2] == ["CCSD", "ENERGY:"]:
            self.metadata["methods"].append("CCSD")
            if not hasattr(self, "ccenergies"):
                self.ccenergies = []
            ccenergy = float(line.split()[2])
            line = next(inputfile)
            if line[8:23] == "CCSD[T] ENERGY:":
                self.metadata["methods"].append("CCSD[T]")
                ccenergy = float(line.split()[2])
                line = next(inputfile)
                if line[8:23] == "CCSD(T) ENERGY:":
                    self.metadata["methods"].append("CCSD(T)")
                    ccenergy = float(line.split()[2])
            self.ccenergies.append(utils.convertor(ccenergy, "hartree", "eV"))

        # Also collect MP2 energies, which are always calculated before CC
        if line[8:23] == "MBPT(2) ENERGY:":
            if not hasattr(self, "mpenergies"):
                self.mpenergies = []
            self.mpenergies.append([])
            mp2energy = float(line.split()[2])
            self.mpenergies[-1].append(utils.convertor(mp2energy, "hartree", "eV"))

        # Extract charge and multiplicity
        if line[1:19] == "CHARGE OF MOLECULE":

            charge = int(round(float(line.split()[-1])))
            self.set_attribute('charge', charge)

            line = next(inputfile)
            mult = int(line.split()[-1])
            self.set_attribute('mult', mult)

        # Electronic transitions (etenergies) for CIS runs and TD-DFT, which
        # have very similar outputs. The outputs EOM look very differentm, though.
        #
        #  ---------------------------------------------------------------------
        #                    CI-SINGLES EXCITATION ENERGIES
        #  STATE       HARTREE        EV      KCAL/MOL       CM-1         NM
        #  ---------------------------------------------------------------------
        #   1A''   0.1677341781     4.5643    105.2548      36813.40     271.64
        #   ...
        if re.match("(CI-SINGLES|TDDFT) EXCITATION ENERGIES", line.strip()):

            if not hasattr(self, "etenergies"):
                self.etenergies = []

            get_etosc = False
            header = next(inputfile).rstrip()
            if header.endswith("OSC. STR."):
                # water_cis_dets.out does not have the oscillator strength
                # in this table...it is extracted from a different section below
                get_etosc = True
                self.etoscs = []

            self.skip_line(inputfile, 'dashes')

            line = next(inputfile)
            broken = line.split()
            while len(broken) > 0:

                # Take hartree value with more numbers, and convert.
                # Note that the values listed after this are also less exact!
                etenergy = float(broken[1])
                self.etenergies.append(utils.convertor(etenergy, "hartree", "wavenumber"))
                if get_etosc:
                    etosc = float(broken[-1])
                    self.etoscs.append(etosc)
                broken = next(inputfile).split()

        # Detect the CI hamiltonian type, if applicable.
        # Should always be detected if CIS is done.
        if line[8:64] == "RESULTS FROM SPIN-ADAPTED ANTISYMMETRIZED PRODUCT (SAPS)":
            self.cihamtyp = "saps"
        if line[8:64] == "RESULTS FROM DETERMINANT BASED ATOMIC ORBITAL CI-SINGLES":
            self.cihamtyp = "dets"

        # etsecs (used only for CIS runs for now)
        if line[1:14] == "EXCITED STATE":
            if not hasattr(self, 'etsecs'):
                self.etsecs = []
            if not hasattr(self, 'etsyms'):
                self.etsyms = []
            statenumber = int(line.split()[2])
            spin = int(float(line.split()[7]))
            if spin == 0:
                sym = "Singlet"
            if spin == 1:
                sym = "Triplet"
            sym += '-' + line.split()[-1]
            self.etsyms.append(sym)
            # skip 5 lines
            for i in range(5):
                line = next(inputfile)
            line = next(inputfile)
            CIScontribs = []
            while line.strip()[0] != "-":
                MOtype = 0
                # alpha/beta are specified for hamtyp=dets
                if self.cihamtyp == "dets":
                    if line.split()[0] == "BETA":
                        MOtype = 1
                fromMO = int(line.split()[-3])-1
                toMO = int(line.split()[-2])-1
                coeff = float(line.split()[-1])
                # With the SAPS hamiltonian, the coefficients are multiplied
                #   by sqrt(2) so that they normalize to 1.
                # With DETS, both alpha and beta excitations are printed.
                # if self.cihamtyp == "saps":
                #    coeff /= numpy.sqrt(2.0)
                CIScontribs.append([(fromMO, MOtype), (toMO, MOtype), coeff])
                line = next(inputfile)
            self.etsecs.append(CIScontribs)

        # etoscs (used only for CIS runs now)
        if line[1:50] == "TRANSITION FROM THE GROUND STATE TO EXCITED STATE":
            if not hasattr(self, "etoscs"):
                self.etoscs = []

            # This was the suggested as a fix in issue #61, and it does allow
            # the parser to finish without crashing. However, it seems that
            # etoscs is shorter in this case than the other transition attributes,
            # so that should be somehow corrected and tested for.
            if "OPTICALLY" in line:
                pass
            else:
                statenumber = int(line.split()[-1])
                # skip 7 lines
                for i in range(8):
                    line = next(inputfile)
                strength = float(line.split()[3])
                self.etoscs.append(strength)

        # TD-DFT for GAMESS-US.
        # The format for excitations has changed a bit between 2007 and 2012.
        # Original format parser was written for:
        #
        #          -------------------
        #          TRIPLET EXCITATIONS
        #          -------------------
        #
        # STATE #   1  ENERGY =    3.027228 EV
        # OSCILLATOR STRENGTH =    0.000000
        #        DRF    COEF       OCC      VIR
        #        ---    ----       ---      ---
        #         35 -1.105383     35  ->   36
        #         69 -0.389181     34  ->   37
        #        103 -0.405078     33  ->   38
        #        137  0.252485     32  ->   39
        #        168 -0.158406     28  ->   40
        #
        # STATE #   2  ENERGY =    4.227763 EV
        # ...
        #
        # Here is the corresponding 2012 version:
        #
        #          -------------------
        #          TRIPLET EXCITATIONS
        #          -------------------
        #
        # STATE #   1  ENERGY =    3.027297 EV
        # OSCILLATOR STRENGTH =    0.000000
        # LAMBDA DIAGNOSTIC   =    0.925 (RYDBERG/CHARGE TRANSFER CHARACTER)
        # SYMMETRY OF STATE   =    A
        #                 EXCITATION  DE-EXCITATION
        #     OCC     VIR  AMPLITUDE      AMPLITUDE
        #      I       A     X(I->A)        Y(A->I)
        #     ---     ---   --------       --------
        #     35      36   -0.929190      -0.176167
        #     34      37   -0.279823      -0.109414
        # ...
        #
        # We discern these two by the presence of the arrow in the old version.
        #
        # The "LET EXCITATIONS" pattern used below catches both
        # singlet and triplet excitations output.
        if line[14:29] == "LET EXCITATIONS":

            self.etenergies = []
            self.etoscs = []
            self.etsecs = []
            etsyms = []

            self.skip_lines(inputfile, ['d', 'b'])

            # Loop while states are still being printed.
            line = next(inputfile)
            while line[1:6] == "STATE":

                self.updateprogress(inputfile, "Excited States")

                etenergy = utils.convertor(float(line.split()[-2]), "eV", "wavenumber")
                etoscs = float(next(inputfile).split()[-1])
                self.etenergies.append(etenergy)
                self.etoscs.append(etoscs)

                # Symmetry is not always present, especially in old versions.
                # Newer versions, on the other hand, can also provide a line
                # with lambda diagnostic and some extra headers.
                line = next(inputfile)
                if "LAMBDA DIAGNOSTIC" in line:
                    line = next(inputfile)
                if "SYMMETRY" in line:
                    etsyms.append(line.split()[-1])
                    line = next(inputfile)
                if "EXCITATION" in line and "DE-EXCITATION" in line:
                    line = next(inputfile)
                if line.count("AMPLITUDE") == 2:
                    line = next(inputfile)

                self.skip_line(inputfile, 'dashes')

                CIScontribs = []
                line = next(inputfile)
                while line.strip():
                    cols = line.split()
                    if "->" in line:
                        i_occ_vir = [2, 4]
                        i_coeff = 1

                    else:
                        i_occ_vir = [0, 1]
                        i_coeff = 2
                    fromMO, toMO = [int(cols[i]) - 1 for i in i_occ_vir]
                    coeff = float(cols[i_coeff])
                    CIScontribs.append([(fromMO, 0), (toMO, 0), coeff])
                    line = next(inputfile)
                self.etsecs.append(CIScontribs)
                line = next(inputfile)

            # The symmetries are not always present.
            if etsyms:
                self.etsyms = etsyms

        # Maximum and RMS gradients.
        if "MAXIMUM GRADIENT" in line or "RMS GRADIENT" in line:

            parts = line.split()

            # Avoid parsing the following...

            ## YOU SHOULD RESTART "OPTIMIZE" RUNS WITH THE COORDINATES
            ## WHOSE ENERGY IS LOWEST.  RESTART "SADPOINT" RUNS WITH THE
            ## COORDINATES WHOSE RMS GRADIENT IS SMALLEST.  THESE ARE NOT
            ## ALWAYS THE LAST POINT COMPUTED!

            if parts[0] not in ["MAXIMUM", "RMS", "(1)"]:
                return

            if not hasattr(self, "geovalues"):
                self.geovalues = []

            # Newer versions (around 2006) have both maximum and RMS on one line:
            #       MAXIMUM GRADIENT =  0.0531540    RMS GRADIENT = 0.0189223
            if len(parts) == 8:
                maximum = float(parts[3])
                rms = float(parts[7])

            # In older versions of GAMESS, this spanned two lines, like this:
            #       MAXIMUM GRADIENT =    0.057578167
            #           RMS GRADIENT =    0.027589766
            if len(parts) == 4:
                maximum = float(parts[3])
                line = next(inputfile)
                parts = line.split()
                rms = float(parts[3])

            # FMO also prints two final one- and two-body gradients (see exam37):
            #   (1) MAXIMUM GRADIENT =  0.0531540    RMS GRADIENT = 0.0189223
            if len(parts) == 9:
                maximum = float(parts[4])
                rms = float(parts[8])

            self.geovalues.append([maximum, rms])

        # This is the input orientation, which is the only data available for
        # SP calcs, but which should be overwritten by the standard orientation
        # values, which is the only information available for all geoopt cycles.
        if line[11:50] == "ATOMIC                      COORDINATES":

            if not hasattr(self, "atomcoords"):
                self.atomcoords = []

            line = next(inputfile)
            atomcoords = []
            atomnos = []
            line = next(inputfile)
            while line.strip():
                temp = line.strip().split()
                atomcoords.append([utils.convertor(float(x), "bohr", "Angstrom") for x in temp[2:5]])
                atomnos.append(int(round(float(temp[1]))))  # Don't use the atom name as this is arbitary
                line = next(inputfile)

            self.set_attribute('atomnos', atomnos)
            self.atomcoords.append(atomcoords)

        if line[12:40] == "EQUILIBRIUM GEOMETRY LOCATED":
            # Prevent extraction of the final geometry twice
            if not hasattr(self, 'optdone'):
                self.optdone = []
            self.optdone.append(len(self.geovalues) - 1)

        # Make sure we always have optdone for geomtry optimization, even if not converged.
        if "GEOMETRY SEARCH IS NOT CONVERGED" in line:
            if not hasattr(self, 'optdone'):
                self.optdone = []

        # This is the standard orientation, which is the only coordinate
        # information available for all geometry optimisation cycles.
        # The input orientation will be overwritten if this is a geometry optimisation
        # We assume that a previous Input Orientation has been found and
        # used to extract the atomnos
        if line[1:29] == "COORDINATES OF ALL ATOMS ARE" and (not hasattr(self, "optdone") or self.optdone == []):

            self.updateprogress(inputfile, "Coordinates")

            if self.firststdorient:
                self.firststdorient = False
                # Wipes out the single input coordinate at the start of the file
                self.atomcoords = []

            self.skip_lines(inputfile, ['line', '-'])

            atomcoords = []
            line = next(inputfile)

            for i in range(self.natom):
                temp = line.strip().split()
                atomcoords.append(list(map(float, temp[2:5])))
                line = next(inputfile)
            self.atomcoords.append(atomcoords)

        # Section with SCF information.
        #
        # The space at the start of the search string is to differentiate from MCSCF.
        # Everything before the search string is stored as the type of SCF.
        # SCF types may include: BLYP, RHF, ROHF, UHF, etc.
        #
        # For example, in exam17 the section looks like this (note that this is GVB):
        #          ------------------------
        #          ROHF-GVB SCF CALCULATION
        #          ------------------------
        # GVB STEP WILL USE    119875 WORDS OF MEMORY.
        #
        #     MAXIT=  30   NPUNCH= 2   SQCDF TOL=1.0000E-05
        #     NUCLEAR ENERGY=        6.1597411978
        #     EXTRAP=T   DAMP=F   SHIFT=F   RSTRCT=F   DIIS=F  SOSCF=F
        #
        # ITER EX     TOTAL ENERGY       E CHANGE        SQCDF       DIIS ERROR
        #   0  0      -38.298939963   -38.298939963   0.131784454   0.000000000
        #   1  1      -38.332044339    -0.033104376   0.026019716   0.000000000
        # ... and will be terminated by a blank line.
        if line.rstrip()[-16:] == " SCF CALCULATION":

            # Remember the type of SCF.
            self.scftype = line.strip()[:-16]

            self.skip_line(inputfile, 'dashes')

            while line[:5] != " ITER":

                self.updateprogress(inputfile, "Attributes")

                # GVB uses SQCDF for checking convergence (for example in exam17).
                if "GVB" in self.scftype and "SQCDF TOL=" in line:
                    scftarget = float(line.split("=")[-1])

                # Normally, however, the density is used as the convergence criterium.
                # Deal with various versions:
                #   (GAMESS VERSION = 12 DEC 2003)
                #     DENSITY MATRIX CONV=  2.00E-05  DFT GRID SWITCH THRESHOLD=  3.00E-04
                #   (GAMESS VERSION = 22 FEB 2006)
                #     DENSITY MATRIX CONV=  1.00E-05
                #   (PC GAMESS version 6.2, Not DFT?)
                #     DENSITY CONV=  1.00E-05
                elif "DENSITY CONV" in line or "DENSITY MATRIX CONV" in line:
                    scftarget = float(line.split()[-1])

                line = next(inputfile)

            if not hasattr(self, "scftargets"):
                self.scftargets = []

            self.scftargets.append([scftarget])

            if not hasattr(self, "scfvalues"):
                self.scfvalues = []

            # Normally the iterations print in 6 columns.
            # For ROHF, however, it is 5 columns, thus this extra parameter.
            if "ROHF" in self.scftype:
                self.scf_valcol = 4
            else:
                self.scf_valcol = 5

            line = next(inputfile)

            # SCF iterations are terminated by a blank line.
            # The first four characters usually contains the step number.
            # However, lines can also contain messages, including:
            #   * * *   INITIATING DIIS PROCEDURE   * * *
            #   CONVERGED TO SWOFF, SO DFT CALCULATION IS NOW SWITCHED ON
            #   DFT CODE IS SWITCHING BACK TO THE FINER GRID
            values = []
            while line.strip():
                try:
                    temp = int(line[0:4])
                except ValueError:
                    pass
                else:
                    values.append([float(line.split()[self.scf_valcol])])
                try:
                    line = next(inputfile)
                except StopIteration:
                    self.logger.warning('File terminated before end of last SCF!')
                    break
            self.scfvalues.append(values)

        # Sometimes, only the first SCF cycle has the banner parsed for above,
        # so we must identify them from the header before the SCF iterations.
        # The example we have for this is the GeoOpt unittest for Firefly8.
        if line[1:8] == "ITER EX":

            # In this case, the convergence targets are not printed, so we assume
            # they do not change.
            self.scftargets.append(self.scftargets[-1])

            values = []
            line = next(inputfile)
            while line.strip():
                try:
                    temp = int(line[0:4])
                except ValueError:
                    pass
                else:
                    values.append([float(line.split()[self.scf_valcol])])
                line = next(inputfile)
            self.scfvalues.append(values)

        # Extract normal coordinate analysis, including vibrational frequencies (vibfreq),
        # IT intensities (vibirs) and displacements (vibdisps).
        #
        # This section typically looks like the following in GAMESS-US:
        #
        # MODES 1 TO 6 ARE TAKEN AS ROTATIONS AND TRANSLATIONS.
        #
        #     FREQUENCIES IN CM**-1, IR INTENSITIES IN DEBYE**2/AMU-ANGSTROM**2,
        #     REDUCED MASSES IN AMU.
        #
        #                          1           2           3           4           5
        #       FREQUENCY:        52.49       41.45       17.61        9.23       10.61
        #    REDUCED MASS:      3.92418     3.77048     5.43419     6.44636     5.50693
        #    IR INTENSITY:      0.00013     0.00001     0.00004     0.00000     0.00003
        #
        # ...or in the case of a numerical Hessian job...
        #
        # MODES 1 TO 5 ARE TAKEN AS ROTATIONS AND TRANSLATIONS.
        #
        #     FREQUENCIES IN CM**-1, IR INTENSITIES IN DEBYE**2/AMU-ANGSTROM**2,
        #     REDUCED MASSES IN AMU.
        #
        #                          1           2           3           4           5
        #       FREQUENCY:         0.05        0.03        0.03       30.89       30.94
        #    REDUCED MASS:      8.50125     8.50137     8.50136     1.06709     1.06709
        #
        # ...whereas PC-GAMESS has...
        #
        # MODES 1 TO 6 ARE TAKEN AS ROTATIONS AND TRANSLATIONS.
        #
        #     FREQUENCIES IN CM**-1, IR INTENSITIES IN DEBYE**2/AMU-ANGSTROM**2
        #
        #                          1           2           3           4           5
        #       FREQUENCY:         5.89        1.46        0.01        0.01        0.01
        #    IR INTENSITY:      0.00000     0.00000     0.00000     0.00000     0.00000
        #
        # If Raman is present we have (for PC-GAMESS)...
        #
        # MODES 1 TO 6 ARE TAKEN AS ROTATIONS AND TRANSLATIONS.
        #
        #     FREQUENCIES IN CM**-1, IR INTENSITIES IN DEBYE**2/AMU-ANGSTROM**2
        #     RAMAN INTENSITIES IN ANGSTROM**4/AMU, DEPOLARIZATIONS ARE DIMENSIONLESS
        #
        #                          1           2           3           4           5
        #       FREQUENCY:         5.89        1.46        0.04        0.03        0.01
        #    IR INTENSITY:      0.00000     0.00000     0.00000     0.00000     0.00000
        # RAMAN INTENSITY:       12.675       1.828       0.000       0.000       0.000
        #  DEPOLARIZATION:        0.750       0.750       0.124       0.009       0.750
        #
        # If GAMESS-US or PC-GAMESS has not reached the stationary point we have
        # and additional warning, repeated twice, like so (see n_water.log for an example):
        #
        #     *******************************************************
        #     * THIS IS NOT A STATIONARY POINT ON THE MOLECULAR PES *
        #     *     THE VIBRATIONAL ANALYSIS IS NOT VALID !!!       *
        #     *******************************************************
        #
        # There can also be additional warnings about the selection of modes, for example:
        #
        # * * * WARNING, MODE 6 HAS BEEN CHOSEN AS A VIBRATION
        #          WHILE MODE12 IS ASSUMED TO BE A TRANSLATION/ROTATION.
        # PLEASE VERIFY THE PROGRAM'S DECISION MANUALLY!
        #
        if "NORMAL COORDINATE ANALYSIS IN THE HARMONIC APPROXIMATION" in line:

            self.vibfreqs = []
            self.vibirs = []
            self.vibdisps = []

            # Need to get to the modes line, which is often preceeded by
            # a list of atomic weights and some possible warnings.
            # Pass the warnings to the logger if they are there.
            while not "MODES" in line:
                self.updateprogress(inputfile, "Frequency Information")

                line = next(inputfile)

                # Typical Atomic Masses section printed in GAMESS
                #               ATOMIC WEIGHTS (AMU)
                #
                # 1     O                15.99491
                # 2     H                 1.00782
                # 3     H                 1.00782
                if "ATOMIC WEIGHTS" in line:
                    atommasses = []
                    self.skip_line(inputfile,['b'])
                    # There is a blank line after ATOMIC WEIGHTS
                    line = next(inputfile)
                    while line.strip():
                        temp = line.strip().split()
                        atommasses.append(float(temp[2]))
                        line = next(inputfile)
                    self.set_attribute('atommasses', atommasses)

                if "THIS IS NOT A STATIONARY POINT" in line:
                    msg = "\n   This is not a stationary point on the molecular PES"
                    msg += "\n   The vibrational analysis is not valid!!!"
                    self.logger.warning(msg)
                if "* * * WARNING, MODE" in line:
                    line1 = line.strip()
                    line2 = next(inputfile).strip()
                    line3 = next(inputfile).strip()
                    self.logger.warning("\n   " + "\n   ".join((line1, line2, line3)))

            # In at least one case (regression zolm_dft3a.log) for older version of GAMESS-US,
            # the header concerning the range of nodes is formatted wrong and can look like so:
            # MODES 9 TO14 ARE TAKEN AS ROTATIONS AND TRANSLATIONS.
            #  ... although it's unclear whether this happens for all two-digit values.
            startrot = int(line.split()[1])
            if len(line.split()[2]) == 2:
                endrot = int(line.split()[3])
            else:
                endrot = int(line.split()[2][2:])

            self.skip_line(inputfile, 'blank')

            # Continue down to the first frequencies
            line = next(inputfile)
            # With GAMESS-US 20 APR 2017 (R1), there are 28 blank spaces,
            # in earlier versions there used to be 26.
            while not line.strip() or not re.search(' {26,}1', line) is not None:
                line = next(inputfile)

            while not "SAYVETZ" in line:
                self.updateprogress(inputfile, "Frequency Information")

                # Note: there may be imaginary frequencies like this (which we make negative):
                #       FREQUENCY:       825.18 I    111.53       12.62       10.70        0.89
                #
                # A note for debuggers: some of these frequencies will be removed later,
                # assumed to be translations or rotations (see startrot/endrot above).
                for col in next(inputfile).split()[1:]:
                    if col == "I":
                        self.vibfreqs[-1] *= -1
                    else:
                        self.vibfreqs.append(float(col))

                line = next(inputfile)

                # Skip the symmetry (appears in newer versions), fixes bug #3476063.
                if line.find("SYMMETRY") >= 0:
                    line = next(inputfile)

                # Skip the reduced mass (not always present).
                if line.find("REDUCED") >= 0:
                    line = next(inputfile)

                # Not present in numerical Hessian calculations.
                if line.find("IR INTENSITY") >= 0:
                    irIntensity = map(float, line.strip().split()[2:])
                    self.vibirs.extend([utils.convertor(x, "Debye^2/amu-Angstrom^2", "km/mol") for x in irIntensity])
                    line = next(inputfile)

                # Read in Raman vibrational intensities if present.
                if line.find("RAMAN") >= 0:
                    if not hasattr(self, "vibramans"):
                        self.vibramans = []
                    ramanIntensity = line.strip().split()
                    self.vibramans.extend(list(map(float, ramanIntensity[2:])))
                    depolar = next(inputfile)
                    line = next(inputfile)

                # This line seems always to be blank.
                assert line.strip() == ''

                # Extract the Cartesian displacement vectors.
                p = [[], [], [], [], []]
                for j in range(self.natom):
                    q = [[], [], [], [], []]
                    for coord in "xyz":
                        line = next(inputfile)[21:]
                        cols = list(map(float, line.split()))
                        for i, val in enumerate(cols):
                            q[i].append(val)
                    for k in range(len(cols)):
                        p[k].append(q[k])
                self.vibdisps.extend(p[:len(cols)])

                # Skip the Sayvetz stuff at the end.
                for j in range(10):
                    line = next(inputfile)

                self.skip_line(inputfile, 'blank')
                line = next(inputfile)

            # Exclude rotations and translations.
            self.vibfreqs = numpy.array(self.vibfreqs[:startrot-1]+self.vibfreqs[endrot:], "d")
            self.vibirs = numpy.array(self.vibirs[:startrot-1]+self.vibirs[endrot:], "d")
            self.vibdisps = numpy.array(self.vibdisps[:startrot-1]+self.vibdisps[endrot:], "d")
            if hasattr(self, "vibramans"):
                self.vibramans = numpy.array(self.vibramans[:startrot-1]+self.vibramans[endrot:], "d")

        if line[5:21] == "ATOMIC BASIS SET":
            self.gbasis = []
            line = next(inputfile)
            while line.find("SHELL") < 0:
                line = next(inputfile)

            self.skip_lines(inputfile, ['blank', 'atomname'])

            # shellcounter stores the shell no of the last shell
            # in the previous set of primitives
            shellcounter = 1
            while line.find("TOTAL NUMBER") < 0:

                self.skip_line(inputfile, 'blank')

                line = next(inputfile)
                shellno = int(line.split()[0])
                shellgap = shellno - shellcounter
                gbasis = []  # Stores basis sets on one atom
                shellsize = 0
                while len(line.split()) != 1 and line.find("TOTAL NUMBER") < 0:
                    shellsize += 1
                    coeff = {}
                    # coefficients and symmetries for a block of rows
                    while line.strip():
                        temp = line.strip().split()
                        sym = temp[1]
                        assert sym in ['S', 'P', 'D', 'F', 'G', 'L']
                        if sym == "L":  # L refers to SP
                            if len(temp) == 6:  # GAMESS US
                                coeff.setdefault("S", []).append((float(temp[3]), float(temp[4])))
                                coeff.setdefault("P", []).append((float(temp[3]), float(temp[5])))
                            else:  # PC GAMESS
                                assert temp[6][-1] == temp[9][-1] == ')'
                                coeff.setdefault("S", []).append((float(temp[3]), float(temp[6][:-1])))
                                coeff.setdefault("P", []).append((float(temp[3]), float(temp[9][:-1])))
                        else:
                            if len(temp) == 5:  # GAMESS US
                                coeff.setdefault(sym, []).append((float(temp[3]), float(temp[4])))
                            else:  # PC GAMESS
                                assert temp[6][-1] == ')'
                                coeff.setdefault(sym, []).append((float(temp[3]), float(temp[6][:-1])))
                        line = next(inputfile)
                    # either a blank or a continuation of the block
                    if sym == "L":
                        gbasis.append(('S', coeff['S']))
                        gbasis.append(('P', coeff['P']))
                    else:
                        gbasis.append((sym, coeff[sym]))
                    line = next(inputfile)
                # either the start of the next block or the start of a new atom or
                # the end of the basis function section

                numtoadd = 1 + (shellgap // shellsize)
                shellcounter = shellno + shellsize
                for x in range(numtoadd):
                    self.gbasis.append(gbasis)

        # The eigenvectors, which also include MO energies and symmetries, follow
        # the *final* report of evalues and the last list of symmetries in the log file:
        #
        #           ------------
        #           EIGENVECTORS
        #           ------------
        #
        #                       1          2          3          4          5
        #                   -10.0162   -10.0161   -10.0039   -10.0039   -10.0029
        #                      BU         AG         BU         AG         AG
        #     1  C  1  S    0.699293   0.699290  -0.027566   0.027799   0.002412
        #     2  C  1  S    0.031569   0.031361   0.004097  -0.004054  -0.000605
        #     3  C  1  X    0.000908   0.000632  -0.004163   0.004132   0.000619
        #     4  C  1  Y   -0.000019   0.000033   0.000668  -0.000651   0.005256
        #     5  C  1  Z    0.000000   0.000000   0.000000   0.000000   0.000000
        #     6  C  2  S   -0.699293   0.699290   0.027566   0.027799   0.002412
        #     7  C  2  S   -0.031569   0.031361  -0.004097  -0.004054  -0.000605
        #     8  C  2  X    0.000908  -0.000632  -0.004163  -0.004132  -0.000619
        #     9  C  2  Y   -0.000019  -0.000033   0.000668   0.000651  -0.005256
        #    10  C  2  Z    0.000000   0.000000   0.000000   0.000000   0.000000
        #    11  C  3  S   -0.018967  -0.019439   0.011799  -0.014884  -0.452328
        #    12  C  3  S   -0.007748  -0.006932   0.000680  -0.000695  -0.024917
        #    13  C  3  X    0.002628   0.002997   0.000018   0.000061  -0.003608
        # ...
        #
        # There are blanks lines between each block.
        #
        # Warning! There are subtle differences between GAMESS-US and PC-GAMES
        # in the formatting of the first four columns. In particular, for F orbitals,
        # PC GAMESS:
        #   19  C   1 YZ   0.000000   0.000000   0.000000   0.000000   0.000000
        #   20  C    XXX   0.000000   0.000000   0.000000   0.000000   0.002249
        #   21  C    YYY   0.000000   0.000000  -0.025555   0.000000   0.000000
        #   22  C    ZZZ   0.000000   0.000000   0.000000   0.002249   0.000000
        #   23  C    XXY   0.000000   0.000000   0.001343   0.000000   0.000000
        # GAMESS US
        #   55  C  1 XYZ   0.000000   0.000000   0.000000   0.000000   0.000000
        #   56  C  1XXXX  -0.000014  -0.000067   0.000000   0.000000   0.000000
        #
        if line.find("EIGENVECTORS") == 10 or line.find("MOLECULAR ORBITALS") == 10:

            # This is the stuff that we can read from these blocks.
            self.moenergies = [[]]
            self.mosyms = [[]]

            if not hasattr(self, "nmo"):
                self.nmo = self.nbasis

            self.mocoeffs = [numpy.zeros((self.nmo, self.nbasis), "d")]

            readatombasis = False
            if not hasattr(self, "atombasis"):
                self.atombasis = []
                self.aonames = []
                for i in range(self.natom):
                    self.atombasis.append([])
                self.aonames = []
                readatombasis = True

            self.skip_line(inputfile, 'dashes')

            for base in range(0, self.nmo, 5):
                self.updateprogress(inputfile, "Coefficients")

                line = next(inputfile)

                # This makes sure that this section does not end prematurely,
                # which happens in regression 2CO.ccsd.aug-cc-pVDZ.out.
                if line.strip() != "":
                    break

                numbers = next(inputfile)  # Eigenvector numbers.

                # This is for regression CdtetraM1B3LYP.
                if "ALPHA SET" in numbers:
                    blank = next(inputfile)
                    numbers = next(inputfile)

                # If not all coefficients are printed, the logfile will go right to
                # the beta section if there is one, so break out in that case.
                if "BETA SET" in numbers:
                    line = numbers
                    break

                # Sometimes there are some blank lines here.
                while not line.strip():
                    line = next(inputfile)

                # Geometry optimizations don't have END OF RHF/DFT
                # CALCULATION, they head right into the next section.
                if "--------" in line:
                    break

                # Eigenvalues for these orbitals (in hartrees).
                try:
                    self.moenergies[0].extend([utils.convertor(float(x), "hartree", "eV") for x in line.split()])
                except:
                    self.logger.warning('MO section found but could not be parsed!')
                    break

                # Orbital symmetries.
                line = next(inputfile)
                if line.strip():
                    self.mosyms[0].extend(list(map(self.normalisesym, line.split())))

                # Now we have nbasis lines. We will use the same method as in normalise_aonames() before.
                p = re.compile("(\d+)\s*([A-Z][A-Z]?)\s*(\d+)\s*([A-Z]+)")
                oldatom = '0'
                i_atom = 0  # counter to keep track of n_atoms > 99
                flag_w = True  # flag necessary to keep from adding 100's at wrong time

                for i in range(self.nbasis):
                    line = next(inputfile)

                    # If line is empty, break (ex. for FMO in exam37 which is a regression).
                    if not line.strip():
                        break

                    # Fill atombasis and aonames only first time around
                    if readatombasis and base == 0:

                        aonames = []
                        start = line[:17].strip()
                        m = p.search(start)

                        if m:
                            g = m.groups()
                            g2 = int(g[2])  # atom index in GAMESS file; changes to 0 after 99

                            # Check if we have moved to a hundred
                            # if so, increment the counter and add it to the parsed value
                            # There will be subsequent 0's as that atoms AO's are parsed
                            # so wait until the next atom is parsed before resetting flag
                            if g2 == 0 and flag_w:
                                i_atom = i_atom + 100
                                flag_w = False  # handle subsequent AO's
                            if g2 != 0:
                                flag_w = True  # reset flag
                            g2 = g2 + i_atom

                            aoname = "%s%i_%s" % (g[1].capitalize(), g2, g[3])
                            oldatom = str(g2)
                            atomno = g2-1
                            orbno = int(g[0])-1
                        else:  # For F orbitals, as shown above
                            g = [x.strip() for x in line.split()]
                            aoname = "%s%s_%s" % (g[1].capitalize(), oldatom, g[2])
                            atomno = int(oldatom)-1
                            orbno = int(g[0])-1

                        self.atombasis[atomno].append(orbno)
                        self.aonames.append(aoname)

                    coeffs = line[15:]  # Strip off the crud at the start.
                    j = 0
                    while j*11+4 < len(coeffs):
                        self.mocoeffs[0][base+j, i] = float(coeffs[j * 11:(j + 1) * 11])
                        j += 1

            # If it's a restricted calc and no more properties, we have:
            #
            #  ...... END OF RHF/DFT CALCULATION ......
            #
            # If there are more properties (such as the density matrix):
            #               --------------
            #
            # If it's an unrestricted calculation, however, we now get the beta orbitals:
            #
            #  ----- BETA SET -----
            #
            #          ------------
            #          EIGENVECTORS
            #          ------------
            #
            #                      1          2          3          4          5
            # ...
            #
            if "BETA SET" not in line:
                line = next(inputfile)
                line = next(inputfile)

            # This can come in between the alpha and beta orbitals (see #130).
            if line.strip() == "LZ VALUE ANALYSIS FOR THE MOS":
                while line.strip():
                    line = next(inputfile)
                line = next(inputfile)

            # Covers label with both dashes and stars (like regression CdtetraM1B3LYP).
            if "BETA SET" in line:
                self.mocoeffs.append(numpy.zeros((self.nmo, self.nbasis), "d"))
                self.moenergies.append([])
                self.mosyms.append([])
                blank = next(inputfile)
                line = next(inputfile)

                # Sometimes EIGENVECTORS is missing, so rely on dashes to signal it.
                if set(line.strip()) == {'-'}:
                    self.skip_lines(inputfile, ['EIGENVECTORS', 'd', 'b'])
                    line = next(inputfile)

                for base in range(0, self.nmo, 5):
                    self.updateprogress(inputfile, "Coefficients")
                    if base != 0:
                        line = next(inputfile)
                        line = next(inputfile)
                    line = next(inputfile)
                    if "properties" in line.lower():
                        break
                    self.moenergies[1].extend([utils.convertor(float(x), "hartree", "eV") for x in line.split()])
                    line = next(inputfile)
                    self.mosyms[1].extend(list(map(self.normalisesym, line.split())))
                    for i in range(self.nbasis):
                        line = next(inputfile)
                        temp = line[15:]  # Strip off the crud at the start
                        j = 0
                        while j * 11 + 4 < len(temp):
                            self.mocoeffs[1][base+j, i] = float(temp[j * 11:(j + 1) * 11])
                            j += 1
                line = next(inputfile)
            self.moenergies = [numpy.array(x, "d") for x in self.moenergies]

        # Natural orbital coefficients and occupation numbers, presently supported only
        # for CIS calculations. Looks the same as eigenvectors, without symmetry labels.
        #
        #          --------------------
        #          CIS NATURAL ORBITALS
        #          --------------------
        #
        #                      1          2          3          4          5
        #
        #                    2.0158     2.0036     2.0000     2.0000     1.0000
        #
        #    1  O  1  S    0.000000  -0.157316   0.999428   0.164938   0.000000
        #    2  O  1  S    0.000000   0.754402   0.004472  -0.581970   0.000000
        # ...
        #
        if line[10:30] == "CIS NATURAL ORBITALS":

            self.nocoeffs = numpy.zeros((self.nmo, self.nbasis), "d")
            self.nooccnos = []

            self.skip_line(inputfile, 'dashes')

            for base in range(0, self.nmo, 5):

                self.skip_lines(inputfile, ['blank', 'numbers'])

                # The eigenvalues that go along with these natural orbitals are
                # their occupation numbers. Sometimes there are blank lines before them.
                line = next(inputfile)
                while not line.strip():
                    line = next(inputfile)
                eigenvalues = map(float, line.split())
                self.nooccnos.extend(eigenvalues)

                # Orbital symemtry labels are normally here for MO coefficients.
                line = next(inputfile)

                # Now we have nbasis lines with the coefficients.
                for i in range(self.nbasis):

                    line = next(inputfile)
                    coeffs = line[15:]
                    j = 0
                    while j*11+4 < len(coeffs):
                        self.nocoeffs[base+j, i] = float(coeffs[j * 11:(j + 1) * 11])
                        j += 1

        # We cannot trust this self.homos until we come to the phrase:
        #   SYMMETRIES FOR INITAL GUESS ORBITALS FOLLOW
        # which either is followed by "ALPHA" or "BOTH" at which point we can say
        # for certain that it is an un/restricted calculations.
        # Note that MCSCF calcs also print this search string, so make sure
        #   that self.homos does not exist yet.
        if line[1:28] == "NUMBER OF OCCUPIED ORBITALS" and not hasattr(self, 'homos'):

            homos = [int(line.split()[-1])-1]
            line = next(inputfile)
            homos.append(int(line.split()[-1])-1)

            self.set_attribute('homos', homos)

        if line.find("SYMMETRIES FOR INITIAL GUESS ORBITALS FOLLOW") >= 0:
            # Not unrestricted, so lop off the second index.
            # In case the search string above was not used (ex. FMO in exam38),
            #   we can try to use the next line which should also contain the
            #   number of occupied orbitals.
            if line.find("BOTH SET(S)") >= 0:
                nextline = next(inputfile)
                if "ORBITALS ARE OCCUPIED" in nextline:
                    homos = int(nextline.split()[0])-1
                    if hasattr(self, "homos"):
                        try:
                            assert self.homos[0] == homos
                        except AssertionError:
                            self.logger.warning("Number of occupied orbitals not consistent. This is normal for ECP and FMO jobs.")
                    else:
                        self.homos = [homos]
                self.homos = numpy.resize(self.homos, [1])

        # Set the total number of atoms, only once.
        # Normally GAMESS print TOTAL NUMBER OF ATOMS, however in some cases
        #   this is slightly different (ex. lower case for FMO in exam37).
        if not hasattr(self, "natom") and "NUMBER OF ATOMS" in line.upper():
            natom = int(line.split()[-1])
            self.set_attribute('natom', natom)

        # The first is from Julien's Example and the second is from Alexander's
        # I think it happens if you use a polar basis function instead of a cartesian one
        if line.find("NUMBER OF CARTESIAN GAUSSIAN BASIS") == 1 or line.find("TOTAL NUMBER OF BASIS FUNCTIONS") == 1:
            nbasis = int(line.strip().split()[-1])
            self.set_attribute('nbasis', nbasis)

        elif line.find("TOTAL NUMBER OF CONTAMINANTS DROPPED") >= 0:
            nmos_dropped = int(line.split()[-1])
            if hasattr(self, "nmo"):
                self.set_attribute('nmo', self.nmo - nmos_dropped)
            else:
                self.set_attribute('nmo', self.nbasis - nmos_dropped)

        # Note that this line is present if ISPHER=1, e.g. for C_bigbasis
        elif line.find("SPHERICAL HARMONICS KEPT IN THE VARIATION SPACE") >= 0:
            nmo = int(line.strip().split()[-1])
            self.set_attribute('nmo', nmo)

        # Note that this line is not always present, so by default
        # NBsUse is set equal to NBasis (see below).
        elif line.find("TOTAL NUMBER OF MOS IN VARIATION SPACE") == 1:
            nmo = int(line.split()[-1])
            self.set_attribute('nmo', nmo)

        elif line.find("OVERLAP MATRIX") == 0 or line.find("OVERLAP MATRIX") == 1:
            # The first is for PC-GAMESS, the second for GAMESS
            # Read 1-electron overlap matrix
            if not hasattr(self, "aooverlaps"):
                self.aooverlaps = numpy.zeros((self.nbasis, self.nbasis), "d")
            else:
                self.logger.info("Reading additional aooverlaps...")
            base = 0
            while base < self.nbasis:
                self.updateprogress(inputfile, "Overlap")

                self.skip_lines(inputfile, ['b', 'basis_fn_number', 'b'])

                for i in range(self.nbasis - base):  # Fewer lines each time
                    line = next(inputfile)
                    temp = line.split()
                    for j in range(4, len(temp)):
                        self.aooverlaps[base+j-4, i+base] = float(temp[j])
                        self.aooverlaps[i+base, base+j-4] = float(temp[j])
                base += 5

        # ECP Pseudopotential information
        if "ECP POTENTIALS" in line:
            if not hasattr(self, "coreelectrons"):
                self.coreelectrons = [0]*self.natom

            self.skip_lines(inputfile, ['d', 'b'])

            header = next(inputfile)
            while header.split()[0] == "PARAMETERS":
                name = header[17:25]
                atomnum = int(header[34:40])
                # The pseudopotnetial is given explicitely
                if header[40:50] == "WITH ZCORE":
                    zcore = int(header[50:55])
                    lmax = int(header[63:67])
                    self.coreelectrons[atomnum-1] = zcore
                # The pseudopotnetial is copied from another atom
                if header[40:55] == "ARE THE SAME AS":
                    atomcopy = int(header[60:])
                    self.coreelectrons[atomnum-1] = self.coreelectrons[atomcopy-1]
                line = next(inputfile)
                while line.split() != []:
                    line = next(inputfile)
                header = next(inputfile)

        # This was used before refactoring the parser, geotargets was set here after parsing.
        #if not hasattr(self, "geotargets"):
        #    opttol = 1e-4
        #    self.geotargets = numpy.array([opttol, 3. / opttol], "d")
        #if hasattr(self,"geovalues"): self.geovalues = numpy.array(self.geovalues, "d")

        # This is quite simple to parse, but some files seem to print certain lines twice,
        # repeating the populations without charges, but not in proper order.
        # The unrestricted calculations are a bit tricky, since GAMESS-US prints populations
        # for both alpha and beta orbitals in the same format and with the same title,
        # but it still prints the charges only at the very end.
        if "TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS" in line:

            if not hasattr(self, "atomcharges"):
                self.atomcharges = {}

            header = next(inputfile)
            line = next(inputfile)

            # It seems that when population are printed twice (without charges),
            # there is a blank line along the way (after the first header),
            # so let's get a flag out of that circumstance.
            doubles_printed = line.strip() == ""
            if doubles_printed:
                title = next(inputfile)
                header = next(inputfile)
                line = next(inputfile)

            # Only go further if the header had five columns, which should
            # be the case when both populations and charges are printed.
            # This is pertinent for both double printing and unrestricted output.
            if not len(header.split()) == 5:
                return
            mulliken, lowdin = [], []
            while line.strip():
                if line.strip() and doubles_printed:
                    line = next(inputfile)
                mulliken.append(float(line.split()[3]))
                lowdin.append(float(line.split()[5]))
                line = next(inputfile)
            self.atomcharges["mulliken"] = mulliken
            self.atomcharges["lowdin"] = lowdin

        #          ---------------------
        #          ELECTROSTATIC MOMENTS
        #          ---------------------
        #
        # POINT   1           X           Y           Z (BOHR)    CHARGE
        #                -0.000000    0.000000    0.000000       -0.00 (A.U.)
        #         DX          DY          DZ         /D/  (DEBYE)
        #     0.000000   -0.000000    0.000000    0.000000
        #
        if line.strip() == "ELECTROSTATIC MOMENTS":

            self.skip_lines(inputfile, ['d', 'b'])
            line = next(inputfile)

            # The old PC-GAMESS prints memory assignment information here.
            if "MEMORY ASSIGNMENT" in line:
                memory_assignment = next(inputfile)
                line = next(inputfile)

            # If something else ever comes up, we should get a signal from this assert.
            assert line.split()[0] == "POINT"

            # We can get the reference point from here, as well as
            # check here that the net charge of the molecule is correct.
            coords_and_charge = next(inputfile)
            assert coords_and_charge.split()[-1] == '(A.U.)'
            reference = numpy.array([float(x) for x in coords_and_charge.split()[:3]])
            reference = utils.convertor(reference, 'bohr', 'Angstrom')
            charge = int(round(float(coords_and_charge.split()[-2])))
            self.set_attribute('charge', charge)

            dipoleheader = next(inputfile)
            assert dipoleheader.split()[:3] == ['DX', 'DY', 'DZ']
            assert dipoleheader.split()[-1] == "(DEBYE)"

            dipoleline = next(inputfile)
            dipole = [float(d) for d in dipoleline.split()[:3]]

            # The dipole is always the first multipole moment to be printed,
            # so if it already exists, we will overwrite all moments since we want
            # to leave just the last printed value (could change in the future).
            if not hasattr(self, 'moments'):
                self.moments = [reference, dipole]
            else:
                try:
                    assert self.moments[1] == dipole
                except AssertionError:
                    self.logger.warning('Overwriting previous multipole moments with new values')
                    self.logger.warning('This could be from post-HF properties or geometry optimization')
                    self.moments = [reference, dipole]

        # Static polarizability from a harmonic frequency calculation
        # with $CPHF/POLAR=.TRUE.
        if line.strip() == 'ALPHA POLARIZABILITY TENSOR (ANGSTROMS**3)':
            if not hasattr(self, 'polarizabilities'):
                self.polarizabilities = []
            polarizability = numpy.zeros(shape=(3, 3))
            self.skip_lines(inputfile, ['d', 'b', 'directions'])
            for i in range(3):
                line = next(inputfile)
                polarizability[i, :i+1] = [float(x) for x in line.split()[1:]]
            polarizability = utils.symmetrize(polarizability, use_triangle='lower')
            # Convert from Angstrom**3 to bohr**3 (a.u.**3).
            volume_convert = numpy.vectorize(lambda x: x * utils.convertor(1, 'Angstrom', 'bohr') ** 3)
            polarizability = volume_convert(polarizability)
            self.polarizabilities.append(polarizability)

        # Static and dynamic polarizability from RUNTYP=TDHF.
        if line.strip() == 'TIME-DEPENDENT HARTREE-FOCK NLO PROPERTIES':
            if not hasattr(self, 'polarizabilities'):
                self.polarizabilities = []
            polarizability = numpy.empty(shape=(3, 3))
            coord_to_idx = {'X': 0, 'Y': 1, 'Z': 2}
            self.skip_lines(inputfile, ['d', 'b', 'dots'])
            line = next(inputfile)
            assert 'ALPHA AT' in line
            self.skip_lines(inputfile, ['dots', 'b'])
            for a in range(3):
                for b in range(3):
                    line = next(inputfile)
                    tokens = line.split()
                    i, j = coord_to_idx[tokens[1][0]], coord_to_idx[tokens[1][1]]
                    polarizability[i, j] = tokens[3]
            self.polarizabilities.append(polarizability)

        if line[:30] == ' ddikick.x: exited gracefully.'\
                or line[:41] == ' EXECUTION OF FIREFLY TERMINATED NORMALLY'\
                or line[:40] == ' EXECUTION OF GAMESS TERMINATED NORMALLY':
            self.metadata['success'] = True