{ "cells": [ { "cell_type": "markdown", "id": "2887ae18", "metadata": {}, "source": [ "# ADLD - HFLD\n", "\n", "This tutorial demonstrates how to use ORCA together with the ORCA python interface (OPI) to perform the **[Atomic Decomposition of London Dispersion energy (ADLD)](https://doi.org/10.1021/acscentsci.5c00356)** analysis to decompose the London dispersion energy of the **[Hartree-Fock plus London Dispersion (HFLD)](https://doi.org/10.1021/acs.jctc.9b00425)** interaction energy within a dimer into atomic contributions. Decomposing the HFLD energy allows a non-empirical estimate of atomic London dispersion contributions. For the quicker but more empirical ADLD(D4) decomposition, see the corresponding notebook. In comparison to ADLD(D4) only one calculation of the whole system is required, since the subsystems consisiting of only one fragment have no contribution.\n", "\n", "In this notebook we will:\n", "1. Import the required python dependencies.\n", "2. Define a working directory.\n", "3. Set up the input structure.\n", "4. Run a HFLD calculation for the supersystem.\n", "7. Process the ADLD results.\n", "8. Visualize the results.\n", "9. Generating the Cube File." ] }, { "cell_type": "markdown", "id": "0e7b75a9", "metadata": {}, "source": [ "## Step 1: Import Dependencies\n", "\n", "We start by importing the modules needed for:\n", "- Interfacing with ORCA input/output\n", "- Numerical calculations and data handling\n", "- Plotting results\n", "\n", "> **Note:** We additionally import modules for visualization/plotting like `py3Dmol`. For this, it might be necessary to install `py3Dmol` into your OPI `venv` (e.g., by activating the `.venv` and using `uv pip install py3Dmol`)." ] }, { "cell_type": "code", "execution_count": 1, "id": "78657414", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "import shutil\n", "\n", "# > pandas and numpy for data handling\n", "import numpy as np\n", "\n", "# > OPI imports for performing ORCA calculations and reading output\n", "from opi.core import Calculator\n", "from opi.output.core import Output\n", "from opi.input.simple_keywords import BasisSet, Dlpno, \\\n", " Scf, AuxBasisSet, Approximation\n", "from opi.input.structures.structure import Structure\n", "from opi.input.blocks import BlockFrag\n", "\n", "# > For visualization of molecules\n", "import py3Dmol" ] }, { "cell_type": "markdown", "id": "28791864", "metadata": {}, "source": [ "## Step 2: Working Directory and Conversion Factor\n", "\n", "We define a subfolder `RUN` in which the actual ORCA calculations will take place. Also, we define a conversion factor, since we want the resulting interaction energies in kcal/mol for better interpretability." ] }, { "cell_type": "code", "execution_count": null, "id": "3e349fdf", "metadata": {}, "outputs": [], "source": [ "# > Calculation is performed in `RUN`\n", "working_dir = Path(\"RUN\")\n", "# > The `working_dir`is automatically (re-)created\n", "shutil.rmtree(working_dir, ignore_errors=True)\n", "working_dir.mkdir()\n", "# > Conversion factor for atomic units to kcal/mol\n", "unit_conversion = 627.509 " ] }, { "cell_type": "markdown", "id": "5f638574", "metadata": {}, "source": [ "## Step 3: Setup the Input Structure\n", "\n", "As an example we will decompose the LD interaction in an cubane dimer. The 3D structure in Cartesian coordinates is defined and visualized:\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "61fac6de", "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", "text/html": [ "
\n", "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# > Define cartesian coordinates in Angstroem as python string \n", "# > In principle, this can be replaced with other neutral NCI dimers\n", "xyz_data = \"\"\"\\\n", "32\n", "cubane dimer\n", "C -2.37438447 -0.01391040 0.09884370\n", "H -1.28799640 -0.03512861 0.15212364\n", "C -4.14778752 -0.93512832 -0.83041188\n", "C -4.12750489 1.22660964 -0.39560995\n", "C -3.20643343 0.25290551 -1.19176459\n", "C -4.23643748 -0.22901780 1.25693113\n", "C -3.31541388 -1.20282062 0.46084844\n", "C -3.29511100 0.95905424 0.89568028\n", "C -5.06903241 0.03775498 -0.03358762\n", "H -4.49182155 -1.69883200 -1.52468842\n", "H -4.45488733 2.20468479 -0.74189213\n", "H -2.79034514 0.44620851 -2.17829474\n", "H -4.65224637 -0.42236391 2.24356374\n", "H -2.98771769 -2.18079631 0.80710684\n", "H -2.95068094 1.72266937 1.58984614\n", "H -6.15558078 0.05880427 -0.08645048\n", "C 2.05734830 0.00538088 -0.06151765\n", "H 0.97096030 0.02659937 -0.11479860\n", "C 2.97807526 -0.96758450 -0.85835285\n", "C 2.99837839 1.19429062 -0.42352221\n", "C 3.91940243 0.22048707 -1.21960352\n", "C 2.88939603 -0.26143458 1.22909153\n", "C 3.81046792 -1.23513946 0.43293829\n", "C 3.83075080 0.92659878 0.86773902\n", "C 4.75199612 -0.04628527 0.07091613\n", "H 2.63364560 -1.73119990 -1.55251858\n", "H 2.67068285 2.17226622 -0.76978149\n", "H 4.33521233 0.41383252 -2.20623583\n", "H 2.47330675 -0.45473693 2.21562140\n", "H 4.13784970 -2.21321450 0.77922136\n", "H 4.17478444 1.69030273 1.56201542\n", "H 5.83854442 -0.06733485 0.12378001\n", "\n", "\\n\n", "\"\"\"\n", "\n", "# > Visualize the input structure\n", "view = py3Dmol.view(width=400, height=400)\n", "view.addModel(xyz_data, 'xyz')\n", "view.setStyle({}, {'stick': {}, 'sphere': {'scale': 0.3}})\n", "view.zoomTo()\n", "view.show()\n", "\n", "# > Write the input structure to a file\n", "with open(working_dir / \"struc.xyz\",\"w\") as f:\n", " f.write(xyz_data)\n", "# > Read structure into object\n", "structure = Structure.from_xyz(working_dir / \"struc.xyz\")" ] }, { "cell_type": "markdown", "id": "951edfa2", "metadata": {}, "source": [ "## Step 4: Run HFLD for the Supramolecular System\n", "\n", "We run a HFLD/def2-TZVP calculation for the complete system. We choose the moderate def2-TZVP basis set, because the HFLD London dispersion energy is dependent on the basis and this is a good compromise for computational costs vs. accuracy. The fragments in the calculation will automatically be defined by ORCA (in this example two fragments, one for each monomer). We define and run three helper functions: `setup_calc`, `run_calc`, and `check_and_parse_output`. \n", "\n", "First, we set up the calculation with `setup_calc`:" ] }, { "cell_type": "code", "execution_count": 4, "id": "2af6539f", "metadata": {}, "outputs": [], "source": [ "def setup_calc(basename : str, working_dir: Path, structure: Structure, ncores: int = 4) -> Calculator:\n", " # > Set up a Calculator object\n", " calc = Calculator(basename=basename, working_dir=working_dir)\n", " # > Assign structure to calculator\n", " calc.structure = structure\n", "\n", " # > Define a simple keyword list for the ORCA calculation.\n", " sk_list = [\n", " Dlpno.HFLD,\n", " BasisSet.DEF2_TZVP,\n", " AuxBasisSet.DEF2_TZVP_C,\n", " Scf.NORMALSCF,\n", " Approximation.RIJCOSX,\n", " Dlpno.ADLD,\n", " Dlpno.LED\n", " ]\n", "\n", " # > Use simple keywords in calculator\n", " calc.input.add_simple_keywords(*sk_list)\n", "\n", " # > Automatically assign fragments by setting the frag block\n", " calc.input.add_blocks(BlockFrag(storefrags=True))\n", "\n", " # Define number of CPUs for the calcualtion\n", " calc.input.ncores = ncores # > 4 CPUs for this ORCA run\n", "\n", " return calc\n", "\n", "calc_supersystem = setup_calc(\"supersystem\", working_dir=working_dir, structure=structure)" ] }, { "cell_type": "markdown", "id": "bd19da89", "metadata": {}, "source": [ "Then, we run the calculation with `run_calc` and obtain the output:" ] }, { "cell_type": "code", "execution_count": 5, "id": "649537a5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running ORCA calculation ... Done\n" ] } ], "source": [ "def run_calc(calc: Calculator) -> Output:\n", " # > Write the ORCA input file\n", " calc.write_input()\n", " # > Run the ORCA calculation\n", " print(\"Running ORCA calculation ...\", end=\"\")\n", " calc.run()\n", " print(\" Done\")\n", "\n", " # > Get the output object\n", " output = calc.get_output()\n", " \n", " return output\n", "\n", "output_supersystem = run_calc(calc_supersystem)" ] }, { "cell_type": "markdown", "id": "2b46cbb5", "metadata": {}, "source": [ "The successful calculation has to be checked and parsed, which is done in `check_and_parse_output`:" ] }, { "cell_type": "code", "execution_count": 6, "id": "da89691a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ORCA terminated normally and the SCF converged\n" ] } ], "source": [ "def check_and_parse_output(output: Output):\n", " # > Check for proper termination of ORCA\n", " status = output.terminated_normally()\n", " if not status:\n", " # > ORCA did not terminate normally\n", " raise RuntimeError(f\"ORCA did not terminate normally, see output file: {output.get_outfile()}\")\n", " else:\n", " # > ORCA did terminate normally so we can parse the output\n", " output.parse()\n", "\n", " # Now check for convergence of the SCF\n", " if output.results_properties.geometries[0].single_point_data.converged:\n", " print(\"ORCA terminated normally and the SCF converged\")\n", " else:\n", " raise RuntimeError(\"SCF DID NOT CONVERGE\")\n", " \n", "check_and_parse_output(output_supersystem)" ] }, { "cell_type": "markdown", "id": "71352834", "metadata": {}, "source": [ "## Step 5: Processing of ADLD Energies\n", "\n", "Now we perform the post-processing for the ADLD. First, we gather atomic contributions of the LD. " ] }, { "cell_type": "code", "execution_count": 7, "id": "5d317a34", "metadata": {}, "outputs": [], "source": [ "# > Path to atomic contributions to the LD in the dimer (product) [Eh]\n", "product_vdw = output_supersystem.results_properties.geometries[0].mdci_adld.adlddispatomic_loewdin\n", "\n", "# > Energies are already relative since reactant contributions are zero\n", "# > Apply conversion factor for result in kcal/mol\n", "atomic_disp: list[float] = []\n", "for i in range(len(product_vdw)):\n", " product_x = float(product_vdw[i][0])\n", " product_x_kcal = (product_x) * unit_conversion\n", " atomic_disp.append(product_x_kcal)" ] }, { "cell_type": "markdown", "id": "956bfcab", "metadata": {}, "source": [ "Now we have an array of atomic contributions to the relative London dispersion interaction energies. A London dispersion density ($\\rho_{disp}$) can be defined for easily visualizaing and analyzing the atomic contributions:" ] }, { "cell_type": "code", "execution_count": 8, "id": "8ab273e9", "metadata": {}, "outputs": [], "source": [ "def rho_disp(r, alpha, epsilons, positions):\n", " \"\"\"\n", " Calculate the dispersion density rho_disp at position r.\n", "\n", " Parameters:\n", " - r: (3,) array-like, the 3D point to evaluate at.\n", " - alpha: float, parameter alpha (typically set to 0.3-0.5).\n", " - epsilons: (N,) array-like, dispersion strengths for each atom A.\n", " - positions: (N, 3) array-like, positions R_A of atoms.\n", "\n", " Returns:\n", " - rho: float, the dispersion density at r.\n", " \"\"\"\n", " # Pre-factor\n", " prefactor = (np.pi / alpha)**(-1.5)\n", " \n", " # Sum over atoms\n", " r = np.array(r)\n", " positions = np.array(positions)\n", " deltas = r - positions # shape (N, 3)\n", " squared_distances = np.sum(deltas**2, axis=1) # shape (N,)\n", " \n", " contributions = epsilons * np.exp(-alpha * squared_distances)\n", " total_sum = np.sum(contributions)\n", " \n", " return prefactor * total_sum" ] }, { "cell_type": "markdown", "id": "feed6238", "metadata": {}, "source": [ "## Step 6: Crude Visualization of ADLD Contributions\n", "\n", "In the original ADLD publication, the `rho_disp` function was evaluated on a grid and the results was written to a `.cube` file for visualization. For the sake of a simple visualization within this notebook, we will first evaluate this function only on the atomic positions and color the atoms by the values at their centers. These plots offer a qualitative analysis of the data and may differ in detail from those published. Further below, we also generate the cube file for the full visualization. We refer the interested reader to the [LDDsuite](https://github.com/bistonigroup/LDDsuite) on GitHub.\n", "\n", "We start by calculating $\\rho_{disp}$ at the atomic positions:" ] }, { "cell_type": "code", "execution_count": 9, "id": "b97f203c", "metadata": {}, "outputs": [], "source": [ "# > Path to the Cartesian coordinates of the dimer\n", "elements = output_supersystem.results_properties.geometries[0].geometry.coordinates.cartesians\n", "\n", "# > List of element symbols\n", "elem_list: list[str] = []\n", "# > List of atomic positions\n", "positions: list[tuple[float,float,float]] = []\n", "\n", "# > Fill the lists\n", "for elem in elements:\n", " elem_list.append(elem[0])\n", " positions.append([elem[1],elem[2],elem[3]])\n", "\n", "# > Calculate the London dispersion density (rho_disp) at atomic positions\n", "rho_disp_atomic: list[float] = []\n", "for i in range(len(positions)):\n", " # > alpha value of 0.3 is selected\n", " rho_disp_atomic.append(rho_disp(positions[i], 0.3, atomic_disp, positions))" ] }, { "cell_type": "markdown", "id": "1bc6d408", "metadata": {}, "source": [ "Next, we define a small helper function for mapping the values of the dispersion density to a color gradient:" ] }, { "cell_type": "code", "execution_count": 10, "id": "28224d51", "metadata": {}, "outputs": [], "source": [ "def value_to_color(value, min_val):\n", " \"\"\"\n", " Maps a value (min_val to 0) to a white-to-red gradient.\n", " Numbers near 0 -> white; near min_val -> red.\n", " \"\"\"\n", " if value >= 0:\n", " return '#FFFFFF'\n", " \n", " # Ensure min_val is negative and avoid division by zero\n", " if min_val == 0:\n", " min_val = -1e-6 # small epsilon to avoid division by zero\n", " \n", " # Normalize: 0 (white) to min_val (red)\n", " norm = min(abs(value / min_val), 1.0)\n", " \n", " r = 255\n", " g = int(255 * (1 - norm))\n", " b = int(255 * (1 - norm))\n", " \n", " return f'#{r:02X}{g:02X}{b:02X}'" ] }, { "cell_type": "markdown", "id": "a3b6865c", "metadata": {}, "source": [ "Finally, we plot $\\rho_{disp}$ by adding colored spheres to the py3Dmol visualization. We require custom vdw radii for the sizes of the spheres to add:" ] }, { "cell_type": "code", "execution_count": 11, "id": "05d83efe", "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", "text/html": [ "
\n", "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# > Custom vdw radii for visualization\n", "vdw_radii = {\n", " 'H': 1.20, 'He': 1.40, 'Li': 1.82, 'Be': 1.53,\n", " 'B': 1.92, 'C': 1.70, 'N': 1.55, 'O': 1.52,\n", " 'F': 1.47\n", "}\n", "\n", "# > Definition of maximum and minimum values for color gradient\n", "min_val = -0.01 # > kcal/mol Bohr³\n", "\n", "# > Setup the viewer\n", "view = py3Dmol.view(width=400, height=400)\n", "view.addModel(xyz_data)\n", "\n", "# > Initial style: white stick representation for all atoms\n", "view.setStyle({}, {'stick': {'color' : 'white'}})\n", "\n", "# Apply spheres as overlays for color-coding\n", "for i, val in enumerate(rho_disp_atomic):\n", " color = value_to_color(val, min_val)\n", " radius = vdw_radii[elem_list[i]] * 0.35 # > arbitrary scaling factor for nice spheres\n", " selection = {'serial': i}\n", " view.addStyle(selection, {'sphere': {'color': color, 'radius': radius}})\n", "\n", "view.zoomTo()\n", "view.show()" ] }, { "cell_type": "markdown", "id": "68b8fc69", "metadata": {}, "source": [ "## Step 7: Generating the Cube File\n", "\n", "For the correct visualization, also the cube file can be generated that might be plotted separately:" ] }, { "cell_type": "code", "execution_count": 12, "id": "69994a25", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Cube file written: rho_disp.cube\n" ] } ], "source": [ "# > Helper dict to obtain atomic numbers\n", "element_dict = {\n", " \"H\": 1,\n", " \"He\": 2,\n", " \"Li\": 3,\n", " \"Be\": 4,\n", " \"B\": 5,\n", " \"C\": 6,\n", " \"N\": 7,\n", " \"O\": 8,\n", " \"F\": 9\n", "}\n", "\n", "# > Helper function for defining bounding box\n", "def get_bounding_box(coords, padding):\n", " coords = np.array(coords)\n", " min_coords = coords.min(axis=0) - padding\n", " max_coords = coords.max(axis=0) + padding\n", " return min_coords, max_coords\n", "\n", "# > Helper function for writing the actual cube file\n", "def write_cube_file(filename, coords, atomic_numbers, values, min_coords, voxel_size, n_points):\n", " with open(filename, 'w') as f:\n", " f.write(\"Cube file generated by adld notebook\\n\")\n", " f.write(\"rho_disp\\n\")\n", " f.write(f\"{len(coords):5d} {min_coords[0]:12.6f} {min_coords[1]:12.6f} {min_coords[2]:12.6f}\\n\")\n", " f.write(f\"{n_points[0]:5d} {voxel_size[0]:12.6f} 0.0 0.0\\n\")\n", " f.write(f\"{n_points[1]:5d} 0.0 {voxel_size[1]:12.6f} 0.0\\n\")\n", " f.write(f\"{n_points[2]:5d} 0.0 0.0 {voxel_size[2]:12.6f}\\n\")\n", " for atom_num, (x, y, z) in zip(atomic_numbers, coords):\n", " f.write(f\"{atom_num:5d} 0.0 {x:12.6f} {y:12.6f} {z:12.6f}\\n\")\n", " flat_values = values.flatten(order='F') # Fortran order for cube files\n", " for i in range(0, len(flat_values), 6):\n", " line = \"\".join(f\"{v:13.5e}\" for v in flat_values[i:i+6])\n", " f.write(line + \"\\n\")\n", "\n", "\n", "# > Define the box size\n", "padding = 3.0 # 3 Å padding around the molecule\n", "min_coords, max_coords = get_bounding_box(positions, padding)\n", "\n", "# > Define the grid\n", "n_points = (40, 40, 40)\n", "x_lin = np.linspace(min_coords[0], max_coords[0], n_points[0])\n", "y_lin = np.linspace(min_coords[1], max_coords[1], n_points[1])\n", "z_lin = np.linspace(min_coords[2], max_coords[2], n_points[2])\n", "voxel_size = (max_coords - min_coords) / (np.array(n_points) - 1)\n", "\n", "# Evaluate rho_disp on the grid\n", "values = np.zeros(n_points)\n", "for i, x in enumerate(x_lin):\n", " for j, y in enumerate(y_lin):\n", " for k, z in enumerate(z_lin):\n", " r = (x, y, z)\n", " # > we set alpha to 0.3\n", " values[i, j, k] = rho_disp(r, 0.3, atomic_disp, positions)\n", "\n", "# > obtain atomic numbers\n", "atomic_numbers = [element_dict[symbol] for symbol in elem_list]\n", "\n", "# Write to cube file\n", "write_cube_file(working_dir / \"rho_disp.cube\", positions, atomic_numbers, values, min_coords, voxel_size, n_points)\n", "print(\"Cube file written: rho_disp.cube\")" ] }, { "cell_type": "markdown", "id": "1bda6268", "metadata": {}, "source": [ "## Summary\n", "\n", "In this notebook we demonstrated how the ORCA Python interface can be employed to decompose the London dispersion energy from HFLD into atomic contributions. We performed a crude visualization within this notebook and generated a cube file for the full analysis." ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all", "executable": "/usr/bin/env python3", "main_language": "python", "notebook_metadata_filter": "-all" }, "kernelspec": { "display_name": ".venv", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }