Problem correctly parsing moenergies
As reported to me by Martin R., a GAMESS file does not correctly parse `moenergies`. After briefly looking at his file, the problem is due to LZ value analysis being printed between the alpha and beta orbitals.

diff --git a/.travis.yml b/.travis.yml
index 218f6d52..bd0a7e00 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -13,4 +13,4 @@ install:
 script:
   - cd test && python testall.py --status --terse
   - cd ../data && bash regression_download.sh
-  - cd ../test && python regression.py
+  - cd ../test && python regression.py --status --traceback
diff --git a/data/regressionfiles.txt b/data/regressionfiles.txt
index 08eaa160..906f4bff 100644
--- a/data/regressionfiles.txt
+++ b/data/regressionfiles.txt
@@ -66,6 +66,7 @@ GAMESS/GAMESS-US2010/dvb_td.log
 GAMESS/GAMESS-US2010/n_water.log
 GAMESS/GAMESS-US2012/dvb_gopt_a_unconverged.out
 GAMESS/GAMESS-US2012/ScriptBicarbNa1B3LYP-B2Vib-good.gamout
+GAMESS/GAMESS-US2013/N-UHF.out
 GAMESS/PCGAMESS/C_bigbasis.out
 GAMESS/PCGAMESS/CN_ir.out
 GAMESS/PCGAMESS/aomixgamess.log
diff --git a/src/cclib/parser/gamessparser.py b/src/cclib/parser/gamessparser.py
index 3ee9c169..f4967ab0 100644
--- a/src/cclib/parser/gamessparser.py
+++ b/src/cclib/parser/gamessparser.py
@@ -153,6 +153,7 @@ class GAMESS(logfileparser.Logfile):
                 if line[8:23] == "CCSD(T) ENERGY:":
                     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"):
@@ -171,7 +172,7 @@ class GAMESS(logfileparser.Logfile):
             mult = int(line.split()[-1])
             self.set_attribute('mult', mult)
 
-        # etenergies (used only for CIS runs now)
+        # etenergies (originally used only for CIS runs, but now also TD-DFT)
         if "EXCITATION ENERGIES" in line and line.find("DONE WITH") < 0:
 
             if not hasattr(self, "etenergies"):
@@ -190,6 +191,7 @@ class GAMESS(logfileparser.Logfile):
             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])
@@ -437,12 +439,13 @@ class GAMESS(logfileparser.Logfile):
             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 == []):
-            # 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
+
             self.updateprogress(inputfile, "Coordinates")
 
             if self.firststdorient:
@@ -496,7 +499,7 @@ class GAMESS(logfileparser.Logfile):
                 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.
+                # 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
@@ -517,14 +520,14 @@ class GAMESS(logfileparser.Logfile):
             if not hasattr(self,"scfvalues"):
                 self.scfvalues = []
 
-            line = next(inputfile)
-
-            # Normally the iteration print in 6 columns.
+            # Normally the iterations print in 6 columns.
             # For ROHF, however, it is 5 columns, thus this extra parameter.
             if "ROHF" in self.scftype:
-                valcol = 4
+                self.scf_valcol = 4
             else:
-                valcol = 5
+                self.scf_valcol = 5
+
+            line = next(inputfile)
 
             # SCF iterations are terminated by a blank line.
             # The first four characters usually contains the step number.
@@ -539,7 +542,28 @@ class GAMESS(logfileparser.Logfile):
                 except ValueError:
                     pass
                 else:
-                    values.append([float(line.split()[valcol])])
+                    values.append([float(line.split()[self.scf_valcol])])
+                line = next(inputfile)
+            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)
 
@@ -773,53 +797,56 @@ class GAMESS(logfileparser.Logfile):
                 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 OBRITALS") == 10:
