geom_opt.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2007 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module geom_opt_oct_m
  use allelectron_oct_m
  use debug_oct_m
  use density_oct_m
  use electrons_oct_m
  use energy_calc_oct_m
  use forces_oct_m
  use global_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use io_function_oct_m
  use ion_dynamics_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use lcao_oct_m
  use loct_oct_m
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use minimizer_oct_m
  use mpi_oct_m
  use multisystem_basic_oct_m
  use namespace_oct_m
  use parser_oct_m
  use pcm_oct_m
  use profiling_oct_m
  use read_coords_oct_m
  use restart_oct_m
  use scf_oct_m
  use species_oct_m
  use states_elec_oct_m
  use states_elec_restart_oct_m
  use symmetries_oct_m
  use symmetrizer_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use v_ks_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public :: geom_opt_run
 
  type geom_opt_t
    private
    integer(int64) :: type
    integer        :: method
    real(real64)   :: step
    real(real64)   :: line_tol
    real(real64)   :: fire_mass
    integer        :: fire_integrator
    real(real64)   :: tolgrad
    real(real64)   :: toldr
    integer        :: max_iter
    integer        :: what2minimize
 
    type(scf_t)                       :: scfv
    type(ions_t),             pointer :: ions
    type(hamiltonian_elec_t), pointer :: hm
    type(electrons_t),        pointer :: syst
    class(mesh_t),            pointer :: mesh
    type(states_elec_t),      pointer :: st
    integer                           :: dim
    integer                           :: periodic_dim
    integer                           :: size
    integer                           :: fixed_atom = 0
 
    real(real64), allocatable         :: cell_force(:, :)
    logical                           :: symmetrize = .false.
    real(real64), allocatable         :: initial_length(:)
    real(real64), allocatable         :: initial_rlattice(:, :)
    real(real64), allocatable         :: inv_initial_rlattice(:, :)
    real(real64)                      :: pressure = m_zero 
 
    logical                           :: poscar_output = .false.
 
    integer                           :: calculator
    logical                           :: ml_preopt = .false.
  end type geom_opt_t
 
  type(geom_opt_t), save :: g_opt
 
  integer, parameter :: &
    MINWHAT_ENERGY = 1, &
    minwhat_forces = 2
 
  integer, parameter ::   &
    GO_IONS   = 1,        &
    go_cell   = 2,        &
    go_volume = 4
 
  integer, parameter ::           &
    GO_CALCULATOR_INTERNAL  = 1, &
    go_calculator_mattersim = 2
 
