From e4c7be43f48be226e2eb5b00a5c4deb93383ce8b Mon Sep 17 00:00:00 2001 From: ddahlbom Date: Thu, 10 Oct 2024 13:10:27 -0400 Subject: [PATCH] Expanded tests for dimer model --- src/EntangledUnits/EntangledSpinWaveTheory.jl | 20 +++++++++++ src/EntangledUnits/EntangledUnits.jl | 7 ++++ src/EntangledUnits/TypesAndAliasing.jl | 2 +- test/test_entangled_units.jl | 34 +++++++++++++++++-- 4 files changed, 59 insertions(+), 4 deletions(-) diff --git a/src/EntangledUnits/EntangledSpinWaveTheory.jl b/src/EntangledUnits/EntangledSpinWaveTheory.jl index 636e5a097..c7d59b72a 100644 --- a/src/EntangledUnits/EntangledSpinWaveTheory.jl +++ b/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 diff --git a/src/EntangledUnits/EntangledUnits.jl b/src/EntangledUnits/EntangledUnits.jl index 2f05db767..7efa0a4bf 100644 --- a/src/EntangledUnits/EntangledUnits.jl +++ b/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}}[] diff --git a/src/EntangledUnits/TypesAndAliasing.jl b/src/EntangledUnits/TypesAndAliasing.jl index ccff2ed9a..b76689112 100644 --- a/src/EntangledUnits/TypesAndAliasing.jl +++ b/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) diff --git a/test/test_entangled_units.jl b/test/test_entangled_units.jl index 762f14c70..0e6b4867a 100644 --- a/test/test_entangled_units.jl +++ b/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