-            # The details returned come from the *final* report of evalues and
-            #   the last list of symmetries in the log file.
-            # Should be followed by lines like this:
-            #           ------------
-            #           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
-            # and so forth... with blanks lines between blocks of 5 orbitals each.
-            # Warning! There are subtle differences between GAMESS-US and PC-GAMES
-            #   in the formatting of the first four columns.
-            #
-            # Watch out 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
-            #
-            # This is fine for GeoOpt and SP, but may be weird for TD and Freq.
 
             # 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 = []
@@ -836,7 +863,9 @@ class GAMESS(logfileparser.Logfile):
                 self.updateprogress(inputfile, "Coefficients")
 
                 line = next(inputfile)
-                # Make sure that this section does not end prematurely - checked by regression test 2CO.ccsd.aug-cc-pVDZ.out.
+
+                # This makes sure that this section does not end prematurely,
+                # which happens in regression 2CO.ccsd.aug-cc-pVDZ.out.
                 if line.strip() != "":
                     break;
                 
@@ -858,9 +887,7 @@ class GAMESS(logfileparser.Logfile):
                 if line.strip():
                     self.mosyms[0].extend(list(map(self.normalisesym, line.split())))
                 
-                # Now we have nbasis lines.
-                # Going to use the same method as for normalise_aonames()
-                # to extract basis set information.
+                # 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
@@ -869,14 +896,16 @@ class GAMESS(logfileparser.Logfile):
                 for i in range(self.nbasis):
                     line = next(inputfile)
 
-                    # If line is empty, break (ex. for FMO in exam37).
+                    # 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
@@ -901,21 +930,27 @@ class GAMESS(logfileparser.Logfile):
                             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
 
             line = next(inputfile)
-            # If it's restricted and no more properties:
+
+            # If it's a restricted calc and no more properties, we have:
+            #
             #  ...... END OF RHF/DFT CALCULATION ......
-            # If there are more properties (DENSITY MATRIX):
+            #
+            # If there are more properties (such as the density matrix):
             #               --------------
             #
-            # If it's unrestricted we have:
+            # If it's an unrestricted calculation, however, we now get the beta orbitals:
             #
             #  ----- BETA SET ----- 
             #
@@ -924,8 +959,16 @@ class GAMESS(logfileparser.Logfile):
             #          ------------
             #
             #                      1          2          3          4          5
-            # ... and so forth.
+            # ...
+            #
             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)
+
             if line[2:22] == "----- BETA SET -----":
                 self.mocoeffs.append(numpy.zeros((self.nmo, self.nbasis), "d"))
                 self.moenergies.append([])
diff --git a/src/cclib/parser/gaussianparser.py b/src/cclib/parser/gaussianparser.py
index f9f3cf1f..be1a0f73 100644
--- a/src/cclib/parser/gaussianparser.py
+++ b/src/cclib/parser/gaussianparser.py
@@ -109,7 +109,87 @@ class Gaussian(logfileparser.Logfile):
     def extract(self, inputfile, line):
         """Extract information from the file object inputfile."""
 
