From 4c102b5e184398b5ca71b80df39ca0429b7ce219 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 11 Feb 2026 10:31:26 +0000
Subject: [PATCH 01/10] Initial plan


From e5bc605ddbc517f2adf9e1b20bbdd11558cfb2cb Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 11 Feb 2026 10:35:23 +0000
Subject: [PATCH 02/10] Add tutorial on how to use 3c composite methods

Create a new documentation page explaining how to use 3c composite
methods (B97-3c, PBEh-3c, HF-3c) with the s-dftd3 Python API,
covering D3 dispersion correction, geometric counter-poise (gCP)
correction, and how to combine both for the full 3c energy and
gradient.

Co-authored-by: awvwgk <28669218+awvwgk@users.noreply.github.com>
---
 doc/tutorial/index.rst            |   1 +
 doc/tutorial/using-3c-methods.rst | 228 ++++++++++++++++++++++++++++++
 2 files changed, 229 insertions(+)
 create mode 100644 doc/tutorial/using-3c-methods.rst

diff --git a/doc/tutorial/index.rst b/doc/tutorial/index.rst
index cc56fd62..19c19b3e 100644
--- a/doc/tutorial/index.rst
+++ b/doc/tutorial/index.rst
@@ -7,3 +7,4 @@ This section contains tutorials on using D3.
 
    first-steps-cli
    first-steps-fortran
