Add AtomsBase conversion support (#76)

Project.toml

@@ -4,16 +4,26 @@ authors = ["Guillaume Fraux <guillaume.fraux@epfl.ch>"]
 version = "0.10.40"
 
 [deps]
+AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
 Chemfiles_jll = "78a364fa-1a3c-552a-b4bb-8fa0f9c1fcca"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+PeriodicTable = "7b2266bf-644c-5ea3-82d8-af4bbd25a884"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [compat]
+AtomsBase = "0.3"
+AtomsBaseTesting = "0.1"
 Chemfiles_jll = "= 0.10.4"
 DocStringExtensions = "0.8.3, 0.9"
+PeriodicTable = "1"
+Unitful = "1"
+UnitfulAtomic = "1"
 julia = "1.6"
 
 [extras]
+AtomsBaseTesting = "ed7c10db-df7e-4efa-a7be-4f4190f7f227"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["AtomsBaseTesting", "Test"]

docs/make.jl

@@ -47,6 +47,7 @@ using Pkg
 Pkg.activate(joinpath(@__DIR__, ".."))
 Pkg.instantiate()
 using Chemfiles
+using AtomsBase
 
 makedocs(
     sitename="Chemfiles.jl",
@@ -67,5 +68,6 @@ makedocs(
         "reference/residue.md",
         "reference/selection.md",
         "reference/misc.md",
+        "reference/atomsbase.md"
     ]
 )

docs/src/index.md

@@ -7,9 +7,13 @@ The Julia interface to chemfiles wraps around the C interface providing a Julian
 API. All the functionalities are in the `Chemfiles` module, which can be
 imported by the `using Chemfiles` expression. The `Chemfiles` module is built
 around the main types of chemfiles: [`Trajectory`](@ref), [`Frame`](@ref),
-[`UnitCell`](@ref), [`Topology`](@ref), [`Residue`](@ref), [`Atom`](@ref), and
+[`UnitCell`](@ref), [`Topology`](@ref), [`Residue`](@ref), [`Chemfiles.Atom`](@ref), and
 [`Selection`](@ref).
 
+Note that the integrates with the [AtomsBase](https://github.com/JuliaMolSim/AtomsBase.jl)
+ecosystem via appropriate conversion routines. See [AtomsBase integration](@ref)
+for more details.
+
 !!! warning
 
     All indexing in chemfiles is 0-based! That means that the first atom in a
@@ -46,5 +50,6 @@ Pages = [
     "reference/residue.md",
     "reference/selection.md",
     "reference/misc.md",
+    "reference/atomsbase.md"
 ]
 ```

docs/src/reference/atom.md

@@ -1,7 +1,7 @@
 # Atom
 
 ```@docs
-Atom
+Chemfiles.Atom
 ```
 
 ```@autodocs

docs/src/reference/atomsbase.md

@@ -0,0 +1,17 @@
+# AtomsBase integration
+
+`Chemfiles` is integrated with the [AtomsBase](https://github.com/JuliaMolSim/AtomsBase.jl)
+by means of appropriate `convert` methods as well as respective constructors,
+which are listed below. Note that these conversion routines copy data unlike the
+majority of functionality in Chemfiles.
+
+If your main interest is in reading structural data from disk, have a look at
+[AtomsIO](https://github.com/mfherbst/AtomsIO.jl). This package provides a uniform
+interface for reading structural data and exposing them in an `AtomsBase`-compatible
+way. Under the hood the package employs a number of parser packages (including Chemfiles)
+to support reading / writing a wide range of file formats.
+
+```@autodocs
+Modules = [Chemfiles, Base, AtomsBase]
+Pages   = ["atomsbase.jl"]
+```

docs/src/tutorials.md

@@ -72,7 +72,7 @@ package
 ```
 
 Everything starts in a [`Topology`](@ref). This is the class that defines the
-atoms and the connectivity in a system. Here, we add three [`Atom`](@ref) and
+atoms and the connectivity in a system. Here, we add three [`Chemfiles.Atom`](@ref) and
 two bonds to create a water molecule.
 
 ```@literalinclude ../examples/generate.jl 7-13

src/Atom.jl

@@ -1,6 +1,7 @@
 # Chemfiles.jl, a modern library for chemistry file reading and writing
 # Copyright (C) Guillaume Fraux and contributors -- BSD license
 
+import AtomsBase: atomic_number
 export mass, set_mass!, charge, set_charge!, name, set_name!, type , set_type!
 export vdw_radius, covalent_radius, atomic_number
 export set_property!, property, properties_count, list_properties

src/Chemfiles.jl

@@ -131,6 +131,8 @@ module Chemfiles
     include("Selection.jl")
     include("Trajectory.jl")
 
+    include("atomsbase.jl")
+
     function __init__()
         if !startswith(version(), "0.10")
             error(

src/atomsbase.jl

@@ -0,0 +1,164 @@
+import AtomsBase
+using AtomsBase: AbstractSystem, FlexibleSystem
+using Unitful
+using UnitfulAtomic
+import PeriodicTable
+
+
+"""
+Construct an AtomsBase `FlexibleSystem` from a Chemfiles Frame.
+"""
+function AtomsBase.FlexibleSystem(frame::Chemfiles.Frame)
+    atoms = map(1:length(frame), frame) do i, atom
+        pos = Chemfiles.positions(frame)[:, i]u"Å"
+        if Chemfiles.has_velocities(frame)
+            velocity = Chemfiles.velocities(frame)[:, i]u"Å/ps"
+        else
+            velocity = zeros(3)u"Å/ps"
+        end
+
+        # Collect atomic properties
+        atprops = Dict(
+            :atomic_mass     => Chemfiles.mass(atom)u"u",
+            :atomic_symbol   => Symbol(Chemfiles.name(atom)),
+            :atomic_number   => Chemfiles.atomic_number(atom),
+            :charge          => Chemfiles.charge(atom)u"e_au",
+            :covalent_radius => Chemfiles.covalent_radius(atom)u"Å",
+            :vdw_radius      => Chemfiles.vdw_radius(atom)*u"Å",
+        )
+        for prop in Chemfiles.list_properties(atom)
+            symbol = Symbol(prop)
+            if !hasfield(AtomsBase.Atom, symbol) && !(symbol in keys(atprops))
+                atprops[symbol] = Chemfiles.property(atom, prop)
+            end
+        end
+
+        AtomsBase.Atom(Chemfiles.atomic_number(atom), pos, velocity; atprops...)
+    end
+
+    # Collect system properties
+    sysprops = Dict{Symbol,Any}()
+    for prop in Chemfiles.list_properties(frame)
+        if hasfield(FlexibleSystem, Symbol(prop))
+            continue
+        elseif prop == "charge"
+            value = Chemfiles.property(frame, "charge")
+            value isa AbstractString && (value = parse(Float64, value))  # Work around a bug
+            sysprops[:charge] = Float64(value)u"e_au"
+        elseif prop == "multiplicity"
+            value = Chemfiles.property(frame, "multiplicity")
+            value isa AbstractString && (value = parse(Float64, value))  # Work around a bug
+            sysprops[:multiplicity] = Int(value)
+        else
+            sysprops[Symbol(prop)] = Chemfiles.property(frame, prop)
+        end
+    end
+
+    # Construct flexible system
+    cell_shape = Chemfiles.shape(Chemfiles.UnitCell(frame))
+    if cell_shape == Chemfiles.Infinite
+        AtomsBase.isolated_system(atoms; sysprops...)
+    else
+        @assert cell_shape in (Chemfiles.Triclinic, Chemfiles.Orthorhombic)
+        box = collect(eachrow(Chemfiles.matrix(Chemfiles.UnitCell(frame))))u"Å"
+        AtomsBase.periodic_system(atoms, box; sysprops...)
+    end
+end
+
+"""
+Construct an AtomsBase `AbstractSystem` from a Chemfiles Frame.
+"""
+AtomsBase.AbstractSystem(frame::Chemfiles.Frame) = FlexibleSystem(frame)
+
+"""
+Convert a Chemfiles Frame to an AtomsBase `AbstractSystem`
+"""
+Base.convert(::Type{AbstractSystem}, frame::Chemfiles.Frame) = AbstractSystem(frame)
+
+"""
+Convert a Chemfiles Frame to an AtomsBase `FlexibleSystem`
+"""
+Base.convert(::Type{FlexibleSystem}, frame::Chemfiles.Frame) = FlexibleSystem(frame)
+
+
+function Base.convert(::Type{Frame}, system::AbstractSystem{D}) where {D}
+    D != 3 && @warn "1D and 2D systems not yet fully supported."
+    frame = Chemfiles.Frame()
+
+    # Cell and boundary conditions
+    if AtomsBase.bounding_box(system) == AtomsBase.infinite_box(D)  # System is infinite
+        cell = Chemfiles.UnitCell(zeros(3, 3))
+        Chemfiles.set_shape!(cell, Chemfiles.Infinite)
+        Chemfiles.set_cell!(frame, cell)
+    else
+        if any(!isequal(AtomsBase.Periodic()), AtomsBase.boundary_conditions(system))
+            @warn("Ignoring specified boundary conditions: Chemfiles only supports " *
+                  "infinite or all-periodic boundary conditions.")
+        end
+
+        box = zeros(3, 3)
+        for i = 1:D
+            box[i, 1:D] = ustrip.(u"Å", AtomsBase.bounding_box(system)[i])
+        end
+        cell = Chemfiles.UnitCell(box)
+        Chemfiles.set_cell!(frame, cell)
+    end
+
+    if any(atom -> !ismissing(AtomsBase.velocity(atom)), system)
+        Chemfiles.add_velocities!(frame)
+    end
+    for atom in system
+        # We are using the atomic_number here, since in AtomsBase the atomic_symbol
+        # can be more elaborate (e.g. D instead of H or "¹⁸O" instead of just "O").
+        # In Chemfiles this is solved using the "name" of an atom ... to which we
+        # map the AtomsBase.atomic_symbol.
+        cf_atom = Chemfiles.Atom(PeriodicTable.elements[atomic_number(atom)].symbol)
+        Chemfiles.set_name!(cf_atom, string(AtomsBase.atomic_symbol(atom)))
+        Chemfiles.set_mass!(cf_atom, ustrip(u"u", AtomsBase.atomic_mass(atom)))
+        @assert Chemfiles.atomic_number(cf_atom) == AtomsBase.atomic_number(atom)
+
+        for (k, v) in pairs(atom)
+            if k in (:atomic_symbol, :atomic_number, :atomic_mass, :velocity, :position)
+                continue  # Dealt with separately
+            elseif k == :charge
+                Chemfiles.set_charge!(cf_atom, ustrip(u"e_au", v))
+            elseif k == :vdw_radius
+                if v != Chemfiles.vdw_radius(cf_atom)u"Å"
+                    @warn "Atom vdw_radius in Chemfiles cannot be mutated"
+                end
+            elseif k == :covalent_radius
+                if v != Chemfiles.covalent_radius(cf_atom)u"Å"
+                    @warn "Atom covalent_radius in Chemfiles cannot be mutated"
+                end
+            elseif v isa Union{Bool, Float64, String, Vector{Float64}}
+                Chemfiles.set_property!(cf_atom, string(k), v)
+            else
+                @warn "Ignoring unsupported property type $(typeof(v)) in Chemfiles for key $k"
+            end
+        end
+
+        pos = convert(Vector{Float64}, ustrip.(u"Å", position(atom)))
+        if ismissing(AtomsBase.velocity(atom))
+            Chemfiles.add_atom!(frame, cf_atom, pos)
+        else
+            vel = convert(Vector{Float64}, ustrip.(u"Å/ps", AtomsBase.velocity(atom)))
+            Chemfiles.add_atom!(frame, cf_atom, pos, vel)
+        end
+    end
+
+    for (k, v) in pairs(system)
+        if k in (:bounding_box, :boundary_conditions)
+            continue  # Already dealt with
+        elseif k in (:charge, )
+            Chemfiles.set_property!(frame, string(k), Float64(ustrip(u"e_au", v)))
+        elseif k in (:multiplicity, )
+            Chemfiles.set_property!(frame, string(k), Float64(v))
+        elseif v isa Union{Bool, Float64, String, Vector{Float64}}
+            Chemfiles.set_property!(frame, string(k), v)
+        else
+            @warn "Ignoring unsupported property type $(typeof(v)) in Chemfiles for key $k"
+        end
+    end
+
+    frame
+end

test/atomsbase.jl

@@ -0,0 +1,92 @@
+using Unitful
+using UnitfulAtomic
+
+@testset "AtomsBase support" begin
+    import AtomsBase
+    using AtomsBase: AbstractSystem, FlexibleSystem
+    using AtomsBaseTesting
+
+    function make_chemfiles_system(D=3; drop_atprop=Symbol[], infinite=false, kwargs...)
+        dropkeys = [:covalent_radius, :vdw_radius]  # Cannot be mutated in Chemfiles
+        data = make_test_system(D; drop_atprop=append!(drop_atprop, dropkeys),
+                                extra_sysprop=(; extra_data=42.0),
+                                cellmatrix=:upper_triangular,
+                                kwargs...)
+        if infinite
+            system = AtomsBase.isolated_system(data.atoms; data.sysprop...)
+        else
+            system = AtomsBase.periodic_system(data.atoms, data.box; data.sysprop...)
+        end
+        merge(data, (; system))
+    end
+
+    @testset "Conversion AtomsBase -> Chemfiles (periodic, velocity)" begin
+        system, atoms, atprop, sysprop, box, bcs = make_chemfiles_system(; infinite=false)
+        frame = convert(Frame, system)
+
+        D = 3
+        cell = Chemfiles.matrix(Chemfiles.UnitCell(frame))
+        for i in 1:D
+            @test cell[i, :] ≈ ustrip.(u"Å", box[i]) atol=1e-12
+        end
+        @test Chemfiles.shape(Chemfiles.UnitCell(frame)) in (Chemfiles.Triclinic,
+                                                             Chemfiles.Orthorhombic)
+        for (i, atom) in enumerate(frame)
+            @test(Chemfiles.positions(frame)[:, i]
+                  ≈ ustrip.(u"Å", atprop.position[i]), atol=1e-12)
+            @test(Chemfiles.velocities(frame)[:, i]
+                  ≈ ustrip.(u"Å/ps", atprop.velocity[i]), atol=1e-12)
+
+            @test Chemfiles.name(atom)            == string(atprop.atomic_symbol[i])
+            @test Chemfiles.atomic_number(atom)   == atprop.atomic_number[i]
+            @test Chemfiles.mass(atom)            == ustrip(u"u", atprop.atomic_mass[i])
+            @test Chemfiles.charge(atom)          == ustrip(u"e_au", atprop.charge[i])
+            @test Chemfiles.list_properties(atom) == ["magnetic_moment"]
+            @test Chemfiles.property(atom, "magnetic_moment") == atprop.magnetic_moment[i]
+
+            if atprop.atomic_number[i] == 1
+                @test Chemfiles.vdw_radius(atom)      == 1.2
+                @test Chemfiles.covalent_radius(atom) == 0.37
+            end
+        end
+
+        @test Chemfiles.list_properties(frame)    == ["charge", "extra_data", "multiplicity"]
+        @test Chemfiles.property(frame, "charge") == ustrip(u"e_au", sysprop.charge)
+        @test Chemfiles.property(frame, "extra_data") == sysprop.extra_data
+        @test Chemfiles.property(frame, "multiplicity") == sysprop.multiplicity
+    end
+
+    @testset "Conversion AtomsBase -> Chemfiles (infinite, no velocity)" begin
+        system = make_chemfiles_system(; infinite=true, drop_atprop=[:velocity]).system
+        frame  = convert(Frame, system)
+        @test Chemfiles.shape(Chemfiles.UnitCell(frame)) == Chemfiles.Infinite
+        @test iszero(Chemfiles.velocities(frame))
+    end
+
+    @testset "Warning about setting invalid data" begin
+        system, atoms, atprop, sysprop, box, bcs = make_test_system()
+        frame  = @test_logs((:warn, r"Atom vdw_radius in Chemfiles cannot be mutated"),
+                            (:warn, r"Atom covalent_radius in Chemfiles cannot be mutated"),
+                            (:warn, r"Ignoring unsupported property type Int[0-9]+.*key extra_data"),
+                            (:warn, r"Ignoring specified boundary conditions:"),
+                            match_mode=:any, convert(Frame, system))
+
+        D = 3
+        cell = Chemfiles.matrix(Chemfiles.UnitCell(frame))
+        for i in 1:D
+            @test cell[i, :] ≈ ustrip.(u"Å", box[i]) atol=1e-12
+        end
+        @test Chemfiles.shape(Chemfiles.UnitCell(frame)) in (Chemfiles.Triclinic,
+                                                             Chemfiles.Orthorhombic)
+    end
+
+    @testset "Conversion AtomsBase -> Chemfiles -> AtomsBase" begin
+        import AtomsBase
+
+        system    = make_chemfiles_system().system
+        frame     = convert(Frame, system)
+        newsystem = convert(AtomsBase.AbstractSystem, frame)
+        test_approx_eq(system, newsystem;
+                       rtol=1e-12, ignore_atprop=[:covalent_radius, :vdw_radius])
+    end
+end

test/runtests.jl

@@ -11,6 +11,7 @@ TESTS = [
     "trajectory.jl",
     "selection.jl",
     "misc.jl",
+    "atomsbase.jl",
 ]
 
 include("utils.jl")