Create and add scripts to Lane-UoM

initiatives/Lane-UoM/Disorder-example.jl

@@ -0,0 +1,81 @@
+using Sunny, GLMakie, LinearAlgebra, Random
+
+latvecs = lattice_vectors(1, 1, 20.00, 90, 90, 120);
+pos = [[0,0,0]]
+cryst = Crystal(latvecs, pos)
+rng = MersenneTwister(1643);
+lenx=30
+leny=30
+ranperp=randn(rng,Float64 , lenx*leny)/6
+ranpar=randn(rng,Float64 , lenx*leny)/6
+gpar=1
+gperp = 1
+infos=[1 => Moment(s=1/2,g =1 ) ]
+dims = (lenx,leny,1)
+sys = System(cryst, infos, :dipole;dims, seed=0)
+J₁=diagm([1,1,1])
+set_exchange!(sys,J₁,Bond(1,1,[-1,0,0]))
+
+
+sys_inhom=to_inhomogeneous(sys)
+for i ∈ 1:lenx*leny
+    sys_inhom.gs[i]= [1.0 0.0 0.0;
+                        0.0 1.0 0.0;
+                        0.0 0.0 1.0+ranperp[i]];
+end
+hsat = 9*1*0.5
+set_field!(sys_inhom, 2.5hsat*[0,0,1.0])
+kT_target =0.0
+Δt = 0.02
+λ = 0.1
+langevin = Langevin(Δt; kT=kT_target,damping= λ)
+
+function anneal!(sys, sampler, kTs,nsweeps)
+    Es = zeros(length(kTs))        # Buffer for saving energy as we proceed
+    for (i, kT) in enumerate(kTs)
+        sampler.kT = kT
+        for j ∈ 1:nsweeps
+            step!(sys, sampler)
+        end                
+        Es[i] = energy(sys)   # Query the energy
+    end
+
+    return Es    # Return the energy values collected during annealing
+end;
+
+kT_upper = 3.0
+kTs = [kT_target + (kT_upper-kT_target) * 0.9^k  for k in 0:100] 
+randomize_spins!(sys_inhom) 
+Es = anneal!(sys_inhom, langevin, kTs,2_000 ) 
+for _ in 1:10_000
+    step!(sys_inhom, langevin)
+end
+for i ∈ 1:100
+    minimize_energy!(sys_inhom;maxiters=10_000)
+end
+
+plot_spins(sys_inhom; color = [S[3] for S ∈ sys_inhom.dipoles])
+print_wrapped_intensities(sys_inhom)
+
+measure = ssf_perp(sys_inhom; )
+swt = SpinWaveTheoryKPM(sys_inhom; measure,tol=0.01)
+
+# Γ—K—M-Γ
+Γ = [0,0,0]
+K = [1/3,1/3,0]
+M = [1/2,0,0]
+
+q_points = [Γ,K,M,Γ]
+
+path = q_space_path(cryst, q_points, 200)
+Emax = 17.5
+σin = 0.025*Emax
+kernel = lorentzian(fwhm=σin)
+
+
+@time begin
+    energies = range(0,Emax,150) 
+    res = intensities(swt, path;energies,kernel)
+    plot_intensities(res)
+end
+

initiatives/Lane-UoM/Okubo-example.jl

@@ -0,0 +1,94 @@
+using Sunny, GLMakie
+
+a, c = 1.0, 10.0
+latvecs = lattice_vectors(a, a, c, 90, 90, 120)
+positions = [[0, 0, 0]]
+cryst = Crystal(latvecs, positions)
+
+J₁ = -1
+J₂ = 0.0
+J₃ = -3J₁
+qmod = 2*acos(0.25*(1+sqrt(1-2*(J₁/J₃))))
+x=2π/(qmod)
+x_rat=rationalize(x;tol=0.001)
+dims = (denominator(x_rat),denominator(x_rat),1)
+spinfos = [1 => Moment(s=1,g=1)]
+sys = System(cryst, spinfos, :dipole;dims, seed=0)
+set_exchange!(sys, J₁, Bond(1, 1, [1,0,0]))
+set_exchange!(sys, J₂, Bond(1, 1, [1,2,0]))
+set_exchange!(sys, J₃, Bond(1, 1, [2,0,0]))
+
+# Parameters from Okubo (close to lattice size 26)
+h = -2.5*J₃
+kT_target = 0.3*J₃
+set_field!(sys, [0, 0, h])
+function anneal!(sys, sampler, kTs,nsweeps)
+    Es = zeros(length(kTs))        # Buffer for saving energy as we proceed
+    for (i, kT) in enumerate(kTs)
+        sampler.kT = kT
+        for j ∈ 1:nsweeps
+            step!(sys, sampler)
+        end                
+        Es[i] = energy(sys)   # Query the energy
+    end
+    return Es    # Return the energy values collected during annealing
+end;
+Δt = 0.02
+λ = 0.1
+
+randomize_spins!(sys)
+kT_upper =10*kT_target  # some large temperature to start with (in meV)
+langevin = Langevin(Δt; kT=kT_upper, λ) # initialize the Langevin integrator
+kTs = [kT_target + (kT_upper-kT_target) * 0.9^k  for k in 0:100] # define a temperature schedule
+Es = anneal!(sys, langevin, kTs,2_000 ) # run the anneal
+langevin.kT = kT_target
+
+for _ in 1:2_000
+    step!(sys, langevin)
+end
+
+minimize_energy!(sys;maxiters=10_000) 
+# you may need to do the simulated anneal more carefully to find the true minimum.
+
+########################################
+# Plot scalar spin chirality 
+function spin_chirality(s₁, s₂, s₃)
+    s₁ ⋅ (s₂ × s₃)
+end
+include(joinpath(pkgdir(Sunny), "examples", "extra", "Plotting","plotting2d.jl"))
+
+begin 
+    fig = Figure()
+    ax = Axis(fig[1,1];)
+    crange = (-0.5,0.5)
+    fig=plot_triangular_plaquettes(spin_chirality, [sys.dipoles];resolution = (1600, 800),fontsize = 48, colorrange = crange  )
+    Colorbar(fig[1,2];label = "Spin chirality", colormap = :RdBu , colorrange = crange);
+    fig
+end
+
+
+plot_spins(sys;color = [s[3] for s ∈ sys.dipoles] )
+energy_per_site(sys) 
+
+
+measure = ssf_perp(sys; )
+swt = SpinWaveTheoryKPM(sys; measure,tol=0.01)
+
+
+# Γ—K—M-Γ
+Γ = [0,0,0]
+K = [1/3,1/3,0]
+M = [1/2,0,0]
+
+q_points = [Γ,K,M,Γ]
+path = q_space_path(cryst, q_points, 150)
+Emax = 17.5
+σin = 0.025*Emax
+kernel = lorentzian(fwhm=σin)
+
+@time begin
+    energies = range(0,Emax,150) 
+    res = intensities(swt, path;energies,kernel)
+    plot_intensities(res)
+end
+

initiatives/Lane-UoM/incommensurate-example.jl

@@ -0,0 +1,82 @@
+using Sunny, GLMakie, LinearAlgebra 
+
+a =  10.6172 
+b =  9.0777
+c = 2.9661     
+latvecs = lattice_vectors(a, b, c, 90, 90, 90) 
+positions = [[0.6594, 0.7597,1/4], [0.6127, 0.4401, 1/4], [0.1005, 0.4165,1/4],
+[0.1603, 0.2027,1/4],[0.4752,0.1172,1/4],[0.7852, 0.5270,1/4],[0.4273, 0.4174,1/4]] 
+types = ["Ca", "Cr", "Cr","O","O","O","O"]
+CaCrO = Crystal(latvecs, positions,62; types, setting = "cab")
+cryst = subcrystal(CaCrO,"Cr")
+dims = (1,1,37)
+spinfos = [1 => Moment(s=3/2, g=2),5 => Moment(s=3/2, g=2)]
+sys = System(cryst,spinfos, :dipole;dims)
+Dyanis=0.2
+Day=0.5
+DMI = dmvec([0,Day,0])
+J2=I(3)*8
+J1 = I(3)*1.
+J22 = J2+DMI 
+J21 = J2+DMI
+J11 = J1
+J12 = J1
+Jb = 0.0
+Ja = -4*J1
+
+set_exchange!(sys, J22, Bond(1, 1, [0, 0, 1]))
+set_exchange!(sys, J21, Bond(5, 5, [0, 0, 1]))
+set_exchange!(sys, J11, Bond(6, 7, [0, 1, 0]))
+set_exchange!(sys, J12, Bond(1, 4, [0, 0, 0]))
+set_exchange!(sys, Jb, Bond(3, 5, [0, 0, 0]))
+set_exchange!(sys, Ja,Bond(4, 5, [0, 0, 0]))
+set_onsite_coupling!(sys, S -> Dyanis*S[2]^2, 1)
+set_onsite_coupling!(sys, S -> Dyanis*S[2]^2, 5)
+randomize_spins!(sys)
+minimize_energy!(sys)
+print_wrapped_intensities(sys)
+plot_spins(sys; color =[ S[3] for S in sys.dipoles ])
+print_wrapped_intensities(sys)
+
+# compare to analytical 
+k=2acos(-norm(Ja)/(4norm(diag(J2))))
+krlu = k/2π
+1-krlu
+rationalize(krlu;tol=0.001)
+k=2acos(-norm(Ja)/(4norm(diag(J2))))
+krlu = k/2π
+1-krlu
+kexp=rationalize(1-krlu;tol=0.001) # k_exp in 1BZ
+rationalize(krlu;tol=0.001) 
+kexp
+round(Float64(kexp),digits=4)  === 0.4595
+
+measure = ssf_perp(sys; )
+swt = SpinWaveTheoryKPM(sys; measure,tol=0.01)
+
+Γ = [0,0,0]
+R = [1/2,1/2,1/2]
+S = [1/2,1/2,0]
+Tp = [0,1/2,1/2]
+U = [1/2,0,1/2]
+X = [1/2,0,0]
+Y = [0,1/2,0]
+Z = [0,0,1/2]
+offset = [0.0,0.0,1.0]
+nqs = 100 
+
+q_points = [Γ,Z,U,X,Γ]
+q_points_shift = [point + offset for point in q_points]
+path = q_space_path(cryst, q_points_shift, 150)
+function br(ω, x, σ)
+    return (1/π) * (σ / ((x - ω)^2 + σ^2))
+end
+Emax = 35
+σin= 0.025*Emax
+kernel = lorentzian(fwhm=σin)
+
+@time begin
+    energies = range(0,Emax,150) 
+    res = intensities(swt, path; energies, kernel)
+    plot_intensities(res)
+end