diff --git a/data/naphthalene-azulene/azulene/azulene.yaml b/data/naphthalene-azulene/azulene/azulene.yaml new file mode 100644 index 0000000..7f4f9c6 --- /dev/null +++ b/data/naphthalene-azulene/azulene/azulene.yaml @@ -0,0 +1,24 @@ +coords: [ + [2.28106300, -0.11553900, 0.00000000], + [1.38735500, 0.95238500, 0.00000000], + [-0.88728900, -0.27262700, 0.00000000], + [0.00000000, 0.92800200, 0.00000000], + [-0.84195500, 2.04842700, 0.00000000], + [-2.16952200, 1.60305000, 0.00000000], + [-0.50265900, -1.60631100, 0.00000000], + [0.78051900, -2.14674200, 0.00000000], + [-2.20504100, 0.20332800, 0.00000000], + [2.00416100, -1.48072100, 0.00000000], + [3.33297500, 0.14999900, 0.00000000], + [1.83292400, 1.94476400, 0.00000000], + [-0.50854400, 3.07638000, 0.00000000], + [-3.03984900, 2.24584600, 0.00000000], + [-1.32037900, -2.32370600, 0.00000000], + [0.83535300, -3.23010800, 0.00000000], + [-3.09025100, -0.41657700, 0.00000000], + [2.87798200, -2.12610800, 0.00000000], +] +charges: [6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1] +charge: 0 +spin: 0 +unit: "angstrom" diff --git a/data/naphthalene-azulene/naphthalene/naphthalene.yaml b/data/naphthalene-azulene/naphthalene/naphthalene.yaml new file mode 100644 index 0000000..244770a --- /dev/null +++ b/data/naphthalene-azulene/naphthalene/naphthalene.yaml @@ -0,0 +1,24 @@ +coords: [ + [-0.70764600, -2.42410900, 0.00000000], + [-1.39804900, -1.24152200 , 0.00000000], + [1.39805100 , 1.24152600 , 0.00000000], + [-0.71093300, -0.00000200, 0.00000000], + [-1.39805300, 1.24152400, 0.00000000], + [-0.70765300, 2.42410600, 0.00000000], + [0.71093300 , -0.00000200 , 0.00000000], + [1.39805100 , -1.24152100 , 0.00000000], + [0.70764900 , 2.42410600 , 0.00000000], + [0.70764900 , -2.42410800 , 0.00000000], + [-1.24307200, -3.36595900, 0.00000000], + [-2.48281300, -1.23717600, 0.00000000], + [2.48281500 , 1.23717300 , 0.00000000], + [-2.48281700, 1.23717000, 0.00000000], + [-1.24307300, 3.36596000, 0.00000000], + [2.48281500 , -1.23717300 , 0.00000000], + [1.24306800 , 3.36596100 , 0.00000000], + [1.24307700 , -3.36595700 , 0.00000000] +] +charges: [6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1] +charge: 0 +spin: 0 +unit: "angstrom" diff --git a/experiment_results/09_azulene_naphthalene/finetune_params.npz b/experiment_results/09_azulene_naphthalene/finetune_params.npz new file mode 100644 index 0000000..54852e7 Binary files /dev/null and b/experiment_results/09_azulene_naphthalene/finetune_params.npz differ diff --git a/experiment_results/09_azulene_naphthalene/molecule_images/azulene-structure.png b/experiment_results/09_azulene_naphthalene/molecule_images/azulene-structure.png new file mode 100755 index 0000000..f537570 Binary files /dev/null and b/experiment_results/09_azulene_naphthalene/molecule_images/azulene-structure.png differ diff --git a/experiment_results/09_azulene_naphthalene/molecule_images/naphthalene-structure.png b/experiment_results/09_azulene_naphthalene/molecule_images/naphthalene-structure.png new file mode 100755 index 0000000..2389011 Binary files /dev/null and b/experiment_results/09_azulene_naphthalene/molecule_images/naphthalene-structure.png differ diff --git a/notebooks/09_azulene_naphthalene/create_figures.ipynb b/notebooks/09_azulene_naphthalene/create_figures.ipynb new file mode 100644 index 0000000..3360693 --- /dev/null +++ b/notebooks/09_azulene_naphthalene/create_figures.ipynb @@ -0,0 +1,410 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "e76c60e1-c1e0-4db5-af6d-a31d92111f07", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pyscf\n", + "import py3Dmol\n", + "from uncertainties import ufloat\n", + "import matplotlib.pyplot as plt\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39b665b4-e3f7-401c-8bfe-123faad9009a", + "metadata": {}, + "outputs": [], + "source": [ + "from oneqmc.analysis.visual import show_mol\n", + "from oneqmc.convert_geo import load_molecules\n", + "from oneqmc.analysis.plot import set_defaults\n", + "\n", + "set_defaults()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a047adae-a23f-4ecf-9428-eaf8759899ef", + "metadata": {}, + "outputs": [], + "source": [ + "with np.load(\n", + " \"../../../data_paper/09_azulene_naphthalene/finetune_params.npz\"\n", + ") as params_npz:\n", + " params = {key: params_npz[key] for key in params_npz.files}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7812b505-fb06-44e4-89bd-859f323bbd5b", + "metadata": {}, + "outputs": [], + "source": [ + "naphthalene = load_molecules(\n", + " \"../../../data/naphthalene-azulene/naphthalene\"\n", + ")[0]\n", + "azulene = load_molecules(\n", + " \"../../../data/naphthalene-azulene/azulene\"\n", + ")[0]" + ] + }, + { + "cell_type": "markdown", + "id": "6be2ff5d-d9da-4f7c-892e-ba496687dc05", + "metadata": {}, + "source": [ + "#### Isosurfaces of envelopes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c98c10d5-5f60-48bd-adee-855e94c1d994", + "metadata": {}, + "outputs": [], + "source": [ + "class CubeDataFormatter:\n", + " def __init__(self, mol, nx=50, ny=50, nz=50):\n", + " self.charges = mol.charges\n", + " self.coords = mol.coords\n", + " margin = 3.0\n", + " extent = np.max(mol.coords, axis=0) - np.min(mol.coords, axis=0) + 2 * margin\n", + " self.box = np.diag(extent)\n", + " self.boxorig = np.min(mol.coords, axis=0) - margin\n", + "\n", + " self.nx = nx\n", + " self.ny = ny\n", + " self.nz = nz\n", + " self.xs = np.linspace(0, 1, nx)\n", + " self.ys = np.linspace(0, 1, ny)\n", + " self.zs = np.linspace(0, 1, nz)\n", + "\n", + " def get_coords(self):\n", + " frac_coords = np.stack(np.meshgrid(self.xs, self.ys, self.zs), axis=-1)\n", + " # permuting x<->y is necessary to match weird ordering of cube format\n", + " return np.einsum(\"yxzi,ij->xyzj\", frac_coords, self.box) + self.boxorig\n", + "\n", + " def __call__(self, field) -> str:\n", + " assert field.ndim == 3\n", + " assert field.shape == (self.nx, self.ny, self.nz)\n", + " comment = \"\"\n", + " string = \"\"\n", + " string += comment + \"\\n\"\n", + " string += \"Created by OneQMC CubeFormatter\\n\"\n", + " string += f\"{len(self.coords):5d}\"\n", + " string += \"{:12.6f}{:12.6f}{:12.6f}\\n\".format(*tuple((self.boxorig).tolist()))\n", + " dx = self.xs[-1] if len(self.xs) == 1 else self.xs[1]\n", + " dy = self.ys[-1] if len(self.ys) == 1 else self.ys[1]\n", + " dz = self.zs[-1] if len(self.zs) == 1 else self.zs[1]\n", + " delta = (self.box.T * np.stack([dx, dy, dz])).T\n", + " string += (\n", + " f\"{self.nx:5d}{delta[0,0]:12.6f}{delta[0,1]:12.6f}{delta[0,2]:12.6f}\\n\"\n", + " )\n", + " string += (\n", + " f\"{self.ny:5d}{delta[1,0]:12.6f}{delta[1,1]:12.6f}{delta[1,2]:12.6f}\\n\"\n", + " )\n", + " string += (\n", + " f\"{self.nz:5d}{delta[2,0]:12.6f}{delta[2,1]:12.6f}{delta[2,2]:12.6f}\\n\"\n", + " )\n", + " for charge, coord in zip(self.charges, self.coords):\n", + " string += \"%5d%12.6f\" % (charge, 0.0)\n", + " string += \"{:12.6f}{:12.6f}{:12.6f}\\n\".format(*tuple((coord).tolist()))\n", + "\n", + " # Sync to CPU if not there already\n", + " field = np.asarray(field)\n", + " for ix in range(self.nx):\n", + " for iy in range(self.ny):\n", + " for iz0, iz1 in pyscf.lib.prange(0, self.nz, 6):\n", + " fmt = \"%13.5E\" * (iz1 - iz0) + \"\\n\"\n", + " string += fmt % tuple(field[ix, iy, iz0:iz1].tolist())\n", + "\n", + " return string\n", + "\n", + "\n", + "def show_isosurface(\n", + " field_data,\n", + " iso_value: float = 0.05,\n", + " view=None,\n", + "):\n", + " if view is None:\n", + " view = py3Dmol.view()\n", + " view.addVolumetricData(\n", + " field_data,\n", + " \"cube\",\n", + " {\n", + " \"isoval\": -iso_value,\n", + " \"smoothness\": 5,\n", + " \"opacity\": 0.8,\n", + " \"volformat\": \"cube\",\n", + " \"color\": \"blue\",\n", + " },\n", + " )\n", + " view.addVolumetricData(\n", + " field_data,\n", + " \"cube\",\n", + " {\n", + " \"isoval\": iso_value,\n", + " \"smoothness\": 5,\n", + " \"opacity\": 0.8,\n", + " \"volformat\": \"cube\",\n", + " \"color\": \"red\",\n", + " },\n", + " )\n", + " return view" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2af67290-e369-49ac-8392-84f2403e6f48", + "metadata": {}, + "outputs": [], + "source": [ + "formatter = CubeDataFormatter(naphthalene, nx=100, nz=100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5149c1cf-348e-4c14-a583-b08eec7deec1", + "metadata": {}, + "outputs": [], + "source": [ + "grid = formatter.get_coords()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6443480-e5fb-4628-8fe5-36cf5acd3bad", + "metadata": {}, + "outputs": [], + "source": [ + "def envelope_fn(x, exponents, centers, coefs):\n", + " r = np.linalg.norm(centers - x[..., None, :], axis=-1)\n", + " exps = np.exp(-exponents * r[..., None])\n", + " return np.einsum(\"...ij,ij->...\", exps, coefs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2cc031a", + "metadata": {}, + "outputs": [], + "source": [ + "orb_idx = 4\n", + "coef = params[\"naphthalene_2_se_envelope_up_feature_selector\"][\n", + " orb_idx, :, :, 0\n", + "]\n", + "field_values = envelope_fn(\n", + " grid,\n", + " params[\"naphthalene_2_exponents\"].squeeze(0),\n", + " naphthalene.coords,\n", + " coef,\n", + ")\n", + "# Easier to visualize if it integrates to 1\n", + "field_values /= np.abs(field_values).sum()\n", + "field_data = formatter(field_values)\n", + "view = show_mol(naphthalene)\n", + "view.rotate(-90, 'z')\n", + "show_isosurface(field_data, view=view, iso_value=1e-5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5982bf5", + "metadata": {}, + "outputs": [], + "source": [ + "orb_idx = 9\n", + "coef = params[\"naphthalene_2_se_envelope_up_feature_selector\"][\n", + " orb_idx, :, :, 0\n", + "]\n", + "field_values = envelope_fn(\n", + " grid,\n", + " params[\"naphthalene_2_exponents\"].squeeze(0),\n", + " naphthalene.coords,\n", + " coef,\n", + ")\n", + "# Easier to visualize if it integrates to 1\n", + "field_values /= np.abs(field_values).sum()\n", + "field_data = formatter(field_values)\n", + "view = show_mol(naphthalene)\n", + "view.rotate(-90, 'z')\n", + "show_isosurface(field_data, view=view, iso_value=1e-5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74e0edb1-c828-49f7-bfdb-26cdc861a056", + "metadata": {}, + "outputs": [], + "source": [ + "orb_idx = 12\n", + "coef = params[\"naphthalene_2_se_envelope_up_feature_selector\"][\n", + " orb_idx, :, :, 0\n", + "]\n", + "field_values = envelope_fn(\n", + " grid,\n", + " params[\"naphthalene_2_exponents\"].squeeze(0),\n", + " naphthalene.coords,\n", + " coef,\n", + ")\n", + "# Easier to visualize if it integrates to 1\n", + "field_values /= np.abs(field_values).sum()\n", + "field_data = formatter(field_values)\n", + "view = show_mol(naphthalene)\n", + "view.rotate(90, 'z')\n", + "show_isosurface(field_data, view=view, iso_value=1e-5)" + ] + }, + { + "cell_type": "markdown", + "id": "5383601a", + "metadata": {}, + "source": [ + "## Structure figures" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b3324e2-5693-4af0-b535-c5160d9f960b", + "metadata": {}, + "outputs": [], + "source": [ + "view = show_mol(naphthalene)\n", + "view.rotate(90, \"z\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d47089d5", + "metadata": {}, + "outputs": [], + "source": [ + "view = show_mol(azulene)\n", + "view.rotate(36, \"z\")" + ] + }, + { + "cell_type": "markdown", + "id": "30f11cbb", + "metadata": {}, + "source": [ + "## Convergence plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6478936b", + "metadata": {}, + "outputs": [], + "source": [ + "%config InlineBackend.figure_format = 'retina'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b1c5446", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.offsetbox import OffsetImage, AnnotationBbox\n", + "from PIL import Image\n", + "\n", + "training_iterations = [48_000, 56_000, 64_000, 72_000, 80_000]\n", + "reaction_energies = [ufloat(31.87, 1.2), ufloat(36.6368, 0.7353), ufloat(37.5258, 0.6017), ufloat(36.6801, 0.5550), ufloat(36.9793, 0.9101)]\n", + "ccsdpt_cbs_reference = 36.82\n", + "\n", + "fig, ax = plt.subplots(figsize=(6, 3.5))\n", + "ax.scatter(\n", + " training_iterations,\n", + " [r.n for r in reaction_energies],\n", + " color=\"#1C909F\",\n", + " label=\"Orbformer\",\n", + ")\n", + "ax.axhline(ccsdpt_cbs_reference, color=\"k\", linestyle=\"--\", label=\"CCSD(T)/CBS\")\n", + "ax.fill_between(\n", + " [40_000, 90_000],\n", + " [ccsdpt_cbs_reference - 1.0],\n", + " [ccsdpt_cbs_reference + 1.0],\n", + " alpha=0.2,\n", + " color=\"gray\",\n", + ")\n", + "ax.set_xlabel(\"Orbformer finetuning iterations\")\n", + "ax.set_ylabel(\"Reaction energy (kcal/mol)\")\n", + "ax.set_ylim(30.5, 41.0)\n", + "ax.set_xlim(44_000, 85_000)\n", + "ax.set_xticks([50_000, 60_000, 70_000, 80_000])\n", + "ax.set_yticks([32, 34, 36, 38, 40])\n", + "ax.legend(loc=(0.02, 0.75))\n", + "\n", + "struct_dir = \"../../experiment_results/09_azulene_naphthalene/molecule_images\"\n", + "naphthalene_img = np.array(Image.open(f\"{struct_dir}/naphthalene-structure.png\"))\n", + "azulene_img = np.array(Image.open(f\"{struct_dir}/azulene-structure.png\"))\n", + "\n", + "x_pos, y_pos = 76_500, 33.0\n", + "gap = 17000\n", + "zoom = 0.13\n", + "im_naph = OffsetImage(naphthalene_img, zoom=zoom)\n", + "ab_naph = AnnotationBbox(im_naph, (x_pos - gap, y_pos), frameon=False)\n", + "ax.add_artist(ab_naph)\n", + "im_azul = OffsetImage(azulene_img, zoom=zoom)\n", + "ab_azul = AnnotationBbox(im_azul, (x_pos, y_pos), frameon=False)\n", + "ax.add_artist(ab_azul)\n", + "\n", + "ax.annotate(\n", + " \"\",\n", + " xy=(x_pos - gap / 2 + 1200, y_pos),\n", + " xytext=(x_pos - gap / 2 - 2200, y_pos),\n", + " arrowprops=dict(arrowstyle=\"->\", lw=1.5, color=\"k\"),\n", + ")\n", + "fig.savefig(\"reaction_energy_convergence.pdf\", dpi=300, bbox_inches=\"tight\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30d88a9a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "env", + "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.13.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}