I am trying to compute the ground-state energy associated to a nonlinear Schrödinger equation without the interaction term and with an external potential given by a 1-D harmonic potential. The expression for the energy involves the square of the absolute value of the gradient of the wave-function and that's where I got stuck.
After doing some research on the internet, I found (at least, I believe) that np.gradient could do the job. However, for my current problem, it hasn't shown any in the sense that when I plot its contribution it returns me a constant function equal to zero. Thus, I think I am doing something wrong.
My code is as follows:
import matplotlib.pyplot as plt
import numpy as np
import h5py as h5
from scipy import integrate
data = h5.File('groundstate.h5', 'r')
phireal = data['3']['phireal']
phiimag = data['3']['phiimag']
lattice = data['3']['y']
time = data['3']['t']
distr = np.power(phireal[:,:],2) + np.power(phiimag[:,:],2)
egradi = np.gradient(phiimag, axis = 0)
egradr = np.gradient(phireal, axis = 0)
V = (1/2) * np.power(lattice,2)
intergy = (1/2) * np.power(egradr,2) + V * distr
energy = integrate.simps(intergy, lattice, 0.1171875)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(time, energy)
plt.ylabel('E')
plt.xlabel('t')
ax.set_ylim(0, 3)
ax.set_xlim(0, 10)
plt.show()
The data stored in the file groundstate.h5 basically refers to the complex and real parts of the wave-function obtained evolving in imaginary time.
Here I provide a link for the input h5 file: https://drive.google.com/open?id=1FPM_sdpfQSOxeEikGuQyO4kfwH88wpcZ
In the plot of the energy, I expect to find for higher values of t (time) the value of 1/2 which is the energy of the ground-state of the 1-D quantum harmonic oscillator. However, I am getting a different value, 1/4 and that's why I believe it is because of the kinetic term, in this case, given by the gradient of the wave-function.
Can someone please tell me if this is the correct way of computing the gradient? Thanks in advance.