diff --git a/README.md b/README.md index 1709039..01f9dc2 100644 --- a/README.md +++ b/README.md @@ -12,46 +12,57 @@ AtomsBase is an abstract interface for representation of atomic geometries in Ju * file I/O with standard formats (.cif, .xyz, ...) * numerical tools: sampling, integration schemes, etc. * automatic differentiation and machine learning systems +* visualization (e.g. plot recipes) -Currently, the design philosophy is to be as lightweight as possible, with only a small set of required function dispatches. We will also provide a couple of standard concrete implementations of the interface that we envision could be broadly applicable, but encourage developers to provide their own implementations as needed in new or existing packages. +Currently, the design philosophy is to be as lightweight as possible, with only a small set of required function dispatches to make adopting the interface into existing packages easy. We also provide a couple of standard concrete implementations of the interface that we envision could be broadly applicable, but encourage developers to provide their own implementations as needed in new or existing packages. ## Overview -AtomsBase defines a few abstract types used for specifying an atomic system. We will describe them briefly here, from the "top down." Users and/or prospective developers may also find `implementation_simple.jl` a useful reference for a simple concrete implementation of the interface. +The main abstract type introduced in AtomsBase `AbstractSystem{D,S}`. The `D` parameter indicates the number of spatial dimensions in the system, and `S` indicates the type identifying an individual species in this system. -**A remark on SoA vs. AoS:** The "struct-of-arrays" (SoA) vs. "array-of-structs" (AoS) is a common design dilemma in representations of systems such as these. We have deliberately designed this interface to be _agnostic_ to how a concrete implementation chooses to structure its data. Some specific notes regarding how implementations might differ for these two paradigms are included below. +The main power of the interface comes from predictable behavior of several core functions to access information about a system. Various categories of such functions are described below. -A way to think about this broadly is that the difference amounts to the ordering of function calls. For example, to get the position of a single particle in an AoS implementation, the explicit funciton chaining would be `position(getindex(sys))` (i.e. extract the single struct representing the particle of interest and query its position), while for SoA, one should prefer `getindex(position(sys))` (extract the array of positions, then index into it for a single particle). The beauty of an abstract interface in Julia is that these details can be, in large part, abstracted away through method dispatch such that the end user sees the same expected behavior irrespective of how things are implemented "under the hood." +### System geometry +Functions that need to be dispatched: +* `bounding_box(::AbstractSystem{D})::SVector{D,SVector{D,<:Unitful.Length}}`: returns `D` vectors of length `D` that describe the "box" in which the system lives +* `boundary_conditions(::AbstractSystem{D})::SVector{D,BoundaryCondition})`: returns a vector of length `D` of `BoundaryCondition` objects to describe what happens at the edges of the box -### System -An object describing a system should be a subtype of `AbstractSystem` and will in general store identifiers, positions, and (if relevant) velocities of the particles that make it up, as well as a set of coordinate bounds. +Functions that will work automatically: +* `get_periodic`: returns a vector of length `D` of `Bool`s for whether each dimension of the system has periodic boundary conditions +* `n_dimensions`: returns `D`, the number of spatial dimensions of the system -An `AbstractSystem` takes several type parameters: `D` (an integer representing the dimensionality of the system), `ET<:AbstractElement`, and `AT<:AbstractParticle` (see below) to indicate what types of particles these are, and requires dispatch of the following functions: -* `bounding_box(::AbstractSystem{D})::SVector{D, SVector{D, <:Unitful.Length}}`: should return a set of basis vectors describing the boundaries of the coordinates in which the system resides -* `boundary_conditions(::AbstractSystem{D})::SVector{D, BoundaryCondition}`: returns the boundary conditions corresponding to each spatial dimension of the system (see below for more on the `BoundaryCondition` type) +### Iteration and Indexing over systems +There is a presumption of the ability to somehow extract an individual component (e.g. a single atom or molecule) of this system, though there are no constraints on the type of this component. To achieve this, an `AbstractSystem` object is expected to implement the Julia interfaces for [iteration](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration) and [indexing](https://docs.julialang.org/en/v1/manual/interfaces/#Indexing) in order to access representations of individual components of the system. Some default dispatches of parts of these interfaces are already included, so the minimal set of functions to dispatch in a concrete implementation is `Base.getindex` and `Base.length`, though it may be desirable to customize additional behavior depending on context. -`AbstractSystem` implements [Julia's indexing interface](https://docs.julialang.org/en/v1/manual/interfaces/#Indexing), thus the following functions must also be dispatched: -* `Base.getindex(::AbstractSystem, ::Int)` -* `Base.size(::AbstractSystem)` +### System state and properties +The only required properties to be specified of the system is the species of each component of the system and the positions and velocities associated with each component. These are accessed through the functions `species`, `position`, and `velocity`, respectively. The default dispatch of these functions onto an `AbstractSystem` object is as a broadcast over it, which will "just work" provided the indexing/iteration interfaces have been implemented (see above) and the functions are defined on individual system components. -By default, the `position` and `velocity` functions dispatch onto `AbstractSystem` objects as a broadcast over the system (e.g. `position(sys::AbstractSystem) = position.(sys)`, invoking the dispatch of `position` onto the `AbstractParticle` type parameter of the system). In an SoA implementation, custom dispatches should probably be included to avoid the construction of the particle objects. +As a concrete example, AtomsBase provides the `StaticAtom` type as this is anticipated to be a commonly needed representation. Its implementation looks as follows: +```julia +struct StaticAtom{D,L<:Unitful.Length} + position::SVector{D,L} + element::Element +end +StaticAtom(position, element) = StaticAtom{length(position)}(position, element) +position(atom::StaticAtom) = atom.position +species(atom::StaticAtom) = atom.element +``` +Note that the default behavior of `velocity` is to return `missing`, so it doesn't need to be explicitly dispatched here. -**A remark on fractional coordinates:** In many contexts, it is desirable to work with fractional coordinates, i.e. unitless multiples of some reference coordinate system (typically this reference is what would be returned by the `box` function above). Because a focus of this interface is interoperability and unitless, non-Cartesian coordinates introduce ambiguity, we've chosen to impose that coordinates returned by functions in the interface be unitful quantities in a Cartesian frame. Of course, a concrete implementation of the interface could (and likely would need to) also dispatch additional functions for working with fractional coordinates. +The two sample implementations provided in this repo are both "composed" of `StaticAtom` objects; refer to them as well as `sandbox/aos_vs_soa.jl` to see how this can work in practice. -### Particles -Particle objects are subtypes of `AbstractParticle`, and also take a type parameter that is `<:AbstractElement` (see below) to indicate how particles are identified. +The `position`, `velocity`, and `species` functions also have indexed signatures to extract a given element directly, as in, for example: -In the SoA case, particle objects would only ever be explicitly constructed when `getindex` is invoked on the system, as a "view" into the system. +```julia +position(sys, i) # position of `i`th particle in `sys` +``` +Currently, this syntax only supports linear indexing. -Particle objects should dispatch methods of the following functions: -* `position(::AbstractParticle)::AbstractVector{<: Unitful.Length}` -* `element(::AbstractParticle)::AbstractElement` +### Struct-of-Arrays vs. Array-of-Structs +The "struct-of-arrays" (SoA) vs. "array-of-structs" (AoS) is a common design dilemma in representations of systems such as these. We have deliberately designed this interface to be _agnostic_ to how a concrete implementation chooses to structure its data. Some specific notes regarding how implementations might differ for these two paradigms are included below. -And, optionally, -* `velocity(::AbstractParticle)::AbstractVector{<: Unitful.Velocity}`, which defaults to returning `nothing` if not dispatched. -### Elements -Subtypes of `AbstractElement` serve as identifiers of particles. As the name suggests, a common case would be a chemical element (e.g. for a DFT simulation). However, it could also contain more detailed isotopic/spin information if that is necessary, or be a molecule (e.g. in MD), or even a much larger-scale object! +A way to think about this broadly is that the difference amounts to the ordering of function calls. For example, to get the position of a single particle in an AoS implementation, the explicit function chaining would be `position(getindex(sys))` (i.e. extract the single struct representing the particle of interest and query its position), while for SoA, one should prefer an implementation like `getindex(position(sys))` (extract the array of positions, then index into it for a single particle). The beauty of an abstract interface in Julia is that these details can be, in large part, abstracted away through method dispatch such that the end user sees the same expected behavior irrespective of how things are implemented "under the hood." -For simulation purposes, the utility of this object would be to serve as a sufficiently specific "index" into a database of simulation parameters (e.g. a pseudopotential library for DFT, or interparticle potential parameters for MD). Because we envision a chemical element being an extremely common case, we provide an explicit subtype in the form of `Element`, which makes use of [PeriodicTable.jl](https://github.com/JuliaPhysics/PeriodicTable.jl) to access information such as atomic numbers and masses. +To demonstrate this, we provide two simple concrete implementations of the interface in `implementation_soa.jl` and `implementation_aos.jl` to show how analogous systems could be constructed within these paradigms (including the `getindex` implementations mentioned above). See also `sandbox/aos_vs_soa.jl` for how they can actually be constructed and queried. ### Boundary Conditions Finally, we include support for defining boundary conditions. Currently included are `Periodic` and `DirichletZero`. There should be one boundary condition specified for each spatial dimension represented. diff --git a/sandbox/aos_vs_soa.jl b/sandbox/aos_vs_soa.jl index 9802ed7..63f26ba 100644 --- a/sandbox/aos_vs_soa.jl +++ b/sandbox/aos_vs_soa.jl @@ -2,13 +2,17 @@ using AtomsBase using StaticArrays using Unitful +import PeriodicTable + +periodic_table = PeriodicTable.elements + box = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]u"m" bcs = [Periodic(), Periodic(), Periodic()] positions = [0.25 0.25 0.25; 0.75 0.75 0.75]u"m" -elements = [ChemicalElement(:C), ChemicalElement(:C)] +elements = [periodic_table[:C], periodic_table[:C]] -atom1 = SimpleAtom(SVector{3}(positions[1,:]),elements[1]) -atom2 = SimpleAtom(SVector{3}(positions[2,:]),elements[2]) +atom1 = StaticAtom(SVector{3}(positions[1,:]),elements[1]) +atom2 = StaticAtom(SVector{3}(positions[2,:]),elements[2]) aos = FlexibleSystem(box, bcs, [atom1, atom2]) soa = FastSystem(box, bcs, positions, elements) diff --git a/src/AtomsBase.jl b/src/AtomsBase.jl index 62ffdb8..c7d150c 100644 --- a/src/AtomsBase.jl +++ b/src/AtomsBase.jl @@ -1,9 +1,8 @@ module AtomsBase include("interface.jl") +include("atoms.jl") include("implementation_soa.jl") include("implementation_aos.jl") -# Write your package code here. - end diff --git a/src/atoms.jl b/src/atoms.jl new file mode 100644 index 0000000..67c5bb4 --- /dev/null +++ b/src/atoms.jl @@ -0,0 +1,35 @@ +# +# Extra stuff only for Systems composed of atoms +# + +export StaticAtom, AbstractAtomicSystem +export atomic_mass, atomic_number, atomic_symbol, atomic_property + +struct StaticAtom{D,L<:Unitful.Length} + position::SVector{D,L} + element::Element +end +StaticAtom(position, element) = StaticAtom{length(position)}(position, element) +position(atom::StaticAtom) = atom.position +species(atom::StaticAtom) = atom.element + +function StaticAtom(position, symbol::Union{Integer,AbstractString,Symbol,AbstractVector}) + StaticAtom(position, elements[symbol]) +end + +function Base.show(io::IO, a::StaticAtom) + print(io, "StaticAtom: $(a.element.symbol)") +end + +const AbstractAtomicSystem{D} = AbstractSystem{D,Element} + +atomic_symbol(a::StaticAtom) = a.element.symbol +atomic_mass(a::StaticAtom) = a.element.atomic_mass +atomic_number(a::StaticAtom) = a.element.number +atomic_property(a::StaticAtom, property::Symbol) = getproperty(a.element, property) + +atomic_symbol(sys::AbstractAtomicSystem) = atomic_symbol.(sys) +atomic_number(sys::AbstractAtomicSystem) = atomic_number.(sys) +atomic_mass(sys::AbstractAtomicSystem) = atomic_mass.(sys) +atomic_property(sys::AbstractAtomicSystem, property::Symbol)::Vector{Any} = + atomic_property.(sys, property) diff --git a/src/implementation_aos.jl b/src/implementation_aos.jl index f818281..45b388b 100644 --- a/src/implementation_aos.jl +++ b/src/implementation_aos.jl @@ -1,32 +1,48 @@ -# Example implementation using as few function definitions as possible # +# Implementation of AtomsBase interface in an array-of-structs style +# + using StaticArrays export FlexibleSystem # TODO Switch order of type arguments? -struct FlexibleSystem{D,ET<:AbstractElement,AT<:AbstractParticle{ET}} <: AbstractSystem{D,ET,AT} - box::SVector{D,<:SVector{D,<:Unitful.Length}} +struct FlexibleSystem{D,S,L<:Unitful.Length,AT} <: AbstractSystem{D,S} + box::SVector{D,<:SVector{D,L}} boundary_conditions::SVector{D,<:BoundaryCondition} particles::Vector{AT} + FlexibleSystem(box, boundary_conditions, particles) = new{ + length(boundary_conditions), + eltype(elements), + eltype(eltype(box)), + eltype(particles), + }( + box, + boundary_conditions, + particles, + ) end # convenience constructor where we don't have to preconstruct all the static stuff... function FlexibleSystem( - box::Vector{Vector{L}}, + box::Vector{<:AbstractVector{L}}, boundary_conditions::Vector{BC}, particles::Vector{AT}, -) where {BC<:BoundaryCondition,L<:Unitful.Length,AT<:AbstractParticle} +) where {BC<:BoundaryCondition,L<:Unitful.Length,AT} D = length(box) if !all(length.(box) .== D) throw(ArgumentError("box must have D vectors of length D")) end - FlexibleSystem(SVector{D,SVector{D,L}}(box), SVector{D,BoundaryCondition}(boundary_conditions), particles) + FlexibleSystem( + SVector{D,SVector{D,L}}(box), + SVector{D,BoundaryCondition}(boundary_conditions), + particles, + ) end bounding_box(sys::FlexibleSystem) = sys.box boundary_conditions(sys::FlexibleSystem) = sys.boundary_conditions -function Base.show(io::IO, ::MIME"text/plain", sys::FlexibleSystem) +function Base.show(io::IO, sys::FlexibleSystem) print(io, "FlexibleSystem with ", length(sys), " particles") end diff --git a/src/implementation_simple.jl b/src/implementation_simple.jl deleted file mode 100644 index 1bb9141..0000000 --- a/src/implementation_simple.jl +++ /dev/null @@ -1,40 +0,0 @@ -# Example implementation using as few function definitions as possible -# -using StaticArrays - -export SimpleAtom, SimpleSystem - -struct SimpleAtom{D, T<:Unitful.Length} <: AbstractAtom - position::SVector{D, T} - element::ChemicalElement -end -SimpleAtom(position, element) = SimpleAtom{length(position), eltype(position)}(position, element) -position(atom::SimpleAtom) = atom.position -element(atom::SimpleAtom) = atom.element - -function SimpleAtom(position, symbol::Union{Integer,AbstractString,Symbol,AbstractVector}) - SimpleAtom(position, ChemicalElement(symbol)) -end - -# TODO Switch order of type arguments? -struct SimpleSystem{D, ET<:AbstractElement, AT<:AbstractParticle{ET}, T<:Unitful.Length} <: AbstractSystem{D,ET,AT} - box::SVector{D, SVector{D, T}} - boundary_conditions::SVector{D, BoundaryCondition} - particles::Vector{AT} -end - -function SimpleSystem(box, boundary_conditions, particles) - D = length(box) - ET = typeof(element(first(particles))) - AT = eltype(particles) - T = eltype(first(box)) - - SimpleSystem{D, ET, AT, T}(box, boundary_conditions, particles) -end - -bounding_box(sys::SimpleSystem) = sys.box -boundary_conditions(sys::SimpleSystem) = sys.boundary_conditions - -Base.size(sys::SimpleSystem) = size(sys.particles) -Base.length(sys::SimpleSystem) = length(sys.particles) -Base.getindex(sys::SimpleSystem, i::Int) = getindex(sys.particles, i) diff --git a/src/implementation_soa.jl b/src/implementation_soa.jl index 1e682e7..fd220d8 100644 --- a/src/implementation_soa.jl +++ b/src/implementation_soa.jl @@ -1,18 +1,19 @@ -# Example implementation using as few function definitions as possible # +# Implementation of AtomsBase interface in a struct-of-arrays style. +# + using StaticArrays export FastSystem -struct FastSystem{D,ET<:AbstractElement,AT<:AbstractParticle{ET},L<:Unitful.Length} <: - AbstractSystem{D,ET,AT} +struct FastSystem{D,S,L<:Unitful.Length} <: AbstractSystem{D,S} box::SVector{D,SVector{D,L}} boundary_conditions::SVector{D,BoundaryCondition} positions::Vector{SVector{D,L}} - elements::Vector{ET} + elements::Vector{S} # janky inner constructor that we need for some reason FastSystem(box, boundary_conditions, positions, elements) = - new{length(boundary_conditions),eltype(elements),SimpleAtom,eltype(eltype(positions))}( + new{length(boundary_conditions),eltype(elements),eltype(eltype(positions))}( box, boundary_conditions, positions, @@ -22,11 +23,11 @@ end # convenience constructor where we don't have to preconstruct all the static stuff... function FastSystem( - box::AbstractVector{Vector{L}}, + box::Vector{<:AbstractVector{L}}, boundary_conditions::AbstractVector{BC}, positions::AbstractMatrix{M}, - elements::AbstractVector{ET}, -) where {L<:Unitful.Length,BC<:BoundaryCondition,M<:Unitful.Length,ET<:AbstractElement} + elements::AbstractVector{S}, +) where {L<:Unitful.Length,BC<:BoundaryCondition,M<:Unitful.Length,S} N = length(elements) D = length(box) if !all(length.(box) .== D) @@ -46,7 +47,7 @@ function FastSystem( ) end -function Base.show(io::IO, ::MIME"text/plain", sys::FastSystem) +function Base.show(io::IO, sys::FastSystem) print(io, "FastSystem with ", length(sys), " particles") end @@ -54,8 +55,14 @@ bounding_box(sys::FastSystem) = sys.box boundary_conditions(sys::FastSystem) = sys.boundary_conditions # Base.size(sys::FastSystem) = size(sys.particles) -Base.length(sys::FastSystem{D,ET,AT}) where {D,ET,AT} = length(sys.elements) +Base.length(sys::FastSystem{D,S}) where {D,S} = length(sys.elements) + +Base.getindex(sys::FastSystem{D,S,L}, i::Int) where {D,S,L} = + StaticAtom{D,L}(sys.positions[i], sys.elements[i]) + +# these dispatches aren't strictly necessary, but they make these functions ~2x faster +position(s::FastSystem) = s.positions +species(s::FastSystem) = s.elements -# first piece of trickiness: can't do a totally abstract dispatch here because we need to know the signature of the constructor for AT -Base.getindex(sys::FastSystem{D,ET,SimpleAtom}, i::Int) where {D,ET} = - SimpleAtom{D}(sys.positions[i], sys.elements[i]) +position(s::FastSystem, i) = s.positions[i] +species(s::FastSystem, i) = s.elements[i] diff --git a/src/interface.jl b/src/interface.jl index 66732ce..fc54c5f 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -4,74 +4,31 @@ using PeriodicTable using StaticArrays import Base.position -export AbstractElement, AbstractParticle, AbstractAtom, AbstractSystem, AbstractAtomicSystem -export ChemicalElement, SimpleAtom +export AbstractSystem export BoundaryCondition, DirichletZero, Periodic -export atomic_mass, - atomic_number, - atomic_symbol, - bounding_box, - element, - position, - velocity, - boundary_conditions, - periodic_dims -export atomic_property, has_atomic_property, atomic_propertynames -export n_dimensions - - -abstract type AbstractElement end -struct ChemicalElement <: AbstractElement - data::PeriodicTable.Element -end - -ChemicalElement(symbol::Union{Symbol,Integer,AbstractString}) = - ChemicalElement(PeriodicTable.elements[symbol]) -Base.show(io::IO, elem::ChemicalElement) = print(io, "Element(", atomic_symbol(elem), ")") - -# These are always only read-only ... and allow look-up into a database -atomic_symbol(el::ChemicalElement) = el.data.symbol -atomic_number(el::ChemicalElement) = el.data.number -atomic_mass(el::ChemicalElement) = el.data.atomic_mass - - - -# -# A distinguishable particle, can be anything associated with coordinate -# information (position, velocity, etc.) -# most importantly: Can have any identifier type -# -# IdType: Type used to identify the particle -# -abstract type AbstractParticle{ET<:AbstractElement} end -velocity(::AbstractParticle)::AbstractVector{<:Unitful.Velocity} = missing -position(::AbstractParticle)::AbstractVector{<:Unitful.Length} = error("Implement me") -(element(::AbstractParticle{ET})::ET) where {ET<:AbstractElement} = error("Implement me") +export species, position, velocity +export bounding_box, boundary_conditions, is_periodic, n_dimensions +""" + velocity(p) -# -# The atom type itself -# - The atom interface is read-only (to allow as simple as possible implementation) -# Writability may be supported in derived or concrete types. -# - The inferface is only in Cartesian coordinates. -# - Has atom-specific defaults (i.e. assumes every entity represents an atom or ion) -# +Return the velocity of a particle `p`. +""" +velocity(p)::Union{Unitful.Velocity,Missing} = missing -const AbstractAtom = AbstractParticle{ChemicalElement} -element(::AbstractAtom)::ChemicalElement = error("Implement me") +""" + position(p) +Return the position of a particle `p`. +""" +position(p)::Unitful.Length = error("Implement me") -# Extracting things ... it might make sense to make some of them writable in concrete -# implementations, therefore these interfaces are forwarded from the Element object. -atomic_symbol(atom::AbstractAtom) = atomic_symbol(element(atom)) -atomic_number(atom::AbstractAtom) = atomic_number(element(atom)) -atomic_mass(atom::AbstractAtom) = atomic_mass(element(atom)) +""" + species(p) -# Custom atomic properties: -atomic_property(::AbstractAtom, ::Symbol, default = missing) = default -has_atomic_property(atom::AbstractAtom, property::Symbol) = - !ismissing(atomic_property(atom, property)) -atomic_propertynames(::AbstractAtom) = Symbol[] +Return the species of a particle `p`. +""" +species(p) = error("Implement me") # # Identifier for boundary conditions per dimension @@ -81,82 +38,80 @@ struct DirichletZero <: BoundaryCondition end # Dirichlet zero boundary (i.e. m struct Periodic <: BoundaryCondition end # Periodic BCs -# -# The system type -# Again readonly. -# +""" + AbstractSystem{D,S} + +A `D`-dimensional system comprised of particles identified by type `S`. +""" +abstract type AbstractSystem{D,S} end -abstract type AbstractSystem{D,ET<:AbstractElement,AT<:AbstractParticle{ET}} end +""" + bounding_box(sys::AbstractSystem{D}) + +Return a vector of length `D` of vectors of length `D` that describe the "box" in which the system `sys` is defined. +""" (bounding_box(::AbstractSystem{D})::SVector{D,SVector{D,<:Unitful.Length}}) where {D} = error("Implement me") + +""" + boundary_conditions(sys::AbstractSystem{D}) + +Return a vector of length `D` of `BoundaryCondition` objects, one for each direction described by `bounding_box(sys)`. +""" (boundary_conditions(::AbstractSystem{D})::SVector{D,BoundaryCondition}) where {D} = error("Implement me") -get_periodic(sys::AbstractSystem) = - [isa(bc, Periodic) for bc in get_boundary_conditions(sys)] +is_periodic(sys::AbstractSystem) = [isa(bc, Periodic) for bc in boundary_conditions(sys)] -# Note: Can't use ndims, because that is ndims(sys) == 1 (because of AbstractVector interface) +# Note: Can't use ndims, because that is ndims(sys) == 1 (because of indexing interface) n_dimensions(::AbstractSystem{D}) where {D} = D - -# indexing interface -Base.getindex(::AbstractSystem, ::Int) = error("Implement me") -Base.size(::AbstractSystem) = error("Implement me") -Base.length(::AbstractSystem) = error("Implement me") +# indexing and iteration interface...need to implement getindex and length, here are default dispatches for others +Base.size(s::AbstractSystem) = (length(s),) Base.setindex!(::AbstractSystem, ::Int) = error("AbstractSystem objects are not mutable.") Base.firstindex(::AbstractSystem) = 1 Base.lastindex(s::AbstractSystem) = length(s) -Base.iterate(S::AbstractSystem, i::Int=1) = (1 <= i <= length(S)) ? (@inbounds S[i], i+1) : nothing - -# iteration interface, needed for default broadcast dispatches below to work -Base.iterate(sys::AbstractSystem{D,ET,AT}, state = firstindex(sys)) where {D,ET,AT} = - state > length(sys) ? nothing : (sys[state], state + 1) +# default to 1D indexing +Base.iterate(S::AbstractSystem, state = firstindex(S)) = + (firstindex(S) <= state <= length(S)) ? (@inbounds S[state], state + 1) : nothing # TODO Support similar, push, ... -# Some implementations might prefer to store data in the System as a flat list and -# expose Atoms as a view. Therefore these functions are needed. Of course this code -# should be autogenerated later on ... +""" + position(sys::AbstractSystem{D}) + position(sys::AbstractSystem, index) + +Return a vector of positions of every particle in the system `sys`. Return type should be a vector of vectors each containing `D` elements that are `<:Unitful.Length`. If an index is passed, return only the position of the particle at that index. +""" position(sys::AbstractSystem) = position.(sys) # in Cartesian coordinates! +position(sys::AbstractSystem, index) = position(sys[index]) + +""" + velocity(sys::AbstractSystem{D}) + velocity(sys::AbstractSystem, index) + +Return a vector of velocities of every particle in the system `sys`. Return type should be a vector of vectors each containing `D` elements that are `<:Unitful.Velocity`. If an index is passed, return only the velocity of the particle at that index. +""" velocity(sys::AbstractSystem) = velocity.(sys) # in Cartesian coordinates! -element(sys::AbstractSystem) = element.(sys) +velocity(sys::AbstractSystem, index) = velocity(sys[index]) -# -# Extra stuff only for Systems composed of atoms -# -const AbstractAtomicSystem{D,AT<:AbstractAtom} = AbstractSystem{D,ChemicalElement,AT} -atomic_symbol(sys::AbstractAtomicSystem) = atomic_symbol.(sys) -atomic_number(sys::AbstractAtomicSystem) = atomic_number.(sys) -atomic_mass(sys::AbstractAtomicSystem) = atomic_mass.(sys) -atomic_property(sys::AbstractAtomicSystem, property::Symbol)::Vector{Any} = - atomic_property.(sys, property) -atomic_propertiesnames(sys::AbstractAtomicSystem) = unique(sort(atomic_propertynames.(sys))) - -struct SimpleAtom{D} <: AbstractAtom - position::SVector{D,<:Unitful.Length} - element::ChemicalElement -end -SimpleAtom(position, element) = SimpleAtom{length(position)}(position, element) -position(atom::SimpleAtom) = atom.position -element(atom::SimpleAtom) = atom.element +""" + species(sys::AbstractSystem{D,S}) + species(sys::AbstractSystem, index) -function SimpleAtom(position, symbol::Union{Integer,AbstractString,Symbol,AbstractVector}) - SimpleAtom(position, ChemicalElement(symbol)) -end +Return a vector of species of every particle in the system `sys`. Return type should be a vector of length `D` containing elements of type `S`. If an index is passed, return only the species of the particle at that index. +""" +species(sys::AbstractSystem) = species.(sys) +(species(sys::AbstractSystem{D,S}, index)::S) where {D,S} = species(sys[index]) # Just to make testing a little easier for now -function Base.show(io::IO, ::MIME"text/plain", part::AbstractParticle) - print(io, "Particle(", element(part), ") @ ", position(part)) -end -function Base.show(io::IO, ::MIME"text/plain", part::AbstractAtom) - print(io, "Atom(", atomic_symbol(part), ") @ ", position(part)) -end function Base.show(io::IO, mime::MIME"text/plain", sys::AbstractSystem) - println(io, "System:") + println(io, "$(string(nameof(typeof(sys)))):") println(io, " BCs: ", boundary_conditions(sys)) println(io, " Box: ", bounding_box(sys)) println(io, " Particles: ") for particle in sys + print(" ") Base.show(io, mime, particle) println(io) end diff --git a/test/runtests.jl b/test/runtests.jl index db0aff9..f7af318 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,43 @@ using AtomsBase using Test +using StaticArrays +using Unitful +using PeriodicTable @testset "AtomsBase.jl" begin - # Write your tests here. + # Basically transcribing aos_vs_soa example for now... + box = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]u"m" + bcs = [Periodic(), Periodic(), DirichletZero()] + positions = [0.25 0.25 0.25; 0.75 0.75 0.75]u"m" + elems = [elements[:C], elements[:C]] + + atom1 = StaticAtom(SVector{3}(positions[1,:]),elems[1]) + atom2 = StaticAtom(SVector{3}(positions[2,:]),elems[2]) + + aos = FlexibleSystem(box, bcs, [atom1, atom2]) + soa = FastSystem(box, bcs, positions, elems) + + @testset "Atoms" begin + @test position(atom1) == [0.25, 0.25, 0.25]u"m" + @test ismissing(velocity(atom1)) + @test species(atom1) == elements[:C] + @test atomic_symbol(atom1) == "C" + @test atomic_number(atom1) == 6 + @test atomic_mass(atom1) == 12.011u"u" + @test atomic_property(atom1, :shells) == [2, 4] + end + + @testset "System" begin + @test bounding_box(aos) == [[1, 0, 0], [0, 1, 0], [0, 0, 1]]u"m" + @test boundary_conditions(aos) == [Periodic(), Periodic(), DirichletZero()] + @test is_periodic(aos) == [1, 1, 0] + @test n_dimensions(aos) == 3 + @test position(aos) == [[0.25, 0.25, 0.25], [0.75, 0.75, 0.75]]u"m" + @test position(aos, 1) == [0.25, 0.25, 0.25]u"m" + @test all(ismissing.(velocity(aos))) + @test species(aos) == [elements[:C], elements[:C]] + @test species(soa, 1) == elements[:C] + @test all(soa .== aos) + end + end