epsi
/
msspec_python3
3
0
Fork 0
Code Issues 6 Pull Requests Releases Wiki Activity

Updated documentation and tutorials. 

Download and install section contains now information about
how to install MsSpec using Docker.
Tutorials have been updated a bit.
This commit is contained in:
Sylvain Tricot 2021-09-24 16:13:14 +02:00
parent 08c6ea7a34
commit b011e4348c
14 changed files with 598 additions and 575 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

148
doc/source/downloads.rst
View File

@ -1,3 +1,7 @@
.. |LINUXSCRIPT| replace:: https://git.ipr.univ-rennes1.fr/epsi/msspec_python3/raw/branch/devel/utils/dockerized/linux/msspec
########################## ##########################
Download and install notes Download and install notes
########################## ##########################
@ -9,36 +13,138 @@ Installation
There are 2 ways to install MsSpec. You can either: There are 2 ways to install MsSpec. You can either:
- Use a Docker image - Use a Docker image. This is, by far, the most straightforward and easy way to work with MsSpec.
- Compile your own version - Compile your own version for your system.
Using a Docker image
--------------------
You first need Docker (Docker Desktop) to be installed on your system. 1. Using a Docker image
This is really straightforward. Please see the `Docker website -----------------------
<https://www.docker.com>`_ for more information.
You first need `Docker <https://www.docker.com>`_ to be installed on your system.
This is really straightforward. Please see the `Docker documentation
<https://docs.docker.com/engine/install/>`_ for more information depending on
your OS.
- pull the image from our repository: * **For Linux**,
::
docker pull iprcnrs/msspec:latest * Download :download:`this script <../../utils/dockerized/linux/msspec>` and
make it executable (with the *chmod u+x msspec* command)
* Place it in a location known to your $PATH (or add this location to your $PATH if needed).
The *~/.local/bin* folder is a good choice.
That's all. MsSpec is installed. Now the command to launch the MsSpec environment depends slighly on the * **For Windows**
OS you are runing:
- For Linux * You need a running X server. Download and install `VcXsrv <https://sourceforge.net/projects/vcxsrv/>`_
::
xhost +
docker run --rm -it -e DISPLAY=$DISPLAY --net=host -v msspec_data:/home/msspec/data msspec
- For Windows * Install the small :download:`MsSpec frontend <../../utils/dockerized/windows/msspec_setup.exe>`
::
docker run --rm -it -e DISPLAY=host.internal.display:0 -v msspec_data:/home/msspec/data msspec * While not necessary, it is also a good idea to use a better terminal application than the default one.
`Windows Terminal <https://www.microsoft.com/fr-fr/p/windows-terminal/9n0dx20hk701#activetab=pivot:overviewtab>`_
may be a good choice.
* **For Mac**
* *To be documented*
Compile your own version Open a terminal in Linux, or a powershell in Windows and type in::
------------------------
msspec
The first time you run the command, it will download the msspec image (it may takes several minutes or half an hour
depending on your internet connexion).
The command will automatically create a new container and start it.
As the command was entered without any argument on the command-line, the help message should be printed on the screen
.. code-block:: text
Usage: 1) msspec -p [PYTHON OPTIONS] SCRIPT [ARGUMENTS...]
2) msspec [-l FILE | -i | -h]
3) msspec [bash | reset]
Form (1) is used to launch a script
Form (2) is used to load a hdf5 data file
Form (3) is used to control the Docker container/image.
List of possible options:
-p Pass every arguments after this option to the msspec
virtual environment Python interpreter.
-i Run the interactive Python interpreter within msspec
virtual environment.
-l Load and display a *.hdf5 data file in a graphical
window.
-v Print the version.
-h Show this help message.
bash This command starts an interactive bash shell in the
MsSpec container.
reset This command removes the MsSpec container (but not the
image). Changes made in the container will be lost and
any new call to msspec will recreate a new fresh container.
2. Compile your own version
---------------------------
To install MsSpec this way, follow the instructions `here <https://git.ipr.univ-rennes1.fr/epsi/msspec_python3/src/branch/devel>`_ To install MsSpec this way, follow the instructions `here <https://git.ipr.univ-rennes1.fr/epsi/msspec_python3/src/branch/devel>`_
***************************
Running your Python scripts
***************************
You can run your MsSpec Python scripts (e.g. *my_script.py*) by typing in::
msspec -p my_script.py
This command is equivalent to activating the Python virtual environment MsSpec is intsalled in
and to run the Python interpreter of that environment (everything that follows the -p option is
passed to the python command).
You can also launch the Interactive Python (iPython) of MsSpec::
msspec -i
Inside this interactive session, you can run a script with the command::
%run my_script.py
You can interact with your filesystem with the classical *cd*, *ls*, *cp*, *rm*... commands.
and you can edit your script with::
%ed my_script.py
.. warning::
**If using the Docker image of MsSpec in Linux**, your home folder on the host machine is bind mounted
in the same location in the container and your UID and GID are also set so that creating files within
your home file hierarchy is totally transparent.
**If using the Docker image of MsSpec in Windows**, the drive containing your "My Documents" folder on the
host machine is bind mounted on the container at the root of the filesystem. For example, if your
"My documents" folder on Windows are in 'D:\\Home\\Bob\\MyDocuments', it will be available in the container as
'/D/Home/Bob/MyDocuments'. It has two consequences:
#. The msspec command will fail if running on another drive than the one where is located "My Documents"
#. You have to specify filenames with the Unix slashes. For example if you want to run the script located
in *.\\results\\my_script.py*, you will have to enter *msspec -p ./results/my_script.py*
**************
Uninstallation
**************
* **Under Linux**, type in::
msspec -u
* **Under Windows**, simply `uninstall the application from the Settings page
<https://support.microsoft.com/en-us/windows/uninstall-or-remove-apps-and-programs-in-windows-10-4b55f974-2cc6-2d2b-d092-5905080eaf98>`_.
A command window will pop-up and you will have to answer 'y' to remove MsSpec.
* **Under Mac OS**, *to be documented.*

