This article discusses using the Time-Dependent Variational Principle (TDVP) to solve the quantum harmonic oscillator. It provides insights into TDVP's application in quantum chemistry, based on a lab seminar homework problem.
Loading content…
Note: This article presents a solution to a problem assigned as homework in our lab's seminar, with the permission of the seminar organizer.
When solving for electronic states in quantum chemistry calculations, variational principles such as:
are often employed. This article provides a concrete example of calculations using TDVP.
The time-dependent Schrödinger equation we will solve is:
Here, the Hamiltonian is defined as:
We express the wave function using time-dependent parameters
This Hamiltonian represents a parabolic potential centered at
The initial conditions are:
In TDVP, we calculate matrix A and C from the partial derivatives of the wave function with respect to its parameters. By utilizing the following relationship, we update the parameters:
where R and I represent the real and imaginary parts of the elements, respectively.
A and C are defined as:
Since we have two parameters in this case, the explicit form is:
Calculating the A and C matrices, we obtain:
Therefore, we have the following system of differential equations:
Solving these equations with the initial conditions given in
Substituting these into
Let's verify that the solution obtained using TDVP,
We add a global phase factor to the wave function:
Substituting this into both sides of
and comparing the coefficients, we obtain the following conditions:
Since they holds for any
By substituting the solutions for
Integrating this equation with the initial condition
Therefore, excluding the global phase factor, the solution
The time evolution of the wave function
The square of the probability amplitude
As can be seen, the probability wave oscillates around
The visualization was generated using the following program.
Note: The plots show the analytical solution, not the result by TDVP's time evolution.
# animation of real and imaginary parts
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.collections import LineCollection
x0 = 0
p0 = 2
xs = np.linspace(-5, 5, 100)
def wf(x, t):
theta_1 = x0 * np.cos(t) + p0 * np.sin(t)
theta_2 = -x0 * np.sin(t) + p0 * np.cos(t)
gamma = ((x0**2+p0**2-1)*t-x0*p0*np.sin(2*t))/2
ret = np.pi**(-1/4) * np.exp(-1/2 * (x - theta_1)**2 + 1j*theta_2*x+1j*gamma)
return ret
fig, ax = plt.subplots()
ax.set_xlim(-5, 5)
ax.set_ylim(-1, 1)
ln_r, = ax.plot([], [], 'r-', label='real')
ln_i, = ax.plot([], [], 'b-', label='imag')
def init():
ln_r.set_data([], [])
ln_i.set_data([], [])
return ln_r, ln_i
def update(frame):
data = [wf(x, frame) for x in xs]
vals_r= [d.real for d in data]
vals_i= [d.imag for d in data]
ln_r.set_data(xs, vals_r)
ln_i.set_data(xs, vals_i)
return ln_r, ln_i
ani = FuncAnimation(
fig, update, frames=np.linspace(0, np.pi*2, 100, endpoint=False),
init_func=init, blit=True, repeat=True
)
ax.legend()
ani.save("wavefunc.gif", writer="pillow", fps=20)
# animation of probability
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.collections import LineCollection
x0 = 0
p0 = 2
xs = np.linspace(-5, 5, 100)
def wf(x, t):
theta_1 = x0 * np.cos(t) + p0 * np.sin(t)
theta_2 = -x0 * np.sin(t) + p0 * np.cos(t)
gamma = ((x0**2+p0**2-1)*t-x0*p0*np.sin(2*t))/2
ret = np.pi**(-1/4) * np.exp(-1/2 * (x - theta_1)**2 + 1j*theta_2*x+1j*gamma)
return ret
fig, ax = plt.subplots()
ax.set_xlim(-5, 5)
ax.set_ylim(0, 1)
ln, = ax.plot([], [], 'r-')
lc = LineCollection([], cmap='twilight', linewidth=2)
ax.add_collection(lc)
norm = plt.Normalize(0, 2 * np.pi)
sm = plt.cm.ScalarMappable(cmap='twilight', norm=norm)
sm.set_array([])
fig.colorbar(sm, ax=ax, orientation='vertical', label="Phase (rad)")
def init():
ln.set_data([], [])
lc.set_segments([])
return ln, lc
def update(frame):
data = [wf(x, frame) for x in xs]
amps = [np.abs(d) for d in data]
phases = [np.angle(d) for d in data]
ln.set_data(xs, amps)
points = np.array([xs, amps]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
lc.set_segments(segments)
lc.set_array(phases)
return ln, lc
ani = FuncAnimation(
fig, update, frames=np.linspace(0, np.pi*2, 100, endpoint=False),
init_func=init, blit=True, repeat=True
)
ani.save("prob.gif", writer="pillow", fps=20)