diff --git a/.gitignore b/.gitignore
index 2bd8a7b..e0cef88 100644
--- a/.gitignore
+++ b/.gitignore
@@ -14,3 +14,5 @@ notebooks/weibel_instablity.ipynb
 notebooks/distributions.ipynb
 notebooks/strong_landau_damping.ipynb
 notebooks/quietstart.ipynb
+*.csv
+*.jld2
diff --git a/Project.toml b/Project.toml
index 7bbc3f5..c57f9c1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,20 +8,24 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-julia = "^1"
 DataFrames = "^0"
 Distributions = "^0"
 FFTW = "^1"
 FastGaussQuadrature = "^0"
 Sobol = "^1"
 StaticArrays = "^0"
+julia = "^1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/docs/src/diagnostics.md b/docs/src/diagnostics.md
index 9718a81..6932ee7 100644
--- a/docs/src/diagnostics.md
+++ b/docs/src/diagnostics.md
@@ -5,3 +5,9 @@ Modules = [GEMPIC]
 Pages   = ["diagnostics.jl"]
 ```
 
+## Spin version
+
+```@autodocs
+Modules = [GEMPIC]
+Pages   = ["diagnostics_spin.jl"]
+```
diff --git a/docs/src/distributions.md b/docs/src/distributions.md
index 589391c..8281e6a 100644
--- a/docs/src/distributions.md
+++ b/docs/src/distributions.md
@@ -5,3 +5,10 @@ Modules = [GEMPIC]
 Pages   = ["distributions.jl"]
 ```
 
+## Spin version
+
+```@autodocs
+Modules = [GEMPIC]
+Pages   = ["distributions_spin.jl"]
+```
+
diff --git a/docs/src/hamiltonian_splitting.md b/docs/src/hamiltonian_splitting.md
index dcde298..ae3b65f 100644
--- a/docs/src/hamiltonian_splitting.md
+++ b/docs/src/hamiltonian_splitting.md
@@ -5,3 +5,10 @@ Modules = [GEMPIC]
 Pages   = ["hamiltonian_splitting.jl"]
 ```
 
+## Spin version
+
+```@autodocs
+Modules = [GEMPIC]
+Pages   = ["hamiltonian_splitting_spin.jl"]
+```
+
diff --git a/docs/src/particle_group.md b/docs/src/particle_group.md
index fe49bfc..02bc898 100644
--- a/docs/src/particle_group.md
+++ b/docs/src/particle_group.md
@@ -5,3 +5,9 @@ Modules = [GEMPIC]
 Pages   = ["particle_group.jl"]
 ```
 
+## Spin version
+
+```@autodocs
+Modules = [GEMPIC]
+Pages   = ["particle_group_spin.jl"]
+```
diff --git a/docs/src/particle_mesh_coupling.md b/docs/src/particle_mesh_coupling.md
index 2de18a3..2bf3645 100644
--- a/docs/src/particle_mesh_coupling.md
+++ b/docs/src/particle_mesh_coupling.md
@@ -5,3 +5,10 @@ Modules = [GEMPIC]
 Pages   = ["particle_mesh_coupling.jl"]
 ```
 
+## Spin version
+
+```@autodocs
+Modules = [GEMPIC]
+Pages   = ["particle_mesh_coupling_spin.jl"]
+```
+
diff --git a/docs/src/particle_sampling.md b/docs/src/particle_sampling.md
index a84bba7..ba6b762 100644
--- a/docs/src/particle_sampling.md
+++ b/docs/src/particle_sampling.md
@@ -5,3 +5,10 @@ Modules = [GEMPIC]
 Pages   = ["particle_sampling.jl"]
 ```
 
+## Spin version
+
+```@autodocs
+Modules = [GEMPIC]
+Pages   = ["particle_sampling_spin.jl"]
+```
+
diff --git a/notebooks/strong_landau_damping_spin.jl b/notebooks/strong_landau_damping_spin.jl
index 262b483..0ef0d33 100755
--- a/notebooks/strong_landau_damping_spin.jl
+++ b/notebooks/strong_landau_damping_spin.jl
@@ -118,16 +118,6 @@ function run( steps)
 
    store = zeros(ComplexF64,nx)
    
-   ss11 = Float64[]
-   ss12 = Float64[]
-   ss13 = Float64[]
-   ss21 = Float64[]
-   ss22 = Float64[]
-   ss23 = Float64[]
-   ss31 = Float64[]
-   ss32 = Float64[]
-   ss33 = Float64[]
-
    electric = zeros(ComplexF64, steps, nx)
    
    @showprogress 1 for j = 1:steps # loop over time
@@ -150,21 +140,13 @@ function run( steps)
        fft!(store)
        electric[j,:] .= store 
        
-       push!(ss11, GEMPIC.get_s1(propagator.particle_group, 1)) 
-       push!(ss12, GEMPIC.get_s2(propagator.particle_group, 1))
-       push!(ss13, GEMPIC.get_s3(propagator.particle_group, 1)) 
-
-       push!(ss21, GEMPIC.get_s1(propagator.particle_group, 100)) 
-       push!(ss22, GEMPIC.get_s2(propagator.particle_group, 100)) 
-       push!(ss23, GEMPIC.get_s3(propagator.particle_group, 100)) 
-
-       push!(ss31, GEMPIC.get_s1(propagator.particle_group, 800)) 
-       push!(ss32, GEMPIC.get_s2(propagator.particle_group, 800)) 
-       push!(ss33, GEMPIC.get_s3(propagator.particle_group, 800)) 
+       if step % 1000 == 0 
+           save( "particles", step, particle_group)
+       end
        
    end
 
-   return thdiag.data
+   thdiag.data
 
 end
 
diff --git a/src/landau_damping.jl b/src/landau_damping.jl
index c967840..47e4825 100644
--- a/src/landau_damping.jl
+++ b/src/landau_damping.jl
@@ -10,7 +10,7 @@ struct LandauDamping
 end
 
 """
-    sample!( LandauDamping(α, kx), pg)
+    sample!( d, pg)
 
 Sampling from a probability distribution to initialize
 a Landau damping
diff --git a/src/particle_group.jl b/src/particle_group.jl
index 3dfab4f..e4ae542 100644
--- a/src/particle_group.jl
+++ b/src/particle_group.jl
@@ -1,3 +1,5 @@
+using Printf, FileIO, JLD2
+
 export ParticleGroup
 
 abstract type AbstractParticleGroup end
@@ -164,9 +166,9 @@ end
 
 Set weights of particle @ i
 """
-@generated function set_weights( p :: ParticleGroup{D,V}, i :: Int64, w :: Float64 ) where {D, V}
+function set_weights( p :: ParticleGroup{D,V}, i :: Int64, w :: Float64 ) where {D, V}
 
-    :(p.particle_array[$D+$V+1, i] = w)
+    p.particle_array[D+V+1, i] = w
     
 end
 
@@ -180,3 +182,13 @@ function set_common_weight( p :: AbstractParticleGroup, x :: Float64 )
     p.common_weight = x
     
 end
+
+function save( file, step, p :: ParticleGroup{D,V}) where {D, V}
+
+    datafile = @sprintf("%s-%06d.jld2", file, step)
+
+    FileIO.save(datafile, Dict("x" => p.particle_array[1:D,:],
+                               "v" => p.particle_array[D+1:D+V,:],
+                               "w" => p.particle_array[D+V+1:D+V+1,:]))
+
+end
diff --git a/src/particle_group_spin.jl b/src/particle_group_spin.jl
index a5a4209..b382044 100755
--- a/src/particle_group_spin.jl
+++ b/src/particle_group_spin.jl
@@ -1,3 +1,6 @@
+using JLD2, FileIO, Printf
+
+
 export SpinParticleGroup
 
 """
@@ -234,3 +237,15 @@ Set weights of particle @ i
     :(p.particle_array[$D+$V+$S+1, i] = w)
     
 end
+
+
+function save( file, step, p :: SpinParticleGroup{D,V,S}) where {D, V, S}
+
+    datafile = @sprintf("%s-%06d.jld2", file, step)
+
+    FileIO.save(datafile, Dict("x" => p.particle_array[1:D,:],
+                               "v" => p.particle_array[D+1:D+V,:], 
+                               "s" => p.particle_array[D+V+1:D+V+S,:],
+                               "w" => p.particle_array[D+V+S+1,:]))
+
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 8208bc9..33dd583 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,3 +10,4 @@ include("test_sampling.jl")
 include("test_particle_mesh_coupling_spline_1d.jl")
 include("test_hamiltonian_splitting.jl")
 include("test_hamiltonian_splitting_boris.jl")
+include("test_save_particles.jl")
diff --git a/test/test_save_particles.jl b/test/test_save_particles.jl
new file mode 100644
index 0000000..298940d
--- /dev/null
+++ b/test/test_save_particles.jl
@@ -0,0 +1,53 @@
+@testset "Save particles data into file" begin
+
+n_particles = 100000
+xmin        = 1.0 :: Float64
+xmax        = 4π + 1.0
+Lx          = xmax - xmin  
+nx          = 64
+
+mesh = Mesh( xmax, xmin, nx)
+
+pg = ParticleGroup{1,2}(n_particles, 1.0, 1.0, 1)
+
+params = ( k = [[0.5]],
+           α = [0.01],
+           σ = [[0.1, 2.0]],
+           μ = [[0.0, 0.0]]
+)
+
+df = CosSumGaussian{1,2}(params...)
+
+n_particles = pg.n_particles
+   
+sampler = ParticleSampler{1, 2}( :sobol, true, n_particles )
+
+sample!( pg, sampler, df, mesh )
+
+GEMPIC.save( "test", 1, pg)
+
+@test isfile("test-000001.jld2")
+
+σ, μ = 0.02, 0.0
+kx, α = 1.004355, 0.001
+xmin, xmax = 0, 2π/kx
+domain = [xmin, xmax, xmax - xmin]
+nx = 64 
+n_particles = 1000
+mesh = Mesh( xmin, xmax, nx)
+spline_degree = 3
+   
+df = SpinCosSumGaussian{1,1,3}([[kx]], [α], [[σ]], [[μ]] )
+   
+mass, charge = 1.0, 1.0
+   
+spg = SpinParticleGroup{1,1,3}( n_particles, mass, charge, 1)   
+sampler = SpinParticleSampler{1,1,3}( :sobol, n_particles)
+   
+sample!(spg, sampler, df, mesh)
+   
+GEMPIC.save( "test", 2, spg)
+
+@test isfile("test-000002.jld2")
+
+end