Source code for haslib.cc2d_ext

# CC_code_2d.py
#
# This script contains a python-function to supply a method for calculating the
# intensities of scattered helium atoms from a two-dimensional grid using the
# close-coupling formalism.
# The example coupling parameters supplied are given in an extra function and
# describe the coupling parameters for a hexagonal-like grid as it is present
# on the (111) surfaces of bismuth and antimony.
# The interaction potential is assumed to be of a Morse-shape but may be
# interchanged by every other potential shape .
#
# Patrick Kraus
# 09.July.2013
# comments translated
# 02.June.2014

# Import needed libraries
import itertools
from numpy import (zeros, sqrt, pi, delete, exp, column_stack,
                   sin, cos, diagflat, dot, eye, array, arange, load)
from numpy.linalg import(norm, inv, solve)
from scipy import constants as con
from .fox_goodwin import fox_goodwin

# ----------------------------------------------------------------------------
# --------------------- Konstanten und Flags ---------------------------------
# ----------------------------------------------------------------------------

A = con.angstrom                        # Angstroem
m = con.m_u*3.0160293                   # Mass of a helium atoms (modified He3)
hb = con.hbar                           # reduced Planck's constant
h2m = hb**2/(2*m)                       # Energy-Momentum-factor
mel = con.eV/1E3                        # milli-electronvolt

# ----------------------------------------------------------------------------
# --------------------- Funktionsdefinition  ---------------------------------
# ----------------------------------------------------------------------------


