The MO diagram's layout: levels, core folding, and symmetry connections

The drawing is the product, and this is the part of it with judgement in
it -- what shares a level, what counts as core, what connects to what.
Kept separate from anything that paints so it can be tested without a
display; the panel will only render the result.

Three decisions, each pinned by a test:

LEVELS.  Orbitals share a level when they share an irrep AND an energy.
Both conditions matter, and the tests cover both directions: two irreps
that happen to be degenerate stay two levels -- an accidental
near-degeneracy is not a degenerate set -- and one irrep at two energies
stays two.  A level sits at the mean of its orbitals, so a set split by
numerical noise is drawn where it belongs.  It keeps its orbital
indices, which is what makes a level clickable through to the viewer.

CORE.  Folded into a counted band rather than drawn, because thirty core
levels squeeze the valence region into nothing.  The cutoff is taken
from the biggest gap in the lower half, judged against the TYPICAL gap
rather than the total span: a span-based rule is scale-dependent and got
the easy case wrong, since an evenly spaced spectrum of five levels has
every gap at exactly a quarter of the span.  A core sits four or more
typical gaps from the valence; even spacing sits at one.  That also
gives the right answer for a semiempirical calculation, which has no
core: nothing is hidden, rather than something arbitrary being hidden.

The count is kept, so the diagram can say how many are folded away.  An
absence the reader cannot see is worse than a cluttered diagram.

CONNECTIONS.  Symmetry is the entire rule: same irrep may combine,
different may not.  A molecular orbital with no fragment partner of its
irrep is non-bonding, and is drawn as such -- water's b2 lone pair is
exactly that case, and is the test for it.

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

Author: andy

include/tdat/MoDiagram.H

@@ -0,0 +1,124 @@
+/**
+ * @file
+ *
+ * The layout of a molecular orbital correlation diagram.
+ *
+ * Deliberately separate from anything that draws: what to put where is
+ * decided here, in terms of levels, labels and connections, and the
+ * panel only paints the result.  That keeps the part with the judgement
+ * in it -- which orbitals share a level, which count as core, what
+ * connects to what -- testable without a display.
+ */
+
+#ifndef MODIAGRAM_HH
+#define MODIAGRAM_HH
+
+#include <string>
+#include <vector>
+using std::string;
+using std::vector;
+
+/** One drawn level: a set of orbitals sharing an energy and a label. */
+struct MoLevel
+{
+   double         energy;          // Hartree
+   string         label;           // the irrep, as the code reports it
+   double         occupancy;       // electrons in the whole level
+   int            degeneracy;      // orbitals drawn as this one level
+   vector<int>    orbitals;        // their indices, for clicking through
+
+   MoLevel() : energy(0.0), occupancy(0.0), degeneracy(0) {}
+};
+
+
+/** One column of the diagram. */
+struct MoColumn
+{
+   string           title;
+   vector<MoLevel>  levels;
+
+   /**
+    * Orbitals folded away below the energy cutoff.
+    *
+    * Drawn as a single band rather than as levels: a molecule with
+    * thirty core orbitals squeezes its valence region into nothing
+    * otherwise.  Kept as a count so the diagram can SAY how many are
+    * hidden -- an absence the reader cannot see is worse than a
+    * cluttered diagram.
+    */
+   int    hiddenCount;
+   double hiddenMaxEnergy;
+
+   MoColumn() : hiddenCount(0), hiddenMaxEnergy(0.0) {}
+};
+
+
+/** A line between a fragment level and a molecular orbital. */
+struct MoConnection
+{
+   int leftLevel;     // index into the left column, or -1
+   int centreLevel;   // index into the centre column
+   int rightLevel;    // index into the right column, or -1
+};
+
+
+class MoDiagram
+{
+  public:
+
+    MoDiagram();
+
+    /**
+     * Group orbitals into drawn levels.
+     *
+     * Orbitals are one level when they share an irrep label AND lie
+     * within tolerance in energy.  The degeneracy that results is
+     * reported, but a caller that knows the irrep's dimension should
+     * prefer that: an accidental near-degeneracy is not a degenerate
+     * set, and grouping by energy alone cannot tell them apart.
+     *
+     * @param energies     one per orbital, Hartree, ascending
+     * @param occupancies  one per orbital
+     * @param labels       one per orbital; may be empty, in which case
+     *                     grouping is by energy alone and every level
+     *                     is labelled unknown
+     * @param tolerance    Hartree
+     */
+    static bool group(const vector<double>& energies,
+                      const vector<double>& occupancies,
+                      const vector<string>& labels,
+                      double tolerance,
+                      vector<MoLevel>& levels);
+
+    /**
+     * Fold levels below a cutoff into the column's hidden band.
+     *
+     * @param cutoff  Hartree; levels below this are hidden
+     */
+    static void hideBelow(MoColumn& column, double cutoff);
+
+    /**
+     * A cutoff taken from the biggest gap in the lower half of the
+     * spectrum, so it adapts to the molecule instead of assuming what
+     * counts as core.
+     *
+     * Returns the lowest energy in the column when no useful gap is
+     * found, which hides nothing.
+     */
+    static double suggestCoreCutoff(const vector<MoLevel>& levels);
+
+    /**
+     * Connect fragment levels to molecular orbitals by matching irrep.
+     *
+     * Symmetry is the whole rule: orbitals of different irreps cannot
+     * mix, and orbitals of the same one must be allowed to.  A centre
+     * level with no fragment partner of its irrep is non-bonding and
+     * gets a connection with -1 on that side.
+     */
+    static void connect(const vector<MoLevel>& left,
+                        const vector<MoLevel>& centre,
+                        const vector<MoLevel>& right,
+                        vector<MoConnection>& connections);
+};
+
+#endif

