# 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)