[docs]def CC_2D_e(Ei, theta_i, G1, G2, potential, z, G_ope=-1, G_clo=0, write_P=False): """ One-dimensional scattering intensity calculation Integrating the CC equations using the Fox-Goodwin integrator The potential functions is assumed to be a Morse function Args: Ei: energy of the impinging helium beam (in joule) theta_i: incident scattering angle of the helium beam with respect to the surface normal (in degrees) G1: reciprocal lattice vector 1 (in m**-1) G2: reciprocal lattice vector 2 (in m**-1) G_ope: number of open channels to be considered in the calculation (-1 for all) G_clo: number of closed channels to be considered in the calculation (-1 for all) potential z: array containing z values write_P: boolean value stating if the wavefunctions should be kept and returned (not fully implemented) Returns: array: Scattering intensities listed by channels following the scheme: :: [[Gx1, Gy1, I1, kz2_1] [Gx2, Gy2, I2, kz2_2] [Gx3, Gy3, I3, kz2_3] [... ... ... ... ]] (only listing open channels) or array containing the values of the channel's wavefunctions, if write_P was set to 'TRUE' (not implemented yet) """ # calculation of help-parameters maxrad = 20 # <-- maximum G-radius of reciprocal lattice vectors to be considered # We now construct a list of reciprocal lattice vectors to be tested for # the relevance of their coupling parameters grid = zeros([(2*maxrad+1)**2, 4]) channels = [] x1 = range(-(maxrad), (maxrad)+1) ki = sqrt(Ei/h2m) Ki = ki*sin(theta_i*pi/180.0) # calculation of the normal energies of each channel # Now we calculate the distance of each channel to Gamma - everything # farther away than G1*maxrad is ignored and cut out of the calculation for run, ind in enumerate(itertools.product(x1, x1)): g1 = G1[0]*ind[0]+G2[0]*ind[1] g2 = G1[1]*ind[0]+G2[1]*ind[1] dist = sqrt(g1**2+g2**2) if dist > maxrad*norm(G1): continue kz2 = ki**2 - (Ki+g1)**2 - g2**2 channels.append([ind[0], ind[1], kz2]) # calculation of the number of needed numerov-steps, # definition of the z-list e_max = max(channels, key=lambda x: x[2])[2] #step = 2*pi / (st_p_wl * sqrt(e_max)) #z = arange(sta_z, end_z, step) NrS = z.size # Number of steps # ----------- Calculation of the couplings -------------------- # Definition of the coupling and the potential # V0 = lambda z: D*( exp(-2*xi*z) - 2*exp(-xi*z) ) # surface averaged interaction potential # Now we sort and choose the 'interesting' channels channels.sort(key=lambda x: -x[2]) # Now everything is sorted, the open channels are atop, the closed afterwards. # We need to know the position of the border # last open channel: last_open = [i for i, c in enumerate(channels) if c[2] > 0][-1] # Now we separate open and closed channels, since they need to be treated # differently chan_op = channels[:last_open+1] chan_cl = channels[last_open+1:] # the condition for open channels if ((G_ope != -1) & (G_ope < last_open + 1)): chan_op = chan_op[-G_ope:] # and the closed ones - finalization if G_clo == 0: channels = chan_op elif (G_clo == -1 | G_clo > len(chan_cl)): channels = chan_op + chan_cl else: channels = chan_op + chan_cl[:G_clo] # set the coupling matrix NrD = len(channels) pot_size = potential.shape[0] Mask = zeros([NrD, NrD]).astype(int) # read-out the final normal energies K2z = [c[2] for c in channels] gmax = (pot_size-1)/4 # This loop calculates the relative distances between the respective channels # and uses this information to deduct the needed coupling factor between # these channels # This information is stored in a Mask-variable, so this constant relationship # doesn't need to be calculated a second time. for ind_l, ind_c in itertools.product(range(NrD), range(NrD)): # Gx and Gy-distance gx_d = int(channels[ind_l][0]-channels[ind_c][0]) gy_d = int(channels[ind_l][1]-channels[ind_c][1]) # This distance-coupling relationship isrecorded in the CCS dictionary if ind_l == ind_c: Mask[ind_l, ind_c] = int(2*gmax+pot_size*2*gmax) else: Mask[ind_l, ind_c] = int(gy_d+2*gmax+(pot_size)*(gx_d+2*gmax)) Mask = Mask.ravel() # This function calculates the dynamic coupling matrix for the integration steps V0 = potential[2*gmax, 2*gmax, :] diagK2z = diagflat(K2z) diagV0 = zeros([NrD, NrD]) def giveme_w(z, i): # The calculated Mask variable is weighed with the exponential decay for the distance w = potential[:, :, i].ravel()[Mask] w.shape = [NrD, NrD] # The only thing left out are the diagonal elements which are dependent # on the surface averaged interaction potential w = (w)/h2m - diagK2z return w # ------------- Numerov Routine ---------------- R_00 = fox_goodwin(z, giveme_w) # ------------ K-Extraction ---------------- # Creation of a S0 (Sine-zero) and the CE (Cosine-Exponential) - Vectors # for the last z-positions S0_M1 = zeros([NrD, 1]) S0_00 = zeros([NrD, 1]) CE_M1 = zeros([NrD, 1]) CE_00 = zeros([NrD, 1]) # Now, depending on the character of the respective K2z-value, everything is set run = len(z)-2 for index, erg in enumerate(K2z): if erg > 0: # If the normal energy is positive, the corresponding Sine-Cosine # terms are calculated and written kz = sqrt(erg) S0_M1[index] = sin(kz*z[run])/sqrt(kz) S0_00[index] = sin(kz*z[run+1])/sqrt(kz) CE_M1[index] = cos(kz*z[run])/sqrt(kz) CE_00[index] = cos(kz*z[run+1])/sqrt(kz) else: # If the normal energy is negative, only the CE-term is set with the exponentials kz = sqrt(abs(erg)) CE_M1[index] = exp(-kz*z[run])/sqrt(kz) CE_00[index] = exp(-kz*z[run+1])/sqrt(kz) S0_M1 = diagflat(S0_M1) S0_00 = diagflat(S0_00) CE_M1 = diagflat(CE_M1) CE_00 = diagflat(CE_00) # After everything is prepared, the matrix equation can be easily solved SS = S0_M1 - dot(R_00, S0_00) CC = dot(R_00, CE_00) - CE_M1 KK = solve(CC, SS) # Here the values for K and E are mixed. We now remove the values from KK # that belong to the closed channels geschl = [i for i, c in enumerate(K2z) if c < 0] K = delete(delete(KK, geschl, 0), geschl, 1) NUM = K.shape[0] S = dot(inv(eye(NUM)-(0+1j)*K), eye(NUM)+(0+1j)*K) S = abs(S)**2 cha_k2z_gt_0 = [c for c in channels if c[2] >= 0] CHA_x, CHA_y, K2z = zip(*cha_k2z_gt_0) spe_c = next(cx for cx, c in enumerate( cha_k2z_gt_0) if c[0] == 0 and c[1] == 0) return column_stack((CHA_x, CHA_y, S[spe_c, :], K2z))
if __name__ == "__main__": print(""" 2D Close-Coupling code using the Fox-Goodwin integrator Patrick Kraus 09.July.2013 Example code performing He-scattering simulation using scattering parameters theta_i = 45.75 degrees E_i = 13 meV in Gamma-M direction of the Sb(111) surface using a plain Morse surface-averaged interaction potential with D = 4.196 meV xi = 0.38E10 m**-1 """) a = 4.3084*A Gx = array([1.68396562e+10, 0.00000000e+00]) Gy = array([8.41982808e+09, 1.45835700e+10]) #h = 0.078 #a = 4.38*A #D = 4.872*mel #xi = 0.218/A gmax = 10 z_list = arange(-8*A, 16*A, 0.01*A) pot = load("test_pot.npy") ERG = CC_2D_e(13*mel, 45.75, Gx, Gy, G_ope=-1, G_clo=2, potential=pot, z=z_list) print(ERG)