Correct sampling

Project.toml

@@ -14,6 +14,15 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
+[weakdeps]
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+
+[extensions]
+NormalModesMakieExt = "GLMakie"
+
+[compat]
+GLMakie = "0.9"
+
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 

ext/NormalModesMakieExt.jl

@@ -0,0 +1,12 @@
+module NormalModesMakieExt
+
+using NormalModes
+using GLMakie
+
+import NormalModes: animate
+
+function animate(nm::NormalDecomposition)
+
+end
+
+end

src/NormalModes.jl

@@ -13,9 +13,13 @@ using UnitfulAtomic
 export NormalDecomposition
 export project, project_per_atom
 export normal_modes, normal_mode, frequencies, wave_number, reduced_masses
-export normal_from_positions
+export spatial_variances, momentum_variances, atom_spatial_variances
 export sample
 
+export animate
+
+const hbar = 1  # Atomic units
+
 function mass_weight_matrix(elements)
     masses = to_atomic_masses(elements)
 
@@ -126,6 +130,7 @@ Return the real space normal modes.
 """
 function normal_modes(nm::NormalDecomposition)
     modes = nm.M * nm.U
+    return modes
     μs = norm.(eachcol(modes))
     return modes ./ reshape(μs, 1, :)
 end
@@ -167,6 +172,14 @@ function wave_number(nm::NormalDecomposition)
     end
 end
 
+spatial_variances(nm::NormalDecomposition) = 1/2 * (hbar ./ nm.ωs)
+momentum_variances(nm::NormalDecomposition) = 1/2 * hbar * nm.ωs
+
+function atom_spatial_variances(nm::NormalDecomposition)
+    X = nm.M * nm.U
+    return abs.(X .* spatial_variances(nm)')
+end
+
 """
     sample(nm::NormalDecomposition[, n_samples])
 
@@ -176,10 +189,9 @@ Return the deviation from the average geometry and the deviation from zero
 momentum.
 """
 function StatsBase.sample(rng::AbstractRNG, nm::NormalDecomposition, n_samples)
-    hbar = 1  # Atomic units
     X = nm.M * nm.U
-    Δx_dist = MvNormal(Diagonal(1/2 * (hbar ./ nm.ωs)))
-    Δp_dist = MvNormal(Diagonal(1/2 * hbar * nm.ωs))
+    Δx_dist = MvNormal(Diagonal(spatial_variances(nm)))
+    Δp_dist = MvNormal(Diagonal(momentum_variances(nm)))
 
     Δx = X * rand(rng, Δx_dist, n_samples)
     Δp = inv(nm.M)^2 * X * rand(rng, Δp_dist, n_samples)
@@ -195,4 +207,6 @@ end
 StatsBase.sample(nm::NormalDecomposition, n_samples) = sample(Random.GLOBAL_RNG, nm, n_samples)
 StatsBase.sample(nm::NormalDecomposition) = sample(Random.GLOBAL_RNG, nm)
 
-end
+function animate() end
+
+end