diff --git a/Project.toml b/Project.toml index eb0f214..71d61f5 100644 --- a/Project.toml +++ b/Project.toml @@ -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" diff --git a/src/atoms.jl b/src/atoms.jl index 4bc933c..7a84e8c 100644 --- a/src/atoms.jl +++ b/src/atoms.jl @@ -1,3 +1,4 @@ +using LinearAlgebra using AtomsBase using Unitful using UnitfulAtomic @@ -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) diff --git a/src/fileio.jl b/src/fileio.jl index af8bd5a..cd98a0e 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 c67e61b..0cddcd5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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