Comity of errors

In theoretical chemistry, some amount of interest is placed in the distinction between adiabatic and diabatic electronic states. Adiabatic states are well-defined as eigenstates of the electronic Hamiltonian, and are what we usually calcluate. However, sometimes they can be inconvenient to use, especially when we consider a large range of geometries. For example, if I pull apart the two atoms in a molecule of lithium fluoride, at short separations, the bond is of ionic character, but at large separations, the state behaves covalently. The concept of diabatic electronic states is introduced to eliminate this inconsistency — in LiF, we would consider two diabatic states, one of which is ionic in character and the other covalent. These potential energy surfaces for these states will cross at some point, which is where the nature of the adiabatic surface changes.

Unfortunately, there is no rigorous definition of diabatic states — the only possible candidate is not achievable in practice, so we frequently find ourselves considering approximately diabatic states.

We can convert a set of diabatic states into a set of adiabatic states by calculating the derivative couplings between them (e.g. ), and using those to construct a matrix Hamiltonian that we diagonalize; its eigenstates are the adiabatic states (or adiabats).

In our research group, we make extensive use of constrained density functional theory, which uses constraints on the molecular charge and/or spin in particular regions of space to find approximately diabatic electronic states. Now, if I wanted to compute a transition dipole moment (which is a physical observable) between two such states, I would be out of luck. Two general constrained states will come from systems with different Hamiltonians (due to the constraint potential), and as such will be non-orthogonal. Since observables are only defined for a set of orthogonal states, it does not make any sense to directly compute a transition dipole moment between two (different) constrained states. However, there is a simple way to make a general collection of constrained states into a collection of orthogonal constrained states — construct an overlap matrix, and use it to orthogonalize the states in question. For our purposes, symmetric orthogonalization makes the most sense, so we end up with coefficient matrices C’=S^{-1/2}C . Now these are some states we can get some transition dipoles from. And we probably will, in the future.

In a previous post, I mentioned a theorem that would speed up computing the overlap between unrelated molecular wavefunctions. I’ve been putting more work into this project, off-and-on, and come up with a bunch of expressions that I’ll want to evaluate; things that look like T_{ij}^{ab}C_{ia}C_{jc}C_{kb}w_{kc}, where T_{ij}^{ab} is the two-electron integral involving orbitals i, j, a, and b; w_{bc} is the matrix element of a single-electron operator (just a multiplicative potential in real-space); and C_{fg} are these overlaps between wavefunctions. Naively applying the Einstein summation convention, this would look to scale as N^6, but we can pre-sum over i and a (a N^4 step), and are then left with sums over j, k, b, and c, which is also only N^4.

But where do these terms come from? w_{kc} is in practice just a vector-matrix-vector product c_k.W.c_c; computing each term is N^2 time, and there are only N^2 of them, which is N^4 work total. So far so good. The two-electron integral T_{ij}^{ab} is just an integral over two variables r1 and r2 of 1/(r1-r2) times the orbitals i, j, a, and b (in real space), where i and a are functions of r1, and j and b are functions of r2. (Incidentally, the code to compute these integrals is the hard part of any atomic-orbital-based computational chemistry code; I am building off of someone else’s base.) Since the basis functions are all gaussians, these integrals can be evaluated in constant time … but only in the atomic orbital basis. We are working in the molecular orbital basis, and the best transformation of these integrals into the molecular orbital basis (each molecular orbital being a weighted sum of the atomic orbitals) requires N^5 time. However, this code is the bottleneck for many things, and is already heavily optimized, so I’ll go ahead and worry about other things.
The C_{ia} are nominally overlaps between two determinantal wavefunctions, with each determinant being M-factorial terms (M being the number of electrons); computing the overlap between a pair of terms is constant time (just lookup the one-electron orbital overlaps), but there are still (M!)^2 terms! M is frequently an order of magnitude less than N, but this is still quite a bit of work. It turns out, though (and I’m not planning to do the derivation here), that if we write things out in terms of permutations, then this collapses down into the determinant of a single M-by-M matrix, the condensation of C_{phi}.S.C_{psi}, where C_{phi} is an M-by-N matrix of orbital coefficients, S is the overlap matrix in the atomic orbital basis, and C_{psi} is the N-by-M matrix of orbital coefficients for psi. The matrix-matrix-matrix product is a (large-constant-factor) N^3 operation, as is the determinant, but there are N^2 of these terms to compute, so the C_{ia} as a whole contribute an N^5 factor. Unfortunately, this is all in code that I am writing, so it’s not terribly optimized.
Now, this whole setup really makes one think that it might be further optimized — S is always the same, and C_{psi} is always the same. Furthermore, the C_{phi} are constructed from a common base, with a single column being changed from that common base each time. (Which column that is can be different, though.) If M was equal to N, then we could use the associativity of the determinant to precompute det(S) and det(C_{psi}), and all of the appropriate minors from C_{phi}, and then do a O(N) expansion over minors for each C_{ia} term. However, the determinant is not associative for non-square matrices, so this speedup doesn’t work. Any linear algebra masters reading are requested to comment on whether such a speedup may or may not be possible in the general case.
For now, I’ll take some relish in the fact that this is not the (only) scaling bottleneck ….