diff --git a/examples/mcpdft/00-simple_mcpdft.py b/examples/mcpdft/00-simple_mcpdft.py
new file mode 100644
index 000000000..6fe893023
--- /dev/null
+++ b/examples/mcpdft/00-simple_mcpdft.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env/python
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M (
+    atom = 'O 0 0 0; O 0 0 1.2',
+    basis = 'ccpvdz',
+    spin = 2)
+
+mf = scf.RHF (mol).run ()
+
+# 1. CASCI density
+
+mc0 = mcpdft.CASCI (mf, 'tPBE', 6, 8).run ()
+
+# 2. CASSCF density
+# Note that the MC-PDFT energy may not be lower, even though
+# E(CASSCF)<=E(CASCI).
+
+mc1 = mcpdft.CASSCF (mf, 'tPBE', 6, 8).run ()
+
+# 3. analyze () does the same thing as CASSCF analyze ()
+
+mc1.verbose = 4
+mc1.analyze ()
+
+# 4. Energy decomposition for additional analysis
+
+e_decomp = mc1.get_energy_decomposition (split_x_c=False)
+print ("e_nuc =",e_decomp[0])
+print ("e_1e =",e_decomp[1])
+print ("e_Coul =",e_decomp[2])
+print ("e_OT =",e_decomp[3])
+print ("e_ncwfn (not included in total energy) =",e_decomp[4])
+print ("e_PDFT - e_MCSCF =", mc1.e_tot - mc1.e_mcscf)
+print ("e_OT - e_ncwfn =", e_decomp[3] - e_decomp[4])
+e_decomp = mc1.get_energy_decomposition (split_x_c=True)
+print ("e_OT (x component) = ",e_decomp[3])
+print ("e_OT (c component) = ",e_decomp[4])
+
+
diff --git a/examples/mcpdft/01-different_functionals.py b/examples/mcpdft/01-different_functionals.py
new file mode 100644
index 000000000..adbb1bda8
--- /dev/null
+++ b/examples/mcpdft/01-different_functionals.py
@@ -0,0 +1,49 @@
+#!/usr/bin/env/python
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M (
+    atom = 'O 0 0 0; O 0 0 1.2',
+    basis = 'ccpvdz',
+    spin = 2)
+
+mf = scf.RHF (mol).run ()
+
+mc = mcpdft.CASCI (mf, 'tPBE', 6, 8).run ()
+
+# 1. Change the functional by setting mc.otxc
+
+mc.otxc = 'ftBLYP'
+mc.kernel ()
+
+# 2. Change the functional and compute energy in one line without
+# reoptimizing the wave function using compute_pdft_energy_
+
+mc.compute_pdft_energy_(otxc='ftPBE')
+
+#
+# The leading "t" or "ft" identifies either a "translated functional"
+# [JCTC 10, 3669 (2014)] or a "fully-translated functional"
+# [JCTC 11, 4077 (2015)]. It can be combined with a general PySCF
+# xc string containing any number of pure LDA or GGA functionals.
+# Meta-GGAs, built-in hybrid functionals ("B3LYP"), and range-separated
+# functionals are not supported.
+# 
+
+# 3. A translated user-defined compound functional
+
+mc.compute_pdft_energy_(otxc="t0.3*B88 + 0.7*SLATER,0.4*VWN5+0.6*LYP")
+
+# 4. A fully-translated functional consisting of "exchange" only
+
+mc.compute_pdft_energy_(otxc="ftPBE,")
+
+# 5. A fully-translated functional consisting of "correlation" only
+
+mc.compute_pdft_energy_(otxc="ft,PBE")
+
+# 6. The sum of 5 and 6 (look at the "Eot" output)
+
+mc.compute_pdft_energy_(otxc="ftPBE")
+
+
diff --git a/examples/mcpdft/02-hybrid_functionals.py b/examples/mcpdft/02-hybrid_functionals.py
new file mode 100644
index 000000000..b8956afd6
--- /dev/null
+++ b/examples/mcpdft/02-hybrid_functionals.py
@@ -0,0 +1,49 @@
+#!/usr/bin/env/python
+
+
+#
+# JPCL 11, 10158 (2020)
+# (but see #5 below)
+#
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M (
+    atom = 'O 0 0 0; O 0 0 1.2',
+    basis = 'ccpvdz',
+    spin = 2)
+
+mf = scf.RHF (mol).run ()
+
+# 1. The only two predefined hybrid functionals as of writing
+
+mc = mcpdft.CASSCF (mf, 'tPBE0', 6, 8).run ()
+mc.compute_pdft_energy_(otxc='ftPBE0')
+
+# 2. Other predefined hybrid functionals are not supported
+
+try:
+    mc.compute_pdft_energy_(otxc='tB3LYP')
+except NotImplementedError as e:
+    print ("otxc='tB3LYP' results in NotImplementedError:")
+    print (str (e))
+
+# 3. Construct a custom hybrid functional by hand
+
+my_otxc = 't.8*B88 + .2*HF, .8*LYP + .2*HF'
+mc.compute_pdft_energy_(otxc=my_otxc)
+
+# 4. Construct the same custom functional using helper function
+
+my_otxc = 't' + mcpdft.hyb ('BLYP', .2)
+mc.compute_pdft_energy_(otxc=my_otxc)
+
+# 5. "lambda-MC-PDFT" of JCTC 16, 2274, 2020.
+
+my_otxc = 't' + mcpdft.hyb('PBE', .5, hyb_type='lambda')
+#       = 't0.5*PBE + 0.5*HF, 0.75*PBE + 0.5*HF'
+mc.compute_pdft_energy_(otxc=my_otxc)
+
+
+
+
diff --git a/examples/mcpdft/03-metaGGA_functionals.py b/examples/mcpdft/03-metaGGA_functionals.py
new file mode 100644
index 000000000..b61c0197e
--- /dev/null
+++ b/examples/mcpdft/03-metaGGA_functionals.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env/python
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M (
+    atom = 'O 0 0 0; O 0 0 1.2',
+    basis = 'ccpvdz',
+    spin = 2)
+
+mf = scf.RHF (mol).run ()
+
+# The translation of Meta-GGAs and hybrid-Meta-GGAs [PNAS, 122, 1, 2025, e2419413121; https://doi.org/10.1073/pnas.2419413121]
+
+# Translated-Meta-GGA
+mc = mcpdft.CASCI(mf, 'tM06L', 6, 8).run ()
+
+# Hybrid-Translated-Meta-GGA
+tM06L0 = 't' + mcpdft.hyb('M06L',0.25, hyb_type='average')
+mc = mcpdft.CASCI(mf, tM06L0, 6, 8).run ()
+
+# MC23: meta-hybrid on-top functional [PNAS, 122, 1, 2025, e2419413121; https://doi.org/10.1073/pnas.2419413121]
+
+# State-Specific
+mc = mcpdft.CASCI(mf, 'MC23', 6, 8)
+mc.kernel()
+
+# State-average
+nroots=2
+mc = mcpdft.CASCI(mf, 'MC23', 6, 8)
+mc.fcisolver.nroots=nroots
+mc.kernel()[0]
diff --git a/examples/mcpdft/11-grid_scheme.py b/examples/mcpdft/11-grid_scheme.py
new file mode 100644
index 000000000..0fc860cfa
--- /dev/null
+++ b/examples/mcpdft/11-grid_scheme.py
@@ -0,0 +1,40 @@
+#!/usr/bin/env python
+
+from pyscf import gto, scf, dft
+from pyscf import mcpdft
+
+# See also pyscf/examples/dft/11-grid-scheme.py
+
+mol = gto.M(
+    verbose = 0,
+    atom = '''
+    o    0    0.       0.
+    h    0    -0.757   0.587
+    h    0    0.757    0.587''',
+    basis = '6-31g')
+mf = scf.RHF (mol).run ()
+mc = mcpdft.CASSCF(mf, 'tLDA', 4, 4).run()
+print('Default grid setup.  E = %.12f' % mc.e_tot)
+
+# See pyscf/dft/radi.py for more radial grid schemes
+mc.grids.radi_method = dft.mura_knowles
+print('radi_method = mura_knowles.  E = %.12f' % mc.kernel()[0])
+
+# All in one command:
+mc = mcpdft.CASSCF (mf, 'tLDA', 4, 4,
+                    grids_attr={'radi_method': dft.gauss_chebyshev})
+print('radi_method = gauss_chebyshev.  E = %.12f' % mc.kernel ()[0])
+
+# Or inline with an already-built mc object:
+e = mc.compute_pdft_energy_(grids_attr={'radi_method': dft.delley})[0]
+print('radi_method = delley.  E = %.12f' % e)
+
+# All grids attributes can be addressed in any of the ways above
+# There is also a shortcut to address grids.level:
+
+mc = mcpdft.CASSCF(mf, 'tLDA', 4, 4, grids_level=4).run()
+print('grids.level = 4.  E = %.12f' % mc.e_tot)
+
+e = mc.compute_pdft_energy_(grids_level=2)[0]
+print('grids.level = 2.  E = %.12f' % e)
+
diff --git a/examples/mcpdft/15-state_average.py b/examples/mcpdft/15-state_average.py
new file mode 100644
index 000000000..886bcaa9d
--- /dev/null
+++ b/examples/mcpdft/15-state_average.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+
+'''
+State average
+
+This works the same way as mcscf.state_average_, but you
+must use the method attribute (mc.state_average, mc.state_average_)
+instead of the function call.
+'''
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M(
+    atom = [
+        ['O', ( 0., 0.    , 0.   )],
+        ['H', ( 0., -0.757, 0.587)],
+        ['H', ( 0., 0.757 , 0.587)],],
+    basis = '6-31g',
+    symmetry = 1)
+
+mf = scf.RHF(mol)
+mf.kernel()
+
+mc = mcpdft.CASSCF(mf, 'tPBE', 4, 4).state_average_([.64,.36]).run (verbose=4)
+
+print ("Average MC-PDFT energy =", mc.e_tot)
+print ("E_PDFT-E_CASSCF for state 0 =", mc.e_states[0]-mc.e_mcscf[0])
+print ("E_PDFT-E_CASSCF for state 1 =", mc.e_states[1]-mc.e_mcscf[1])
+
+
diff --git a/examples/mcpdft/16-multi_state.py b/examples/mcpdft/16-multi_state.py
new file mode 100644
index 000000000..6ec8c4a44
--- /dev/null
+++ b/examples/mcpdft/16-multi_state.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+
+'''
+multi-state
+
+A "quasidegenerate" extension of MC-PDFT for states close in energy. Currently
+only "compressed multi-state" (CMS-) [JCTC 16, 7444 (2020)] and "extended multi-
+state" (XMS-) [Faraday Discuss 224, 348-372 (2020)] is supported.
+'''
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M(
+    atom = [
+        ['Li', ( 0., 0.    , 0.   )],
+        ['H', ( 0., 0., 3)]
+    ], basis = 'sto-3g',
+    symmetry = 0 # symmetry enforcement is not recommended for MS-PDFT
+    )
+
+mf = scf.RHF(mol)
+mf.kernel()
+
+mc = mcpdft.CASSCF(mf, 'tpbe', 2, 2)
+mc.fix_spin_(ss=0) # often necessary!
+mc_sa = mc.state_average ([.5, .5]).run ()
+
+# For CMS-PDFT
+mc_cms = mc.multi_state([.5, .5], "cms").run(verbose=4)
+
+# For XMS-PDFT
+mc_xms = mc.multi_state([.5, .5], "xms").run(verbose=4)
+
+print ('{:>21s} {:>12s} {:>12s}'.format ('state 0','state 1', 'gap'))
+fmt_str = '{:>9s} {:12.9f} {:12.9f} {:12.9f}'
+print (fmt_str.format ('CASSCF', mc_sa.e_mcscf[0], mc_sa.e_mcscf[1],
+    mc_sa.e_mcscf[1]-mc_sa.e_mcscf[0]))
+print (fmt_str.format ('MC-PDFT', mc_sa.e_states[0], mc_sa.e_states[1],
+    mc_sa.e_states[1]-mc_sa.e_states[0]))
+print (fmt_str.format ('CMS-PDFT', mc_cms.e_states[0], mc_cms.e_states[1],
+    mc_cms.e_states[1]-mc_cms.e_states[0]))
+print (fmt_str.format ('XMS-PDFT', mc_xms.e_states[0], mc_xms.e_states[1],
+    mc_xms.e_states[1]-mc_xms.e_states[0]))
diff --git a/examples/mcpdft/41-state_average.py b/examples/mcpdft/41-state_average.py
new file mode 100644
index 000000000..693e22422
--- /dev/null
+++ b/examples/mcpdft/41-state_average.py
@@ -0,0 +1,80 @@
+#!/usr/bin/env python
+
+'''
+State average over states of different spins and/or spatial symmetry
+
+In MC-PDFT, you must use the method attribute rather than the function:
+
+NO  : mc = mcscf.state_average_mix (mc, [solver1, solver2], weights)
+NO  : mcscf.state_average_mix_(mc, [solver1, solver2], weights)
+YES : mc = mc.state_average_mix ([solver1, solver2], weights)
+YES : mc.state_average_mix_([solver1, solver2], weights)
+'''
+
+import numpy as np
+from pyscf import gto, scf, mcpdft, fci
+
+r = 1.8
+mol = gto.Mole()
+mol.atom = [
+    ['C', ( 0., 0.    , -r/2   )],
+    ['C', ( 0., 0.    ,  r/2)],]
+mol.basis = 'cc-pvdz'
+mol.unit = 'B'
+mol.symmetry = True
+mol.build()
+mf = scf.RHF(mol)
+mf.irrep_nelec = {'A1g': 4, 'E1gx': 0, 'E1gy': 0, 'A1u': 4,
+                  'E1uy': 2, 'E1ux': 2, 'E2gx': 0, 'E2gy': 0, 'E2uy': 0, 'E2ux': 0}
+ehf = mf.kernel()
+#mf.analyze()
+
+#
+# state-average over 1 triplet + 2 singlets
+# Note direct_spin1 solver is called here because the CI solver will take
+# spin-mix solution as initial guess which may break the spin symmetry
+# required by direct_spin0 solver
+#
+weights = np.ones(3)/3
+solver1 = fci.direct_spin1_symm.FCI(mol)
+solver1.spin = 2
+solver1 = fci.addons.fix_spin(solver1, shift=.2, ss=2)
+solver1.nroots = 1
+solver2 = fci.direct_spin0_symm.FCI(mol)
+solver2.spin = 0
+solver2.nroots = 2
+
+mc = mcpdft.CASSCF(mf, 'tPBE', 8, 8)
+mc.state_average_mix_([solver1, solver2], weights)
+
+# Mute warning msgs
+mc.check_sanity = lambda *args: None
+
+mc.verbose = 4
+mc.kernel()
+
+
+#
+# Example 2: Mix FCI wavefunctions with different symmetry irreps
+#
+mol = gto.Mole()
+mol.build(atom='''
+ O     0.    0.000    0.1174
+ H     0.    0.757   -0.4696
+ H     0.   -0.757   -0.4696
+''', symmetry=True, basis='631g')
+
+mf = scf.RHF(mol).run()
+
+weights = [.5, .5]
+solver1 = fci.direct_spin1_symm.FCI(mol)
+solver1.wfnsym= 'A1'
+solver1.spin = 0
+solver2 = fci.direct_spin1_symm.FCI(mol)
+solver2.wfnsym= 'A2'
+solver2.spin = 0
+
+mc = mcpdft.CASSCF(mf, 'tPBE', 4, 4)
+mc.state_average_mix_([solver1, solver2], weights)
+mc.verbose = 4
+mc.kernel()
diff --git a/examples/mcpdft/42-linearized-pdft.py b/examples/mcpdft/42-linearized-pdft.py
new file mode 100644
index 000000000..51a51c6a0
--- /dev/null
+++ b/examples/mcpdft/42-linearized-pdft.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+'''
+linearized pair-density functional theory
+
+A multi-states extension of MC-PDFT for states close in energy. Generates an
+effective L-PDFT Hamiltonian through a Taylor expansion of the MC-PDFT energy
+expression [J Chem Theory Comput: 10.1021/acs.jctc.3c00207].
+'''
+
+from pyscf import gto, scf, mcpdft
+
+mol = gto.M(
+    atom = [
+        ['Li', ( 0., 0., 0.)],
+        ['H', ( 0., 0., 3)]
+    ], basis = 'sto-3g',
+    symmetry = 0
+    )
+
+mf = scf.RHF(mol)
+mf.kernel()
+
+mc = mcpdft.CASSCF(mf, 'tpbe', 2, 2)
+mc.fix_spin_(ss=0) # often necessary!
+mc_sa = mc.state_average ([.5, .5]).run ()
+# Note, this operates in the same way as "CMS-PDFT" but constructs and
+# diagonalizes a different matrix.
+mc_ms = mc.multi_state ([.5, .5], "lin").run (verbose=4)
+
+print ('{:>21s} {:>12s} {:>12s}'.format ('state 0','state 1', 'gap'))
+fmt_str = '{:>9s} {:12.9f} {:12.9f} {:12.9f}'
+print (fmt_str.format ('CASSCF', mc_sa.e_mcscf[0], mc_sa.e_mcscf[1],
+    mc_sa.e_mcscf[1]-mc_sa.e_mcscf[0]))
+print (fmt_str.format ('MC-PDFT', mc_sa.e_states[0], mc_sa.e_states[1],
+    mc_sa.e_states[1]-mc_sa.e_states[0]))
+print (fmt_str.format ('L-PDFT', mc_ms.e_states[0], mc_ms.e_states[1],
+    mc_ms.e_states[1]-mc_ms.e_states[0]))
diff --git a/examples/mcpdft/43-multi_state_mix.py b/examples/mcpdft/43-multi_state_mix.py
new file mode 100644
index 000000000..235f1064c
--- /dev/null
+++ b/examples/mcpdft/43-multi_state_mix.py
@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+
+
+'''
+Perform multi-state PDFT averaging over states of different spins and/or
+spatial symmetry
+
+The mcpdft.multi_state function maybe not generate the right spin or spatial
+symmetry as one needs.  This example shows how to put states with different
+spins and spatial symmetry in a state-average solver using the method
+multi_state_mix.
+'''
+
+import numpy as np
+from pyscf import gto, scf, mcpdft, fci
+
+mol = gto.M(atom='''
+   O     0.    0.000    0.1174
+   H     0.    0.757   -0.4696
+   H     0.   -0.757   -0.4696
+      ''', symmetry=True, basis="6-31G", verbose=3)
+
+mf = scf.RHF(mol).run()
+
+#
+# state-average over 1 triplet + 2 singlets
+# Note direct_spin1 solver is called here because the CI solver will take
+# spin-mix solution as initial guess which may break the spin symmetry
+# required by direct_spin0 solver
+#
+weights = np.ones(3) / 3
+solver1 = fci.direct_spin1_symm.FCI(mol)
+solver1.spin = 2
+solver1 = fci.addons.fix_spin(solver1, shift=.2, ss=2)
+solver1.nroots = 1
+solver2 = fci.direct_spin0_symm.FCI(mol)
+solver2.spin = 0
+solver2.nroots = 2
+
+mc = mcpdft.CASSCF(mf, "tPBE", 4, 4, grids_level=1)
+# Currently, only the Linearized PDFT method is available for
+# multi_state_mix
+mc = mc.multi_state_mix([solver1, solver2], weights, "lin")
+mc.run()
+
+
+#
+# Example 2: Mix FCI wavefunctions with different symmetry irreps
+#
+mol = gto.Mole()
+mol.build(atom='''
+ O     0.    0.000    0.1174
+ H     0.    0.757   -0.4696
+ H     0.   -0.757   -0.4696
+''', symmetry=True, basis='631g', verbose=3)
+
+mf = scf.RHF(mol).run()
+
+# Also possible to construct 2 solvers of different wfnsym, but
+# of the same spin symmetry
+weights = [.5, .5]
+solver1 = fci.direct_spin1_symm.FCI(mol)
+solver1.wfnsym= 'A1'
+solver1.spin = 0
+solver2 = fci.direct_spin1_symm.FCI(mol)
+solver2.wfnsym= 'A2'
+solver2.spin = 2
+
+mc = mcpdft.CASSCF(mf, "tPBE", 4, 4, grids_level=1)
+mc = mc.multi_state_mix([solver1, solver2], weights, "lin")
+mc.kernel()
diff --git a/source/modules.rst b/source/modules.rst
index a0307d406..f35808780 100644
--- a/source/modules.rst
+++ b/source/modules.rst
@@ -27,6 +27,7 @@ API documentation
    modules/lib.rst
    modules/lo.rst
    modules/mcscf.rst
