diff --git a/Project.toml b/Project.toml index ffa22fe..6c1f8de 100644 --- a/Project.toml +++ b/Project.toml @@ -7,16 +7,14 @@ version = "0.10.41" 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" +AtomsBase = "0.5" +AtomsBaseTesting = "0.4" Chemfiles_jll = "= 0.10.4" DocStringExtensions = "0.8.3, 0.9" -PeriodicTable = "1" Unitful = "1" UnitfulAtomic = "1" julia = "1.6" diff --git a/src/atomsbase.jl b/src/atomsbase.jl index ec4b0ff..aef34d7 100644 --- a/src/atomsbase.jl +++ b/src/atomsbase.jl @@ -1,8 +1,7 @@ import AtomsBase -using AtomsBase: AbstractSystem, FlexibleSystem +using AtomsBase: AbstractSystem, FlexibleSystem, ChemicalSpecies, PeriodicCell, IsolatedCell using Unitful using UnitfulAtomic -import PeriodicTable """ @@ -11,17 +10,22 @@ 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"Å" + velocity_arg = () if Chemfiles.has_velocities(frame) - velocity = Chemfiles.velocities(frame)[:, i]u"Å/ps" - else - velocity = zeros(3)u"Å/ps" + velocity_arg = (Chemfiles.velocities(frame)[:, i]u"Å/ps", ) + end + + species = ChemicalSpecies(Chemfiles.atomic_number(atom); + atom_name=Symbol(Chemfiles.name(atom))) + if Symbol(species) != Symbol(Chemfiles.type(atom)) + @warn("Non-standard atom type $(Chemfiles.type(atom)) " * + "for atom $i taken as $(Symbol(species)) in AtomsBase.") 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), + :mass => Chemfiles.mass(atom)u"u", + :species => ChemicalSpecies(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"Å", @@ -33,7 +37,7 @@ function AtomsBase.FlexibleSystem(frame::Chemfiles.Frame) end end - AtomsBase.Atom(Chemfiles.atomic_number(atom), pos, velocity; atprops...) + AtomsBase.Atom(Chemfiles.atomic_number(atom), pos, velocity_arg...; atprops...) end # Collect system properties @@ -57,12 +61,14 @@ function AtomsBase.FlexibleSystem(frame::Chemfiles.Frame) # Construct flexible system cell_shape = Chemfiles.shape(Chemfiles.UnitCell(frame)) if cell_shape == Chemfiles.Infinite - AtomsBase.isolated_system(atoms; sysprops...) + cell = IsolatedCell(3) else @assert cell_shape in (Chemfiles.Triclinic, Chemfiles.Orthorhombic) - box = collect(eachrow(Chemfiles.matrix(Chemfiles.UnitCell(frame))))u"Å" - AtomsBase.periodic_system(atoms, box; sysprops...) + cell_vectors = collect(eachrow(Chemfiles.matrix(Chemfiles.UnitCell(frame))))u"Å" + cell = PeriodicCell(; cell_vectors, periodicity=(true, true, true)) end + + AtomsBase.FlexibleSystem(atoms, cell; sysprops...) end """ @@ -86,19 +92,19 @@ function Base.convert(::Type{Frame}, system::AbstractSystem{D}) where {D} frame = Chemfiles.Frame() # Cell and boundary conditions - if AtomsBase.bounding_box(system) == AtomsBase.infinite_box(D) # System is infinite + if AtomsBase.cell(system) isa IsolatedCell 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)) + if !all(AtomsBase.periodicity(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]) + box[i, 1:D] = ustrip.(u"Å", AtomsBase.cell_vectors(system)[i]) end cell = Chemfiles.UnitCell(box) Chemfiles.set_cell!(frame, cell) @@ -108,17 +114,20 @@ function Base.convert(::Type{Frame}, system::AbstractSystem{D}) where {D} 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))) + # We are using the symbol of the symbol from the element here + # instead of the AtomsBase.atomic_symbol because the latter may have + # isotope information attached (e.g. C13), which Chemfiles cannot parse. + cf_atom = Chemfiles.Atom(string(AtomsBase.element(atom).symbol)) + Chemfiles.set_name!(cf_atom, string(AtomsBase.atom_name(atom))) + Chemfiles.set_mass!(cf_atom, ustrip(u"u", AtomsBase.mass(atom))) + + if string(AtomsBase.atomic_symbol(atom)) != string(AtomsBase.element(atom).symbol) + @warn "Custom neutron count or custom atomic symbol not supported." + end @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) + if k in (:species, :mass, :velocity, :position) continue # Dealt with separately elseif k == :charge Chemfiles.set_charge!(cf_atom, ustrip(u"e_au", v)) @@ -147,7 +156,7 @@ function Base.convert(::Type{Frame}, system::AbstractSystem{D}) where {D} end for (k, v) in pairs(system) - if k in (:bounding_box, :boundary_conditions) + if k in (:cell_vectors, :periodicity) continue # Already dealt with elseif k in (:charge, ) Chemfiles.set_property!(frame, string(k), Float64(ustrip(u"e_au", v))) diff --git a/test/atomsbase.jl b/test/atomsbase.jl index b049013..ea93e64 100644 --- a/test/atomsbase.jl +++ b/test/atomsbase.jl @@ -4,6 +4,8 @@ using UnitfulAtomic @testset "AtomsBase support" begin import AtomsBase using AtomsBase: AbstractSystem, FlexibleSystem + using AtomsBase: PeriodicCell, IsolatedCell + using AtomsBase: atomic_symbol using AtomsBaseTesting function make_chemfiles_system(D=3; drop_atprop=Symbol[], infinite=false, kwargs...) @@ -13,21 +15,23 @@ using UnitfulAtomic cellmatrix=:upper_triangular, kwargs...) if infinite - system = AtomsBase.isolated_system(data.atoms; data.sysprop...) + cell = IsolatedCell(3) else - system = AtomsBase.periodic_system(data.atoms, data.box; data.sysprop...) + cell = PeriodicCell(; data.cell_vectors, periodicity=(true, true, true)) end + system = AtomsBase.FlexibleSystem(data.atoms, cell; data.sysprop...) merge(data, (; system)) end @testset "Conversion AtomsBase -> Chemfiles (periodic, velocity)" begin - system, atoms, atprop, sysprop, box, bcs = make_chemfiles_system(; infinite=false) + data = make_chemfiles_system(; infinite=false) + system, atoms, cell, cell_vectors, periodicity, atprop, sysprop = data 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 + @test cell[i, :] ≈ ustrip.(u"Å", cell_vectors[i]) atol=1e-12 end @test Chemfiles.shape(Chemfiles.UnitCell(frame)) in (Chemfiles.Triclinic, Chemfiles.Orthorhombic) @@ -37,14 +41,15 @@ using UnitfulAtomic @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]) + species_i = atprop.species[i] + @test Chemfiles.name(atom) == string(atomic_symbol(species_i)) + @test Chemfiles.atomic_number(atom) == atomic_number(species_i) + @test Chemfiles.mass(atom) == ustrip(u"u", atprop.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 + if atomic_number(atprop.species[i]) == 1 @test Chemfiles.vdw_radius(atom) == 1.2 @test Chemfiles.covalent_radius(atom) == 0.37 end @@ -64,7 +69,8 @@ using UnitfulAtomic end @testset "Warning about setting invalid data" begin - system, atoms, atprop, sysprop, box, bcs = make_test_system() + data = make_test_system() + system, atoms, cell, cell_vectors, periodicity, atprop, sysprop = data 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"), @@ -74,7 +80,7 @@ using UnitfulAtomic D = 3 cell = Chemfiles.matrix(Chemfiles.UnitCell(frame)) for i in 1:D - @test cell[i, :] ≈ ustrip.(u"Å", box[i]) atol=1e-12 + @test cell[i, :] ≈ ustrip.(u"Å", cell_vectors[i]) atol=1e-12 end @test Chemfiles.shape(Chemfiles.UnitCell(frame)) in (Chemfiles.Triclinic, Chemfiles.Orthorhombic) @@ -89,4 +95,12 @@ using UnitfulAtomic test_approx_eq(system, newsystem; rtol=1e-12, ignore_atprop=[:covalent_radius, :vdw_radius]) end + + @testset "Make sure the water example can be parsed without error" begin + import AtomsBase + traj = Chemfiles.Trajectory(joinpath(@__DIR__, "data", "water.xyz")) + frame = Chemfiles.read_step(traj, 0) + sys = convert(AtomsBase.FlexibleSystem, frame) + @test length(sys) == length(frame) + end end