From b4ad424fd9ef90f135f6737376dba2cc65d7cece Mon Sep 17 00:00:00 2001
From: "Michael F. Herbst" <info@michael-herbst.com>
Date: Fri, 9 Aug 2024 09:22:43 +0200
Subject: [PATCH 1/7] Not working ideas for printing

---
 src/optimization.jl | 37 +++++++++++++++++++++++++++++++++++++
 1 file changed, 37 insertions(+)

diff --git a/src/optimization.jl b/src/optimization.jl
index 072cf89..6311aff 100644
--- a/src/optimization.jl
+++ b/src/optimization.jl
@@ -104,6 +104,43 @@ function minimize_energy!(system, calculator, solver=Autoselect(); kwargs...)
     _minimize_energy!(system, calculator, solver; kwargs...)
 end
 
+struct GeoOptTabularPrint
+end
+function (cb::GeoOptTabularPrint)(optim_state, geoopt_state)
+    show_force = true  # TODO Set to false if have no forces
+    show_time  = true
+
+    if optim_state.iter == 0
+        label_force = show_force ? ("   max(Force)", "   ----------") : ("", "")
+        label_time  = show_time ? ("   Δtime", "   ------") : ("", "")
+        println("n     Energy            log10(ΔE)", label_force[1], label_time[1])
+        println("---   ---------------   ---------", label_force[2], label_time[2])
+    end
+    return false
+
+    fstr = 
+
+    tstr = " "^9
+    if show_time
+        deltatime_ns = 0
+        tstr = @sprintf "   % 6s" TimerOutputs.prettytime(deltatime_ns)
+    end
+
+
+    # iteration
+    # change in energy (logscale)
+    # current force
+    # current virial
+    # maximal atomic displacement
+    # trace of lattice deformation matrix
+    # timing
+    # other flags (e.g. 
+
+
+
+    return false
+end
+
 # Function that does all the work. The idea is that calculator implementations
 # can provide more specific methods for minimize_energy! calling the function below.
 # This allows to adjust (based on the calculator type) the default parameters,

From 258a0c2f5b529780556cf55be3bed547bc1d6860 Mon Sep 17 00:00:00 2001
From: "Michael F. Herbst" <info@michael-herbst.com>
Date: Fri, 9 Aug 2024 09:37:22 +0200
Subject: [PATCH 2/7] Clarify wording

---
 src/optimization.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/optimization.jl b/src/optimization.jl
index 6311aff..c8cbc14 100644
--- a/src/optimization.jl
+++ b/src/optimization.jl
@@ -92,7 +92,7 @@ can also be employed here.
 - `tol_energy`: Tolerance in the energy to stop the minimisation (all `tol_*` need to be satisfied)
 - `tol_force`:  Tolerance in the force  to stop the minimisation (all `tol_*` need to be satisfied)
 - `tol_virial`: Tolerance in the virial to stop the minimisation (all `tol_*` need to be satisfied)
-- `maxstep`: Maximal step size (in AU) to be taken in a single optimisation step
+- `maxstep`: Maximal step size (in AU or length units) to be taken in a single optimisation step
   (not supported for all `solver`s)
 - `callback`: A custom callback, which obtains the pair `(optimization_state, geoopt_state)` and is
   expected to return `false` (continue iterating) or `true` (halt iterations).

From 88e1215379239b383c5ddff7aaedca688cdf7ee7 Mon Sep 17 00:00:00 2001
From: "Michael F. Herbst" <info@michael-herbst.com>
Date: Fri, 9 Aug 2024 22:12:21 +0200
Subject: [PATCH 3/7] Add printing routines

---
 Project.toml                        |  7 ++-
 docs/src/examples/aluminium_dftk.md | 16 +++--
 docs/src/examples/tial_lj.md        |  2 +-
 docs/src/index.md                   |  4 +-
 src/GeometryOptimization.jl         |  4 ++
 src/clamping_updating_positions.jl  |  1 -
 src/optimization.jl                 | 92 ++++++++++++-----------------
 src/printing.jl                     | 58 ++++++++++++++++++
 8 files changed, 120 insertions(+), 64 deletions(-)
 create mode 100644 src/printing.jl

diff --git a/Project.toml b/Project.toml
index 218e3c5..a9adae4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "GeometryOptimization"
 uuid = "673bf261-a53d-43b9-876f-d3c1fc8329c2"
 authors = ["JuliaMolSim community"]
-version = "0.0.2"
+version = "0.1.0"
 
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
@@ -11,7 +11,10 @@ LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
@@ -23,8 +26,10 @@ DocStringExtensions = "0.9"
 LineSearches = "7"
 Optimization = "3"
 OptimizationOptimJL = "0.3"
+PrettyTables = "2"
 StaticArrays = "1"
 TestItemRunner = "0.2"
+TimerOutputs = "0.5.12"
 Unitful = "1"
 UnitfulAtomic = "1"
 julia = "1.10"
diff --git a/docs/src/examples/aluminium_dftk.md b/docs/src/examples/aluminium_dftk.md
index 14e5ae5..1166716 100644
--- a/docs/src/examples/aluminium_dftk.md
+++ b/docs/src/examples/aluminium_dftk.md
@@ -20,8 +20,8 @@ using DFTK
 
 model_kwargs = (; functionals=[:lda_x, :lda_c_pw], temperature=1e-3)
 basis_kwargs = (; kgrid=(3, 3, 3), Ecut=10.0)
