How to do Lowdin symmetric orthonormalisation?

Question

Löwdin symmetric orthonormalisation seems to be a common practice in quantum chemistry. I come from a different background though and have to understand it and possibly implement in a computer code. Here is the definition of the transformation;

$$ |{\phi'}\rangle = S^{-1/2}|\phi\rangle, $$

where $S$ is the overlap matrix between the non-orthogonal basis,

$$ S = \left< \phi|\phi \right> . $$

1. My first question is how to obtain $S^{-1/2}$ analytically?

2. Are there any parallel algorithms available to compute it numerically?

Answer 1

You obtain the diagonalization of $S$ as defined $$ S = U \cdot s \cdot U^{\dagger} $$ (which holds for the case of symmetric $S$, which it is) and apply the power needed $$ S^{-\frac{1}{2}} = U \cdot s^{-\frac{1}{2}} \cdot U^{\dagger} $$

For more background, try the application section of Wikipedia's article on diagonalizable matrices. As far as understanding it goes, I'm afraid not much beyond doing some model calculations on 2D or 3D vectors can be done.

I do not have a strong background in parallel computing, but all operations involved are well-known and your linear algebra framework (BLAS, anyone?) will document how to run any of this in parallel and what storage may be needed for it. Note that, typically, diagonalization is not the most expensive step of a self-consistent field calculation (Hartree-Fock or DFT) anyway, so parallelization may not be critical for this step.

Answer 2

Your first question has already been answered, but in words: to find $f(A)$ for some matrix $A$, you diagonalize it to obtain the eigenvalues $a$ and eigenvectors $U$, apply $f$ to the diagonalized matrix (the eigenvalues), then back-transform $f(a)$ using the eigenvectors to the original non-diagonal basis.

For the second question, yes. Most revolve around calling the LAPACK routine dsyev (real double-precision symmetric eigen decomposition) or its variants, ssyev for real single-precision, {c,z}heev for complex single- and double-precision Hermitian decomposition. The vast majority of quantum chemistry uses entirely real coefficients and double-precision.

Assume I have already compute the overlap matrix $S$ by some method, usually by calling a program's integral engine. Here is a sample implementation in Python using NumPy:

print("Overlap matrix")
print(S)

lam_s, l_s = np.linalg.eigh(S)
lam_s = lam_s * np.eye(len(lam_s))
lam_sqrt_inv = np.sqrt(np.linalg.inv(lam_s))
symm_orthog = np.dot(l_s, np.dot(lam_sqrt_inv, l_s.T))

print("Symmetric orthogonalization matrix")
print(symm_orthog)

From the documentation of numpy.linalg.eigh:

The eigenvalues/eigenvectors are computed using LAPACK routines _syevd, _heevd

Here is a sample implementation in C++ using Armadillo:

S.print("Overlap matrix");

arma::vec lam_s_vec;
arma::mat l_s;
arma::eig_sym(lam_s_vec, l_s, S);
arma::mat lam_s_mat = arma::diagmat(lam_s_vec);
arma::mat lam_sqrt_inv = arma::sqrt(arma::inv(lam_s_mat));
arma::mat symm_orthog = l_s * lam_sqrt_inv * l_s.t();

symm_orthog.print("Symmetric orthogonalization matrix");

Both working examples can be found here; check the Makefile for how to run.

In the Armadillo source tree, the files include/armadillo_bits/{def,wrapper}_lapack.hpp contain more information about which LAPACK routines are called for which types.

Regarding parallelization, your BLAS + LAPACK implementation (MKL, OpenBLAS, ATLAS, ...) is most likely threaded and can be controlled by <something>_NUM_THREADS=4, where <something> might be OMP, MKL, OpenBLAS, or possibly something else, but you should check the documentation. This means that as long as your environment is set up properly, math library calls with NumPy or C++ template libraries like Armadillo or Eigen will run in parallel without explicit OpenMP annotations or MPI code. For distributed parallelization (MPI), there is ScaLAPACK, which shares a similar interface to regular LAPACK.