diff --git a/msspecbook/Activity05/Activity05.ipynb b/msspecbook/Activity05/Activity05.ipynb index ca3f110..a7fe237 100644 --- a/msspecbook/Activity05/Activity05.ipynb +++ b/msspecbook/Activity05/Activity05.ipynb @@ -878,7 +878,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.3" + "version": "3.11.13" } }, "nbformat": 4, diff --git a/msspecbook/Activity06/Activity06.ipynb b/msspecbook/Activity06/Activity06.ipynb index f114ebd..a95699d 100644 --- a/msspecbook/Activity06/Activity06.ipynb +++ b/msspecbook/Activity06/Activity06.ipynb @@ -133,7 +133,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.3" + "version": "3.11.13" } }, "nbformat": 4, diff --git a/msspecbook/Activity07/Activity07.ipynb b/msspecbook/Activity07/Activity07.ipynb index 9f612ce..1e99041 100644 --- a/msspecbook/Activity07/Activity07.ipynb +++ b/msspecbook/Activity07/Activity07.ipynb @@ -103,7 +103,8 @@ ":::{figure-md} filters-fig\n", "\"path\n", "\n", - "Some examples of scattering paths with `forward_scattering`, `backward_scattering` and `distance` filters selected. The accepted forward angle is 45°, the accepted backscattering angle is 20° and the threshold distance is $6a_0$ where $a_0$ is the lattice parameter. Note that the yellow path is rejected but if the `off_cone_events` option is set to a value > 1, then it could have been accepted." + "Some examples of scattering paths with `forward_scattering`, `backward_scattering` and `distance` filters selected. The accepted forward angle is 45°, the accepted backscattering angle is 20° and the threshold distance is $6a_0$ where $a_0$ is the lattice parameter. Note that the yellow path is rejected but if the `off_cone_events` option is set to a value > 1, then it could have been accepted.\n", + ":::" ] }, { @@ -113,7 +114,9 @@ "source": [ "## Application to a deep plane in a Si(001) sample\n", "\n", - "The following script will compute contribution of all the planes of a Si(001) substrate to get the total photoelectron intensity of a Si(2s) polar scan. \n", + "The following script will compute the contribution of a Si(2p) atom in the 4{sup}`th` plane of a Si(001) cluster at scattering order 3.\n", + "\n", + "Taking into account all scattering paths took 15 minutes to compute.\n", "\n", "(msd-paper)=\n", ":::{seealso}\n", @@ -126,14 +129,46 @@ "\n", ":::{tab-item} Quiz\n", "\n", - "The script is almost completed, try to define path filtering options and compare results with and without filtering for emitter in plane n° 3 at scattering order 2.\n", + "The following script is almost completed, try to define path filtering options (no backscattering, accept all paths with forward angles < 40° and reject paths longer than the diameter of the cluster).\n", "\n", - "Compute the contribution of plane n° 7\n", + "```{literalinclude} Si001.py\n", + ":lineno-match:\n", + ":emphasize-lines: 37-41\n", + "```\n", + "\n", + "1. How long was your calculation ?\n", + "2. How does it compare to the calculation with **all** scattering paths up to order 3 ?\n", + "3. What is the proportion of scattering paths of order 3 that were actually taken into account ?\n", "\n", ":::\n", "\n", "::::" ] + }, + { + "cell_type": "markdown", + "id": "19fbd486-b0c1-450c-a00d-79984945aefd", + "metadata": {}, + "source": [ + "```{toggle}\n", + "The calculation took few seconds and the result is very close to the calculation with all scattering paths.\n", + "\n", + "Only 0.01% of 3{sup}`rd` order paths were actually taken into account\n", + "\n", + ":::{figure-md} si-fig\n", + "\"Si\n", + "\n", + "Si(2p) polar scan (contribution of an emitter in the 4{sup}`th` plane with all 7 114 945 scattering paths taken into account (orange curve), and for only 1525 filtered paths (blue curve).\n", + "\n", + ":::\n", + "\n", + ":::{literalinclude} Si001_completed.py\n", + ":lineno-match:\n", + ":emphasize-lines: 37-41\n", + ":::\n", + "\n", + "``` " + ] } ], "metadata": { @@ -152,7 +187,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.3" + "version": "3.11.13" } }, "nbformat": 4, diff --git a/msspecbook/Activity07/Si001.py b/msspecbook/Activity07/Si001.py index dce39d8..af37f6f 100644 --- a/msspecbook/Activity07/Si001.py +++ b/msspecbook/Activity07/Si001.py @@ -1,165 +1,56 @@ # coding: utf8 -import logging import numpy as np -import os -import sys - -from ase.build import bulk, diamond111 +from ase.build import bulk from msspec.calculator import MSSPEC, XRaySource from msspec.iodata import Data -from msspec.utils import hemispherical_cluster, get_atom_index, cut_sphere, cut_plane - -logging.basicConfig(level=logging.INFO) - -def create_clusters(nplanes=6, direction='100', a=5.43, c=5.43, radius=30): - clusters = [] - Si = bulk('Si', a=a, cubic=True) - - # Get scaled positions and cell - p = Si.get_scaled_positions() - cell = Si.get_cell() - - # Resize - covera = c / a - cell[2,:] *= covera - p[:,2] *= covera - Si.set_scaled_positions(p) - Si.set_cell(cell) - - if direction in ('001', '010', '100'): - pass - elif direction in ('111'): - Si.rotate([1,1,1], [1,1,0], rotate_cell=True) - Si.rotate([1,1,0], [1,0,0], rotate_cell=True) - Si.rotate([1,0,0], [0,0,1], rotate_cell=True) - - for iplane in range(nplanes): - #for iplane in (nplanes-1,): - logging.info(f'Building cluster #{iplane:d}/{nplanes-1:d}') - cluster = hemispherical_cluster(Si, #planes=max(iplane+1, 4), - diameter=70, - emitter_plane=iplane, - shape = 'cylindrical', - ) - - cluster = cut_sphere(cluster, radius = radius) - cluster = cut_plane(cluster, z = -a/4) - - for atom in cluster: - atom.set('mean_square_vibration', 0.006) - atom.set('mt_radius', 1.1) - cluster.absorber = get_atom_index(cluster, 0, 0, 0) - cluster.info.update(emitter_plane=iplane) - clusters.append(cluster) - return clusters - -def compute_polar_scan(cluster, data=None, phi=0, folder='calc', RA_cutoff=2, - max_ndif=6, level='2p', kinetic_energy=1382.28): - emitter_plane = cluster.info['emitter_plane'] - - calc = MSSPEC(spectroscopy='PED', algorithm='expansion', - folder=folder) - - calc.muffintin_parameters.interstitial_potential = 0. - - #calc.tmatrix_parameters.exchange_correlation = 'x_alpha_real' - calc.tmatrix_parameters.exchange_correlation = 'hedin_lundqvist_complex' - - calc.source_parameters.energy = XRaySource.AL_KALPHA - calc.source_parameters.theta = -54.7 - calc.source_parameters.phi = 90 - calc.spectroscopy_parameters.final_state = 1 - calc.detector_parameters.angular_acceptance = 1.0 - calc.detector_parameters.average_sampling = 'high' - - calc.calculation_parameters.scattering_order = min(max_ndif, - max(1, emitter_plane)) - #calc.tmatrix_parameters.lmax_mode ='imposed' - #calc.tmatrix_parameters.lmaxt = 10 - calc.tmatrix_parameters.tl_threshold = 1e-4 - #calc.calculation_parameters.scattering_order = 6 - #calc.calculation_parameters.mean_free_path = 10. - calc.calculation_parameters.vibrational_damping = 'averaged_tl' - calc.calculation_parameters.RA_cutoff = RA_cutoff - my_filters = ('forward_scattering', 'distance_cutoff') - calc.calculation_parameters.path_filtering = my_filters - #calc.calculation_parameters.path_filtering = 'forward_scattering' - calc.calculation_parameters.distance = 30 - calc.calculation_parameters.off_cone_events = 0 - calc.calculation_parameters.RA_cutoff_damping = 1 - [a.set('forward_angle', 40) for a in cluster] - - calc.phagen_parameters.noproto = '.false.' - - calc.set_atoms(cluster) - - data = calc.get_theta_scan(level=level, - theta=np.arange(-30., 80., 0.5), - kinetic_energy = kinetic_energy, - phi=phi, data=data) - - return data - -def sum_planes(data): - cs = None - ds = None - for dset in data: - theta = dset.theta - phi = dset.phi - try: - cs += dset.cross_section - ds += dset.direct_signal - except: - cs = dset.cross_section.copy() - ds = dset.direct_signal.copy() - - dset = data.add_dset('Substrate signal') - dset.add_columns(theta=theta, phi=phi, cross_section=cs, - direct_signal=ds, chi=(cs-ds)/ds) - - phi_values = np.unique(phi) - - view = dset.add_view('Substrate signal', - title=r'Si(2p)', - xlabel=r'$\Theta (\degree$)', - ylabel='Signal (a.u.)') - - for phi_value in phi_values: - condition = f'phi == {phi_value:.1f}' - legend = '$\Phi = $' + f'{phi_value:.1f}' + '$^\circ$' - view.select('theta', 'cross_section', where=condition, legend=legend) - view.set_plot_options(autoscale=True) +from msspec.utils import hemispherical_cluster, get_atom_index -# Create the list of clusters -RA = 2 #int(sys.argv[1]) -NDIF = 5 #int(sys.argv[2]) -RAD = 20. #float(sys.argv[3]) -PHI = float(sys.argv[1]) -KE = float(sys.argv[2]) -#basename = f'RA{RA:d}_NDIF{NDIF:d}_RAD{RAD:.1f}_PHI{PHI:.1f}' -basename = f'PHI{PHI:.1f}_KE{KE:.2f}' -data = Data(basename) -clusters = create_clusters(nplanes=15, radius=RAD) +# Create the cluster +a = 5.43 +Si = bulk('Si', a=a, cubic=True) +cluster = hemispherical_cluster(Si, + diameter=30, planes=4, + emitter_plane=3, + shape = 'cylindrical', + ) +for atom in cluster: + atom.set('mean_square_vibration', 0.006) + atom.set('mt_radius', 1.1) +cluster.emitter = get_atom_index(cluster, 0, 0, 0) -if 'view' in sys.argv: - for cluster in clusters: - cluster.edit() - exit(0) +# Create a calculator and set parameters +calc = MSSPEC(spectroscopy='PED', algorithm='expansion') -for cluster in clusters: -#for cluster in (clusters[-1],): - data = compute_polar_scan(cluster, data, phi=PHI, - folder='calc_' + basename, - RA_cutoff=RA, - max_ndif=NDIF, - kinetic_energy=KE - ) +calc.source_parameters.energy = XRaySource.AL_KALPHA +calc.source_parameters.theta = -54.7 +calc.source_parameters.phi = 90 +calc.spectroscopy_parameters.final_state = 1 -sum_planes(data) -data.save('simulation_' + basename + '.hdf5', append=False) -data.export('simulation_' + basename) -#data.view() +calc.calculation_parameters.scattering_order = 3 +calc.tmatrix_parameters.tl_threshold = 1e-4 +calc.calculation_parameters.vibrational_damping = 'averaged_tl' +calc.calculation_parameters.RA_cutoff = 2 + +# Define path filtering options such that you only +# accept scattering paths with a forward cone <= 40° +# and whose length are <= cluster diameter +# +# +calc.set_atoms(cluster) + +# Compute and add previous data for comparison +data = calc.get_theta_scan(level='2p', + theta=np.arange(-30., 80., 0.5), + phi=0, + kinetic_energy=1382.28) +no_filters = Data.load('path_filtering.hdf5') +data[0].add_columns(**{'no_filters': no_filters[0].cross_section}) +view = data[0].views[0] +view.select('theta', 'cross_section', index=0, legend="With path filtering") +view.select('theta', 'no_filters', legend="Without path filtering") + +data.view() diff --git a/msspecbook/Activity07/Si001_completed.py b/msspecbook/Activity07/Si001_completed.py new file mode 100644 index 0000000..81caa94 --- /dev/null +++ b/msspecbook/Activity07/Si001_completed.py @@ -0,0 +1,56 @@ +# coding: utf8 + +import numpy as np +from ase.build import bulk + +from msspec.calculator import MSSPEC, XRaySource +from msspec.iodata import Data +from msspec.utils import hemispherical_cluster, get_atom_index + + +# Create the cluster +a = 5.43 +Si = bulk('Si', a=a, cubic=True) +cluster = hemispherical_cluster(Si, + diameter=30, planes=4, + emitter_plane=3, + shape = 'cylindrical', + ) +for atom in cluster: + atom.set('mean_square_vibration', 0.006) + atom.set('mt_radius', 1.1) +cluster.emitter = get_atom_index(cluster, 0, 0, 0) + +# Create a calculator and set parameters +calc = MSSPEC(spectroscopy='PED', algorithm='expansion') + +calc.source_parameters.energy = XRaySource.AL_KALPHA +calc.source_parameters.theta = -54.7 +calc.source_parameters.phi = 90 +calc.spectroscopy_parameters.final_state = 1 + +calc.calculation_parameters.scattering_order = 3 +calc.tmatrix_parameters.tl_threshold = 1e-4 +calc.calculation_parameters.vibrational_damping = 'averaged_tl' +calc.calculation_parameters.RA_cutoff = 2 + +my_filters = ('forward_scattering', 'distance_cutoff') +calc.calculation_parameters.path_filtering = my_filters +calc.calculation_parameters.distance = 30 +calc.calculation_parameters.off_cone_events = 0 +[a.set('forward_angle', 40) for a in cluster] + +calc.set_atoms(cluster) + +# Compute and add previous data for comparison +data = calc.get_theta_scan(level='2p', + theta=np.arange(-30., 80., 0.5), + phi=0, + kinetic_energy=1382.28) +no_filters = Data.load('path_filtering.hdf5') +data[0].add_columns(**{'no_filters': no_filters[0].cross_section}) +view = data[0].views[0] +view.select('theta', 'cross_section', index=0, legend="With path filtering") +view.select('theta', 'no_filters', legend="Without path filtering") + +data.view() diff --git a/msspecbook/Activity07/nbpaths.py b/msspecbook/Activity07/nbpaths.py index e90e9f4..448532f 100644 --- a/msspecbook/Activity07/nbpaths.py +++ b/msspecbook/Activity07/nbpaths.py @@ -11,6 +11,7 @@ natoms = [] allt = [] for emitter_plane in planes: cluster = hemispherical_cluster(Si, diameter=40, emitter_plane=emitter_plane, planes=emitter_plane+1) + cluster.edit() N = len(cluster) npaths = np.sum([(N-1)**i for i in range(emitter_plane + 1)]) t = npaths * 1e-6 @@ -54,4 +55,4 @@ plt.ylim(1e-6, 1e20) plt.xlabel("Number of atoms for emitter in plane 1 to 9") plt.ylabel("Computation time (s)") plt.savefig('my_plot.png', transparent=True) -plt.show() \ No newline at end of file +plt.show() diff --git a/msspecbook/Activity07/path_filtering.hdf5 b/msspecbook/Activity07/path_filtering.hdf5 new file mode 100644 index 0000000..893077e Binary files /dev/null and b/msspecbook/Activity07/path_filtering.hdf5 differ diff --git a/msspecbook/Activity07/results.png b/msspecbook/Activity07/results.png new file mode 100644 index 0000000..2dbf754 Binary files /dev/null and b/msspecbook/Activity07/results.png differ