contains
 
  ! ---------------------------------------------------------
  subroutine mattersim_paths(mattersim_share, mattersim_python, namespace)
    character(len=*), intent(out) :: mattersim_share, mattersim_python
    type(namespace_t), intent(in) :: namespace
 
    logical :: python_exists
 
    ! In build-tree runs conf%share points to build/share, but the
    ! MatterSim virtual environment is created when installing Octopus.
    mattersim_share = trim(conf%share)
    mattersim_python = trim(mattersim_share) // "/mattersim_env/bin/python"
    inquire(file=trim(mattersim_python), exist=python_exists)
    if (python_exists) return
 
    mattersim_share = share_dir
    mattersim_python = trim(mattersim_share) // "/mattersim_env/bin/python"
    inquire(file=trim(mattersim_python), exist=python_exists)
    if (python_exists) return
 
    message(1) = 'Unable to find the MatterSim Python environment.'
    message(2) = 'Please install Octopus with MatterSim support before using the MatterSim calculator.'
    call messages_fatal(2, namespace=namespace)
  end subroutine mattersim_paths
 
  ! ---------------------------------------------------------
  subroutine geom_opt_run(system, from_scratch)
    class(*),        intent(inout) :: system
    logical,         intent(inout) :: from_scratch
 
    push_sub(geom_opt_run)
 
    select type (system)
    class is (multisystem_basic_t)
      message(1) = "CalculationMode = go not implemented for multi-system calculations"
      call messages_fatal(1, namespace=system%namespace)
    type is (electrons_t)
      call geom_opt_run_legacy(system, from_scratch)
    end select
 
    pop_sub(geom_opt_run)
  end subroutine geom_opt_run
 
  ! ---------------------------------------------------------
  subroutine geom_opt_run_legacy(sys, fromscratch)
    type(electrons_t),   target, intent(inout) :: sys
    logical,                     intent(inout) :: fromscratch
 
    integer :: ierr
    real(real64), allocatable :: coords(:)
    real(real64) :: energy
 
    real(real64), allocatable :: mass(:)
    integer :: iatom, imass
    type(restart_t) :: restart_load
    logical :: known_lower_bound
    logical :: default_val
 
    push_sub(geom_opt_run_legacy)
 
    ! For 1D-periodic systems the ion-ion interaction (hence the forces) and the stress
    ! tensor are not implemented.
    if (sys%space%periodic_dim == 1) then
      message(1) = "Geometry optimization is not allowed for systems periodic in 1D, "
      message(2) = "as the ion-ion interaction and the stress tensor are not available."
      call messages_fatal(2, namespace=sys%namespace)
    end if
 
    if (sys%space%periodic_dim == 2) then
      call messages_experimental("Geometry optimization for 2D-periodic systems", namespace=sys%namespace)
    end if
 
 
    if (sys%hm%pcm%run_pcm) then
      call messages_not_implemented("PCM for CalculationMode /= gs or td", namespace=sys%namespace)
    end if
 
    if (sys%kpoints%use_symmetries) then
      call messages_experimental("KPoints symmetries with CalculationMode = go", namespace=sys%namespace)
    end if
 
    g_opt%symmetrize = sys%kpoints%use_symmetries .or. sys%st%symmetrize_density
 
    call init_(fromscratch)
 
    ! load wavefunctions
    if (.not. fromscratch) then
      call restart_load%init(sys%namespace, restart_gs, restart_type_load, sys%mc, ierr, mesh=sys%gr)
      if (ierr == 0) then
        call states_elec_load(restart_load, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
          fixed_occ=.false., ierr=ierr)
      end if
      call restart_load%end()
      if (ierr /= 0) then
        message(1) = "Unable to read wavefunctions: Starting from scratch."
        call messages_warning(1, namespace=sys%namespace)
        fromscratch = .true.
      end if
    end if
 
    call scf_init(g_opt%scfv, sys%namespace, sys%gr, sys%ions, sys%st, sys%mc, sys%hm, sys%space)
 
    if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0) then
      if (.not. g_opt%scfv%calc_stress) then
        message(1) = "In order to optimize the cell, one needs to set SCFCalculateStress = yes."
        call messages_fatal(1, namespace=sys%namespace)
      end if
    end if
 
    if (fromscratch) then
      call lcao_run(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, sys%hm, &
        lmm_r = g_opt%scfv%lmm_r, known_lower_bound=known_lower_bound)
      call scf_set_lower_bound_is_known(g_opt%scfv, known_lower_bound)
    else
      ! setup Hamiltonian
      message(1) = 'Info: Setting up Hamiltonian.'
      call messages_info(1, namespace=sys%namespace)
      call v_ks_h_setup(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, sys%hm)
      call scf_set_lower_bound_is_known(g_opt%scfv, .true.)
    end if
 
    !Initial point
    safe_allocate(coords(1:g_opt%size))
    call to_coords(g_opt, coords)
 
    if (sys%st%pack_states .and. sys%hm%apply_packed()) call sys%st%pack()
 
    !Minimize
    select case (g_opt%method)
    case (minmethod_nmsimplex)
      call minimize_multidim_nograd(g_opt%method, g_opt%size, coords, g_opt%step,&
        g_opt%toldr, g_opt%max_iter, &
        calc_point_ng, write_iter_info_ng, energy, ierr)
 
    case (minmethod_fire)
 
      safe_allocate(mass(1:g_opt%size))
      mass = g_opt%fire_mass
      imass = 1
      do iatom = 1, sys%ions%natoms
        if (g_opt%fixed_atom == iatom) cycle
        if (g_opt%ions%fixed(iatom)) cycle
        if (g_opt%fire_mass <= m_zero) mass(imass:imass + 2) = sys%ions%mass(iatom)
        imass = imass + g_opt%dim
      end do
 
      !TODO: add variable to use Euler integrator
      call minimize_fire(g_opt%size, g_opt%ions%space%dim, coords, g_opt%step, g_opt%tolgrad, &
        g_opt%max_iter, calc_point, write_iter_info, energy, ierr, mass, integrator=g_opt%fire_integrator)
 
      safe_deallocate_a(mass)
 
    case default
      call minimize_multidim(g_opt%method, g_opt%size, coords, g_opt%step ,&
        g_opt%line_tol , g_opt%tolgrad, g_opt%toldr, g_opt%max_iter, &
        calc_point, write_iter_info, energy, ierr)
 
    end select
 
    if (ierr == 1025) then
      ! not a GSL error, set by our minimize routines, so we must handle it separately
      message(1) = "Reached maximum number of iterations allowed by GOMaxIter."
      call messages_info(1, namespace=sys%namespace)
    else if (ierr /= 0 .and. g_opt%method /= minmethod_fire) then
      message(1) = "Error occurred during the GSL minimization procedure:"
      call loct_strerror(ierr, message(2))
      call messages_fatal(2, namespace=sys%namespace)
    end if
 
    if (sys%st%pack_states .and. sys%hm%apply_packed()) call sys%st%unpack()
 
 
    ! print out geometry
    message(1) = "Writing final coordinates to min.xyz"
    call messages_info(1, namespace=sys%namespace)
    call from_coords(g_opt, coords)
    call g_opt%ions%write_xyz('./min')
 
    safe_deallocate_a(coords)
    call scf_end(g_opt%scfv)
    ! Because g_opt has the "save" attribute, we need to explicitly empty the criteria list here, or there will be a memory leak.
    call g_opt%scfv%criterion_list%empty()
    call end_()
 
    pop_sub(geom_opt_run_legacy)
  contains
 
    ! ---------------------------------------------------------
    subroutine init_(fromscratch)
      logical,  intent(inout) :: fromscratch
 
      logical :: center, does_exist
      integer :: iter, iatom, idir
      character(len=100) :: filename
      real(real64) :: default_toldr
      real(real64) :: default_step
      type(read_coords_info) :: xyz
 
      push_sub(geom_opt_run_legacy.init_)
 
      if (sys%space%is_periodic()) then
        call messages_experimental('Geometry optimization for periodic systems', namespace=sys%namespace)
      end if
 
      !%Variable GOType
      !%Type flag
      !%Default ions
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% This variable defines which parameters are allowed to change during the optimization.
      !% Multiple options can be chosen e.g. “ion_positions + cell_shape”.
      !% Only one type of lattice vectors relaxation is possible.
      !%Option ion_positions 1
      !% Relax position of ions based on the forces acting on the ions.
      !%Option cell_shape 2
      !% Relax cell shape. This changes lattice vector lengths and directions
      !% based on the stress acting on the lattice vectors.
      !% See for instance Wentzcovitch, PRB 44, 2358 (1991).
      !%Option cell_volume 4
      !% Relax cell volume. Only allow for rescaling the lengths of lattice vectors.
      !% This is a simplication of the option cell_shape, where only a diagonal strain is allowed.
      !%End
 
      call parse_variable(sys%namespace, 'GOType', go_ions, g_opt%type)
      if (.not. varinfo_valid_option('GOType', g_opt%type, is_flag=.true.)) then
        call messages_input_error(sys%namespace, 'GOType')
      end if
 
      write(message(1),'(a)') 'Input: [GOType = '
      if (bitand(g_opt%type, go_ions) /= 0) then
        write(message(1),'(a,1x,a)') trim(message(1)), 'ion_positions'
      end if
      if (bitand(g_opt%type, go_cell) /= 0) then
        if (len_trim(message(1)) > 16) then
          write(message(1),'(a,1x,a)') trim(message(1)), '+'
        end if
        write(message(1),'(a,1x,a)') trim(message(1)), 'cell_shape'
      end if
      if (bitand(g_opt%type, go_volume) /= 0) then
        if (len_trim(message(1)) > 16) then
          write(message(1),'(a,1x,a)') trim(message(1)), '+'
        end if
        write(message(1),'(a,1x,a)') trim(message(1)), 'cell_volume'
      end if
      write(message(1),'(2a)') trim(message(1)), ']'
      call messages_info(1, namespace=sys%namespace)
 
      if (bitand(g_opt%type, go_volume) /= 0 .and. bitand(g_opt%type, go_cell) /= 0) then
        message(1) = "Cell and volume optimization cannot be used simultaneously."
        call messages_fatal(1, namespace=sys%namespace)
      end if
 
 
      if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0) then
        if (parse_is_defined(sys%namespace, 'TDMomentumTransfer') .or. &
          parse_is_defined(sys%namespace, 'TDReducedMomentumTransfer')) then
          call messages_not_implemented("Cell dynamics with TDMomentumTransfer and TDReducedMomentumTransfer")
        end if
      end if
 
      do iatom = 1, sys%ions%natoms
        select type(spec=>sys%ions%atom(iatom)%species)
        class is(allelectron_t)
          write(message(1),'(a)') "Geometry optimization for all-electron potential is not implemented."
          call messages_fatal(1)
        end select
      end do
 
 
      call states_elec_allocate_wfns(sys%st, sys%gr, packed=.true.)
 
      ! shortcuts
      g_opt%mesh   => sys%gr
      g_opt%ions   => sys%ions
      g_opt%st     => sys%st
      g_opt%hm     => sys%hm
      g_opt%syst   => sys
      g_opt%dim    =  sys%space%dim
      g_opt%periodic_dim    =  sys%space%periodic_dim
 
      g_opt%size = 0
      ! Ion dyamics
      if (bitand(g_opt%type, go_ions) /= 0) then
        g_opt%size = g_opt%dim * g_opt%ions%natoms
      end if
 
      ! Cell dynamics
      if (bitand(g_opt%type, go_cell) /= 0) then
        g_opt%size = g_opt%size + (g_opt%periodic_dim +1) * g_opt%periodic_dim / 2
        safe_allocate(g_opt%cell_force(1:g_opt%periodic_dim, 1:g_opt%periodic_dim))
      end if
 
      ! Volume dynamics
      if (bitand(g_opt%type, go_volume) /= 0) then
        g_opt%size = g_opt%size + g_opt%periodic_dim
        safe_allocate(g_opt%cell_force(1:g_opt%periodic_dim, 1:1))
        ! Store the length of the original lattic vectors, to work with reduced lengthes
        safe_allocate(g_opt%initial_length(1:g_opt%periodic_dim))
        do idir = 1, g_opt%periodic_dim
          g_opt%initial_length(idir) = norm2(g_opt%ions%latt%rlattice(1:g_opt%periodic_dim, idir))
        end do
      end if
 
      ! Store the initial lattice vectors and the inverse matrix
      if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0) then
        safe_allocate(g_opt%initial_rlattice(1:g_opt%periodic_dim, 1:g_opt%periodic_dim))
        g_opt%initial_rlattice(1:g_opt%periodic_dim, 1:g_opt%periodic_dim) &
          = g_opt%ions%latt%rlattice(1:g_opt%periodic_dim, 1:g_opt%periodic_dim)
        safe_allocate(g_opt%inv_initial_rlattice(1:g_opt%periodic_dim, 1:g_opt%periodic_dim))
        g_opt%inv_initial_rlattice(:, :) = g_opt%initial_rlattice(:, :)
        call lalg_inverse(g_opt%periodic_dim, g_opt%inv_initial_rlattice, 'dir')
      end if
 
      if(g_opt%ions%space%is_periodic()) then
        call parse_variable(sys%namespace, 'HydrostaticPressure', m_zero, g_opt%pressure)
      end if
 
      !%Variable GOCenter
      !%Type logical
      !%Default no
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% (Experimental) If set to yes, Octopus centers the geometry at
      !% every optimization step. It also reduces the degrees of
      !% freedom of the optimization by using the translational
      !% invariance.
      !%End
      call parse_variable(sys%namespace, 'GOCenter', .false.,  center)
 
      if (center  .and. bitand(g_opt%type, go_ions) /= 0) then
        g_opt%fixed_atom = 1
        g_opt%size = g_opt%size  - g_opt%dim
        call messages_experimental('GOCenter', namespace=sys%namespace)
      end if
 
      !Check if atoms are allowed to move and redifine g_opt%size
      do iatom = 1, g_opt%ions%natoms
        if (g_opt%ions%fixed(iatom) .and. bitand(g_opt%type, go_ions) /= 0) then
          g_opt%size = g_opt%size  - g_opt%dim
        end if
      end do
 
      assert(g_opt%size > 0)
 
      !%Variable GOMethod
      !%Type integer
      !%Default fire
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% Method by which the minimization is performed. For more information see the
      !% <a href=http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html>
      !% GSL documentation</a>.
      !%Option steep 1
      !% Simple steepest descent.
      !%Option steep_native -1
      !% (Experimental) Non-gsl implementation of steepest descent.
      !%Option cg_fr 2
      !% Fletcher-Reeves conjugate-gradient algorithm. The
      !% conjugate-gradient algorithm proceeds as a succession of line
      !% minimizations. The sequence of search directions is used to build
      !% up an approximation to the curvature of the function in the
      !% neighborhood of the minimum.
      !%Option cg_pr 3
      !% Polak-Ribiere conjugate-gradient algorithm.
      !%Option cg_bfgs 4
      !% Vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) conjugate-gradient algorithm.
      !% It is a quasi-Newton method which builds up an approximation to the second
      !% derivatives of the function <i>f</i> using the difference between successive gradient
      !% vectors.  By combining the first and second derivatives, the algorithm is able
      !% to take Newton-type steps towards the function minimum, assuming quadratic
      !% behavior in that region.
      !%Option cg_bfgs2 5
      !% The bfgs2 version of this minimizer is the most efficient version available,
      !% and is a faithful implementation of the line minimization scheme described in
      !% Fletcher, <i>Practical Methods of Optimization</i>, Algorithms 2.6.2 and 2.6.4.
      !%Option simplex 6
      !% This is experimental, and in fact, <b>not</b> recommended unless you just want to
      !% fool around. It is the Nead-Melder simplex algorithm, as implemented in the
      !% GNU Scientific Library (GSL). It does not make use of the gradients (<i>i.e.</i>, the
      !% forces) which makes it less efficient than other schemes. It is included here
      !% for completeness, since it is free.
      !%Option fire 8
      !% The FIRE algorithm. See also <tt>GOFireMass</tt> and <tt>GOFireIntegrator</tt>.
      !% Ref: E. Bitzek, P. Koskinen, F. Gahler, M. Moseler, and P. Gumbsch, <i>Phys. Rev. Lett.</i> <b>97</b>, 170201 (2006).
      !%End
      call parse_variable(sys%namespace, 'GOMethod', minmethod_fire, g_opt%method)
      if (.not. varinfo_valid_option('GOMethod', g_opt%method)) call messages_input_error(sys%namespace, 'GOMethod')
 
 
      call messages_print_var_option("GOMethod", g_opt%method, namespace=sys%namespace)
 
      !%Variable GOTolerance
      !%Type float
      !%Default 0.001 H/b (0.051 eV/A)
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% Convergence criterion, for stopping the minimization. In
      !% units of force; minimization is stopped when all forces on
      !% ions are smaller than this criterion, or the
      !% <tt>GOMinimumMove</tt> is satisfied. If <tt>GOTolerance < 0</tt>,
      !% this criterion is ignored.
      !%End
      call parse_variable(sys%namespace, 'GOTolerance', 0.001_real64, g_opt%tolgrad, units_inp%force)
 
      !%Variable GOMinimumMove
      !%Type float
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% Convergence criterion, for stopping the minimization. In
      !% units of length; minimization is stopped when the coordinates
      !% of all species change less than <tt>GOMinimumMove</tt>, or the
      !% <tt>GOTolerance</tt> criterion is satisfied.
      !% If <tt>GOMinimumMove < 0</tt>, this criterion is ignored.
      !% Default is -1, except 0.001 b with <tt>GOMethod = simplex</tt>.
      !% Note that if you use <tt>GOMethod = simplex</tt>,
      !% then you must supply a non-zero <tt>GOMinimumMove</tt>.
      !%End
      if (g_opt%method == minmethod_nmsimplex) then
        default_toldr = 0.001_real64
      else
        default_toldr = -m_one
      end if
      call parse_variable(sys%namespace, 'GOMinimumMove', default_toldr, g_opt%toldr, units_inp%length)
 
      if (g_opt%method == minmethod_nmsimplex .and. g_opt%toldr <= m_zero) call messages_input_error(sys%namespace, 'GOMinimumMove')
 
      !%Variable GOStep
      !%Type float
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% Initial step for the geometry optimizer. The default is 0.5.
      !% WARNING: in some weird units.
      !% For the FIRE minimizer, default value is 0.1 fs,
      !% and corresponds to the initial time-step for the MD.
      !%End
      if (g_opt%method /= minmethod_fire) then
        default_step = m_half
        call parse_variable(sys%namespace, 'GOStep', default_step, g_opt%step)
      else
        default_step = 0.1_real64*unit_femtosecond%factor
        call parse_variable(sys%namespace, 'GOStep', default_step, g_opt%step, unit = units_inp%time)
      end if
 
      !%Variable GOLineTol
      !%Type float
      !%Default 0.1
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% Tolerance for line-minimization. Applies only to GSL methods
      !% that use the forces.
      !% WARNING: in some weird units.
      !%End
      call parse_variable(sys%namespace, 'GOLineTol', 0.1_real64, g_opt%line_tol)
 
      !%Variable GOMaxIter
      !%Type integer
      !%Default 200
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% Even if the convergence criterion is not satisfied, the minimization will stop
      !% after this number of iterations.
      !%End
      call parse_variable(sys%namespace, 'GOMaxIter', 200, g_opt%max_iter)
      if (g_opt%max_iter <= 0) then
        message(1) = "GOMaxIter has to be larger than 0"
        call messages_fatal(1, namespace=sys%namespace)
      end if
 
      !%Variable GOFireMass
      !%Type float
      !%Default 1.0 amu
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% The Fire algorithm (<tt>GOMethod = fire</tt>) assumes that all degrees of freedom
      !% are comparable. All the velocities should be on the same
      !% scale,  which  for  heteronuclear  systems  can  be  roughly
      !% achieved by setting all the atom masses equal, to the value
      !% specified by this variable.
      !% By default the mass of a proton is selected (1 amu).
      !% However, a selection of <tt>GOFireMass = 0.01</tt> can, in manys systems,
      !% speed up the geometry optimization procedure.
      !% If <tt>GOFireMass</tt> <= 0, the masses of each
      !% species will be used.
      !%End
      call parse_variable(sys%namespace, 'GOFireMass', m_one*unit_amu%factor, g_opt%fire_mass, unit = unit_amu)
 
      !%Variable GOFireIntegrator
      !%Type integer
      !%Default verlet
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% The Fire algorithm (<tt>GOMethod = fire</tt>) uses a molecular dynamics
      !% integrator to compute new geometries and velocities.
      !% Currently, two integrator schemes can be selected
      !%Option euler 0
      !% The Explicit Euler method.
      !%Option verlet 1
      !% The Velocity Verlet algorithm.
      !%Option semi_implicit_euler 2
      !% Semi-implicit Euler integration, see J. Guénolé, et al. Computational Materials Science 175 (2020) 109584.
      !%End
      call parse_variable(sys%namespace, 'GOFireIntegrator', option__gofireintegrator__verlet, g_opt%fire_integrator)
 
      call messages_obsolete_variable(sys%namespace, 'GOWhat2Minimize', 'GOObjective')
 
      !%Variable GOObjective
      !%Type integer
      !%Default minimize_energy
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% This rather esoteric option allows one to choose which
      !% objective function to minimize during a geometry
      !% minimization. The use of this variable may lead to
      !% inconsistencies, so please make sure you know what you are
      !% doing.
      !%Option minimize_energy 1
      !% Use the total energy as objective function.
      !%Option minimize_forces 2
      !% Use <math>\sqrt{\sum_i \left| f_i \right|^2}</math> as objective function.
      !% Note that in this case one still uses the forces as the gradient of the objective function.
      !% This is, of course, inconsistent, and may lead to very strange behavior.
      !%End
      call parse_variable(sys%namespace, 'GOObjective', minwhat_energy, g_opt%what2minimize)
      if (.not. varinfo_valid_option('GOObjective', g_opt%what2minimize)) call messages_input_error(sys%namespace, 'GOObjective')
      call messages_print_var_option("GOObjective", g_opt%what2minimize, namespace=sys%namespace)
 
 
      !%Variable XYZGOConstrains
      !%Type string
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% <tt>Octopus</tt> will try to read the coordinate-dependent constrains from the XYZ file
      !% specified by the variable <tt>XYZGOConstrains</tt>.
      !% Note: It is important for the contrains to maintain the ordering
      !% in which the atoms were defined in the coordinates specifications.
      !% Moreover, constrains impose fixed absolute coordinates, therefore
      !% constrains are not compatible with GOCenter = yes
      !%End
 
      !%Variable XSFGOConstrains
      !%Type string
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% Like <tt>XYZGOConstrains</tt> but in XCrySDen format, as in <tt>XSFCoordinates</tt>.
      !%End
 
      !%Variable PDBGOConstrains
      !%Type string
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% Like <tt>XYZGOConstrains</tt> but in PDB format, as in <tt>PDBCoordinates</tt>.
      !%End
 
      !%Variable GOConstrains
      !%Type block
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% If <tt>XYZGOConstrains</tt>, <tt>PDBConstrains</tt>, and <tt>XSFGOConstrains</tt>
      !% are not present, <tt>Octopus</tt> will try to fetch the geometry optimization
      !% contrains from this block. If this block is not present, <tt>Octopus</tt>
      !% will not set any constrains. The format of this block can be
      !% illustrated by this example:
      !%
      !% <tt>%GOConstrains
      !% <br>&nbsp;&nbsp;'C'  |      1 | 0 | 0
      !% <br>&nbsp;&nbsp;'O'  | &nbsp;1 | 0 | 0
      !% <br>%</tt>
      !%
      !% Coordinates with a constrain value of 0 will be optimized, while
      !% coordinates with a constrain different from zero will be kept fixed. So,
      !% in this example the x coordinates of both atoms will remain fixed and the
      !% distance between the two atoms along the x axis will be constant.
      !%
      !% Note: It is important for the constrains to maintain the ordering
      !% in which the atoms were defined in the coordinates specifications.
      !% Moreover, constrains impose fixed absolute coordinates, therefore
      !% constrains are not compatible with GOCenter = yes
      !%End
 
      call read_coords_init(xyz)
      call read_coords_read('GOConstrains', xyz, g_opt%ions%space, sys%namespace)
      if (xyz%source /= read_coords_err) then
        !Sanity check
        if (g_opt%ions%natoms /= xyz%n) then
          write(message(1), '(a,i4,a,i4)') 'I need exactly ', g_opt%ions%natoms, ' constrains, but I found ', xyz%n
          call messages_fatal(1, namespace=sys%namespace)
        end if
        ! copy information and adjust units
        do iatom = 1, g_opt%ions%natoms
          where(abs(xyz%atom(iatom)%x) <= m_epsilon)
            g_opt%ions%atom(iatom)%c = m_zero
          elsewhere
            g_opt%ions%atom(iatom)%c = m_one
          end where
        end do
 
        call read_coords_end(xyz)
 
 
        if (g_opt%fixed_atom > 0) then
          call messages_not_implemented("GOCenter with constrains", namespace=sys%namespace)
        end if
      else
        do iatom = 1, g_opt%ions%natoms
          g_opt%ions%atom(iatom)%c = m_zero
        end do
      end if
 
      !%Variable GOCalculator
      !%Type integer
      !%Default internal
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% Defines what calculator is used for performing geometry optimization.
      !% By default, Octopus is used for computing the forces and energy.
      !% However, it is possible to use (if selected at compile time) a machine-learned
      !% universal force field to do the optimization of the geometry instead of Octopus.
      !%
      !%Option internal 1
      !% Usual case:  forces/stress are only computed by Octopus
      !%Option mattersim 2
      !% MatterSim deep learning atomistic model.
      !%End
      call parse_variable(sys%namespace, 'GOCalculator', go_calculator_internal, g_opt%calculator)
      if (.not. varinfo_valid_option('GOCalculator', g_opt%calculator)) then
        call messages_input_error(sys%namespace, 'GOCalculator')
      end if
