Fixed wrong size of X,Y,Z ranges for cube file (#19)

The axes in the cube format are the ones for each voxel, not the
axes for the whole volume.

src/volume/cube.jl

@@ -27,18 +27,18 @@ function read_cube(filename::AbstractString)
     # to Å unit
     origin .*= Bohr
 
-    # number of voxels per spanning vector
+    # number of voxels in domain
     n_voxels = zeros(Int, 3)
-    span_vectors = zeros(Float64, 3, 3)
+    voxel_vectors = zeros(Float64, 3, 3)
     for i in 1:3
         line = split(strip(readline(io)))
         n_v = parse(Int, line[1])
         n_voxels[i] = n_v
         # bohr unit
-        span_vectors[:, i] = parse.(Float64, line[2:4])
+        voxel_vectors[:, i] = parse.(Float64, line[2:4])
     end
     # to Å unit
-    span_vectors .*= Bohr
+    voxel_vectors .*= Bohr
 
     atom_positions = zeros(Float64, 3, n_atoms)
     atom_numbers = zeros(Int, n_atoms)
@@ -53,10 +53,10 @@ function read_cube(filename::AbstractString)
     atom_positions .*= Bohr
 
     n_x, n_y, n_z = n_voxels
-    # fractional w.r.t. span_vectors
-    X = range(0, 1, n_x)
-    Y = range(0, 1, n_y)
-    Z = range(0, 1, n_z)
+    # fractional w.r.t. voxel_vectors
+    X = range(0, n_x-1, n_x)
+    Y = range(0, n_y-1, n_y)
+    Z = range(0, n_z-1, n_z)
 
     W = zeros(Float64, n_x, n_y, n_z)
     # 6 columns per line
@@ -84,7 +84,7 @@ function read_cube(filename::AbstractString)
     end
 
     close(io)
-    return (; atom_positions, atom_numbers, origin, span_vectors, X, Y, Z, W)
+    return (; atom_positions, atom_numbers, origin, voxel_vectors, X, Y, Z, W)
 end
 
 """
@@ -96,20 +96,20 @@ Write `cube` file.
 - `atom_positions`: `3 * n_atoms`, Å, cartesian coordinates
 - `atom_numbers`: `n_atoms`, atomic numbers
 - `origin`: `3`, Å, origin of the grid
-- `span_vectors`: `3 * 3`, Å, each column is a spanning vector
+- `voxel_vectors`: `3 * 3`, Å, each column is a voxel vector
 - `W`: `nx * ny * nz`, volumetric data
 """
 function write_cube(
     filename::AbstractString,
     atom_positions::AbstractMatrix{T},
     atom_numbers::AbstractVector{Int},
     origin::AbstractVector{T},
-    span_vectors::AbstractMatrix{T},
+    voxel_vectors::AbstractMatrix{T},
     W::AbstractArray{T,3},
 ) where {T<:Real}
     n_atoms = length(atom_numbers)
     size(atom_positions, 2) == n_atoms || error("incompatible n_atoms")
-    size(span_vectors) == (3, 3) || error("incompatible span_vectors")
+    size(voxel_vectors) == (3, 3) || error("incompatible voxel_vectors")
     length(origin) == 3 || error("origin must be 3-vector")
 
     @info "Writing cube file: " filename
@@ -127,7 +127,7 @@ function write_cube(
     for i in 1:3
         # number of voxels
         n_v = n_xyz[i]
-        ax = span_vectors[:, i] ./ (n_v * Bohr)
+        ax = voxel_vectors[:, i] ./ Bohr
         @printf(io, "%d %12.6f %12.6f %12.6f\n", n_v, ax...)
     end
 

test/volume/cube.jl

@@ -0,0 +1,27 @@
+@testitem "read/write cube" begin
+    using LazyArtifacts
+    using WannierIO: Bohr
+    cube = read_cube(artifact"Si2_valence/reference/Si2_valence_00001.cube")
+    @test cube.origin ≈ [-7.10461, -9.47281, -9.47281]*Bohr 
+    @test cube.voxel_vectors ≈
+        transpose([0.0 0.39470 0.39470; 0.39470 0.00000 0.39470; 0.39470 0.39470 0.00000])*Bohr 
+    @test cube.X ≈ range(0, 19, 20)
+    @test cube.Y ≈ cube.X
+    @test cube.Z ≈ cube.X
+    @test size(cube.W) == (20,20,20)
+    @test size(cube.atom_positions) == (3,8)
+    @test cube.atom_positions[:,1] ≈ transpose([0.00000 -5.13111 -5.13111])*Bohr
+
+    tmpfile = tempname(; cleanup=true)
+    write_cube(tmpfile, cube.atom_positions, cube.atom_numbers, cube.origin, cube.voxel_vectors, cube.W)
+    cube2 = read_cube(tmpfile)
+
+    @test cube.atom_positions ≈ cube2.atom_positions
+    @test cube.atom_numbers ≈ cube2.atom_numbers
+    @test cube.origin ≈ cube2.origin
+    @test cube.voxel_vectors ≈ cube2.voxel_vectors
+    @test cube.X ≈ cube2.X
+    @test cube.Y ≈ cube2.Y
+    @test cube.Z ≈ cube2.Z
+    @test cube.W ≈ cube2.W
+end