diff --git a/Project.toml b/Project.toml index 14e2362..560240d 100644 --- a/Project.toml +++ b/Project.toml @@ -3,11 +3,24 @@ uuid = "f5cc8831-eeb7-4288-8d9f-d6c1ddb77004" authors = ["Christoph Ortner and contributors"] version = "0.0.1-dev" +[deps] +AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a" +JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" + [compat] -julia = "1.6.7" +AtomsBase = "0.3.5" +JSON = "0.21.4" +LinearAlgebra = "1.9" +StaticArrays = "1.9" +Unitful = "1.19" +julia = "1.9" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +JuLIP = "945c410c-986d-556a-acb1-167a618e0462" [targets] -test = ["Test"] +test = ["Test", "JuLIP"] diff --git a/README.md b/README.md index 46f0c9e..aed0fa5 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,47 @@ # AtomsBuilder -[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaMolSim.github.io/AtomsBuilder.jl/stable/) -[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaMolSim.github.io/AtomsBuilder.jl/dev/) -[![Build Status](https://github.com/JuliaMolSim/AtomsBuilder.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaMolSim/AtomsBuilder.jl/actions/workflows/CI.yml?query=branch%3Amain) + +[![Build Status](https://github.com/JuliaMolSim/AtomsBuilder.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaMolSim/AtomsBuilder.jl/actions/workflows/CI.yml?query=branch%3Amain) + +This package provides utilities to build atomic structures. At the moment the functionality is very limited - see examples below. + + +## Preliminary Documentation + +Currently there are just two exported functions: +* `bulk` +* `rattle!` +In addition we overload +* `repeat` (with alias `*`) + +```julia +using AtomsBuilder + +# generate a diamond cubic bulk Si unit cell +at1 = bulk(:Si) +@show length(at1) + +# generate a minimal cubic Si cell (diamond cubic) +at2 = bulk(:Si, cubic=true) +@show length(at2) + +# repeat the cell 3 times in each coordinate direction +at3 = at2 * 3 +@show length(at3) + +# repeat the unit cell in only one direction +at4 = at2 * (3, 1, 1) +@show length(at3) + +# create a bulk supercell and then rattle the atoms +at5 = rattle!( bulk(:Si, cubic=true) * 3 ) +``` + +See `?bulk` and `?rattle!` for more information. + +## Contributions + +The current version of the package is essentially a copy-paste of a subset of functionality from an older package that is no longer developed. Contributions to expand the capabilities, improve the implementation, or entirely replace it are very welcome. Some packages that contain overlapping functionalities that could replace or add to `AtomsBuilder.jl` include +* [`Electrum.jl`](https://github.com/brainandforce/Electrum.jl) +* [`AtomsToolbox.jl`](https://github.com/rashidrafeek/AtomsToolbox.jl) diff --git a/data/asedata.json b/data/asedata.json new file mode 100644 index 0000000..e8cac54 --- /dev/null +++ b/data/asedata.json @@ -0,0 +1,801 @@ +{ + "symbols": [ + "X", + "H", + "He", + "Li", + "Be", + "B", + "C", + "N", + "O", + "F", + "Ne", + "Na", + "Mg", + "Al", + "Si", + "P", + "S", + "Cl", + "Ar", + "K", + "Ca", + "Sc", + "Ti", + "V", + "Cr", + "Mn", + "Fe", + "Co", + "Ni", + "Cu", + "Zn", + "Ga", + "Ge", + "As", + "Se", + "Br", + "Kr", + "Rb", + "Sr", + "Y", + "Zr", + "Nb", + "Mo", + "Tc", + "Ru", + "Rh", + "Pd", + "Ag", + "Cd", + "In", + "Sn", + "Sb", + "Te", + "I", + "Xe", + "Cs", + "Ba", + "La", + "Ce", + "Pr", + "Nd", + "Pm", + "Sm", + "Eu", + "Gd", + "Tb", + "Dy", + "Ho", + "Er", + "Tm", + "Yb", + "Lu", + "Hf", + "Ta", + "W", + "Re", + "Os", + "Ir", + "Pt", + "Au", + "Hg", + "Tl", + "Pb", + "Bi", + "Po", + "At", + "Rn", + "Fr", + "Ra", + "Ac", + "Th", + "Pa", + "U", + "Np", + "Pu", + "Am", + "Cm", + "Bk", + "Cf", + "Es", + "Fm", + "Md", + "No", + "Lr", + "Rf", + "Db", + "Sg", + "Bh", + "Hs", + "Mt", + "Ds", + "Rg", + "Cn", + "Nh", + "Fl", + "Mc", + "Lv", + "Ts", + "Og" + ], + "refstates": [ + null, + { + "symmetry": "diatom", + "d": 0.74 + }, + { + "symmetry": "atom" + }, + { + "symmetry": "bcc", + "a": 3.49 + }, + { + "symmetry": "hcp", + "c/a": 1.567, + "a": 2.29 + }, + { + "symmetry": "tetragonal", + "c/a": 0.576, + "a": 8.73 + }, + { + "symmetry": "diamond", + "a": 3.57 + }, + { + "symmetry": "diatom", + "d": 1.1 + }, + { + "symmetry": "diatom", + "d": 1.21 + }, + { + "symmetry": "diatom", + "d": 1.42 + }, + { + "symmetry": "fcc", + "a": 4.43 + }, + { + "symmetry": "bcc", + "a": 4.23 + }, + { + "symmetry": "hcp", + "c/a": 1.624, + "a": 3.21 + }, + { + "symmetry": "fcc", + "a": 4.05 + }, + { + "symmetry": "diamond", + "a": 5.43 + }, + { + "symmetry": "cubic", + "a": 7.17 + }, + { + "symmetry": "orthorhombic", + "b/a": 1.229, + "c/a": 2.339, + "a": 10.47 + }, + { + "symmetry": "orthorhombic", + "b/a": 0.718, + "c/a": 1.324, + "a": 6.24 + }, + { + "symmetry": "fcc", + "a": 5.26 + }, + { + "symmetry": "bcc", + "a": 5.23 + }, + { + "symmetry": "fcc", + "a": 5.58 + }, + { + "symmetry": "hcp", + "c/a": 1.594, + "a": 3.31 + }, + { + "symmetry": "hcp", + "c/a": 1.588, + "a": 2.95 + }, + { + "symmetry": "bcc", + "a": 3.02 + }, + { + "symmetry": "bcc", + "a": 2.88 + }, + { + "symmetry": "cubic", + "a": 8.89 + }, + { + "symmetry": "bcc", + "a": 2.87 + }, + { + "symmetry": "hcp", + "c/a": 1.622, + "a": 2.51 + }, + { + "symmetry": "fcc", + "a": 3.52 + }, + { + "symmetry": "fcc", + "a": 3.61 + }, + { + "symmetry": "hcp", + "c/a": 1.856, + "a": 2.66 + }, + { + "symmetry": "orthorhombic", + "b/a": 1.001, + "c/a": 1.695, + "a": 4.51 + }, + { + "symmetry": "diamond", + "a": 5.66 + }, + { + "symmetry": "rhombohedral", + "alpha": 54.1, + "a": 4.13 + }, + { + "symmetry": "hcp", + "c/a": 1.136, + "a": 4.36 + }, + { + "symmetry": "orthorhombic", + "b/a": 0.672, + "c/a": 1.307, + "a": 6.67 + }, + { + "symmetry": "fcc", + "a": 5.72 + }, + { + "symmetry": "bcc", + "a": 5.59 + }, + { + "symmetry": "fcc", + "a": 6.08 + }, + { + "symmetry": "hcp", + "c/a": 1.571, + "a": 3.65 + }, + { + "symmetry": "hcp", + "c/a": 1.593, + "a": 3.23 + }, + { + "symmetry": "bcc", + "a": 3.3 + }, + { + "symmetry": "bcc", + "a": 3.15 + }, + { + "symmetry": "hcp", + "c/a": 1.604, + "a": 2.74 + }, + { + "symmetry": "hcp", + "c/a": 1.584, + "a": 2.7 + }, + { + "symmetry": "fcc", + "a": 3.8 + }, + { + "symmetry": "fcc", + "a": 3.89 + }, + { + "symmetry": "fcc", + "a": 4.09 + }, + { + "symmetry": "hcp", + "c/a": 1.886, + "a": 2.98 + }, + { + "symmetry": "tetragonal", + "c/a": 1.076, + "a": 4.59 + }, + { + "symmetry": "tetragonal", + "c/a": 0.546, + "a": 5.82 + }, + { + "symmetry": "rhombohedral", + "alpha": 57.6, + "a": 4.51 + }, + { + "symmetry": "hcp", + "c/a": 1.33, + "a": 4.45 + }, + { + "symmetry": "orthorhombic", + "b/a": 0.659, + "c/a": 1.347, + "a": 7.27 + }, + { + "symmetry": "fcc", + "a": 6.2 + }, + { + "symmetry": "bcc", + "a": 6.05 + }, + { + "symmetry": "bcc", + "a": 5.02 + }, + { + "symmetry": "hcp", + "c/a": 1.619, + "a": 3.75 + }, + { + "symmetry": "fcc", + "a": 5.16 + }, + { + "symmetry": "hcp", + "c/a": 1.614, + "a": 3.67 + }, + { + "symmetry": "hcp", + "c/a": 1.614, + "a": 3.66 + }, + null, + { + "symmetry": "rhombohedral", + "alpha": 23.13, + "a": 9.0 + }, + { + "symmetry": "bcc", + "a": 4.61 + }, + { + "symmetry": "hcp", + "c/a": 1.588, + "a": 3.64 + }, + { + "symmetry": "hcp", + "c/a": 1.581, + "a": 3.6 + }, + { + "symmetry": "hcp", + "c/a": 1.573, + "a": 3.59 + }, + { + "symmetry": "hcp", + "c/a": 1.57, + "a": 3.58 + }, + { + "symmetry": "hcp", + "c/a": 1.57, + "a": 3.56 + }, + { + "symmetry": "hcp", + "c/a": 1.57, + "a": 3.54 + }, + { + "symmetry": "fcc", + "a": 5.49 + }, + { + "symmetry": "hcp", + "c/a": 1.585, + "a": 3.51 + }, + { + "symmetry": "hcp", + "c/a": 1.582, + "a": 3.2 + }, + { + "symmetry": "bcc", + "a": 3.31 + }, + { + "symmetry": "bcc", + "a": 3.16 + }, + { + "symmetry": "hcp", + "c/a": 1.615, + "a": 2.76 + }, + { + "symmetry": "hcp", + "c/a": 1.579, + "a": 2.74 + }, + { + "symmetry": "fcc", + "a": 3.84 + }, + { + "symmetry": "fcc", + "a": 3.92 + }, + { + "symmetry": "fcc", + "a": 4.08 + }, + { + "symmetry": "rhombohedral", + "alpha": 70.45, + "a": 2.99 + }, + { + "symmetry": "hcp", + "c/a": 1.599, + "a": 3.46 + }, + { + "symmetry": "fcc", + "a": 4.95 + }, + { + "symmetry": "rhombohedral", + "alpha": 57.14, + "a": 4.75 + }, + { + "symmetry": "sc", + "a": 3.35 + }, + null, + null, + null, + null, + { + "symmetry": "fcc", + "a": 5.31 + }, + { + "symmetry": "fcc", + "a": 5.08 + }, + { + "symmetry": "tetragonal", + "c/a": 0.825, + "a": 3.92 + }, + { + "symmetry": "orthorhombic", + "b/a": 1.736, + "c/a": 2.056, + "a": 2.85 + }, + { + "symmetry": "orthorhombic", + "b/a": 1.035, + "c/a": 1.411, + "a": 4.72 + }, + { + "symmetry": "monoclinic" + }, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null + ], + "rnn": [ + -1.0, + -1.0, + -1.0, + 3.02243, + 2.22873, + -1.0, + 1.54586, + -1.0, + -1.0, + -1.0, + 3.13248, + 3.66329, + 3.19823, + 2.86378, + 2.35126, + -1.0, + -1.0, + -1.0, + 3.71938, + 4.52931, + 3.94566, + 3.25752, + 2.89607, + 2.6154, + 2.49415, + -1.0, + 2.48549, + 2.49875, + 2.48902, + 2.55266, + 2.66, + -1.0, + 2.45085, + -1.0, + 3.53122, + -1.0, + 4.04465, + 4.84108, + 4.29921, + 3.55822, + 3.17748, + 2.85788, + 2.72798, + 2.70767, + 2.64627, + 2.68701, + 2.75065, + 2.89207, + 2.98, + -1.0, + -1.0, + -1.0, + 3.91893, + -1.0, + 4.38406, + 5.23945, + 4.34745, + 3.72861, + 3.64867, + 3.6416, + 3.63168, + -1.0, + -1.0, + 3.99238, + 3.57345, + 3.524, + 3.50263, + 3.48854, + 3.46905, + 3.44956, + 3.88202, + 3.44157, + 3.13374, + 2.86654, + 2.73664, + 2.73976, + 2.67994, + 2.71529, + 2.77186, + 2.885, + -1.0, + 3.41215, + 3.50018, + -1.0, + 3.35, + -1.0, + -1.0, + -1.0, + -1.0, + 3.75474, + 3.5921, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0 + ], + "masses": [ + 1.0, + 1.008, + 4.002602, + 6.94, + 9.0121831, + 10.81, + 12.011, + 14.007, + 15.999, + 18.998403163, + 20.1797, + 22.98976928, + 24.305, + 26.9815385, + 28.085, + 30.973761998, + 32.06, + 35.45, + 39.948, + 39.0983, + 40.078, + 44.955908, + 47.867, + 50.9415, + 51.9961, + 54.938044, + 55.845, + 58.933194, + 58.6934, + 63.546, + 65.38, + 69.723, + 72.63, + 74.921595, + 78.971, + 79.904, + 83.798, + 85.4678, + 87.62, + 88.90584, + 91.224, + 92.90637, + 95.95, + 97.90721, + 101.07, + 102.9055, + 106.42, + 107.8682, + 112.414, + 114.818, + 118.71, + 121.76, + 127.6, + 126.90447, + 131.293, + 132.90545196, + 137.327, + 138.90547, + 140.116, + 140.90766, + 144.242, + 144.91276, + 150.36, + 151.964, + 157.25, + 158.92535, + 162.5, + 164.93033, + 167.259, + 168.93422, + 173.054, + 174.9668, + 178.49, + 180.94788, + 183.84, + 186.207, + 190.23, + 192.217, + 195.084, + 196.966569, + 200.592, + 204.38, + 207.2, + 208.9804, + 208.98243, + 209.98715, + 222.01758, + 223.01974, + 226.02541, + 227.02775, + 232.0377, + 231.03588, + 238.02891, + 237.04817, + 244.06421, + 243.06138, + 247.07035, + 247.07031, + 251.07959, + 252.083, + 257.09511, + 258.09843, + 259.101, + 262.11, + 267.122, + 268.126, + 271.134, + 270.133, + 269.1338, + 278.156, + 281.165, + 281.166, + 285.177, + 286.182, + 289.19, + 289.194, + 293.204, + 293.208, + 294.214 + ] +} diff --git a/data/import_asedata.jl b/data/import_asedata.jl new file mode 100644 index 0000000..eba4a18 --- /dev/null +++ b/data/import_asedata.jl @@ -0,0 +1,92 @@ +# +# ====================================================================== +# Import some chemistry and materials science data tables from +# ASE. Without this data, JuLIP can do very little! +# ====================================================================== +# +using PyCall, JSON +@pyimport ase.data as ase_data +asedata = Dict( + :symbols => ase_data.chemical_symbols, + :masses => ase_data.atomic_masses, + :refstates => ase_data.reference_states + ) + + +write(@__DIR__() * "/asedata.json", JSON.json(asedata, 0)) + +# ====================================================================== + + +# NOTE: +# some other data that we could consider adding +# asedata.atomic_numbers +# ase_data.atomic_names +# ase_data.covalent_radii +# ase_data.ground_state_magnetic_moments +# ase_data.vdw_radii +# can we get some more stuff like electron affinity somewhere? + + +# function rnn_old(species::Symbol) +# X = positions(bulk(species) * 2) +# return minimum( norm(X[n]-X[m]) for n = 1:length(X) for m = n+1:length(X) ) +# end +# +# _rnn = fill(-1.0, length(_symbols)) +# for n = 2:length(_symbols) +# z = n-1 +# try +# _rnn[n] = rnn_old(JuLIP.Chemistry.chemical_symbol(z)) +# catch +# end +# end + + +# ====================================================================== + + +# get access to the atomic numbers +ase_emt = pyimport("ase.calculators.emt") +ase_data = pyimport("ase.data") + +function emt_default_parameters() + # import everything we need from ASE + Bohr = ase_emt.Bohr + # get_atomic_numbers(id::AbstractString) = ase_data.atomic_numbers[id] + ase_parameters = ase_emt.parameters + beta = ase_emt.beta + # convert ASE style Dict to a new Dict where parameters + # are stored in a type instead of tuple. + params = Dict{String, Dict{String,Float64}}() + maxseq = maximum([par[2] for par in values(ase_parameters)]) * Bohr # ✓ + rc = beta * maxseq * 0.5 * (sqrt(3) + sqrt(4)) # ✓ + rr = rc * 2 * sqrt(4) / (sqrt(3) + sqrt(4)) # ✓ + acut = log(9999.0) / (rr - rc) # ✓ + + for (key, p) in ase_parameters # p is a tuple of parameters, key a species + s0 = p[2] * Bohr + eta2 = p[4] / Bohr + kappa = p[5] / Bohr + x = eta2 * beta * s0 + gamma1 = 0.0 + gamma2 = 0.0 + for (i, n) in enumerate([12; 6; 24]) + r = s0 * beta * sqrt(i) + x = n / (12.0 * (1.0 + exp(acut * (r - rc)))) + gamma1 += x * exp(-eta2 * (r - beta * s0)) + gamma2 += x * exp(-kappa / beta * (r - beta * s0)) + end + params[key] = Dict{String, Float64}( + "E0" => p[1], "s0" => s0, "V0" => p[3], "η2" => eta2, + "κ" => kappa, "λ" => p[6] / Bohr, "n0" => p[7] / Bohr^3, + "γ1" => gamma1, "γ2" => gamma2 ) + par = params[key] + par["Cpair"] = -0.5 * par["V0"] / par["γ2"] / par["n0"] + end + + return Dict("params" => params, "acut" => acut, "rc" => rc, "beta" => beta) +end + +emt_data = emt_default_parameters() +write(@__DIR__() * "/emt.json", JSON.json(emt_data, 0)) diff --git a/src/AtomsBuilder.jl b/src/AtomsBuilder.jl index 3d58e66..0171037 100644 --- a/src/AtomsBuilder.jl +++ b/src/AtomsBuilder.jl @@ -1,5 +1,15 @@ module AtomsBuilder -# Write your package code here. +import AtomsBase +using StaticArrays: SMatrix, SVector + +const Mat3{T} = SMatrix{3, 3, T} +const Vec3{T} = SVector{3, T} + +include("chemistry.jl") + +include("utils.jl") + +include("bulk.jl") end diff --git a/src/bulk.jl b/src/bulk.jl new file mode 100644 index 0000000..dacdd59 --- /dev/null +++ b/src/bulk.jl @@ -0,0 +1,76 @@ +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), + :bcc => ( [ [0.0,0.0,0.0], ], + [-1 1 1; 1 -1 1; 1 1 -1], 0.5), + :diamond => ( [ [0.0, 0.0, 0.0], [0.5, 0.5, 0.5] ], + [0 1 1; 1 0 1; 1 1 0], 0.5) +) + +const _cubic_cells = Dict( # (positions, factor of a) + :fcc => ( [ [0 0 0], [0 1 1], [1 0 1], [1 1 0] ], 0.5 ), + :bcc => ( [ [0 0 0], [1 1 1] ], 0.5 ), + :diamond => ( [ [0 0 0], [1 1 1], [0 2 2], [1 3 3], [2 0 2], + [3 1 3], [2 2 0], [3 3 1] ], 0.25 ) +) + +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)] + C = Matrix(1.0I, 3, 3) / scal + else + X, C, scal = _unit_cells[symmetry(sym)] + end + a === nothing && (a = 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 + + +function _bulk_hcp(sym::Symbol; a=nothing, c=nothing, T=Float64) + D = Chemistry.refstate(sym) + a === nothing && (a = D["a"]) + c === nothing && (c = a * D["c/a"]) + return [ a * Vec3{T}(0.0, 0.0, 0.0), + a * Vec3{T}(0.0, 1 / sqrt(3), c / a / 2) ], + a * Mat3{T}( [1.0, -1/2, 0.0, 0.0, sqrt(3)/2, 0.0, 0.0, 0.0, c/a] ) +end + + +_convert_pbc(pbc::NTuple{3, Bool}) = pbc +_convert_pbc(pbc::Bool) = (pbc, pbc, pbc) + +""" +`bulk(sym)` : generates a `FlexibleSystem` unit cell for a bulk +crystal structure. If `sym` is a chemical symbol then the phase and +lattice constant are taken from a database that is consistent with ASE. +If `sym` is one of `[:fcc, :bcc, :diamond, :hcp]` then one needs to +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) + 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) + 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 + # it again ... + # (x !== nothing || y !== nothing || z !== nothing) && rotate!(at, x=x, y=y, z=z) + return at +end diff --git a/src/chemistry.jl b/src/chemistry.jl index ebd3cf3..468fa5d 100644 --- a/src/chemistry.jl +++ b/src/chemistry.jl @@ -1,6 +1,7 @@ module Chemistry using JSON +using Unitful # -------------------------------------------------------------------- # Load data and prepare it a bit ... @@ -9,10 +10,18 @@ const ase_data_path = joinpath(@__DIR__(), "..", "data", "asedata.json") const ase_data = JSON.parsefile(ase_data_path) -const _rnn = [Float64(d) for d in ase_data["rnn"]] -const _masses = [Float64(m) for m in ase_data["masses"]] +const _rnn = [Float64(d) * u"Å" for d in ase_data["rnn"]] +const _masses = [Float64(m) * u"u" for m in ase_data["masses"]] const _refstates = Dict{String, Any}[ d == nothing ? Dict{String, Any}() : d for d in ase_data["refstates"]] +for D in _refstates + if haskey(D, "a") + D["a"] *= u"Å" + end + if haskey(D, "d") + D["d"] *= u"Å" + end +end const _symbols = Dict{Int, Symbol}() const _numbers = Dict{Symbol, Int}() @@ -28,6 +37,7 @@ end atomic_number(sym::Symbol) = _numbers[sym] chemical_symbol(z::Integer) = _symbols[z] +chemical_symbol(sym::Symbol) = sym atomic_mass(z::Integer) = _masses[z+1] atomic_mass(sym::Symbol) = atomic_mass(atomic_number(sym)) diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..e241962 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,130 @@ + +using AtomsBase +using AtomsBase: Atom, FlexibleSystem, Periodic +using Unitful: unit, ustrip +using LinearAlgebra: norm + +export rattle! + +""" +Helper function to convert construct a FlexibleSystem from + a list of positions, elements, cell matrix and pbc tuple +""" +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 ] + bc = [ (pbc[i] ? Periodic() : nothing) for i = 1:3 ] + bb = [cell[i, :] for i = 1:3] + return FlexibleSystem(atoms; + bounding_box = bb, + boundary_conditions = bc) +end + +_set_position(x::Atom, 𝐫) = Atom(; position = 𝐫, + velocity = x.velocity, + atomic_mass = x.atomic_mass, + atomic_number = x.atomic_number, + atomic_symbol = x.atomic_symbol) + + +""" +``` +repeat(at, n::NTuple{3}) +repeat(at, n::Integer) +``` + +Takes a structure and repeats it n_j times +into the j-th cell-vector direction. For example, +``` +at = repeat(bulk(:C), (3,2,4)) +``` +creates 3 x 2 x 4 unit cells of carbon. + +The same can be achieved by `*`: +``` +at = bulk(:) * (3, 2, 4) +``` +""" +function Base.repeat(at::FlexibleSystem, n::NTuple{3}) + c1, c2, c3 = at.bounding_box + + particles = eltype(at.particles)[] + for a in CartesianIndices( (1:n[1], 1:n[2], 1:n[3]) ) + b = c1 * (a[1]-1) + c2 * (a[2]-1) + c3 * (a[3]-1) + for i in 1:length(at) + p_i = at.particles[i] + p_new = _set_position(p_i, b + p_i.position) + push!(particles, p_new) + end + end + + bb = [c1 * n[1], c2 * n[2], c3 * n[3]] + return FlexibleSystem(particles, bb, at.boundary_conditions) +end + +Base.repeat(at::FlexibleSystem, n::Integer) = repeat(at, (n,n,n)) + +import Base.* +*(at::FlexibleSystem, n) = repeat(at, n) +*(n, at::FlexibleSystem) = repeat(at, n) + + + + +""" +`rattle!(at, r::Float64) -> at` + +Randomly perturbs the atom positions within a ball of radius `r`. The perturbation +is uniform in angular component, and uniform in radial component. (Note this is +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])) + ui = randn(Vec3{T}) + p_new = _set_position(p, 𝐫i + r * ui / norm(ui) * unit(𝐫i[1])) + at.particles[i] = p_new + end + return at +end + + + +# union(at1::Atoms, at2::Atoms) = +# Atoms( X = union(at1.X, at2.X), +# P = union(at1.P, at2.P), +# M = union(at1.M, at2.M), +# Z = union(at1.Z, at2.Z), +# 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) ) + + +# """ +# `deleteat!(at::Atoms, n) -> at`: + +# returns the same atoms object `at`, but with the atom(s) specified by `n` +# removed. +# """ +# function Base.deleteat!(at::Atoms, n) +# deleteat!(at.X, n) +# deleteat!(at.P, n) +# deleteat!(at.M, n) +# deleteat!(at.Z, n) +# update_data!(at, Inf) +# JuLIP.reset_clamp!(at) +# return at +# end + + diff --git a/test/runtests.jl b/test/runtests.jl index a5c327a..590db00 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,4 +3,7 @@ using Test @testset "AtomsBuilder.jl" begin # Write your tests here. + @testset "bulk" begin include("test_bulk.jl"); end end + + diff --git a/test/test_bulk.jl b/test/test_bulk.jl new file mode 100644 index 0000000..10d06e6 --- /dev/null +++ b/test/test_bulk.jl @@ -0,0 +1,42 @@ + +using AtomsBuilder, Test, AtomsBase +import JuLIP + +## + +@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 + +for sym in [:Si, :Ge, :W, :Ti] + at = bulk(sym) + at_f = FlexibleSystem(JuLIP.bulk(sym)) + @test _compare_at(at, at_f) + + at = bulk(sym, cubic=true) + at_f = FlexibleSystem(JuLIP.bulk(sym, cubic=true)) + @test _compare_at(at, at_f) + + 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) + + at = bulk(sym, cubic=true) * nn + at_f = FlexibleSystem(JuLIP.bulk(sym, cubic=true) * nn) + @test _compare_at(at, at_f) + end +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