Repository: chemosim-lab/ProLIF

Issue:
Tutorials improvements
The info is currently spread out between different notebooks which makes running a typical analysis from begin to end a bit too tedious.

New proposed structure: organize tutorials in sections:

### Section 1: Generating the fingerprint
This section should also include tips on how to speed things up (restricting to binding site residues and so on)
- MD trajectory:
  - ligand-protein
  - protein-protein
- Docking
  - Requirement on explicit hydrogens and/or bond orders/charges
- PDB
  - Structure preparation
- Troubleshooting

### Section 2: Analysis
- Data analysis with pandas
- Visualization


Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/chemosim-lab/ProLIF/pull/150?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chemosim-lab) Report
Patch coverage: **`75.00%`** and project coverage change: **`-0.14%`** :warning:
> Comparison is base [(`5acfa79`)](https://app.codecov.io/gh/chemosim-lab/ProLIF/commit/5acfa794489ed860bef6752411972693a897dd80?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chemosim-lab) 95.80% compared to head [(`77dcbe4`)](https://app.codecov.io/gh/chemosim-lab/ProLIF/pull/150?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chemosim-lab) 95.67%.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master     #150      +/-   ##
==========================================
- Coverage   95.80%   95.67%   -0.14%     
==========================================
  Files          19       19              
  Lines        1358     1363       +5     
==========================================
+ Hits         1301     1304       +3     
- Misses         57       59       +2     
```


| [Files Changed](https://app.codecov.io/gh/chemosim-lab/ProLIF/pull/150?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chemosim-lab) | Coverage Δ | |
|---|---|---|
| [prolif/parallel.py](https://app.codecov.io/gh/chemosim-lab/ProLIF/pull/150?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chemosim-lab#diff-cHJvbGlmL3BhcmFsbGVsLnB5) | `78.04% <0.00%> (ø)` | |
| [prolif/plotting/barcode.py](https://app.codecov.io/gh/chemosim-lab/ProLIF/pull/150?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chemosim-lab#diff-cHJvbGlmL3Bsb3R0aW5nL2JhcmNvZGUucHk=) | `78.88% <50.00%> (-1.58%)` | :arrow_down: |
| [prolif/fingerprint.py](https://app.codecov.io/gh/chemosim-lab/ProLIF/pull/150?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chemosim-lab#diff-cHJvbGlmL2ZpbmdlcnByaW50LnB5) | `98.24% <100.00%> (+0.02%)` | :arrow_up: |


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/chemosim-lab/ProLIF/pull/150?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chemosim-lab).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chemosim-lab).


Patch:
diff --git a/.binder/postBuild b/.binder/postBuild
deleted file mode 100644
index b6fa644..0000000
--- a/.binder/postBuild
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/bin/bash
-
-set -e
-
-pip install .
\ No newline at end of file
diff --git a/.readthedocs.yml b/.readthedocs.yml
index 97baeb0..b6c81c8 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -10,4 +10,4 @@ python:
   - method: pip
     path: .
 conda:
-  environment: .binder/environment.yml
+  environment: environment.yml
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1051bfd..aadd8a5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- Added a `display_residues` function to quickly visualize the residues in a `Molecule`
+  object.
 - Added a `Complex3D` class for plotting interactions in 3D. Added the corresponding
   `Fingerprint.plot_3d` method to generate the plot directly from an FP object.
 - Added a `Barcode` class for plotting interactions. Added the corresponding
diff --git a/README.rst b/README.rst
index 0633ddf..68cfb73 100644
--- a/README.rst
+++ b/README.rst
@@ -6,8 +6,6 @@ ProLIF
 
     * - **Documentation**
       - |docs|
-    * - **Tutorials**
-      - |binder|
     * - **CI**
       - |tests| |codecov| |codeql|
     * - **Builds**
@@ -25,32 +23,36 @@ interaction fingerprints for complexes made of ligands, protein, DNA or RNA mole
 extracted from molecular dynamics trajectories, docking simulations and experimental
 structures.
 
-You can try it out prior to any installation on `Binder <https://mybinder.org/v2/gh/chemosim-lab/ProLIF/HEAD?filepath=docs%2Fnotebooks>`_.
-
 Documentation
 -------------
 
-The installation instructions, documentation and tutorials can be found online on `ReadTheDocs <https://prolif.readthedocs.io/en/latest/>`_.
+The installation instructions, documentation and tutorials can be found online on
+`ReadTheDocs <https://prolif.readthedocs.io>`_.
 
 Issues
 ------
 
-If you have found a bug, please open an issue on the `GitHub Issues <https://github.com/chemosim-lab/ProLIF/issues>`_ page.
+If you have found a bug, please open an issue on the
+`GitHub Issues <https://github.com/chemosim-lab/ProLIF/issues>`_ page.
 
 Discussion
 ----------
 
-If you have questions on how to use ProLIF, or if you want to give feedback or share ideas and new features, please head to the `GitHub Discussions <https://github.com/chemosim-lab/ProLIF/discussions>`_ page.
+If you have questions on how to use ProLIF, or if you want to give feedback or share
+ideas and new features, please head to the
+`GitHub Discussions <https://github.com/chemosim-lab/ProLIF/discussions>`_ page.
 
 Citing ProLIF
 -------------
 
-Please refer to the `citation page <https://prolif.readthedocs.io/en/latest/source/citation.html>`_ on the documentation.
+Please refer to the `citation page <https://prolif.readthedocs.io/en/latest/source/citation.html>`_
+on the documentation.
 
 License
 -------
 
-Unless otherwise noted, all files in this directory and all subdirectories are distributed under the Apache License, Version 2.0 ::
+Unless otherwise noted, all files in this directory and all subdirectories are
+distributed under the Apache License, Version 2.0 ::
 
     Copyright 2017-2023 Cédric BOUYSSET
 
@@ -87,7 +89,7 @@ Unless otherwise noted, all files in this directory and all subdirectories are d
     :target: https://codecov.io/gh/chemosim-lab/ProLIF
 
 .. |docs| image:: https://readthedocs.org/projects/prolif/badge/?version=latest
-    :target: https://prolif.readthedocs.io/en/latest/?badge=latest
+    :target: https://prolif.readthedocs.io/
     :alt: Documentation Status
 
 .. |codeql| image:: https://github.com/chemosim-lab/ProLIF/workflows/CodeQL/badge.svg?branch=master
@@ -98,10 +100,6 @@ Unless otherwise noted, all files in this directory and all subdirectories are d
     :target: http://www.apache.org/licenses/LICENSE-2.0
     :alt: License
 
-.. |binder| image:: https://mybinder.org/badge_logo.svg
-    :target: https://mybinder.org/v2/gh/chemosim-lab/ProLIF/HEAD?filepath=docs%2Fnotebooks
-    :alt: Try it on binder
-
 .. |mdanalysis| image:: https://img.shields.io/badge/Powered%20by-MDAnalysis-orange.svg?logoWidth=16&logo=data:image/x-icon;base64,AAABAAEAEBAAAAEAIAAoBAAAFgAAACgAAAAQAAAAIAAAAAEAIAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAJD+XwCY/fEAkf3uAJf97wGT/a+HfHaoiIWE7n9/f+6Hh4fvgICAjwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACT/yYAlP//AJ///wCg//8JjvOchXly1oaGhv+Ghob/j4+P/39/f3IAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAJH8aQCY/8wAkv2kfY+elJ6al/yVlZX7iIiI8H9/f7h/f38UAAAAAAAAAAAAAAAAAAAAAAAAAAB/f38egYF/noqAebF8gYaagnx3oFpUUtZpaWr/WFhY8zo6OmT///8BAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgICAn46Ojv+Hh4b/jouJ/4iGhfcAAADnAAAA/wAAAP8AAADIAAAAAwCj/zIAnf2VAJD/PAAAAAAAAAAAAAAAAICAgNGHh4f/gICA/4SEhP+Xl5f/AwMD/wAAAP8AAAD/AAAA/wAAAB8Aov9/ALr//wCS/Z0AAAAAAAAAAAAAAACBgYGOjo6O/4mJif+Pj4//iYmJ/wAAAOAAAAD+AAAA/wAAAP8AAABhAP7+FgCi/38Axf4fAAAAAAAAAAAAAAAAiIiID4GBgYKCgoKogoB+fYSEgZhgYGDZXl5e/m9vb/9ISEjpEBAQxw8AAFQAAAAAAAAANQAAADcAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAjo6Mb5iYmP+cnJz/jY2N95CQkO4pKSn/AAAA7gAAAP0AAAD7AAAAhgAAAAEAAAAAAAAAAACL/gsAkv2uAJX/QQAAAAB9fX3egoKC/4CAgP+NjY3/c3Nz+wAAAP8AAAD/AAAA/wAAAPUAAAAcAAAAAAAAAAAAnP4NAJL9rgCR/0YAAAAAfX19w4ODg/98fHz/i4uL/4qKivwAAAD/AAAA/wAAAP8AAAD1AAAAGwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAALGxsVyqqqr/mpqa/6mpqf9KSUn/AAAA5QAAAPkAAAD5AAAAhQAAAAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADkUFBSuZ2dn/3V1df8uLi7bAAAATgBGfyQAAAA2AAAAMwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB0AAADoAAAA/wAAAP8AAAD/AAAAWgC3/2AAnv3eAJ/+dgAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA9AAAA/wAAAP8AAAD/AAAA/wAKDzEAnP3WAKn//wCS/OgAf/8MAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAIQAAANwAAADtAAAA7QAAAMAAABUMAJn9gwCe/e0Aj/2LAP//AQAAAAAAAAAA
     :alt: Powered by MDAnalysis
     :target: https://www.mdanalysis.org
diff --git a/docs/conf.py b/docs/conf.py
index 941d761..452cd9e 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -31,7 +31,7 @@ author = "Cédric Bouysset"
 # -- General configuration ---------------------------------------------------
 
 github_doc_root = "https://github.com/chemosim-lab/ProLIF/tree/master/docs/"
-needs_sphinx = "4.1.2"
+needs_sphinx = "5.3.0"
 
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
@@ -43,17 +43,24 @@ extensions = [
     "sphinx.ext.viewcode",
     "sphinx.ext.napoleon",
     "sphinx.ext.autosectionlabel",
-    "sphinx_rtd_theme",
     "recommonmark",
     "IPython.sphinxext.ipython_console_highlighting",
     "IPython.sphinxext.ipython_directive",
-    "nbsphinx",
+    "myst_nb",
+    "sphinx_copybutton",
 ]
+myst_enable_extensions = [
+    "colon_fence",
+]
+# nb_execution_raise_on_error = True
 
 mathjax_path = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js?config=TeX-AMS-MML_HTMLorMML"
 autosectionlabel_prefix_document = True
 napoleon_google_docstring = False
 
+# copybutton
+copybutton_exclude = ".linenos, .gp"
+
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ["_templates"]
 
@@ -62,15 +69,37 @@ templates_path = ["_templates"]
 # This pattern also affects html_static_path and html_extra_path.
 exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"]
 
-source_suffix = [".rst", ".md"]
+source_suffix = {
+    ".rst": "restructuredtext",
+    ".md": "myst-nb",
+    ".ipynb": "myst-nb",
+    ".myst": "myst-nb",
+}
 
 # -- Options for HTML output -------------------------------------------------
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
 #
-html_theme = "sphinx_rtd_theme"
+html_theme = "sphinx_book_theme"
 pygments_style = "sphinx"
+html_theme_options = {
+    "repository_url": "https://github.com/chemosim-lab/ProLIF",
+    "path_to_docs": "docs",
+    "use_source_button": True,
+    "use_download_button": True,
+    "use_repository_button": True,
+    "use_issues_button": True,
+    "launch_buttons": {"colab_url": "https://colab.research.google.com"},
+    "icon_links": [
+        {
+            "name": "GitHub",
+            "url": "https://github.com/chemosim-lab/ProLIF",
+            "icon": "fa-brands fa-square-github",
+            "type": "fontawesome",
+        }
+    ],
+}
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
@@ -83,7 +112,6 @@ intersphinx_mapping = {
     "https://docs.mdanalysis.org/stable/": None,
     "https://www.rdkit.org/docs/": None,
     "https://pandas.pydata.org/docs/": None,
-    "https://tqdm.github.io/docs/": None,
     "https://multiprocess.readthedocs.io/en/latest/": None,
 }
 
diff --git a/docs/index.rst b/docs/index.rst
index b68485f..10a09e0 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -4,7 +4,7 @@ Welcome to ProLIF's documentation!
 .. include:: ../README.rst
 
 .. toctree::
-   :caption: User guide
+   :caption: Information
    :maxdepth: 2
 
    source/installation
@@ -12,14 +12,10 @@ Welcome to ProLIF's documentation!
    source/changelog
 
 .. toctree::
-   :caption: Tutorials
+   :caption: User Guide
    :maxdepth: 2
-   :numbered:
 
-   notebooks/quickstart
-   notebooks/how-to
-   notebooks/protein-protein_interactions
-   notebooks/visualisation
+   source/tutorials
 
 .. toctree::
    :caption: Documentation
diff --git a/docs/notebooks/advanced.ipynb b/docs/notebooks/advanced.ipynb
new file mode 100644
index 0000000..7cc34a6
--- /dev/null
+++ b/docs/notebooks/advanced.ipynb
@@ -0,0 +1,408 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Advanced usage\n",
+    "\n",
+    "## Prerequisite\n",
+    "\n",
+    "We'll setup the protein and ligand files for the different usecases here:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import MDAnalysis as mda\n",
+    "import prolif as plf\n",
+    "\n",
+    "# load protein\n",
+    "protein_file = str(plf.datafiles.datapath / \"vina\" / \"rec.pdb\")\n",
+    "u = mda.Universe(protein_file)\n",
+    "protein_mol = plf.Molecule.from_mda(u)\n",
+    "\n",
+    "# load docking poses\n",
+    "poses_path = str(plf.datafiles.datapath / \"vina\" / \"vina_output.sdf\")\n",
+    "pose_iterable = plf.sdf_supplier(poses_path)\n",
+    "\n",
+    "# load 1 ligand from poses\n",
+    "ligand_mol = pose_iterable[0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Interactions\n",
+    "\n",
+    "This section shows how to modify existing interactions or defined new ones.\n",
+    "\n",
+    "### Modifying interaction parameters\n",
+    "\n",
+    "In ProLIF, each interaction is defined as a Python class. You can get a list of all the available interactions (including base abstract classes) with the following code snippet:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plf.Fingerprint.list_available(show_hidden=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "When you create an interaction fingerprint with `prolif.Fingerprint`, those classes are\n",
+    "created with default parameters and attached to the new fingerprint object as methods:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp = plf.Fingerprint([\"Hydrophobic\", \"HBDonor\", \"HBAcceptor\"])\n",
+    "fp.hydrophobic"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "These method will yield all occurences for that specific interaction between 2 residues (1 for the ligand and 1 for the protein).\n",
+    "\n",
+    "In this example, we'll reparametrize the hydrophobic interaction with a shorter\n",
+    "distance and see how this affects the number of occurences of an interaction for a\n",
+    "given pair of residues. \n",
+    "\n",
+    ":::{note}\n",
+    "To know which parameters are available for interactions, see the\n",
+    "{mod}`prolif.interactions` module.\n",
+    ":::\n",
+    "\n",
+    "Let's start with a test case. With the default parameters, TYR109 is interacting with our ligand in 4 different occasions:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp = plf.Fingerprint()\n",
+    "hydrophobic_interactions_tyr109 = list(\n",
+    "    fp.hydrophobic(ligand_mol[0], protein_mol[\"TYR109.A\"])\n",
+    ")\n",
+    "len(hydrophobic_interactions_tyr109)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we'll simply change the distance threshold to `4.0` instead of the default `4.5`. To do this we simply provide the `parameters` argument with a dictionnary mapping the name of the interaction to reconfigure with updated parameters.\n",
+    "\n",
+    "This modification will affect all the interaction analysis code run by the fingerprint object, i.e. when using\n",
+    "`run` or `run_from_iterable`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp = plf.Fingerprint(parameters={\"Hydrophobic\": {\"distance\": 4.0}})\n",
+    "hydrophobic_interactions_tyr109 = list(\n",
+    "    fp.hydrophobic(ligand_mol[0], protein_mol[\"TYR109.A\"])\n",
+    ")\n",
+    "len(hydrophobic_interactions_tyr109)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As you can see, by reducing the distance threshold the interaction between our ligand\n",
+    "and TYR109 is now ignored."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Repurposing an existing interaction\n",
+    "\n",
+    "In case you want to reuse an existing class for a different type of interaction, the\n",
+    "easiest way is to define an interaction class that inherits one of the classes listed\n",
+    "in the {mod}`prolif.interactions` module, and just update its `__init__` method with the\n",
+    "appropriate parameters.\n",
+    "\n",
+    "There are some generic interactions, like the {class}`~prolif.interactions.base.Distance`\n",
+    "class, if you just need to define two chemical moieties within a certain distance. Both\n",
+    "the Hydrophobic, Ionic, and Metallic interactions inherit from this class!\n",
+    "\n",
+    "In most cases, defining an interaction only based on a distance is not enough and\n",
+    "requires one or two angles constraints as well. For this purpose, the\n",
+    "{class}`~prolif.interactions.base.SingleAngle` and {class}`~prolif.interactions.base.DoubleAngle`\n",
+    "interactions can be used.\n",
+    "\n",
+    "Here we'll define a `C-H...O` HBond interaction by reusing the existing HBond classes:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class CHOAcceptor(plf.interactions.HBAcceptor):\n",
+    "    def __init__(\n",
+    "        self,\n",
+    "        acceptor=\"O\",\n",
+    "        donor=\"[#6]-[H]\",\n",
+    "        distance=3.5,\n",
+    "        DHA_angle=(90, 180),\n",
+    "    ):\n",
+    "        super().__init__(\n",
+    "            acceptor=acceptor, donor=donor, distance=distance, DHA_angle=DHA_angle\n",
+    "        )\n",
+    "\n",
+    "\n",
+    "# create inverse interaction as well\n",
+    "CHODonor = CHOAcceptor.invert_role(\n",
+    "    \"CHODonor\",\n",
+    "    \"C-H...O Hbond interaction between a ligand (donor) and a residue (acceptor)\",\n",
+    ")\n",
+    "\n",
+    "# calculate both classical and weak hbonds\n",
+    "fp = plf.Fingerprint([\"HBAcceptor\", \"CHOAcceptor\", \"HBDonor\", \"CHODonor\"])\n",
+    "fp.run_from_iterable(pose_iterable, protein_mol)\n",
+    "# show dataframe\n",
+    "df = fp.to_dataframe()\n",
+    "df"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can also display these new interactions:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from prolif.plotting.complex3d import Complex3D\n",
+    "\n",
+    "# assign colors for the new interactions on the 3D plot\n",
+    "Complex3D.COLORS.update(\n",
+    "    {\n",
+    "        \"CHODonor\": \"red\",\n",
+    "        \"CHOAcceptor\": \"blue\",\n",
+    "    }\n",
+    ")\n",
+    "# show specific docking pose\n",
+    "pose_index = 4\n",
+    "fp.plot_3d(pose_iterable[pose_index], protein_mol, frame=pose_index)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Defining a custom interaction\n",
+    "\n",
+    ":::{important}\n",
+    "Before you dive into this section, make sure that there isn't already an interaction\n",
+    "that could just be **repurposed** to do what you want!\n",
+    "\n",
+    "For this, the best is to check the {mod}`prolif.interactions` module.\n",
+    ":::\n",
+    "\n",
+    "With that being said, there are a few rules that you must respect when writing your own\n",
+    "interaction:\n",
+    "\n",
+    "#### Inherit the ProLIF Interaction class\n",
+    "  \n",
+    "The {class}`~prolif.interactions.base.Interaction` class is the base class that provides\n",
+    "some functionalities to automatically register the interactions in `Fingerprint` objects.\n",
+    "\n",
+    "#### Naming convention\n",
+    "\n",
+    "For non-symmetrical interactions, like hydrogen bonds or salt-bridges, the convention\n",
+    "used here is to name the class after the role of the ligand.\n",
+    "\n",
+    "For example, the class `HBDonor` detects if a ligand acts as a hydrogen bond donor, and\n",
+    "the class `Cationic` detects if a ligand acts as a cation.\n",
+    "\n",
+    "#### Define a `detect` method\n",
+    "\n",
+    "This method takes exactly two positional arguments: a ligand Residue and a protein\n",
+    "Residue (in this order).\n",
+    "\n",
+    "#### Return value for the `detect` method\n",
+    "\n",
+    "You must yield a dictionary containing some basic metadata about the interaction when it\n",
+    "is detected. To help with this process, the `metadata` method should be used (see\n",
+    "example below) for which the arguments are listed here:\n",
+    "\n",
+    "* the input residues (`lig_res` and `prot_res` arguments, of type `rdkit.Chem.Mol`),\n",
+    "* the indices of atoms responsible for the interaction, (`lig_indices` and\n",
+    "  `prot_indices` arguments, of type `tuple[int, ...]`),\n",
+    "* any other relevant metric (distances or angles), named as you want. Distances should\n",
+    "  be in angstroms, and preferably named `distance`, and angles should be in degrees.\n",
+    "\n",
+    ":::{note}\n",
+    "You don't have to return anything if no interaction is detected for a pair of residues.\n",
+    ":::\n",
+    "\n",
+    "Here's an example implementing a close-contact interaction using numpy:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from scipy.spatial import distance_matrix\n",
+    "\n",
+    "\n",
+    "class CloseContact(plf.interactions.Interaction):\n",
+    "    def __init__(self, contact_threshold=2.0):\n",
+    "        self.contact_threshold = contact_threshold\n",
+    "\n",
+    "    def detect(self, ligand_residue, protein_residue):\n",
+    "        # distance matrix between atoms of both residues\n",
+    "        dist_matrix = distance_matrix(ligand_residue.xyz, protein_residue.xyz)\n",
+    "        below_threshold = dist_matrix < self.contact_threshold\n",
+    "        for ligand_indices in np.argwhere(below_threshold.any(axis=1)):\n",
+    "            ligand_index = int(ligand_indices[0])\n",
+    "            for protein_indices in np.argwhere(below_threshold[ligand_index]):\n",
+    "                protein_index = int(protein_indices[0])\n",
+    "                # yield dict with metadata on the interaction\n",
+    "                # required arguments: input residues, and tuple of indices of atoms\n",
+    "                #                     responsible for the interaction\n",
+    "                # optional arguments: any additional `key=value` pair (e.g. distance)\n",
+    "                yield self.metadata(\n",
+    "                    lig_res=ligand_residue,\n",
+    "                    prot_res=protein_residue,\n",
+    "                    lig_indices=(ligand_index,),\n",
+    "                    prot_indices=(protein_index,),\n",
+    "                    distance=dist_matrix[ligand_index, protein_index],\n",
+    "                )\n",
+    "\n",
+    "\n",
+    "# run analysis\n",
+    "fp = plf.Fingerprint([\"CloseContact\"])\n",
+    "fp.run_from_iterable([ligand_mol], protein_mol)\n",
+    "\n",
+    "# show results\n",
+    "df = fp.to_dataframe()\n",
+    "df"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# assign colors for the new interactions on the 3D plot\n",
+    "Complex3D.COLORS.update(\n",
+    "    {\n",
+    "        \"CloseContact\": \"brown\",\n",
+    "    }\n",
+    ")\n",
+    "# display\n",
+    "fp.plot_3d(ligand_mol, protein_mol, frame=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Fingerprint generation\n",
+    "\n",
+    "This section contains information about modifying some aspects of the fingerprint\n",
+    "generation.\n",
+    "\n",
+    "### Ignoring the protein backbone\n",
+    "\n",
+    "In some cases, you might want to dismiss backbone interactions. You cannot simply remove\n",
+    "the backbone from your protein input file(s), as it will either result in charges added\n",
+    "on the side-chains' end, or you would need to add dummy atoms at the end, but these\n",
+    "could also result in artifacts during the interaction detection.\n",
+    "\n",
+    "One workaround is to use a substructure search (SMARTS) to delete the backbone atoms\n",
+    "after the structures has been parsed by MDAnalysis."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from rdkit import Chem\n",
+    "from rdkit.Chem.AllChem import DeleteSubstructs\n",
+    "\n",
+    "# SMARTS for backbone\n",
+    "backbone_smarts = \"[C^2](=O)-[C;X4](-[H])-[N;+0]-[H]\"\n",
+    "backbone_query = Chem.MolFromSmarts(backbone_smarts)\n",
+    "\n",
+    "\n",
+    "def make_mol_and_strip_backbone(atomgroup, **kwargs):\n",
+    "    mol = atomgroup.convert_to.rdkit(**kwargs)\n",
+    "    mol = DeleteSubstructs(mol, backbone_query)\n",
+    "    return plf.Molecule(mol)\n",
+    "\n",
+    "\n",
+    "# patch the `from_mda` method with our modified version\n",
+    "plf.Molecule.from_mda = make_mol_and_strip_backbone"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can then prepare your system and run the analysis as you normally would."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "prolif",
+   "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.9.13"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/notebooks/docking.ipynb b/docs/notebooks/docking.ipynb
new file mode 100644
index 0000000..c4cee2b
--- /dev/null
+++ b/docs/notebooks/docking.ipynb
@@ -0,0 +1,827 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Docking\n",
+    "\n",
+    "This tutorial showcases how to use ProLIF to generate an interaction fingerprint for interactions between a protein and different docking poses.\n",
+    "\n",
+    "ProLIF currently provides file readers for docking poses for MOL2, SDF and PDBQT files which rely on RDKit or MDAnalysis (or both).\n",
+    "\n",
+    ":::{important}\n",
+    "For convenience, the different files for this tutorial are included with the ProLIF installation, and we'll access these files through `plf.datafiles.datapath` which is a {class}`pathlib.Path` object pointing to the tutorials data directory. This makes it easier to manipulate paths to files, match filenames using wildcards...etc. in a Pythonic way, but you can also use plain strings, e.g. `\"/home/user/prolif/data/vina/rec.pdb\"` instead of `plf.datafiles.datapath / \"vina\" / \"rec.pdb\"` if you prefer.\n",
+    "\n",
+    "**Remember to replace any reference to `plf.datafiles.datapath` with the actual paths to your inputs outside of this tutorial.**\n",
+    ":::\n",
+    "\n",
+    ":::{tip}\n",
+    "At the top of the page you can find links to either download this notebook or run it in\n",
+    "Google Colab. You can install the dependencies for the tutorials with the command:\n",
+    "\n",
+    "```shell\n",
+    "pip install prolif[tutorials]\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Protein preparation\n",
+    "\n",
+    "Let's start by importing MDAnalysis and ProLIF to read our protein for this tutorial.\n",
+    "\n",
+    "In this tutorial we'll be using a PDB file but you can also use the MOL2 format as shown later on.\n",
+    "\n",
+    ":::{important}\n",
+    "No matter which file format you chose to use for the protein, it **must contain explicit hydrogens**.\n",
+    ":::\n",
+    "\n",
+    "### PDB file\n",
+    "\n",
+    "We have 2 options to read the protein file: MDAnalysis (preferred) or RDKit (faster but risky, see below).\n",
+    "\n",
+    ":::{important}\n",
+    "While RDKit is going to be much faster at parsing the file, it can only recognize standard residue names, so some protonated residues may be incorrectly parsed.\n",
+    "\n",
+    "For example, the protonated histidine residue HSE347 in our protein is not correcltly parsed which removes aromaticity on the ring, meaning that ProLIF will not be able to match this residue for pi-stacking interactions.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import prolif as plf\n",
+    "from rdkit import Chem\n",
+    "\n",
+    "protein_file = str(plf.datafiles.datapath / \"vina\" / \"rec.pdb\")\n",
+    "\n",
+    "rdkit_prot = Chem.MolFromPDBFile(protein_file, removeHs=False)\n",
+    "protein_mol = plf.Molecule(rdkit_prot)\n",
+    "# histidine HSE347 not recognized by RDKit\n",
+    "plf.display_residues(protein_mol, slice(260, 263))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Unless you know for sure that all the residue names in your PDB file are standard, it is thus recommanded that you use MDAnalysis for parsing the protein."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import MDAnalysis as mda\n",
+    "import prolif as plf\n",
+    "\n",
+    "u = mda.Universe(protein_file)\n",
+    "protein_mol = plf.Molecule.from_mda(u)\n",
+    "plf.display_residues(protein_mol, slice(260, 263))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{admonition} Troubleshooting\n",
+    "\n",
+    "In some cases, some atomic clashes may be incorrectly classified as bonds and will prevent the conversion of the MDAnalysis molecule to RDKit through ProLIF. Since MDAnalysis uses van der Waals radii for bond detection, one can modify the default radii that are used:\n",
+    "\n",
+    "```python\n",
+    "u.atoms.guess_bonds(vdwradii={\"H\": 1.05, \"O\": 1.48})\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### MOL2 file\n",
+    "\n",
+    "You could also use a MOL2 file for the protein, here's a code snippet to guide you:\n",
+    "\n",
+    "```python\n",
+    "u = mda.Universe(\"protein.mol2\")\n",
+    "# add \"elements\" category\n",
+    "elements = mda.topology.guessers.guess_types(u.atoms.names)\n",
+    "u.add_TopologyAttr(\"elements\", elements)\n",
+    "# create protein mol\n",
+    "protein_mol = plf.Molecule.from_mda(u)\n",
+    "```\n",
+    "\n",
+    ":::{admonition} Troubleshooting\n",
+    "\n",
+    "While doing so, you may run into one of these errors:\n",
+    "- **`RDKit ERROR: Can't kekulize mol. Unkekulized atoms`**\n",
+    "- **`RDKit ERROR: non-ring atom marked aromatic`**\n",
+    "\n",
+    "This usually happens when some of the bonds in the MOL2 file are unconventional. For example in MOE, charged histidines are represented part with aromatic bonds and part with single and double bonds, presumably to capture the different charged resonance structures in a single one. A practical workaround for this is to redefine problematic bonds as single bonds in the `Universe` object:\n",
+    "\n",
+    "```python\n",
+    "u = mda.Universe(\"protein.mol2\")\n",
+    "# replace aromatic bonds with single bonds\n",
+    "for i, bond_order in enumerate(u._topology.bonds.order):\n",
+    "    # you may need to replace double bonds (\"2\") as well\n",
+    "    if bond_order == \"ar\":\n",
+    "        u._topology.bonds.order[i] = 1\n",
+    "# clear the bond cache, just in case\n",
+    "u._topology.bonds._cache.pop(\"bd\", None)\n",
+    "# infer bond orders again\n",
+    "protein_mol = plf.Molecule.from_mda(u)\n",
+    "```\n",
+    ":::\n",
+    "\n",
+    "The parsing of residue info in MOL2 files can sometimes fail and lead to the residue index being appended to the residue name and number. You can fix this with the following snippet:\n",
+    "\n",
+    "```python\n",
+    "import numpy as np\n",
+    "\n",
+    "u = mda.Universe(\"protein.mol2\")\n",
+    "resids = [\n",
+    "    plf.ResidueId.from_string(x) for x in u.residues.resnames\n",
+    "]\n",
+    "u.residues.resnames = np.array([x.name for x in resids], dtype=object)\n",
+    "u.residues.resids = np.array([x.number for x in resids], dtype=np.uint32)\n",
+    "u.residues.resnums = u.residues.resids\n",
+    "protein_mol = plf.Molecule.from_mda(u)\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Docking poses preparation\n",
+    "\n",
+    "Loading docking poses is done through one of the `supplier` functions available in ProLIF. These will read your input file with either RDKit or MDAnalysis and handle any additional preparation.\n",
+    "\n",
+    "### SDF format\n",
+    "\n",
+    "The SDF format is the easiest way for ProLIF to parse docking poses:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# load ligands\n",
+    "poses_path = str(plf.datafiles.datapath / \"vina\" / \"vina_output.sdf\")\n",
+    "pose_iterable = plf.sdf_supplier(poses_path)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### MOL2 format\n",
+    "\n",
+    "MOL2 is another format that should be easy for ProLIF to parse:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# load ligands\n",
+    "poses_path = str(plf.datafiles.datapath / \"vina\" / \"vina_output.mol2\")\n",
+    "pose_iterable = plf.mol2_supplier(poses_path)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### PDBQT format\n",
+    "\n",
+    "The typical use case here is getting the IFP from AutoDock Vina's output. It requires a few additional steps and informations compared to other formats like MOL2, since the PDBQT format gets rid of most hydrogen atoms and doesn't contain bond order information, which are both needed for ProLIF to work.\n",
+    "\n",
+    ":::{important}\n",
+    "Consider using [Meeko](https://github.com/forlilab/Meeko) to prepare your inputs for Vina as it contains some utilities to [convert the output poses](https://github.com/forlilab/Meeko#examples-using-the-command-line-scripts) to the SDF format which is a lot safer than the solution proposed here.\n",
+    ":::\n",
+    "\n",
+    ":::{tip}\n",
+    "Do not use OpenBabel to convert your PDBQT poses to SDF, it will not be able to guess the bond orders and charges correctly.\n",
+    ":::\n",
+    "\n",
+    "The prerequisites for a successfull usage of ProLIF in this case is having external files that contain bond orders and formal charges for your ligand (like SMILES, SDF or MOL2).\n",
+    "\n",
+    ":::{note}\n",
+    "Please note that your PDBQT input must have a single model per file (this is required by `MDAnalysis`). Splitting a multi-model file can be done using the `vina_split` command-line tool that comes with AutoDock Vina: `vina_split --input vina_output.pdbqt`\n",
+    ":::\n",
+    "\n",
+    "Let's start by loading our \"template\" file with bond orders. It can be a SMILES string, MOL2, SDF file or anything supported by RDKit."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from rdkit import Chem\n",
+    "\n",
+    "template = Chem.MolFromSmiles(\n",
+    "    \"C[NH+]1CC(C(=O)NC2(C)OC3(O)C4CCCN4C(=O)C(Cc4ccccc4)N3C2=O)C=C2c3cccc4[nH]cc(c34)CC21\"\n",
+    ")\n",
+    "template"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, we'll use the PDBQT supplier which loads each file from a list of paths, and assigns bond orders and charges using the template."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# load list of ligands\n",
+    "pdbqt_files = sorted(plf.datafiles.datapath.glob(\"vina/*.pdbqt\"))\n",
+    "pdbqt_files"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pose_iterable = plf.pdbqt_supplier(pdbqt_files, template)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Fingerprint generation\n",
+    "\n",
+    "We can now generate a fingerprint. By default, ProLIF will calculate the following interactions: Hydrophobic, HBDonor, HBAcceptor, PiStacking, Anionic, Cationic, CationPi, PiCation, VdWContact.\n",
+    "You can list all interactions that are available with the following command:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plf.Fingerprint.list_available()"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "The default fingerprint will only keep track of the first group of atoms that satisfied the constraints per interaction type and residue pair.\n",
+    "\n",
+    "If you want to keep track of all possible interactions to generate a count-fingerprint (e.g. when there are two atoms in the ligand that make an HBond-donor interaction with residue X), use `plf.Fingerprint(count=True)`.\n",
+    "This is also quite useful for visualization purposes as you can then display the atom pair that has the shortest distance which will look more accurate.\n",
+    "This fingerprint type is however a bit slower to compute.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# use default interactions\n",
+    "fp = plf.Fingerprint()\n",
+    "# run on your poses\n",
+    "fp.run_from_iterable(pose_iterable, protein_mol)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "The `run_from_iterable` method will automatically select residues that are close to the ligand (6.0 Å) when computing the fingerprint. You can modify the 6.0 Å cutoff by specifying `plf.Fingerprint(vicinity_cutoff=7.0)`, but this is only useful if you decide to change the distance parameters for an interaction class (see in the advanced section of the tutorials).\n",
+    "\n",
+    "Alternatively, you can pass a list of residues like so:\n",
+    "\n",
+    "```python\n",
+    "fp.run_from_iterable(<other parameters>, residues=[\"TYR38.A\", \"ASP129.A\"])\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can save the fingerprint object with `fp.to_pickle` and reload it later with `Fingerprint.from_pickle`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp.to_pickle(\"fingerprint.pkl\")\n",
+    "fp = plf.Fingerprint.from_pickle(\"fingerprint.pkl\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Analysis\n",
+    "\n",
+    "Once the execution is done, you can access the results through `fp.ifp` which is a nested dictionary:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pose_index = 0\n",
+    "ligand_residue = \"LIG1.G\"\n",
+    "protein_residue = \"ASP129.A\"\n",
+    "\n",
+    "fp.ifp[pose_index][(ligand_residue, protein_residue)]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "While this contains all the details about the different interactions that were detected, it's not the easiest thing to digest.\n",
+    "\n",
+    "The best way to analyse our results is to export the interaction fingerprint to a Pandas DataFrame. You can read more about pandas in their\n",
+    "[user_guide](https://pandas.pydata.org/docs/user_guide/index.html)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = fp.to_dataframe(index_col=\"Pose\")\n",
+    "# show only the 5 first poses\n",
+    "df.head(5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can change the type used for the bits if need be:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "df_uint8 = fp.to_dataframe(dtype=np.uint8)\n",
+    "df_uint8.head(5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here are some common pandas snippets to extract useful information from the fingerprint table.\n",
+    "\n",
+    ":::{important}\n",
+    "Make sure to remove the `.head(5)` at the end of the commands to display the results for all the frames.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# hide an interaction type (Hydrophobic)\n",
+    "df.drop(\"Hydrophobic\", level=\"interaction\", axis=1).head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# show only one protein residue (ASP129.A)\n",
+    "df.xs(\"ASP129.A\", level=\"protein\", axis=1).head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# show only an interaction type (PiStacking)\n",
+    "df.xs(\"PiStacking\", level=\"interaction\", axis=1).head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of poses where each interaction is present\n",
+    "(df.mean().sort_values(ascending=False).to_frame(name=\"%\").T * 100)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# same but we regroup all interaction types\n",
+    "(\n",
+    "    df.T.groupby(level=[\"ligand\", \"protein\"])\n",
+    "    .sum()\n",
+    "    .T.astype(bool)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of poses where PiStacking interactions are present, by residue\n",
+    "(\n",
+    "    df.xs(\"PiStacking\", level=\"interaction\", axis=1)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of poses where interactions with PHE330 occur, by interaction type\n",
+    "(\n",
+    "    df.xs(\"PHE330.B\", level=\"protein\", axis=1)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of poses where each interaction type is present\n",
+    "(\n",
+    "    df.T.groupby(level=\"interaction\")\n",
+    "    .sum()\n",
+    "    .T.astype(bool)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 10 residues most frequently interacting with the ligand\n",
+    "(\n",
+    "    df.T.groupby(level=[\"ligand\", \"protein\"])\n",
+    "    .sum()\n",
+    "    .T.astype(bool)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .head(10)\n",
+    "    .to_frame(\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can compute a Tanimoto similarity between poses:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Tanimoto similarity between the first frame and the rest\n",
+    "from rdkit import DataStructs\n",
+    "\n",
+    "bitvectors = fp.to_bitvectors()\n",
+    "tanimoto_sims = DataStructs.BulkTanimotoSimilarity(bitvectors[0], bitvectors)\n",
+    "tanimoto_sims"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Comparing docking poses with a reference ligand\n",
+    "\n",
+    "If you have a reference ligand that you wish to use for comparison with your docking poses, you can do so the following way:\n",
+    "\n",
+    ":::{important}\n",
+    "Just like for the protein, you reference ligand must contain explicit hydrogens. If the reference comes from a PDB file, you could use PyMOL to add hydrogens to it and export the prepared reference.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# load the reference\n",
+    "ref = mda.Universe(plf.datafiles.datapath / \"vina\" / \"lig.pdb\")\n",
+    "ref_mol = plf.Molecule.from_mda(ref)\n",
+    "\n",
+    "# generate IFP for the reference\n",
+    "fp_ref = plf.Fingerprint(fp.interactions)\n",
+    "fp_ref.run_from_iterable([ref_mol], protein_mol)\n",
+    "df_ref = fp_ref.to_dataframe(index_col=\"Pose\")\n",
+    "\n",
+    "# set the \"pose index\" to -1\n",
+    "df_ref.rename(index={0: -1}, inplace=True)\n",
+    "# set the ligand name to be the same as poses\n",
+    "df_ref.rename(columns={str(ref_mol[0].resid): df.columns.levels[0][0]}, inplace=True)\n",
+    "df_ref"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# concatenate both dataframes\n",
+    "import pandas as pd\n",
+    "\n",
+    "df_ref_poses = (\n",
+    "    pd.concat([df_ref, df])\n",
+    "    .fillna(False)\n",
+    "    .sort_index(\n",
+    "        axis=1,\n",
+    "        level=1,\n",
+    "        key=lambda index: [plf.ResidueId.from_string(x) for x in index],\n",
+    "    )\n",
+    ")\n",
+    "df_ref_poses"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can now calculate the Tanimoto similarity between our reference ligand and the docking poses:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bitvectors = plf.to_bitvectors(df_ref_poses)\n",
+    "tanimoto_sims = DataStructs.BulkTanimotoSimilarity(bitvectors[0], bitvectors[1:])\n",
+    "tanimoto_sims"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualisation\n",
+    "\n",
+    "There are a few different options builtin when it comes to visualisation.\n",
+    "\n",
+    "You can start by plotting the interactions over time:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib ipympl\n",
+    "\n",
+    "fp.plot_barcode(xlabel=\"Pose\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "If you uncomment `%matplotlib ipympl` at the top of the above cell, you should be able to see an interactive version of the plot.\n",
+    ":::\n",
+    "\n",
+    "You can also display the interactions in a 2D interactive diagram:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp.plot_lignetwork(pose_iterable[0])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This diagram is interactive, you can:\n",
+    "- zoom and pan,\n",
+    "- move residues around,\n",
+    "- click on the legend to display or hide types of residues or interactions,\n",
+    "- hover an interaction line to display the distance.\n",
+    "\n",
+    ":::{note}\n",
+    "It is not possible to export it as an image, but you can always take a screenshot.\n",
+    ":::\n",
+    "\n",
+    "You can generate 2 types of diagram with this function, controlled by the `kind` argument:\n",
+    "- `frame`: shows a single specific docking pose (specified with `frame`, corresponds to the `Pose` index in the dataframe).\n",
+    "- `aggregate` (default): the interactions from all poses are grouped and displayed. An optional `threshold` parameter controls the minimum frequency required for an interaction to be displayed (default `0.3`, meaning that interactions occuring in less than 30% of poses will be hidden). The width of interactions is linked to the frequency."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Showing a specific pose:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pose_index = 0\n",
+    "fp.plot_lignetwork(pose_iterable[pose_index], kind=\"frame\", frame=pose_index)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Up to now we've only been using the default fingerprint generator, but you can optionally enable the `count` parameter to enumerate all occurences of an interaction (the default fingerprint generator will stop at the first occurence), and then display all of them by specifying `display_all=True`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp_count = plf.Fingerprint(count=True)\n",
+    "fp_count.run_from_iterable(pose_iterable, protein_mol)\n",
+    "fp_count.plot_lignetwork(\n",
+    "    pose_iterable[pose_index], kind=\"frame\", frame=pose_index, display_all=True\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can also visualize this information in 3D:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "view = fp_count.plot_3d(\n",
+    "    pose_iterable[pose_index], protein_mol, frame=pose_index, display_all=False\n",
+    ")\n",
+    "view"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As in the lignetwork plot, you can hover atoms and interactions to display more information.\n",
+    "\n",
+    "The advantage of using a count fingerprint in that case is that it will automatically select the interaction occurence with the shortest distance for a more intuitive visualization.\n",
+    "\n",
+    "Once you're satisfied with the orientation, you can export the view as a PNG image with the following snippet:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from IPython.display import Javascript\n",
+    "\n",
+    "Javascript(\n",
+    "    \"\"\"\n",
+    "    var png = viewer_%s.pngURI()\n",
+    "    var a = document.createElement('a')\n",
+    "    a.href = png\n",
+    "    a.download = \"prolif-3d.png\"\n",
+    "    a.click()\n",
+    "    a.remove()\n",
+    "\"\"\"\n",
+    "    % view.uniqueid\n",
+    ")"
+   ]
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae"
+  },
+  "kernelspec": {
+   "display_name": "Python 3.8.5 ('prolif')",
+   "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.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/notebooks/how-to.ipynb b/docs/notebooks/how-to.ipynb
deleted file mode 100644
index 66e6449..0000000
--- a/docs/notebooks/how-to.ipynb
+++ /dev/null
@@ -1,646 +0,0 @@
-{
- "cells": [
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# How-to\n",
-    "\n",
-    "This notebook serves as a practical guide to common questions users might have.\n",
-    "\n",
-    "**Table of content**\n",
-    "\n",
-    "- [Changing the parameters for an interaction](#Changing-the-parameters-for-an-interaction)\n",
-    "- [Writing your own interaction](#Writing-your-own-interaction)\n",
-    "- [Working with docking poses instead of MD simulations](#Working-with-docking-poses-instead-of-MD-simulations)\n",
-    "- [Using PDBQT files](#Using-PDBQT-files)\n",
-    "\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import MDAnalysis as mda\n",
-    "import prolif as plf\n",
-    "import pandas as pd"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)\n",
-    "lig = u.select_atoms(\"resname LIG\")\n",
-    "# restrict to 7.0 angstroms around the ligand\n",
-    "prot = u.select_atoms(\"protein and byres around 7 group ligand\", ligand=lig)\n",
-    "lmol = plf.Molecule.from_mda(lig)\n",
-    "pmol = plf.Molecule.from_mda(prot)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Changing the parameters for an interaction"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "You can list all the available interactions as follow:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plf.Fingerprint.list_available(show_hidden=True)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "In this example, we'll reparametrize the hydrophobic interaction with a shorter distance. To know which parameters are available for interactions, see the [prolif.interactions](../source/modules/interaction-fingerprint.html#detecting-interactions-between-residues-prolif-interactions-interactions) module.\n",
-    "\n",
-    "Let's start with a test case: with the default parameters, Y109 is interacting with our ligand in 4 different occasions."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fp = plf.Fingerprint()\n",
-    "hydrophobic_interactions = fp.hydrophobic(lmol[\"LIG1.G\"], pmol[\"TYR109.A\"])\n",
-    "len(list(hydrophobic_interactions))"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Now we'll simply change the distance threshold to `4.0` instead of the default `4.5`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fp = plf.Fingerprint(parameters={\"Hydrophobic\": {\"distance\": 4.0}})\n",
-    "hydrophobic_interactions = fp.hydrophobic(lmol[\"LIG1.G\"], pmol[\"TYR109.A\"])\n",
-    "len(list(hydrophobic_interactions))"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "You can then use `fp.run` and other methods as usual."
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Repurposing an existing interaction"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "In case you want to reuse an existing class for a different type of interaction, the easiest way is to \n",
-    "define an interaction class that inherits one of the classes listed in the `prolif.interactions` module, and just update its `__init__` method with the appropriate parameters:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "class ShorterHydrophobic(plf.interactions.Hydrophobic):\n",
-    "    def __init__(self, distance=4.0):\n",
-    "        super().__init__(distance=distance)\n",
-    "\n",
-    "\n",
-    "fp = plf.Fingerprint([\"Hydrophobic\", \"ShorterHydrophobic\"])\n",
-    "hydrophobic_interactions = fp.hydrophobic(lmol[\"LIG1.G\"], pmol[\"TYR109.A\"])\n",
-    "print(len(list(hydrophobic_interactions)))  # 4\n",
-    "shorterhydrophobic_interactions = fp.shorterhydrophobic(\n",
-    "    lmol[\"LIG1.G\"], pmol[\"TYR109.A\"]\n",
-    ")\n",
-    "print(len(list(shorterhydrophobic_interactions)))  # 0"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Writing your own interaction"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Before you dive into this section, make sure that there isn't already an interaction that could just be **repurposed** to do what you want!\n",
-    "\n",
-    "For this, the best is to check the section of the documentation corresponding to the `prolif.interactions` module. There are some generic interactions, like the `Distance` class, if you just need to define two chemical moieties within a certain distance. Both the Hydrophobic, Ionic, and Metallic interactions inherit from this class!\n",
-    "In most cases, defining an interaction only based on a distance is not enough and requires one or two angles constraints as well. For this purpose, the `SingleAngle` and `DoubleAngle` interactions can be used."
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "With that being said, there are a few rules that you must respect when writing your own interaction:\n",
-    "\n",
-    "* **Inherit the ProLIF Interaction class**\n",
-    "  \n",
-    "This `prolif.interactions.Interaction` class in the base class that provides some functionalities to automatically register the interactions in `Fingerprint` objects.\n",
-    "\n",
-    "* **Naming convention**\n",
-    "\n",
-    "For non-symmetrical interactions, like hydrogen bonds or salt-bridges, the convention used here is to name the class after the role of the ligand. For example, the class `HBDonor` detects if a ligand acts as a hydrogen bond donor, and the class `Cationic` detects if a ligand acts as a cation.\n",
-    "\n",
-    "* **Define a `detect` method**\n",
-    "\n",
-    "This method takes exactly two positional arguments: a ligand Residue and a protein Residue (in this order).\n",
-    "\n",
-    "* **Return value for the `detect` method**\n",
-    "\n",
-    "You must yield a dictionary containing some basic metadata about the interaction when it is detected. To help with this process, the `metadata` method should be used (see example below) for which the arguments are listed here:\n",
-    "\n",
-    "  * the input residues (`lig_res` and `prot_res` arguments, of type `rdkit.Chem.Mol`),\n",
-    "  * the indices of atoms responsible for the interaction, (`lig_indices` and `prot_indices` arguments, of type `tuple[int, ...]`),\n",
-    "  * any other relevant metric (distances or angles), named as you want. Distances should be in angstroms, and preferably named `distance`, and angles should be in degrees.\n",
-    "  \n",
-    "\n",
-    "Note that you don't have to return anything if no interaction is detected for a pair of residues."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import numpy as np\n",
-    "from scipy.spatial import distance_matrix\n",
-    "\n",
-    "\n",
-    "class CloseContact(plf.interactions.Interaction):\n",
-    "    def __init__(self, contact_threshold=2.0):\n",
-    "        self.contact_threshold = contact_threshold\n",
-    "\n",
-    "    def detect(self, ligand_residue, protein_residue):\n",
-    "        # distance matrix between atoms of both residues\n",
-    "        dist_matrix = distance_matrix(ligand_residue.xyz, protein_residue.xyz)\n",
-    "        below_threshold = dist_matrix < self.contact_threshold\n",
-    "        for ligand_indices in np.argwhere(below_threshold.any(axis=1)):\n",
-    "            ligand_index = int(ligand_indices[0])\n",
-    "            for protein_indices in np.argwhere(below_threshold[ligand_index]):\n",
-    "                protein_index = int(protein_indices[0])\n",
-    "                # yield dict with metadata on the interaction\n",
-    "                # required arguments: input residues, and tuple of indices of atoms\n",
-    "                #                     responsible for the interaction\n",
-    "                # optional arguments: any additional `key=value` pair (e.g. distance)\n",
-    "                yield self.metadata(\n",
-    "                    lig_res=ligand_residue,\n",
-    "                    prot_res=protein_residue,\n",
-    "                    lig_indices=(ligand_index,),\n",
-    "                    prot_indices=(protein_index,),\n",
-    "                    distance=dist_matrix[ligand_index, protein_index],\n",
-    "                )\n",
-    "\n",
-    "\n",
-    "fp = plf.Fingerprint([\"CloseContact\"])\n",
-    "len(list(fp.closecontact(lmol[\"LIG1.G\"], pmol[\"ASP129.A\"])))"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "If you want to access metadata about the interaction, you can do one of the following:"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "* Recommended: use the `run` or `run_from_iterable` methods and then interact with the `fp.ifp` attribute:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fp = plf.Fingerprint()\n",
-    "fp.run(u.trajectory[0:1], lig, prot)\n",
-    "# check the interactions between the ligand and ASP129 on the first frame\n",
-    "frame_number, ligand_residue, protein_residue = 0, \"LIG1.G\", \"ASP129.A\"\n",
-    "fp.ifp[frame_number][(ligand_residue, protein_residue)]"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "* For custom workflows: use the `metadata=True` argument when calling the `generate` method:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "ifp = fp.generate(lmol, pmol, metadata=True)\n",
-    "# check the interaction between the ligand and ASP129\n",
-    "ifp[(\"LIG1.G\", \"ASP129.A\")]"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "* For debugging: use the `metadata=True` argument when calling the interaction directly through the fingerprint object:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "list(fp.closecontact(lmol[\"LIG1.G\"], pmol[\"ASP129.A\"], metadata=True))"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "* For debugging: use the `metadata` method instead of `bitvector`:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fp.metadata(lmol[\"LIG1.G\"], pmol[\"ASP129.A\"])"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Working with docking poses instead of MD simulations\n",
-    "\n",
-    "ProLIF currently provides file readers for MOL2, SDF and PDBQT files. The API is slightly different compared to the quickstart example but the end result is the same.\n",
-    "\n",
-    "Please note that this part of the tutorial is only suitable for interactions between one protein and several ligands, or in more general terms, between one molecule with multiple residues and one molecule with a single residue. This is not suitable for protein-protein or DNA-protein interactions.\n",
-    "\n",
-    "Let's start by loading the protein. Here I'm using a PDB file but you can use any format supported by MDAnalysis as long as it contains explicit hydrogens.\n",
-    "\n",
-    "Note that for this tutorial, we're using example files that come with the package. These files are accessed through the `plf.datafiles.datapath` variable which holds a `pathlib.Path` object. This makes it easier to manipulate paths to file, match filenames using wildcards...etc. in a Pythonic way, but you can also use plain strings, i.e. `\"/home/cedric/projects/ProLIF/prolif/data/vina/rec.pdb\"` instead of `plf.datafiles.datapath / \"vina\" / \"rec.pdb\"` if you prefer."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# load protein\n",
-    "prot = mda.Universe(plf.datafiles.datapath / \"vina\" / \"rec.pdb\")\n",
-    "prot = plf.Molecule.from_mda(prot)\n",
-    "prot.n_residues"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Using an SDF file"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# load ligands\n",
-    "path = str(plf.datafiles.datapath / \"vina\" / \"vina_output.sdf\")\n",
-    "lig_suppl = plf.sdf_supplier(path)\n",
-    "# generate fingerprint\n",
-    "fp = plf.Fingerprint()\n",
-    "fp.run_from_iterable(lig_suppl, prot)\n",
-    "df = fp.to_dataframe()\n",
-    "df"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "If you want to calculate the Tanimoto similarity between your docked poses and a reference ligand, here's how to do it.\n",
-    "\n",
-    "We first need to generate the interaction fingerprint for the reference, and concatenate it to the previous one"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# load the reference\n",
-    "ref = mda.Universe(plf.datafiles.datapath / \"vina\" / \"lig.pdb\")\n",
-    "ref = plf.Molecule.from_mda(ref)\n",
-    "# generate IFP\n",
-    "fp.run_from_iterable([ref], prot)\n",
-    "df0 = fp.to_dataframe()\n",
-    "df0.rename({0: \"ref\"}, inplace=True)\n",
-    "# drop the ligand level on both dataframes\n",
-    "df0.columns = df0.columns.droplevel(0)\n",
-    "df.columns = df.columns.droplevel(0)\n",
-    "# concatenate and sort columns\n",
-    "df = (\n",
-    "    pd.concat([df0, df])\n",
-    "    .fillna(False)\n",
-    "    .sort_index(\n",
-    "        axis=1, level=0, key=lambda index: [plf.ResidueId.from_string(x) for x in index]\n",
-    "    )\n",
-    ")\n",
-    "df"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Lastly, we can convert the dataframe to a list of RDKit bitvectors to finally compute the Tanimoto similarity between our reference pose and the docking poses generated by Vina:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from rdkit import DataStructs\n",
-    "\n",
-    "bvs = plf.to_bitvectors(df)\n",
-    "for i, bv in enumerate(bvs[1:]):\n",
-    "    tc = DataStructs.TanimotoSimilarity(bvs[0], bv)\n",
-    "    print(f\"{i}: {tc:.3f}\")"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Interestingly, the best scored docking pose (#0) isn't the most similar to the reference (#5)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Using a MOL2 file\n",
-    "\n",
-    "The input mol2 file can contain multiple ligands in different conformations."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# load ligands\n",
-    "path = str(plf.datafiles.datapath / \"vina\" / \"vina_output.mol2\")\n",
-    "lig_suppl = plf.mol2_supplier(path)\n",
-    "# generate fingerprint\n",
-    "fp = plf.Fingerprint()\n",
-    "fp.run_from_iterable(lig_suppl, prot)\n",
-    "df = fp.to_dataframe()\n",
-    "df"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "If your protein is also a MOL2 file, here's a code snippet to guide you:\n",
-    "\n",
-    "```python\n",
-    "u = mda.Universe(\"protein.mol2\")\n",
-    "# add \"elements\" category\n",
-    "elements = mda.topology.guessers.guess_types(u.atoms.names)\n",
-    "u.add_TopologyAttr(\"elements\", elements)\n",
-    "# create protein mol and run\n",
-    "prot = plf.Molecule.from_mda(u)\n",
-    "fp = plf.Fingerprint()\n",
-    "suppl = plf.mol2_supplier(\"ligands.mol2\")\n",
-    "fp.run_from_iterable(suppl, prot)\n",
-    "df = fp.to_dataframe()\n",
-    "df\n",
-    "```"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "\n",
-    "While doing so, you may run into these errors:\n",
-    "\n",
-    "* **`RDKit ERROR: [17:50:37] Can't kekulize mol.  Unkekulized atoms: 1123`**\n",
-    "  This usually happens when some of the bonds in the MOL2 file are unconventional. For example, charged histidines are represented part with aromatic bonds and part with single and double bonds in MOE, presumably to capture the different charged resonance structures in a single one. A practical workaround for this is to redefine problematic bonds as single bonds in the Universe:\n",
-    "\n",
-    "```python\n",
-    "u = mda.Universe(\"protein.mol2\")\n",
-    "# replace aromatic bonds with single bonds\n",
-    "for i, bond_order in enumerate(u._topology.bonds.order):\n",
-    "    # you may need to replace double bonds (\"2\") as well\n",
-    "    if bond_order == \"ar\":\n",
-    "        u._topology.bonds.order[i] = 1\n",
-    "# clear the bond cache, just in case\n",
-    "u._topology.bonds._cache.pop(\"bd\", None)\n",
-    "# infer bond orders again\n",
-    "prot = plf.Molecule.from_mda(u)\n",
-    "```"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "\n",
-    "* **`RDKit ERROR: [17:50:37] non-ring atom 33 marked aromatic`**\n",
-    "  This is very similar to the previous error. You can use the same fix or, as an alternative, use RDKit:\n",
-    "  \n",
-    "```python\n",
-    "from rdkit import Chem\n",
-    "mol = Chem.MolFromMol2File(\"protein.mol2\", removeHs=False)\n",
-    "# assign residue info (needed for fingerprint generation)\n",
-    "u = mda.Universe(\"protein.mol2\")\n",
-    "for atom, resname in zip(mol.GetAtoms(), u.atoms.resnames):\n",
-    "    resid = plf.ResidueId.from_string(resname)\n",
-    "    mi = Chem.AtomPDBResidueInfo()\n",
-    "    mi.SetResidueNumber(resid.number)\n",
-    "    mi.SetResidueName(resid.name)\n",
-    "    atom.SetMonomerInfo(mi)\n",
-    "prot = plf.Molecule(mol)\n",
-    "```"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "\n",
-    "* **Residue naming** in the resulting dataframe may be different from what was expected as the residue index is appended to the residue name and number:\n",
-    "\n",
-    "```python\n",
-    "import numpy as np\n",
-    "u = mda.Universe(\"protein.mol2\")\n",
-    "resids = [plf.ResidueId.from_string(x) for x in u.residues.resnames]\n",
-    "u.residues.resnames = np.array([x.name for x in resids], dtype=object)\n",
-    "u.residues.resids = np.array([x.number for x in resids], dtype=np.uint32)\n",
-    "u.residues.resnums = u.residues.resids\n",
-    "prot = plf.Molecule.from_mda(u)\n",
-    "```"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Using PDBQT files\n",
-    "\n",
-    "The typical use case here is getting the IFP from AutoDock Vina's output. It requires a few additional steps and informations compared to other formats like MOL2, since the PDBQT format gets rid of most hydrogen atoms and doesn't contain bond order information.\n",
-    "\n",
-    "The prerequisites for a successfull usage of ProLIF in this case is having external files that contain bond orders and formal charges for your ligand (like SMILES, SDF or MOL2), or at least a file with explicit hydrogen atoms. \n",
-    "\n",
-    "Please note that your PDBQT input must have a single model per file (this is required by MDAnalysis). Splitting a multi-model file can be done using the `vina_split` command-line tool that comes with AutoDock Vina: `vina_split --input vina_output.pdbqt`\n",
-    "\n",
-    "Let's start by loading our \"template\" file with bond orders. It can be a SMILES string, MOL2, SDF file or anything supported by RDKit."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from rdkit import Chem\n",
-    "\n",
-    "template = Chem.MolFromSmiles(\n",
-    "    \"C[NH+]1CC(C(=O)NC2(C)OC3(O)C4CCCN4C(=O)C(Cc4ccccc4)N3C2=O)C=C2c3cccc4[nH]cc(c34)CC21\"\n",
-    ")\n",
-    "template"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Next, we'll use the PDBQT supplier which loads each file from a list of paths, and assigns bond orders and charges using the template. The template and PDBQT file must have the exact same atoms, **even hydrogens**, otherwise no match will be found. Since PDBQT files partially keep the hydrogen atoms, we have the choice between:\n",
-    "\n",
-    "* Manually selecting where to add the hydrogens on the template, do the matching, then add the remaining hydrogens (not covered here)\n",
-    "* Or just remove the hydrogens from the PDBQT file, do the matching, then re-add the original hydrogens.\n",
-    "\n",
-    "For the protein, there's usually no need to load the PDBQT that was used by Vina. The original file that was used to generate the PDBQT can be used directly, but **it must contain explicit hydrogen atoms**:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# load list of ligands\n",
-    "pdbqt_files = sorted(plf.datafiles.datapath.glob(\"vina/*.pdbqt\"))\n",
-    "lig_suppl = plf.pdbqt_supplier(pdbqt_files, template)\n",
-    "# generate fingerprint\n",
-    "fp = plf.Fingerprint()\n",
-    "fp.run_from_iterable(lig_suppl, prot)\n",
-    "df = fp.to_dataframe()\n",
-    "df"
-   ]
-  }
- ],
- "metadata": {
-  "interpreter": {
-   "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae"
-  },
-  "kernelspec": {
-   "display_name": "Python 3.8.13 ('prolif')",
-   "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.9.13"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/notebooks/md-ligand-protein.ipynb b/docs/notebooks/md-ligand-protein.ipynb
new file mode 100644
index 0000000..03671c1
--- /dev/null
+++ b/docs/notebooks/md-ligand-protein.ipynb
@@ -0,0 +1,701 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Ligand-protein MD\n",
+    "\n",
+    "This tutorial showcases how to use ProLIF to generate an interaction fingerprint for a ligand-protein complex from an MD simulation.\n",
+    "\n",
+    "ProLIF uses MDAnalysis to process MD simulations, as it supports many [file formats](https://userguide.mdanalysis.org/stable/formats/index.html) and is inter-operable with RDKit (which ProLIF uses under the hood to implement the different types of interactions). To learn more on how to use MDAnalysis, you can find their user guide [here](https://userguide.mdanalysis.org/stable/index.html).\n",
+    "\n",
+    "We will use 3 objects from the MDAnalysis library:\n",
+    "\n",
+    "- The `Universe` which bundles the atoms and bonds of your system with the coordinates in your trajectory.\n",
+    "- The `AtomGroup` which is a collection of atoms that you can define by applying a [selection](https://userguide.mdanalysis.org/stable/selections.html) on the `Universe`.\n",
+    "- The trajectory (or most often a subset of it) to know which frames to process.\n",
+    "\n",
+    ":::{important}\n",
+    "For convenience, the topology and trajectory files for this tutorial are included with the ProLIF installation, and you can access the path to both files through `prolif.datafiles.TOP` and `prolif.datafiles.TRAJ` respectively. **Remember to switch these with the actual paths to your inputs outside of this tutorial.**\n",
+    ":::\n",
+    "\n",
+    ":::{tip}\n",
+    "At the top of the page you can find links to either download this notebook or run it in\n",
+    "Google Colab. You can install the dependencies for the tutorials with the command:\n",
+    "\n",
+    "```shell\n",
+    "pip install prolif[tutorials]\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Preparation\n",
+    "\n",
+    "Let's start by importing MDAnalysis and ProLIF to read our tutorial files, and create selections for the ligand and protein:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import MDAnalysis as mda\n",
+    "import prolif as plf\n",
+    "\n",
+    "# load topology and trajectory\n",
+    "u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)\n",
+    "\n",
+    "# create selections for the ligand and protein\n",
+    "ligand_selection = u.select_atoms(\"resname LIG\")\n",
+    "protein_selection = u.select_atoms(\"protein\")\n",
+    "ligand_selection, protein_selection"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "MDAnalysis should automatically recognize the file type that you're using from its extension. \n",
+    "\n",
+    ":::{note}\n",
+    "Click [here](https://userguide.mdanalysis.org/stable/examples/quickstart.html) to learn more about loading files with MDAnalysis,\n",
+    "and [here](https://userguide.mdanalysis.org/stable/selections.html) to learn more about their atom selection language.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In our case the protein is of reasonable size, but if you're working with a very large system, to save some time and memory it may be wise to restrict the protein selection to a sphere around the ligand:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "protein_selection = u.select_atoms(\n",
+    "    \"protein and byres around 20.0 group ligand\", ligand=ligand_selection\n",
+    ")\n",
+    "protein_selection"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Depending on your system, it may be of interest to also include water molecules in the protein selection. There are none in this tutorial example but something like this could be used:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "protein_selection = u.select_atoms(\n",
+    "    \"(protein or resname WAT) and byres around 20.0 group ligand\",\n",
+    "    ligand=ligand_selection,\n",
+    ")\n",
+    "protein_selection"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, lets make sure that our ligand was correctly read by MDAnalysis.\n",
+    "\n",
+    ":::{important}\n",
+    "This next step is crucial if you're loading a structure from a file that doesn't\n",
+    "explicitely contain bond orders and formal charges, such as a PDB file or most MD\n",
+    "trajectory files. MDAnalysis will infer those from the atoms connectivity, which\n",
+    "**requires all atoms including hydrogens to be present in the input file**.\n",
+    ":::\n",
+    "\n",
+    "Since ProLIF molecules are built on top of RDKit, we can use RDKit functions to display molecules. Let's have a quick look at our ligand:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create a molecule from the MDAnalysis selection\n",
+    "ligand_mol = plf.Molecule.from_mda(ligand_selection)\n",
+    "# display\n",
+    "plf.display_residues(ligand_mol, size=(400, 200))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can do the same for the residues in the protein (only showing the first 20 to keep the notebook short):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "protein_mol = plf.Molecule.from_mda(protein_selection)\n",
+    "plf.display_residues(protein_mol, slice(20))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{admonition} Troubleshooting\n",
+    "\n",
+    "If one of the two molecules was not processed correctly, it might be because your input file does not contain bond information. If that's the case, please add this snippet right after creating your selections:\n",
+    "\n",
+    "```python\n",
+    "ligand_selection.guess_bonds()\n",
+    "protein_selection.guess_bonds()\n",
+    "```\n",
+    "\n",
+    "In some cases, some atomic clashes may be incorrectly classified as bonds and will prevent the conversion of MDAnalysis molecules to RDKit. Since MDAnalysis uses van der Waals radii for bond detection, one can modify the default radii that are used:\n",
+    "\n",
+    "```python\n",
+    "protein_selection.guess_bonds(vdwradii={\"H\": 1.05, \"O\": 1.48})\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Fingerprint generation\n",
+    "\n",
+    "Everything looks good, we can now generate a fingerprint. By default, ProLIF will calculate the following interactions: Hydrophobic, HBDonor, HBAcceptor, PiStacking, Anionic, Cationic, CationPi, PiCation, VdWContact.\n",
+    "You can list all interactions that are available with the following command:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plf.Fingerprint.list_available()"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "The default fingerprint will only keep track of the first group of atoms that satisfied the constraints per interaction type and residue pair.\n",
+    "\n",
+    "If you want to keep track of all possible interactions to generate a count-fingerprint (e.g. when there are two atoms in the ligand that make an HBond-donor interaction with residue X), use `plf.Fingerprint(count=True)`.\n",
+    "This is also quite useful for visualization purposes as you can then display the atom pair that has the shortest distance which will look more accurate.\n",
+    "This fingerprint type is however a bit slower to compute.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# use default interactions\n",
+    "fp = plf.Fingerprint()\n",
+    "# run on a slice of the trajectory frames: from begining to end with a step of 10\n",
+    "fp.run(u.trajectory[::10], ligand_selection, protein_selection)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "The `run` method will automatically select residues that are close to the ligand (6.0 Å) when computing the fingerprint. You can modify the 6.0 Å cutoff by specifying `plf.Fingerprint(vicinity_cutoff=7.0)`, but this is only useful if you decide to change the distance parameters for an interaction class (see in the advanced section of the tutorials).\n",
+    "\n",
+    "Alternatively, you can pass a list of residues like so:\n",
+    "\n",
+    "```python\n",
+    "fp.run(<other parameters>, residues=[\"TYR38.A\", \"ASP129.A\"])\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can save the fingerprint object with `fp.to_pickle` and reload it later with `Fingerprint.from_pickle`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp.to_pickle(\"fingerprint.pkl\")\n",
+    "fp = plf.Fingerprint.from_pickle(\"fingerprint.pkl\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Analysis\n",
+    "\n",
+    "Once the execution is done, you can access the results through `fp.ifp` which is a nested dictionary:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "frame_number = 0\n",
+    "ligand_residue = \"LIG1.G\"\n",
+    "protein_residue = \"ASP129.A\"\n",
+    "\n",
+    "fp.ifp[frame_number][(ligand_residue, protein_residue)]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "While this contains all the details about the different interactions that were detected, it's not the easiest thing to digest.\n",
+    "\n",
+    "The best way to analyse our results is to export the interaction fingerprint to a Pandas DataFrame. You can read more about pandas in their\n",
+    "[user_guide](https://pandas.pydata.org/docs/user_guide/index.html)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = fp.to_dataframe()\n",
+    "# show only the 10 first frames\n",
+    "df.head(10)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can change the type used for the bits if need be:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "df_uint8 = fp.to_dataframe(dtype=np.uint8)\n",
+    "df_uint8.head(5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here are some common pandas snippets to extract useful information from the fingerprint table.\n",
+    "\n",
+    ":::{important}\n",
+    "Make sure to remove the `.head(5)` at the end of the commands to display the results for all the frames.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# hide an interaction type (Hydrophobic)\n",
+    "df.drop(\"Hydrophobic\", level=\"interaction\", axis=1).head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# show only one protein residue (ASP129.A)\n",
+    "df.xs(\"ASP129.A\", level=\"protein\", axis=1).head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# show only an interaction type (PiStacking)\n",
+    "df.xs(\"PiStacking\", level=\"interaction\", axis=1).head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of the trajectory where each interaction is present\n",
+    "(df.mean().sort_values(ascending=False).to_frame(name=\"%\").T * 100)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# same but we regroup all interaction types\n",
+    "(\n",
+    "    df.T.groupby(level=[\"ligand\", \"protein\"])\n",
+    "    .sum()\n",
+    "    .T.astype(bool)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of the trajectory where PiStacking interactions are present, by residue\n",
+    "(\n",
+    "    df.xs(\"PiStacking\", level=\"interaction\", axis=1)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of the trajectory where interactions with SER212 occur, by interaction type\n",
+    "(\n",
+    "    df.xs(\"SER212.A\", level=\"protein\", axis=1)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of the trajectory where each interaction type is present\n",
+    "(\n",
+    "    df.T.groupby(level=\"interaction\")\n",
+    "    .sum()\n",
+    "    .T.astype(bool)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 10 residues most frequently interacting with the ligand\n",
+    "(\n",
+    "    df.T.groupby(level=[\"ligand\", \"protein\"])\n",
+    "    .sum()\n",
+    "    .T.astype(bool)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .head(10)\n",
+    "    .to_frame(\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can compute a Tanimoto similarity between frames:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Tanimoto similarity between the first frame and the rest\n",
+    "from rdkit import DataStructs\n",
+    "\n",
+    "bitvectors = fp.to_bitvectors()\n",
+    "tanimoto_sims = DataStructs.BulkTanimotoSimilarity(bitvectors[0], bitvectors)\n",
+    "tanimoto_sims"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To compare binding modes in your trajectory, it's possible to compute the entire similarity matrix:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pandas as pd\n",
+    "import seaborn as sns\n",
+    "from matplotlib import pyplot as plt\n",
+    "\n",
+    "# Tanimoto similarity matrix\n",
+    "bitvectors = fp.to_bitvectors()\n",
+    "similarity_matrix = []\n",
+    "for bv in bitvectors:\n",
+    "    similarity_matrix.append(DataStructs.BulkTanimotoSimilarity(bv, bitvectors))\n",
+    "similarity_matrix = pd.DataFrame(similarity_matrix, index=df.index, columns=df.index)\n",
+    "\n",
+    "# display heatmap\n",
+    "fig, ax = plt.subplots(figsize=(3, 3), dpi=200)\n",
+    "colormap = sns.diverging_palette(\n",
+    "    300, 145, s=90, l=80, sep=30, center=\"dark\", as_cmap=True\n",
+    ")\n",
+    "sns.heatmap(\n",
+    "    similarity_matrix,\n",
+    "    ax=ax,\n",
+    "    square=True,\n",
+    "    cmap=colormap,\n",
+    "    vmin=0,\n",
+    "    vmax=1,\n",
+    "    center=0.5,\n",
+    "    xticklabels=5,\n",
+    "    yticklabels=5,\n",
+    ")\n",
+    "ax.invert_yaxis()\n",
+    "plt.yticks(rotation=\"horizontal\")\n",
+    "fig.patch.set_facecolor(\"white\")"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualisation\n",
+    "\n",
+    "There are a few different options builtin when it comes to visualisation.\n",
+    "\n",
+    "You can start by plotting the interactions over time:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib ipympl\n",
+    "\n",
+    "fp.plot_barcode()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "If you uncomment `%matplotlib ipympl` at the top of the above cell, you should be able to see an interactive version of the plot.\n",
+    ":::\n",
+    "\n",
+    "You can also display the interactions in a 2D interactive diagram:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp.plot_lignetwork(ligand_mol)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This diagram is interactive, you can:\n",
+    "- zoom and pan,\n",
+    "- move residues around,\n",
+    "- click on the legend to display or hide types of residues or interactions,\n",
+    "- hover an interaction line to display the distance.\n",
+    "\n",
+    ":::{note}\n",
+    "It is not possible to export it as an image, but you can always take a screenshot.\n",
+    ":::\n",
+    "\n",
+    "You can generate 2 types of diagram with this function, controlled by the `kind` argument:\n",
+    "- `frame`: shows a single specific frame (specified with `frame`, corresponds to the frame number in the simulation).\n",
+    "- `aggregate` (default): the interactions from all frames are grouped and displayed. An optional `threshold` parameter controls the minimum frequency required for an interaction to be displayed (default `0.3`, meaning that interactions occuring in less than 30% of frames will be hidden). The width of interactions is linked to the frequency."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Some more examples:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp.plot_lignetwork(ligand_mol, threshold=0.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Up to now we've only been using the default fingerprint generator, but you can optionally enable the `count` parameter to enumerate all occurences of an interaction (the default fingerprint generator will stop at the first occurence), and then display all of them by specifying `display_all=True`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp_count = plf.Fingerprint(count=True)\n",
+    "fp_count.run(u.trajectory[0:1], ligand_selection, protein_selection)\n",
+    "fp_count.plot_lignetwork(ligand_mol, kind=\"frame\", frame=0, display_all=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can also visualize this information in 3D:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "view = fp_count.plot_3d(ligand_mol, protein_mol, frame=0, display_all=False)\n",
+    "view"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As in the lignetwork plot, you can hover atoms and interactions to display more information.\n",
+    "\n",
+    "The advantage of using a count fingerprint in that case is that it will automatically select the interaction occurence with the shortest distance for a more intuitive visualization.\n",
+    "\n",
+    "Once you're satisfied with the orientation, you can export the view as a PNG image with the following snippet:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from IPython.display import Javascript\n",
+    "\n",
+    "Javascript(\n",
+    "    \"\"\"\n",
+    "    var png = viewer_%s.pngURI()\n",
+    "    var a = document.createElement('a')\n",
+    "    a.href = png\n",
+    "    a.download = \"prolif-3d.png\"\n",
+    "    a.click()\n",
+    "    a.remove()\n",
+    "\"\"\"\n",
+    "    % view.uniqueid\n",
+    ")"
+   ]
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae"
+  },
+  "kernelspec": {
+   "display_name": "Python 3.8.5 ('prolif')",
+   "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.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/notebooks/md-protein-protein.ipynb b/docs/notebooks/md-protein-protein.ipynb
new file mode 100644
index 0000000..842767e
--- /dev/null
+++ b/docs/notebooks/md-protein-protein.ipynb
@@ -0,0 +1,830 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Protein-protein MD\n",
+    "\n",
+    "This tutorial showcases how to use ProLIF to generate an interaction fingerprint for protein-protein interactions in an MD simulation.\n",
+    "\n",
+    "ProLIF uses MDAnalysis to process MD simulations, as it supports many [file formats](https://userguide.mdanalysis.org/stable/formats/index.html) and is inter-operable with RDKit (which ProLIF uses under the hood to implement the different types of interactions). To learn more on how to use MDAnalysis, you can find their user guide [here](https://userguide.mdanalysis.org/stable/index.html).\n",
+    "\n",
+    "We will use 3 objects from the MDAnalysis library:\n",
+    "\n",
+    "- The `Universe` which bundles the atoms and bonds of your system with the coordinates in your trajectory.\n",
+    "- The `AtomGroup` which is a collection of atoms that you can define by applying a [selection](https://userguide.mdanalysis.org/stable/selections.html) on the `Universe`.\n",
+    "- The trajectory (or most often a subset of it) to know which frames to process.\n",
+    "\n",
+    ":::{important}\n",
+    "For convenience, the topology and trajectory files for this tutorial are included with the ProLIF installation, and you can access the path to both files through `prolif.datafiles.TOP` and `prolif.datafiles.TRAJ` respectively. **Remember to switch these with the actual paths to your inputs outside of this tutorial.**\n",
+    ":::\n",
+    "\n",
+    ":::{tip}\n",
+    "At the top of the page you can find links to either download this notebook or run it in\n",
+    "Google Colab. You can install the dependencies for the tutorials with the command:\n",
+    "\n",
+    "```shell\n",
+    "pip install prolif[tutorials]\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Preparation\n",
+    "\n",
+    "Let's start by importing MDAnalysis and ProLIF to read our tutorial files, and create selections for the protein components. To keep the size of tutorial files short, we are going to reuse the same simulation as for the ligand-protein system, and decompose our protein in 2 virtual segments: one of the seven transmembrane domains of the GPCR, and the rest of the GPCR. In a real-world scenario this could be for example a peptide and a protein."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import MDAnalysis as mda\n",
+    "import prolif as plf\n",
+    "\n",
+    "# load topology and trajectory\n",
+    "u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)\n",
+    "\n",
+    "# create selections for both protein components\n",
+    "small_protein_selection = u.select_atoms(\"resid 119:152\")\n",
+    "large_protein_selection = u.select_atoms(\n",
+    "    \"protein and not group peptide\", peptide=small_protein_selection\n",
+    ")\n",
+    "small_protein_selection, large_protein_selection"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "MDAnalysis should automatically recognize the file type that you're using from its extension. \n",
+    "\n",
+    ":::{note}\n",
+    "Click [here](https://userguide.mdanalysis.org/stable/examples/quickstart.html) to learn more about loading files with MDAnalysis,\n",
+    "and [here](https://userguide.mdanalysis.org/stable/selections.html) to learn more about their atom selection language.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In our case the protein is of reasonable size, but if you're working with a very large system, to save some time and memory it may be wise to restrict the protein selection to a sphere around the peptide:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "large_protein_selection = u.select_atoms(\n",
+    "    \"protein and byres around 20.0 group peptide\", peptide=small_protein_selection\n",
+    ")\n",
+    "large_protein_selection"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Depending on your system, it may be of interest to also include water molecules in the protein selection. There are none in this tutorial example but something like this could be used:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "large_protein_selection = u.select_atoms(\n",
+    "    \"(protein or resname WAT) and byres around 20.0 group peptide\",\n",
+    "    peptide=small_protein_selection,\n",
+    ")\n",
+    "large_protein_selection"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, lets make sure that both components were correctly read by MDAnalysis.\n",
+    "\n",
+    ":::{important}\n",
+    "This next step is crucial if you're loading a structure from a file that doesn't\n",
+    "explicitely contain bond orders and formal charges, such as a PDB file or most MD\n",
+    "trajectory files. MDAnalysis will infer those from the atoms connectivity, which\n",
+    "**requires all atoms including hydrogens to be present in the input file**.\n",
+    ":::\n",
+    "\n",
+    "Since ProLIF molecules are built on top of RDKit, we can use RDKit functions to display molecules. Let's have a quick look at our protein selections. We'll only show the first 20 residue to keep the notebook short, this should be enough to detect any compatibility problem between ProLIF/MDAnalysis and your input if any."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create a molecule from the MDAnalysis selection\n",
+    "small_protein_mol = plf.Molecule.from_mda(small_protein_selection)\n",
+    "\n",
+    "# display (remove `slice(20)` to show all residues)\n",
+    "plf.display_residues(small_protein_mol, slice(20))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can do the same for the residues in the other protein selection:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create a molecule from the MDAnalysis selection\n",
+    "large_protein_mol = plf.Molecule.from_mda(large_protein_selection)\n",
+    "\n",
+    "# display\n",
+    "plf.display_residues(large_protein_mol, slice(20))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{admonition} Troubleshooting\n",
+    "\n",
+    "If one of the two molecules was not processed correctly, it might be because your input file does not contain bond information. If that's the case, please add this snippet right after creating your selections:\n",
+    "\n",
+    "```python\n",
+    "small_protein_selection.guess_bonds()\n",
+    "large_protein_selection.guess_bonds()\n",
+    "```\n",
+    "\n",
+    "In some cases, some atomic clashes may be incorrectly classified as bonds and will prevent the conversion of MDAnalysis molecules to RDKit. Since MDAnalysis uses van der Waals radii for bond detection, one can modify the default radii that are used:\n",
+    "\n",
+    "```python\n",
+    "small_protein_selection.guess_bonds(vdwradii={\"H\": 1.05, \"O\": 1.48})\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Fingerprint generation\n",
+    "\n",
+    "Everything looks good, we can now generate a fingerprint. By default, ProLIF will calculate the following interactions: Hydrophobic, HBDonor, HBAcceptor, PiStacking, Anionic, Cationic, CationPi, PiCation, VdWContact.\n",
+    "You can list all interactions that are available with the following command:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plf.Fingerprint.list_available()"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "The default fingerprint will only keep track of the first group of atoms that satisfied the constraints per interaction type and residue pair.\n",
+    "\n",
+    "If you want to keep track of all possible interactions to generate a count-fingerprint (e.g. when there are two atoms in the ligand that make an HBond-donor interaction with residue X), use `plf.Fingerprint(count=True)`.\n",
+    "This is also quite useful for visualization purposes as you can then display the atom pair that has the shortest distance which will look more accurate.\n",
+    "This fingerprint type is however a bit slower to compute.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ignore VdWContact and Hydrophobic interactions\n",
+    "fp = plf.Fingerprint(\n",
+    "    [\n",
+    "        \"HBDonor\",\n",
+    "        \"HBAcceptor\",\n",
+    "        \"PiStacking\",\n",
+    "        \"PiCation\",\n",
+    "        \"CationPi\",\n",
+    "        \"Anionic\",\n",
+    "        \"Cationic\",\n",
+    "    ]\n",
+    ")\n",
+    "# run on a slice of the trajectory frames: from begining to end with a step of 10\n",
+    "fp.run(u.trajectory[::10], small_protein_selection, large_protein_selection)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "The `run` method will automatically select residues that are close to the ligand (6.0 Å) when computing the fingerprint. You can modify the 6.0 Å cutoff by specifying `plf.Fingerprint(vicinity_cutoff=7.0)`, but this is only useful if you decide to change the distance parameters for an interaction class (see in the advanced section of the tutorials).\n",
+    "\n",
+    "Alternatively, you can pass a list of residues like so:\n",
+    "\n",
+    "```python\n",
+    "fp.run(<other parameters>, residues=[\"TYR38.A\", \"ASP129.A\"])\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can save the fingerprint object with `fp.to_pickle` and reload it later with `Fingerprint.from_pickle`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp.to_pickle(\"fingerprint.pkl\")\n",
+    "fp = plf.Fingerprint.from_pickle(\"fingerprint.pkl\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Analysis\n",
+    "\n",
+    "Once the execution is done, you can access the results through `fp.ifp` which is a nested dictionary:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "frame_number = 0\n",
+    "residues = (\"ASP129.A\", \"TYR359.B\")\n",
+    "\n",
+    "fp.ifp[frame_number][residues]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{note}\n",
+    "Internally, ProLIF uses `ligand` and `protein` as a naming convention for components, but it uses whatever was passed as the first selection in the `fp.run` call as the `ligand` and the second object as the `protein`.\n",
+    ":::\n",
+    "\n",
+    "While this contains all the details about the different interactions that were detected, it's not the easiest thing to digest.\n",
+    "\n",
+    "The best way to analyse our results is to export the interaction fingerprint to a Pandas DataFrame. You can read more about pandas in their\n",
+    "[user_guide](https://pandas.pydata.org/docs/user_guide/index.html)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = fp.to_dataframe()\n",
+    "# show only the 10 first frames\n",
+    "df.head(10)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can change the type used for the bits if need be:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "df_uint8 = fp.to_dataframe(dtype=np.uint8)\n",
+    "df_uint8.head(5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here are some common pandas snippets to extract useful information from the fingerprint table.\n",
+    "\n",
+    ":::{important}\n",
+    "Make sure to remove the `.head(5)` at the end of the commands to display the results for all the frames.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# hide an interaction type (HBAcceptor)\n",
+    "df.drop(\"HBAcceptor\", level=\"interaction\", axis=1).head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# show only one protein residue (LYS191.A)\n",
+    "df.xs(\"LYS191.A\", level=\"protein\", axis=1).head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# show only an interaction type (PiStacking)\n",
+    "df.xs(\"PiStacking\", level=\"interaction\", axis=1).head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of the trajectory where each interaction is present\n",
+    "(df.mean().sort_values(ascending=False).to_frame(name=\"%\").T * 100)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# same but we regroup all interaction types\n",
+    "(\n",
+    "    df.T.groupby(level=[\"ligand\", \"protein\"])\n",
+    "    .sum()\n",
+    "    .T.astype(bool)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of the trajectory where PiStacking interactions are present, by residue\n",
+    "(\n",
+    "    df.xs(\"PiStacking\", level=\"interaction\", axis=1)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of the trajectory where interactions with LYS191 occur, by interaction type\n",
+    "(\n",
+    "    df.xs(\"LYS191.A\", level=\"protein\", axis=1)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# percentage of the trajectory where each interaction type is present\n",
+    "(\n",
+    "    df.T.groupby(level=\"interaction\")\n",
+    "    .sum()\n",
+    "    .T.astype(bool)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .to_frame(name=\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 10 residue pairs most frequently interacting\n",
+    "(\n",
+    "    df.T.groupby(level=[\"ligand\", \"protein\"])\n",
+    "    .sum()\n",
+    "    .T.astype(bool)\n",
+    "    .mean()\n",
+    "    .sort_values(ascending=False)\n",
+    "    .head(10)\n",
+    "    .to_frame(\"%\")\n",
+    "    .T\n",
+    "    * 100\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can compute a Tanimoto similarity between frames:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Tanimoto similarity between the first frame and the rest\n",
+    "from rdkit import DataStructs\n",
+    "\n",
+    "bitvectors = fp.to_bitvectors()\n",
+    "tanimoto_sims = DataStructs.BulkTanimotoSimilarity(bitvectors[0], bitvectors)\n",
+    "tanimoto_sims"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To compare binding modes in your trajectory, it's possible to compute the entire similarity matrix:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pandas as pd\n",
+    "import seaborn as sns\n",
+    "from matplotlib import pyplot as plt\n",
+    "\n",
+    "# Tanimoto similarity matrix\n",
+    "bitvectors = fp.to_bitvectors()\n",
+    "similarity_matrix = []\n",
+    "for bv in bitvectors:\n",
+    "    similarity_matrix.append(DataStructs.BulkTanimotoSimilarity(bv, bitvectors))\n",
+    "similarity_matrix = pd.DataFrame(similarity_matrix, index=df.index, columns=df.index)\n",
+    "\n",
+    "# display heatmap\n",
+    "fig, ax = plt.subplots(figsize=(3, 3), dpi=200)\n",
+    "colormap = sns.diverging_palette(\n",
+    "    300, 145, s=90, l=80, sep=30, center=\"dark\", as_cmap=True\n",
+    ")\n",
+    "sns.heatmap(\n",
+    "    similarity_matrix,\n",
+    "    ax=ax,\n",
+    "    square=True,\n",
+    "    cmap=colormap,\n",
+    "    vmin=0,\n",
+    "    vmax=1,\n",
+    "    center=0.5,\n",
+    "    xticklabels=5,\n",
+    "    yticklabels=5,\n",
+    ")\n",
+    "ax.invert_yaxis()\n",
+    "plt.yticks(rotation=\"horizontal\")\n",
+    "fig.patch.set_facecolor(\"white\")"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualisation\n",
+    "\n",
+    "There are a few different options builtin when it comes to visualisation.\n",
+    "\n",
+    "You can start by plotting the interactions over time:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib ipympl\n",
+    "\n",
+    "fp.plot_barcode()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "If you uncomment `%matplotlib ipympl` at the top of the above cell, you should be able to see an interactive version of the plot.\n",
+    ":::\n",
+    "\n",
+    "You can also display the interactions in a 2D interactive diagram if you have a small peptide (not the case here so the plot below won't be readable but that's ok):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp.plot_lignetwork(small_protein_mol)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This diagram is interactive, you can:\n",
+    "- zoom and pan,\n",
+    "- move residues around,\n",
+    "- click on the legend to display or hide types of residues or interactions,\n",
+    "- hover an interaction line to display the distance.\n",
+    "\n",
+    ":::{note}\n",
+    "It is not possible to export it as an image, but you can always take a screenshot.\n",
+    ":::\n",
+    "\n",
+    "You can generate 2 types of diagram with this function, controlled by the `kind` argument:\n",
+    "- `frame`: shows a single specific frame (specified with `frame`, corresponds to the frame number in the simulation).\n",
+    "- `aggregate` (default): the interactions from all frames are grouped and displayed. An optional `threshold` parameter controls the minimum frequency required for an interaction to be displayed (default `0.3`, meaning that interactions occuring in less than 30% of frames will be hidden). The width of interactions is linked to the frequency.\n",
+    "\n",
+    ":::{tip}\n",
+    "Make sure to check the Ligand-Protein tutorial for more examples on this plot.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can also visualize this information in 3D.\n",
+    "\n",
+    "Up to now we've only been using the default fingerprint generator, but you can optionally enable the `count` parameter to enumerate all occurences of an interaction (the default fingerprint generator will stop at the first occurence), and then ProLIF will automatically display the occurence with the smallest distance. You can show all of them with `display_all=True`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp_count = plf.Fingerprint(count=True)\n",
+    "fp_count.run(u.trajectory[0:1], small_protein_selection, large_protein_selection)\n",
+    "view = fp_count.plot_3d(\n",
+    "    small_protein_mol, large_protein_mol, frame=0, display_all=False\n",
+    ")\n",
+    "view"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As in the lignetwork plot, you can hover atoms and interactions to display more information.\n",
+    "\n",
+    "Once you're satisfied with the orientation, you can export the view as a PNG image with the following snippet:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from IPython.display import Javascript\n",
+    "\n",
+    "Javascript(\n",
+    "    \"\"\"\n",
+    "    var png = viewer_%s.pngURI()\n",
+    "    var a = document.createElement('a')\n",
+    "    a.href = png\n",
+    "    a.download = \"prolif-3d.png\"\n",
+    "    a.click()\n",
+    "    a.remove()\n",
+    "\"\"\"\n",
+    "    % view.uniqueid\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Another powerfull visualization option for protein-protein interaction is to use networks to abstract away atomic details.\n",
+    "\n",
+    "`networkx` is a great library for working with graphs in Python, but since the drawing options are quickly limited we will also use `pyvis` to create interactive plots. The following code snippet will convert each residue in our selections to a node, each interaction to an edge, and the occurence of each interaction between residues will be used to control the weight and thickness of each edge."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from html import escape\n",
+    "import networkx as nx\n",
+    "import pandas as pd\n",
+    "from IPython.display import HTML\n",
+    "from matplotlib import colormaps, colors\n",
+    "from pyvis.network import Network\n",
+    "\n",
+    "\n",
+    "def make_graph(\n",
+    "    values: pd.Series,\n",
+    "    df: pd.DataFrame,\n",
+    "    node_color=[\"#FFB2AC\", \"#ACD0FF\"],\n",
+    "    node_shape=\"dot\",\n",
+    "    edge_color=\"#a9a9a9\",\n",
+    "    width_multiplier=1,\n",
+    ") -> nx.Graph:\n",
+    "    \"\"\"Convert a pandas DataFrame to a NetworkX object\n",
+    "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    values : pandas.Series\n",
+    "        Series with 'ligand' and 'protein' levels, and a unique value for\n",
+    "        each lig-prot residue pair that will be used to set the width and weigth\n",
+    "        of each edge. For example:\n",
+    "\n",
+    "            ligand  protein\n",
+    "            LIG1.G  ALA216.A    0.66\n",
+    "                    ALA343.B    0.10\n",
+    "\n",
+    "    df : pandas.DataFrame\n",
+    "        DataFrame obtained from the fp.to_dataframe() method\n",
+    "        Used to label each edge with the type of interaction\n",
+    "\n",
+    "    node_color : list\n",
+    "        Colors for the ligand and protein residues, respectively\n",
+    "\n",
+    "    node_shape : str\n",
+    "        One of ellipse, circle, database, box, text or image, circularImage,\n",
+    "        diamond, dot, star, triangle, triangleDown, square, icon.\n",
+    "\n",
+    "    edge_color : str\n",
+    "        Color of the edge between nodes\n",
+    "\n",
+    "    width_multiplier : int or float\n",
+    "        Each edge's width is defined as `width_multiplier * value`\n",
+    "    \"\"\"\n",
+    "    lig_res = values.index.get_level_values(\"ligand\").unique().tolist()\n",
+    "    prot_res = values.index.get_level_values(\"protein\").unique().tolist()\n",
+    "\n",
+    "    G = nx.Graph()\n",
+    "    # add nodes\n",
+    "    # https://pyvis.readthedocs.io/en/latest/documentation.html#pyvis.network.Network.add_node\n",
+    "    for res in lig_res:\n",
+    "        G.add_node(\n",
+    "            res, title=res, shape=node_shape, color=node_color[0], dtype=\"ligand\"\n",
+    "        )\n",
+    "    for res in prot_res:\n",
+    "        G.add_node(\n",
+    "            res, title=res, shape=node_shape, color=node_color[1], dtype=\"protein\"\n",
+    "        )\n",
+    "\n",
+    "    for resids, value in values.items():\n",
+    "        label = \"{} - {}<br>{}\".format(\n",
+    "            *resids,\n",
+    "            \"<br>\".join(\n",
+    "                [\n",
+    "                    f\"{k}: {v}\"\n",
+    "                    for k, v in (\n",
+    "                        df.xs(resids, level=[\"ligand\", \"protein\"], axis=1)\n",
+    "                        .sum()\n",
+    "                        .to_dict()\n",
+    "                        .items()\n",
+    "                    )\n",
+    "                ]\n",
+    "            ),\n",
+    "        )\n",
+    "        # https://pyvis.readthedocs.io/en/latest/documentation.html#pyvis.network.Network.add_edge\n",
+    "        G.add_edge(\n",
+    "            *resids,\n",
+    "            title=label,\n",
+    "            color=edge_color,\n",
+    "            weight=value,\n",
+    "            width=value * width_multiplier,\n",
+    "        )\n",
+    "\n",
+    "    return G"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data = df.T.groupby(level=[\"ligand\", \"protein\"], sort=False).sum().T.astype(bool).mean()\n",
+    "\n",
+    "G = make_graph(data, df, width_multiplier=8)\n",
+    "\n",
+    "# color each node based on its degree\n",
+    "max_nbr = len(max(G.adj.values(), key=lambda x: len(x)))\n",
+    "blues = colormaps.get_cmap(\"Blues\")\n",
+    "reds = colormaps.get_cmap(\"Reds\")\n",
+    "for n, d in G.nodes(data=True):\n",
+    "    n_neighbors = len(G.adj[n])\n",
+    "    # show the smaller domain in red and the larger one in blue\n",
+    "    palette = reds if d[\"dtype\"] == \"ligand\" else blues\n",
+    "    d[\"color\"] = colors.to_hex(palette(n_neighbors / max_nbr))\n",
+    "\n",
+    "# convert to pyvis network\n",
+    "width, height = (700, 700)\n",
+    "net = Network(width=f\"{width}px\", height=f\"{height}px\", notebook=True, heading=\"\")\n",
+    "net.from_nx(G)\n",
+    "\n",
+    "html_doc = net.generate_html(notebook=True)\n",
+    "iframe = (\n",
+    "    f'<iframe width=\"{width+25}px\" height=\"{height+25}px\" frameborder=\"0\" '\n",
+    "    'srcdoc=\"{html_doc}\"></iframe>'\n",
+    ")\n",
+    "HTML(iframe.format(html_doc=escape(html_doc)))"
+   ]
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae"
+  },
+  "kernelspec": {
+   "display_name": "Python 3.8.5 ('prolif')",
+   "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.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/notebooks/pdb.ipynb b/docs/notebooks/pdb.ipynb
new file mode 100644
index 0000000..3eb5766
--- /dev/null
+++ b/docs/notebooks/pdb.ipynb
@@ -0,0 +1,493 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# PDB\n",
+    "\n",
+    "This tutorial showcases how to use ProLIF to generate an interaction fingerprint for interactions between a protein and a ligand in a PDB file. There are several preparation steps which you need to follow before generating an interaction fingerprint with ProLIF.\n",
+    "\n",
+    ":::{important}\n",
+    "For convenience, the different files for this tutorial are included with the ProLIF installation, and we'll access these files through `plf.datafiles.datapath` which is a {class}`pathlib.Path` object pointing to the tutorials data directory. This makes it easier to manipulate paths to files, match filenames using wildcards...etc. in a Pythonic way, but you can also use plain strings, e.g. `\"/home/user/prolif/data/vina/rec.pdb\"` instead of `plf.datafiles.datapath / \"vina\" / \"rec.pdb\"` if you prefer.\n",
+    "\n",
+    "**Remember to replace any reference to `plf.datafiles.datapath` with the actual paths to your inputs outside of this tutorial.**\n",
+    ":::\n",
+    "\n",
+    ":::{tip}\n",
+    "At the top of the page you can find links to either download this notebook or run it in\n",
+    "Google Colab. You can install the dependencies for the tutorials with the command:\n",
+    "\n",
+    "```shell\n",
+    "pip install prolif[tutorials]\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Prerequisite\n",
+    "\n",
+    "The first step is to split your PDB file in two: one file for the protein and one file for the ligand.\n",
+    "There are different ways to do that, from the command line with `grep`, or from a GUI such as PyMOL.\n",
+    "\n",
+    ":::{important}\n",
+    "Make sure that each file produced either has all bonds explicitely stated, or none of them.\n",
+    ":::\n",
+    "\n",
+    "## Protein preparation\n",
+    "\n",
+    "In order to generate a valid molecule (with bond orders and charges) from a PDB file and correctly detect HBond interactions, **ProLIF requires the protein file to contain explicit hydrogens.**\n",
+    "\n",
+    ":::{note}\n",
+    "Having only polar hydrogens explicit may work if you use RDKit for reading the file, although they are some caveats around protonable residues such as histidines which RDKit may not recognize correctly due to naming conventions, we therefore encourage users to have **all hydrogens explicit** in the file and use MDAnalysis instead.\n",
+    ":::\n",
+    "\n",
+    "In order to protonate residues in your protein file, there are several strategies available, via command-line tools that you can install locally, or webservers, here are a few suggestions:\n",
+    "- [PypKa](https://pypka.org/) available through a [webserver](https://pypka.org/run-pypka/) or as a [CLI tool](https://pypka.readthedocs.io/),\n",
+    "- [H++](http://newbiophysics.cs.vt.edu/H++/) webserver\n",
+    "- [PDB2PQR](https://server.poissonboltzmann.org/) available through a [webserver](https://server.poissonboltzmann.org/pdb2pqr) or as a [CLI tool](https://pdb2pqr.readthedocs.io/). Make sure to then convert the output PQR file to a PDB file.\n",
+    "\n",
+    "Once this step is done, and you've verified that the prepared PDB file contains explicit hydrogens (and optionally explicit bonds), we can finally use that as input.\n",
+    "We have 2 options to read the protein file: MDAnalysis (preferred) or RDKit (faster but risky, see below).\n",
+    "\n",
+    "### With MDAnalysis\n",
+    "\n",
+    "This is the recommended way of reading your prepared protein PDB file:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import MDAnalysis as mda\n",
+    "import prolif as plf\n",
+    "\n",
+    "protein_file = str(plf.datafiles.datapath / \"vina\" / \"rec.pdb\")\n",
+    "\n",
+    "u = mda.Universe(protein_file)\n",
+    "protein_mol = plf.Molecule.from_mda(u)\n",
+    "protein_mol.n_residues"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can have a quick look at each residue individually using the following command (we'll only check residues 260 to 262 as there's an histidine in there which could be problematic):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# check that residues 260 to 262 were processed correctly\n",
+    "# remove the `slice(260, 263)` part to show all residues\n",
+    "plf.display_residues(protein_mol, slice(260, 263))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{admonition} Troubleshooting\n",
+    "\n",
+    "**`Explicit valence for atom <...> is greater than permitted`**\n",
+    "\n",
+    "If your file did not contain explicit bonds, MDAnalysis will automatically run its connectivity-perception algorithm. This could result in some atomic clashes that may be incorrectly classified as bonds and will prevent the conversion of the MDAnalysis molecule to RDKit through ProLIF. Since MDAnalysis uses van der Waals radii for bond detection, one can modify the default radii that are used:\n",
+    "\n",
+    "```python\n",
+    "u.atoms.guess_bonds(vdwradii={\"H\": 1.05, \"O\": 1.48})\n",
+    "```\n",
+    "\n",
+    "Alternatively you could use some GUI tools like PyMOL to write the bonds (CONECT record) in your file and then reread it here.\n",
+    "\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### With RDKit\n",
+    "\n",
+    ":::{important}\n",
+    "While RDKit is going to be much faster at parsing the file, it can only recognize standard residue names, so some protonated residues may be incorrectly parsed.\n",
+    "\n",
+    "For example, the protonated histidine residue HSE347 in our protein is not correcltly parsed which removes aromaticity on the ring, meaning that ProLIF will not be able to match this residue for pi-stacking interactions.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from rdkit import Chem\n",
+    "\n",
+    "rdkit_prot = Chem.MolFromPDBFile(protein_file, removeHs=False)\n",
+    "rdkit_protein_mol = plf.Molecule(rdkit_prot)\n",
+    "# histidine HSE347 not recognized by RDKit\n",
+    "plf.display_residues(rdkit_protein_mol, slice(260, 263))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Ligand preparation\n",
+    "\n",
+    "As for the protein structure, we'll need our prepared ligand file to contain explicit hydrogens.\n",
+    "There are different ways to do this.\n",
+    "\n",
+    "One could load the ligand file (or the whole structure) in PyMOL, select the ligand, add explicit hydrogens, and export the file as SDF, then use the following snippet:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# read SDF\n",
+    "ligand_file = str(plf.datafiles.datapath / \"vina\" / \"vina_output.sdf\")\n",
+    "ligand_mol = plf.sdf_supplier(ligand_file)[0]\n",
+    "# display ligand\n",
+    "plf.display_residues(ligand_mol, size=(400, 200))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If you only have a PDB file with hydrogens, you could use this code snippet instead:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# load PDB with explicit hydrogens\n",
+    "ligand_file = str(plf.datafiles.datapath / \"vina\" / \"lig.pdb\")\n",
+    "u = mda.Universe(ligand_file)\n",
+    "ligand_mol = plf.Molecule.from_mda(u)\n",
+    "# display ligand\n",
+    "plf.display_residues(ligand_mol, size=(400, 200))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Alternatively, PDB files from the Protein Data Bank usually provide the ligand as a SMILES string.\n",
+    "\n",
+    "We can use the following snippet to combine the SMILES string and the ligand PDB file to prepare a valid molecule for ProLIF:\n",
+    "\n",
+    "```python\n",
+    "from rdkit import Chem\n",
+    "from rdkit.Chem.AllChem import AssignBondOrdersFromTemplate\n",
+    "\n",
+    "pdb_ligand = Chem.MolFromPDBFile(str(plf.datafiles.datapath / \"vina\" / \"lig.pdb\"))\n",
+    "smiles_mol = Chem.MolFromSmiles(\"C[NH+]1CC(C(=O)NC2(C)OC3(O)C4CCCN4C(=O)C(Cc4ccccc4)N3C2=O)C=C2c3cccc4[nH]cc(c34)CC21\")\n",
+    "mol = AssignBondOrdersFromTemplate(smiles_mol, pdb_ligand)\n",
+    "mol_h = Chem.AddHs(mol)\n",
+    "ligand_mol = plf.Molecule.from_rdkit(mol_h)\n",
+    "```"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Fingerprint generation\n",
+    "\n",
+    "We can now generate a fingerprint. By default, ProLIF will calculate the following interactions: Hydrophobic, HBDonor, HBAcceptor, PiStacking, Anionic, Cationic, CationPi, PiCation, VdWContact.\n",
+    "You can list all interactions that are available with the following command:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plf.Fingerprint.list_available()"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "The default fingerprint will only keep track of the first group of atoms that satisfied the constraints per interaction type and residue pair.\n",
+    "\n",
+    "If you want to keep track of all possible interactions to generate a count-fingerprint (e.g. when there are two atoms in the ligand that make an HBond-donor interaction with residue X), use `plf.Fingerprint(count=True)`.\n",
+    "This is also quite useful for visualization purposes as you can then display the atom pair that has the shortest distance which will look more accurate.\n",
+    "This fingerprint type is however a bit slower to compute.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# use default interactions\n",
+    "fp = plf.Fingerprint()\n",
+    "# run on your poses\n",
+    "fp.run_from_iterable([ligand_mol], protein_mol)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{tip}\n",
+    "The `run_from_iterable` method will automatically select residues that are close to the ligand (6.0 Å) when computing the fingerprint. You can modify the 6.0 Å cutoff by specifying `plf.Fingerprint(vicinity_cutoff=7.0)`, but this is only useful if you decide to change the distance parameters for an interaction class (see in the advanced section of the tutorials).\n",
+    "\n",
+    "Alternatively, you can pass a list of residues like so:\n",
+    "\n",
+    "```python\n",
+    "fp.run_from_iterable(<other parameters>, residues=[\"TYR38.A\", \"ASP129.A\"])\n",
+    "```\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can save the fingerprint object with `fp.to_pickle` and reload it later with `Fingerprint.from_pickle`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp.to_pickle(\"fingerprint.pkl\")\n",
+    "fp = plf.Fingerprint.from_pickle(\"fingerprint.pkl\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Analysis\n",
+    "\n",
+    "Once the execution is done, you can access the results through `fp.ifp` which is a nested dictionary:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ligand_residue = \"LIG1.G\"\n",
+    "protein_residue = \"ASP129.A\"\n",
+    "\n",
+    "fp.ifp[0][(ligand_residue, protein_residue)]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "While this contains all the details about the different interactions that were detected, it's not the easiest thing to digest.\n",
+    "\n",
+    "The best way to analyse our results is to export the interaction fingerprint to a Pandas DataFrame. You can read more about pandas in their\n",
+    "[user_guide](https://pandas.pydata.org/docs/user_guide/index.html)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = fp.to_dataframe()\n",
+    "df.T"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here are some common pandas snippets to extract useful information from the fingerprint table."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# hide an interaction type (Hydrophobic)\n",
+    "df.drop(\"Hydrophobic\", level=\"interaction\", axis=1).T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# show only one protein residue (ASP129.A)\n",
+    "df.xs(\"ASP129.A\", level=\"protein\", axis=1).T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# show only an interaction type (PiStacking)\n",
+    "df.xs(\"PiStacking\", level=\"interaction\", axis=1).T"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualisation\n",
+    "\n",
+    "There are a few different options builtin when it comes to visualisation.\n",
+    "\n",
+    "You can display the interactions in a 2D interactive diagram:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp.plot_lignetwork(ligand_mol, kind=\"frame\", frame=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This diagram is interactive, you can:\n",
+    "- zoom and pan,\n",
+    "- move residues around,\n",
+    "- click on the legend to display or hide types of residues or interactions,\n",
+    "- hover an interaction line to display the distance.\n",
+    "\n",
+    ":::{note}\n",
+    "It is not possible to export it as an image, but you can always take a screenshot.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Up to now we've only been using the default fingerprint generator, but you can optionally enable the `count` parameter to enumerate all occurences of an interaction (the default fingerprint generator will stop at the first occurence), and then display all of them by specifying `display_all=True`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fp_count = plf.Fingerprint(count=True)\n",
+    "fp_count.run_from_iterable([ligand_mol], protein_mol)\n",
+    "fp_count.plot_lignetwork(ligand_mol, kind=\"frame\", frame=0, display_all=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can also visualize this information in 3D:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "view = fp_count.plot_3d(ligand_mol, protein_mol, frame=0, display_all=False)\n",
+    "view"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As in the lignetwork plot, you can hover atoms and interactions to display more information.\n",
+    "\n",
+    "The advantage of using a count fingerprint in that case is that it will automatically select the interaction occurence with the shortest distance for a more intuitive visualization.\n",
+    "\n",
+    "Once you're satisfied with the orientation, you can export the view as a PNG image with the following snippet:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from IPython.display import Javascript\n",
+    "\n",
+    "Javascript(\n",
+    "    \"\"\"\n",
+    "    var png = viewer_%s.pngURI()\n",
+    "    var a = document.createElement('a')\n",
+    "    a.href = png\n",
+    "    a.download = \"prolif-3d.png\"\n",
+    "    a.click()\n",
+    "    a.remove()\n",
+    "\"\"\"\n",
+    "    % view.uniqueid\n",
+    ")"
+   ]
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae"
+  },
+  "kernelspec": {
+   "display_name": "Python 3.8.5 ('prolif')",
+   "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.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/notebooks/protein-protein_interactions.ipynb b/docs/notebooks/protein-protein_interactions.ipynb
deleted file mode 100644
index 1039185..0000000
--- a/docs/notebooks/protein-protein_interactions.ipynb
+++ /dev/null
@@ -1,236 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Protein-protein interactions\n",
-    "\n",
-    "This notebooks shows how to compute a fingerprint for protein-protein interactions.\n",
-    "\n",
-    "Here we will investigate the interactions in a G-protein coupled receptor (GPCR) between a particular helix (called TM3) and the rest of the protein.\n",
-    "\n",
-    "This can obviously be applied to proteins that don't belong to the same chain/segment, as long as you can figure out an appropriate [MDAnalysis selection](https://docs.mdanalysis.org/stable/documentation_pages/selections.html)\n",
-    "\n",
-    "There is also an example at the end of this tutorial for generating an IFP of PPI without considering the backbone."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2020-03-12T18:24:55.324128Z",
-     "start_time": "2020-03-12T18:24:54.314383Z"
-    },
-    "tags": []
-   },
-   "outputs": [],
-   "source": [
-    "import MDAnalysis as mda\n",
-    "import prolif as plf"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# load traj\n",
-    "u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)\n",
-    "tm3 = u.select_atoms(\"resid 119:152\")\n",
-    "prot = u.select_atoms(\"protein and not group tm3\", tm3=tm3)\n",
-    "tm3, prot"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": true,
-    "tags": []
-   },
-   "outputs": [],
-   "source": [
-    "# prot-prot interactions\n",
-    "fp = plf.Fingerprint(\n",
-    "    [\n",
-    "        \"HBDonor\",\n",
-    "        \"HBAcceptor\",\n",
-    "        \"PiStacking\",\n",
-    "        \"PiCation\",\n",
-    "        \"CationPi\",\n",
-    "        \"Anionic\",\n",
-    "        \"Cationic\",\n",
-    "    ]\n",
-    ")\n",
-    "fp.run(u.trajectory[::10], tm3, prot)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "df = fp.to_dataframe()\n",
-    "df.head()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# show interactions for a specific ligand residue\n",
-    "df.xs(\"ARG147.A\", level=\"ligand\", axis=1).head(5)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# same for a protein residue\n",
-    "df.xs(\"GLU309.B\", level=\"protein\", axis=1).head(5)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# display a specific type of interaction\n",
-    "df.xs(\"Cationic\", level=\"interaction\", axis=1).head(5)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# calculate the occurence of each interaction on the trajectory\n",
-    "occ = df.mean()\n",
-    "# restrict to the frequent ones\n",
-    "occ.loc[occ > 0.3]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# regroup all interactions together and do the same\n",
-    "g = df.groupby(level=[\"ligand\", \"protein\"], axis=1).sum().astype(bool).mean()\n",
-    "g.loc[g > 0.3]"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Note: Add `%matplotlib ipympl` at the top of the cell for an interactive visualization:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fp.plot_barcode()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Ignoring backbone interactions\n",
-    "\n",
-    "In some cases, you might want to dismiss backbone interactions. While it might be tempting to just modify the MDAnalysis selection with `\"protein and not backbone\"`, this won't work as expected and will lead to adding a charges where the backbone was bonding with the sidechain.  \n",
-    "However there is a temporary workaround (which will be directly included in the code in the near future):"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from rdkit import Chem\n",
-    "from rdkit.Chem import AllChem\n",
-    "from tqdm.auto import tqdm\n",
-    "\n",
-    "# remove backbone\n",
-    "backbone = Chem.MolFromSmarts(\"[C^2](=O)-[C;X4](-[H])-[N;+0]\")\n",
-    "fix_h = Chem.MolFromSmarts(\"[H&D0]\")\n",
-    "\n",
-    "\n",
-    "def remove_backbone(atomgroup):\n",
-    "    mol = plf.Molecule.from_mda(atomgroup)\n",
-    "    mol = AllChem.DeleteSubstructs(mol, backbone)\n",
-    "    mol = AllChem.DeleteSubstructs(mol, fix_h)\n",
-    "    return plf.Molecule(mol)\n",
-    "\n",
-    "\n",
-    "# generate IFP\n",
-    "ifp = {}\n",
-    "for ts in tqdm(u.trajectory[::10]):\n",
-    "    tm3_mol = remove_backbone(tm3)\n",
-    "    prot_mol = remove_backbone(prot)\n",
-    "    ifp[ts.frame] = fp.generate(tm3_mol, prot_mol, metadata=True)\n",
-    "df = plf.to_dataframe(ifp, fp.interactions.keys())\n",
-    "df.head()"
-   ]
-  }
- ],
- "metadata": {
-  "hide_input": false,
-  "interpreter": {
-   "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae"
-  },
-  "kernelspec": {
-   "display_name": "Python 3.8.5 ('prolif')",
-   "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.9.13"
-  },
-  "toc": {
-   "base_numbering": 1,
-   "nav_menu": {},
-   "number_sections": true,
-   "sideBar": true,
-   "skip_h1_title": false,
-   "title_cell": "Table of Contents",
-   "title_sidebar": "Contents",
-   "toc_cell": false,
-   "toc_position": {
-    "height": "calc(100% - 180px)",
-    "left": "10px",
-    "top": "150px",
-    "width": "275px"
-   },
-   "toc_section_display": true,
-   "toc_window_display": true
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/notebooks/quickstart.ipynb b/docs/notebooks/quickstart.ipynb
deleted file mode 100644
index 15cd2dc..0000000
--- a/docs/notebooks/quickstart.ipynb
+++ /dev/null
@@ -1,327 +0,0 @@
-{
- "cells": [
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Quickstart\n",
-    "\n",
-    "This is a very short guide on how to use ProLIF to generate an interaction fingerprint for a ligand-protein complex from an MD simulation.\n",
-    "\n",
-    "Let's start by importing MDAnalysis and ProLIF to read our input files:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import MDAnalysis as mda\n",
-    "import prolif as plf\n",
-    "\n",
-    "# load trajectory\n",
-    "u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)\n",
-    "# create selections for the ligand and protein\n",
-    "lig = u.atoms.select_atoms(\"resname LIG\")\n",
-    "prot = u.atoms.select_atoms(\"protein and byres around 20 group ligand\", ligand=lig)\n",
-    "lig, prot"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "MDAnalysis should automatically recognize the file type that you're using from its extension. Click [here](https://userguide.mdanalysis.org/stable/examples/quickstart.html) to learn more about loading files with MDAnalysis, and [here](https://userguide.mdanalysis.org/stable/selections.html) to learn more about their atom selection language.\n",
-    "\n",
-    "Next, lets make sure that our ligand was correctly read by MDAnalysis.\n",
-    "\n",
-    "⚠ This next step is crucial if you're loading a structure from a file that doesn't explicitely contain bond orders and formal charges, such as a PDB file or most MD trajectory files. MDAnalysis will infer those from the atoms connectivity, which **requires all atoms including hydrogens to be present in the input file**.\n",
-    "\n",
-    "ProLIF molecules are built on top of RDKit and are compatible with its drawing code. Let's have a quick look at our ligand:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from rdkit import Chem\n",
-    "from rdkit.Chem import Draw\n",
-    "\n",
-    "# create a molecule from the MDAnalysis selection\n",
-    "lmol = plf.Molecule.from_mda(lig)\n",
-    "# cleanup before drawing\n",
-    "mol = Chem.RemoveHs(lmol)\n",
-    "mol.RemoveAllConformers()\n",
-    "Draw.MolToImage(mol, size=(400, 200))"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We can do the same for the residues in the protein (I'll only show the first 20 to keep the notebook short):"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "pmol = plf.Molecule.from_mda(prot)\n",
-    "frags = []\n",
-    "# to show all residues, simply use `for res in pmol:`\n",
-    "for i in range(20):\n",
-    "    res = pmol[i]\n",
-    "    mol = Chem.RemoveHs(res)\n",
-    "    mol.RemoveAllConformers()\n",
-    "    frags.append(mol)\n",
-    "Draw.MolsToGridImage(\n",
-    "    frags,\n",
-    "    legends=[str(res.resid) for res in pmol],\n",
-    "    subImgSize=(200, 140),\n",
-    "    molsPerRow=4,\n",
-    "    maxMols=prot.n_residues,\n",
-    ")"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "If only one of the two molecules was not read correctly, it might be because your input file doesn't contain bond information (typically the case for PDB files) so MDAnalysis has to guess it. This is usually done automatically, but when it's done on two different atom selections it will fail, so one has to explicitely tell MDAnalysis to do bond detection for the whole system.  \n",
-    "To fix this, simply add `guess_bonds=True` when creating the Universe: `u = mda.Universe(\"file.pdb\", guess_bonds=True)`.\n",
-    "\n",
-    "In some cases, some atomic clashes may be incorrectly classified as bonds and will prevent the conversion of MDAnalysis molecules to RDKit, but workarounds exist:\n",
-    "\n",
-    "* Since MDAnalysis uses van der Waals radii for bond detection, one can modify the default radii that are used:\n",
-    "  `u = mda.Universe(\"file.pdb\", guess_bonds=True, vdwradii={\"H\": 1.05, \"O\": 1.48})`\n",
-    "* If you have separate PDB files for your ligand and protein, you can also directly use RDKit:\n",
-    "\n",
-    "```python\n",
-    "mol = Chem.MolFromPDBFile(\"prot.pdb\", removeHs=False)\n",
-    "prot = plf.Molecule(mol)\n",
-    "mol = Chem.MolFromPDBFile(\"lig.pdb\", removeHs=False)\n",
-    "lig = plf.Molecule(mol)\n",
-    "fp = plf.Fingerprint()\n",
-    "fp.run_from_iterable([lig], prot)\n",
-    "```"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Everything looks good, we can now generate a fingerprint. By default, ProLIF will calculate the following interactions: Hydrophobic, HBDonor, HBAcceptor, PiStacking, Anionic, Cationic, CationPi, PiCation, VdWContact.\n",
-    "You can list all interactions that are available with the following command:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plf.Fingerprint.list_available()"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The default fingerprint will only keep track of the first group of atoms that satisfied the constraints per interaction type and residue pair. If you want to keep track of all possible interactions to generate a count-fingerprint (e.g. \"there are two atoms in the ligand that make an HBond-donor interaction with residue X\"), use `plf.Fingerprint(count=True)`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# use defaults\n",
-    "fp = plf.Fingerprint()\n",
-    "# run on a slice of frames from begining to end with a step of 10\n",
-    "fp.run(u.trajectory[::10], lig, prot)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The `run` method will automatically select residues that are close to the ligand (6.0 Å) when computing the fingerprint. Alternatively, you can pass a list of residues like so:\n",
-    "\n",
-    "```python\n",
-    "fp.run(..., residues=[\"TYR38.A\", \"ASP129.A\"])\n",
-    "```\n",
-    "Or simply use `fp.run(..., residues=\"all\")` to use all residues in the `prot` selection.\n",
-    "\n",
-    "The best way to access our results is to export our interaction fingerprints to a Pandas DataFrame. By default, the resulting DataFrame only keeps track of residues and interaction types that were seen in at least one of the frames in your trajectory. If needed, you can access the full results with `fp.to_dataframe(drop_empty=False)`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "df = fp.to_dataframe()\n",
-    "# show only the 10 first frames\n",
-    "df.head(10)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "If the generated DataFrame is empty and you are using a PDB file, please manually call `u.atoms.guess_bonds()` right after creating the Universe with MDAnalysis, and rerun the code up to this point. This problem typically happens when the PDB file is missing all or some of the `CONECT` records to indicate bonds, which leads to structures with atoms but no bonds. By design, those structures cannot match any of the patterns that are used to define interactions, hence no interaction is detected, resulting in an empty dataframe."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# drop the ligand residue column since there's only a single ligand residue\n",
-    "df = df.droplevel(\"ligand\", axis=1)\n",
-    "df.head(5)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# show all pi-stacking interactions\n",
-    "df.xs(\"PiStacking\", level=\"interaction\", axis=1).head(5)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# show all interactions with a specific protein residue\n",
-    "df.xs(\"ASP129.A\", level=\"protein\", axis=1).head(5)\n",
-    "# or more simply\n",
-    "df[\"ASP129.A\"].head(5)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Here's a simple example to plot the interactions over time. Add `%matplotlib ipympl` at the top of the cell for an interactive visualization:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fp.plot_barcode()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# calculate the occurence of each interaction on the trajectory\n",
-    "occ = df.mean()\n",
-    "# restrict to the frequent ones\n",
-    "occ.loc[occ > 0.3]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# regroup all interactions together and do the same\n",
-    "g = df.groupby(level=[\"protein\"], axis=1).sum().astype(bool).mean()\n",
-    "g.loc[g > 0.3]"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "You can also compute a Tanimoto similarity between each frame:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from rdkit import DataStructs\n",
-    "\n",
-    "bvs = fp.to_bitvectors()\n",
-    "tanimoto_sims = DataStructs.BulkTanimotoSimilarity(bvs[0], bvs)\n",
-    "tanimoto_sims"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "You can save the fingerprint object with `fp.to_pickle` and reload it later with `Fingerprint.from_pickle`:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fp.to_pickle(\"fingerprint.pkl\")\n",
-    "fp = plf.Fingerprint.from_pickle(\"fingerprint.pkl\")"
-   ]
-  }
- ],
- "metadata": {
-  "interpreter": {
-   "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae"
-  },
-  "kernelspec": {
-   "display_name": "Python 3.8.5 ('prolif')",
-   "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.9.13"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/notebooks/visualisation.ipynb b/docs/notebooks/visualisation.ipynb
deleted file mode 100644
index 242cce9..0000000
--- a/docs/notebooks/visualisation.ipynb
+++ /dev/null
@@ -1,524 +0,0 @@
-{
- "cells": [
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Visualisation\n",
-    "\n",
-    "This notebook showcases different ways of visualizing ligand-protein and protein-protein interactions, either with atomistic details or simply at the residue level.\n",
-    "\n",
-    "By default, ProLIF will only retain the first group of atoms that satisfy the constraints (SMARTS patterns and distance, angles...etc) for a given interaction and pair of residues. This means that the some of the visualization features may highlight atoms that are not the \"optimal\" one.\n",
-    "To mitigate this, we can use a count-fingerprint (`prolif.Fingerprint(count=True)`) to keep track of all occurences of an interaction between a pair of residues, and then restrict ourselves to only display the one with the shortest distance. The downside is that this count-fingerprint can be slightly slower to generate."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import MDAnalysis as mda\n",
-    "import prolif as plf\n",
-    "\n",
-    "# load topology\n",
-    "u = mda.Universe(plf.datafiles.TOP, plf.datafiles.TRAJ)\n",
-    "lig = u.select_atoms(\"resname LIG\")\n",
-    "# restrict protein to 10 angstroms around the ligand\n",
-    "prot = u.select_atoms(\"protein and byres around 10 group ligand\", ligand=lig)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# create RDKit-like molecules for visualisation\n",
-    "lmol = plf.Molecule.from_mda(lig)\n",
-    "pmol = plf.Molecule.from_mda(prot)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## py3Dmol (3Dmol.js)\n",
-    "\n",
-    "With py3Dmol we can easily display the interactions.\n",
-    "\n",
-    "For interactions involving a ring (pi-cation, pi-stacking...etc.) ProLIF returns the indices of the ring atoms, but for visualisation having the centroid of the ring looks nicer. We'll start by writing a function to find the centroid, given the indices of the ring atoms."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Now let's run the fingerprint analysis (we'll just do 1 frame)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fp = plf.Fingerprint(count=True)\n",
-    "fp.run(u.trajectory[0:1], lig, prot)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Finally, the actual visualisation code. The API of py3Dmol is exactly the same as the GLViewer class of 3Dmol.js, for which the documentation can be found [here](https://3dmol.csb.pitt.edu/doc/$3Dmol.GLViewer.html)."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "fp.plot_3d(lmol, pmol, frame=0, display_all=True)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Ligand Interaction Network (LigPlot)\n",
-    "\n",
-    "Protein-ligand interactions are typically represented with the ligand in atomic details, residues as nodes, and interactions as edges. Such diagram can be easily displayed by calling ProLIF's builtin class `prolif.plotting.network.LigNetwork`.\n",
-    "This diagram is interactive and allows moving around the residues, as well as clicking the legend to toggle the display of specific residues types or interactions.\n",
-    "LigNetwork can generate two kinds of depictions:\n",
-    "\n",
-    "- Based on a single specific frame\n",
-    "- By aggregating results from several frames\n",
-    "\n",
-    "In the latter case, the frequency with which an interaction is seen will control the width of the corresponding edge. You can hide the least frequent interactions by using a threshold, *i.e.* `threshold=0.3` will hide interactions that occur in less than 30% of frames."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# run the analysis for a couple of frames\n",
-    "fp = plf.Fingerprint(count=True)\n",
-    "fp.run(u.trajectory[::10], lig, prot)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# display\n",
-    "fp.plot_lignetwork(\n",
-    "    lmol,\n",
-    "    # replace with `kind=\"frame\", frame=0,` for the other depiction\n",
-    "    kind=\"aggregate\",\n",
-    "    threshold=0.3,\n",
-    ")"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "You can further customize the diagram by changing the colors in `LigNetwork.COLORS` or the residues types in `LigNetwork.RESIDUE_TYPES`. Type `help(LigNetwork)` for more details.  \n",
-    "The diagram can be saved as an HTML file by calling `net.save(\"output.html\")`. It is not currently possible to export it as an image, so please make a screenshot instead.  \n",
-    "You can combine both saving and displaying the diagram with `net.show(\"output.html\")`."
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## NetworkX and pyvis\n",
-    "\n",
-    "NetworkX is a great library for working with graphs, but the drawing options are quickly limited so we will use networkx to create a graph, and pyvis to create interactive plots. The following code snippet will calculate the IFP, each residue (ligand or protein) is converted to a node, each interaction to an edge, and the occurence of each interaction between residues will be used to control the weight and thickness of each edge."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import networkx as nx\n",
-    "from pyvis.network import Network\n",
-    "from matplotlib import cm, colors\n",
-    "from IPython.display import IFrame"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "# get lig-prot interactions and distance between residues\n",
-    "\n",
-    "fp = plf.Fingerprint()\n",
-    "fp.run(u.trajectory[::10], lig, prot)\n",
-    "df = fp.to_dataframe()\n",
-    "df.head()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def make_graph(\n",
-    "    values,\n",
-    "    df=None,\n",
-    "    node_color=[\"#FFB2AC\", \"#ACD0FF\"],\n",
-    "    node_shape=\"dot\",\n",
-    "    edge_color=\"#a9a9a9\",\n",
-    "    width_multiplier=1,\n",
-    "):\n",
-    "    \"\"\"Convert a pandas DataFrame to a NetworkX object\n",
-    "\n",
-    "    Parameters\n",
-    "    ----------\n",
-    "    values : pandas.Series\n",
-    "        Series with 'ligand' and 'protein' levels, and a unique value for\n",
-    "        each lig-prot residue pair that will be used to set the width and weigth\n",
-    "        of each edge. For example:\n",
-    "\n",
-    "            ligand  protein\n",
-    "            LIG1.G  ALA216.A    0.66\n",
-    "                    ALA343.B    0.10\n",
-    "\n",
-    "    df : pandas.DataFrame\n",
-    "        DataFrame obtained from the fp.to_dataframe() method\n",
-    "        Used to label each edge with the type of interaction\n",
-    "\n",
-    "    node_color : list\n",
-    "        Colors for the ligand and protein residues, respectively\n",
-    "\n",
-    "    node_shape : str\n",
-    "        One of ellipse, circle, database, box, text or image, circularImage,\n",
-    "        diamond, dot, star, triangle, triangleDown, square, icon.\n",
-    "\n",
-    "    edge_color : str\n",
-    "        Color of the edge between nodes\n",
-    "\n",
-    "    width_multiplier : int or float\n",
-    "        Each edge's width is defined as `width_multiplier * value`\n",
-    "    \"\"\"\n",
-    "    lig_res = values.index.get_level_values(\"ligand\").unique().tolist()\n",
-    "    prot_res = values.index.get_level_values(\"protein\").unique().tolist()\n",
-    "\n",
-    "    G = nx.Graph()\n",
-    "    # add nodes\n",
-    "    # https://pyvis.readthedocs.io/en/latest/documentation.html#pyvis.network.Network.add_node\n",
-    "    for res in lig_res:\n",
-    "        G.add_node(\n",
-    "            res, title=res, shape=node_shape, color=node_color[0], dtype=\"ligand\"\n",
-    "        )\n",
-    "    for res in prot_res:\n",
-    "        G.add_node(\n",
-    "            res, title=res, shape=node_shape, color=node_color[1], dtype=\"protein\"\n",
-    "        )\n",
-    "\n",
-    "    for resids, value in values.items():\n",
-    "        label = \"{} - {}<br>{}\".format(\n",
-    "            *resids,\n",
-    "            \"<br>\".join(\n",
-    "                [\n",
-    "                    f\"{k}: {v}\"\n",
-    "                    for k, v in (\n",
-    "                        df.xs(resids, level=[\"ligand\", \"protein\"], axis=1)\n",
-    "                        .sum()\n",
-    "                        .to_dict()\n",
-    "                        .items()\n",
-    "                    )\n",
-    "                ]\n",
-    "            ),\n",
-    "        )\n",
-    "        # https://pyvis.readthedocs.io/en/latest/documentation.html#pyvis.network.Network.add_edge\n",
-    "        G.add_edge(\n",
-    "            *resids,\n",
-    "            title=label,\n",
-    "            color=edge_color,\n",
-    "            weight=value,\n",
-    "            width=value * width_multiplier,\n",
-    "        )\n",
-    "\n",
-    "    return G"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Regrouping all interactions\n",
-    "\n",
-    "We will regroup all interactions as if they were equivalent."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "data = df.groupby(level=[\"ligand\", \"protein\"], axis=1).sum().astype(bool).mean()\n",
-    "\n",
-    "G = make_graph(data, df, width_multiplier=3)\n",
-    "\n",
-    "# display graph\n",
-    "net = Network(width=600, height=500, notebook=True, heading=\"\")\n",
-    "net.from_nx(G)\n",
-    "net.write_html(\"lig-prot_graph.html\")\n",
-    "IFrame(\"lig-prot_graph.html\", width=610, height=510)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Only plotting a specific interaction\n",
-    "\n",
-    "We can also plot a specific type of interaction."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "data = df.xs(\"Hydrophobic\", level=\"interaction\", axis=1).mean()\n",
-    "\n",
-    "G = make_graph(data, df, width_multiplier=3)\n",
-    "\n",
-    "# display graph\n",
-    "net = Network(width=600, height=500, notebook=True, heading=\"\")\n",
-    "net.from_nx(G)\n",
-    "net.write_html(\"lig-prot_hydrophobic_graph.html\")\n",
-    "IFrame(\"lig-prot_hydrophobic_graph.html\", width=610, height=510)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Protein-protein interaction\n",
-    "\n",
-    "This kind of \"residue-level\" visualisation is especially suitable for protein-protein interactions. Here we'll show the interactions between one helix of our G-Protein coupled receptor (transmembrane helix 3, or TM3) in red and the rest of the protein in blue."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "tm3 = u.select_atoms(\"resid 119:152\")\n",
-    "prot = u.select_atoms(\"protein and not group tm3\", tm3=tm3)\n",
-    "fp = plf.Fingerprint()\n",
-    "fp.run(u.trajectory[::25], tm3, prot)\n",
-    "df = fp.to_dataframe()\n",
-    "df.head()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "data = (\n",
-    "    df.groupby(level=[\"ligand\", \"protein\"], axis=1, sort=False)\n",
-    "    .sum()\n",
-    "    .astype(bool)\n",
-    "    .mean()\n",
-    ")\n",
-    "\n",
-    "G = make_graph(data, df, width_multiplier=8)\n",
-    "\n",
-    "# color each node based on its degree\n",
-    "max_nbr = len(max(G.adj.values(), key=lambda x: len(x)))\n",
-    "blues = cm.get_cmap(\"Blues\", max_nbr)\n",
-    "reds = cm.get_cmap(\"Reds\", max_nbr)\n",
-    "for n, d in G.nodes(data=True):\n",
-    "    n_neighbors = len(G.adj[n])\n",
-    "    # show TM3 in red and the rest of the protein in blue\n",
-    "    palette = reds if d[\"dtype\"] == \"ligand\" else blues\n",
-    "    d[\"color\"] = colors.to_hex(palette(n_neighbors / max_nbr))\n",
-    "\n",
-    "# convert to pyvis network\n",
-    "net = Network(width=640, height=500, notebook=True, heading=\"\")\n",
-    "net.from_nx(G)\n",
-    "net.write_html(\"prot-prot_graph.html\")\n",
-    "IFrame(\"prot-prot_graph.html\", width=650, height=510)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Residue interaction network\n",
-    "\n",
-    "Another possible application is the visualisation of the residue interaction network of the whole protein. Since this protein is a GPCR, the graph will mostly display the HBond interactions reponsible for the secondary structure of the protein (7 alpha-helices). It would also show hydrophobic interactions between neighbor residues, so I'm simply going to disable it in the Fingerprint."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# prot = u.select_atoms(\"protein\")\n",
-    "# fp = plf.Fingerprint(\n",
-    "#     [\n",
-    "#         \"HBDonor\",\n",
-    "#         \"HBAcceptor\",\n",
-    "#         \"PiStacking\",\n",
-    "#         \"Anionic\",\n",
-    "#         \"Cationic\",\n",
-    "#         \"CationPi\",\n",
-    "#         \"PiCation\",\n",
-    "#     ]\n",
-    "# )\n",
-    "# fp.run(u.trajectory[::25], prot, prot, n_jobs=1)\n",
-    "# df = fp.to_dataframe()\n",
-    "# df.head()"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "To hide most of the HBond interactions responsible for the alpha-helix structuration, I will show how to do it on the pandas DataFrame for simplicity, but ideally you should copy-paste the source code inside the `fp.run` method and add the condition shown below before calculating the bitvector for a residue pair, then use the custom function instead of `fp.run`. This would make the analysis faster and more memory efficient."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# # remove interactions between residues i and i±4 or less\n",
-    "# mask = []\n",
-    "# for l, p, interaction in df.columns:\n",
-    "#     lr = plf.ResidueId.from_string(l)\n",
-    "#     pr = plf.ResidueId.from_string(p)\n",
-    "#     if (pr == lr) or (\n",
-    "#         abs(pr.number - lr.number) <= 4\n",
-    "#         and interaction in [\"HBDonor\", \"HBAcceptor\", \"Hydrophobic\"]\n",
-    "#     ):\n",
-    "#         mask.append(False)\n",
-    "#     else:\n",
-    "#         mask.append(True)\n",
-    "# df = df[df.columns[mask]]\n",
-    "# df.head()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# data = (\n",
-    "#     df.groupby(level=[\"ligand\", \"protein\"], axis=1, sort=False)\n",
-    "#     .sum()\n",
-    "#     .astype(bool)\n",
-    "#     .mean()\n",
-    "# )\n",
-    "\n",
-    "# G = make_graph(data, df, width_multiplier=5)\n",
-    "\n",
-    "# # color each node based on its degree\n",
-    "# max_nbr = len(max(G.adj.values(), key=lambda x: len(x)))\n",
-    "# palette = cm.get_cmap(\"YlGnBu\", max_nbr)\n",
-    "# for n, d in G.nodes(data=True):\n",
-    "#     n_neighbors = len(G.adj[n])\n",
-    "#     d[\"color\"] = colors.to_hex(palette(n_neighbors / max_nbr))\n",
-    "\n",
-    "# # convert to pyvis network\n",
-    "# net = Network(width=640, height=500, notebook=True, heading=\"\")\n",
-    "# net.from_nx(G)\n",
-    "\n",
-    "# # use specific layout\n",
-    "# layout = nx.circular_layout(G)\n",
-    "# for node in net.nodes:\n",
-    "#     node[\"x\"] = layout[node[\"id\"]][0] * 1000\n",
-    "#     node[\"y\"] = layout[node[\"id\"]][1] * 1000\n",
-    "# net.toggle_physics(False)\n",
-    "\n",
-    "# net.write_html(\"residue-network_graph.html\")\n",
-    "# IFrame(\"residue-network_graph.html\", width=650, height=510)"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "End of this notebook. If you have other suggestions for displaying interaction fingerprints, please create a new [Discussion](https://github.com/chemosim-lab/ProLIF/discussions) on GitHub 👍"
-   ]
-  },
-  {
-   "attachments": {},
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "*List of files to be automatically copied to the docs:*\n",
-    "\n",
-    "* [lig-prot_graph.html](lig-prot_graph.html)\n",
-    "* [lig-prot_hydrophobic_graph.html](lig-prot_hydrophobic_graph.html)\n",
-    "* [prot-prot_graph.html](prot-prot_graph.html)\n",
-    "* [residue-network_graph.html](residue-network_graph.html)"
-   ]
-  }
- ],
- "metadata": {
-  "interpreter": {
-   "hash": "8787e9fc73b27535744a25d17e74686c0add9df598b8e27ca04412fce7f0c7ae"
-  },
-  "kernelspec": {
-   "display_name": "Python 3.8.5 ('prolif')",
-   "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.9.13"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/docs/source/api.rst b/docs/source/api.rst
index 3c8e32c..f49b68a 100644
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -20,8 +20,8 @@ Interaction fingerprint
 -----------------------
 
 An interaction fingerprint decomposes the interactions between two molecules
-in a binary vector. The interactions are detected by looking up for predefined
-molecular patterns (like SMARTS queries) that satisfy geometrical
+in a binary vector. The interactions are detected by looking up predefined
+molecular patterns (using SMARTS queries in ProLIF) that satisfy geometrical
 constraints (distance, angle, dihedral...).
 To learn more on how these interactions are defined and how to use prolif to
 extract fingerprints, please refer to the
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
index 2df4269..da5212b 100644
--- a/docs/source/installation.rst
+++ b/docs/source/installation.rst
@@ -30,5 +30,9 @@ Alternatively, you can install ProLIF with pip instead of conda::
 
     pip install rdkit prolif
 
+To run the tutorials, you can install the optional dependencies with::
+
+    pip install prolif[tutorials]
+
 
 .. _conda: https://docs.conda.io/projects/conda/en/latest/user-guide/index.html
diff --git a/docs/source/modules/plotting.rst b/docs/source/modules/plotting.rst
index b0b9d35..60e2b74 100644
--- a/docs/source/modules/plotting.rst
+++ b/docs/source/modules/plotting.rst
@@ -1,8 +1,10 @@
 Plotting
 ========
 
-.. automodule:: prolif.plotting.network
+.. automodule:: prolif.plotting.residues
 
 .. automodule:: prolif.plotting.barcode
 
-.. automodule:: prolif.plotting.complex3d
\ No newline at end of file
+.. automodule:: prolif.plotting.complex3d
+
+.. automodule:: prolif.plotting.network
diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst
new file mode 100644
index 0000000..8583aed
--- /dev/null
+++ b/docs/source/tutorials.rst
@@ -0,0 +1,65 @@
+Tutorials
+=========
+
+In this section you will learn how to generate an interaction fingerprint for different
+scenarios, how to extract information from the fingerprint data, and how to visualize
+it.
+
+You can get a very brief overview of how ProLIF is implemented in the
+:ref:`source/api:Summary` section of this documentation, with links to the different
+modules for a more in-depth explanation.
+
+All the examples here showcase a transmembrane protein (GPCR family) in complex with a
+small molecule. For convenience, the different files used for this tutorial are included
+with the ProLIF installation, and you can access these through either:
+
+- ``prolif.datafiles.TOP`` and ``prolif.datafiles.TRAJ`` for the topology (PDB file) and
+  trajectory (XTC file) used in the molecular dynamics tutorials.
+- ``prolif.datafiles.datapath``, which points to the directory containing all the other
+  tutorial files. This is a :class:`pathlib.Path` object which offers a convenient way
+  to navigate filesystems as shown in their
+  `basic use <https://docs.python.org/3/library/pathlib.html#basic-use>`__ section.
+
+.. warning::
+    Outside of the tutorials, remember to switch any reference to ``prolif.datafiles``
+    with the actual paths to your inputs.
+
+.. tip::
+    At the top of each tutorial's page you can find links to either download the
+    notebook or run it in Google Colab. You can install the dependencies for the
+    tutorials with the command::
+      
+      pip install prolif[tutorials]
+
+
+Molecular dynamics
+------------------
+
+There are two tutorial notebooks for MD simulations, depending on the type of components
+that are being analyzed:
+
+- :ref:`notebooks/md-ligand-protein:Ligand-protein MD`
+- :ref:`notebooks/md-protein-protein:Protein-protein MD`
+
+Docking
+-------
+
+Follow this tutorial for docking of a ligand with a protein:
+
+:ref:`notebooks/docking:Docking`
+
+PDB file
+--------
+
+This tutorial showcases how to use ProLIF from a PDB file:
+
+:ref:`notebooks/pdb:PDB`
+
+Advanced usage
+--------------
+
+For more advanced usecases such as modifying interaction parameters, defining your own
+interactions, ignoring backbone interactions and such, you can have a look at this
+tutorial:
+
+:ref:`notebooks/advanced:Advanced usage`
diff --git a/.binder/environment.yml b/environment.yml
similarity index 59%
rename from .binder/environment.yml
rename to environment.yml
index db8e919..253f86f 100644
--- a/.binder/environment.yml
+++ b/environment.yml
@@ -9,17 +9,21 @@ dependencies:
   - ipywidgets
   - pandoc
   - pip
-  - nbsphinx==0.8.7
+  - myst-nb<1.0
+  - mdit-py-plugins
   - networkx
   - recommonmark
-  - matplotlib==3.5.2
-  - sphinx==4.1.2
+  - matplotlib>=3.7
+  - seaborn<1.0
+  - sphinx==5.3.0
+  - sphinx-copybutton
   - tqdm
-  - python==3.9.*
-  - mdanalysis==2.2.*
+  - python==3.11.*
+  - mdanalysis
   - rdkit
   - dill
   - multiprocess
   - pip:
     - pyvis
     - py3Dmol
+    - sphinx-book-theme
diff --git a/prolif/__init__.py b/prolif/__init__.py
index ddf957f..0e69964 100644
--- a/prolif/__init__.py
+++ b/prolif/__init__.py
@@ -2,5 +2,6 @@ from prolif import datafiles
 from prolif._version import __version__
 from prolif.fingerprint import Fingerprint
 from prolif.molecule import Molecule, mol2_supplier, pdbqt_supplier, sdf_supplier
+from prolif.plotting.residues import display_residues
 from prolif.residue import ResidueId
 from prolif.utils import get_residues_near_ligand, to_bitvectors, to_dataframe
diff --git a/prolif/fingerprint.py b/prolif/fingerprint.py
index 56deedc..88ac7ae 100644
--- a/prolif/fingerprint.py
+++ b/prolif/fingerprint.py
@@ -24,6 +24,7 @@ Calculate a Protein-Ligand Interaction Fingerprint --- :mod:`prolif.fingerprint`
     fp = prolif.Fingerprint.from_pickle("fingerprint.pkl")
 
 """
+import warnings
 from collections.abc import Sized
 from functools import wraps
 from typing import Literal, Optional, Tuple
@@ -479,7 +480,7 @@ class Fingerprint:
         for ts in iterator:
             lig_mol = Molecule.from_mda(lig, **converter_kwargs[0])
             prot_mol = Molecule.from_mda(prot, **converter_kwargs[1])
-            ifp[ts.frame] = self.generate(
+            ifp[int(ts.frame)] = self.generate(
                 lig_mol, prot_mol, residues=residues, metadata=True
             )
         self.ifp = ifp
@@ -809,10 +810,15 @@ class Fingerprint:
             Switched to dill instead of pickle
 
         """
-        if isinstance(path_or_bytes, bytes):
-            return dill.loads(path_or_bytes)
-        with open(path_or_bytes, "rb") as f:
-            return dill.load(f)
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                "ignore",
+                r"The .+ interaction has been superseded by a new class",  # pragma: no cover
+            )
+            if isinstance(path_or_bytes, bytes):
+                return dill.loads(path_or_bytes)
+            with open(path_or_bytes, "rb") as f:
+                return dill.load(f)
 
     def plot_lignetwork(
         self,
diff --git a/prolif/parallel.py b/prolif/parallel.py
index cb455ab..3705087 100644
--- a/prolif/parallel.py
+++ b/prolif/parallel.py
@@ -143,7 +143,7 @@ class TrajectoryPool:
             data = cls.fp.generate(
                 lig_mol, prot_mol, residues=cls.residues, metadata=True
             )
-            ifp[ts.frame] = data
+            ifp[int(ts.frame)] = data
             with cls.tracker.get_lock():
                 cls.tracker.value += 1
         return ifp
diff --git a/prolif/plotting/barcode.py b/prolif/plotting/barcode.py
index dde1c6b..2d4f9a2 100644
--- a/prolif/plotting/barcode.py
+++ b/prolif/plotting/barcode.py
@@ -179,7 +179,11 @@ class Barcode:
 
         # legend
         values: List[int] = np.unique(self.df.values).tolist()
-        values.pop(values.index(0))  # remove None color
+        try:
+            values.pop(values.index(0))  # remove None color
+        except ValueError:
+            # 0 not in values (e.g. plotting a single frame)
+            pass
         legend_colors = {
             self.inv_color_mapper[value]: im.cmap(value) for value in values
         }
diff --git a/prolif/plotting/complex3d.py b/prolif/plotting/complex3d.py
index 50f158a..75802ae 100644
--- a/prolif/plotting/complex3d.py
+++ b/prolif/plotting/complex3d.py
@@ -49,6 +49,11 @@ class Complex3D:
         interactions.
     PROTEIN_STYLE : Dict[str, Dict] = {"cartoon": {"style": "edged"}}
         Style object passed to ``3Dmol.js`` for the entire protein.
+    PEPTIDE_STYLE : Dict[str, Dict] = "cartoon": {"style": "edged", "colorscheme": "cyanCarbon"}
+        Style object passed to ``3Dmol.js`` for the ligand as a peptide if appropriate.
+    PEPTIDE_THRESHOLD : int = 2
+        Ligands with this number of residues or more will be displayed using
+        ``PEPTIDE_STYLE`` in addition to the ``LIGAND_STYLE``.
     LIGAND_DISPLAYED_ATOM : Dict[str, int]
         Which atom should be used to display an atom-to-atom interaction for the ligand.
         Refers to the order defined in the SMARTS pattern used in interaction
@@ -73,6 +78,10 @@ class Complex3D:
     LIGAND_STYLE: ClassVar[Dict] = {"stick": {"colorscheme": "cyanCarbon"}}
     RESIDUES_STYLE: ClassVar[Dict] = {"stick": {}}
     PROTEIN_STYLE: ClassVar[Dict] = {"cartoon": {"style": "edged"}}
+    PEPTIDE_STYLE: ClassVar[Dict] = {
+        "cartoon": {"style": "edged", "colorscheme": "cyanCarbon"}
+    }
+    PEPTIDE_THRESHOLD: ClassVar[int] = 5
     LIGAND_DISPLAYED_ATOM = {
         "HBDonor": 1,
         "XBDonor": 1,
@@ -273,5 +282,13 @@ class Complex3D:
         model = v.getModel()
         model.setStyle({}, self.PROTEIN_STYLE)
 
+        # do the same for ligand if multiple residues
+        if self.lig_mol.n_residues >= self.PEPTIDE_THRESHOLD:
+            mol = Chem.RemoveAllHs(self.lig_mol)
+            pdb = Chem.MolToPDBBlock(mol, flavor=0x20 | 0x10)
+            v.addModel(pdb, "pdb")
+            model = v.getModel()
+            model.setStyle({}, self.PEPTIDE_STYLE)
+
         v.zoomTo({"model": list(models.values())})
         return v
diff --git a/prolif/plotting/residues.py b/prolif/plotting/residues.py
new file mode 100644
index 0000000..7ae68b5
--- /dev/null
+++ b/prolif/plotting/residues.py
@@ -0,0 +1,64 @@
+"""
+Plot residues --- :mod:`prolif.plotting.residues`
+=================================================
+
+.. versionadded:: 2.0.0
+
+.. autofunction:: display_residues
+
+"""
+from typing import Any, Optional, Tuple
+
+from rdkit import Chem
+from rdkit.Chem import Draw
+
+from prolif.molecule import Molecule
+
+
+def display_residues(
+    mol: Molecule,
+    residues_slice: Optional[slice] = None,
+    *,
+    size: Tuple[int, int] = (200, 140),
+    mols_per_row: int = 4,
+    use_svg: bool = True,
+) -> Any:
+    """Display a grid image of the residues in the molecule. The hydrogens are stripped
+    and the 3D coordinates removed for a clearer visualisation.
+
+    Parameters
+    ----------
+    mol: prolif.Molecule
+        The molecule to show residues from.
+    residues_slice: Optional[slice] = None
+        Optionally, a slice of residues to display, e.g. ``slice(20)`` for the first 20
+        residues, or ``slice(<start>, <stop>, <step>)`` for a more complex selection.
+    size: Tuple[int, int] = (200, 140)
+        Size of each residue image.
+    mols_per_row: int = 4
+        Number of residues displayed per row.
+    use_svg: bool = True
+        Generate an SVG or PNG image.
+    """
+    frags = []
+    residues_iterable = (
+        mol if residues_slice is None else mol.residues.select(residues_slice).values()
+    )
+    ipython_kwargs = (
+        {"maxMols": mol.n_residues} if hasattr(Chem.Mol, "_repr_svg_") else {}
+    )
+
+    for residue in residues_iterable:
+        resmol = Chem.RemoveHs(residue)
+        resmol.RemoveAllConformers()
+        resmol.SetProp("_Name", str(residue.resid))
+        frags.append(resmol)
+
+    return Draw.MolsToGridImage(
+        frags,
+        legends=[mol.GetProp("_Name") for mol in frags],
+        subImgSize=size,
+        molsPerRow=mols_per_row,
+        useSVG=use_svg,
+        **ipython_kwargs,
+    )
diff --git a/pyproject.toml b/pyproject.toml
index 049809d..107c1a5 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -46,6 +46,14 @@ plots = [
     "py3Dmol",
     "matplotlib",
 ]
+tutorials = [
+    "rdkit",
+    "matplotlib>=3.7",
+    "seaborn<1.0",
+    "networkx",
+    "pyvis",
+    "prolif[plots]",
+]
 build = [
     "build",
 ]


Test patch:
diff --git a/tests/plotting/test_residues.py b/tests/plotting/test_residues.py
new file mode 100644
index 0000000..a889519
--- /dev/null
+++ b/tests/plotting/test_residues.py
@@ -0,0 +1,28 @@
+import pytest
+from PIL.PngImagePlugin import PngImageFile
+
+import prolif as plf
+
+
+@pytest.mark.parametrize(
+    "kwargs, expected_size",
+    [
+        ({"residues_slice": slice(10)}, "width='800px' height='420px'"),
+        ({"residues_slice": slice(10, 20)}, "width='800px' height='420px'"),
+        ({"residues_slice": slice(0, 20, 2)}, "width='800px' height='420px'"),
+        ({"size": (100, 100)}, "width='400px' height='1400px'"),
+        ({"mols_per_row": 2}, "width='400px' height='3920px'"),
+    ],
+)
+def test_display_residues_svg(
+    protein_mol: plf.Molecule, kwargs: dict, expected_size: str
+) -> None:
+    img = plf.display_residues(protein_mol, **kwargs)
+    assert "<svg" in img
+    assert expected_size in img
+
+
+def test_display_residues_png(protein_mol: plf.Molecule) -> None:
+    img = plf.display_residues(protein_mol, slice(10), use_svg=False)
+    assert isinstance(img, PngImageFile)
+    assert img.size == (800, 420)