diff --git a/doc/source/derivs/derivs.ipynb b/doc/source/derivs/derivs.ipynb
index aa39bc31243defd2b4cff0dd1a92e88f426da8fa..408013f39e0e14a6c6e2c1484124166f435c5fe7 100644
--- a/doc/source/derivs/derivs.ipynb
+++ b/doc/source/derivs/derivs.ipynb
@@ -5,7 +5,26 @@
    "id": "1afc4517",
    "metadata": {},
    "source": [
-    "# Thermodynamic Derivatives"
+    "# Thermodynamic Derivatives\n",
+    "\n",
+    "## Helmholtz energy derivatives\n",
+    "Thermodynamic derivatives are at the very heart of teqp.  All models are defined in the form $\\alpha^r(T, \\rho, z)$, where $\\rho$ is the molar density, and z are mole fractions.  There are exceptions for models for which the independent variables are in simulation units (Lennard-Jones and its ilk).\n",
+    "\n",
+    "Thereofore, to obtain the residual pressure, it is obtained as a derivative:\n",
+    "\n",
+    "$$ p^r = \\rho R T \\left( \\rho \\left(\\frac{\\partial \\alpha^r}{\\partial \\rho}\\right)_{T}\\right)$$\n",
+    "\n",
+    "and other residual thermodynamic properties are defined likewise.\n",
+    "\n",
+    "We can define the concise derivative\n",
+    "\n",
+    "$$ \\Lambda^r_{xy} = (1/T)^x(\\rho)^y\\left(\\frac{\\partial^{x+y}(\\alpha^r)}{\\partial (1/T)^x\\partial \\rho^y}\\right)$$\n",
+    "\n",
+    "so we can re-write the derivative above as \n",
+    "\n",
+    "$$ p^r = \\rho R T \\Lambda^r_{01}$$\n",
+    "\n",
+    "Similar definitions apply for all the other residual thermodynamic properties"
    ]
   },
   {
@@ -74,6 +93,161 @@
     "z = np.array([1.0])\n",
     "model.get_Ar01(300,300,z)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "498a2350",
+   "metadata": {},
+   "source": [
+    "And there are additional methods to obtain all the derivatives up to a given order:\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "3ae755e1",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[-0.06966138343515413,\n",
+       " -0.06836660379313926,\n",
+       " 0.0025357822532378147,\n",
+       " -0.00015701162203571184,\n",
+       " 1.6818628788290574e-05,\n",
+       " -2.2305940927885907e-06,\n",
+       " 3.8259258513417917e-07]"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "model.get_Ar06n(300,300,z) # derivatives 00, 01, 02, ... 06"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cd644a3e",
+   "metadata": {},
+   "source": [
+    "But more derivatives are slower than fewer:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "4e0609ae",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1.93 µs ± 32.9 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n",
+      "2.13 µs ± 36.7 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit model.get_Ar01(300,300,z)\n",
+    "%timeit model.get_Ar04n(300,300,z)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d05e9070",
+   "metadata": {},
+   "source": [
+    "Note: calling overhead is usually on the order of 1 microsecond"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "502e9830",
+   "metadata": {},
+   "source": [
+    "## Virial coefficients\n",
+    "\n",
+    "Virial coefficients represent the thermodynamics of the interaction of two-, three-, ... bodies interacting with each other.  They can be obtained rigorously if the potential energy surface of interaction is fully known.  In general, such a surface can only be constructed for small rigid molecules.  Many simple thermodynamic models do a poor job of predicting the thermodynamics captured by the virial coefficients.  \n",
+    "\n",
+    "The i-th virial coefficient is defined by \n",
+    "$$\n",
+    "\tB_i = \\frac{(\\alpha^r)^{(i-1)}}{(i-2)!}\n",
+    "$$\n",
+    "with the concise derivative term\n",
+    "$$\n",
+    "\t(\\alpha^r)^{(i)} = \\lim_{\\rho \\to 0} \\left(\\frac{\\partial^{i}\\alpha^r}{\\partial \\rho^{i}}\\right)_{T,\\vec{x}}\n",
+    "$$\n",
+    "\n",
+    "teqp supports the virial coefficient directly, there is the ``get_B2vir`` for the second virial coefficient:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "720e120c",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "-0.00023661263734465424"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "model.get_B2vir(300, z)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2fba4369",
+   "metadata": {},
+   "source": [
+    "And the ``get_Bnvir`` method that allows for the calculation of higher virial coefficients:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "92d01992",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{2: -0.0002366126373446542,\n",
+       " 3: 3.001768410777936e-08,\n",
+       " 4: -3.2409760373816355e-12,\n",
+       " 5: 3.9617816466337214e-16,\n",
+       " 6: -4.552923983836698e-20,\n",
+       " 7: 5.3759278511184914e-24}"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "model.get_Bnvir(7, 300, z)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d50998b1",
+   "metadata": {},
+   "source": [
+    "The ``get_Bnvir`` function was implemented because when doing automatic differentiation, all the intermediate derivatives are also retained."
+   ]
   }
  ],
  "metadata": {