physical chemistry - How to do Lowdin symmetric orthonormalisation?

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

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 _NUM_THREADS=4, where 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.