Handling of complex arithmetic

[jjgoings]

It would be good to toggle on or off real vs complex arithmetic.

Right now we (for the most part) treat everything as complex-valued, because the goal at the time was to allow for real-time electronic dynamics. For most users this is just additional computational overhead, so it would be good to implement the ability to switch between real and complex, and include safeguards to force complex if doing RT-TDSCF.

[pwborthwick]

Hi Josh, I used the other approach using real for everything then switching to complex for the RT-TDSCF (still a work in progress!). The only place I had to coerse to real was returning from the diis extrapolation routine. (and in CI where eigh returns complex values). This seems OK as both our initial F and P are real on starting the RT_TDSCF. It shouldn't be too difficult to implement a switch and going real might reduce time with big molecules/basis sets.

[jjgoings]

@pwborthwick That makes sense. It would be best to toggle real/complex for the overall calculation, because there are some cases we want everything complex. Singlet O2 is a great example where the lowest energy Slater determinant is complex RHF. And, in the future, should we implement variational relativistic methods like X2C, complex is necessary as well. I am surprised that eigh returns complex. It should always be real for eigenvalues, and I would think that if the input matrix is real and Hermitian, the eigensolver would return real eigenvectors as well.

[pwborthwick]

I've just checked and its 'eigvals' in the random phase approximation that despite having a real Hamiltonian returns (a few) complex eigenvalues. I'm going to re-write my code now to switch between real/complex. I agree that its a pity about the angular momentum commutators as they would provide a self-contained test for the 'E' auxilliary function.

[jjgoings]

@pwborthwick Unless you do some tricks to reduce the dimension, the full RPA response matrix is non-Hermitian, and can in principle give complex eigenvalues. See, for example, Eq(31) here. So it's not necessarily something to worry about, although RPA will give imaginary eigenvalues if the underlying SCF wave function is unstable. For a stable wave function all the eigenvalues will be real.

[pwborthwick]

Yes the RPA Hamiltonian is not symmetric I know this because I had to write a QR algorithm in a BASIC Hartree-Fock program I wrote because the Jacobi doesn't do non-symmetric. The hamiltonian is water in STO-3G, it's totally real but 'eigvals' return a few eigenvalues with a complex part ~1e-15i. So the complex is in the noise and not anything to worry about as you say. I've decided not to write the main SCF in complex rather than real, I will look at O2 singlet later and maybe try to incorporate some sort of stability check. I have UHF coded. I had, of course, read your article. Thanks.

[AlexB67]

It would be good to toggle on or off real vs complex arithmetic.

Right now we (for the most part) treat everything as complex-valued, because the goal at the time was to allow for real-time electronic dynamics. For most users this is just additional computational overhead, so it would be good to implement the ability to switch between real and complex, and include safeguards to force complex if doing RT-TDSCF.

Hi I learned a lot from your code getting back into theochem. I think there is a lot of optimization that could be done before worrying about complex arithmetic. I think there is an order of magnitude in there for the ints .. easily. even changing the loop structure for the ERI ints already gives me a 10 - 20% increase for h2o cc-pvdz. around 2.3 secs with old code down to 1.9 secs with a small change

for t in range(l1+l2+1): l12 = E(l1,l2,t,A[0]-B[0],a,b) for u in range(m1+m2+1): m12 = E(m1,m2,u,A[1]-B[1],a,b) for v in range(n1+n2+1): n12 = E(n1,n2,v,A[2]-B[2],a,b) for tau in range(l3+l4+1): l34 = E(l3,l4,tau,C[0]-D[0],c,d) for nu in range(m3+m4+1): m34 = E(m3,m4,nu ,C[1]-D[1],c,d) for phi in range(n3+n4+1): val += l12 m12 n12 l34 m34 \ E(n3,n4,phi,C[2]-D[2],c,d) \ pow(-1,tau+nu+phi) * \ R(t+tau,u+nu,v+phi,0,\ alpha,Px-Qx,Py-Qy,Pz-Qz,RPQ)

but more importantly when computing the polynomial terms E(l1, l2) term and R(t+ tau etc. They currently are recomputed inside a N^6. loop. The E terms can be precomputed as arrays with N^2 operations and the Hermite coulomb R(.. terms with N^4, This will give you an almighty boost. Also, the gradient ints can gain a 25% boost using translational invariance, since centers A + B +-C = -D so you compute 3 and get the fourth for free. similar relationships exist for the 1e ints. Where you calculate A and B, you only need B or A. Not aware if complex arithmatic has any bearing on this mind you.

I have C++ versions fort this if you want, but I hope you get the gist. I've attached some stuff for the future perhaps. I would add a shell based structure with the OS scheme would be must faster still which is what I ended up doing.

Cheers. Your blogs were great, not having access to a UNI library they saved my day early on. ints.zip

[jjgoings]

Hey thanks @AlexB67!

These comments and suggestions are great. You are right, with regards to complex, all the integrals are done in real arithmetic. The complex comes afterward when we build Fock, do SCF, etc. So complex arithmetic and the integrals are really two separate issues.

For the integrals, I am split on re-structuring. On the one hand, removing the explicit recursion and doing (for example) pre-computed arrays like you mention would be a smart change for efficiency. On the other hand, I like keeping the structure for pedagogical purposes. This project has been, and I think will continue to be, an experiment for a readable end-to-end implementation of integrals for quantum chemistry -- hence the clunky name "McMurchie-Davidson". I never intended it to be a standalone electronic structure package, but hey, here we are :) I have often toyed with the idea of interfacing an external integral library like LibInt, but I think that defeats the educational merit of the project. If efficiency is desired, PySCF or Psi4 (+NumPy) are probably the way to go, and not MMD.

Kind of thinking out loud here. If I find time (unlikely), I might try some of these changes. But I'm happy to take contributions via pull requests! I think a few modest changes to the integral library are OK if they can strike a balance between readability and efficiency.

[AlexB67]

I agree, for ease of reading it very easy to follow. I never got to a chance to implement them in your code but it was of great help to me to make those changes in C++. with what I mentioned h2o cc-pVQZ went from a minutes regime, to a seconds regime. however as my attached files show mine ended up not nearly as easy to follow ... C++ aside

In any case, just to throw the ideas out there. I can give it a crack for you some day. The translational invariance tricks would be very easy and not detract from readability but only to shorten the code.

as an example (shameless plug) See how I approached it here.

https://github.com/AlexB67/hfscf/blob/master/src/gradient/hfscf_grad_ints.cpp

Cheers.

[jjgoings]

Excellent @AlexB67! Impressive timing, though I'm not surprised. I'll gladly take any help I can get :)