Skip to content

Commit

Permalink
Update to AtomsBase 0.5 (#77)
Browse files Browse the repository at this point in the history
  • Loading branch information
@mfherbst
mfherbst authored Nov 27, 2024
1 parent 45574b9 commit 764a9d3
Show file tree
Hide file tree
Showing 3 changed files with 59 additions and 38 deletions.
6 changes: 2 additions & 4 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
57 changes: 33 additions & 24 deletions src/atomsbase.jl
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
import AtomsBase
using AtomsBase: AbstractSystem, FlexibleSystem
using AtomsBase: AbstractSystem, FlexibleSystem, ChemicalSpecies, PeriodicCell, IsolatedCell
using Unitful
using UnitfulAtomic
import PeriodicTable


"""
Expand All @@ -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"Å",
Expand All @@ -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
Expand All @@ -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

"""
Expand All @@ -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)
Expand All @@ -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))
Expand Down Expand Up @@ -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)))
Expand Down
34 changes: 24 additions & 10 deletions test/atomsbase.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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...)
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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"),
Expand All @@ -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)
Expand All @@ -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

0 comments on commit 764a9d3

Please sign in to comment.