diff --git a/Project.toml b/Project.toml index 98b8ceb..11547d8 100644 --- a/Project.toml +++ b/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"] diff --git a/examples/Al_supercell.jl b/examples/Al_supercell.jl index 69a4483..b906cee 100644 --- a/examples/Al_supercell.jl +++ b/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,7 +30,7 @@ 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; @@ -40,9 +38,9 @@ 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) @@ -50,11 +48,11 @@ energy_true = AtomsCalculators.potential_energy(al_supercell, calculator) 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) diff --git a/examples/H2.jl b/examples/H2.jl index c3aea99..65fbfe1 100644 --- a/examples/H2.jl +++ b/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]) diff --git a/examples/H2_LJ.jl b/examples/H2_LJ.jl index 1e7e864..6343b37 100644 --- a/examples/H2_LJ.jl +++ b/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])).") diff --git a/src/atomsbase_interface.jl b/src/atomsbase_interface.jl index d98b900..22decca 100644 --- a/src/atomsbase_interface.jl +++ b/src/atomsbase_interface.jl @@ -1,19 +1,16 @@ # -# 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)] @@ -21,12 +18,12 @@ function update_positions(system, positions::AbstractVector{<:AbstractVector{<:U 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, :)) @@ -34,40 +31,28 @@ function update_optimizable_positions(system, positions::AbstractVector{<:Unitfu 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 diff --git a/src/optimization.jl b/src/optimization.jl index a551f09..b69e5eb 100644 --- a/src/optimization.jl +++ b/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 diff --git a/test/dummy_calculator.jl b/test/dummy_calculator.jl index c32f077..ba29c0e 100644 --- a/test/dummy_calculator.jl +++ b/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 diff --git a/test/geometry_optimization.jl b/test/geometry_optimization.jl index b8bf797..da6c99d 100644 --- a/test/geometry_optimization.jl +++ b/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()