Add AtomsBase conversion support

Project.toml

@@ -4,16 +4,22 @@ 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]
+AtomsBaseTesting = "0.1"
 Chemfiles_jll = "= 0.10.4"
 DocStringExtensions = "0.8.3, 0.9"
 julia = "1.6"
 
 [extras]
+AtomsBaseTesting = "ed7c10db-df7e-4efa-a7be-4f4190f7f227"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["AtomsBaseTesting", "Test"]

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,150 @@
+import AtomsBase
+using AtomsBase: AbstractSystem, FlexibleSystem
+using Unitful
+using UnitfulAtomic
+import PeriodicTable
+
+
+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
+AtomsBase.AbstractSystem(frame::Chemfiles.Frame) = FlexibleSystem(frame)
+
+Base.convert(::Type{AbstractSystem}, frame::Chemfiles.Frame) = AbstractSystem(frame)
+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 Int64.*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")