{ "cells": [ { "cell_type": "markdown", "id": "17f297b4", "metadata": {}, "source": [ "## Initial Imports" ] }, { "cell_type": "code", "execution_count": 1, "id": "57c123de", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PLAMS working folder: /path/plams/examples/BasisSetBenchmark/plams_workdir\n" ] } ], "source": [ "import sys\n", "import multiprocessing\n", "from scm.plams import JobRunner, config, from_smiles, Settings, AMSJob, init\n", "import numpy as np\n", "\n", "# this line is not required in AMS2025+\n", "init();" ] }, { "cell_type": "markdown", "id": "99e4d81e", "metadata": {}, "source": [ "## Set Up Job Runner\n", "Set up job runner, running as many jobs as possible in parallel." ] }, { "cell_type": "code", "execution_count": 2, "id": "d9855c38", "metadata": {}, "outputs": [], "source": [ "config.default_jobrunner = JobRunner(parallel=True, maxjobs=multiprocessing.cpu_count())\n", "config.job.runscript.nproc = 1" ] }, { "cell_type": "markdown", "id": "7606c8bd", "metadata": {}, "source": [ "## Set Up Molecules\n", "Create the molecules we want to use in our benchmark from SMILES." ] }, { "cell_type": "code", "execution_count": 3, "id": "4c9dd47a", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Methane Atoms: \n", " 1 C 0.000000 0.000000 0.000000\n", " 2 H 0.538912 0.762358 -0.599295\n", " 3 H 0.731244 -0.596616 0.583182\n", " 4 H -0.567129 -0.670302 -0.678108\n", " 5 H -0.703028 0.504560 0.694220\n", " Bonds: \n", " (1)--1.0--(2)\n", " (1)--1.0--(3)\n", " (1)--1.0--(4)\n", " (1)--1.0--(5)\n", "\n", "Ethane Atoms: \n", " 1 C -0.757196 -0.040522 0.044605\n", " 2 C 0.757196 0.040522 -0.044605\n", " 3 H -1.205222 0.185290 -0.945970\n", " 4 H -1.130281 0.694397 0.788688\n", " 5 H -1.061719 -1.061491 0.357407\n", " 6 H 1.205222 -0.185290 0.945971\n", " 7 H 1.130281 -0.694396 -0.788689\n", " 8 H 1.061719 1.061491 -0.357406\n", " Bonds: \n", " (1)--1.0--(2)\n", " (1)--1.0--(3)\n", " (1)--1.0--(4)\n", " (1)--1.0--(5)\n", " (2)--1.0--(6)\n", " (2)--1.0--(7)\n", " (2)--1.0--(8)\n", "\n", "Ethylene Atoms: \n", " 1 C 0.664485 0.027988 -0.023685\n", " 2 C -0.664485 -0.027988 0.023685\n", " 3 H 1.253433 -0.878614 0.070299\n", " 4 H 1.167038 0.980564 -0.156575\n", " 5 H -1.253433 0.878614 -0.070299\n", " 6 H -1.167038 -0.980564 0.156575\n", " Bonds: \n", " (1)--2.0--(2)\n", " (1)--1.0--(3)\n", " (1)--1.0--(4)\n", " (2)--1.0--(5)\n", " (2)--1.0--(6)\n", "\n", "Acetylene Atoms: \n", " 1 C -0.587409 0.175060 -0.002211\n", " 2 C 0.587409 -0.094463 0.002211\n", " 3 H -1.618985 0.411721 -0.006095\n", " 4 H 1.618985 -0.331124 0.006094\n", " Bonds: \n", " (1)--3.0--(2)\n", " (1)--1.0--(3)\n", " (2)--1.0--(4)\n", "\n" ] } ], "source": [ "# The molecules we want to use in our benchmark:\n", "mol_smiles = {\"Methane\": \"C\", \"Ethane\": \"C-C\", \"Ethylene\": \"C=C\", \"Acetylene\": \"C#C\"}\n", "molecules = {}\n", "for name, smiles in mol_smiles.items():\n", " # Compute 10 conformers, optimize with UFF and pick the lowest in energy.\n", " molecules[name] = from_smiles(smiles, nconfs=10, forcefield=\"uff\")[0]\n", " print(name, molecules[name])" ] }, { "cell_type": "markdown", "id": "2164f3d9", "metadata": {}, "source": [ "## Initialize Calculation Settings\n", "Set up the settings which are common across jobs. The basis type is added later for each job." ] }, { "cell_type": "code", "execution_count": 4, "id": "8596ead9", "metadata": {}, "outputs": [], "source": [ "common_settings = Settings()\n", "common_settings.input.ams.Task = \"SinglePoint\"\n", "common_settings.input.ams.System.Symmetrize = \"Yes\"\n", "common_settings.input.adf.Basis.Core = \"None\"" ] }, { "cell_type": "code", "execution_count": 5, "id": "ea59fcef", "metadata": {}, "outputs": [], "source": [ "basis = [\"QZ4P\", \"TZ2P\", \"TZP\", \"DZP\", \"DZ\", \"SZ\"]\n", "reference_basis = \"QZ4P\"" ] }, { "cell_type": "markdown", "id": "83f31b12", "metadata": {}, "source": [ "## Run Calculations" ] }, { "cell_type": "code", "execution_count": 6, "id": "46585822", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[25.02|15:23:48] JOB Methane_QZ4P STARTED\n", "[25.02|15:23:48] JOB Ethane_QZ4P STARTED\n", "[25.02|15:23:48] JOB Ethylene_QZ4P STARTED\n", "[25.02|15:23:48] JOB Acetylene_QZ4P STARTED\n", "[25.02|15:23:48] JOB Methane_TZ2P STARTED\n", "[25.02|15:23:48] JOB Ethane_TZ2P STARTED\n", "[25.02|15:23:48] JOB Ethylene_TZ2P STARTED\n", "[25.02|15:23:48] JOB Acetylene_TZ2P STARTED\n", "[25.02|15:23:48] JOB Methane_TZP STARTED\n", "[25.02|15:23:48] JOB Ethane_TZP STARTED\n", "[25.02|15:23:48] JOB Ethylene_TZP STARTED\n", "[25.02|15:23:48] JOB Methane_QZ4P RUNNING\n", "[25.02|15:23:48] JOB Acetylene_TZP STARTED\n", "[25.02|15:23:48] JOB Methane_DZP STARTED\n", "[25.02|15:23:48] JOB Ethane_DZP STARTED\n", "[25.02|15:23:48] JOB Ethylene_DZP STARTED\n", "[25.02|15:23:48] JOB Acetylene_DZP STARTED\n", "[25.02|15:23:48] JOB Ethane_QZ4P RUNNING\n", "[25.02|15:23:48] JOB Methane_DZ STARTED\n", "[25.02|15:23:48] JOB Ethane_DZ STARTED\n", "[25.02|15:23:48] JOB Ethylene_DZ STARTED\n", "[25.02|15:23:48] JOB Acetylene_DZ STARTED\n", "[25.02|15:23:48] JOB Methane_SZ STARTED\n", "[25.02|15:23:48] JOB Acetylene_QZ4P RUNNING\n", "[25.02|15:23:48] JOB Ethane_SZ STARTED\n", "[25.02|15:23:48] JOB Ethylene_SZ STARTED\n", "[25.02|15:23:48] JOB Ethylene_QZ4P RUNNING\n", "[25.02|15:23:48] JOB Acetylene_SZ STARTED\n", "[25.02|15:23:48] JOB Methane_TZ2P RUNNING\n" ] } ], "source": [ "results = {}\n", "jobs = []\n", "for bas in basis:\n", " for name, molecule in molecules.items():\n", " settings = common_settings.copy()\n", " settings.input.adf.Basis.Type = bas\n", " job = AMSJob(name=name + \"_\" + bas, molecule=molecule, settings=settings)\n", " jobs.append(job)\n", " results[(name, bas)] = job.run()" ] }, { "cell_type": "markdown", "id": "5bbe7c66", "metadata": {}, "source": [ "## Results\n", "Extract the energy from each calculation. Calculate the average absolute error in bond energy per atom for each basis set." ] }, { "cell_type": "code", "execution_count": 7, "id": "de91b7a9-72d4-4abe-b95e-77d6174ee967", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[25.02|15:23:48] Waiting for job Methane_QZ4P to finish\n", "[25.02|15:23:48] JOB Ethylene_TZ2P RUNNING\n", "[25.02|15:23:48] JOB Acetylene_TZ2P RUNNING\n", "[25.02|15:23:48] JOB Ethane_TZ2P RUNNING\n", "[25.02|15:23:48] JOB Acetylene_TZP RUNNING\n", "[25.02|15:23:48] JOB Acetylene_DZ RUNNING\n", "[25.02|15:23:48] JOB Ethane_DZP RUNNING\n", "[25.02|15:23:48] JOB Ethylene_DZP RUNNING\n", "[25.02|15:23:48] JOB Methane_DZP RUNNING\n", "[25.02|15:23:48] JOB Methane_SZ RUNNING\n", "[25.02|15:23:48] JOB Ethane_TZP RUNNING\n", "[25.02|15:23:48] JOB Ethylene_TZP RUNNING\n", "[25.02|15:23:48] JOB Ethylene_SZ RUNNING\n", "[25.02|15:23:48] JOB Ethylene_DZ RUNNING\n", "[25.02|15:23:48] JOB Methane_DZ RUNNING\n", "[25.02|15:23:48] JOB Acetylene_DZP RUNNING\n", "[25.02|15:23:48] JOB Acetylene_SZ RUNNING\n", "[25.02|15:23:48] JOB Ethane_SZ RUNNING\n", "[25.02|15:23:48] JOB Ethane_DZ RUNNING\n", "[25.02|15:23:48] JOB Methane_TZP RUNNING\n", "[25.02|15:23:50] JOB Methane_QZ4P FINISHED\n", "[25.02|15:23:50] JOB Methane_QZ4P SUCCESSFUL\n", "[25.02|15:23:50] Waiting for job Ethane_QZ4P to finish\n", "[25.02|15:23:50] JOB Methane_TZ2P FINISHED\n", "[25.02|15:23:50] JOB Methane_TZ2P SUCCESSFUL\n", "[25.02|15:23:51] JOB Acetylene_QZ4P FINISHED\n", "[25.02|15:23:51] JOB Acetylene_QZ4P SUCCESSFUL\n", "[25.02|15:23:51] JOB Ethylene_TZ2P FINISHED\n", "[25.02|15:23:51] JOB Ethylene_TZ2P SUCCESSFUL\n", "[25.02|15:23:51] JOB Acetylene_TZP FINISHED\n", "[25.02|15:23:51] JOB Acetylene_TZP SUCCESSFUL\n", "[25.02|15:23:51] JOB Ethylene_QZ4P FINISHED\n", "[25.02|15:23:51] JOB Ethylene_QZ4P SUCCESSFUL\n", "[25.02|15:23:51] JOB Acetylene_DZ FINISHED\n", "[25.02|15:23:51] JOB Acetylene_DZ SUCCESSFUL\n", "[25.02|15:23:51] JOB Ethylene_DZP FINISHED\n", "[25.02|15:23:51] JOB Ethylene_DZP SUCCESSFUL\n", "[25.02|15:23:51] JOB Acetylene_TZ2P FINISHED\n", "[25.02|15:23:51] JOB Acetylene_TZ2P SUCCESSFUL\n", "[25.02|15:23:52] JOB Ethane_TZP FINISHED\n", "[25.02|15:23:52] JOB Ethane_TZP SUCCESSFUL\n", "[25.02|15:23:52] JOB Ethane_TZ2P FINISHED\n", "[25.02|15:23:52] JOB Ethane_TZ2P SUCCESSFUL\n", "[25.02|15:23:52] JOB Ethane_QZ4P FINISHED\n", "[25.02|15:23:52] JOB Ethane_QZ4P SUCCESSFUL\n", "[25.02|15:23:52] Waiting for job Methane_TZP to finish\n", "[25.02|15:23:52] JOB Ethylene_TZP FINISHED\n", "[25.02|15:23:52] JOB Ethylene_TZP SUCCESSFUL\n", "[25.02|15:23:52] JOB Ethylene_DZ FINISHED\n", "[25.02|15:23:52] JOB Ethylene_DZ SUCCESSFUL\n", "[25.02|15:23:52] JOB Acetylene_SZ FINISHED\n", "[25.02|15:23:52] JOB Acetylene_SZ SUCCESSFUL\n", "[25.02|15:23:53] JOB Methane_DZ FINISHED\n", "[25.02|15:23:53] JOB Methane_DZ SUCCESSFUL\n", "[25.02|15:23:53] JOB Methane_DZP FINISHED\n", "[25.02|15:23:53] JOB Methane_DZP SUCCESSFUL\n", "[25.02|15:23:53] JOB Methane_SZ FINISHED\n", "[25.02|15:23:53] JOB Methane_SZ SUCCESSFUL\n", "[25.02|15:23:53] JOB Ethylene_SZ FINISHED\n", "[25.02|15:23:53] JOB Ethylene_SZ SUCCESSFUL\n", "[25.02|15:23:53] JOB Ethane_DZ FINISHED\n", "[25.02|15:23:53] JOB Ethane_DZ SUCCESSFUL\n", "[25.02|15:23:53] JOB Acetylene_DZP FINISHED\n", "[25.02|15:23:53] JOB Acetylene_DZP SUCCESSFUL\n", "[25.02|15:23:54] JOB Methane_TZP FINISHED\n", "[25.02|15:23:54] JOB Methane_TZP SUCCESSFUL\n", "[25.02|15:23:54] Waiting for job Ethane_DZP to finish\n", "[25.02|15:23:54] JOB Ethane_SZ FINISHED\n", "[25.02|15:23:54] JOB Ethane_SZ SUCCESSFUL\n", "[25.02|15:23:54] JOB Ethane_DZP FINISHED\n", "[25.02|15:23:54] JOB Ethane_DZP SUCCESSFUL\n" ] }, { "data": { "text/markdown": [ "| Formula | Smiles | Basis | NAtoms | Energy [kcal/mol] | Average Error [kcal/mol] |\n", "|---------|--------|-------|--------|-------------------|--------------------------|\n", "| C2H2 | C#C | DZ | 4 | -537.10 | 4.91 |\n", "| C2H2 | C#C | DZP | 4 | -550.65 | 1.53 |\n", "| C2H2 | C#C | TZP | 4 | -552.96 | 0.95 |\n", "| C2H2 | C#C | TZ2P | 4 | -555.67 | 0.27 |\n", "| C2H2 | C#C | QZ4P | 4 | -556.76 | 0.00 |\n", "| C2H2 | C#C | SZ | 4 | -647.50 | 22.69 |\n", "| CH4 | C | DZ | 5 | -560.93 | 2.34 |\n", "| CH4 | C | DZP | 5 | -569.12 | 0.70 |\n", "| CH4 | C | TZP | 5 | -571.04 | 0.32 |\n", "| CH4 | C | TZ2P | 5 | -572.11 | 0.10 |\n", "| CH4 | C | QZ4P | 5 | -572.63 | 0.00 |\n", "| CH4 | C | SZ | 5 | -723.55 | 30.18 |\n", "| C2H4 | C=C | DZ | 6 | -750.17 | 3.37 |\n", "| C2H4 | C=C | DZP | 6 | -764.41 | 1.00 |\n", "| C2H4 | C=C | TZP | 6 | -767.33 | 0.51 |\n", "| C2H4 | C=C | TZ2P | 6 | -769.43 | 0.16 |\n", "| C2H4 | C=C | QZ4P | 6 | -770.41 | 0.00 |\n", "| C2H4 | C=C | SZ | 6 | -934.66 | 27.37 |\n", "| C2H6 | CC | SZ | 8 | -1216.91 | 30.49 |\n", "| C2H6 | CC | DZ | 8 | -951.17 | 2.73 |\n", "| C2H6 | CC | DZP | 8 | -966.09 | 0.87 |\n", "| C2H6 | CC | TZP | 8 | -970.08 | 0.37 |\n", "| C2H6 | CC | TZ2P | 8 | -971.88 | 0.14 |\n", "| C2H6 | CC | QZ4P | 8 | -973.02 | 0.00 |" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "try:\n", " # For AMS2025+ can use JobAnalysis class to perform results analysis\n", " from scm.plams import JobAnalysis\n", "\n", " ja = (\n", " JobAnalysis(jobs=jobs, standard_fields=[\"Formula\", \"Smiles\"])\n", " .add_settings_field((\"Input\", \"ADF\", \"Basis\", \"Type\"), display_name=\"Basis\")\n", " .add_field(\"NAtoms\", lambda j: len(j.molecule))\n", " .add_field(\n", " \"Energy\", lambda j: j.results.get_energy(unit=\"kcal/mol\"), display_name=\"Energy [kcal/mol]\", fmt=\".2f\"\n", " )\n", " .sort_jobs([\"NAtoms\", \"Energy\"])\n", " )\n", "\n", " ref_ja = ja.copy().filter_jobs(lambda data: data[\"InputAdfBasisType\"] == \"QZ4P\")\n", "\n", " ref_energies = {f: e for f, e in zip(ref_ja.Formula, ref_ja.Energy)}\n", "\n", " def get_average_error(job):\n", " return abs(job.results.get_energy(unit=\"kcal/mol\") - ref_energies[job.molecule.get_formula()]) / len(\n", " job.molecule\n", " )\n", "\n", " ja.add_field(\"AvErr\", get_average_error, display_name=\"Average Error [kcal/mol]\", fmt=\".2f\")\n", "\n", " # Pretty-print if running in a notebook\n", " if \"ipykernel\" in sys.modules:\n", " ja.display_table()\n", " else:\n", " print(ja.to_table())\n", "\n", "except ImportError:\n", "\n", " average_errors = {}\n", " for bas in basis:\n", " if bas != reference_basis:\n", " errors = []\n", " for name, molecule in molecules.items():\n", " reference_energy = results[(name, reference_basis)].get_energy(unit=\"kcal/mol\")\n", " energy = results[(name, bas)].get_energy(unit=\"kcal/mol\")\n", " errors.append(abs(energy - reference_energy) / len(molecule))\n", " print(\"Energy for {} using {} basis set: {} [kcal/mol]\".format(name, bas, energy))\n", " average_errors[bas] = sum(errors) / len(errors)" ] }, { "cell_type": "code", "execution_count": 8, "id": "ccac8e42-a5cc-4184-8d2e-21cba07f0b43", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "== Results ==\n", "Average absolute error in bond energy per atom\n", "Error for basis set TZ2P: 0.170 [kcal/mol]\n", "Error for basis set TZP : 0.537 [kcal/mol]\n", "Error for basis set DZP : 1.024 [kcal/mol]\n", "Error for basis set DZ : 3.339 [kcal/mol]\n", "Error for basis set SZ : 27.683 [kcal/mol]\n" ] } ], "source": [ "print(\"== Results ==\")\n", "print(\"Average absolute error in bond energy per atom\")\n", "for bas in basis:\n", " if bas != reference_basis:\n", " if ja:\n", " av = np.average(ja.copy().filter_jobs(lambda data: data[\"InputAdfBasisType\"] == bas).AvErr)\n", " else:\n", " av = average_errors[bas]\n", " print(\"Error for basis set {:<4}: {:>10.3f} [kcal/mol]\".format(bas, av))" ] }, { "cell_type": "code", "execution_count": null, "id": "a6ae4a50-3ded-4147-8f22-cb155586c232", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all", "executable": "/usr/bin/env plams", "main_language": "python", "notebook_metadata_filter": "-all" }, "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.8.16" } }, "nbformat": 4, "nbformat_minor": 5 }