Update to AtomsBase 0.5

.github/workflows/CI.yml

@@ -19,7 +19,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - '1.9'
           - 'nightly'
         os:
           - ubuntu-latest
@@ -44,6 +44,7 @@ jobs:
     runs-on: ubuntu-latest
     permissions:
       contents: write
+      statuses: write
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2

CondaPkg.toml

@@ -1,4 +1,4 @@
 channels = ["conda-forge"]
 
 [deps]
-ase = ">=3.22,<4"
+ase = ">=3.23,<3.24"

Project.toml

@@ -1,7 +1,7 @@
 name = "ASEconvert"
 uuid = "3da9722f-58c2-4165-81be-b4d7253e8fd2"
 authors = ["Michael F. Herbst <[email protected]>"]
-version = "0.1.7"
+version = "0.1.8"
 
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
@@ -13,15 +13,15 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [compat]
-AtomsBaseTesting = "0.1"
-AtomsBase = "0.3"
+AtomsBaseTesting = "0.4"
+AtomsBase = "0.5"
 AtomsCalculators = "0.2"
 CondaPkg = "0.2"
 PeriodicTable = "1"
 PythonCall = "0.9"
 Unitful = "1"
 UnitfulAtomic = "1"
-julia = "1.6"
+julia = "1.9"
 
 [extras]
 AtomsBaseTesting = "ed7c10db-df7e-4efa-a7be-4f4190f7f227"

docs/Project.toml

@@ -3,6 +3,6 @@ ASEconvert = "3da9722f-58c2-4165-81be-b4d7253e8fd2"
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
 AtomsBuilder = "f5cc8831-eeb7-4288-8d9f-d6c1ddb77004"
 AtomsCalculators = "a3e0e189-c65a-42c1-833c-339540406eb1"
-AtomsIO = "1692102d-eeb4-4df9-807b-c9517f998d44"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ExtXYZ = "352459e4-ddd7-4360-8937-99dcb397b478"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"

docs/make.jl

@@ -12,20 +12,23 @@ using Documenter
 DocMeta.setdocmeta!(ASEconvert, :DocTestSetup, :(using ASEconvert); recursive=true)
 
 makedocs(;
-         modules=[ASEconvert],
-         authors="Michael F. Herbst <[email protected]> and contributors",
-         repo="https://github.com/mfherbst/ASEconvert.jl/blob/{commit}{path}#{line}",
-         sitename="ASEconvert.jl",
-         format=Documenter.HTML(;
-                                prettyurls=get(ENV, "CI", "false") == "true",
-                                canonical="https://mfherbst.github.io/ASEconvert.jl",
-                                edit_link="master",
-                                assets=String[]),
-         pages=[
-             "Home" => "index.md",
-             "apireference.md",
-         ])
+    modules=[ASEconvert],
+    authors="Michael F. Herbst <[email protected]> and contributors",
+    repo="https://github.com/mfherbst/ASEconvert.jl/blob/{commit}{path}#{line}",
+    sitename="ASEconvert.jl",
+    format=Documenter.HTML(;
+        prettyurls=get(ENV, "CI", "false") == "true",
+        canonical="https://mfherbst.github.io/ASEconvert.jl",
+        edit_link="master",
+        assets=String[],
+    ),
+    pages=[
+        "Home" => "index.md",
+        "apireference.md",
+    ],
+)
 
 deploydocs(;
-           repo="github.com/mfherbst/ASEconvert.jl",
-           devbranch="master")
+    repo="github.com/mfherbst/ASEconvert.jl",
+    devbranch="master",
+)

docs/src/index.md

@@ -54,18 +54,16 @@ scfres.energies
 ```
 
 ## Conversion from AtomsBase to ASE
-
 ```@example extxyz
 using ASEconvert
-using AtomsIO
+using ExtXYZ
 
 # Read an extxyz file using AtomsIO.jl.
