Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update to AtomsBase 0.5 #77

Merged
merged 7 commits into from
Nov 27, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 2 additions & 4 deletions Project.toml
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
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
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