# inelastic_cc_2d.py
#
# Supplies functions to perform two-dimensional inelastic close-coupling
# calculations for scattering He atoms from corrugated surfaces.
#
# Patrick Kraus
# 06.02.14
# If this script is executed on its own, import neccesary libraries
#from pylab import *
#import sys
#from haslib import *
# Now import the libraries, which certainly are ONLY needed here
from haslib.science_basics import *
from multiprocessing import Pool
from scipy.special import iv
from numpy.polynomial.legendre import Legendre
from numpy import *
from numpy.linalg import norm
# -----------------------------------------------------------------------------------
# ----------------------- Constants and Flags ---------------------------------------
# -----------------------------------------------------------------------------------
# -----------------------------------------------------------------------------------
# ----------------------- Function definition ---------------------------------------
# -----------------------------------------------------------------------------------
[docs]def iCC_2D(h,AP,a,D,xi,Ei,theta_i,G1,G2,Gope=-1,Gclo=0,phon=[0,0,0],floq=1,sta_z=-7.0*A,end_z=20.0*A,st_p_wl=150,CC_border=1E-7,write_P=False,verbose=False,m=m_He3,maxrad=20):
"""
Two-dimensional calculation of scattering intensities,
taking inelastic contributions into account.
The CC-calculation is done using a Fox-Goodwin integrator
The potential function is assumed to be a Corrugated Morse shape
Patrick Kraus
06.Feb.2014
email: patrick.kraus@tugraz.at
Args:
h: Corrugation Parameter
AP: Inelastic coupling Parameter
a: 1D lattice constant (in meter)
D: Potential depth (in joule)
xi: Potential stiffnes of the central morse potential (in m**-1)
Ei: Incident helium energy (in joule)
theta_i: Incident angle of the scattering (in degrees)
G1: Reciprocal lattice vector 1 (in m**-1)
G2: Reciprocal lattice vector 2 (in m**-1)
Gope: Number of open channels to be considered - (Standard Value "-1" - all the channels)
Gclo: Number of closed channels to be considered - (Standard Value "0")
phon: Array carrying the phonon values array([Kval_x, Kval_y, Eval])
floq: Number of Floquet Blocks on both sides - default value is 1
sta_z: starting value of the numerov interation (in meter) - (Standard Value "7.0*A")
end_z: terminal value of the numerov interation (in meter) - (Standard Value "20.0*A")
st_p_wl: Number of steps to be performed in the interation within the length of the shortest calculated wavelength - (Standard Value "150")
CC_border: Minimum value for the coupling parameters to be considered - (Standard Value "1E-7")
write_P: Boolean Flag to determine if the wavefunctions will be written and returned (Standard Value - "False")
verbose: Boolean Flag to signal if Parameter Information should be displayed - (Standard Value "False")
Returns:
array: Scattering Intensities sorted by channels and Floquet Blocks (basically by energy)
::
[[Gx1 , Gx2 , Gx3 , ...],
[Gy1 , Gy2 , Gy3 , ...],
[F1 , F2 , F3 , ...],
[Kx1 , Kx2 , Kx3 , ...],
[Ky1 , Ky2 , Ky3 , ...],
[Kz2_1, Kz2_2, Kz2_3, ...],
[I1 , I2 , I3 , ...]]
Only considering open channels
or array, holding the wavefunctions of the open channels in its rows,
if the "write_P" Flag is set.
References:
.. [1] A.S. Sanz and S. Miret-Artes - Selective Adsorption Resonances: Quantum And Stochastic Approaches, Phys. Rep. 451 (2007), 37-154
.. [2] S. Miret-Artes - Multiphonon contributions to diffracted intensities and resonance profiles in atom-surface scattering, Surf. Sci. 366 (1996), L681-L684
.. [3] S. Miret-Artes - Resonant inelastic scattering of atoms from surfaces, Surf. Sci. 339 (1995), 205-220
"""
# ###############################################################################
# ----------------------- Channel calculation -----------------------------------
# ###############################################################################
# Phonon properties
phon_K = array([phon[0],phon[1]])
phon_E = phon[2]
# Kinematical properties
ki = sqrt(Ei/h2m(m))
Ki = array([ki*sin(theta_i*pi/180.0)])
# Determine Number of Floquet Blocks considered (including the 0 Block)
Nr_fl = 2*floq+1
# Calculation of help-parameters for determining a maximum radius
maxrad2 = 4+int(round(h,1)*10) # <-- maximal distance for the G-vectors to be considered
if maxrad < maxrad2:
maxrad = maxrad2
# Now we produce a list of G-vectors to be checked for relevance
grid = zeros([Nr_fl*(2*maxrad+3)**2,5])
run = 0
for ind1 in arange(2*maxrad+3)-(maxrad+1):
for ind2 in arange(2*maxrad+3)-(maxrad+1):
for flo in arange(-floq,floq+1,1):
grid[run,0] = G1[0]*ind1+G2[0]*ind2+flo*phon_K[0]
grid[run,1] = G1[1]*ind1+G2[1]*ind2+flo*phon_K[1]
grid[run,2] = ind1
grid[run,3] = ind2
grid[run,4] = flo
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)
# Produce a container with the chosen values
container = zeros([sum(flag),7])
container[:,0:5] = grid[flag,:]
Ki_x_new = Ki + container[:,0]
Ki_y_new = container[:,1]
# Calculate final energy
container[:,5] = Ki_x_new
container[:,6] = ki**2 - Ki_x_new**2 - Ki_y_new**2 - container[:,4]*phon_E/h2m(m)
# Now sort the channels kz2-wise
I = argsort(container[:,6])[::-1]
container = container[I,:]
# seperate into open and closed channels
chan_open = container[container[:,6]>=0,:]
chan_clos = container[container[:,6]<0,:]
# Now choose only the channels within the range given by the starting parameters
if Gope==-1 or Gope>chan_open.shape[0]:
pass
else:
choice = -1*Gope
chan_open = chan_open[choice:,:]
if Gclo==-1 or Gclo>chan_clos.shape[0]:
pass
else:
chan_clos = chan_clos[:Gclo,:]
# Read some simulation parameters resulting from those limits
Nr_open = chan_open.shape[0]
Nr_clos = chan_clos.shape[0]
Nr_C = Nr_open+Nr_clos
# Find specular channel
spe_c = where( (chan_open[:,2] == 0) &
(chan_open[:,3] == 0) &
(chan_open[:,4] == 0))[0][0]
# Combine closed and open containers again
channels = vstack([chan_open,chan_clos])
# Calculate the numerov step
step = 2*pi/(st_p_wl*sqrt(channels[:,6].max()))
# Definition of the z-Step list
z = arange(sta_z,end_z,step)
NrS = z.size # <-- Number of steps
# Time for the first information piece for the user
if verbose:
print("""
#-------------------------------------------------
iCC_2D v0.5 - 11.Feb.2014 - Patrick Kraus
including elastic Channel integration
#-------------------------------------------------
Channel Information:
Total Number of Channels --> {0}
Number of Floquet Blocks --> {1}
Number of reciprocal Lattice Positions --> {2}
Open channels --> {3}
Closed channels --> {4}
Integration Information:
Number of Steps --> {5}
Number of Steps per minimum wavelength --> {6}
Length of Step (A) --> {7}
""".format(Nr_C,Nr_fl,floor(Nr_C/Nr_fl).astype(int),Nr_open,Nr_clos,NrS,st_p_wl,step/A))
# ###############################################################################
# ----------------------- Coupling calculation ----------------------------------
# ###############################################################################
# Now a function to calculate the matrix for the system (at z==0)
# This matrix can stay constant throughout the calculation. To determine
# the z-dependence we only need to apply a multiplicative factor (for CMP)
Mas_0 = zeros([Nr_C,Nr_C]) # <-- will hold markers for specular contributions
Mas_F0 =zeros([Nr_C,Nr_C]) # <-- will hold markers for inelastic specular contributions
Mas_K = zeros([Nr_C,Nr_C]) # <-- will hold markers for K2z values on specular contributions
Mas_E = zeros([Nr_C,Nr_C]) # <-- will hold the asymptotic energy of spec. contrib.
for ind_l in range(Nr_C):
for ind_c in range(Nr_C):
if (channels[ind_l,2]==channels[ind_c,2])&(channels[ind_l,3]==channels[ind_c,3]):
if (channels[ind_l,4]==channels[ind_c,4]):
# Elastic specular contrib - mark Mas_0 and Mas_K
Mas_0[ind_l,ind_c] = 1
Mas_K[ind_l,ind_c] = channels[ind_l,6]
elif abs(channels[ind_l,4]-channels[ind_c,4])==1:
# Inelastic specular contribution - mark Mas_F0
Mas_F0[ind_l,ind_c] = 1
else:
Mas_E[ind_l,ind_c] = couplings(channels[ind_l,2]-channels[ind_c,2],channels[ind_l,3]-channels[ind_c,3],channels[ind_l,4]-channels[ind_c,4],h,a,xi,AP,CC_border)
# define a function to supply the coupling matrix for every z
def giveme_w(z):
# Define specular contributions
SPE = Mas_0*centr_pot(z,D,xi)/h2m(m) - Mas_K
# Define inelastic specular contributions (ignored a.t.m.)
ISP = 0*Mas_F0*centr_pot_dev(z,D,xi)/h2m(m)
# Weigh couplings with the distance
COU = Mas_E*D*exp(-2*xi*z)/h2m(m)
# return complete coupling matrix
return SPE+ISP+COU
K2z = channels[:,6]
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# ------------- Numerov Routine ----------------
# Preparation of used containers
R_00 = zeros([Nr_C,Nr_C]) # Filling it up with zeros - seemed to work
W_M1 = giveme_w(z[0]) # The first w-minus-one matrix, reset every run with a new one
W_00 = giveme_w(z[1]) # The first W-zero Matrix of the potential we are currently using - will be changed every run
I = eye(Nr_C) # Unity matrix - to save time
R_dict = {} # A dictionary to save the R-matrices - if needed
# Start the main loop
for ru in range(NrS-2):
# The loop really starts at value 1 - since the first two values have
# already be set
run = ru+1
# Calculate the w-plus-one-matrix for this position
W_P1 = giveme_w(z[run+1])
# I split the parts of the equation into blocks
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
#-----
R_00 = linalg.solve((AA-BB),CC)
# if the flag was set, the R-matrix is saved in the dictionary
if write_P:
R_dict[z[run]] = R_00
# And shift the used values for the next step
W_M1 = W_00
W_00 = W_P1
# ------------ K-Extraction ----------------
# Calculation of the S0 (Sinus-Zero) and the CE (Cosinus-Exponential) - Vectors
# for the last two z-positions
S0_M1 = zeros([Nr_C,1])
S0_00 = zeros([Nr_C,1])
CE_M1 = zeros([Nr_C,1])
CE_00 = zeros([Nr_C,1])
# Now, depending on the character of K2z - the values are set:
for index in range(Nr_C):
if K2z[index]>0:
# If the normal energy is positive, the corresponding sin and cos
# values are calculated and set
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 values are set for
# the use inside the decaying (hopefully) 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)
# Now everything should be prepared and we can solve the simple matrix
# equation
SS = S0_M1 - dot(R_00,S0_00)
CC = dot(R_00,CE_00) - CE_M1
KK = linalg.solve(CC,SS)
# Here the values of K and E are still mixed
# The values belonging to closed channels in KK are removed
geschl = where(K2z<0)[0]
K = delete(delete(KK,geschl,0),geschl,1)
NUM = K.shape[0]
S = linalg.solve((eye(NUM)-(0+1j)*K),eye(NUM)+(0+1j)*K)
S = abs(S)**2
CHA_x = channels[K2z>0,2].reshape([-1,1])
CHA_y = channels[K2z>0,3].reshape([-1,1])
FLO = channels[K2z>0,4].reshape([-1,1])
Kval_x = channels[K2z>0,5].reshape([-1,1])
Kval_y = channels[K2z>0,1].reshape([-1,1])
K2z = K2z[K2z>0].reshape([-1,1])
# Return the prepared values
return hstack([CHA_x,CHA_y,FLO,Kval_x,Kval_y,K2z,S[spe_c,:].reshape([-1,1])])
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
[docs]def two_D_el_coupling(a,h,xi,G1,G2,div=1.0):
"""
Calculate coupling factor
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
Results:
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)
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
# Define a function for the inter- and intrablock couplings
[docs]def couplings(dGx,dGy,dF,h,a,xi,B,CC_border):
# Now check for multi-phonon interaction
if abs(dF)>1:
return 0
# for all the other cases we can calculate
V_0 = two_D_el_coupling(a,h,xi,0,0)
V_G = two_D_el_coupling(a,h,xi,dGx,dGy,V_0)
# Now we only need to check if the interaction is elastic or inelastic
if dF==0:
# Elastic case
if V_G > CC_border:
return V_G
else:
return 0
else:
# inelastic case
if V_G > CC_border:
return -2*B*xi*V_G
else:
return 0
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
[docs]def centr_pot(z,D,xi):
# Define central interaction potential (Morse?)
return D*( exp(-2*xi*z) - 2*exp(-xi*z) )
[docs]def centr_pot_dev(z,D,xi):
return 2*D*xi*(exp(-xi*z) - exp(-2*xi*z))
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
[docs]def Temp_dependent_B(T,Qc,Td=155.0,a=4.3084*A,mas=1):
# function to calculate the temperature dependent inelastic coupling parameter
return sqrt( ( pi*T*384*hb**2 )/( mas*kb*Td**2 ) )/(Qc*a)
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
# Define a function to extract the elastic channel intensities
[docs]def get_iCC_elastics(phon_E,h,AP,a,D,xi,Ei,theta_i,G1,G2,G_ope,G_clos,floq=1,sta_z=-7.0*A,end_z=20.0*A,st_p_wl=150,CC_border=1E-7,write_P=False,verbose=False,only_int=True):
# The additional Parameter phon_E states the phonon energy considered
# Perform inelastic CC
phon = array([0,0,phon_E])
ZWE = iCC_2D(h,AP,a,D,xi,Ei,theta_i,G1,G2,G_ope,G_clos,phon,floq,sta_z,end_z,st_p_wl,CC_border,write_P,verbose)
ZWE = ZWE[ZWE[:,2].astype(int)==0,:]
if only_int==True:
return ZWE[:,6].reshape([1,-1])
else:
ZWE2 = zeros([3,ZWE.shape[0]])
ZWE2[0,:] = ZWE[:,0].transpose()
ZWE2[1,:] = ZWE[:,1].transpose()
ZWE2[2,:] = ZWE[:,6].transpose()
ZWE = ZWE2
return ZWE
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
[docs]def get_iCC_specs(phon_E,h,AP,a,D,xi,Ei,theta_i,G1,G2,G_ope,G_clos,floq=1,sta_z=-7.0*A,end_z=20.0*A,st_p_wl=150,CC_border=1E-7,write_P=False,verbose=False,only_int=True):
"""
return array of specular intensites for every floque block
"""
# The additional Parameter phon_E states the phonon energy considered
# Perform inelastic CC
phon = array([0,0,phon_E])
iccResult = iCC_2D(h,AP,a,D,xi,Ei,theta_i,G1,G2,G_ope,G_clos,phon,floq,sta_z,end_z,st_p_wl,CC_border,write_P,verbose)
specResult = iccResult[(iccResult[:,0]==0)&(iccResult[:,1]==0),:]
if only_int==True:
return specResult[:,6]
else:
return specResult
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# Define a function to integrate the elastic intensities up to the debye Frequency
[docs]def integr_el_intensity(h,AP,a,D,xi,Ei,theta_i,G1,G2,G_ope,G_clo,floq=1,sta_z=-7.0*A,end_z=20.0*A,st_p_wl=150,CC_border=1E-7,write_P=False,verbose=False,wD=10*mel,n_leg=10):
# The additional Parameter wD states the Debye energy
# First calculate the roots of the Legendre polynomial of n_th order
LP = Legendre.basis(n_leg)
LP_D = LP.deriv()
x_i = LP.roots()
# Now the corresponding weights
w_i = 2.0 / ( ( 1 - x_i**2 )*( LP_D(x_i)**2 ) )
# Now transform the Legendre roots to the desired interval
x_i_mod = wD*0.5*x_i+ wD*0.5
# Loop over the get_iCC_elastics function to obtain all the elastic
# intensities for the different phonon modes
run = 0
for E_now in x_i_mod:
if run==0:
ARR_mod = get_iCC_elastics(E_now,h,AP,a,D,xi,Ei,theta_i,G1,G2,G_ope,G_clo,floq,sta_z,end_z,st_p_wl,CC_border,verbose=False,only_int=False)
ARR = ARR_mod[2,:]
run = run+1
else:
ARR = vstack([ARR,get_iCC_elastics(E_now,h,AP,a,D,xi,Ei,theta_i,G1,G2,G_ope,G_clo,floq,sta_z,end_z,st_p_wl,CC_border,verbose=False,only_int=True)])
# Now multiply with the corresponding weights and the Debye distribution
rho = 3*x_i_mod**2/wD**3
INT = sum(ARR*tile(reshape(w_i,[-1,1]),[1,ARR.shape[1]])*tile(reshape(rho,[-1,1]),[1,ARR.shape[1]]),0)*wD*0.5
#INT = INT/sum(INT)
return vstack([ARR_mod[0,:],ARR_mod[1,:],INT])
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
[docs]def integrate_GL(a,b,f,n=10,cpus=None,args=[]):
"""
integrate f from a to b using n net points
Method: Gauss-Legendre quadrature
f can return a 2d array and the integration will be per element
"""
LP = Legendre.basis(n)
LP_D = LP.deriv()
# calculate mesh points
x_i = LP.roots()
# calculate weights
w_i = 2.0 / ((1 - x_i**2) * (LP_D(x_i)**2))
# Now transform the Legendre roots to the desired interval
x_i_mod = ((b - a) * x_i + (a + b)) / 2.0
#create parameter list
params = [[x]+args for x in x_i_mod]
# evaluate functions at the mesh points
if cpus:
p = Pool(cpus)
points = dstack(p.map(f,params))
p.close()
else:
points = dstack(map(f,params))
# multiply every layer [:,:,i] by w[i] and sum along the 3rd axis
result = sum(points * w_i, 2)
# multiply by factor due to interval change
result = result * (b - a) / 2
return result
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
[docs]def calc_gauss_legendre(wD,n):
"""
calculate positions and width for gauss legendre integration
wD Integrate from -wD to wD
n order of Legendre polynomial to use
return array of phonon ernergies
"""
LP = Legendre.basis(n)
LP_D = LP.deriv()
x_i = LP.roots()
# Now the corresponding weights (Gauss Legendre quadrature)
w_i = 2.0 / ((1 - x_i**2) * (LP_D(x_i)**2))
# Now transform the Legendre roots to the desired interval
x_i_mod = wD * 0.5 * x_i + wD * 0.5
return x_i_mod,w_i
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# Define a function to integrate the elastic intensities up to the debye Frequency
[docs]def integrate_spec_int(h,AP,a,D,xi,Ei,theta_i,G1,G2,G_ope,G_clo,floq=1,sta_z=-7.0*A,end_z=20.0*A,st_p_wl=150,CC_border=1E-7,write_P=False,verbose=False,wD=10*mel,n_leg=10):
phon_energies,w = calc_gauss_legendre(wD,n_leg)
# Loop over the get_iCC_elastics function to obtain all the elastic
# intensities for the different phonon modes
intensities = zeros([len(phon_energies),floq*2+1])
for ind,phon_energy in enumerate(phon_energies):
intensities[ind,:] = get_iCC_specs(phon_energy,h,AP,a,D,xi,Ei,theta_i,
G1,G2,G_ope,G_clo,floq,sta_z,end_z,st_p_wl,CC_border,verbose=False)
# Now multiply with the corresponding weights and the Debye distribution
rho = 3.0 * phon_energies**2 / wD**3
result = sum(transpose(intensities)*(rho*w),1)*wD*0.5
return result
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
if __name__=="__main__":
a = 4.3084*A
Sb111 = Crystal(a,a,120) # Define Sb structure
Sb111.rotate_crystal(1,0) # Define scattering direction
Gx = Sb111.b1[0:2]
Gy = Sb111.b2[0:2]
for i in range(10):
ERG = iCC_2D(0.1,0,4.308*A,4.196*mel,0.38/A,13*mel,45.75,Gx,Gy,Gope=-1,Gclo=2,phon=array([0,0,0]),floq=0,sta_z=-7.0*A,end_z=20.0*A,st_p_wl=100,CC_border=1E-4,verbose=True)
print(ERG)