-system = load_system("Mn3Si.extxyz")
+system = ExtXYZ.Atoms(ExtXYZ.read_frame("Mn3Si.extxyz"))
 ```
-
-This example uses [AtomsIO](https://github.com/mfherbst/AtomsIO.jl)
+This example uses [ExtXYZ](https://github.com/libAtoms/ExtXYZ.jl)
 to read the extended XYZ file file `Mn3Si.extxyz`. The data is returned
-as a subtype of `AtomsBase.AbstractSystem`
+as a subtype of [`AtomsBase.AbstractSystem`](https://juliamolsim.github.io/AtomsBase.jl/)
 (in this case an `ExtXYZ.Atoms` from [ExtXYZ](https://github.com/libAtoms/ExtXYZ.jl)).
 We can thus directly convert this system to an `ase.Atoms` using [`convert_ase`](@ref)
 and write it again as an ASE json file
@@ -74,6 +72,9 @@ and write it again as an ASE json file
 ase.io.write("out.json", convert_ase(system));
 ```
 
+For a more convenient and feature-rich way of reading and writing atomic
+structures in julia see [AtomsIO](https://github.com/mfherbst/AtomIO.jl).
+
 ## Employing ASE calculators in Julia
 
 ```@example calculators

src/ASEconvert.jl

@@ -5,8 +5,8 @@ using Unitful
 using UnitfulAtomic
 
 export ase
-export pyconvert, pytuple  # Reexport from PythonCall
-export AbstractSystem      # Reexport from AtomsBase
+export pyconvert       # Reexport from PythonCall
+export AbstractSystem  # Reexport from AtomsBase
 
 export ASEcalculator
 export convert_ase

src/ase_conversions.jl

@@ -1,30 +1,53 @@
 import PeriodicTable
 
+# Some general notes:
+#    - ASE has a cell object, but does neither make the distinction
+#      between isolated and periodic cells nor do cells store
+#      the boundary conditions like we do in AtomsBase. So we will
+#      *not* support conversion to and from cells for now.
+#    - Isotope handling is not yet fully implemented in AtomsBase,
+#      so this will be ignored in this interface.
+#    - Atomic numbers and atomic symbols cannot differ in ASE,
+#      so this is not supported in the interface.
+
+# TODO One could maybe also implement this as a wrapper, where data
+#      is not copied, but updated in-place in the python object.
 
 # For ASE units, see https://wiki.fysik.dtu.dk/ase/ase/units.html
 # In particular note that uTime = u"Å" * sqrt(u"u" / u"eV") and thus
 const uVelocity = sqrt(u"eV" / u"u")
 
 
 function ase_to_system(S::Type{<:AbstractSystem}, ase_atoms::Py)
-    box = [pyconvert(Vector, ase_atoms.cell[i])u"Å" for i = 0:2]
+    cell_vectors = [pyconvert(Vector, ase_atoms.cell[i])u"Å" for i = 0:2]
+    periodicity  = pyconvert(Vector, ase_atoms.pbc)
+
+    # Convention in ASE isolated cells have zero (undefined) cell vectors
+    if all(iszero, cell_vectors)
+        if any(periodicity)
+            @warn "Ignoring ASE pbc settings which ASE cell vectors are zero (undefined)"
+        end
+        cϵll = IsolatedCell(3, typeof(1.0u"Å"))
+    else
+        cϵll = PeriodicCell(; cell_vectors, periodicity)
+    end
 
     atnums     = pyconvert(Vector, ase_atoms.get_atomic_numbers())
-    atsyms     = pyconvert(Vector, ase_atoms.symbols)
     atmasses   = pyconvert(Vector, ase_atoms.get_masses())
     positions  = pyconvert(Matrix, ase_atoms.get_positions())
     velocities = pyconvert(Matrix, ase_atoms.get_velocities())
     magmoms    = pyconvert(Vector, ase_atoms.get_initial_magnetic_moments())
     charges    = pyconvert(Vector, ase_atoms.get_initial_charges())
     ase_info   = pyconvert(Dict{String,Any}, ase_atoms.info)
 
-    atoms = map(1:length(atnums)) do i
-        AtomsBase.Atom(atnums[i], positions[i, :]u"Å", velocities[i, :] * uVelocity;
-                       atomic_symbol=Symbol(atsyms[i]),
-                       atomic_number=atnums[i],
-                       atomic_mass=atmasses[i]u"u",
-                       magnetic_moment=magmoms[i],
-                       charge=charges[i]u"e_au")
+    particles = map(1:length(atnums)) do i
+        Atom(atnums[i],
+             positions[i, :]u"Å",
+             velocities[i, :] * uVelocity;
+             species=ChemicalSpecies(atnums[i]),
+             mass=atmasses[i]u"u",
+             magnetic_moment=magmoms[i],
+             charge=charges[i]u"e_au")
     end
 
     # Parse extra data in info struct
@@ -37,8 +60,7 @@
         end
     end
 
-    bcs = [p ? Periodic() : DirichletZero() for p in pyconvert(Vector, ase_atoms.pbc)]
-    PythonCall.pyconvert_return(atomic_system(atoms, box, bcs; info...))
+    PythonCall.pyconvert_return(FlexibleSystem(particles, cϵll; info...))
 end
 
 """
@@ -50,33 +72,36 @@
 """
 function convert_ase(system::AbstractSystem{D}) where {D}
     D != 3 && @warn "1D and 2D systems not yet fully supported."
-
     n_atoms = length(system)
-    pbc     = map(isequal(Periodic()), boundary_conditions(system))
-    numbers = atomic_number(system)
-    masses  = ustrip.(u"u", atomic_mass(system))
-
-    symbols_match = [
-        PeriodicTable.elements[atnum].symbol == string(atomic_symbol(system, i))
-        for (i, atnum) in enumerate(numbers)
-    ]
-    if !all(symbols_match)
-        @warn("Mismatch between atomic numbers and atomic symbols, which is not " *
-              "supported in ASE. Atomic numbers take preference.")
+
+    ase_cell = zeros(3, 3)
+    if cell(system) isa IsolatedCell
+        pbc = [false, false, false]
+    else
+        pbc = periodicity(system)
+        for (i, v) in enumerate(cell_vectors(system))
+            ase_cell[i, 1:D] = ustrip.(u"Å", v)
+        end
     end
 
-    cell = zeros(3, 3)
-    for (i, v) in enumerate(bounding_box(system))
-        cell[i, 1:D] = ustrip.(u"Å", v)
+    for at = 1:n_atoms
+        spec = species(system, at)
+        if spec.n_neutrons ≥ 0
+            @warn("Atom $at has a non-default neutron count. This information " *
+                  "is dropped when converting to ASE.")
+        end
     end
 
+    numbers = [atomic_number(species(system, at)) for at = 1:n_atoms]
+    masses  = [ustrip(u"u", mass(system, at))     for at = 1:n_atoms]
+
     positions = zeros(n_atoms, 3)
     for at = 1:n_atoms
         positions[at, 1:D] = ustrip.(u"Å", position(system, at))
     end
 
     velocities = nothing
-    if !ismissing(velocity(system))
+    if n_atoms > 0 && !ismissing(velocity(system, 1))
         velocities = zeros(n_atoms, 3)
         for at = 1:n_atoms
             velocities[at, 1:D] = ustrip.(uVelocity, velocity(system, at))
@@ -90,7 +115,7 @@
     charges = nothing
     magmoms = nothing
     for key in atomkeys(system)
-        if key in (:position, :velocity, :atomic_symbol, :atomic_number, :atomic_mass)
+        if key in (:position, :velocity, :species, :mass)
             continue  # Already dealt with
         elseif key == :charge
             charges = ustrip.(u"e_au", system[:, :charge])
@@ -104,7 +129,7 @@
     # Map extra system properties
     info = Dict{String, Any}()
     for (k, v) in pairs(system)
-        if k in (:bounding_box, :boundary_conditions)
+        if k in (:cell_vectors, :periodicity)
             continue
         elseif k in (:charge, )
             info[string(k)] = ustrip(u"e_au", v)
@@ -117,9 +142,7 @@
     end
 
     ase.Atoms(; positions, numbers, masses, magmoms, charges,
-              cell, pbc, velocities, info)
+              cell=ase_cell, pbc, velocities, info)
 end
 
 # TODO Could have a convert_ase(Vector{AbstractSystem}) to make an ASE trajectory
-# TODO Could have a convert_ase(Vector{Vector{Unitful}}) to make an ASE cell
-# TODO Could have a way to make an ASE calculator from an InteratomicPotential object

src/atoms_calculators.jl

@@ -46,7 +46,7 @@ end
     fdata = calc.calculator.get_forces(convert_ase(system))
     # We need to convert from Python row major to Julia column major
     # and from there to Vector with SVector{3,Float64} element.
-    # We do reallocation in the end to ensure we have the correct memory allignement
+    # We do reallocation in the end to ensure we have the correct memory alignment
     tmp  = pyconvert(Array, fdata)
     FT   = AtomsCalculators.promote_force_type(system, calc)
     tmp2 = reinterpret(FT, tmp')