imaginary time propagation to find ground state wavefunction

Question

I understand the basic idea of imaginary time propagation method:

The wavefunction $\psi(x,t)$ as a superposition of energy eigenstates $\phi_m(x)$: $$ \psi(x,t)=\sum_m \phi_m(x)e^{-iE_mt/\hbar} $$ In imaginary time, ${t}\rightarrow{-it}$, $$ \psi(x,-it)=\sum_m \phi_m(x)e^{-E_mt/\hbar}=e^{-E_ot/\hbar}\bigg(\phi_o+\phi_1e^{-(E_1-E_o)t/\hbar}+\phi_1e^{-(E_2-E_o)t/\hbar}+........\bigg) $$ As time $t$ increases the terms with greater exponents will decay more rapidly, leaving behind the ground state.

But how do i make a program or subroutine to find the ground state wavefunction using imaginary time propagation for a basic one dimensional problem, preferably in Fortran ?

Ref to 7.12 in an introduction to computational physics by Pang.

The macroscopic state of Bose–Einstein condensate is described by the Gross–Pitaevskii equation: $$ i\hbar\frac{\partial \Psi(x,t)}{\partial t}=\frac{-\hbar^{2}}{2m}\frac{\partial^2 \Psi(x,t)}{\partial x^2}+V_{ext}\Psi(x,t)+g|\Psi(x,t)|^2\Psi(x,t) $$ where $\psi(x,t)$ is the macroscopic wavefunction, $m$ and $g$ are two system-related parameters, and $V_{ext}(x)$ is the external potential. Develop a program that solves this equation numerically. Consider various potentials, such as a parabolic or a square well.

Can I apply Finite-difference method or Crank-Nickelson method, similar to simple diffusion equation, by substituting for $\frac{\partial \Psi(x,t)}{\partial t}, \frac{\partial^2 \Psi(x,t)}{\partial x^2}, |\Psi(x,t)|^2\Psi(x,t), \Psi(x,t)$, after inputting the imaginary time to get rid of '$i$' ?

or how do I apply imaginary time evolution to Schrodinger equation of very basic 1D system ?

Note: It'd be very helpful if someone can give me advice on where to start.

Answer 1

It's a bit heavy, but not all of it must be read to get a useable method : Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow, Bao and Du.

I believe -- it's been a while since I worked on it -- that you have to move to imaginary time (as you said, by using $-it$ instead of $t$). Then you want to propagate along time the equation obtained (something like : $\frac{\partial \Psi}{\partial t} = \frac{1}{2} \Delta \Psi - V(x,y)\Psi - \beta |\Psi|^2 \Psi$ ), which can be done as follows.

Each time-step of this propagation can then be separated in two "sub-steps" (basically, it means to take into account two different parts of the righ hand side of the above equation), using "Time splitting spectral method" (TSSP) described p. 10 in the above reference.

You can then solve these "sub-timesteps" on $|\Psi|^2 = \rho$. One is easier to solve in Fourier space, the other is exactly solvable.

For each iteration, the energy can then be computed and you just need to set a convergence criterion to get your ground state.