-        # Number of atoms.
+        # This block contains some general information as well as coordinates,
+        # which could be parsed in the future:
+        #
+        # Symbolic Z-matrix:
+        # Charge =  0 Multiplicity = 1
+        # C                     0.73465   0.        0. 
+        # C                     1.93465   0.        0. 
+        # C 
+        # ...
+        #
+        # It also lists fragments, if there are any, which is potentially valuable:
+        #
+        # Symbolic Z-matrix:
+        # Charge =  0 Multiplicity = 1 in supermolecule
+        # Charge =  0 Multiplicity = 1 in fragment      1.
+        # Charge =  0 Multiplicity = 1 in fragment      2.
+        # B(Fragment=1)         0.06457  -0.0279    0.01364 
+        # H(Fragment=1)         0.03117  -0.02317   1.21604 
+        # ...
+        #
+        # Note, however, that currently we only parse information for the whole system
+        # or supermolecule as Gaussian calls it.
+        if line.strip() == "Symbolic Z-matrix:":
+
+            self.updateprogress(inputfile, "Symbolic Z-matrix", self.fupdate)
+
+            line = inputfile.next()
+            while line.split()[0] == 'Charge':
+
+                # For the supermolecule, we can parse the charge and multicplicity.
+                regex = ".*=(.*)Mul.*=\s*-?(\d+).*"
+                match = re.match(regex, line)
+                assert match, "Something unusual about the line: '%s'" % line
+                
+                self.set_attribute('charge', int(match.groups()[0]))
+                self.set_attribute('mult', int(match.groups()[1]))
+
+                if line.split()[-2] == "fragment":
+                    self.nfragments = int(line.split()[-1].strip('.'))
+
+                if line.strip()[-13:] == "model system.":
+                    self.nmodels = getattr(self, 'nmodels', 0) + 1
+
+                line = inputfile.next()
+
+            # The remaining part will allow us to get the atom count.
+            # When coordinates are given, there is a blank line at the end, but if
+            # there is a Z-matrix here, there will also be variables and we need to
+            # stop at those to get the right atom count.
+            # Also, in older versions there is bo blank line (G98 regressions),
+            # so we need to watch out for leaving the link.
+            natom = 0
+            while line.split() and not "Variables" in line and not "Leave Link" in line:
+                natom += 1
+                line = inputfile.next()
+            self.set_attribute('natom', natom)
+
+        # Continuing from above, there is not always a symbolic matrix, for example
+        # if the Z-matrix was in the input file. In such cases, try to match the
+        # line and get at the charge and multiplicity.
+        # 
+        #   Charge =  0 Multiplicity = 1 in supermolecule
+        #   Charge =  0 Multiplicity = 1 in fragment  1.
+        #   Charge =  0 Multiplicity = 1 in fragment  2.
+        if line[1:7] == 'Charge' and line.find("Multiplicity") >= 0:
+
+            self.updateprogress(inputfile, "Charge and Multiplicity", self.fupdate)
+
+            if line.split()[-1] == "supermolecule" or not "fragment" in line:
+
+                regex = ".*=(.*)Mul.*=\s*-?(\d+).*"
+                match = re.match(regex, line)
+                assert match, "Something unusual about the line: '%s'" % line
+                
+                self.set_attribute('charge', int(match.groups()[0]))
+                self.set_attribute('mult', int(match.groups()[1]))
+
+            if line.split()[-2] == "fragment":
+                self.nfragments = int(line.split()[-1].strip('.'))
+
+        # Number of atoms is also explicitely printed after the above.
         if line[1:8] == "NAtoms=":
 
             self.updateprogress(inputfile, "Attributes", self.fupdate)
@@ -213,6 +293,21 @@ class Gaussian(logfileparser.Logfile):
             self.set_attribute('natom', len(atomnos))
             self.set_attribute('atomnos', atomnos)
 
+        # This is a bit of a hack for regression Gaussian09/BH3_fragment_guess.pop_minimal.log
+        # to skip output for all fragments, assuming the supermolecule is always printed first.
+        # Eventually we want to make this more general, or even better parse the output for
+        # all fragment, but that will happen in a newer version of cclib.
+        if line[1:16] == "Fragment guess:" and getattr(self, 'nfragments', 0) > 1:
+            if not "full" in line:
+                inputfile.file.seek(0, 2)
+
+        # Another hack for regression Gaussian03/ortho_prod_prod_freq.log, which is an ONIOM job.
+        # Basically for now we stop parsing after the output for the real system, because
+        # currently we don't support changes in system size or fragments in cclib. When we do,
+        # we will want to parse the model systems, too, and that is what nmodels could track.
+        if "ONIOM: generating point" in line and line.strip()[-13:] == 'model system.' and getattr(self, 'nmodels', 0) > 0:
+            inputfile.file.seek(0,2)
+
         # With the gfinput keyword, the atomic basis set functios are:
         #
         # AO basis set in the form of general basis input (Overlap normalization):
