From b011e4348c7b4c6e8ef6a57b9472f531566b6610 Mon Sep 17 00:00:00 2001 From: sylvain tricot Date: Fri, 24 Sep 2021 16:13:14 +0200 Subject: [PATCH] Updated documentation and tutorials. Download and install section contains now information about how to install MsSpec using Docker. Tutorials have been updated a bit. --- doc/source/downloads.rst | 148 +++++++-- doc/source/modules/config.rst | 6 - doc/source/modules/looper.rst | 6 + doc/source/modules/version.rst | 6 + doc/source/tutorials/AlN/AlN.py | 311 ++++++++---------- doc/source/tutorials/RhO/RhO.py | 79 ++--- doc/source/tutorials/RhO/sf.py | 34 +- doc/source/tutorials/RhO/tuto_ped_RhO.rst | 14 +- .../tutorials/copper/{Cu_simple.py => Cu.py} | 23 +- .../tutorials/copper/tuto_ped_copper.rst | 91 ++--- doc/source/tutorials/index.rst | 16 + doc/source/tutorials/nickel_chain/Ni_chain.py | 155 ++++----- .../nickel_chain/tuto_ped_nickel_chain.rst | 2 +- .../tutorials/temperature/Cu_temperature.py | 282 ++++++++-------- 14 files changed, 598 insertions(+), 575 deletions(-) delete mode 100644 doc/source/modules/config.rst create mode 100644 doc/source/modules/looper.rst create mode 100644 doc/source/modules/version.rst rename doc/source/tutorials/copper/{Cu_simple.py => Cu.py} (51%) create mode 100644 doc/source/tutorials/index.rst diff --git a/doc/source/downloads.rst b/doc/source/downloads.rst index a287be1..2fb51da 100644 --- a/doc/source/downloads.rst +++ b/doc/source/downloads.rst @@ -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 ########################## @@ -9,36 +13,138 @@ Installation There are 2 ways to install MsSpec. You can either: - - Use a Docker image - - Compile your own version + - Use a Docker image. This is, by far, the most straightforward and easy way to work with MsSpec. + - Compile your own version for your system. -Using a Docker image --------------------- -You first need Docker (Docker Desktop) to be installed on your system. -This is really straightforward. Please see the `Docker website -`_ for more information. +1. Using a Docker image +----------------------- + +You first need `Docker `_ to be installed on your system. +This is really straightforward. Please see the `Docker documentation +`_ for more information depending on +your OS. -- pull the image from our repository: - :: - docker pull iprcnrs/msspec:latest +* **For Linux**, + + * 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 -OS you are runing: +* **For Windows** -- For Linux - :: - xhost + - docker run --rm -it -e DISPLAY=$DISPLAY --net=host -v msspec_data:/home/msspec/data msspec + * You need a running X server. Download and install `VcXsrv `_ -- For Windows - :: - docker run --rm -it -e DISPLAY=host.internal.display:0 -v msspec_data:/home/msspec/data msspec + * Install the small :download:`MsSpec frontend <../../utils/dockerized/windows/msspec_setup.exe>` + + * While not necessary, it is also a good idea to use a better terminal application than the default one. + `Windows Terminal `_ + 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 `_ + + +*************************** +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 + `_. + A command window will pop-up and you will have to answer 'y' to remove MsSpec. + +* **Under Mac OS**, *to be documented.* diff --git a/doc/source/modules/config.rst b/doc/source/modules/config.rst deleted file mode 100644 index 7ae9e81..0000000 --- a/doc/source/modules/config.rst +++ /dev/null @@ -1,6 +0,0 @@ -:orphan: - -.. automodule:: config - :members: - :undoc-members: - :show-inheritance: \ No newline at end of file diff --git a/doc/source/modules/looper.rst b/doc/source/modules/looper.rst new file mode 100644 index 0000000..ad049c4 --- /dev/null +++ b/doc/source/modules/looper.rst @@ -0,0 +1,6 @@ +:orphan: + +.. automodule:: looper + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/source/modules/version.rst b/doc/source/modules/version.rst new file mode 100644 index 0000000..9b0ced6 --- /dev/null +++ b/doc/source/modules/version.rst @@ -0,0 +1,6 @@ +:orphan: + +.. automodule:: version + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/source/tutorials/AlN/AlN.py b/doc/source/tutorials/AlN/AlN.py index bcc5fe5..7322380 100644 --- a/doc/source/tutorials/AlN/AlN.py +++ b/doc/source/tutorials/AlN/AlN.py @@ -1,181 +1,130 @@ -# 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() +# coding: utf8 + +from ase.build import bulk +import numpy as np +from msspec.calculator import MSSPEC, XRaySource +from msspec.utils import hemispherical_cluster, get_atom_index + +def create_clusters(nplanes=6): + def get_AlN_tags_planes(side, emitter): + AlN = bulk('AlN', crystalstructure='wurtzite', a=3.11, c=4.975) + [atom.set('tag', i) for i, atom in enumerate(AlN)] + if side == 'Al': + AlN.rotate([0,0,1],[0,0,-1]) + Al_planes = range(0, nplanes, 2) + N_planes = range(1, nplanes, 2) + else: + N_planes = range(0, nplanes, 2) + Al_planes = range(1, nplanes, 2) + if emitter == 'Al': + tags = [0, 2] + planes = Al_planes + else: + tags = [1, 3] + planes = N_planes + return AlN, tags, planes + + clusters = [] + for side in ('Al', 'N'): + for emitter in ('Al', 'N'): + AlN, tags, planes = get_AlN_tags_planes(side, emitter) + for emitter_tag in tags: + for emitter_plane in planes: + cluster = hemispherical_cluster(AlN, + emitter_tag=emitter_tag, + emitter_plane=emitter_plane, + planes=emitter_plane+2) + cluster.absorber = get_atom_index(cluster, 0, 0, 0) + cluster.info.update({ + 'emitter_plane': emitter_plane, + 'emitter_tag' : emitter_tag, + 'emitter' : emitter, + 'side' : side, + }) + clusters.append(cluster) + print("Added cluster {}-side, emitter {}(tag {:d}) in " + "plane #{:d}".format(side, emitter, emitter_tag, + emitter_plane)) + return clusters + + +def compute(clusters, theta=np.arange(-20., 80., 1.), phi=0.): + data = None + for ic, cluster in enumerate(clusters): + # Retrieve info from cluster object + side = cluster.info['side'] + emitter = cluster.info['emitter'] + plane = cluster.info['emitter_plane'] + tag = cluster.info['emitter_tag'] + + # Set the level and the kinetic energy + if emitter == 'Al': + level = '2p' + ke = 1407. + elif emitter == 'N': + level = '1s' + ke = 1083. + + calc = MSSPEC(spectroscopy='PED', algorithm='expansion') + + 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.scattering_order = max(1, min(4, plane)) + calc.calculation_parameters.path_filtering = 'forward_scattering' + calc.calculation_parameters.off_cone_events = 1 + [a.set('forward_angle', 30.) 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 = "\'{}\' side - {}({}) tag #{:d}, plane #{:d}".format( + side, emitter, level, tag, plane) + + return data + + +def analysis(data): + tmp_data = {} + for dset in data: + info = dset.get_cluster().info + side = info['side'] + emitter = info['emitter'] + try: + key = '{}_{}'.format(side, emitter) + tmp_data[key] += dset.cross_section + except KeyError: + tmp_data[key] = dset.cross_section.copy() + + tmp_data['theta'] = dset.theta.copy() + tmp_data['Al_side'] = tmp_data['Al_Al'] / tmp_data['Al_N'] + tmp_data['N_side'] = tmp_data['N_Al'] / tmp_data['N_N'] + + # now add all columns + substrate_dset = data.add_dset('Total substrate signal') + substrate_dset.add_columns(**tmp_data) + + view = substrate_dset.add_view('Ratios', + title=r'Al(2p)/N(1s) ratios on both polar ' + r'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) + + return data + + +clusters = create_clusters() +data = compute(clusters) +data = analysis(data) +data.view() diff --git a/doc/source/tutorials/RhO/RhO.py b/doc/source/tutorials/RhO/RhO.py index e719c43..90c8e46 100644 --- a/doc/source/tutorials/RhO/RhO.py +++ b/doc/source/tutorials/RhO/RhO.py @@ -1,54 +1,25 @@ -# -*- encoding: utf-8 -*- -# 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 ase.build import fcc111, add_adsorbate -from ase.visualize import view -from msspec.iodata import cols2matrix - -from matplotlib import pyplot as plt -import numpy as np -import sys - -data = None -all_z = np.arange(1.10, 1.50, 0.02) -all_z=(1.1,) -for zi, z0 in enumerate(all_z): - # construct the cluster - cluster = fcc111('Rh', size = (2,2,1)) - add_adsorbate(cluster, 'O', z0, position = 'fcc') - cluster.pop(3) - #cluster.rotate('z',np.pi/3.) - #view(cluster) - - 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() +# coding: utf8 + +from msspec.calculator import MSSPEC +from ase.build import fcc111, add_adsorbate +import numpy as np + +data = None +all_z = np.arange(1.10, 1.65, 0.05) +for zi, z0 in enumerate(all_z): + # construct the cluster + cluster = fcc111('Rh', size = (2,2,1)) + cluster.pop(3) + add_adsorbate(cluster, 'O', z0, position = 'fcc') + cluster.absorber = len(cluster) - 1 + + # Define a calculator + calc = MSSPEC(spectroscopy='PED', algorithm='inversion') + calc.set_atoms(cluster) + + # Compute + data = calc.get_theta_phi_scan(level='1s', kinetic_energy=723, data=data) + dset = data[-1] + dset.title = "{:d}) z = {:.2f} angstroms".format(zi, z0) + +data.view() diff --git a/doc/source/tutorials/RhO/sf.py b/doc/source/tutorials/RhO/sf.py index 02e2fca..e42c786 100644 --- a/doc/source/tutorials/RhO/sf.py +++ b/doc/source/tutorials/RhO/sf.py @@ -1,31 +1,17 @@ -# -*- encoding: utf-8 -*- -# 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 * +# coding: utf8 +from msspec.calculator import MSSPEC 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 -a0 = 6.0 -symbols = ('Rh', 'O') -ke = 723. -level = '1s' +# Create the calculator +calc = MSSPEC(spectroscopy = 'PED') +calc.set_atoms(cluster) +cluster.absorber = 0 -data = None -for symbol in symbols: - cluster = Atoms(symbol*2, positions = [(0,0,0), (0,0,a0)]) - [a.set('mt_radius', 1.5) for a in cluster] +# compute +data = calc.get_scattering_factors(level='1s', kinetic_energy=723) - # 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() diff --git a/doc/source/tutorials/RhO/tuto_ped_RhO.rst b/doc/source/tutorials/RhO/tuto_ped_RhO.rst index 5247145..0b55f03 100644 --- a/doc/source/tutorials/RhO/tuto_ped_RhO.rst +++ b/doc/source/tutorials/RhO/tuto_ped_RhO.rst @@ -82,13 +82,13 @@ Here is the script for the computation. (:download:`download `) .. literalinclude:: RhO.py :linenos: -.. note:: - After runing this script, you will get 20 images in your folder. You can merge them in one animated gif image - like this: - - .. code-block:: bash - - convert -delay 50 -loop 0 image*.png animation.gif +.. .. note:: +.. After runing this script, you will get 20 images in your folder. You can merge them in one animated gif image +.. like this: +.. +.. .. code-block:: bash +.. +.. convert -delay 50 -loop 0 image*.png animation.gif .. seealso:: diff --git a/doc/source/tutorials/copper/Cu_simple.py b/doc/source/tutorials/copper/Cu.py similarity index 51% rename from doc/source/tutorials/copper/Cu_simple.py rename to doc/source/tutorials/copper/Cu.py index 9801b03..7bb5e48 100644 --- a/doc/source/tutorials/copper/Cu_simple.py +++ b/doc/source/tutorials/copper/Cu.py @@ -1,7 +1,7 @@ -# coding: utf-8 +# coding: utf8 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.visualize import view @@ -9,20 +9,15 @@ from ase.visualize import view a0 = 3.6 # The lattice parameter in angstroms # Create the copper cubic cell -cluster = bulk('Cu', a=a0, cubic=True) -# Repeat the cell many times along x, y and z -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) +copper = bulk('Cu', a=a0, cubic=True) +cluster = hemispherical_cluster(copper, planes=3, emitter_plane=2) # 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 calc = MSSPEC(spectroscopy='PED') + # Set the cluster to use for the calculation calc.set_atoms(cluster) @@ -31,4 +26,6 @@ data = calc.get_theta_scan(level='2p3/2') # Show the results data.view() -#calc.shutdown() + +# Clean temporary files +calc.shutdown() diff --git a/doc/source/tutorials/copper/tuto_ped_copper.rst b/doc/source/tutorials/copper/tuto_ped_copper.rst index ed7037c..0c564f3 100644 --- a/doc/source/tutorials/copper/tuto_ped_copper.rst +++ b/doc/source/tutorials/copper/tuto_ped_copper.rst @@ -16,7 +16,7 @@ Begin by typing: :linenos: :lineno-start: 1 - # coding: utf-8 + # coding: utf8 from msspec.calculator import MSSPEC 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.visualize import view - a0 = 3.6 - cluster = bulk('Cu', a = a0, cubic = True) + a0 = 3.6 # The lattice parameter in angstroms + + # Create the copper cubic cell + copper = bulk('Cu', a=a0, cubic=True) view(cluster) 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. -The creation of the cluster is done line 10. The :py:func:`bulk` needs one -argument which are the chemical species ('Cu' in our example). We also pass 2 -keyword (optional) arguments here: +The creation of the cluster starts on line 12. We create first the copper cubic +cell with the The :py:func:`bulk` function. It needs one argument which are the +chemical species ('Cu' in our example). We also pass 2 keyword (optional) arguments +here: * The lattice parameter *a* in units of angströms. * The *cubic* keyword, to obtain a cubic cell rather than the fully reduced one which is not cubic 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 - 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 like this one: @@ -80,34 +90,38 @@ like this one: Obviously, multiple scattering calculations need more atoms to be accurate. Due to the forward focusing effect in photo-electron diffraction, the best suited geometry for the cluster is hemispherical. Obtaining such a cluster is a -straightforward process: - - 1. Repeat the previous cell in all 3 directions - 2. center the structure - 3. Keep only atoms within a given radius from center - 4. Keep only atoms within one halh of the sphere. +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 +cell here). Modify your script like this and run it again. -.. literalinclude:: Cu_simple.py +.. literalinclude:: Cu.py :linenos: :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 -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 :align: center :width: 60% 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 -will absorbe the light. This atom is called the *absorber*. +Once your cluster is built the next step is to define which atom in the cluster +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 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 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 @@ -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 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 -4 arguments: the cluster to look the index for, and the x, y and z coordinates. +program (select an atom with the left mouse button and look at its index in the +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 -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: - :lineno-start: 22 - :lines: 22-23 + :lineno-start: 15 + :lines: 15-16 That's all for the cluster part. We now need to create a calculator for that object. This is a 2 lines procedure: -.. literalinclude:: Cu_simple.py +.. literalinclude:: Cu.py :linenos: - :lineno-start: 24 - :lines: 24-28 + :lineno-start: 18 + :lines: 18-22 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*. @@ -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, store the results in the data object and display it in a graphical window. -.. literalinclude:: Cu_simple.py +.. literalinclude:: Cu.py :linenos: - :lineno-start: 29 - :lines: 29-33 + :lineno-start: 24 + :lines: 24-28 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. .. code-block:: python - :linenos: #For a polar scan between 0° and 60° with 100 points 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: .. code-block:: python - :linenos: # To set the kinetic energy... 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 -` +` -.. literalinclude:: Cu_simple.py +.. literalinclude:: Cu.py :linenos: .. seealso:: diff --git a/doc/source/tutorials/index.rst b/doc/source/tutorials/index.rst new file mode 100644 index 0000000..7efd61e --- /dev/null +++ b/doc/source/tutorials/index.rst @@ -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 diff --git a/doc/source/tutorials/nickel_chain/Ni_chain.py b/doc/source/tutorials/nickel_chain/Ni_chain.py index 321ebd2..8b5bfc5 100644 --- a/doc/source/tutorials/nickel_chain/Ni_chain.py +++ b/doc/source/tutorials/nickel_chain/Ni_chain.py @@ -1,77 +1,78 @@ -# coding: utf-8 - -# import all we need and start by msspec -from msspec.calculator import MSSPEC - -# we will build a simple atomic chain -from ase import Atom, Atoms - -# we need some numpy functions -import numpy as np - - -symbol = 'Ni' # The kind of atom for the chain -orders = (1, 5) # We will run the calculation for single scattering - # and multiple scattering (5th diffusion order) -chain_lengths = (2,3,5) # We will run the calculation for differnt lengths - # of the atomic chain -a = 3.5 # The distance bewteen 2 atoms - -# Define an empty variable to store all the results -all_data = None - -# 2 for nested loops over the chain length and the order of diffusion -for chain_length in chain_lengths: - for order in orders: - # We build the atomic chain by - # 1- stacking each atom one by one along the z axis - chain = Atoms([Atom(symbol, position = (0., 0., i*a)) for i in - range(chain_length)]) - # 2- rotating the chain by 45 degrees with respect to the y axis - chain.rotate('y', np.radians(45.)) - # 3- setting a custom Muffin-tin radius of 1.5 angstroms for all - # atoms (needed if you want to enlarge the distance between - # the atoms while keeping the radius constant) - [atom.set('mt_radius', 1.5) for atom in chain] - # 4- defining the absorber to be the first atom in the chain at - # x = y = z = 0 - chain.absorber = 0 - - # We define a new PED calculator - calc = MSSPEC(spectroscopy = 'PED') - calc.set_atoms(chain) - # Here is how to tweak the scattering order - calc.calculation_parameters.scattering_order = order - # This line below is where we actually run the calculation - 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 - t = "order {:d}, n = {:d}".format(order, chain_length) # A useful title - dset = all_data[-1] # get the last dataset - dset.title = t # change its title - # get its last view (there is only one defined for each dataset) - 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 -# maximum of the data -min_cs = max_cs = 0 -for dset in all_data: - 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: - v = dset.views()[-1] - v.set_plot_options(ylim=[min_cs, max_cs]) - -# Pop up the graphical window -all_data.view() -# You can end your script with the line below to remove the temporary -# folder needed for the calculation -calc.shutdown() +# coding: utf-8 + +# import all we need and start by msspec +from msspec.calculator import MSSPEC + +# we will build a simple atomic chain +from ase import Atom, Atoms + +# we need some numpy functions +import numpy as np + + +symbol = 'Ni' # The kind of atom for the chain +orders = (1, 5) # We will run the calculation for single scattering + # and multiple scattering (5th diffusion order) +chain_lengths = (2,3,5) # We will run the calculation for differnt lengths + # of the atomic chain +a = 3.499 * np.sqrt(2)/2 # The distance bewteen 2 atoms + +# Define an empty variable to store all the results +all_data = None + +# 2 for nested loops over the chain length and the order of diffusion +for chain_length in chain_lengths: + for order in orders: + # We build the atomic chain by + # 1- stacking each atom one by one along the z axis + chain = Atoms([Atom(symbol, position = (0., 0., i*a)) for i in + range(chain_length)]) + # 2- rotating the chain by 45 degrees with respect to the y axis + #chain.rotate('y', np.radians(45.)) + chain.rotate(45., 'y') + # 3- setting a custom Muffin-tin radius of 1.5 angstroms for all + # atoms (needed if you want to enlarge the distance between + # the atoms while keeping the radius constant) + #[atom.set('mt_radius', 1.5) for atom in chain] + # 4- defining the absorber to be the first atom in the chain at + # x = y = z = 0 + chain.absorber = 0 + + # We define a new PED calculator + calc = MSSPEC(spectroscopy = 'PED') + calc.set_atoms(chain) + # Here is how to tweak the scattering order + calc.calculation_parameters.scattering_order = order + # This line below is where we actually run the calculation + 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 + t = "order {:d}, n = {:d}".format(order, chain_length) # A useful title + dset = all_data[-1] # get the last dataset + dset.title = t # change its title + # get its last view (there is only one defined for each dataset) + 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 +# maximum of the data +min_cs = max_cs = 0 +for dset in all_data: + 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: + v = dset.views()[-1] + v.set_plot_options(ylim=[min_cs, max_cs]) + +# Pop up the graphical window +all_data.view() +# You can end your script with the line below to remove the temporary +# folder needed for the calculation +calc.shutdown() diff --git a/doc/source/tutorials/nickel_chain/tuto_ped_nickel_chain.rst b/doc/source/tutorials/nickel_chain/tuto_ped_nickel_chain.rst index 27f04b6..affaf1a 100644 --- a/doc/source/tutorials/nickel_chain/tuto_ped_nickel_chain.rst +++ b/doc/source/tutorials/nickel_chain/tuto_ped_nickel_chain.rst @@ -30,7 +30,7 @@ trajectory for this electron which in turn lowers the signal. :align: center :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. diff --git a/doc/source/tutorials/temperature/Cu_temperature.py b/doc/source/tutorials/temperature/Cu_temperature.py index 4120366..c159c0d 100644 --- a/doc/source/tutorials/temperature/Cu_temperature.py +++ b/doc/source/tutorials/temperature/Cu_temperature.py @@ -1,153 +1,129 @@ -# coding: utf-8 - -from msspec.calculator import MSSPEC, XRaySource -from msspec.utils import * -from msspec.iodata import Data -from ase.build import bulk - -import numpy as np -import sys - -a0 = 3.6 # The lattice parameter in angstroms -phi = np.arange(0, 100) # An array of phi angles -all_T = np.arange(300, 1000, 100) # All temperatures used for the calculation -all_theta = np.array([45, 83]) # 2 polar angles, 83° is grazing detection -eps = 0.01 # a useful small value - -DATA_FNAME = 'all_data.hdf5' # Where to store all the data -ANALYSIS_FNAME = 'analysis.hdf5' - - -def compute(filename): - """Will compute a phi scan for an emitter in the first, second and third plane - of a copper substrate for various polar angles and temperatures. - In a second step (outside this function), intensities from all these emitters - are added to get the signal of a substrate.""" - calc = MSSPEC(spectroscopy='PED') - calc.source_parameters.energy = XRaySource.AL_KALPHA - calc.muffintin_parameters.interstitial_potential = 14.1 - - calc.calculation_parameters.vibrational_damping = 'averaged_tl' - calc.calculation_parameters.use_debye_model = True - calc.calculation_parameters.debye_temperature = 343 - calc.calculation_parameters.vibration_scaling = 1.2 - - calc.detector_parameters.average_sampling = 'high' - calc.detector_parameters.angular_acceptance = 5.7 - - filters = ['forward_scattering', 'backward_scattering'] - calc.calculation_parameters.path_filtering = filters - - calc.calculation_parameters.RA_cutoff = 3 - - # You can also choose a real potential and manually define the - # electron mean free path - # - # calc.tmatrix_parameters.exchange_correlation = 'hedin_lundqvist_real' - # calc.calculation_parameters.mean_free_path = 10. - - data = None # init an empty data object - for temperature in all_T: - for plane in range(3): - # We create a cylindrical cluster here with one plane below the emitter - # and 0, 1 or to planes above the emitter - cluster = bulk('Cu', a=a0, cubic=True) - cluster = cluster.repeat((20, 20, 8)) - center_cluster(cluster) - cluster = cut_cylinder(cluster, radius=1.5 * a0 + eps) - cluster = cut_plane(cluster, z=-( a0/2 + eps)) - cluster = cut_plane(cluster, z=(plane) * a0 / 2 + eps) - cluster.absorber = get_atom_index(cluster, 0, 0, 0) - - calc.calculation_parameters.temperature = temperature - # the scattering order depends on the number of planes above - # the emitter to speed up calculations - calc.calculation_parameters.scattering_order = 1 + plane - - calc.set_atoms(cluster) - # Here starts the calculation - data = calc.get_phi_scan(level='2p3/2', theta=all_theta, phi=phi, - kinetic_energy=560, data=data) - - data.save(filename) - -def process_data(datafile, outputfile): - """Will create another Data file with some phi-scans corresponding to 3 - planes at different temperatures and the anisotropy curve for 45° and - grazing detection. - """ - def get_signal(datafile, T=300, theta=45): - data = Data.load(datafile) - total = None - for dset in data: - p = {_['group'] + '.' + _['name']: _['value'] for _ in dset.parameters()} - temperature = np.float(p['CalculationParameters.temperature']) - if temperature != T: continue - i = np.where(dset.plane == -1) - j = np.where(dset.theta[i] == theta) - signal = dset.cross_section[i][j] - try: - total += signal - except: - total = signal - phi = dset.phi[i][j] - return phi, total - - analysis = Data('Temperature tutorial') - scans_dset = analysis.add_dset('Phi scans') - scans_view = scans_dset.add_view('Figure 1', - title=r'Cu(2p) Azimuthal scan for $\theta = 83\degree$', - xlabel=r'$\Phi (\degree$)', - ylabel='Signal (a.u.)') - anisotropy_dset = analysis.add_dset('Anisotropies') - anisotropy_view = anisotropy_dset.add_view('Figure 2', - title='Relative anisotropies for Cu(2p)', - marker='o', - xlabel='T (K)', - ylabel=r'$\frac{\Delta I / I_{max}(T)}{\Delta I_{300} / I_{max}(300)} (\%)$') - for theta in all_theta: - for T in all_T: - PHI, SIGNAL = get_signal(datafile, T=T, theta=theta) - for phi, signal in zip(PHI, SIGNAL): - scans_dset.add_row(phi=phi, signal=signal, theta=theta, temperature=T) - - middleT = all_T[np.abs(all_T - np.mean(all_T)).argmin()] - if theta == 83 and T in [np.min(all_T), middleT, np.max(all_T)]: - scans_view.select('phi', 'signal', - where='temperature == {:f} and theta == {:f}'.format(T, theta), - legend='{:.0f} K'.format(T)) - - PHI, SIGNAL0 = get_signal(datafile, T=np.min(all_T), theta=theta) - Imax = np.max(SIGNAL) - Imax0 = np.max(SIGNAL0) - dI = Imax - np.min(SIGNAL) - dI0 = Imax0 - np.min(SIGNAL0) - ani = (dI / Imax) / (dI0 / Imax0) - anisotropy_dset.add_row(temperature=T, anisotropy=ani, theta=theta) - - anisotropy_view.select('temperature', 'anisotropy', - where='theta == {:f}'.format(theta), - 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() - - - - - +# coding: utf8 + +from ase.build import bulk +import numpy as np +from msspec.calculator import MSSPEC, XRaySource +from msspec.utils import hemispherical_cluster, get_atom_index + + +def create_clusters(nplanes=3): + copper = bulk('Cu', a=3.6, cubic=True) + clusters = [] + for emitter_plane in range(nplanes): + cluster = hemispherical_cluster(copper, + emitter_plane=emitter_plane, + planes=emitter_plane+2, + shape='cylindrical') + cluster.absorber = get_atom_index(cluster, 0, 0, 0) + cluster.info.update({ + 'emitter_plane': emitter_plane, + }) + clusters.append(cluster) + return clusters + +def compute(clusters, all_theta=[45., 83.], + all_T=np.arange(300., 1000., 400.)): + data = None + for ic, cluster in enumerate(clusters): + # Retrieve info from cluster object + plane = cluster.info['emitter_plane'] + + calc = MSSPEC(spectroscopy='PED', algorithm='expansion') + calc.source_parameters.energy = XRaySource.AL_KALPHA + calc.muffintin_parameters.interstitial_potential = 14.1 + + calc.calculation_parameters.vibrational_damping = 'averaged_tl' + calc.calculation_parameters.use_debye_model = True + calc.calculation_parameters.debye_temperature = 343 + calc.calculation_parameters.vibration_scaling = 1.2 + + calc.detector_parameters.average_sampling = 'high' + calc.detector_parameters.angular_acceptance = 5.7 + + for atom in cluster: + atom.set('forward_angle', 30) + atom.set('backward_angle', 30) + filters = ['forward_scattering', 'backward_scattering'] + calc.calculation_parameters.path_filtering = filters + + calc.calculation_parameters.RA_cutoff = 2 + + for T in all_T: + for theta in all_theta: + calc.calculation_parameters.temperature = T + calc.calculation_parameters.scattering_order = min(1 + plane, 3) + calc.set_atoms(cluster) + data = calc.get_phi_scan(level='2p3/2', theta=theta, + phi=np.linspace(0, 100), + kinetic_energy=560, data=data) + dset = data[-1] + dset.title = "plane #{:d}, T={:f}K, theta={:f}°".format(plane, + T, + theta) + + dset.add_parameter(group='misc', name='plane', value=plane, unit='') + dset.add_parameter(group='misc', name='T', value=T, unit='') + dset.add_parameter(group='misc', name='theta', value=theta, unit='') + return data + + +def analysis(data): + all_plane = [] + all_T = [] + all_theta = [] + for dset in data: + plane = dset.get_parameter('misc', 'plane')['value'] + T = dset.get_parameter('misc', 'T')['value'] + theta = dset.get_parameter('misc', 'theta')['value'] + cs = dset.cross_section.copy() + phi = dset.phi.copy() + + if plane not in all_plane: all_plane.append(plane) + if T not in all_T: all_T.append(T) + if theta not in all_theta: all_theta.append(theta) + + def get_anisotropy(theta, T): + cs = None + for dset in data: + try: + _T = dset.get_parameter('misc', 'T')['value'] + _theta = dset.get_parameter('misc', 'theta')['value'] + _cs = dset.cross_section.copy() + phi = dset.phi.copy() + except: + continue + if _theta == theta and _T == T: + try: + cs += _cs + except: + cs = _cs + Imax = np.max(cs) + Imin = np.min(cs) + return (Imax - Imin)/Imax + + # create a substrate dataset for each T and theta + anisotropy_dset = data.add_dset("all") + anisotropy_view = anisotropy_dset.add_view('Anisotropies', + title='Relative anisotropies for Cu(2p)', + marker='o', + xlabel='T (K)', + ylabel=r'$\frac{\Delta I / I_{max}(T)}{\Delta I_{300}' + r'/ I_{max}(300)} (\%)$') + + for theta in all_theta: + for T in all_T: + A = get_anisotropy(theta, T) + A0 = get_anisotropy(theta, np.min(all_T)) + + anisotropy_dset.add_row(temperature=T, theta=theta, anisotropy=A/A0) + + anisotropy_view.select('temperature', 'anisotropy', + where='theta == {:f}'.format(theta), + legend=r'$\theta = {:.0f} \degree$'.format(theta)) + return data + + +clusters = create_clusters() +data = compute(clusters) +data = analysis(data) +data.view()