-scf_kwargs   = (; tol=1e-6, mixing=KerkerMixing())
-calc = DFTKCalculator(; model_kwargs, basis_kwargs, scf_kwargs, verbose=true)
+scf_kwargs   = (; mixing=KerkerMixing())
+calc = DFTKCalculator(; model_kwargs, basis_kwargs, scf_kwargs)
 nothing
 ```
 
@@ -42,18 +42,24 @@ nothing
 We perform the structure optimisation using the LBFGS solver
 from Optim with solver parameters adapted for our geometry optimisation setting.
 This is selected by passing the [GeometryOptimization.OptimLBFGS](@ref)
-solver as the third argument.
+solver as the third argument. The `verbosity=2` flag makes sure we get
+output from both the geometry optimisation as well as the inner SCF solver.
 
 ```@example dftk-aluminium
 using GeometryOptimization
 GO = GeometryOptimization
 
 results = minimize_energy!(system, calc, GO.OptimLBFGS();
-                           tol_force=1e-4u"eV/Å",
-                           show_trace=true)
+                           tol_force=1e-4u"eV/Å", verbosity=2)
 nothing
 ```
 
+!!! tip "Automatically adapted calculator parameters"
+    Some calculators (such as DFTK) are able to adapt to the keyword arguments
+    and parameters passed to `minimize_energy!`. In this case the SCF tolerance
+    is automatically adapted according to the convergence parameters
+    (here `tol_force`) passed to `minimize_energy!`.
+
 The final energy is
 ```@example dftk-aluminium
 results.energy
diff --git a/docs/src/examples/tial_lj.md b/docs/src/examples/tial_lj.md
index 9c60414..19525da 100644
--- a/docs/src/examples/tial_lj.md
+++ b/docs/src/examples/tial_lj.md
@@ -26,7 +26,7 @@ Minimise energy:
 ```julia
 ## TODO: Should run as @example once EmpiricalPotentials is compatible
 
-results = minimize_energy!(system, calc, GO.OptimCG(); maxiters=10, show_trace=true)
+results = minimize_energy!(system, calc, GO.OptimCG(); maxiters=10, verbosity=1)
 results.energy
 ```
 
diff --git a/docs/src/index.md b/docs/src/index.md
index 1b567dd..25fb893 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -27,8 +27,8 @@ system = isolated_system([:H => [0, 0, 0.0]u"bohr",
                           :H => [0, 0, 1.9]u"bohr"])
 calc = LennardJones(-1.17u"hartree", 0.743u"angstrom", 1, 1, 0.6u"nm")
 
-# Run the geometry optimisation
-results = minimize_energy!(system, calc)
+# Run the geometry optimisation (using verbosity=1 to print the progress)
+results = minimize_energy!(system, calc; verbosity=1)
 
 # Inspect the results
 optimised_system = results.system
diff --git a/src/GeometryOptimization.jl b/src/GeometryOptimization.jl
index 47a5188..1426486 100644
--- a/src/GeometryOptimization.jl
+++ b/src/GeometryOptimization.jl
@@ -24,5 +24,9 @@ $(DOCSTRING)
 
 include("clamping_updating_positions.jl")
 include("optimization.jl")
+include("printing.jl")
+
+export minimize_energy!
+export GeoOptDefaultPrint
 
 end
diff --git a/src/clamping_updating_positions.jl b/src/clamping_updating_positions.jl
index c8cf187..2f5d587 100644
--- a/src/clamping_updating_positions.jl
+++ b/src/clamping_updating_positions.jl
@@ -2,7 +2,6 @@
 # Interface between AtomsBase.jl and GeometryOptimization.jl that provides 
 # utility functions for manipulating systems.
 #
-export clamp_atoms
 
 """
 Convert the input system to a kind of system we can work with, i.e. one where
diff --git a/src/optimization.jl b/src/optimization.jl
index c8cbc14..18c7100 100644
--- a/src/optimization.jl
+++ b/src/optimization.jl
@@ -3,21 +3,25 @@
 # We always work in cartesian coordinates.
 #
 
-export minimize_energy!
-
 using AtomsCalculators: Energy, Forces, calculate
 
 mutable struct GeometryOptimizationState
-    calc_state  # Reference to the current calculator state
-    energy      # Most recent energy value
-    forces      # Most recent force value
-    virial      # Most recent virial value
+    calc_state          # Reference to the most recent calculator state
+    energy              # Current energy value
+    forces              # Current force value
+    virial              # Current virial value
+    energy_change       # Change in energy since previous step
+    start_time::UInt64  # Time when the calculation started from time_ns()
 end
-function GeometryOptimizationState(system, calculator)
-    GeometryOptimizationState(AC.get_state(calculator),
-                              AC.zero_energy(system, calculator),
-                              AC.zero_forces(system, calculator),
-                              AC.zero_virial(system, calculator))
+function GeometryOptimizationState(system, calculator; start_time=time_ns())
+    GeometryOptimizationState(
+        AC.get_state(calculator),
+        AC.zero_energy(system, calculator),
+        AC.zero_forces(system, calculator),
+        AC.zero_virial(system, calculator),
+        AC.zero_energy(system, calculator),
+        start_time,
+    )
 end
 
 """