+   using-3c-methods
diff --git a/doc/tutorial/using-3c-methods.rst b/doc/tutorial/using-3c-methods.rst
new file mode 100644
index 00000000..945e2d5c
--- /dev/null
+++ b/doc/tutorial/using-3c-methods.rst
@@ -0,0 +1,228 @@
+Using 3c composite methods
+==========================
+
+In this tutorial we will learn how to use the 3c family of composite electronic structure methods with the ``dftd3`` Python API.
+Composite methods like B97-3c\ :footcite:`brandenburg2018` and PBEh-3c\ :footcite:`grimme2015` combine a density functional with a specific basis set, a D3 dispersion correction, and a geometric counter-poise (gCP) correction\ :footcite:`kruse2012` to correct for basis set superposition error (BSSE).
+
+
+Background on 3c methods
+------------------------
+
+The 3c methods are designed as efficient and accurate composite electronic structure methods.
+The name "3c" refers to the three corrections applied on top of a base density functional:
+
+1. **D3 dispersion correction** – Accounts for London dispersion interactions using the DFT-D3 model with specifically fitted damping parameters.
+2. **Geometric counter-poise (gCP) correction** – Corrects the basis set superposition error that arises from using small basis sets.
+3. **Short-range basis set (SRB) correction** – In some methods an additional short-range bond length correction is applied to improve covalent bond lengths.
+
+The total energy of a 3c method is computed as
+
+.. math::
+
+   E_\text{3c} = E_\text{DFT/basis} + E_\text{D3} + E_\text{gCP}
+
+where :math:`E_\text{DFT/basis}` is the SCF energy computed with the corresponding density functional and basis set, :math:`E_\text{D3}` is the DFT-D3 dispersion energy, and :math:`E_\text{gCP}` is the geometric counter-poise correction energy.
+
+
+Computing the D3 contribution
+------------------------------
+
+The D3 dispersion correction is the first component we need.
+Each 3c method uses a specific damping scheme and parameters.
+For example, PBEh-3c uses rational (BJ) damping with ``s8=0.0``, while B97-3c uses rational damping with its own set of parameters.
+
+Using the Python API we can easily obtain the D3 energy for a given molecule.
+In this example we will use a caffeine molecule:
+
+.. code-block:: python
+
+   import numpy as np
+   from dftd3.interface import DispersionModel, RationalDampingParam
+
+   # Caffeine molecule (atomic numbers and positions in Bohr)
+   numbers = np.array(
+       [6, 7, 6, 7, 6, 6, 6, 8, 7, 6, 8, 7, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
+   )
+   positions = np.array([
+       [+2.02799738646442, +0.09231312124713, -0.14310895950963],
+       [+4.75011007621000, +0.02373496014051, -0.14324124033844],
+       [+6.33434307654413, +2.07098865582721, -0.14235306905930],
+       [+8.72860718071825, +1.38002919517619, -0.14265542523943],
+       [+8.65318821103610, -1.19324866489847, -0.14231527453678],
+       [+6.23857175648671, -2.08353643730276, -0.14218299370797],
+       [+5.63266886875962, -4.69950321056008, -0.13940509630299],
+       [+3.44931709749015, -5.48092386085491, -0.14318454855466],
+       [+7.77508917214943, -5.24064112783473, -0.13206210149840],
+       [+9.68504443463980, -3.43480556543577, -0.13376503914194],
+       [+11.7766597740572, -2.85589272667498, -0.13347836327959],
+       [+0.71292821912498, -1.81184541295565, -0.14404507712755],
+       [-1.07804988915028, -0.36933811262178, -0.14399838668498],
+       [+9.84554065797340, -6.86700842661498, -0.13277505395063],
+       [-2.48328028863736, -1.73067674389689, -0.14259442617502],
+       [-1.34385752710948, +0.57786478045678, +1.58564153988000],
+       [-1.34425478091498, +0.57834894571524, -1.87640258713638],
+       [+8.96081504094682, -8.37942090821983, +1.00589803206426],
+       [+11.5017680845878, -6.28787412600376, -0.13181456387625],
+       [+8.96116647088646, -8.37865627213056, -1.27402975539366],
+       [+1.52389863680920, +4.42491864930852, +1.59083729873498],
+       [+1.52355529826621, +4.42437909527069, -1.87861372682396],
+       [-0.43240687424107, +5.48666820368743, -0.14223727178700],
+       [+3.59762505724975, +5.19469189498972, -0.14192321862387],
+   ])
+
+   # PBEh-3c D3 parameters (rational damping)
+   model = DispersionModel(numbers, positions)
+   res = model.get_dispersion(RationalDampingParam(method="pbeh3c"), grad=False)
+   print(f"D3 dispersion energy: {res['energy']:16.10f} Hartree")
+
+Since the D3 parameters are already stored in the library for the 3c methods, we can simply pass the method name (e.g. ``"pbeh3c"`` or ``"b973c"``) to ``RationalDampingParam``.
+
+.. note::
+
+   B97-3c also has parameters for zero damping.
+   When using it as part of the 3c composite method, the rational damping (BJ) parameters should be used.
+
+
+Computing the gCP contribution
+-------------------------------
+
+The geometric counter-poise (gCP) correction accounts for the basis set superposition error (BSSE) that occurs with small or medium-sized basis sets.
+Using the Python API, we can compute the gCP energy using the ``GeometricCounterpoise`` class.
+
+.. code-block:: python
+
+   from dftd3.interface import GeometricCounterpoise
+
+   # Create gCP object for a 3c method
+   # For B97-3c, the method name alone is sufficient
+   gcp = GeometricCounterpoise(
+       numbers,
+       positions,
+       method="b973c",
+   )
+
+   res = gcp.get_counterpoise(grad=False)
+   print(f"gCP energy: {res['energy']:16.10f} Hartree")
+
+For methods like PBEh-3c which use a named basis set, both the method and basis set name can be specified:
+
+.. code-block:: python
+
+   gcp = GeometricCounterpoise(
+       numbers,
+       positions,
+       method="pbeh3c",
+   )
+
+   res = gcp.get_counterpoise(grad=False)
+   print(f"gCP energy: {res['energy']:16.10f} Hartree")
+
+The ``GeometricCounterpoise`` class loads the appropriate gCP parameters internally, including the short-range bond (SRB) correction where applicable.
+
+
+Putting it all together
+-----------------------
+
+To assemble the full 3c composite energy, we combine both the D3 dispersion correction and the gCP correction with the SCF energy from a DFT calculation.
+
+Here is a complete example computing both correction terms for the B97-3c method:
+
+.. code-block:: python
+
+   import numpy as np
+   from dftd3.interface import (
+       DispersionModel,
+       RationalDampingParam,
+       GeometricCounterpoise,
+   )
+
+   numbers = np.array(
+       [6, 7, 6, 7, 6, 6, 6, 8, 7, 6, 8, 7, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
+   )
+   positions = np.array([
+       [+2.02799738646442, +0.09231312124713, -0.14310895950963],
+       [+4.75011007621000, +0.02373496014051, -0.14324124033844],
+       [+6.33434307654413, +2.07098865582721, -0.14235306905930],
+       [+8.72860718071825, +1.38002919517619, -0.14265542523943],
+       [+8.65318821103610, -1.19324866489847, -0.14231527453678],
+       [+6.23857175648671, -2.08353643730276, -0.14218299370797],
+       [+5.63266886875962, -4.69950321056008, -0.13940509630299],
+       [+3.44931709749015, -5.48092386085491, -0.14318454855466],
+       [+7.77508917214943, -5.24064112783473, -0.13206210149840],
+       [+9.68504443463980, -3.43480556543577, -0.13376503914194],
+       [+11.7766597740572, -2.85589272667498, -0.13347836327959],
+       [+0.71292821912498, -1.81184541295565, -0.14404507712755],
+       [-1.07804988915028, -0.36933811262178, -0.14399838668498],
+       [+9.84554065797340, -6.86700842661498, -0.13277505395063],
+       [-2.48328028863736, -1.73067674389689, -0.14259442617502],
+       [-1.34385752710948, +0.57786478045678, +1.58564153988000],
+       [-1.34425478091498, +0.57834894571524, -1.87640258713638],
+       [+8.96081504094682, -8.37942090821983, +1.00589803206426],
+       [+11.5017680845878, -6.28787412600376, -0.13181456387625],
+       [+8.96116647088646, -8.37865627213056, -1.27402975539366],
+       [+1.52389863680920, +4.42491864930852, +1.59083729873498],
+       [+1.52355529826621, +4.42437909527069, -1.87861372682396],
+       [-0.43240687424107, +5.48666820368743, -0.14223727178700],
+       [+3.59762505724975, +5.19469189498972, -0.14192321862387],
+   ])
+
+   # 1. Compute D3 dispersion correction
+   model = DispersionModel(numbers, positions)
+   d3_res = model.get_dispersion(RationalDampingParam(method="b973c"), grad=False)
+   d3_energy = d3_res["energy"]
+
+   # 2. Compute gCP correction
+   gcp = GeometricCounterpoise(numbers, positions, method="b973c")
+   gcp_res = gcp.get_counterpoise(grad=False)
+   gcp_energy = gcp_res["energy"]
+
+   print(f"D3 dispersion energy:  {d3_energy:16.10f} Hartree")
+   print(f"gCP correction energy: {gcp_energy:16.10f} Hartree")
+   print(f"Total correction:      {d3_energy + gcp_energy:16.10f} Hartree")
+
+.. note::
+
+   The SCF energy :math:`E_\text{DFT/basis}` must be computed separately using a quantum chemistry program (such as PySCF, ORCA, or Turbomole) with the correct functional and basis set.
+   The total 3c composite energy is then:
+
+   :math:`E_\text{3c} = E_\text{SCF} + E_\text{D3} + E_\text{gCP}`
+
+
+Including gradients
+-------------------
+
+For geometry optimizations, both the D3 and gCP gradients need to be included.
+To request gradients, pass ``grad=True``:
+
+.. code-block:: python
+
+   # D3 gradient
+   d3_res = model.get_dispersion(RationalDampingParam(method="b973c"), grad=True)
+   d3_gradient = d3_res["gradient"]
+
+   # gCP gradient
+   gcp_res = gcp.get_counterpoise(grad=True)
+   gcp_gradient = gcp_res["gradient"]
+
+   # Total correction gradient
+   total_gradient = d3_gradient + gcp_gradient
+
+The gradient arrays have the shape ``(natoms, 3)`` and are in atomic units (Hartree/Bohr).
+
+
+Summary
+-------
+
+In this tutorial we learned how to
+
+- compute the D3 dispersion energy for 3c composite methods using stored parameters
+- compute the geometric counter-poise (gCP) correction for basis set superposition error
+- combine both corrections to obtain the total correction energy and gradient for 3c methods
+
+The available 3c methods with built-in parameter support include B97-3c, PBEh-3c, and HF-3c.
+
+
+Literature
+----------
+
+.. footbibliography::

From abd4b5b515453a69546be78e4d1803df59ae63da Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 11 Feb 2026 10:48:15 +0000
Subject: [PATCH 03/10] Convert 3c methods tutorial to Jupyter notebook with
 myst_nb

Replace the RST tutorial with a Jupyter notebook that uses MyST
Markdown syntax for rich content (math, notes, citations). Add
myst_nb Sphinx extension for rendering notebooks in documentation,
with nb_execution_mode=off since the library is not available during
doc builds.

Co-authored-by: awvwgk <28669218+awvwgk@users.noreply.github.com>
---
 doc/conf.py                         |   4 +
 doc/requirements.txt                |   1 +
 doc/tutorial/using-3c-methods.ipynb | 328 ++++++++++++++++++++++++++++
 doc/tutorial/using-3c-methods.rst   | 228 -------------------
 4 files changed, 333 insertions(+), 228 deletions(-)
 create mode 100644 doc/tutorial/using-3c-methods.ipynb
 delete mode 100644 doc/tutorial/using-3c-methods.rst

diff --git a/doc/conf.py b/doc/conf.py
index a1288084..5af12d9e 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -13,6 +13,7 @@
 release = version
 
 extensions = [
+    "myst_nb",
     "sphinx_design",
     "sphinx_copybutton",
     "sphinx.ext.autosummary",
@@ -25,6 +26,9 @@
     "sphinxcontrib.bibtex",
 ]
 
+nb_execution_mode = "off"
+myst_enable_extensions = ["dollarmath"]
+
 html_theme = "sphinx_book_theme"
 html_title = "Simple DFT-D3"
 
diff --git a/doc/requirements.txt b/doc/requirements.txt
index 512fb2a7..29b471ef 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -2,6 +2,7 @@ sphinx-book-theme
 sphinx-copybutton
 sphinx-design
 sphinxcontrib-bibtex
+myst-nb
 numpy
 ase
 cffi
diff --git a/doc/tutorial/using-3c-methods.ipynb b/doc/tutorial/using-3c-methods.ipynb
new file mode 100644
index 00000000..4daba8b4
--- /dev/null
+++ b/doc/tutorial/using-3c-methods.ipynb
@@ -0,0 +1,328 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Using 3c composite methods\n",
+    "\n",
+    "In this tutorial we will learn how to use the 3c family of composite electronic structure methods with the `dftd3` Python API.\n",
+    "Composite methods like B97-3c {footcite}`brandenburg2018` and PBEh-3c {footcite}`grimme2015` combine a density functional with a specific basis set, a D3 dispersion correction, and a geometric counter-poise (gCP) correction {footcite}`kruse2012` to correct for basis set superposition error (BSSE)."
+   ],
+   "id": "cell-0"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Background on 3c methods\n",
+    "\n",
+    "The 3c methods are designed as efficient and accurate composite electronic structure methods.\n",
+    "The name \"3c\" refers to the three corrections applied on top of a base density functional:\n",
+    "\n",
+    "1. **D3 dispersion correction** \u2013 Accounts for London dispersion interactions using the DFT-D3 model with specifically fitted damping parameters.\n",
+    "2. **Geometric counter-poise (gCP) correction** \u2013 Corrects the basis set superposition error that arises from using small basis sets.\n",
+    "3. **Short-range basis set (SRB) correction** \u2013 In some methods an additional short-range bond length correction is applied to improve covalent bond lengths.\n",
+    "\n",
+    "The total energy of a 3c method is computed as\n",
+    "\n",
+    "$$\n",
+    "E_\\text{3c} = E_\\text{DFT/basis} + E_\\text{D3} + E_\\text{gCP}\n",
+    "$$\n",
+    "\n",
+    "where $E_\\text{DFT/basis}$ is the SCF energy computed with the corresponding density functional and basis set, $E_\\text{D3}$ is the DFT-D3 dispersion energy, and $E_\\text{gCP}$ is the geometric counter-poise correction energy."
+   ],
+   "id": "cell-1"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Computing the D3 contribution\n",
+    "\n",
+    "The D3 dispersion correction is the first component we need.\n",
+    "Each 3c method uses a specific damping scheme and parameters.\n",
+    "For example, PBEh-3c uses rational (BJ) damping with `s8=0.0`, while B97-3c uses rational damping with its own set of parameters.\n",
+    "\n",
+    "Using the Python API we can easily obtain the D3 energy for a given molecule.\n",
+    "In this example we will use a caffeine molecule:"
+   ],
+   "id": "cell-2"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from dftd3.interface import DispersionModel, RationalDampingParam\n",
+    "\n",
+    "# Caffeine molecule (atomic numbers and positions in Bohr)\n",
+    "numbers = np.array(\n",
+    "    [6, 7, 6, 7, 6, 6, 6, 8, 7, 6, 8, 7, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
+    ")\n",
+    "positions = np.array([\n",
+    "    [+2.02799738646442, +0.09231312124713, -0.14310895950963],\n",
+    "    [+4.75011007621000, +0.02373496014051, -0.14324124033844],\n",
+    "    [+6.33434307654413, +2.07098865582721, -0.14235306905930],\n",
+    "    [+8.72860718071825, +1.38002919517619, -0.14265542523943],\n",
+    "    [+8.65318821103610, -1.19324866489847, -0.14231527453678],\n",
+    "    [+6.23857175648671, -2.08353643730276, -0.14218299370797],\n",
+    "    [+5.63266886875962, -4.69950321056008, -0.13940509630299],\n",
+    "    [+3.44931709749015, -5.48092386085491, -0.14318454855466],\n",
+    "    [+7.77508917214943, -5.24064112783473, -0.13206210149840],\n",
+    "    [+9.68504443463980, -3.43480556543577, -0.13376503914194],\n",
+    "    [+11.7766597740572, -2.85589272667498, -0.13347836327959],\n",
+    "    [+0.71292821912498, -1.81184541295565, -0.14404507712755],\n",
+    "    [-1.07804988915028, -0.36933811262178, -0.14399838668498],\n",
+    "    [+9.84554065797340, -6.86700842661498, -0.13277505395063],\n",
+    "    [-2.48328028863736, -1.73067674389689, -0.14259442617502],\n",
+    "    [-1.34385752710948, +0.57786478045678, +1.58564153988000],\n",
+    "    [-1.34425478091498, +0.57834894571524, -1.87640258713638],\n",
+    "    [+8.96081504094682, -8.37942090821983, +1.00589803206426],\n",
+    "    [+11.5017680845878, -6.28787412600376, -0.13181456387625],\n",
+    "    [+8.96116647088646, -8.37865627213056, -1.27402975539366],\n",
+    "    [+1.52389863680920, +4.42491864930852, +1.59083729873498],\n",
+    "    [+1.52355529826621, +4.42437909527069, -1.87861372682396],\n",
+    "    [-0.43240687424107, +5.48666820368743, -0.14223727178700],\n",
+    "    [+3.59762505724975, +5.19469189498972, -0.14192321862387],\n",
+    "])\n",
+    "\n",
+    "# PBEh-3c D3 parameters (rational damping)\n",
+    "model = DispersionModel(numbers, positions)\n",
+    "res = model.get_dispersion(RationalDampingParam(method=\"pbeh3c\"), grad=False)\n",
+    "print(f\"D3 dispersion energy: {res['energy']:16.10f} Hartree\")"
+   ],
+   "id": "cell-3"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Since the D3 parameters are already stored in the library for the 3c methods, we can simply pass the method name (e.g. `\"pbeh3c\"` or `\"b973c\"`) to `RationalDampingParam`.\n",
+    "\n",
+    "```{note}\n",
+    "B97-3c also has parameters for zero damping.\n",
+    "When using it as part of the 3c composite method, the rational damping (BJ) parameters should be used.\n",
+    "```"
+   ],
+   "id": "cell-4"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Computing the gCP contribution\n",
+    "\n",
+    "The geometric counter-poise (gCP) correction accounts for the basis set superposition error (BSSE) that occurs with small or medium-sized basis sets.\n",
+    "Using the Python API, we can compute the gCP energy using the `GeometricCounterpoise` class."
+   ],
+   "id": "cell-5"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from dftd3.interface import GeometricCounterpoise\n",
+    "\n",
+    "# Create gCP object for a 3c method\n",
+    "# For B97-3c, the method name alone is sufficient\n",
+    "gcp = GeometricCounterpoise(\n",
+    "    numbers,\n",
+    "    positions,\n",
+    "    method=\"b973c\",\n",
+    ")\n",
+    "\n",
+    "res = gcp.get_counterpoise(grad=False)\n",
+    "print(f\"gCP energy: {res['energy']:16.10f} Hartree\")"
+   ],
+   "id": "cell-6"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For methods like PBEh-3c which use a named basis set, both the method and basis set name can be specified:"
+   ],
+   "id": "cell-7"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gcp = GeometricCounterpoise(\n",
+    "    numbers,\n",
+    "    positions,\n",
+    "    method=\"pbeh3c\",\n",
+    ")\n",
+    "\n",
+    "res = gcp.get_counterpoise(grad=False)\n",
+    "print(f\"gCP energy: {res['energy']:16.10f} Hartree\")"
+   ],
+   "id": "cell-8"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `GeometricCounterpoise` class loads the appropriate gCP parameters internally, including the short-range bond (SRB) correction where applicable."
+   ],
+   "id": "cell-9"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Putting it all together\n",
+    "\n",
+    "To assemble the full 3c composite energy, we combine both the D3 dispersion correction and the gCP correction with the SCF energy from a DFT calculation.\n",
+    "\n",
+    "Here is a complete example computing both correction terms for the B97-3c method:"
+   ],
+   "id": "cell-10"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from dftd3.interface import (\n",
+    "    DispersionModel,\n",
+    "    RationalDampingParam,\n",
+    "    GeometricCounterpoise,\n",
+    ")\n",
+    "\n",
+    "numbers = np.array(\n",
+    "    [6, 7, 6, 7, 6, 6, 6, 8, 7, 6, 8, 7, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
+    ")\n",
+    "positions = np.array([\n",
+    "    [+2.02799738646442, +0.09231312124713, -0.14310895950963],\n",
+    "    [+4.75011007621000, +0.02373496014051, -0.14324124033844],\n",
+    "    [+6.33434307654413, +2.07098865582721, -0.14235306905930],\n",
+    "    [+8.72860718071825, +1.38002919517619, -0.14265542523943],\n",
+    "    [+8.65318821103610, -1.19324866489847, -0.14231527453678],\n",
+    "    [+6.23857175648671, -2.08353643730276, -0.14218299370797],\n",
+    "    [+5.63266886875962, -4.69950321056008, -0.13940509630299],\n",
+    "    [+3.44931709749015, -5.48092386085491, -0.14318454855466],\n",
+    "    [+7.77508917214943, -5.24064112783473, -0.13206210149840],\n",
+    "    [+9.68504443463980, -3.43480556543577, -0.13376503914194],\n",
+    "    [+11.7766597740572, -2.85589272667498, -0.13347836327959],\n",
+    "    [+0.71292821912498, -1.81184541295565, -0.14404507712755],\n",
+    "    [-1.07804988915028, -0.36933811262178, -0.14399838668498],\n",
+    "    [+9.84554065797340, -6.86700842661498, -0.13277505395063],\n",
+    "    [-2.48328028863736, -1.73067674389689, -0.14259442617502],\n",
+    "    [-1.34385752710948, +0.57786478045678, +1.58564153988000],\n",
+    "    [-1.34425478091498, +0.57834894571524, -1.87640258713638],\n",
+    "    [+8.96081504094682, -8.37942090821983, +1.00589803206426],\n",
+    "    [+11.5017680845878, -6.28787412600376, -0.13181456387625],\n",
+    "    [+8.96116647088646, -8.37865627213056, -1.27402975539366],\n",
+    "    [+1.52389863680920, +4.42491864930852, +1.59083729873498],\n",
+    "    [+1.52355529826621, +4.42437909527069, -1.87861372682396],\n",
+    "    [-0.43240687424107, +5.48666820368743, -0.14223727178700],\n",
+    "    [+3.59762505724975, +5.19469189498972, -0.14192321862387],\n",
+    "])\n",
+    "\n",
+    "# 1. Compute D3 dispersion correction\n",
+    "model = DispersionModel(numbers, positions)\n",
+    "d3_res = model.get_dispersion(RationalDampingParam(method=\"b973c\"), grad=False)\n",
+    "d3_energy = d3_res[\"energy\"]\n",
+    "\n",
+    "# 2. Compute gCP correction\n",
+    "gcp = GeometricCounterpoise(numbers, positions, method=\"b973c\")\n",
+    "gcp_res = gcp.get_counterpoise(grad=False)\n",
+    "gcp_energy = gcp_res[\"energy\"]\n",
+    "\n",
+    "print(f\"D3 dispersion energy:  {d3_energy:16.10f} Hartree\")\n",
+    "print(f\"gCP correction energy: {gcp_energy:16.10f} Hartree\")\n",
+    "print(f\"Total correction:      {d3_energy + gcp_energy:16.10f} Hartree\")"
+   ],
+   "id": "cell-11"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "```{note}\n",
+    "The SCF energy $E_\\text{DFT/basis}$ must be computed separately using a quantum chemistry program (such as PySCF, ORCA, or Turbomole) with the correct functional and basis set.\n",
+    "The total 3c composite energy is then:\n",
+    "\n",
+    "$E_\\text{3c} = E_\\text{SCF} + E_\\text{D3} + E_\\text{gCP}$\n",
+    "```"
+   ],
+   "id": "cell-12"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Including gradients\n",
+    "\n",
+    "For geometry optimizations, both the D3 and gCP gradients need to be included.\n",
+    "To request gradients, pass `grad=True`:"
+   ],
+   "id": "cell-13"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# D3 gradient\n",
+    "d3_res = model.get_dispersion(RationalDampingParam(method=\"b973c\"), grad=True)\n",
+    "d3_gradient = d3_res[\"gradient\"]\n",
+    "\n",
+    "# gCP gradient\n",
+    "gcp_res = gcp.get_counterpoise(grad=True)\n",
+    "gcp_gradient = gcp_res[\"gradient\"]\n",
+    "\n",
+    "# Total correction gradient\n",
+    "total_gradient = d3_gradient + gcp_gradient"
+   ],
+   "id": "cell-14"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The gradient arrays have the shape `(natoms, 3)` and are in atomic units (Hartree/Bohr).\n",
+    "\n",
+    "## Summary\n",
+    "\n",
+    "In this tutorial we learned how to\n",
+    "\n",
+    "- compute the D3 dispersion energy for 3c composite methods using stored parameters\n",
+    "- compute the geometric counter-poise (gCP) correction for basis set superposition error\n",
+    "- combine both corrections to obtain the total correction energy and gradient for 3c methods\n",
+    "\n",
+    "The available 3c methods with built-in parameter support include B97-3c, PBEh-3c, and HF-3c.\n",
+    "\n",
+    "## Literature\n",
+    "\n",
+    "```{footbibliography}\n",
+    "```"
+   ],
+   "id": "cell-15"
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
\ No newline at end of file
diff --git a/doc/tutorial/using-3c-methods.rst b/doc/tutorial/using-3c-methods.rst
deleted file mode 100644
index 945e2d5c..00000000
--- a/doc/tutorial/using-3c-methods.rst
+++ /dev/null
@@ -1,228 +0,0 @@
-Using 3c composite methods
-==========================
-
-In this tutorial we will learn how to use the 3c family of composite electronic structure methods with the ``dftd3`` Python API.
-Composite methods like B97-3c\ :footcite:`brandenburg2018` and PBEh-3c\ :footcite:`grimme2015` combine a density functional with a specific basis set, a D3 dispersion correction, and a geometric counter-poise (gCP) correction\ :footcite:`kruse2012` to correct for basis set superposition error (BSSE).
-
-
-Background on 3c methods
-------------------------
-
-The 3c methods are designed as efficient and accurate composite electronic structure methods.
-The name "3c" refers to the three corrections applied on top of a base density functional:
-
-1. **D3 dispersion correction** – Accounts for London dispersion interactions using the DFT-D3 model with specifically fitted damping parameters.
-2. **Geometric counter-poise (gCP) correction** – Corrects the basis set superposition error that arises from using small basis sets.
-3. **Short-range basis set (SRB) correction** – In some methods an additional short-range bond length correction is applied to improve covalent bond lengths.
-
-The total energy of a 3c method is computed as
-
-.. math::
-
-   E_\text{3c} = E_\text{DFT/basis} + E_\text{D3} + E_\text{gCP}
-
-where :math:`E_\text{DFT/basis}` is the SCF energy computed with the corresponding density functional and basis set, :math:`E_\text{D3}` is the DFT-D3 dispersion energy, and :math:`E_\text{gCP}` is the geometric counter-poise correction energy.
-
-
-Computing the D3 contribution
-------------------------------
-
-The D3 dispersion correction is the first component we need.
-Each 3c method uses a specific damping scheme and parameters.
-For example, PBEh-3c uses rational (BJ) damping with ``s8=0.0``, while B97-3c uses rational damping with its own set of parameters.
-
-Using the Python API we can easily obtain the D3 energy for a given molecule.
-In this example we will use a caffeine molecule:
-
-.. code-block:: python
-
-   import numpy as np
-   from dftd3.interface import DispersionModel, RationalDampingParam
-
-   # Caffeine molecule (atomic numbers and positions in Bohr)
-   numbers = np.array(
-       [6, 7, 6, 7, 6, 6, 6, 8, 7, 6, 8, 7, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
-   )
-   positions = np.array([
-       [+2.02799738646442, +0.09231312124713, -0.14310895950963],
-       [+4.75011007621000, +0.02373496014051, -0.14324124033844],
-       [+6.33434307654413, +2.07098865582721, -0.14235306905930],
-       [+8.72860718071825, +1.38002919517619, -0.14265542523943],
-       [+8.65318821103610, -1.19324866489847, -0.14231527453678],
-       [+6.23857175648671, -2.08353643730276, -0.14218299370797],
-       [+5.63266886875962, -4.69950321056008, -0.13940509630299],
-       [+3.44931709749015, -5.48092386085491, -0.14318454855466],
-       [+7.77508917214943, -5.24064112783473, -0.13206210149840],
-       [+9.68504443463980, -3.43480556543577, -0.13376503914194],
-       [+11.7766597740572, -2.85589272667498, -0.13347836327959],
-       [+0.71292821912498, -1.81184541295565, -0.14404507712755],
-       [-1.07804988915028, -0.36933811262178, -0.14399838668498],
-       [+9.84554065797340, -6.86700842661498, -0.13277505395063],
-       [-2.48328028863736, -1.73067674389689, -0.14259442617502],
-       [-1.34385752710948, +0.57786478045678, +1.58564153988000],
-       [-1.34425478091498, +0.57834894571524, -1.87640258713638],
-       [+8.96081504094682, -8.37942090821983, +1.00589803206426],
-       [+11.5017680845878, -6.28787412600376, -0.13181456387625],
-       [+8.96116647088646, -8.37865627213056, -1.27402975539366],
-       [+1.52389863680920, +4.42491864930852, +1.59083729873498],
-       [+1.52355529826621, +4.42437909527069, -1.87861372682396],
-       [-0.43240687424107, +5.48666820368743, -0.14223727178700],
-       [+3.59762505724975, +5.19469189498972, -0.14192321862387],
-   ])
-
-   # PBEh-3c D3 parameters (rational damping)
-   model = DispersionModel(numbers, positions)
-   res = model.get_dispersion(RationalDampingParam(method="pbeh3c"), grad=False)
-   print(f"D3 dispersion energy: {res['energy']:16.10f} Hartree")
-
-Since the D3 parameters are already stored in the library for the 3c methods, we can simply pass the method name (e.g. ``"pbeh3c"`` or ``"b973c"``) to ``RationalDampingParam``.
-
-.. note::
-
-   B97-3c also has parameters for zero damping.
-   When using it as part of the 3c composite method, the rational damping (BJ) parameters should be used.
-
-
-Computing the gCP contribution
--------------------------------
-
-The geometric counter-poise (gCP) correction accounts for the basis set superposition error (BSSE) that occurs with small or medium-sized basis sets.
-Using the Python API, we can compute the gCP energy using the ``GeometricCounterpoise`` class.
-
-.. code-block:: python
-
-   from dftd3.interface import GeometricCounterpoise
-
-   # Create gCP object for a 3c method
-   # For B97-3c, the method name alone is sufficient
-   gcp = GeometricCounterpoise(
-       numbers,
-       positions,
-       method="b973c",
-   )
-
-   res = gcp.get_counterpoise(grad=False)
-   print(f"gCP energy: {res['energy']:16.10f} Hartree")
-
-For methods like PBEh-3c which use a named basis set, both the method and basis set name can be specified:
-
-.. code-block:: python
-
-   gcp = GeometricCounterpoise(
-       numbers,
-       positions,
-       method="pbeh3c",
-   )
-
-   res = gcp.get_counterpoise(grad=False)
-   print(f"gCP energy: {res['energy']:16.10f} Hartree")
-
-The ``GeometricCounterpoise`` class loads the appropriate gCP parameters internally, including the short-range bond (SRB) correction where applicable.
-
-
-Putting it all together
------------------------
-
-To assemble the full 3c composite energy, we combine both the D3 dispersion correction and the gCP correction with the SCF energy from a DFT calculation.
-
-Here is a complete example computing both correction terms for the B97-3c method:
-
-.. code-block:: python
-
-   import numpy as np
-   from dftd3.interface import (
-       DispersionModel,
-       RationalDampingParam,
-       GeometricCounterpoise,
-   )
-
-   numbers = np.array(
-       [6, 7, 6, 7, 6, 6, 6, 8, 7, 6, 8, 7, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
-   )
-   positions = np.array([
-       [+2.02799738646442, +0.09231312124713, -0.14310895950963],
-       [+4.75011007621000, +0.02373496014051, -0.14324124033844],
-       [+6.33434307654413, +2.07098865582721, -0.14235306905930],
-       [+8.72860718071825, +1.38002919517619, -0.14265542523943],
-       [+8.65318821103610, -1.19324866489847, -0.14231527453678],
-       [+6.23857175648671, -2.08353643730276, -0.14218299370797],
-       [+5.63266886875962, -4.69950321056008, -0.13940509630299],
-       [+3.44931709749015, -5.48092386085491, -0.14318454855466],
-       [+7.77508917214943, -5.24064112783473, -0.13206210149840],
-       [+9.68504443463980, -3.43480556543577, -0.13376503914194],
-       [+11.7766597740572, -2.85589272667498, -0.13347836327959],
-       [+0.71292821912498, -1.81184541295565, -0.14404507712755],
-       [-1.07804988915028, -0.36933811262178, -0.14399838668498],
-       [+9.84554065797340, -6.86700842661498, -0.13277505395063],
-       [-2.48328028863736, -1.73067674389689, -0.14259442617502],
-       [-1.34385752710948, +0.57786478045678, +1.58564153988000],
-       [-1.34425478091498, +0.57834894571524, -1.87640258713638],
-       [+8.96081504094682, -8.37942090821983, +1.00589803206426],
-       [+11.5017680845878, -6.28787412600376, -0.13181456387625],
-       [+8.96116647088646, -8.37865627213056, -1.27402975539366],
-       [+1.52389863680920, +4.42491864930852, +1.59083729873498],
-       [+1.52355529826621, +4.42437909527069, -1.87861372682396],
-       [-0.43240687424107, +5.48666820368743, -0.14223727178700],
-       [+3.59762505724975, +5.19469189498972, -0.14192321862387],
-   ])
-
-   # 1. Compute D3 dispersion correction
-   model = DispersionModel(numbers, positions)
-   d3_res = model.get_dispersion(RationalDampingParam(method="b973c"), grad=False)
-   d3_energy = d3_res["energy"]
-
-   # 2. Compute gCP correction
-   gcp = GeometricCounterpoise(numbers, positions, method="b973c")
-   gcp_res = gcp.get_counterpoise(grad=False)
-   gcp_energy = gcp_res["energy"]
-
-   print(f"D3 dispersion energy:  {d3_energy:16.10f} Hartree")
-   print(f"gCP correction energy: {gcp_energy:16.10f} Hartree")
-   print(f"Total correction:      {d3_energy + gcp_energy:16.10f} Hartree")
-
-.. note::
-
-   The SCF energy :math:`E_\text{DFT/basis}` must be computed separately using a quantum chemistry program (such as PySCF, ORCA, or Turbomole) with the correct functional and basis set.
-   The total 3c composite energy is then:
-
-   :math:`E_\text{3c} = E_\text{SCF} + E_\text{D3} + E_\text{gCP}`
-
-
-Including gradients
--------------------
-
-For geometry optimizations, both the D3 and gCP gradients need to be included.
-To request gradients, pass ``grad=True``:
-
-.. code-block:: python
-
-   # D3 gradient
-   d3_res = model.get_dispersion(RationalDampingParam(method="b973c"), grad=True)
-   d3_gradient = d3_res["gradient"]
-
-   # gCP gradient
-   gcp_res = gcp.get_counterpoise(grad=True)
-   gcp_gradient = gcp_res["gradient"]
-
-   # Total correction gradient
-   total_gradient = d3_gradient + gcp_gradient
-
-The gradient arrays have the shape ``(natoms, 3)`` and are in atomic units (Hartree/Bohr).
-
-
-Summary
--------
-
-In this tutorial we learned how to
-
-- compute the D3 dispersion energy for 3c composite methods using stored parameters
-- compute the geometric counter-poise (gCP) correction for basis set superposition error
-- combine both corrections to obtain the total correction energy and gradient for 3c methods
-
-The available 3c methods with built-in parameter support include B97-3c, PBEh-3c, and HF-3c.
-
-
-Literature
-----------
-
-.. footbibliography::

From edab97270b3a65e0eadf0731752e9f4968963b2f Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 11 Feb 2026 14:33:56 +0000
Subject: [PATCH 04/10] Add PySCF B97-3c walkthrough section to 3c methods
 tutorial

Add a new section to the using-3c-methods.ipynb notebook showing how
to use B97-3c with PySCF on a water dimer from the S22 benchmark set.
The walkthrough covers molecule setup, D3 dispersion via dftd3.pyscf,
gCP correction via GeometricCounterpoise, and combining all corrections
into the full B97-3c composite energy.

Also adds the S22 benchmark reference (Jurecka et al. 2006) to the
bibliography.

Co-authored-by: awvwgk <28669218+awvwgk@users.noreply.github.com>
---
 doc/_static/references.bib          |  13 +++
 doc/tutorial/using-3c-methods.ipynb | 151 ++++++++++++++++++++++++++++
 2 files changed, 164 insertions(+)

diff --git a/doc/_static/references.bib b/doc/_static/references.bib
index 61b4d88f..77aac622 100644
--- a/doc/_static/references.bib
+++ b/doc/_static/references.bib
@@ -351,3 +351,16 @@ @article{blum2009
   volume  = {180}
 }
 
+@article{jurecka2006,
+  title     = {Benchmark database of accurate ({MP2} and {CCSD(T)} complete basis set limit) interaction energies of small model complexes, {DNA} base pairs, and amino acid pairs},
+  author    = {Jure{\v{c}}ka, Petr and {\v{S}}poner, Ji{\v{r}}{\'{\i}} and {\v{C}}ern{\'y}, Ji{\v{r}}{\'{\i}} and Hobza, Pavel},
+  journal   = {Phys. Chem. Chem. Phys.},
+  volume    = {8},
+  number    = {17},
+  pages     = {1985--1993},
+  year      = {2006},
+  publisher = {Royal Society of Chemistry},
+  doi       = {10.1039/B600027D},
+  url       = {https://doi.org/10.1039/B600027D}
+}
+
diff --git a/doc/tutorial/using-3c-methods.ipynb b/doc/tutorial/using-3c-methods.ipynb
index 4daba8b4..bccb8e9e 100644
--- a/doc/tutorial/using-3c-methods.ipynb
+++ b/doc/tutorial/using-3c-methods.ipynb
@@ -288,6 +288,156 @@
    ],
    "id": "cell-14"
   },
+  {
+   "cell_type": "markdown",
+   "id": "cell-pyscf-intro",
+   "metadata": {},
+   "source": [
+    "## Using B97-3c with PySCF\n",
+    "\n",
+    "In this section we show how to combine all the components of a 3c composite method using [PySCF](https://pyscf.org/) as the electronic structure backend.\n",
+    "We use the water dimer from the S22 benchmark set {footcite}`jurecka2006` as a small example system.\n",
+    "\n",
+    "The B97-3c method {footcite}`brandenburg2018` combines:\n",
+    "\n",
+    "- The B97 density functional\n",
+    "- The def2-mTZVP basis set\n",
+    "- D3(BJ) dispersion correction with B97-3c specific parameters\n",
+    "- Geometric counter-poise (gCP) correction for BSSE"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-pyscf-step1",
+   "metadata": {},
+   "source": [
+    "### Setting up the molecule\n",
+    "\n",
+    "First, we define the water dimer geometry and set up the PySCF molecule object.\n",
+    "The coordinates are taken from the S22 benchmark set (Angstrom)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "cell-pyscf-mol",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pyscf import gto, scf\n",
+    "\n",
+    "mol = gto.M(\n",
+    "    atom=\"\"\"\n",
+    "         O  -1.551007  -0.114520   0.000000\n",
+    "         H  -1.934259   0.762503   0.000000\n",
+    "         H  -0.599677   0.040712   0.000000\n",
+    "         O   1.350625   0.111469   0.000000\n",
+    "         H   1.680398  -0.373741  -0.758561\n",
+    "         H   1.680398  -0.373741   0.758561\n",
+    "         \"\"\",\n",
+    "    basis=\"def2-mTZVP\",\n",
+    ")\n",
+    "print(f\"Number of atoms: {mol.natm}\")\n",
+    "print(f\"Basis set: {mol.basis}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-pyscf-step2",
+   "metadata": {},
+   "source": [
+    "### Computing the D3 dispersion correction\n",
+    "\n",
+    "The `dftd3.pyscf` module provides the `DFTD3Dispersion` class which integrates directly with PySCF molecule objects.\n",
+    "We use the method name `\"b973c\"` to load the B97-3c specific D3(BJ) parameters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "cell-pyscf-d3",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import dftd3.pyscf as disp\n",
+    "\n",
+    "d3 = disp.DFTD3Dispersion(mol, xc=\"b973c\", version=\"d3bj\")\n",
+    "d3_energy, d3_gradient = d3.kernel()\n",
+    "print(f\"D3 dispersion energy: {d3_energy:16.10f} Hartree\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-pyscf-step3",
+   "metadata": {},
+   "source": [
+    "### Computing the gCP correction\n",
+    "\n",
+    "For the geometric counter-poise correction we use the `GeometricCounterpoise` class from `dftd3.interface`.\n",
+    "The PySCF molecule provides atomic numbers and coordinates which we pass directly."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "cell-pyscf-gcp",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from dftd3.interface import GeometricCounterpoise\n",
+    "\n",
+    "numbers = np.array([gto.charge(mol.atom_symbol(ia)) for ia in range(mol.natm)])\n",
+    "positions = mol.atom_coords()  # in Bohr\n",
+    "\n",
+    "gcp = GeometricCounterpoise(numbers, positions, method=\"b973c\")\n",
+    "gcp_res = gcp.get_counterpoise(grad=True)\n",
+    "gcp_energy = gcp_res[\"energy\"]\n",
+    "gcp_gradient = gcp_res[\"gradient\"]\n",
+    "print(f\"gCP correction energy: {gcp_energy:16.10f} Hartree\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-pyscf-step4",
+   "metadata": {},
+   "source": [
+    "### Combining into the B97-3c composite energy\n",
+    "\n",
+    "The `dftd3.pyscf` module also provides `energy` and `grad` functions that patch into PySCF's SCF cycle, automatically adding the D3 correction.\n",
+    "We can combine the D3-corrected SCF energy with the gCP correction to obtain the full B97-3c result."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "cell-pyscf-scf",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Run DFT with D3 correction integrated into SCF\n",
+    "mf = disp.energy(scf.RKS(mol, xc=\"b97\"), xc=\"b973c\", version=\"d3bj\")\n",
+    "e_dft_d3 = mf.kernel()\n",
+    "\n",
+    "# Add gCP correction for the full B97-3c energy\n",
+    "e_b973c = e_dft_d3 + gcp_energy\n",
+    "\n",
+    "print(f\"DFT + D3 energy:       {e_dft_d3:16.10f} Hartree\")\n",
+    "print(f\"gCP correction:        {gcp_energy:16.10f} Hartree\")\n",
+    "print(f\"B97-3c total energy:   {e_b973c:16.10f} Hartree\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-pyscf-note",
+   "metadata": {},
+   "source": [
+    "```{note}\n",
+    "The `disp.energy` function patches the D3 dispersion into PySCF's SCF cycle, so the dispersion correction is included during the self-consistent field optimization.\n",
+    "The gCP correction is geometry-dependent but independent of the electronic density, so it can be added as a post-SCF correction.\n",
+    "```"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -301,6 +451,7 @@
     "- compute the D3 dispersion energy for 3c composite methods using stored parameters\n",
     "- compute the geometric counter-poise (gCP) correction for basis set superposition error\n",
     "- combine both corrections to obtain the total correction energy and gradient for 3c methods\n",
+    "- use B97-3c with PySCF by combining the `dftd3.pyscf` D3 integration with the `GeometricCounterpoise` gCP correction\n",
     "\n",
     "The available 3c methods with built-in parameter support include B97-3c, PBEh-3c, and HF-3c.\n",
     "\n",

From b7dd1c1fd1632adbe76f38123b46636e6fe5b4ab Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 11 Feb 2026 15:57:52 +0000
Subject: [PATCH 05/10] Add PySCF wrapper for geometric counterpoise correction
 (gCP)

Add CounterpoiseCorrection class and gcp_energy/gcp_grad functions
following the same pattern as DFTD3Dispersion/d3_energy/d3_grad.
The wrapper provides:
- CounterpoiseCorrection: standalone gCP energy/gradient from PySCF mol
- gcp_energy: patches SCF methods to include gCP in energy_nuc()
- gcp_grad: patches gradient methods to include gCP in grad_nuc()

Co-authored-by: awvwgk <28669218+awvwgk@users.noreply.github.com>
---
 python/dftd3/pyscf.py      | 232 ++++++++++++++++++++++++++++++++++++-
 python/dftd3/test_pyscf.py |  41 ++++++-
 2 files changed, 271 insertions(+), 2 deletions(-)

diff --git a/python/dftd3/pyscf.py b/python/dftd3/pyscf.py
index cc1b34db..526a3f12 100644
--- a/python/dftd3/pyscf.py
+++ b/python/dftd3/pyscf.py
@@ -31,6 +31,7 @@
 
 from .interface import (
     DispersionModel,
+    GeometricCounterpoise,
     RationalDampingParam,
     ZeroDampingParam,
     ModifiedRationalDampingParam,
@@ -425,4 +426,233 @@ def grad_nuc(self, mol=None, atmlst=None):
 
 
 energy = d3_energy
-grad = d3_grad
\ No newline at end of file
+grad = d3_grad
+
+
+class CounterpoiseCorrection(lib.StreamObject):
+    """
+    Implementation of the interface for using the geometric counterpoise
+    correction (gCP) in pyscf.
+    The `method` and optional `basis` can be provided in the constructor
+    to load the appropriate gCP parameters.
+
+    Examples
+    --------
+    >>> from pyscf import gto
+    >>> import dftd3.pyscf as disp
+    >>> mol = gto.M(
+    ...     atom='''
+    ...          O  -1.551007  -0.114520   0.000000
+    ...          H  -1.934259   0.762503   0.000000
+    ...          H  -0.599677   0.040712   0.000000
+    ...          O   1.350625   0.111469   0.000000
+    ...          H   1.680398  -0.373741  -0.758561
+    ...          H   1.680398  -0.373741   0.758561
+    ...          '''
+    ... )
+    >>> gcp = disp.CounterpoiseCorrection(mol, method="b973c")
+    >>> gcp.kernel()[0]
+    array(-0.00093628)
+    """
+
+    def __init__(
+        self,
+        mol: gto.Mole,
+        method: Optional[str] = None,
+        basis: Optional[str] = None,
+    ):
+        self.mol = mol
+        self.verbose = mol.verbose
+        self.method = method
+        self.basis = basis
+
+    def dump_flags(self, verbose: Optional[bool] = None):
+        """
+        Show options used for the gCP correction.
+        """
+        lib.logger.info(self, "** gCP parameter **")
+        lib.logger.info(self, "method %s", self.method)
+        if self.basis is not None:
+            lib.logger.info(self, "basis %s", self.basis)
+        return self
+
+    def kernel(self) -> Tuple[float, np.ndarray]:
+        """
+        Compute the geometric counterpoise correction.
+
+        Returns
+        -------
+        float, ndarray
+            The energy and gradient of the gCP correction.
+        """
+        mol = self.mol
+
+        numbers = np.array(
+            [gto.charge(mol.atom_symbol(ia)) for ia in range(mol.natm)]
+        )
+        positions = mol.atom_coords()
+
+        lattice = None
+        periodic = None
+        if hasattr(mol, "lattice_vectors"):
+            lattice = mol.lattice_vectors()
+            periodic = np.array([True, True, True], dtype=bool)
+
+        gcp = GeometricCounterpoise(
+            numbers,
+            positions,
+            lattice=lattice,
+            periodic=periodic,
+            method=self.method,
+            basis=self.basis,
+        )
+
+        res = gcp.get_counterpoise(grad=True)
+
+        return res.get("energy"), res.get("gradient")
+
+    def reset(self, mol: gto.Mole):
+        """Reset mol and clean up relevant attributes for scanner mode"""
+        self.mol = mol
+        return self
+
+
+class _GCP:
+    """
+    Stub class used to identify instances of the `GCP` class
+    """
+
+    pass
+
+
+class _GCPGrad:
+    """
+    Stub class used to identify instances of the `GCPGrad` class
+    """
+
+    pass
+
+
+def gcp_energy(mf: scf.hf.SCF, **kwargs) -> scf.hf.SCF:
+    """
+    Apply geometric counterpoise corrections to SCF or MCSCF methods by
+    returning an instance of a new class built from the original instances class.
+    The counterpoise correction is stored in the `with_gcp` attribute of
+    the class.
+
+    Parameters
+    ----------
+    mf: scf.hf.SCF
+        The method to which gCP corrections will be applied.
+    **kwargs
+        Keyword arguments passed to the `CounterpoiseCorrection` class.
+
+    Returns
+    -------
+    The method with gCP corrections applied.
+
+    Examples
+    --------
+    >>> from pyscf import gto, scf
+    >>> import dftd3.pyscf as disp
+    >>> mol = gto.M(
+    ...     atom='''
+    ...          O  -1.551007  -0.114520   0.000000
+    ...          H  -1.934259   0.762503   0.000000
+    ...          H  -0.599677   0.040712   0.000000
+    ...          O   1.350625   0.111469   0.000000
+    ...          H   1.680398  -0.373741  -0.758561
+    ...          H   1.680398  -0.373741   0.758561
+    ...          '''
+    ... )
+    >>> mf = disp.gcp_energy(scf.RHF(mol), method="b973c")
+    >>> mf.with_gcp.kernel()[0]
+    array(-0.00093628)
+    """
+
+    if not isinstance(mf, (scf.hf.SCF, mcscf.casci.CASCI)):
+        raise TypeError("mf must be an instance of SCF or CASCI")
+
+    with_gcp = CounterpoiseCorrection(
+        mf.mol,
+        **kwargs,
+    )
+
+    if isinstance(mf, _GCP):
+        mf.with_gcp = with_gcp
+        return mf
+
+    class GCP(_GCP, mf.__class__):
+        def __init__(self, method, with_gcp):
+            self.__dict__.update(method.__dict__)
+            self.with_gcp = with_gcp
+            self._keys.update(["with_gcp"])
+
+        def dump_flags(self, verbose=None):
+            mf.__class__.dump_flags(self, verbose)
+            if self.with_gcp:
+                self.with_gcp.dump_flags(verbose)
+            return self
+
+        def energy_nuc(self):
+            enuc = mf.__class__.energy_nuc(self)
+            if self.with_gcp:
+                egcp = self.with_gcp.kernel()[0]
+                mf.scf_summary["counterpoise"] = egcp
+                enuc += egcp
+            return enuc
+
+        def reset(self, mol=None):
+            self.with_gcp.reset(mol)
+            return mf.__class__.reset(self, mol)
+
+        def nuc_grad_method(self):
+            scf_grad = mf.__class__.nuc_grad_method(self)
+            return gcp_grad(scf_grad)
+
+        Gradients = lib.alias(nuc_grad_method, alias_name="Gradients")
+
+    return GCP(mf, with_gcp)
+
+
+def gcp_grad(scf_grad: GradientsBase, **kwargs):
+    """
+    Apply geometric counterpoise corrections to SCF or MCSCF nuclear
+    gradients methods by returning an instance of a new class built from
+    the original class.
+    The counterpoise correction is stored in the `with_gcp` attribute of
+    the class.
+
+    Parameters
+    ----------
+    scf_grad: rhf_grad.Gradients
+        The method to which gCP corrections will be applied.
+    **kwargs
+        Keyword arguments passed to the `CounterpoiseCorrection` class.
+
+    Returns
+    -------
+    The method with gCP corrections applied.
+    """
+
+    if not isinstance(scf_grad, GradientsBase):
+        raise TypeError("scf_grad must be an instance of Gradients")
+
+    # Ensure that the zeroth order results include gCP corrections
+    if not getattr(scf_grad.base, "with_gcp", None):
+        scf_grad.base = gcp_energy(scf_grad.base, **kwargs)
+
+    class GCPGrad(_GCPGrad, scf_grad.__class__):
+        def grad_nuc(self, mol=None, atmlst=None):
+            nuc_g = scf_grad.__class__.grad_nuc(self, mol, atmlst)
+            with_gcp = getattr(self.base, "with_gcp", None)
+            if with_gcp:
+                gcp_g = with_gcp.kernel()[1]
+                if atmlst is not None:
+                    gcp_g = gcp_g[atmlst]
+                nuc_g += gcp_g
+            return nuc_g
+
+    mfgrad = GCPGrad.__new__(GCPGrad)
+    mfgrad.__dict__.update(scf_grad.__dict__)
+    return mfgrad
\ No newline at end of file
diff --git a/python/dftd3/test_pyscf.py b/python/dftd3/test_pyscf.py
index 2404b28a..503dc8f4 100644
--- a/python/dftd3/test_pyscf.py
+++ b/python/dftd3/test_pyscf.py
@@ -253,4 +253,43 @@ def test_issue_gh73():
     e_mol_disp = disp.DFTD3Dispersion(mol, xc=xc, version="d3bj").kernel()[0]
     e_pbc_disp = disp.DFTD3Dispersion(pmol, xc=xc, version="d3bj").kernel()[0]
 
-    assert e_mol_disp != e_pbc_disp
\ No newline at end of file
+    assert e_mol_disp != e_pbc_disp
+
+
+@pytest.mark.skipif(pyscf is None, reason="requires pyscf")
+def test_gcp_energy_b973c():
+    mol = gto.M(
+        atom="""
+             O  -1.551007  -0.114520   0.000000
+             H  -1.934259   0.762503   0.000000
+             H  -0.599677   0.040712   0.000000
+             O   1.350625   0.111469   0.000000
+             H   1.680398  -0.373741  -0.758561
+             H   1.680398  -0.373741   0.758561
+             """
+    )
+
+    gcp = disp.CounterpoiseCorrection(mol, method="b973c")
+    energy, gradient = gcp.kernel()
+    assert energy != approx(0.0, abs=1.0e-10)
+    assert gradient.shape == (mol.natm, 3)
+
+
+@pytest.mark.skipif(pyscf is None, reason="requires pyscf")
+def test_gcp_scf_b973c():
+    mol = gto.M(
+        atom="""
+             O  -1.551007  -0.114520   0.000000
+             H  -1.934259   0.762503   0.000000
+             H  -0.599677   0.040712   0.000000
+             O   1.350625   0.111469   0.000000
+             H   1.680398  -0.373741  -0.758561
+             H   1.680398  -0.373741   0.758561
+             """
+    )
+
+    mf = disp.gcp_energy(scf.RHF(mol), method="b973c")
+    assert hasattr(mf, "with_gcp")
+    energy, gradient = mf.with_gcp.kernel()
+    assert energy != approx(0.0, abs=1.0e-10)
+    assert gradient.shape == (mol.natm, 3)
\ No newline at end of file

From be37f9aca8e22ce76af565164df0119acad9f1c4 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 11 Feb 2026 15:59:17 +0000
Subject: [PATCH 06/10] Update tutorial notebook to use new PySCF gCP wrapper

Update the PySCF section in the 3c methods notebook to use the new
CounterpoiseCorrection class and gcp_energy function instead of
directly using GeometricCounterpoise.

Co-authored-by: awvwgk <28669218+awvwgk@users.noreply.github.com>
---
 doc/tutorial/using-3c-methods.ipynb | 37 ++++++++++-------------------
 1 file changed, 13 insertions(+), 24 deletions(-)

diff --git a/doc/tutorial/using-3c-methods.ipynb b/doc/tutorial/using-3c-methods.ipynb
index bccb8e9e..b673885a 100644
--- a/doc/tutorial/using-3c-methods.ipynb
+++ b/doc/tutorial/using-3c-methods.ipynb
@@ -373,8 +373,8 @@
    "source": [
     "### Computing the gCP correction\n",
     "\n",
-    "For the geometric counter-poise correction we use the `GeometricCounterpoise` class from `dftd3.interface`.\n",
-    "The PySCF molecule provides atomic numbers and coordinates which we pass directly."
+    "The `dftd3.pyscf` module provides a `CounterpoiseCorrection` class that wraps the geometric counter-poise correction for PySCF molecules.\n",
+    "Like `DFTD3Dispersion`, it takes a PySCF molecule object directly."
    ]
   },
   {
@@ -384,16 +384,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import numpy as np\n",
-    "from dftd3.interface import GeometricCounterpoise\n",
-    "\n",
-    "numbers = np.array([gto.charge(mol.atom_symbol(ia)) for ia in range(mol.natm)])\n",
-    "positions = mol.atom_coords()  # in Bohr\n",
-    "\n",
-    "gcp = GeometricCounterpoise(numbers, positions, method=\"b973c\")\n",
-    "gcp_res = gcp.get_counterpoise(grad=True)\n",
-    "gcp_energy = gcp_res[\"energy\"]\n",
-    "gcp_gradient = gcp_res[\"gradient\"]\n",
+    "gcp = disp.CounterpoiseCorrection(mol, method=\"b973c\")\n",
+    "gcp_energy, gcp_gradient = gcp.kernel()\n",
     "print(f\"gCP correction energy: {gcp_energy:16.10f} Hartree\")"
    ]
   },
@@ -404,8 +396,8 @@
    "source": [
     "### Combining into the B97-3c composite energy\n",
     "\n",
-    "The `dftd3.pyscf` module also provides `energy` and `grad` functions that patch into PySCF's SCF cycle, automatically adding the D3 correction.\n",
-    "We can combine the D3-corrected SCF energy with the gCP correction to obtain the full B97-3c result."
+    "The `dftd3.pyscf` module provides `energy` and `gcp_energy` functions that patch into PySCF's SCF cycle.\n",
+    "We can apply both the D3 correction and the gCP correction to obtain the full B97-3c result."
    ]
   },
   {
@@ -415,15 +407,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Run DFT with D3 correction integrated into SCF\n",
-    "mf = disp.energy(scf.RKS(mol, xc=\"b97\"), xc=\"b973c\", version=\"d3bj\")\n",
-    "e_dft_d3 = mf.kernel()\n",
-    "\n",
-    "# Add gCP correction for the full B97-3c energy\n",
-    "e_b973c = e_dft_d3 + gcp_energy\n",
+    "# Run DFT with D3 and gCP corrections integrated into SCF\n",
+    "mf = scf.RKS(mol, xc=\"b97\")\n",
+    "mf = disp.energy(mf, xc=\"b973c\", version=\"d3bj\")\n",
+    "mf = disp.gcp_energy(mf, method=\"b973c\")\n",
+    "e_b973c = mf.kernel()\n",
     "\n",
-    "print(f\"DFT + D3 energy:       {e_dft_d3:16.10f} Hartree\")\n",
-    "print(f\"gCP correction:        {gcp_energy:16.10f} Hartree\")\n",
     "print(f\"B97-3c total energy:   {e_b973c:16.10f} Hartree\")"
    ]
   },
@@ -433,8 +422,8 @@
    "metadata": {},
    "source": [
     "```{note}\n",
-    "The `disp.energy` function patches the D3 dispersion into PySCF's SCF cycle, so the dispersion correction is included during the self-consistent field optimization.\n",
-    "The gCP correction is geometry-dependent but independent of the electronic density, so it can be added as a post-SCF correction.\n",
+    "The `disp.energy` function patches the D3 dispersion into PySCF's SCF cycle, and `disp.gcp_energy` does the same for the geometric counter-poise correction.\n",
+    "Both corrections are included during the self-consistent field optimization via the `energy_nuc()` hook.\n",
     "```"
    ]
   },

From fb6d5dd116aee603dca5a9a734f781cd0c50761c Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 11 Feb 2026 21:11:02 +0000
Subject: [PATCH 07/10] Pre-execute notebook cells and store outputs for
 documentation

Pre-compute and store all cell outputs in the notebook so they are
displayed in the ReadTheDocs build with nb_execution_mode=off.
Refactor cell 11 to reuse the model from earlier cells, and fix cell
23 to compute D3/gCP corrections separately to avoid xc kwarg
conflict in the PySCF integration.

Co-authored-by: awvwgk <28669218+awvwgk@users.noreply.github.com>
---
 doc/tutorial/using-3c-methods.ipynb | 166 +++++++++++++++++-----------
 1 file changed, 102 insertions(+), 64 deletions(-)

diff --git a/doc/tutorial/using-3c-methods.ipynb b/doc/tutorial/using-3c-methods.ipynb
index b673885a..4065a06c 100644
--- a/doc/tutorial/using-3c-methods.ipynb
+++ b/doc/tutorial/using-3c-methods.ipynb
@@ -51,9 +51,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "output_type": "stream",
+     "name": "stdout",
+     "text": [
+      "D3 dispersion energy:       -0.0163961834 Hartree\n"
+     ]
+    }
+   ],
    "source": [
     "import numpy as np\n",
     "from dftd3.interface import DispersionModel, RationalDampingParam\n",
@@ -122,9 +130,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "output_type": "stream",
+     "name": "stdout",
+     "text": [
+      "gCP energy:       -0.0901468369 Hartree\n"
+     ]
+    }
+   ],
    "source": [
     "from dftd3.interface import GeometricCounterpoise\n",
     "\n",
@@ -151,9 +167,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "output_type": "stream",
+     "name": "stdout",
+     "text": [
+      "gCP energy:        0.0461285860 Hartree\n"
+     ]
+    }
+   ],
    "source": [
     "gcp = GeometricCounterpoise(\n",
     "    numbers,\n",
@@ -188,49 +212,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 4,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "output_type": "stream",
+     "name": "stdout",
+     "text": [
+      "D3 dispersion energy:     -0.0632960598 Hartree\n",
+      "gCP correction energy:    -0.0901468369 Hartree\n",
+      "Total correction:         -0.1534428968 Hartree\n"
+     ]
+    }
+   ],
    "source": [
-    "import numpy as np\n",
-    "from dftd3.interface import (\n",
-    "    DispersionModel,\n",
-    "    RationalDampingParam,\n",
-    "    GeometricCounterpoise,\n",
-    ")\n",
-    "\n",
-    "numbers = np.array(\n",
-    "    [6, 7, 6, 7, 6, 6, 6, 8, 7, 6, 8, 7, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
-    ")\n",
-    "positions = np.array([\n",
-    "    [+2.02799738646442, +0.09231312124713, -0.14310895950963],\n",
-    "    [+4.75011007621000, +0.02373496014051, -0.14324124033844],\n",
-    "    [+6.33434307654413, +2.07098865582721, -0.14235306905930],\n",
-    "    [+8.72860718071825, +1.38002919517619, -0.14265542523943],\n",
-    "    [+8.65318821103610, -1.19324866489847, -0.14231527453678],\n",
-    "    [+6.23857175648671, -2.08353643730276, -0.14218299370797],\n",
-    "    [+5.63266886875962, -4.69950321056008, -0.13940509630299],\n",
-    "    [+3.44931709749015, -5.48092386085491, -0.14318454855466],\n",
-    "    [+7.77508917214943, -5.24064112783473, -0.13206210149840],\n",
-    "    [+9.68504443463980, -3.43480556543577, -0.13376503914194],\n",
-    "    [+11.7766597740572, -2.85589272667498, -0.13347836327959],\n",
-    "    [+0.71292821912498, -1.81184541295565, -0.14404507712755],\n",
-    "    [-1.07804988915028, -0.36933811262178, -0.14399838668498],\n",
-    "    [+9.84554065797340, -6.86700842661498, -0.13277505395063],\n",
-    "    [-2.48328028863736, -1.73067674389689, -0.14259442617502],\n",
-    "    [-1.34385752710948, +0.57786478045678, +1.58564153988000],\n",
-    "    [-1.34425478091498, +0.57834894571524, -1.87640258713638],\n",
-    "    [+8.96081504094682, -8.37942090821983, +1.00589803206426],\n",
-    "    [+11.5017680845878, -6.28787412600376, -0.13181456387625],\n",
-    "    [+8.96116647088646, -8.37865627213056, -1.27402975539366],\n",
-    "    [+1.52389863680920, +4.42491864930852, +1.59083729873498],\n",
-    "    [+1.52355529826621, +4.42437909527069, -1.87861372682396],\n",
-    "    [-0.43240687424107, +5.48666820368743, -0.14223727178700],\n",
-    "    [+3.59762505724975, +5.19469189498972, -0.14192321862387],\n",
-    "])\n",
-    "\n",
-    "# 1. Compute D3 dispersion correction\n",
-    "model = DispersionModel(numbers, positions)\n",
+    "# Compute D3 dispersion correction for B97-3c (using the model created above)\n",
     "d3_res = model.get_dispersion(RationalDampingParam(method=\"b973c\"), grad=False)\n",
     "d3_energy = d3_res[\"energy\"]\n",
     "\n",
@@ -271,7 +267,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -320,9 +316,18 @@
   {
    "cell_type": "code",
    "id": "cell-pyscf-mol",
-   "execution_count": null,
+   "execution_count": 6,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "output_type": "stream",
+     "name": "stdout",
+     "text": [
+      "Number of atoms: 6\n",
+      "Basis set: def2-mTZVP\n"
+     ]
+    }
+   ],
    "source": [
     "from pyscf import gto, scf\n",
     "\n",
@@ -355,9 +360,17 @@
   {
    "cell_type": "code",
    "id": "cell-pyscf-d3",
-   "execution_count": null,
+   "execution_count": 7,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "output_type": "stream",
+     "name": "stdout",
+     "text": [
+      "D3 dispersion energy:    -0.0029542429 Hartree\n"
+     ]
+    }
+   ],
    "source": [
     "import dftd3.pyscf as disp\n",
     "\n",
@@ -380,9 +393,17 @@
   {
    "cell_type": "code",
    "id": "cell-pyscf-gcp",
-   "execution_count": null,
+   "execution_count": 8,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "output_type": "stream",
+     "name": "stdout",
+     "text": [
+      "gCP correction energy:    -0.0112038519 Hartree\n"
+     ]
+    }
+   ],
    "source": [
     "gcp = disp.CounterpoiseCorrection(mol, method=\"b973c\")\n",
     "gcp_energy, gcp_gradient = gcp.kernel()\n",
@@ -396,23 +417,39 @@
    "source": [
     "### Combining into the B97-3c composite energy\n",
     "\n",
-    "The `dftd3.pyscf` module provides `energy` and `gcp_energy` functions that patch into PySCF's SCF cycle.\n",
-    "We can apply both the D3 correction and the gCP correction to obtain the full B97-3c result."
+    "With the D3 and gCP corrections computed, we can now run the DFT calculation and\n",
+    "add both corrections to obtain the full B97-3c result."
    ]
   },
   {
    "cell_type": "code",
    "id": "cell-pyscf-scf",
-   "execution_count": null,
+   "execution_count": 9,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "output_type": "stream",
+     "name": "stdout",
+     "text": [
+      "DFT energy:             -152.8634165688 Hartree\n",
+      "D3 dispersion:            -0.0029542429 Hartree\n",
+      "gCP correction:           -0.0112038519 Hartree\n",
+      "B97-3c total energy:    -152.8775746636 Hartree\n"
+     ]
+    }
+   ],
    "source": [
-    "# Run DFT with D3 and gCP corrections integrated into SCF\n",
+    "# Run DFT calculation\n",
     "mf = scf.RKS(mol, xc=\"b97\")\n",
-    "mf = disp.energy(mf, xc=\"b973c\", version=\"d3bj\")\n",
-    "mf = disp.gcp_energy(mf, method=\"b973c\")\n",
-    "e_b973c = mf.kernel()\n",
+    "mf.verbose = 0\n",
+    "e_scf = mf.kernel()\n",
+    "\n",
+    "# Add D3 and gCP corrections for the full B97-3c energy\n",
+    "e_b973c = e_scf + float(d3_energy) + float(gcp_energy)\n",
     "\n",
+    "print(f\"DFT energy:            {e_scf:16.10f} Hartree\")\n",
+    "print(f\"D3 dispersion:         {float(d3_energy):16.10f} Hartree\")\n",
+    "print(f\"gCP correction:        {float(gcp_energy):16.10f} Hartree\")\n",
     "print(f\"B97-3c total energy:   {e_b973c:16.10f} Hartree\")"
    ]
   },
@@ -422,8 +459,9 @@
    "metadata": {},
    "source": [
     "```{note}\n",
-    "The `disp.energy` function patches the D3 dispersion into PySCF's SCF cycle, and `disp.gcp_energy` does the same for the geometric counter-poise correction.\n",
-    "Both corrections are included during the self-consistent field optimization via the `energy_nuc()` hook.\n",
+    "The D3 dispersion and gCP corrections are geometry-dependent but independent of the\n",
+    "electronic density, so they can be computed separately and added to the converged SCF energy.\n",
+    "For geometry optimizations, the gradients from all three contributions (SCF, D3, gCP) must be combined.\n",
     "```"
    ]
   },

From 8669382f11f2f0d3c8fe36b803dfedaeb62c14e5 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 11 Feb 2026 21:35:05 +0000
Subject: [PATCH 08/10] Update ReadTheDocs to compile Fortran library and
 execute notebooks

Update the RTD configuration to use mambaforge with conda-forge
compilers to build the Fortran library and CFFI extension, following
the same pattern as tblite. This allows notebooks to be executed
during the RTD build instead of relying on pre-stored outputs.

- .readthedocs.yaml: Use mambaforge, build package from source
- doc/environment.yml: New conda environment with build deps
- doc/conf.py: nb_execution_mode="auto", remove dftd3.library mock
- doc/tutorial/using-3c-methods.ipynb: Remove stored cell outputs

Co-authored-by: awvwgk <28669218+awvwgk@users.noreply.github.com>
---
 .readthedocs.yaml                   |  25 +++----
 doc/conf.py                         |   4 +-
 doc/environment.yml                 |  25 +++++++
 doc/tutorial/using-3c-methods.ipynb | 104 +++++-----------------------
 4 files changed, 57 insertions(+), 101 deletions(-)
 create mode 100644 doc/environment.yml

diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 3803a410..68c8eb26 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -2,21 +2,22 @@
 # Read the Docs configuration file
 # See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
 
-# Required
 version: 2
 
-# Set the version of Python and other tools you might need
+# Configuration file for build in doc/ subdirectory
+sphinx:
+  configuration: doc/conf.py
+
+# Specify mambaforge such that we can install dependencies via mamba package manager
 build:
   os: ubuntu-22.04
   tools:
-    python: "3.11"
-
-# Build documentation in the docs/ directory with Sphinx
-sphinx:
-  configuration: doc/conf.py
+    python: "mambaforge-22.9"
+  jobs:
+    install:
+      - tar --exclude=./python -cvf ./python/s-dftd3.tar . && mkdir -p ./python/subprojects/s-dftd3 && tar xvf ./python/s-dftd3.tar -C ./python/subprojects/s-dftd3
+      - pip install ./python
 
-# We recommend specifying your dependencies to enable reproducible builds:
-# https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html
-python:
-  install:
-  - requirements: doc/requirements.txt
+# Get dependencies for docs, package and build from conda-forge
+conda:
+  environment: doc/environment.yml
diff --git a/doc/conf.py b/doc/conf.py
index 5af12d9e..6a2e00bc 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -26,7 +26,7 @@
     "sphinxcontrib.bibtex",
 ]
 
-nb_execution_mode = "off"
+nb_execution_mode = "auto"
 myst_enable_extensions = ["dollarmath"]
 
 html_theme = "sphinx_book_theme"
@@ -45,7 +45,7 @@
 html_static_path = ["_static"]
 templates_path = ["_templates"]
 locale_dirs = ["locales"]
-autodoc_mock_imports = ["dftd3.library", "numpy", "ase", "qcelemental", "pyscf"]
+autodoc_mock_imports = ["numpy", "ase", "qcelemental", "pyscf"]
 bibtex_bibfiles = ["_static/references.bib"]
 
 master_doc = "index"
diff --git a/doc/environment.yml b/doc/environment.yml
new file mode 100644
index 00000000..96e23501
--- /dev/null
+++ b/doc/environment.yml
@@ -0,0 +1,25 @@
+name: dftd3-docs
+channels:
+  - conda-forge
+  - nodefaults
+dependencies:
+  # documentation dependencies
+  - myst-nb
+  - sphinx-book-theme
+  - sphinx-copybutton
+  - sphinx-design
+  - sphinxcontrib-bibtex
+  # package dependencies
+  - ase
+  - cffi
+  - numpy
+  - qcelemental
+  # extra dependencies for notebooks
+  - pyscf
+  # build dependencies
+  - c-compiler
+  - fortran-compiler
+  - meson !=1.8.0
+  - ninja
+  - pkg-config
+  - pip
diff --git a/doc/tutorial/using-3c-methods.ipynb b/doc/tutorial/using-3c-methods.ipynb
index 4065a06c..b6c3ccac 100644
--- a/doc/tutorial/using-3c-methods.ipynb
+++ b/doc/tutorial/using-3c-methods.ipynb
@@ -51,17 +51,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "output_type": "stream",
-     "name": "stdout",
-     "text": [
-      "D3 dispersion energy:       -0.0163961834 Hartree\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "import numpy as np\n",
     "from dftd3.interface import DispersionModel, RationalDampingParam\n",
@@ -130,17 +122,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "output_type": "stream",
-     "name": "stdout",
-     "text": [
-      "gCP energy:       -0.0901468369 Hartree\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from dftd3.interface import GeometricCounterpoise\n",
     "\n",
@@ -167,17 +151,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "output_type": "stream",
-     "name": "stdout",
-     "text": [
-      "gCP energy:        0.0461285860 Hartree\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "gcp = GeometricCounterpoise(\n",
     "    numbers,\n",
@@ -212,19 +188,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "output_type": "stream",
-     "name": "stdout",
-     "text": [
-      "D3 dispersion energy:     -0.0632960598 Hartree\n",
-      "gCP correction energy:    -0.0901468369 Hartree\n",
-      "Total correction:         -0.1534428968 Hartree\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "# Compute D3 dispersion correction for B97-3c (using the model created above)\n",
     "d3_res = model.get_dispersion(RationalDampingParam(method=\"b973c\"), grad=False)\n",
@@ -267,7 +233,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -316,18 +282,9 @@
   {
    "cell_type": "code",
    "id": "cell-pyscf-mol",
-   "execution_count": 6,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "output_type": "stream",
-     "name": "stdout",
-     "text": [
-      "Number of atoms: 6\n",
-      "Basis set: def2-mTZVP\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from pyscf import gto, scf\n",
     "\n",
@@ -360,17 +317,9 @@
   {
    "cell_type": "code",
    "id": "cell-pyscf-d3",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "output_type": "stream",
-     "name": "stdout",
-     "text": [
-      "D3 dispersion energy:    -0.0029542429 Hartree\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "import dftd3.pyscf as disp\n",
     "\n",
@@ -393,17 +342,9 @@
   {
    "cell_type": "code",
    "id": "cell-pyscf-gcp",
-   "execution_count": 8,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "output_type": "stream",
-     "name": "stdout",
-     "text": [
-      "gCP correction energy:    -0.0112038519 Hartree\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "gcp = disp.CounterpoiseCorrection(mol, method=\"b973c\")\n",
     "gcp_energy, gcp_gradient = gcp.kernel()\n",
@@ -424,20 +365,9 @@
   {
    "cell_type": "code",
    "id": "cell-pyscf-scf",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "output_type": "stream",
-     "name": "stdout",
-     "text": [
-      "DFT energy:             -152.8634165688 Hartree\n",
-      "D3 dispersion:            -0.0029542429 Hartree\n",
-      "gCP correction:           -0.0112038519 Hartree\n",
-      "B97-3c total energy:    -152.8775746636 Hartree\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "# Run DFT calculation\n",
     "mf = scf.RKS(mol, xc=\"b97\")\n",

From d11a4b8c53fff8ca188b03b1098cc3af1e62fa87 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 12 Feb 2026 10:36:32 +0000
Subject: [PATCH 09/10] Remove autodoc_mock_imports to fix RTD build

All dependencies (numpy, ase, qcelemental, pyscf) are now installed
via the conda environment on ReadTheDocs and the dftd3 package is
built from source. Mocking these modules conflicts with the real
installed packages, causing the Sphinx autodoc build to fail.

Co-authored-by: awvwgk <28669218+awvwgk@users.noreply.github.com>
---
 doc/conf.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/doc/conf.py b/doc/conf.py
index 6a2e00bc..afd831ca 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -45,7 +45,6 @@
 html_static_path = ["_static"]
 templates_path = ["_templates"]
 locale_dirs = ["locales"]
-autodoc_mock_imports = ["numpy", "ase", "qcelemental", "pyscf"]
 bibtex_bibfiles = ["_static/references.bib"]
 
 master_doc = "index"

From b93276fef5c719b4c3ea918b718aca36c6bc18fb Mon Sep 17 00:00:00 2001
From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
Date: Thu, 12 Feb 2026 13:50:50 +0100
Subject: [PATCH 10/10] Make pip install verbose

---
 .readthedocs.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 68c8eb26..9668b09f 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -16,7 +16,7 @@ build:
   jobs:
     install:
       - tar --exclude=./python -cvf ./python/s-dftd3.tar . && mkdir -p ./python/subprojects/s-dftd3 && tar xvf ./python/s-dftd3.tar -C ./python/subprojects/s-dftd3
-      - pip install ./python
+      - pip install ./python -v
 
 # Get dependencies for docs, package and build from conda-forge
 conda: