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.

Sign up for GitHub

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

Jump to bottom

Fix for issue #37 - missing cell in extxyz files #39

Closed
wants to merge 2 commits into from
Closed

Fix for issue #37 - missing cell in extxyz files #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d8b3842
Fix for issue #37
jameskermode Sep 20, 2023
00c4c68
Merge branch 'master' into issue37
jameskermode Sep 20, 2023
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
1 change: 1 addition & 0 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ version = "0.1.14-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"
Expand Down
19 changes: 19 additions & 0 deletions src/atoms.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
using LinearAlgebra
using AtomsBase
using Unitful
using UnitfulAtomic
Expand Down Expand Up @@ -166,6 +167,24 @@ function Atoms(dict::Dict{String, Any})
end
end

system_data = _dict_remap_fwd(dict["info"])

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
read_dict(dict::Dict{String,Any}) = Atoms(dict)
Expand Down
10 changes: 6 additions & 4 deletions src/fileio.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
211 changes: 209 additions & 2 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,213 @@
using Test

@testset "ExtXYZ.jl" begin
include("dict.jl")
include("atomsbase.jl")
test_xyz(frame) = """2
cutoff=5.50000000 pbc=[T, T, T] nneightol=1.20000000 config_type=isolated_atom matrix=[[$frame.00000000, 1.00000000, 2.00000000], [3.00000000, 4.00000000, 5.00000000], [6.00000000, 7.00000000, 8.00000000]] dft_energy=-158.$(frame)4496821 Lattice="1.00000000 2.00000000 3.00000000 4.00000000 5.00000000 6.00000000 7.00000000 $(frame).00000000 9.00000000" gap_energy=-157.72725320 Properties=species:S:1:pos:R:3:n_neighb:I:1:dft_force:R:3:gap_force:R:3:map_shift:I:3
Si 10.00000000 11.00000000 $frame.00000000 0 0.10000000 0.20000000 0.30000000 0.1$(frame)000000 0.22000000 0.330000000 -1 -2 0
Si 13.00000000 14.00000000 $(frame+1).00000000 0 0.10000000 0.20000000 0.30000000 0.2$(frame)000000 0.22000000 0.330000000 -1 -2 0
"""

infile = "test.xyz"
open(infile, "w") do io
for frame=1:10
print(io, test_xyz(frame))
end
end
outfile = "dump.xyz"

try
@testset "read" begin
seq1 = read_frames(infile)

@test all([s["cell"][3,2] for s in seq1] .== 1.0:10.0)
@test all([s["arrays"]["pos"][3,1] for s in seq1] .== 1.0:10.0)
@test all([s["arrays"]["pos"][3,2] for s in seq1] .== 2.0:11.0)
@test all([s["info"]["matrix"][1,1] for s in seq1] .== 1.0:10.0)
@test all([s["info"]["matrix"][2,1] for s in seq1] .== 3.0)

@test all([s["info"]["dft_energy"] for s in seq1] .== [parse(Float64, "-158.$(frame)4496821") for frame=1:10])
@test all([s["arrays"]["gap_force"][1,1] for s in seq1] .== [parse(Float64, "0.1$(frame)") for frame=1:10])
@test all([s["arrays"]["gap_force"][1,2] for s in seq1] .== [parse(Float64, "0.2$(frame)") for frame=1:10])

seq2 = read_frames(infile, 4:10)
@test all(seq1[4:10] .== seq2)

frame4 = read_frame(infile, 4)
@test frame4 == seq1[4]

frame1 = read_frame(infile)
@test frame1 == seq1[1]

f = open(infile, "r")
seq3 = read_frames(f)
close(f)

@test all(seq1 .== seq3)
end

@testset "convert" begin
nat, info, arrays = ExtXYZ.cfopen(infile, "r") do fp
ExtXYZ.read_frame_dicts(fp)
end
cinfo = convert(Ptr{ExtXYZ.DictEntry}, info)
carrays = convert(Ptr{ExtXYZ.DictEntry}, arrays)

ninfo = convert(Dict{String,Any}, cinfo)
narrays = convert(Dict{String,Any}, carrays)

@test ninfo == info
@test narrays == arrays
end

