Source code for haslib.icc_2d

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