Remove arbitrary-dimension lattices and systems

docs/src/examples.md

@@ -180,7 +180,7 @@ meas_rate = 40     # Number of timesteps between snapshots of LLD to input to FF
                    # The maximum frequency we resolve is set by 2π/(meas_rate * Δt)
 dyn_meas = 400     # Total number of frequencies we'd like to resolve
 dynsf = dynamic_structure_factor(
-    sys, sampler; therm_samples=10, dynΔt=Δt, meas_rate=meas_rate,
+    sys, sampler; nsamples=10, dynΔt=Δt, meas_rate=meas_rate,
     dyn_meas=dyn_meas, bz_size=(1,1,2), thermalize=10, verbose=true,
     reduce_basis=true, dipole_factor=false,
 )
@@ -325,7 +325,7 @@ meas_rate = convert(Int, div(2π, (2 * target_max_ω * Δt)))
 sampler = MetropolisSampler(system, kT, 500)
 println("Starting structure factor measurement...")
 dynsf = dynamic_structure_factor(
-    system, sampler; therm_samples=15, thermalize=15,
+    system, sampler; nsamples=15, thermalize=15,
     bz_size=(2,0,0), reduce_basis=true, dipole_factor=true,
     dynΔt=Δt, meas_rate=meas_rate, dyn_meas=1000, verbose=true, 
 )
@@ -496,7 +496,6 @@ To discover the symmetry allowed couplings for all bonds up to a certain
 distance, we can use the function [`print_bond_table`](@ref).
 
 ```
-lat_vecs = lattice_vectors(1, 1, 1, 90, 90, 90)
 crystal = Sunny.diamond_crystal()
 print_bond_table(crystal, 1.0)
 ```
@@ -548,7 +547,7 @@ find the ones starting from atom 2 we can use
 ```
 julia> all_symmetry_related_bonds_for_atom(crystal, 2, Bond(2, 3, [0, 0, 0]))
 
-4-element Vector{Bond{3}}:
+4-element Vector{Bond}:
  Bond(2, 7, [0, 0, -1])
  Bond(2, 8, [0, 0, -1])
  Bond(2, 3, [0, 0, 0])
@@ -563,8 +562,8 @@ julia> print_bond(crystal, Bond(1, 6, [1,-1,0]))
 Bond(1, 6, [1, -1, 0])
 Distance 1.225, coordination 24
 Connects [0, 0, 0] to [1, -0.5, 0.5]
-Allowed exchange matrix: |   A  D+E -D-E |
-                         | D-E    B    C |
-                         |-D+E    C    B |
-Allowed DM vector: [0 E E]
+Allowed exchange matrix: |   A  D-E -D+E |
+                         | D+E    B    C |
+                         |-D-E    C    B |
+Allowed DM vector: [0 -E -E]
 ```

docs/src/index.md

@@ -14,10 +14,10 @@ This is the documentation for [Sunny.jl](https://github.com/MagSims/Sunny.jl). T
 - Structure factor calculations
 - Automated symmetry analysis and bond equivalency class discovery
 - Interactive plotting using [Makie.jl](https://github.com/JuliaPlots/Makie.jl)
+- Parallel tempering
 
 **Planned Features**
 - CPU parallelization of dynamics and Monte Carlo simulations
-- Parallel tempering
 - Generalized ``\mathrm{SU(N)}`` models
 - GPU acceleration
 - Accelerated local MC updates in the presence of long-range dipole interactions

examples/PT_afm_diamond.jl

@@ -12,7 +12,7 @@ crystal = Sunny.diamond_crystal()
 # interactions -- units of K  
 J = 28.28
 interactions = [
-    heisenberg(J, Bond{3}(1, 3, [0,0,0])),
+    heisenberg(J, Bond(1, 3, [0,0,0])),
 ]
 
 # spin system -- setup to use K

examples/reproduce_testcases.jl

@@ -32,7 +32,7 @@ function test_diamond_heisenberg_sf()
                        # The maximum frequency we resolve is set by 2π/(meas_rate * Δt)
     dyn_meas = 400     # Total number of frequencies we'd like to resolve
     dynsf = dynamic_structure_factor(
-        sys, sampler; therm_samples=10, dynΔt=Δt, meas_rate=meas_rate,
+        sys, sampler; nsamples=10, dynΔt=Δt, meas_rate=meas_rate,
         dyn_meas=dyn_meas, bz_size=(1,1,2), thermalize=10, verbose=true,
         reduce_basis=true, dipole_factor=false,
     )
@@ -174,7 +174,7 @@ function test_FeI2_MC()
     println("Starting structure factor measurement...")
     dynsf = dynamic_structure_factor(
         system, sampler; bz_size=(2,0,0), thermalize=15,
-        therm_samples=15, dipole_factor=true, dyn_meas=1000,
+        nsamples=15, dipole_factor=true, dyn_meas=1000,
         meas_rate=meas_rate, verbose=true,
     )
     S = dynsf.sfactor
@@ -248,79 +248,20 @@ function plot_ET_data(Ts, Es, Eerrors)
     p
 end
 
-function test_gtensor_sf()
-    # An orthorhombic crystal, to allow z-component of g-tensor to be distinct
-    crystal = Crystal(
-        [1. 0. 0.; 0. 1. 0.; 0. 0. 2.],
-        [[0., 0., 0.]],
-    )
-
-    # A ferromagnetic XXZ model, with Jz > Jx = Jy
-    Jx, Jy, Jz = -0.2, -0.2, -1
-    interactions = [
-        exchange(diagm([Jx, Jy, Jz]), Bond(1, 1, [1, 0, 0]))
-        exchange(diagm([Jx, Jy, Jz]), Bond(1, 1, [0, 0, 1]))
-    ]
-    # Only z-component of spin contributes to magnetic moment -- no signal anywhere
-    # sites_info = [
-    #     SiteInfo(1, 1, diagm([0, 0, 1]))
-    # ]
-    # Only xy-component of spin contributes to magnetic moment -- excitation signal in Sxx/Syy
-    site_info = [
-        SiteInfo(1, 1, diagm([1, 1, 0]))
-    ]
-    sys = SpinSystem(crystal, interactions, (8, 8, 8), site_info)
-    rand!(sys)
-
-    Δt = 0.02
-    kT = 1.
-    α  = 0.1
-    nsteps = 20000
-    sampler = LangevinSampler(sys, kT, α, Δt, nsteps)
-
-    meas_rate = 20     # Number of timesteps between snapshots of LLD to input to FFT
-                       # The maximum frequency we resolve is set by 2π/(meas_rate * Δt)
-    dyn_meas = 800     # Total number of frequencies we'd like to resolve
-    dynsf = dynamic_structure_factor(
-        sys, sampler; therm_samples=10, dynΔt=Δt, meas_rate=meas_rate,
-        dyn_meas=dyn_meas, bz_size=(1,1,2), thermalize=10, verbose=true,
-        reduce_basis=true, dipole_factor=false,
-    )
-
-    # Sxx and Syy are symmetry-equivalent, but Szz is not
-    Sxx = real.(dynsf.sfactor[1, 1, :, :, :, :] .+ dynsf.sfactor[2, 2, :, :, :, :])
-    Szz = real.(dynsf.sfactor[3, 3, :, :, :, :])
-
-    # Calculate the maximum ω present in our FFT. Since the time gap between
-    #  our snapshots is meas_rate * Δt, the maximum frequency we resolve
-    #  is 2π / (meas_rate * Δt)
-    # This is implicitly in the same units as the units you use to define
-    #  the interactions in the Hamiltonian. Here, we defined our interactions
-    #  in K, but we want to see ω in units of meV (to compare to a baseline
-    #  linear spin-wave solution we have).
-    maxω = 2π / (meas_rate * Δt)
-    p = plot_many_cuts_afmdiamond(Sxx, 1.0, 1.; maxω=maxω, chopω=5.0)
-
-    display(p)
-    return Sxx, Szz
-end
-
 #= Binned Statistics Routines =#
 
 """
 Calculates the average and binned standard deviation of a set of data.
 The number of bins used is equal to the length of the preallocated `bins` vector
 passed to the function.
 """
-function binned_statistics(data::AbstractVector{T}; nbins::Int=10)::Tuple{T,T} where {T<:Number}
-
+function binned_statistics(data::AbstractVector{T}; nbins::Int=10)::Tuple{T,T} where {T<:Number}    
     bins = zeros(T, nbins)
     avg, stdev = binned_statistics(data, bins)
     return avg, stdev
 end
 
 function binned_statistics(data::AbstractVector{T}, bins::Vector{T})::Tuple{T,T} where {T<:Number}
-
     N = length(data)
     n = length(bins)
     @assert length(data)%length(bins) == 0

src/Ewald.jl

@@ -17,7 +17,7 @@ Specifically, computes:
       - \frac{π}{2Vη^2}Q^2 - \frac{η}{\sqrt{π}} \sum_{i} q_i^2
 ```
 """
-function ewald_sum_monopole(lattice::Lattice{3}, charges::Array{Float64, 4}; η=1.0, extent=5) :: Float64
+function ewald_sum_monopole(lattice::Lattice, charges::Array{Float64, 4}; η=1.0, extent=5) :: Float64
     extent_idxs = CartesianIndices((-extent:extent, -extent:extent, -extent:extent))
 
     recip = gen_reciprocal(lattice)
@@ -80,7 +80,7 @@ function ewald_sum_monopole(lattice::Lattice{3}, charges::Array{Float64, 4}; η=
     return 0.5 * real_space_sum + 2π / vol * real(recip_space_sum) + tot_charge_term + charge_square_term
 end
 
-function direct_sum_monopole(lattice::Lattice{3}, charges::Array{Float64, 4}; s=0.0, extent=5) :: Float64
+function direct_sum_monopole(lattice::Lattice, charges::Array{Float64, 4}; s=0.0, extent=5) :: Float64
     extent_idxs = CartesianIndices((-extent:extent, -extent:extent, -extent:extent))
 
     # Vectors spanning the axes of the entire system
@@ -128,7 +128,7 @@ end
 Performs ewald summation to calculate the potential energy of a 
 system of dipoles with periodic boundary conditions.
 """
-function ewald_sum_dipole(lattice::Lattice{3}, spins::Array{Vec3, 4}; extent=2, η=1.0) :: Float64
+function ewald_sum_dipole(lattice::Lattice, spins::Array{Vec3, 4}; extent=2, η=1.0) :: Float64
     extent_idxs = CartesianIndices((-extent:extent, -extent:extent, -extent:extent))
 
     # Vectors spanning the axes of the entire system
@@ -204,7 +204,7 @@ function ewald_sum_dipole(lattice::Lattice{3}, spins::Array{Vec3, 4}; extent=2,
 end
 
 "Precompute the dipole interaction matrix, not yet in ± compressed form."
-function precompute_monopole_ewald(lattice::Lattice{3}; extent=10, η=1.0) :: OffsetArray{5, Float64} where {D}
+function precompute_monopole_ewald(lattice::Lattice; extent=10, η=1.0) :: OffsetArray{5, Float64} where {D}
     nb = nbasis(lattice)
     A = zeros(Float64, nb, nb, map(n->2*(n-1)+1, lattice.size)...)
     A = OffsetArray(A, 1:nb, 1:nb, map(n->-(n-1):n-1, lattice.size)...)
@@ -293,7 +293,7 @@ function contract_monopole(charges::Array{Float64, 4}, A::OffsetArray{Float64})
 end
 
 "Precompute the dipole interaction matrix, in ± compressed form."
-function precompute_dipole_ewald(lattice::Lattice{3}; extent=3, η=1.0) :: OffsetArray{Mat3, 5}
+function precompute_dipole_ewald(lattice::Lattice; extent=3, η=1.0) :: OffsetArray{Mat3, 5}
     nb = nbasis(lattice)
     A = zeros(Mat3, nb, nb, lattice.size...)
     A = OffsetArray(A, 1:nb, 1:nb, map(n->0:n-1, lattice.size)...)

src/Hamiltonian.jl

@@ -1,30 +1,12 @@
 # Functions associated with HamiltonianCPU, which maintains the actual internal
 # interaction types and orchestrates energy/field calculations.
 
-
 function validate_and_clean_interactions(ints::Vector{<:Interaction}, crystal::Crystal, latsize::Vector{Int64})
-    D = dimension(crystal)
-
-    # Now that we know dimension D, we can convert every OnSiteQuadratic to
-    # QuadraticInteraction
-    ints = map(ints) do int
-        if isa(int, OnSiteQuadratic)
-            return QuadraticInteraction(int.J, Bond{D}(int.site, int.site, zeros(D)), int.label)
-        else
-            return int
-        end
-    end
-
     # Validate all interactions
     for int in ints
         if isa(int, QuadraticInteraction)
             b = int.bond
 
-            # Verify that the dimension is correct
-            if length(b.n) != D
-                error("Interaction $(repr(MIME("text/plain"), int)) inconsistent with crystal dimension $D.")
-            end
-
             # Verify that both basis sites indexed actually exist
             if !(1 <= b.i <= nbasis(crystal)) || !(1 <= b.j <= nbasis(crystal))
                 error("Provided interaction $(repr(MIME("text/plain"), int)) indexes a non-existent basis site.")
@@ -52,11 +34,6 @@ function validate_and_clean_interactions(ints::Vector{<:Interaction}, crystal::C
                 println("Distance-violating interaction: $int.")
                 error("Interaction wraps system.")
             end
-
-        elseif isa(int, DipoleDipole)
-            if D != 3
-                error("Dipole-dipole interactions require three dimensions.")
-            end
         end
     end
 
@@ -65,24 +42,24 @@ end
 
 
 """
-    HamiltonianCPU{D}
+    HamiltonianCPU
 
 Stores and orchestrates the types that perform the actual implementations
 of all interactions internally.
 """
-struct HamiltonianCPU{D}
+struct HamiltonianCPU
     ext_field   :: Union{Nothing, ExternalFieldCPU}
-    heisenbergs :: Vector{HeisenbergCPU{D}}
-    diag_coups  :: Vector{DiagonalCouplingCPU{D}}
-    gen_coups   :: Vector{GeneralCouplingCPU{D}}
+    heisenbergs :: Vector{HeisenbergCPU}
+    diag_coups  :: Vector{DiagonalCouplingCPU}
+    gen_coups   :: Vector{GeneralCouplingCPU}
     dipole_int  :: Union{Nothing, DipoleRealCPU, DipoleFourierCPU}
     spin_mags   :: Vector{Float64}
 end
 
 """
     HamiltonianCPU(ints::Vector{<:Interaction}, crystal, latsize, sites_info::Vector{SiteInfo})
 
-Construct a `HamiltonianCPU{3}` from a list of interactions, converting
+Construct a `HamiltonianCPU` from a list of interactions, converting
 each of the interactions into the proper backend type specialized
 for the given `crystal` and `latsize`.
 
@@ -92,9 +69,9 @@ function HamiltonianCPU(ints::Vector{<:Interaction}, crystal::Crystal,
                         latsize::Vector{Int64}, sites_info::Vector{SiteInfo};
                         μB=BOHR_MAGNETON::Float64, μ0=VACUUM_PERM::Float64)
     ext_field   = nothing
-    heisenbergs = Vector{HeisenbergCPU{3}}()
-    diag_coups  = Vector{DiagonalCouplingCPU{3}}()
-    gen_coups   = Vector{GeneralCouplingCPU{3}}()
+    heisenbergs = Vector{HeisenbergCPU}()
+    diag_coups  = Vector{DiagonalCouplingCPU}()
+    gen_coups   = Vector{GeneralCouplingCPU}()
     dipole_int  = nothing
     spin_mags   = [site.S for site in sites_info]
 
@@ -128,7 +105,7 @@ function HamiltonianCPU(ints::Vector{<:Interaction}, crystal::Crystal,
         end
     end
 
-    return HamiltonianCPU{3}(
+    return HamiltonianCPU(
         ext_field, heisenbergs, diag_coups, gen_coups, dipole_int, spin_mags
     )
 end
@@ -167,7 +144,7 @@ Note that all `_accum_neggrad!` functions should return _just_ the
 this function. Likewise, all code which utilizes local fields should
 be calling _this_ function, not the `_accum_neggrad!`'s directly.
 """
-function field!(B::Array{Vec3}, spins::Array{Vec3}, ℋ::HamiltonianCPU{D}) where {D}
+function field!(B::Array{Vec3}, spins::Array{Vec3}, ℋ::HamiltonianCPU)
     fill!(B, SA[0.0, 0.0, 0.0])
     if !isnothing(ℋ.ext_field)
         _accum_neggrad!(B, ℋ.ext_field)