Implemented suggested fixes for PR.

Project.toml

@@ -9,12 +9,10 @@ AtomsCalculators = "a3e0e189-c65a-42c1-833c-339540406eb1"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [compat]
-julia = "1.9"
 AtomsBase = "0.3"
 AtomsCalculators = "0.1"
 Optimization = "3.20"
@@ -23,9 +21,18 @@ StaticArrays = "1.8"
 TestItemRunner = "0.2"
 Unitful = "1.19"
 UnitfulAtomic = "1.0"
+julia = "1.9"
+DFTK = "0.6"
+ASEconvert = "0.1"
+EmpiricalPotentials = "0.0.1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
+DFTK = "acf6eb54-70d9-11e9-0013-234b7a5f5337"
+ASEconvert = "3da9722f-58c2-4165-81be-b4d7253e8fd2"
+EmpiricalPotentials = "38527215-9240-4c91-a638-d4250620c9e2"
 
 [targets]
-test = ["Test"]
+examples = ["DFTK", "ASEconvert", "EmpiricalPotentials"]
+test = ["Test", "TestItemRunner"]

examples/Al_supercell.jl

@@ -1,6 +1,5 @@
 #= Test Geometry Optimization on an aluminium supercell.
 =#
-using Printf
 using LinearAlgebra
 using DFTK
 using ASEconvert
@@ -13,13 +12,12 @@ using OptimizationOptimJL
 
 using GeometryOptimization
 
-# Get PseudoDojo pseudopotential.
-psp_upf  = load_psp(artifact"pd_nc_sr_lda_standard_0.4.1_upf/Al.upf");
 
 function build_al_supercell(rep=1)
+    pseudodojo_psp = artifact"pd_nc_sr_lda_standard_0.4.1_upf/Al.upf"
     a = 7.65339 # true lattice constant.
     lattice = a * Matrix(I, 3, 3)
-    Al = ElementPsp(:Al; psp=psp_upf)
+    Al = ElementPsp(:Al; psp=load_psp(pseudodojo_psp))
     atoms     = [Al, Al, Al, Al]
     positions = [[0.0, 0.0, 0.0], [0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]]
     unit_cell = periodic_system(lattice, atoms, positions)
@@ -32,29 +30,29 @@ function build_al_supercell(rep=1)
 
     # Unfortunately right now the conversion to ASE drops the pseudopotential information,
     # so we need to reattach it:
-    supercell = attach_psp(supercell; Al=artifact"pd_nc_sr_lda_standard_0.4.1_upf/Al.upf")
+    supercell = attach_psp(supercell; Al=pseudodojo_psp)
     return supercell
 end;
 
 al_supercell = build_al_supercell(1)
 
 # Create a simple calculator for the model.
 model_kwargs = (; functionals = [:lda_x, :lda_c_pw], temperature = 1e-4)
-basis_kwargs = (; kgrid = [8, 8, 8], Ecut = 30.0)
+basis_kwargs = (; kgrid = [6, 6, 6], Ecut = 30.0)
 scf_kwargs = (; tol = 1e-6)
-calculator = DFTKCalculator(al_supercell; model_kwargs, basis_kwargs, scf_kwargs, verbose_scf=true)
+calculator = DFTKCalculator(al_supercell; model_kwargs, basis_kwargs, scf_kwargs, verbose=true)
 
 energy_true = AtomsCalculators.potential_energy(al_supercell, calculator)
 
 # Starting position is a random perturbation of the equilibrium one.
 Random.seed!(1234)
 x0 = vcat(position(al_supercell)...)
 σ = 0.5u"angstrom"; x0_pert = x0 + σ * rand(Float64, size(x0))
-al_supercell = update_optimizable_coordinates(al_supercell, x0_pert)
+al_supercell = update_not_clamped_positions(al_supercell, x0_pert)
 energy_pert = AtomsCalculators.potential_energy(al_supercell, calculator)
 
-@printf "Initial guess distance (norm) from true parameters %.3e bohrs.\n" austrip(norm(x0 - x0_pert))
-@printf "Initial regret %.3e.\n" energy_pert - energy_true
+println("Initial guess distance (norm) from true parameters $(norm(x0 - x0_pert)).")
+println("Initial regret $(energy_pert - energy_true).")
 
 optim_options = (f_tol=1e-6, iterations=6, show_trace=true)
 

examples/H2.jl

@@ -18,7 +18,7 @@ positions = [[0, 0, 0], [0, 0, .16]]
 system = periodic_system(lattice, atoms, positions)
 
 # Set everything to optimizable.
-# system = clamp_atoms(system, [1])
+system = clamp_atoms(system, [1])
 
 # Create a simple calculator for the model.
 model_kwargs = (; functionals = [:lda_x, :lda_c_pw])
