Source code for haslib.cc2d_new

# 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 matplotlib.pylab import *
from matplotlib.mlab import find
from scipy.special import iv
from scipy import constants as con
from haslib.fox_goodwin import *
from .science_basics import *
from numpy import*
from numpy.linalg import *

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


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

[docs]def CC_2D(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,verbose=False,maxrad=20,m=m_He3): #def CC_2D(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,verbose=False,maxrad=20,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, 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]) channels = [] x1 = range(-(maxrad+1),(maxrad+1)+1) ki = sqrt(Ei/h2m(m)) 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 = np.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 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 channel in channels: ggg = calc_coupl_2d_hex(a,h,xi,channel[0],channel[1],V_0) if logical_and(abs(ggg) > CC_border, ggg > 0): CCS[(channel[0],channel[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 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) Mask = zeros([NrD,NrD]) # read-out the final normal energies K2z = [c[2] for c in channels] 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, 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 (gx_d,gy_d) in CCS.keys(): Mask[ind_l,ind_c] = CCS[(gx_d,gy_d)] # This function calculates the dynamic coupling matrix for the integration steps 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 = Mask*D*exp(-2*xi*z) # The only thing left out are the diagonal elements which are dependent # on the surface averaged interaction potential fill_diagonal(diagV0,V0(z)) w = (w+diagV0)/h2m(m) - 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 = linalg.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))
kmax = 10 k = arange(-1*kmax,kmax+1,1)
[docs]def calc_coupl_2d_hex(a,h,xi,G1,G2,div=1.0): """ Calculate coupling factors between (0,0) and (G1,G2) 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
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) ERG = CC_2D(h,a,D,xi,13*mel,45.75,Gx,Gy,z=z_list,G_ope=-1,G_clo=2) print (ERG)