diff --git a/include/chemist/density/decomposable_density.hpp b/include/chemist/density/decomposable_density.hpp
index 1e89d049d..9243cba65 100644
--- a/include/chemist/density/decomposable_density.hpp
+++ b/include/chemist/density/decomposable_density.hpp
@@ -15,7 +15,7 @@
  */
 
 #pragma once
-#include <chemist/density/density.hpp>
+#include <chemist/density/density_class.hpp>
 #include <chemist/quantum_mechanics/wavefunction/transformed.hpp>
 
 namespace chemist {
@@ -172,4 +172,6 @@ bool operator!=(const DecomposableDensity<LHSParticles...>& lhs,
 /// Type of the one-Electron density
 using Decomposable1EDensity = DecomposableDensity<Electron>;
 
+extern template class DecomposableDensity<Electron>;
+
 } // namespace chemist
diff --git a/include/chemist/density/density.hpp b/include/chemist/density/density.hpp
index 3efdbb99d..af7d138e7 100644
--- a/include/chemist/density/density.hpp
+++ b/include/chemist/density/density.hpp
@@ -1,5 +1,5 @@
 /*
- * Copyright 2022 NWChemEx-Project
+ * Copyright 2024 NWChemEx-Project
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
  * you may not use this file except in compliance with the License.
@@ -15,180 +15,6 @@
  */
 
 #pragma once
-#include <chemist/electron/electron.hpp>
-#include <chemist/quantum_mechanics/wavefunction/aos.hpp>
-#include <tensorwrapper/tensorwrapper.hpp>
-
-namespace chemist {
-
-/** @brief Class describing a density.
- *
- *  In Electronic structure theory we place a lot of emphasis on the
- *  one-Electron density; however, there are many other densities. This class
- *  is templated on the types and number of particles that the density
- *  describes. In turn this allows for this class to model densities involving
- *  an arbitrary number of arbitrary particles.
- *
- *  @note At the moment this class is basically a strong type which wraps a
- *        tensor and an AOs object.
- *
- *  @note If we want to support Particles... != Electron... we will need to add
- *        TMP to work out the basis type.
- *
- *  @tparam Particles The types of the particles described by this density. At
- *                    the moment the types are used purely for strong typing
- *                    purposes. Expected to be unqualified types.
- *
- */
-template<typename... Particles>
-class Density {
-public:
-    /// Type used to represent the numeric form of the density
-    using value_type = tensorwrapper::Tensor;
-
-    /// Type of a mutable reference to the density matrix
-    using reference = value_type&;
-
-    /// Type of a read-only reference to the density matrix
-    using const_reference = const value_type&;
-
-    /// Type of the basis set for value_type
-    using basis_type = wavefunction::AOs;
-
-    /// Type of a mutable reference to an object of type basis_type
-    using basis_reference = basis_type&;
-
-    /// Type of a read-only reference to an object of basis_type
-    using const_basis_reference = const basis_type&;
-
-    /// Default ctor
-    explicit Density() = default;
-
-    /** @brief Creates a new density object with the provided values.
-     *
-     *  @param[in] rho The tensor representation of the density matrix.
-     *  @param[in] orbs The basis set the density matrix is expressed in.
-     *
-     *  @throw None No throw guarantee.
-     */
-    Density(value_type rho, basis_type orbs) :
-      m_density_(std::move(rho)), m_orbs_(std::move(orbs)) {}
-
-    /** @brief Returns a mutable reference to the density matrix.
-     *
-     *  @return A mutable reference to the tensor representation of the density.
-     *
-     *  @throw None No throw guarantee.
-     */
-    reference value() noexcept { return m_density_; }
-
-    /** @brief Provides read-only access to the density matrix.
-     *
-     *  Same as the non-const version except that the return is read-only. See
-     *  description of non-const version for more details.
-     *
-     *  @return The density value.
-     *
-     *  @throw None No throw guarantee.
-     */
-    const_reference value() const noexcept { return m_density_; }
-
-    /** @brief Returns a mutable reference to the orbital basis set for the
-     *         density.
-     *
-     *  @return A mutable reference to the basis set for *this.
-     *
-     *  @throw None No throw guarantee.
-     */
-    basis_reference basis_set() noexcept { return m_orbs_; }
-
-    /** @brief Provides read-only access to the orbital basis set used to define
-     *         the density matrix.
-     *
-     *  Same as the non-const version except that the return is read-only. See
-     *  description of the non-const version for more details.
-     *
-     *  @return The orbital basis the density is formed in.
-     *
-     *  @throw ??? Throws if the derived class's implementation of `basis_set_`
-     *             throws. Same throw guarantee.
-     */
-    const_basis_reference basis_set() const { return basis_set_(); }
-
-protected:
-    /// The density tensor
-    value_type m_density_;
-
-    /// For assigning density value in derived types
-    explicit Density(value_type rho) : m_density_(rho) {}
-
-    /** @brief To be overridden by the derived class to implement basis_set().
-     *
-     *  The derived class is responsible for ensuring this function returns the
-     *  basis set.
-     *
-     *  @return The orbital basis the density is formed in.
-     *
-     *  @throw None No throw guarantee.
-     */
-    virtual const_basis_reference basis_set_() const { return m_orbs_; }
-
-private:
-    /// The orbital space used to make the density
-    basis_type m_orbs_;
-};
-
-/** @brief Compares two Density instances for equality.
- *
- *  @relates Density
- *
- *  Two Density instances are equal if they describe the same particles and they
- *  have the same numeric value.
- *
- *  @tparam LHSParticles The particle types for @p lhs.
- *  @tparam RHSParticles The particle types for @p rhs.
- *
- *  @param[in] lhs The density on the left side of the equality operator.
- *  @param[in] rhs The density on the right side of the equality operator.
- *
- *  @return True if @p lhs has the same type and numeric value as @p rhs. False
- *          otherwise.
- *
- *  @throw None No throw guarantee.
- */
-template<typename... LHSParticles, typename... RHSParticles>
-bool operator==(const Density<LHSParticles...>& lhs,
-                const Density<RHSParticles...>& rhs) {
-    if constexpr(std::is_same_v<decltype(lhs), decltype(rhs)>) {
-        return std::tie(lhs.basis_set(), lhs.value()) ==
-               std::tie(rhs.basis_set(), rhs.value());
-    } else {
-        return false;
-    }
-}
-
-/** @brief Determines if two Density instances are different.
- *
- *  @relates Density
- *
- *  Two Densities are different if the equality comparison returns false. This
- *  indicates that the densities either have different types and/or different
- *  numeric values.
- *
- *  @tparam LHSParticles The particle types for @p lhs.
- *  @tparam RHSParticles The particle types for @p rhs.
- *
- *  @param[in] lhs The density on the left side of the not-equal operator.
- *  @param[in] rhs The density on the right side of the not-equal operator.
- *
- *  @return False if @p lhs and @p rhs are equal and true otherwise.
- *
- *  @throw None No throw guarantee.
- */
-template<typename... LHSParticles, typename... RHSParticles>
-bool operator!=(const Density<LHSParticles...>& lhs,
-                const Density<RHSParticles...>& rhs) {
-    return !(lhs == rhs);
-}
-
-} // namespace chemist
+#include <chemist/density/decomposable_density.hpp>
+#include <chemist/density/density_class.hpp>
+#include <chemist/traits/density_traits.hpp>
\ No newline at end of file
diff --git a/include/chemist/density/density_class.hpp b/include/chemist/density/density_class.hpp
new file mode 100644
index 000000000..9a6d73080
--- /dev/null
+++ b/include/chemist/density/density_class.hpp
@@ -0,0 +1,196 @@
+/*
+ * Copyright 2022 NWChemEx-Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+#include <chemist/electron/electron.hpp>
+#include <chemist/quantum_mechanics/wavefunction/aos.hpp>
+#include <tensorwrapper/tensorwrapper.hpp>
+
+namespace chemist {
+
+/** @brief Class describing a density.
+ *
+ *  In Electronic structure theory we place a lot of emphasis on the
+ *  one-Electron density; however, there are many other densities. This class
+ *  is templated on the types and number of particles that the density
+ *  describes. In turn this allows for this class to model densities involving
+ *  an arbitrary number of arbitrary particles.
+ *
+ *  @note At the moment this class is basically a strong type which wraps a
+ *        tensor and an AOs object.
+ *
+ *  @note If we want to support Particles... != Electron... we will need to add
+ *        TMP to work out the basis type.
+ *
+ *  @tparam Particles The types of the particles described by this density. At
+ *                    the moment the types are used purely for strong typing
+ *                    purposes. Expected to be unqualified types.
+ *
+ */
+template<typename... Particles>
+class Density {
+public:
+    /// Type used to represent the numeric form of the density
+    using value_type = tensorwrapper::Tensor;
+
+    /// Type of a mutable reference to the density matrix
+    using reference = value_type&;
+
+    /// Type of a read-only reference to the density matrix
+    using const_reference = const value_type&;
+
+    /// Type of the basis set for value_type
+    using basis_type = wavefunction::AOs;
+
+    /// Type of a mutable reference to an object of type basis_type
+    using basis_reference = basis_type&;
+
+    /// Type of a read-only reference to an object of basis_type
+    using const_basis_reference = const basis_type&;
+
+    /// Default ctor
+    explicit Density() = default;
+
+    /** @brief Creates a new density object with the provided values.
+     *
+     *  @param[in] rho The tensor representation of the density matrix.
+     *  @param[in] orbs The basis set the density matrix is expressed in.
+     *
+     *  @throw None No throw guarantee.
+     */
+    Density(value_type rho, basis_type orbs) :
+      m_density_(std::move(rho)), m_orbs_(std::move(orbs)) {}
+
+    /** @brief Returns a mutable reference to the density matrix.
+     *
+     *  @return A mutable reference to the tensor representation of the density.
+     *
+     *  @throw None No throw guarantee.
+     */
+    reference value() noexcept { return m_density_; }
+
+    /** @brief Provides read-only access to the density matrix.
+     *
+     *  Same as the non-const version except that the return is read-only. See
+     *  description of non-const version for more details.
+     *
+     *  @return The density value.
+     *
+     *  @throw None No throw guarantee.
+     */
+    const_reference value() const noexcept { return m_density_; }
+
+    /** @brief Returns a mutable reference to the orbital basis set for the
+     *         density.
+     *
+     *  @return A mutable reference to the basis set for *this.
+     *
+     *  @throw None No throw guarantee.
+     */
+    basis_reference basis_set() noexcept { return m_orbs_; }
+
+    /** @brief Provides read-only access to the orbital basis set used to define
+     *         the density matrix.
+     *
+     *  Same as the non-const version except that the return is read-only. See
+     *  description of the non-const version for more details.
+     *
+     *  @return The orbital basis the density is formed in.
+     *
+     *  @throw ??? Throws if the derived class's implementation of `basis_set_`
+     *             throws. Same throw guarantee.
+     */
+    const_basis_reference basis_set() const { return basis_set_(); }
+
+protected:
+    /// The density tensor
+    value_type m_density_;
+
+    /// For assigning density value in derived types
+    explicit Density(value_type rho) : m_density_(rho) {}
+
+    /** @brief To be overridden by the derived class to implement basis_set().
+     *
+     *  The derived class is responsible for ensuring this function returns the
+     *  basis set.
+     *
+     *  @return The orbital basis the density is formed in.
+     *
+     *  @throw None No throw guarantee.
+     */
+    virtual const_basis_reference basis_set_() const { return m_orbs_; }
+
+private:
+    /// The orbital space used to make the density
+    basis_type m_orbs_;
+};
+
+/** @brief Compares two Density instances for equality.
+ *
+ *  @relates Density
+ *
+ *  Two Density instances are equal if they describe the same particles and they
+ *  have the same numeric value.
+ *
+ *  @tparam LHSParticles The particle types for @p lhs.
+ *  @tparam RHSParticles The particle types for @p rhs.
+ *
+ *  @param[in] lhs The density on the left side of the equality operator.
+ *  @param[in] rhs The density on the right side of the equality operator.
+ *
+ *  @return True if @p lhs has the same type and numeric value as @p rhs. False
+ *          otherwise.
+ *
+ *  @throw None No throw guarantee.
+ */
+template<typename... LHSParticles, typename... RHSParticles>
+bool operator==(const Density<LHSParticles...>& lhs,
+                const Density<RHSParticles...>& rhs) {
+    if constexpr(std::is_same_v<decltype(lhs), decltype(rhs)>) {
+        return std::tie(lhs.basis_set(), lhs.value()) ==
+               std::tie(rhs.basis_set(), rhs.value());
+    } else {
+        return false;
+    }
+}
+
+/** @brief Determines if two Density instances are different.
+ *
+ *  @relates Density
+ *
+ *  Two Densities are different if the equality comparison returns false. This
+ *  indicates that the densities either have different types and/or different
+ *  numeric values.
+ *
+ *  @tparam LHSParticles The particle types for @p lhs.
+ *  @tparam RHSParticles The particle types for @p rhs.
+ *
+ *  @param[in] lhs The density on the left side of the not-equal operator.
+ *  @param[in] rhs The density on the right side of the not-equal operator.
+ *
+ *  @return False if @p lhs and @p rhs are equal and true otherwise.
+ *
+ *  @throw None No throw guarantee.
+ */
+template<typename... LHSParticles, typename... RHSParticles>
+bool operator!=(const Density<LHSParticles...>& lhs,
+                const Density<RHSParticles...>& rhs) {
+    return !(lhs == rhs);
+}
+
+extern template class Density<Electron>;
+
+} // namespace chemist
diff --git a/include/chemist/quantum_mechanics/operator/coulomb.hpp b/include/chemist/quantum_mechanics/operator/coulomb.hpp
index 0f882626b..a768d5e1e 100644
--- a/include/chemist/quantum_mechanics/operator/coulomb.hpp
+++ b/include/chemist/quantum_mechanics/operator/coulomb.hpp
@@ -113,6 +113,10 @@ class Coulomb : public detail_::OperatorImpl<Coulomb<LHSParticle, RHSParticle>,
 
 extern template class Coulomb<Electron, Electron>;
 extern template class Coulomb<ManyElectrons, ManyElectrons>;