6
doc/source/modules/config.rst
View File

@ -1,6 +0,0 @@
:orphan:
.. automodule:: config
:members:
:undoc-members:
:show-inheritance:

6
doc/source/modules/looper.rst Normal file
View File

@ -0,0 +1,6 @@
:orphan:
.. automodule:: looper
:members:
:undoc-members:
:show-inheritance:

6
doc/source/modules/version.rst Normal file
View File

@ -0,0 +1,6 @@
:orphan:
.. automodule:: version
:members:
:undoc-members:
:show-inheritance:

311
doc/source/tutorials/AlN/AlN.py
View File

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

79
doc/source/tutorials/RhO/RhO.py
View File

@ -1,54 +1,25 @@
# -*- encoding: utf-8 -*- # coding: utf8
# vim: set fdm=indent ts=4 sw=4 sts=4 et ai tw=80 cc=+0 mouse=a nu : #
from msspec.calculator import MSSPEC
from msspec.calculator import MSSPEC, XRaySource from ase.build import fcc111, add_adsorbate
from msspec.utils import * import numpy as np
from ase.build import fcc111, add_adsorbate data = None
from ase.visualize import view all_z = np.arange(1.10, 1.65, 0.05)
from msspec.iodata import cols2matrix for zi, z0 in enumerate(all_z):
# construct the cluster
from matplotlib import pyplot as plt cluster = fcc111('Rh', size = (2,2,1))
import numpy as np cluster.pop(3)
import sys add_adsorbate(cluster, 'O', z0, position = 'fcc')
cluster.absorber = len(cluster) - 1
data = None
all_z = np.arange(1.10, 1.50, 0.02) # Define a calculator
all_z=(1.1,) calc = MSSPEC(spectroscopy='PED', algorithm='inversion')
for zi, z0 in enumerate(all_z): calc.set_atoms(cluster)
# construct the cluster
cluster = fcc111('Rh', size = (2,2,1)) # Compute
add_adsorbate(cluster, 'O', z0, position = 'fcc') data = calc.get_theta_phi_scan(level='1s', kinetic_energy=723, data=data)
cluster.pop(3) dset = data[-1]
#cluster.rotate('z',np.pi/3.) dset.title = "{:d}) z = {:.2f} angstroms".format(zi, z0)
#view(cluster)
data.view()
cluster.absorber = len(cluster) - 1
calc = MSSPEC(spectroscopy = 'PED', folder = './RhO_z')
calc.set_atoms(cluster)
calc.calculation_parameters.scattering_order = 3
calc.calculation_parameters.RA_cutoff = 1
calc.source_parameters.energy = XRaySource.AL_KALPHA
# compute
level = '1s'
ke = 723.
data = calc.get_theta_phi_scan(level=level, kinetic_energy=ke, data=data)
# OPTIONAL, this will create an image of the data that you can combine
# in an animated gif
dset = data[-1]
theta, phi, Xsec = cols2matrix(dset.theta, dset.phi, dset.cross_section)
X, Y = np.meshgrid(np.radians(phi), 2*np.tan(np.radians(theta/2.)))
fig = plt.figure()
ax = fig.add_subplot(111, projection='polar')
im = ax.pcolormesh(X, Y, Xsec)
theta_ticks = np.arange(0, 91, 15)
plt.yticks(2 * np.tan(np.radians(theta_ticks/2.)), theta_ticks)
plt.title(r"$z_0 = {:.2f} \AA$".format(z0))
plt.savefig('image{:03d}.png'.format(zi))
data.view()

34
doc/source/tutorials/RhO/sf.py
View File

