Errors in multipole expansion

[jagot]

When deriving energy expressions for open-shell configurations, with p electrons present, the correct terms are present, but their prefactors seem to wrong (normalization of spherical harmonics?) and also their signs:

julia> using AtomicLevels

julia> using AngularMomentumAlgebra

julia> cfgs = spin_configurations([c"1s 2p"])
12-element Array{Configuration{SpinOrbital{Orbital{Int64}}},1}:
 1s₀α 2p₋₁α
 1s₀α 2p₋₁β
 1s₀α 2p₀α
 1s₀α 2p₀β
 1s₀α 2p₁α
 1s₀α 2p₁β
 1s₀β 2p₋₁α
 1s₀β 2p₋₁β
 1s₀β 2p₀α
 1s₀β 2p₀β
 1s₀β 2p₁α
 1s₀β 2p₁β

julia> H = OneBodyHamiltonian() + CoulombInteraction()
ĥ + ĝ

julia> E = Matrix(H, cfgs)
12×12 Array{EnergyExpressions.NBodyMatrixElement,2}:
 (1s₀α|1s₀α) + (2p₋₁α|2p₋₁α) + 0.9999999999999998G¹(1s₀α,2p₋₁α) + 1.7320508075688772F⁰(1s₀α,2p₋₁α)       0                                                               …  0
 0                                                                                                       (1s₀α|1s₀α) + (2p₋₁β|2p₋₁β) + 1.7320508075688772F⁰(1s₀α,2p₋₁β)     0
 (2p₀α|2p₋₁α) + 0.9999999999999998R¹(1s₀α,2p₀α;2p₋₁α,1s₀α) + 1.7320508075688772R⁰(1s₀α,2p₀α;1s₀α,2p₋₁α)  0                                                                  0
 0                                                                                                       (2p₀β|2p₋₁β) + 1.7320508075688772R⁰(1s₀α,2p₀β;1s₀α,2p₋₁β)          0
 (2p₁α|2p₋₁α) + 0.9999999999999998R¹(1s₀α,2p₁α;2p₋₁α,1s₀α) + 1.7320508075688772R⁰(1s₀α,2p₁α;1s₀α,2p₋₁α)  0                                                                  0
 0                                                                                                       (2p₁β|2p₋₁β) + 1.7320508075688772R⁰(1s₀α,2p₁β;1s₀α,2p₋₁β)       …  0
 0                                                                                                       0.9999999999999998R¹(1s₀β,2p₋₁α;2p₋₁β,1s₀α)                        0
 0                                                                                                       0                                                                  (2p₋₁β|2p₁β) + 0.9999999999999998R¹(1s₀β,2p₋₁β;2p₁β,1s₀β) + 1.7320508075688772R⁰(1s₀β,2p₋₁β;1s₀β,2p₁β)
 0                                                                                                       0.9999999999999998R¹(1s₀β,2p₀α;2p₋₁β,1s₀α)                         0
 0                                                                                                       0                                                                  (2p₀β|2p₁β) + 0.9999999999999998R¹(1s₀β,2p₀β;2p₁β,1s₀β) + 1.7320508075688772R⁰(1s₀β,2p₀β;1s₀β,2p₁β)
 0                                                                                                       0.9999999999999998R¹(1s₀β,2p₁α;2p₋₁β,1s₀α)                      …  0
 0                                                                                                       0                                                                  (1s₀β|1s₀β) + (2p₁β|2p₁β) + 0.9999999999999998G¹(1s₀β,2p₁β) + 1.7320508075688772F⁰(1s₀β,2p₁β)

julia> o = cfgs[1].orbitals[1]
1s₀α

julia> co = Conjugate(o)
1s₀α†

