Completed Activity 7

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,

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,

msspecbook/Activity07/Activity07.ipynb

@@ -103,7 +103,8 @@
     ":::{figure-md} filters-fig\n",
     "<img src=\"filters.jpg\" alt=\"path filtering\" width=\"600px\" align=\"center\">\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} <i class=\"fa-solid fa-circle-question\"></i> 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",
+    "<img src=\"results.png\" alt=\"Si polar scan\" width=\"600px\" align=\"center\">\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,

msspecbook/Activity07/Si001.py

@@ -1,165 +1,56 @@
 # coding: utf8
 
-import logging
 import numpy as np
-import os
+from ase.build import bulk
-import sys
-
-from ase.build import bulk, diamond111
 
 from msspec.calculator import MSSPEC, XRaySource
 from msspec.iodata import Data
-from msspec.utils import hemispherical_cluster, get_atom_index, cut_sphere, cut_plane
+from msspec.utils import hemispherical_cluster, get_atom_index
 
-logging.basicConfig(level=logging.INFO)
 
-def create_clusters(nplanes=6, direction='100', a=5.43, c=5.43, radius=30):
+# Create the cluster
-    clusters = []
+a = 5.43
-    Si = bulk('Si', a=a, cubic=True)
+Si = bulk('Si', a=a, cubic=True)
-
+cluster = hemispherical_cluster(Si,
-    # Get scaled positions and cell
+                                diameter=30, planes=4,
-    p = Si.get_scaled_positions()
+                                emitter_plane=3,
-    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',
                                 )
-
+for atom in cluster:
-        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.emitter = 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,
+# Create a calculator and set parameters
-                       max_ndif=6, level='2p', kinetic_energy=1382.28):
+calc = MSSPEC(spectroscopy='PED', algorithm='expansion')
-    emitter_plane = cluster.info['emitter_plane']
 
-    calc = MSSPEC(spectroscopy='PED', algorithm='expansion',
+calc.source_parameters.energy = XRaySource.AL_KALPHA
-                  folder=folder)
+calc.source_parameters.theta  = -54.7
+calc.source_parameters.phi    = 90
+calc.spectroscopy_parameters.final_state = 1
 
-    calc.muffintin_parameters.interstitial_potential = 0.
+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
  
-    #calc.tmatrix_parameters.exchange_correlation = 'x_alpha_real'
+# Define path filtering options such that you only
-    calc.tmatrix_parameters.exchange_correlation = 'hedin_lundqvist_complex'
+# accept scattering paths with a forward cone <= 40°
+# and whose length are <= cluster diameter
+#
+#
 
-    calc.source_parameters.energy = XRaySource.AL_KALPHA
+calc.set_atoms(cluster)
-    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,
+# Compute and add previous data for comparison
-                                                       max(1, emitter_plane))
+data = calc.get_theta_scan(level='2p',
-    #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=0,
-                               phi=phi, data=data)
+                            kinetic_energy=1382.28)
-
+no_filters = Data.load('path_filtering.hdf5')
-    return data
+data[0].add_columns(**{'no_filters': no_filters[0].cross_section})
-
+view = data[0].views[0]
-def sum_planes(data):
+view.select('theta', 'cross_section', index=0, legend="With path filtering")
-    cs = None
+view.select('theta', 'no_filters', legend="Without path filtering")
-    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)
-
-
-# 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)
-
-if 'view' in sys.argv:
-    for cluster in clusters:
-        cluster.edit()
-    exit(0)
-
-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
-                              )
-
-sum_planes(data)
-data.save('simulation_' + basename + '.hdf5', append=False)
-data.export('simulation_' + basename)
-#data.view()
 
+data.view()

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

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