diff --git a/src/EntangledUnits/TypesAndAliasing.jl b/src/EntangledUnits/TypesAndAliasing.jl
index b76689112..dc41f9221 100644
--- a/src/EntangledUnits/TypesAndAliasing.jl
+++ b/src/EntangledUnits/TypesAndAliasing.jl
@@ -149,7 +149,29 @@ function set_field_at!(::EntangledSystem, _, _)
 end
 
 
-set_dipole!(esys::EntangledSystem, dipole, site; kwargs...) = error("Setting dipoles of an EntangledSystem not well defined.") 
+function set_dipole!(::EntangledSystem, dipole, site)
+    error("`set_dipole!` operation for `EntangledSystem` not well defined. Consider using `set_coherent!` to set the state of each entangled unit.")
+end
+
+# Find the unique coherent state corresponding to a set of fully-polarized
+# dipoles on each site inside a specified entangled unit.
+function coherent_state_from_dipoles(esys::EntangledSystem, dipoles, unit)
+    (; sys_origin, contraction_info) = esys
+
+    original_sites = [id.site for id in contraction_info.inverse[unit]]
+    Ns = Ns_in_units(sys_origin, contraction_info)[unit]
+
+    @assert length(dipoles) == length(original_sites) == length(Ns) "Invalid number of dipoles for specified unit."
+
+    coherents = []
+    for (dipole, N) in zip(dipoles, Ns)
+        S = spin_matrices((N-1)/2)
+        coherent = eigvecs(S' * dipole)[:,N] # Retrieve highest-weight eigenvector
+        push!(coherents, coherent)
+    end
+
+    return kron(coherents...)
+end
 
 # Sets the coherent state of a specified unit. The `site` refers to the
 # contracted lattice (i.e., to a "unit"). The function then updates all dipoles
diff --git a/test/test_entangled_units.jl b/test/test_entangled_units.jl
index ea89ca637..6b8837c96 100644
--- a/test/test_entangled_units.jl
+++ b/test/test_entangled_units.jl
@@ -55,6 +55,13 @@ end
     field_offset = 2*(Sl[3] + Su[3]) # 2 for g-factor
     @test onsite_operator ≈ onsite_ref + field_offset 
 
+    # Test apparatus for setting coherent states from dipoles specification
+    dipoles = [[0, 1/2, 0], [0, -1/2, 0]] # Dipoles specifying a dimer state
+    cs = Sunny.coherent_state_from_dipoles(esys, dipoles, 1)
+    set_coherent!(esys, cs, CartesianIndex(1,1,1,1))
+    @test esys.sys_origin.dipoles[1,1,1,1][2] ≈ 1/2
+    @test esys.sys_origin.dipoles[1,1,1,2][2] ≈ -1/2
+
     # Test external field works in action
     set_field!(esys, [0, 0, 10])
     randomize_spins!(esys)