+   modules/mcpdft.rst
    modules/md.rst
    modules/mp.rst
    modules/mrpt.rst
diff --git a/source/modules/mcpdft.rst b/source/modules/mcpdft.rst
new file mode 100644
index 000000000..15ad31c78
--- /dev/null
+++ b/source/modules/mcpdft.rst
@@ -0,0 +1,72 @@
+.. _mcpdft:
+
+:mod:`mcpdft` --- Multi-configuration pair-density functional theory
+********************************************************************
+
+.. module:: mcpdft
+   :synopsis: Multi-configuration pair-density functional theory
+
+The :mod:`mcpdft` module implements multi-configuration pair-density functional theory.
+
+
+Examples
+========
+
+* :source:`examples/mcpdft/00-simple_mcpdft.py`
+* :source:`examples/mcpdft/01-different_functionals.py`
+* :source:`examples/mcpdft/02-hybrid_functionals.py`
+* :source:`examples/mcpdft/03-metaGGA_functionals.py`
+* :source:`examples/mcpdft/11-grid_scheme.py`
+* :source:`examples/mcpdft/15-state_average.py`
+* :source:`examples/mcpdft/16-multi_state.py`
+* :source:`examples/mcpdft/41-state_average.py`
+* :source:`examples/mcpdft/42-linearized-pdft.py`
+* :source:`examples/mcpdft/43-multi_state_mix.py`
+
+Program reference
+=================
+
+.. automodule:: pyscf.mcpdft
+   :members:
+
+.. automodule:: pyscf.mcpdft.mcpdft
+   :members:
+
+.. automodule:: pyscf.mcpdft._dms
+   :members:
+
+.. automodule:: pyscf.mcpdft.otpd
+   :members:
+
+.. automodule:: pyscf.mcpdft.otfnal
+   :members:
+
+.. automodule:: pyscf.mcpdft._libxc
+   :members:
+
+.. automodule:: pyscf.mcpdft.tfnal_derivs
+   :members:
+
+.. automodule:: pyscf.mcpdft.pdft_eff
+   :members:
+
+.. automodule:: pyscf.mcpdft.pdft_veff
+   :members:
+
+.. automodule:: pyscf.mcpdft.pdft_feff
+   :members:
+
+.. automodule:: pyscf.mcpdft.mspdft
+   :members:
+
+.. automodule:: pyscf.mcpdft.cmspdft
+   :members:
+
+.. automodule:: pyscf.mcpdft.xmspdft
+   :members:
+
+.. automodule:: pyscf.mcpdft.lpdft
+   :members:
+
+.. automodule:: pyscf.mcpdft.chkfile
+   :members:
diff --git a/source/user/index.md b/source/user/index.md
index 51401aca0..8018b5187 100644
--- a/source/user/index.md
+++ b/source/user/index.md
@@ -34,6 +34,7 @@ cc.rst
 adc.rst
 agf2.rst
 mcscf.rst
+mcpdft.rst
 mrpt.rst
 tddft.rst
 ```
diff --git a/source/user/mcpdft.rst b/source/user/mcpdft.rst
new file mode 100644
index 000000000..88d9785e8
--- /dev/null
+++ b/source/user/mcpdft.rst
@@ -0,0 +1,114 @@
+.. _user_mcpdft:
+
+Multi-configuration pair-density functional theory (MC-PDFT)
+************************************************************
+
+*Modules*: :py:mod:`pyscf.mcpdft`
+
+
+Introduction
+------------
+
+Multiconfiguration pair-density functional theory (MC-PDFT) refers to methods in which the total electronic energy is obtained or derived from a functional of the on-top pair density as well as the total electron density, and these densities are in turn obtained from a multiconfigurational self-consistent field (MCSCF) wave function or wave functions of some sort.
+The on-top pair density is the probability of two electrons existing simultaneously at one point in space.
+The non-classical part of the MC-PDFT energy is called the "on-top energy" and is directly analogous to the exchange-correlation energy in standard density functional theory.
+In practice, the on-top energy functional is always a generalization or a "translation" of an exchange-correlation functional from Kohn-Sham DFT.
+MC-PDFT furnishes a practical and effective generalization of KS-DFT functionals from single Slater determinants to general multiconfigurational wave functions because the pair density, unlike the spin density, can simultaneously encode both strong electron correlation effects and spin symmetry constraints.
+PySCF implements the MC-PDFT protocol of Reference :cite:`Manni2014`, in which the molecular orbital and/or configuration interaction coefficients are chosen to minimize the expectation value of the molecular Hamiltonian in the underlying MCSCF wave function, not the MC-PDFT energy itself.
+
+The relationship of :mod:`mcpdft` to :mod:`mcscf` is roughly similar to the relationship between :mod:`dft` and :mod:`scf`, in that instances of MC-PDFT methods are also instances of MCSCF methods.
+:mod:`mcpdft` thus inherits :mod:`mcscf`'s compatibility with density fitting (:ref:`user_df`) and x2c (:ref:`user_x2c`) as well as CASSCF's ``state_average_`` and ``state_average_mix_`` extensions.
+As in :mod:`dft`, the MC-PDFT object has a ``grids`` attribute which stores information about the quadrature grid.
+The on-top functional is identified with a string stored in the ``otxc`` attribute of the MC-PDFT object, which is passed as the required second argument to the method constructor; see below.
+
+CASCI-PDFT
+""""""""""
+.. code-block:: python
+
+  from pyscf import gto, mp, mcpdft
+  mol = gto.M(
+      atom = 'O 0 0 0; O 0 0 1.2',
+      basis = 'ccpvdz',
+      spin = 2)
+  myhf = mol.RHF().run()
+  # Use MP2 natural orbitals to define the active space for the single-point CAS-CI calculation
+  mymp = mp.UMP2(myhf).run()
+
+  noons, natorbs = mcscf.addons.make_natural_orbitals(mymp)
+  ncas, nelecas = (6,8)
+  otfnal = 'tPBE'
+  mycas = mcpdft.CASCI(myhf, otfnal, ncas, nelecas)
+  mycas.kernel(natorbs)
+
+
+CASSCF-PDFT
+"""""""""""
+.. code-block:: python
+
+  import pyscf
+  from pyscf import mcpdft
+  mol = pyscf.M(
+      atom = 'O 0 0 0; O 0 0 1.2',
+      basis = 'ccpvdz',
+      spin = 2)
+  myhf = mol.RHF().run()
+  ncas, nelecas = (6,(5,3))
+  otfnal = 'tPBE'
+  # For convenience, mycas.grids.level can be set at construction using an optional kwarg
+  mycas = mcpdft.CASSCF(myhf, otfnal, ncas, nelecas, grids_level=6).run()
+
+When state-averaged CASSCF wave functions are used, the MC-PDFT energy is computed separately for each state.
+When either state-averaged CASSCF or multiple roots with CASCI are calculated, the ordinal label of each state is as assigned by the underlying MCSCF calculation, regardless of the relative MC-PDFT energies.
+In other words, state 2 may sometimes have a lower MC-PDFT energy than state 1.
+
+Available on-top functionals
+----------------------------
+
+"Translated" on-top functionals corresponding to any pure LDA or GGA exchange-correlation functional :cite:`Manni2014`, as well as meta-GGA functionals :cite:`Bao2025` which depend on the kinetic energy, are available.
+These are specified by prepending the letter ``t`` to the name of the base DFT functional, as in ``tPBE`` or ``tBLYP``.
+"Fully-translated" on-top functionals :cite:`Carlson2015`  corresponding to pure LDA or GGA functionals are also available and are specified by prepending the letters ``ft`` to the name of the base functional.
+
+With the exception of PBE0, for which the translated and fully-translated generalizations are implemented precisely as stipulated in Reference :cite:`Pandharkar2020`, prepending ``t`` or ``ft`` to the names of hybrid DFT functionals in this way will not work, because there are multiple forms of hybrid on-top functional in the literature (for instance, compare References :cite:`Pandharkar2020` and :cite:`Mostafanejad2020`).
+However, it is still possible to specify a variety of global hybrid on-top functionals:
+
+.. code-block:: python
+
+  from pyscf import mcpdft
+
+  # Construct a custom hybrid functional with 20% CASSCF energy (JPCL 11, 10158, 2020) 
+  
+  my_otxc = 't' + mcpdft.hyb ('BLYP', .2)
+  mc = mcpdft.CASSCF (mf, my_otxc, 6, 8).run ()
+  
+  # "lambda-MC-PDFT" of JCTC 16, 2274, 2020.
+  
+  my_otxc = 't' + mcpdft.hyb('PBE', .5, hyb_type='lambda')
+  mc.run (otxc=my_otxc)
+
+
+Multi-state extension methods
+-----------------------------
+
+MC-PDFT evaluates the total energy of a single electronic state through a nonlinear functional expression involving density variables.
+For approximate functionals evaluated at approximate densities, it cannot be guaranteed that potential energy surfaces computed in this way for states of the same symmetry only cross along `3N-8`-dimensional conical intersection seams, which is a requirement for qualitatively accurate non-adiabatic molecular dynamics simulations.
+Multi-state methods are extensions of MC-PDFT which address this difficulty by defining the electronic energies of a set of states as eigenvalues of a small Hamiltonian matrix, of which the elements are related in one way or another to the MC-PDFT energy functional expression.
+In PySCF, multi-state MC-PDFT methods are implemented as child classes of state-averaged CASSCF methods.
+
+The preferred multi-state MC-PDFT method is L-PDFT :cite:`Hennefarth2023`, which is performed in the following way:
+
+.. code-block:: python
+
+  from pyscf import mcpdft
+
+  mc = mcpdft.CASSCF(mf, 'tPBE', 2, 2)
+  mc.fix_spin_(ss=0) # often necessary!
+  mc_ms = mc.multi_state ([.5, .5], method="lin")
+  mc_ms.kernel () # Performs SA-CASSCF orbital optimization followed by L-PDFT
+
+Note that the first argument of ``multi_state`` is a sequence of weights, exactly as in ``state_average_``.
+The second argument stipulates the type of multi-state MC-PDFT.
+XMS-PDFT :cite:`Bao2020_XMS` and CMS-PDFT :cite:`Bao2020_CMS` are available by specifying ``method="XMS"`` and ``method="CMS"``, respectively, in place of ``method="lin"``.
+For L-PDFT only, it is also possible to carry out calculations involving multiple different point group symmetries or spins using ``multi_state_mix``, similar to ``state_average_mix_``.
+
+
+
diff --git a/source/user/ref.bib b/source/user/ref.bib
index 64a795e78..94d85e101 100644
--- a/source/user/ref.bib
+++ b/source/user/ref.bib
@@ -898,3 +898,93 @@ @article{VanLenthe2006_JCC_926
   doi = {https://doi.org/10.1002/jcc.20393},
   year = {2006}
 }
+
+@article{Manni2014,
+  author = {Manni, Giovanni Li and Carlson, Rebecca K and Luo, Sijie and Ma, Dongxia and Olsen, Jeppe and Truhlar, Donald G and Gagliardi, Laura},
+  doi = {10.1021/ct500483t},
+  journal = {J. Chem. Theory Comput.},
+  keywords = {MC-PDFT},
+  pages = {3669},
+  title = {{Multiconfiguration Pair-Density Functional Theory}},
+  volume = {10},
+  year = {2014}
+}
+
+@article{Carlson2015,
+  author = {Carlson, Rebecca K and Truhlar, Donald G and Gagliardi, Laura},
+  doi = {10.1021/acs.jctc.5b00609},
+  journal = {J. Chem. Theory Comput.},
+  keywords = {MC-PDFT},
+  pages = {4077--4085},
+  title = {{Multiconfiguration Pair-Density Functional Theory: A Fully Translated Gradient Approximation and Its Performance for Transition Metal Dimers and the Spectroscopy of Re2Cl82−}},
+  volume = {11},
+  year = {2015}
+}
+
+@article{Pandharkar2020,
+  author = {Pandharkar, Riddhish and Hermes, Matthew R. and Truhlar, Donald G and Gagliardi, Laura},
+  doi = {10.1021/acs.jpclett.0c02956},
+  journal = {J. Phys. Chem. Lett.},
+  pages = {10158--10163},
+  title = {{A New Mixing of Nonlocal Exchange and Nonlocal Correlation with Multiconfiguration Pair-Density Functional Theory}},
+  volume = {11},
+  year = {2020}
+}
+
+@article{Mostafanejad2020,
+  author = {Mostafanejad, Mohammad and {Dante Liebenthal}, Marcus and {Eugene DePrince III}, A},
+  doi = {10.1021/acs.jctc.9b01178},
+  journal = {J. Chem. Theory Comput.},
+  keywords = {DFT functional,MC-PDFT,Mar20 lit search},
+  number = {4}, 
+  pages = {2274--2283},
+  title = {{Global Hybrid Multiconfiguration Pair-Density Functional Theory}},
+  volume = {16},
+  year = {2020}
+}
+
+@article{Bao2025,
+  author = {Bao, Jie J. and Zhang, Dayou and Zhang, Shaoting and Gagliardi, Laura and Truhlar, Donald G.},
+  doi = {10.1073/pnas.2419413121},
+  journal = {Proc. Natl. Acad. Sci. USA},
+  number = {1},
+  pages = {e2419413121},
+  title = {{A hybrid meta on-­top functional for multiconfiguration pair-density functional theory}},
+  volume = {122},
+  year = {2025}
+}
+
+@article{Bao2020_CMS,
+  author = {Bao, Jie J. and Zhou, Chen and Truhlar, Donald G.},
+  doi = {10.1021/acs.jctc.0c00908},
+  issn = {15499626},
+  journal = {J. Chem. Theory. Comput.},
+  number = {12},
+  pages = {7444--7452},
+  title = {{Compressed-State Multistate Pair-Density Functional Theory}},
+  volume = {16},
+  year = {2020}
+}
+
+@article{Bao2020_XMS,
+  author = {Bao, Jie J. and Zhou, Chen and Varga, Zoltan and Kanchanakungwankul, Siriluk and Gagliardi, Laura and Truhlar, Donald G.},
+  doi = {10.1039/d0fd00037j},
+  journal = {Faraday Discuss.},
+  pages = {348--372},
+  title = {{Multi-state pair-density functional theory}},
+  volume = {224},
+  year = {2020}
+}
+
+@article{Hennefarth2023,
+  author = {Hennefarth, Matthew R. and Hermes, Matthew R. and Truhlar, Donald G. and Gagliardi, Laura},
+  doi = {10.1021/acs.jctc.3c00207},
+  journal = {J. Chem. Theory Comput.},
+  number = {11},
+  pages = {3172--3183},
+  title = {{Linearized Pair-Density Functional Theory}},
+  volume = {19},
+  year = {2023}
+}
+
+