@@ -48,7 +52,6 @@ function Optimization.OptimizationProblem(system, calculator, geoopt_state; kwar
         new_system = update_not_clamped_positions(system, x * u"bohr")
         res = calculate(Energy(), new_system, calculator, ps, geoopt_state.calc_state)
         geoopt_state.calc_state = res.state
-        geoopt_state.energy = res.energy
         austrip(res.energy), geoopt_state
     end
     g! = function(G::AbstractVector{<:Real}, x::AbstractVector{<:Real}, ps)
@@ -89,11 +92,14 @@ can also be employed here.
 
 ## Keyword arguments:
 - `maxiters`: Maximal number of iterations
+- `maxtime`:  Maximal allowed runtime (in seconds)
 - `tol_energy`: Tolerance in the energy to stop the minimisation (all `tol_*` need to be satisfied)
 - `tol_force`:  Tolerance in the force  to stop the minimisation (all `tol_*` need to be satisfied)
 - `tol_virial`: Tolerance in the virial to stop the minimisation (all `tol_*` need to be satisfied)
 - `maxstep`: Maximal step size (in AU or length units) to be taken in a single optimisation step
   (not supported for all `solver`s)
+- `verbosity`: Printing level. The idea is that `0` is silent, `1` displays the optimisation
+  progress and `≥ 2` starts displaying things from the calculator as well (e.g SCF iterations).
 - `callback`: A custom callback, which obtains the pair `(optimization_state, geoopt_state)` and is
   expected to return `false` (continue iterating) or `true` (halt iterations).
 - `kwargs`: All other keyword arguments are passed to the call to `solve`. Note, that
@@ -104,43 +110,6 @@ function minimize_energy!(system, calculator, solver=Autoselect(); kwargs...)
     _minimize_energy!(system, calculator, solver; kwargs...)
 end
 
-struct GeoOptTabularPrint
-end
-function (cb::GeoOptTabularPrint)(optim_state, geoopt_state)
-    show_force = true  # TODO Set to false if have no forces
-    show_time  = true
-
-    if optim_state.iter == 0
-        label_force = show_force ? ("   max(Force)", "   ----------") : ("", "")
-        label_time  = show_time ? ("   Δtime", "   ------") : ("", "")
-        println("n     Energy            log10(ΔE)", label_force[1], label_time[1])
-        println("---   ---------------   ---------", label_force[2], label_time[2])
-    end
-    return false
-
-    fstr = 
-
-    tstr = " "^9
-    if show_time
-        deltatime_ns = 0
-        tstr = @sprintf "   % 6s" TimerOutputs.prettytime(deltatime_ns)
-    end
-
-
-    # iteration
-    # change in energy (logscale)
-    # current force
-    # current virial
-    # maximal atomic displacement
-    # trace of lattice deformation matrix
-    # timing
-    # other flags (e.g. 
-
-
-
-    return false
-end
-
 # Function that does all the work. The idea is that calculator implementations
 # can provide more specific methods for minimize_energy! calling the function below.
 # This allows to adjust (based on the calculator type) the default parameters,
@@ -148,11 +117,13 @@ end
 # by calling this function the actual minimisation is started off.
 function _minimize_energy!(system, calculator, solver;
                            maxiters=100,
+                           maxtime=60*60*24*365,   # 1 year
                            tol_energy=Inf*u"eV",
                            tol_force=1e-4u"eV/Å",  # VASP default
                            tol_virial=1e-6u"eV",   # TODO How reasonable ?
                            maxstep=0.8u"bohr",
-                           callback=(x,y) -> false,
+                           verbosity=0,
+                           callback=default_printing_callback(verbosity),
                            kwargs...)
     solver = setup_solver(system, calculator, solver; maxstep)
     system = convert_to_updatable(system)
@@ -161,26 +132,39 @@ function _minimize_energy!(system, calculator, solver;
     problem = OptimizationProblem(system, calculator, geoopt_state; sense=Optimization.MinSense)
     converged = false
 
-    Eold = AC.zero_energy(system, calculator)
+    Eold = 0  # Atomic units
     function inner_callback(optim_state, ::Any, geoopt_state)
+        geoopt_state.energy        = optim_state.objective * u"hartree"
+        geoopt_state.energy_change = (optim_state.objective - Eold) * u"hartree"
+
         halt = callback(optim_state, geoopt_state)
         halt && return true
 
-        energy_converged = austrip(abs(geoopt_state.energy - Eold))    < austrip(tol_energy)
+        energy_converged = optim_state.objective - Eold < austrip(tol_energy)
         force_converged  = austrip(maximum(norm, geoopt_state.forces)) < austrip(tol_force)
         virial_converged = austrip(maximum(abs,  geoopt_state.virial)) < austrip(tol_virial)
 
-        Eold = geoopt_state.energy
+        Eold = optim_state.objective
         converged = energy_converged && force_converged && virial_converged
         return converged
     end
 
-    optimres = solve(problem, solver; maxiters, callback=inner_callback, kwargs...)
+    optimres = solve(problem, solver; maxiters, maxtime, callback=inner_callback, kwargs...)
     (; system=update_not_clamped_positions(system, optimres.u * u"bohr"), converged,
        energy=optimres.objective, geoopt_state.forces, geoopt_state.virial,
        state=geoopt_state.calc_state, optimres.stats, optimres.alg, optimres)
 end
 
+function default_printing_callback(verbosity::Integer)
+    if verbosity <= 0
+        return (os, gs) -> false
+    elseif verbosity == 1
+        return GeoOptDefaultPrint(; always_show_header=false)
+    elseif verbosity ≥  2
+        return GeoOptDefaultPrint(; always_show_header=true)
+    end
+end
+
 # Default setup_solver function just passes things through
 setup_solver(system, calculator, solver::Any; kwargs...) = solver
 
diff --git a/src/printing.jl b/src/printing.jl
new file mode 100644
index 0000000..0649e32
--- /dev/null
+++ b/src/printing.jl
@@ -0,0 +1,58 @@
+import TimerOutputs
+using PrettyTables
+using Printf
+
+"""
+Callback producing a convergence table summarising the geometry
+optimisation convergence. If `always_show_header=true` the
+header is shown in each iteration. This is helpful if the calculator
+produces output as well.
+"""
+struct GeoOptDefaultPrint{F}
+    always_show_header::Bool
+    callback_extra_fields::F
+    prev_time::Ref{UInt64}
+end
+function GeoOptDefaultPrint(; always_show_header=false, callback_extra_fields=nothing)
+    GeoOptDefaultPrint(always_show_header, callback_extra_fields, Ref{UInt64}(0))
+end
+
+format_log8(value) = (value < 0 ? " " : "+") * (@sprintf "%8.2f" log10(abs(value)))
+
+function (cb::GeoOptDefaultPrint)(optim_state, geoopt_state)
+    # If first iteration clear a potentially cached previous time
+    optim_state.iter ≤ 0 && (cb.prev_time[] = 0)
+
+    runtime_ns = time_ns() - geoopt_state.start_time
+    tstr = @sprintf "% 6s" TimerOutputs.prettytime(runtime_ns - cb.prev_time[])
+    cb.prev_time[] = runtime_ns
+
+    Estr  = (@sprintf "%+15.12f" round(austrip(geoopt_state.energy), sigdigits=13))[1:15]
+    logΔE = optim_state.iter < 1 ? "" : format_log8(austrip(geoopt_state.energy_change))
+
+    maxforce = austrip(maximum(norm, geoopt_state.forces))
+    fstr = iszero(maxforce) ? "" : round(maxforce, sigdigits=8)
+
+    data = [  # Header, width, value
+        ("n",           3, optim_state.iter),
+        ("Energy",     15, Estr),
+        ("log10(ΔE)",   9, logΔE),
+        ("max(Force)", 10, fstr),
+        # TODO Maximal atomic displacement
+        # TODO Current virial, trace of lattice deformation matrix
+        ("Δtime",       6, tstr),
+    ]
+    if !isnothing(cb.callback_extra_fields)
+        append!(data, cb.additional_fields(optim_state, geoopt_state))
+    end
+
+    title = iszero(optim_state.iter) ? "Geometry optimisation convergence (in atomic units)" : ""
+    show_header = optim_state.iter == 0 || cb.always_show_header
+    pretty_table(reshape(getindex.(data, 3), 1, length(data));
+                 header=Symbol.(getindex.(data, 1)), title,
+                 columns_width=getindex.(data, 2),
+                 show_header, hlines=show_header ? [0, 1] : :none)
+
+    flush(stdout)
+    return false
+end

From a6d36880e9176d828295e67a8d0065b911a02554 Mon Sep 17 00:00:00 2001
From: "Michael F. Herbst" <info@michael-herbst.com>
Date: Sat, 10 Aug 2024 14:28:06 +0200
Subject: [PATCH 4/7] Printing and tests

---
 Project.toml                        |  4 +++-
 docs/src/examples/aluminium_dftk.md |  4 ++--
 docs/src/examples/other_solvers.md  |  2 +-
 src/clamping_updating_positions.jl  |  1 +
 src/optimization.jl                 | 26 ++++++++++-----------
 src/printing.jl                     | 18 +++++++++++----
 test/interface_test.jl              | 20 +++++++++++-----
 test/minimization.jl                | 36 +++++++++++++++++++++++++++++
 test/testcalculators.jl             | 19 ---------------
 9 files changed, 83 insertions(+), 47 deletions(-)
 create mode 100644 test/minimization.jl
 delete mode 100644 test/testcalculators.jl

diff --git a/Project.toml b/Project.toml
index a9adae4..6294249 100644
--- a/Project.toml
+++ b/Project.toml
@@ -35,9 +35,11 @@ UnitfulAtomic = "1"
 julia = "1.10"
 
 [extras]
+ASEconvert = "3da9722f-58c2-4165-81be-b4d7253e8fd2"
 AtomsBuilder = "f5cc8831-eeb7-4288-8d9f-d6c1ddb77004"
+PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["AtomsBuilder", "Test", "TestItemRunner"]
+test = ["ASEconvert", "AtomsBuilder", "PythonCall", "Test", "TestItemRunner"]
diff --git a/docs/src/examples/aluminium_dftk.md b/docs/src/examples/aluminium_dftk.md
index 1166716..abe2985 100644
--- a/docs/src/examples/aluminium_dftk.md
+++ b/docs/src/examples/aluminium_dftk.md
@@ -50,7 +50,7 @@ using GeometryOptimization
 GO = GeometryOptimization
 
 results = minimize_energy!(system, calc, GO.OptimLBFGS();
-                           tol_force=1e-4u"eV/Å", verbosity=2)
+                           tol_forces=1e-4u"eV/Å", verbosity=2)
 nothing
 ```
 
@@ -58,7 +58,7 @@ nothing
     Some calculators (such as DFTK) are able to adapt to the keyword arguments
     and parameters passed to `minimize_energy!`. In this case the SCF tolerance
     is automatically adapted according to the convergence parameters
-    (here `tol_force`) passed to `minimize_energy!`.
+    (here `tol_forces`) passed to `minimize_energy!`.
 
 The final energy is
 ```@example dftk-aluminium
diff --git a/docs/src/examples/other_solvers.md b/docs/src/examples/other_solvers.md
index 300199f..eadfd5f 100644
--- a/docs/src/examples/other_solvers.md
+++ b/docs/src/examples/other_solvers.md
@@ -46,7 +46,7 @@ using OptimizationNLopt
 solver = NLopt.LD_TNEWTON()
 
 results = minimize_energy!(system, calc, solver;
-                           tol_force=1e-4u"eV/Å", maxeval=100)
+                           tol_forces=1e-4u"eV/Å", maxeval=100)
 nothing
 ```
 
diff --git a/src/clamping_updating_positions.jl b/src/clamping_updating_positions.jl
index 2f5d587..77262fa 100644
--- a/src/clamping_updating_positions.jl
+++ b/src/clamping_updating_positions.jl
@@ -47,6 +47,7 @@ indices corresponding to their positions in the system.
     in the future.
 """
 function clamp_atoms(system, clamped_indexes::Union{AbstractVector{<:Integer},Nothing})
+    system = convert_to_updatable(system)
     clamped = falses(length(system))
     clamped[clamped_indexes] .= true
     particles = [Atom(atom; clamped=msk) for (atom, msk) in zip(system, clamped)]
diff --git a/src/optimization.jl b/src/optimization.jl
index 18c7100..cb9a205 100644
--- a/src/optimization.jl
+++ b/src/optimization.jl
@@ -94,14 +94,16 @@ can also be employed here.
 - `maxiters`: Maximal number of iterations
 - `maxtime`:  Maximal allowed runtime (in seconds)
 - `tol_energy`: Tolerance in the energy to stop the minimisation (all `tol_*` need to be satisfied)
-- `tol_force`:  Tolerance in the force  to stop the minimisation (all `tol_*` need to be satisfied)
+- `tol_forces`:  Tolerance in the force  to stop the minimisation (all `tol_*` need to be satisfied)
 - `tol_virial`: Tolerance in the virial to stop the minimisation (all `tol_*` need to be satisfied)
 - `maxstep`: Maximal step size (in AU or length units) to be taken in a single optimisation step
   (not supported for all `solver`s)
 - `verbosity`: Printing level. The idea is that `0` is silent, `1` displays the optimisation
   progress and `≥ 2` starts displaying things from the calculator as well (e.g SCF iterations).
 - `callback`: A custom callback, which obtains the pair `(optimization_state, geoopt_state)` and is
-  expected to return `false` (continue iterating) or `true` (halt iterations).
+  expected to return `false` (continue iterating) or `true` (halt iterations). Note that
+  specifying this overwrites the default printing callback. The calculation thus becomes
+  silent unless a [`GeoOptDefaultPrint`](@ref) is included in the callback.
 - `kwargs`: All other keyword arguments are passed to the call to `solve`. Note, that
   if special `kwargs` should be passed to the `Optimization.OptimizationProblem` the user
   needs to setup the problem manually (e.g. `OptimizationProblem(system, calculator)`)
@@ -116,13 +118,13 @@ end
 # do some additional calculator-specific setup (e.g. callbacks) and so on. Then
 # by calling this function the actual minimisation is started off.
 function _minimize_energy!(system, calculator, solver;
-                           maxiters=100,
-                           maxtime=60*60*24*365,   # 1 year
+                           maxiters::Integer=100,
+                           maxtime::Integer=60*60*24*365,  # 1 year
                            tol_energy=Inf*u"eV",
-                           tol_force=1e-4u"eV/Å",  # VASP default
-                           tol_virial=1e-6u"eV",   # TODO How reasonable ?
+                           tol_forces=1e-4u"eV/Å",  # VASP default
+                           tol_virial=1e-6u"eV",    # TODO How reasonable ?
                            maxstep=0.8u"bohr",
-                           verbosity=0,
+                           verbosity::Integer=0,
                            callback=default_printing_callback(verbosity),
                            kwargs...)
     solver = setup_solver(system, calculator, solver; maxstep)
@@ -141,7 +143,7 @@ function _minimize_energy!(system, calculator, solver;
         halt && return true
 
         energy_converged = optim_state.objective - Eold < austrip(tol_energy)
-        force_converged  = austrip(maximum(norm, geoopt_state.forces)) < austrip(tol_force)
+        force_converged  = austrip(maximum(norm, geoopt_state.forces)) < austrip(tol_forces)
         virial_converged = austrip(maximum(abs,  geoopt_state.virial)) < austrip(tol_virial)
 
         Eold = optim_state.objective
@@ -151,17 +153,15 @@ function _minimize_energy!(system, calculator, solver;
 
     optimres = solve(problem, solver; maxiters, maxtime, callback=inner_callback, kwargs...)
     (; system=update_not_clamped_positions(system, optimres.u * u"bohr"), converged,
-       energy=optimres.objective, geoopt_state.forces, geoopt_state.virial,
+       energy=optimres.objective * u"hartree", geoopt_state.forces, geoopt_state.virial,
        state=geoopt_state.calc_state, optimres.stats, optimres.alg, optimres)
 end
 
 function default_printing_callback(verbosity::Integer)
     if verbosity <= 0
         return (os, gs) -> false
-    elseif verbosity == 1
-        return GeoOptDefaultPrint(; always_show_header=false)
-    elseif verbosity ≥  2
-        return GeoOptDefaultPrint(; always_show_header=true)
+    else
+        return GeoOptDefaultPrint(; always_show_header=verbosity > 1)
     end
 end
 
diff --git a/src/printing.jl b/src/printing.jl
index 0649e32..a61b56e 100644
--- a/src/printing.jl
+++ b/src/printing.jl
@@ -43,15 +43,23 @@ function (cb::GeoOptDefaultPrint)(optim_state, geoopt_state)
         ("Δtime",       6, tstr),
     ]
     if !isnothing(cb.callback_extra_fields)
-        append!(data, cb.additional_fields(optim_state, geoopt_state))
+        append!(data, cb.callback_extra_fields(geoopt_state.calc_state))
     end
 
+    if cb.always_show_header
+        hlines = :all
+        show_header = true
+    elseif optim_state.iter == 0
+        hlines = [0, 1]
+        show_header = true
+    else
+        hlines = :none
+        show_header = false
+    end
     title = iszero(optim_state.iter) ? "Geometry optimisation convergence (in atomic units)" : ""
-    show_header = optim_state.iter == 0 || cb.always_show_header
     pretty_table(reshape(getindex.(data, 3), 1, length(data));
-                 header=Symbol.(getindex.(data, 1)), title,
-                 columns_width=getindex.(data, 2),
-                 show_header, hlines=show_header ? [0, 1] : :none)
+                 header=Symbol.(getindex.(data, 1)), columns_width=getindex.(data, 2),
+                 title, show_header, hlines)
 
     flush(stdout)
     return false
diff --git a/test/interface_test.jl b/test/interface_test.jl
index dde767b..173c26f 100644
--- a/test/interface_test.jl
+++ b/test/interface_test.jl
@@ -1,21 +1,29 @@
-@testitem "Test GeometryOptimization with AtomsCalculators interface" #=
-        =# setup=[TestCalculators] begin
+@testitem "Test GeometryOptimization with AtomsCalculators interface" begin
     # This is a fake test. It does not really test anything,
     # just that the interface code compiles and works properly.
 
     using AtomsBuilder
+    using AtomsCalculators
     using GeometryOptimization
     using LinearAlgebra
     using Unitful
     using UnitfulAtomic
+    AC = AtomsCalculators
     GO = GeometryOptimization
 
-    system = rattle!(bulk(:Al; cubic=true), 0.1u"Å")
-    system = clamp_atoms(system, [1])
+    struct DummyCalc end
+    AC.energy_unit(::DummyCalc) = u"hartree"
+    AC.length_unit(::DummyCalc) = u"bohr"
+    AC.@generate_interface function AC.potential_energy(system, calc::DummyCalc; kwargs...)
+        AC.zero_energy(system, calc)
+    end
+    AC.@generate_interface function AC.forces(system, calc::DummyCalc; kwargs...)
+        AC.zero_forces(system, calc)
+    end
 
-    calculator = TestCalculators.DummyCalc()
+    system = GO.clamp_atoms(bulk(:Al; cubic=true), [1])
     for solver in (GO.Autoselect(), GO.OptimCG(), GO.OptimLBFGS(), GO.OptimSD())
-        results = minimize_energy!(system, calculator, solver)
+        results = minimize_energy!(system, DummyCalc(), solver; verbosity=1)
         @test norm(position(results.system[1]) - position(system[1])) < 1e-5u"bohr"
     end
 end
diff --git a/test/minimization.jl b/test/minimization.jl
new file mode 100644
index 0000000..63f1dd0
--- /dev/null
+++ b/test/minimization.jl
@@ -0,0 +1,36 @@
+@testitem "Test aluminium EMT minimisation" begin
+    using ASEconvert
+    using AtomsBuilder
+    using AtomsCalculators
+    using GeometryOptimization
+    using LinearAlgebra
+    using PythonCall
+    using Unitful
+    using UnitfulAtomic
+    AC = AtomsCalculators
+    GO = GeometryOptimization
+
+    ase_emt = pyimport("ase.calculators.emt")
+    calculator = ASEcalculator(ase_emt.EMT())
+
+    # Get a minimised reference
+    aluminium_final = minimize_energy!(bulk(:Al; cubic=true), calculator;
+                                       tol_forces=1e-10).system
+    aluminium_init  = rattle!(bulk(:Al; cubic=true), 0.2u"Å")
+    energy_final = AC.potential_energy(aluminium_final, calculator)
+    energy_init  = AC.potential_energy(aluminium_init,  calculator)
+    @assert energy_final ≤ energy_init
+
+    for solver in (GO.Autoselect(), GO.OptimCG(), GO.OptimLBFGS(), GO.OptimSD())
+        results = minimize_energy!(aluminium_init, calculator, solver;
+                                   tol_forces=1e-9, maxiters=500)
+        @test results.energy ≤ energy_init
+        @test austrip(abs(results.energy - energy_final)) < 1e-12
+
+        delta = position(results.system, 1) - position(aluminium_final, 1)
+        for i in 2:length(aluminium_final)
+            @test norm(  position(results.system, i) - delta
+                       - position(aluminium_final, i)) < 1e-6u"bohr"
+        end
+    end
+end
diff --git a/test/testcalculators.jl b/test/testcalculators.jl
deleted file mode 100644
index adbf300..0000000
--- a/test/testcalculators.jl
+++ /dev/null
@@ -1,19 +0,0 @@
-# Create a dummy AtomsCalculator to test the geometry optimization interfcae.
-@testsetup module TestCalculators
-    using AtomsCalculators
-    using Unitful
-    using UnitfulAtomic
-    AC = AtomsCalculators
-
-    struct DummyCalc end
-    AC.energy_unit(::DummyCalc) = u"hartree"
-    AC.length_unit(::DummyCalc) = u"bohr"
-
-    AC.@generate_interface function AC.potential_energy(system, calc::DummyCalc; kwargs...)
-        AC.zero_energy(system, calc)
-    end
-
-    AC.@generate_interface function AC.forces(system, calc::DummyCalc; kwargs...)
-        AC.zero_forces(system, calc)
-    end
-end

From b321e2144b0e75b20f32ae953133b36841186caf Mon Sep 17 00:00:00 2001
From: "Michael F. Herbst" <info@michael-herbst.com>
Date: Sat, 10 Aug 2024 16:37:51 +0200
Subject: [PATCH 5/7] Simplify

---
 src/GeometryOptimization.jl |  2 +-
 src/optimization.jl         | 13 ++++---------
 src/printing.jl             | 28 ++++++++++++++++------------
 3 files changed, 21 insertions(+), 22 deletions(-)

diff --git a/src/GeometryOptimization.jl b/src/GeometryOptimization.jl
index 1426486..28750d3 100644
--- a/src/GeometryOptimization.jl
+++ b/src/GeometryOptimization.jl
@@ -27,6 +27,6 @@ include("optimization.jl")
 include("printing.jl")
 
 export minimize_energy!
-export GeoOptDefaultPrint
+export GeoOptDefaultCallback
 
 end
diff --git a/src/optimization.jl b/src/optimization.jl
index cb9a205..dfcc7ed 100644
--- a/src/optimization.jl
+++ b/src/optimization.jl
@@ -6,6 +6,7 @@
 using AtomsCalculators: Energy, Forces, calculate
 
 mutable struct GeometryOptimizationState
+    calculator
     calc_state          # Reference to the most recent calculator state
     energy              # Current energy value
     forces              # Current force value
@@ -15,6 +16,7 @@ mutable struct GeometryOptimizationState
 end
 function GeometryOptimizationState(system, calculator; start_time=time_ns())
     GeometryOptimizationState(
+        calculator,
         AC.get_state(calculator),
         AC.zero_energy(system, calculator),
         AC.zero_forces(system, calculator),
@@ -125,7 +127,7 @@ function _minimize_energy!(system, calculator, solver;
                            tol_virial=1e-6u"eV",    # TODO How reasonable ?
                            maxstep=0.8u"bohr",
                            verbosity::Integer=0,
-                           callback=default_printing_callback(verbosity),
+                           callback=GeoOptDefaultCallback(verbosity),
                            kwargs...)
     solver = setup_solver(system, calculator, solver; maxstep)
     system = convert_to_updatable(system)
@@ -152,19 +154,12 @@ function _minimize_energy!(system, calculator, solver;
     end
 
     optimres = solve(problem, solver; maxiters, maxtime, callback=inner_callback, kwargs...)
+    converged || @warn "Geometry optimisation not converged."
     (; system=update_not_clamped_positions(system, optimres.u * u"bohr"), converged,
        energy=optimres.objective * u"hartree", geoopt_state.forces, geoopt_state.virial,
        state=geoopt_state.calc_state, optimres.stats, optimres.alg, optimres)
 end
 
-function default_printing_callback(verbosity::Integer)
-    if verbosity <= 0
-        return (os, gs) -> false
-    else
-        return GeoOptDefaultPrint(; always_show_header=verbosity > 1)
-    end
-end
-
 # Default setup_solver function just passes things through
 setup_solver(system, calculator, solver::Any; kwargs...) = solver
 
diff --git a/src/printing.jl b/src/printing.jl
index a61b56e..e9e96f1 100644
--- a/src/printing.jl
+++ b/src/printing.jl
@@ -8,21 +8,26 @@ optimisation convergence. If `always_show_header=true` the
 header is shown in each iteration. This is helpful if the calculator
 produces output as well.
 """
-struct GeoOptDefaultPrint{F}
+struct GeoOptDefaultCallback
     always_show_header::Bool
-    callback_extra_fields::F
     prev_time::Ref{UInt64}
 end
-function GeoOptDefaultPrint(; always_show_header=false, callback_extra_fields=nothing)
-    GeoOptDefaultPrint(always_show_header, callback_extra_fields, Ref{UInt64}(0))
+function GeoOptDefaultCallback(; always_show_header=false)
+    GeoOptDefaultCallback(always_show_header, Ref{UInt64}(0))
+end
+function GeoOptDefaultCallback(verbosity::Integer=1; kwargs...)
+    if verbosity <= 0
+        return (os, gs) -> false  # No printing => no callback
+    else
+        return GeoOptDefaultCallback(; always_show_header=verbosity > 1; kwargs...)
+    end
 end
 
 format_log8(value) = (value < 0 ? " " : "+") * (@sprintf "%8.2f" log10(abs(value)))
 
-function (cb::GeoOptDefaultPrint)(optim_state, geoopt_state)
+function (cb::GeoOptDefaultCallback)(optim_state, geoopt_state)
     # If first iteration clear a potentially cached previous time
     optim_state.iter ≤ 0 && (cb.prev_time[] = 0)
-
     runtime_ns = time_ns() - geoopt_state.start_time
     tstr = @sprintf "% 6s" TimerOutputs.prettytime(runtime_ns - cb.prev_time[])
     cb.prev_time[] = runtime_ns
@@ -33,7 +38,7 @@ function (cb::GeoOptDefaultPrint)(optim_state, geoopt_state)
     maxforce = austrip(maximum(norm, geoopt_state.forces))
     fstr = iszero(maxforce) ? "" : round(maxforce, sigdigits=8)
 
-    data = [  # Header, width, value
+    fields = [  # Header, width, value
         ("n",           3, optim_state.iter),
         ("Energy",     15, Estr),
         ("log10(ΔE)",   9, logΔE),
@@ -41,10 +46,9 @@ function (cb::GeoOptDefaultPrint)(optim_state, geoopt_state)
         # TODO Maximal atomic displacement
         # TODO Current virial, trace of lattice deformation matrix
         ("Δtime",       6, tstr),
+        # TODO Would be nice to have some simple way to add in
+        #      a few calculator-specific things (e.g. total number of SCF iterations)
     ]
-    if !isnothing(cb.callback_extra_fields)
-        append!(data, cb.callback_extra_fields(geoopt_state.calc_state))
-    end
 
     if cb.always_show_header
         hlines = :all
@@ -57,8 +61,8 @@ function (cb::GeoOptDefaultPrint)(optim_state, geoopt_state)
         show_header = false
     end
     title = iszero(optim_state.iter) ? "Geometry optimisation convergence (in atomic units)" : ""
-    pretty_table(reshape(getindex.(data, 3), 1, length(data));
-                 header=Symbol.(getindex.(data, 1)), columns_width=getindex.(data, 2),
+    pretty_table(reshape(getindex.(fields, 3), 1, length(fields));
+                 header=Symbol.(getindex.(fields, 1)), columns_width=getindex.(fields, 2),
                  title, show_header, hlines)
 
     flush(stdout)

From 05c898b1c50dae792b92cb08d9a45a998495ec65 Mon Sep 17 00:00:00 2001
From: "Michael F. Herbst" <info@michael-herbst.com>
Date: Sat, 10 Aug 2024 18:36:04 +0200
Subject: [PATCH 6/7] up

---
 src/GeometryOptimization.jl       |  2 +-
 src/{printing.jl => callbacks.jl} | 19 +++++++++----------
 2 files changed, 10 insertions(+), 11 deletions(-)
 rename src/{printing.jl => callbacks.jl} (82%)

diff --git a/src/GeometryOptimization.jl b/src/GeometryOptimization.jl
index 28750d3..c776f91 100644
--- a/src/GeometryOptimization.jl
+++ b/src/GeometryOptimization.jl
@@ -24,7 +24,7 @@ $(DOCSTRING)
 
 include("clamping_updating_positions.jl")
 include("optimization.jl")
-include("printing.jl")
+include("callbacks.jl")
 
 export minimize_energy!
 export GeoOptDefaultCallback
diff --git a/src/printing.jl b/src/callbacks.jl
similarity index 82%
rename from src/printing.jl
rename to src/callbacks.jl
index e9e96f1..fa385ce 100644
--- a/src/printing.jl
+++ b/src/callbacks.jl
@@ -9,23 +9,19 @@ header is shown in each iteration. This is helpful if the calculator
 produces output as well.
 """
 struct GeoOptDefaultCallback
+    verbosity::Int
     always_show_header::Bool
     prev_time::Ref{UInt64}
 end
-function GeoOptDefaultCallback(; always_show_header=false)
-    GeoOptDefaultCallback(always_show_header, Ref{UInt64}(0))
-end
-function GeoOptDefaultCallback(verbosity::Integer=1; kwargs...)
-    if verbosity <= 0
-        return (os, gs) -> false  # No printing => no callback
-    else
-        return GeoOptDefaultCallback(; always_show_header=verbosity > 1; kwargs...)
-    end
+function GeoOptDefaultCallback(verbosity=1; always_show_header=verbosity > 1)
+    GeoOptDefaultCallback(verbosity, always_show_header, Ref{UInt64}(0))
 end
 
 format_log8(value) = (value < 0 ? " " : "+") * (@sprintf "%8.2f" log10(abs(value)))
 
 function (cb::GeoOptDefaultCallback)(optim_state, geoopt_state)
+    cb.verbosity ≤ 0 && return false  # No printing, just continue iterations
+
     # If first iteration clear a potentially cached previous time
     optim_state.iter ≤ 0 && (cb.prev_time[] = 0)
     runtime_ns = time_ns() - geoopt_state.start_time
@@ -61,9 +57,12 @@ function (cb::GeoOptDefaultCallback)(optim_state, geoopt_state)
         show_header = false
     end
     title = iszero(optim_state.iter) ? "Geometry optimisation convergence (in atomic units)" : ""
-    pretty_table(reshape(getindex.(fields, 3), 1, length(fields));
+
+    cb.always_show_header && println(stdout)
+    pretty_table(stdout, reshape(getindex.(fields, 3), 1, length(fields));
                  header=Symbol.(getindex.(fields, 1)), columns_width=getindex.(fields, 2),
                  title, show_header, hlines)
+    cb.always_show_header && println(stdout)
 
     flush(stdout)
     return false

From 05b43040c738ce7f3e8f2f947e6efd3a47bf46ca Mon Sep 17 00:00:00 2001
From: "Michael F. Herbst" <info@michael-herbst.com>
Date: Sat, 10 Aug 2024 19:20:51 +0200
Subject: [PATCH 7/7] up

---
 src/optimization.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/optimization.jl b/src/optimization.jl
index dfcc7ed..8ea0c83 100644
--- a/src/optimization.jl
+++ b/src/optimization.jl
@@ -102,10 +102,10 @@ can also be employed here.
   (not supported for all `solver`s)
 - `verbosity`: Printing level. The idea is that `0` is silent, `1` displays the optimisation
   progress and `≥ 2` starts displaying things from the calculator as well (e.g SCF iterations).
-- `callback`: A custom callback, which obtains the pair `(optimization_state, geoopt_state)` and is
-  expected to return `false` (continue iterating) or `true` (halt iterations). Note that
-  specifying this overwrites the default printing callback. The calculation thus becomes
-  silent unless a [`GeoOptDefaultPrint`](@ref) is included in the callback.
+- `callback`: A custom callback, which obtains the pair `(optimization_state, geoopt_state)`
+  and is expected to return `false` (continue iterating) or `true` (halt iterations). Note
+  that specifying this overwrites the default printing callback. The calculation thus becomes
+  silent unless a [`GeoOptDefaultCallback`](@ref) is included in the callback.
 - `kwargs`: All other keyword arguments are passed to the call to `solve`. Note, that
   if special `kwargs` should be passed to the `Optimization.OptimizationProblem` the user
   needs to setup the problem manually (e.g. `OptimizationProblem(system, calculator)`)