+extern template class Coulomb<Electron, chemist::Density<Electron>>;
+extern template class Coulomb<ManyElectrons, chemist::Density<Electron>>;
+extern template class Coulomb<Electron, DecomposableDensity<Electron>>;
+extern template class Coulomb<ManyElectrons, DecomposableDensity<Electron>>;
 extern template class Coulomb<Electron, Nuclei>;
 extern template class Coulomb<ManyElectrons, Nuclei>;
 extern template class Coulomb<Nuclei, Nuclei>;
diff --git a/include/chemist/quantum_mechanics/operator/detail_/operator_impl.hpp b/include/chemist/quantum_mechanics/operator/detail_/operator_impl.hpp
index b92134933..63d2a346f 100644
--- a/include/chemist/quantum_mechanics/operator/detail_/operator_impl.hpp
+++ b/include/chemist/quantum_mechanics/operator/detail_/operator_impl.hpp
@@ -15,6 +15,7 @@
  */
 
 #pragma once
+#include <chemist/density/density.hpp>
 #include <chemist/dsl/dsl.hpp>
 #include <chemist/electron/electron.hpp>
 #include <chemist/nucleus/nucleus.hpp>
diff --git a/include/chemist/quantum_mechanics/operator/exchange.hpp b/include/chemist/quantum_mechanics/operator/exchange.hpp
index a7b558e28..7e760d224 100644
--- a/include/chemist/quantum_mechanics/operator/exchange.hpp
+++ b/include/chemist/quantum_mechanics/operator/exchange.hpp
@@ -112,10 +112,9 @@ class Exchange
     }
 };
 