@ -1,31 +1,17 @@
# -*- encoding: utf-8 -*- # coding: utf8
# vim: set fdm=indent ts=4 sw=4 sts=4 et ai tw=80 cc=+0 mouse=a nu : #
from msspec.calculator import MSSPEC, XRaySource
from msspec.utils import *
from msspec.calculator import MSSPEC
from ase import Atoms from ase import Atoms
import numpy as np # Create an atomic chain O-Rh
cluster = Atoms(['O', 'Rh'], positions = [(0,0,0), (0,0,4.)])
# Defining global variables # Create the calculator
a0 = 6.0 calc = MSSPEC(spectroscopy = 'PED')
symbols = ('Rh', 'O') calc.set_atoms(cluster)
ke = 723. cluster.absorber = 0
level = '1s'
data = None # compute
for symbol in symbols: data = calc.get_scattering_factors(level='1s', kinetic_energy=723)
cluster = Atoms(symbol*2, positions = [(0,0,0), (0,0,a0)])
[a.set('mt_radius', 1.5) for a in cluster]
# Create the calculator
calc = MSSPEC(spectroscopy = 'PED')
calc.source_parameters.energy = XRaySource.AL_KALPHA
calc.set_atoms(cluster)
cluster.absorber = 0
# compute
data = calc.get_scattering_factors(level=level, kinetic_energy=ke, data=data)
data.view() data.view()

14
doc/source/tutorials/RhO/tuto_ped_RhO.rst
View File

@ -82,13 +82,13 @@ Here is the script for the computation. (:download:`download <RhO.py>`)
.. literalinclude:: RhO.py .. literalinclude:: RhO.py
:linenos: :linenos:
.. note:: .. .. note::
After runing this script, you will get 20 images in your folder. You can merge them in one animated gif image .. After runing this script, you will get 20 images in your folder. You can merge them in one animated gif image
like this: .. like this:
..
.. code-block:: bash .. .. code-block:: bash
..
convert -delay 50 -loop 0 image*.png animation.gif .. convert -delay 50 -loop 0 image*.png animation.gif
.. seealso:: .. seealso::

23
doc/source/tutorials/copper/Cu_simple.py → doc/source/tutorials/copper/Cu.py
View File

@ -1,7 +1,7 @@
# coding: utf-8 # coding: utf8
from msspec.calculator import MSSPEC from msspec.calculator import MSSPEC
from msspec.utils import * from msspec.utils import hemispherical_cluster, get_atom_index
from ase.build import bulk from ase.build import bulk
from ase.visualize import view from ase.visualize import view
@ -9,20 +9,15 @@ from ase.visualize import view
a0 = 3.6 # The lattice parameter in angstroms a0 = 3.6 # The lattice parameter in angstroms
# Create the copper cubic cell # Create the copper cubic cell
cluster = bulk('Cu', a=a0, cubic=True) copper = bulk('Cu', a=a0, cubic=True)
# Repeat the cell many times along x, y and z cluster = hemispherical_cluster(copper, planes=3, emitter_plane=2)
cluster = cluster.repeat((4, 4, 4))
# Put the center of the structure at the origin
center_cluster(cluster)
# Keep atoms inside a given radius
cluster = cut_sphere(cluster, radius=a0 + .01)
# Keep only atoms below the plane z <= 0
cluster = cut_plane(cluster, z=0.1)
# Set the absorber (the deepest atom centered in the xy-plane) # Set the absorber (the deepest atom centered in the xy-plane)
cluster.absorber = get_atom_index(cluster, 0, 0, -a0) cluster.absorber = get_atom_index(cluster, 0, 0, 0)
# Create a calculator for the PhotoElectron Diffration # Create a calculator for the PhotoElectron Diffration
calc = MSSPEC(spectroscopy='PED') calc = MSSPEC(spectroscopy='PED')
# Set the cluster to use for the calculation # Set the cluster to use for the calculation
calc.set_atoms(cluster) calc.set_atoms(cluster)
@ -31,4 +26,6 @@ data = calc.get_theta_scan(level='2p3/2')
# Show the results # Show the results
data.view() data.view()
#calc.shutdown()
# Clean temporary files
calc.shutdown()

91
doc/source/tutorials/copper/tuto_ped_copper.rst
View File