@@ -591,21 +686,6 @@ class Gaussian(logfileparser.Logfile):
             if not hasattr(self, "scannames"):
                 self.scannames = colmnames.split()[1:-1]
 
-        # Charge and multiplicity.
-        # If counterpoise correction is used, multiple lines match.
-        # The first one contains charge/multiplicity of the whole molecule.:
-        #   Charge =  0 Multiplicity = 1 in supermolecule
-        #   Charge =  0 Multiplicity = 1 in fragment  1.
-        #   Charge =  0 Multiplicity = 1 in fragment  2.
-        if line[1:7] == 'Charge' and line.find("Multiplicity")>=0:
-
-            regex = ".*=(.*)Mul.*=\s*-?(\d+).*"
-            match = re.match(regex, line)
-            assert match, "Something unusual about the line: '%s'" % line
-            
-            self.set_attribute('charge', int(match.groups()[0]))
-            self.set_attribute('mult', int(match.groups()[1]))
-
         # Orbital symmetries.
         if line[1:20] == 'Orbital symmetries:' and not hasattr(self, "mosyms"):
 
diff --git a/src/cclib/parser/psiparser.py b/src/cclib/parser/psiparser.py
index 05f2034f..263041e2 100644
--- a/src/cclib/parser/psiparser.py
+++ b/src/cclib/parser/psiparser.py
@@ -637,8 +637,8 @@ class Psi(logfileparser.Logfile):
                 self.ccenergies = []
             self.ccenergies.append(ccsd_energy)
 
-        # The geometry convergence targets and values are printed in a table, but the legends
-        # described do not seems to be working. Still, probably exact slicing of the line needs
+        # The geometry convergence targets and values are printed in a table, with the legends
+        # describing the convergence annotation. Probably exact slicing of the line needs
         # to be done in order to extract the numbers correctly.
         #
         #  ==> Convergence Check <==
@@ -663,20 +663,25 @@ class Psi(logfileparser.Logfile):
             criteria = next(inputfile)
             geotargets = []
             for istart in starts:
-                if line[istart:istart+9].strip():
-                    geotargets.append(float(line[istart:istart+9]))
+                if criteria[istart:istart+9].strip():
+                    geotargets.append(float(criteria[istart:istart+9]))
+                else:
+                    geotargets.append(numpy.nan)
 
             self.skip_line(inputfile, 'dashes')
 
             values = next(inputfile)
             geovalues = []
             for istart in starts:
-                if line[istart:istart+9].strip():
-                    geovalues.append(float(line[istart:istart+9]))
+                if values[istart:istart+9].strip():
+                    geovalues.append(float(values[istart:istart+9]))
 
+            # This assertion may be too restrictive, but we haven't seen the geotargets change.
+            # If such an example comes up, update the value since we're interested in the last ones.
             if not hasattr(self, 'geotargets'):
-                self.geotargets = []
-            self.geotargets.append(geotargets)
+                self.geotargets = geotargets
+            else:
+                assert self.geotargets == geotargets
 
             if not hasattr(self, 'geovalues'):
                 self.geovalues = []
diff --git a/test/regression.py b/test/regression.py
index 34787d59..7a900c69 100644
--- a/test/regression.py
+++ b/test/regression.py
@@ -75,6 +75,10 @@ def testGAMESS_GAMESS_US2012_dvb_gopt_a_unconverged_out(logfile):
     """An unconverged geometry optimization to test for empty optdone (see #103 for details)."""
     assert hasattr(logfile.data, 'optdone') and logfile.data.optdone == []
 
