I posted a similar question yesterday but I deleted it since I think that I had to reformulate it after some insights.

I want to calculate $$ \exp(-i\Delta t\,\mathcal{H}) = V\,\mathrm{diag}(\{\exp(-i\Delta t\,\epsilon_k\})\,V^{-1} ~~~~~~~~~~~~~(*) $$ where $\mathcal{H}$ is a hermitian matrix (the hamilton matrix) and $V$, $\epsilon_k$ are its eigenvectors and -values.

Since I'm implementing this in a given DFT package, I make use of the existing eigensolver routine. All numerics are done by ScaLAPACK. It first computes the Cholesky factorization of the overlap matrix $S_{ij}=\langle\phi_i|\phi_j\rangle$: $$ S = U^TU $$ where $U$ is the upper triangular matrix. In the next step, $U$ is inverted for further usage. It then solves the eigenvalue problem for the transformed hamilton matrix $$ \mathcal{H}' = U^{-\dagger}\mathcal{H}\,U^{-1} $$ which yields the eigenvectors $V'$. This is done since we have a generalized eigenvalue problem because a non-orthogonal basis set is used. The eigenvectors are afterwards backtransformed via $$ V = U^{-1} V'. $$

I tried to calculate $(*)$ by inverting $V$ which needs a LU decomposition since $V$ is a non-hermitian matrix. This turned out to be unstable after some iterations although it initially produced correct results.

I found a different way to circumvent inverting $V$ which uses the not backtransformed eigenvectors $V'$. Here, the exponential is explicitly calculated via $$ (U^{-1}V')\,\mathrm{diag}(\{\exp(-i\Delta t\,\epsilon_k\})(U^{-1}V')^{-1}=U^{-1}V'\,\mathrm{diag}(\{\exp(-i\Delta t\,\epsilon_k\})V'^\dagger U $$ where $V^{-1}=V'^{\dagger}$ can be applied since $V'$ is unitary. Unfortunately, this yields oscillating solutions around the right solution (e.g. total energy).

I'm not sure why this isn't working and I'm quite new to the field. I tried $S^{-1}$ and $S$ in above formula instead of $U^{-1}$ and $U$ but it didn't help, also I tried a lot more modifications but I just can't figure out where the problem could be.

I would really appreciate some help. Thanks and best regards!


(Heavily edited because I made a big mistake in the original post)

We write the Time Dependent Kohn-Sham equations as:

$Hc = i\hbar S\frac{\delta c}{\delta t} $

where $c$ are the coefficients of your wavefunction. This is also the same as:

$S^{-1}Hc = i\hbar\frac{\delta c}{\delta t} $

and so we can update the wavefunctions as:

$c(t+\Delta t) = exp(-i\hbar S^{-1}H\Delta t)c(t)$.

So now we want to compute $exp(-i\hbar S^{-1}H\Delta t)$, through the eigendecomposition as you said.

$S^{-1}H = X^{-1}DX$

Let's write that in terms of a familiar problem:

$S^{-1}HX = DX$

$HX = DSX$

Using the cholesky decomposition, as you suggested, we can turn this into the standard eigenvalue problem.

$S = LL^{T}$

$(L^{-1}HL^{-T})(L^{T}X) = D(L^{T}X)$

$L^T X = X'$

$L^{-1}HL^{-T} = H'$

and then just solve

$H'X' = DX'$.

Now let's work our way back.

$X = L^{-T}X'$

$X^{-1} = X'^{-1}L^{T} = L^{T}X'^{T}$

(since $X'$ is an orthogonal matrix).

So now we know D, X, and $X^{-1}$, and thus can compute:

$exp(-i\hbar S^{-1}H\Delta t) = exp(-i\hbar X^{-1}DX\Delta t) = X^{-1}exp(-i\hbar D\Delta t)X$.

Let's make sure I'm correct with some source code:

import scipy.linalg
import numpy.random
import numpy

# Setup
N = 10
H = 0.01*numpy.random.rand(N,N)
H = H+H.T # symmetric
S = numpy.random.rand(N,N)
S = S+S.T # symmetric
S = S.dot(S) # positive definite
Sinv = scipy.linalg.inv(S)
I = numpy.identity(N)

# Reference answer
ref = scipy.linalg.expm(Sinv.dot(H))

# Our solution
## first our cholesky factors
L = scipy.linalg.cholesky(S,lower=True)
Linv = scipy.linalg.solve_triangular(L, I, lower=True)

## now solve the normal eigenvalue problem
Hprime = Linv.dot(H).dot(Linv.T)
D, Xprime = scipy.linalg.eigh(Hprime)
D = numpy.diag(D)

## back subsite for X and X^{-1}
X = Linv.T.dot(Xprime)
Xinv = Xprime.T.dot(L.T)

our_answer = X.dot(scipy.linalg.expm(D)).dot(Xinv)

print(scipy.linalg.norm(ref - our_answer))
| cite | improve this answer | |
  • $\begingroup$ Thanks for your remark. So, if I understand you right, with this scheme, I need to solve the eigenvalue problem for $S^{-1}H$ to get the respective eigenvectors $V$. I wanted to avoid a modification of the eigensolver but this seems straightforward. So the transformed hamiltonian $U^{-\dagger}HU^{-1}:=H'$ is again modified to $S^{-1}H'$ before the EVP is solved or is this not necessary anymore because we handled the overlap already? $\endgroup$ – Lukk Sep 28 '17 at 8:46
  • $\begingroup$ Thanks a lot, KobeGote, I got it working now. I really appreciate it! $\endgroup$ – Lukk Sep 29 '17 at 8:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.