from numpy import * # constants a0 = 0.529177 # in Anstroem Ha = 27.21138 # in eV # define grid startx=7.5 endx = 10.5 num_grid_points = 200 mass = 131.29*1836 # Xenon plot_num_states = 5 scale_wavefunctions = 0.0001 # generate grid x = linspace(startx, endx, num_grid_points) dx = (endx-startx)/float(num_grid_points) # define potentials def HomogeneousPotential(grid,constant=0): return [constant for i in grid] def HarmonicPotential(grid, alpha=1., beta=0., gamma=0.): return [ alpha*x**2 + beta*x+gamma for x in grid] def LJpotential(grid,alpha,sigma): return [ 4*alpha*( (sigma/x)**12 - (sigma/x)**6. ) for x in grid] # define Operators def OpV(potential): return diag(potential) def OpT(M): """ Kinetic energy (1/(2M)*Laplacian) approximated by (-1/2*\delta_{i\pm1}+\delta_ii)/(M*dx**2) """ T = diag([1./(dx**2*M) for i in range(num_grid_points)]) for i in range(num_grid_points-1): T[i,i+1] = -0.5/M/dx**2 T[i+1,i] = -0.5/M/dx**2 return T # change potentials here #potential = HarmonicPotential(x) #potential = HomogeneousPotential(x) potential = LJpotential(x, 0.02/Ha, 3.98/a0) # generate Hamiltonian and solve it H = OpT(mass) + OpV(potential) energies, wf = linalg.eigh(H) # calculate expectation values expR = [] fluctR = [] for s in range(30): R = sum(wf[:,s]**2*x) expR.append(R) fluctR.append(sqrt(sum(wf[:,s]**2*x**2)-R**2)) # gnuplot output f = open("plotfile.gnu", 'w') scaledWF = scale_wavefunctions*wf print >>f, """ #set terminal postscript enh solid color "Helvetica" 22 lw 2 #set output "solved.ps" set yrange[:%f] # set xrange[:] plot """ % ceil(energies[plot_num_states]), for i in range(plot_num_states): print >>f, energies[i], " lw 2 lt 2 t '', ", print >>f, "\\" print >>f, """ "-" t "potential" with lines lt 1 lw 2,\\""" for i in range(plot_num_states-1): print >>f, """ "" t "" with lines lw 1 lt 3,""", print >>f, """ "" t "" with lines lw 1 lt 3,""", print >>f, """ "" t "" with lines lw 2 lt 5,""", print >>f, """ "" t "" with lines lw 1 lt 7,""", print >>f, """ "" t "" with lines lw 1 lt 7\n""" for (i,v) in enumerate(potential): print >>f, x[i], v print >>f, "e" for i in range(plot_num_states): for (j,v) in enumerate(scaledWF[:,i]): print >>f, x[j],v+energies[i] print >>f, "e" showR = 6 for i in range(showR): print >>f, expR[i], energies[i] print >>f, "e" for i in range(showR): print >>f, expR[i]-fluctR[i], energies[i] print >>f, "e" for i in range(showR): print >>f, expR[i]+fluctR[i], energies[i] f.close()