+def testGAMESS_GAMESS_US2013_N_UHF_out(logfile):
+    """An UHF job that has an LZ value analysis between the alpha and beta orbitals."""
+    assert len(logfile.data.moenergies) == 2
+
 def testGAMESS_WinGAMESS_dvb_td_trplet_2007_03_24_r1_out(logfile):
     """Do some basic checks for this old unit test that was failing.
 
@@ -729,17 +733,19 @@ def main(which=[], traceback=False, status=False):
         print("\n".join(orphaned_tests))
 
     if status and errors > 0:
-        os.exit(1)
+        sys.exit(1)
 
 if __name__=="__main__":
 
     # If 'test' is passed as the first argument, do a doctest on this module.
     # Otherwise, any arguments are used to limit the test to the packages/parsers
     # passed as arguments. No arguments implies all parsers.
+    # In general, it would be best to replace this trickery by argparse magic.
     if len(sys.argv) == 2 and sys.argv[1] == "test":
         import doctest
         doctest.testmod()
     else:
-        traceback = "traceback" in sys.argv or "tb" in sys.argv
-        status = "status" in sys.argv or "--status" in sys.argv
-        main(sys.argv[1:], traceback, status)
+        traceback = "--traceback" in sys.argv
+        status = "--status" in sys.argv
+        which = [arg for arg in sys.argv[1:] if not arg in ["--status", "--traceback"]]
+        main(which, traceback, status)
diff --git a/test/testSP.py b/test/testSP.py
index c13c4a69..e7dfe74e 100644
--- a/test/testSP.py
+++ b/test/testSP.py
@@ -143,36 +143,46 @@ class GenericSPTest(bettertest.TestCase):
         # is at the origin, so we now what to expect.
         reference = self.data.moments[0]
         self.assertEquals(len(reference), 3)
-        self.assertInside(sum(reference), 0.0, 0.0001)
+        for x in reference:
+            self.assertInside(x, 0.0, 0.001)
 
         # Length and value of dipole moment should always be correct (zero for this test).
         dipole = self.data.moments[1]
         self.assertEquals(len(dipole), 3)
-        self.assertInside(sum(dipole), 0.0, 0.0001)
+        for d in dipole:
+            self.assertInside(d, 0.0, 0.001)
 
-        # If the quadrupole is there, we can expect roughly -50B for the XX moment.
+        # If the quadrupole is there, we can expect roughly -50B for the XX moment,
+        # -50B for the YY moment and and -60B for the ZZ moment.
         if len(self.data.moments) > 2:
             quadrupole = self.data.moments[2]
             self.assertEquals(len(quadrupole), 6)
             self.assertInside(quadrupole[0], -50, 5)
+            self.assertInside(quadrupole[3], -50, 5)
+            self.assertInside(quadrupole[5], -60, 5)
 
         # If the octupole is there, it should have 10 components and be zero.
         if len(self.data.moments) > 3:
             octupole = self.data.moments[3]
             self.assertEquals(len(octupole), 10)
-            self.assertInside(sum(octupole), 0.0, 0.0001)
+            for m in octupole:
+                self.assertInside(m, 0.0, 0.001)
 
-        # Hexadecapole should have 15 elements and an XX componenet of around -1900 Debye*ang^2.
+        # Hexadecapole should have 15 elements and an XXXX componenet of around -1900 Debye*ang^2,
+        # YYYY component of -330B and ZZZZ component of -50B.
         if len(self.data.moments) > 4:
             hexadecapole = self.data.moments[4]
             self.assertEquals(len(hexadecapole), 15)
             self.assertInside(hexadecapole[0], -1900, 100)
+            self.assertInside(hexadecapole[10], -330, 10)
+            self.assertInside(hexadecapole[14], -50, 5)
 
         # The are 21 unique 32-pole moments, and all are zero in this test case.
         if len(self.data.moments) > 5:
             moment32 = self.data.moments[5]
             self.assertEquals(len(moment32), 21)
-            self.assertInside(sum(moment32), 0.0, 0.0001)
+            for m in moment32:
+                self.assertInside(m, 0.0, 0.001)
 
 
 
@@ -315,7 +325,4 @@ class QChemSPTest(GenericSPTest):
 if __name__=="__main__":
 
     from testall import testall
-    parser = None
-    if len(sys.argv)==2:
-        parser = sys.argv[1]
-    testall(parser, modules=["SP"])
+    testall(modules=["SP"])
diff --git a/test/testTD.py b/test/testTD.py
index 2737de01..f3b513ae 100644
--- a/test/testTD.py
+++ b/test/testTD.py
@@ -21,7 +21,11 @@ class GenericTDTest(bettertest.TestCase):
 
     def testenergies(self):
         """Is the l_max reasonable?"""
+
         self.assertEqual(len(self.data.etenergies), self.number)
+
+        # Note that if all oscillator strengths are zero (like for triplets)
+        # then this will simply pick out the first energy.
         idx_lambdamax = [i for i, x in enumerate(self.data.etoscs)
                          if x==max(self.data.etoscs)][0]
         self.assertInside(self.data.etenergies[idx_lambdamax], self.expected_l_max, 5000)
@@ -95,17 +99,11 @@ class OrcaTDDFTTest(GenericTDTest):
         self.assertInside(max(self.data.etoscs), 1.0, 0.1)
 
 
-class GenericTDTesttrp(GenericTDTest):
+class GenericTDDFTtrpTest(GenericTDTest):
     """Generic time-dependent HF/DFT (triplet) unittest"""
 
     number = 5
-
-    def testenergies(self):
-        """Is the l_max reasonable?"""
-        self.assertEqual(len(self.data.etenergies), self.number)
-        idx_lambdamax = [i for i, x in enumerate(self.data.etoscs)
-                         if x==max(self.data.etoscs)][0]
-        self.assertInside(self.data.etenergies[idx_lambdamax], 24500, 100)
+    expected_l_max = 24500
     
     def testoscs(self):
         """Triplet excitations should be disallowed."""
@@ -113,16 +111,6 @@ class GenericTDTesttrp(GenericTDTest):
         self.assertInside(max(self.data.etoscs), 0.0, 0.01)
 
 
-class GAMESSUSTDDFTtrpTest(GenericTDTesttrp):
-    """Customized TD DFT (restricted) triplet unittest"""
-
-    #old_tests = ["GAMESS/WinGAMESS/dvb_td_trplet_2007.03.24.r1.out.gz"]
-
-    def testsymsnumber(self):
-        """Is the length of etsyms correct? PASS"""
-        pass
-
-
 if __name__=="__main__":
 
     from testall import testall
diff --git a/test/testdata b/test/testdata
index 141d96de..9ae4ee91 100644
--- a/test/testdata
+++ b/test/testdata
@@ -165,9 +165,9 @@ SPun      QChem       QChemSPunTest         basicQChem4.2         dvb_sp_un.out
 TD        ADF         ADFTDDFTTest          basicADF2007.01       dvb_td.adfout
 TD        ADF         ADFTDDFTTest          basicADF2013.01       dvb_td.adfout
 TD        GAMESS      GenericTDTest         basicFirefly8.0       dvb_td.out
-TD        GAMESS      GenericTDTest         basicFirefly8.0       dvb_td_trplet.out
+TD        GAMESS      GenericTDDFTtrpTest   basicFirefly8.0       dvb_td_trplet.out
 TD        GAMESS      GAMESSUSTDDFTTest     basicGAMESS-US2012    dvb_td.out
-TD        GAMESS      GAMESSUSTDDFTtrpTest  basicGAMESS-US2012    dvb_td_trplet.out
+TD        GAMESS      GenericTDDFTtrpTest   basicGAMESS-US2012    dvb_td_trplet.out
 TD        Gaussian    GaussianTDDFTTest     basicGaussian03       dvb_td.out
 TD        Gaussian    GaussianTDDFTTest     basicGaussian09       dvb_td.out
 TD        Jaguar      GenericTDTest         basicJaguar8.3        dvb_td.out
['test/testSP.py', 'data/regressionfiles.txt', 'test/regression.py', 'test/testTD.py', 'src/cclib/parser/gaussianparser.py', 'src/cclib/parser/psiparser.py', 'src/cclib/parser/gamessparser.py', '.travis.yml', 'test/testdata']