SOFT Method

The goal of this section is to introduce the grid-based Split Operator Fourier Transform (SOFT) method for step-wise propagation of density matrices from an inital high-temperature $T_i=1/(k_B \epsilon)$ to a final lower-temperature $T_f=1/(k_B \beta)$, where $\beta \equiv n \epsilon$. For simplicitly, we write all expressions in mass-weighted coordinates, so all degrees have the same mass m.

Boltzmann-operator matrix-elements $\rho({\bf x},{\bf x}';\beta) \equiv \langle {\bf x} \vert e^{-\beta \hat{H}_0} \vert {\bf x}' \rangle$ are obtained by solving the Bloch equation,

$$\{\frac{\partial}{\partial \beta} -\frac{1}{2 m} \nabla_{\bf x}^2 + V_0({\bf x}) \} \rho({\bf x},{\bf x}';\beta) = 0,$$ (530)

subject to the initial condition given by the high-temperature approximation,

$$\rho({\bf x},{\bf x}';\epsilon) = \left(\frac{m}{2 \pi \epsilon \... ...\bf x})+V_0({\bf x}')]} e^{-\frac{m}{2 \hbar^2 \epsilon} ({\bf x}-{\bf x}')^2},$$ (531)

where $\epsilon$ defines a sufficiently high temperature $T_i=1/(k_B \epsilon)$. Equation (530) is formally integrated as follows,

$$\begin{split} \rho({\bf x},{\bf x}';\beta) = \int d{\bf x}''... ...\beta-\epsilon) \rho({\bf x}'',{\bf x}';\epsilon), \end{split}$$ (532)

where the propagator $\rho({\bf x},{\bf x}'';\beta-\epsilon) \equiv \langle {\bf x} \vert e^{- (\beta-\epsilon) \hat{H}_0} \vert {\bf x}'' \rangle$ is imaginary-time sliced by repeatedly inserting the resolution of identity,

$$\hat{{\bf 1}} = \int d{\bf x}_{j} \vert {\bf x}_j \rangle \langle {\bf x}_j \vert,$$ (533)

yielding,

$$\begin{split} \langle {\bf x} \vert e^{- (\beta-\epsilon) \h... ...i/\hbar) \hat{H}_0 \tau } \mid {\bf x}'' \rangle, \end{split}$$ (534)

where $\tau \equiv - i \hbar (\beta-\epsilon)/s$ is a sufficiently thin imaginary-time slice. Each finite-time propagator, introduced by Eq. (534), is approximated for sufficiently small imaginary-time slices $\tau$ by the Trotter expansion to second-order accuracy,

$$\begin{split} e^{-(i/\hbar) \hspace{.1cm} \hat{H}_0 \tau} \a... ...i/\hbar) \hspace{.1cm} V_0(\hat{{\bf x}}) \tau/2}, \end{split}$$ (535)