-extern template class Exchange<Electron, Electron>;
-extern template class Exchange<ManyElectrons, ManyElectrons>;
-extern template class Exchange<Electron, Nuclei>;
-extern template class Exchange<ManyElectrons, Nuclei>;
-extern template class Exchange<Nuclei, Nuclei>;
+extern template class Exchange<Electron, chemist::Density<Electron>>;
+extern template class Exchange<ManyElectrons, chemist::Density<Electron>>;
+extern template class Exchange<Electron, DecomposableDensity<Electron>>;
+extern template class Exchange<ManyElectrons, DecomposableDensity<Electron>>;
 
 } // namespace chemist::qm_operator
\ No newline at end of file
diff --git a/include/chemist/quantum_mechanics/operator/operator_visitor.hpp b/include/chemist/quantum_mechanics/operator/operator_visitor.hpp
index 4519d707b..739f3d005 100644
--- a/include/chemist/quantum_mechanics/operator/operator_visitor.hpp
+++ b/include/chemist/quantum_mechanics/operator/operator_visitor.hpp
@@ -15,6 +15,7 @@
  */
 
 #pragma once
+#include <chemist/density/density.hpp>
 #include <chemist/electron/electron.hpp>
 #include <chemist/nucleus/nucleus.hpp>
 #include <chemist/quantum_mechanics/operator/operator_fwd.hpp>
