Attempt to recreate Hammet constants from Peter Ertl's paper: https://chemistry-europe.onlinelibrary.wiley.com/doi/pdf/10.1002/cmtd.202200041

ertl_hammet_constants.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f99032df-1cc9-4440-9536-35063dd717be",
   "metadata": {},
   "source": [
    "Attempt to recreate Hammet constants from Peter Ertl's paper: https://chemistry-europe.onlinelibrary.wiley.com/doi/pdf/10.1002/cmtd.202200041"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "id": "4d47869b-442f-4904-aa97-00827a9170f8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from xtb.libxtb import VERBOSITY_MINIMAL, VERBOSITY_MUTED\n",
    "from xtb.interface import Calculator, Param\n",
    "import numpy as np\n",
    "\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem\n",
    "from rdkit.Chem import rdDistGeom\n",
    "\n",
    "mol = Chem.MolFromSmiles(\"c1ccccc1(C)\")\n",
    "mol = Chem.AddHs(mol, addCoords=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "id": "dfcd852f-900c-4a95-9826-dd71571165f6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for i, atom in enumerate(mol.GetAtoms()):\n",
    "    atom.SetProp(\"atomLabel\", f\"{i}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "id": "9a210dbd-3685-4f91-8aa4-5d10c8876d22",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAACWCAIAAADCEh9HAAAABmJLR0QA/wD/AP+gvaeTAAAejUlEQVR4nO3deVzM+R8H8M8cnZQp03ToIHTIsXRSRPhVjta1aR2LsuknMh2LZJOjRSXVuu/83KyfEIskZSVyFjXjalPSPR3TMc3x++PrN5LW1cx8RvN+/tX3U495vxyPd5/v9fmQRCIRAgAA8LXIuAMAAMC3DdooAAB0CLRRAADoEGijAADQIdBGAQCgQ6CNdlo1NTUfPoZRX1+fn59fVlaGJZK0lZaWfvhHrq2tbWhowJIHKAhoo53Q7du3vby86HQ6k8lsPR4ZGclgMCwtLfX09Dw8PGpra3EllCyhUHju3LmxY8fq6eklJiaKx1+8eLFkyRIDAwMnJyeM8UCnB220s9m6devkyZO1tLQMDQ1bT81OnDixcuXKiIiImpqaq1ev3rhxY+HChRhzSpC3t3d4eLijoyOZ/O7/8927d+3t7aurqx0cHDBmA4qABI/fdzJNTU0qKiokEsne3t7e3j4hIYEYHz58uEgkunHjBnEYFRW1YsWKN2/e0Ol0fGElg8fjKSsrI4SoVOqePXvmzp2LEBIKhUKhkEqlBgcHX7t27d69e5hTgs4LZqOdjaqqKolE+nD80aNH7u7u4kM3NzeBQJCbmyvDaNJC9NA2yGQylUqVfRiggKCNKoTa2tra2lp9fX3xCPF1cXExvlAAdBLQRhWCUChECLWepRKXEQUCAbZMAHQW0EYVAo1GU1dXLy8vF48QzzwZGhriCwVAJwFtVFFYWlpeu3ZNfJiWlkYikSwtLTFGAqBzgDaqKHx8fFJSUs6ePYsQevXqVXR09Lhx41pfLQUAfCUR6Fy4XK6WlpaWlhaVSlVRUdHS0goKChKJRAKBYN68eQghAwMDCoUyYMCA4uJi3GEl4/r1623+V+/atUskEg0ePLj1YP/+/XEnBZ0TPDfa2QiFwvv377ce0dHRMTY2Jr5msVh5eXl6enq2trYUCgVHQEkSCAQUCqW+vp7FYrUeNzExodPpeXl5rV8DVVNT69evn8wzgs4P2qji4nK5GzZsEIlE69atw53lKy1YsKCysjI2Nlb8ewIA2YM2qrju379vY2NDpVJzcnLMzMxwx/liDx48sLGxIZPJOTk55ubmuOMAxQW3mBTX4MGDfXx8eDxeUFAQ7ixfg8lkCgSCgIAA6KEAL5iNKrSysjIzM7OampoLFy60flVU/p04cWL69OkMBoPFYtFoNNxxgEKD2ahCYzAYK1euRAgFBQW1tLTgjvO5Ghsbly1bhhBat24d9FCAHbRRRUecFOfn52/duhV3ls8VHR1dUFDw3XffeXt7484CAJzUA4TOnz8/ceJEGo3GZrN1dHRwx/mE4uJic3NzLpeblpbm7OyMOw4AMBsFCE2YMMHd3Z3D4YSHh+PO8mlLly7lcrmenp7QQ4GcgNkoQAih/Pz8gQMHCoXCu3fvDho0CHecf5SZmeno6KiqqvrkyZOePXvijgMAQjAbBQQLC4uFCxcKBILAwEDcWf6RUChkMpkikSgkJAR6KJAfMBsFb1VXV5uZmVVUVJw6dWrq1Km447Rj//793t7ePXr0YLFYXbp0wR0HgLdgNgre0tLSWrNmDUIoKChIDnckrqurCwsLQwht3LgReiiQK9BGwTu+vr6DBg0qLCzcvHkz7ixtRUZGlpSUDB06dMaMGbizAPAeOKkH77l27ZqLi4u6unpeXp78rPfx4sULKyur5ubmW7du2dnZ4Y4DwHtgNgreM2rUqKlTpzY0NBBvN8mJ4ODgpqamOXPmQA8Fcghmo6Ctly9f9uvXr7m5OT093cnJCXcclJqaOnr06K5du7JYLAMDA9xxEEKovLw8MTGxb9++33///YffPXXqVHZ2trm5ObFONuj8MC0XDeTaihUrEELW1tYCgQBvEj6fP2DAAITQ+vXr8SYhvHz50tfXV01NjUwme3h4fPgDLBZLTU3NwMDA3d1d9vEAFnBSD9oRFhZmZGR09+7dQ4cO4U2yY8eOnJwcU1NTJpOJNwmhsrJSX1//r7/+mjJlyoffFYlECxcuHDduHLxhpVCgjYJ2qKurE0viHzx4EGMMoVC4fft2hFBMTIyqqirGJGLW1tYRERGDBw8mkUgffnf37t23b9+Oi4uTfTCAERV3ACCnZs+eTaFQPD09MWYgk8k3b948cuTI5MmTMcb4TCUlJcuXL1+3bp2hoSHuLECmoI2C9pFIpJkzZyKEysvLt23bdu/ePRqN5urq6uXlRSZL/STm8OHDSUlJenp6CQkJfn5+0i4nEf7+/qampv7+/riDAFmDNgo+prCw0NHR0crKauTIka9fv/bx8cnKyoqPj5d2UT8/PxMTkzZbnMqzU6dOnT17NisrqxPstwq+FDzwBD7m9u3bp0+f3rBhA3EYHh4eGxtbU1Mj1Wbh4eHB4/Hs7OyOHj369OlT6RXqCE9Pz+bm5qSkJOLQyMiotrZW/Fhrbm5uS0vLsGHDzp49iy8jkBGYjYKPsbOza/3EO4/HU1JSkupJ/aFDh1JSUnJycvDe3fpSERERAoFAfFhTU9PY2PhNXNIFHQdtFHxaZWUli8W6efPmtm3bYmNjSSTSpUuX0tLSJPX5VCp17dq1RKHg4ODw8PDevXtL6sMl7uDBg01NTS9evODz+bt27TIyMnJ3d/fx8Wn9M2lpaRwOBx6/VxDQRsGnZWRkMJnM0tLSIUOGjBo1CiF07dq1jRs3SurzlZSUiDYaEBCgo6MTHBwsqU+WhoMHD3I4HIQQlUrdtWvXsGHDPtxUtVevXnV1dTjSAQygjYJPmzRp0qRJk+rq6ubOnTtq1Cg2m+3m5tatWzdJfT5xpfXixYvHjh1LT09XUlKS1CdLQ0pKilAo5PP5FArln64RR0ZGyjgVwAhuMYEvkJ6e7uzsnJmZ6eDgIPEPt7e3z87OFnfnxsZGHo/XrVu3/Px8BoMh8XIdsXfv3vnz53t7e+/duxd3FoAfzEbBx1RWVmZlZY0bN444fPXqFUJIQ0NDGrXi4+Pr6+vFh4mJiampqYmJibATPZBz0EbBx6SkpMycOfPnn392cHAoLCyMi4tzc3OzsrKSRq02M9yMjIxbt26NGTNGGrUAkCBoo+Bjpk+f3qNHj+3bt2/btq1r167Lli0LCAiQTWlra2uhUCibWgB0BLRR0D6BQJCenj5q1CgnJycsq456eHh4eHhUVFTs378/JCSk3aVAAJAHsMITaN/evXtdXFwWLFiAMYNQKBw+fPjSpUuPHTuGMQYAHwdtFLSjtrZ21apVCCG8lybJZHJISAhCaOnSpVwuF2MSAD4C2ihoR0RExJs3bxwdHadNm0aMXL9+fefOna3fd5QSNps9ceLErKws4nDevHm2trZFRUXR0dHSLg3A14E2Ctp69uzZ1q1byWRyfHw8cUWSz+cvWrTIz89vz5490q6emJh4/vz5JUuWEE80i2NERUX9/fff0q4OwFeANgraWrJkCY/Hmz9/vrW1NTGyY8eO3NxcU1PTOXPmSLt6aGiogYFBVlbWf/7zH2Jk6NCh06dPb2xsXL58ubSrA/A18G4FBeTN5cuXEUKampolJSXESFVVFZ1ORwj997//lU2GxMREhJCurm5NTQ0x8urVqy5duiCErl+/LpsMH0fMyr29vXEHAXIBZqPgHT6fHxgYiBBatWqVnp4eMRgeHl5RUeHi4jJp0iTZxJg9e7a9vX1paWlUVBQxYmhoSNxrYjKZMrg+C8AXgTYK3klISHj8+HGfPn3EO2Hk5eXt3LmTQqHIcps2EokUFxdHIpFiYmKePXtGDC5btoxYD//AgQMySwLA54A2Ct4qLy8nVquLj49XUVEhBgMDA1taWvz8/IjN4mXGwcFh1qxZzc3N4uuhampqxCL8K1asqKmpkWUYAD4O2ih469dff+VwOGPGjBEvRJKUlHTp0iUtLa2IiAjZ54mKitLU1Pzjjz+uXLlCjHh5eY0YMaKsrAyWoQNyBdooQAihhw8f7tmzh0qlik/eeTze0qVLEUKrV68mbjHJmJ6eHhEgMDCQz+cTg3FxccQjUGw2W/aRAGgXtFGA0P9v3SxevFi8elNcXBybzba0tMS4v3FISEifPn0eP34sfl518ODB8+bN4/F4xB0nAOQC7kcFAH4nTpxACGlra1dWVhIjpaWlxPLJf/75J95sJ0+eJLJVVFS0yXbx4kVcqeCBJ9AazEYVXVNTE3Hu/Ntvv2lraxODoaGhNTU1Hh4erq6uWNOhadOmjR07tqqqat26dcQIg8EICwtDCAUFBbW0tGBNBwBCcFIPYmJiCgoKrKysxHtbEg8VKSsry8lr7Js3b6ZSqVu2bMnNzSVGlixZYmZmlpeXt337drzZAEDQRhVccXExscHnli1bqFQqQkgkEi1ZskQoFDKZTDMzM9wBEULIyspq/vz54lcDEELiFh8REVFRUYE1HQDQRhXb8uXL6+vrp02bNnLkSGLk6NGjGRkZDAZjxYoVWKO9JzIyUltbOyUlJTk5mRjx8PBwc3Orrq4mFvQDACNoo4rr1q1bhw8fVlVVFb9z2djYSHTP9evXS3D/5I7T1tZeuXIlQojJZDY3NxODsbGxSkpKO3fufPToEdZ0QNFRsDxZDeQBlUqtqqpydXWdMmUKMRIZGZmUlDR48OBt27bJ26Ydtra2p0+fZrFYNBpt2LBhCCEdHZ3y8nJjY+Pvv/9exk1fKBR279595MiRlpaWsqwL5BPsU4+ZSCR69OgRh8OxsrKSwVPuRUVFL1++1NfX7927d5tGWVRUZGFhweVyr1+/PmLECGkn+QopKSljx47V0NBgs9nEyilCoZBMJjc1NT1+/Lj1Tw4YMEBZWVl6SUQiUX5+/ps3b4yMjPr06SO9QoTa2trKyspevXq1HiwtLWWz2To6On379qVQKNLOAD4G7/NWCo7FYg0YMIBCoaiqqlKp1MjISOnVqqmpmTp1KoVCodPpZDJ52LBh4qXwCAsXLkQIeXl5SS9DxxEvqjKZzNaDN2/eRAjRaDSt/6urq5Nehlu3bpmbmyOEunfvTqVSw8LCpFfr9evXq1at0tLSQgi1tLQQg01NTb6+vlQqlU6nU6nU/v375+fnSy8D+CRoo9gIhcJBgwb17dv31atXLS0tCQkJCKGUlBQplQsKCjI3N2exWCKR6NmzZwwGY/bs2a1/gMvlrl69+u+//5ZSAInIz88PCQnhcDitB0+ePKmsrCwUCmUQoLi4mEajubu7v379WiQScTic8vJyKdWqr683MjKaNWsWsbGguI1u3rxZT0/vzp07IpHozZs35ubmI0eOlFIG8DmgjWKTn5+PEEpMTBSPWFpatmltEtTc3FxcXCw+XLRokYmJiZRqyVhCQoLM/izLly/X1NQUv+4lbcTvhiNHjrRuowKBoLCwUPwzGzduVFZWFn8XyB7sU49NQ0MDQkhXV1c8Mnjw4OfPnyOEqqurhUKhpAqRSCRtbW1lZWUDAwPxoEAg6Nq1q6RK4FVSUsJgMDIzM4uKikxMTGxtbYnLvpWVlRKsQqFQaDRaamrqqFGjNDQ0Hj58WFtba2Njo6amhhDicrlNTU0SLKeiotK1a9d2b/SRyWQjIyPxoUAgUFNTg8ujOOHu44qrublZV1fXzc2trq6Oz+dfvXq1V69ePXv2FIlE4pXnJYJEIn1Y2tDQMCAgAMefW/KCg4NVVVUdHBzs7OyoVKqLi4vEmxpCyMzMTCQS6erqDhs2rGfPnlpaWioqKnQ6ndjXZO7cuZIt1/q8pM1stDWhUGhjYzN58mSZ/W2DD8FsFBtlZeUjR454eXlpa2tTKBRTU1MGg1FdXY0Q0tbW5vF4kipEJrd9OnjZsmVcLpd4ErMTiImJWbt2LTErzM7OHjFiRHR0dGhoqHiJAIkg7vPw+XwOh3Ps2DF7e/v6+vpJkybNnj37xYsXXbt2lWw5DQ2Nz/mxzZs35+TkHDp0SIKlwRfD3ccVXWNj482bN1kslkAgcHZ2dnV1lXbF1atXq6urp6enS7sQLi4uLmPGjJHSh/fv33/BggXiQ2IBqmfPnkmpHOGfZqN79uxRVlY+deqUVKuDT4K3mDBTVVUdOnSomZlZSUnJrVu3xowZI71aXC535syZCQkJycnJw4cPl14hGWtzHbmqqopGo0mpVp8+fe7evSs+LCkpQQjJ/iozscJAQEDA4cOHp06dKuPqoC3cfVyhFRcXV1VViUQiFovl5ORkYGDQ5lEeCWKz2f379ye6QNX/CQQCKZWTmcbGRltb2/Dw8JcvX3I4nNWrV5NIpPPnz0up3B9//IEQ2r59O4/Hy8/PNzU1dXFxkVItkUjE5XKrqqp2796NECorK6uqqhIKhSUlJcOHD9fT00tNTRX/U8KdeoygjeJELFlE7MBuZ2f35MkT6dWysbH58JeonD8l+pn27dsnvnOtq6u7b98+qZb79ddflZSUiAWxHB0di4qKpFfL19e3zT9ZQ0OD+OXd1jIyMqQXA3wcvAyKWUFBQWlpqb6+vrGxsVQLFRUVfXjbytjYmM/nq6qqSrW09PD5fPH6fgUFBXw+39TUVAaP/lRXV7PZbDqd3rt3b6kWqqioqK2tbT3Sq1ev0tJS4mm51gwMDL7df8dvHbRRxZWXl+fn52dubr5r1y7cWb7SjBkzEEKbNm3S19eXZd2srKxjx445ODhMnz5dlnWBfIJbTIpLSUkpKytr79692dnZuLN8jZs3bx47duzMmTMSfDjsM+Xm5sbFxV2+fFnGdYF8gjaquPr06ePv708sdP/NnZSIYy9btszExAR3HKDQoI0qtFWrVunp6f3111/E84/fkH379t25c8fQ0BB2WgbYQRtVaJqamqtXr0YIhYSEfHjXQm7V1dWFh4cjhKKjo4nnHADACNqoops/f76Njc2rV682bdqEO8vnWrNmTUlJybBhw+AOD5AH0EYVHZlMjouLI5FI69evLywsxB3n054/f/7777+LY+OOAwC0UYCQo6PjtGnTxPvZyTliVztvb29bW1vcWQBACNooIGzatEldXf3IkSMZGRm4s3zM1atXz58/r6GhsWbNGtxZAHgL2ihACCEjI6Pg4GCRSMRkMiW4YrRk8fl8JpOJEAoPD5fx8/YAfAS0UfBWaGiosbHxvXv3EhMTcWdp39atW3Nzc3v37r148WLcWQB4B9ooeEtNTe23335DCIWGhrZ5j1seVFVVrV27FiG0efNmFRUV3HEAeAfaKHhnxowZTk5OpaWl69evx52lrZUrV1ZWVo4ePXrixIm4swDwHmij4B0SiRQfH08mk2NjY58+fYo7zjuPHz/evXs3lUqNi4vDnQWAtqCNgvcMGTLkp59+4vF4S5cuxZ3lncDAQD6f7+/v379/f9xZAGgL2ihoa8OGDZqammfOnLl06RLuLAghdPr06StXrmhra//666+4swDQDmijoC1dXd3Q0FCEUFBQUEtLC94wPB5v+fLlCKG1a9d2794dbxgA2gVtFLQjKCiob9++T548wb6ic0xMzNOnT/v16/fhdhoAyAnYpx60Q1lZOSoqavLkyefOnfP39//wB2pqaiT7lH63bt3I5HZ+qaempiKEEhISiM1CAJBD8F8TtG/SpEnnz593d3dv97v9+/cvKiqSYLmcnJx2bx9duXIlNTV19OjREqwFgGRBGwX/aPz48cQXDx8+jIuLI5PJe/fuJUZoNFp9fb0EaxH70LW0tCQkJJw7d665udnOzi4sLIzBYEAPBXIO2ij4mNevX8+ZMycjI0NHR6f1RiM5OTnSKOfn53f48OHg4GA6nR4fH3/58uW7d++qq6tLoxYAkgK3mMDH0Ol0X1/f8vLydq+QStaTJ0/2798fGxsbGRkZGBj4559/stnsffv2SbsuAB0EbRR8jLKy8g8//KChoSGDWlevXqVQKLNmzSIOLSws7O3tU1JSZFD6S6mrq+vq6nbr1g13ECAX4KQefI2QkJCqqioJfuDatWuLiooYDIampqZ4sG/fvo8fP5ZgFUn58ccff/zxR9wpgLyANgq+xvHjxyV7pz4oKKi6urp1D0UIaWpqSrZZSwSPxzt69Oj27dtHjBgRFRVFDAqFwvj4+P3795eUlJiYmCxatGju3LlYYwLZgTYKvsamTZske6fe0NCQwWC0aZqVlZV6enoSrNJx6enpXl5eXbp04fF4LBZLPB4YGJicnLx27VojI6OLFy/OmzdPXV3d09MTY1QgM9BGwdeQRoMwNjauqKgoKytjMBjESF5enoWFhcQLdUS/fv2uXLliZWXl6enZ3NwsHh83bty///1vIq2Tk1NSUtLFixehjSoIuMUE5IWbmxuFQtm5cydxmJmZ+eDBA3lbXZROp1tZWX047urqKu74HA6nuLi4Z8+eMk0G8IHZKPiE+Pj4J0+ePHjwgMPhLFiwQENDIyYmRhqFjI2Nly9fHhERwWaz6XT6wYMHnZ2dv62d6JOTk+/cuXP8+HEnJ6egoCDccYCMwGwUfIKGhoaWltaoUaMWLVqkpaVFo9GkV2vNmjWnTp2iUCilpaXh4eGXLl0i3m76VrDZ7EePHhELDjQ2NuKOA2SE1PrVFADEeDzeqlWrmEymrq4u3hiHDh2aM2eOvPVT4tpoUlLSh9/icrlDhgyxtrY+cuSI7IMB2YPZKGjf77//vmHDBg8PD7wxJkyY4OPj8229y9SlS5eRI0fevn0bdxAgI9BGQTvKysqIbTgjIiLwJpk/fz5CKCwsjMPh4E3ycdHR0WlpacTXLS0tt27d6tWrF9ZEQHagjYJ2hIWF1dTUTJgw4Z8WypMZT09PZ2fn8vJyoq3Lg2nTpo0dOzY9PT0zM3Ps2LErVqxACFVWVrq6uk6ZMmX+/PkDBw4sKCiIjIzEnRTICFwbBW3dv3/f1taWQqE8evTI3Nwcdxz04MEDGxsbMpmck5MjD3kSEhJa3z4yNTX94YcfEELZ2dnJycmVlZU9e/b08vIyMDDAlxHIFLRR0Jazs3N6enpISEh0dDTuLG/5+vru3r17/Pjx58+fx50FgLagjYL3HD9+3MvLi8FgsNls+VnBqKyszNzcnMPhXLhwAft1BgDagGuj4J3GxkZiG85169bJTw9FCDEYjJUrVyL52KwUgDagjYJ3oqKiCgoKvvvuO29vb9xZ2goICDA3N8/Pz9+yZQvuLAC8B07qwVtFRUUWFhZcLjctLc3Z2Rl3nHYkJydPmDBBU1OTxWLJ28pPQJHBbBS8tXTpUi6XO336dPnsoQih8ePHu7u719bWYn+aFYDWYDYKEEIoMzPT0dFRVVX1yZMn8rw0UX5+/sCBAwUCwe3bt62trXHHAQAhmI0ChJBQKGQymSKR6JdffpHnHooQsrCw8Pf3FwfGHQcAhGA2ChBC+/bt8/Hx6dGjB4vF6tKlC+44n8DhcMzMzMrLy0+ePDlt2jTccQCA2ajCq6urI54lioqKkv8eihCi0Whr1qxBCAUFBTU0NOCOAwC0UYW3bt26kpKSoUOHfkNbXf7888+DBg169epVbGws7iwAwEm9Ynv+/LmVlVVLS0tmZqadnR3uOF/gxo0bI0aMUFNTy8vLMzY2xh0HKDSYjSq0oKCg5ubmuXPnfls9FCHk5OQ0ZcqUhoYG4ooEABjBbFRxXb16dcyYMRoaGiwWS19fH3ecL1ZYWGhpadnY2Jienu7k5IQ7DlBcsKWd4rK3tw8NDdXR0fkWeyhCyNjY+Jdffqmurm53q04AZAZmo4rl7t27BQUF2traDg4OampquONIQENDQ3JycptBfX19mJ8CmYE2qigaGhomTJiQnZ09ZMiQ58+fi0Si5OTkQYMG4c7VUcXFxSNGjGg98vr1a1dX1zNnzuCKBBQNnNQril27dt28eTM/P79nz54tLS3Dhw8PDAxMTU3FnaujevTo8fz5c/FhdXW1sbGxp6cnxkhA0cCdekVRUFBAo9FMTEwQQkpKSv369SssLMQdSvJ27NjRrVs3YlcPAGQD2qiicHFxKS0t9fX1LS0tLSwsvHDhwk8//YQ7lIS1tLRs27ZtyZIlSkpKuLMABQLXRhXI4sWLt2zZoqSkpKmpOXLkyJMnT5JIJNyhJOnAgQMBAQGFhYU0Gg13FqBAYDaqKJKSkvbt23f8+PEzZ864uLicPn160aJFuENJWFxcnI+PD/RQIGMwG1UU5ubmw4cP37NnD3G4bds2f3//vLw8CwsLvMEk5dKlS+PHj2ez2aamprizAMUCs1FFweFwtLW1xYfEo07l5eX4EknYpk2bpk6dCj0UyB60UUXh5uZ24MCB7OxshFB1dXVkZKSurm6nWUA+Nzc3JSUlMDAQdxCgiKCNKoqEhIR//etfDg4OGhoadDq9pKTkwoUL6urquHNJRkxMjKOjo4ODA+4gQBHBtVHF0tTUVFhYqKGh8Y2+Rw+AHII2CgAAHQIn9QAA0CHQRgEAoEOgjQIAQIdAGwUAgA6BNgoAAB3yP7gVX1qX/TN0AAAA6HpUWHRyZGtpdFBLTCByZGtpdCAyMDIyLjA5LjQAAHice79v7T0GIOBlgABGIOaH4gZGtoQMIM3MQiRDA8xgU0gA0izMjAwaYBPJpbkZGDOYGJkSmJgzmJhZElhYM5hY2RhYGRIY2BkYORiYOBmYuRhYuBnYeBjYeBnY+BhEGIHSLMxMjOL7oD4BA/6zZ3xsH7ot2w/izJq5cz+MffYMz4EX2dxqSOL2EHGwehjb/uqf20tBbGNjYYfTmZ+XwsQfuqkdgJmDYOcceKU0VQ1J3AFJvQPEnGKH89evLYGZycBwYD8SG2yvGAC9T0tVWDx9sQAAATh6VFh0TU9MIHJka2l0IDIwMjIuMDkuNAAAeJyNk0tOxDAMhvc9hS8wkZ1H2yyn7QgQmlaCwh3Yc39hZ+jYFZCSh5Q4X+3kt9uAtJfp+eMT7s1PTQOAlZFzhveAiM0VZAHD5eFphnE9D5tlXN7m9RUoyUDpe/a8LtfNQjACOkzsEOFEzueMoRdTafqpZ5DP4zEZhPTOf5vv5z/AuHNZ8ZjsJStcW7j+Zj1VIncMeheOwR4ejUfvUu6xa38BM4P66BpJKGhw3aE+RDunNZ/e3rMGBgbDfySiyKRqRC79SSZLVsDLPO1q71aNwzJPWo2Rp9eaizyDFlaUqeUjPWmR8AZarQXi2WnGo5g0r8S91+xxVMgmRcSRyGaiXI2M4MTRyBthiYNQMPoVQzQyFUOyctjHy377jXndfAGz+MLkLeDLOgAAAMh6VFh0U01JTEVTIHJka2l0IDIwMjIuMDkuNAAAeJxdjbEOgzAMRH+lAwNIiWU7SZvgsQvduiOGigmpqKhl5OMLoaWki8/vTj7XVdNSm9dVU3znOWocO5vmeZhybeCkUIlm4KgE1ikCDgHtx2ZwwaNXgoC7DMHvI7+cM5hFzI/mQrducYkVOunQW8n6XKfft7hQWRAWI4RihUicHIVYaHasoJyExGfqNj766/MxlAjd69IP967tRqDSpOhS9CmGFOmvi/iPbcLTG2sfY8aCXFp4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<rdkit.Chem.rdchem.Mol at 0x12f5b4ac0>"
      ]
     },
     "execution_count": 179,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5177659b-adc2-414b-9862-8284d905b9cb",
   "metadata": {},
   "source": [
    "Atom positions that map to Peter's ring positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "id": "8af05ab7-9144-4d1e-87b1-4d3728bf4ad8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "class ChargeMap:\n",
    "    ATTACHMENT = 5\n",
    "    ORTHO = 0\n",
    "    META = 1\n",
    "    PARA = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "id": "07329c7f-cab8-4fbe-a9b2-597cf62b30a6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def atoms_and_pos_from_rdkit_mol(mol: Chem.Mol):\n",
    "    atoms = np.zeros(mol.GetNumAtoms())\n",
    "    pos = np.zeros((mol.GetNumAtoms(), 3))\n",
    "    \n",
    "    for i, atom in enumerate(mol.GetAtoms()):\n",
    "        positions = mol.GetConformer().GetAtomPosition(i)\n",
    "        atoms[i] = atom.GetAtomicNum()\n",
    "        pos[i] = (positions.x, positions.y, positions.z)\n",
    "    return atoms, pos\n",
    "\n",
    "def calculate_conformers(mol: Chem.Mol):\n",
    "    etkdg = rdDistGeom.ETKDGv3()\n",
    "    etkdg.randomSeed = 42\n",
    "    etkdg.verbose = False\n",
    "    etkdg.numThreads = 8\n",
    "    etkdg.useRandomCoords = True\n",
    "\n",
    "    num_conformers = 100\n",
    "\n",
    "    Chem.AssignStereochemistryFrom3D(mol)\n",
    "    conformation_ids = rdDistGeom.EmbedMultipleConfs(mol, numConfs=num_conformers, params=etkdg)\n",
    "    \n",
    "    return mol\n",
    "\n",
    "def calculate_charges(param, mol: Chem.Mol):\n",
    "    mol = calculate_conformers(mol)\n",
    "    \n",
    "    atoms, pos = atoms_and_pos_from_rdkit_mol(mol)\n",
    "    calc = Calculator(param, atoms, pos)\n",
    "    calc.set_verbosity(VERBOSITY_MUTED)\n",
    "    res = calc.singlepoint()\n",
    "    partial_charges = res.get_charges()\n",
    "    \n",
    "    \n",
    "    return partial_charges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "id": "d40a729e-51d9-4695-88d9-1ea321cef9db",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def hammet_constants_func(q_attachment: float, q_ortho: float, q_meta: float, q_para: float):\n",
    "    sigma_meta = 1.9381 + 0.9273 * q_attachment + 43.9107 * q_meta + 8.6774 * q_para \n",
    "    sigma_para = 2.0849 + 0.2074 * q_attachment + 28.4679 * q_meta + 28.9006 * q_para\n",
    "    \n",
    "    \n",
    "    return sigma_meta, sigma_para\n",
    "\n",
    "def calculate_hammet_constants(param, mol: Chem.Mol):\n",
    "    partial_charges = calculate_charges(param, mol)\n",
    "    \n",
    "    q_attachment=partial_charges[ChargeMap.ATTACHMENT]\n",
    "    q_ortho=partial_charges[ChargeMap.ORTHO]\n",
    "    q_meta=partial_charges[ChargeMap.META]\n",
    "    q_para=partial_charges[ChargeMap.PARA]\n",
    "    \n",
    "    print(f\"Charges:\\n{q_attachment =  :.3f}\\n{q_ortho =  :.3f}\\n{q_meta =  :.3f}\\n{q_para =  :.3f}\")\n",
    "    sigma_meta, sigma_para = hammet_constants(\n",
    "        q_attachment=q_attachment, \n",
    "        q_ortho=q_ortho,\n",
    "        q_meta=q_meta,\n",
    "        q_para=q_para\n",
    "    )\n",
    "    return (sigma_meta, sigma_para)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "id": "25b6e3f0-8980-4f0a-919b-d396712cd13e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Charges:\n",
      "q_attachment =  0.386\n",
      "q_ortho =  0.212\n",
      "q_meta =  0.209\n",
      "q_para =  0.071\n",
      "GFN1xTB: (12.07879763040109, 10.154674008151897)\n"
     ]
    }
   ],
   "source": [
    "print(f\"GFN1xTB: {calculate_hammet_constants(Param.GFN1xTB, mol)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "id": "5bdfb6f3-84e1-44ed-864b-0aca08edea41",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Charges:\n",
      "q_attachment =  0.240\n",
      "q_ortho =  -0.153\n",
      "q_meta =  -0.060\n",
      "q_para =  -0.146\n",
      "GFN2xTB: (-1.7473574765373785, -3.8024574359170535)\n"
     ]
    }
   ],
   "source": [
    "print(f\"GFN2xTB: {calculate_hammet_constants(Param.GFN2xTB, mol)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c340b834-6fd9-444d-90a9-d0044ab77ca9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}