Changed kwargs in minimize enery and fixed bug in position updating.

examples/Al_supercell.jl

@@ -1,5 +1,3 @@
-#= Test Geometry Optimization on an aluminium supercell.
-=#
 using LinearAlgebra
 using DFTK
 using ASEconvert

examples/Al_variable_cell.jl

@@ -43,4 +43,4 @@ scf_kwargs = (; tol = 1e-3)
 calculator = DFTKCalculator(; model_kwargs, basis_kwargs, scf_kwargs, verbose=true)
 
 optim_options = (f_tol=1e-6, iterations=6, show_trace=true)
-results = minimize_energy!(al_supercell, calculator; procedure="vc_relax", optim_options...)
+results = minimize_energy!(al_supercell, calculator; relaxation="unit_cell", optim_options...)

examples/H2.jl

@@ -24,9 +24,9 @@ basis_kwargs = (; kgrid = [1, 1, 1], Ecut = 10.0)
 scf_kwargs = (; tol = 1e-7)
 calculator = DFTKCalculator(; model_kwargs, basis_kwargs, scf_kwargs, verbose=true)
 
-solver = OptimizationOptimJL.LBFGS()
+optimizermizationOptimJL.LBFGS()
 optim_options = (f_tol=1e-32, iterations=20, show_trace=true)
 
-results = minimize_energy!(system, calculator; solver=solver, optim_options...)
+results = minimize_energy!(system, calculator; optimizer, optim_options...)
 println(results)
 @printf "Bond length: %3f bohrs.\n" norm(results.minimizer[1:end])

examples/H2_LJ.jl

@@ -1,5 +1,3 @@
-#= Test Geometry Optimization on an aluminium supercell.
-=#
 using LinearAlgebra
 using EmpiricalPotentials
 using Unitful
@@ -16,8 +14,8 @@ system = periodic_system(atoms, bounding_box)
 
 lj = LennardJones(-1.17u"hartree", 0.743u"angstrom", 1, 1, 0.6u"nm")
 
-solver = OptimizationOptimJL.LBFGS()
+optimizer = OptimizationOptimJL.LBFGS()
 optim_options = (f_tol=1e-6, iterations=100, show_trace=false)
 
-results = minimize_energy!(system, lj; solver=solver, optim_options...)
+results = minimize_energy!(system, lj; optimizer, optim_options...)
 println("Bond length: $(norm(results.minimizer[1:3] - results.minimizer[4:end])).")

examples/TiAl_EmpiricalPotentials.jl

@@ -19,8 +19,8 @@ particles = map(data) do atom
 end
 system = AbstractSystem(data; particles)
 
-solver = OptimizationOptimJL.LBFGS()
+optimizer = OptimizationOptimJL.LBFGS()
 optim_options = (f_tol=1e-8, g_tol=1e-8, iterations=10, show_trace=true)
 
-results = minimize_energy!(system, lj; solver, optim_options...)
+results = minimize_energy!(system, lj; optimizer, optim_options...)
 println(results)

src/atomsbase_interface.jl

@@ -32,8 +32,9 @@ function update_positions(system, positions::ComponentVector)
     # Compatibility with Optimization.jl (which is not able to handle units), 
     # requires us to artificially attach units to strain upstream. We here remove them.
     deformation_tensor = I + voigt_to_full(austrip.(positions.strain))
-    # TODO: Do we want to apply the strain to the atoms too?
-    particles = [Atom(atom; position = deformation_tensor * position) for (atom, position)
+
+    # auconvert needed here since multiplication drops the type of the array (elements still typed)
+    particles = [Atom(atom; position = auconvert.(deformation_tensor * position)) for (atom, position)
                  in zip(system, collect.(positions.atoms))]
 
     bbox = eachcol(deformation_tensor * bbox_to_matrix(bounding_box(system)))

src/optimization.jl

@@ -19,7 +19,7 @@ By default we work in cartesian coordinates.
 Note that internally, when optimizing the cartesian positions, atomic units 
 are used.
 """
-function Optimization.OptimizationFunction(system, calculator; pressure=0.0, kwargs...)
+function Optimization.OptimizationFunction(system, calculator; kwargs...)
     mask = not_clamped_mask(system)  # mask is assumed not to change during optim.
 
     # TODO: Note that this function will dispatch appropriately when called with 
@@ -60,19 +60,33 @@ function Optimization.OptimizationFunction(system, calculator; pressure=0.0, kwa
     OptimizationFunction(f; grad=g!)
 end
 
-function minimize_energy!(system, calculator; pressure=0.0, procedure="relax", 
-                          solver=Optim.LBFGS(), kwargs...)
+@doc raw"""
+    minimize_energy!(system, calculator; [relaxation, optimizer])
+
+Minimize the system's energy with respect to it's internal degrees of freedom 
+(geometry optimization). Degrees of freedom can be either the atomic coordinates 
+only, or the atomic coordinates together with the unit cell vectors.
+
+Keyword arguments:
+- `relaxation`: Defines if only atoms or atoms and unit cell vectors are optimized.
+                Must be one of `["atoms", "unit_cell"]`.
+- `optimizer`: Optimizer to use. Use one from Optimization.jl. Defaults to LBFGS.
+
+Returns a named tuple (; optimized_system, optimization_data).
+
+"""
+function minimize_energy!(system, calculator; relaxation="atoms", optimizer=Optim.LBFGS(), kwargs...)
     # Use current system parameters as starting positions.
-    if procedure == "relax"
+    if relaxation == "atoms"
         x0 = austrip.(not_clamped_positions(system))
-    elseif procedure == "vc_relax"
+    elseif relaxation == "unit_cell"
         x0 = ComponentVector(atoms = austrip.(reduce(vcat, position(system))), strain = zeros(6))
     else
-        print("Error: unknown optimization procedure. Please use one of [`relax`, `vc_relax`].")
+        print("Error: unknown relaxation type. Please use one of [`atoms`, `unit_cell`].")
     end
-    f_opt = OptimizationFunction(system, calculator; pressure)
+    f_opt = OptimizationFunction(system, calculator)
     problem = OptimizationProblem(f_opt, x0, nothing)  # Last argument needed in Optimization.jl.
-    optimization_data = solve(problem, solver; kwargs...)
+    optimization_data = solve(problem, optimizer; kwargs...)
 
     # Return final structure.
     optimized_system = update_not_clamped_positions(system, optimization_data.u * u"bohr")

test/geometry_optimization.jl

@@ -15,9 +15,9 @@
 
     calculator = TestCalculators.DummyCalculator()
 
-    solver = OptimizationOptimJL.LBFGS()
+    optimizer = OptimizationOptimJL.LBFGS()
     optim_options = (f_tol=1e-6, iterations=4, show_trace=false)
 
-    results = minimize_energy!(system, calculator; solver=solver, optim_options...)
+    results = minimize_energy!(system, calculator; optimizer, optim_options...)
     @test isapprox(results.u[1], 0; atol=1e-5)
 end