minimizer.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 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 minimizer_oct_m
  use debug_oct_m
  use global_oct_m
#ifndef NDEBUG
  use, intrinsic :: ieee_exceptions
#endif
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use profiling_oct_m
  use messages_oct_m
 
  implicit none
 
  private
  public ::                    &
    loct_1dminimize,           &
    minimize_fire,             &
    minimize_multidim,         &
    minimize_multidim_nograd,  &
    minimize_multidim_nlopt
 
 
  integer, public, parameter ::      &
    MINMETHOD_STEEPEST_DESCENT =  1, &
    minmethod_fr_cg            =  2, &
    minmethod_pr_cg            =  3, &
    minmethod_bfgs             =  4, &
    minmethod_bfgs2            =  5, &
    minmethod_nmsimplex        =  6, &
    minmethod_sd_native        = -1, &
    minmethod_nlopt_bobyqa     =  7, &
    minmethod_fire             =  8, &
    minmethod_nlopt_lbfgs      =  9
 
  abstract interface
    subroutine minimizer_function_i(n, x, val)
      import real64
      implicit none
      integer :: n
      real(real64) :: x(n)
      real(real64) :: val
    end subroutine minimizer_function_i
    subroutine minimizer_with_grad_i(n, x, val, getgrad, grad)
      import real64
      implicit none
      integer, intent(in)    :: n
      real(real64), intent(in)    :: x(n)
      real(real64), intent(inout) :: val
      integer, intent(in)    :: getgrad
      real(real64), intent(inout) :: grad(n)
    end subroutine minimizer_with_grad_i
    subroutine info_i(iter, n, val, maxdr, maxgrad, x)
      import real64
      implicit none
      integer, intent(in) :: iter
      integer, intent(in) :: n
      real(real64), intent(in) :: val
      real(real64), intent(in) :: maxdr
      real(real64), intent(in) :: maxgrad
      real(real64), intent(in) :: x(n)
    end subroutine info_i
    subroutine info_no_grad_i(iter, n, val, maxdr, x)
      import real64
      implicit none
      integer, intent(in) :: iter
      integer, intent(in) :: n
      real(real64), intent(in) :: val
      real(real64), intent(in) :: maxdr
      real(real64), intent(in) :: x(n)
    end subroutine info_no_grad_i
  end interface
 
 
  interface loct_1dminimize
    subroutine oct_1dminimize(a, b, m, f, status)
      import real64
      implicit none
      real(real64), intent(inout) :: a, b, m
      interface
        subroutine f(x, fx)
          import real64
          implicit none
          real(real64), intent(in)  :: x
          real(real64), intent(out) :: fx
        end subroutine f
      end interface
      integer, intent(out) :: status
    end subroutine oct_1dminimize
  end interface loct_1dminimize
 
  interface loct_minimize
    integer function oct_minimize(method, dim, x, step, line_tol, &
      tolgrad, toldr, maxiter, f, write_iter_info, minimum)
      import minimizer_with_grad_i, info_i, real64
      implicit none
      integer, intent(in)              :: method
      integer, intent(in)              :: dim
      real(real64), intent(inout)      :: x
      real(real64), intent(in)         :: step
      real(real64), intent(in)         :: line_tol
      real(real64), intent(in)         :: tolgrad
      real(real64), intent(in)         :: toldr
      integer, intent(in)              :: maxiter
      procedure(minimizer_with_grad_i) :: f
      procedure(info_i)                :: write_iter_info
      real(real64), intent(out)        :: minimum
    end function oct_minimize
  end interface loct_minimize
 
  interface loct_minimize_direct
    function oct_minimize_direct(method, dim, x, step, toldr, maxiter, f, write_iter_info, minimum)
      import minimizer_function_i, info_no_grad_i, real64
      implicit none
      integer                         :: oct_minimize_direct
      integer, intent(in)             :: method
      integer, intent(in)             :: dim
      real(real64), intent(inout)     :: x
      real(real64), intent(in)        :: step
      real(real64), intent(in)        :: toldr
      integer, intent(in)             :: maxiter
      procedure(minimizer_function_i) :: f
      procedure(info_no_grad_i)       :: write_iter_info
      real(real64), intent(out)       :: minimum
    end function oct_minimize_direct
  end interface loct_minimize_direct
 