julia> diff(E, co)
12×12 SparseArrays.SparseMatrixCSC{EnergyExpressions.LinearCombinationEquation,Int64} with 27 stored entries:
  [1 ,  1]  =  ĥ|1s₀α⟩ + Y¹(2p₋₁α,1s₀α)|2p₋₁α⟩0.9999999999999998 + Y⁰(2p₋₁α,2p₋₁α)|1s₀α⟩1.7320508075688772
  [3 ,  1]  =  Y¹(2p₀α,1s₀α)|2p₋₁α⟩0.9999999999999998 + Y⁰(2p₀α,2p₋₁α)|1s₀α⟩1.7320508075688772
  [5 ,  1]  =  Y¹(2p₁α,1s₀α)|2p₋₁α⟩0.9999999999999998 + Y⁰(2p₁α,2p₋₁α)|1s₀α⟩1.7320508075688772
  [2 ,  2]  =  ĥ|1s₀α⟩ + Y⁰(2p₋₁β,2p₋₁β)|1s₀α⟩1.7320508075688772
  [4 ,  2]  =  Y⁰(2p₀β,2p₋₁β)|1s₀α⟩1.7320508075688772
  [6 ,  2]  =  Y⁰(2p₁β,2p₋₁β)|1s₀α⟩1.7320508075688772
  [1 ,  3]  =  Y¹(2p₋₁α,1s₀α)|2p₀α⟩0.9999999999999998 + Y⁰(2p₋₁α,2p₀α)|1s₀α⟩1.7320508075688772
  [3 ,  3]  =  ĥ|1s₀α⟩ + Y¹(2p₀α,1s₀α)|2p₀α⟩0.9999999999999998 + Y⁰(2p₀α,2p₀α)|1s₀α⟩1.7320508075688772
  [5 ,  3]  =  Y¹(2p₁α,1s₀α)|2p₀α⟩0.9999999999999998 + Y⁰(2p₁α,2p₀α)|1s₀α⟩1.7320508075688772
  [2 ,  4]  =  Y⁰(2p₋₁β,2p₀β)|1s₀α⟩1.7320508075688772
  [4 ,  4]  =  ĥ|1s₀α⟩ + Y⁰(2p₀β,2p₀β)|1s₀α⟩1.7320508075688772
  [6 ,  4]  =  Y⁰(2p₁β,2p₀β)|1s₀α⟩1.7320508075688772
  [1 ,  5]  =  Y¹(2p₋₁α,1s₀α)|2p₁α⟩0.9999999999999998 + Y⁰(2p₋₁α,2p₁α)|1s₀α⟩1.7320508075688772
  [3 ,  5]  =  Y¹(2p₀α,1s₀α)|2p₁α⟩0.9999999999999998 + Y⁰(2p₀α,2p₁α)|1s₀α⟩1.7320508075688772
  [5 ,  5]  =  ĥ|1s₀α⟩ + Y¹(2p₁α,1s₀α)|2p₁α⟩0.9999999999999998 + Y⁰(2p₁α,2p₁α)|1s₀α⟩1.7320508075688772
  [2 ,  6]  =  Y⁰(2p₋₁β,2p₁β)|1s₀α⟩1.7320508075688772
  [4 ,  6]  =  Y⁰(2p₀β,2p₁β)|1s₀α⟩1.7320508075688772
  [6 ,  6]  =  ĥ|1s₀α⟩ + Y⁰(2p₁β,2p₁β)|1s₀α⟩1.7320508075688772
  [2 ,  7]  =  Y¹(2p₋₁β,1s₀β)|2p₋₁α⟩0.9999999999999998
  [4 ,  7]  =  Y¹(2p₀β,1s₀β)|2p₋₁α⟩0.9999999999999998
  [6 ,  7]  =  Y¹(2p₁β,1s₀β)|2p₋₁α⟩0.9999999999999998
  [2 ,  9]  =  Y¹(2p₋₁β,1s₀β)|2p₀α⟩0.9999999999999998
  [4 ,  9]  =  Y¹(2p₀β,1s₀β)|2p₀α⟩0.9999999999999998
  [6 ,  9]  =  Y¹(2p₁β,1s₀β)|2p₀α⟩0.9999999999999998
  [2 , 11]  =  Y¹(2p₋₁β,1s₀β)|2p₁α⟩0.9999999999999998
  [4 , 11]  =  Y¹(2p₀β,1s₀β)|2p₁α⟩0.9999999999999998
  [6 , 11]  =  Y¹(2p₁β,1s₀β)|2p₁α⟩0.9999999999999998

C.f. Eqs. (3.11–12) of

• Froese Fischer, C., Brage, T., & Jönsson, P. (1997). Computational
  Atomic Structure : An Mchf Approach. Bristol, UK Philadelphia, Penn:
  Institute of Physics Publ.