Integrate with DofManager from GeomOpt.jl (#24)

Author: mfherbst

Project.toml

@@ -19,7 +19,6 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [compat]
-ASEconvert = "0.1.8"
 AtomsBase = "0.5"
 AtomsBuilder = "0.2.2"
 AtomsCalculators = "0.2.3"
@@ -36,12 +35,10 @@ UnitfulAtomic = "1"
 julia = "1.10"
 
 [extras]
-ASEconvert = "3da9722f-58c2-4165-81be-b4d7253e8fd2"
 AtomsBuilder = "f5cc8831-eeb7-4288-8d9f-d6c1ddb77004"
 EmpiricalPotentials = "38527215-9240-4c91-a638-d4250620c9e2"
-PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["ASEconvert", "AtomsBuilder", "EmpiricalPotentials", "PythonCall", "Test", "TestItemRunner"]
+test = ["AtomsBuilder", "EmpiricalPotentials", "Test", "TestItemRunner"]

docs/make.jl

@@ -26,6 +26,7 @@ makedocs(;
         "Home" => "index.md",
         "Examples" => [
             "examples/aluminium_dftk.md",
+            "examples/variablecell.md",
             "examples/other_solvers.md",
             "examples/tial_lj.md",
         ],

docs/src/examples/variablecell.md

@@ -0,0 +1,38 @@
+# Silicon variable cell relaxation
+
+In this example we optimise both atomic positions and lattice geometry
+of a rattled silicon structure using a Stillinger-Weber potential.
+
+We build an initial structure by first rattling the positions
+and than the lattice geometry:
+
+```@example silicon
+using AtomsBase
+using AtomsBuilder
+using LinearAlgebra
+using Unitful
+using UnitfulAtomic
+
+silicon_posrattle = rattle!(bulk(:Si, cubic=true) * (2, 2, 2), 0.1u"Å")
+F = I + 1e-3randn(3, 3)
+new_cell_vectors = tuple([F * v for v in cell_vectors(silicon_posrattle)]...)
+silicon_rattle = AbstractSystem(silicon_posrattle; cell_vectors=new_cell_vectors)
+```
+
+and optimise with `variablecell=true`
+
+```@example silicon
+using EmpiricalPotentials
+using GeometryOptimization
+
+sw = StillingerWeber()
+silicon = minimize_energy!(silicon_posrattle, sw;
+                           variablecell=true, verbosity=1, tol_virial=1e-4u"hartree").system
+nothing  # hide
+```
+
+Final structure:
+
+```@example silicon
+silicon
+```

src/GeometryOptimization.jl

@@ -1,18 +1,20 @@
 module GeometryOptimization
 
+using AtomsBase
+using AtomsCalculators
 using DocStringExtensions
 using LinearAlgebra
-using StaticArrays
 using Optimization
-using AtomsBase
-using AtomsCalculators
+using StaticArrays
 using Unitful
 using UnitfulAtomic
 
 # Make sure Optim is always available
 using OptimizationOptimJL
 using LineSearches
 
+# Useful shortcuts
+using AtomsCalculators: Energy, Forces, Virial
 AC = AtomsCalculators
 
 @template METHODS =
@@ -22,7 +24,7 @@ $(TYPEDSIGNATURES)
 $(DOCSTRING)
 """
 
-include("clamping_updating_positions.jl")
+include("dof_management.jl")
 include("optimization.jl")
 include("callbacks.jl")
 

src/callbacks.jl

@@ -11,10 +11,12 @@ produces output as well.
 struct GeoOptDefaultCallback
     verbosity::Int
     always_show_header::Bool
+    show_virial::Bool
     prev_time::Ref{UInt64}
 end
-function GeoOptDefaultCallback(verbosity=1; always_show_header=verbosity > 1)
-    GeoOptDefaultCallback(verbosity, always_show_header, Ref{UInt64}(0))
+function GeoOptDefaultCallback(verbosity=1;
+                               show_virial=true, always_show_header=verbosity > 1)
+    GeoOptDefaultCallback(verbosity, always_show_header, show_virial, Ref{UInt64}(0))
 end
 
 format_log8(value) = (value < 0 ? " " : "+") * (@sprintf "%8.2f" log10(abs(value)))
@@ -29,7 +31,11 @@ function (cb::GeoOptDefaultCallback)(optim_state, geoopt_state)
     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))
+    if iszero(geoopt_state.energy_change) && optim_state.iter < 1
+        logΔE = ""
+    else
+        logΔE = format_log8(austrip(geoopt_state.energy_change))
+    end
 
     maxforce = austrip(maximum(norm, geoopt_state.forces))
     fstr = iszero(maxforce) ? "" : round(maxforce, sigdigits=8)
@@ -38,13 +44,20 @@ function (cb::GeoOptDefaultCallback)(optim_state, geoopt_state)
         ("n",           3, optim_state.iter),
         ("Energy",     15, Estr),
         ("log10(ΔE)",   9, logΔE),
-        ("max(Force)", 10, fstr),
+        ("max(Force)", 11, fstr),
         # 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 cb.show_virial
+        maxvirial = austrip(maximum(abs, geoopt_state.virial))
+        pressure  = -austrip(tr(geoopt_state.virial)) / 3
+        vstr = iszero(maxvirial) ? "" : round(maxvirial, sigdigits=8)
+        pstr = iszero(pressure)  ? "" : round(pressure,  sigdigits=2)
+        push!(fields, ("max(Virial)", 11, vstr))
+        push!(fields, ("Pressure",     8, pstr))
+    end
+    push!(fields, ("Δ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 cb.always_show_header
         hlines = :all