and implemented according to a grid-based representation of the initial condition $\rho({\bf x}_j,{\bf x'}_k;{\epsilon})$ as follows:

• Step [1]: Multiply the initial state $\rho({\bf x}_j,{\bf x'}_k;{\epsilon})$ by the operator $e^{-i V_0({\bf x}_j) \tau/2}$ (i.e., the first half of the potential energy part of the propagator):

  $$\tilde{\rho}({\bf x}_j,{\bf x'}_k;{\epsilon}) \equiv e^{-\frac{i}{\hbar} V_0({\bf x}_j) \tau/2} \rho({\bf x}_j,{\bf x'}_k;{\epsilon}).$$ (536)

  
  

• Step [2]: Fourier transform $\tilde{\rho}({\bf x}_j,{\bf x'}_k;{\epsilon})$ with respect to ${\bf x}_j$ and multiply the resulting function by the kinetic energy part of the Trotter expansion,

  $$\tilde{\rho}({\bf p}_j,{\bf x'}_k;{\epsilon}) = e^{- \frac{i}{\hb... ..._j e^{- \frac{i}{\hbar} p_j x_j} \tilde{\rho}({\bf x}_j,{\bf x'}_k;{\epsilon}),$$ (537)

  
  
  where $\Delta x = x_j-x_{j-1}$ is the grid spacing.
• Step [3]: Inverse Fourier transform $\tilde{\rho}({\bf p}_j,{\bf x'}_k;{\epsilon})$ with respect to $p_j$ back to the coordinate representation and multiply the resulting function by the second half of the potential energy part of the propagator, obtaining the imaginary-time evolved Boltzmann-operator matrix elements at $2 \epsilon$ as follows:

  $$\rho({\bf x}_j,{\bf x'}_k;{2 \epsilon}) = e^{- \frac{i}{\hbar} V_... ...um_j e^{\frac{i}{\hbar} p_j x_j} \tilde{\rho}({\bf p}_j,{\bf x'}_k;{\epsilon}),$$ (538)

  
  
  where $\Delta p = p_j-p_{j-1}$.
• Step [3] Go to [1], substituting $\rho({\bf x}_j,{\bf x'}_k;{\epsilon})$ by $\rho({\bf x}_j,{\bf x'}_k;{2 \epsilon})$.

After $n$ iterations of the SOFT algorithm described by steps [1]--[3], the intial state $\rho({\bf x}_j,{\bf x'}_k;{\epsilon})$ is evolved to the final state $\rho({\bf x}_j,{\bf x'}_k;{\beta})$. Note that the SOFT method can also be implemented for propagating wave-packets $\Psi({\bf x};t)$ in real time, exploiting the similarity between the Bloch equation introduced by Eq. (530) and the Schrödinger equation,

$$\{-i \hbar \frac{\partial}{\partial t} -\frac{1}{2 m} \nabla_{\bf x}^2 + V_0({\bf x}) \} \Psi({\bf x};t) = 0.$$ (539)

The fortran program wp.f, which can be downloaded from here, implements the grid-based SOFT method and propagates two wavepacket components on two 1-dimensional poential energy surfaces. As given, the program propagate Gaussian wavepackets on a harmonic potential and tunneling through a double well, although the subroutine ``Hamil'' can be easily modified to model any arbitrary potential including 2 coupled 1-dimensional surfaces.


Exercise(optional)
            Item (A): Generalize the program wp.f to perform 1-dimensional wavepacket propagation of a state $\Psi(x,x';t)$, which depends parametrically on $x'$.
            Item (B): Make the variable substitution $\beta = i t/ \hbar$ and use your program to propagate the density matrix of a particle in a harmonic potential from a high-temperature $T_i$ to a final temperature $T_f$.
            Item (C): Compare the density of states $P(x;\beta)= Z^{-1} \rho(x,x;\beta)$, obtained in (B) at $\beta_i = 1/(k_B T_i)$ and $\beta_f = 1/(k_B T_f)$, to the corresponding analytic expressions given by Eq. (523) at $T_i$ and $T_f$, respectively.
            Item (D): Compare the density of states $P(x;\beta)= Z^{-1} \rho(x,x;\beta)$, obtained in (B) at $\beta_i = 1/(k_B T_i)$ and $\beta_f = 1/(k_B T_f)$, to the corresponding classical expression $P_c(x;\beta) = Z^{-1} exp(-\beta V(x))$.
            Item (E): Repeat items (B)--(D) for the double-well potential and analyze the importance of quantum effects, such as tunneling, at high and low temperature.


It is important to note that a problem requiring $O(l)$ grid points for an accurate propagation of the state in 1-dimension, requires $O(l^N)$ points for the solution of a similar problem in $N$-dimensions. Therefore, the applicability of the grid-based SOFT method is limited to systems with very few degrees of freedom since both the storage and manipulation of multidimensional grids is prohibited for other than very small values of $l$ and $N$. This problem, however, can be partially overcome by using compact coherent state representations as implemented in the MP/SOFT approach [Chen, X.; Wu, Y.; Batista V.S. J. Chem. Phys. 122, 64102 (2005); Wu, Y.; Batista V.S. J. Chem. Phys. 121, 1676 (2004)].