Update documentation with new IO functionality (#301)

* Update documentation with new IO functionality

* Reduce documentation dependencies

* Update documentation.yml

* Add back StaticArrays for doctest

* Fix errors in doc build

.github/workflows/documentation.yml

@@ -32,7 +32,7 @@ jobs:
           Pkg.Registry.add(RegistrySpec(url="https://github.com/JuliaMolSim/MolSim"))
           Pkg.Registry.add(RegistrySpec(url="https://github.com/ACEsuit/ACEregistry"))
           Pkg.Registry.add(RegistrySpec(name="General"))
-        shell: julia {0}
+        shell: julia --color=yes {0}
 
       - name: Install ase
         run: python3 -m pip install ase
@@ -42,12 +42,12 @@ jobs:
           using Pkg
           Pkg.develop([(;path=pwd()), (;name="NQCModels")])
           Pkg.instantiate()
-        shell: julia --project=docs/ {0}
+        shell: julia --color=yes --project=docs/ {0}
 
       - name: Build and deploy
+        run: julia --color=yes --project=docs/ docs/make.jl
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
           GKSwstype: "100" # https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988
-        run: julia --project=docs/ docs/make.jl
 

docs/Project.toml

@@ -1,6 +1,7 @@
 [deps]
+ASEconvert = "3da9722f-58c2-4165-81be-b4d7253e8fd2"
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+AtomsIO = "1692102d-eeb4-4df9-807b-c9517f998d44"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 CubeLDFAModel = "56146587-2a70-4d59-809c-88d5ba27a8d3"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -9,7 +10,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-JuLIP = "945c410c-986d-556a-acb1-167a618e0462"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NNInterfaces = "07253886-9aba-4a5d-b3d5-29b8b100a664"
 NQCBase = "78c76ebc-5665-4934-b512-82d81b5cbfb7"
@@ -18,13 +18,12 @@ NQCDynamics = "36248dfb-79eb-4f4d-ab9c-e29ea5f33e14"
 NQCModels = "c814dc9f-a51f-4eaf-877f-82eda4edad48"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 RingPolymerArrays = "81cd853e-d334-476c-b156-f95acc8a32dc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
-UnitfulRecipes = "42071c24-d89e-48dd-8a24-8a12d9b8861f"
 
 [compat]
 Documenter = "0.27"

docs/src/NQCModels/analyticmodels.md

@@ -2,14 +2,6 @@
 
 This page plots many of the analytic models included in `NQCDynamics`.
 
-!!! tip 
-    To produce the plots we use two of Julia's plotting options [Plots](http://docs.juliaplots.org/latest/)
-    and [Makie](https://makie.juliaplots.org/dev/).
-    Plots has a mature recipe system that allows us to define custom plots for the 1D models
-    but we use Makie to produce the more complex images.
-    Each of these has their pros and cons and if you are interested in producing plots
-    using Julia you should visit their documentation to decide which is best for you.
-
 ## [`AdiabaticModels`](@ref NQCModels.AdiabaticModels)
 These models are used for classical dynamics and provide a single potential energy surface.
 
@@ -22,8 +14,7 @@ plot(-10:0.1:10, Harmonic())
 ### [`DiatomicHarmonic`](@ref)
 
 ```@example
-using NQCModels
-using Plots
+using NQCModels, Plots
 
 model = DiatomicHarmonic(r₀=10.0)
 f(x,y) = potential(model, [x y 0])
@@ -35,26 +26,23 @@ ylabel!("y coordinate /a₀")
 ### [`DarlingHollowayElbow`](@ref)
 
 ```@example
-using NQCModels
+using NQCModels, Plots
 using NQCBase: eV_to_au
-using CairoMakie
 
 model = DarlingHollowayElbow()
 V(x,z) = potential(model, [x, z])
 
 x = range(-0.5, 3.5, length=200)
 z = range(-0.5, 4.5, length=200)
 
-f = Figure()
-ax = Axis(f[1,1], xlabel="Bond length /a₀", ylabel="Surface molecule distance /a₀")
+plot(
+    xlabel="Bond length /a₀",
+    ylabel="Surface molecule distance /a₀",
+    xlims=(-0.5, 3.5),
+    ylims=(-0.5, 4.5)
+)
 
-levels = eV_to_au.([-0.1, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1])
-contourf!(ax, x, z, V, levels=levels)
-contour!(ax, x, z, V, levels=levels, color=:black)
-Colorbar(f[1,2], limits=(-0.1, 1.1))
-xlims!(-0.5, 3.5)
-ylims!(-0.5, 4.5)
-f
+contourf!(x, z, V)
 ```
 
 ## [`DiabaticModels`](@ref NQCModels.DiabaticModels)
@@ -93,22 +81,15 @@ plot(-5:0.1:5, DoubleWell(); coupling=true)
 ```
 ### [`GatesHollowayElbow`](@ref)
 ```@example
-using NQCModels
-using CairoMakie
+using NQCModels, Plots
 
 model = GatesHollowayElbow()
 v1(x,z) = potential(model, [x z])[1,1]
 v2(x,z) = potential(model, [x z])[2,2]
-coupling(x,z) = potential(model, [x z])[1,2]
 
 x = range(-0.5, 4.0, length=200)
 z = range(-0.5, 4.0, length=200)
 
-f = Figure()
-ax = Axis(f[1,1], xlabel="x coordinate", ylabel="z coordinate", limits=(-0.5, 4.0, -0.5, 4.0))
-
-contour!(ax, x, z, coupling, color=:black, levels=10, label="V12")
-contour!(ax, x, z, v1, color=:blue, levels=0:0.01:0.1, label="V11")
-contour!(ax, x, z, v2, color=:red, levels=0:0.01:0.1, label="V22")
-f
+contour(x, z, v1, color=:blue, levels=0:0.01:0.1, label="V11", colorbar=false)
+contour!(x, z, v2, color=:red, levels=0:0.01:0.1, label="V22")
 ```

docs/src/NQCModels/neuralnetworkmodels.md

@@ -3,7 +3,7 @@
 Using the [ASE interface](@ref) we can directly use models trained using
 [SchNetPack](https://github.com/atomistic-machine-learning/schnetpack).
 
-!!! danger
+!!! warning
 
     The examples on this page do not run during the documentation build due to `schnetpack`
     causing segfaults when installed in the build environment.
@@ -35,7 +35,7 @@ h2.set_calculator(calc)
 ```
 
 We can obtain the energies and forces from `ase` directly in the usual way, converting
-them to atomic units using [UnifulAtomic](https://github.com/sostock/UnitfulAtomic.jl).
+them to atomic units using [UnitfulAtomic](https://github.com/sostock/UnitfulAtomic.jl).
 ```julia-repl
 using Unitful, UnitfulAtomic;
 austrip(h2.get_total_energy() * u"eV")

docs/src/atoms.md

@@ -26,54 +26,59 @@ auno = [au; no]
 Atoms(auno)
 ```
 
-## [Reading and writing atomic structures](@id reading-and-writing)
+# [Manipulating atomic structures with AtomsBase.jl](@id atoms-base)
 
-When using a complex system however, it is likely more effective to read structures directly
-from a file.
-We provide two ways to do this, either using the 
-[ExtXYZ.jl](https://github.com/libAtoms/ExtXYZ.jl) package which works for xyz files.
-Or instead there is a conversion to and from the `ase.Atoms` type which can be used
-when [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) is loaded.
+[AtomsBase](https://github.com/JuliaMolSim/AtomsBase.jl) provides a convenient format for
+representing atomic geometries, facilitating interoperability between a collection of
+different packages.
 
-First we can use `ase` to build a system and write it to an `.xyz` file.
-```@example atoms
-using PyCall
+When working with NQCDynamics, the most useful packages are [AtomsIO](https://github.com/mfherbst/AtomsIO.jl)
+for reading and writing structures and trajectories, and [ASEconvert](https://github.com/mfherbst/ASEconvert.jl)
+for working with [ASE](https://wiki.fysik.dtu.dk/ase/index.html) from within Julia.
 
-build = pyimport("ase.build")
+## Using ASE with ASEconvert.jl
 
-slab = build.fcc100("Al", size=(2, 2, 3))
-build.add_adsorbate(slab, "Au", 1.7, "hollow")
-slab.center(axis=2, vacuum=4.0)
+This example shows how ASEconvert can be used to build a structure, then convert
+from the ASE format into an AtomsBase compatible system:
 
-slab.write("slab.xyz")
-```
+```julia
+using ASEconvert
 
-Now we can read it in with the [`read_extxyz`](@ref NQCBase.read_extxyz)
-function.
-```@repl atoms
-atoms, positions, cell = read_extxyz("slab.xyz")
-atoms
-positions
-cell
+# Make a silicon supercell using ASE
+atoms_ase = ase.build.bulk("Si") * pytuple((4, 1, 1))
+
+# Convert to an AtomsBase-compatible structure
+atoms_ab = pyconvert(AbstractSystem, atoms_ase)
 ```
 
-Similarly, we can write the file with
-[`write_extxyz`](@ref NQCBase.write_extxyz):
-```@repl atoms
-write_extxyz("out.xyz", atoms, positions, cell)
+It is currently not possible to use an AtomsBase system directly with NQCDynamics, but can
+be quickly converted to the correct format:
+
+```julia
+using NQCDynamics
+atoms_nqcd = Atoms(atoms_ab)
+r = Position(atoms_ab)
+v = Velocity(atoms_ab)
+c = Cell(atoms_ab)
 ```
-Both of these functions also work with trajectories such that the positions will be a vector
-of configurations, whilst the atoms and cell will remain unchanged.
 
-If not using `.xyz` files, we can directly use the IO capability of ase to read or the write
-the files.
-This can be done by using the conversions between our data types and the `ase.Atoms` object.
-```@repl atoms
-atoms = Atoms([:H, :H, :C])
-ase_atoms = NQCBase.convert_to_ase_atoms(atoms, rand(3, 3))
-NQCBase.convert_from_ase_atoms(ase_atoms)
+## Saving and loading with AtomsIO.jl
+
+After running a simulation it often desirable to save the trajectory in a standard format for visualization.
+For this, convert the NQCDynamics output into the AtomsBase format,
+then use AtomsIO to write the file in your chosen format.
+
+```julia
+using AtomsIO
+
+system = System(atoms_nqcd, r, v, c)
+
+AtomsIO.save_system("Si.xyz", system) # Save a single image
+
+trajectory = Trajectory(atoms_nqcd, [r, r, r, r], [v, v, v, v], c)
+AtomsIO.save_trajectory("Si.xyz", trajectory) # Save a trajectory
 ```
-These conversions work both ways such that you can read any file format using
-ase then convert the `ase.Atoms` object to our types afterwards.
-Then at the end when you are finished, you can convert them back and write your output
-with ase.
+
+AtomsIO also provides `load_system` and `load_trajectory` which can be converted to the
+NQCDynamics format as above to initialise simulations.
+Refer to [AtomsIO](https://mfherbst.github.io/AtomsIO.jl/stable/) for more information.

docs/src/dynamicssimulations/dynamicsmethods/ehrenfest.md

@@ -40,9 +40,8 @@ final_velocities = [r[:OutputVelocity][end] for r in solution]
 momenta = reduce(vcat, final_velocities*atoms.masses[1])
 ```
 
-Resulting momenta can be plotted by using `StatsPlots` package.
 ```@example ehrenfest
-using StatsPlots
-plot(density(momenta))
+using Plots
+histogram(momenta)
 xlims!(-20,20)
 ```

docs/src/dynamicssimulations/dynamicsmethods/nrpmd.md

@@ -122,7 +122,7 @@ We can check the distribution by plotting the phasespace diagram for each of the
 in our distribution:
 
 ```@example nrpmd
-using CairoMakie
+using Plots
 
 nuclear = distribution.nuclear
 flat_position = reduce(vcat, (nuclear.position[i][:] for i in 1:length(nuclear)))
@@ -175,8 +175,7 @@ the figure from the paper we were attempting to reproduce. Nice!
 ensemble = run_dynamics(sim, (0.0, 30.0), distribution;
     trajectories=100, output, reduction=MeanReduction(), dt=0.1)
 
-plt = lines(0:0.1:30, [p[1,1]-p[2,1] for p in ensemble[:PopulationCorrelationFunction]])
-plt.axis.xlabel = "Time"
-plt.axis.ylabel = "Population difference"
-plt
+plot(0:0.1:30, [p[1,1]-p[2,1] for p in ensemble[:PopulationCorrelationFunction]])
+xlabel!("Time")
+ylabel!("Population difference")
 ```

docs/src/dynamicssimulations/dynamicsmethods/rpmd.md

@@ -85,38 +85,29 @@ This animation shows the cyclic nature of the ring polymer, and how every bead i
 to its two neighbours.
 
 ```@example rpmd
-using CairoMakie
+using Plots
 
 rs = traj[:OutputPosition]
 
-index = Observable(1)
-xs = @lift(rs[$index][1,1,:])
-ys = @lift(rs[$index][2,1,:])
-close_loop_x = @lift([rs[$index][1,1,end], rs[$index][1,1,begin]])
-close_loop_y = @lift([rs[$index][2,1,end], rs[$index][2,1,begin]])
-fig = scatter(xs, ys, axis = (title = @lift("t = $(round(Int, dt*($index-1)))"),))
-lines!(xs, ys)
-lines!(close_loop_x, close_loop_y)
-xlims!(-3, 3)
-ylims!(-3, 3)
-
 timestamps = 1:length(traj[:OutputPosition])
-filepath = "../../assets/figures/rpmd.mp4" # hide
-record(fig, filepath, timestamps;
-        framerate = 30) do i
-    index[] = i
+@gif for i in timestamps
+    xs = rs[i][1,1,:]
+    ys = rs[i][2,1,:]
+    close_loop_x = [rs[i][1,1,end], rs[i][1,1,begin]]
+    close_loop_y = [rs[i][2,1,end], rs[i][2,1,begin]]
+
+    plot(
+        xlims=(-3, 3),
+        ylims=(-3, 3),
+        legend=false
+    )
+
+    plot!(xs, ys, color=:black)
+    scatter!(xs, ys)
+    plot!(close_loop_x, close_loop_y)
 end
-nothing # hide
 ```
 
-![rpmd fig](../../assets/figures/rpmd.mp4)
-
-!!! note
-
-    We have used Makie's animation features to produce this animation. If you want
-    information how Makie works, take a look at the
-    [Makie documentation](https://docs.makie.org/stable/documentation/animation/).
-
 Since this package is focused on nonadiabatic dynamics, you won't see much adiabatic RPMD
-elsewhere in the documentation but it's useful to understand how the original adiabatic
-version works before moving onto the nonadiabatic extensions.
+elsewhere in the documentation, but it's useful to understand how the original adiabatic
+version works before moving on to the nonadiabatic extensions.