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 accel_oct_m
  use allelectron_oct_m
  use debug_oct_m
  use density_oct_m
  use electrons_oct_m
  use energy_calc_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
    type(restart_t)                   :: restart_dump
 
    real(real64), allocatable         :: cell_force(:,:)
    logical                           :: symmetrize = .false.
    real(real64), allocatable         :: initial_length(:)
    real(real64), allocatable         :: initial_rlattice(:,:)
    real(real64), allocatable         :: inv_initial_rlattice(:,:)
  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
 
contains
 
  ! ---------------------------------------------------------
  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
 
    push_sub(geom_opt_run_legacy)
 
    if (sys%space%periodic_dim == 1 .or. sys%space%periodic_dim == 2) then
      message(1) = "Geometry optimization not allowed for systems periodic in 1D and 2D, "
      message(2) = "as in those cases the total energy is not correct."
      call messages_fatal(2, 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
 
    call init_(fromscratch)
 
    ! load wavefunctions
    if (.not. fromscratch) then
      call restart_init(restart_load, 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, ierr)
      end if
      call restart_end(restart_load)
      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 SCFCalculeStress = 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)
    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)
    end if
 
    g_opt%symmetrize = sys%kpoints%use_symmetries .or. sys%st%symmetrize_density
 
    !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, real(g_opt%step, real64) ,&
        real(g_opt%toldr, real64) , 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 + 3
      end do
 
      !TODO: add variable to use Euler integrator
      call minimize_fire(g_opt%size, g_opt%ions%space%dim, coords, real(g_opt%step, real64) , real(g_opt%tolgrad, real64) , &
        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, real(g_opt%step, real64) ,&
        real(g_opt%line_tol, real64) , real(g_opt%tolgrad, real64) , real(g_opt%toldr, real64) , 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" atribute, 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
      !%Option cell_shape 2
      !% Relax cell shape. This changes lattice vector lengths and directions.
      !%Option cell_volume 4
      !% Relax cell volume. Only allow for rescaling the lengths of lattice vectors.
      !%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,x,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,x,a)') trim(message(1)), '+'
        end if
        write(message(1),'(a,x,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,x,a)') trim(message(1)), '+'
        end if
        write(message(1),'(a,x,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
 
      if (bitand(g_opt%type, go_cell) /= 0 .or. bitand(g_opt%type, go_volume) /= 0) then
        if (accel_is_enabled()) then
          message(1) = "Cell dynamics not supported on GPUs."
          call messages_fatal(1, namespace=sys%namespace)
        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
 
      assert(g_opt%size > 0)
 
      !%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) 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)) then
          g_opt%size = g_opt%size  - g_opt%dim
        end if
      end do
 
      !%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 verlet 1
      !% The Velocity Verlet algorithm.
      !%Option euler 0
      !% The Euler method.
      !%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
 
      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)
          iter = iter + 1
        else
          exit
        end if
        ! TODO: clean forces directory
      end do
 
      call restart_init(g_opt%restart_dump, 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 restart_end(g_opt%restart_dump)
 
      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
    real(real64), dimension(g_opt%periodic_dim, g_opt%periodic_dim) :: stress, scaling
 
 
    push_sub(calc_point)
 
    assert(size == g_opt%size)
 
    call from_coords(g_opt, coords)
 
    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%syst%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
 
    call g_opt%ions%write_xyz('./work-geom', append = .true.)
 
    call scf_mix_clear(g_opt%scfv)
 
    ! 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 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%restart_dump)
 
    call scf_print_mem_use(g_opt%syst%namespace)
 
    ! 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 right stretch tensor is defined as U = (1+\epsilon) = h h_0^{-1},
    ! where \epsilon is the infinitesimal strain tensor
    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)
      g_opt%cell_force = -stress * g_opt%ions%latt%rcell_volume
 
      ! We employ here the expression 6.a of Wentzcovitch, PRB 44, 2358 (1991)
      scaling = matmul(g_opt%ions%latt%rlattice(1:g_opt%periodic_dim,1:g_opt%periodic_dim), g_opt%inv_initial_rlattice)
      call lalg_inverse(g_opt%periodic_dim, scaling, 'dir')
      g_opt%cell_force = matmul(g_opt%cell_force, transpose(scaling))
 
      ! Thanks to the polar decomposition, any tensor can be expressed as a rotation times a symmetric matric
      ! In order to suppress any rotation, we only keep the symmetric part
      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) = -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
      if (debug%info) 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)
 
        write(message(1),'(a,3a,a)') '  Cell force 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%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)
      end if
    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 io_mkdir('geom', g_opt%ions%namespace)
    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 (mpi_grp_is_root(mpi_world)) 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 = "f16.10,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 = "f16.10,1x", print_units = .true., new_line = .true.)
      end if
 
      call messages_write("  Max dr    = ")
      call messages_write(maxdx, units = units_out%length, fmt = "f16.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
      strain = matmul(gopt%ions%latt%rlattice, 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)
 
    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
      gopt%ions%latt%rlattice = matmul(strain, g_opt%initial_rlattice)
 
      call symmetrize_lattice_vectors(g_opt, g_opt%dim, gopt%ions%latt%rlattice)
 
      call gopt%ions%update_lattice_vectors(gopt%ions%latt, gopt%symmetrize)
    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
 
      if (gopt%symmetrize) then
        call symmetrize_lattice_vectors(g_opt, g_opt%dim, gopt%ions%latt%rlattice)
      end if
      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
 
    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 symmetrize_lattice_vectors(g_opt, size, rlattice)
    type(geom_opt_t), intent(inout) :: g_opt
    integer,          intent(in)    :: size
    real(real64),     intent(inout) :: rlattice(size,size)
 
    real(real64) :: strain(size,size)
    real(real64) :: tol = 1.0e-14_real64
 
    push_sub(symmetrize_lattice_vectors)
 
    ! Remove too small elements
    call dzero_small_elements_matrix(g_opt%periodic_dim, rlattice, tol)
 
    ! Compute strain as rlattice * initial_rlattice^{-1}
    strain = matmul(rlattice, g_opt%inv_initial_rlattice)
 
    ! Remove too small elements
    call dzero_small_elements_matrix(g_opt%periodic_dim, strain, tol)
 
    if (g_opt%symmetrize) then
      ! Symmetrize the strain tensor
      call dsymmetrize_tensor_cart(g_opt%syst%gr%symm, strain, use_non_symmorphic=.true.)
    else ! The tensor should be at least symmetric, as we forbid rotations
      strain = m_half * (strain + transpose(strain))
    end if
 
    ! Remove too small elements
    call dzero_small_elements_matrix(g_opt%periodic_dim, strain, tol)
 
    ! Get the symmetrized lattice vectors
    rlattice = matmul(strain, g_opt%initial_rlattice)
 
    ! Remove too small elements
    call dzero_small_elements_matrix(g_opt%periodic_dim, rlattice, tol)
 
    pop_sub(symmetrize_lattice_vectors)
  end subroutine symmetrize_lattice_vectors
 
end module geom_opt_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: