diff --git a/Project.toml b/Project.toml
index 3936b61..f94ab9b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "AtomsBase"
 uuid = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
 authors = ["JuliaMolSim community"]
-version = "0.1.1"
+version = "0.2.0"
 
 [deps]
 PeriodicTable = "7b2266bf-644c-5ea3-82d8-af4bbd25a884"
diff --git a/README.md b/README.md
index 01f9dc2..94aed9a 100644
--- a/README.md
+++ b/README.md
@@ -14,12 +14,25 @@ AtomsBase is an abstract interface for representation of atomic geometries in Ju
 * 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 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.
+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 flexible implementations of the interface](#atomic-systems)
+that we envision to be broadly applicable.
+If features beyond these are required we
+we encourage developers to open PRs or provide their own implementations.
 
 ## Overview
-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.
+The main abstract type introduced in AtomsBase `AbstractSystem{D}`. The `D`
+parameter indicates the number of spatial dimensions in the system.
+Contained inside the system are species, which may have an arbitrary type,
+accessible via the `species_type(system)` function.
+While AtomsBase provides some default species types (e.g. `Atom` and `AtomView`
+for standard atoms) in principle no constraints are made on this species type.
 
-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.
+The main power of the interface comes from predictable behavior of several core
+functions to access information about a system and the species.
+Various categories of such functions are described below.
 
 ### System geometry
 Functions that need to be dispatched:
@@ -27,42 +40,144 @@ Functions that need to be dispatched:
 * `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
 
 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
+* `periodicity`: 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
 
 ### 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.
+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.
 
 ### 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.
-
-As a concrete example, AtomsBase provides the `StaticAtom` type as this is anticipated to be a commonly needed representation. Its implementation looks as follows:
+The only required properties to be specified of the system is the species
+and implementations of standard functions accessing the properties of the species,
+currently `position`, `velocity`, `atomic_symbol`, `atomic_mass`, `atomic_number`,
+`n_dimensions`, `element`. Based on these methods respective equivalent methods acting
+on an `AbstractSystem` will be automatically available, e.g. using the iteration
+interface of the `AbstractSystem` (see above). Most of the property accessors on the
+`AbstractSystem` also have indexed signatures to extract data from a particular species
+directly, for example:
 ```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
+position(sys, i) # position of `i`th particle in `sys`
 ```
-Note that the default behavior of `velocity` is to return `missing`, so it doesn't need to be explicitly dispatched here.
+Currently, this syntax only supports linear indexing.
 
-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.
+To simplify working with `AtomsBase`, default implementations for systems
+composed of atoms are provided, [see below](#atomic-systems)
 
-The `position`, `velocity`, and `species` functions also have indexed signatures to extract a given element directly, as in, for example:
+### 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.
+
+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 concrete implementations see the section on [atomic systems](#atomic-systems) below.
 
+### 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.
+
+## Atomic systems
+Since we anticipate atomic systems to be a commonly needed representation,
+`AtomsBase` provides two flexible implementations for this setting.
+One implementation follows the struct-of-arrays approach introducing the `AtomView`
+type to conveniently expose atomic data.
+The more flexible implementation is based on an array-of-structs approach
+and can be easily customised, e.g. by adding custom properties or by swapping
+the underlying `Atom` struct by a custom one.
+In both cases the respective datastructures can be used either fully
+or in parts in downstream packages and we hope these to develop into universally
+useful types within the Julia ecosystem over time.
+
+### Struct of Arrays / FastSystem
+The file [src/fast_system.jl](src/fast_system.jl) contains an implementation of
+AtomsBase based on the struct-of-arrays approach. All species data is stored
+as plain arrays, but for convenience indexing of individual atoms is supported
+by a light-weight [`AtomView`](src/atomview.jl). See the implementation files
+as well as the tests for how these can be used.
+
+### Atoms and FlexibleSystem
+A flexible implementation of the interface is provided by the
+[`FlexibleSystem`]( src/flexible_system.jl) and the [`Atom`]( src/atom.jl) structs
+for representing atomic systems.
+
+An `Atom` object can be constructed
+just by passing an identifier (e.g. symbol like `:C`, atomic number like `6`) and a vector
+of positions as
 ```julia
-position(sys, i) # position of `i`th particle in `sys`
+atom = Atoms(:C, [0, 1, 2.]u"bohr")
 ```
-Currently, this syntax only supports linear indexing.
-
-### 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.
+This automatically fills the atom with standard data such as the atomic mass. Such data
+can be accessed using the `AtomsBase` interface functions
+such as `atomic_mass(atom)`, `position(atom)`, `velocity(atom)`, `atomic_mass(atom)`, etc.
+See [src/atom.jl](src/atom.jl) for details.
+
+Custom properties can be easily attached to an `Atom` by supplying arbitrary
+keyword arguments upon construction. For example to attach a pseudopotential
+for using the structure with [DFTK](https://dftk.org), construct the atom as
+```julia
+atom = Atoms(:C, [0, 1, 2.]u"bohr", pseudopotential="hgh/lda/c-q4")
+```
+which will make the pseudopotential identifier available as `atom.pseudopotential`.
+Updating an atomic property proceeds similarly. E.g.
+```julia
+newatom = Atoms(;atom=atom, atomic_mass=13u"u")
+```
+makes a new carbon atom with all properties identical to `atom` (including custom ones),
+but setting the `atomic_mass` to 13 units.
 
-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."
+Once the atoms are constructed these can be assembled into a system.
+For example to place a hydrogen molecule into a cubic box of `10Å` and periodic
+boundary conditions, use:
+```julia
+bounding_box = [[10.0, 0.0, 0.0], [0.0, 10.0, 0.0], [0.0, 0.0, 10.0]]u"Å"
+boundary_conditions = [Periodic(), Periodic(), Periodic()]
+hydrogen = FlexibleSystem([Atom(:H, [0, 0, 1.]u"bohr"),
+                           Atom(:H, [0, 0, 3.]u"bohr")],
+                          bounding_box, boundary_conditions)
+```
+An update constructor is supported as well (see [src/flexible_system.jl](src/flexible_system.jl)).
 
-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.
+Oftentimes more convenient are the functions
+`atomic_system`, `isolated_system`, `periodic_system`,
+which cover some standard atomic system setups:
+```julia
+# Same hydrogen system with periodic BCs:
+bounding_box = [[10.0, 0.0, 0.0], [0.0, 10.0, 0.0], [0.0, 0.0, 10.0]]u"Å"
+hydrogen = periodic_system([:H => [0, 0, 1.]u"bohr",
+                            :H => [0, 0, 3.]u"bohr"],
+                           bounding_box)
+
+# Silicon unit cell using fractional positions
+# (Common for solid-state simulations)
+bounding_box = 10.26 / 2 * [[0, 0, 1], [1, 0, 1], [1, 1, 0]]u"bohr"
+silicon = periodic_system([:Si =>  ones(3)/8,
+                           :Si => -ones(3)/8],
+                           bounding_box, fractional=true)
+
+# Isolated H2 molecule in vacuum (Infinite box and zero dirichlet BCs)
+# (Standard setup for molecular simulations)
+hydrogen = isolated_system([:H => [0, 0, 1.]u"bohr",
+                            :H => [0, 0, 3.]u"bohr"])
 
-### 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/Fe_afm.pwi b/sandbox/Fe_afm.pwi
deleted file mode 100644
index 013703e..0000000
--- a/sandbox/Fe_afm.pwi
+++ /dev/null
@@ -1,44 +0,0 @@
-! Slightly modified from
-! https://gitlab.com/QEF/material-for-ljubljana-qe-summer-school/-/raw/master/Day-1/example5.Fe/pw.fe_afm.scf.in
-
- &CONTROL
-    prefix='fe',
-
-    !pseudo_dir = 'directory with pseudopotentials',
-    !outdir = 'temporary directory for large files'
-    !verbosity = 'high',
- /
-
- &SYSTEM
-    ibrav = 1,
-    celldm(1) = 5.42,
-    nat = 2,
-    ntyp = 2,
-    ecutwfc = 25.0,
-    ecutrho = 200.0,
-
-    occupations='smearing',
-    smearing='mv',
-    degauss=0.01,
-
-    nspin=2,
-    starting_magnetization(1) =  0.6
-    starting_magnetization(2) = -0.6
- /
-
- &ELECTRONS
- /
-
-ATOMIC_SPECIES
-# the second field, atomic mass, is not actually used
-# except for MD calculations
-   Fe1  1.  Fe.pbe-nd-rrkjus.UPF
-   Fe2  1.  Fe.pbe-nd-rrkjus.UPF
-
-ATOMIC_POSITIONS crystal
-   Fe1 0.0  0.0  0.0
-   Fe2 0.5  0.5  0.0
-! this is a comment that the code will ignore
-
-K_POINTS automatic
-   8 8 8   1 1 1
diff --git a/sandbox/Manifest.toml b/sandbox/Manifest.toml
deleted file mode 100644
index 2972e47..0000000
--- a/sandbox/Manifest.toml
+++ /dev/null
@@ -1,88 +0,0 @@
-# This file is machine-generated - editing it directly is not advised
-
-[[AtomsBase]]
-deps = ["PeriodicTable", "StaticArrays", "Unitful", "UnitfulAtomic"]
-path = ".."
-uuid = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
-version = "0.1.0"
-
-[[Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-
-[[ConstructionBase]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
-uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-version = "1.3.0"
-
-[[Dates]]
-deps = ["Printf"]
-uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
-
-[[InteractiveUtils]]
-deps = ["Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-
-[[Libdl]]
-uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-
-[[LinearAlgebra]]
-deps = ["Libdl"]
-uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-
-[[Logging]]
-uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-
-[[Markdown]]
-deps = ["Base64"]
-uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-
-[[PeriodicTable]]
-deps = ["Base64", "Test", "Unitful"]
-git-tree-sha1 = "067ef5738ef3157f29749c32b64f0ff9b3d07ab2"
-uuid = "7b2266bf-644c-5ea3-82d8-af4bbd25a884"
-version = "1.1.0"
-
-[[Printf]]
-deps = ["Unicode"]
-uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-
-[[Random]]
-deps = ["Serialization"]
-uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-
-[[Serialization]]
-uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-
-[[SparseArrays]]
-deps = ["LinearAlgebra", "Random"]
-uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-
-[[StaticArrays]]
-deps = ["LinearAlgebra", "Random", "Statistics"]
-git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
-uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.2.12"
-
-[[Statistics]]
-deps = ["LinearAlgebra", "SparseArrays"]
-uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
-[[Test]]
-deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
-uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[[Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
-
-[[Unitful]]
-deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"]
-git-tree-sha1 = "a981a8ef8714cba2fd9780b22fd7a469e7aaf56d"
-uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
-version = "1.9.0"
-
-[[UnitfulAtomic]]
-deps = ["Unitful"]
-git-tree-sha1 = "903be579194534af1c4b4778d1ace676ca042238"
-uuid = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
-version = "1.0.0"
diff --git a/sandbox/Project.toml b/sandbox/Project.toml
deleted file mode 100644
index f51b2c0..0000000
--- a/sandbox/Project.toml
+++ /dev/null
@@ -1,2 +0,0 @@
-[deps]
-AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
diff --git a/sandbox/aos_vs_soa.jl b/sandbox/aos_vs_soa.jl
deleted file mode 100644
index 63f26ba..0000000
--- a/sandbox/aos_vs_soa.jl
+++ /dev/null
@@ -1,21 +0,0 @@
-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 = [periodic_table[:C], periodic_table[:C]]
-
-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)
-
-# And now we can ask questions like...
-soa .== aos
diff --git a/sandbox/example.jl b/sandbox/example.jl
deleted file mode 100644
index 8a4677f..0000000
--- a/sandbox/example.jl
+++ /dev/null
@@ -1,11 +0,0 @@
-include("load_using_ase.jl")
-include("supercell.jl")
-
-iron  = load_using_ase("./Fe_afm.pwi")
-super = make_supercell(iron, (2, 1, 1))
-
-display(iron)
-println()
-println()
-println()
-display(super)
diff --git a/sandbox/load_using_ase.jl b/sandbox/load_using_ase.jl
deleted file mode 100644
index 8fa1fe6..0000000
--- a/sandbox/load_using_ase.jl
+++ /dev/null
@@ -1,22 +0,0 @@
-using PyCall
-using Unitful
-using UnitfulAtomic
-using AtomsBase
-
-# Convert ase Atoms to 3D SimpleAtomicSystem ... could do this via the interface lazily later
-function ase_to_simple(atoms)
-    box = [vec * 1u"Å" for vec in atoms.cell]  # Check this picks up the vectors in the right order
-    boundary_conditions = [(isperiodic ? Periodic() : DirichletZero()) for isperiodic in atoms.pbc]
-
-    pos = atoms.positions * 1u"Å"
-    symbols = atoms.get_chemical_symbols()
-    simple_atoms = SimpleAtom.(eachrow(pos), Symbol.(symbols))
-
-    return SimpleSystem(box, boundary_conditions, simple_atoms)
-end
-
-load_using_ase(filename; kwargs...) = ase_to_simple(pyimport("ase.io").read(filename; kwargs...))
-
-# Field for fractional
-# Field for pseudos
-# Field for magnetic moments
diff --git a/sandbox/supercell.jl b/sandbox/supercell.jl
deleted file mode 100644
index fdf1c39..0000000
--- a/sandbox/supercell.jl
+++ /dev/null
@@ -1,17 +0,0 @@
-using AtomsBase
-
-function make_supercell(system::AbstractSystem{<: AbstractAtom}, repeat::Tuple{<:Integer, <:Integer, <:Integer})
-    # TODO This destroys structure as it only works for SimpleAtomicSystem objects
-    #      ... should ideally be made that it works for all AbstractSystem objects
-
-    newbox = [r * system.box[i] for (i, r) in enumerate(repeat)]
-    @assert n_dimensions(system) == 3
-    newparticles = SimpleAtom{3}[]
-    for part in system
-        for i in 1:repeat[1], j in 1:repeat[2], k in 1:repeat[3]
-            shift = (i-1) * get_box(system)[1] + (j-1) * get_box(system)[2] + (k-1) * get_box(system)[3]
-            push!(newparticles, SimpleAtom{3}(shift + get_position(part), get_element(part)))
-        end
-    end
-    SimpleAtomicSystem{3}(newbox, get_boundary_conditions(system), newparticles)
-end
diff --git a/src/AtomsBase.jl b/src/AtomsBase.jl
index c7d150c..74c7fad 100644
--- a/src/AtomsBase.jl
+++ b/src/AtomsBase.jl
@@ -1,8 +1,13 @@
 module AtomsBase
+using Unitful
+using UnitfulAtomic
+import PeriodicTable: elements
+using StaticArrays
 
 include("interface.jl")
-include("atoms.jl")
-include("implementation_soa.jl")
-include("implementation_aos.jl")
+include("flexible_system.jl")
+include("atom.jl")
+include("atomview.jl")
+include("fast_system.jl")
 
 end
diff --git a/src/atom.jl b/src/atom.jl
new file mode 100644
index 0000000..816eb69
--- /dev/null
+++ b/src/atom.jl
@@ -0,0 +1,81 @@
+#
+# A simple and flexible atom implementation
+#
+export Atom, atomic_system, periodic_system, isolated_system
+
+struct Atom{D, L<:Unitful.Length, V<:Unitful.Velocity, M<:Unitful.Mass}
+    position::SVector{D, L}
+    velocity::SVector{D, V}
+    atomic_symbol::Symbol
+    atomic_number::Int
+    atomic_mass::M
+    data::Dict{Symbol, Any}  # Store arbitrary data about the atom.
+end
+velocity(atom::Atom)      = atom.velocity
+position(atom::Atom)      = atom.position
+atomic_mass(atom::Atom)   = atom.atomic_mass
+atomic_symbol(atom::Atom) = atom.atomic_symbol
+atomic_number(atom::Atom) = atom.atomic_number
+element(atom::Atom)       = elements[atomic_symbol(atom)]
+n_dimensions(atom::Atom{D}) where {D} = D
+
+Base.hasproperty(at::Atom, x::Symbol) = hasfield(Atom, x) || haskey(at.data, x)
+Base.getproperty(at::Atom, x::Symbol) = hasfield(Atom, x) ? getfield(at, x) : getindex(at.data, x)
+function Base.propertynames(at::Atom, private::Bool=false)
+    if private
+        (fieldnames(Atom)..., keys(at.data)...)
+    else
+        (filter(!isequal(:data), fieldnames(Atom))..., keys(at.data)...)
+    end
+end
+
+function Atom(identifier::Union{Symbol,AbstractString,Integer},
+              position::AbstractVector{L},
+              velocity::AbstractVector{V}=zeros(3)u"bohr/s";
+              atomic_symbol=Symbol(elements[identifier].symbol),
+              atomic_number=elements[identifier].number,
+              atomic_mass::M=elements[identifier].atomic_mass,
+              kwargs...) where {L <: Unitful.Length, V <: Unitful.Velocity, M <: Unitful.Mass}
+    Atom{length(position), L, V, M}(position, velocity, atomic_symbol,
+                                    atomic_number, atomic_mass, Dict(kwargs...))
+end
+Atom(id_pos::Pair; kwargs...)  = Atom(id_pos...; kwargs...)
+
+# Update constructor: Amend any atom by extra data.
+function Atom(;atom, kwargs...)
+    extra = atom isa Atom ? atom.data : (; )
+    Atom(atomic_symbol(atom), position(atom), velocity(atom);
+         atomic_symbol=atomic_symbol(atom),
+         atomic_number=atomic_number(atom),
+         atomic_mass=atomic_mass(atom),
+         extra..., kwargs...)
+end
+Atom(atom::Atom; kwargs...) = Atom(; atom, kwargs...)
+
+#
+# Special high-level functions to construct atomic systems
+#
+atomic_system(atoms::AbstractVector{<:Atom}, box, bcs) = FlexibleSystem(atoms, box, bcs)
+atomic_system(atoms::AbstractVector, box, bcs) = atomic_system(Atom.(atoms), box, bcs)
+
+function isolated_system(atoms::AbstractVector{<:Atom})
+    # Use dummy box and boundary conditions
+    D = n_dimensions(first(atoms))
+    atomic_system(atoms, infinite_box(D), fill(DirichletZero(), D))
+end
+isolated_system(atoms::AbstractVector) = isolated_system(Atom.(atoms))
+
+function periodic_system(atoms::AbstractVector,
+                         box::AbstractVector{<:AbstractVector},
+                         boundary_conditions=fill(Periodic(), length(box));
+                         fractional=false, kwargs...)
+    lattice = hcat(box...)
+    !fractional && return atomic_system(atoms, box, boundary_conditions)
+
+    parse_fractional(atom::Atom) = atom
+    function parse_fractional(atom::Pair)
+        id, pos_fractional = atom
+        Atom(id, lattice * pos_fractional)
+    end
+    atomic_system(parse_fractional.(atoms), box, boundary_conditions)
+end
diff --git a/src/atoms.jl b/src/atoms.jl
deleted file mode 100644
index fd3115c..0000000
--- a/src/atoms.jl
+++ /dev/null
@@ -1,36 +0,0 @@
-#
-# 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
-velocity(::StaticAtom) = missing
-
-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/atomview.jl b/src/atomview.jl
new file mode 100644
index 0000000..bd937bf
--- /dev/null
+++ b/src/atomview.jl
@@ -0,0 +1,15 @@
+#
+# A simple view datastructure for atoms of struct of array systems
+#
+export AtomView
+struct AtomView{FS <: AbstractSystem}
+    system::FS
+    index::Int
+end
+velocity(v::AtomView)      = velocity(v.system, v.index)
+position(v::AtomView)      = position(v.system, v.index)
+atomic_mass(v::AtomView)   = atomic_mass(v.system, v.index)
+atomic_symbol(v::AtomView) = atomic_symbol(v.system, v.index)
+atomic_number(v::AtomView) = atomic_number(v.system, v.index)
+element(v::AtomView)       = elements[atomic_symbol(v)]
+n_dimensions(v::AtomView)  = n_dimensions(v.system)
diff --git a/src/fast_system.jl b/src/fast_system.jl
new file mode 100644
index 0000000..8cf91f2
--- /dev/null
+++ b/src/fast_system.jl
@@ -0,0 +1,62 @@
+#
+# Implementation of AtomsBase interface in a struct-of-arrays style.
+#
+export FastSystem
+
+struct FastSystem{D, L <: Unitful.Length, M <: Unitful.Mass} <: AbstractSystem{D}
+    box::SVector{D, SVector{D, L}}
+    boundary_conditions::SVector{D, BoundaryCondition}
+    positions::Vector{SVector{D, L}}
+    atomic_symbols::Vector{Symbol}
+    atomic_numbers::Vector{Int}
+    atomic_masses::Vector{M}
+end
+
+# Constructor to fetch the types
+function FastSystem(box, boundary_conditions, positions, atomic_symbols, atomic_numbers, atomic_masses)
+    FastSystem{length(box),eltype(eltype(positions)),eltype(atomic_masses)}(
+        box, boundary_conditions, positions, atomic_symbols, atomic_numbers, atomic_masses
+    )
+end
+
+# Constructor to take data from another system
+function FastSystem(system::AbstractSystem)
+    FastSystem(bounding_box(system), boundary_conditions(system), position(system),
+               atomic_symbol(system), atomic_number(system), atomic_mass(system))
+end
+
+# Convenience constructor where we don't have to preconstruct all the static stuff...
+function FastSystem(particles, box, boundary_conditions)
+    D = length(box)
+    if !all(length.(box) .== D)
+        throw(ArgumentError("Box must have D vectors of length D=$D."))
+    end
+    if length(boundary_conditions) != D
+        throw(ArgumentError("Boundary conditions should be of length D=$D."))
+    end
+    if !all(n_dimensions.(particles) .== D)
+        throw(ArgumentError("Particles must have positions of length D=$D."))
+    end
+    FastSystem(box, boundary_conditions, position.(particles), atomic_symbol.(particles),
+               atomic_number.(particles), atomic_mass.(particles))
+end
+
+bounding_box(sys::FastSystem)        = sys.box
+boundary_conditions(sys::FastSystem) = sys.boundary_conditions
+Base.length(sys::FastSystem)         = length(sys.positions)
+Base.size(sys::FastSystem)           = size(sys.positions)
+
+species_type(sys::FS) where {FS <: FastSystem} = AtomView{FS}
+Base.getindex(sys::FastSystem, index::Int)     = AtomView(sys, index)
+
+position(s::FastSystem)       = s.positions
+atomic_symbol(s::FastSystem)  = s.atomic_symbols
+atomic_number(s::FastSystem)  = s.atomic_numbers
+atomic_mass(s::FastSystem)    = s.atomic_masses
+velocity(s::FastSystem)       = missing
+
+position(s::FastSystem, i)      = s.positions[i]
+atomic_symbol(s::FastSystem, i) = s.atomic_symbols[i]
+atomic_number(s::FastSystem, i) = s.atomic_numbers[i]
+atomic_mass(s::FastSystem, i)   = s.atomic_masses[i]
+velocity(s::FastSystem, i)      = missing
diff --git a/src/flexible_system.jl b/src/flexible_system.jl
new file mode 100644
index 0000000..4646d9d
--- /dev/null
+++ b/src/flexible_system.jl
@@ -0,0 +1,39 @@
+#
+# Implementation of AtomsBase interface in an array-of-structs style
+#
+export FlexibleSystem
+
+struct FlexibleSystem{D, S, L<:Unitful.Length} <: AbstractSystem{D}
+    particles::AbstractVector{S}
+    box::SVector{D, SVector{D, L}}
+    boundary_conditions::SVector{D, BoundaryCondition}
+end
+
+function FlexibleSystem(
+    particles::AbstractVector{S},
+    box::AbstractVector{<:AbstractVector{L}},
+    boundary_conditions::AbstractVector{BC}
+) where {BC<:BoundaryCondition, L<:Unitful.Length, S}
+    D = length(box)
+    if !all(length.(box) .== D)
+        throw(ArgumentError("Box must have D vectors of length D"))
+    end
+    FlexibleSystem{D, S, L}(particles, box, boundary_conditions)
+end
+
+# Update constructor
+function FlexibleSystem(system::AbstractSystem;
+                        particles=nothing, atoms=nothing,
+                        box=bounding_box(system),
+                        boundary_conditions=boundary_conditions(system))
+    particles = something(particles, atoms, collect(system))
+    FlexibleSystem(particles, box, boundary_conditions)
+end
+
+bounding_box(sys::FlexibleSystem)        = sys.box
+boundary_conditions(sys::FlexibleSystem) = sys.boundary_conditions
+species_type(sys::FlexibleSystem{D, S, L}) where {D, S, L} = S
+
+Base.size(sys::FlexibleSystem)   = size(sys.particles)
+Base.length(sys::FlexibleSystem) = length(sys.particles)
+Base.getindex(sys::FlexibleSystem, i::Integer) = getindex(sys.particles, i)
diff --git a/src/implementation_aos.jl b/src/implementation_aos.jl
deleted file mode 100644
index 45b388b..0000000
--- a/src/implementation_aos.jl
+++ /dev/null
@@ -1,51 +0,0 @@
-#
-# Implementation of AtomsBase interface in an array-of-structs style
-#
-
-using StaticArrays
-
-export FlexibleSystem
-
-# TODO Switch order of type arguments?
-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{<:AbstractVector{L}},
-    boundary_conditions::Vector{BC},
-    particles::Vector{AT},
-) 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,
-    )
-end
-
-bounding_box(sys::FlexibleSystem) = sys.box
-boundary_conditions(sys::FlexibleSystem) = sys.boundary_conditions
-
-function Base.show(io::IO, sys::FlexibleSystem)
-    print(io, "FlexibleSystem with ", length(sys), " particles")
-end
-
-Base.size(sys::FlexibleSystem) = size(sys.particles)
-Base.length(sys::FlexibleSystem) = length(sys.particles)
-Base.getindex(sys::FlexibleSystem, i::Int) = getindex(sys.particles, i)
diff --git a/src/implementation_soa.jl b/src/implementation_soa.jl
deleted file mode 100644
index fd220d8..0000000
--- a/src/implementation_soa.jl
+++ /dev/null
@@ -1,68 +0,0 @@
-#
-# Implementation of AtomsBase interface in a struct-of-arrays style.
-#
-
-using StaticArrays
-
-export FastSystem
-
-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{S}
-    # janky inner constructor that we need for some reason
-    FastSystem(box, boundary_conditions, positions, elements) =
-        new{length(boundary_conditions),eltype(elements),eltype(eltype(positions))}(
-            box,
-            boundary_conditions,
-            positions,
-            elements,
-        )
-end
-
-# convenience constructor where we don't have to preconstruct all the static stuff...
-function FastSystem(
-    box::Vector{<:AbstractVector{L}},
-    boundary_conditions::AbstractVector{BC},
-    positions::AbstractMatrix{M},
-    elements::AbstractVector{S},
-) where {L<:Unitful.Length,BC<:BoundaryCondition,M<:Unitful.Length,S}
-    N = length(elements)
-    D = length(box)
-    if !all(length.(box) .== D)
-        throw(ArgumentError("box must have D vectors of length D"))
-    end
-    if !(size(positions, 1) == N)
-        throw(ArgumentError("box must have D vectors of length D"))
-    end
-    if !(size(positions, 2) == D)
-        throw(ArgumentError("box must have D vectors of length D"))
-    end
-    FastSystem(
-        SVector{D,SVector{D,L}}(box),
-        SVector{D,BoundaryCondition}(boundary_conditions),
-        Vector{SVector{D,eltype(positions)}}([positions[i, :] for i = 1:N]),
-        elements,
-    )
-end
-
-function Base.show(io::IO, sys::FastSystem)
-    print(io, "FastSystem with ", length(sys), " particles")
-end
-
-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,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
-
-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 c9ee32e..cdf3c08 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -1,34 +1,9 @@
-using Unitful
-using UnitfulAtomic
-using PeriodicTable
-using StaticArrays
 import Base.position
 
 export AbstractSystem
-export BoundaryCondition, DirichletZero, Periodic
-export species, position, velocity
-export bounding_box, boundary_conditions, periodicity, n_dimensions
-
-"""
-    velocity(p)
-
-Return the velocity of a particle `p`.
-"""
-function velocity end
-
-"""
-    position(p)
-
-Return the position of a particle `p`.
-"""
-function position end
-
-"""
-    species(p)
-
-Return the species of a particle `p`.
-"""
-function species end
+export BoundaryCondition, DirichletZero, Periodic, infinite_box
+export bounding_box, boundary_conditions, periodicity, n_dimensions, species_type
+export position, velocity, element, atomic_mass, atomic_number, atomic_symbol
 
 #
 # Identifier for boundary conditions per dimension
@@ -37,13 +12,21 @@ abstract type BoundaryCondition end
 struct DirichletZero <: BoundaryCondition end  # Dirichlet zero boundary (i.e. molecular context)
 struct Periodic <: BoundaryCondition end  # Periodic BCs
 
+infinite_box(::Val{1}) = [[Inf]]u"bohr"
+infinite_box(::Val{2}) = [[Inf, 0], [0, Inf]]u"bohr"
+infinite_box(::Val{3}) = [[Inf, 0, 0], [0, Inf, 0], [0, 0, Inf]]u"bohr"
+infinite_box(dim::Int) = infinite_box(Val(dim))
+
 
+#
+# Abstract system
+#
 """
-    AbstractSystem{D,S}
+    AbstractSystem{D}
 
-A `D`-dimensional system comprised of particles identified by type `S`.
+A `D`-dimensional system.
 """
-abstract type AbstractSystem{D,S} end
+abstract type AbstractSystem{D} end
 
 """
     bounding_box(sys::AbstractSystem{D})
@@ -59,45 +42,102 @@ Return a vector of length `D` of `BoundaryCondition` objects, one for each direc
 """
 function boundary_conditions end
 
+"""
+    species_type(::AbstractSystem)
+
+Return the type used to represent a species or atom.
+"""
+function species_type end
+
+"""Return vector indicating whether the system is periodic along a dimension."""
 periodicity(sys::AbstractSystem) = [isa(bc, Periodic) for bc in boundary_conditions(sys)]
 
-# Note: Can't use ndims, because that is ndims(sys) == 1 (because of indexing interface)
+"""
+    n_dimensions(::AbstractSystem)
+    n_dimensions(atom)
+
+Return number of dimensions.
+"""
 n_dimensions(::AbstractSystem{D}) where {D} = D
+# Note: Can't use ndims, because that is ndims(sys) == 1 (because of indexing interface)
 
 # 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)
 # default to 1D indexing
-Base.iterate(S::AbstractSystem, state = firstindex(S)) =
-    (firstindex(S) <= state <= length(S)) ? (@inbounds S[state], state + 1) : nothing
+Base.iterate(sys::AbstractSystem, state = firstindex(sys)) =
+    (firstindex(sys) <= state <= length(sys)) ? (@inbounds sys[state], state + 1) : nothing
 
 # TODO Support similar, push, ...
 
+#
+# Species property accessors from systems and species
+#
+"""The element corresponding to a species/atom (or missing)."""
+function element end
+
+
 """
     position(sys::AbstractSystem{D})
     position(sys::AbstractSystem, index)
+    position(species)
 
-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.
+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 or the action is on a `species`,
+return only the position of the referenced `species` / species on that index.
 """
-position(sys::AbstractSystem) = position.(sys)    # in Cartesian coordinates!
+position(sys::AbstractSystem)        = position.(sys)    # in Cartesian coordinates!
 position(sys::AbstractSystem, index) = position(sys[index])
 
+
 """
     velocity(sys::AbstractSystem{D})
     velocity(sys::AbstractSystem, index)
+    velocity(species)
 
-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.
+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 or the action is on a `species`,
+return only the velocity of the referenced `species`. Returned value of the function
+may be `missing`.
 """
-velocity(sys::AbstractSystem) = velocity.(sys)    # in Cartesian coordinates!
+velocity(sys::AbstractSystem)        = velocity.(sys)    # in Cartesian coordinates!
 velocity(sys::AbstractSystem, index) = velocity(sys[index])
 
+
+"""
+    atomic_mass(sys::AbstractSystem)
+    atomic_mass(sys::AbstractSystem, i)
+    atomic_mass(species)
+
+Vector of atomic masses in the system `sys` or the atomic mass of a particular `species` /
+the `i`th species in `sys`. The elements are `<: Unitful.Mass`.
+"""
+atomic_mass(sys::AbstractSystem)        = atomic_mass.(sys)
+atomic_mass(sys::AbstractSystem, index) = atomic_mass(sys[index])
+
+
+"""
+    atomic_symbol(sys::AbstractSystem)
+    atomic_symbol(sys::AbstractSystem, i)
+    atomic_symbol(species)
+
+Vector of atomic symbols in the system `sys` or the atomic symbol of a particular `species` /
+the `i`th species in `sys`.
+"""
+atomic_symbol(sys::AbstractSystem)        = atomic_symbol.(sys)
+atomic_symbol(sys::AbstractSystem, index) = atomic_symbol(sys[index])
+
+
 """
-    species(sys::AbstractSystem{D,S})
-    species(sys::AbstractSystem, index)
+    atomic_number(sys::AbstractSystem)
+    atomic_number(sys::AbstractSystem, i)
+    atomic_number(species)
 
-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.
+Vector of atomic numbers in the system `sys` or the atomic number of a particular `species` /
+the `i`th species in `sys`.
 """
-species(sys::AbstractSystem) = species.(sys)
-(species(sys::AbstractSystem{D,S}, index)::S) where {D,S} = species(sys[index])
+atomic_number(sys::AbstractSystem)        = atomic_number.(sys)
+atomic_number(sys::AbstractSystem, index) = atomic_number(sys[index])
diff --git a/test/atom.jl b/test/atom.jl
new file mode 100644
index 0000000..841c1bb
--- /dev/null
+++ b/test/atom.jl
@@ -0,0 +1,83 @@
+using AtomsBase
+using Unitful
+using Test
+
+
+@testset "atomic systems" begin
+    @testset "Atom construction" begin
+        at = Atom(:Si, zeros(3) * u"m", extradata=42)
+
+        @test n_dimensions(at) == 3
+        @test position(at) == zeros(3) * u"m"
+        @test velocity(at) == zeros(3) * u"bohr/s"
+        @test at.atomic_symbol == :Si
+        @test at.atomic_number == 14
+        @test hasproperty(at, :atomic_mass)
+        @test hasproperty(at, :atomic_symbol)
+        @test hasproperty(at, :atomic_number)
+        @test hasproperty(at, :extradata)
+        @test at.extradata == 42
+
+        @test propertynames(at) == (:position, :velocity, :atomic_symbol,
+                                    :atomic_number, :atomic_mass, :extradata)
+
+        # Test update constructor
+        newatom = Atom(at; extradata=43, atomic_number=15)
+        @test propertynames(at) == propertynames(newatom)
+        @test newatom.extradata == 43
+        @test newatom.atomic_number == 15
+
+        newatom = Atom(; atom=at, extradata=43, atomic_number=15)
+        @test propertynames(at) == propertynames(newatom)
+        @test newatom.extradata == 43
+        @test newatom.atomic_number == 15
+    end
+
+    @testset "flexible atomic systems" begin
+        box = [[10, 0.0, 0.0], [0.0, 5, 0.0], [0.0, 0.0, 7]]u"Å"
+        bcs = [Periodic(), DirichletZero(), DirichletZero()]
+        atoms = [:Si => [0.0, 1.0, 1.5]u"Å",
+                 :C  => [0.0, 0.8, 1.7]u"Å",
+                 Atom(:H, zeros(3) * u"Å", ones(3) * u"bohr/s")]
+        system = atomic_system(atoms, box, bcs)
+        @test length(system) == 3
+        @test atomic_symbol(system) == [:Si, :C, :H]
+        @test boundary_conditions(system) == [Periodic(), DirichletZero(), DirichletZero()]
+        @test position(system) == [[0.0, 1.0, 1.5], [0.0, 0.8, 1.7], [0.0, 0.0, 0.0]]u"Å"
+        @test velocity(system) == [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]u"bohr/s"
+
+        # Test update constructor
+        newatoms  = [system[1], system[2]]
+        newsystem = FlexibleSystem(system; atoms=newatoms,
+                                   boundary_conditions=[Periodic(), Periodic(), Periodic()])
+        @test length(newsystem) == 2
+        @test atomic_symbol(newsystem) == [:Si, :C]
+        @test boundary_conditions(newsystem) == [Periodic(), Periodic(), Periodic()]
+        @test position(newsystem) == [[0.0, 1.0, 1.5], [0.0, 0.8, 1.7]]u"Å"
+    end
+
+    @testset "isolated_system" begin
+        system = isolated_system([:Si => [0.0, 1.0, 1.5]u"Å",
+                                  :C  => [0.0, 0.8, 1.7]u"Å",
+                                  Atom(:H, zeros(3) * u"Å")])
+        @test length(system) == 3
+        @test atomic_symbol(system) == [:Si, :C, :H]
+        @test boundary_conditions(system) == [DirichletZero(), DirichletZero(), DirichletZero()]
+        @test position(system) == [[0.0, 1.0, 1.5], [0.0, 0.8, 1.7], [0.0, 0.0, 0.0]]u"Å"
+        @test velocity(system) == [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]u"bohr/s"
+        @test bounding_box(system) == infinite_box(3)
+    end
+
+    @testset "periodic_system" begin
+        box = [[10, 0.0, 0.0], [0.0, 5, 0.0], [0.0, 0.0, 7]]u"Å"
+        atoms = [:Si => [0.0, -0.125, 0.0],
+                 :C  => [0.125, 0.0, 0.0],
+                 Atom(:H, zeros(3) * u"Å")]
+        system = periodic_system(atoms, box; fractional=true)
+
+        @test length(system) == 3
+        @test atomic_symbol(system) == [:Si, :C, :H]
+        @test boundary_conditions(system) == [Periodic(), Periodic(), Periodic()]
+        @test position(system) == [[0.0, -0.625, 0.0], [1.25, 0.0, 0.0], [0.0, 0.0, 0.0]]u"Å"
+    end
+end
diff --git a/test/fast_system.jl b/test/fast_system.jl
new file mode 100644
index 0000000..350d80a
--- /dev/null
+++ b/test/fast_system.jl
@@ -0,0 +1,28 @@
+using AtomsBase
+using Test
+using Unitful
+using PeriodicTable
+
+@testset "Fast system" begin
+    box = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]u"m"
+    bcs = [Periodic(), Periodic(), DirichletZero()]
+    atoms = Atom.([:C => [0.25, 0.25, 0.25]u"m",
+                   :C => [0.75, 0.75, 0.75]u"m"])
+    system = FastSystem(atoms, box, bcs)
+
+    @test length(system) == 2
+    @test size(system)   == (2, )
+    @test atomic_mass(system) == [12.011, 12.011]u"u"
+    @test boundary_conditions(system) == bcs
+    @test bounding_box(system) == box
+
+
+    # Test AtomView
+    for method in (position, atomic_mass, atomic_symbol, atomic_number)
+        @test method(system[1]) == method(system, 1)
+        @test method(system[2]) == method(system, 2)
+    end
+    @test ismissing(velocity(system[1]))
+    @test n_dimensions(system[1]) == n_dimensions(system)
+    @test element(system[1])      == elements[:C]
+end
diff --git a/test/interface.jl b/test/interface.jl
new file mode 100644
index 0000000..269a437
--- /dev/null
+++ b/test/interface.jl
@@ -0,0 +1,49 @@
+using AtomsBase
+using Test
+using Unitful
+using UnitfulAtomic
+using PeriodicTable
+
+@testset "Interface" begin
+    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"
+    elements = [:C, :C]
+    atoms = [Atom(elements[i], positions[i]) for i in 1:2]
+
+    @testset "Atoms" begin
+        @test position(atoms[1]) == [0.25, 0.25, 0.25]u"m"
+        @test velocity(atoms[1]) == [0.0, 0.0, 0.0]u"bohr/s"
+        @test element(atoms[1]) == PeriodicTable.elements[:C]
+        @test atomic_symbol(atoms[1]) == :C
+        @test atomic_number(atoms[1]) == 6
+        @test atomic_mass(atoms[1])   == 12.011u"u"
+    end
+
+    @testset "System" begin
+        flexible = FlexibleSystem(atoms, box, bcs)
+        fast     = FastSystem(flexible)
+        @test length(flexible) == 2
+        @test size(flexible)   == (2, )
+
+        @test bounding_box(flexible) == [[1, 0, 0], [0, 1, 0], [0, 0, 1]]u"m"
+        @test boundary_conditions(flexible) == [Periodic(), Periodic(), DirichletZero()]
+        @test periodicity(flexible) == [1, 1, 0]
+        @test n_dimensions(flexible) == 3
+        @test position(flexible) == [[0.25, 0.25, 0.25], [0.75, 0.75, 0.75]]u"m"
+        @test position(flexible, 1) == [0.25, 0.25, 0.25]u"m"
+        @test velocity(flexible)[1] == [0.0, 0.0, 0.0]u"bohr/s"
+        @test velocity(flexible)[2] == [0.0, 0.0, 0.0]u"bohr/s"
+        @test atomic_mass(flexible)   == [12.011, 12.011]u"u"
+        @test atomic_number(fast) == [6, 6]
+        @test atomic_number(fast, 1) == 6
+        @test ismissing(velocity(fast, 2))
+        @test atomic_symbol(flexible, 2) == :C
+        @test atomic_number(flexible, 2) == 6
+        @test atomic_mass(flexible, 1)   == 12.011u"u"
+
+        @test ismissing(velocity(fast))
+        @test all(position(fast)      .== position(flexible))
+        @test all(atomic_symbol(fast) .== atomic_symbol(flexible))
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 4f95e11..d571164 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,43 +1,7 @@
-using AtomsBase
 using Test
-using StaticArrays
-using Unitful
-using PeriodicTable
 
 @testset "AtomsBase.jl" begin
-    # 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 periodicity(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
-
+    include("interface.jl")
+    include("fast_system.jl")
+    include("atom.jl")
 end