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 1 commit
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.4"
AtomsBaseTesting = "0.2"
Chemfiles_jll = "= 0.10.4"
DocStringExtensions = "0.8.3, 0.9"
PeriodicTable = "1"
Unitful = "1"
UnitfulAtomic = "1"
julia = "1.6"
Expand Down
45 changes: 24 additions & 21 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, ChemialSpecies, PeriodicCell, IsolatedCell
using Unitful
using UnitfulAtomic
import PeriodicTable


"""
Expand All @@ -11,17 +10,20 @@ 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 = AtomsBase.ChemicalSpecies(Chemfiles.atomic_number(atom))
mfherbst marked this conversation as resolved.
Show resolved Hide resolved
if Symbol(species) != Symbol(Chemfiles.name(atom))
@warn "Ignoring non-standard atom name $(Chemfiles.name(atom)) for atom $i."
mfherbst marked this conversation as resolved.
Show resolved Hide resolved
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 => AtomsBase.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 +35,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 +59,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,12 +90,12 @@ 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
Expand All @@ -109,16 +113,15 @@ function Base.convert(::Type{Frame}, system::AbstractSystem{D}) where {D}
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)
# usually has things like annotations for the isotopes etc., which in
# Chemfiles is mapped to the "name" of an atom
cf_atom = Chemfiles.Atom(AtomsBase.element(AtomsBase.species(atom)).symbol)
Chemfiles.set_name!(cf_atom, string(AtomsBase.atomic_symbol(atom)))
Chemfiles.set_mass!(cf_atom, ustrip(u"u", AtomsBase.atomic_mass(atom)))
Chemfiles.set_mass!(cf_atom, ustrip(u"u", AtomsBase.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)
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 +150,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 (:bounding_box, :periodicity)
continue # Already dealt with
elseif k in (:charge, )
Chemfiles.set_property!(frame, string(k), Float64(ustrip(u"e_au", v)))
Expand Down
16 changes: 10 additions & 6 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, atomic_species
using AtomsBaseTesting

function make_chemfiles_system(D=3; drop_atprop=Symbol[], infinite=false, kwargs...)
Expand All @@ -13,10 +15,11 @@ 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(; cell_vectors=data.box, periodicity=(true, true, true))
end
system = AtomsBase.FlexibleSystem(data.atoms, cell; data.sysprop...)
merge(data, (; system))
end

Expand All @@ -37,14 +40,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 Down
Loading