diff --git a/Project.toml b/Project.toml index 2277a9c..d58d34a 100644 --- a/Project.toml +++ b/Project.toml @@ -4,16 +4,26 @@ authors = ["Guillaume Fraux "] 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"] diff --git a/docs/make.jl b/docs/make.jl index 22edad2..222fc30 100644 --- a/docs/make.jl +++ b/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" ] ) diff --git a/docs/src/index.md b/docs/src/index.md index 34bfe63..03b782d 100644 --- a/docs/src/index.md +++ b/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" ] ``` diff --git a/docs/src/reference/atom.md b/docs/src/reference/atom.md index d1c856c..920e398 100644 --- a/docs/src/reference/atom.md +++ b/docs/src/reference/atom.md @@ -1,7 +1,7 @@ # Atom ```@docs -Atom +Chemfiles.Atom ``` ```@autodocs diff --git a/docs/src/reference/atomsbase.md b/docs/src/reference/atomsbase.md new file mode 100644 index 0000000..5039ef3 --- /dev/null +++ b/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"] +``` diff --git a/docs/src/tutorials.md b/docs/src/tutorials.md index 887061d..5590007 100644 --- a/docs/src/tutorials.md +++ b/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 diff --git a/src/Atom.jl b/src/Atom.jl index 5475a24..57095f1 100644 --- a/src/Atom.jl +++ b/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 diff --git a/src/Chemfiles.jl b/src/Chemfiles.jl index 8843f95..ec8ec58 100644 --- a/src/Chemfiles.jl +++ b/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( diff --git a/src/atomsbase.jl b/src/atomsbase.jl new file mode 100644 index 0000000..ec4b0ff --- /dev/null +++ b/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 diff --git a/test/atomsbase.jl b/test/atomsbase.jl new file mode 100644 index 0000000..b049013 --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index aae40d4..ca0a6d4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,7 @@ TESTS = [ "trajectory.jl", "selection.jl", "misc.jl", + "atomsbase.jl", ] include("utils.jl")