@@ -37,6 +38,11 @@ namespace chemist::qm_operator {
     OVERLOADS(T<Electron, Nuclei>);             \
     OVERLOADS(T<ManyElectrons, Nuclei>);        \
     OVERLOADS(T<Nuclei, Nuclei>)
+#define DENSITY_OVERLOADS(T)                                 \
+    OVERLOADS(T<Electron, chemist::Density<Electron>>);      \
+    OVERLOADS(T<ManyElectrons, chemist::Density<Electron>>); \
+    OVERLOADS(T<Electron, DecomposableDensity<Electron>>);   \
+    OVERLOADS(T<ManyElectrons, DecomposableDensity<Electron>>)
 
 /** @brief Base class for all operator visitors.
  *
@@ -78,7 +84,8 @@ class OperatorVisitor {
     OVERLOADS(Density<wavefunction::CMOs, Electron>);
     ONE_PARTICLE_OVERLOADS(Kinetic);
     TWO_PARTICLE_OVERLOADS(Coulomb);
-    TWO_PARTICLE_OVERLOADS(Exchange);
+    DENSITY_OVERLOADS(Coulomb);
+    DENSITY_OVERLOADS(Exchange);
     TWO_PARTICLE_OVERLOADS(ExchangeCorrelation);
 
 private:
@@ -86,6 +93,7 @@ class OperatorVisitor {
     bool m_throw_ = true;
 };
 
+#undef DENSITY_OVERLOADS
 #undef TWO_PARTICLE_OVERLOADS
 #undef ONE_PARTICLE_OVERLOADS
 #undef OVERLOADS
diff --git a/include/chemist/quantum_mechanics/operator/typedefs.hpp b/include/chemist/quantum_mechanics/operator/typedefs.hpp
index 91380bdc8..6e2dc0f97 100644
--- a/include/chemist/quantum_mechanics/operator/typedefs.hpp
+++ b/include/chemist/quantum_mechanics/operator/typedefs.hpp
@@ -15,6 +15,7 @@
  */
 
 #pragma once
+#include <chemist/density/density.hpp>
 #include <chemist/electron/electron.hpp>
 #include <chemist/nucleus/nucleus.hpp>
 #include <chemist/quantum_mechanics/operator/operator_fwd.hpp>
@@ -53,11 +54,11 @@ using v_en_type = Coulomb<Electron, Nuclei>;
 using V_en_type = Coulomb<ManyElectrons, Nuclei>;
 using V_nn_type = Coulomb<Nuclei, Nuclei>;
 
-using k_ee_type = Exchange<Electron, Electron>;
-using K_ee_type = Exchange<ManyElectrons, ManyElectrons>;
-using k_en_type = Exchange<Electron, Nuclei>;
-using K_en_type = Exchange<ManyElectrons, Nuclei>;
-using K_nn_type = Exchange<Nuclei, Nuclei>;
+using j_e_type = Coulomb<Electron, DecomposableDensity<Electron>>;
+using J_e_type = Coulomb<ManyElectrons, DecomposableDensity<Electron>>;
+
+using k_e_type = Exchange<Electron, DecomposableDensity<Electron>>;
+using K_e_type = Exchange<ManyElectrons, DecomposableDensity<Electron>>;
 
 using xc_ee_type = ExchangeCorrelation<Electron, Electron>;
 using XC_ee_type = ExchangeCorrelation<ManyElectrons, ManyElectrons>;
diff --git a/include/chemist/traits/density_traits.hpp b/include/chemist/traits/density_traits.hpp
new file mode 100644
index 000000000..9d15db5be
--- /dev/null
+++ b/include/chemist/traits/density_traits.hpp
@@ -0,0 +1,64 @@
+/*
+ * Copyright 2024 NWChemEx-Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+#include <chemist/traits/chemist_class_traits.hpp>
+
+namespace chemist {
+
+template<typename... ParticleTypes>
+class Density;
+
+template<typename... ParticleTypes>
+class DecomposableDensity;
+
+template<typename... ParticleTypes>
+struct ChemistClassTraits<Density<ParticleTypes...>> {
+    using value_type      = Density<ParticleTypes...>;
+    using reference       = value_type&;
+    using const_reference = const value_type&;
+    using view_type       = value_type;
+    using const_view_type = value_type;
+};
+
+template<typename... ParticleTypes>
+struct ChemistClassTraits<const Density<ParticleTypes...>> {
+    using value_type      = Density<ParticleTypes...>;
+    using reference       = const value_type&;
+    using const_reference = const value_type&;
+    using view_type       = value_type;
+    using const_view_type = value_type;
+};
+
+template<typename... ParticleTypes>
+struct ChemistClassTraits<DecomposableDensity<ParticleTypes...>> {
+    using value_type      = DecomposableDensity<ParticleTypes...>;
+    using reference       = value_type&;
+    using const_reference = const value_type&;
+    using view_type       = value_type;
+    using const_view_type = value_type;
+};
+
+template<typename... ParticleTypes>
+struct ChemistClassTraits<const DecomposableDensity<ParticleTypes...>> {
+    using value_type      = DecomposableDensity<ParticleTypes...>;
+    using reference       = const value_type&;
+    using const_reference = const value_type&;
+    using view_type       = value_type;
+    using const_view_type = value_type;
+};
+
+} // namespace chemist
\ No newline at end of file
diff --git a/src/chemist/density/decomposable_density.cpp b/src/chemist/density/decomposable_density.cpp
new file mode 100644
index 000000000..17dd231f3
--- /dev/null
+++ b/src/chemist/density/decomposable_density.cpp
@@ -0,0 +1,23 @@
+/*
+ * Copyright 2024 NWChemEx-Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <chemist/density/decomposable_density.hpp>
+
+namespace chemist {
+
+template class DecomposableDensity<Electron>;
+
+}
\ No newline at end of file
diff --git a/src/chemist/density/density_class.cpp b/src/chemist/density/density_class.cpp
new file mode 100644
index 000000000..14f760863
--- /dev/null
+++ b/src/chemist/density/density_class.cpp
@@ -0,0 +1,23 @@
+/*
+ * Copyright 2024 NWChemEx-Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <chemist/density/density_class.hpp>
+
+namespace chemist {
+
+template class Density<Electron>;
+
+}
\ No newline at end of file
diff --git a/src/chemist/quantum_mechanics/operator/coulomb.cpp b/src/chemist/quantum_mechanics/operator/coulomb.cpp
index b0159871c..0942435f7 100644
--- a/src/chemist/quantum_mechanics/operator/coulomb.cpp
+++ b/src/chemist/quantum_mechanics/operator/coulomb.cpp
@@ -31,6 +31,10 @@ COULOMB::Coulomb(lhs_value_type p0, rhs_value_type p1) noexcept :
 
 template class Coulomb<Electron, Electron>;
 template class Coulomb<ManyElectrons, ManyElectrons>;
+template class Coulomb<Electron, chemist::Density<Electron>>;
+template class Coulomb<ManyElectrons, chemist::Density<Electron>>;
+template class Coulomb<Electron, DecomposableDensity<Electron>>;
+template class Coulomb<ManyElectrons, DecomposableDensity<Electron>>;
 template class Coulomb<Electron, Nuclei>;
 template class Coulomb<ManyElectrons, Nuclei>;
 template class Coulomb<Nuclei, Nuclei>;
diff --git a/src/chemist/quantum_mechanics/operator/exchange.cpp b/src/chemist/quantum_mechanics/operator/exchange.cpp
index 36a799c13..0eb4a14bf 100644
--- a/src/chemist/quantum_mechanics/operator/exchange.cpp
+++ b/src/chemist/quantum_mechanics/operator/exchange.cpp
@@ -29,10 +29,9 @@ EXCHANGE::Exchange(lhs_value_type p0, rhs_value_type p1) noexcept :
 #undef EXCHANGE
 #undef TPARAMS
 
-template class Exchange<Electron, Electron>;
-template class Exchange<ManyElectrons, ManyElectrons>;
-template class Exchange<Electron, Nuclei>;
-template class Exchange<ManyElectrons, Nuclei>;
-template class Exchange<Nuclei, Nuclei>;
+template class Exchange<Electron, chemist::Density<Electron>>;
+template class Exchange<ManyElectrons, chemist::Density<Electron>>;
+template class Exchange<Electron, DecomposableDensity<Electron>>;
+template class Exchange<ManyElectrons, DecomposableDensity<Electron>>;
 
 } // namespace chemist::qm_operator
\ No newline at end of file
diff --git a/src/chemist/quantum_mechanics/operator/operator_visitor.cpp b/src/chemist/quantum_mechanics/operator/operator_visitor.cpp
index 937935e46..288ea27b1 100644
--- a/src/chemist/quantum_mechanics/operator/operator_visitor.cpp
+++ b/src/chemist/quantum_mechanics/operator/operator_visitor.cpp
@@ -46,14 +46,21 @@ namespace chemist::qm_operator {
     OVERLOADS(T<Electron, Nuclei>);             \
     OVERLOADS(T<ManyElectrons, Nuclei>);        \
     OVERLOADS(T<Nuclei, Nuclei>)
+#define DENSITY_OVERLOADS(T)                                 \
+    OVERLOADS(T<Electron, chemist::Density<Electron>>);      \
+    OVERLOADS(T<ManyElectrons, chemist::Density<Electron>>); \
+    OVERLOADS(T<Electron, DecomposableDensity<Electron>>);   \
+    OVERLOADS(T<ManyElectrons, DecomposableDensity<Electron>>)
 
 OVERLOADS(Density<wavefunction::MOs, Electron>);
 OVERLOADS(Density<wavefunction::CMOs, Electron>);
 ONE_PARTICLE_OVERLOADS(Kinetic);
 TWO_PARTICLE_OVERLOADS(Coulomb);
-TWO_PARTICLE_OVERLOADS(Exchange);
+DENSITY_OVERLOADS(Coulomb);
+DENSITY_OVERLOADS(Exchange);
 TWO_PARTICLE_OVERLOADS(ExchangeCorrelation);
 
+#undef DENSITY_OVERLOADS
 #undef TWO_PARTICLE_OVERLOADS
 #undef ONE_PARTICLE_OVERLOADS
 #undef OVERLOADS
diff --git a/tests/cxx/unit_tests/quantum_mechanics/operator/coulomb.cpp b/tests/cxx/unit_tests/quantum_mechanics/operator/coulomb.cpp
index 783e1fff7..6186b11be 100644
--- a/tests/cxx/unit_tests/quantum_mechanics/operator/coulomb.cpp
+++ b/tests/cxx/unit_tests/quantum_mechanics/operator/coulomb.cpp
@@ -23,6 +23,10 @@ using namespace chemist::qm_operator;
 using types2test =
   std::tuple<std::pair<Electron, Electron>,
              std::pair<ManyElectrons, ManyElectrons>,
+             std::pair<Electron, chemist::Density<Electron>>,
+             std::pair<ManyElectrons, chemist::Density<Electron>>,
+             std::pair<Electron, DecomposableDensity<Electron>>,
+             std::pair<ManyElectrons, DecomposableDensity<Electron>>,
              std::pair<Electron, Nuclei>, std::pair<ManyElectrons, Nuclei>,
              std::pair<Nuclei, Nuclei>>;
 
diff --git a/tests/cxx/unit_tests/quantum_mechanics/operator/exchange.cpp b/tests/cxx/unit_tests/quantum_mechanics/operator/exchange.cpp
index 2577d09da..23ea5a162 100644
--- a/tests/cxx/unit_tests/quantum_mechanics/operator/exchange.cpp
+++ b/tests/cxx/unit_tests/quantum_mechanics/operator/exchange.cpp
@@ -21,10 +21,10 @@ using namespace chemist;
 using namespace chemist::qm_operator;
 
 using types2test =
-  std::tuple<std::pair<Electron, Electron>,
-             std::pair<ManyElectrons, ManyElectrons>,
-             std::pair<Electron, Nuclei>, std::pair<ManyElectrons, Nuclei>,
-             std::pair<Nuclei, Nuclei>>;
+  std::tuple<std::pair<Electron, chemist::Density<Electron>>,
+             std::pair<ManyElectrons, chemist::Density<Electron>>,
+             std::pair<Electron, DecomposableDensity<Electron>>,
+             std::pair<ManyElectrons, DecomposableDensity<Electron>>>;
 
 TEMPLATE_LIST_TEST_CASE("Exchange", "", types2test) {
     auto defaulted_values = test_chemist::defaulted_particles();
diff --git a/tests/cxx/unit_tests/quantum_mechanics/test_qm.hpp b/tests/cxx/unit_tests/quantum_mechanics/test_qm.hpp
index 6bbeb03f0..c27658c23 100644
--- a/tests/cxx/unit_tests/quantum_mechanics/test_qm.hpp
+++ b/tests/cxx/unit_tests/quantum_mechanics/test_qm.hpp
@@ -23,6 +23,7 @@
 #pragma once
 #include "../catch.hpp"
 #include <chemist/basis_set/ao_basis_set.hpp>
+#include <chemist/density/density.hpp>
 #include <chemist/electron/electron.hpp>
 #include <chemist/nucleus/nucleus.hpp>
 #include <chemist/quantum_mechanics/wavefunction/wavefunction.hpp>
@@ -83,7 +84,9 @@ inline auto h2_mos() {
 /// Creates a defaulted instance of each particle type we want to test
 inline auto defaulted_particles() {
     return std::make_tuple(chemist::Electron{}, chemist::ManyElectrons{},
-                           chemist::Nucleus{}, chemist::Nuclei{});
+                           chemist::Nucleus{}, chemist::Nuclei{},
+                           chemist::Density<chemist::Electron>{},
+                           chemist::DecomposableDensity<chemist::Electron>{});
 }
 
 /// Creates a non-default instance of each particle type we want to test
@@ -91,8 +94,12 @@ inline auto non_defaulted_particles() {
     auto r = h2_coords();
     chemist::Nucleus h0("H", 1ul, 1837.15264648179, r[0], r[1], r[2]);
     chemist::Nucleus h1("H", 1ul, 1837.15264648179, r[3], r[4], r[5]);
+    using density_type = chemist::Density<chemist::Electron>;
+    typename density_type::value_type rho{{1.0, 0.0}, {0.0, 1.0}};
+    density_type rho_e(rho, h2_aos());
+    chemist::DecomposableDensity<chemist::Electron> rho_e2(rho, h2_mos());
     return std::make_tuple(chemist::Electron{}, chemist::ManyElectrons{3}, h0,
-                           chemist::Nuclei{h0, h1});
+                           chemist::Nuclei{h0, h1}, rho_e, rho_e2);
 }
 
 } // namespace test_chemist