Source code for haslib.morse

#from pylab import *
from scipy.special import iv
from scipy import optimize
import itertools
from matplotlib import animation
from matplotlib.colors import LogNorm
from .science_basics import *
from numpy import array,sqrt,zeros

[docs]class Morse: def __init__(self,crystal,h,xi,D): self.crystal = crystal self.a = crystal.a1[0] self.h = h self.xi = xi self.D = D
[docs] def calc_morse_potential(self, a, h, xi, D, z_list, gmax): cc00 = self.calc_coupl_2d_hex(a, h, xi, 0, 0) gg_range = xrange(-2 * gmax , 2 * gmax + 1) pot = zeros([len(gg_range), len(gg_range), len(z_list)]) for i, gg in enumerate(itertools.product(gg_range, gg_range)): gg_x = gg[0] gg_y = gg[1] for z_ind, z in enumerate(z_list): if gg_x == 0 and gg_y == 0: pot[2 * gmax, 2 * gmax, z_ind] = D * (exp(-2 * xi * z) - 2 * exp(-xi * z)) else: cc = self.calc_coupl_2d_hex(a, h, xi, gg_x, gg_y, cc00) pot[gg_x + 2 * gmax, gg_y + 2 * gmax, z_ind] = cc * D * exp(-2 * xi * z) return pot
[docs] def calc_coupl_2d_hex(self, a, h, xi, G1, G2, div=1.0): """ Calculate coupling factor Args: a: Edge length of the primitive cell hexagon (in meters) h: Corrugation parameter of the cosine corrugation function xi: Stiffnes parameter of the corrugated morse parameter (in m**-1) G1: First index of G-Vector to the channel G2: Second index of G-Vector to the channel div: divison parameter - The solution will be divided by this number returns: scalar: Coupling factor from (0,0) to (G1,G2) """ kmax = 10 k = arange(-1 * kmax, kmax + 1, 1) alpha = 2 * xi * h * a / 3.0 qu = iv(k, alpha) * (iv(k + G2, alpha) * iv(k - G1, alpha) + iv(k - G2, alpha) * iv(k + G1, alpha)) return sum(qu) / div
[docs] def show(self, pot): fig = figure() a=pot quad = pcolormesh(a[:, :, 7], norm=LogNorm(vmin=10E-100, vmax=a.max())) l = len(a[:, 0, 0]) xlim([0, l]) ylim([0, l]) title("Potential") colorbar() def animate(i): title("z-index:{}".format(i)) quad.set_array(a[:, :, i].ravel()) anim = animation.FuncAnimation(fig, animate, frames=len(a[0, 0, :]), blit=False) show()
[docs] @staticmethod def boundstate(alpha, D, xi, m): omg = xi * sqrt(2 * D / m) gam = 2 * D / (hb * omg) E_alpha = -D + hb * omg * (alpha + 0.5) * (1 - (alpha + 0.5) / (2 * gam)) return E_alpha
[docs] @staticmethod def fit(bse,index,m): # Definiere Parameter param_0 = array([5.0 * meV, 0.2 / A]) def func(a,b,c): return Morse.boundstate(a,b,c,m) popt, pcov = optimize.curve_fit(func, index, bse, p0=param_0) D,xi = popt omg = xi * sqrt(2 * D / m) gam = 2 * D / (hb * omg) print ('maxbound = ' + str(gam + 0.5)) print ('D = ' + str(D / meV)[0:5] + '+/-' + str(sqrt(pcov[0][0]) / meV)[0:5] + ' meV') print ('xi = ' + str(xi * A)[0:5] + '+/-' + str(sqrt(pcov[1][1]) * A)[0:5] + ' A**-1') for i in range(int(gam + 0.5)): print (Morse.boundstate(i, D, xi,m) / meV) cont_E = zeros(len(index)) for n in range(len(index)): cont_E[n] = Morse.boundstate(n, D, xi,m) print (sqrt(sum((bse / meV - cont_E / meV) ** 2)) / len(index))