trying to move to strict use of interface

Project.toml

@@ -12,10 +12,10 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 AtomsBase = "0.3.5"
 JSON = "0.21.4"
-LinearAlgebra = "1.9"
+LinearAlgebra = "1.9.0, 1.10.0"
 StaticArrays = "1.9"
 Unitful = "1.19"
-julia = "1.9"
+julia = "1.9.0, 1.10.0"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/bulk.jl

@@ -1,11 +1,8 @@
-export bulk 
+export bulk
 
 using Unitful: unit 
 using LinearAlgebra: I
 
-using AtomsBuilder.Chemistry: symmetry, lattice_constant, atomic_mass, 
-                 atomic_number, chemical_symbol, refstate
-
 const _unit_cells = Dict(     # (positions, cell matrix, factor of a)
    :fcc => ( [ [0 0 0], ],
              [0 1 1; 1 0 1; 1 1 0],  0.5),
@@ -26,12 +23,12 @@ const _simple_structures = [:fcc, :bcc, :diamond]
 
 function _simple_bulk(sym::Symbol, cubic::Bool; a=nothing, T = Float64)
    if cubic
-      X, scal = _cubic_cells[symmetry(sym)]
+      X, scal = _cubic_cells[Chemistry.symmetry(sym)]
       C = Matrix(1.0I, 3, 3) / scal
    else
-      X, C, scal = _unit_cells[symmetry(sym)]
+      X, C, scal = _unit_cells[Chemistry.symmetry(sym)]
    end
-   a === nothing && (a = lattice_constant(sym))
+   a === nothing && (a = Chemistry.lattice_constant(sym))
    TU = typeof( one(T) * unit(a) )
    return [ Vec3{TU}(x * a * scal)  for x in X ], Mat3{TU}(C * a * scal)
 end
@@ -59,14 +56,14 @@ specify the kwargs `a` or `c` to determine the lattice constants.
 """
 function bulk(sym::Symbol; T=Float64, cubic = false, pbc = (true,true,true), 
                            a=nothing, c=nothing) # , x=nothing, y=nothing, z=nothing)
-   symm = symmetry(sym)
+   symm = Chemistry.symmetry(sym)
    if symm in _simple_structures
       X, C = _simple_bulk(sym, cubic; a=a)
    elseif symm == :hcp
       X, C = _bulk_hcp(sym; a=a, c=c)  # cubic parameter is irrelevant for hcp
    end
-   m = atomic_mass(sym)
-   Z = atomic_number(sym)
+   m = Chemistry.atomic_mass(sym)
+   Z = Chemistry.atomic_number(sym)
    nat = length(X)
    at = _flexible_system( X, fill(Z, nat), C, _convert_pbc(pbc))
    # I don't remember what this does so I'm commenting it out until I understand 

src/chemistry.jl

@@ -1,7 +1,6 @@
 module Chemistry
 
-using JSON 
-using Unitful 
+using JSON, Unitful 
 
 # --------------------------------------------------------------------
 # Load data and prepare it a bit ... 
@@ -32,7 +31,8 @@ end
 
 # --------------------------------------------------------------------
 #  some convenient accessor and conversion functions 
-
+#  these will not be exported and only accessed via Chemistry.atomic_number etc 
+#  in order to avoid confusion with the AtomsBase interface 
 
 atomic_number(sym::Symbol) = _numbers[sym]
 
@@ -42,9 +42,6 @@ chemical_symbol(sym::Symbol) = sym
 atomic_mass(z::Integer) = _masses[z+1]
 atomic_mass(sym::Symbol) = atomic_mass(atomic_number(sym))
 
-rnn(z::Integer) = _rnn[z+1]
-rnn(sym::Symbol) = rnn(atomic_number(sym))
-
 symmetry(z::Integer) = Symbol(_refstates[z+1]["symmetry"])
 symmetry(sym::Symbol) = symmetry(atomic_number(sym))
 
@@ -54,4 +51,11 @@ lattice_constant(sym::Symbol) = lattice_constant(atomic_number(sym))
 refstate(z::Integer) = _refstates[z+1]
 refstate(sym::Symbol) = refstate(atomic_number(sym))
 
+# The `rnn` could be exported though since it is quite useful for a lot of 
+# modelling scripts. I'll keep it here for now, unexported. Something to 
+# think about. 
+
+rnn(z::Integer) = _rnn[z+1]
+rnn(sym::Symbol) = rnn(atomic_number(sym))
+
 end

src/utils.jl

@@ -4,7 +4,9 @@ using AtomsBase: Atom, FlexibleSystem, Periodic
 using Unitful: unit, ustrip 
 using LinearAlgebra: norm 
 
-export rattle!
+export rattle!, 
+       set_positions, 
+       set_elements
 
 """
 Helper function to convert construct a FlexibleSystem from 
@@ -13,40 +15,37 @@ Helper function to convert construct a FlexibleSystem from
 function _flexible_system(positions, elements, cell, pbc)
    Nat = length(positions)
    syms = Chemistry.chemical_symbol.(elements)
-   atoms = [ Atom(; position = positions[i],
-                    atomic_symbol = syms[i]) for i in 1:Nat ]
+   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]
+   bb = [cell[i, :] for i = 1:3]   # IS THIS A BUG? SHOULD IT BE cell[:, i] ??? 
    return FlexibleSystem(atoms; 
                          bounding_box = bb, 
                          boundary_conditions = bc)
 end
 
