Expanded tests for dimer model

src/EntangledUnits/EntangledSpinWaveTheory.jl

@@ -255,6 +255,26 @@ function excitations!(T, tmp, swt::EntangledSpinWaveTheory, q)
     end
 end
 
+# No changes
+function excitations(swt::EntangledSpinWaveTheory, q)
+    L = nbands(swt)
+    T = zeros(ComplexF64, 2L, 2L)
+    H = zeros(ComplexF64, 2L, 2L)
+    energies = excitations!(T, copy(H), swt, q)
+    return (energies, T)
+end
+
+# No changes
+function dispersion(swt::EntangledSpinWaveTheory, qpts)
+    L = nbands(swt)
+    qpts = convert(AbstractQPoints, qpts)
+    disp = zeros(L, length(qpts.qs))
+    for (iq, q) in enumerate(qpts.qs)
+        view(disp, :, iq) .= view(excitations(swt, q)[1], 1:L)
+    end
+    return reshape(disp, L, size(qpts.qs)...)
+end
+
 # Only change is no Ewald
 function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::EntangledSpinWaveTheory, q_reshaped::Vec3)
     (; sys) = swt

src/EntangledUnits/EntangledUnits.jl

@@ -258,6 +258,13 @@ function entangle_system(sys::System{M}, units) where M
     spin_infos = [i => Moment(s=(N-1)/2, g=1.0) for (i, N) in enumerate(Ns_contracted)]  # TODO: Decisions about g-factor 
     sys_entangled = System(contracted_crystal, spin_infos, :SUN; dims)
 
+    # Transfer rng from origin system to entangled system
+    sys_entangled.rng.s0 = sys.rng.s0
+    sys_entangled.rng.s1 = sys.rng.s1
+    sys_entangled.rng.s2 = sys.rng.s2
+    sys_entangled.rng.s3 = sys.rng.s3
+    sys_entangled.rng.s4 = sys.rng.s4
+
     # TODO: Extend to inhomogenous systems
     # For each contracted site, scan original interactions and reconstruct as necessary.
     new_pair_data = Tuple{Bond, Matrix{ComplexF64}}[]

src/EntangledUnits/TypesAndAliasing.jl

@@ -80,7 +80,7 @@ function EntangledSystem(sys::System{N}, units) where {N}
     # contracted system) in terms of the atom index (of the uncontracted
     # system).
     dipole_operators_origin = all_dipole_observables(sys_origin; apply_g=false) 
-    (; observables, source_idcs)  = observables_to_product_space(dipole_operators_origin, sys_origin, contraction_info)
+    (; observables, source_idcs) = observables_to_product_space(dipole_operators_origin, sys_origin, contraction_info)
 
     esys = EntangledSystem(sys_entangled, sys_origin, contraction_info, observables, source_idcs)
 

test/test_entangled_units.jl

@@ -33,7 +33,7 @@ end
     positions = [[0, 0, 0], [0.0, 0.5, 0.0]] 
 
     crystal = Crystal(latvecs, positions, 1; types = ["A", "B"])
-    sys = System(crystal, [1 => Moment(s=1/2, g=2), 2 => Moment(s=1/2, g=2)], :SUN)
+    sys = System(crystal, [1 => Moment(s=1/2, g=2), 2 => Moment(s=1/2, g=2)], :SUN; seed=1)
     set_exchange!(sys, J, Bond(1, 2, [0, 0, 0]))
     set_exchange!(sys, J′, Bond(1, 1, [1, 0, 0]))
     set_exchange!(sys, J′, Bond(2, 2, [1, 0, 0]))  # Needed because we broke the symmetry equivalence of the two sites
@@ -65,11 +65,39 @@ end
     bond_ref = J′*((Sl2' * Sl1) .+ (Su2' * Su1))
     @test bond_operator ≈ bond_ref
 
-    # Test dispersion against analytical formula
+
+    # Test dispersion against analytical formula for antisymmetric channel.
+    qpts = q_space_path(crystal, [[0, 1, 0], [1, 1, 0]], 5)
+    qs = qpts.qs
+
     ω_ref(q, J, J′) = J*sqrt(1 + 2(J′/J) * cos(2π*q))
+    ωs_analytical = ω_ref.([q[1] for q in qs], J, J′)
+
+    for unit in Sunny.eachunit(esys)
+        set_coherent!(esys, [0, 1/√2, -1/√2, 0], unit)
+    end
+    swt = SpinWaveTheory(esys; measure=Sunny.empty_measurespec(sys), regularization=0.0)
+    disp = dispersion(swt, qptss)
+    ωs_numerical = disp[1,:]
 
-    qs = [[0, 1.0, 0], []]
+    @test all(both -> isapprox(both[1], both[2]; atol=1e-12), zip(ωs_analytical, ωs_numerical))
 
+
+    # Test classical dynamics and perform golden test.
+    esys = repeat_periodically(esys, (8, 1, 1))
+    energies = range(0, 2, 10)
+    dt = 0.1
+    measure = ssf_trace(esys)
+    sc = SampledCorrelations(esys; dt, energies, measure)
+    integrator = Langevin(dt ; damping=0.4, kT=0.05)
+
+    for _ in 1:100
+        step!(esys, integrator)
+    end
+    add_sample!(sc, esys)
+    res = intensities(sc, qpts; energies, kT=0.05)
+    intensities_ref = [-0.0021805338918195185 -0.00323580926131731 -0.003475460861305824 -0.003235809261317303 -0.0021805338918195185; 0.015394114132126614 0.05504458917464571 0.006363628894456872 0.024858210899731406 0.015394114132126614; -0.062439105112575104 -0.24430343196998291 -0.1124375404098702 -0.02581168756492061 -0.062439105112575104; 0.392194235986946 2.8241763179059833 4.85842932490045 0.4259092656302079 0.392194235986946; 9.761448823359627 18.38878057796522 15.41704091322518 2.3935803890294105 9.761448823359627; 19.84319371801515 22.70779362388382 11.689296089047344 2.883234285853617 19.84319371801515; 9.747007990151367 5.296141862559427 0.8103824155070479 0.6877505519539169 9.747007990151367; -0.42870193438878 -0.9643857517663812 -0.0903901605398698 -0.22212424683216814 -0.42870193438878; 0.14524976514187082 0.2766893084879453 0.20093635005032084 0.032539493064106036 0.14524976514187082; -0.11486352541237474 -0.2600088555205098 -0.0036428642723143265 -0.26000885552050884 -0.11486352541237474]
+    @test isapprox(res.data, intensities_ref)
 end
 
 # @testitem "Ba3Mn2O8 Dispersion and Golden Test" begin