contains
 
  subroutine minimize_multidim_nograd(method, dim, x, step, toldr, maxiter, f, write_iter_info, minimum, ierr)
    integer, intent(in)    :: method
    integer, intent(in)    :: dim
    real(real64), intent(inout) :: x(:)
    real(real64), intent(in)    :: step
    real(real64), intent(in)    :: toldr
    integer, intent(in)         :: maxiter
    procedure(minimizer_function_i) :: f
    procedure(info_no_grad_i)   :: write_iter_info
    real(real64), intent(out)   :: minimum
    integer, intent(out)   :: ierr
 
    push_sub(minimize_multidim_nograd)
 
    assert(ubound(x, dim = 1) >= dim)
 
    ierr = 0
    select case (method)
    case (minmethod_nmsimplex)
      ierr = loct_minimize_direct(method, dim, x(1), step, toldr, maxiter, f, write_iter_info, minimum)
    end select
 
    pop_sub(minimize_multidim_nograd)
 
  end subroutine minimize_multidim_nograd
 
 
  subroutine minimize_multidim_nlopt(ierr, method, dim, x, step, toldr, maxiter, f, minimum, lb, ub)
    integer,                  intent(out)   :: ierr
    integer,                  intent(in)    :: method
    integer,                  intent(in)    :: dim
    real(real64), contiguous, intent(inout) :: x(:)
    real(real64),             intent(in)    :: step
    real(real64),             intent(in)    :: toldr
    integer,                  intent(in)    :: maxiter
    interface
      subroutine f(val, n, x, grad, need_gradient, f_data)
        use iso_c_binding
        real(c_double), intent(out) :: val
        integer(c_int), intent(in)  :: n
        real(c_double), intent(in)  :: x(*)
        real(c_double), intent(out) :: grad(*)
        integer(c_int), intent(in)  :: need_gradient
        type(c_ptr),    intent(in)  :: f_data
      end subroutine f
    end interface
    real(real64), intent(out)   :: minimum
    real(real64), intent(in), optional :: lb(:), ub(:)
#if defined(HAVE_NLOPT)
 
    interface
      subroutine nlo_create(opt, alg, n)
        use iso_c_binding
        type(c_ptr),    intent(out) :: opt
        integer(c_int), intent(in)  :: alg
        integer(c_int), intent(in)  :: n
      end subroutine nlo_create
 
      subroutine nlo_set_lower_bounds(ret, opt, lower_bounds)
        use iso_c_binding
        integer(c_int), intent(out)   :: ret
        type(c_ptr),    intent(inout) :: opt
        real(c_double), intent(in)    :: lower_bounds(*)
      end subroutine nlo_set_lower_bounds
 
      subroutine nlo_set_upper_bounds(ret, opt, upper_bounds)
        use iso_c_binding
        integer(c_int), intent(out)   :: ret
        type(c_ptr),    intent(inout) :: opt
        real(c_double), intent(in)    :: upper_bounds(*)
      end subroutine nlo_set_upper_bounds
 
      subroutine nlo_set_min_objective(ret, opt, f, f_data)
        use iso_c_binding
        integer(c_int), intent(out)   :: ret
        type(c_ptr),    intent(inout) :: opt
        interface
          subroutine f(val, n, x, grad, need_gradient, f_data)
            use iso_c_binding
            real(c_double), intent(out) :: val
            integer(c_int), intent(in)  :: n
            real(c_double), intent(in)  :: x(*)
            real(c_double), intent(out) :: grad(*)
            integer(c_int), intent(in)  :: need_gradient
            type(c_ptr),    intent(in)  :: f_data
          end subroutine f
        end interface
        type(c_ptr),    intent(in)    :: f_data
      end subroutine nlo_set_min_objective
 
      subroutine nlo_set_xtol_abs1(ret, opt, xtol_abs)
        use iso_c_binding
        integer(c_int), intent(out)   :: ret
        type(c_ptr),    intent(inout) :: opt
        real(c_double), intent(in)    :: xtol_abs
      end subroutine nlo_set_xtol_abs1
 
      subroutine nlo_set_initial_step1(ret, opt, initial_step1)
        use iso_c_binding
        integer(c_int), intent(out)   :: ret
        type(c_ptr),    intent(inout) :: opt
        real(c_double), intent(in)    :: initial_step1
      end subroutine nlo_set_initial_step1
 
      subroutine nlo_set_maxeval(ret, opt, maxeval)
        use iso_c_binding
        integer(c_int), intent(out)   :: ret
        type(c_ptr),    intent(inout) :: opt
        integer(c_int), intent(in)    :: maxeval
      end subroutine nlo_set_maxeval
 
      subroutine nlo_optimize(ret, opt, x, optf)
        use iso_c_binding
        integer(c_int), intent(out)   :: ret
        type(c_ptr),    intent(inout) :: opt
        real(c_double), intent(inout) :: x(*)
        real(c_double), intent(out)   :: optf
      end subroutine nlo_optimize
 
      subroutine nlo_destroy(opt)
        use iso_c_binding
        type(c_ptr), intent(inout) :: opt
      end subroutine nlo_destroy
 
      integer(c_int) function nlopt_algorithm_from_string(name) bind(c)
        use iso_c_binding
        character(kind=c_char), intent(in) :: name(*)
      end function nlopt_algorithm_from_string
    end interface
 
    type(c_ptr) :: opt
    integer :: ires
    integer(c_int) :: algorithm
