private static IEnumerable ToTwoBodySpinOrbitalTerms()

in Chemistry/src/DataModel/OrbitalIntegral/OrbitalIntegralExtensions.cs [124:231]
81 lines of code
28 McCabe index (conditional complexity)

        private static IEnumerable<(HermitianFermionTerm, double)> ToTwoBodySpinOrbitalTerms(
            this OrbitalIntegral orbitalIntegral,
            int nOrbitals,
            IndexConvention indexConvention)
        {
            // Two-electron orbital integral symmetries
            // ijkl = lkji = jilk = klij = ikjl = ljki = kilj = jlik.
            var pqrsSpinOrbitals = orbitalIntegral.EnumerateOrbitalSymmetries().EnumerateSpinOrbitals();
            var coefficient = orbitalIntegral.Coefficient;


            // We only need to see one of these.
            // Now iterate over pqrsArray
            foreach (var pqrs in pqrsSpinOrbitals)
            {
                var p = pqrs[0];
                var q = pqrs[1];
                var r = pqrs[2];
                var s = pqrs[3];

                var pInt = p.ToInt(indexConvention, nOrbitals);
                var qInt = q.ToInt(indexConvention, nOrbitals);
                var rInt = r.ToInt(indexConvention, nOrbitals);
                var sInt = s.ToInt(indexConvention, nOrbitals);

                // Only consider terms on the lower diagonal due to Hermitian symmetry.

                // For terms with two different orbital indices, possibilities are
                // PPQQ (QQ = 0), PQPQ, QPPQ (p<q), PQQP, QPQP (p<q), QQPP (PP=0)
                // Hence, if we only count PQQP, and PQPQ, we need to double the coefficient.
                // iU jU jU iU | iU jD jD iD | iD jU jU iD | iD jD jD iD
                if (pInt == sInt && qInt == rInt && pInt < qInt)
                {   // PQQP
                    yield return (new HermitianFermionTerm(new[] { pInt, qInt, rInt, sInt }.ToLadderSequence()), 1.0 * coefficient );
                }
                else if (pInt == rInt && qInt == sInt && pInt < qInt)
                {
                    // iU jU iU jU | iD jD iD jD
                    // PQPQ
                    yield return (new HermitianFermionTerm(new[] { pInt, qInt, sInt, rInt }.ToLadderSequence()), -1.0 * coefficient );
                }
                else if (qInt == rInt && pInt < sInt && rInt != sInt && pInt != qInt)
                {
                    // PQQR
                    // For any distinct pqr, [i;j;j;k] generates PQQR ~ RQQP ~ QPRQ ~ QRPQ. We only need to record one.
                    if (rInt < sInt)
                    {
                        if (pInt < qInt)
                        {
                            yield return (new HermitianFermionTerm(new[]  { pInt, qInt, sInt, rInt }.ToLadderSequence()), -2.0 * coefficient );
                        }
                        else
                        {
                            yield return (new HermitianFermionTerm(new[]  { qInt, pInt, sInt, rInt }.ToLadderSequence()), 2.0 * coefficient );
                        }

                    }
                    else
                    {
                        if (pInt < qInt)
                        {
                            yield return (new HermitianFermionTerm(new[]  { pInt, qInt, rInt, sInt }.ToLadderSequence()), 2.0 * coefficient );
                        }
                        else
                        {
                            yield return (new HermitianFermionTerm(new[]  { qInt, pInt, rInt, sInt }.ToLadderSequence()), -2.0 * coefficient );
                        }
                    }
                }
                else if (qInt == sInt && pInt < rInt && rInt != sInt && pInt != sInt)
                {
                    // PQRQ
                    // For any distinct pqr, [i;j;k;j] generates {p, q, r, q}, {q, r, q, p}, {q, p, q, r}, {r, q, p, q}. We only need to record one.
                    if (pInt < qInt)
                    {
                        if (rInt > qInt)
                        {
                            yield return (new HermitianFermionTerm(new[]  { pInt, qInt, rInt, sInt }.ToLadderSequence()), 2.0 * coefficient );
                        }
                        else
                        {
                            yield return (new HermitianFermionTerm(new[]  { pInt, qInt, sInt, rInt }.ToLadderSequence()), -2.0 * coefficient );
                        }
                    }
                    else
                    {
                        yield return (new HermitianFermionTerm(new[]  { qInt, pInt, rInt, sInt }.ToLadderSequence()), -2.0 * coefficient );
                    }
                }
                else if (pInt < qInt && pInt < rInt && pInt < sInt && qInt != rInt && qInt != sInt && rInt != sInt)
                {
                    // PQRS
                    // For any distinct pqrs, [i;j;k;l] generates 
                    // {{p, q, r, s}<->{s, r, q, p}<->{q, p, s, r}<->{r, s, p, q}, 
                    // {1,2,3,4}<->{4,3,2,1}<->{2,1,4,3}<->{3,4,1,2}
                    // {p, r, q, s}<->{s, q, r, p}<->{r, p, s, q}<->{q, s, p, r}}
                    // 1324, 4231, 3142, 2413
                    if (rInt < sInt)
                    {
                        yield return (new HermitianFermionTerm(new[]  { pInt, qInt, sInt, rInt }.ToLadderSequence()), -2.0 * coefficient );
                    }
                    else
                    {
                        yield return (new HermitianFermionTerm(new[]  { pInt, qInt, rInt, sInt }.ToLadderSequence()), 2.0 * coefficient );
                    }
                }
            }
        }