Merge pull request #9 from JuliaMolSim/co/ab04

Update to [email protected]

Project.toml

@@ -1,6 +1,6 @@
 name = "AtomsBuilder"
 uuid = "f5cc8831-eeb7-4288-8d9f-d6c1ddb77004"
-version = "0.1.0"
+version = "0.2.0-dev"
 
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
@@ -11,7 +11,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
-AtomsBase = "0.3.5"
+AtomsBase = "0.4"
 JSON = "0.21.4"
 LinearAlgebra = "1.9.0, 1.10.0"
 Random = "1.9.0, 1.10.0"
@@ -20,8 +20,7 @@ Unitful = "1.19"
 julia = "1.9.0, 1.10.0"
 
 [extras]
-JuLIP = "945c410c-986d-556a-acb1-167a618e0462"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "JuLIP"]
+test = ["Test", ]

data/generate_testdata.jl

@@ -0,0 +1,35 @@
+
+# this requires to be run in an environment that has JuLIP and JSON
+# installed. 
+using JuLIP, JSON 
+
+make_dict(sym, cubic, pbc, nn, at) = 
+   Dict("sym" => sym, "cubic" => cubic, "pbc" => pbc, "nn" => nn, 
+        "sys" => write_dict(at))
+
+test_systems = []    
+
+for sym in [:Si, :Ge, :W, :Ti]
+   at = bulk(sym)
+   push!(test_systems, make_dict(sym, false, true, 1, at))
+
+   at = bulk(sym, cubic=true) 
+   push!(test_systems, make_dict(sym, true, true, 1, at))
+end
+
+for ntest = 1:30
+   sym = rand([:Si, :Ge, :W, :Ti]) 
+   cubic = rand(Bool)
+   pbc = tuple(rand(Bool, 3)...)
+   nn = (rand(1:3), rand(2:4), 2)
+
+   at = bulk(sym, cubic=cubic, pbc = pbc) * nn 
+   push!(test_systems, make_dict(sym, cubic, pbc, nn, at))
+end
+
+open(@__DIR__() * "/test_systems.json", "w") do f
+   JSON.print(f, test_systems)
+end
+
+# read the file back 
+# test_systems2 = JSON.parsefile("test_systems.json")

src/bulk.jl

@@ -46,6 +46,7 @@ end
 
 _convert_pbc(pbc::NTuple{3, Bool}) = pbc
 _convert_pbc(pbc::Bool) = (pbc, pbc, pbc)
+_convert_pbc(pbc::AbstractVector) = tuple(pbc...)
 
 """
 `bulk(sym)` : generates a `FlexibleSystem` unit cell for a bulk 

src/utils.jl

@@ -1,13 +1,13 @@
 
 using AtomsBase 
-using AtomsBase: Atom, FlexibleSystem, Periodic
+using AtomsBase: Atom, FlexibleSystem
 using Unitful: unit, ustrip, Quantity
 using LinearAlgebra: norm 
 
 export rattle!, 
        set_positions, 
        set_elements, 
-         randz!
+       randz!
 
 """
 Helper function to convert construct a FlexibleSystem from 
@@ -17,29 +17,29 @@ function _flexible_system(positions, elements, cell, pbc)
    Nat = length(positions)
    syms = Chemistry.chemical_symbol.(elements)
    atoms = [ Atom(syms[i], positions[i]) for i in 1:Nat ]