@ -16,7 +16,7 @@ Begin by typing:
:linenos: :linenos:
:lineno-start: 1 :lineno-start: 1
# coding: utf-8 # coding: utf8
from msspec.calculator import MSSPEC from msspec.calculator import MSSPEC
from msspec.utils import * from msspec.utils import *
@ -45,27 +45,37 @@ This is the job of the *ase* module, so load the module
from ase.build import bulk from ase.build import bulk
from ase.visualize import view from ase.visualize import view
a0 = 3.6 a0 = 3.6 # The lattice parameter in angstroms
cluster = bulk('Cu', a = a0, cubic = True)
# Create the copper cubic cell
copper = bulk('Cu', a=a0, cubic=True)
view(cluster) view(cluster)
In line 6 we load the :py:func:`bulk` function to create our atomic object and, In line 6 we load the :py:func:`bulk` function to create our atomic object and,
in line 7, we load the :py:func:`view` function to actually view our cluster. in line 7, we load the :py:func:`view` function to actually view our cluster.
The creation of the cluster is done line 10. The :py:func:`bulk` needs one The creation of the cluster starts on line 12. We create first the copper cubic
argument which are the chemical species ('Cu' in our example). We also pass 2 cell with the The :py:func:`bulk` function. It needs one argument which are the
keyword (optional) arguments here: chemical species ('Cu' in our example). We also pass 2 keyword (optional) arguments
here:
* The lattice parameter *a* in units of angströms. * The lattice parameter *a* in units of angströms.
* The *cubic* keyword, to obtain a cubic cell rather than the fully * The *cubic* keyword, to obtain a cubic cell rather than the fully
reduced one which is not cubic reduced one which is not cubic
From now on, you can save your script as 'Cu.py' and open a terminal window in From now on, you can save your script as 'Cu.py' and open a terminal window in
the same directory as this file. Launch your script using python: the same directory as this file. Launch your script using 'python' or 'msspec -p'
(depending on your installation):
.. code-block:: bash .. code-block:: bash
python2 Cu.py python Cu.py
or
.. code-block:: bash
msspec -p Cu.py
and a graphical window (the ase-gui) should open with a cubic cell of copper and a graphical window (the ase-gui) should open with a cubic cell of copper
like this one: like this one:
@ -80,34 +90,38 @@ like this one:
Obviously, multiple scattering calculations need more atoms to be accurate. Due Obviously, multiple scattering calculations need more atoms to be accurate. Due
to the forward focusing effect in photo-electron diffraction, the best suited to the forward focusing effect in photo-electron diffraction, the best suited
geometry for the cluster is hemispherical. Obtaining such a cluster is a geometry for the cluster is hemispherical. Obtaining such a cluster is a
straightforward process: straightforward process thanks to the :py:func:`utils.hemispherical_cluster` function.
This function will basically create a cluster based on a pattern (the cubic copper
1. Repeat the previous cell in all 3 directions cell here).
2. center the structure
3. Keep only atoms within a given radius from center
4. Keep only atoms within one halh of the sphere.
Modify your script like this and run it again. Modify your script like this and run it again.
.. literalinclude:: Cu_simple.py .. literalinclude:: Cu.py
:linenos: :linenos:
:lineno-start: 1 :lineno-start: 1
:lines: 1-21 :lines: 1-13
Don't forget to add the line to view the cluster at the end of the script and run Don't forget to add the line to view the cluster at the end of the script and run
the script again. the script again. The :py:func:`hemispherical_cluster` works in 3 simple steps:
#. Repeat the given *pattern* in all 3 directions
#. Center this new set of atoms and cut a sphere from the center
#. Remove the upper half of the created 'sphere'.
To get more information about how to use this function, have a look at the :ref:`hemispherical_cluster_faq` section in the :ref:`faq`.
.. figure:: Cu_fig2.png .. figure:: Cu_fig2.png
:align: center :align: center
:width: 60% :width: 60%
Figure 2. The different steps to output a cluster. Figure 2. The different steps to output a cluster.
a) After repeat, b) after cut_sphere, c) after cut_plane a) After repeat, b) after cutting a sphere, c) final cluster
Once you a cluster is built the next step is to define which atom in the cluster Once your cluster is built the next step is to define which atom in the cluster
will absorbe the light. This atom is called the *absorber*. will absorb the light. This atom is called the *absorber* (or the *emitter* since
it emits the photoelectron).
To specify which atom is the absorber, you need to understand that the cluster To specify which atom is the absorber, you need to understand that the cluster
object is like a list of atoms. Each member of this list is an atom with its object is like a list of atoms. Each member of this list is an atom with its
@ -115,7 +129,7 @@ own position. You need to locate the index of the atom in the list that you want
it to be the absorber. Then, put that number in the *cluster.absorber* attribute it to be the absorber. Then, put that number in the *cluster.absorber* attribute
For example, suppose that you want the first atom of the list to be the For example, suppose that you want the first atom of the list to be the
absorber. You should write: absorber. You would write:
.. code-block:: python .. code-block:: python
@ -123,24 +137,27 @@ absorber. You should write:
To find what is the index of the atom you'd like to be the absorber, you can To find what is the index of the atom you'd like to be the absorber, you can
either get it while you are visualizing the structure within the ase-gui either get it while you are visualizing the structure within the ase-gui
program. Or, you can use :py:func:`get_atom_index` function. This function takes program (select an atom with the left mouse button and look at its index in the
4 arguments: the cluster to look the index for, and the x, y and z coordinates. status line). Or, you can use :py:func:`utils.get_atom_index` function. This function
takes 4 arguments: the cluster to look the index for, and the x, y and z coordinates.
It will return the index of the closest atom to these coordinates. In our It will return the index of the closest atom to these coordinates. In our
example, to get the deepest atom centered in the xy-plane, we write: example, as we used the :py:func:`utils.hemispherical_cluster` function to create our
cluster, the *emitter* (*absorber*) is always located at the origin, so defining it
is straightforward:
.. literalinclude:: Cu_simple.py .. literalinclude:: Cu.py
:linenos: :linenos:
:lineno-start: 22 :lineno-start: 15
:lines: 22-23 :lines: 15-16
That's all for the cluster part. We now need to create a calculator for that That's all for the cluster part. We now need to create a calculator for that
object. This is a 2 lines procedure: object. This is a 2 lines procedure:
.. literalinclude:: Cu_simple.py .. literalinclude:: Cu.py
:linenos: :linenos:
:lineno-start: 24 :lineno-start: 18
:lines: 24-28 :lines: 18-22
When creating a new calculator, you must choose the kind of spectroscopy you When creating a new calculator, you must choose the kind of spectroscopy you
will work with. In this example we choose 'PED' for *PhotoElectron Diffraction*. will work with. In this example we choose 'PED' for *PhotoElectron Diffraction*.
@ -153,14 +170,14 @@ Other types of spectroscopies are:
Now that everything is ready you can actually perform a calculation. The 2 lines Now that everything is ready you can actually perform a calculation. The lines
below will produce a polar scan of the Cu(2p3/2) level with default parameters, below will produce a polar scan of the Cu(2p3/2) level with default parameters,
store the results in the data object and display it in a graphical window. store the results in the data object and display it in a graphical window.
.. literalinclude:: Cu_simple.py .. literalinclude:: Cu.py
:linenos: :linenos:
:lineno-start: 29 :lineno-start: 24
:lines: 29-33 :lines: 24-28
running this script will produce the following output running this script will produce the following output
@ -177,7 +194,6 @@ By default, the program computes for :math:`\theta` angles in the -70°..+70°
range. This can be changed by using the *angles* keyword. range. This can be changed by using the *angles* keyword.
.. code-block:: python .. code-block:: python
:linenos:
#For a polar scan between 0° and 60° with 100 points #For a polar scan between 0° and 60° with 100 points
import numpy as np import numpy as np
@ -201,7 +217,6 @@ work function of the sample, arbitrary set to 4.5 eV.
Of course, you can choose any kinetic energy you'd like: Of course, you can choose any kinetic energy you'd like:
.. code-block:: python .. code-block:: python
:linenos:
# To set the kinetic energy... # To set the kinetic energy...
data = calc.get_theta_scan(level = '2p3/2', kinetic_energy = 300) data = calc.get_theta_scan(level = '2p3/2', kinetic_energy = 300)
@ -209,9 +224,9 @@ Of course, you can choose any kinetic energy you'd like:
Below is the full code of this example. You can download it :download:`here Below is the full code of this example. You can download it :download:`here
<Cu_simple.py>` <Cu.py>`
.. literalinclude:: Cu_simple.py .. literalinclude:: Cu.py
:linenos: :linenos:
.. seealso:: .. seealso::