@@ -31,4 +31,4 @@ optim_options = (f_tol=1e-32, iterations=20, show_trace=true)
 
 results = optimize_geometry(system, calculator; solver=solver, optim_options...)
 println(results)
-@printf "Bond length: %3f bohrs.\n" norm(results.minimizer[1:3] - results.minimizer[4:end])
+@printf "Bond length: %3f bohrs.\n" norm(results.minimizer[1:end])

examples/H2_LJ.jl

@@ -1,9 +1,7 @@
 #= Test Geometry Optimization on an aluminium supercell.
 =#
-using Printf
 using LinearAlgebra
 using EmpiricalPotentials
-using DFTK
 using Unitful
 using UnitfulAtomic
 using OptimizationOptimJL
@@ -23,4 +21,4 @@ solver = OptimizationOptimJL.LBFGS()
 optim_options = (f_tol=1e-6, iterations=100, show_trace=false)
 
 results = optimize_geometry(system, lj; solver=solver, optim_options...)
-@printf "Bond length: %3f bohrs.\n" norm(results.minimizer[1:3] - results.minimizer[4:end])
+println("Bond length: $(norm(results.minimizer[1:3] - results.minimizer[4:end])).")

src/atomsbase_interface.jl

@@ -1,73 +1,58 @@
 #
-# Interface between AtomsBase and GeomOpt that provides 
+# Interface between AtomsBase.jl and GeometryOptimization.jl that provides 
 # utility functions for manipulating systems.
 #
-# This interface helps defining which particles in the system are considered 
-# optimizable. Note that by default all particles are assumed to be optimizable. 
-# The user can call clamp_atoms to fix atoms whose position should not be optimized.
-#
 # IMPORTANT: Note that we always work in cartesian coordinates.
 #
-export update_positions, update_optimizable_positions, clamp_atoms
+export update_positions, update_not_clamped_positions, clamp_atoms
 
 
 @doc raw"""
 Creates a new system based on ``system`` but with atoms positions updated 
-to the ones specified in `positions`, using cartesian coordinates.
+to the ones provided.
+
 """
 function update_positions(system, positions::AbstractVector{<:AbstractVector{<:Unitful.Length}})
     particles = [Atom(atom; position) for (atom, position) in zip(system, positions)]
     AbstractSystem(system; particles)
 end
 
 @doc raw"""
-Creates a new system based on ``system`` with the optimizable coordinates are
-updated to the ones provided (in the order in which they appear in the system),
-cartesian coordinates version..
+Creates a new system based on ``system`` where the non clamped positions are
+updated to the ones provided (in the order in which they appear in the system).
+
 """
-function update_optimizable_positions(system, positions::AbstractVector{<:Unitful.Length})
-    mask = get_optimizable_mask(system)
+function update_not_clamped_positions(system, positions::AbstractVector{<:Unitful.Length})
+    mask = not_clamped_mask(system)
     new_positions = deepcopy(position(system))
     new_positions[mask] = reinterpret(reshape, SVector{3, eltype(positions)},
                                       reshape(positions, 3, :))
     update_positions(system, new_positions)
 end
 
 @doc raw"""
-Sets the mask defining which coordinates of the system can be optimized. 
-The mask is a vector of booleans, specifying which of the atoms can be optimized.
-
-By default (when no mask is specified), all particles are assumed optimizable.
+Returns a mask for selecting the not clamped atoms in the system.
 
 """
-function set_optimizable_mask(system, mask::AbstractVector{<:Bool})
-    particles = [Atom(atom; optimizable=m) for (atom, m) in zip(system, mask)]
-    AbstractSystem(system, particles=particles)
+function not_clamped_mask(system)
+    # If flag not set, the atom is considered not clamped.
+    [haskey(a, :clamped) ? !a[:clamped] : true for a in system]
 end
 
-@doc raw"""
-    get_optimizable_mask(system) -> AbstractVector{<:Bool}
-
-Returns the optimizable mask of the system (see documentation for `set_optimizable_mask`.
-"""
-function get_optimizable_mask(system)
-    # If flag not set, the atom is considered to be optimizable.
-    [haskey(a, :optimizable) ? a[:optimizable] : true for a in system]
-end
-
-function get_optimizable_positions(system)
-    mask = get_optimizable_mask(system)
-    collect(Iterators.flatten(system[mask, :position]))
+function not_clamped_positions(system)
+    mask = not_clamped_mask(system)
+    Iterators.flatten(system[mask, :position])
 end
 
 @doc raw"""
-    Clamp given atoms if the system. Clamped atoms are fixed and their positions 
+    Clamp given atoms in the system. Clamped atoms are fixed and their positions 
     will not be optimized. The atoms to be clamped should be given as a list of 
-    indies corresponding to their positions in the system.
+    indices corresponding to their positions in the system.
 
     """
 function clamp_atoms(system, clamped_indexes::Union{AbstractVector{<:Integer},Nothing})
