From 26b16d6da7b28b890a7d6993cea6510f3251b707 Mon Sep 17 00:00:00 2001
From: ddahlbom <david.dahlbom@gmail.com>
Date: Tue, 24 Sep 2024 16:12:05 -0400
Subject: [PATCH] Fix external field initialization and updating

---
 src/EntangledUnits/EntangledUnits.jl   |  9 +++++++++
 src/EntangledUnits/TypesAndAliasing.jl | 13 ++++++-------
 2 files changed, 15 insertions(+), 7 deletions(-)

diff --git a/src/EntangledUnits/EntangledUnits.jl b/src/EntangledUnits/EntangledUnits.jl
index 0ebf76b12..2f05db767 100644
--- a/src/EntangledUnits/EntangledUnits.jl
+++ b/src/EntangledUnits/EntangledUnits.jl
@@ -244,6 +244,10 @@ function entangle_system(sys::System{M}, units) where M
     # Construct contracted crystal
     contracted_crystal, contraction_info = contract_crystal(sys.crystal, units)
 
+    # Make sure we have a uniform external field
+    @assert allequal(@view sys.extfield[:,:,:,:]) "`EntangledSystems` requires a uniform applied field." 
+    B = sys.extfield[1,1,1,1]
+
     # Determine Ns for local Hilbert spaces (all must be equal). (TODO: Determine if alternative behavior preferable in mixed case.)
     Ns_unit = Ns_in_units(sys, contraction_info)
     Ns_contracted = map(Ns -> prod(Ns), Ns_unit)
@@ -267,9 +271,14 @@ function entangle_system(sys::System{M}, units) where M
         # Pair interactions that become within-unit interactions
         original_interactions = sys.interactions_union[relevant_sites] 
         for (site, interaction) in zip(relevant_sites, original_interactions)
+            # Onsite anisotropy portion
             onsite_original = interaction.onsite
             unit_index = contraction_info.forward[site][2]
             unit_operator += local_op_to_product_space(onsite_original, unit_index, Ns)
+
+            # Zeeman portion
+            S = [local_op_to_product_space(S, unit_index, Ns) for S in spin_matrices((Ns[unit_index]-1)/2)]
+            unit_operator += Hermitian((sys.gs[1, 1, 1, site] * B)' * S)
         end
 
         # Sort all PairCouplings in couplings that will be within a unit and couplings that will be between units
diff --git a/src/EntangledUnits/TypesAndAliasing.jl b/src/EntangledUnits/TypesAndAliasing.jl
index d2d6a0341..ccff2ed9a 100644
--- a/src/EntangledUnits/TypesAndAliasing.jl
+++ b/src/EntangledUnits/TypesAndAliasing.jl
@@ -69,15 +69,14 @@ function EntangledSystem(sys::System{N}, units) where {N}
     for atom in axes(sys.coherents, 4)
         @assert allequal(@view sys.gs[:,:,:,atom]) "`EntangledSystem` require g-factors be uniform across unit cells" 
     end
-    @assert allequal(@view sys.extfield[:,:,:,:]) "`EntangledSystems` requires a uniform applied field." 
 
     # Generate pair of contracted and uncontracted systems
     (; sys_entangled, contraction_info) = entangle_system(sys, units)
     sys_origin = clone_system(sys)
 
-    # Generate observable field. This observable field has as many entres as the
-    # uncontracted system but contains operators in the local product spaces of
-    # the contracted system. `source_idcs` provides the unit index (of the
+    # Generate observable field. This observable field has as many entries as
+    # the uncontracted system but contains operators in the local product spaces
+    # of the contracted system. `source_idcs` provides the unit index (of the
     # contracted system) in terms of the atom index (of the uncontracted
     # system).
     dipole_operators_origin = all_dipole_observables(sys_origin; apply_g=false) 
@@ -87,7 +86,6 @@ function EntangledSystem(sys::System{N}, units) where {N}
 
     # Coordinate sys_entangled and sys_origin
     set_expected_dipoles_of_entangled_system!(esys)
-    set_field!(esys, sys.extfield[1,1,1,1]) # Note external field checked to be uniform
 
     return esys
 end
@@ -141,12 +139,13 @@ function set_field!(esys::EntangledSystem, B)
         unit = source_idcs[1, 1, 1, atom]
         S = dipole_operators[:, 1, 1, 1, atom]
         ΔB = sys_origin.gs[1, 1, 1, atom]' * (B - B_old) 
-        sys.interactions_union[unit].onsite -= Hermitian(ΔB' * S)
+        sys.interactions_union[unit].onsite += Hermitian(ΔB' * S)
     end
 end
 
+# TODO: Actually, we could give a well-defined meaning to this procedure. Implement this.
 function set_field_at!(::EntangledSystem, _, _)
-    error("`EntangledSystem`s do not support inhomogenous external fields. Use `set_field!(sys, B).")
+    error("`EntangledSystem`s do not support inhomogenous external fields. Use `set_field!(sys, B) to set a uniform field.")
 end