diff --git a/.gitignore b/.gitignore index 5a16984..8d0be91 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ /Manifest.toml /docs/Manifest.toml /docs/build/ +*~ +.*.swp diff --git a/src/AtomsBuilder.jl b/src/AtomsBuilder.jl index 0171037..f837c34 100644 --- a/src/AtomsBuilder.jl +++ b/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 diff --git a/src/bulk.jl b/src/bulk.jl index b6ab040..d64a3af 100644 --- a/src/bulk.jl +++ b/src/bulk.jl @@ -1,11 +1,10 @@ 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] ], @@ -13,34 +12,34 @@ const _unit_cells = Dict( # (positions, cell matrix, factor of a) ) 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 diff --git a/src/chemistry.jl b/src/chemistry.jl index 06539c5..a02eae0 100644 --- a/src/chemistry.jl +++ b/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 \ No newline at end of file +end diff --git a/src/utils.jl b/src/utils.jl index 4530529..5ff27b9 100644 --- a/src/utils.jl +++ b/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 \ No newline at end of file +end diff --git a/test/test_bulk.jl b/test/test_bulk.jl index 3b171eb..8051d39 100644 --- a/test/test_bulk.jl +++ b/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 diff --git a/test/test_utils.jl b/test/test_utils.jl index e1a0c83..ec32ab3 100644 --- a/test/test_utils.jl +++ b/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