MOPAC: read orbitals from the GRAPHF file, and teach tests/e2e File= types

MOPAC's printed EIGENVECTORS block is the weaker source for orbitals.
Its GRAPHF file (<job>.mgf) carries the same data better:

  precision   Fortran doubles, against four printed decimals
  occupancy   one field per orbital, so mopac.orbocc no longer has to
              reconstruct it from atomic orbital populations plus the
              open-shell count passed in argv
  shape       free-form, not the fixed-width column groups that merged
              two coefficients into one for ORCA (#108)

ai.mopac now asks for GRAPHF alongside VECTORS -- they serve different
readers, VECTORS the person reading the output and GRAPHF ECCE -- and a
new File= parse type reads it, the same mechanism Gaussian uses for
fort.7 and GAMESS-UK for ftn058.

Both entries emit the same keys on purpose.  The file is delivered
after the matched blocks, so it wins where it exists, and the printed
block stays as the fallback for decks generated before GRAPHF was
requested.  That matters: a calculation run an hour ago would otherwise
lose its orbitals.

The .mgf also records the Slater exponents of the valence basis, which
is what a 3-D orbital renderer for a semiempirical code would need.
Nothing reads them yet, so they are deliberately not emitted rather
than shipped unused.

tests/e2e previously skipped File= parse types entirely, silently --
the codes' most important data for three of them.  It now decodes the
jmFILE announcements and runs the scripts as the client does.  Two
traps in that path, both now handled and documented:

  importDir   MsgSendFile resolves "$idir$file", concatenated with no
              separator, and idir defaults to "-".  Unset, every File=
              type looks for "-fort.7" and is skipped with no error.
  deck name   MOPAC names its outputs after its input, so a fixture
              called ch4_mos.mop writes ch4_mos.mgf while mopac.desc
              names mopac.mgf -- correct, since ECCE always writes
              mopac.mop.  The fixtures are renamed to match, and the
              case now asserts the GRAPHF file really was the source by
              the precision of the energies: eight decimals from the
              file, six from the printed block.

That assertion was checked both ways -- pointing File= at a name that
does not exist makes it fail with "6 decimals; 6 means the fallback
ran", and restoring it passes.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFMUNbw2UpRPmgsqviAowb

Author: andy

eccejobmonitor.propbuf

@@ -325,8 +325,16 @@ sub KeywordLine {
   #  all, and treating absent as "0" silently gave those calculations no
   #  orbitals -- which is every calculation anyone already had.  Only an
   #  explicit "0" turns it off.
+  #  GRAPHF as well as VECTORS.  They serve different readers: VECTORS
+  #  puts the eigenvectors in the output file where a person can read
+  #  them, GRAPHF writes <job>.mgf, which ECCE parses instead because it
+  #  carries full-precision coefficients, the occupancies and the
+  #  symmetry labels in one machine-readable place -- see mopac.mgf.
+  #  GRAPHF also records the Slater exponents, which is what a 3-D
+  #  orbital renderer for a semiempirical code would need.
   my $printMOs = $AbiDict{"ES.Theory.SCF.PrintMOs"};
-  push(@kw, "VECTORS", "ALLVEC") if (!defined($printMOs) || $printMOs ne "0");
+  push(@kw, "VECTORS", "ALLVEC", "GRAPHF")
+    if (!defined($printMOs) || $printMOs ne "0");
 
   # PRECISE: "TIGHTER CRITERIA TO BE USED" (binary's own banner text).
   push(@kw, "PRECISE") if ($AbiDict{"ES.Theory.SCF.Precise"} eq "1");

scripts/parsers/ai.mopac

@@ -187,6 +187,27 @@ Frequency=last
 End=HEATS OF FORMATION ARE RELATIVE
 [END]
 
+#  MO data from MOPAC's GRAPHF file rather than its printed output.
+#
+#  A File= parse type, the same mechanism Gaussian uses for fort.7 and
+#  GAMESS-UK for ftn058: eccejobmonitor ships the named file whole and
+#  the script parses it, instead of matching Begin/End in the job's
+#  stdout.
+#
+#  This entry and the printed-eigenvector one above emit the SAME keys,
+#  deliberately.  The .mgf is the better source -- full-precision
+#  coefficients, occupancies and symmetry labels all in one place -- but
+#  it only exists when the deck asked for GRAPHF, which decks generated
+#  before that did not.  Both are kept so an older calculation still
+#  gets its orbitals from the printed block, and a new one gets the
+#  better data.  Whichever is delivered later wins, which is what the
+#  test fixture pins.
+[MOFILE]
+Script=mopac.mgf
+File=mopac.mgf
+Begin=mopac.mgf
+[END]
+
 [ORBOCC][ORBOCCBETA]
 Script=mopac.orbocc
 Begin=ATOMIC ORBITAL ELECTRON POPULATIONS

scripts/parsers/mopac.desc

@@ -0,0 +1,178 @@
+#!/usr/bin/env perl
+################################################################################
+#
+# Purpose:
+# Parses MOPAC's GRAPHF file (<job>.mgf) into MO/ORBENG/ORBOCC/ORBSYM and
+# the UHF beta counterparts.
+#
+# WHY THIS RATHER THAN THE PRINTED EIGENVECTORS.  mopac.mo reads the
+# human-readable "EIGENVECTORS" block out of the main output, which works
+# but is the weaker source in three ways:
+#
+#   precision   the printed block carries four decimals; this file carries
+#               Fortran D-notation doubles.
+#   occupancy   the printed block has none, so mopac.orbocc has to
+#               reconstruct it from the atomic orbital populations plus
+#               the open-shell count ECCE passes in argv.  Here it is
+#               simply the second field of each ORBITAL line.
+#   shape       the printed block is fixed-width and column-grouped, the
+#               exact shape that merged two coefficients into one for
+#               ORCA (issue #108).  This file is free-form.
+#
+# Real MOPAC 23.1.2 shape:
+#
+#       3 MOPAC-Graphical data Version 23.1.2
+#      8    0.0000000  -0.0000030   0.1149079  -0.6491      <- Z x y z charge
+#      1    0.0000000   0.7594689  -0.4644394   0.3246
+#      1   -0.0000000  -0.7594570  -0.4644300   0.3246
+#     5.9723090  2.3490170  0.0000000                       <- zeta s p d
+#     1.2602370  0.0000000  0.0000000
+#     1.2602370  0.0000000  0.0000000
+#    ORBITAL 2 1a1    -29.6190                              <- occ sym E(eV)
+#    0.85155170D+00-0.28614321D-16 ...                      <- coefficients
+#
+# The atom count on line 1 gives the size of BOTH leading blocks: one
+# coordinate line and one exponent line per atom.  The exponent lines are
+# the Slater exponents of the valence basis, which is what a future
+# 3-D orbital renderer would need; nothing reads them yet, so they are
+# deliberately not emitted as a property rather than shipped unused.
+#
+# Energies are eV here and Hartree in ECCE (MoPanel hard-codes the column
+# label), so they are converted.
+#
+# The symmetry label is a single token such as "1a1", "2b2" or "1a''" --
+# an index within the irrep followed by the irrep itself.  Only the irrep
+# is kept: MoPanel::fillGraph() groups the orbital plot by this string,
+# so keeping the index would put every orbital in a group of its own.
+#
+################################################################################
+
+$| = 1;
+
+use constant EV_PER_HARTREE => 27.211386245988;
+
+my @lines;
+while (<STDIN>) { chomp; push(@lines, $_); }
+
+exit(0) if (!@lines);
+
+#  Line 1 begins with the atom count.
+my $natoms = 0;
+$natoms = $1 if ($lines[0] =~ /^\s*(\d+)\s/);
+exit(0) if ($natoms <= 0);
+
+#  Coordinates then exponents, one line each per atom.  The number of
+#  basis functions is not stated anywhere, so it is taken from the
+#  coefficient count of the first orbital instead.
+my $idx = 1 + 2 * $natoms;
+
+my (@occ, @eng, @sym, @coeff);
+while ($idx <= $#lines) {
+  my $line = $lines[$idx];
+  if ($line =~ /^\s*ORBITAL\s+(\S+)\s+(\S+)\s+(-?[\d.]+)\s*$/) {
+    my ($o, $s, $e) = ($1, $2, $3);
+    #  "1a1" -> "a1"; "2a''" -> "a''".  A label that is nothing but the
+    #  irrep is left alone.
+    $s =~ s/^\d+//;
+    push(@occ, $o);
+    push(@eng, $e / EV_PER_HARTREE);
+    push(@sym, $s);
+    push(@coeff, []);
+    $idx++;
+    #  Coefficients run until the next line that is not purely numbers.
+    #
+    #  Terminating on a whitelist of what a coefficient line looks like,
+    #  rather than a blacklist of what follows, because more follows than
+    #  is obvious: after the last orbital MOPAC writes an
+    #  "INVERSE_MATRIX[6x6]=" block of the same D-notation numbers, and a
+    #  loop that stopped only at "ORBITAL" or "Keywords:" swallowed all
+    #  of it into the final orbital -- 27 coefficients where the other
+    #  five had 6.  That then failed the rectangular check below and the
+    #  script emitted nothing at all, silently.
+    while ($idx <= $#lines &&
+           $lines[$idx] =~ /^[-+\s\d.DdEe]+$/ &&
+           $lines[$idx] =~ /\d[DdEe][-+]?\d/) {
+      my $text = $lines[$idx];
+      #  Fortran D exponents, and values packed without separators:
+      #  "-0.63614535D-01-0.45803202D-14" is two numbers, not one.
+      while ($text =~ /(-?\d*\.\d+[DdEe][-+]?\d+)/g) {
+        my $v = $1;
+        $v =~ s/[Dd]/e/;
+        push(@{$coeff[-1]}, $v + 0);
+      }
+      $idx++;
+    }
+    next;
+  }
+  $idx++;
+}
+
+my $nmo = scalar(@eng);
+exit(0) if ($nmo == 0);
+
+my $nbas = scalar(@{$coeff[0]});
+exit(0) if ($nbas == 0);
+#  Every orbital must have the same number of coefficients, or the table
+#  is not rectangular and PropTable rejects it on load with "input vector
+#  length does not match rows*columns" -- which presents as the Orbitals
+#  panel refusing to compute, not as a parse error.
+foreach my $c (@coeff) {
+  exit(0) if (scalar(@$c) != $nbas);
+}
+
+#  UHF writes the alpha set then the beta set into the same file, so twice
+#  the basis size means both are present.
+my $half = ($nmo == 2 * $nbas) ? $nbas : $nmo;
+
+sub emitSet {
+  my ($from, $count, $moKey, $engKey, $occKey, $symKey) = @_;
+  return if ($count <= 0);
+
+  print "key: $engKey\n";
+  print "size:\n$count\n";
+  print "values:\n";
+  printf("%.8f ", $eng[$from + $_]) foreach (0 .. $count - 1);
+  print "\nunits:\nHartree\n";
+  print "END\n";
+
+  print "key: $occKey\n";
+  print "size:\n$count\n";
+  print "values:\n";
+  printf("%.4f ", $occ[$from + $_]) foreach (0 .. $count - 1);
+  print "\nunits:\nelectrons\n";
+  print "END\n";
+
+  print "key: $symKey\n";
+  print "size:\n$count\n";
+  print "rowlabels:\n";
+  print " " . ($_ + 1) . " " foreach (0 .. $count - 1);
+  print "\nvalues:\n";
+  print join(" ", @sym[$from .. $from + $count - 1]), "\n";
+  print "END\n";
+
+  print "key: $moKey\n";
+  print "size:\n$count $nbas\n";
+  print "rowlabels:\n";
+  for my $m (1 .. $count) {
+    print "MO-$m ";
+    print "\n" if (($m % 8) == 0);
+  }
+  print "\ncolumnlabels:\n";
+  for my $b (1 .. $nbas) {
+    print "BasFun-$b ";
+    print "\n" if (($b % 6) == 0);
+  }
+  print "\nvalues:\n";
+  for my $m (0 .. $count - 1) {
+    my $row = $coeff[$from + $m];
+    printf("%.8f ", $_) foreach (@$row);
+    print "\n";
+  }
+  print "END\n";
+}
+
+&emitSet(0, $half, "MO", "ORBENG", "ORBOCC", "ORBSYM");
+&emitSet($half, $nmo - $half,
+         "MOBETA", "ORBENGBETA", "ORBOCCBETA", "ORBSYMBETA")
+  if ($nmo > $half);
+exit(0);

scripts/parsers/mopac.mgf

@@ -75,6 +75,39 @@ properties first — a coefficient matrix exceeds `msgLength` and arrives
 split, so treating each packet as a message truncates exactly the
 properties most worth testing.
 
+## File= parse types
+
+Some codes do not put their best data in the job output at all. Gaussian
+writes MO coefficients to `fort.7`, GAMESS-UK to `ftn058`, MOPAC to the
+GRAPHF file. Those are `File=` parse types, and they work differently:
+eccejobmonitor never matches them against the output. It sends a
+**jmFILE** message naming the file and the parse type, and the client
+fetches that file and runs the script on it. The message carries no
+content.
+
+`pipeline.unpack_files()` decodes those announcements and `run_parsers()`
+acts on them, after the matched blocks -- which is the monitor's own
+order, `JobFilesGet()` running after `JobOutputGet()`. That ordering is
+load-bearing where both kinds of entry emit the same keys: MOPAC's
+GRAPHF file is the better source for orbitals, but only newer decks
+request it, so the printed-eigenvector entry remains as a fallback and
+the file wins when present.
+
+Two traps found wiring this up:
+
+- **`importDir`**. In anything but live mode `MsgSendFile` resolves the
+  name against `$idir$file` -- concatenated, no separator -- and `idir`
+  defaults to `"-"`. Left unset, every `File=` parse type silently looks
+  for `-fort.7`, finds nothing, and is skipped with no error. The
+  harness passes the work directory with a trailing separator.
+- **The deck's name**. MOPAC derives every output name from the input
+  name, so a fixture called `ch4_mos.mop` writes `ch4_mos.mgf` while
+  `mopac.desc` names `mopac.mgf` -- correct, because ECCE always writes
+  `mopac.mop`. The fixtures are therefore named the way ECCE names them,
+  and `mopac-ch4-mos` asserts the GRAPHF file really was the source by
+  checking the precision of the energies: eight decimals from the file,
+  six from the printed block.
+
 ## Assertion style
 
 Loose about exact numbers, strict about structure. Pinning an energy to

tests/e2e/README.md

@@ -64,6 +64,22 @@ def expect_mopac_ch4(props, report):
     report.check(occ[:4] == [2.0] * 4 and occ[4:] == [0.0] * 4,
                  'four filled levels then virtuals')
 
+    #  THE GRAPHF FILE MUST BE THE SOURCE, not the printed eigenvector
+    #  block.  Both parse types emit these keys and the File= one is
+    #  delivered last, so it wins -- but only if MOPAC actually wrote it
+    #  and mopac.desc's File= name matched, neither of which is visible
+    #  from the values alone.
+    #
+    #  The tell is precision: mopac.mo prints six decimals because that
+    #  is all the printed block has, mopac.mgf prints eight because the
+    #  GRAPHF file carries Fortran doubles.  A value that stops at six
+    #  means the fallback ran and the better source was missed.
+    decimals = max(len(t.split('.')[1]) for t in
+                   sect(orbeng[-1], 'values').split() if '.' in t)
+    report.check(decimals >= 8,
+                 'ORBENG came from the GRAPHF file, not the printed block '
+                 '(%d decimals; 6 means the fallback ran)' % decimals)
+
     #  Energies must be in Hartree, not the eV MOPAC prints.  MoPanel
     #  hard-codes its column label to "Energy Hartree", so a missing
     #  conversion shows numbers 27x too large with no error at all.
@@ -178,20 +194,20 @@ def expect_nwchem_h2o(props, report):
         name='mopac-ch4-mos',
         code='mopac',
         desc='mopac.desc',
-        deck='mopac/ch4_mos.mop',
+        deck='mopac/mos/mopac.mop',
         #  ". <runtype> <theory category> <theory name> <open shells>",
         #  as DavCalculation::getParseScriptArgs builds it.
         parse_args=('.', 'Energy', 'SE', 'PM7', '0'),
-        output='ch4_mos.out',
+        output='mopac.out',
         expect=expect_mopac_ch4,
     ),
     dict(
         name='mopac-ch4-thermo',
         code='mopac',
         desc='mopac.desc',
-        deck='mopac/ch4_thermo.mop',
+        deck='mopac/thermo/mopac.mop',
         parse_args=('.', 'Vibration', 'SE', 'PM7', '0'),
-        output='ch4_thermo.out',
+        output='mopac.out',
         expect=expect_mopac_ch4_thermo,
     ),
     dict(

tests/e2e/cases.py

@@ -1,4 +1,4 @@
-PM7 VECTORS ALLVEC CHARGE=0
+PM7 VECTORS ALLVEC GRAPHF CHARGE=0
 ECCE MOPAC calculation
 
 C   0.000000 1  0.000000 1  0.000000 1