16
doc/source/tutorials/index.rst Normal file
View File

@ -0,0 +1,16 @@
.. _tutorials:
####################
Learn from tutorials
####################
...About Photodiffraction
=========================
.. toctree::
copper/tuto_ped_copper
nickel_chain/tuto_ped_nickel_chain
RhO/tuto_ped_RhO
temperature/tuto_ped_temperature
AlN/tuto_cluster

155
doc/source/tutorials/nickel_chain/Ni_chain.py
View File

@ -1,77 +1,78 @@
# coding: utf-8 # coding: utf-8
# import all we need and start by msspec # import all we need and start by msspec
from msspec.calculator import MSSPEC from msspec.calculator import MSSPEC
# we will build a simple atomic chain # we will build a simple atomic chain
from ase import Atom, Atoms from ase import Atom, Atoms
# we need some numpy functions # we need some numpy functions
import numpy as np import numpy as np
symbol = 'Ni' # The kind of atom for the chain symbol = 'Ni' # The kind of atom for the chain
orders = (1, 5) # We will run the calculation for single scattering orders = (1, 5) # We will run the calculation for single scattering
# and multiple scattering (5th diffusion order) # and multiple scattering (5th diffusion order)
chain_lengths = (2,3,5) # We will run the calculation for differnt lengths chain_lengths = (2,3,5) # We will run the calculation for differnt lengths
# of the atomic chain # of the atomic chain
a = 3.5 # The distance bewteen 2 atoms a = 3.499 * np.sqrt(2)/2 # The distance bewteen 2 atoms
# Define an empty variable to store all the results # Define an empty variable to store all the results
all_data = None all_data = None
# 2 for nested loops over the chain length and the order of diffusion # 2 for nested loops over the chain length and the order of diffusion
for chain_length in chain_lengths: for chain_length in chain_lengths:
for order in orders: for order in orders:
# We build the atomic chain by # We build the atomic chain by
# 1- stacking each atom one by one along the z axis # 1- stacking each atom one by one along the z axis
chain = Atoms([Atom(symbol, position = (0., 0., i*a)) for i in chain = Atoms([Atom(symbol, position = (0., 0., i*a)) for i in
range(chain_length)]) range(chain_length)])
# 2- rotating the chain by 45 degrees with respect to the y axis # 2- rotating the chain by 45 degrees with respect to the y axis
chain.rotate('y', np.radians(45.)) #chain.rotate('y', np.radians(45.))
# 3- setting a custom Muffin-tin radius of 1.5 angstroms for all chain.rotate(45., 'y')
# atoms (needed if you want to enlarge the distance between # 3- setting a custom Muffin-tin radius of 1.5 angstroms for all
# the atoms while keeping the radius constant) # atoms (needed if you want to enlarge the distance between
[atom.set('mt_radius', 1.5) for atom in chain] # the atoms while keeping the radius constant)
# 4- defining the absorber to be the first atom in the chain at #[atom.set('mt_radius', 1.5) for atom in chain]
# x = y = z = 0 # 4- defining the absorber to be the first atom in the chain at
chain.absorber = 0 # x = y = z = 0
chain.absorber = 0
# We define a new PED calculator
calc = MSSPEC(spectroscopy = 'PED') # We define a new PED calculator
calc.set_atoms(chain) calc = MSSPEC(spectroscopy = 'PED')
# Here is how to tweak the scattering order calc.set_atoms(chain)
calc.calculation_parameters.scattering_order = order # Here is how to tweak the scattering order
# This line below is where we actually run the calculation calc.calculation_parameters.scattering_order = order
all_data = calc.get_theta_scan(level='3s', kinetic_energy=1000., # This line below is where we actually run the calculation
theta=np.arange(0., 80.), data=all_data) all_data = calc.get_theta_scan(level='3s', #kinetic_energy=1000.,
theta=np.arange(0., 80.), data=all_data)
# OPTIONAL, to improve the display of the data we will change the dataset
# default title as well as the plot title # OPTIONAL, to improve the display of the data we will change the dataset
t = "order {:d}, n = {:d}".format(order, chain_length) # A useful title # default title as well as the plot title
dset = all_data[-1] # get the last dataset t = "order {:d}, n = {:d}".format(order, chain_length) # A useful title
dset.title = t # change its title dset = all_data[-1] # get the last dataset
# get its last view (there is only one defined for each dataset) dset.title = t # change its title
v = dset.views()[-1] # get its last view (there is only one defined for each dataset)
v.set_plot_options(title=t) # change the title of the figure v = dset.views()[-1]
v.set_plot_options(title=t) # change the title of the figure
# OPTIONAL, set the same scale for all plots
# 1. iterate over all computed cross_sections to find the absolute minimum and # OPTIONAL, set the same scale for all plots
# maximum of the data # 1. iterate over all computed cross_sections to find the absolute minimum and
min_cs = max_cs = 0 # maximum of the data
for dset in all_data: min_cs = max_cs = 0
min_cs = min(min_cs, np.min(dset.cross_section)) for dset in all_data:
max_cs = max(max_cs, np.max(dset.cross_section)) min_cs = min(min_cs, np.min(dset.cross_section))
max_cs = max(max_cs, np.max(dset.cross_section))
# 2. for each view in each dataset, change the scale accordingly
for dset in all_data: # 2. for each view in each dataset, change the scale accordingly
v = dset.views()[-1] for dset in all_data:
v.set_plot_options(ylim=[min_cs, max_cs]) v = dset.views()[-1]
v.set_plot_options(ylim=[min_cs, max_cs])
# Pop up the graphical window
all_data.view() # Pop up the graphical window
# You can end your script with the line below to remove the temporary all_data.view()
# folder needed for the calculation # You can end your script with the line below to remove the temporary
calc.shutdown() # folder needed for the calculation
calc.shutdown()

