{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7fb27b941602401d91542211134fc71a",
   "metadata": {},
   "source": [
    "## Diffusion coefficient temperature dependence\n",
    "\n",
    "The diffusion coefficient often follows an Arrhenius relation with temperature. This example runs Lennard-Jones argon simulations at several temperatures in parallel, extracts the diffusion coefficient from each trajectory, and fits ``ln(D)`` versus ``1/T`` to estimate the prefactor and activation barrier.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9a63283cbaf04dbcab1f6479b197f3a8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scm.plams import Settings\n",
    "\n",
    "\n",
    "def get_lennard_jones_settings():\n",
    "    s = Settings()\n",
    "    s.input.lennardjones.eps = 3.785e-4\n",
    "    s.input.lennardjones.rmin = 3.81637\n",
    "    s.input.lennardjones.cutoff = 8.0\n",
    "    s.runscript.nproc = 1\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dd0d8092fe74a7c96281538738b07e2",
   "metadata": {},
   "source": [
    "Use a 50-atom argon box and evaluate the diffusion coefficient at 50, 100, 200, and 300 K.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "72eea5119410473aa328ad9291626812",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAADICAYAAACtWK6eAABbYElEQVR4Ae19CXxcVb3/79x7Z0kmkz2ZmSRNuqUre2nZadlBEBBpVBR8ooCoiPD+bD4kCVBABRR9qMBzeQrqa12eyhMQpCkqlaUt3fclaTJbZs3sM/fe8//9bmZi2iZplplJWuZ8mt479557zu/8zvn+fmf5nd8BKIQCBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECB0bFATaqWFmOxDlnGPgokz2CRvx+lJ+OLlpbW9sRedCX+Dy7GQ1Bzhj4MMTXhUe55sCQDSPXmabT1/IeqbFj46GoOW+keSrv8ZTNEe1mpHocT8GzLbTGK4iOKOh4CjOWb0h73H333ZeVlJTESktLo6lUShjue0EQmCzLVXTNxMGCMkVRekVRTGSeTfSKSRYhHRLmkwEj3goBTDc00bQz3yO9DMuulTWRSIAkSVSUZDwep3yyEoqKikg1q5nEotEoKy4u5niNZEsbYvqUfIZPmayO2+tAw8tHCZcvXy6edNJJZ2BDsRoMhtfvvffe8Ej5khTBCjZlKjoTlxpaLBbLGu1Go7Ec09QjGDnSpWVDDQ2BmM2GYFZVtRwT1xow3jOdTidjvrFMuSZ6RXoNCMQaTHOAbsyHIxi9dJ1o+vQ9CatjSWhhPcp79uwR6+vrO9vb2weEx2h5IY024kTjtba2CtjYz0dJvUuv11OFkQg9WqXR+6xJ8RHKEBjh3fHwirRXVsoxnNDCusy60EKwDyW0RlUO0tioqZMoHGaUl5cv6evrexY/HLMwygtACBwopc9PJpMHvvGNb/Q8/vjji7Br4RpVSfMQCRtP1rTRcORmu089XD7YUId6RV2voZ6P59kxI7Sw3dVhucM2m+01FMx6LOzUA0gGHIjk/Q888EAnEkk1NWZVN56aHO03yMTsiNeRMxwyD1T7I3/1IX07XqG1atUqoaWlRXnkkUdmYHdwJgrl1ahNTkSBPC4JMewAORv1Mhgc99xzT+fKlStFTHfIhpKN/AppHD8cIKE11j/S0gSOb3/7200okJs++OCDDhp34NiM2t24Qs4Acjg46DcO0qeU5hgXxwofTUkOUPsiMBA4cBJhBmqMt1Aga+0NNci4hXJOADIUOIj4KcnZAlHHPAeGAge1t+HGfWMpcNYBUgDHWNhfiDtRDgwHDkoXATJuzZGhK6uzWAVwZNh65HXlchCXL8fnW3EMthAnKvDK2qfWZMWRVE/tJyOBI1uUZw0gBXAMXyUoxhhbBQqsOjQObwWhAJJDeTLaX/kAB9GSFYAUwDF8tbamQXDvlfXXL5vNP3mCBUqNet7d61GfY+2udwogGZ53w73JFzgo/wkDpACO4aoR57PT4Lh6ie3JUxrYv58xUw+V00wAJgalzr4bd7dbb2CtzlV8JYisBTVMIRyVA/kEBxEzkUE6o3UNNBtZitNqnZl1DhoY0SLPUH+Yn/Yc56UZrqwfEYfeE1E5DLlOfzDpWveposJyQqWJ3VVTAkpRaakMs2Yqyrz5KUOzRaouF5/56WfLymE5kKFUtmljqL0kAilamQir6Z5nPY/B5c35fb7BQQUaL0BwDYdxNAJbggZ329F8ZD8lRos0R1ncUek92mQlcJ46dnhcTGLCsw5ExwghV+lT4x74I+GwfPkCbJCcGYulFiynEE0xnuJ6CUSTKBjMOhBMiqlIZ5tebTyFit3RulQcSqjQM2oYQ/2NFB/p4bhIL9MYB61M1AvonqFSw7RgGAE2VHrpcuElZ2FUgoHoPnydIx9LB6Mi7jDWaN/cd999V6ERWAS1x260gNWPxvIV45LmUJ1O5wkIEjt+H6RnpFEwDwKPF4GTtUaMeRzSbSGrYMwnilrukOeHlW9UP5FWbGcoko9iplJaXvnA2bOLHr16AahnzzFLC+ZZQFeOZkG+XqV7v1d4c4f/gs/+NLxmVJmOIdJFi+uqrpnH719cD6dNKwVIKOrfVu+Xn/rCT7z5MP4cA6Wjikr1xgcvAo4GHAR4qh+0/TsV29V+/CYwqtwGRaKMxxwQHGVojXtOVVXVFrR30Y0GHJRJBiB2u/0U3LuQRIAEMvtBaLUT30cnsuo5uCBklo10WZFOEe810BG3MPjxWWpw3Anel2GaBqyMDLAF/B3AvMO450X+w5/+fHr3vq0/vmCWKp/ZJAlzLHpWaxaUVDIhru9RkzuKL726xNLcS0JiuLJj0nrKA+keyAPjxvC5ZzDtxN/qkmrlL2+8XDZTXf+rTy/WzT91hhHAUoTREuDbGXx3pfvCz3kqT4lJGPfwb7G7PNCjIHrw/aQLLRRmjMDwxBNPNGJ5Z9IK+WjAQWXD+Mi2iQFkXIN0Gj9g5XfhuKNrMJNHe//0008XYwXvufPOO7O26Wm0eVO8NOM0DTCW74aKixVIrU9CbZVpvED88fl8AmpKdtklF3Y+9Z3NZ7x9QLw1klDhgC8OFSaQIljyngD7cu01l78d712LKqVURWExVBaAaRej6bYZX6qZ/SrYgKn2D4kvijITigU14N7V3jhXmF9WJCXAZpXUuXUgxENyJd+zZElow22/g7OerhGc+phqGtCkyBM9ArBmEAiBAIugK8XyDJTtkAzH+GMooUW8wuBHKT9k4bELz7ELb0JzdSPev4bgyAotoyV9XAChxLHC9NTQ8Jb+Rku0FnfFihUGbAwm7FcewhQa4FPambCqpUX4xG9+o2A+mUfAVy+V2jqWkRnBvx4OvD305vAGlH5LDUv7dpj3hyZy9F/Ro0eB23CP2MFAjH1qkwMsZiPsjCTZMxt39ayE11tG8TkERhMpHUd/6hzLooQsqRFZ1AEYBKYzabMAgmBUa0RH84p7vtI5THrbh3me88fDCS2qI2wnFSgjpqE1uJqJl3OC0hmMGyAkXaihUePNNLijEZ2Ji9KCI0DUw1Ul/h5IohUnEPAXSTjdj26ynnVKvSj2+tgudsGaHoA1MDjuwEeTd6MB/7DsM8/oCh6381G8PNYNjWUAXX56hrNMOPAccTV9IA3i3UiBuiIkNMrKFpaE4v6S7gBn+zwpsOz3QS1+yHgCFFcM9vnFMkyLwaoWAZb3G/Nl0qU0MveZ62gEUSbuaK7DCKURhRZ2x0lIDnT/RpNPtuKMGyDZImCodDINp+Us6+IL5wo/XTpTXDizVoBkoxpyPG572vaA42GcvgQ2cuMaKumJPGNoLoKNCq1F0EzksIY9VOvNPMtcyeQaAd/lJ5nS0qKBY6CLMwxhmW+P6E4NER9paqc8/MYi6zubncJlFnMSO3Je/YzePjDhoKMvrhg2O+DNZSjYVq7EWce0Jh2UVia/gUdTQRDR9mrsig3QlM+bqQgQoa0d+H/V1zdXm9XXTqoTK2Y2mBV9cw3odYq5xOdu3f8oT7EHnSvytcDWitrsYez/twwyF8FGjoPxEaX/4fVIYNAkNE2A4f3RwHH496P5rTVwnJZ4YKdLXQIgVrhCSagrS0GxnhkO+NTtm5zxJ/uV/phoH03ex2WcKQeQpUux4a0BuVlSvlxaJFUIICTF0go9n9YEvIgpgpwCkzF8xz1XV30PlnvD6co+QvJlsbZQUWmNiX330w2nn17HS/x9sV2M+XowUxomjyXvscQdTxFoS4Hgsds3VNfVXbTNqbRhZ+6cciOP6kX22sYDyYeiUX8QiSag5pqW8dA/5b6ZlH7dSFxYtqZfsiUU4bRAlKu9ERDDwRSw3ijWfFSAPkVMyqx4Rl1Z0SrspoyU1kTfUVePAHDZIttpz3y6/h/nz4R3T50uvXnBqabt9m9av0PvsKtHNFCDmyphACQup+OavQH97HU7Yd7abfZbolGPA4ksgGMMNTXlNEg7SkCkX40mlc17vdK5u3tlpWyvX5yfksFcDDwa7GObHNzzpR/tc1M5eQuk+/ZjKPXoojLq6j1pstTWlrBXT6oXauZMM6nFzdUARWA2+Xu/dmAFtrb/cN6Vr67e6MjWYmkgwTsGwc4AIeJ6NLdftUoTPgXNMQZG5lQCj4GOwVG1CqR5lt29nL19gOve2ReV/7qxF/74TzdbvSPKzQYwdz1m++2eb9ReRGbkJOkHJ5CNe+zqidiweHUFuxFtqWoEYEmduUzArp6gzJitQnmVYjZJt37/xroqMjSkrl428s1iGgQSbdyDDGUIDrovgGOMDJ5yGiRdqYLX5epgtbZb3+lk/4ndLH0VWo2c1sDZ1QslOHG2qRqqpOtUf+Q65xO2W633O14Yx6B5VKxSOSvvi4PqjXIh1CdDlTcGoqxgVw89v3FmKNcrJ2JCHenuXi4G3qOic4RIBVCMwJyjvcq65D1ahqN839+PdjteUOJw0rqD/MFwnMfnVAu8rroEYEYTV+cvlIXmBl5WpXtu/T3VzTSjlB4PjDKLkaOtWdP/Ppbi8n4vCHt6Fb55fxB61x+A2Lt7eHivl29xKPIOP9tPMbeuOv6kM+2CXE1WwGiOPwU15MgVmKW3U1GDZIpGINHZvfadporaP6KdQXtcZjyJDhlBZ0JvuqUSGMOy0WyQQJ+8BuM+iX/a+CWTwASvlD8ojP92r0dtfbuT6VSIpnr7EqyiCLjCQffOfvWvK37v7KIu3mHrIhPMetI/xxVg7JexQ6eiSQDhlJ7Gl0mnME8ETGWAEAs0TYJTvQ00m2UPqup+dxzK97ihKIwWHmE/xPyyejDIzRS5g/7LXtDy9rtcW0RL3RfePwhPe6OsrKlCBZMes07wf+xwsk+nG9Lx1I2hsRQuTAJ851N1Z82oVi+Zboat63bDW6zd2XscCoMRW8yxABBQRPkDNOzzbbSz8rKimJyUHWJdl8glNcn3eWThAxd/lUrZuy3r3RwNJB6X/Se1tbUdW53iR/d6WDlmtR4tkl/BK+2x0BoU3h8PQdMcp8+1VV+5kP18SZNwxeLpBtCVCNDUkLCfu8Dysbntrnc/TJpkqgOEJLMY9XgcMYvt2+93C48nVS70BOOAQxGGS9nCAY/y3NOvuP9Jkq2l/dAuQZZaLIFEdLvd+/D6TCZNkrDpfvlxoz2w4ZOpFL/0VP7iPIt02SyLXtbNsALYzLw85KnTfeD80z++aDkR2ly9vK0/boYfY7lS720Vmu0sX6lt5prS/JvqACG+08yQ0OtyPMEtdQfWHuBoRgH1ZiOL4KsX1m51PkoMx75xLhmt0UB04PQv1K4Bvoofd2sK2vhCX2KdL0nCxVxVZZCKRF5Vy3h9DQhuKWWyBGuNB+XPIIieRjBR25Hxb0xhQPvQtDMJGZwAoC3HBMwxJZSnyMcCQIgVJMUZdnV+jdeVds0iVo0D2HHT0PglGSU8hkA0qJnZrTF8d6xEpa4iSDpoCieYgGM+7grKrNYZAkGHzaQPz2NBy4aIAjUTKBAJMvXGi+qqzm7g82SV7WAtdu8E0sv5p8cKQIgRJGEGLGKpNml1GCUPSfdCmDgHNAmOtjWd9iDnO3tBrTZHmSj0MFuPBwQ5xt2+pLDVzfdQVm1jyy8zK8YfvLbuO2iyc9N8i1RZVgy+T5xU99JrH8Qf/MyLvhBpFFIqY0s6t7GPJYAQJwgMmqQjLqZXh+l5IUycA9pYK+xz7ghK1t+8e1BqEfEsqWCkT7CYGZdE0G1xqO7X9wt/wKzIDEdpH2WeOD7UtoFce4bth2idfdspjXqomY3TgeX6SnM8cMelzG9B8/8byAQfppj7o2MNIFQlU0rCjLKNHCvRCCRg0Eu3H/Qp5lCcXYGaBCpwN7Gq8K3dQeXmf25zebDBj2U9RMQ1IqW0ynq6US/cpsddvLqiEpFPaxRUaxkX7V1KbXW05b7ptm+yFsd67DKPdRtBTnl7LAIkpwz5kCdOwod1d3f78PqReJV1cU+fWKJD9xcBj3MtPqMt0gwbvAYkvB91YDhTjFYJaiDOWCiG9jshBUQzOobQzrURVE+KX4iJre9oy+pi76jpGy5iASDDcWYKPifJvQxn0vCPGqhmNUCD3iyTqoGE0uzzOt/LpI1jPXjoIRi3xYCscm93gAn7PaqyvTsMBsNBqOhxo9l2HzicsnAwwA9SXjlYy8oUYVzXvAKko20Z2fbAu0pcNCYS2lhiXFR/CD8icJDkxn5/BhDalY+tuzNazmW6sTQpogXs+qAPgYG8M49HcyU6mV41bbAHo453D4o2gy6e7It5RCt6m8Qes36TXent2M/exEzRhntceYyGjnHFyRdA0rMYa7R58zueOd9fWe5NTcF9FONi4jg+Yq1oTt+2rP/LVWgBoG3nHSYhlCkaOL5wgWXhhc3Cv82rgVOrimGDL6j+DB1gb80RSIiabMwQEtgEv39fsNJiuXaTQ/1zKClWze5VoLKYidGE2nfAx297db2zt6V/z0o28hyGk2N/nC+AaLY930LbnpNq1Uv8nfcb9/QYf8j+B8I5rNyxcyM/X5Dm5O24rRj/BsJwfGhNa46lJ9b92wl18Px5s3W6hgZ0xVUqXFTTG75zS6vlVgTJz6ba4HagYP03pEUEn8v1bnV13dk7nHLLfo/QVKznnZ1eWJkMOXalyzmlwEGk5xogmgOqOput6qbT2c8XN7Arls0thm1uB1zUqN702RNtd7AHHb8frnH08/a4+l8DB/H96RvqP7Wkkc+sN3Hm7lP+xtrdfx1i0VPb1fhUia262sy/O7NK1JWWmVIwE2eAakxqUVW3rinS893Vd9teBnB4h/h+0piHhwWJX1oAbBlRsBCNH/unbwWPx74LnzxKjzMhIwQyv6fSNacAQSYJtONvWS08Nc8iXjGrRpJheh0kSjlUVsTr1YOel9beaV0Ebc4dqIfHMnU4lXg4WloywqL6htPY/541XThnyQw9CNU6qA8noGuF7WnGHP9+mLAQEFFKhZFdrBNYWUIGRQaDDoylqEHKRQiGlJLS3jKRRS7GQfSvNbON7HSLRlumIeNltMGqQW/T4NU0CT4W0HsSdy8FtmbNuMc2g1LP3W0uAcJ+o23ztFVLAlyrcpxJF40ir6hmqoh1a46lKhOxovLS+CewctuwYdCAkBh4XIa0sFDProKnZtcI59jKpKRQbxVgVg3XRfzCNEP33bu/YXkZu0urB43NSONg4KXuCFPtfRy60Plb6V43SCE0RfMHIIQ7G51RjogB6Njav4hK95MVaKBNM2t3fqT+jKUzlIvn1QiSnJDfZsz9uvauv45VDTyDupiTRe/R8s0lQMhjtGKsgLmhBCtGBz2MbHsaHSGsXNzzFE8yCCtk2zMwU3I0Yo/h92lh0VQOPHm1rHA01NfroLiMqeXo91CvlwWTV2VS9GYs4+pB5dQEhsLkNzt9OmGLA2SzIaImUt3MUoqeLdUE32pPiVs9+jfpm47xzTINyu7QW9IEC7cB27rgCEd5h0ZM/2qlXgAC4IZzbE+d3QR3XzDHCOUW3DzDZfB/V3qu7Wv2L6Eg1AA0ZAJT8GEuAaJVrpSCnc4QejF3C4aakrhq2mQXfGEFyiu56peT0gYHCxBfOui/4zigsMCQEvsSwFwhLnT7U7y+OwDFZAiYijDwJpkjwmZo0eigz/5APBT6env3GWpsP1rbKX4xIavg6ItBpQm9iYIg7HarP/rOnx370t0ajefpbydyydhODaR3tEkAGnO0Y4+hucl63Ywa8e7maqaUW6tUmGvFHTMJKK/ovu0LK3gX+w/HY4M05ERozMu3uQQIVbIYDjs8RSbr8+/1CF9jTOGRVEgNxZLcXqHTO/u4+x+d/BcYjy3LwV4OajTLsIFh+oBXlVQ/3Y82DAw02/C7NkynDXfgjs8sm3jBoEkX8UYSrm0uVmo1JxQRbcinu4JgFFJyIBI3bLSr2sJcRz/NGVo1sPT2Om7nzNrTERdu2eaCYvTsElVV5YU1mx2Ppvv3oy3W0eKhEujfUfjgtdb5Z84QrX3B+E48usWOhGDxB8A7ZDp6id0sCtSf1nFuLNFBlQVUMaUyf1CV9H03t3626WlY3okuUY+NkEuAEAeoknEvh/MuACusjrMv7fdzPaAUNBUpWyIJ6fNvbrS7qSGPtfEejb0ZidreT4MWfSySK/29ovWVMREMKuCVj59WATo743KN9esf2NlvsC8i+aJRdaczxorQLejeXsX/Xo/4A8yHdQyiGX9rAKHn/Q6wG57qAsGIUdDcvzvb5v4aAC7Ew3cunA6/OH8GXHJCnSiJuqKQvdn6Y/YA+gDrLz/RlKELbwEWpJ1WRBKs3BMG5gghJlwRqDzgAVGHJiW9CcDntZ7uFDoUADol6ahg0xKe5P9yDZAMIxmBRLXZnnunC2wnzVswY11fw++Dm3/p1xg1Rsl+NJ4R59sxzQevrpt7SoN61TQz1kZceBn3HuwcTQNPg0P96uW1Fy2bJV3dWKYusJrh1X12+Dnty05L7EMayNFowvc0xy/09jp/ixu/Pv/PTt62z8umoSGgYpT4m10+uH/z3p7drRiHaB8iPcpPYqw7huCKybjd6PnbFukYWzfmTUtDpE2P6ExDtLNqMlQbkn86rUE668QmIy+bjd0ks2g293m/5voWOjq613H3UIIGZQdyHffqJ9XdO9ziOdMrU4pB7xPnhmJgMnAlEYuKW5ywe4uMxvTjFzKURV5DrgGSKQxVruBzOnZg49px3z3fsSdd+9St9l9mX3NQA+PAPn+B9clFjexLOFDUl1VJEI0lH+9aYX0UPSE+PFIDz4Bj+dm2W05tEGhhDqrrcGHOBBeXV4Xu3tFuvZYx53vjrGRNo9Ie94aGhl9tc6unxJNKCCKuLcQoDRyHao8M/7Qr0kbmHuilDluj1tVblzVrBJx2JWAqtbXJy0qNwlmiACldSZnEG6cLapVRFTsF1WSKfPWVL9Y9R4Imw6dBBFIdo794ZcWuXmG5UcdMoXhc6fImoERP2hIAx0uPr8EF0lW1x87ETL4AQnxUr79+Obq/XAWbNu0owfMChSeeGFJSDuL52G4zA8VTmq2fnV2r+1pzjcDL6qrwZNlaKE6FpeKDjvZtD9XuxinHXw0lBTE3TXpXWCwn4Jbe5xrLBbWkolSB5mkClOoUU3lXXXnI9WLrcsB92WjZio11HIFAIqLFbAyvZCGr7W1vw6EsJjeU5qAo/eBBzfLFy+rmnj+dn9+AWjGVYm+NVitqiYzwH61JwBqMILAz8PQr1R/lLBxRWHEwAdjJEtDmAUT0AiAa5OkYayfNbh2WnAb+3t7ePVUWy0c/6BG+bw8KC2tL0FOrjveEY8pDf93k+g3lMJJZzWFpTvpPIZ8U4LHRxEQFj9Kio8RI4mQ1LMDpSIBFOlEU7pAEVWECurIyV0jcWi+pddMUKC1RTUbxNi3T5Yf2oekZeZanq0knflyHexWTClNVwYjTseWiWlqNctCslBaLs8+wWU8mCU6O1Sj+OAJ1t6iB0fe0/4GPCA7sktD7L11su2XpDNh88Rz98+edWPT8koXS5q5HrA8hsvDkpSMa7JjIQsmuhZSi+jr9IOz1qnxLZx/0ftAJiff2QHSXh293KGybh/tGSFgDCXrFXG2z1Z26x6cse3ufsuwP76TmITh+QloHv6U4x0zIpwbJNVPS+xR6y0IJPr03DIK9T4Hp6BW++KAPdT9aD+M6TCyZBkbbEOSkG0kixXWuEC3MqdDtjELzHhcIZTif7w6xAKbrjrOKIb4e6yMSEASUEYOmFbHrs3BW3SdqStnzCywir2molGF2DZTIUamkx9m+s9WCxzG4tJX0tEnHiGkO85IaLtMx9XedPqH9n52sWIBYyh9OsmoTCgo8w3DdQXn9Pb+u3Nja6hzJg4wGknXr1qUwPY2jDG9oe3R7DmYqhylL1h4fTwChBocSqivQF7Pu2NkrnN1YkZSLdvTqZgUiUIRTjdFIVNjgAFyCxu26R3YRqDY1rZbiynt7PIKwxa7KRl0fjyCqqJEocgLWd6V873bLmzBiXkyz+7UiriVK6peNOpELOCUEpZUSt03DrllcFsIRoQTPS8Ei/Q96B9Hop/KNI2gNm9wb1dZab1vfLTyLZyqWzuhVodQAEIxz+95e8WaAbUlMW9O0I+ShgS0TD4lSj9Xt0ccTQKi+SFipssy/h0eNnWPEfdThREQ50BvDBTnQOUMctjvUpyiiNn1LN4cGrZGUFRW95gol/vm3A9KZETx2oQcBVl4EIp5ZAlu74akfvdnrXIYSMR996bZ26n41GSOJZKU3wsHZpwozPVgeexC1YlIEv4xHp7CZFy+aWcrYviABF5kwXqCkQeJ8Ea1u393m5Mv3ob8GPI4y1HmQvwTg7EV2jXZHIdFwVA15KPun3q/jDSBUIWgx6kJvS3Ulazvh8a4Aq7WY0ZGuyN3BqPLAXza4VlOcYSSa1rA6cb0CZ5mu7PIrT3oj7KpNds7QNDsWl9Vvv73F8Wwrfp8PcFBzwcaO/zrjwbjVh7NDML0yqRajemwORcEgKUoyEhZw3/i+N9bt6+PYx6fxCH03gaCBJG11uyKTDkmeh7Dc7cfYGCJD/3ivxxtAiA/9FYxTqSaL5WWvUzpBcmIXwatsAXC5qXEfpZIJJJl92TdDRUXZ/q5SbB8Krv52o3MoXGOhXPIXKEuc2uXPbnXCecU6LsaSEbmLtAje+2IMDxRi38covGXbUbs+o6Va4yFG1hzl0Uc4iCcvJvT8QxWOR4BQBVJFihGXi06hepMeUCvTBoqahTE9GTFoIMEYAvjxTD/wa5HTpzTlu9vQrxVdrv+BWlvp2wfYs91BprPhlgGBqSmcjn30L+udv06v7WSTNuLh8ewob8QGkHl5vAKEykeNhXChDSixxQ8eKGrGePQW/xEYhgr0XEsjE2GYbtlQ32b7GTVWwYPnpVRV1b21Kc5P3uHmZk+U/R1CDm3RDgceFKcQssyB4xkgxKpMIx9gG61d0PiBRrIUhlkw7H/Z/38GH4OfTca9BhIvnpeCmdOfFjLTwJnfhWt2OXC8A+QQbmFLZ7TD8Uzc1H3HRwXj/25Q46ylf1wxVVBwCMFH/tBAgo8J3vQ3WCseGbvwZMIcONp89kQywP0KWiVOJI1sfUtdKm36s/U664PtLXzHBdPVXS/cADuc37Q+SOCgGSC8pPVKtrLNSToEEur6kZFioVuVExb/K9FcAISlTTDwoFpcy6Lz7cjF/SQGNHEg37AC+ob90UKr7pHTGg2NtgWW6rJT6hstJ5c/4noS94OTucbKrM0CTWJpC1lnkwO5AEjGx5N4+5WNFWT6QH9pCZ1N2kebloAmDmpxme3kYoNwG5qWy3qTWVUbm7gyZ56KM0Nyeanuq+vur5k9yXSOtjyFeHnkQDbHINrpRNiV0T1yfd19ZzbCDTMqUrWPXVz396BP/j65tVm5rYU0STanIkfDKhICKm5KuhC9gqi45RUicRBKY9iximOXKiFwPZ4Y448J52K8PR39s16FrstoOPshiJM1gKCGoM344g3n1f33aQ3ip86YqQdzYymAGa4pdQSu2fuwdfmsh1aRuXO+Aw0xAO2rDnYHJOGAT+U7e8JQVNQNZXYDyKEAs7tkhoarHorXm/1zDinZQjhGOZAVgOBmGxHBoUyz2T5aaWKfwg0yKX15hajOns6gTJ8SDHv0le6eFSuXL/hjyyrN2C2f7CJtwPD06De7A2ro3S4B+1dx1RF0CzVogMhAEfZ5FM87Xfw9RFJeDBDzWfhCXhPjQFbGINpmG6RDEGBRErsxoSRjKVkQBAV3iKo63IOuVyVJnM2qA6cTuQI6r54Y2WP6mjQICzudvXGZfwkdI3g9aPRnFJK8WEhwe0BRXlwHn//xm24XHSxJEwtjSr0Q+bjmQFY0CNrpaCEqg9AdRNM6nwK7u/pgPrq4N5glUHx+3uVVBXuoPx5tFMhzwB4gCG9b+e9Pmsbuuf5kofLU6UaOm6/Z2WoCWk6KXvy4F15GTVgAR54rZrTZKYrSPzu6qoU+yds4NisAQYK1hqWo8NbuXviP9WauSGJYCeLeTfRCrsiKon/ngLrjufXl7+PBmyCKRiWVQqcceQq42iygzyb1TEX40cxq8aQas14xTKvDLqAFhFhANeg777jJa/kAvRr+ZBQr63mievKzoen6mox/XZzooKnwSaEq5JD0jQtkzYIa7X203ZNt2pg35/RkCyCEaMHX6/iLKNie+0eneBtusIG9niSUFTHRH1FD+/zKl7Zt2zYJygMGvBqqavKqeIrzhIrOUHXFwEzl6LNJ4ILkpP25n8Ey/GSCm44mpf3kIFNtUfVwk36aqp8MkDz8zH977r3zy/Oeu8lygbVE72bs4FYsM08baOZU62cLIFRHRCidZ/5FgLp33o7BJ7Y4uVish4OOIPynw+Fav3y55rQhb+qRiKKgcbBUEQIxEB1BzjrxaJK6Th+UYzNgahQUd4K8Guq0yG3a/x/m/3AYRkIaXbC31H307EY4zVbCAyE/e5W1j85tUraY14rtqR0115cvrv7qkuRvvrFovq66ogRk3/fqV+84qN6Pzr7X5xq02QYI8RV9YNl/itef2gdxqhWlT1vbSjwwnvif10A0CdA3qy9QtGv9Zidbai1NpJBK3Qz0algkyUowmhI3OXkHUTXUVlx6/iEJZB4EVZWVpTedYVx53gzhsrNwul5XIUIsnEg6n7B9hd3veAHRk/ODNjUjTLSbO+cE21fm1LJnTq1Hx5TNNSpYS0ScBrpkPnhO+cPna9G7jNvNqbuluUHKfi1lZRbrMLKIxzRLJRIoqKB4r61mHxYvnz8RlWtk9If50DYXhzd2M33HrgR7dWuYvbotof/FeynXz99j30WCGHYrct6vzWfBx5IXjdWQUXxhnf5hPK7isumVuqSuwSLDgrly0cJpekuD8flt37Auxcao4pHNuZyJZCu1erAVGyT4hknPVb2hSOFWm6A0NzOwNSXKp5tqFjbCpzRgtOWOlmxqkMF1oXWj2lE/TpGgjZE8HudbNTW2yz/o5q0HfDC3sghdEAn8b7u97AG73e5FWkm9EcCneqAZHSHLA+j0WK2mBN1ffZKOq5CZXgJzpaDW1qHLRKMseAKiIMRuQeaswRm/XAbqZ6h6M5uGY0Yzdo2ZL8qF+pAMYh8OY8OyhKfjck9MuBCJ+C60oVDLUVvLFUByybzxpq0tGKIT6NcwgdecZY0VeAANPvMH0wlOdXBooFiOvr+Eh0HNxQC6XzLoDMEYN7rDIPYEFF5v7wODEd2aJEO4sVKmRjld49e/PNCn2ZfVi0aKTlX87rCQ3Oflhu2OOJRsczCbrw+YEuMRXxz2+5jWi+9owy40AiqrFKQT+zABhIqc6f6pEOzyEyIewm5gWtPlhMFpPk/0ooH3X6BYoH/6M4HLzmyA0yr0cqC7F/6Etm77JjhgJd5gPrqIJ5L07XDxEqs5js4unGKTOwB6dKPXF4rrt7qH9EA/0fId/j3RIkQiLnewxNqxvke8pkSfSsaTfqHeGeQ6xsVOv8rWdam/pg9/kEPzoA8bQIifWvcPr9qAFMGRd2DQ2GwZNgD8U/H8Pmr8OPgZoIt+Dg7U3cBx8SLdY8sd9587nS+rMAbq8FiReU0WIxhrddDoi7ajrdtX2EPOFyc4gEZJ3BlPpqxPbrCLz4pMhUA8oticMXQfCno849z/gVP8AZ7QK6HNmkBrRnSq1TI0MyIaBxOdhXsNsCKDW/b1ylXRpHDufh+HWrNChwNAX1S99Zd/d3UQL7EOM3WahWwPTeLDCJAMB7JdoZl0R7y2pics2g/rEmDDHmomRgMxg9mGG5fa/7KkSTp/8SwDGBtM+HVKjdlDaqK2iRc3JcvwwMfn1t5rXYe+q3ZMQJNQQyMP9D8EwVL7t/3C3Xt9grm6GP3yCnyrO6x+bu1W5+6BAq4auNMW77QB878eTfROA4gTTYRaW1uvePnPr97z1t5OtVgvKIEEe93ncryT4eVEMxrp+w8zQEbiS07eaRWKTfvaxXVzrjmJ33KyFU6tLWHBaFz9FXqM/83hDRt9BYsMvaHPmRFe3lAunV9mZEldVaWoNjcxpudCEd8q+A54QVoyK2Wo8RSXdCaycd4jaVScqne1VVY2vLA5xucwtPMIB5zkaBtHyOhI+7q6i09v4Nc2lYG5RM/fWP4N568RHCls0Zq2ozhZChpIMC3pyUfaf3HB5ZfvyaSb1hw51/4FgGQ4nusrdgUexu7cWQutixc1sVfOmaGram4qwkUHXJ+MRq87+Jjth+zrDjrDb2C1ura2v9uiqOwTSTzXMJISREXViXqhCFSsOVZcCvI+B6h7PUyMcDWUytp0JzVM0efr7sEr/WU0BNx4vu0Hp9QJty9t1oHZioN3Sb3pze+Lt+3er14LT+FR1ASS7K5JEC3QsXZtOQoMae7cRez559fRMRA561ZRfplAo/9CyAMHeBu1saVSsY79oqmSVZWa9ElobFDUBQvweIZGuWF60e3v319zKUMQZbYok9kRkRZNgrknAHDAJ0NnVxCULQdB2HoQUj0hsPs5uLY51F0H4sK73f3nPbbRRxMPWncLkxHRP7CeVMNMPH+wrly4vaFMkItr0IH23GYF5s9Oli6sOMdaDU9RHPTpmpM2ZTQalTWoTZ977n0Zc8kLOIiFBQ1CXMh90Faedabt8/HktZnxJMcdKbgNQEJ7sKJS3NmIex2L9GqJQTgfSfkL/Os4Z236MqXwP2x3w9L6MiUhCUHmDcZYBSqRvmgCtjpVMLrj+t29yv7NPdJP02OZbDUgrQuzcCG6rMaZIty8cDOCQE1xJEtXLPHSKlys4AJzelVR57+69ZqmMtbSGRhmPJV7LucghwJAcsDUIZIk2QqiHkq8UTy/r4/Dfjyzr3pvL5gS6Cw91geyNwU9fUcIX02Km4oMP3QGE8tX75XO8uK52Xt642A24iE8ATy1qRdH1QJf2+ljX9q4u9vXhs4pMKus9s0z5w9GU6zMHUK3rEGVNTnCULfbBaIRc3NFBV+EMw9L9a+ua6XF58dBKAAkP5VIDZZVGMWNOBMU2OgQqiqK47KiusRpPV6c10+pzqAsvdfNdxM5aA+mda3wVruSM+2ZMyuu6PSyu/AgzE+WGnkJdtXeweOkX+zsE3ugr/u9dFzsoWUXHEQPpqmBDg8nfm+7m53bVC7LeskvzA3FWZlRlePxuLSuh7u2bNBHjiftQWUvAIS4kPtADV1wOBzRGovl3g097CcKB8kRioOlhDFJZPpdbmXt914t+XUrxjvMHoy+Zfv2aSv+bV40sgWoQbmNRwSlA/X9H+Ig5AIc6SxIkzEwSk8e8MqfWb1XrAnGk4BW0WA2MJ0fHWDsdvGvr0Egt7RoEwXZ6uJlijhp1wJA8sd60iK4HcD1UzzlNvJuJ//anl6YU1oEEZz0+d3bB6Ptq1sdyjJcOGwjQyca7LbgJqV+LaKBBJ9SF0bGM5XD5Egb7wEH8jT6V3MIDspGA7j34EG7zWa7fLtLfQK7iefWmEBvMqjbAlF1xTvbXb8lcOPGtOMGHFTwAkCIC/kLGkjwlNuVmOVKB6A9WPq8cyLhAmzlQwTEiNZAqZHSDI62eDgJjrS1biJqwfVIw6WBMsuMfW59cSpycDv+VvO1LoF55TUUAJJXdmuZUUMj6Y/XLr+AzV/Fpr/87Lq5l82D/4cLcDPqsNsVTamrH/ud85vPrztiAY6AMlmB8qbxCE8EXfuJCOreXX89HJPnDxL9RwsFgByNQ7l5r3VDaK8M7Xs4obn+pLk18NezpotVC2YUA9TiEEOJX7jCKFyysMZ+KbwCSWyZ2V6AG2/JCOAUtCk37N7xSdBm/RTk4f8CQPLA5OGyoAM6GSwVzzHs/H5DpVRlNhoSUGeT1Nk2EKI+pTq17/xrIrXL8Vz3F3GFneqKulhTJWSAkhN6qMu2LA3C9dAnJIHmk/MfCgDJP88zOTIyPQHYW66q7ORokvNwStADN+C0VjE6k0jg+oYeFxTFT+AHL6LV76R1rTKNdVmG8jbtfJWc0ZMZz+CQTANhU+tXo3cKP9MmJTIk5Ot6TACEuhfUPLC/m7NKyRfDB+ejFaYUwBdFg6cAeYFJQM1eD1TjybpMjkLKhZuC/Oi8lUJuNygNJuuQ+1aamUIgZxqr9hJ/oEYbsBk75IMJ/siA48Gr6+YuaoSrGkpk7k99+q2fJh4ITjDpcX0+1QGi7aKjQ28QIprBHPoNy9sJs+Pi6Og/InxoziT8RTs3brQL51WbkklV9eiaXAFuFGTZF1EMH9j5HyjJtg6tz5/Tbs3hpGvgQCm+/DzLjGsXSFcvruelEvC+7v38l6zd2ZvtRcEMGG9aanvqlAa4Y1mzXldeVYS2aA612f21//oVwO3Qhg2hf1bvcHJz8nsqA4Q4oR2lsHxRRdkFS8xoz9Tlx2dKRsrkhCN5T3SNzNWaO7a72GoVxAp7XxJspSlmkJhhv0dZ+587DT/UGs6a/K4vZHh85SLLkrOmSX8+e4ZUNb1BD2ASoNYaubdruvVaNNF/L1uaJOPF5MRm26dm1Yh3N9cIvNyGBpHNFihWIsL0Hvut+x+x/JN9A9eRcKPWCBvMslqDUxUgGjguPbOh8oqZypOLG9hVM6sU8fqn63d6/Wr7/HbHa9mqmKxyc+yJkUbADUq9G6urbRdvdagPdnr54rJiHtEx+PU72+PfAXDE26mLmUepiXmxtnbg7SXWGjxT5U/NNayqsqw4BTOnMag2qiavo05NOF56+66Gk6GtO87bJj7DRhMWlK8k8Nv1IleZqFPBXCHhoir2HNCYsy8EBl34Zozzs3w695uKAGEovdh3v9tYPtOcevX0RmnxoplFYGwsAzAqZ1V2+V7d1V57KWt1v55PSYIVk6uggcTj0RbgrnNBTQmakcQxs8yMVb7BQeWk5Rml2sgv1EmsNkW7XCXcdGuuALWyDDeLx2VzqbdZ1MXOx3Hha2nz/AmtoGuAxPWhSIIVob0ZswcVYbo7AuaDPkQfGkGiWYsrzGZfvGhmKWP7gogmWoIhUOU0aHPZOc1hjInTLjoaFFoqUrdWlYiLjRJLSJVVXG1u5ursuSmxqQrKTNLjrUuXSvmUJGMsxlijayDBj3CmpjeMFS/T5iD8PRnggKX9+dJU1TRvBFQ0ioSe3hioBzwg7HGiUutjkbCq9kS4lQpK+9LpOpGACVBblANx3DOMFso7XSllw04P2N/ZC/51+xRHZ5htdcO+N9bt69N6D3kAB5VnymmQZctAO7weN+ZfHklwNZgQJDnBmDHBQBWYBEkJcFp0dmntXpQk4Mv2QHEilTzBbzMDcM2UZA1uDppgeuP+fE16elXl7E10lCBsdYJSYozwlNzNaktFzuU43+mUxY12tpkyWZaOP+4MB32IGyef3uZkVxlFpgsnokqnJw6480RyhzlscbFnMCpvQYcRgz7J6e2UAwhqD0Ylxl10ERqS7/fKvPaAH5pR6esMHBK+AN/vU9Uut76/QWmxs8sjbYWbDAZpapW8juCVoVbLbi7DppbzbsOwOf/rhabRgh7HRoNo/cvbndKlSUUGd19IxQOScLJEkLbYlZVP/p9jI/EKeTOh7lU6W0pD8Lpcq5ml7vN4oNHj3UFWazHjhjCd6JMl84rfdWxbmW//zlMOIMgkrcmnVOVX293ClbUmVeYQUjzBGJQbuZyQVcPaA/D3Z9Z0Bsk9vzYF/K+KnfAdjn9o3l/BdA8JpNbb8EkbAoYsaNMm6VOhMSNVOQlUNi4K7NPdAeVb/ih8fItTKDXreZ8Kyn+/vj7xIGpvlQYCWQwaMNGY8ycmi+Vlv0s6YYtLgWkzT+7++lfvVE9cermQb//OUxEg1LUQLnS7f93BrEvX7Jdu6Y1y2FmRApOBia6Qat/eK/w7xuFb+2c+slY/aXCoD1xVe9EZM6QvnWxlpUU69aDTpT6HZ4e8Qxm1D8rtOOreDSrVwK0GfrTe9eCTm4sqG77hjvA5Ea+wC0Bz5qCNkgdiZ++GQCJGXC43Xt+kZH/9y0+ZXn5v7XwUXGpbW78Apef5CFMRIFRujkJaBZfzVlW1/i0QZddusrNyg8Tf2rJH/QH6fsJhHDBiWLaYpM3Do+b42JK6/zfPJnz73JkSVDWW4Lw/Oh7pDd+45xHbTb/8gO34+ElQiVqs77RHnO+j9ESM5Gc2JVvlHEc6pCOEWNrDCf2gvSioRYn3udKgmvl8RwcIOBbjv169TV9WVpaNbtyYiz9lAZIuCfP2On+B9/Q3EHLQKAWySDWb65tLi/gTFjPworJSBd14MKgt5olt+6TVbxx88fyZgjDLogMJvd14n6n/W8QXupa19x1PEwUDPB50QyDQxgd4pQmEnFrvUreZXKyi8KN8NQ8vT+wdRE2eb/M2GzDOchGTCMTaUQqZqU+UYvQ8awHT1fhgMvFrUB+I8RRTUqpe4jqTSHMoj/43nrqDUvTkJhMYZjWo4smzlMpTqs8rKS/5v+dutRUDqn0kiITrpAfqJq5Gy18aM9H6RHqNIht0kcYgoGRNax9OFHVZ0/6HhZdub6y468yGIlox9wmz8UhY8m6a/zBVNchgTmjTnShRKOSmctb0Z4fudczOEFN7girstYfhtP0u2Pp3GbZuD7Izz9FxpjcxXNkVeH0ZCAcFuaI6dOappeGzWTu8oTXE4f3r9meQ4/9bEcTU7URWHcKnqTpWIjAvQ5pxmzFpJRG7rMm7Lq/75NLZ8MBJNrX+o/Mh/PWPW377as/qJ3daT0vlmH1DJn8sAGRIwrP5EPGhaSSVKa/u8woPbbRzRSdGVB07yDq2y2oQj7XGBTPBg2dTlOG5ADjRCRBIgJrgiivGFiMtb2RjsWwiZdLAgcD4xLl1cz4yj99yEro1tZrRrWn0SLempO3oyGvKb7Jm4wgch4FZufAk6xebLcIPFzVK0DCzDH1TFFWZk8G7LxTXnvtu77mfI3rb2tqQ+dntQVC6w4VJBcjAKaooQaiBLaPp1TwWfhBTtD623+3+p1SDh5AeEG8LJRTwRGLgDjPRHRFgh0tVK0xhMpiEmk49bvgLQ5dbFvf3ie9ROujtXAPZoDTzd4uN7WHUHOTW9CQbe+VsdGs6m9yaVupwQSl23cEn0K3p/f1uTemIOW1qfJBzBeqOoRY8ROvkkngNzJjfly6zLFw2Q7h5ViU/CWfx33/q7+zLpQauGozoLbthmqBMr+ai1yHXBfYtafH/6Nzvo/+6ZdAhtB+mIXNJ66QAhM68Rgf2Q5qtT2J3QGvgeMDOFzmzdv8txj611cltxRJs9Ee5bW2nOFdgSTkQ8UN1CUN3gky3qUf+4Lf74G2Shi15bGCHNwjeRguZS6US/a5fNKJbU7PJQG5NRbWxigt+NzQkO29/74Ha3+MZIn/Fb5WlC2pKPnemsLihDCIXf8f1PoEjXyDJaI4rF9vQebfw7DmzJF1dczFs2C1f7I+EIZzElpFA5YGuhES0oIAUnkisMrVUdTZQuZfBmsOLn9PfkwIQvV6vdl3+dKj9g5cuOq1J/FijGcxmY069hI+GiRkNwDxu56P4wWPd0FCOc/4+dHUzb4+H/6wvzs5AF6CAjtsgkVK37vEon92wyx296iqtu5L5fjR5ZTMOYhVUcmvKedqtqYJuTcW0W9NYXAajXi3WaceVvf7Vy223njeTPbSkUayvLhUg+Gz9Zo9b/gKu87ybB5Bomo5mC8uL+HPTqxgzl5tT0DxdkHHbvbtvo0huLHY4YlCy3QHVdJpUMsSDriTb6RdJy0NHNjk3irTyCpBM/7FPLU/J/zH7ZyfOlm5YNscAJeQlXFBueuMZ8ba1e5SPQaXTi1KRFhgmo9FRniL+YYV0+/olnmPHUoBzP6iqP98e5Oegn4JNXrfjzxiHBo5ZXY/B9MYaUBsjwWm3pkgfujWNQxW6NS1B527k1jTlTbGYKvjq6+svnFbBnjupToTGWRUK1JchhyMnGne4X9ndWnMGtPXuxcLnrLtFs4W4rqEWmZSr9aLIEgrDQweNOmClcPJZRWCp3wOvb4/yEkNSiSc9rK7bx0VQhAMBxtYF524HsENHB9ZEHrVIXgHS0dFBA0N15QuPL7v+xNQN0ysl8hIOMBfP+JVSSnlp1zkLI97H2dfhC3xhppFS9ec9kLTSGl57f9dJQMUug7fnTXxOfwTeXHszpGxGE2jscIRbU1npd2uqh5TqQLemr+2Qu/EE5a8X42YLSadXoKpWUmdOAyEeSOpDsUroTv4HlulzOBunDd5Hk/GY4yATKaiqUNMbAe5Ao5VuVxTm7HKCvq4Ynr2hFJY/FWV/2gbSQq+CXVkVEkkFfIr50ZOWf209vLQad1Z2KO3E/DyFvAKkds0aks4Q7vNdoiomRea4Q01v0ryEcx06EUQv4ejz9SP3XF1lFlq8oUkcjxCZGq10g0FrhHgVUH4xrGeOtB0xnarFzP9/ROehbk1VdGvaF8fjysitqaDf55Zf+c6r7I+nNMPz5GS6N8zFpiAeJuKLoQ7EvRbYWGMym6ORvvyQcme1NMQ3SlAF+a97PNJ9W+zAi/VhHkt0Q9VugddJCX7PhYJy38vqj7sDbG6JAbr6wurPu5073/rWihUnZ5WYUSaWV4CsSjNI5lIVzpYKB4OKOs0eAmupC5geAYJzprhZxhSNGPQaJ/MnKEbDLiJJwUqeioEA/C+3pl38a+j1fU4ZujXFI2F/t3ZLvK2hwaR6o0pkr4cXNboSzLzDxRr7IriimOCxvijf6QUvFQz3/OeS66SZBZ/b/brErK/9vVO6LCbLYA9GAI9zYAmZ4cm1/Ht+j+NuPxGTDg888JWqsrIGrUlknuXrmleAYKG0LlY0qa7f4YZzZzhlWSf6hTnoJbzUoMoJ9BK+wQ57tsjYk87z1GO+GJ7DfDSQDLg11Y65VnFnYjeqCTyyoNsPVTWW5zbYdQ+V6JVkIhUUOnGThUFQRReu8Wzu1vZa5JC8gaSpoQtFRsMnD/pSDwci8MnNDs5MelQkCnv+H5sdj+J7He42UHFDHPkoBoOhmoA1KSGvAEFbflwxXQUXLD335Y0bN97y5i65KIgbPw6kvYR7I5zWGx7HgZzcUquNQSaFKcdwpgQSmmAYOOY6bVhIjZKXmZNPO4LsjL/uES/rDAJY8cRYMrd0hvjTf3rX8VeMM+TUOz7PZtBowSMdApjoVz2VlQ/uP2hGwang5G53jIYXSFIKcaGBgy6t9F9/wLdtjBKgaJmHubzmFSArV67E/QMMbrzxxi6Xy3XZ5gOdDzpD5O6m30u4P6J5Cf8NFlgzHsxlwY/jtDPSVmtIZISZLmvmCIXLlaq6q/Egn1NRe8iRuPiGz2d/B1sbtU0CWL4CZoc9Ch/O5QLuO8cfaTBn6D2EDkWJiwQKxh7WgKGZ9izX9qPkFCh5BUimxDt27Ch5BwOC5bJQwjJjj+v49xKeKfs4rwylqhawFY+2QRwej35TO+Rer/2PeKU/LdBUNuufrcs8yseV6CEwEE1aoQaBmR4NhEQiIZbOOiuAH7AXPj+/Ihp3yqyFgKV9pymdgchZvhGynN6oktPpdOrdd9+Nmy1wYIZewtGF/lYspZrek5FPKTYqeiczEvEE80fJ2f+HQGE0PhsnTRnQiEvRSjptHa05yRhnetn4jGjK/A2VHnv88cf9B/7v25946RbL+vOnB3ffeHrxnt6nbT95+67SStIqxJOhPszGs0nRIER4aWkpAYEYo1U2FjKn+wwoz2MtkGTHdRiSssIPbmjEPUMoUfqd52mNYgza5PCi98/GrTn88dT6TcKBtMri+baPLT9FenTJdCPMnoVMsBYjoZHPFemFE/52P7sM2oJB0rAT4MewBZ80gAyiqKAxBjEjc9uKoKBFyi9cXNdy8Wx48LR6ud5aJoLz6bp3Orvha4zZd5HkzEWjyNAw2deVK2mM0WQs0SefqCzi3Gg0yGCzSUqzDcSQL2kKJxdbeuVbWHvwW4gPasva1ohs0j1eVZ1NGgppHcYBrauJA+bT51o/PrtK+J/TpulObJ5TXWk+eVql5bSaK+bOEN566w5rDX2Wy+7FYWTl+yftKsB2nzKhYbUFDRhZIIZun+LA0PE9GjHifUJQcQfCRzTC2nIzwTAVNEi+GT/l88tIziJ98jE06uMGPUpOC0rOOXUghn2psnDC0ujzfw61x7dwPJITyTkVmITowKBTeiNJuRM3Nu9Eb/dlO5xgwwVOpsR43Btn3XTKL4VVuRmHFACicXdK/fcvyYmwQIsR5o+B1BhFd7V4RjpKUBElJw+mhMuR6m8BSc72KUV/toghfOAERWcwGLX+4/1u9tEyYyqRlH1SoyPI9Wjca+/jAjqv+y/KEMcqOQkFgOSErRNLVJOcZTrFE07yAz6uouRkZTtcUBeKotFzjMd8MTINwVVybBi5NQ2ZWEEm/jWxgptLiv+9O6QsenWnXNcZUHGBE4etqD4dAeWHv3jLsboVx2vpvewTz/GwFAoAOYwhU+Bnv+QMdgbDBtsf19uFm8qMMkpOvzTNGeRGUQVHiAvvd3NNck4BenNJgjaBs3fv3oP33nvvJ5//8U8XOYL6jxYZoI8r8NKeLsdvCByoQHM20VMASC6rd/xpaxUuSezuAx5l3l9kcUmnHyWnNjOOXsH61B/+8m/O3+VSco6f9Ox/iWtmxTNnNgSCPs93gwDfzeSQnsXLGTgonwJAMtyeWlfSIsxut3tnV1Ze0ukz3IzmmxeY9EIwofJfHDzoeD3XkpP8BeBsWn9owenW/jWrSeESrZkFgzEdZi6hc150Ktg/4MBe1pBmKdkksgCQbHIzu2lpINmj2SppUvO7meRzLDk1sxat8Q0a+NLq/SSYo2SKDKIoEj/k5StXshammZkNvMvlTQEgueTuxNPWQILJaOYmmkTHRptDyYmKQjNpgW+3WBef2sDqq4rUwBdWuP6B4Egd7wuTQ1VXASBDcWVqPdMkJ5F0tKlMjMjaWkFso0EreqEfoyslTXOcPtdWfcV89vMzZwhXnIb+qYpNDN78Qf3Wzm7l3/BMwvcnW5Pku2oKAMk3x3OXHyoWXHluB7l9HHlgw9fMVi49lb+4sE66rNmqV4pnWwBqS3hp1LNwhupa+fLtZYvW1QXD7z+Hlh64XXZZGx1CpdnTjSPHY+OTAkCOjXo6GpVa1+i8edW2TyzW3b24AU6rK2GyzPkr1//e8ey6ddg90npmwzZmbXyhL7HON+iEi9F+RRZ0xSKvtjDeWAOCU5cqsfTNmO6UrznhtuDPBohBJB7vGqUAkIHazumNdt475bAcjw3QJH32siN7OnX2tJpZS2bp/obHNthOnmFEyY9/auLSPxvZx34w034JrEyDZGiJTwADSQdN4SQTPWGu9vbJrAb34jI9OtYIxAQFHTsc7IMl37/BarxwDrPhZqv96/ez/8vFeenZY83EUyoAZOI8HDGF1vRC1uCVXtoNl61zvnHgTideMUup1DazWrCVFUkJsnhV59jwUI+AUqvuO/+msOVGxlw/1nbh9W9SOpxmGufgNBV09gRUZZeH8ZruKBPFHmbt8QBPRcX9jji6GBFuvWi+Tpw/3YAnDqP3w5qYfXNDds9LP5ywyf5dAEhuawAnf7RVXuGpT1lOn1sNpgNO9T3Wop1kO2x/Zwwksd9oW2qbypNq8spYkvOILOmBGRnTl6D/KRUE0aAqgvBpTPPHI5wKTIttQtjn3BHQWd9456B4mSik5GA0gEMQNA1DT6t4hjvMqBLEhhoTnptSz6HKqJZ67HX1SedLP/+M5VRoc0U5ni+ZZe04BlbkJmrB3D03fIXW/l1//MpFliXPfKbugwubde8snWt487MXGrc7HrfdQWMC6r9PNHtN9JcqAjp/FjWvit4keA54gW3pAmGXAxR3HPAwTFpkA2jT/h/uPy0pPJPwM7h6/8prO0D6w1aQfrNJFZ5fqya3uRkU6QTgUpEIZVWSWl2nh+JKpaLU0Dzbxk7XBuu5dDo3HNU5fn68a5D+vdwo1vIs2VhbO/DvV9Y3NFTwP50+Tayd11SiGmdVoaM8paGk1/O9/Y9aFfag8wcT7G5Roxagb1ZfoGjnxk1O4bxacxL1iEc33RUEk6TIgahs+MDOO6gdrcKpX7oOEygtlj6T8COxytoz3SHxCnz4j4TMly2dLdzfG1EVhzchleLxw0IcnZD0RlgojFqLqxoAjzYNPUy+U/rxcQsQ2nRE3Q/S+RTy5QWD8sK93iJD10XNZvXWCpNYKzIhqSur0POmGXi2G8gCV1hJUeT//fSzTT+B5Z20/YcmaCcQ1shcrblju5OtxjXFCmcoCXWlKSjSM8NBv/Leenv8SdJoLWgycpRMNJBQnJDP/U+80B9UWixBdOj2wDYngLkozGW5i1WXSuicJ6rscKR0b3cJQYp3PIbjBiAECLQf0s42x4PmaeDav5f79say3RvUOGvRHKhRHU64OR6tISxbA+oajJRS4Ny+GOde7P5EQwqY/YiFBBMhBCylsAqHmjIhgOMTXKHWxg+9vb0bq6vrLtpsV9oO+OCcCiOPorHja2u3iPcB+INvrBt1uTNQFZcir7AczOdyvSdB/3npKTVFM1y8EreF43npuvXdyl+e+j/nOg2A/fvnj8aeY+r9cQEQqpx2OnxnkO3Q1z5S13LeDPiPk21yww1zhdC919W9vrYT7r/uWbt3gg1yqApmrag12pahhG4D+N4VoOOvQLLWyh37fAB7emVesc8P82QFivBcm1QoqG5z4dHWXbp+YzuC7MSCBhKPx74Bk7nGVdZUDsEkuna3R9PJUg6Zhj/anPodOyAO6AM6L/2AV/k5npd+yRYnw7NIUJGo8sp/7GR3UNrI/+MyHPMASYNDfeCa2otOq5O+vGQalL24gYdRMl+9ZLoeGqaX4XyksbJMDn/hfNG3cOXyhougrZvsirK1Cqw1vnbsUuEfQH9DSdyJTy0W/sJON7+h3AiiqkZSLtxUXWHE1QkV9P88wP/432s6g2Q1q534NPHmpYEEk2EQ7AwQUVk6rpmANTA2iVbZ5qE7ZVskwR0QduxIk03ZjRWA6U+n9uWYBgh1q0hzXINnm8+tEb598XwdlKAX5D1v9sHZ01Su0xdxta6eCbNrOYQ9qapI8qxmX+gq1g6r8HgFXAEDlLITCjTEQawt0j16vR3ByT4ysxw4ng/0l99vSb5w50uujtpa221rO9l/ort/XVOFijNB6Ecwyv+xzgX/TjmnzwicEBGDPiaQUDjcq2L/0/H/r4GEPg97NVBowEgLJ3p3XIKDynssA0RzT1pRYTmhohiesOK0f7m1Ql3jN7GNXV5Y0iCJgRguoCkCdmlUtn9jn07sSqlFOvaVDfdatrEW11aOXTMES6ZRET/GEtI+YhfpbrnI8cZ5M6Xzz5ylB70V+1CSfOnnSkIt88yWj17+I8fztRbL2xt62DU7ewWJc3VfT4/zV5gRuajJleTNRYPNpElT0xrdKJzGy7ux8HlS4x6zAMGZov7TiorEjwsCnm2ugJySdZLFUgYhWYR9Hg7T8aQl9n4PvLJlO8w0hIUzm428Zq75/HgiscH9hO3LeLDlCyj+tSPMxloLqL0E6hqdMNvxwPQq4fxKk5DUW2pEmF+P8jQhm/X7l8zxex/FVnWbp9e1RVFhiy+dSVrt5AocYy3KWOMf96AYzJAJL1QNTiyv92v6c0soqt4ZAtXRh57h9/fBXAjDLJTif92pwGZ7Unnijy48sCcCp8wsgZr5DQxOnqsYFzZKNdOMz2/7hnUpNlY1bYIxJvIXLOjvVuBGjUtkhatxGXf0CLiIZiwT1eJy7EgZVVEnXHTrVbaiX30cxEXovwZBrbn7pEkCzCwjkceUbyFyfjlwzGoQxIfWwLDL0rHXI359E55trpfCqqLI7M6zGb9jFWO/Wg8iHisMpzeKaKmHbdZUDmq1TcQDJ2ShzCcKQuwWZHcaamNkfHt//ECCq/Y+dAgaUHj9wSDYDN1o3oGn+DqjgKAtisgmCccZUSRWxTnTAijGyGacHZSqDR4xJVSLq/t9gOV8m+1gEo9ZgGAhtHUOOq1ItPSfbR5JqnDQH4Fp5QK77iQRfvqu8germX2kNwy6bn8K6nqCeOArAiUZYuBPcTweb7rGjK1jb7iID0bfRhKwb6sTzseFOVlgXmGWLwxlBp6MJhIG7N2te+nVPX1ZnKnSyP2w/IfShKFQkWHZf/qvE1+uvBv3uvC2Nur15E3QHMsAoXZCjMKjx/Bsc27t9seEWzY5uKlIp/QYRPaY1+38H51o24drDo2W0gSXRJfU5AqAniWVYDiu3+KC9yiRjg7NJmqsfWvKG+345BV7PdI1qI4qfJEUrXlAiYEZHEHVv9ejnf3Ctqa7Y5RXIYyKA+kJEOCPXFd/42nBho9YS/VF33zI+lvG2n9BKWjgyQNQjgeAaA21/2zzhqcOgmAE6AoiD6nB44q1+tQHdt1/CkwFfySi2hxRMOiYfp9X8a3vEn9AcTpogW8N3o0tUPoCrmDvqa6ruwhPxmqzB+DMimI059Oxdais7tvT6dhAcT4Msz1jY93IsXH6GHGAx09dbn3u7FnCLWfP1IGxBmflmXiN68m6T/7gHfvHjrK/ZeQMxvD2WAdIpqgEEhwva+Yk2jFe11+vuc7nvS7XD3itteLtA+y2fV6hARswrgrzza6A+oV1ux2703P5Y9UemXw1EHrs6RVsqChDayXsAqA1H4YJpp3J40N1zaxtndhsfdBaJt1iNQtJoxVnB+fZUOTFlNrizo/cHJfvYcy9YoT9LVnj2fECEGIIjUkY3eAsETk4GBjMkXapqan57taYfnEipUTkPtf7GE1tzY50x9y0Lhpe/EEc3OC0GDpPYMAKmgNZMsaQnh0UBGBXqypXk6okgr5E5OYq7FYlQND14iaX4CdvXbToW7B8Ha0l5TRMKYCQxF2IhobLqc+e9spxAQ7MxsABaqxDBRG7QmF8sZpe0jrEQ7j+gVp8vJrj8Dwy6fT3nfu3tQ5Hy+HfFn4P4gAKFeTbbH0oESl29XF20C/DNOyvVhU7sN5SAqCHd1cIKpyGMO4KgxBGpurMGa8nDSBFRUWHF2pYiTuICbi/4/DPBnF3+Fu1DWc/Ojo6tHWfZcuW0W/ehszFq6Z1Bn9K7wb/Hst9SwsTFrhBS3NbLVrZr+QZ8BySDLkQwUD50M2486NEhgm5SneY7Mb/mGb5ahYAW4ZJ3PYysOfX7UmEEtYPtrqFhfXlqZQkenWzcXayRKekovG4fpMT1v/x7Z2hfMwOThZA1FgshmMGahwcli9vQbORVcq/nWc7bekc9uCZ01Chmpgqy+zVugd6vqNFwsgY8Kw+qvdxBWqEWmNds2YNtKOoSocjGuegd5k4Y7kOdO3oo1HQe0T+Y8lshLi5SveIChin0AISTshrPni/PpWHFlR37dI/vcstX6kThfJALAn7epNgMjB9l5/Dbq+worW1Vdja0YECs+MIegbxRBOoK1asGCnOoOhH3uYVINRYqGChUKgOXUlWPvbYY+FU6iGppuaCZDLUe4YhuemJcqPBnNAbwaHTQSAYuvRHn2tc6m6+83kJIpKqCgoeAOrAvdZZqXxFUZgkScV0zbBGQLuVVCrVi/nglrnRB0yH+3T10Q0/u3vRmY3JxWUGWd7Ua3mj69LVey6Vf1KJi5NaYnhiK10Fg8EQimIYfQ4jx6ysrFTj8fgAXzBpKpOMjTA28pdHf0v1lj4UkdIfyCPz5SiEQCbq4VctrRc+P+2Sky387LpSVS7Ti/9rvvvgVpzoWI+Rz4snqn+8v5fNspaKuL1F3d3pjdyztyfyDnS0U1pq+8gCU0v/29/+dhLr9PC8R/V7oGGMKnY6EjbycqPROOOBBx7YgIxD/mhd76MmQeDA787H+Pvuv//+Lvrdhl+x9pfFi0+1b77xdHHuJQuKk7ZF08XEtGpIupxK0d4D+vXbdJ/pKL35j2b1gOSKo2VulgLSQhKmCpMT0r5ftZSTSRwaYoMfbTY8HhcSRqMaee8X31tiCVx56jTctl2GkFY4BMLl3/81u+3H1eCpSCLAKU0CIeabQD5M1Jp4MIlmTLOcofBAiJMjW7QxIxM1yTs40kTukW4d5lGBfxpvUFAxfKYg75xjEVp45jkzmUzq5s17ijybftt6xXzxstMbdVCOZzB6gonUml3C476T739lVrnIPcGgfuO2nQsikZhQNe2kv1x12amRPe+9V8KQ30crC9KF7FBTXq/3BKTzH9/85jeDR/vm8Pd50yAZcGDD23/PPfd0ISGkXgHlgKoz1S5UVGlmFM02YuSVg+tR4BaDoaKM46qbWqXzzr7vvvtyNSDzH86U8fw+vbnq/92wWH/lkiZDqulEK4NZNbwmGhBhe+eXp+2645fV7XLHeNKd+DektR8S2tvb0s16/CmiNqJucTEGDSDYALXEUFv134wy6UTCI5rNlan9H7zReuVc42VnzZDkGSdYQJ1ezWsjAWl6o+OhrY7fd73tu/Dvlbix/ozTFgTC4UhdbWmJsmP9erPOaDykGztUtgQO1NZJbG/NCA4RwdGH8UghaLQP9c1Qz/ICkMPA0Um/ERwkATQNhtYfDPdLwMEghz04S1G52w3l0TgeNxaChDsBe4MC7d1g625F8/znYXy6cqjSY5oTaTXUh27DdKu+/1KJqIveg6vpPIGHg/GicgFt79FesVgWqwKSX0w92Np669UfrXtZ/JP9qoHKxe/HVFlDF6H/aUtL/zju9isbZ17UqHz0BBsvB1ndd9vDzl8hr2X0SYoKa2yNY4j8aEaRGlo2QtGZC6wtRp2omE1GodRSJcDMBoCYWS5NhYS4c/PlFz/U8WPK6NZbby1ubm7mt99114GxZPytb33LKstyDK9v4ndjBgfllXOAjAAOyl8DSSQwf5u3aNeejXZhfmVxPJlS3Lppdh83MBkNARXp7QOwE+Ny3LV6VLVKiY4hTGTQr2WDGhAbeZ0uXMGL3GGV9QRkaETzYqPJgyszOA7wyxz3pVS0P96uomNpOL193UAZSINmI9D0OK37XL2k7ubTauSnz5qpK6urR5mi5/Dysw239biVa6HN4UVCswGScTW0QeUUMAGVl8wxhRPhcjToFDt9KW7r7oMigxtHTVGmulPMGdfXo/CivOCJJ54woDbQp3+PmD9O9qBzihblkUcemYHjjukIrLfS341LGOUUIEcBR4ZnOPW6RlaU2q9usrNXUoqg7w4kwVKaYrgPWr/fo/7qh687X6JG0DL1nAIQ0xk06SK+aNKxww2m+rKkatC5pRmePjAIKSUUjglo80UAh1W4xpMpdBavqI1BNVfXN1cW8/+aVSOyympzCprRtL9EUEtqe86peq/3RdQeV+DIYcTGNUqaxtXQBqVNMywMqhJhb5Tbd7jY7LqyhKqXXNJ0DxqTIs/6wnFxkwN2npUe22I7wjPSjSTMsK1jIYYZ86Z7JgoOyptQczThpMUaBMuAQBpEw6hvcwaQUYKDCKUuh+B2u9+wWCznbrTDffu8/PTyInArXPnfDTud38T3apaELeWX7SBAZ2dcrrE+iBpwpU5QhWA8ouxxxWh7rX6vB3dK9TIqA+TCaBGnRLWNY5VG5eoitNCUOUPLfpOOl1UDrzBwIRhWioy+i1++3ToPjy/YzlHQsMnfCfgvnjmElZKAZMaixDNmRDu5Lr/i2+mQVxCQENGjAmQaHGoaHDMQHG+lu/ETqu+cAGQM4MgQr3W1XC7XO/jgOhdYTADoORlDNkReJpMcXTWA9/Y6V6Gbhq+8fUB86ICP1VaXcNAJfGt3UHn4ve0u0iCapM8RDRCThSpvGFQ8v5C5e2PQ1OkBFjAwcEVZJM5YkLNGzHt7jrTYWIt1KM86xYf2+1ltFY79kWebsav6hX9u693Tt3y5iDZDFHfEkCtwUKZZB8g4wJEpPEkK7G4RJlwREh1pg0MCz6ikSCahSbgSjWjZ63rWZLGs8rqkhcyuKOFA/VqAdZmjByhO1gOueWpBlpW9ezySsN2pyiZDEGKJFFSYBFVJRPn7B1Kptw6KWygieZfv/2LS/x/gWYnVurLPJZ4QSalhuc+5DilTW0mgTDI4iENZBcgEwJGprUzlYWcT+GCDw0yEKXwl2sWIy4UjTaA/7B84NZBnya0PJTlUIAnLcCbotz3B5Nf/tl+cGZdT3BHo49hNFVKqgFuP+fef+3N3Tz5MM4YicIRnGs/CTmcvxllN8UgwPvQQgmMU3cBcao4MzRMGCE11YmI8C+DI0ETXqa4xBtM6+F5rrPiANCFZFat5ADnxinV2dgaqq+sv3+VWnkGzzEvqSplUpOchtF7+8V8/cN6VnuTICCAib6qEI3g2VcBBDMoGQKiCtBXy9CLg4HWOqVIJ+aSD+EGVns+ggcTj6dmNmX4kWmmdj7Om1ohXxLHPQTu+pA1IUxEcGR6NiWf50BwZwjRJl/kxlmvGfokWqMh8pACOsXAvJ3GpkWljODznY3vE27OaIThIc6CKp3fHRaB1DpqdyvZs1XDMGTdAcClfISSfcsopy3DOuRPNRzpXrlwp0uowLcyM9W84AgvPx8QB0hIZoGj2WKPprowph8mNzGgRMF/goKKOq4tFeznQ2hVd7CinoGHfdkS0nRIj4un6IQo0/jok0EJWNkN6jHdEkkcxUyEiMoQcQeMRiQ3zYLgFuWGi5/oxCXNOK+TpRcCsrHMcjehxMe+uu+4qQg3y5aqqqneR2P1ouq1HsGQqZMQ8yYgMG5Gmucj0G7tmWA8siQs7gRE/HMNLAjCmOdDnJtNvMrDDa+QoDWvUuWD6FHdUZR51ooWIQ3Lg8ccfr0CDw1o0HQmjYJ6D9biGullDRs7yw3EBhLpWPT09TbNnz1awsY9VC5nRBLkcy6EVkEymsdB4YjGf8L6FDG8QrAYEYg2mOdCAyRwbweilaybeRK5YYXjOslxF10w6hHTMuxfzTmSeTfSKSRZhw5AwnwzdeCsEMN3QRNPOfD/VhRYaG5agML0EeeHctGnTO7S5LkN7rq8DlZvrjKZI+qS9skIKdX1QK5ky2imTKKYv4G7JrPEVJ0DKMU09gpGjptaywYZC3dvsFKSf8CkrtKgr39fXVxwOh4uefvrp17Ds2Sx3ptqGvY67IrHSxvXtOPvUwxZguBfIyKFe5ZW5QxFQeHYEB0YUWoPqcVLqbshWdEQRCg8mjQPjFURjIfhYEFrEh3xrj7HwsBC3wIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHCgwIECBwocKHDgQ8+B/w/2DAjEMdFZ/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<PIL.Image.Image image mode=RGBA size=200x200>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scm.plams import Molecule, Atom, packmol, preoptimize, view\n",
    "\n",
    "T_list = [50, 100, 200, 300]\n",
    "density = 1.0  # g/cm^3\n",
    "max_correlation_time_fs = 6000\n",
    "start_time_fit_fs = 3000\n",
    "\n",
    "Ar = Molecule()\n",
    "Ar.add_atom(Atom(symbol=\"Ar\", coords=(0, 0, 0)))\n",
    "\n",
    "mol = packmol(Ar, n_molecules=50, density=density)\n",
    "mol = preoptimize(mol, settings=get_lennard_jones_settings(), maxiterations=50)\n",
    "view(mol, height=200, width=200, direction=\"tilt_z\", padding=-2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8edb47106e1a46a883d545849b8ab81b",
   "metadata": {},
   "source": [
    "If you have multiple cores available, a parallel ``JobRunner`` can start the trajectories concurrently.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "10185d26023b46108eb7d9f57d49d2b3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scm.plams import config, JobRunner\n",
    "\n",
    "maxjobs = 4\n",
    "config.default_jobrunner = JobRunner(parallel=True, maxjobs=maxjobs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8763a12b2bbd4a93a75aff182afb95dc",
   "metadata": {},
   "source": [
    "To keep PLAMS from blocking too early, the production jobs use ``use_prerun=True`` and the MSD results are only accessed after all jobs have been submitted.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "06a3dcb1-1e50-474f-ab5b-08198569b6ed",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[23.03|14:56:29] JOB eq_T_50 STARTED\n",
      "[23.03|14:56:29] JOB prod_T_50 STARTED\n",
      "[23.03|14:56:29] JOB msd_T_50 STARTED\n",
      "[23.03|14:56:29] JOB eq_T_100 STARTED\n",
      "[23.03|14:56:29] Waiting for job eq_T_50 to finish\n",
      "[23.03|14:56:29] JOB prod_T_100 STARTED\n",
      "[23.03|14:56:29] JOB msd_T_100 STARTED\n",
      "[23.03|14:56:29] JOB eq_T_200 STARTED\n",
      "[23.03|14:56:29] JOB eq_T_50 RUNNING\n",
      "[23.03|14:56:29] Waiting for job prod_T_50 to finish\n",
      "[23.03|14:56:29] JOB prod_T_200 STARTED\n",
      "[23.03|14:56:29] JOB msd_T_200 STARTED\n",
      "[23.03|14:56:29] JOB eq_T_300 STARTED\n",
      "[23.03|14:56:29] Waiting for job eq_T_100 to finish\n",
      "[23.03|14:56:29] JOB prod_T_300 STARTED\n",
      "[23.03|14:56:29] JOB msd_T_300 STARTED\n",
      "[23.03|14:56:29] Waiting for job prod_T_100 to finish\n",
      "[23.03|14:56:29] Waiting for job msd_T_50 to finish\n",
      "[23.03|14:56:29] JOB eq_T_100 RUNNING\n",
      "[23.03|14:56:29] Waiting for job eq_T_200 to finish\n",
      "[23.03|14:56:29] Waiting for job prod_T_200 to finish\n",
      "[23.03|14:56:29] Waiting for job eq_T_300 to finish\n",
      "[23.03|14:56:29] Waiting for job prod_T_300 to finish\n",
      "[23.03|14:56:29] JOB eq_T_200 RUNNING\n",
      "[23.03|14:56:29] JOB eq_T_300 RUNNING\n",
      "[23.03|14:56:31] JOB eq_T_50 FINISHED\n",
      "[23.03|14:56:31] JOB eq_T_50 SUCCESSFUL\n",
      "[23.03|14:56:31] JOB prod_T_50 RUNNING\n",
      "[23.03|14:56:31] JOB eq_T_100 FINISHED\n",
      "[23.03|14:56:31] JOB eq_T_100 SUCCESSFUL\n",
      "[23.03|14:56:31] JOB prod_T_100 RUNNING\n",
      "[23.03|14:56:31] JOB eq_T_200 FINISHED\n",
      "[23.03|14:56:31] JOB eq_T_300 FINISHED\n",
      "[23.03|14:56:31] JOB eq_T_200 SUCCESSFUL\n",
      "[23.03|14:56:31] JOB eq_T_300 SUCCESSFUL\n",
      "[23.03|14:56:31] JOB prod_T_200 RUNNING\n",
      "[23.03|14:56:31] JOB prod_T_300 RUNNING\n",
      "[23.03|14:56:42] JOB prod_T_50 FINISHED\n",
      "[23.03|14:56:43] JOB prod_T_100 FINISHED\n",
      "[23.03|14:56:43] JOB prod_T_50 SUCCESSFUL\n",
      "[23.03|14:56:43] JOB msd_T_50 RUNNING\n",
      "[23.03|14:56:43] JOB prod_T_300 FINISHED\n",
      "[23.03|14:56:43] JOB prod_T_200 FINISHED\n",
      "[23.03|14:56:43] JOB prod_T_100 SUCCESSFUL\n",
      "[23.03|14:56:43] JOB msd_T_100 RUNNING\n",
      "[23.03|14:56:43] JOB msd_T_50 FINISHED\n",
      "[23.03|14:56:43] JOB msd_T_100 FINISHED\n",
      "[23.03|14:56:43] JOB prod_T_300 SUCCESSFUL\n",
      "[23.03|14:56:43] JOB prod_T_200 SUCCESSFUL\n",
      "[23.03|14:56:43] JOB msd_T_300 RUNNING\n",
      "[23.03|14:56:43] JOB msd_T_200 RUNNING\n",
      "[23.03|14:56:44] JOB msd_T_300 FINISHED\n",
      "[23.03|14:56:44] JOB msd_T_200 FINISHED\n",
      "[23.03|14:56:44] JOB msd_T_50 SUCCESSFUL\n",
      "[23.03|14:56:44] JOB msd_T_100 SUCCESSFUL\n",
      "[23.03|14:56:44] Temperature 50 K: D = 1.73e-09 m^2 s^-1\n",
      "[23.03|14:56:44] Temperature 100 K: D = 6.03e-09 m^2 s^-1\n",
      "[23.03|14:56:44] JOB msd_T_200 SUCCESSFUL\n",
      "[23.03|14:56:44] JOB msd_T_300 SUCCESSFUL\n",
      "[23.03|14:56:44] Waiting for job msd_T_200 to finish\n",
      "[23.03|14:56:44] Temperature 200 K: D = 1.19e-08 m^2 s^-1\n",
      "[23.03|14:56:44] Temperature 300 K: D = 1.75e-08 m^2 s^-1\n"
     ]
    }
   ],
   "source": [
    "from scm.plams import AMSNVTJob, AMSMSDJob, log\n",
    "\n",
    "msd_jobs = []\n",
    "DT_list = []\n",
    "for T in T_list:\n",
    "    eq_job = AMSNVTJob(\n",
    "        name=f\"eq_T_{T}\",\n",
    "        settings=get_lennard_jones_settings(),\n",
    "        molecule=mol,\n",
    "        temperature=T,\n",
    "        nsteps=10000,\n",
    "        samplingfreq=500,\n",
    "        thermostat=\"Berendsen\",\n",
    "        tau=100,\n",
    "        timestep=1.0,\n",
    "        writebonds=False,\n",
    "        writecharges=False,\n",
    "        writemolecules=False,\n",
    "        writevelocities=False,\n",
    "    )\n",
    "    eq_job.run()\n",
    "    prod_job = AMSNVTJob.restart_from(\n",
    "        eq_job,\n",
    "        name=f\"prod_T_{T}\",\n",
    "        nsteps=50000,\n",
    "        samplingfreq=25,\n",
    "        temperature=T,\n",
    "        thermostat=\"NHC\",\n",
    "        use_prerun=True,\n",
    "    )\n",
    "    prod_job.run()\n",
    "    msd_job = AMSMSDJob(\n",
    "        prod_job, name=f\"msd_T_{T}\", max_correlation_time_fs=max_correlation_time_fs\n",
    "    )\n",
    "    msd_job.run()\n",
    "    msd_jobs.append(msd_job)\n",
    "\n",
    "\n",
    "for T, msd_job in zip(T_list, msd_jobs):\n",
    "    D = msd_job.results.get_diffusion_coefficient(start_time_fit_fs=start_time_fit_fs)\n",
    "    DT_list.append(D)\n",
    "    log(f\"Temperature {T} K: D = {D:.2e} m^2 s^-1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7cdc8c89c7104fffa095e18ddfef8986",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[23.03|14:56:44] D0: 2.48e-08 m²/s. Barrier: 0.012 eV\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAADgCAYAAAAt6jC5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0PElEQVR4nO3deZxN9f/A8dd7zAyyb8mWpULGjGGGhrJUZKkspY32pPqStCj1/VbqV33TivgmKVIqpc23KPlGKGJGdhWJTKkmIUUY8/798TkzXde9Y+aaO3eW9/PxuI+595zPOed97r1z3vdsn7eoKsYYY0x+RUU6AGOMMcWTJRBjjDEhsQRijDEmJJZAjDHGhMQSiDHGmJBYAjHGGBMSSyDGGGNCYgnEGJNDROaIyFWRjsMUD5ZASigR2SIi+0Rkj4jsEpHPReRGEYnyaSMiMlpEdniP0SIieZx/ioh8LCK/iUiGiLwpInXyMN0pIvKXiLziN/xmEflORH4XkVQROSP/a50zr+oi8o6I/CkiW0VkQLiW5Tffq0XkdRFZIiI3+I27W0QeOYZ5dxGRLBH5w3v8ICIPHHvUh1PVnqr6UkHPN5uIDPA+kz9F5F0RqZ5L20QRSRORvd7fRJ9xZ4rIfBHZLSJb/KY7XkReE5EfvfGfichpIcRazvvfOSvAuKdFZGZ+51niqKo9SuAD2AJ09Z5XAXoD3wFTfNrcAHwN1AfqAeuBG/M4/57ARUBl4DjgReDDPEw3F1gEvOIz7DTgTyAJEOAmIAMoE+K6vwbMACoCZwC7gbhwLMtvueL9rQWk+Y1bDJxxDPPuAqT7vG4MpAN9Q5xftH/sQNQxxBedhzZxwB6gk/fZvAq8HqRtLLAVuBUoCwzzXsd649sBVwCDgS1+0zYBbgPqAGW8Nr8CFUNYr+eAqX7DygA/Aecf63emuD8iHoA9wvTB+iQQn2HtgCygpff6c2Cwz/jrgKUhLq8NsOcobS4F3gBG+SWQS4BlPq8rAArU8V6XBZ4Avgd+BiYC5YMsowJwAGjqM+xl4NG8LCvA/BYAD3nv1R/Af4EawHTgd2A50MinfRQwzXfjAlQDfvE2PDWB94FdwG+4ZHrUDbd/AvGGvQHc4/N6LLDNiysN6OgzbhQwE3jFGz/IW7eHgc+AfcDJ3rBBPtNdC2wAdgIfAQ19xikwBNgIfJeHdXgEeNXn9UneZ1UpQNtzgB/wkrI37Hugh1+7rvglkCDL/h1ICjIuChgJfAvs8N7X6t64Drikd5xP+17e53nUpFnSH3YIqxRR1WW4X60dvUFxwCqfJqu8YQCIyGr/wz+56ASsCzZSRCoDD+J+GfqbA5QRkdNEpAxuo7US9ysP4FGgKZCI28jVA+4LsqimQKaqfuMzzHe9jrasQC7F/dqth9voLQGmANVxG9f7vXWMASbj9sT+6zN9d+B/qnoIuB33GdQCagP34DbE+SIipwCnA0t9Bi/HvUfVcb/u3xSRcj7j++CSSFVcAoS/f8VXwv3C911GHy++C7x4F+H27nz1xe3VtfCm2ZXLIcHDvm+q+i1esg/SdrV6W2zPany+n3nlHfqKBTYFaXIzbj06A3VxyXKCF+PnwHbce5DtClwizMxvLCVNdKQDMIXuR9wGBtxhhN0+43YDFUVE1EnIywxFJAG3Qe+TS7P/A15Q1fQAp1n2AG/hDvMI7td5T1VV75zMYCBBVX/zlvcIbgN5d4DlVMT92vS1G7eBzHVZucQ+xdvYISJzgBaqOs97/aa3buASSVfgeBHpqapXeMPPBWZ7zw/iDq00VNVNuI1yXtUVkV24X8yVgHe89QBAVX3PKz0pIv8CmvH3RnuJqr7rPd/nfQ5TVTUn8ft9NjcC/1bVDd64R4B7RKShqmYnm39nfy5eDFVzid//+waHfzahtg3K++HyMvCAqvrPL9uNwFBVTfemGQV8LyJXeEliGnAl8Io3vz645F3q2R5I6VMPd+gE3CGZyj7jKgN/HGVjehgRORn3q/4WVQ24MfR+AXYFng4ym+uAa3C/LmOBy4H3RaQu7pfvcUCa9+t2F/ChNzz7qqHsE8sDA6xT9nrtycOygvnZ5/m+AK8rAqjqv1T1RFU9Lzt5eBctdPNiBngc90t4rohsFpGRuSzX34+qWlVVK+P2IvYBOSe8ReQOEdngnTjehTv3VdNn+m0B5hloWLaGwFif9/03XNKtl8fp/R3tswm1bUAiUh53yHGpqv47l6YNgXd81nMDcAi3hwguAZ3pfUf6A9+q6pd5jaMkswRSiohIW9w/f/av1nVAK58mrcjlMFSA+TUE5gH/p6ov59K0C9AI96vuJ+AO4EIRWeGNTwTeV9VvVDVLVT/EHTbogDv5uQ93Eryq96iiqtkb7Z6qWtF7TAe+AaK9QzyB1iu3ZYVDW2CrqmZ48e5R1dtVtQnuwobbROTs/M7U+zX9KnA+gIh0BO4ELgaqeXsCu3Eb/JzJAs0ql8VsA27wed+rqmp577BOXqb3d9j3TUSa4M5vfROkbYIcvkuUQB6/nyJSFngXd7jwhtxbsw23F+q7nuVU9QcAb29rEe7HxhX4JO3SzhJIKSAilUXkPOB13MnrNd6oabgNWD3v19XtwNQ8zrMe8AkwXlUnHqX5JNy5g0TvMRH4AHduANyx+3NFpIk43XDHxdeqahbwPPC0iByfvWwR6U4Aqvon8DbwoIhUEJHTcYccshNc0GXlZb1D0Au3rnixnyciJ3sbxt24X7pZ+Z2piFTEnZvJ3qBWAjJxV5RFi8h9HPkLPr8mAneLSJy3zCoictExzG86cL6IdBSRCrhzYm+raqC9igW492aYiJQVkaHe8E+8WKK88zsx7qWUE5FYb1wM7lzPPuAq7zt0tPV82PtBhIjU8s7/+HoJGIo7dDUd40T6LL49wvPAXYW1D7fLvxt34ncIPper4n6dPoY7NPGb99z3qpd1wMAg878f9+vzD9+Hz/h7gDlBph3F4VdhCW5j8r0X7wbgCp/x5XBX8GzGnd/YAAzLZd2r4359/unNc0BelxVgXgs4/Kqkh/C5rBN3aG5TLtOnAsk+r2/1Pps/cb+O7/UZNwefq6r85tMFl2iy3+sduMR0sje+DO5S6t9xe1R3cvil3Ie954HWLcj6XgGs8ea7DXjRZ5xmL99n2B/4XP0VYD0GeO/9n8B7eFc7BVp/oDXuarJ9wAqgtd/7oX6PBd64zt7rvRz+/QwYF+6H9G24S9r34K7GesSvTUVvHgG/06X1kX3dujGmgIlIbeBLoJ7aP5opgewQljHhUwW43ZKHKalsD8QYY0xIbA/EGGNMSCyBGGOMCUmpuhO9Zs2a2qhRo0iHYYwxxUpaWtqvqlrLf3ipSiCNGjUiNTU10mEYY0yxIiJbAw23Q1jGGGNCEpEEIiIXicg6cQVykn2GDxSRlT6PLPEpIuPTbpS4gjrZ7XoV6goYY4yJ2CGstbjukZ/zHaiuL6PpACISD7yrqiuDzONpVX0inEEaY4wJLiIJRP/uHjq3Zpfh+m4yxkTAwYMHSU9P56+//op0KKaQlCtXjvr16xMTE5On9kX5JPol5F5fYqiIXInra+h2Vd0ZtkjWvQPpqdBhGFSqffT2xpQA6enpVKpUiUaNGh3tx54pAVSVHTt2kJ6eTuPGjfM0TdjOgYjIPBFZG+CRW1LInvY0YK+qBush9Vn+7t11O/BkLvMaLCKpIpKakZERwpoAP62Fpf+BsQkw+07Y/UNo8zGmGPnrr7+oUaOGJY9SQkSoUaNGvvY4w7YHoqpdj2HySzmydKbvvHMK+ojI87ga08HaTsJ1J05ycnJo/bacfS8kDoBFT0HqC5D6IrQeCGfcCtUahTRLY4oDSx6lS34/7yJ3Ga9Xwe1icjn/ISJ1fF72I3y1HP5W4yToOwFuXgGtL4eVr8K4NvDuP2DHt2FfvDGl0cMPP0xcXBwJCQkkJibyxRdfADBo0CDWr19fIMto1KgRv/76a65tHnnkkcNed+gQrvpjBaNixYoAbNmyhVdffTVsy4nUZbz9RCQdaA98ICIf+YzuBGxT1c1+00z2ueT3MRFZIyKrgTNxNRYKR7WGcP4YGLYS2l0Pa9+C8ckw8zr4ZUOhhWFMSbdkyRLef/99VqxYwerVq5k3bx4NGjQAYPLkybRo0aLQYvFPIJ9//nmQlqHLzMws8HmWyASiqu+oan1VLauqtVW1u8+4BaqaEmCaQaqa6j2/QlXjVTVBVXur6vbCjB+AKvWg52i4ZTW0HwJfz4H/pMCMK2D7qkIPx5iiIG3rTibM30Ta1mO/pmX79u3UrFmTsmXLAlCzZk3q1nWl67t06ZLTq0TFihUZMWIEcXFxdO3alWXLltGlSxeaNGnCrFmzAJg6dSpDhw7Nmfd5553HggULjlhm3759SUpKIi4ujkmTJgEwcuRI9u3bR2JiIgMHDsxZJrgTzyNGjKBly5bEx8czY8YMABYsWECXLl3o378/zZs3Z+DAgQTq+bxLly4MHz6c5ORkxo4dS1paGp07dyYpKYnu3buzfbvbtI0bN44WLVqQkJDApZdeCsCoUaN44om/72Ro2bIlW7ZsOWz+I0eOZNGiRSQmJvL000+zbt062rVrR2JiIgkJCWzcuDEfn0gAka5oVZiPpKQkDZs/flWd96DqI/VV76+sOv1i1W2p4VueMWG2fv36fLVP3fKbNvvXbG088n1t9q/Zmrrlt2Na/p49e7RVq1Z6yimn6E033aQLFizIGde5c2ddvny5qqoCOnv2bFVV7du3r3br1k0PHDigK1eu1FatWqmq6pQpU3TIkCE505977rk6f/58VVVt2LChZmRkqKrqjh07VFV17969GhcXp7/++quqqlaoUOGw2LJfz5w5U7t27aqZmZn6008/aYMGDfTHH3/U+fPna+XKlXXbtm166NAhTUlJ0UWLFh2xjp07d9abbrpJVVUPHDig7du3119++UVVVV9//XW95pprVFW1Tp06+tdff6mq6s6dO1VV9f7779fHH388Z15xcXH63XffHRbf/Pnz9dxzz81pM3ToUH3llVdUVXX//v26d+/eI2IK9LkDqRpgm1rkzoEUWxVquJPtw9fAmf+E75fC5LPg5X6wteB3d40papZu3sGBzCyyFA5mZrF0845jml/FihVJS0tj0qRJ1KpVi0suuYSpU6ce0S42NpYePXoAEB8fT+fOnYmJiSE+Pv6IX+RHM27cOFq1akVKSgrbtm076i/0xYsXc9lll1GmTBlq165N586dWb58OQDt2rWjfv36REVFkZiYGDSWSy65BICvv/6atWvX0q1bNxITE3nooYdIT08HICEhgYEDB/LKK68QHR36tU/t27fnkUceYfTo0WzdupXy5cuHPC8ogifRi73yVaHznXDrWug6Cravhik9Ycq5sHkBWAEvU0KlNKlBbHQUZQRioqNIaVLjmOdZpkwZunTpwgMPPMD48eN56623jmgTExOTc/VQVFRUziGvqKionPMK0dHRZGVl5UwT6FLVBQsWMG/ePJYsWcKqVato3br1Md1EmR1H9noEO8dRoUIFwB0NiouLY+XKlaxcuZI1a9Ywd+5cAD744AOGDBnCihUraNu2LZmZmXlaJ38DBgxg1qxZlC9fnl69evHJJ5+EvH5gCSR8ylZyl/kOXwPd/w07NsG0PvDCOfDNXEskpsRJaliN6YNSuO2cZkwflEJSw2rHNL+vv/76sD2AlStX0rBhw5Dm1ahRI1auXElWVhbbtm1j2bJlR7TZvXs31apV47jjjuOrr75i6dKlOeNiYmI4ePDgEdN07NiRGTNmcOjQITIyMli4cCHt2rULKcZmzZqRkZHBkiVLANcTwLp163JiPvPMMxk9ejS7d+/mjz/+oFGjRqxYsQKAFStW8N133x0xz0qVKrFnz56c15s3b6ZJkyYMGzaMPn36sHr16pBizVaU70QvGWKPg/b/gORrYeUrsHgMvHoR1EmETiOgWS+IsjxuSoakhtWOOXFk++OPP7j55pvZtWsX0dHRnHzyyTkntvPr9NNPp3HjxrRo0YJTTz2VNm3aHNGmR48eTJw4kVNPPZVmzZqRkvL3tTyDBw8mISGBNm3aMH369Jzh/fr1Y8mSJbRq1QoR4bHHHuOEE07gq6++yneMsbGxzJw5k2HDhrF7924yMzMZPnw4TZs25fLLL2f37t2oKsOGDaNq1apceOGFTJs2jbi4OE477TSaNm16xDwTEhIoU6YMrVq14uqrr2b//v28/PLLxMTEcMIJJ3DPPffkO05fpaomenJyska8HkjmAVj9urspced3cHwcdLodWvSFqDKRjc0YHxs2bODUU0+NdBimkAX63EUkTVWT/dvaT9/CFh0Lba6EoanQbxJkHYSZ18KE02Dla3Co4K8FN8aYcLAEEilloqHVJfCPpXDRVIguC+/eCOOTIO0lt6dijDFFmCWQSIsqA3H94IZFcOmrUL4a/HcYjGsNy56Hg9aVtjGmaLIEUlRERUHzc+H6+TDwLXen++w7YGwrWDIBDvwZ6QiNMeYwlkCKGhE4pStc+xFcOQtqngIf3QNjEtyJ9/17jj4PY4wpBJZAiioRaNIZrn7fJZO6ifC/B+DplrBgNOwLX/0sY4zJC0sgxcGJKXD5W3D9J9CwAyx4xO2R/O9B+DN/3UUUZGd3xoRbdqeFviZOnMi0adMKNY5FixYRFxdHYmIiP/zwA/379wfczY2zZ88u1FiKEruRsDiplwSXvQY/rYGFT7hDWkufdTcp5qHcbtrWnQycvJQDmVnERkcVyN3CxhS2G2+8Mazzz+ko0OcG3+nTp3P33Xdz+eWXAzBz5kzAJZDU1FR69eoV1piKKtsDKY5OiIeLX3KXADc/L8/ldgu6sztjIsG3G/MuXbpw11130a5dO5o2bcqiRYsAOHToECNGjKBt27YkJCTw3HPPAe7u9rPPPps2bdoQHx/Pe++9B7i6Gc2aNePKK6+kZcuWbNu2LWd5kydP5o033uDee+9l4MCBbNmyhZYtW3LgwAHuu+8+ZsyYQWJiYk5X7qWJ7YEUZ8c3hwufhy4jYfHRy+1md3Z3MDOrwDq7M6XEnJFuz7cgnRAPPR895tlkZmaybNkyZs+ezQMPPMC8efN44YUXqFKlCsuXL2f//v2cfvrpnHPOOTRo0IB33nmHypUr8+uvv5KSkkLv3r0B2LhxIy+99NJhXZiAq364ePFizjvvPPr375/Tq25sbCwPPvggqampjB8//pjXoziyBFIS1DgJ+kyATnfCZ2Phy5dhxcvQ6lI44zaoeTLwd2d3SzfvIKVJDTt8ZUqECy64AICkpKScjfvcuXNZvXp1zqGm3bt3s3HjRurXr88999zDwoULiYqK4ocffuDnn38GoGHDhkckD5M7SyAlSbWGcN5T0OkO+GwcpE2BVa9B3AVu2PGnFmhnd6YUKYA9hXDJ7jbdt8t0VeWZZ56he/fuh7WdOnUqGRkZpKWlERMTQ6NGjXK6Qc/uVt3kXcTOgYjIRSKyTkSyfGqdIyIxIvKSV/N8g4jcHWT6xiLyhYhsEpEZIhJbeNEXcZXrun/44Wugw80+5XYvt3K7plTo3r07zz77bE4X7N988w1//vknu3fv5vjjjycmJob58+ezdevWY1qOf3fppU0kT6KvBS4AFvoNvwgoq6rxQBJwg4g0CjD9aOBpVT0Z2AlcF8ZYi6eKx0O3B11xq04jYPOn8FwnePUSSI9wr8TG5MHevXupX79+zuOpp57K03SDBg2iRYsWtGnThpYtW3LDDTeQmZnJwIEDSU1NJT4+nmnTptG8efNjiu/MM89k/fr1pfYkesS7cxeRBcAdqprqvb4MGAD0A6oAS4AUVf3NZxoBMoATVDVTRNoDo1S1u//8fRWJ7twjad8u17/W0gnuRsQmZ7rqiQ07RDoyUwRZd+6lU3Hvzn0m8CewHfgeeMI3eXhqALtUNbvv83SgXuGFWEyVrwqdR8DwtW7P5Oe1Vm7XGBOysCYQEZknImsDPPrkMlk74BBQF2gM3C4iTY4hhsEikioiqRkZGaHOpmQpWxFOvwVuWQ09HoXfvvXK7XazcrvGmDwL61VYqto1hMkGAB+q6kHgFxH5DEgGNvu02QFUFZFoby+kPhDwDjpVnQRMAncIK4R4Sq7Y4yDlJncn+5e+5XZbeeV2z7Vyu8aYoIri1uF74CwAEakApACHFRhWd+JmPtDfG3QV8F4hxliyRJeFttfBsBXufpL9e9wVWxNPh7VvQdahSEdoIiTS50hN4crv5x3Jy3j7iUg60B74QEQ+8kZNACqKyDpgOTBFVVd708wWkbpeu7uA20RkE+6cyAuFuwYlUJkYaH05DFkOFzzvEoeV2y21ypUrx44dOyyJlBKqyo4dOyhXrlyep4n4VViFqdRfhZVfWVmwYZbruPHnNa5rlDNug1aXudrupkQ7ePAg6enpOTfamZKvXLly1K9fn5iYmMOGB7sKyxKIOTpVdzPiwsfgxy+hcn04Yzi0vgJi8v5rxRhTPBWny3hNUSMCzXu5cruXvwVV6nvldhPg8/FWbteYUsoSiMk7ETi5K1z7IVz1X6jVDOb+8+9yu3/9HukIjTGFyBKIyT8RaNzJJZFr5/5dbndMPCx41MrtGlNKWAIxx+bE07xyu/Oh4emw4N8hl9s1xhQvlkBMwajXBi57FW78DE46yx3SGtMSPvon7Pk50tEZY8LAEogpWCe0dOV2h3wBp57vyu2OiYfZI2B3eqSjM8YUIEsgpsClbd3JhLVlSEsaDUNTIeFiV2p3bCL89xbYuSXSIRpjCoBVJDQFKm3rTgZOXsqBzCxio6OYPiiFpD7jXbfxi8f8XW434RLoeHtOuV1jTPFjeyCmQC3dvIMDmVlkKRzMzGLpZu9EetUTXbndW1bDaTfAundgQluYeR38vD6yQRtjQmIJxBSolCY1iI2OooxATHQUKU1qHN6gch3o8W+v3O4w+OZDeLa9lds1phiyrkxMgUvbupOlm3eQ0qQGSQ2r5d5472+w9Fn44jnYvxtO6e4Od9U/otcEY0yEWF9YWAIp0v7aDcsmwRKfcrudRkCj0yMdmTGlnvWFZYq2clVcwvAttzu1F0zpBd/OtyqJxhRBlkBM0ZJdbnf4GugxGn77Dl7u65Xb/cgSiTFFiCUQUzTFlIeUG+GWlXDe0+5u9lcvhuc6wfpZrlaJMSaiLIGYoi26rKvZnl1u98Af8MYVrtzumplWbteYCLIEYoqHw8rtTgbNgreugwntYOWrcOhgpCM0ptSxBGKKlzLRkHAR3LQELnoJosvDuzfBM0mQNhUyD0Q6QmNKjYgkEBG5SETWiUiWiCT7DI8RkZdEZI2IbBCRu4NMP1VEvhORld4jsdCCN0VDVBTE9YUbF8Flr8NxNVw/W+MS4YtJcHBfpCM0psSL1B7IWuACYKHf8IuAsqoaDyQBN4hIoyDzGKGqid5jZdgiNUWbCDTrCdd/Ape/7bpMmTMCxraCz5+xcrvGhFFEEoiqblDVrwONAiqISDRQHjgAWJ1Uc3QicPLZcM0cuOp9qNUc5v7LdSW/6Ekrt2tMGBS1cyAzgT+B7cD3wBOq+luQtg+LyGoReVpEygaboYgMFpFUEUnNyMgIQ8imSBGBxh3hqllw3cdQt42rjjimJcz/t5XbNaYAhS2BiMg8EVkb4NEnl8naAYeAukBj4HYRaRKg3d1Ac6AtUB24K9gMVXWSqiaranKtWrVCXyFT/DRoB5fPhMELoFFH+PRReDoe5o2CP3+NdHTGFHthqweiql1DmGwA8KGqHgR+EZHPgGRgs9+8t3tP94vIFOCOYwrWlGx1W8Ol0+Gnte5w1uIxrvPG5Guhw81Q6YRIR2hMsVTUDmF9D5wFICIVgBTgK/9GIlLH+ytAX9xJeWNyd0JLuGgKDFkGp/Z2vQCPSYAP7rByu8aEIFKX8fYTkXSgPfCBiHzkjZoAVBSRdcByYIqqrvammS0idb1200VkDbAGqAk8VLhrYIq1Wk3hgufg5lRodYm7f2RsIsy62fW9ZYzJE+vO3Zhd2+CzMa7Ublamq+He8XaoeUqkIzOmSLDu3I0JpmoDOPdJuGUVnHYjrHsXxreFN6+xcrvG5MISiDHZKteBHo+4ruTPGA4b57pyu68PhB9XRjo6Y4ocSyDG+KtYC7qOcomk813w3SKY1BmmXwTblkc6OmOKDEsgxgRzXHU48x64dQ2cdS+kp8ILXeGl3rBlcaSjMybiLIEYczTlqkCnO9weyTkPwS8bYOq58GJP+PYTq5JoSi1LIMbkVdmK7sbD4auh52Owcwu83A8md4WvP7REYkqdXC/j9To17InrNgRgA+5O8cxCiK3A2WW8pkBl7nfFrBY/Bbu+hxPiodMIaH6+627emBIi35fxikg9YB1wO65vqnrACGCdzw19xpRe0WUh+Rq4eQX0fRYO7IU3roRnO1i5XVMqBN0DEZGpwEpVHeM3fBiQpKpXhT26AmZ7ICassg7Bundg4ROQsQGqn+RuSEy42JXkNaaYCrYHklsC+UpVmwcZ97WqNivgGMPOEogpFFlZ8NX7sPAx+GmNK3J1xq2QONDttRhTzIRyJ3puNUH3HntIxpRQUVHQojfcsAgumwEVasH7t8K41q4XYCu3a0qI3LpzryIiFwQYLkDlMMVjTMkhAs16QNPusHk+fPo4zLnTHeLqcLPrTr5sxUhHaUzIcjuENSW3CVX1mrBEFEZ2CMtE3JbP3KGtzQugfHVoPwTaXe/uNTGmiMr3OZCSyBKIKTK2LXN7Ihs/csnjtBvd47jqkY7MmCNYb7zGFCUN2sHAN2Dwp1653dEwJh4+vh/+yIh0dMbkiSUQYyKpbqIrt3vT5+5cyWdjXSL58B7Y81OkozMmV5ZAjImwtK07mbC+LGltn4ShyyGuH3wx0Su3e7sreGVMEZTbVVg5RKQD0Mi3vapOC1NMxpQaaVt3MnDyUg5kZhEbHcX0QSkk9XsWOt8Ji5+GtJfcI/Eydy9J9SaRDtmYHEfdAxGRl4EngDOAtt7jiJMp+SEiF4nIOhHJEpFkn+GxIjJFRNaIyCoR6RJk+uoi8rGIbPT+VjuWeIyJlKWbd3AgM4sshYOZWSzdvMONqN4Yeo+DYV+67lJWzYBnkuHtGyDjm8gGbYwnL3sgyUALLdjLtdYCFwDP+Q2/HkBV40XkeGCOiLRV1Sy/diOB/6nqoyIy0nt9VwHGZ0yhSGlSg9joKA5mZhETHUVKkxqHN6jaAHo97rpE+fwZSH0RVs+AuL6u48bacRGJ2xjIWwJZC5wAbC+oharqBgAR8R/VAvjEa/OLiOzCJbBlfu36AF285y8BC7AEYoqhpIbVmD4ohaWbd5DSpAZJDYPsTFc6Abo/7A5jLZkAy553/W41P8/VKqnbunADN4a8JZCawHoRWQbszx6oqr3DEM8qoLeIvAY0AJK8v/4JpLaqZie0n4DaYYjFmEKR1LBa8MThr0JN6Hq/u5P9i+fgi2ddv1snd3PnTRq0C2+wxvjISwIZFcqMRWQebs/F3z9V9b0gk70InAqkAluBz4Fc+8RWVRWRoIfXRGQwMBjgxBNPzEPkxhQDx1WHM+92d7Ivf97tlbzQDRp3gk53QqMzXFcqxoRRRO9EF5EFwB2qGvD2cBH5HBikquv9hn8NdFHV7SJSB1iQl96B7U50U2Id+NOdH/lsHPz5C5zY3p0jOeksSyTmmIVSUGqPiPwe4LFHRH4PU5DHiUgF73k3INM/eXhmAdn1SK4Cgu3RGFM6xFbwKbf7uKuQ+MoFMPls+HqOlds1YRGRPRAR6Qc8A9QCduEKV3UXkUbAR0AW8ANwnapu9aaZDExU1VQRqQG8AZyIO9R1sar+drTl2h6IKTUy98Oq12DRU7BrK9SOdyfbT+1t5XZNvllnilgCMaXQoYOw5k1Y9CTs2AS1mkPHO6DlBRBVJtLRmWLCOlM0pjQqEwOJA2DIMrjwBUDg7UEwvi18+YpLMMaEyBKIMaVBVBmI7+86bbzkFXfO5L0h8Ewbd/I9c//R52GMH0sgxpQmUVFw6vlww0IY8AZUON6V2x2bCEsnWrldky+WQIwpjURc9/GD5sEV77q+tz68y/UA/Nk42P9HpCM0xYAlEGNKMxE46Uy4ZjZcPdv1rfXxva4mycLH4a/dkY7QFGGWQIwxTqPT4cp34bp5UL8tfPIQPB0PnzwMe496lbwphSyBGGMO16CtK7d7w0Jo0gkWPmbldk1AlkCMMYHVaeWu2LppiV+53bvh9wLrnNsUY5ZAjDG5q90C+r/oU273ORibAO/f5rpMMaWWJRBjTN7UPAX6PQs3p7mbE1dMg3Gt4b2h8NvmSEdnIsASiDEmf6o3hvPHwi0rIflaWP0GPJMEbw+2cruljCUQY0xoqtR35XaHr4aUf8CG/8KEdvDm1fDT2khHZwqBJRBjzLHJLrc7fI0rubtxHkw8HV4bAD9+GenoTBhZAjHGFIzscru3roEud8PWxTCpC7zSH77/ItLRmTCwBGKMKVjlq0GXkTB8LZx9P/y4Al48B146H75baMWtShBLIMaY8ChXGTre5g5tnfMwZHztksiLPWDTPEskJYAlEGNMeMVWgA5D4ZbV0OsJ2J0Or1wIz58FX822RFKMWQIxxhSOmHLQ7noY9qW7DHjvDnj9MpjYEda9C1lZkY7Q5JMlEGNM4YqOhaSr4eYV0HciZO6DN6+C/6S4e0oOZUY6QpNHEUkgInKRiKwTkSwRSfYZHisiU0RkjYisEpEuQaYfJSI/iMhK79GrsGI3xhSQMtGQeJkrt9v/RVc18e3rYYKV2y0uIrUHsha4AFjoN/x6AFWNB7oBT4pIsBifVtVE7zE7fKEaY8Iqqgy0vBBu/Mx13li2kiu3O64NLH/Byu0WYRFJIKq6QVW/DjCqBfCJ1+YXYBeQHKCdMaakyS63O/hTGPAmVKoNH9wGY1vB0mfhwN5IR2j8FLVzIKuA3iISLSKNgSSgQZC2Q0VktYi8KCLVgs1QRAaLSKqIpGZkWC0DY4o8EWh6Dlz3MVz5HlQ/CT4c6XoA/mysldstQsKWQERknoisDfDok8tkLwLpQCowBvgcOBSg3bPASUAisB14MtgMVXWSqiaranKtWrVCXBtjTKETgSZd4JoP4Jo5ULslfHwfjGkJn1q53aIgOlwzVtWuIUyTCdya/VpEPgeO6N5TVX/2afM88H6IYRpjioOGHVy53fRUV6t9/kPw+TNw2mDXkeNx1SMdYalUpA5hichxIlLBe94NyFTV9QHa1fF52Q93Ut4YU9LVT4YBM3zK7T7uldu9z8rtRkCkLuPtJyLpQHvgAxH5yBt1PLBCRDYAdwFX+Ewz2eeS38e8S31XA2fis9dijCkFssvt/mMpNO3h9kbGxMOckfD7j5GOrtQQLUXdCCQnJ2tqamqkwzDGFLRfN8Hip2DV6+6y4NZXwBnDoeqJkY6sRBCRNFU94orYInUIyxhjQlLzZOj7Hxi2wq/c7hDY8W2koyuxLIEYY0qOao28crurIPk6WDMTxid75XYD3XpmjoUlEGNMyVOlHvR6zPUA3H6IV273NHjjKiu3W4AsgRhjSq5KteGch1xxq463wab//V1u94cVkY6u2LMEYowp+SrUgLPv88rt3gNbP4Pnz3R1SazcbsgsgRhjSo/y1aDLXa5K4tn3w49funK7U8+zcrshsARijCl9fMvtdn8Eft34d7ndjVZuN68sgRhjSq/YCu4k+y2r/i63O/1Cd3jLyu0elSUQY4w5rNzuONi30yu3ewase8fK7QZhCcQYY7JFx0LSVTA0Dfo954pZvXm1K7e7aoaV2/VjCcQYY/yViYZWl8KQL6D/FIiKhncGu3K7K16GzAORjrBIsARijDHBRJWBlhfAjYvhkumu3O6sofBMG1g+udSX27UEYowxRxMVBaee58rtDpwJlerAB7eX+nK7lkCMMSavROCUbnDd3CPL7S4eA/v3RDrCQmUJxBhj8su/3O4J8TDvfleT5NPHYd+uSEdYKCyBGGPMsWjYAa54BwZ9Ag1SXLndMfHwyUOw97dIRxdWlkCMMaYg1E+CAa/DDYvc3snCx+HpljD3Xvjjl0hHFxaWQIwxpiDVSYBLXnbldpv3giXjYUxCiSy3G7EEIiKPi8hXIrJaRN4Rkao+4+4WkU0i8rWIdA8yfWMR+cJrN0NEYgsteGOMOZrjT4ULJ8OQ5dDyQlj+vLtq6/1bYdf3kY6uQERyD+RjoKWqJgDfAHcDiEgL4FIgDugB/EdEygSYfjTwtKqeDOwEriuUqI0xJj9qngx9J8DNKyBxoLsRcVxreLf4l9uNWAJR1bmqmt0vwFKgvve8D/C6qu5X1e+ATUA732lFRICzgJneoJeAvmEP2hhjQlWtIZw/xnXc2HYQrPXK7b51fbEtt1tUzoFcC8zxntcDtvmMS/eG+aoB7PJJQIHaGGNM0VOlHvQc/Xe53a8+8MrtXgk/rYl0dPkS1gQiIvNEZG2ARx+fNv8EMoHpYYphsIikikhqRkZGOBZhjDH5l1Nudw10vB2+ne96/33tMvghLdLR5Ul0OGeuql1zGy8iVwPnAWer5nS8/wPQwKdZfW+Yrx1AVRGJ9vZCArXJjmESMAkgOTnZOvc3xhQtFWrA2fdCh5th2SRYMgGePwtOOhs63wknpkQ6wqAieRVWD+BOoLeq+nYkMwu4VETKikhj4BRgme+0XrKZD/T3Bl0FvBf+qI0xJkzKV3UJ49a10HUUbF8FL3Z35XY3f1oki1tF8hzIeKAS8LGIrBSRiQCqug54A1gPfAgMUdVDACIyW0TqetPfBdwmIptw50ReKOwVMMaYAle2EpxxKwxfDd3/7crtTuvtksnGj4tUIhEtQsGEW3JysqampkY6DGOMybuDf8GXL7vOGn9Ph7qtodMIaNrT9RJcCEQkTVWT/YcXlauwjDHGBOJbbrf3M1653QHwXEdY+zZkHYpYaJZAjDGmOIiOhTZXHl5ud+Y1ES23awnEGGOKkyPK7ca4crvjk2HFtEItt2sJxBhjiiPfcruXvgrlqsCsm/8ut3vwr5ymaVt3MmH+JtK27izQEOwkujHGlASqsGkefPoYpC9zZXc7DGPF8X0ZMHUVBzKziI2OYvqgFJIaVsvXrO0kujHGlGSHldudBTVOho/upvnrHbhG36W87uNgZhZLN+8osEVaAjHGmJJEBJp0hqvfh2s+5ODx8dwV/Tody6wlJjqKlCY1CmxRYe3KxBhjTAQ1bE+V6//L+hWLiN9Vm+tPqpXvw1e5sQRijDElXIs2HWkRhvnaISxjjDEhsQRijDEmJJZAjDHGhMQSiDHGmJBYAjHGGBOSUnUnuohkAFtDnLwm8GsBhlNQLK78sbjyx+LKn6IaFxxbbA1VtZb/wFKVQI6FiKQGupU/0iyu/LG48sfiyp+iGheEJzY7hGWMMSYklkCMMcaExBJI3k2KdABBWFz5Y3Hlj8WVP0U1LghDbHYOxBhjTEhsD8QYY0xILIH4EZEeIvK1iGwSkZEBxpcVkRne+C9EpFEhxNRAROaLyHoRWScitwRo00VEdovISu9xX7jj8pa7RUTWeMs8olqXOOO892u1iLQphJia+bwPK0XkdxEZ7temUN4vEXlRRH4RkbU+w6qLyMcistH7G7B7VBG5ymuzUUSuKoS4HheRr7zP6R0RqRpk2lw/8zDENUpEfvD5rHoFmTbX/90wxDXDJ6YtIrIyyLThfL8CbhsK7TumqvbwHkAZ4FugCRALrAJa+LX5BzDRe34pMKMQ4qoDtPGeVwK+CRBXF+D9CLxnW4CauYzvBcwBBEgBvojAZ/oT7jr2Qn+/gE5AG2Ctz7DHgJHe85HA6ADTVQc2e3+rec+rhTmuc4Bo7/noQHHl5TMPQ1yjgDvy8Dnn+r9b0HH5jX8SuC8C71fAbUNhfcdsD+Rw7YBNqrpZVQ8ArwN9/Nr0AV7yns8EzhYRCWdQqrpdVVd4z/cAG4B64VxmAeoDTFNnKVBVROoU4vLPBr5V1VBvID0mqroQ+M1vsO936CWgb4BJuwMfq+pvqroT+BjoEc64VHWuqmZ6L5cC9QtqeccSVx7l5X83LHF5//8XA68V1PLyKpdtQ6F8xyyBHK4esM3ndTpHbqhz2nj/bLuBgivxdRTeIbPWwBcBRrcXkVUiMkdE4gopJAXmikiaiAwOMD4v72k4XUrwf+xIvF8AtVV1u/f8J6B2gDaRft+uxe05BnK0zzwchnqH1l4Mcjgmku9XR+BnVd0YZHyhvF9+24ZC+Y5ZAilGRKQi8BYwXFV/9xu9AneYphXwDPBuIYV1hqq2AXoCQ0SkUyEt96hEJBboDbwZYHSk3q/DqDuWUKQuhRSRfwKZwPQgTQr7M38WOAlIBLbjDhcVJZeR+95H2N+v3LYN4fyOWQI53A9AA5/X9b1hAduISDRQBSi4KvVBiEgM7gsyXVXf9h+vqr+r6h/e89lAjIjUDHdcqvqD9/cX4B3coQRfeXlPw6UnsEJVf/YfEan3y/Nz9mE87+8vAdpE5H0TkauB84CB3obnCHn4zAuUqv6sqodUNQt4PsjyIvV+RQMXADOCtQn3+xVk21Ao3zFLIIdbDpwiIo29X6+XArP82swCsq9W6A98EuwfraB4x1hfADao6lNB2pyQfS5GRNrhPtuwJjYRqSAilbKf407CrvVrNgu4UpwUYLfPrnW4Bf1lGIn3y4fvd+gq4L0AbT4CzhGRat4hm3O8YWEjIj2AO4Heqro3SJu8fOYFHZfvObN+QZaXl//dcOgKfKWq6YFGhvv9ymXbUDjfsXBcGVCcH7irhr7BXdHxT2/Yg7h/KoByuEMim4BlQJNCiOkM3C7oamCl9+gF3Ajc6LUZCqzDXX2yFOhQCHE18Za3ylt29vvlG5cAE7z3cw2QXEifYwVcQqjiM6zQ3y9cAtsOHMQdY74Od87sf8BGYB5Q3WubDEz2mfZa73u2CbimEOLahDsmnv0dy77asC4wO7fPPMxxvex9d1bjNox1/OPyXh/xvxvOuLzhU7O/Uz5tC/P9CrZtKJTvmN2JbowxJiR2CMsYY0xILIEYY4wJiSUQY4wxIbEEYowxJiSWQIwxxoTEEogxxpiQWAIxxk+grru94QG7yPZukgzYZX1u3WWLyEgR+adPl+CHfJ4PCxBXX/G6nRfXxfkd3vNyXjyjgqxPJREZLSLzfYa9LiKnHMPbZIwlEGMCmErgXklHAv9T1VNwN2ll15zoCZziPQbj+m5CRKoD9wOn4bqvuN+vI8DuwCRVTVTVRGBf9nNVHRdg+XcC//Ed4N11/RaQpqqjAq2Mqu5R1btwN3Vme9abnzEhi450AMYUNaq6UAIXCuuDqyMCrovsBcBd+HRZDywVkewu67vgdZcNICLZ3WW/JiKVgVhVzchLTCLSFNivqr/6DI7G9cG0UVVHeu3icF24ZEtV1XcDzHIRMFVEovXvLtyNyRdLIMbkXbAusoN1i51bd9ldcXsxeXU6rgdhX3fiEtTw7AGqug74l//EIjIe11fULao6VlWzRGQT0ApIy0ccxuSwQ1jGhMDb2ziWfoB6ELzeRiB1AP+9lcVAB2/vJFeqOlRV66nqWJ/Bv+D6bTImJJZAjMm7YF1kB+sWO7fustvhOuPMq324jjx9LQSGA3MktCqP5bz5GhMSSyDG5F2wLrKDdVkfsLts7zzFV6p6KB/L3gCc7D9QVd8CngA+FJGq+VyfpoS5K3ZTslkCMcaPiLwGLAGaiUi6iFznjXoU6CYiG3HnMB71hs8GNuO6xH4e+AeAd/L8/3C1KpYDD3rDegIf5jOshUDr7BomvlT1WVyholki4r+XEmwda+Ou+vopn3EYk8O6czemkHlXY12p+SysJSJjgf+q6rwCiOFW4HdVfeFY52VKL9sDMaaQqWq3/CYPzyPAcQUUxi7cpcjGhMz2QIwxxoTE9kCMMcaExBKIMcaYkFgCMcYYExJLIMYYY0JiCcQYY0xI/h9AZLcq33bOGQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scm.plams.tools.plot import linear_fit_extrapolate_to_0\n",
    "\n",
    "\n",
    "inverted_T_list = 1.0 / np.array(T_list)\n",
    "ln_D_list = np.log(DT_list)\n",
    "fit_x, fit_y, slope, intercept = linear_fit_extrapolate_to_0(inverted_T_list, ln_D_list)\n",
    "\n",
    "D0 = np.exp(intercept)\n",
    "boltzmann_eV_per_K = 8.617e-5\n",
    "barrier = -boltzmann_eV_per_K * slope\n",
    "title = f\"D0: {D0:.2e} m²/s. Barrier: {barrier:.3f} eV\"\n",
    "log(title)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 3))\n",
    "ax.plot(inverted_T_list * 1000, ln_D_list, \".\", label=\"Simulation results\")\n",
    "ax.plot(fit_x * 1000, fit_y, label=\"Linear fit\")\n",
    "ax.set_title(title)\n",
    "ax.legend()\n",
    "ax.set_xlabel(\"1000/T (K⁻¹)\")\n",
    "ax.set_ylabel(\"ln D\")\n",
    "ax;"
   ]
  }
 ],
 "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}