HEAVY28.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "9ccb80e9-e19e-44c6-aa0b-a70a32c38156",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymolpro\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "c2281a20-b900-4964-bb52-d328d378797b",
   "metadata": {},
   "outputs": [],
   "source": [
    "backend = 'local' \n",
    "project_name = 'HEAVY28_methodtest'\n",
    "parallel = 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "30498226-59eb-4577-9ff7-2da521cd1385",
   "metadata": {},
   "outputs": [],
   "source": [
    "methods = {\"HF\": \"HF;counterpoise\",\n",
    "\"MP2\": \"hf;MP2;counterpoise\",\n",
    "\"MP3\": \"hf;MP3;counterpoise\",\n",
    "\"MP4\": \"hf;MP4;counterpoise\",\n",
    "\"PAPT2\":\"hf;PAPT;MP2;counterpoise\",\n",
    "\"PAPT3\":\"hf;PAPT;MP3;counterpoise\",\n",
    "\"PAPT4\":\"hf;PAPT;MP4;counterpoise\",\n",
    "\"CCSD\":\"hf;CCSD;counterpoise\",\n",
    "\"CCSD(T)\":\"hf;CCSD(T);counterpoise\" }\n",
    "bases = ['aug-cc-pVDZ-pp'],['aug-cc-pVDZ-pp']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "72a3df68-9aed-43f6-ba7e-5feb8a68e39c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GMTKN55 HEAVY28 (closed shell only) (maximum number of electrons 70):\n",
      "\n",
      "Molecules:\n",
      "h2o: {'geometry': '           3\\n  \\n    O          0.00000000      0.00000000      0.39219533\\n    H         -0.75610000      0.00000000     -0.19609767\\n    H          0.75610000      0.00000000     -0.19609767\\n', 'description': 'h2o'}\n",
      "teh2_nh3: {'geometry': '           7\\n  \\n    TE        -0.27801386     -1.82151771      0.00000000\\n    H          0.21816614     -3.38866571      0.00000000\\n    H          1.27654414     -1.33346771      0.00000000\\n    N         -0.27801386      1.35373929      0.00000000\\n    H          0.15335914      1.77217829      0.81226000\\n    H          0.15335914      1.77217829     -0.81226000\\n    H         -1.24540086      1.64555529      0.00000000\\n', 'description': 'teh2_nh3'}\n",
      "sbh3: {'geometry': '           4\\n  \\n    SB         0.00000000      0.00000000      0.70962975\\n    H         -0.69522101      1.20415811     -0.23654319\\n    H         -0.69522101     -1.20415811     -0.23654319\\n    H          1.39044201      0.00000000     -0.23654319\\n', 'description': 'sbh3'}\n",
      "sbh3_h2o: {'geometry': '           7\\n  \\n    SB         0.07940814     -1.16740071      0.00000000\\n    H          0.92542714     -1.96223271      1.21137400\\n    H          0.92542714     -1.96223271     -1.21137400\\n    H         -1.13926886     -2.32204371      0.00000000\\n    O          0.07940814      2.66198329      0.00000000\\n    H         -0.06916886      1.71249229      0.00000000\\n    H         -0.80123286      3.03943429      0.00000000\\n', 'description': 'sbh3_h2o'}\n",
      "teh2: {'geometry': '           3\\n  \\n    TE         0.00000000      0.00000000      0.76954933\\n    H         -1.15536600      0.00000000     -0.38477467\\n    H          1.15536600      0.00000000     -0.38477467\\n', 'description': 'teh2'}\n",
      "sbh3_nh3: {'geometry': '           8\\n  \\n    SB         0.36106250     -1.68885413      0.00000000\\n    H         -0.68440450     -1.17992612      1.21267800\\n    H         -0.34318450     -3.23078712      0.00000000\\n    H         -0.68440450     -1.17992612     -1.21267800\\n    N          0.36106250      1.53738488      0.00000000\\n    H         -0.13718250      1.87593687      0.81122600\\n    H          1.26423350      1.99023488      0.00000000\\n    H         -0.13718250      1.87593687     -0.81122600\\n', 'description': 'sbh3_nh3'}\n",
      "nh3: {'geometry': '           4\\n  \\n    N         -0.00000025      0.28324650      0.00000000\\n    H          0.93623275     -0.09441550      0.00000000\\n    H         -0.46811625     -0.09441550      0.81080100\\n    H         -0.46811625     -0.09441550     -0.81080100\\n', 'description': 'nh3'}\n",
      "teh2_h2o: {'geometry': '           6\\n  \\n    TE         0.00000000      0.00000000      1.62022367\\n    H         -1.14729200      0.00000000      2.78676467\\n    H          1.14729200      0.00000000      2.78676467\\n    O          0.00000000      0.00000000     -2.00527833\\n    H         -0.75610200      0.00000000     -2.59423733\\n    H          0.75610200      0.00000000     -2.59423733\\n', 'description': 'teh2_h2o'}\n",
      "\n",
      "Reactions:\n",
      "16: sbh3_h2o -> sbh3 + h2o, energy = 0.0025019542608000003\n",
      "21: sbh3_nh3 -> sbh3 + nh3, energy = 0.0045258280896\n",
      "23: teh2_h2o -> teh2 + h2o, energy = 0.0010836489792000002\n",
      "28: teh2_nh3 -> teh2 + nh3, energy = 0.005338564824\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "db = pymolpro.database.load(\"GMTKN55_HEAVY28\")\n",
    "small=db.subset(open_shell=False,max_electrons=70)\n",
    "print(small)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "b1187ea2-8d36-482f-b3a6-909ecd7fa1e5",
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "unhashable type: 'list'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn [57], line 5\u001b[0m\n\u001b[1;32m      3\u001b[0m     results[method] \u001b[38;5;241m=\u001b[39m {}\n\u001b[1;32m      4\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m basis \u001b[38;5;129;01min\u001b[39;00m bases:\n\u001b[0;32m----> 5\u001b[0m         \u001b[43mresults\u001b[49m\u001b[43m[\u001b[49m\u001b[43mmethod\u001b[49m\u001b[43m]\u001b[49m\u001b[43m[\u001b[49m\u001b[43mbasis\u001b[49m\u001b[43m]\u001b[49m \u001b[38;5;241m=\u001b[39m pymolpro\u001b[38;5;241m.\u001b[39mdatabase\u001b[38;5;241m.\u001b[39mrun(small, methods[method], basis, location\u001b[38;5;241m=\u001b[39mproject_name,\n\u001b[1;32m      6\u001b[0m                                                        backend\u001b[38;5;241m=\u001b[39mbackend,\n\u001b[1;32m      7\u001b[0m                                                        preamble\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcore,small\u001b[39m\u001b[38;5;124m\"\u001b[39m, parallel\u001b[38;5;241m=\u001b[39mparallel)\n\u001b[1;32m      8\u001b[0m         \u001b[38;5;28;01mif\u001b[39;00m results[method][basis]\u001b[38;5;241m.\u001b[39mfailed: \u001b[38;5;28mprint\u001b[39m(method, basis, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mfailed\u001b[39m\u001b[38;5;124m'\u001b[39m, results[method][basis]\u001b[38;5;241m.\u001b[39mproject_directory)\n\u001b[1;32m     10\u001b[0m \u001b[38;5;28mprint\u001b[39m(results)\n",
      "\u001b[0;31mTypeError\u001b[0m: unhashable type: 'list'"
     ]
    }
   ],
   "source": [
    "results = {}\n",
    "for method in methods:\n",
    "    results[method] = {}\n",
    "    for basis in bases:\n",
    "        results[method][basis] = pymolpro.database.run(small, methods[method], basis, location=project_name,\n",
    "                                                       backend=backend,\n",
    "                                                       preamble=\"core,small\", parallel=parallel)\n",
    "        if results[method][basis].failed: print(method, basis, 'failed', results[method][basis].project_directory)\n",
    "        \n",
    "print(results)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cc27890-5b33-4627-8c6f-b80eae1be616",
   "metadata": {},
   "outputs": [],
   "source": [
    "analyse=pymolpro.database.analyse([result['aug-cc-pVDZ-pp'] for result in results.values()],small)\n",
    "analyse['reaction statistics']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f76f3df-f492-4dfa-8808-06b7c1a286fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas \n",
    "\n",
    "methods_pruned = [method for method in methods if (method != 'HF') and (method != 'CCSD(T)')]\n",
    "bases_pruned = ['aug-cc-pVDZ-pp','aug-cc-pVDZ-pp']\n",
    "fig, panes = plt.subplots(nrows=1, ncols=len(bases_pruned), sharey=True, figsize=(18, 6))\n",
    "\n",
    "for pane in range(len(bases_pruned)):\n",
    "    data = []\n",
    "    for method in methods_pruned:\n",
    "        data.append(\n",
    "            pymolpro.database.analyse(results[method][bases_pruned[pane]],\n",
    "                                      results['CCSD(T)']['aug-cc-pVDZ-pp'],'kJ/mol')['reaction energy deviations'].to_numpy()[:, 0]\n",
    "        )\n",
    "    panes[pane].violinplot(data, showmeans=True, showextrema=True, vert=True, bw_method='silverman')\n",
    "    panes[pane].set_xticks(range(1, len(methods_pruned) + 1), labels=methods_pruned, rotation=-90)\n",
    "    panes[pane].set_title(bases_pruned[pane])\n",
    "panes[0].set_ylabel('Error / kJ/mol')\n",
    "plt.savefig(project_name + \".violin.pdf\")\n",
    "\n",
    "df = pd.DataFrame(data)\n",
    "df.to_excel(\"PAPT.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "608037de-a623-41f9-98dd-1a5fb001aac6",
   "metadata": {},
   "outputs": [],
   "source": [
    "extrapolate=pymolpro.database.basis_extrapolate(results[method].values(), results[\"HF\"].values())\n",
    "extrapolate"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:persistent] *",
   "language": "python",
   "name": "conda-env-persistent-py"
  },
  "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}