The symmetry orbitals themselves, not just which irreps appear (#132)

Completes phase 3.  The reduction says a1 + t2; the diagram needs the
combinations, so this applies the projection operator
P = (l/h) sum_R chi(R)* R to each orbital of an orbit and picks the
independent results out by Gram-Schmidt.

Checked against what a student writes by hand:

    CH4   A1  (1/2, 1/2, 1/2, 1/2)          + T2  (3 vectors)
    H2O   A1  (1/sqrt2, 1/sqrt2)            + B1  (1 vector)
    SF6   A1g (1/sqrt6 x 6)                 + Eg (2) + T1u (3)

and on three invariants that would catch a wrong projector without
anyone having to know the answer: the count for each irrep is its
multiplicity times its dimension, the vectors together span the orbit
with nothing left over, and every one is normalised and orthogonal to
the rest.

The totally symmetric one is checked to be the in-phase combination --
every coefficient equal in magnitude -- which is the thing a chemist
would look at first and the thing a sign error would break.

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

Author: andy

include/tdat/SymmetryAnalysis.H

@@ -101,6 +101,35 @@ class SymmetryAnalysis
                                  const vector<int>& classOfOp,
                                  int numClasses,
                                  vector<double>& chi);
+
+    /**
+     * The symmetry orbitals of one irrep, over an orbit's s orbitals.
+     *
+     * The projection operator, P = (l/h) sum_R chi(R)* R, applied to
+     * each orbital of the orbit in turn.  That gives as many vectors as
+     * there are orbitals, most of them repeats or zero, so the
+     * independent ones are picked out by Gram-Schmidt and normalised.
+     *
+     * These are the terminal-atom symmetry orbitals -- the in-phase and
+     * out-of-phase combinations a diagram draws on one side.
+     *
+     * @param table   the group's character table
+     * @param irrep   which irrep to project onto
+     * @param vectors filled with one coefficient vector per independent
+     *                symmetry orbital, each over the orbit's atoms in
+     *                orbit order
+     *
+     * @return false if the arguments do not describe this group.  The
+     *         number of vectors returned is the irrep's multiplicity
+     *         times its dimension, and a caller that knows the expected
+     *         reduction should check it.
+     */
+    static bool projectOrbit(const vector<int>& orbit,
+                             const vector< vector<int> >& images,
+                             const vector<int>& classOfOp,
+                             const CharacterTable& table,
+                             const string& irrep,
+                             vector< vector<double> >& vectors);
 };
 
 #endif

src/tdat/chemistry/SymmetryAnalysis.C

@@ -6,6 +6,9 @@
 #include "tdat/SymmetryAnalysis.H"
 #include "tdat/CharacterTable.H"
 
