# 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