import numpy as np
from scipy.linalg import eigh
from scipy.optimize import minimize_scalar
from scipy.integrate import quad
import matplotlib.pyplot as plt
def compute_V(n_basis):
"""Compute the matrix <k|x|l> for an electron in a box from -1 to 1 in a unit external field, in a.u."""
# initialize V to a zero matrix
V = np.zeros((n_basis,n_basis))
# Because Python is zero-indexed, our V matrix will be shifted by 1. I'll
# make this explicit by making the counters km1 and lm1 (k minus 1 and l minus 1)
for km1 in range(n_basis):
for lm1 in range(n_basis):
if (km1 + lm1) % 2 == 1:
# The matrix element is zero unless (km1 + lm1) is odd, which means that (km1 + lm1) mod 2 = 1.
# Either km1 is even or km1 is odd. If km1 is odd, then the km1 corresponds to a sine and
# lm1 is even, and corresponds to a cosine. If km1 is even and lm1 is odd, then the roles of the
# sine and cosine are reversed, and one needs to multiply the first term below by -1. The
# factor -1**lm1 achieves this switching.
V[km1,lm1] = 4. / np.pi**2 * (-1**((km1 - lm1 - 1)/2) / (km1-lm1)**2 * -1**(lm1)
+ -1**((km1 + lm1 - 1)/2) / (km1+lm1+2)**2)
return V
def energy_pt2(k,F,n_basis):
"""kth excited state energy in a.u. for an electron in a box of length 2 in the field F estimated with 2nd-order PT.
Parameters
----------
k : scalar, int
k = 0 is the ground state and k = 1 is the first excite state.
F : scalar
the external field strength
n_basis : scalar, int
the number of terms to include in the sum over states in the second order pert. th. correction.
Returns
-------
energy_pt2 : scalar
The estimated energy of the kth-excited state of the particle in a box of length 2 in field F.
"""
# It makes no sense for n_basis to be less than k.
assert(k < n_basis), "The excitation level of interest should be smaller than n_basis"
# Energy of the kth-excited state in a.u.
energy = (np.pi**2 / 8.) * np.array([(k + 1)**2 for k in range(n_basis)])
V = compute_V(n_basis)
der2 = 0
for j in range(n_basis):
if j != k:
der2 += 2*V[j,k]**2/(energy[k]-energy[j])
return energy[k] + der2 * F**2 / 2
def energy_variational(F):
"""ground state energy for a electron0in-a-Box in a box of length 2 in the field F estimated with the var. principle.
The variational wavefunction ansatz is psi(x) = (1+cx)cos(pi*x/2) where c is a variational parameter.
"""
gamma = 1/3 - 2/np.pi**2
def func(c):
return (np.pi**2/8*(1+gamma*c**2)+c**2/2 - 2*c*gamma*F) / ( 1 + gamma * c**2)
res = minimize_scalar(func,(0,1))
return res.fun
def energy_basis(F,n_basis):
"""Eigenenergies in a.u. of an electron in a box of length 2 in the field F estimated by basis-set expansion.
n_basis basis functions from the F=0 case are used.
Parameters
----------
F : scalar
the external field strength
n_basis : scalar, int
the number of terms to include in the sum over states in the second order pert. th. correction.
Returns
-------
energy_basis_exp : array_like
list of n_basis eigenenergies
"""
energy = (np.pi**2 / 8.) * np.array([(k + 1)**2 for k in range(n_basis)])
V = compute_V(n_basis)
# assign Hamiltonian to the potential matrix, times the field strength:
h = F*V
np.fill_diagonal(h,energy)
# solve Hc = Ec to get eigenvalues E
e_vals = eigh(h, None, eigvals_only=True)
return e_vals
print("Energy of these models vs. reference values:")
print("Energy of the unperturbed ground state (field = 0):", np.pi**2 / 8.)
print(" ")
print("Field value: ", 1./16)
print("Exact Energy of the ground state:", 10.3685/8 - 1./16)
print("Energy of the ground state estimated with 2nd-order perturbation theory:", energy_pt2(0,1./16,50))
print("Energy of the ground state estimated with the variational principle:", energy_variational(1./16))
print("Energy of the ground state estimated with basis set expansion:", energy_basis(1./16,50)[0])
print(" ")
print("Field value: ", 25./4)
print("Exact Energy of the ground state:", 32.2505/8 - 25./8)
print("Energy of the ground state estimated with 2nd-order perturbation theory:", energy_pt2(0,25./8,50))
print("Energy of the ground state estimated with the variational principle:", energy_variational(25./8))
print("Energy of the ground state estimated with basis set expansion:", energy_basis(25./8,50)[0])
print(" ")
print(" ")
print("Energy of the unperturbed first excited state (field = 0):", np.pi**2 * 2**2 / 8.)
print(" ")
print("Field value: ", 1./16)
print("Exact Energy of the first excited state:", 39.9787/8 - 1./16)
print("Energy of the first excited state estimated with 2nd-order perturbation theory:", energy_pt2(1,1./16,50))
print("Energy of the first excited state estimated with basis set expansion:", energy_basis(1./16,50)[1])
print(" ")
print("Field value: ", 25./4)
print("Exact Energy of the first excited state:", 65.177/8 - 25./8)
print("Energy of the first excited state estimated with 2nd-order perturbation theory:", energy_pt2(1,25./4,50))
print("Energy of the first excited state estimated with basis set expansion:", energy_basis(25./4,50)[1])
print(" ")
print(" ")
print("Energy of the unperturbed second excited state (field = 0):", np.pi**2 * 3**2 / 8.)
print(" ")
print("Field value: ", 1./16)
print("Exact Energy of the second excited state:", 89.3266/8 - 1./16)
print("Energy of the second excited state estimated with 2nd-order perturbation theory:", energy_pt2(2,1./16,50))
print("Energy of the first excited state estimated with basis set expansion:", energy_basis(1./16,50)[2])
print(" ")
print("Field value: ", 25./4)
print("Exact Energy of the ground state:", 114.309/8 - 25./8)
print("Energy of the second excited state estimated with 2nd-order perturbation theory:", energy_pt2(2,25./4,50))
print("Energy of the second excited state estimated with basis set expansion:", energy_basis(25./4,50)[2])