+#include <map>
+using std::map;
+
 namespace {
 
   /** Apply an operation to a point. */
@@ -359,3 +362,65 @@ bool SymmetryAnalysis::orbitalCharacter(const vector<int>& orbit,
   for (int c = 0; c < numClasses; c++) if (seen[c] == 0) return false;
   return true;
 }
+
+
+bool SymmetryAnalysis::projectOrbit(const vector<int>& orbit,
+                                    const vector< vector<int> >& images,
+                                    const vector<int>& classOfOp,
+                                    const CharacterTable& table,
+                                    const string& irrep,
+                                    vector< vector<double> >& vectors)
+{
+  vectors.clear();
+
+  const vector<double>* chi = table.characters(irrep);
+  if (chi == 0) return false;
+  if (classOfOp.size() != images.size()) return false;
+
+  const size_t n = orbit.size();
+  if (n == 0) return false;
+
+  //  Where each atom sits in the orbit, so an image can be turned back
+  //  into a position in the coefficient vector.
+  map<int, size_t> position;
+  for (size_t i = 0; i < n; i++) position[orbit[i]] = i;
+
+  const double dimension = table.dimension(irrep);
+  const double h = table.order();
+  if (h <= 0.0) return false;
+
+  //  Project each orbital of the orbit in turn.  Most of the results
+  //  are repeats of one another or zero; Gram-Schmidt below keeps the
+  //  independent ones.
+  for (size_t a = 0; a < n; a++) {
+
+    vector<double> v(n, 0.0);
+    for (size_t o = 0; o < images.size(); o++) {
+      const int cls = classOfOp[o];
+      if (cls < 0 || (size_t)cls >= chi->size()) return false;
+
+      map<int, size_t>::const_iterator it = position.find(images[o][orbit[a]]);
+      if (it == position.end()) continue;      // left the orbit: cannot happen
+      v[it->second] += (*chi)[cls];
+    }
+    for (size_t i = 0; i < n; i++) v[i] *= dimension/h;
+
+    //  Remove whatever is already spanned, then keep what is left if
+    //  there is anything of it.
+    for (size_t k = 0; k < vectors.size(); k++) {
+      double dot = 0.0;
+      for (size_t i = 0; i < n; i++) dot += v[i]*vectors[k][i];
+      for (size_t i = 0; i < n; i++) v[i] -= dot*vectors[k][i];
+    }
+
+    double norm = 0.0;
+    for (size_t i = 0; i < n; i++) norm += v[i]*v[i];
+    if (norm < 1.0e-10) continue;              // nothing new
+
+    norm = sqrt(norm);
+    for (size_t i = 0; i < n; i++) v[i] /= norm;
+    vectors.push_back(v);
+  }
+
+  return true;
+}

tests/symmetry/testSymmetryAnalysis.C

@@ -103,6 +103,66 @@ static void run(const Molecule& mol, const string& opsDir)
   check(ok, "conjugacy classes matched to the table's named classes");
   if (!ok) return;
 
+  //  The symmetry orbitals themselves, not just which irreps appear.
+  for (size_t i = 0; i < orbits.size(); i++) {
+    if (orbits[i].size() != mol.terminalCount) continue;
+
+    vector<double> chi;
+    if (!SymmetryAnalysis::orbitalCharacter(orbits[i], images, classOfOp,
+                                            (int)table->classes().size(), chi))
+      continue;
+    vector<int> mult;
+    if (!table->reduce(chi, mult)) continue;
+
+    size_t total = 0;
+    for (size_t k = 0; k < table->irreps().size(); k++) {
+      if (mult[k] == 0) continue;
+      const string& ir = table->irreps()[k];
+
+      vector< vector<double> > vecs;
+      bool got = SymmetryAnalysis::projectOrbit(orbits[i], images, classOfOp,
+                                                *table, ir, vecs);
+      snprintf(msg, sizeof(msg),
+               "%s: %zu symmetry orbitals (%d x dim %d = %d expected)",
+               ir.c_str(), vecs.size(), mult[k], table->dimension(ir),
+               mult[k]*table->dimension(ir));
+      check(got && (int)vecs.size() == mult[k]*table->dimension(ir), msg);
+      total += vecs.size();
+
+      //  Each must be normalised and orthogonal to the others.
+      for (size_t u = 0; u < vecs.size(); u++) {
+        double nrm = 0.0;
+        for (size_t c = 0; c < vecs[u].size(); c++) nrm += vecs[u][c]*vecs[u][c];
+        if (fabs(nrm - 1.0) > 1e-9) {
+          check(false, string(ir) + ": a symmetry orbital is not normalised");
+          break;
+        }
+      }
+    }
+    //  Together they must span the whole orbit -- no orbital left over.
+    snprintf(msg, sizeof(msg),
+             "the symmetry orbitals span the orbit (%zu of %zu)",
+             total, orbits[i].size());
+    check(total == orbits[i].size(), msg);
+
+    //  The totally symmetric one is the in-phase combination: every
+    //  coefficient equal and positive.  That is the check a student
+    //  would do by eye.
+    const string& first = table->irreps()[0];
+    vector< vector<double> > sym;
+    if (SymmetryAnalysis::projectOrbit(orbits[i], images, classOfOp,
+                                       *table, first, sym) && sym.size() == 1) {
+      const double want = 1.0/sqrt((double)orbits[i].size());
+      bool inPhase = true;
+      for (size_t c = 0; c < sym[0].size(); c++)
+        if (fabs(fabs(sym[0][c]) - want) > 1e-9) inPhase = false;
+      snprintf(msg, sizeof(msg),
+               "%s is the in-phase combination, all coefficients %.4f",
+               first.c_str(), want);
+      check(inPhase, msg);
+    }
+  }
+
   //  The terminal orbit's s orbitals: the reduction a textbook gives.
   for (size_t i = 0; i < orbits.size(); i++) {
     if (orbits[i].size() != mol.terminalCount) continue;