psi4
ADC Eigenvectors
Hi, I've checked the ADC eigenvalues against pyscf and there is agreement. However when looking at the eigenvectors there seems to be a problem. It looks like the Davidson routine is returning the full expanded sub-space rather than the eigenvectors that correspond to the returned eigenvalues. For IP the returned vector v_ip[:,0] is just a column of zeroes. Could someone check this please. Peter
Hi Peter, thanks for spotting this - I have addressed it in #120. You were correct that too many vectors were being returned from davidson - though v_ip[:, 0] should still have been correct in the old code, there would just be too many columns in v_ip. In the IP example v_ip[8, 0] has non-zero weight, which corresponds to the 8th (0-based) spin-orbital being strongly coupled to the IP - the rest of the vector being mostly zeros indicates that most of the character of this excitation is in that particular orbital.
I have also added printing of the quasiparticle weight, which may be useful to you - this is equal to the "spectral factor" in PySCF for unrestricted ADC, and equal to half the "spectral factor" for restricted ADC.
Hi @obackhouse , Thank you very much for you speedy reply. I should say that I don't have psi4/psi4numpy loaded on my computer so am providing 3rd-party integrals to the psi4numpy code. For water 6-31g IP, I get the following energies
IP (Ha) IP (eV)
0 0.37045812 10.08068273
1 0.37045812 10.08068273
2 0.45411512 12.35710670
3 0.45411512 12.35710670
4 0.60030477 16.33513128
and from pySCF
ADC calculation summary
adc(2) root 0 | Energy (Eh) = 0.3704581157 | Energy (eV) = 10.08068397 | Spec factors = 1.77529061
adc(2) root 1 | Energy (Eh) = 0.4541151172 | Energy (eV) = 12.35710810 | Spec factors = 1.80479914
adc(2) root 2 | Energy (Eh) = 0.6003047718 | Energy (eV) = 16.33513327 | Spec factors = 1.85441321
So the excitation energies are in excellent agreement (the difference in eV values just down to precision of conversion factors used). However, the spectroscopic factors calculated via psi4numpy (I've implemented your changes) are
1.7903528
1.7903528
1.8276638
1.8276638
1.8845276
and my v_ip[:10,0] is
[ 0.00000000e+00 1.05820404e-18 1.47466953e-16 -2.35402293e-17
3.60822483e-16 2.63178130e-13 -3.06005221e-15 -6.07153217e-15
-2.72004641e-15 9.46137631e-01]
This value for spin orbital 9 corresponds to the pySCF value
1h block:
i U(i)
------------------
5 0.9461
but v_ip[:10,2] is
[ 3.29288537e-04 -1.69256154e-04 -4.87733530e-02 2.50697768e-02
3.08000153e-15 3.18646708e-16 8.48807362e-01 -4.36291741e-01
-6.30435633e-14 -1.02999206e-18]
The largest value is spin orbital 6 but pySCF has
1h block:
i U(i)
------------------
4 0.9544
I appreciate that I'm comparing a spin treatment with a spatial one but this looks wrong. If you get different values running genuine psi4numpy code then I guess I've done something wrong in my implementation. Do the vectors you're getting match psi4? Once again I appreciate you looking at this. Peter PS I see you're at KCL, is that Kings' College London?
I have tried different combinations of integrals and Davidson solver from pyscf and psi4 with the psi4numpy ADC equations, and all are in good agreement with the eigenvectors from pyscf's ADC implementation:
pyscf integrals, psi4numpy equations, psi4numpy solver:
0.45411530211638756
[3.7046e-04 5.4884e-02 3.3517e-06 9.5437e-01 1.0092e-14]
pyscf integrals, psi4numpy equations, pyscf solver:
0.45411530211640455
[3.7047e-04 5.4884e-02 3.3517e-06 9.5437e-01 9.3286e-15]
psi4 integrals, psi4numpy equations, psi4numpy solver:
0.45411541019518575
[3.7046e-04 5.4884e-02 3.3517e-06 9.5437e-01 5.9368e-13]
psi4 integrals, psi4numpy equations, pyscf solver:
-0.4541154101951957
[3.7046e-04 5.4884e-02 3.3517e-06 9.5437e-01 3.6885e-14]
pyscf.adc.radc.RADC:
0.4541153021163934
[-3.7046e-04 5.4884e-02 -3.3517e-06 -9.5437e-01 -4.6132e-16]
Note that I have also summed over spin to remove the degeneracy due to the spin-orbital formulation. So I am satisfied with the values printed by psi4numpy for this system. Can you check if you are mixing up the eigenvector indices in your implementation - I define the vectors as (ADC space, excitation space).
PS I see you're at KCL, is that Kings' College London?
I am indeed :smile:
Hi, I'll have to check my spin-summation it's probably why I get incorrect results. I'm happy from your results that all is OK with psi4numpy now. I was an undergraduate(maths)/graduate(theoretical physics) at Kings' in the late 60's early '70's, we were the first to go into the Strand building from Surrey Street. Thanks for all your efforts. Peter
I'll have to check my spin-summation it's probably why I get incorrect results. I'm happy from your results that all is OK with psi4numpy now.
@pwborthwick no worries - let me know if I can help any more.
I was an undergraduate(maths)/graduate(theoretical physics) at Kings' in the late 60's early '70's, we were the first to go into the Strand building from Surrey Street.
Ah wow, small world! We're still going strong on the Strand.
Hi @obackhouse, Can I just conform that if I do an IP calculation in 6-31g that for the eigenvalue 0.4541153021163 the eigenvectors are v_ip[:,2] and v_ip[:,3] and they split into [nocc] (10) and [nocc, nocc, none] (10x10x16=1600) blocks and you have spin-summed over the first block to get the results above? If I've got this right I can figure out the rest. Peter
@pwborthwick yes, v_ip has shape (nocc + nocc*nocc*nvir, nroots). I projected that into the singles space for a particular root i, i.e. v_ip[:nocc, i], and then took the norms over the degenerate pairs i.e. sqrt(v_ip[:nocc, i][::2]**2 + v_ip[:nocc, i][1::2]**2). For a restricted reference, that should be the same for i=2 and i=3 which have the eigenvalue 0.454115(...). Note that this isn't as simple for the doubles (nocc*nocc*nvir) part.
I've uploaded the script I used to get the above results here - should run fine if you comment out the psi4 parts.
Hi, I've tried your code and it's fine for 6-31g, but out of interest I tried sto-3g (given that on my computer it's the only basis small enough to run EE computations) and it doesn't agree with pyscf.
pyscf integrals, psi4numpy equations, psi4numpy solver:
-0.3777435472744927
[0. 0. 0. 1. 0.]
pyscf integrals, psi4numpy equations, pyscf solver:
-0.3777435472745175
[2.4511e-04 4.8115e-02 2.7979e-06 9.6912e-01 3.0290e-15]
psi4numpy solver is converged OK. As this a almost certainly a problem with the Davidson solver and the point of the code is really to demonstrate ADC you might want to ignore this? Peter
Yes, this look like the Davidson solver in psi4numpy is converging such that the character of the ionisation potential is completely due to that MO and no other excitation - but still getting the excitation correct... the Davidson code I wrote for psi4numpy is very simple, the solver in pyscf is quite a bit more robust. I definitely won't ignore this though because it's a really interesting example, I'll take a closer look at it.
... (given that on my computer it's the only basis small enough to run EE computations) ...
Maybe try again after merge of #121 - I just realised that the initialisation of guess vectors in the psi4numpy ADC codes were scaling horrifically with system size. My bad!
Hi @obackhouse, Yes that's the same guess as pyscf and much better. As for the solver try...
de = np.linalg.norm(theta[:k] - theta_old)
if de < tol:
conv = True
b = np.dot(b, alpha)
break
else:
if b.shape[1] >= (k * nvec_per_root):
b = np.dot(b, alpha)
theta = theta_old
else:
b = np.concatenate([b, np.column_stack(b_new)], axis=1)
b = b[:, :k]
return theta, b
Peter
@pwborthwick good catch! That is another bug. Seems like I wrote a really buggy solver...
No it's a nice compact solver that does the job really well. Thanks for this ADC code I was ploughing through Hodecker's thesis on ADC and trying my own implementation and your psi4numpy presentation was an excellent accompaniment to the theory. Peter