-_set_position(x::Atom, 𝐫) = Atom(; position = 𝐫, 
+_set_position(x::Atom, 𝐫) = Atom(atomic_number(x), 𝐫; 
                                    velocity = x.velocity,
-                                   atomic_mass = x.atomic_mass,
-                                   atomic_number = x.atomic_number, 
-                                   atomic_symbol = x.atomic_symbol)
+                                   atomic_mass = x.atomic_mass)
 
-function _get_positions(at::FlexibleSystem)
-   [ x.position for x in at.particles ]
-end
 
-function _set_positions(at::FlexibleSystem, 
-                        X::AbstractVector{SVector{3, T}}) where {T}
-   particles = [ _set_position(at.particles[i], X[i]) 
+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, at.bounding_box, at.boundary_conditions)
+   return FlexibleSystem(particles, 
+                         bounding_box(at), 
+                         boundary_conditions(at))
 end
 
-function _get_atomic_numbers(at::FlexibleSystem)
-   [ x.atomic_number for x in at.particles ]
-end
 
-function _set_atomic_numbers(at::FlexibleSystem, Z::AbstractVector{<: Integer})
-   particles = [ Atom(Z[i], x.position, x.velocity)
+function set_elements(at::FlexibleSystem, Z::AbstractVector)
+   @assert length(Z) == length(at)                       
+   particles = [ Atom(Z[i], position(x), velocity(x))
                  for (i, x) in enumerate(at.particles) ]
-   return FlexibleSystem(particles, at.bounding_box, at.boundary_conditions)
+   return FlexibleSystem(particles, 
+                         bounding_box(at), 
+                         boundary_conditions(at))
 end
 
 
@@ -69,7 +68,7 @@ at = bulk(:) * (3, 2, 4)
 ```
 """
 function Base.repeat(at::FlexibleSystem, n::NTuple{3})
-   c1, c2, c3 = at.bounding_box
+   c1, c2, c3 = bounding_box(at)
 
    particles = eltype(at.particles)[] 
    for a in CartesianIndices( (1:n[1], 1:n[2], 1:n[3]) )
@@ -82,7 +81,7 @@ function Base.repeat(at::FlexibleSystem, n::NTuple{3})
    end
 
    bb = [c1 * n[1], c2 * n[2], c3 * n[3]]
-   return FlexibleSystem(particles, bb, at.boundary_conditions)
+   return FlexibleSystem(particles, bb, boundary_conditions(at))
 end
 
 Base.repeat(at::FlexibleSystem, n::Integer) = repeat(at, (n,n,n))
@@ -93,7 +92,6 @@ import Base.*
 
 
 
-
 """
 `rattle!(at, r::Float64) -> at`
 
@@ -104,10 +102,10 @@ not the same as choosing them uniform in cartesian coordinates!).
 function rattle!(at::FlexibleSystem, r::AbstractFloat)
    for i = 1:length(at.particles)
       p = at.particles[i]
-      𝐫i = p.position 
-      T = typeof(ustrip(𝐫i[1]))
+      𝐫ᵢ = p.position 
+      T = typeof(ustrip(𝐫ᵢ[1]))
       ui = randn(Vec3{T})
-      p_new = _set_position(p, 𝐫i + r * ui / norm(ui) * unit(𝐫i[1]))
+      p_new = _set_position(p, 𝐫ᵢ + r * ui / norm(ui) * unit(𝐫ᵢ[1]))
       at.particles[i] = p_new
    end
    return at
@@ -123,13 +121,6 @@ end
 #           cell = cell(at1),
 #           pbc = pbc(at1) )
 
-# append(at::Atoms, X::AbstractVector{<:JVec}) =
-#    Atoms( X = union(at.X, X),
-#           P = union(at.P, zeros(JVecF, length(X))),
-#           M = union(at.M, zeros(length(X))),
-#           Z = union(at.Z, zeros(Int, length(X))),
-#           cell = cell(at),
-#           pbc = pbc(at) )
 
 
 # """

test/test_bulk.jl

@@ -1,5 +1,5 @@
 
-using AtomsBuilder, Test, AtomsBase
+using AtomsBuilder, Test, AtomsBase, Unitful
 import JuLIP
 
 ##
@@ -38,5 +38,19 @@ end
 
 ##
 
-# not sure how to write a test for this: 
-rattle!(bulk(:C, cubic=true) * (2,3,4), 0.1)
+# not sure how to write a test for this, but at least it should execute
+sys1 = rattle!(bulk(:C, cubic=true) * (2,3,4), 0.1)
+sys2 = rattle!(bulk(:C) * (2,3,4), 0.1)
+
+X = position(sys1)
+Xnew = [ x .+ 0.01u"Å" for x in X ]
+sys3 = set_positions(sys1, Xnew)
+@test position(sys3) == Xnew
+
+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)
+