#ifndef HAVE_MATTERSIM
      if (g_opt%calculator == go_calculator_mattersim) then
        message(1) = 'Octopus was not compiled with support for MatterSim'
        call messages_fatal(1, namespace=sys%namespace)
      end if
#endif
 
      !%Variable GOMLPreoptimize
      !%Type logical
      !%Default yes (if compiled with MatterSim)
      !%Section Calculation Modes::Geometry Optimization
      !%Description
      !% If true, a machine-learned force field is used to pre-optimize the geometry,
      !% once, before the actual optimization starts.
      !% Currently, only uses MatterSim. This is not compatible with <tt>GOCalculator=mattersim</tt>.
      !% Requires Octopus to be compiled with MatterSim; the default is no otherwise.
      !% The default is also no for <tt>GOCalculator=mattersim</tt>, and when relaxing
      !% the cell of a system which is not periodic along all three directions.
      !% The pre-optimization keeps the coordinates fixed by <tt>GOConstrains</tt> and
      !% <tt>GOCenter</tt>, as well as the atoms of a fixed species.
      !%End
#ifdef HAVE_MATTERSIM
      default_val = .true.
      if (g_opt%calculator == go_calculator_mattersim) default_val = .false.
      !TODO: Add support for mixed periodicity
      if (sys%space%periodic_dim /= 3 .and. &
        (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0)) default_val = .false.