-   bc =  [ (pbc[i] ? Periodic() : nothing) for i = 1:3 ]
-   bb = [cell[i, :] for i = 1:3]   # IS THIS A BUG? SHOULD IT BE cell[:, i] ??? 
+   bc = pbc isa Bool ? (pbc, pbc, pbc) : tuple(pbc...)
+   bb = tuple([cell[i, :] for i = 1:3]...)
    return FlexibleSystem(atoms; 
                          bounding_box = bb, 
-                         boundary_conditions = bc)
+                         periodicity = bc)
 end
 
 _set_position(x::Atom, 𝐫) = Atom(atomic_number(x), 𝐫; 
-                                   velocity = x.velocity,
-                                   atomic_mass = x.atomic_mass)
+                                   velocity = velocity(x),
+                                   mass = mass(x))
 
 _set_element(x::Atom, Z) = Atom(Z, position(x); 
-                                velocity = x.velocity,
-                                atomic_mass = x.atomic_mass)
+                                velocity = velocity(x),
+                                mass = mass(x))
 
 function set_positions(at::FlexibleSystem, 
                        X::AbstractVector{SVector{3, T}}) where {T}
    @assert length(X) == length(at)                       
    particles = [ _set_position(at.particles[i], X[i])
                  for i in 1:length(at) ] 
    return FlexibleSystem(particles, 
-                         bounding_box(at), 
-                         boundary_conditions(at))
+                         bounding_box = bounding_box(at), 
+                         periodicity = periodicity(at))
 end
 
 
@@ -48,8 +48,8 @@ function set_elements(at::FlexibleSystem, Z::AbstractVector)
    particles = [ Atom(Z[i], position(x), velocity(x))
                  for (i, x) in enumerate(at.particles) ]
    return FlexibleSystem(particles, 
-                         bounding_box(at), 
-                         boundary_conditions(at))
+                         bounding_box = bounding_box(at), 
+                         periodicity = periodicity(at))
 end
 
 
@@ -84,8 +84,9 @@ function Base.repeat(at::FlexibleSystem, n::NTuple{3})
       end
    end
 
-   bb = [c1 * n[1], c2 * n[2], c3 * n[3]]
-   return FlexibleSystem(particles, bb, boundary_conditions(at))
+   bb = (c1 * n[1], c2 * n[2], c3 * n[3])
+   return FlexibleSystem(particles; bounding_box = bb, 
+                           periodicity = periodicity(at))
 end
 
 Base.repeat(at::FlexibleSystem, n::Integer) = repeat(at, (n,n,n))
@@ -120,7 +121,7 @@ function rattle!(at::FlexibleSystem, r::Quantity)
 end
 
 rattle!(sys::FlexibleSystem, r::AbstractFloat) = 
-      rattle!(sys, r * unit(position(sys)[1][1]))
+      rattle!(sys, r * unit(position(sys, 1)[1]))
 
 
 """
@@ -168,11 +169,11 @@ union(sys1::FlexibleSystem, sys2::FlexibleSystem)
 takes the union of two particle systems provided their cells are identical. 
 """
 function union(sys1::FlexibleSystem, sys2::FlexibleSystem) 
