From 288377c3bbc293a18f04bdf6b3d6a22e21b3afe2 Mon Sep 17 00:00:00 2001 From: JuliaMolSim Date: Wed, 22 May 2024 09:50:06 -0700 Subject: [PATCH] trying to move to strict use of interface --- Project.toml | 4 ++-- src/bulk.jl | 17 ++++++-------- src/chemistry.jl | 16 ++++++++----- src/utils.jl | 59 ++++++++++++++++++++--------------------------- test/test_bulk.jl | 20 +++++++++++++--- 5 files changed, 61 insertions(+), 55 deletions(-) diff --git a/Project.toml b/Project.toml index 016b2ea..bc7a6dd 100644 --- a/Project.toml +++ b/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" diff --git a/src/bulk.jl b/src/bulk.jl index dacdd59..2111df3 100644 --- a/src/bulk.jl +++ b/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 diff --git a/src/chemistry.jl b/src/chemistry.jl index 468fa5d..06539c5 100644 --- a/src/chemistry.jl +++ b/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 \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index 964ca70..5734a9a 100644 --- a/src/utils.jl +++ b/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) ) # """ diff --git a/test/test_bulk.jl b/test/test_bulk.jl index 10d06e6..dc3cf12 100644 --- a/test/test_bulk.jl +++ b/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) \ No newline at end of file +# 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) +