Boundstate and Continuum WF about fac HOT 2 OPEN

all in atomic units

from fac.

Sorry to disturb old thread. About bound state orthogonality. Is the ortohgonality condiction:
∫(P_{n1κ1}(r)*P_{n2κ2}(r)+Q_{n1κ1}(r)Q_{n2κ2}(r))dr=δ_{n1n2}*δ_{κ1κ2}
correct?(Sorry bad at LATEX) If it is, then the following results confuses me.
I tried to calculate ∫(P_{n1=3,κ1=1}(r)*P_{n2=3,κ2=-2}(r))dr of Si1+:1s2_2s2_2p5_3s2_3p2. If the above orthogonality condition is correct, then this should be approximately 0? Since Q functions are generally much smaller than P functions? Here's my code:

from pfac import fac
fac.SetAtom('Si')
fac.Config('g','1s2 2s2 2p5 3s2 3p2')
fac.ConfigEnergy(0)
fac.OptimizeRadial('g')
fac.ConfigEnergy(1)
fac.Structure('Si1+.lev.b', ['g'])
fac.MemENTable('Si1+.lev.b')
fac.PrintTable('Si1+.lev.b','Si1+.lev',1)
fac.WaveFuncTable("Si1+(3p-).wf",3,1)
fac.WaveFuncTable("Si1+(3p+).wf",3,2)

According to another issue, the data I need is the 5th column. Then the integral returns 0.999486, which is too much larger than "approximately 0", but it is very close to 1. I guess the orthogonality condition above may be wrong, may I know what is the correct orthogonality condition? Or does there even exist such a orthogonality condition for bound states that looks like what Guadaloupe mentioned? Or is the program somehow acting faulty?
Sorry for the disturbance.

from fac.