From d8b38421b71b133c3fff493e12dc8fe939d558d5 Mon Sep 17 00:00:00 2001 From: James Kermode Date: Wed, 20 Sep 2023 12:05:12 +0100 Subject: [PATCH] Fix for issue #37 --- Project.toml | 3 ++- src/atoms.jl | 20 ++++++++++++++++--- src/fileio.jl | 10 ++++++---- test/runtests.jl | 51 ++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 76 insertions(+), 8 deletions(-) diff --git a/Project.toml b/Project.toml index c77dfdb..af23fa5 100644 --- a/Project.toml +++ b/Project.toml @@ -5,17 +5,18 @@ version = "0.1.9-DEV" [deps] AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PeriodicTable = "7b2266bf-644c-5ea3-82d8-af4bbd25a884" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" extxyz_jll = "6ecdc6fc-93a8-5528-aee3-ac7ae1c60be7" [compat] -julia = "1" AtomsBase = "0.2" PeriodicTable = "1" StaticArrays = "1.5" Unitful = "1" +julia = "1" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/atoms.jl b/src/atoms.jl index d8c0481..5fba5db 100644 --- a/src/atoms.jl +++ b/src/atoms.jl @@ -1,3 +1,4 @@ +using LinearAlgebra using AtomsBase using PeriodicTable using StaticArrays @@ -99,12 +100,25 @@ function Atoms(dict::Dict{String}{Any}) else atom_data[:atomic_masses] .* u"u" end - + system_data = _dict_remap_fwd(dict["info"]) - system_data[:box] = [dict["cell"][i, :] for i in 1:D ].*u"Å" # lattice vectors are rows from cell matrix - pbc = get(dict, "pbc", [true, true, true]) # default to periodic in all directions + + pbc = dict["pbc"] system_data[:boundary_conditions] =[p ? Periodic() : DirichletZero() for p in pbc] + if haskey(dict, "cell") + box = [dict["cell"][i, :] for i in 1:D ].*u"Å" # lattice vectors are rows from cell matrix + else + if (pbc == [false, false, false]) + extents = [extrema(getindex.(atom_data[:positions], dim)) for dim=1:3] + sizes = [extents[dim][2] - extents[dim][1] for dim=1:3] + box = diagm(sizes) + else + error("No cell specified but pbc = $pbc. This is inconsistent.") + end + end + system_data[:box] = box + Atoms(NamedTuple(atom_data), NamedTuple(system_data)) end diff --git a/src/fileio.jl b/src/fileio.jl index 67c4879..8318f68 100644 --- a/src/fileio.jl +++ b/src/fileio.jl @@ -269,15 +269,17 @@ function read_frame(fp::Ptr{Cvoid}; verbose=false) dict = Dict{String, Any}() dict["N_atoms"] = nat # number of atoms + # cell is transpose of the stored lattice + lattice = extract_lattice!(info) + if (!isnothing(lattice)) dict["cell"] = permutedims(lattice, (2, 1)) end + # periodic boundary conditions + # default is "FFF" if no cell, "TTT" if cell present + dict["pbc"] = isnothing(lattice) ? [false, false, false] : [true, true, true] if "pbc" in keys(info) dict["pbc"] = pop!(info, "pbc") end - # cell is transpose of the stored lattice - lattice = extract_lattice!(info) - if (!isnothing(lattice)) dict["cell"] = permutedims(lattice, (2, 1)) end - delete!(info, "Properties") # everything else stays in info and arrays diff --git a/test/runtests.jl b/test/runtests.jl index 973b8a1..4cf83c7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -161,6 +161,57 @@ Si 13.00000000 14.00000000 $(frame+1).00000000 0 @test(ismissing(velocity(sys2[1]))) @test(ismissing(velocity(sys2, 1))) end + + @testset "AtomsBase missing cell" begin + fname = tempname() + try + # We do this twice, first with pbc = "F F F" and no cell - should work and give a bounding box + # then again with a pbc = "T T T" but with cell missing - should raise an error + for pbc in ["F F F", "T T T"] + + text = """13 + Properties=species:S:1:pos:R:3:forces:R:3 energy=-66.79083251953125 pbc="$pbc" + C -5.13553286 0.00000000 0.00000000 0.00006186 -0.00000000 0.00000000 + H -5.72485781 -0.91726500 0.00000000 -0.00001280 0.00000756 0.00000000 + H -5.72485781 0.91726500 0.00000000 -0.00001280 -0.00000756 -0.00000000 + C -3.82464290 0.00000000 0.00000000 -0.00004186 0.00000000 0.00000000 + C -2.55226111 0.00000000 0.00000000 0.00003090 -0.00000000 -0.00000000 + C -1.27494299 0.00000000 0.00000000 -0.00009426 0.00000000 -0.00000000 + C 0.00000000 0.00000000 0.00000000 0.00000000 -0.00000000 -0.00000000 + C 1.27494299 0.00000000 0.00000000 0.00009426 0.00000000 0.00000000 + C 2.55226111 0.00000000 0.00000000 -0.00003090 -0.00000000 0.00000000 + C 3.82464290 0.00000000 0.00000000 0.00004186 0.00000000 -0.00000000 + H 5.72485781 0.00000000 0.91726500 0.00001280 -0.00000000 -0.00000756 + H 5.72485781 0.00000000 -0.91726500 0.00001280 0.00000000 0.00000756 + C 5.13553286 0.00000000 0.00000000 -0.00006186 -0.00000000 -0.00000000 + """ + open(fname, "w") do io + print(io, text) + end + + if pbc == "F F F" + atoms = ExtXYZ.load(fname) + pos = position(atoms) + box = bounding_box(atoms) + frac_pos = inv(box) * hcat(pos...) + @test maximum(frac_pos) <= 0.5 + @test minimum(frac_pos) >= -0.5 + else + try + ExtXYZ.load(fname) + catch e + buf = IOBuffer() + showerror(buf, e) + message = String(take!(buf)) + @test message == "No cell specified but pbc = Bool[1, 1, 1]. This is inconsistent." + end + end + end + + finally + isfile(fname) && rm(fname) + end + end end try