{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Inferring topology from atom coordinates\n\n.. currentmodule:: biotite.structure\n\nQuantum chemistry programs typically report the result of a geometry optimization as a\nbare list of elements and coordinates, for example in the *XYZ* format.\nWhich atoms are bonded to each other, let alone the bond orders and formal charges, is\nnot part of such an output.\n\nThis information is required for a lot of downstream tasks though, so this example\ndemonstrates how to reconstruct the bond graph (i.e. the topology) from coordinates\nalone.\nFinally, the assigned bond orders are compared with the *Wiberg bond indices* from the\nDFT calculation.\n\nAs input the DFT-optimized geometry (B3LYP-D3BJ/DZVP) of the antibiotic *ciprofloxacin*\nis fetched from the *MolSSI QCArchive* :footcite:`Smith2021`, more precisely from the\n*FDA Optimization Dataset 1* of the *Open Force Field Initiative*.\nCiprofloxacin is an instructive test case, as it is a zwitterion:\nThe piperazine nitrogen is protonated and the carboxyl group is deprotonated, so the\nformal charges must be placed correctly, although the total charge of the molecule is\nzero.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Code source: Patrick Kunzmann\n# License: BSD 3 clause\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport requests\nimport biotite\nimport biotite.interface.pymol as pymol_interface\nimport biotite.structure as struc\n\nQCARCHIVE_URL = \"https://api.qcarchive.molssi.org/api/v1/\"\n# The QCArchive ID of the molecule to fetch:\n# This is the final geometry of the optimization record '6095247', i.e. the first\n# conformer of ciprofloxacin in the 'FDA Optimization Dataset 1'\nMOLECULE_ID = 3954680\n# The QCArchive ID of the final single point calculation of that optimization,\n# which provides the Wiberg bond indices\nRECORD_ID = 6357374\n# QCArchive reports coordinates in atomic units\nBOHR_TO_ANGSTROM = 0.529177210903"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The *QCArchive* REST API returns the molecule as JSON, which contains the element\nsymbols and a flat list of coordinates, among other data.\nOnly these two fields are used to build the :class:`AtomArray`, exactly as if the\ngeometry was read from an *XYZ* file.\nThe total charge of the molecule is taken from the record as well, as it is required\nas constraint for the bond order assignment later.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "response = requests.get(QCARCHIVE_URL + f\"molecules/{MOLECULE_ID}\")\nresponse.raise_for_status()\nrecord = response.json()\n\nmolecule = struc.AtomArray(len(record[\"symbols\"]))\nmolecule.coord = np.array(record[\"geometry\"]).reshape(-1, 3) * BOHR_TO_ANGSTROM\nmolecule.element = np.char.upper(record[\"symbols\"])\nmolecule.res_name[:] = \"UNL\"\nmolecule.res_id[:] = 1\nmolecule.hetero[:] = True\nmolecule.atom_name = struc.create_atom_names(molecule)\ntotal_charge = int(record[\"molecular_charge\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The topology is reconstructed in two steps:\n:func:`connect_via_distances()` determines which atoms are bonded to each other based\non their distances and covalent radii, and :func:`infer_bond_types()` assigns bond\norders and formal charges to these bonds, so that each atom satisfies its valence and\nthe formal charges sum up to the total charge of the molecule.\nThe formal charges are stored as ``charge`` annotation for use in the following steps.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "molecule.bonds = struc.connect_via_distances(molecule)\nmolecule.bonds, charges = struc.infer_bond_types(molecule, total_charge)\nmolecule.set_annotation(\"charge\", charges)\n\nprint(\"Number of bonds:\", molecule.bonds.get_bond_count())\nbond_types = molecule.bonds.as_array()[:, 2]\nfor bond_type in np.unique(bond_types):\n    count = np.count_nonzero(bond_types == bond_type)\n    print(f\"{struc.BondType(bond_type).name}: {count}\")\nprint(\"Charged atoms:\")\nfor i in np.where(molecule.charge != 0)[0]:\n    print(f\"    {molecule.atom_name[i]}: {molecule.charge[i]:+d}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the visualization the molecule is rendered with *PyMOL*.\nThe heavy atoms are labeled with their names to be able to identify them in the\nfollowing analysis, and the formal charges are indicated by ``+``/``-`` signs.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pymol_obj = pymol_interface.PyMOLObject.from_structure(molecule)\npymol_obj.show(\"sticks\")\npymol_obj.color(\"white\", molecule.element == \"H\")\n# Lighter carbon atoms give a better contrast to the labels\npymol_obj.color((0.7, 0.7, 0.7), molecule.element == \"C\")\npymol_interface.cmd.set(\"stick_radius\", 0.15)\n# Draw double bonds as two parallel sticks\npymol_interface.cmd.set(\"valence\", 1)\npymol_interface.cmd.set(\"valence_size\", 0.1)\n# Label the heavy atoms with their names\nfor i in np.where(molecule.element != \"H\")[0]:\n    pymol_obj.label(i, molecule.atom_name[i])\npymol_interface.cmd.set(\"label_size\", 14)\n# Offset the labels from their atoms, so that they do not hide each other\npymol_interface.cmd.set(\"label_position\", (1.5, -1.5, 4))\n# As an atom can only have a single label, the formal charges are indicated by\n# labeled pseudoatoms at the positions of the charged atoms\nfor i in np.where(molecule.charge != 0)[0]:\n    pymol_interface.cmd.pseudoatom(\n        \"charges\",\n        pos=molecule.coord[i].tolist(),\n        label=\"+\" if molecule.charge[i] > 0 else \"-\",\n    )\npymol_interface.cmd.set(\"label_size\", 35, \"charges\")\npymol_interface.cmd.set(\"label_position\", (0.4, 0.4, 4), \"charges\")\npymol_interface.cmd.orient()\npymol_interface.show((1500, 1000))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "How well do the assigned integer bond orders reflect the electronic structure?\nThe DFT calculation provides the *Wiberg bond indices*, a quantum mechanical measure\nof the bond order between each pair of atoms.\nThey are stored in the *QCArchive* record of the final single point calculation of\nthe geometry optimization.\nIn the following plot the Wiberg index of each atom pair is compared with its assigned\nbond order, where pairs without a bond count as bond order zero.\nBonds whose Wiberg index deviates strongly from the assigned bond order are labeled.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "response = requests.get(QCARCHIVE_URL + f\"records/singlepoint/{RECORD_ID}\")\nresponse.raise_for_status()\nproperties = response.json()[\"properties\"]\nn_atoms = molecule.array_length()\nwiberg_indices = np.array(properties[\"wiberg lowdin indices\"]).reshape(n_atoms, n_atoms)\n\n# Consider each atom pair only once\natom_i, atom_j = np.triu_indices(n_atoms, k=1)\npair_wiberg_indices = wiberg_indices[atom_i, atom_j]\n# The values of the single, double and triple bond types correspond to the bond order,\n# non-bonded pairs are marked with -1 in the bond type matrix\nbond_orders = molecule.bonds.bond_type_matrix()[atom_i, atom_j].astype(int)\nbond_orders[bond_orders == -1] = 0\nhydrogen_mask = (molecule.element[atom_i] == \"H\") | (molecule.element[atom_j] == \"H\")\n\nfig, ax = plt.subplots(figsize=(8.0, 5.0))\n# Spread the points horizontally to reduce overlap\nx = bond_orders + np.random.default_rng(0).uniform(-0.15, 0.15, len(bond_orders))\nfor mask, color, label in [\n    (bond_orders == 0, \"gray\", \"No bond\"),\n    (\n        (bond_orders != 0) & hydrogen_mask,\n        biotite.colors[\"dimgreen\"],\n        \"Bond to hydrogen\",\n    ),\n    (\n        (bond_orders != 0) & ~hydrogen_mask,\n        biotite.colors[\"dimorange\"],\n        \"Heavy atom bond\",\n    ),\n]:\n    ax.scatter(x[mask], pair_wiberg_indices[mask], s=20, color=color, label=label)\nfor i in np.where(np.abs(pair_wiberg_indices - bond_orders) > 0.5)[0]:\n    # Put the label on the side of the point that faces the plot center\n    on_left_side = bond_orders[i] > 1\n    ax.annotate(\n        f\"{molecule.atom_name[atom_i[i]]}-{molecule.atom_name[atom_j[i]]}\",\n        (x[i], pair_wiberg_indices[i]),\n        xytext=(-8 if on_left_side else 8, 0),\n        textcoords=\"offset points\",\n        ha=\"right\" if on_left_side else \"left\",\n        va=\"center\",\n    )\nax.set_xticks([0, 1, 2])\nax.set_xticklabels([\"none\", \"single\", \"double\"])\nax.set_xlabel(\"Assigned bond order\")\nax.set_ylabel(\"Wiberg bond index\")\nax.legend(loc=\"upper left\")\nfig.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The Wiberg indices of all non-bonded pairs are far below the ones of bonded pairs,\nwhich confirms the connectivity found via atom distances.\nBonds to hydrogen and most single bonds between heavy atoms have Wiberg indices close\nto one and the carbonyl double bonds are close to two, as expected.\n\nThe labeled bonds are those, where the assigned integer bond order is furthest from\nthe Wiberg index.\nAll of them are cases of *resonance*:\nThe molecule cannot be described properly by a single Lewis structure, but only by\nthe superposition of multiple ones, so its electrons are delocalized over several\nbonds.\nThe assigned bond orders and formal charges represent only one of these resonance\nstructures, while the Wiberg indices reflect the actual delocalized electrons.\n\n## References\n\n.. footbibliography::\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}