Source code for haslib.CC

# 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
from numpy import (arange, argsort, array, cos, delete, diagflat, dot, exp,
                   eye, hstack, logical_and, ones, pi, sin, sqrt, vstack,
                   where, zeros)
from numpy.linalg import inv, norm, solve
from scipy.special import iv

from .science_basics import A, h2m, m_He3

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


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

[docs]def CC_2D_old(h, a, D, xi, Ei, theta_i, G1, G2, G_ope=-1, G_clo=0, sta_z=-7.0*A, end_z=20.0*A, st_p_wl=150, CC_border=1E-7, write_P=False, maxrad=20, verbose=False, m=m_He3): """ 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: h: Corrugation parameter a: 1D lattice constant (in meter) D: depth of the surface averaged interaction potential (in joule) xi: stiffness parameter of the surface averaged interaction potential (in m**-1) 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) sta_z: z-position where the integration should start (in meter) end_z: z-position where the integration should end (in meter) st_p_wl: number of iteration steps to be done within the shortest occuring wavelength CC_border: lowest value for a coupling factor to be considered 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 an array containing the values of the channel's wavefunctions is returned 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+3)**2, 4]) run = 0 for ind1 in arange(2*maxrad+3)-(maxrad+1): for ind2 in arange(2*maxrad+3)-(maxrad+1): grid[run, 0] = G1[0]*ind1+G2[0]*ind2 grid[run, 1] = G1[1]*ind1+G2[1]*ind2 grid[run, 2] = ind1 grid[run, 3] = ind2 run += 1 # Now determine the distances of the channels - ignore everything farther # away than maxrad dist = sqrt(grid[:, 0]**2+grid[:, 1]**2) flag = dist <= maxrad*norm(G1) grid = grid[flag, :] # calculation of the normal energies of each channel NrC = grid.shape[0] channels = zeros([NrC, 3]) ki = sqrt(Ei/h2m(m)) Ki = array([ki*sin(theta_i*pi/180.0)]) for Gn in range(grid.shape[0]): G_scat = grid[Gn, [0, 1]] kz2 = ki**2 - (Ki+G_scat[0])**2 - G_scat[1]**2 channels[Gn, 0] = grid[Gn, 2] channels[Gn, 1] = grid[Gn, 3] channels[Gn, 2] = kz2 # Now we calculate the distance of each channel to Gamma - everything # farther away than G1*maxrad is ignored and cut out of the calculation Cx = G1[0]*channels[:, 0]+G2[0]*channels[:, 1] Cy = G1[1]*channels[:, 0]+G2[1]*channels[:, 1] dist = sqrt(Cx**2+Cy**2) flag = dist <= maxrad*norm(G1) # the container for the coupling calculation container_neu = channels[flag, :] channels = container_neu # calculation of the number of needed numerov-steps, # definition of the z-list step = 2*pi/(st_p_wl*sqrt(channels[:, 2].max())) # print(step) z = arange(sta_z, end_z, step) NrS = z.size # Number of steps # ----------- Calculation of the couplings -------------------- # Definition of the coupling and the potential # surface averaged interaction potential def V0(z): return D*(exp(-2*xi*z) - 2*exp(-xi*z)) CCS = {} # The CCS dicitonary contains all coupling parameters to be considered # We now calculate a dictionary containing all the coupling parameters # for the different channels. All the values too small to be considered # should be included with a zero-value # the function calc_coupl_2d_hex is included in the same script and should # be exchanged if a different geometry is considered # Calculate zero-value: V_0 = calc_coupl_2d_hex(a, h, xi, 0, 0) for run in range(container_neu.shape[0]): ggg = calc_coupl_2d_hex(a, h, xi, int( container_neu[run, 0]), int(container_neu[run, 1]), V_0) if logical_and(abs(ggg) > CC_border, ggg > 0): CCS[(int(container_neu[run, 0]), int(container_neu[run, 1]))] = ggg # Set (0,0) channel to 0, so we don't run into problems with the coupling matrix CCS[(0, 0)] = 0 # Now we sort and choose the 'interesting' channels I = argsort(-1*channels[:, 2]) channels = channels[I, :] # looking for the specular channels spe_c_x = where(channels[:, 0].astype(int) == 0)[0] spe_c_y = where(channels[:, 1].astype(int) == 0)[0] for cx in spe_c_x: if cx in spe_c_y: spe_c = cx # Now everything is sorted, the open channels are atop, the closed afterwards. # We need to know the position of the border # last open channel: LastOpen = where(channels[:, 2] > 0)[0].max() # Now we separate open and closed channels, since they need to be treated # differently chan_op = channels[:LastOpen+1, :] chan_cl = channels[LastOpen+1:, :] # the condition for open channels if ((G_ope != -1) & (G_ope < LastOpen+1)): chan_op = chan_op[-1*G_ope:, :] # and the closed ones - finalization if G_clo == 0: channels = chan_op elif (G_clo == -1 | G_clo > chan_cl.shape[0]): channels = vstack([chan_op, chan_cl]) else: channels = vstack([chan_op, chan_cl[:G_clo, :]]) # set the coupling matrix NrD = channels.shape[0] Mask = zeros([NrD, NrD]) # read-out the final normal energies K2z = channels[:, 2] if verbose: print(""" #------------------------------------------------- CC_2D v0.5 - 11.Feb.2014 - Patrick Kraus including elastic Channel integration #------------------------------------------------- Channel Information: Total Number of Channels --> {0} Open channels --> {1} Closed channels --> {2} Integration Information: Number of Steps --> {3} Number of Steps per minimum wavelength --> {4} Length of Step (A) --> {5} """.format(NrD, len(chan_op), NrD-len(chan_op), NrS, st_p_wl, step/A)) # 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 in range(NrD): for ind_c in 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 (gx_d, gy_d) in CCS.keys(): Mask[ind_l, ind_c] = CCS[(gx_d, gy_d)] # print Mask # This function calculates the dynamic coupling matrix for the integration steps def giveme_w(z): # The calculated Mask variable is weighed with the exponential decay for the distance w = Mask*D*exp(-2*xi*z) # The only thing left out are the diagonal elements which are dependent # on the surface averaged interaction potential w = (w+diagflat(V0(z)*ones(NrD)))/h2m(m) - diagflat(K2z) return w # ------------- Numerov Routine ---------------- # preparation of the needed containers R_00 = zeros([NrD, NrD]) # arbitrarily chosen zero-container # the first w-minus-one matrix, they will be passed over from every previous run from now on W_M1 = giveme_w(z[0]) # the first w-zero matrix of the actual position we sit on - will be passed to w-monus-one in the next step W_00 = giveme_w(z[1]) # a stored unitary matrix- so it won't be recalculated again I = eye(NrD) # a dictionary to store the R-matrices, if asked for (not implemented) R_dict = {} # start the iteration for ru in range(NrS-2): # The loop needs to be started with the value 2 - since the first two # values have been set from the beginning on run = ru+1 # calculate the w-plus-one matrix W_P1 = giveme_w(z[run+1]) # performing a Numerov step AA = 2*I+10.0*step**2*W_00/12.0 BB = dot(I-step**2*W_M1/12.0, R_00) CC = I-step**2*W_P1/12.0 # ----- AB = inv(AA-BB) # ----- R_00 = dot(AB, CC) # if the flag was set, the R-matrix is stored in the R-dictionary if write_P: R_dict[z[run]] = R_00 # preparations for the next step W_M1 = W_00 W_00 = W_P1 # ------------ 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 for index in range(NrD): if K2z[index] > 0: # If the normal energy is positive, the corresponding Sine-Cosine # terms are calculated and written kz = sqrt(K2z[index]) 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(K2z[index])) 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 = where(K2z < 0)[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_x = channels[K2z > 0, 0] CHA_y = channels[K2z > 0, 1] spe_c_x = where(CHA_x.astype(int) == 0)[0] spe_c_y = where(CHA_y.astype(int) == 0)[0] for cx in spe_c_x: if cx in spe_c_y: spe_c = cx K2z = K2z[K2z > 0] # Look for the specular channel return hstack([CHA_x.reshape(-1, 1), CHA_y.reshape(-1, 1), S[spe_c, :].reshape(-1, 1), K2z.reshape(-1, 1)])
[docs]def calc_coupl_2d_hex(a, h, xi, G1, G2, div=1.0): """ Calculate coupling factors 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))/div return sum(qu)
# 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]) # ERG = CC_2D(0.1,4.3084*A,4.196*mel,0.38/A,13*mel,45.75,Gx,Gy,G_ope=-1,G_clo=2) # print ERG