2
doc/source/tutorials/nickel_chain/tuto_ped_nickel_chain.rst
View File

@ -30,7 +30,7 @@ trajectory for this electron which in turn lowers the signal.
:align: center :align: center
:width: 80% :width: 80%
Figure 4. Polar scan of a Ni chain of 2-5 atoms for single and mutliple (5th order) Figure 1. Polar scan of a Ni chain of 2-5 atoms for single and mutliple (5th order)
scattering. scattering.

282
doc/source/tutorials/temperature/Cu_temperature.py
View File

@ -1,153 +1,129 @@
# coding: utf-8 # coding: utf8
from msspec.calculator import MSSPEC, XRaySource from ase.build import bulk
from msspec.utils import * import numpy as np
from msspec.iodata import Data from msspec.calculator import MSSPEC, XRaySource
from ase.build import bulk from msspec.utils import hemispherical_cluster, get_atom_index
import numpy as np
import sys def create_clusters(nplanes=3):
copper = bulk('Cu', a=3.6, cubic=True)
a0 = 3.6 # The lattice parameter in angstroms clusters = []
phi = np.arange(0, 100) # An array of phi angles for emitter_plane in range(nplanes):
all_T = np.arange(300, 1000, 100) # All temperatures used for the calculation cluster = hemispherical_cluster(copper,
all_theta = np.array([45, 83]) # 2 polar angles, 83° is grazing detection emitter_plane=emitter_plane,
eps = 0.01 # a useful small value planes=emitter_plane+2,
shape='cylindrical')
DATA_FNAME = 'all_data.hdf5' # Where to store all the data cluster.absorber = get_atom_index(cluster, 0, 0, 0)
ANALYSIS_FNAME = 'analysis.hdf5' cluster.info.update({
'emitter_plane': emitter_plane,
})
def compute(filename): clusters.append(cluster)
"""Will compute a phi scan for an emitter in the first, second and third plane return clusters
of a copper substrate for various polar angles and temperatures.
In a second step (outside this function), intensities from all these emitters def compute(clusters, all_theta=[45., 83.],
are added to get the signal of a substrate.""" all_T=np.arange(300., 1000., 400.)):
calc = MSSPEC(spectroscopy='PED') data = None
calc.source_parameters.energy = XRaySource.AL_KALPHA for ic, cluster in enumerate(clusters):
calc.muffintin_parameters.interstitial_potential = 14.1 # Retrieve info from cluster object
plane = cluster.info['emitter_plane']
calc.calculation_parameters.vibrational_damping = 'averaged_tl'
calc.calculation_parameters.use_debye_model = True calc = MSSPEC(spectroscopy='PED', algorithm='expansion')
calc.calculation_parameters.debye_temperature = 343 calc.source_parameters.energy = XRaySource.AL_KALPHA
calc.calculation_parameters.vibration_scaling = 1.2 calc.muffintin_parameters.interstitial_potential = 14.1
calc.detector_parameters.average_sampling = 'high' calc.calculation_parameters.vibrational_damping = 'averaged_tl'
calc.detector_parameters.angular_acceptance = 5.7 calc.calculation_parameters.use_debye_model = True
calc.calculation_parameters.debye_temperature = 343
filters = ['forward_scattering', 'backward_scattering'] calc.calculation_parameters.vibration_scaling = 1.2
calc.calculation_parameters.path_filtering = filters
calc.detector_parameters.average_sampling = 'high'
calc.calculation_parameters.RA_cutoff = 3 calc.detector_parameters.angular_acceptance = 5.7
# You can also choose a real potential and manually define the for atom in cluster:
# electron mean free path atom.set('forward_angle', 30)
# atom.set('backward_angle', 30)
# calc.tmatrix_parameters.exchange_correlation = 'hedin_lundqvist_real' filters = ['forward_scattering', 'backward_scattering']
# calc.calculation_parameters.mean_free_path = 10. calc.calculation_parameters.path_filtering = filters
data = None # init an empty data object calc.calculation_parameters.RA_cutoff = 2
for temperature in all_T:
for plane in range(3): for T in all_T:
# We create a cylindrical cluster here with one plane below the emitter for theta in all_theta:
# and 0, 1 or to planes above the emitter calc.calculation_parameters.temperature = T
cluster = bulk('Cu', a=a0, cubic=True) calc.calculation_parameters.scattering_order = min(1 + plane, 3)
cluster = cluster.repeat((20, 20, 8)) calc.set_atoms(cluster)
center_cluster(cluster) data = calc.get_phi_scan(level='2p3/2', theta=theta,
cluster = cut_cylinder(cluster, radius=1.5 * a0 + eps) phi=np.linspace(0, 100),
cluster = cut_plane(cluster, z=-( a0/2 + eps)) kinetic_energy=560, data=data)
cluster = cut_plane(cluster, z=(plane) * a0 / 2 + eps) dset = data[-1]
cluster.absorber = get_atom_index(cluster, 0, 0, 0) dset.title = "plane #{:d}, T={:f}K, theta={:f}°".format(plane,
T,
calc.calculation_parameters.temperature = temperature theta)
# the scattering order depends on the number of planes above
# the emitter to speed up calculations dset.add_parameter(group='misc', name='plane', value=plane, unit='')
calc.calculation_parameters.scattering_order = 1 + plane dset.add_parameter(group='misc', name='T', value=T, unit='')
dset.add_parameter(group='misc', name='theta', value=theta, unit='')
calc.set_atoms(cluster) return data
# Here starts the calculation
data = calc.get_phi_scan(level='2p3/2', theta=all_theta, phi=phi,
kinetic_energy=560, data=data) def analysis(data):
all_plane = []
data.save(filename) all_T = []
all_theta = []
def process_data(datafile, outputfile): for dset in data:
"""Will create another Data file with some phi-scans corresponding to 3 plane = dset.get_parameter('misc', 'plane')['value']
planes at different temperatures and the anisotropy curve for 45° and T = dset.get_parameter('misc', 'T')['value']
grazing detection. theta = dset.get_parameter('misc', 'theta')['value']
""" cs = dset.cross_section.copy()
def get_signal(datafile, T=300, theta=45): phi = dset.phi.copy()
data = Data.load(datafile)
total = None if plane not in all_plane: all_plane.append(plane)
for dset in data: if T not in all_T: all_T.append(T)
p = {_['group'] + '.' + _['name']: _['value'] for _ in dset.parameters()} if theta not in all_theta: all_theta.append(theta)
temperature = np.float(p['CalculationParameters.temperature'])
if temperature != T: continue def get_anisotropy(theta, T):
i = np.where(dset.plane == -1) cs = None
j = np.where(dset.theta[i] == theta) for dset in data:
signal = dset.cross_section[i][j] try:
try: _T = dset.get_parameter('misc', 'T')['value']
total += signal _theta = dset.get_parameter('misc', 'theta')['value']
except: _cs = dset.cross_section.copy()
total = signal phi = dset.phi.copy()
phi = dset.phi[i][j] except:
return phi, total continue
if _theta == theta and _T == T:
analysis = Data('Temperature tutorial') try:
scans_dset = analysis.add_dset('Phi scans') cs += _cs
scans_view = scans_dset.add_view('Figure 1', except:
title=r'Cu(2p) Azimuthal scan for $\theta = 83\degree$', cs = _cs
xlabel=r'$\Phi (\degree$)', Imax = np.max(cs)
ylabel='Signal (a.u.)') Imin = np.min(cs)
anisotropy_dset = analysis.add_dset('Anisotropies') return (Imax - Imin)/Imax
anisotropy_view = anisotropy_dset.add_view('Figure 2',
title='Relative anisotropies for Cu(2p)', # create a substrate dataset for each T and theta
marker='o', anisotropy_dset = data.add_dset("all")
xlabel='T (K)', anisotropy_view = anisotropy_dset.add_view('Anisotropies',
ylabel=r'$\frac{\Delta I / I_{max}(T)}{\Delta I_{300} / I_{max}(300)} (\%)$') title='Relative anisotropies for Cu(2p)',
for theta in all_theta: marker='o',
for T in all_T: xlabel='T (K)',
PHI, SIGNAL = get_signal(datafile, T=T, theta=theta) ylabel=r'$\frac{\Delta I / I_{max}(T)}{\Delta I_{300}'
for phi, signal in zip(PHI, SIGNAL): r'/ I_{max}(300)} (\%)$')
scans_dset.add_row(phi=phi, signal=signal, theta=theta, temperature=T)
for theta in all_theta:
middleT = all_T[np.abs(all_T - np.mean(all_T)).argmin()] for T in all_T:
if theta == 83 and T in [np.min(all_T), middleT, np.max(all_T)]: A = get_anisotropy(theta, T)
scans_view.select('phi', 'signal', A0 = get_anisotropy(theta, np.min(all_T))
where='temperature == {:f} and theta == {:f}'.format(T, theta),
legend='{:.0f} K'.format(T)) anisotropy_dset.add_row(temperature=T, theta=theta, anisotropy=A/A0)
PHI, SIGNAL0 = get_signal(datafile, T=np.min(all_T), theta=theta) anisotropy_view.select('temperature', 'anisotropy',
Imax = np.max(SIGNAL) where='theta == {:f}'.format(theta),
Imax0 = np.max(SIGNAL0) legend=r'$\theta = {:.0f} \degree$'.format(theta))
dI = Imax - np.min(SIGNAL) return data
dI0 = Imax0 - np.min(SIGNAL0)
ani = (dI / Imax) / (dI0 / Imax0)
anisotropy_dset.add_row(temperature=T, anisotropy=ani, theta=theta) clusters = create_clusters()
data = compute(clusters)
anisotropy_view.select('temperature', 'anisotropy', data = analysis(data)
where='theta == {:f}'.format(theta), data.view()
legend=r'$\theta = {:.0f} \degree$'.format(theta))
analysis.save(outputfile)
# A convenient way to run the script, just specify one or more of these calc, analyse or
# view keywords as arguments
# ... to calculate all the data, analyse the data and view the results, run
# python Cu_temperature.py calc analyse view
# ... to just view the results, simply run
# python Cu_temperature.py view
if 'calc' in sys.argv:
compute(DATA_FNAME)
if 'analyse' in sys.argv:
process_data(DATA_FNAME, ANALYSIS_FNAME)
if 'view' in sys.argv:
data = Data.load(ANALYSIS_FNAME)
data.view()