Integrate with AtomsBase 0.3 (#14)

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.1.2
+current_version = 0.1.3
 commit = True
 tag = False
 

Project.toml

@@ -1,7 +1,7 @@
 name = "AtomsIO"
 uuid = "1692102d-eeb4-4df9-807b-c9517f998d44"
 authors = ["Michael F. Herbst <[email protected]> and contributors"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 ASEconvert = "3da9722f-58c2-4165-81be-b4d7253e8fd2"
@@ -17,10 +17,10 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [compat]
-ASEconvert = "0.1.2"
-AtomsBase = "0.2"
-Chemfiles = "0.10"
-ExtXYZ = "0.1.8"
+ASEconvert = "0.1.4"
+AtomsBase = "0.3"
+Chemfiles = "0.10.3"
+ExtXYZ = "0.1.11"
 PeriodicTable = "1"
 Reexport = "1"
 StaticArrays = "1"

README.md

@@ -22,10 +22,12 @@ and supports all file formats any of these packages support.
 The best-matching backend for reading / writing is automatically chosen.
 
 Amongst others AtomsIO supports the following formats
+
   - [Crystallographic Information Framework](https://www.iucr.org/resources/cif) (CIF) files
   - [Quantum Espresso](https://www.quantum-espresso.org/Doc/INPUT_PW.html) / [ABINIT](https://docs.abinit.org/variables/) / [VASP](https://www.vasp.at/wiki/) input files
   - ASE / [Gromacs](http://manual.gromacs.org/archive/5.0.7/online/trj.html) / [LAMMPS](https://lammps.sandia.gov/doc/dump.html) / [Amber](http://ambermd.org/netcdf/nctraj.xhtml) trajectory files
   - [XYZ](https://openbabel.org/wiki/XYZ) and [extxyz](https://github.com/libAtoms/extxyz#extended-xyz-specification-and-parsing-tools) files
+
 For more details see the [documentation](https://mfherbst.github.io/AtomsIO.jl/stable).
 
 ## Usage example

docs/src/index.md

@@ -25,13 +25,12 @@ The best-matching backend for reading / writing is automatically chosen.
     reading a file. AtomsIO does this by default.
 
 !!! tip
-    For the reasons mentioned above the authors of AtomsIO recommend using
-    extended XYZ format via the [`ExtxyzParser`](@ref) for long-term storage.
+    For the reasons mentioned above a good long-term storage format
+    (in the eyes of the AtomsIO authors) is the extended XYZ format
+    via the [`ExtxyzParser`](@ref).
     This format and parser has a
     [well-documented specification](https://github.com/libAtoms/extxyz#extended-xyz-specifcation)
     and moreover leads to a human-readable plain-text file.
-    If you have a different opinion we are happy to hear about it.
-    Please open an issue for discussion.
 
 ## Usage example
 

src/chemfiles.jl

@@ -63,12 +63,15 @@ function parse_chemfiles(frame::Chemfiles.Frame)
     # TODO Should be moved to Chemfiles as a
     #      convert(AbstractSystem, ::Chemfile.Frame)
 
+    # Regarding the unit conventions, see:
+    # https://chemfiles.org/chemfiles/latest/overview.html#units
+    #
     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"Å/s"
+            velocity = Chemfiles.velocities(frame)[:, i]u"Å/ps"
         else
-            velocity = zeros(3)u"Å/s"
+            velocity = zeros(3)u"Å/ps"
         end
 
         # Collect atomic properties
@@ -158,52 +161,46 @@ function convert_chemfiles(system::AbstractSystem{D}) where {D}
         Chemfiles.set_mass!(cf_atom, ustrip(u"u", atomic_mass(atom)))
         @assert Chemfiles.atomic_number(cf_atom) == atomic_number(atom)
 
-        if atom isa Atom
-            # TODO not a good idea to directly access the field
-            # TODO Implement and make use of a property interface on the atom level
-            for (k, v) in atom.data
-                if 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 atom key $k")
+        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(velocity(atom))
             Chemfiles.add_atom!(frame, cf_atom, pos)
         else
-            vel = convert(Vector{Float64}, ustrip.(u"Å/s", velocity(atom)))
+            vel = convert(Vector{Float64}, ustrip.(u"Å/ps", velocity(atom)))
             Chemfiles.add_atom!(frame, cf_atom, pos, vel)
         end
     end
 
-    if system isa FlexibleSystem
-        # TODO not a good idea to directly access the field
-        # TODO Implement and make use of a property interface on the system level
-        for (k, v) in system.data
-            if k == :charge
-                Chemfiles.set_property!(frame, string(k), Float64(ustrip(u"e_au", v)))
-            elseif k == :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 system key $k")
-            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
 

src/extxyz.jl

@@ -10,229 +10,21 @@ function supports_parsing(::ExtxyzParser, file; save, trajectory)
     ext in (".xyz", ".extxyz")
 end
 
-
 function load_system(::ExtxyzParser, file::AbstractString, index=nothing)
-    isfile(file) || throw(ArgumentError("File $file not found"))
     frame = (isnothing(index) ? last(ExtXYZ.read_frames(file))
                               : only(ExtXYZ.read_frames(file, index)))
-    parse_extxyz(frame)
+    ExtXYZ.Atoms(frame)
 end
 
 function save_system(::ExtxyzParser, file::AbstractString, system::AbstractSystem)
-    ExtXYZ.write_frame(file, convert_extxyz(system))
+    ExtXYZ.write_frame(file, ExtXYZ.write_dict(system))
 end
 
 function load_trajectory(::ExtxyzParser, file::AbstractString)
-    isfile(file) || throw(ArgumentError("File $file not found"))
-    parse_extxyz.(ExtXYZ.read_frames(file))
+    ExtXYZ.Atoms.(ExtXYZ.read_frames(file))
 end
 
 function save_trajectory(::ExtxyzParser, file::AbstractString,
                          systems::AbstractVector{<:AbstractSystem})
-    ExtXYZ.write_frames(file, convert_extxyz.(systems))
-end
-
-#
-# Convert to/from ExtXYZ.Atoms
-#
-
-function parse_extxyz(dict::Dict{String, Any})
-    # TODO Move this upstream to replace the current implementation of ExtXYZ.Atoms
-    arrays = dict["arrays"]
-    info   = dict["info"]
-
-    if haskey(arrays, "Z")
-        Z = Int.(arrays["Z"])
-    elseif haskey(arrays, "species")
-        Z = [PeriodicTable.elements[Symbol(spec)].number for spec in arrays["species"]]
-    else
-        error("Cannot determine atomic numbers. Either 'Z' or 'species' must " *
-              "be present in arrays")
-    end
-    @assert length(Z) == dict["N_atoms"]
-
-    atom_data = Dict{Symbol, Any}(
-        :positions      => collect(eachcol(arrays["pos"]))u"Å",
-        :atomic_numbers => Z,
-    )
-    if haskey(arrays, "species")
-        atom_data[:atomic_symbols] = Symbol.(arrays["species"])
-    else
-        atom_data[:atomic_symbols] = [Symbol(element(num).symbol) for num in Z]
-    end
-    if haskey(arrays, "mass")
-        atom_data[:atomic_masses] = arrays["mass"]u"u"
-    else
-        atom_data[:atomic_masses] = [PeriodicTable.elements[num].atomic_mass for num in Z]
-    end
-    if haskey(arrays, "velocities")
-        atom_data[:velocities] = collect(eachcol(arrays["velocities"]))u"Å/s"
-    end
-
-    for key in keys(arrays)
-        key in ("mass", "species", "Z", "pos", "velocities") && continue  # Already done
-        if key in ("vdw_radius", "covalent_radius")  # Add length unit
-            atom_data[Symbol(key)] = arrays[key] * u"Å"
-        elseif key in ("charge", )  # Add charge unit
-            atom_data[Symbol(key)] = arrays[key] * u"e_au"
-        else
-            atom_data[Symbol(key)] = arrays[key]
-        end
-    end
-
-    system_data = Dict{Symbol, Any}(
-        :box => collect(eachrow(dict["cell"]))u"Å",
-    )
-    if haskey(dict, "pbc")
-        system_data[:boundary_conditions] = [p ? Periodic() : DirichletZero()
-                                             for p in dict["pbc"]]
-    else
-        @warn "'pbc' not contained in 'info' dict. Defaulting to all-periodic boundary. "
-        system_data[:boundary_conditions] = fill(Periodic(), 3)
-    end
-
-    for key in keys(info)
-        if key in ("charge", )
-            system_data[Symbol(key)] = info[key] * u"e_au"
-        else
-            system_data[Symbol(key)] = info[key]
-        end
-    end
-
-    extxyz_atoms = ExtXYZ.Atoms(NamedTuple(atom_data), NamedTuple(system_data))
-
-    # TODO Once ExtXYZ.Atoms supports atom and system property access,
-    #      this is no longer needed
-    atoms = map(1:length(extxyz_atoms), extxyz_atoms) do i, atom
-        atprops = Dict{Symbol,Any}()
-        for k in keys(atom_data)
-            if !(k in (:atomic_number, :atomic_symbol, :atomic_mass, :position, :velocity))
-                atprops[k] = atom_data[k][i]
-            end
-        end
-        Atom(; atom, atprops...)
-    end
-    extra = Dict{Symbol,Any}()
-    for k in keys(system_data)
-        if !(k in (:box, :boundary_conditions))
-            extra[k] = system_data[k]
-        end
-    end
-    FlexibleSystem(atoms, bounding_box(extxyz_atoms),
-                   boundary_conditions(extxyz_atoms); extra...)
-end
-
-
-function convert_extxyz(system::AbstractSystem)
-    dict  = ExtXYZ.write_dict(make_atoms(system))
-
-    # Strip off units from the AtomsBase-supported extra properties
-    # (the ones where units are attached back on in parse_extxyz
-    for k in keys(dict["arrays"])
-        if k in ("charge", "vdw_radius", "covalent_radius")
-            dict["arrays"][k] = austrip.(dict["arrays"][k])
-        end
-    end
-    for k in keys(dict["info"])
-        if k in ("charge", )
-            dict["info"][k] = austrip.(dict["info"][k])
-        end
-    end
-
-    # Fix velocities (which write_dict currently gets wrong)
-    # TODO Fix upstream
-    velocities = map(system) do atom
-        vel = zeros(3)
-        if !ismissing(velocity(atom))
-            vel[1:n_dimensions(system)] = ustrip.(u"Å/s", velocity(atom))
-        end
-        vel
-    end
-    dict["arrays"]["velocities"] = hcat(velocities...)
-
-    dict
-end
-
-function make_atoms(system::AbstractSystem{D}) where {D}
-    D != 3 && @warn "1D and 2D systems not yet fully supported."
-
-    # Types supported losslessly in ExtXYZ
-    ExtxyzType = Union{Integer, AbstractFloat, AbstractString}
-
-    n_atoms = length(system)
-    atomic_symbols = [Symbol(PeriodicTable.elements[atomic_number(at)].symbol)
-                      for at in system]
-    if atomic_symbols != atomic_symbol(system)
-        @warn("Mismatch between atomic numbers and atomic symbols, which is not supported " *
-              "in ExtXYZ. Atomic numbers take preference.")
-    end
-    atom_data = Dict{Symbol,Any}(
-        :atomic_symbols => atomic_symbols,
-        :atomic_numbers => atomic_number(system),
-        :atomic_masses  => atomic_mass(system)
-    )
-    atom_data[:positions] = map(1:n_atoms) do at
-        pos = zeros(3)u"Å"
-        pos[1:D] = position(system, at)
-        pos
-    end
-    atom_data[:velocities] = map(1:n_atoms) do at
-        vel = zeros(3)u"Å/s"
-        if !ismissing(velocity(system)) && !ismissing(velocity(system, at))
-            vel[1:D] = velocity(system, at)
-        end
-        vel
-    end
-
-    # Extract extra atomic properties
-    if first(system) isa Atom
-        # TODO not a good idea to directly access the field
-        # TODO Implement and make use of a property interface on the atom level
-        for (k, v) in system[1].data
-            atoms_base_keys = (:charge, :covalent_radius, :vdw_radius,
-                               :magnetic_moment, :pseudopotential)
-            if k in atoms_base_keys || v isa ExtxyzType
-                # These are either Unitful quantities, which are uniformly supported
-                # across all of AtomsBase or the value has a type that Extxyz can write
-                # losslessly, so we can just write them no matter the value
-                atom_data[k] = [at.data[k] for at in system]
-            elseif v isa Quantity || (v isa AbstractArray && eltype(v) <: Quantity)
-                @warn "Unitful quantity $k is not yet supported in convert_extxyz."
-            else
-                @warn "Writing quantities of type $(typeof(v)) is not supported in convert_extxyz."
-            end
-        end
-    end
-
-    box = map(1:3) do i
-        v = zeros(3)u"Å"
-        i ≤ D && (v[1:D] = bounding_box(system)[i])
-        v
-    end
-    system_data = Dict{Symbol,Any}(
-        :box => box,
-        :boundary_conditions => boundary_conditions(system)
-    )
-
-    # Extract extra system properties
-    if system isa FlexibleSystem
-        # TODO not a good idea to directly access the field
-        # TODO Implement and make use of a property interface on the system level
-        for (k, v) in system.data
-            atoms_base_keys = (:charge, :multiplicity)
-            if k in atoms_base_keys || v isa ExtxyzType
-                # These are either Unitful quantities, which are uniformly supported
-                # across all of AtomsBase or the value has a type that Extxyz can write
-                # losslessly, so we can just write them no matter the value
-                system_data[k] = v
-            elseif v isa Quantity || (v isa AbstractArray && eltype(v) <: Quantity)
-                @warn "Unitful quantity $k is not yet supported in convert_extxyz."
-            else
-                @warn "Writing quantities of type $(typeof(v)) is not supported in convert_extxyz."
-            end
-        end
-    end
-
-    ExtXYZ.Atoms(NamedTuple(atom_data), NamedTuple(system_data))
+    ExtXYZ.write_frames(file, ExtXYZ.write_dict.(systems))
 end
-make_atoms(system::ExtXYZ.Atoms) = system