{ "cells": [ { "cell_type": "markdown", "id": "d962cb38", "metadata": {}, "source": [ "# From an MD morphology to a state file\n", "\n", "VOTCA-XTP works on molecular morphologies from molecular dynamics (MD). The first step of every multiscale workflow is the same: the atoms of the MD snapshot are grouped into the molecules (*segments*) that the electronic-structure calculations work with, each segment is assigned a QM geometry and classical multipoles, and pairs of neighbouring segments are listed. The result is a *state file* (`state.hdf5`) that all later calculators read and extend.\n", "\n", "This tutorial does this step for an amorphous film of 1000 thiophene molecules simulated with LAMMPS. The tutorials *Charge transport in amorphous thiophene* and *Excitations of thiophene in its environment* start from the same morphology.\n", "\n", "**You will:** write a mapping file, map the MD snapshot, check the mapping, and build the neighbour list.\n", "\n", "## Requirements\n", "* VOTCA installed and the environment activated with `VOTCARC.bash` (see the [installation instructions](https://github.com/votca/votca/blob/master/share/sphinx/INSTALL.rst)).\n", "* Python with `numpy`, `h5py` and `matplotlib`.\n", "\n", "The VOTCA programs are command-line tools. Below they are called through a small helper `run()`, which prints only the last lines of their output, if any; in a terminal you would call them directly. At the end of their output, the programs list the references to cite for the methods they used." ] }, { "cell_type": "code", "execution_count": 1, "id": "4377d984", "metadata": { "execution": { "iopub.execute_input": "2026-11-10T20:59:41.346320Z", "iopub.status.busy": "2026-11-10T20:59:41.346131Z", "iopub.status.idle": "2026-11-10T20:59:41.696271Z", "shell.execute_reply": "2026-11-10T20:59:41.695736Z" } }, "outputs": [], "source": [ "import subprocess\n", "import numpy as np\n", "import h5py\n", "import matplotlib.pyplot as plt\n", "\n", "\n", "def run(command, show=0, output=False):\n", " \"\"\"Run a VOTCA command; print the last `show` lines of its output (0: none).\"\"\"\n", " result = subprocess.run(command, shell=True, text=True,\n", " stdout=subprocess.PIPE, stderr=subprocess.STDOUT)\n", " if result.returncode != 0:\n", " print(result.stdout)\n", " raise RuntimeError(f\"'{command}' failed\")\n", " if show:\n", " print(\"\\n\".join(result.stdout.splitlines()[-show:]))\n", " return result.stdout if output else None\n", "\n", "run(\"rm -f state.hdf5 skeleton.xml *.pdb\")" ] }, { "cell_type": "markdown", "id": "69810321", "metadata": {}, "source": [ "## The input\n", "\n", "The morphology consists of two LAMMPS files in `MD_FILES`:\n", "\n", "* `thiophene.data`, the *topology*: atoms with their types, masses and charges, and the bonds, which define the molecules;\n", "* `traj.dump`, the *trajectory*: atom positions (here a single snapshot at time step 10000).\n", "\n", "For GROMACS the topology is the run input file (`topol.tpr`) and the trajectory a `.gro`, `.trr` or `.xtc` file; everything below works the same way.\n", "\n", "In addition, each kind of molecule needs data from separate gas-phase calculations: optimized QM geometries for the neutral, anionic and cationic molecule (`QC_FILES/thiophene_{n,e,h}.xyz`) and atomic multipoles with polarizabilities for the same states (`MP_FILES/thiophene_{n,e,h}.mps`).\n", "\n", "## The mapping file\n", "\n", "The mapping file tells `xtp_map` which MD atoms form a segment and how they correspond to the atoms of the QM geometry and of the multipole file. The names it must use are those of the MD topology. `xtp_map -m` writes a skeleton with these names:" ] }, { "cell_type": "code", "execution_count": 2, "id": "531b5c32", "metadata": { "execution": { "iopub.execute_input": "2026-11-10T20:59:41.698952Z", "iopub.status.busy": "2026-11-10T20:59:41.698321Z", "iopub.status.idle": "2026-11-10T20:59:41.847063Z", "shell.execute_reply": "2026-11-10T20:59:41.846498Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\t\n", "\t\t\n", "\t\t\tC4H4S1\n", "\t\t\t\n", "\t\t\t\t\n", "\t\t\t\t\tUPTOYOU_BUTUNIQUE\n", "\t\t\t\t\tXYZFILE_GROUNDSTATE\n", "\t\t\t\t\tMPSFILE_GROUNDSTATE\n", "\t\t\t\t\tWANTTOMAPTOMDGEOMETRY\n", "\t\t\t\t\tREORG1_hole\n", "\t\t\t\t\tREORG2_hole\n", "\t\t\t\t\tREORG3_hole\n", "\t\t\t\t\t\n", "\t\t\t\t\t\t\n", "\t\t\t\t\t\t\tUPTOYOU_BUTUNIQUE\n", "\t\t\t\t\t\t\t 0:C63:0 0:C63:1 0:S23:2 0:C63:3 0:C63:4 0:H20:5 0:H20:6 0:H20:7 0:H20:8\n", "\t\t\t\t\t\t\tIDS of QMATOMS i.e 0:C 1:H 2:C\n", "\t\t\t\t\t\t\tIDS of MPOLES i.e 0:C 1:H 2:C\n", "\t\t\t\t\t\t\tweights for mapping(often atomic mass) i.e. 12 1 12\n", "\t\t\t\t\t\t\tIDs of up to 3 qmatoms or mpoles i.e. 0 1 2\n", "\t\t\t\t\t\t\n", "\t\t\t\t\t\n", "\t\t\t\t\n", "\t\t\t\n", "\t\t\n", "\t\n", "\n", "\n" ] } ], "source": [ "run(\"xtp_map -t MD_FILES/thiophene.data -c MD_FILES/traj.dump -s skeleton.xml -m\")\n", "print(open(\"skeleton.xml\").read())" ] }, { "cell_type": "markdown", "id": "7ad41de9", "metadata": {}, "source": [ "The LAMMPS reader names the molecule by its composition (`C4H4S1`) and every atom by its element (recognized from the mass) followed by its LAMMPS atom type (`C63`, `S23`, `H20`). An MD atom is written as `residue:name:index`, where the index counts the atoms within the molecule. With GROMACS the molecule and atom names of the topology are used.\n", "\n", "The complete mapping file `system.xml` fills in the placeholders:" ] }, { "cell_type": "code", "execution_count": 3, "id": "8bf106df", "metadata": { "execution": { "iopub.execute_input": "2026-11-10T20:59:41.849225Z", "iopub.status.busy": "2026-11-10T20:59:41.848819Z", "iopub.status.idle": "2026-11-10T20:59:41.852532Z", "shell.execute_reply": "2026-11-10T20:59:41.852020Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", " \n", " \n", " C4H4S1\n", " \n", " \n", " thiophene\n", " QC_FILES/thiophene_n.xyz\n", " QC_FILES/thiophene_e.xyz\n", " QC_FILES/thiophene_h.xyz\n", " MP_FILES/thiophene_n.mps\n", " MP_FILES/thiophene_e.mps\n", " MP_FILES/thiophene_h.mps\n", " 0\n", " 8.7140\n", " 0.1908\n", " 0.5711\n", " 1.5972\n", " 0.1940\n", " 0.1879\n", " \n", " \n", " thiophene\n", " 0:C63:0 0:C63:1 0:S23:2 0:C63:3 0:C63:4 0:H20:5 0:H20:6 0:H20:7 0:H20:8\n", " 0:C 1:C 2:S 3:C 4:C 5:H 6:H 7:H 8:H \n", " 0:C 1:C 2:S 3:C 4:C 5:H 6:H 7:H 8:H\n", " 12 12 32 12 12 1 1 1 1\n", " 1 2 3\n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "\n" ] } ], "source": [ "print(open(\"system.xml\").read())" ] }, { "cell_type": "markdown", "id": "9b560bc6", "metadata": {}, "source": [ "* A **segment** is the unit of the electronic-structure calculations; here one segment per molecule. A large molecule can be split into several segments, e.g. one per conjugated unit of a polymer.\n", "* `qmcoords_*` and `multipoles_*` give the QM geometry and the multipoles for each state: `n` neutral, `e` electron (anion), `h` hole (cation).\n", "* A **fragment** is a rigid part of a segment. `mdatoms`, `qmatoms` and `mpoles` list corresponding atoms in the same order; the QM geometry is placed onto the MD atoms of each fragment by a rigid-body fit, using the `weights` (here the masses) and the local frame defined by the atoms in `localframe`. With `map2md` set to 1, the QM atoms are put directly onto the MD positions instead.\n", "* `U_xX_nN_*`, `U_nX_nN_*` and `U_xN_xX_*` are internal energies (in eV) of the molecule in the charged state from gas-phase calculations; they enter the site energies and reorganization energies of the charge-transport tutorial.\n", "\n", "## Mapping the snapshot\n", "\n", "`xtp_map` reads the topology and trajectory, applies the mapping, and writes the state file:" ] }, { "cell_type": "code", "execution_count": 4, "id": "f993c2b0", "metadata": { "execution": { "iopub.execute_input": "2026-11-10T20:59:41.854984Z", "iopub.status.busy": "2026-11-10T20:59:41.854481Z", "iopub.status.idle": "2026-11-10T20:59:42.375150Z", "shell.execute_reply": "2026-11-10T20:59:42.374589Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Creating statefile state.hdf5\n", "Frame with id 10000 was not in statefile state.hdf5 ,adding it now.\n", "Wrote MD topology (step = 10000, time = 0) to state.hdf5\n", "... . \n" ] } ], "source": [ "run(\"xtp_map -t MD_FILES/thiophene.data -c MD_FILES/traj.dump -s system.xml -f state.hdf5\", show=4)" ] }, { "cell_type": "markdown", "id": "a9fd5bcb", "metadata": {}, "source": [ "The state file is an HDF5 file with one group per frame; each segment stores its atoms with the MD positions (in bohr):" ] }, { "cell_type": "code", "execution_count": 5, "id": "4feb4b5c", "metadata": { "execution": { "iopub.execute_input": "2026-11-10T20:59:42.377651Z", "iopub.status.busy": "2026-11-10T20:59:42.377102Z", "iopub.status.idle": "2026-11-10T20:59:42.385649Z", "shell.execute_reply": "2026-11-10T20:59:42.385068Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "box diagonal (nm): [5.15643722 5.15643722 5.15643722]\n", "segments: 1000\n", "atoms of segment 0: C63 C63 S23 C63 C63 H20 H20 H20 H20\n" ] } ], "source": [ "with h5py.File(\"state.hdf5\", \"r\") as f:\n", " frame = f[\"frame_10000\"]\n", " print(\"box diagonal (nm):\", np.diag(frame[\"box\"][()]) * 0.0529177)\n", " segments = frame[\"segments\"]\n", " print(\"segments:\", len(segments))\n", " atoms = segments[\"segment0\"][\"atoms\"][()].ravel()\n", " print(\"atoms of segment 0:\", \" \".join(a.decode() for a in atoms[\"name\"]))" ] }, { "cell_type": "markdown", "id": "0a317d42", "metadata": {}, "source": [ "## Checking the mapping\n", "\n", "`mapchecker` writes the segments as they are in the MD (`md_segments_*.pdb`), with the QM geometries mapped onto them (`qm_segments_*.pdb`), and with the multipoles of the charged states (`mp_segments_*.pdb`). Viewing them in e.g. VMD on top of each other is the quickest check that the atoms correspond. Here we compare the mapped QM geometries with the MD geometry numerically: since the QM geometry is fitted as a rigid body, the deviation measures how far the molecules in the MD differ from the optimized gas-phase geometry." ] }, { "cell_type": "code", "execution_count": 6, "id": "0d080b13", "metadata": { "execution": { "iopub.execute_input": "2026-11-10T20:59:42.387973Z", "iopub.status.busy": "2026-11-10T20:59:42.387566Z", "iopub.status.idle": "2026-11-10T20:59:43.081100Z", "shell.execute_reply": "2026-11-10T20:59:43.080581Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD between MD and mapped QM geometry: mean 0.180 A, max 0.372 A\n" ] } ], "source": [ "run(\"xtp_run -e mapchecker -c map_file=system.xml -f state.hdf5\")\n", "\n", "\n", "def read_pdb(filename):\n", " return np.array([[float(line[30:38]), float(line[38:46]), float(line[46:54])]\n", " for line in open(filename) if line.startswith((\"ATOM\", \"HETATM\"))])\n", "\n", "\n", "md = read_pdb(\"md_segments_step_10000.pdb\").reshape(-1, 9, 3)\n", "qm = read_pdb(\"qm_segments_n_step_10000.pdb\").reshape(-1, 9, 3)\n", "rmsd = np.sqrt(np.mean(np.sum((md - qm)**2, axis=2), axis=1))\n", "print(f\"RMSD between MD and mapped QM geometry: mean {rmsd.mean():.3f} A, max {rmsd.max():.3f} A\")" ] }, { "cell_type": "markdown", "id": "79a6424f", "metadata": {}, "source": [ "A wrong atom order in the mapping file shows up as a large deviation (several Å) for every molecule.\n", "\n", "## Neighbour list\n", "\n", "Charge or energy transfer is computed for pairs of segments that are close to each other. `neighborlist` lists all pairs with at least one pair of atoms closer than the cutoff (here 0.6 nm):" ] }, { "cell_type": "code", "execution_count": 7, "id": "ee94a292", "metadata": { "execution": { "iopub.execute_input": "2026-11-10T20:59:43.083132Z", "iopub.status.busy": "2026-11-10T20:59:43.082753Z", "iopub.status.idle": "2026-11-10T20:59:43.715931Z", "shell.execute_reply": "2026-11-10T20:59:43.715349Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10151 pairs, 20.3 neighbours per molecule\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAANJpJREFUeJzt3Xtc1GXe//H3yElAQMBkIAm0kDynWCZtq60o65pW3t260cF2bcuHRpKaq+uaWi2obeq2pa0uqR2M7ruik1biVnhaW7XcPJsGHgqWzQxQWVC8fn/8cm5HoQQZBi9fz8djHo/m+l7fmc9nhpp31/c733EYY4wAAAAs1czbBQAAAHgSYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGq+3i6gKTh16pS+/vprhYSEyOFweLscAABwHowxKi8vV0xMjJo1q339hrAj6euvv1ZsbKy3ywAAAPVw8OBBtWnTptbthB1JISEh0vcvVmhoqLfLAQAA56GsrEyxsbGuz/HaEHYk16Gr0NBQwg4AABeZHzsFhROUAQCA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFbz9XYBQFMXP2n5ec8tnDnIo7UAAOqOlR0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKvx1XNckurydXIAwMWNlR0AAGA1wg4AALAaYQcAAFjNq2Fn+vTpcjgcbjen0+nabozR9OnTFRMTo8DAQPXt21fbt293e4zKykqlp6erVatWCg4O1pAhQ3To0CEvdAMAAJoir6/sdOrUSUVFRa7b1q1bXdtmz56tOXPm6JlnntHGjRvldDrVv39/lZeXu+ZkZGQoNzdXOTk5Wrt2rY4ePaqbb75Z1dXVXuoIAAA0JV7/Npavr6/bas5pxhjNmzdPU6ZM0dChQyVJS5cuVVRUlJYtW6YHHnhApaWlys7O1osvvqiUlBRJ0ksvvaTY2FitWrVKqampjd4PAABoWry+svPFF18oJiZGbdu21S9/+Ut9+eWXkqSCggIVFxdrwIABrrkBAQHq06eP1q9fL0navHmzTpw44TYnJiZGnTt3ds2pSWVlpcrKytxuAADATl4NO7169dILL7ygDz74QIsWLVJxcbGSk5N1+PBhFRcXS5KioqLc9omKinJtKy4ulr+/v8LDw2udU5OsrCyFhYW5brGxsR7pDwAAeJ9Xw87AgQP1X//1X+rSpYtSUlK0fPn/v9Db0qVLXXMcDofbPsaYc8bO9mNzJk+erNLSUtft4MGDF9wLAABomrx+GOtMwcHB6tKli7744gvXeTxnr9CUlJS4VnucTqeqqqp05MiRWufUJCAgQKGhoW43AABgpyYVdiorK7Vz505FR0erbdu2cjqdysvLc22vqqpSfn6+kpOTJUlJSUny8/Nzm1NUVKRt27a55gAAgEubV7+NNWHCBA0ePFhXXHGFSkpK9MQTT6isrEwjRoyQw+FQRkaGMjMzlZCQoISEBGVmZiooKEhpaWmSpLCwMI0cOVLjx49XZGSkIiIiNGHCBNdhMaCx1fU3twpnDvJYLQCA/8+rYefQoUO644479M033+iyyy7T9ddfrw0bNiguLk6SNHHiRFVUVGj06NE6cuSIevXqpZUrVyokJMT1GHPnzpWvr6+GDRumiooK9evXT0uWLJGPj48XOwMAAE2FwxhjvF2Et5WVlSksLEylpaWcv3OJaCq/es7KDgDU3/l+fjepc3YAAAAaGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Xy9XQDQUOInLfd2CQCAJoiVHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1bjODuBFdbk2UOHMQR6tBQBsxcoOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNV9vFwDg/MRPWn7ecwtnDvJoLQBwMWkyKztZWVlyOBzKyMhwjRljNH36dMXExCgwMFB9+/bV9u3b3farrKxUenq6WrVqpeDgYA0ZMkSHDh3yQgcAAKApahJhZ+PGjVq4cKG6du3qNj579mzNmTNHzzzzjDZu3Cin06n+/furvLzcNScjI0O5ubnKycnR2rVrdfToUd18882qrq72QicAAKCp8XrYOXr0qO68804tWrRI4eHhrnFjjObNm6cpU6Zo6NCh6ty5s5YuXarjx49r2bJlkqTS0lJlZ2frqaeeUkpKirp3766XXnpJW7du1apVq7zYFQAAaCq8HnbGjBmjQYMGKSUlxW28oKBAxcXFGjBggGssICBAffr00fr16yVJmzdv1okTJ9zmxMTEqHPnzq45NamsrFRZWZnbDQAA2MmrJyjn5OTo008/1caNG8/ZVlxcLEmKiopyG4+KitL+/ftdc/z9/d1WhE7POb1/TbKysjRjxowG6gIAADRlXgs7Bw8e1NixY7Vy5Uo1b9681nkOh8PtvjHmnLGz/dicyZMna9y4ca77ZWVlio2NrVP98Ly6fPsIAIDaeO0w1ubNm1VSUqKkpCT5+vrK19dX+fn5evrpp+Xr6+ta0Tl7haakpMS1zel0qqqqSkeOHKl1Tk0CAgIUGhrqdgMAAHbyWtjp16+ftm7dqi1btrhuPXv21J133qktW7aoXbt2cjqdysvLc+1TVVWl/Px8JScnS5KSkpLk5+fnNqeoqEjbtm1zzQEAAJc2rx3GCgkJUefOnd3GgoODFRkZ6RrPyMhQZmamEhISlJCQoMzMTAUFBSktLU2SFBYWppEjR2r8+PGKjIxURESEJkyYoC5dupxzwjMAALg0NekrKE+cOFEVFRUaPXq0jhw5ol69emnlypUKCQlxzZk7d658fX01bNgwVVRUqF+/flqyZIl8fHy8WjsAAGgaHMYY4+0ivK2srExhYWEqLS3l/J0mhBOU64+fiwBwKTjfz2+vX2cHAADAkwg7AADAaoQdAABgNcIOAACwGmEHAABYrUl/9RwALlV1+TYi374DfhgrOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKv5ersAAEDjip+0/LznFs4c5NFagMbAyg4AALAaKzsAcJGry0oNcCliZQcAAFiNlR0AnMOBWvG3ARuwsgMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAq/HbWICF+BVsAPg/rOwAAACrEXYAAIDVCDsAAMBqhB0AAGA1TlBGo+LEWQBAY2NlBwAAWI2wAwAArMZhLABoJBzGBbyDlR0AAGA1wg4AALAaYQcAAFjNq2FnwYIF6tq1q0JDQxUaGqrevXvrvffec203xmj69OmKiYlRYGCg+vbtq+3bt7s9RmVlpdLT09WqVSsFBwdryJAhOnTokBe6AQAATZFXw06bNm00c+ZMbdq0SZs2bdLPfvYz3XLLLa5AM3v2bM2ZM0fPPPOMNm7cKKfTqf79+6u8vNz1GBkZGcrNzVVOTo7Wrl2ro0eP6uabb1Z1dbUXOwMAAE2FV8PO4MGD9Ytf/ELt27dX+/bt9Yc//EEtWrTQhg0bZIzRvHnzNGXKFA0dOlSdO3fW0qVLdfz4cS1btkySVFpaquzsbD311FNKSUlR9+7d9dJLL2nr1q1atWpVrc9bWVmpsrIytxsAALBTkzlnp7q6Wjk5OTp27Jh69+6tgoICFRcXa8CAAa45AQEB6tOnj9avXy9J2rx5s06cOOE2JyYmRp07d3bNqUlWVpbCwsJct9jYWA93BwAAvMXr19nZunWrevfurf/85z9q0aKFcnNz1bFjR1dYiYqKcpsfFRWl/fv3S5KKi4vl7++v8PDwc+YUFxfX+pyTJ0/WuHHjXPfLysoIPABwgepyHaHCmYM8WgtwJq+HncTERG3ZskXfffedXn/9dY0YMUL5+fmu7Q6Hw22+MeacsbP92JyAgAAFBAQ0QPUAAKCp8/phLH9/f1111VXq2bOnsrKy1K1bN/3pT3+S0+mUvl+9OVNJSYlrtcfpdKqqqkpHjhypdQ4AALi0eT3snM0Yo8rKSrVt21ZOp1N5eXmubVVVVcrPz1dycrIkKSkpSX5+fm5zioqKtG3bNtccAABwaavXYayDBw/K4XCoTZs2kqR//OMfWrZsmTp27Kj777//vB/nd7/7nQYOHKjY2FiVl5crJydHH3/8sd5//305HA5lZGQoMzNTCQkJSkhIUGZmpoKCgpSWliZJCgsL08iRIzV+/HhFRkYqIiJCEyZMUJcuXZSSklKf1gAAgGXqFXbS0tJ0//336+6771ZxcbH69++vTp066aWXXlJxcbEeffTR83qcf/3rX7r77rtVVFSksLAwde3aVe+//7769+8vSZo4caIqKio0evRoHTlyRL169dLKlSsVEhLieoy5c+fK19dXw4YNU0VFhfr166clS5bIx8enPq0B+BGchArgYuMwxpi67hQeHq4NGzYoMTFRTz/9tF599VWtW7dOK1eu1KhRo/Tll196ploPKSsrU1hYmEpLSxUaGurtcqzGrz5fWgg77vj7/z/8baAhnO/nd71Wdk6cOOH6NtOqVas0ZMgQSdLVV1+toqKi+tYMAB7hqZDBBzZwcahX2OnUqZOee+45DRo0SHl5eXr88cclSV9//bUiIyMbukYAFykOeQFoCuoVdmbNmqXbbrtNTz75pEaMGKFu3bpJkt5++21dd911DV0jAMAydV1tIwzjQtQ57Bhj1LZtW+3fv1/V1dVuVy++//77FRQU1NA1AkCTxDk4wMWhztfZMcYoISFB//rXv875mYb4+Hi1bt26IesDAAC4IHVe2WnWrJkSEhJ0+PBhJSQkeKYqAPgRrKoAOF/1uoLy7Nmz9cgjj2jbtm0NXxEAAEADqtcJynfddZeOHz+ubt26yd/fX4GBgW7bv/3224aqDwAA4ILUK+zMmzev4SsBAADwgHqFnREjRjR8JQAAAB5w3mGnrKzMdSnmsrKyH5zLTy4AAICm4rzDTnh4uIqKitS6dWu1bNlSDofjnDnGGDkcDlVXVzd0nQAAAPVy3mHnww8/VEREhCTpo48+8mRNAAAADea8w06fPn1q/GcAAICmrF4nKJ92/PhxHThwQFVVVW7jXbt2vdC6AAAAGkS9ws6///1v/epXv9J7771X43bO2QEAAE1Fva6gnJGRoSNHjmjDhg0KDAzU+++/r6VLlyohIUFvv/12w1cJAABQT/Va2fnwww/11ltv6dprr1WzZs0UFxen/v37KzQ0VFlZWRo0aFDDVwoAAFAP9VrZOXbsmOvXzSMiIvTvf/9bktSlSxd9+umnDVshAADABajXyk5iYqJ2796t+Ph4XXPNNfrLX/6i+Ph4Pffcc4qOjm74KgFYj18xB+Ap9Qo7GRkZKioqkiRNmzZNqampeumll+Tv76+lS5c2dI0AAAD1Vq+wc+edd7r++ZprrlFhYaF27dqlK664Qq1atWrI+gAAAC5Ivc7ZkaTs7Gx17txZzZs3V3h4uO655x69+eabDVsdAADABarXys7UqVM1d+5cpaenq3fv3pKkv//973r44YdVWFioJ554oqHrBAAAqJd6hZ0FCxZo0aJFuuOOO1xjQ4YMUdeuXZWenk7YAQAATUa9DmNVV1erZ8+e54wnJSXp5MmTDVEXAABAg6hX2Lnrrru0YMGCc8YXLlzodvIyAACAt9X7h0Czs7O1cuVKXX/99ZKkDRs26ODBg7rnnns0btw417w5c+Y0TKUAAAD1UK+ws23bNvXo0UOStG/fPknSZZddpssuu0zbtm1zzXM4HA1VJwAAQL3UK+x89NFHDV8JAACAB9T7OjsAAAAXg3qfswMAQFNUl99ZK5w5yKO1oGkg7AAAmjx+KBYXgsNYAADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrcZ0dXDCufwEAaMpY2QEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYzathJysrS9dee61CQkLUunVr3Xrrrdq9e7fbHGOMpk+frpiYGAUGBqpv377avn2725zKykqlp6erVatWCg4O1pAhQ3To0KFG7gYAADRFXg07+fn5GjNmjDZs2KC8vDydPHlSAwYM0LFjx1xzZs+erTlz5uiZZ57Rxo0b5XQ61b9/f5WXl7vmZGRkKDc3Vzk5OVq7dq2OHj2qm2++WdXV1V7qDAAANBUOY4zxdhGn/fvf/1br1q2Vn5+vn/70pzLGKCYmRhkZGfrtb38rfb+KExUVpVmzZumBBx5QaWmpLrvsMr344osaPny4JOnrr79WbGysVqxYodTU1B993rKyMoWFham0tFShoaEe79M28ZOWe7sEAKiXwpmDvF0CLsD5fn43qXN2SktLJUkRERGSpIKCAhUXF2vAgAGuOQEBAerTp4/Wr18vSdq8ebNOnDjhNicmJkadO3d2zTlbZWWlysrK3G4AAMBOTSbsGGM0btw4/eQnP1Hnzp0lScXFxZKkqKgot7lRUVGubcXFxfL391d4eHitc86WlZWlsLAw1y02NtZDXQEAAG9rMmHnwQcf1Oeff65XXnnlnG0Oh8PtvjHmnLGz/dCcyZMnq7S01HU7ePDgBVYPAACaqiYRdtLT0/X222/ro48+Ups2bVzjTqdTOmOF57SSkhLXao/T6VRVVZWOHDlS65yzBQQEKDQ01O0GAADs5OvNJzfGKD09Xbm5ufr444/Vtm1bt+1t27aV0+lUXl6eunfvLkmqqqpSfn6+Zs2aJUlKSkqSn5+f8vLyNGzYMElSUVGRtm3bptmzZ3uhKwDAxaIuX7DgZOaLl1fDzpgxY7Rs2TK99dZbCgkJca3ghIWFKTAwUA6HQxkZGcrMzFRCQoISEhKUmZmpoKAgpaWlueaOHDlS48ePV2RkpCIiIjRhwgR16dJFKSkp3mwPAAA0AV4NOwsWLJAk9e3b12188eLFuvfeeyVJEydOVEVFhUaPHq0jR46oV69eWrlypUJCQlzz586dK19fXw0bNkwVFRXq16+flixZIh8fn0buCAAANDVN6jo73sJ1di4M19kBcCngMFbTc1FeZwcAAKChEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWM3X2wUAAHAxiJ+0vE7zC2cO8lgtqBtWdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALCar7cLAADARvGTlp/33MKZgzxay6WOlR0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAq/l6uwA0TfGTlnu7BAAAGgQrOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAVvNq2Fm9erUGDx6smJgYORwOvfnmm27bjTGaPn26YmJiFBgYqL59+2r79u1ucyorK5Wenq5WrVopODhYQ4YM0aFDhxq5EwAA0FR5NewcO3ZM3bp10zPPPFPj9tmzZ2vOnDl65plntHHjRjmdTvXv31/l5eWuORkZGcrNzVVOTo7Wrl2ro0eP6uabb1Z1dXUjdgIAAJoqr15UcODAgRo4cGCN24wxmjdvnqZMmaKhQ4dKkpYuXaqoqCgtW7ZMDzzwgEpLS5Wdna0XX3xRKSkpkqSXXnpJsbGxWrVqlVJTU2t87MrKSlVWVrrul5WVeaQ/AADgfU32nJ2CggIVFxdrwIABrrGAgAD16dNH69evlyRt3rxZJ06ccJsTExOjzp07u+bUJCsrS2FhYa5bbGysh7sBAADe0mTDTnFxsSQpKirKbTwqKsq1rbi4WP7+/goPD691Tk0mT56s0tJS1+3gwYMe6QEAAHhfk/9tLIfD4XbfGHPO2Nl+bE5AQIACAgIarEYAANB0NdmVHafTKZ2xwnNaSUmJa7XH6XSqqqpKR44cqXUOAAC4tDXZsNO2bVs5nU7l5eW5xqqqqpSfn6/k5GRJUlJSkvz8/NzmFBUVadu2ba45AADg0ubVw1hHjx7V3r17XfcLCgq0ZcsWRURE6IorrlBGRoYyMzOVkJCghIQEZWZmKigoSGlpaZKksLAwjRw5UuPHj1dkZKQiIiI0YcIEdenSxfXtLAAAcGnzatjZtGmTbrrpJtf9cePGSZJGjBihJUuWaOLEiaqoqNDo0aN15MgR9erVSytXrlRISIhrn7lz58rX11fDhg1TRUWF+vXrpyVLlsjHx8crPQEAUFfxk5af99zCmYM8WouNHMYY4+0ivK2srExhYWEqLS1VaGiot8tpEuryLx4AoPEQdv7P+X5+N9lzdgAAABoCYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNWa/G9joWHwVXIAwKWKlR0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAar7eLgAAAJy/+EnL6zS/cOYgj9VysWBlBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKv5ersA1F/8pOXeLgEA0MTV5bOicOYgj9biLazsAAAAqxF2AACA1Qg7AADAapyzAwAAJIvP72FlBwAAWI2VnSaGb1gBANCwWNkBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1a76NNX/+fD355JMqKipSp06dNG/ePN14443eLotvVwEA4GVWrOy8+uqrysjI0JQpU/TZZ5/pxhtv1MCBA3XgwAFvlwYAALzMYYwx3i7iQvXq1Us9evTQggULXGMdOnTQrbfeqqysrB/dv6ysTGFhYSotLVVoaGiD1sbKDgDgUuepqy2f7+f3RX8Yq6qqSps3b9akSZPcxgcMGKD169fXuE9lZaUqKytd90tLS6XvX7SGdqryeIM/JgAAFxNPfL6e+bg/tm5z0Yedb775RtXV1YqKinIbj4qKUnFxcY37ZGVlacaMGeeMx8bGeqxOAAAuVWHzPPv45eXlCgsLq3X7RR92TnM4HG73jTHnjJ02efJkjRs3znX/1KlT+vbbbxUZGVnrPk1FWVmZYmNjdfDgwQY/5NaUXYp90zM92+xS7JueG75nY4zKy8sVExPzg/Mu+rDTqlUr+fj4nLOKU1JScs5qz2kBAQEKCAhwG2vZsqVH62xooaGhl8y/LGe6FPum50vDpdizLtG+6blh/dCKzmkX/bex/P39lZSUpLy8PLfxvLw8JScne60uAADQNFz0KzuSNG7cON19993q2bOnevfurYULF+rAgQMaNWqUt0sDAABeZkXYGT58uA4fPqzHHntMRUVF6ty5s1asWKG4uDhvl9bgAgICNG3atHMOw9nuUuybni8Nl2LPukT7pmfvseI6OwAAALW56M/ZAQAA+CGEHQAAYDXCDgAAsBphBwAAWI2w0wTNnz9fbdu2VfPmzZWUlKQ1a9ac137r1q2Tr6+vrrnmGo/X2NDq0vPHH38sh8Nxzm3Xrl2NWnNDqOt7XVlZqSlTpiguLk4BAQG68sor9fzzzzdavQ2hLj3fe++9Nb7XnTp1atSaL1Rd3+eXX35Z3bp1U1BQkKKjo/WrX/1Khw8fbrR6G0pd+3722WfVoUMHBQYGKjExUS+88EKj1doQVq9ercGDBysmJkYOh0Nvvvnmj+6Tn5+vpKQkNW/eXO3atdNzzz3XKLU2lLr2XFRUpLS0NCUmJqpZs2bKyMhonEINmpScnBzj5+dnFi1aZHbs2GHGjh1rgoODzf79+39wv++++860a9fODBgwwHTr1q3R6m0Ide35o48+MpLM7t27TVFRket28uTJRq/9QtTnvR4yZIjp1auXycvLMwUFBeaTTz4x69ata9S6L0Rde/7uu+/c3uODBw+aiIgIM23atEavvb7q2vOaNWtMs2bNzJ/+9Cfz5ZdfmjVr1phOnTqZW2+9tdFrvxB17Xv+/PkmJCTE5OTkmH379plXXnnFtGjRwrz99tuNXnt9rVixwkyZMsW8/vrrRpLJzc39wflffvmlCQoKMmPHjjU7duwwixYtMn5+fua1115rtJovVF17LigoMA899JBZunSpueaaa8zYsWMbpU7CThNz3XXXmVGjRrmNXX311WbSpEk/uN/w4cPN73//ezNt2rSLLuzUtefTYefIkSONVKFn1LXv9957z4SFhZnDhw83UoUNr75/36fl5uYah8NhCgsLPVRhw6trz08++aRp166d29jTTz9t2rRp49E6G1pd++7du7eZMGGC29jYsWPNDTfc4NE6PeV8PvgnTpxorr76arexBx54wFx//fUers4zzqfnM/Xp06fRwg6HsZqQqqoqbd68WQMGDHAbHzBggNavX1/rfosXL9a+ffs0bdq0RqiyYdW3Z0nq3r27oqOj1a9fP3300UcerrRh1afvt99+Wz179tTs2bN1+eWXq3379powYYIqKioaqeoLcyHv9WnZ2dlKSUm5aC4YWp+ek5OTdejQIa1YsULGGP3rX//Sa6+9pkGDBjVS1ReuPn1XVlaqefPmbmOBgYH6xz/+oRMnTni0Xm/5+9//fs5rlJqaqk2bNlnbs7cQdpqQb775RtXV1ef8gGlUVNQ5P3R62hdffKFJkybp5Zdflq/vxXdB7Pr0HB0drYULF+r111/XG2+8ocTERPXr10+rV69upKovXH36/vLLL7V27Vpt27ZNubm5mjdvnl577TWNGTOmkaq+MPXp+UxFRUV67733dN9993mwyoZVn56Tk5P18ssva/jw4fL395fT6VTLli315z//uZGqvnD16Ts1NVV//etftXnzZhljtGnTJj3//PM6ceKEvvnmm0aqvHEVFxfX+BqdPHnS2p695eL7dLwEOBwOt/vGmHPGJKm6ulppaWmaMWOG2rdv34gVNrzz7VmSEhMTlZiY6Lrfu3dvHTx4UH/84x/105/+1OO1NqS69H3q1Ck5HA69/PLLrl/5nTNnjm6//XY9++yzCgwMbJSaL1Rdej7TkiVL1LJlS916660erM4z6tLzjh079NBDD+nRRx9VamqqioqK9Mgjj2jUqFHKzs5upIobRl36njp1qoqLi3X99dfLGKOoqCjde++9mj17tnx8fBqp4sZX02tU0zguDCs7TUirVq3k4+Nzzv/5lJSUnJP+Jam8vFybNm3Sgw8+KF9fX/n6+uqxxx7TP//5T/n6+urDDz9sxOrrp6491+b666/XF1984YEKPaM+fUdHR+vyyy93BR1J6tChg4wxOnTokMdrvlAX8l4bY/T888/r7rvvlr+/v4crbTj16TkrK0s33HCDHnnkEXXt2lWpqamaP3++nn/+eRUVFTVS5RemPn0HBgbq+eef1/Hjx1VYWKgDBw4oPj5eISEhatWqVSNV3ricTmeNr5Gvr68iIyO9VpeNCDtNiL+/v5KSkpSXl+c2npeXp+Tk5HPmh4aGauvWrdqyZYvrNmrUKCUmJmrLli3q1atXI1ZfP3XtuTafffaZoqOjPVChZ9Sn7xtuuEFff/21jh496hrbs2ePmjVrpjZt2ni85gt1Ie91fn6+9u7dq5EjR3q4yoZVn56PHz+uZs3c/9N8emXjYvkpwwt5r/38/NSmTRv5+PgoJydHN9988zmvhy169+59zmu0cuVK9ezZU35+fl6ry0qNcho0ztvpr2tmZ2ebHTt2mIyMDBMcHOz69smkSZPM3XffXev+F+O3sera89y5c01ubq7Zs2eP2bZtm5k0aZKRZF5//XUvdlF3de27vLzctGnTxtx+++1m+/btJj8/3yQkJJj77rvPi13UTX3/vu+66y7Tq1cvL1R84era8+LFi42vr6+ZP3++2bdvn1m7dq3p2bOnue6667zYRd3Vte/du3ebF1980ezZs8d88sknZvjw4SYiIsIUFBR4sYu6KS8vN5999pn57LPPjCQzZ84c89lnn7m+bn92z6e/ev7www+bHTt2mOzs7Ivuq+d17dkY45qflJRk0tLSzGeffWa2b9/u0ToJO03Qs88+a+Li4oy/v7/p0aOHyc/Pd20bMWKE6dOnT637Xoxhx9Sx51mzZpkrr7zSNG/e3ISHh5uf/OQnZvny5V6q/MLU9b3euXOnSUlJMYGBgaZNmzZm3Lhx5vjx416ovP7q2vN3331nAgMDzcKFC71QbcOoa89PP/206dixowkMDDTR0dHmzjvvNIcOHfJC5RemLn3v2LHDXHPNNSYwMNCEhoaaW265xezatctLldfP6ctinH0bMWKEMbW81x9//LHp3r278ff3N/Hx8WbBggVeqr5+6tNzTfPj4uI8WqfDXCzrogAAAPVg54FQAACA7xF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAqBdjjO6//35FRETI4XBoy5Yt3i6pwcTHx2vevHmu+w6HQ2+++aZXa/KUqqoqXXXVVVq3bp3Hn+vaa6/VG2+84fHnAc5G2AEuQmd/GHvD+++/ryVLlujdd99VUVGROnfu7NV6PKmoqEgDBw48r7kXWzBauHCh4uLidMMNN3j8uaZOnapJkybp1KlTHn8u4EyEHcBS1dXVHv1Q2bdvn6Kjo5WcnCyn0ylfX1+PPZe3OZ1OBQQEeLsMj/jzn/+s++67r1Gea9CgQSotLdUHH3zQKM8HnEbYARrYqVOnNGvWLF111VUKCAjQFVdcoT/84Q+u7V999ZWGDx+u8PBwRUZG6pZbblFhYaFr+7333qtbb71Vf/zjHxUdHa3IyEiNGTNGJ06ckCT17dtX+/fv18MPPyyHwyGHwyFJWrJkiVq2bKl3331XHTt2VEBAgPbv36+qqipNnDhRl19+uYKDg9WrVy99/PHHP9pHfn6+rrvuOgUEBCg6OlqTJk3SyZMnXTWmp6frwIEDcjgcio+Pr/ExTveSmZmpqKgotWzZUjNmzNDJkyf1yCOPKCIiQm3atNHzzz/vtt9vf/tbtW/fXkFBQWrXrp2mTp3q6l+S/vnPf+qmm25SSEiIQkNDlZSUpE2bNkmS9u/fr8GDBys8PFzBwcHq1KmTVqxYUWufJSUlGjx4sAIDA9W2bVu9/PLL58w5c7WmqqpKDz74oKKjo9W8eXPFx8crKytL+n7FTZJuu+02t9dl3759uuWWWxQVFaUWLVro2muv1apVq9yeIz4+XpmZmfr1r3+tkJAQXXHFFVq4cKHbnEOHDumXv/ylIiIiFBwcrJ49e+qTTz5xbX/nnXeUlJSk5s2bq127dq7Xujaffvqp9u7dq0GDBrnGCgsL5XA49MYbb+imm25SUFCQunXrpr///e+uOWf+rSUmJiooKEi33367jh07pqVLlyo+Pl7h4eFKT09XdXW1az8fHx/94he/0CuvvFJrTYBHePRnRoFL0MSJE014eLhZsmSJ2bt3r1mzZo1ZtGiRMcaYY8eOmYSEBPPrX//afP7552bHjh0mLS3NJCYmmsrKSmO+/5Xg0NBQM2rUKLNz507zzjvvmKCgINevfh8+fNi0adPGPPbYY6aoqMgUFRUZY4xZvHix8fPzM8nJyWbdunVm165d5ujRoyYtLc0kJyeb1atXm71795onn3zSBAQEmD179tTaw6FDh0xQUJAZPXq02blzp8nNzTWtWrUy06ZNM+b7XyJ/7LHHTJs2bUxRUZEpKSmp8XFGjBhhQkJCzJgxY8yuXbtMdna2kWRSU1PNH/7wB7Nnzx7z+OOPGz8/P3PgwAHXfo8//rhZt26dKSgoMG+//baJiooys2bNcm3v1KmTueuuu8zOnTvNnj17zP/8z/+YLVu2GGOMGTRokOnfv7/5/PPPzb59+8w777zj9mvbZxs4cKDp3LmzWb9+vdm0aZNJTk42gYGBZu7cua45kkxubq4xxpgnn3zSxMbGmtWrV5vCwkKzZs0as2zZMmOMMSUlJUaSWbx4sdvrsmXLFvPcc8+Zzz//3OzZs8dMmTLFNG/e3Ozfv9/1HHFxcSYiIsI8++yz5osvvjBZWVmmWbNmZufOncYYY8rLy027du3MjTfeaNasWWO++OIL8+qrr5r169cbY4x5//33TWhoqFmyZInZt2+fWblypYmPjzfTp0+vtfe5c+eaq6++2m2soKDASDJXX321effdd83u3bvN7bffbuLi4syJEyeMOeNvrX///ubTTz81+fn5JjIy0gwYMMAMGzbMbN++3bzzzjvG39/f5OTkuD3+/PnzTXx8fK01AZ5A2AEaUFlZmQkICHCFm7NlZ2ebxMREc+rUKddYZWWlCQwMNB988IEx3weEuLg4c/LkSdec//7v/zbDhw933Y+Li3P7MDbffwBJcn3oG2PM3r17jcPhMF999ZXb3H79+pnJkyfX2sfvfve7c+p89tlnTYsWLUx1dbUx339QxsXF/eDrcbqX0/sYY0xiYqK58cYbXfdPnjxpgoODzSuvvFLr48yePdskJSW57oeEhJglS5bUOLdLly4/+AF/pt27dxtJZsOGDa6xnTt3Gkm1hp309HTzs5/9zO21OdOZc39Ix44dzZ///GfX/bi4OHPXXXe57p86dcq0bt3aLFiwwBhjzF/+8hcTEhJiDh8+XOPj3XjjjSYzM9Nt7MUXXzTR0dG11jB27Fjzs5/9zG3sdNj561//6hrbvn27keQKXqf/1vbu3eua88ADD5igoCBTXl7uGktNTTUPPPCA2+O/9dZbplmzZm5/E4CncRgLaEA7d+5UZWWl+vXrV+P2zZs3a+/evQoJCVGLFi3UokULRURE6D//+Y/27dvnmtepUyf5+Pi47kdHR6ukpORHn9/f319du3Z13f/0009ljFH79u1dz9eiRQvl5+e7nu/M8VGjRrn66N27t+sQmSTdcMMNOnr0qA4dOnTO8x44cMDtcTIzM916adbs//5TExUVpS5durju+/j4KDIy0q2/1157TT/5yU/kdDrVokULTZ06VQcOHHBtHzdunO677z6lpKRo5syZbq/dQw89pCeeeEI33HCDpk2bps8//7zW12vnzp3y9fVVz549XWNXX321WrZsWes+9957r7Zs2aLExEQ99NBDWrlyZa1zTzt27JgmTpyojh07qmXLlmrRooV27drl1pMkt/fO4XDI6XS6XpctW7aoe/fuioiIqPE5Nm/erMcee8ztffjNb36joqIiHT9+vMZ9Kioq1Lx58xq3nVlLdHS09P0hv9OCgoJ05ZVXuu5HRUUpPj5eLVq0cBs7++82MDBQp06dUmVlZS2vFtDw7D2jEPCCwMDAH9x+6tQpJSUl1XheyGWXXeb6Zz8/P7dtDofjvE42DgwMdAsop06dko+PjzZv3uwWnvR9yNH3H6KnhYaGSt9/rfzMxzk9drqWs8XExLg9zpkfyDX18kP9bdiwQb/85S81Y8YMpaamKiwsTDk5OXrqqadc86dPn660tDQtX75c7733nqZNm6acnBzddtttuu+++5Samqrly5dr5cqVysrK0lNPPaX09PRz6v6hnmrTo0cPFRQU6L333tOqVas0bNgwpaSk6LXXXqt1n0ceeUQffPCB/vjHP+qqq65SYGCgbr/9dlVVVbnN+6HX5Xz+tmbMmKGhQ4ees622QNOqVStt3bq1xm1n1nL69Tnzb7Cu7+tp3377rYKCgn60H6AhEXaABpSQkKDAwED97W9/q/EbLj169NCrr76q1q1bu4JFffj7+7ud+Fmb7t27q7q6WiUlJbrxxhtrnHPVVVedM9axY0e9/vrrbqFn/fr1CgkJ0eWXX37OfF9f3xofpz7WrVunuLg4TZkyxTW2f//+c+a1b99e7du318MPP6w77rhDixcv1m233SZJio2N1ahRozRq1ChNnjxZixYtqjHsdOjQQSdPntSmTZt03XXXSZJ2796t77777gdrDA0N1fDhwzV8+HDdfvvt+vnPf65vv/1WERER8vPzO+e9WbNmje69915XfUePHnU7Kf18dO3aVX/9619dz3O2Hj16aPfu3XV6H7p3764FCxbUGG49Zdu2berRo0ejPBdwGoexgAbUvHlz/fa3v9XEiRP1wgsvaN++fdqwYYOys7MlSXfeeadatWqlW265RWvWrFFBQYHy8/M1duzYGg8P1SY+Pl6rV6/WV199pW+++abWee3bt9edd96pe+65R2+88YYKCgq0ceNGzZo16we/oTR69GgdPHhQ6enp2rVrl9566y1NmzZN48aNczsk5QlXXXWVDhw4oJycHO3bt09PP/20cnNzXdsrKir04IMP6uOPP9b+/fu1bt06bdy4UR06dJAkZWRk6IMPPlBBQYE+/fRTffjhh65tZ0tMTNTPf/5z/eY3v9Enn3yizZs367777vvBVYe5c+cqJydHu3bt0p49e/S///u/cjqdrkNf8fHx+tvf/qbi4mIdOXLE1dMbb7yhLVu26J///KfS0tLqfFmAO+64Q06nU7feeqvWrVunL7/8Uq+//rrrW1KPPvqoXnjhBU2fPl3bt2/Xzp079eqrr+r3v/99rY9500036dixY9q+fXudarkQa9as0YABAxrt+QARdoCGN3XqVI0fP16PPvqoOnTooOHDh7vOWwgKCtLq1at1xRVXaOjQoerQoYN+/etfq6Kiok4rPY899pgKCwt15ZVXuh3+qsnixYt1zz33aPz48UpMTNSQIUP0ySefKDY2ttZ9Lr/8cq1YsUL/+Mc/1K1bN40aNUojR478wQ/OhnLLLbfo4Ycf1oMPPqhrrrlG69ev19SpU13bfXx8dPjwYd1zzz1q3769hg0bpoEDB2rGjBnS99cXGjNmjDp06KCf//znSkxM1Pz582t9vsWLFys2NlZ9+vTR0KFDdf/996t169a1zm/RooVmzZqlnj176tprr1VhYaFWrFjhCoFPPfWU8vLyFBsbq+7du0vfB6Tw8HAlJydr8ODBSk1NrfPqhr+/v1auXKnWrVvrF7/4hbp06aKZM2e6Dk+mpqbq3XffVV5enq699lpdf/31mjNnjuLi4mp9zMjISA0dOrTGw6qe8NVXX2n9+vX61a9+1SjPB5zmMKcPWgMALjlbt25VSkqK68R5T3rkkUdUWlp6zvWDAE9jZQcALmFdunTR7Nmz63wOUX20bt1ajz/+uMefBzgbKzsAAMBqrOwAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKv9P6J740QxfBSlAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "run(\"xtp_run -e neighborlist -c constant=0.6 -f state.hdf5\", show=0)\n", "with h5py.File(\"state.hdf5\", \"r\") as f:\n", " pairs = f[\"frame_10000\"][\"neighborlist\"][\"pairs\"][()]\n", "distance = np.sqrt(pairs[\"delta_Rx\"]**2 + pairs[\"delta_Ry\"]**2 + pairs[\"delta_Rz\"]**2) * 0.0529177\n", "print(f\"{len(pairs)} pairs, {2 * len(pairs) / 1000:.1f} neighbours per molecule\")\n", "plt.hist(distance, bins=40)\n", "plt.xlabel(\"centre-of-mass distance (nm)\")\n", "plt.ylabel(\"pairs\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b01bda5f", "metadata": {}, "source": [ "The cutoff acts on the closest atoms, so the centre-of-mass distances of the pairs extend beyond it. Whether the cutoff is large enough can be checked later from the couplings: those of the most distant pairs should be negligible.\n", "\n", "## Summary\n", "\n", "* The mapping file connects the MD topology to the QM geometries and multipoles of each segment; `xtp_map -m` writes a skeleton with the names of the MD topology.\n", "* `xtp_map` writes the state file; `mapchecker` checks the mapping; `neighborlist` adds the pairs.\n", "* The state file is the input for all further steps: site energies, couplings and transport (*Charge transport in amorphous thiophene*) or excited states in the environment (*Excitations of thiophene in its environment*)." ] } ], "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.14.8" } }, "nbformat": 4, "nbformat_minor": 5 }