#ifndef NDEBUG
    logical :: halting_mode(size(ieee_all))
#endif
 
#ifndef NDEBUG
    ! NLopt raises benign floating-point exceptions when handling unbounded
    ! optimizations (BOBYQA), so suspend trapping for builds that enable it
    call ieee_get_halting_mode(ieee_all, halting_mode)
    call ieee_set_halting_mode(ieee_all, .false.)
#endif
 
    ! Look up the algorithm by name: the numerical values of the nlopt_algorithm
    ! enum are not stable across NLopt versions (they changed in NLopt 2.9.0).
    select case (method)
    case (minmethod_nlopt_bobyqa)
      algorithm = nlopt_algorithm_from_string("LN_BOBYQA"//c_null_char)
    case (minmethod_nlopt_lbfgs)
      algorithm = nlopt_algorithm_from_string("LD_LBFGS"//c_null_char)
    case default
      algorithm = -1_c_int
    end select
    if (algorithm < 0_c_int) then
      message(1) = "The requested algorithm is not available in the linked NLopt library."
      call messages_fatal(1)
    end if
 
    opt = c_null_ptr
    call nlo_create(opt, algorithm, dim)
 
    if (present(lb)) then
      call nlo_set_lower_bounds(ires, opt, lb)
    end if
    if (present(ub)) then
      call nlo_set_upper_bounds(ires, opt, ub)
    end if
 
    call nlo_set_min_objective(ires, opt, f, c_null_ptr)
    ! This would set an inequality constraint (TODO)
    ! call nlo_add_inequality_constraint(ires, opt, myconstraint, d1, 1.0e-8_real64)
 
    call nlo_set_xtol_abs1(ires, opt, toldr)
    call nlo_set_initial_step1(ires, opt, step)
    call nlo_set_maxeval(ires, opt, maxiter)
 
    call nlo_optimize(ires, opt, x, minimum)
    ierr = ires
    call nlo_destroy(opt)
#ifndef NDEBUG
    call ieee_set_halting_mode(ieee_all, halting_mode)
#endif
#else
    ierr = 0
    minimum = -m_huge
#endif
  end subroutine minimize_multidim_nlopt
 
 
  !----------------------------------------------
  subroutine minimize_multidim(method, dim, x, step, line_tol, tolgrad, toldr, maxiter, f, write_iter_info, minimum, ierr)
    integer, intent(in)         :: method
    integer, intent(in)         :: dim
    real(real64), intent(inout) :: x(:)
    real(real64), intent(in)    :: step
    real(real64), intent(in)    :: line_tol
    real(real64), intent(in)    :: tolgrad
    real(real64), intent(in)    :: toldr
    integer, intent(in)         :: maxiter
    procedure(minimizer_with_grad_i) :: f
    procedure(info_i)           :: write_iter_info
    real(real64), intent(out)   :: minimum
    integer, intent(out)        :: ierr
 
    push_sub(minimize_multidim)
 
    assert(ubound(x, dim = 1) >= dim)
 
    select case (method)
    case (minmethod_sd_native)
      call minimize_sd(dim, x, step, maxiter, f, write_iter_info, minimum, ierr)
 
    case default
      ierr = loct_minimize(method, dim, x(1), step, line_tol, tolgrad, toldr, maxiter, f, write_iter_info, minimum)
 
    end select
 
    pop_sub(minimize_multidim)
 
  end subroutine minimize_multidim
 
  !----------------------------------------------
 
  subroutine minimize_sd(dim, x, step, maxiter, f, write_iter_info, minimum, ierr)
    integer, intent(in)         :: dim
    real(real64), intent(inout) :: x(:)
    real(real64), intent(in)    :: step
    integer, intent(in)         :: maxiter
    procedure(minimizer_with_grad_i) :: f
    procedure(info_i)           :: write_iter_info
    real(real64), intent(out)   :: minimum
    integer,      intent(out)   :: ierr
 
    integer :: iter
    real(real64), allocatable :: grad(:)
    real(real64) :: step2, maxgrad
 
    push_sub(minimize_sd)
 
    safe_allocate(grad(1:dim))
 
    step2 = step*10.0_real64
    do iter = 1, maxiter
      call f(dim, x, minimum, 1, grad)
 
      maxgrad = maxval(abs(grad))
 
      call write_iter_info(iter, dim, minimum, maxgrad*step2, maxgrad, x)
 
      x(1:dim) = x(1:dim) - step2*grad(1:dim)
 
      step2 = step2*0.99_real64
    end do
    ierr = 0
 
    pop_sub(minimize_sd)
  end subroutine minimize_sd
 
  !----------------------------------------------
  subroutine minimize_fire(dim, space_dim, x, step, tolgrad, maxiter, gradf, write_iter_info, en, ierr, mass, integrator)
    integer,      intent(in)    :: dim
    integer,      intent(in)    :: space_dim
    real(real64), intent(inout) :: x(:)
    real(real64), intent(in)    :: step
    real(real64), intent(in)    :: tolgrad
    integer,      intent(in)    :: maxiter
    procedure(minimizer_with_grad_i) :: gradf
    procedure(info_i)           :: write_iter_info
    real(real64), intent(out)   :: en
    integer,      intent(out)   :: ierr
    real(real64), intent(in)    :: mass(:)
    integer,      intent(in)    :: integrator
 
    real(real64), allocatable :: grad(:)
 
    integer :: p_times, n_times, iter, ia
    real(real64) :: dt, alpha, p_value, dt_max, dt_min, mix
    real(real64), allocatable :: grad_atoms(:), vel(:), dr_i(:), x_old(:), dr_atoms(:)
 
    ! Parameters from the original paper
    integer, parameter :: n_min = 5, n_decrease_max = 500
    real(real64), parameter :: f_alpha = 0.99_real64
    real(real64), parameter :: f_inc = 1.1_real64
    real(real64), parameter :: f_dec = 0.5_real64
    real(real64), parameter :: alpha_start = 0.25_real64
 
    real(real64), parameter :: tol_zero_forces = 1e-10_real64
 
    push_sub(minimize_fire)
 
    ! dim must be a multiple of space_dim
    assert(mod(dim,space_dim) == 0)
    assert(size(x) == dim)
 
    safe_allocate(grad_atoms(1:dim/space_dim))
    safe_allocate(grad(1:dim))
    safe_allocate(vel(1:dim))
    safe_allocate(dr_atoms(1:dim/space_dim))
    safe_allocate(x_old(1:dim))
    safe_allocate(dr_i(1:dim))
 
    ! Initial values
    ierr = 0
    vel = m_zero
    dt = step
    alpha = alpha_start
    dr_i = m_zero
    dt_max = 10.0_real64 * dt
    dt_min = 0.02_real64 * dt
    p_times = 0
    n_times = 0
    call gradf(dim, x, en, 1, grad)
    grad = -grad
 
    do iter = 1, maxiter
 
      ! Perform step F1: compute P = F.v
      p_value = dot_product(grad, vel)
 
      ! Perform step F3
      if (p_value > m_zero) then
        p_times = p_times + 1
        n_times = 0
        if (p_times > n_min) then
          dt = min(dt * f_inc , dt_max)
          alpha = alpha * f_alpha
        end if
      else ! or step F4
        p_times = 0
        n_times = n_times + 1
 
        if (n_times > n_decrease_max) exit ! If we are stuck, we exit
 
        if (iter >= n_min) then ! No decrease for the first few steps
          dt = max(dt * f_dec, dt_min)
          alpha = alpha_start
        end if
        x = x - m_half * vel * dt  ! Correct uphill motion
        vel = m_zero
      end if
 
      ! Store old values
      x_old = x
 
      ! Perform the MD step: get new gradients from the new positions
      ! Step F2: v-> (1-\alpha)v + \alpha \hat{F}.|v| is done inside the MD step
      ! We follow the integration schemes of the FIRE 2.0 paper
 
      ! Prevent division by zero in the mixing
      mix = lalg_nrm2(dim,grad)
      if (mix > tol_zero_forces) then
        mix = m_one / mix
      else
        mix = m_zero
      end if
 
      select case (integrator)
      case (option__gofireintegrator__verlet) ! Velocity Verlet, see Algorithm 6
        ! Velocity Verlet - update velocities 1
        vel = vel + m_half * grad * dt / mass
 
        ! Mixing
        if (p_value > m_zero) then
          vel = (m_one - alpha) * vel + alpha * mix * grad * lalg_nrm2(dim, vel)
        end if
 
        ! Get the new position
        dr_i = vel * dt
        call prevent_large_changes(dim, space_dim, dr_i, dr_atoms)
        x = x_old + dr_i
 
        ! Get the new gradient
        call gradf(dim, x, en, 1, grad)
        grad = -grad
 
        ! Velocity Verlet - update velocities 2
        vel = vel + m_half * grad * dt / mass
 
      case (option__gofireintegrator__euler) ! Explicit Euler method, see Algorithm 3 in Fire2.0 paper
        ! Mixing
        if (p_value > m_zero) then
          vel = (m_one - alpha) * vel + alpha * mix * grad * lalg_nrm2(dim, vel)
        end if
 
        ! Get the new position
        dr_i = vel * dt
        call prevent_large_changes(dim, space_dim, dr_i, dr_atoms)
        x = x_old + dr_i
 
        ! Get the new gradient
        call gradf(dim, x, en, 1, grad)
        grad = -grad
 
        ! Velocity Verlet - update velocities
        vel = vel + grad * dt / mass
 
 
      case (option__gofireintegrator__semi_implicit_euler) ! Semi-implicit Euler method, see Algorithm 4
        ! Velocity Verlet - update velocities
        vel = vel + grad * dt / mass
 
        ! Mixing
        if (p_value > m_zero) then
          vel = (m_one - alpha) * vel + alpha * mix * grad * lalg_nrm2(dim, vel)
        end if
 
        ! Get the new position
        dr_i = vel * dt
        call prevent_large_changes(dim, space_dim, dr_i, dr_atoms)
        x = x_old + dr_i
 
        ! Get the new gradient
        call gradf(dim, x, en, 1, grad)
        grad = -grad
      end select
 
      ! Get the norms of the gradients for each atom
      do ia = 0, dim/space_dim - 1
        grad_atoms(ia+1) = norm2(grad(space_dim*ia+1:space_dim*ia+space_dim))
      end do
 
      ! Output for each iteration step
      call write_iter_info(iter, dim, en, maxval(dr_atoms), maxval(abs(grad_atoms)), x)
 
      ! Check convergence
      if (maxval(abs(grad_atoms(1:))) < tolgrad) exit
    end do
 
    ierr = -iter
 
    safe_deallocate_a(dr_atoms)
    safe_deallocate_a(x_old)
    safe_deallocate_a(dr_i)
    safe_deallocate_a(vel)
    safe_deallocate_a(grad)
    safe_deallocate_a(grad_atoms)
 
    pop_sub(minimize_fire)
  end subroutine minimize_fire
 
  ! TODO: we should instead change the velocities,
  ! see S. Echeverri Restrepo and P. Andric, Computational Materials Science 218 (2023) 111978
  subroutine prevent_large_changes(dim, space_dim, dr_i, dr_atoms)
    integer,      intent(in)    :: dim
    integer,      intent(in)    :: space_dim
    real(real64), intent(inout) :: dr_i(:), dr_atoms(:)
 
    integer :: i1, i2, ia
 
    real(real64), parameter :: maxmove = 0.2_real64 * p_ang
 
    ! Get the norms of the displaments for each atoms
    do ia = 0, dim/space_dim - 1
      i1 = space_dim*ia+1
      i2 = space_dim*ia+space_dim
      dr_atoms(ia+1) = norm2(dr_i(i1:i2))
      ! Rescale the displacement to avoid too large changes
      if (dr_atoms(ia+1) > maxmove) then
        dr_i(i1:i2) = maxmove * dr_i(i1:i2) / dr_atoms(ia+1)
        dr_atoms(ia+1) = maxmove
      end if
    end do
  end subroutine prevent_large_changes
 
end module minimizer_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: