diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 3803a410..9668b09f 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 -v
 
-# 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/_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/conf.py b/doc/conf.py
index a1288084..24070b5a 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -1,7 +1,6 @@
 import os
 import sys
 
-sys.path.insert(0, os.path.join(os.path.abspath(".."), "python"))
 import dftd3
 
 
@@ -13,6 +12,7 @@
 release = version
 
 extensions = [
+    "myst_nb",
     "sphinx_design",
     "sphinx_copybutton",
     "sphinx.ext.autosummary",
@@ -25,6 +25,9 @@
     "sphinxcontrib.bibtex",
 ]
 
+nb_execution_mode = "auto"
+myst_enable_extensions = ["dollarmath"]
+
 html_theme = "sphinx_book_theme"
 html_title = "Simple DFT-D3"
 
@@ -41,7 +44,6 @@
 html_static_path = ["_static"]
 templates_path = ["_templates"]
 locale_dirs = ["locales"]
-autodoc_mock_imports = ["dftd3.library", "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/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/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.ipynb b/doc/tutorial/using-3c-methods.ipynb
new file mode 100644
index 00000000..b6c3ccac
--- /dev/null
+++ b/doc/tutorial/using-3c-methods.ipynb
@@ -0,0 +1,436 @@
+{
+ "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": [
+    "# 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",
+    "# 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",
+   "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",
+    "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "cell-pyscf-gcp",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gcp = disp.CounterpoiseCorrection(mol, method=\"b973c\")\n",
+    "gcp_energy, gcp_gradient = gcp.kernel()\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",
+    "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,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Run DFT calculation\n",
+    "mf = scf.RKS(mol, xc=\"b97\")\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\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-pyscf-note",
+   "metadata": {},
+   "source": [
+    "```{note}\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",
+    "```"
+   ]
+  },
+  {
+   "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",
+    "- 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",
+    "## 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/python/dftd3/pyscf.py b/python/dftd3/pyscf.py
index d1dc2ff7..e4a6d703 100644
--- a/python/dftd3/pyscf.py
+++ b/python/dftd3/pyscf.py
@@ -31,6 +31,7 @@
 
 from .interface import (
     DispersionModel,
+    GeometricCounterpoise,
     RationalDampingParam,
     ZeroDampingParam,
     ModifiedRationalDampingParam,
@@ -443,4 +444,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