Some refactoring (#16)

.gitignore

@@ -1,3 +1,5 @@
 /Manifest.toml
 /docs/Manifest.toml
 /docs/build/
+*~
+.*.swp

src/AtomsBuilder.jl

@@ -3,13 +3,11 @@ module AtomsBuilder
 import AtomsBase
 using StaticArrays: SMatrix, SVector
 
-const Mat3{T} = SMatrix{3, 3, T} 
+const Mat3{T} = SMatrix{3, 3, T}
 const Vec3{T} = SVector{3, T}
 
 include("chemistry.jl")
-
 include("utils.jl")
-
 include("bulk.jl")
 
 end

src/bulk.jl

@@ -1,46 +1,45 @@
 export bulk
 
-using Unitful: unit 
+using Unitful: unit
 using LinearAlgebra: I
 
 const _unit_cells = Dict(     # (positions, cell matrix, factor of a)
-   :fcc => ( [ [0 0 0], ],
-             [0 1 1; 1 0 1; 1 1 0],  0.5),
+   :fcc => ( [ [0 0 0], ], [0 1 1; 1 0 1; 1 1 0],  0.5),
    :bcc => ( [ [0.0,0.0,0.0], ],
              [-1 1 1; 1 -1 1; 1 1 -1], 0.5),
    :diamond => ( [ [0.0, 0.0, 0.0], [0.5, 0.5, 0.5] ],
                  [0 1 1; 1 0 1; 1 1 0], 0.5)
 )
 
 const _cubic_cells = Dict(   # (positions, factor of a)
-   :fcc => ( [ [0 0 0], [0 1 1], [1 0 1], [1 1 0] ], 0.5 ),
-   :bcc => ( [ [0 0 0], [1 1 1] ], 0.5 ),
+   :fcc     => ( [ [0 0 0], [0 1 1], [1 0 1], [1 1 0] ], 0.5 ),
+   :bcc     => ( [ [0 0 0], [1 1 1] ], 0.5 ),
    :diamond => ( [ [0 0 0], [1 1 1], [0 2 2], [1 3 3], [2 0 2],
                    [3 1 3], [2 2 0], [3 3 1] ], 0.25 )
 )
 
-const _simple_structures = [:fcc, :bcc, :diamond]
-
-function _simple_bulk(sym::Symbol, cubic::Bool; a=nothing, T = Float64)
+function _simple_bulk(sym::Symbol, cubic::Bool; a=nothing, T=Float64)
    if cubic
       X, scal = _cubic_cells[Chemistry.symmetry(sym)]
       C = Matrix(1.0I, 3, 3) / scal
    else
       X, C, scal = _unit_cells[Chemistry.symmetry(sym)]
    end
-   a === nothing && (a = Chemistry.lattice_constant(sym))
+   a = @something a Chemistry.lattice_constant(sym)
    TU = typeof( one(T) * unit(a) )
-   return [ Vec3{TU}(x * a * scal)  for x in X ], Mat3{TU}(C * a * scal)
+   [ Vec3{TU}(x * a * scal)  for x in X ], Mat3{TU}(C * a * scal)
 end
 
 
 function _bulk_hcp(sym::Symbol; a=nothing, c=nothing, T=Float64)
    D = Chemistry.refstate(sym)
-   a === nothing && (a = D["a"])
-   c === nothing && (c = a * D["c/a"])
-   return [ a * Vec3{T}(0.0, 0.0, 0.0), 
-            a * Vec3{T}(0.0, 1 / sqrt(3), c / a / 2) ],
-         a * Mat3{T}( [1.0, -1/2,  0.0, 0.0, sqrt(3)/2, 0.0, 0.0, 0.0, c/a] )
+   a = @something a   D["a"]
+   c = @something c a*D["c/a"]
+   [a * Vec3{T}(0.0,         0.0,       0.0),
+    a * Vec3{T}(0.0, 1 / sqrt(3), c / a / 2) ],
+    a * Mat3{T}([1,      -1/2,   0,
+                 0, sqrt(3)/2,   0,
+                 0,         0, c/a])
 end
 
 
@@ -49,36 +48,37 @@ _convert_pbc(pbc::Bool) = (pbc, pbc, pbc)
 _convert_pbc(pbc::AbstractVector) = tuple(pbc...)
 
 """
-`bulk(sym)` : generates a `FlexibleSystem` unit cell for a bulk 
-crystal structure. If `sym` is a chemical symbol then the phase and  
-lattice constant are taken from a database that is consistent with ASE. 
-If `sym` is one of `[:fcc, :bcc, :diamond, :hcp]` then one needs to 
-specify the kwargs `a` or `c` to determine the lattice constants. 
+`bulk(sym)` : generates a `FlexibleSystem` unit cell for a bulk
+crystal structure. If `sym` is a chemical symbol then the phase and
+lattice constant are taken from a database that is consistent with ASE.
+Optional alternative values can be chosen via the kwargs
+`a`, `b` or `c` to specify the lattice constants.
 """
-function bulk(sym::Symbol; T=Float64, cubic = false, pbc = (true,true,true), 
-                           a=nothing, c=nothing) # , x=nothing, y=nothing, z=nothing)
-   symm = try
-      Chemistry.symmetry(sym)
-   catch e
-      if e isa KeyError
-         throw(ArgumentError("No symmetry information known for element $sym"))
-      else
-         rethrow()
-      end
-   end
-   if symm in _simple_structures
-      X, C = _simple_bulk(sym, cubic; a=a)
-   elseif symm == :hcp
-      X, C = _bulk_hcp(sym; a=a, c=c)  # cubic parameter is irrelevant for hcp
-   else
-      throw(ArgumentError("Currently bulk not immplemented for symmetry $symm"))
-   end
-   m = Chemistry.atomic_mass(sym)
-   Z = Chemistry.atomic_number(sym)
-   nat = length(X)
-   at = _flexible_system( X, fill(Z, nat), C, _convert_pbc(pbc))
-   # I don't remember what this does so I'm commenting it out until I understand 
-   # it again ...                         
-   # (x !== nothing || y !== nothing || z !== nothing) && rotate!(at, x=x, y=y, z=z)
-   return at
+function bulk(sym::Symbol; cubic=false, pbc=(true, true, true),
+                           a=nothing, c=nothing, b=nothing, T=Float64)
+                           # , x=nothing, y=nothing, z=nothing)
+    symm = try
+        Chemistry.symmetry(sym)
+    catch e
+        if e isa KeyError
+            throw(ArgumentError("No symmetry information known for element $sym"))
+        else
+            rethrow()
+        end
+    end
+    if symm in (:fcc, :bcc, :diamond)
+        X, C = _simple_bulk(sym, cubic; a, T)
+    elseif symm == :hcp
+        cubic && throw(ArgumentError("cubic=true cannot be selected for $symm lattices"))
+        X, C = _bulk_hcp(sym; a, c, T)
+    else
+        throw(ArgumentError("Currently bulk not immplemented for symmetry $symm"))
+    end
+    Z = Chemistry.atomic_number(sym)
+    nat = length(X)
+    at = _flexible_system( X, fill(Z, nat), C, _convert_pbc(pbc))
+    # I don't remember what this does so I'm commenting it out until I understand 
+    # it again ...                         
+    # (x !== nothing || y !== nothing || z !== nothing) && rotate!(at, x=x, y=y, z=z)
+    return at
 end

src/chemistry.jl

@@ -1,9 +1,9 @@
 module Chemistry
 
-using JSON, Unitful 
+using JSON, Unitful
 
 # --------------------------------------------------------------------
-# Load data and prepare it a bit ... 
+# Load data and prepare it a bit ...
 
 const ase_data_path = joinpath(@__DIR__(), "..", "data", "asedata.json")
 
@@ -13,7 +13,7 @@ const _rnn = [Float64(d) * u"Å" for d in ase_data["rnn"]]
 const _masses = [Float64(m) * u"u" for m in ase_data["masses"]]
 const _refstates = Dict{String, Any}[ d == nothing ? Dict{String, Any}() : d
                                       for d in ase_data["refstates"]]
-for D in _refstates 
+for D in _refstates
    if haskey(D, "a")
       D["a"] *= u"Å"
    end
@@ -30,9 +30,9 @@ for (n, sym) in enumerate(Symbol.(ase_data["symbols"]))
 end
 
 # --------------------------------------------------------------------
-#  some convenient accessor and conversion functions 
-#  these will not be exported and only accessed via Chemistry.atomic_number etc 
-#  in order to avoid confusion with the AtomsBase interface 
+#  some convenient accessor and conversion functions
+#  these will not be exported and only accessed via Chemistry.atomic_number etc
+#  in order to avoid confusion with the AtomsBase interface
 
 atomic_number(sym::Symbol) = _numbers[sym]
 
@@ -51,11 +51,11 @@ lattice_constant(sym::Symbol) = lattice_constant(atomic_number(sym))
 refstate(z::Integer) = _refstates[z+1]
 refstate(sym::Symbol) = refstate(atomic_number(sym))
 
-# The `rnn` could be exported though since it is quite useful for a lot of 
-# modelling scripts. I'll keep it here for now, unexported. Something to 
-# think about. 
+# The `rnn` could be exported though since it is quite useful for a lot of
+# modelling scripts. I'll keep it here for now, unexported. Something to
+# think about.
 
 rnn(z::Integer) = _rnn[z+1]
 rnn(sym::Symbol) = rnn(atomic_number(sym))
 
-end
+end

src/utils.jl

@@ -17,11 +17,9 @@ function _flexible_system(positions, elements, cell, pbc)
    Nat = length(positions)
    syms = Chemistry.chemical_symbol.(elements)
    atoms = [ Atom(syms[i], positions[i]) for i in 1:Nat ]
-   bc = pbc isa Bool ? (pbc, pbc, pbc) : tuple(pbc...)
-   bb = tuple([cell[i, :] for i = 1:3]...)
-   return FlexibleSystem(atoms; 
-                         cell_vectors = bb, 
-                         periodicity = bc)
+   return FlexibleSystem(atoms;
+                         cell_vectors=tuple([cell[i, :] for i = 1:3]...),
+                         periodicity=pbc)
 end
 
 _set_position(x::Atom, 𝐫) = Atom(atomic_number(x), 𝐫; 
@@ -188,4 +186,4 @@ removed.
 function Base.deleteat!(sys::FlexibleSystem, n)
    deleteat!(sys.particles, n)
    return sys
-end 
+end 

test/test_bulk.jl

@@ -4,56 +4,63 @@ using Test, AtomsBase, Unitful, Random, JSON
 
 ##
 
+@testset "Testing `bulk` and `repeat` against JuLIP reference data" begin
+    test_systems = JSON.parsefile(joinpath(@__DIR__(), "..", "data", "test_systems.json"))
 
-@info("Testing `bulk` and `repeat` against JuLIP reference data")
+    _ustripvecvec(X) = [ ustrip.(x) for x in X ]
 
-test_systems = JSON.parsefile(joinpath(@__DIR__(), "..", "data", "test_systems.json"))
+    compare_system(sys_f, sys_j) = (
+          all( ustrip.( hcat(cell_vectors(sys_f)...) ) .≈ hcat(sys_j["cell"]...)' )  && 
+          all( AtomsBuilder._convert_pbc(sys_j["pbc"]) .== periodicity(sys_f) ) && 
+          all( atomic_number(sys_f, :) .== sys_j["Z"] ) && 
+          all( _ustripvecvec(position(sys_f, :)) .≈ sys_j["X"] )    )
 
-_ustripvecvec(X) = [ ustrip.(x) for x in X ]
+    for D in test_systems
+        @testset "Testing $(D["sym"])" begin
+            if D["sym"] in ("Ti", )
+                cubic = false
+            else
+                cubic = D["cubic"]
+            end
 
-compare_system(sys_f, sys_j) = (
-      all( ustrip.( hcat(cell_vectors(sys_f)...) ) .≈ hcat(sys_j["cell"]...)' )  && 
-      all( AtomsBuilder._convert_pbc(sys_j["pbc"]) .== periodicity(sys_f) ) && 
-      all( atomic_number(sys_f, :) .== sys_j["Z"] ) && 
-      all( _ustripvecvec(position(sys_f, :)) .≈ sys_j["X"] )    )
+            sys_f = bulk(Symbol(D["sym"]); cubic, pbc=D["pbc"])
+            nn = D["nn"] isa Integer ? D["nn"] : tuple(D["nn"]...)
+            if nn != 1
+               sys_f = sys_f * nn
+            end
+            @test compare_system(sys_f, D["sys"])
+        end
+    end
+end
 
-for D in test_systems 
-   nn = D["nn"] isa Integer ? D["nn"] : tuple(D["nn"]...)
-   sys_f = bulk( Symbol(D["sym"]); 
-                 cubic = D["cubic"], 
-                 pbc = D["pbc"] ) 
-   if nn != 1 
-      sys_f = sys_f * nn 
-   end 
-   @test compare_system(sys_f, D["sys"])
+@testset "Test argument errors are raised" begin
+    @test_throws ArgumentError bulk(:Og)  # Case where no bulk structure is known
+    @test_throws ArgumentError bulk(:Og)  # Case where no bulk structure is known
+    @test_throws ArgumentError bulk(:Ti; cubic=true)  # Can't make hcp crystals cubic
 end
 
-##
+@testset "Execute code" begin
+    # not sure how to write a test for this, but at least it should execute
+    sys0 = rattle!(bulk(:C, cubic=true) * (2,3,4), 0.1u"Å")
+    sys1 = rattle!(bulk(:C, cubic=true) * (2,3,4), 0.1)
+    sys2 = rattle!(bulk(:C) * (2,3,4), 0.1)
+    rattle!(bulk(:C) * (2,3,4), 0.01u"nm")
+
+    X = position(sys1, :)
+    Xnew = [ x .+ 0.01u"Å" for x in X ]
+    sys3 = set_positions(sys1, Xnew)
+    @test position(sys3, :) == Xnew
 
-# not sure how to write a test for this, but at least it should execute
-sys0 = rattle!(bulk(:C, cubic=true) * (2,3,4), 0.1u"Å")
-sys1 = rattle!(bulk(:C, cubic=true) * (2,3,4), 0.1)
-sys2 = rattle!(bulk(:C) * (2,3,4), 0.1)
-rattle!(bulk(:C) * (2,3,4), 0.01u"nm")
-
-X = position(sys1, :)
-Xnew = [ x .+ 0.01u"Å" for x in X ]
-sys3 = set_positions(sys1, Xnew)
-@test position(sys3, :) == Xnew
-
-Z = atomic_number(sys1, :)
-@test all(Z .== AtomsBuilder.Chemistry.atomic_number(:C))
-zO = AtomsBuilder.Chemistry.atomic_number(:O)
-Znew = copy(Z); Znew[3:5:end] .= zO
-sys4 = set_elements(sys3, Znew)
-@test all(atomic_number(sys4, :) .== Znew)
-
-pold = deepcopy(sys4.particles) 
-deleteat!(sys4, 1:5)
-@test position.(pold[6:end]) == position(sys4, :)
-@test mass.(pold[6:end]) == mass(sys4, :)
-@test atomic_number.(pold[6:end]) == atomic_number(sys4, :)
-
-
-@test_throws ArgumentError bulk(:Og)  # Case where no bulk structure is known
-@test_throws ArgumentError bulk(:B)   # Tetragonal, which is currently not implemented
+    Z = atomic_number(sys1, :)
+    @test all(Z .== AtomsBuilder.Chemistry.atomic_number(:C))
+    zO = AtomsBuilder.Chemistry.atomic_number(:O)
+    Znew = copy(Z); Znew[3:5:end] .= zO
+    sys4 = set_elements(sys3, Znew)
+    @test all(atomic_number(sys4, :) .== Znew)
+
+    pold = deepcopy(sys4.particles) 
+    deleteat!(sys4, 1:5)
+    @test position.(pold[6:end]) == position(sys4, :)
+    @test mass.(pold[6:end]) == mass(sys4, :)
+    @test atomic_number.(pold[6:end]) == atomic_number(sys4, :)
+end

test/test_utils.jl

@@ -4,7 +4,7 @@ using AtomsBuilder, Test, AtomsBase, Unitful, Random
 ##
 
 Random.seed!(1234)
-sys = bulk(:Ti, cubic=true) * 5
+sys = bulk(:Ti, cubic=false) * 5
 sys = randz!(sys, [ :Ti => 0.5, :O => 0.5 ])
 Z = atomic_number(sys, :)
 @test count( Z .== 8 ) / length(Z) > 0.4