Extend AtomsBaseTesting and AtomsBase test routines (#121)

Author: mfherbst

lib/AtomsBaseTesting/Project.toml

@@ -1,7 +1,7 @@
 name = "AtomsBaseTesting"
 uuid = "ed7c10db-df7e-4efa-a7be-4f4190f7f227"
 authors = ["JuliaMolSim community"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"

lib/AtomsBaseTesting/src/AtomsBaseTesting.jl

@@ -6,7 +6,7 @@ using LinearAlgebra
 using Unitful
 using UnitfulAtomic
 
-using AtomsBase: AbstractSystem 
+using AtomsBase: AbstractSystem
 
 export test_approx_eq
 export make_test_system
@@ -17,12 +17,10 @@ properties can be ignored during the comparison using the respective kwargs.
 """
 function test_approx_eq(s::AbstractSystem, t::AbstractSystem;
                         rtol=1e-14, ignore_atprop=Symbol[], ignore_sysprop=Symbol[],
-                        common_only=false)
+                        common_only=false, quiet=false)
     rnorm(a, b) = (ustrip(norm(a)) < rtol ? norm(a - b) / 1unit(norm(a))
                                           : norm(a - b) / norm(a))
 
-    _isinfinite(s) = any(isinf, reduce(vcat, bounding_box(s)))
-
     for method in (length, size, periodicity, )
         @test method(s) == method(t)
     end
@@ -32,7 +30,7 @@ function test_approx_eq(s::AbstractSystem, t::AbstractSystem;
         @test rnorm(method(s, 1), method(t, 1)) < rtol
     end
 
-    # TODO: add element_symbol back in 
+    # TODO: add element_symbol back in
     for method in (species, atomic_symbol, atomic_number, )
         @test method(s, :) == method(t, :)
         @test method(s, 1) == method(t, 1)
@@ -46,6 +44,7 @@ function test_approx_eq(s::AbstractSystem, t::AbstractSystem;
         end
     end
 
+    # test properties of atoms
     if common_only
         test_atprop = [k for k in atomkeys(s) if hasatomkey(t, k)]
     else
@@ -56,7 +55,7 @@ function test_approx_eq(s::AbstractSystem, t::AbstractSystem;
         prop in ignore_atprop && continue
         prop in (:velocity, :position) && continue
         if hasatomkey(s, prop) != hasatomkey(t, prop)
-            println("hashatomkey mismatch for $prop")
+            quiet || println("hashatomkey mismatch for $prop")
             @test hasatomkey(s, prop) == hasatomkey(t, prop)
             continue
         end
@@ -72,6 +71,12 @@ function test_approx_eq(s::AbstractSystem, t::AbstractSystem;
         end
     end
 
+    # Test some things on cell objects
+    @test typeof(cell(s))       == typeof(cell(t))
+    @test periodicity(cell(s))  == periodicity(cell(t))
+    @test n_dimensions(cell(s)) == n_dimensions(cell(t))
+
+    # test properties of systems
     if common_only
         test_sysprop = [k for k in keys(s) if haskey(t, k)]
     else
@@ -81,15 +86,15 @@ function test_approx_eq(s::AbstractSystem, t::AbstractSystem;
     for prop in test_sysprop
         prop in ignore_sysprop && continue
         if haskey(s, prop) != haskey(t, prop)
-            println("haskey mismatch for $prop")
+            quiet || println("haskey mismatch for $prop")
             @test haskey(s, prop) == haskey(t, prop)
             continue
         end
         (haskey(s, prop) && haskey(t, prop)) || continue
 
         if s[prop] isa Quantity
             @test rnorm(s[prop], t[prop]) < rtol
-        elseif prop in (:bounding_box, ) && !(_isinfinite(s))
+        elseif prop in (:bounding_box, ) && (cell(s) isa PeriodicCell)
             @test maximum(map(rnorm, s[prop], t[prop])) < rtol
         else
             @test s[prop] == t[prop]
@@ -104,23 +109,44 @@ Extra atomic or system properties can be specified using `extra_atprop` and `ext
 and specific standard keys can be ignored using `drop_atprop` and `drop_sysprop`.
 """
 function make_test_system(D=3; drop_atprop=Symbol[], drop_sysprop=Symbol[],
-                          extra_atprop=(; ), extra_sysprop=(; ), cellmatrix=:full, 
+                          extra_atprop=(; ), extra_sysprop=(; ), cellmatrix=:full,
                           n_atoms = 5, )
     @assert D == 3
 
-    # Generate some random data to store in Atoms
+    if cellmatrix == :lower_triangular
+        box = ([1.54732, -0.807289, -0.500870]u"Å",
+               [    0.0, 0.4654985, 0.5615117]u"Å",
+               [    0.0,       0.0, 0.7928950]u"Å")
+    elseif cellmatrix == :upper_triangular
+        box = ([1.54732, 0.0, 0.0]u"Å",
+               [-0.807289, 0.4654985, 0.0]u"Å",
+               [-0.500870, 0.5615117, 0.7928950]u"Å")
+    elseif cellmatrix == :diagonal
+        box = ([1.54732, 0.0, 0.0]u"Å",
+              [0.0, 0.4654985, 0.0]u"Å",
+              [0.0, 0.0, 0.7928950]u"Å")
+    else
+        box = ([1.50304, 0.850344, 0.717239]u"Å",
+               [ 0.36113, 1.008144, 0.814712]u"Å",
+               [ 0.06828, 0.381122, 2.129081]u"Å")
+    end
+
+    # Generate some random data to store in atoms and system
     atprop = Dict{Symbol,Any}(
         :position        => [randn(3) for _ = 1:n_atoms]u"Å",
         :velocity        => [randn(3) for _ = 1:n_atoms] * 10^6*u"m/s",
-        #                   Note: reasonable velocity range in au
-        :atomic_symbol   => [:H, :H, :C, :N, :He],
+        #                   Note to above: Reasonable velocity range in au
+        :species         => ChemicalSpecies.([:H, :H, :C, :N, :He]),
         :charge          => [2, 1, 3.0, -1.0, 0.0]u"e_au",
-        :atomic_mass     => 10rand(n_atoms)u"u",
+        :mass            => 10rand(n_atoms)u"u",
         :vdw_radius      => randn(n_atoms)u"Å",
         :covalent_radius => randn(n_atoms)u"Å",
         :magnetic_moment => [0.0, 0.0, 1.0, -1.0, 0.0],
     )
     sysprop = Dict{Symbol,Any}(
+        :bounding_box => box,
+        :periodicity  => (true, true, false),
+        #
         :extra_data   => 42,
         :charge       => -1u"e_au",
         :multiplicity => 2,
@@ -136,36 +162,23 @@ function make_test_system(D=3; drop_atprop=Symbol[], drop_sysprop=Symbol[],
     atprop  = merge(atprop,  pairs(extra_atprop))
 
     atoms = map(1:n_atoms) do i
-        atargs = Dict(k => v[i] for (k, v) in pairs(atprop)
-                      if !(k in (:position, :velocity)))
+        atargs = Dict(k => v[i] for (k, v) in pairs(atprop) if !(k in (:position, :velocity)))
         if haskey(atprop, :velocity)
-            Atom(atprop[:atomic_symbol][i], atprop[:position][i], atprop[:velocity][i];
-                 atargs...)
+            Atom(atprop[:species][i], atprop[:position][i], atprop[:velocity][i]; atargs...)
         else
-            Atom(atprop[:atomic_symbol][i], atprop[:position][i]; atargs...)
+            Atom(atprop[:species][i], atprop[:position][i]; atargs...)
         end
     end
-    if cellmatrix == :lower_triangular
-        box = ([1.54732, -0.807289, -0.500870]u"Å",
-               [    0.0, 0.4654985, 0.5615117]u"Å",
-               [    0.0,       0.0, 0.7928950]u"Å")
-    elseif cellmatrix == :upper_triangular
-        box = ([1.54732, 0.0, 0.0]u"Å",
-               [-0.807289, 0.4654985, 0.0]u"Å",
-               [-0.500870, 0.5615117, 0.7928950]u"Å")
-    elseif cellmatrix == :diagonal
-        box = ([1.54732, 0.0, 0.0]u"Å",
-              [0.0, 0.4654985, 0.0]u"Å",
-              [0.0, 0.0, 0.7928950]u"Å")
-    else
-        box = ([1.50304, 0.850344, 0.717239]u"Å",
-               [ 0.36113, 1.008144, 0.814712]u"Å",
-               [ 0.06828, 0.381122, 2.129081]u"Å")
-    end
-    pbcs = (true, true, false)
-    system = atomic_system(atoms, box, pbcs; sysprop...)
+    cell = PeriodicCell(; cell_vectors=sysprop[:bounding_box],
+                          periodicity=sysprop[:periodicity])
+
+    sysargs = Dict(k => v for (k, v) in pairs(sysprop)
+                   if !(k in (:bounding_box, :periodicity)))
+    system = FlexibleSystem(atoms, cell; sysargs...)
 
-    (; system, atoms, atprop=NamedTuple(atprop), sysprop=NamedTuple(sysprop), box, pbcs)
+    (; system, atoms, cell,
+     bounding_box=sysprop[:bounding_box], periodicity=sysprop[:periodicity],
+       atprop=NamedTuple(atprop), sysprop=NamedTuple(sysprop))
 end
 
 end

lib/AtomsBaseTesting/test/runtests.jl

@@ -9,6 +9,16 @@ include("testmacros.jl")
 
 @testset "AtomsBaseTesting.jl" begin
     @testset "make_test_system" begin
+        let case = make_test_system()
+            # Data that is delivered agrees with the constructed system
+            # TODO Could test more here
+            @test sort(collect(keys(case.atprop)))  == sort(collect(atomkeys(case.system)))
+            @test sort(collect(keys(case.atprop)))  == sort(collect(keys(case.atoms[1])))
+            @test sort(collect(keys(case.sysprop))) == sort(collect(keys(case.system)))
+            @test case.bounding_box  == bounding_box(case.system)
+            @test case.periodicity   == periodicity(case.system)
+        end
+
         let case = make_test_system(; cellmatrix=:full)
             box = reduce(hcat, bounding_box(case.system))
             @test UpperTriangular(box) != box
@@ -33,7 +43,7 @@ include("testmacros.jl")
             @test UpperTriangular(box) == box
             @test LowerTriangular(box) == box
         end
- 
+
         @test  hasatomkey(make_test_system().system,                            :vdw_radius)
         @test !hasatomkey(make_test_system(; drop_atprop=[:vdw_radius]).system, :vdw_radius)
 
@@ -52,9 +62,9 @@ include("testmacros.jl")
         # once we require Julia 1.7
         case = make_test_system()
         system = case.system
-        atoms = case.atoms
-        box = case.box
-        bcs = case.pbcs
+        atoms  = case.atoms
+        box    = case.bounding_box
+        bcs    = case.periodicity
         sysprop = case.sysprop
         # end simplify
 
@@ -71,17 +81,17 @@ include("testmacros.jl")
         # once we require Julia 1.7
         case = make_test_system()
         system = case.system
-        atoms = case.atoms
-        box = case.box
-        bcs = case.pbcs
+        atoms  = case.atoms
+        box    = case.bounding_box
+        bcs    = case.periodicity
         sysprop = case.sysprop
         # end simplify
 
         sysprop_dict = Dict(pairs(sysprop))
         pop!(sysprop_dict, :multiplicity)
         system_edit = atomic_system(atoms, box, bcs; sysprop_dict...)
 
-        @testfail test_approx_eq(system, system_edit)
+        @testfail test_approx_eq(system, system_edit, quiet=true)
         @testpass test_approx_eq(system, system_edit; ignore_sysprop=[:multiplicity])
         @testpass test_approx_eq(system, system_edit; common_only=true)
     end

src/AtomsBase.jl

@@ -5,12 +5,12 @@ using UnitfulAtomic
 using StaticArrays
 using Requires
 
-export Atom, FlexibleSystem, FastSystem, AbstractSystem 
+export Atom, FlexibleSystem, FastSystem, AbstractSystem
 
-# Main Interface specification and inline docs 
+# Main Interface specification and inline docs
 include("interface.jl")
 
-# utilities useful to share across implementations 
+# utilities useful to share across implementations
 include("utils/cells.jl")
 include("utils/chemspecies.jl")
 include("utils/properties.jl")
@@ -19,15 +19,15 @@ include("utils/show.jl")
 include("utils/atomview.jl")
 
 
-# prototype implementations 
+# prototype implementations
 include("implementation/atom.jl")
 include("implementation/flexible_system.jl")
 include("implementation/fast_system.jl")
 include("implementation/utils.jl")
 
 
-# TODO: 
-#  - this should be converted to an extension 
+# TODO:
+#  - this should be converted to an extension
 #  - should work for AbstractSystem
 function __init__()
     @require AtomsView="ee286e10-dd2d-4ff2-afcb-0a3cd50c8041" begin

src/implementation/atom.jl

@@ -16,10 +16,10 @@ struct Atom{D, L<:Unitful.Length, V<:Unitful.Velocity, M<:Unitful.Mass}
     data::Dict{Symbol, Any}  # Store arbitrary data about the atom.
 end
 
-velocity(atom::Atom)      = atom.velocity
-position(atom::Atom)      = atom.position
-mass(atom::Atom)   = atom.mass
-species(atom::Atom)       = atom.species
+velocity(atom::Atom) = atom.velocity
+position(atom::Atom) = atom.position
+mass(atom::Atom)     = atom.mass
+species(atom::Atom)  = atom.species
 
 n_dimensions(::Atom{D}) where {D} = D
 
@@ -68,7 +68,7 @@ function _default_velocity(position::AbstractVector{L}) where {L <: Unitful.Leng
     elseif uL == u"bohr" 
         return zeros(TFL, length(position))u"nm/s"
     elseif uL == u"m" 
-        return zeros(TFL, length(position))u"m/s"    
+        return zeros(TFL, length(position))u"m/s"
     end 
     @warn("Cannot infer default velocity for position with unit $(unit(position[1]))")
     return zeros(TFL, length(position)) * (uL / u"s")

src/implementation/fast_system.jl

@@ -27,7 +27,7 @@ function FastSystem(cϵll::TCELL, positions::AbstractVector{<: SVector{D, L}},
                     species::AbstractVector{S}, masses) where {TCELL, D, L, S} 
     if D != n_dimensions(cϵll)
         throw(ArgumentError("Cell dimension D=$(n_dimensions(cϵll)) does not match particle dimension D=$D."))
-    end                    
+    end
     FastSystem{D, TCELL, L, typeof(masses), S}(
         cϵll, positions, species, masses)
 end
@@ -91,8 +91,6 @@ Base.getindex(system::FastSystem, ::Colon, x::Symbol) = getfield(system, x)
 position(s::FastSystem, ::Colon) = s.position
 position(sys::FastSystem, i::Union{Integer, AbstractVector}) = sys.position[i]
 
-velocity(::FastSystem, args...) = missing
-
 mass(s::FastSystem, ::Colon) = s.mass
 mass(sys::FastSystem, i::Union{Integer, AbstractVector}) = sys.mass[i]
 

src/implementation/flexible_system.jl

@@ -116,7 +116,7 @@ velocity(sys::FlexibleSystem, i::Integer) =
         velocity(sys.particles[i])
 
 velocity(sys::FlexibleSystem, i::Union{AbstractVector, Colon}) = 
-        [ velocity(x) for x in sys.particles[i] ]         
+        [ velocity(x) for x in sys.particles[i] ]
 
 mass(sys::FlexibleSystem, i::Integer) = 
         mass(sys.particles[i])
@@ -128,4 +128,4 @@ species(sys::FlexibleSystem, i::Integer) =
         species(sys.particles[i])
 
 species(sys::FlexibleSystem, i::Union{AbstractVector, Colon}) = 
-        [ species(x) for x in sys.particles[i] ]                         
+        [ species(x) for x in sys.particles[i] ]                         

src/implementation/utils.jl

@@ -22,11 +22,11 @@ julia> hydrogen = atomic_system([:H => [0, 0, 1.]u"bohr",
                                   bounding_box, pbcs)
 ```
 """
-atomic_system(atoms::AbstractVector{<:Atom}, box, bcs; kwargs...) = 
-      FlexibleSystem(atoms, box, bcs; kwargs...)
+atomic_system(atoms::AbstractVector{<:Atom}, box, pbcs; kwargs...) = 
+      FlexibleSystem(atoms, box, pbcs; kwargs...)
 
-atomic_system(atoms::AbstractVector, box, bcs; kwargs...) = 
-      FlexibleSystem(convert.(Atom, atoms), box, bcs; kwargs...)
+atomic_system(atoms::AbstractVector, box, pbcs; kwargs...) = 
+      FlexibleSystem(convert.(Atom, atoms), box, pbcs; kwargs...)
 
 
 """