@testset "write" begin
seq1 = read_frames(infile)
write_frames(outfile, seq1)
seq2 = read_frames(outfile)
@test all(seq1 .≈ seq2) # use custom isapprox(), since expect some loss of precision on round-trip
end

@testset "iread" begin
seq1 = read_frames(infile)
seq2 = collect(iread_frames(infile))
@test all(seq1 .≈ seq2)
end

@testset "iwrite" begin
seq1 = read_frames(infile)
ch = Channel()
job = @async write_frames(outfile, ch)
for frame in seq1
put!(ch, frame)
end
close(ch)
wait(job)
seq2 = read_frames(outfile)
@test all(seq1 .≈ seq2)
end

@testset "missingfile" begin
try
read_frame("bla.extxyz")
catch e
buf = IOBuffer()
showerror(buf, e)
message = String(take!(buf))
@test message == "file bla.extxyz cannot be opened for reading"
end
end

@testset "AtomsBase" begin
seq1 = ExtXYZ.load(infile)
ExtXYZ.save(outfile, seq1)
seq2 = ExtXYZ.load(outfile)
@test all(seq1 .≈ seq2)
@test all(repr(seq1) == repr(seq2))

@testset "AtomsBase unit_tests" begin
at = ExtXYZ.load(infile, 1)
@test at ≈ seq1[1]
@test all(Test_Units.(vcat(seq1, seq2)))
atom_data = Dict(:positions => [[1,2,3],[2,3,4]].*10^(-12)*u"m", :velocities => [[3,4,5],[4,5,6]]*u"Å/(Å*sqrt(u/eV))", :atomic_numbers => [25, 25], :atomic_symbols => [:Mn, :Mn])
system_data = Dict(:box => [[3,3,3],[4,4,4],[5,5,5]]*u"Å", :boundary_conditions => [Periodic(), Periodic(), Periodic()],
:dft_energy => 10*u"eV", :Momentum => 11*u"u*Å/(Å*sqrt(u/eV))", :Force => 12*u"eV/Å", :Stress => [13,14]*u"eV/Å^3", :conv_vec => [10,11,12].*10^(-12)*u"m", :conv_val => 13*10^(-12)*u"m")
sys1 = Atoms(NamedTuple(atom_data), NamedTuple(system_data))
ExtXYZ.save(outfile, sys1)
sys2 = ExtXYZ.load(outfile)
@test(position(sys2) == [[0.010,0.020,0.030],[0.020,0.030,0.040]]*u"Å")
@test(velocity(sys2) == [[3,4,5],[4,5,6]]*u"Å/(Å*sqrt(u/eV))")
@test(sys2.system_data[:dft_energy] == 10)
@test(sys2.system_data[:Momentum] == 11)
@test(sys2.system_data[:Force] == 12)
@test(sys2.system_data[:Stress] == [13,14])
@test(sys2.system_data[:conv_vec] == [0.10, 0.11, 0.12])
@test(sys2.system_data[:conv_val] == 0.13)
#Test without velocities
atom_data = Dict(:positions => [[1,2,3],[2,3,4]].*10^(-12)*u"m", :atomic_numbers => [25, 25], :atomic_symbols => [:Mn, :Mn])
sys1 = Atoms(NamedTuple(atom_data), NamedTuple(system_data))
ExtXYZ.save(outfile, sys1)
sys2 = ExtXYZ.load(outfile)
@test(ismissing(velocity(sys2)))
@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
ase_io = pyimport("ase.io")

@testset "ASE" begin
seq = read_frames(infile)
ase_seq = ase_io.read(infile * "@:")

for (frame, ase_atoms) in zip(seq, ase_seq)
@test frame["N_atoms"] ≈ length(ase_atoms)
@test frame["arrays"]["pos"] ≈ ase_atoms.positions'
@test frame["arrays"]["gap_force"] ≈ ase_atoms.arrays["gap_force"]'
@test frame["arrays"]["n_neighb"] ≈ ase_atoms.arrays["n_neighb"]
@test frame["cell"] ≈ ase_atoms.cell.array
end
end
catch e
if e isa PyError
println("ASE not installed, skipping ASE comparison tests")
else
throw(e)
end
end

finally
rm(infile, force=true)
rm(outfile, force=true)
end
end
Loading