src/tdat/chemistry/MoDiagram.C

@@ -0,0 +1,156 @@
+#include <cmath>
+#include <cstddef>
+
+#include "tdat/MoDiagram.H"
+
+MoDiagram::MoDiagram()
+{
+}
+
+
+bool MoDiagram::group(const vector<double>& energies,
+                      const vector<double>& occupancies,
+                      const vector<string>& labels,
+                      double tolerance,
+                      vector<MoLevel>& levels)
+{
+  levels.clear();
+
+  const size_t n = energies.size();
+  if (n == 0) return false;
+  if (!occupancies.empty() && occupancies.size() != n) return false;
+  if (!labels.empty() && labels.size() != n) return false;
+
+  for (size_t i = 0; i < n; i++) {
+
+    const string label = labels.empty() ? string("?") : labels[i];
+
+    //  Same level as the one before it?  Both conditions matter: two
+    //  orbitals of DIFFERENT irreps can be accidentally degenerate and
+    //  are not one level, and two of the same irrep at different
+    //  energies are two levels.
+    bool extend = false;
+    if (!levels.empty()) {
+      MoLevel& last = levels.back();
+      extend = (last.label == label)
+               && (fabs(energies[i] - last.energy) <= tolerance);
+    }
+
+    if (extend) {
+      MoLevel& last = levels.back();
+      //  The level sits at the mean of its orbitals, so a set split by
+      //  numerical noise is drawn where it belongs rather than at
+      //  whichever of them happened to come first.
+      last.energy = (last.energy*last.degeneracy + energies[i])
+                    / (last.degeneracy + 1);
+      last.degeneracy += 1;
+      last.occupancy += occupancies.empty() ? 0.0 : occupancies[i];
+      last.orbitals.push_back((int)i);
+    } else {
+      MoLevel level;
+      level.energy = energies[i];
+      level.label = label;
+      level.occupancy = occupancies.empty() ? 0.0 : occupancies[i];
+      level.degeneracy = 1;
+      level.orbitals.push_back((int)i);
+      levels.push_back(level);
+    }
+  }
+  return true;
+}
+
+
+void MoDiagram::hideBelow(MoColumn& column, double cutoff)
+{
+  vector<MoLevel> kept;
+  column.hiddenCount = 0;
+  column.hiddenMaxEnergy = 0.0;
+  bool anyHidden = false;
+
+  for (size_t i = 0; i < column.levels.size(); i++) {
+    if (column.levels[i].energy < cutoff) {
+      column.hiddenCount += column.levels[i].degeneracy;
+      if (!anyHidden || column.levels[i].energy > column.hiddenMaxEnergy) {
+        column.hiddenMaxEnergy = column.levels[i].energy;
+      }
+      anyHidden = true;
+    } else {
+      kept.push_back(column.levels[i]);
+    }
+  }
+  column.levels = kept;
+}
+
+
+double MoDiagram::suggestCoreCutoff(const vector<MoLevel>& levels)
+{
+  if (levels.size() < 3) return levels.empty() ? -1.0e30
+                                               : levels.front().energy;
+
+  //  Every gap, and the typical one.
+  vector<double> gaps;
+  for (size_t i = 1; i < levels.size(); i++) {
+    gaps.push_back(levels[i].energy - levels[i-1].energy);
+  }
+  vector<double> sorted(gaps);
+  for (size_t i = 0; i < sorted.size(); i++) {
+    for (size_t j = i+1; j < sorted.size(); j++) {
+      if (sorted[j] < sorted[i]) {
+        const double t = sorted[i]; sorted[i] = sorted[j]; sorted[j] = t;
+      }
+    }
+  }
+  const double median = sorted[sorted.size()/2];
+
+  //  The biggest gap in the lower half, where a core/valence split
+  //  would be.
+  const size_t limit = levels.size()/2;
+  double bestGap = 0.0;
+  size_t bestAt = 0;
+  for (size_t i = 1; i <= limit && i < levels.size(); i++) {
+    if (gaps[i-1] > bestGap) { bestGap = gaps[i-1]; bestAt = i; }
+  }
+
+  //  Judged against the TYPICAL gap, not against the total span.  A
+  //  span-based rule is scale-dependent and gets the easy case wrong:
+  //  an evenly spaced spectrum of five levels has every gap at exactly
+  //  a quarter of the span, so any threshold near that either hides an
+  //  orbital for no reason or misses a real core.  A core sits four or
+  //  more typical gaps away from the valence; even spacing sits at one.
+  //
+  //  This also gives the right answer for a semiempirical calculation,
+  //  which has no core at all: its gaps are comparable, so nothing is
+  //  hidden rather than something arbitrary being hidden.
+  if (median <= 0.0 || bestGap < 4.0*median) return levels.front().energy;
+
+  return 0.5*(levels[bestAt].energy + levels[bestAt-1].energy);
+}
+
+
+void MoDiagram::connect(const vector<MoLevel>& left,
+                        const vector<MoLevel>& centre,
+                        const vector<MoLevel>& right,
+                        vector<MoConnection>& connections)
+{
+  connections.clear();
+
+  for (size_t c = 0; c < centre.size(); c++) {
+
+    MoConnection link;
+    link.leftLevel = -1;
+    link.centreLevel = (int)c;
+    link.rightLevel = -1;
+
+    //  Symmetry is the entire rule: only orbitals of the same irrep may
+    //  combine.  A centre level with no partner on a side is
+    //  non-bonding from that side, which is a result worth drawing --
+    //  water's b2 lone pair is exactly this.
+    for (size_t l = 0; l < left.size(); l++) {
+      if (left[l].label == centre[c].label) { link.leftLevel = (int)l; break; }
+    }
+    for (size_t r = 0; r < right.size(); r++) {
+      if (right[r].label == centre[c].label) { link.rightLevel = (int)r; break; }
+    }
+    connections.push_back(link);
+  }
+}

tests/symmetry/run_tests.py

@@ -529,6 +529,11 @@ def main():
         print("FAILED  the orbit/reduction analysis")
         return 1
 
+    if standalone("testMoDiagram",
+                  ["src/tdat/chemistry/MoDiagram.C"]) != 0:
+        print("FAILED  the diagram layout")
+        return 1
+
     print("PASSED")
     return 0