{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# solving rate equations for the multi-photon multiple ionization dynamics due to intense x-ray radiation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "CFEL 2022, Ludger Inhester" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T08:43:54.468239Z", "start_time": "2022-05-31T08:43:50.188704Z" } }, "outputs": [], "source": [ "import numpy as np\n", "from scipy.integrate import solve_ivp\n", "from matplotlib import pyplot as plt\n", "AU2SEC = 2.418884326500e-17\n", "FS2AU = 1e-15 / AU2SEC\n", "AU2FS = 1 / FS2AU\n", "AU2M = 0.5291772083e-10\n", "M2AU = 1. / AU2M" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T08:43:54.488688Z", "start_time": "2022-05-31T08:43:54.471183Z" } }, "outputs": [], "source": [ "def solveREQ(pIni,\n", " processData,\n", " configurations,\n", " fluence,\n", " pulseEnv,\n", " tStart=-150 * FS2AU,\n", " tEnd=150 * FS2AU,\n", " dt=0.1 * FS2AU,\n", " method='LSODA'):\n", " \"\"\" \n", " solves the rate equation and returns results for the population\n", " Mandatory arguments: \n", " pIni: initial population, array of shape(nconfigurations)\n", " processData: process data, dictionary obtained from readProcessData\n", " configurations: array of configurations, array of shape (nconfigurations, number of occupations)\n", " fluence: fluence value, scalar (number of photons per area in atomic units)\n", " pulseEnv: pulseEnvelop, function f(t), t in atomic units \n", " optional arguments:\n", " tStart: time to start integration (i.e., where p=pIni) in atomic units, default -100*FSAU\n", " tEnd: time to end integration in atomic units, default 100*FSAU\n", " dt: timestep in atomic units for which to provide results, default 0.1*FSAU.\n", " The timestep is also used as maximal timestep for the integrator. \n", " It must be equal or smaller than the pulse duration.\n", " method: one of 'RK23', 'RK45', 'BDF', 'LSODA', 'Radau', 'DOP853'\n", " returns:\n", " res: return results object of scipy.integrate.solve_ivp\n", " in particular\n", " res['t'] time values\n", " res['y'] population values at different times\n", " \"\"\"\n", " confDict = {str(c): i for i, c in enumerate(configurations)}\n", " gammas = [{\n", " 'i': confDict[str(proc['iconfig'])],\n", " 'f': confDict[str(proc['fconfig'])],\n", " 'rate': proc['rate']\n", " } for proc in processData if proc['type'] == 'A' or proc['type'] == 'F']\n", " sigmas = [{\n", " 'i': confDict[str(proc['iconfig'])],\n", " 'f': confDict[str(proc['fconfig'])],\n", " 'pcs': proc['pcs']\n", " } for proc in processData if proc['type'] == 'P' or proc['type'] == 'V']\n", "\n", " def REQ(t, p):\n", " \"\"\" \n", " returns dp/dt (the right-hand-side of the diff. eq.)\n", " \"\"\"\n", " pDot = np.zeros(p.shape)\n", " Ipulse = pulseEnv(t) * fluence\n", " for s in sigmas:\n", " pDot[s['i']] += -s['pcs'] * Ipulse * p[s['i']]\n", " pDot[s['f']] += +s['pcs'] * Ipulse * p[s['i']]\n", " for g in gammas:\n", " pDot[g['i']] += -g['rate'] * p[g['i']]\n", " pDot[g['f']] += +g['rate'] * p[g['i']]\n", " return (pDot)\n", "\n", " def REQjac(t, p):\n", " \"\"\" \n", " returns d/dpj dp_i/dt (jacobian of the right-hand-side of the diff. eq.)\n", " \"\"\"\n", " pDot = REQ(t, p)\n", " pDotjac = np.zeros((p.shape[0], p.shape[0]))\n", " Ipulse = pulseEnv(t) * fluence\n", " for s in sigmas:\n", " pDotjac[s['i'], s['i']] += -s['pcs'] * Ipulse\n", " pDotjac[s['f'], s['i']] += +s['pcs'] * Ipulse\n", " for g in gammas:\n", " pDotjac[g['i'], g['i']] += -g['rate']\n", " pDotjac[g['f'], g['i']] += +g['rate']\n", " return (pDotjac)\n", "\n", " times = np.arange(tStart, tEnd, dt)\n", " if method == 'RK45' or method == 'RK23' or method == 'Radau' or method == 'DOP853':\n", " res = solve_ivp(REQ, (tStart - dt, tEnd + dt),\n", " pIni,\n", " t_eval=times,\n", " method=method,\n", " max_step=dt)\n", " elif method == 'LSODA' or method == 'BDF':\n", " res = solve_ivp(REQ, (tStart - dt, tEnd + dt),\n", " pIni,\n", " t_eval=times,\n", " method='LSODA',\n", " jac=REQjac,\n", " max_step=dt)\n", " else:\n", " raise ValueError(f\"method={method} unknown\")\n", " if not res.success:\n", " print(res.message)\n", " return res" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T08:43:54.643959Z", "start_time": "2022-05-31T08:43:54.490269Z" } }, "outputs": [], "source": [ "def readProcessData(filename):\n", " \"\"\"\n", " reads in and returns process data generated by xatom from file filename\n", " \"\"\"\n", " processData = []\n", " with open(filename, 'r') as f:\n", " lines = f.readlines()\n", " for line in lines:\n", " line = line.rstrip()\n", " els = line.split()\n", " if len(els) == 0:\n", " continue\n", " configs = line.split('#')[1].split('->')\n", " if (els[0] == 'P' or els[0] == 'V'):\n", " processData.append({\n", " \"type\":\n", " els[0],\n", " \"i\":\n", " els[1],\n", " \"f\":\n", " els[2],\n", " \"pcs\":\n", " float(els[3]),\n", " \"energy\":\n", " els[4],\n", " \"iconfig\":\n", " np.array([int(n) for n in configs[0].split()]),\n", " \"fconfig\":\n", " np.array([int(n) for n in configs[1].split()]),\n", " })\n", " else:\n", " processData.append({\n", " \"type\":\n", " els[0],\n", " \"i\":\n", " els[1],\n", " \"f\":\n", " els[2],\n", " \"rate\":\n", " float(els[3]),\n", " \"energy\":\n", " els[4],\n", " \"iconfig\":\n", " np.array([int(n) for n in configs[0].split()]),\n", " \"fconfig\":\n", " np.array([int(n) for n in configs[1].split()]),\n", " })\n", " return processData" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T08:43:54.768006Z", "start_time": "2022-05-31T08:43:54.647046Z" } }, "outputs": [], "source": [ "def gaussian(x, sigma=10 * FS2AU, x0=0.):\n", " \"\"\"\n", " Gaussian function:\n", " x0 is the center\n", " sigma is the stdev\n", " \"\"\"\n", " return 1. / np.sqrt(2 * np.pi) * 1 / sigma * np.exp(-0.5 *\n", " ((x - x0) / sigma)**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initial 3-configuration example" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T08:43:54.963549Z", "start_time": "2022-05-31T08:43:54.771246Z" } }, "outputs": [], "source": [ "pIni = np.array([1., 0., 0., 0.])\n", "sigma = 1.\n", "Gamma = 3 / 2\n", "processData = [\n", " {\n", " \"iconfig\": np.array([1, 1]),\n", " \"fconfig\": np.array([0, 1]),\n", " \"pcs\": sigma,\n", " \"type\": \"P\"\n", " },\n", " {\n", " \"iconfig\": np.array([0, 1]),\n", " \"fconfig\": np.array([0, 0]),\n", " \"pcs\": sigma,\n", " \"type\": \"P\"\n", " },\n", " {\n", " \"iconfig\": np.array([0, 1]),\n", " \"fconfig\": np.array([1, 0]),\n", " \"rate\": Gamma,\n", " \"type\": \"A\"\n", " },\n", "]\n", "configurations = [\n", " np.array([1, 1]),\n", " np.array([0, 1]),\n", " np.array([0, 0]),\n", " np.array([1, 0])\n", "]\n", "fluence = 1.\n", "pulseEnv = lambda t: 1. if t > 0 else 0.\n", "\n", "res = solveREQ(pIni,\n", " processData,\n", " configurations,\n", " fluence,\n", " pulseEnv,\n", " tStart=0,\n", " tEnd=10,\n", " dt=0.1,\n", " method='LSODA')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T08:43:55.758204Z", "start_time": "2022-05-31T08:43:54.966649Z" } }, "outputs": [], "source": [ "times = res['t']\n", "p = res['y']\n", "# plot populations and compare to analytical solution\n", "fig, ax = plt.subplots()\n", "ax.plot(times, p[0, :], label=r'$p_1$')\n", "ax.plot(times, p[1, :], label=r'$p_2$')\n", "ax.plot(times, p[2, :], label=r'$p_3$')\n", "ax.plot(times, p[3, :], label=r'$p_4$')\n", "\n", "ax.set_xlabel('time')\n", "ax.set_ylabel('population')\n", "ax.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3a put your code in the cells below" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-09T13:27:35.229312Z", "start_time": "2022-05-09T13:27:35.219365Z" } }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-09T13:27:38.267176Z", "start_time": "2022-05-09T13:27:38.099613Z" }, "scrolled": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Rate Equation Solving with Calculated Atomic Data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# read an preprocess the process data (rates and cross sections calculated by xatom)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T08:43:55.843318Z", "start_time": "2022-05-31T08:43:55.762042Z" } }, "outputs": [], "source": [ "# read in process data: neon, 2000eV\n", "processData = readProcessData('results_Ne_2000.dat')\n", "\n", "# get a list of involved electronic configurations\n", "configurations = np.unique(np.array([p['iconfig'] for p in processData] +\n", " [p['fconfig'] for p in processData]),\n", " axis=0)\n", "Nc = configurations.shape[0]\n", "\n", "# compute number of electrons for each configuration\n", "Nelectron = np.array([np.sum(c) for c in configurations])\n", "\n", "# initial population we start with the configuration that has the maximum number of electrons\n", "startCfgIdx = np.argmax(Nelectron)\n", "pIni = np.zeros(Nc)\n", "pIni[startCfgIdx] = 1.\n", "\n", "# Assume that nuclear charge is identical to maximum number of electron,\n", "# i.e., the configuration with the most electrons is neutral\n", "nuclearCharge = np.max(Nelectron)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# solve rate equations" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T09:05:50.527859Z", "start_time": "2022-05-31T09:05:49.309425Z" } }, "outputs": [], "source": [ "%%time\n", "\n", "# set the pulse parameters:\n", "# the fluence number of photons per micrometer square\n", "nPhotonsPerSqmuM = 1e12\n", "# the pulse width (FWHM)\n", "fwhmPulse = 25 * FS2AU\n", "# convert fluence into atomic units and pulse width FWHM into stdev (atomic units)\n", "fluence = nPhotonsPerSqmuM / (1e-6 * M2AU)**2\n", "sigmaPulse = fwhmPulse / (2. * np.sqrt(2. * np.log(2.)))\n", "# pulse envelope: simple Gaussian with stddev 25fs and t0 = 0fs\n", "pulseEnvelope = lambda t: gaussian(t, sigma=sigmaPulse, x0=0.0 * FS2AU)\n", "\n", "res = solveREQ(pIni,\n", " processData,\n", " configurations,\n", " fluence,\n", " pulseEnvelope,\n", " method='LSODA')\n", "\n", "# p is population, times is timepoints\n", "p = res['y']\n", "times = res['t']\n", "\n", "# compute the average charge per time\n", "avgCharge = nuclearCharge - np.dot(Nelectron, p)\n", "\n", "# compute the population in a given charge state per time\n", "chargePops = np.zeros((nuclearCharge + 1, p.shape[1]))\n", "for q in range(chargePops.shape[0]):\n", " idx = np.where(nuclearCharge - Nelectron == q)[0]\n", " chargePops[q, :] = np.sum(p[idx, :], axis=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# plot/analyze the results" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T09:07:21.453185Z", "start_time": "2022-05-31T09:07:21.242570Z" } }, "outputs": [], "source": [ "# plot average charge\n", "fig, ax = plt.subplots()\n", "ax2 = ax.twinx()\n", "\n", "ax.plot(times / FS2AU, avgCharge, label='avg. charge')\n", "ax2.plot(times / FS2AU,\n", " fluence * pulseEnvelope(times) * 1 / AU2FS * 1 / ((1e6 * AU2M)**2),\n", " label='pulse',\n", " c='black',\n", " alpha=0.3)\n", "\n", "ax.set_xlabel('time (fs)')\n", "ax.set_ylabel('avg. charge')\n", "ax2.set_ylabel(r'intensity (n. photons / $\\mathrm{\\mu m}^2$ fs)')\n", "ax.legend()\n", "ax2.legend()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T09:07:23.538517Z", "start_time": "2022-05-31T09:07:23.250195Z" } }, "outputs": [], "source": [ "# plot charge population\n", "fig, ax = plt.subplots()\n", "ax2 = ax.twinx()\n", "for q in range(chargePops.shape[0]):\n", " ax.plot(times / FS2AU, chargePops[q, :], label=f'q={q}')\n", "\n", "ax2.plot(times / FS2AU,\n", " fluence * pulseEnvelope(times) * 1 / AU2FS * 1 / ((1e6 * AU2M)**2),\n", " label='pulse',\n", " c='black',\n", " alpha=0.3)\n", "ax.legend()\n", "ax.set_xlabel('time (fs)')\n", "ax.set_ylabel('population')\n", "ax2.set_ylabel(r'intensity (n. photons / $\\mathrm{\\mu m}^2$ fs)')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-05-31T09:07:26.007685Z", "start_time": "2022-05-31T09:07:25.979855Z" } }, "outputs": [], "source": [ "# ratio of PP vs PA processes:\n", "# i.e. either A or P starting from configuration 1,2,6\n", "\n", "# determine the singly core ionized configuration\n", "cfgP = [\n", " p['fconfig'] for p in processData\n", " if np.all(p['iconfig'] == np.array(configurations[startCfgIdx]))\n", " and p['type'] == 'P'\n", "][0]\n", "print(\n", " f\"The singly core ionized configuration has electronic occupation numbers {cfgP}\"\n", ")\n", "\n", "# index of singly core ionized configuration\n", "idxP = np.argmax(np.all(configurations == cfgP, axis=1))\n", "\n", "# find the total photoionization c.s. and the Auger-Meitner decay rate for configuration idxP\n", "sigmaPP = np.sum([\n", " p['pcs'] for p in processData\n", " if np.all(p['iconfig'] == configurations[idxP]) and p['type'] == 'P'\n", "])\n", "GammaPA = np.sum([\n", " p['rate'] for p in processData\n", " if np.all(p['iconfig'] == configurations[idxP]) and p['type'] == 'A'\n", "])\n", "\n", "# probability for PP is\n", "# int dt sigma * F(t) p_P(t)\n", "# probability for PA is\n", "# int dt Gamma * p_P(t)\n", "dt = times[1] - times[0]\n", "probPP = np.sum(sigmaPP * fluence * pulseEnvelope(times) * p[idxP, :]) * dt\n", "probPA = np.sum(GammaPA * p[idxP, :]) * dt\n", "\n", "print(\n", " f\"The probability for PP is {probPP}.\\nThe probability for PA is {probPA}.\\nRatio PP/PA is {probPP/probPA}\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3b, 3c, 3d in the cells below" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }