cofactor sometimes return zero overlap for one case, but non-zero for the conjugate

[jagot]

In some instances, the generated energy expression is not Hermitian. I tracked this down to cofactor returning zero determinantal overlap between Slater determinants a and b, but non-zero for the conjugate case b and a:

using EnergyExpressions
import EnergyExpressions: nonzero_minors, cofactor

using Combinatorics

function test_case(N,nn)
    orbitals = 1:N
    excite = i -> i+100
    virtuals = map(excite, orbitals)

    configurations = [SlaterDeterminant(vcat(1:i-1,excite(i),i+1:N))
                      for i=1:N]
    overlaps = map(((a,b),) -> OrbitalOverlap(a,b),
                   reduce(vcat, [collect(combinations(virtuals,2)),
                                 [[a,a] for a in virtuals]]))

    for (i,j) = [(N-nn,N),(N,N-nn)]
        println(repeat("-", 100))
        a = configurations[i]
        b = configurations[j]
        S = overlap_matrix(a, b, overlaps)
        @info "Overlap matrix" S
        ks,ls = nonzero_minors(1, S)
        for (k,l) in zip(ks,ls)
            @show k,l
            Dkl = cofactor(k, l, S)
            @show Dkl
        end
    end
end

E.g. this example fails:

julia> test_case(9,3)
----------------------------------------------------------------------------------------------------
┌ Info: Overlap matrix
│   S =
│    9×9 SparseArrays.SparseMatrixCSC{EnergyExpressions.NBodyTerm,Int64} with 8 stored entries:
│      [1, 1]  =  1
│      [2, 2]  =  1
│      [3, 3]  =  1
│      [4, 4]  =  1
│      [5, 5]  =  1
│      [7, 7]  =  1
│      [8, 8]  =  1
└      [6, 9]  =  ⟨106|109⟩
(k, l) = ([9], [6])
Dkl = - ⟨106|109⟩
----------------------------------------------------------------------------------------------------
┌ Info: Overlap matrix
│   S =
│    9×9 SparseArrays.SparseMatrixCSC{EnergyExpressions.NBodyTerm,Int64} with 8 stored entries:
│      [1, 1]  =  1
│      [2, 2]  =  1
│      [3, 3]  =  1
│      [4, 4]  =  1
│      [5, 5]  =  1
│      [9, 6]  =  ⟨109|106⟩
│      [7, 7]  =  1
└      [8, 8]  =  1
(k, l) = ([6], [9])
Dkl = 0

whereas this very similar case works as expected:

julia> test_case(9,2)
----------------------------------------------------------------------------------------------------
┌ Info: Overlap matrix
│   S =
│    9×9 SparseArrays.SparseMatrixCSC{EnergyExpressions.NBodyTerm,Int64} with 8 stored entries:
│      [1, 1]  =  1
│      [2, 2]  =  1
│      [3, 3]  =  1
│      [4, 4]  =  1
│      [5, 5]  =  1
│      [6, 6]  =  1
│      [8, 8]  =  1
└      [7, 9]  =  ⟨107|109⟩
(k, l) = ([9], [7])
Dkl = - ⟨107|109⟩
----------------------------------------------------------------------------------------------------
┌ Info: Overlap matrix
│   S =
│    9×9 SparseArrays.SparseMatrixCSC{EnergyExpressions.NBodyTerm,Int64} with 8 stored entries:
│      [1, 1]  =  1
│      [2, 2]  =  1
│      [3, 3]  =  1
│      [4, 4]  =  1
│      [5, 5]  =  1
│      [6, 6]  =  1
│      [9, 7]  =  ⟨109|107⟩
└      [8, 8]  =  1
(k, l) = ([7], [9])
Dkl = - ⟨109|107⟩