#else
      default_val = .false.
#endif
      call parse_variable(sys%namespace, 'GOMLPreoptimize', default_val, g_opt%ml_preopt)
#ifndef HAVE_MATTERSIM
      if (g_opt%ml_preopt) then
        message(1) = 'Octopus was not compiled with support for MatterSim'
        call messages_fatal(1, namespace=sys%namespace)
      end if
#endif
 
      call io_mkdir('geom', sys%namespace)
 
      if (g_opt%ml_preopt) then
        call ml_preopt(g_opt)
      end if
 
 
      call io_rm("geom/optimization.log", sys%namespace)
 
      call io_rm("work-geom.xyz", sys%namespace)
 
      if (.not. fromscratch) then
        inquire(file = './last.xyz', exist = does_exist)
        if (.not. does_exist) fromscratch = .true.
      end if
 
      if (.not. fromscratch) call g_opt%ions%read_xyz('./last')
 
      ! clean out old geom/go.XXXX.xyz files. must be consistent with write_iter_info
      iter = 1
      do
        write(filename, '(a,i4.4,a)') "geom/go.", iter, ".xyz"
        inquire(file = trim(filename), exist = does_exist)
        if (does_exist) then
          call io_rm(trim(filename), sys%namespace)
          if (sys%space%periodic_dim > 0) then
            write(filename, '(a,i4.4,a)') "geom/go.", iter, ".xyz_red"
            call io_rm(trim(filename), sys%namespace)
            write(filename, '(a,i4.4,a)') "geom/go.", iter, ".xsf"
            call io_rm(trim(filename), sys%namespace)
          end if
          iter = iter + 1
        else
          exit
        end if
        ! TODO: clean forces directory
      end do
 
      call g_opt%scfv%restart_dump%init(sys%namespace, restart_gs, restart_type_dump, sys%mc, ierr, mesh=sys%gr)
 
      pop_sub(geom_opt_run_legacy.init_)
    end subroutine init_
 
 
    ! ---------------------------------------------------------
    subroutine end_()
      push_sub(geom_opt_run_legacy.end_)
 
      call states_elec_deallocate_wfns(sys%st)
 
      call g_opt%scfv%restart_dump%end()
 
      nullify(g_opt%mesh)
      nullify(g_opt%ions)
      nullify(g_opt%st)
      nullify(g_opt%hm)
      nullify(g_opt%syst)
 
      safe_deallocate_a(g_opt%cell_force)
 
      pop_sub(geom_opt_run_legacy.end_)
    end subroutine end_
 
  end subroutine geom_opt_run_legacy
 
 
  ! ---------------------------------------------------------
  subroutine calc_point(size, coords, objective, getgrad, df)
    integer,          intent(in)    :: size
    real(real64),     intent(in)    :: coords(size)
    real(real64),     intent(inout) :: objective
    integer,          intent(in)    :: getgrad
    real(real64),     intent(inout) :: df(size)
 
    integer :: iatom, idir, jdir, iunit, command_status, command_exit
    real(real64), dimension(g_opt%periodic_dim, g_opt%periodic_dim) :: stress, strain, right_stretch, inv_right_stretch, &
      rotation, sym_stress
    real(real64) :: det
    character(len=MAX_PATH_LEN) :: cmd, mattersim_share, mattersim_python
 
    push_sub(calc_point)
 
    assert(size == g_opt%size)
    command_status = 0
    command_exit = 0
 
    call from_coords(g_opt, coords)
 
    if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0 ) then
      call ion_dynamics_box_update(g_opt%syst%namespace, g_opt%syst%gr, g_opt%syst%space, g_opt%ions%latt)
    end if
 
    if (g_opt%fixed_atom /= 0) then
      call g_opt%ions%translate(g_opt%ions%center())
    end if
 
    ! When the system is periodic in some directions, the atoms might have moved to a an adjacent cell,
    ! so we need to move them back to the original cell
    call g_opt%ions%fold_atoms_into_cell()
 
    ! Some atoms might have moved outside the simulation box. We stop if this happens.
    do iatom = 1, g_opt%ions%natoms
      if (.not. g_opt%syst%gr%box%contains_point(g_opt%ions%pos(:, iatom))) then
        if (g_opt%syst%space%periodic_dim /= g_opt%syst%space%dim) then
          ! FIXME: This could fail for partial periodicity systems
          ! because contains_point is too strict with atoms close to
          ! the upper boundary to the cell.
          write(message(1), '(a,i5,a)') "Atom ", iatom, " has moved outside the box during the geometry optimization."
          call messages_fatal(1, namespace=g_opt%syst%namespace)
        end if
      end if
    end do
 
    if (g_opt%calculator==go_calculator_internal) then
      call g_opt%ions%write_xyz('./work-geom', append = .true.)
    else
      if (g_opt%ions%space%is_periodic()) then
        call write_extended_xyz_file('.', 'work-geom', g_opt%ions%space, g_opt%ions%latt, g_opt%ions%pos, g_opt%ions%atom, &
          g_opt%syst%gr%box, g_opt%syst%namespace)
      else
        call g_opt%ions%write_xyz('./work-geom', append = .false.)
      end if
    end if
 
    ! Update lattice vectors and regenerate grid
    if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0 ) then
      call electrons_lattice_vectors_update(g_opt%syst%namespace, g_opt%syst%gr, &
        g_opt%syst%space, g_opt%syst%hm%psolver, g_opt%syst%hm%kpoints, &
        g_opt%syst%mc, g_opt%syst%st%qtot, g_opt%ions%latt)
    end if
 
    ! This is already set to zero for periodic systems
    if (g_opt%ions%force_total_enforce) then
      call forces_set_total_to_zero(g_opt%ions, g_opt%ions%tot_force)
    end if
 
    call scf_print_mem_use(g_opt%syst%namespace)
 
    ! Select the calculator: Octopus or a ML force field
    select case (g_opt%calculator)
    case (go_calculator_internal)
 
      call scf_mix_clear(g_opt%scfv)
      call hamiltonian_elec_epot_generate(g_opt%hm, g_opt%syst%namespace, g_opt%syst%space, g_opt%syst%gr, &
        g_opt%ions, g_opt%syst%ext_partners, g_opt%st)
      call density_calc(g_opt%st, g_opt%syst%gr, g_opt%st%rho)
      call v_ks_calc(g_opt%syst%ks, g_opt%syst%namespace, g_opt%syst%space, g_opt%hm, g_opt%st, &
        g_opt%ions,g_opt%syst%ext_partners, calc_eigenval = .true.)
      call energy_calc_total(g_opt%syst%namespace, g_opt%syst%space, g_opt%hm, g_opt%syst%gr, g_opt%st, g_opt%syst%ext_partners)
 
      ! do SCF calculation
      call scf_run(g_opt%scfv, g_opt%syst%namespace, g_opt%syst%space, g_opt%syst%mc, g_opt%syst%gr, &
        g_opt%ions, g_opt%syst%ext_partners, &
        g_opt%st, g_opt%syst%ks, g_opt%hm, outp = g_opt%syst%outp, verbosity = verb_compact, restart_dump=g_opt%scfv%restart_dump)
 
      call scf_print_mem_use(g_opt%syst%namespace)
 
    case (go_calculator_mattersim)
      call mattersim_paths(mattersim_share, mattersim_python, g_opt%syst%namespace)
      if (mpi_world%is_root()) then
        ! Execute a system call for the current geometry
        ! -I isolates the venv python from any external PYTHONPATH/PYTHONHOME
        cmd = trim(mattersim_python) // " -I " // trim(mattersim_share) &
          // "/scripts/run_mattersim.py work-geom.xyz ml_forces.xyz "  // trim(mattersim_share)
        call execute_command_line(trim(cmd), exitstat=command_exit, cmdstat=command_status)
        if (command_status == 0) command_status = command_exit
      end if
      call mpi_world%bcast(command_status, 1, mpi_integer, 0)
      if (command_status /= 0) then
        write(message(1), '(a,i0)') 'MatterSim force calculation failed with exit status ', command_status
        call messages_fatal(1, namespace=g_opt%syst%namespace)
      end if
      call mpi_world%barrier()
 
      ! Read the result from MatterSim
      iunit = io_open('ml_forces.xyz', global_namespace, action='read')
      read(iunit, *) g_opt%hm%energy%total
      do iatom = 1, g_opt%ions%natoms
        read(iunit,*) g_opt%ions%tot_force(:, iatom)
      end do
      read(iunit, *) ! skip comment line
      if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0) then
        do idir = 1, 3
          read(iunit,*) g_opt%syst%st%stress_tensors%total(:, idir)
        end do
      end if
      call io_close(iunit)
 
    end select
 
    ! Convert stress into cell force
    ! This is the Parrinello-Rahman equation of motion of the cell,
    ! see Parrinello and Rahman, J. Appl. Pys. 52, 7182 (1981), Eq. 2.10.
    ! Here we use a slightly different definition, in which we evolve the right stretch tensor
    ! instead of h, because this is a symmetric matrix.
    ! The deformation tensor is defined as F = (1+\epsilon) = h h_0^{-1},
    ! where \epsilon is the infinitesimal strain tensor
    ! Using the polar decomposition F = R U, we can request to have \ddot{U} that remains symmetric
    ! in order to kill all the rotation
    if (bitand(g_opt%type, go_cell) /= 0) then
      stress = -g_opt%syst%st%stress_tensors%total(1:g_opt%periodic_dim, 1:g_opt%periodic_dim)
 
      ! Replace \sigma = sym(R^T \sigma R)
      ! We first compute R from the right polar decomposition R = F U^-1
      strain = matmul(g_opt%ions%latt%rlattice(1:g_opt%periodic_dim,1:g_opt%periodic_dim), g_opt%inv_initial_rlattice)
      right_stretch = lalg_remove_rotation(g_opt%periodic_dim, strain)
      inv_right_stretch = right_stretch
      call lalg_inverse(g_opt%periodic_dim, inv_right_stretch, 'dir')
      rotation = matmul(strain, inv_right_stretch)
      ! In some cases, R can contain a reflection. To guarantee that U is SPD, we multiply by the determinant of R
      det = lalg_determinant(g_opt%periodic_dim, rotation, preserve_mat=.true.)
 
      ! Compute the symmetrized stress
      sym_stress = matmul(transpose(rotation), matmul(stress, rotation))
      sym_stress = m_half*(sym_stress + transpose(sym_stress)) * det
 
      ! The external hydrostatic pressure enters as a full -P on each diagonal component of the
      ! stress, such that at equilibrium the internal stress balances the pressure (\sigma = -P I).
      do idir = 1, g_opt%periodic_dim
        sym_stress(idir, idir) = sym_stress(idir, idir) - g_opt%pressure
      end do
      g_opt%cell_force = sym_stress * g_opt%ions%latt%rcell_volume
 
      g_opt%cell_force = matmul(g_opt%cell_force, inv_right_stretch)
      ! The cell force is the gradient with respect to the symmetric right stretch tensor U,
      ! so we symmetrize it.
      g_opt%cell_force = m_half*(g_opt%cell_force + transpose(g_opt%cell_force))
    end if
 
    ! Convert stress into cell force
    if (bitand(g_opt%type, go_volume) /= 0) then
      stress = g_opt%syst%st%stress_tensors%total(1:g_opt%periodic_dim, 1:g_opt%periodic_dim)
      do idir = 1, g_opt%periodic_dim
        g_opt%cell_force(idir, 1) = -(g_opt%pressure + stress(idir, idir)) * g_opt%ions%latt%rcell_volume
      end do
    end if
 
    if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0) then
      write(message(1),'(a,3a,a)') '  Stress tensor [', trim(units_abbrev(units_out%length)), ']'
      do idir = 1, g_opt%periodic_dim
        write(message(1+idir),'(9e18.10)') (units_from_atomic(units_out%length, g_opt%syst%st%stress_tensors%total(jdir, idir)), &
          jdir = 1, g_opt%periodic_dim)
      end do
      call messages_info(1+g_opt%periodic_dim, namespace=g_opt%ions%namespace, debug_only=.true.)
      write(message(1),'(a,3a,a)') '  Cell force tensor [', trim(units_abbrev(units_out%length)), ']'
      do idir = 1, ubound(g_opt%cell_force, 2)
        write(message(1+idir),'(9e18.10)') (units_from_atomic(units_out%length, g_opt%cell_force(jdir, idir)), &
          jdir = 1, g_opt%periodic_dim)
      end do
      call messages_info(1+g_opt%periodic_dim, namespace=g_opt%ions%namespace, debug_only=.true.)
    end if
 
 
    ! store results
    if (getgrad == 1) call to_grad(g_opt, df)
 
    if (g_opt%what2minimize == minwhat_forces) then
      objective = m_zero
      do iatom = 1, g_opt%ions%natoms
        if (g_opt%ions%fixed(iatom)) cycle
        objective = objective + sum(g_opt%ions%tot_force(:, iatom)**2)
      end do
      if (bitand(g_opt%type, go_cell) /= 0) then
        do idir = 1, g_opt%periodic_dim
          objective = objective + sum(g_opt%cell_force(:, idir)**2)
        end do
      end if
      if (bitand(g_opt%type, go_volume) /= 0) then
        objective = objective + sum(g_opt%cell_force(:,1)**2)
      end if
      objective = sqrt(objective)
    else
      objective = g_opt%hm%energy%total
    end if
 
    pop_sub(calc_point)
  end subroutine calc_point
 
 
  ! ---------------------------------------------------------
  subroutine calc_point_ng(size, coords, objective)
    integer     :: size
    real(real64)    :: coords(size)
    real(real64)    :: objective
 
    integer :: getgrad
    real(real64), allocatable :: df(:)
 
    push_sub(calc_point_ng)
 
    assert(size == g_opt%size)
 
    getgrad = 0
    safe_allocate(df(1:size))
    df = m_zero
 
    call calc_point(size, coords, objective, getgrad, df)
    safe_deallocate_a(df)
 
    pop_sub(calc_point_ng)
  end subroutine calc_point_ng
 
 
  ! ---------------------------------------------------------
  subroutine write_iter_info(geom_iter, size, energy, maxdx, maxdf, coords)
    integer,     intent(in) :: geom_iter
    integer,     intent(in) :: size
    real(real64),    intent(in) :: energy, maxdx, maxdf
    real(real64),    intent(in) :: coords(size)
 
    character(len=256) :: c_geom_iter, title, c_forces_iter
    integer :: iunit
 
    push_sub(write_iter_info)
 
    write(c_geom_iter, '(a,i4.4)') "go.", geom_iter
    write(title, '(f16.10,2x,a)') units_from_atomic(units_out%energy, energy), trim(units_abbrev(units_out%energy))
    call g_opt%ions%write_xyz('geom/'//trim(c_geom_iter), comment = trim(title))
    call g_opt%ions%write_xyz('./last')
 
    if(g_opt%periodic_dim > 0) then
      call g_opt%ions%write_xyz('geom/'//trim(c_geom_iter), comment = 'Reduced coordinates', reduce_coordinates = .true.)
      call write_xsf_geometry_file('geom', trim(c_geom_iter), g_opt%ions%space, g_opt%ions%latt, &
        g_opt%ions%pos, g_opt%ions%atom, g_opt%syst%gr, g_opt%syst%namespace)
    end if
 
    if (g_opt%syst%outp%what(option__output__forces)) then
      write(c_forces_iter, '(a,i4.4)') "forces.", geom_iter
      if (bitand(g_opt%syst%outp%how(option__output__forces), option__outputformat__bild) /= 0) then
        call g_opt%ions%write_bild_forces_file('forces', trim(c_forces_iter))
      else
        call write_xsf_geometry_file('forces', trim(c_forces_iter), g_opt%ions%space, g_opt%ions%latt, &
          g_opt%ions%pos, g_opt%ions%atom, g_opt%syst%gr, g_opt%syst%namespace, total_forces=g_opt%ions%tot_force)
      end if
    end if
 
    if (g_opt%syst%st%system_grp%is_root()) then
      iunit = io_open(trim('geom/optimization.log'), g_opt%syst%namespace, &
        action = 'write', position = 'append')
 
      if (geom_iter == 1) then
        if (bitand(g_opt%type, go_cell) /= 0) then
          write(iunit, '(a10,5(5x,a20),a)') '#     iter','energy [' // trim(units_abbrev(units_out%energy)) // ']', &
            'max_force [' // trim(units_abbrev(units_out%force)) // ']',&
            ' max_dr [' // trim(units_abbrev(units_out%length)) // ']', &
            ' a, b, c ['// trim(units_abbrev(units_out%length)) // ']', &
            ' volume ['// trim(units_abbrev(units_out%length**3)) // ']',&
            ' alpha, beta, gamma [degrees]'
        else
          write(iunit, '(a10,3(5x,a20))') '#     iter','energy [' // trim(units_abbrev(units_out%energy)) // ']', &
            'max_force [' // trim(units_abbrev(units_out%force)) // ']',&
            ' max_dr [' // trim(units_abbrev(units_out%length)) // ']'
        end if
      end if
 
      if (bitand(g_opt%type, go_cell) /= 0) then
        write(iunit, '(i10,10f25.15)')  geom_iter, units_from_atomic(units_out%energy, energy), &
          units_from_atomic(units_out%force,maxdf), &
          units_from_atomic(units_out%length,maxdx), &
          units_from_atomic(units_out%length,norm2(g_opt%ions%latt%rlattice(1:3, 1))),&
          units_from_atomic(units_out%length,norm2(g_opt%ions%latt%rlattice(1:3, 2))),&
          units_from_atomic(units_out%length,norm2(g_opt%ions%latt%rlattice(1:3, 3))),&
          units_from_atomic(units_out%length**3, g_opt%ions%latt%rcell_volume), &
          g_opt%ions%latt%alpha, g_opt%ions%latt%beta, g_opt%ions%latt%gamma
      else
        write(iunit, '(i10,3f25.15)')  geom_iter, units_from_atomic(units_out%energy, energy), &
          units_from_atomic(units_out%force,maxdf), &
          units_from_atomic(units_out%length,maxdx)
      end if
 
      call io_close(iunit)
    end if
 
    call messages_new_line()
    call messages_new_line()
 
    call messages_write("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++", new_line = .true.)
 
    call messages_write("+++++++++++++++++++++ MINIMIZATION ITER #:")
    call messages_write(geom_iter, fmt = "I5")
    call messages_write(" ++++++++++++++++++++++", new_line = .true.)
 
    call messages_write("  Energy    = ")
    call messages_write(energy, units = units_out%energy, fmt = "es15.8,1x", print_units = .true., new_line = .true.)
 
    if (g_opt%periodic_dim == 0) then
      if (maxdf > m_zero) then
        call messages_write("  Max force = ")
        call messages_write(maxdf, units = units_out%force, fmt = "f18.10,1x", print_units = .true., new_line = .true.)
      end if
 
      call messages_write("  Max dr    = ")
      call messages_write(maxdx, units = units_out%length, fmt = "f18.10,1x", print_units = .true., new_line = .true.)
    else
      if (maxdf > m_zero) then
        call messages_write("  Max reduced force = ")
        call messages_write(maxdf, fmt = "f16.10,1x", print_units = .false., new_line = .true.)
      end if
 
      call messages_write("  Max reduced dr    = ")
      call messages_write(maxdx, fmt = "f16.10,1x", print_units = .false., new_line = .true.)
    end if
 
    call messages_write("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++", new_line = .true.)
    call messages_write("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++", new_line = .true.)
    call messages_new_line()
    call messages_new_line()
    call messages_info()
 
    pop_sub(write_iter_info)
  end subroutine write_iter_info
 
  ! ---------------------------------------------------------
  subroutine to_coords(gopt, coords)
    type(geom_opt_t), intent(in)  :: gopt
    real(real64),         intent(out) :: coords(:)
 
    integer :: iatom, idir, jdir, icoord
    real(real64) :: tmp_pos(gopt%dim), strain(g_opt%periodic_dim,g_opt%periodic_dim)
 
    push_sub(to_coords)
 
    icoord = 1
    ! Ion dynamics
    if (bitand(g_opt%type, go_ions) /= 0) then
      do iatom = 1, gopt%ions%natoms
        if (gopt%fixed_atom == iatom) cycle
        if (gopt%ions%fixed(iatom)) cycle
        tmp_pos = gopt%ions%pos(1:gopt%dim, iatom)
        if (gopt%fixed_atom > 0) tmp_pos = tmp_pos - gopt%ions%pos(1:gopt%dim, gopt%fixed_atom)
        tmp_pos = gopt%ions%latt%cart_to_red(tmp_pos)
        do idir = 1, gopt%dim
          coords(icoord) = tmp_pos(idir)
          icoord = icoord + 1
        end do
      end do
    end if
 
    ! Cell dynamics
    if (bitand(g_opt%type, go_cell) /= 0) then
      ! We compute the change in the right stretch tensor U, defined as h = U h_0.
      ! Only the periodic block of the lattice is relaxed; the non-periodic (e.g. vacuum)
      ! direction is left untouched, which is essential for 2D systems.
      strain = matmul(gopt%ions%latt%rlattice(1:g_opt%periodic_dim, 1:g_opt%periodic_dim), g_opt%inv_initial_rlattice)
      do idir = 1, g_opt%periodic_dim
        do jdir = idir, g_opt%periodic_dim
          coords(icoord) = strain(idir, jdir)
          icoord = icoord + 1
        end do
      end do
    end if
 
    ! Volume dynamics
    if (bitand(g_opt%type, go_volume) /= 0) then
      do idir = 1, g_opt%periodic_dim
        coords(icoord) = norm2(gopt%ions%latt%rlattice(1:g_opt%periodic_dim, idir))/g_opt%initial_length(idir)
        icoord = icoord + 1
      end do
    end if
 
 
    pop_sub(to_coords)
  end subroutine to_coords
 
  ! ---------------------------------------------------------
  subroutine to_grad(gopt, grad)
    type(geom_opt_t), intent(in)  :: gopt
    real(real64),      intent(out) :: grad(:)
 
    integer :: iatom, idir, jdir, icoord
    real(real64) :: tmp_force(1:gopt%dim)
 
    push_sub(to_grad)
 
    icoord = 1
    ! Ion dynamics
    if (bitand(g_opt%type, go_ions) /= 0) then
      do iatom = 1, gopt%ions%natoms
        if (gopt%fixed_atom == iatom) cycle
        if (gopt%ions%fixed(iatom)) cycle
        do idir = 1, gopt%dim
          if (abs(gopt%ions%atom(iatom)%c(idir)) <= m_epsilon) then
            tmp_force(idir) = -gopt%ions%tot_force(idir, iatom)
          else
            tmp_force(idir) = m_zero
          end if
          if (gopt%fixed_atom > 0) then
            tmp_force(idir) = tmp_force(idir) + gopt%ions%tot_force(idir, gopt%fixed_atom)
          end if
        end do
        tmp_force = gopt%ions%latt%cart_to_red(tmp_force)
        do idir = 1, gopt%dim
          grad(icoord) = tmp_force(idir)
          icoord = icoord + 1
        end do
      end do
    end if
 
    ! Cell dynamics
    if (bitand(g_opt%type, go_cell) /= 0) then
      do idir = 1, g_opt%periodic_dim
        do jdir = idir, g_opt%periodic_dim
          grad(icoord) = -g_opt%cell_force(idir, jdir)
          icoord = icoord + 1
        end do
      end do
    end if
 
    ! Volume dynamics
    if (bitand(g_opt%type, go_volume) /= 0) then
      do idir = 1, g_opt%periodic_dim
        grad(icoord) = -g_opt%cell_force(idir, 1)
        icoord = icoord + 1
      end do
    end if
 
 
    pop_sub(to_grad)
  end subroutine to_grad
 
  ! ---------------------------------------------------------
  subroutine from_coords(gopt, coords)
    type(geom_opt_t), intent(inout) :: gopt
    real(real64),         intent(in)    :: coords(:)
 
    integer :: iatom, idir, jdir, icoord
    real(real64) :: tmp_pos(gopt%dim, gopt%ions%natoms), strain(g_opt%periodic_dim,g_opt%periodic_dim)
    real(real64) :: eigenvalues(g_opt%periodic_dim)
 
    push_sub(from_coords)
 
    ! Get the new reduced atomic coordinates
    tmp_pos = m_zero
    icoord = 1
    ! Ion dynamics
    if (bitand(g_opt%type, go_ions) /= 0) then
      do iatom = 1, gopt%ions%natoms
        if (gopt%fixed_atom == iatom) cycle
        if (gopt%ions%fixed(iatom)) cycle
        do idir = 1, gopt%dim
          tmp_pos(idir, iatom) = coords(icoord)
          icoord = icoord + 1
        end do
      end do
    else
      do iatom = 1, gopt%ions%natoms
        tmp_pos(:, iatom) = gopt%ions%latt%cart_to_red(gopt%ions%pos(:, iatom))
      end do
    end if
 
    ! Updating the lattice vectors
    if (bitand(g_opt%type, go_cell) /= 0) then
      do idir = 1, g_opt%periodic_dim
        do jdir = idir, g_opt%periodic_dim
          strain(idir, jdir) = coords(icoord)
          icoord = icoord + 1
        end do
      end do
      call upper_triangular_to_hermitian(g_opt%periodic_dim, strain)
 
      ! Get the new lattice vectors from the new right stretch tensor U
      ! The strain tensor is defined as A = U * A_0.
      ! Only the periodic block is updated; the non-periodic (e.g. vacuum) lattice
      ! direction is preserved, which is essential for 2D systems (otherwise the
      ! vacuum direction would be corrupted/collapse).
      gopt%ions%latt%rlattice(1:g_opt%periodic_dim, 1:g_opt%periodic_dim) = matmul(strain, g_opt%initial_rlattice)
 
      ! Sanity check: the right stretch tensor U must remain symmetric positive-definite.
      ! Note: lalg_eigensolve overwrites strain with the eigenvectors, but strain is no longer
      ! needed below, so this is safe.
      call lalg_eigensolve(g_opt%periodic_dim, strain, eigenvalues)
      assert(all(eigenvalues > m_zero))
    end if
 
    ! Updating the lattice vectors
    if (bitand(g_opt%type, go_volume) /= 0) then
      do idir = 1, g_opt%periodic_dim
        gopt%ions%latt%rlattice(1:g_opt%periodic_dim, idir) = coords(icoord) &
          * gopt%initial_rlattice(1:g_opt%periodic_dim, idir)
        icoord = icoord + 1
      end do
    end if
 
    ! Symmetrize and update the lattice vectors
    if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0) then
      call g_opt%syst%gr%symmetrizer%symmetrize_lattice_vectors(g_opt%periodic_dim, g_opt%initial_rlattice, &
        gopt%ions%latt%rlattice(1:g_opt%periodic_dim, 1:g_opt%periodic_dim), gopt%symmetrize)
      call gopt%ions%update_lattice_vectors(gopt%ions%latt, gopt%symmetrize)
    end if
 
    ! Ion dynamics
    if (bitand(g_opt%type, go_ions) /= 0) then
      ! To Cartesian coordinates
      do iatom = 1, gopt%ions%natoms
        if (gopt%fixed_atom == iatom) cycle
        if (gopt%ions%fixed(iatom)) cycle
        tmp_pos(:, iatom) = gopt%ions%latt%red_to_cart(tmp_pos(:, iatom))
        do idir = 1, gopt%dim
          if (abs(gopt%ions%atom(iatom)%c(idir)) <= m_epsilon) then
            gopt%ions%pos(idir, iatom) = tmp_pos(idir, iatom)
          end if
        end do
        if (gopt%fixed_atom > 0) then
          gopt%ions%pos(:, iatom) = gopt%ions%pos(:, iatom) + gopt%ions%pos(:, gopt%fixed_atom)
        end if
      end do
    else
      do iatom = 1, gopt%ions%natoms
        gopt%ions%pos(:, iatom) = gopt%ions%latt%red_to_cart(tmp_pos(:, iatom))
      end do
    end if
 
    if (gopt%symmetrize) then
      call gopt%ions%symmetrize_atomic_coord()
    end if
 
    if (bitand(g_opt%type, go_cell) /= 0) then
      call gopt%ions%print_spacegroup()
    end if
 
    pop_sub(from_coords)
  end subroutine from_coords
 
  ! ---------------------------------------------------------
  subroutine write_iter_info_ng(geom_iter, size, energy, maxdx, coords)
    integer,     intent(in) :: geom_iter
    integer,     intent(in) :: size
    real(real64),    intent(in) :: energy, maxdx
    real(real64),    intent(in) :: coords(size)
 
    push_sub(write_iter_info_ng)
    call write_iter_info(geom_iter, size, energy, maxdx, -m_one, coords)
 
    pop_sub(write_iter_info_ng)
  end subroutine write_iter_info_ng
 
 
  ! ---------------------------------------------------------
  subroutine ml_preopt(g_opt)
    type(geom_opt_t), intent(inout) :: g_opt
 
    integer :: iunit, iatom, idir, command_status, command_exit
    logical :: has_constrains
    character(len=MAX_PATH_LEN) :: cmd, mattersim_share, mattersim_python
    character(len=64) :: constrains_arg
 
    push_sub(ml_preopt)
    command_status = 0
    command_exit = 0
 
    call io_rm('ml_cell.xyz', g_opt%syst%namespace)
 
    if (g_opt%ions%space%is_periodic()) then
      call write_extended_xyz_file('.', 'work-geom', g_opt%ions%space, g_opt%ions%latt, g_opt%ions%pos, g_opt%ions%atom, &
        g_opt%syst%gr%box, g_opt%syst%namespace)
    else
      call g_opt%ions%write_xyz('./work-geom', append = .false.)
    end if
 
    ! MatterSim relaxes every atom, so the fixed coordinates (GOConstrains, GOCenter
    ! and fixed species) have to be passed explicitly to the pre-optimizer.
    has_constrains = any(g_opt%ions%fixed) .or. g_opt%fixed_atom > 0
    do iatom = 1, g_opt%ions%natoms
      has_constrains = has_constrains .or. any(g_opt%ions%atom(iatom)%c /= 0)
    end do
 
    constrains_arg = ''
    if (has_constrains) then
      constrains_arg = ' --constrains ml_constrains.dat'
      if (mpi_world%is_root()) then
        iunit = io_open('ml_constrains.dat', global_namespace, action='write')
        do iatom = 1, g_opt%ions%natoms
          if (g_opt%ions%fixed(iatom) .or. iatom == g_opt%fixed_atom) then
            write(iunit, '(3(1x,i1))') (1, idir = 1, g_opt%ions%space%dim)
          else
            write(iunit, '(3(1x,i1))') (g_opt%ions%atom(iatom)%c(idir), idir = 1, g_opt%ions%space%dim)
          end if
        end do
        call io_close(iunit)
      end if
    end if
 
    call mattersim_paths(mattersim_share, mattersim_python, g_opt%syst%namespace)
    if (mpi_world%is_root()) then
      ! Execute a system call for the current geometry
      ! -I isolates the venv python from any external PYTHONPATH/PYTHONHOME
      write(cmd, '(a,1x,es24.16,1x,i0,1x,l1,a)') &
        trim(mattersim_python) // " -I " // trim(mattersim_share) &
        // "/scripts/run_mattersim_preopt.py work-geom.xyz ml_cell.xyz "  // trim(mattersim_share), &
        g_opt%pressure, g_opt%type, g_opt%symmetrize, trim(constrains_arg)
      call execute_command_line(trim(cmd), exitstat=command_exit, cmdstat=command_status)
      if (command_status == 0) command_status = command_exit
    end if
 
    call mpi_world%bcast(command_status, 1, mpi_integer, 0)
    if (command_status /= 0) then
      write(message(1), '(a,i0)') 'MatterSim pre-optimization failed with exit status ', command_status
      call messages_fatal(1, namespace=g_opt%syst%namespace)
    end if
 
    call mpi_world%barrier()
 
    ! Read the result from MatterSim
    iunit = io_open('ml_cell.xyz', global_namespace, action='read')
    do iatom = 1, g_opt%ions%natoms
      read(iunit,*) g_opt%ions%pos(:, iatom)
    end do
    read(iunit, *) ! skip comment line
    if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0) then
      do idir = 1, 3
        read(iunit,*) g_opt%ions%latt%rlattice(:, idir)
      end do
    end if
    call io_close(iunit)
 
    call g_opt%ions%update_lattice_vectors(g_opt%ions%latt, g_opt%symmetrize)
 
    if (debug%info) then
      call g_opt%ions%print_spacegroup()
    end if
 
    call g_opt%ions%fold_atoms_into_cell()
 
    ! Update lattice vectors and regenerate grid
    if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0 ) then
      call electrons_lattice_vectors_update(g_opt%syst%namespace, g_opt%syst%gr, &
        g_opt%syst%space, g_opt%syst%hm%psolver, g_opt%syst%hm%kpoints, &
        g_opt%syst%mc, g_opt%syst%st%qtot, g_opt%ions%latt)
    end if
 
    call hamiltonian_elec_epot_generate(g_opt%hm, g_opt%syst%namespace, g_opt%syst%space, g_opt%syst%gr, &
      g_opt%ions, g_opt%syst%ext_partners, g_opt%st)
 
    call g_opt%ions%print_spacegroup()
 
    pop_sub(ml_preopt)
  end subroutine ml_preopt
 
end module geom_opt_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: