Update Activity06

msspecbook/Activity06/Activity06.ipynb

@@ -17,11 +17,11 @@
     "\n",
     "For each azimutal scan, they looked at the anisotropy of the signal, that is:\n",
     "\n",
-    "$\\frac{\\Delta I}{I_{max}}$\n",
+    "$\\frac{I_{max} - I_{min}}{I_{max}}=\\frac{\\Delta I}{I_{max}}$\n",
     "\n",
     "This value is representative of how clear are the *modulations* of the signal. As it was shown by their experiments, this anisotropy decreases when the temperature is increased due to the increased disorder in the structure coming from thermal agitation. They also showed that this variation in anisotropy is more pronounced for grazing incidence angles. This is related to the fact that surface atoms are expected to vibrate more than bulk ones. They also proposed single scattering calculations that reproduced well these results.\n",
     "\n",
-    "We propose here to reproduce this kind of calculation to introduce the parameters that control the vibrational damping.\n",
+    "We will reproduce this kind of calculations to introduce the parameters that control the vibrational damping.\n",
     "\n",
     ":::{seealso}\n",
     "based on this paper from R. Trehan and C.S. Fadley\n",
@@ -30,19 +30,38 @@
     "\n",
     "### The script\n",
     "\n",
-    "Let's start by downloading this [python script ](./Cu_temperature.py). \n",
+    "Since we want to distinguish between bulk (low polar angles, $\\theta = 45°$ in the paper) and surface effects (large polar angles, $\\theta = 83°$ in the paper), we need to compute scans for different emitters and different depths in the cluster.\n",
-    "\n",
-    "Since we want to distinguish between bulk (low polar angles) and surface effects (large polar angles), we need to compute scans for different emitters and different depths in the cluster.\n",
     "\n",
     "The script contains 3 functions:\n",
-    "1. The `create_clusters` function will build 3 cluster: one with the emitter on the surface, one in the subsurface and one n the 3{sup}`rd` plane.\n",
+    "1. The `create_clusters` function will build clusters with increasing number of planes, with the emitter being in the deepest plane for each cluster.\n",
-    "2. The function `compute` will compute the azimuthal scan"
+    "2. The function `compute` will compute the azimuthal scan for a given cluster.\n",
+    "3. The `analysis` function will sum the intensity from each plane for a given temperature. This will be the 'substrate total signal' (more on this in the following activity). The function will also compute the anisotropy"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "69b8c17d-ab66-4092-a492-005f05d80495",
+   "metadata": {},
+   "source": [
+    "::::{tab-set}\n",
+    "\n",
+    ":::{tab-item} <i class=\"fa-solid fa-circle-question\"></i> Quiz\n",
+    "Complete the hilighted lines in the following script\n",
+    "\n",
+    "```{literalinclude} Cu_temperature.py\n",
+    ":lineno-match:\n",
+    ":emphasize-lines:1,2\n",
+    "```\n",
+    "\n",
+    ":::\n",
+    "\n",
+    "::::"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "b8de98a0-12f6-49f4-99d9-db561e6738bf",
+   "id": "cfc8af7b-fec1-449c-9cf7-4d9e2ebe026a",
    "metadata": {},
    "outputs": [],
    "source": []

msspecbook/Activity06/Cu_temperature.py

@@ -0,0 +1,130 @@
+# 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.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()