msspec_python3/doc/source/tutorials/AlN/AlN.py

182 lines
8.3 KiB
Python

# coding: utf-8
from msspec.utils import *
from ase.build import bulk
from ase.visualize import view
import numpy as np
from msspec.calculator import MSSPEC, XRaySource
from msspec.iodata import Data
from itertools import product
DATA = None
def AlN_cluster(emitter_tag, emitter_plane, diameter=0, planes=0, term_anion=True, tetra_down=True):
"""
This function is a kind of overload of the hemispherical_cluster function with specific attributes
to the AlN structure
:param emitter_tag: An integer that allows to identify the kind of atom to use as the emitter
:param emitter_plane: The plane where the emitter is. 0 means the surface.
:param diameter: The diameter of the cluster (in Angstroms).
:param planes: The total number of planes.
:param term_anion: True if the surface plane is anion terminated, False otherwise
:param tetra_down: The orientation of the tetrahedral
:return:
"""
# create the initial cluster of AlN
cluster = bulk('AlN', crystalstructure='wurtzite', a=3.11, c=4.975)
# tag each atom in the unit cell because they are all in a different chemical/geometrical environment
# (0 and 2 for Aluminum's atoms and 1 and 3 for Nitride's atoms)
[atom.set('tag', i) for i, atom in enumerate(cluster)]
# change the orientation of the tetrahedron
if not tetra_down:
cluster.rotate(180, 'y')
# From this base pattern, creat an hemispherically shaped cluster
cluster = hemispherical_cluster(cluster, emitter_tag=emitter_tag, emitter_plane=emitter_plane, diameter=diameter,
planes=planes)
# Depending on the number of planes above the emitter, the termination may not be the one you wish,
# we test this and raise an exception in such a case
# Get the termination by finding the kind of one atom located at the topmost z coordinate
termination = cluster[np.argsort(cluster.get_positions()[:, 2])[-1]].symbol
# test if the combination of parameters is possible
assert (termination == 'N' and term_anion) or (termination == 'Al' and not term_anion), (
"This termination isn't compatible with your others parameters, you must change the tag of "
"your emitter, the plane of your emitter or your termination")
return cluster
def create_clusters(side='Al', emitter='Al', diameter=15, planes=6):
clusters = []
if side == 'Al':
term_anion = tetra_down = False
elif side == 'N':
term_anion = tetra_down = True
if emitter == 'Al':
tags = (0, 2)
level = '2p'
ke = 1407.
elif emitter == 'N':
tags = (1, 3)
level = '1s'
ke = 1083.
for emitter_tag, emitter_plane in product(tags, range(0, planes)):
nplanes = emitter_plane + 2
try:
cluster = AlN_cluster(emitter_tag, emitter_plane, diameter=diameter, planes=nplanes, term_anion=term_anion,
tetra_down=tetra_down)
except AssertionError:
continue
cluster.absorber = get_atom_index(cluster, 0, 0, 0)
cluster.info.update({'emitter_plane': emitter_plane,
'emitter_tag': emitter_tag,
'emitter': emitter,
'side': side,
'level': level,
'ke': ke})
clusters.append(cluster)
return clusters
def compute_polar_scans(clusters, theta=np.arange(-20., 80., 1.), phi=0., data=DATA):
for ic, cluster in enumerate(clusters):
# select the spectroscopy of the calculation and create a new folder for each calculation
side, emitter, tag, plane, level, ke = [cluster.info[k] for k in ('side', 'emitter', 'emitter_tag',
'emitter_plane', 'level', 'ke')]
calc = MSSPEC(spectroscopy='PED', folder='calc{:0>3d}_S{}_E{}_T{:d}_P{:d}'.format(ic, side, emitter, tag,
plane))
calc.calculation_parameters.scattering_order = max(1, min(4, plane))
calc.source_parameters.energy = XRaySource.AL_KALPHA
calc.source_parameters.theta = -35
calc.detector_parameters.angular_acceptance = 4.
calc.detector_parameters.average_sampling = 'medium'
calc.calculation_parameters.path_filtering = 'forward_scattering'
calc.calculation_parameters.off_cone_events = 1
[a.set('forward_angle', 30.) for a in cluster]
calc.calculation_parameters.vibrational_damping = 'averaged_tl'
[a.set('mean_square_vibration', 0.006) for a in cluster]
calc.set_atoms(cluster)
data = calc.get_theta_scan(level=level, theta=theta, phi=phi, kinetic_energy=ke, data=data)
dset = data[-1]
dset.title = 'Polar scan {emitter:s}({level:s} tag {tag:d}) in plane #{plane:d}'.format(emitter=emitter,
level=level, tag=tag,
plane=plane)
for name, value in cluster.info.items():
dset.add_parameter(group='AlnTuto', name=name, value=str(value), unit="")
#data.save('all_polar_scans.hdf5', append=True)
data.save('all_polar_scans.hdf5')
def analysis(filename='all_polar_scans.hdf5'):
data=Data.load(filename)
# create datasets to store the sum of all emitter
tmp_data = {}
for dset in data:
# retrieve some info
side = dset.get_parameter(group='AlnTuto', name='side')['value']
emitter = dset.get_parameter(group='AlnTuto', name='emitter')['value']
try:
key = '{}_{}'.format(side, emitter)
tmp_data[key] += dset.cross_section
except KeyError:
tmp_data[key] = dset.cross_section.copy()
# get the index of theta = 0;
it0 = np.where(dset.theta == 0)[0][0]
for key, cs in tmp_data.items():
tmp_data[key + '_norm'] = cs / cs[it0]
tmp_data['Al_side'] = tmp_data['Al_Al_norm'] / tmp_data['Al_N_norm']
tmp_data['N_side'] = tmp_data['N_Al_norm'] / tmp_data['N_N_norm']
# now add all columns
substrate_dset = data.add_dset('Total substrate signal')
substrate_dset.add_columns(theta=dset.theta.copy())
substrate_dset.add_columns(**tmp_data)
view = substrate_dset.add_view('Al side',
title=r'AlN Polar scan in the (10$\bar{1}$0) azimuthal plane - Al side polarity',
xlabel=r'$\Theta (\degree$)',
ylabel='Signal (a.u.)')
view.select('theta', 'Al_Al_norm', legend='Al 2p')
view.select('theta', 'Al_N_norm', legend='N 1s')
view.set_plot_options(autoscale=True)
view = substrate_dset.add_view('N side',
title=r'AlN Polar scan in the (10$\bar{1}$0) azimuthal plane - N side polarity',
xlabel=r'$\Theta (\degree$)',
ylabel='Signal (a.u.)')
view.select('theta', 'N_Al_norm', legend='Al 2p')
view.select('theta', 'N_N_norm', legend='N 1s')
view.set_plot_options(autoscale=True)
view = substrate_dset.add_view('Ratios',
title=r'Al(2p)/N(1s) ratios on both polar sides of AlN in the (10$\bar{1}$0) '
r'azimuthal plane',
xlabel=r'$\Theta (\degree$)',
ylabel='Intenisty ratio')
view.select('theta', 'Al_side', legend='Al side', where="theta >= 0 and theta <=70")
view.select('theta', 'N_side', legend='N side', where="theta >= 0 and theta <=70")
view.set_plot_options(autoscale=True)
data.save('analysis.hdf5')
data.view()
DIAMETER = 10
PLANES = 4
clusters = create_clusters(side='Al', emitter='Al', diameter=DIAMETER, planes=PLANES) + \
create_clusters(side='Al', emitter='N', diameter=DIAMETER, planes=PLANES) + \
create_clusters(side='N', emitter='Al', diameter=DIAMETER, planes=PLANES) + \
create_clusters(side='N', emitter='N', diameter=DIAMETER, planes=PLANES)
compute_polar_scans(clusters, phi=0.)
analysis()