Basic implementation of geometry optimization interface.

examples/Al_supercell.jl

@@ -32,7 +32,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="hgh/lda/al-q3")
+    supercell = attach_psp(supercell; Al=artifact"pd_nc_sr_lda_standard_0.4.1_upf/Al.upf")
     return supercell
 end;
 

src/atomsbase_interface.jl

@@ -8,73 +8,66 @@
 #
 # IMPORTANT: Note that we always work in cartesian coordinates.
 #
-export update_positions, update_optimizable_coordinates, clamp_atoms
+export update_positions, update_optimizable_positions, clamp_atoms
 
 
 @doc raw"""
-    update_postions(system::AbstractSystem, positions::Vector{<:AbstractVector{<:Real}}) where {L <: Unitful.Length}
-
-Creates a new system based on ``system`` but with atoms positions updated to the ones specified in `positions`, 
-using cartesian coordinates.
+Creates a new system based on ``system`` but with atoms positions updated 
+to the ones specified in `positions`, using cartesian coordinates.
 """
-function update_positions(system::AbstractSystem, positions::AbstractVector{<:AbstractVector{L}}) where {L <: Unitful.Length}
-    particles = [
-                 Atom(atom, position=SVector{3,L}(position)) 
-                 for (atom, position) in zip(system.particles, positions)]
-    AbstractSystem(system, particles=particles)
+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"""
-    update_optimizable_coordinates(system::AbstractSystem, positions::Vector{<:AbstractVector{<:L}}) where {L <: Unitful.Length}
-
 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 system.particles), cartesian coordinates version..
+updated to the ones provided (in the order in which they appear in the system),
+cartesian coordinates version..
 """
-function update_optimizable_coordinates(system::AbstractSystem, positions::AbstractVector{<:L}) where {L <: Unitful.Length}
+function update_optimizable_positions(system, positions::AbstractVector{<:Unitful.Length})
     mask = get_optimizable_mask(system)
-    new_positions = Vector.(position(system)) # make mutable.
-    new_positions[mask] = [eachcol(reshape(positions, 3, :))...]
+    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"""
-    set_optimizable_mask(system::AbstractSystem, mask::AbstractVector{<:Bool})
-
 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.
 
 """
-function set_optimizable_mask(system::AbstractSystem, mask::AbstractVector{<:Bool})
-    particles = [Atom(atom; optimizable=m) for (atom, m) in zip(system.particles, mask)]
+function set_optimizable_mask(system, mask::AbstractVector{<:Bool})
+    particles = [Atom(atom; optimizable=m) for (atom, m) in zip(system, mask)]
     AbstractSystem(system, particles=particles)
 end
 
 @doc raw"""
-    get_optimizable_mask(system::AbstractSystem) -> AbstractVector{<:Bool}
+    get_optimizable_mask(system) -> AbstractVector{<:Bool}
 
 Returns the optimizable mask of the system (see documentation for `set_optimizable_mask`.
 """
-function get_optimizable_mask(system::AbstractSystem)
+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.particles]
+    [haskey(a, :optimizable) ? a[:optimizable] : true for a in system]
 end
 
-function get_optimizable_coordinates(system::AbstractSystem)
+function get_optimizable_positions(system)
     mask = get_optimizable_mask(system)
-    return collect(Iterators.flatten(system[mask, :position]))
+    collect(Iterators.flatten(system[mask, :position]))
 end
 
 @doc raw"""
     Clamp given atoms if 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 system.particles.
+    indies corresponding to their positions in the system.
 
     """
-function clamp_atoms(system::AbstractSystem, clamped_indexes::Union{AbstractVector{<:Integer},Nothing})
-    mask = trues(length(system.particles))
+function clamp_atoms(system, clamped_indexes::Union{AbstractVector{<:Integer},Nothing})
+    mask = trues(length(system))
     mask[clamped_indexes] .= false
-    clamped_system = set_optimizable_mask(system, mask)
-    return clamped_system
+    set_optimizable_mask(system, mask)  # Clamp by setting the mask.
 end

src/optimization.jl

@@ -25,15 +25,13 @@ function Optimization.OptimizationFunction(system::AbstractSystem, calculator; k
     mask = get_optimizable_mask(system) # mask is assumed not to change during optim.
 
     f = function(x::AbstractVector{<:Real}, p)
-        x = 1u"bohr" .* x # Work in atomic units.
-        new_system = update_optimizable_coordinates(system, x)
+        new_system = update_optimizable_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)
-        x = 1u"bohr" .* x # Work in atomic units.
-        new_system = update_optimizable_coordinates(system, x)
+        new_system = update_optimizable_positions(system, x * u"bohr")
         energy = AtomsCalculators.potential_energy(new_system, calculator; kwargs...)
 
         forces = AtomsCalculators.forces(new_system, calculator; kwargs...)
@@ -48,7 +46,7 @@ end
 
 function optimize_geometry(system::AbstractSystem, calculator; solver=Optim.LBFGS(), kwargs...)
     # Use current system parameters as starting positions.
-    x0 = Vector(austrip.(get_optimizable_coordinates(system))) # Optim modifies x0 in-place, so need a mutable type.
+    x0 = Vector(austrip.(get_optimizable_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.
     solve(problem, solver; kwargs...)