-    mask = trues(length(system))
-    mask[clamped_indexes] .= false
-    set_optimizable_mask(system, mask)  # Clamp by setting the mask.
+    clamped = falses(length(system))
+    clamped[clamped_indexes] .= true
+    particles = [Atom(atom; clamped=m) for (atom, m) in zip(system, clamped)]
+    AbstractSystem(system, particles=particles)
 end

src/optimization.jl

@@ -1,15 +1,6 @@
 #=
-For the moment, optimization is delegated to Optim.jl (hardcoded dependency). 
-In order for Optim.jl to use gradients that are returned directly upon function 
-evaluation, one has to use the fix described at: 
-https://github.com/JuliaNLSolvers/Optim.jl/blob/master/docs/src/user/tipsandtricks.md
-
-This is the reason for the cumbersome implementation of optimization with gradients.
-
-Note that by default all particles in the system are assumed optimizable.
-
-IMPORTANT: Note that we always work in cartesian coordinates.
-
+# Note that by default all particles in the system are assumed optimizable.
+# IMPORTANT: Note that we always work in cartesian coordinates.
 =#
 
 export optimize_geometry
@@ -22,16 +13,16 @@ export optimize_geometry
 
 """
 function Optimization.OptimizationFunction(system::AbstractSystem, calculator; kwargs...)
-    mask = get_optimizable_mask(system) # mask is assumed not to change during optim.
+    mask = not_clamped_mask(system)  # mask is assumed not to change during optim.
 
     f = function(x::AbstractVector{<:Real}, p)
-        new_system = update_optimizable_positions(system, x * u"bohr")
+        new_system = update_not_clamped_positions(system, x * u"bohr")
         energy = AtomsCalculators.potential_energy(new_system, calculator; kwargs...)
         austrip(energy)
     end
 
     g! = function(G::AbstractVector{<:Real}, x::AbstractVector{<:Real}, p)
-        new_system = update_optimizable_positions(system, x * u"bohr")
+        new_system = update_not_clamped_positions(system, x * u"bohr")
         energy = AtomsCalculators.potential_energy(new_system, calculator; kwargs...)
 
         forces = AtomsCalculators.forces(new_system, calculator; kwargs...)
@@ -46,8 +37,8 @@ end
 
 function optimize_geometry(system::AbstractSystem, calculator; solver=Optim.LBFGS(), kwargs...)
     # Use current system parameters as starting positions.
-    x0 = Vector(austrip.(get_optimizable_positions(system))) # Optim modifies x0 in-place, so need a mutable type.
+    x0 = austrip.(not_clamped_positions(system)) # Optim modifies x0 in-place, so need a mutable type.
     f_opt = OptimizationFunction(system, calculator)
-    problem = OptimizationProblem(f_opt, x0, nothing) # Last argument needed in Optimization.jl.
+    problem = OptimizationProblem(f_opt, x0, nothing)  # Last argument needed in Optimization.jl.
     solve(problem, solver; kwargs...)
 end

test/dummy_calculator.jl

@@ -4,14 +4,15 @@
     using Unitful
     using UnitfulAtomic
 
-    struct DummyCalculator
-    end
+    struct DummyCalculator end
 
-    AtomsCalculators.@generate_interface function AtomsCalculators.potential_energy(system, calculator::DummyCalculator; kwargs...)
-        return 0.0u"eV"
+    AtomsCalculators.@generate_interface function AtomsCalculators.potential_energy(
+            system, calculator::DummyCalculator; kwargs...)
+        0.0u"eV"
     end
 
-    AtomsCalculators.@generate_interface function AtomsCalculators.forces(system, calculator::DummyCalculator; kwargs...)
-        return AtomsCalculators.zero_forces(system, calculator)
+    AtomsCalculators.@generate_interface function AtomsCalculators.forces(
+            system, calculator::DummyCalculator; kwargs...)
+        AtomsCalculators.zero_forces(system, calculator)
     end
 end

test/geometry_optimization.jl

@@ -11,6 +11,7 @@
     bounding_box = 10.0u"angstrom" .* [[1, 0, 0.], [0., 1, 0], [0., 0, 1]]
     atoms = [:H => [0, 0, 0.0]u"bohr", :H => [0, 0, 1.9]u"bohr"]
     system = periodic_system(atoms, bounding_box)
+    system = clamp_atoms(system, [1])
 
     calculator = TestCalculators.DummyCalculator()