-   @assert boundary_conditions(sys1) == boundary_conditions(sys2)
+   @assert periodicity(sys1) == periodicity(sys2)
    @assert bounding_box(sys1) == bounding_box(sys2)
    return FlexibleSystem(union(sys1.particles, sys2.particles),  
-                        bounding_box(at), 
-                        boundary_conditions(at) )
+                        bounding_box = bounding_box(at), 
+                        periodicit = periodicity(at) )
 end
 
 """

test/test_bulk.jl

@@ -1,40 +1,33 @@
 
-using AtomsBuilder, Test, AtomsBase, Unitful, Random
-import JuLIP
+using AtomsBuilder
+using Test, AtomsBase, Unitful, Random, JSON
 
 ##
 
-@info("Testing `bulk` and `repeat` against JuLIP reference implementation")
 
-function _compare_particle(x1, x2)
-   return (x1.position ≈ x2.position) && (x1.atomic_symbol == x2.atomic_symbol)
-end
-
-function _compare_at(at, at_f)
-   return ( at.boundary_conditions == at_f.boundary_conditions &&
-            at.bounding_box == at_f.bounding_box &&
-            all(_compare_particle.(at.particles, at_f.particles)) )
-end
+@info("Testing `bulk` and `repeat` against JuLIP reference data")
 
-for sym in [:Si, :Ge, :W, :Ti]
-   at = bulk(sym)
-   at_f = FlexibleSystem(JuLIP.bulk(sym))
-   @test _compare_at(at, at_f)
+test_systems = JSON.parsefile(joinpath(@__DIR__(), "..", "data", "test_systems.json"))
 
-   at = bulk(sym, cubic=true) 
-   at_f = FlexibleSystem(JuLIP.bulk(sym, cubic=true))
-   @test _compare_at(at, at_f)
+_ustripvecvec(X) = [ ustrip.(x) for x in X ]
 
-   for nn in (rand(1:3), rand(2:4), 2)
-      at = bulk(sym) * nn 
-      at_f = FlexibleSystem(JuLIP.bulk(sym) * nn)
-      @test _compare_at(at, at_f)
+compare_system(sys_f, sys_j) = (
+      all( ustrip.( hcat(bounding_box(sys_f)...) ) .≈ hcat(sys_j["cell"]...)' )  && 
+      all( AtomsBuilder._convert_pbc(sys_j["pbc"]) .== periodicity(sys_f) ) && 
+      all( atomic_number(sys_f, :) .== sys_j["Z"] ) && 
+      all( _ustripvecvec(position(sys_f, :)) .≈ sys_j["X"] )    )
 
-      at = bulk(sym, cubic=true) * nn
-      at_f = FlexibleSystem(JuLIP.bulk(sym, cubic=true) * nn)
-      @test _compare_at(at, at_f)
-   end
+for D in test_systems 
+   nn = D["nn"] isa Integer ? D["nn"] : tuple(D["nn"]...)
+   sys_f = bulk( Symbol(D["sym"]); 
+                 cubic = D["cubic"], 
+                 pbc = D["pbc"] ) 
+   if nn != 1 
+      sys_f = sys_f * nn 
+   end 
+   @test compare_system(sys_f, D["sys"])
 end
+
 
 ##
 
@@ -43,23 +36,24 @@ sys0 = rattle!(bulk(:C, cubic=true) * (2,3,4), 0.1u"Å")
 sys1 = rattle!(bulk(:C, cubic=true) * (2,3,4), 0.1)
 sys2 = rattle!(bulk(:C) * (2,3,4), 0.1)
 
-X = position(sys1)
+X = position(sys1, :)
 Xnew = [ x .+ 0.01u"Å" for x in X ]
 sys3 = set_positions(sys1, Xnew)
-@test position(sys3) == Xnew
+@test position(sys3, :) == Xnew
 
-Z = atomic_number(sys1)
+Z = atomic_number(sys1, :)
 @test all(Z .== AtomsBuilder.Chemistry.atomic_number(:C))
 zO = AtomsBuilder.Chemistry.atomic_number(:O)
 Znew = copy(Z); Znew[3:5:end] .= zO
 sys4 = set_elements(sys3, Znew)
-@test all(atomic_number(sys4) .== Znew)
+@test all(atomic_number(sys4, :) .== Znew)
 
 pold = deepcopy(sys4.particles) 
 deleteat!(sys4, 1:5)
-@test position.(pold[6:end]) == position(sys4)
-@test atomic_mass.(pold[6:end]) == atomic_mass(sys4)
-@test atomic_number.(pold[6:end]) == atomic_number(sys4)
+@test position.(pold[6:end]) == position(sys4, :)
+@test mass.(pold[6:end]) == mass(sys4, :)
+@test atomic_number.(pold[6:end]) == atomic_number(sys4, :)
+
 
 
 

test/test_utils.jl

@@ -6,7 +6,7 @@ using AtomsBuilder, Test, AtomsBase, Unitful, Random
 Random.seed!(1234)
 sys = bulk(:Ti, cubic=true) * 5
 sys = randz!(sys, [ :Ti => 0.5, :O => 0.5 ])
-Z = atomic_number(sys)
-@test count( Z .== 8 ) / length(Z) > 0.45
-@test count( Z .== 22 ) / length(Z) > 0.45
+Z = atomic_number(sys, :)
+@test count( Z .== 8 ) / length(Z) > 0.4
+@test count( Z .== 22 ) / length(Z) > 0.4