cg_utils.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Utilities for Geometry optimization using  Conjugate Gradients
!> \author Teodoro Laino [teo]
!>      10.2005
! **************************************************************************************************
MODULE cg_utils
   USE cp_external_control, ONLY: external_control
   USE dimer_types, ONLY: dimer_env_type
   USE dimer_utils, ONLY: dimer_thrs, &
                          rotate_dimer
   USE global_types, ONLY: global_environment_type
   USE gopt_f_types, ONLY: gopt_f_type
   USE gopt_param_types, ONLY: gopt_param_type
   USE input_constants, ONLY: default_cell_method_id, &
                              default_minimization_method_id, &
                              default_shellcore_method_id, &
                              default_ts_method_id, &
                              ls_2pnt, &
                              ls_fit, &
                              ls_gold
   USE kinds, ONLY: dp
   USE mathconstants, ONLY: pi
   USE memory_utilities, ONLY: reallocate
   USE message_passing, ONLY: mp_para_env_type
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
 
! **************************************************************************************************
!> \brief evaluate the potential energy and its gradients using an array
!>      with same dimension as the particle_set
!> \param gopt_env the geometry optimization environment
!> \param x the position where the function should be evaluated
!> \param f the function value
!> \param gradient the value of its gradient
!> \par History
!>      none
!> \author Teodoro Laino [tlaino] - University of Zurich - 01.2008
! **************************************************************************************************
INTERFACE
 
   SUBROUTINE cp_eval_at(gopt_env, x, f, gradient, master, &
                         final_evaluation, para_env)
 
      USE message_passing, ONLY: mp_para_env_type
      USE gopt_f_types, ONLY: gopt_f_type
      USE kinds, ONLY: dp
 
      TYPE(gopt_f_type), POINTER               :: gopt_env
      REAL(KIND=dp), DIMENSION(:), POINTER     :: x
      REAL(KIND=dp), INTENT(out), OPTIONAL     :: f
      REAL(KIND=dp), DIMENSION(:), OPTIONAL, &
         POINTER                                :: gradient
      INTEGER, INTENT(IN)                      :: master
      LOGICAL, INTENT(IN), OPTIONAL            :: final_evaluation
      TYPE(mp_para_env_type), POINTER          :: para_env
 
   END SUBROUTINE cp_eval_at
 
END INTERFACE
 
   PUBLIC :: cg_linmin, get_conjugate_direction
   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cg_utils'
 
CONTAINS
 
! **************************************************************************************************
!> \brief Main driver for line minimization routines for CG
!> \param gopt_env ...
!> \param xvec ...
!> \param xi ...
!> \param g ...
!> \param opt_energy ...
!> \param output_unit ...
!> \param gopt_param ...
!> \param globenv ...
!> \par History
!>      10.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   RECURSIVE SUBROUTINE cg_linmin(gopt_env, xvec, xi, g, opt_energy, output_unit, gopt_param, &
                                  globenv)
 
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(kind=dp), DIMENSION(:), POINTER               :: xvec, xi, g
      REAL(kind=dp), INTENT(INOUT)                       :: opt_energy
      INTEGER                                            :: output_unit
      TYPE(gopt_param_type), POINTER                     :: gopt_param
      TYPE(global_environment_type), POINTER             :: globenv
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'cg_linmin'
 
      INTEGER                                            :: handle
 
      CALL timeset(routinen, handle)
      gopt_env%do_line_search = .true.
      SELECT CASE (gopt_env%type_id)
      CASE (default_minimization_method_id)
         SELECT CASE (gopt_param%cg_ls%type_id)
         CASE (ls_2pnt)
            CALL linmin_2pnt(gopt_env, xvec, xi, g, opt_energy, gopt_param, use_only_grad=gopt_param%cg_ls%grad_only, &
                             output_unit=output_unit)
         CASE (ls_fit)
            CALL linmin_fit(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brack_limit, &
                            gopt_param%cg_ls%initial_step, output_unit, gopt_param, globenv)
         CASE (ls_gold)
            CALL linmin_gold(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brent_tol, &
                             gopt_param%cg_ls%brent_max_iter, gopt_param%cg_ls%brack_limit, &
                             gopt_param%cg_ls%initial_step, output_unit, globenv)
         CASE DEFAULT
            cpabort("Line Search type not yet implemented in CG.")
         END SELECT
      CASE (default_ts_method_id)
         SELECT CASE (gopt_param%cg_ls%type_id)
         CASE (ls_2pnt)
            IF (gopt_env%dimer_rotation) THEN
               CALL rotmin_2pnt(gopt_env, gopt_env%dimer_env, xvec, xi, opt_energy)
            ELSE
               CALL tslmin_2pnt(gopt_env, gopt_env%dimer_env, xvec, xi, opt_energy, gopt_param, &
                                output_unit)
            END IF
         CASE DEFAULT
            cpabort("Line Search type not yet implemented in CG for TS search.")
         END SELECT
      CASE (default_cell_method_id)
         SELECT CASE (gopt_param%cg_ls%type_id)
         CASE (ls_2pnt)
            CALL linmin_2pnt(gopt_env, xvec, xi, g, opt_energy, gopt_param, use_only_grad=.true., &
                             output_unit=output_unit)
         CASE (ls_gold)
            CALL linmin_gold(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brent_tol, &
                             gopt_param%cg_ls%brent_max_iter, gopt_param%cg_ls%brack_limit, &
                             gopt_param%cg_ls%initial_step, output_unit, globenv)
         CASE DEFAULT
            cpabort("Line Search type not yet implemented in CG for cell optimization.")
         END SELECT
      CASE (default_shellcore_method_id)
         SELECT CASE (gopt_param%cg_ls%type_id)
         CASE (ls_2pnt)
            CALL linmin_2pnt(gopt_env, xvec, xi, g, opt_energy, gopt_param, use_only_grad=.true., &
                             output_unit=output_unit)
         CASE DEFAULT
            cpabort("Line Search type not yet implemented in CG for shellcore optimization.")
         END SELECT
 
      END SELECT
      gopt_env%do_line_search = .false.
      CALL timestop(handle)
 
   END SUBROUTINE cg_linmin
 
! **************************************************************************************************
!> \brief Line search subroutine based on 2 points (using gradients and energies
!>        or only gradients)
!> \param gopt_env ...
!> \param x0 ...
!> \param ls_vec ...
!> \param g ...
!> \param opt_energy ...
!> \param gopt_param ...
!> \param use_only_grad ...
!> \param output_unit ...
!> \author Teodoro Laino - created [tlaino] - 03.2008
! **************************************************************************************************
   RECURSIVE SUBROUTINE linmin_2pnt(gopt_env, x0, ls_vec, g, opt_energy, gopt_param, use_only_grad, &
                                    output_unit)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(kind=dp), DIMENSION(:), POINTER               :: x0, ls_vec, g
      REAL(kind=dp), INTENT(INOUT)                       :: opt_energy
      TYPE(gopt_param_type), POINTER                     :: gopt_param
      LOGICAL, INTENT(IN), OPTIONAL                      :: use_only_grad
      INTEGER, INTENT(IN)                                :: output_unit
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'linmin_2pnt'
 
      INTEGER                                            :: handle
      LOGICAL                                            :: my_use_only_grad, &
                                                            save_consistent_energy_force
      REAL(kind=dp)                                      :: a, b, c, dx, dx_min, dx_min_save, &
                                                            dx_thrs, norm_grad1, norm_grad2, &
                                                            norm_ls_vec, opt_energy2, x_grad_zero
      REAL(kind=dp), DIMENSION(:), POINTER               :: gradient2, ls_norm
 
      CALL timeset(routinen, handle)
      norm_ls_vec = sqrt(dot_product(ls_vec, ls_vec))
      my_use_only_grad = .false.
      IF (PRESENT(use_only_grad)) my_use_only_grad = use_only_grad
      IF (norm_ls_vec /= 0.0_dp) THEN
         ALLOCATE (ls_norm(SIZE(ls_vec)))
         ALLOCATE (gradient2(SIZE(ls_vec)))
         ls_norm = ls_vec/norm_ls_vec
         dx = norm_ls_vec
         dx_thrs = gopt_param%cg_ls%max_step
 
         x0 = x0 + dx*ls_norm
         ![NB] don't need consistent energies and forces if using only gradient
         save_consistent_energy_force = gopt_env%require_consistent_energy_force
         gopt_env%require_consistent_energy_force = .NOT. my_use_only_grad
         CALL cp_eval_at(gopt_env, x0, opt_energy2, gradient2, master=gopt_env%force_env%para_env%mepos, &
                         final_evaluation=.false., para_env=gopt_env%force_env%para_env)
         gopt_env%require_consistent_energy_force = save_consistent_energy_force
 
         norm_grad1 = -dot_product(g, ls_norm)
         norm_grad2 = dot_product(gradient2, ls_norm)
         IF (my_use_only_grad) THEN
            ! a*x+b=y
            ! per x=0; b=norm_grad1
            b = norm_grad1
            ! per x=dx; a*dx+b=norm_grad2
            a = (norm_grad2 - b)/dx
            x_grad_zero = -b/a
            dx_min = x_grad_zero
         ELSE
            ! ax**2+b*x+c=y
            ! per x=0 ; c=opt_energy
            c = opt_energy
            ! per x=dx;          a*dx**2 + b*dx + c = opt_energy2
            ! per x=dx;        2*a*dx    + b        = norm_grad2
            !
            !                  - a*dx**2        + c = (opt_energy2-norm_grad2*dx)
            !                    a*dx**2            = c - (opt_energy2-norm_grad2*dx)
            a = (c - (opt_energy2 - norm_grad2*dx))/dx**2
            b = norm_grad2 - 2.0_dp*a*dx
            dx_min = 0.0_dp
            IF (a /= 0.0_dp) dx_min = -b/(2.0_dp*a)
            opt_energy = opt_energy2
         END IF
         dx_min_save = dx_min
         ! In case the solution is larger than the maximum threshold let's assume the maximum allowed
         ! step length
         IF (abs(dx_min) > dx_thrs) dx_min = sign(1.0_dp, dx_min)*dx_thrs
         x0 = x0 + (dx_min - dx)*ls_norm
 
         ! Print out LS info
         IF (output_unit > 0) THEN
            WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
            WRITE (unit=output_unit, fmt="(T2,A,T31,A,T78,A)") &
               "***", "2PNT LINE SEARCH INFO", "***"
            WRITE (unit=output_unit, fmt="(T2,A,T78,A)") "***", "***"
            WRITE (unit=output_unit, fmt="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
               "***", "DX (EVALUATED)=", dx, "DX (THRESHOLD)=", dx_thrs, "***"
            WRITE (unit=output_unit, fmt="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
               "***", "DX (FITTED   )=", dx_min_save, "DX (ACCEPTED )=", dx_min, "***"
            WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
         END IF
         DEALLOCATE (ls_norm)
         DEALLOCATE (gradient2)
      ELSE
         ! Do Nothing, since.. if the effective force is 0 means that we are already
         ! in the saddle point..
      END IF
      CALL timestop(handle)
   END SUBROUTINE linmin_2pnt
 
! **************************************************************************************************
!> \brief Translational minimization for the Dimer Method - 2pnt LS
!> \param gopt_env ...
!> \param dimer_env ...
!> \param x0 ...
!> \param tls_vec ...
!> \param opt_energy ...
!> \param gopt_param ...
!> \param output_unit ...
!> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
! **************************************************************************************************
   SUBROUTINE tslmin_2pnt(gopt_env, dimer_env, x0, tls_vec, opt_energy, gopt_param, output_unit)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      TYPE(dimer_env_type), POINTER                      :: dimer_env
      REAL(kind=dp), DIMENSION(:), POINTER               :: x0, tls_vec
      REAL(kind=dp), INTENT(INOUT)                       :: opt_energy
      TYPE(gopt_param_type), POINTER                     :: gopt_param
      INTEGER, INTENT(IN)                                :: output_unit
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'tslmin_2pnt'
 
      INTEGER                                            :: handle
      REAL(kind=dp)                                      :: dx, dx_min, dx_min_acc, dx_min_save, &
                                                            dx_thrs, norm_tls_vec, opt_energy2
      REAL(kind=dp), DIMENSION(:), POINTER               :: tls_norm
 
      CALL timeset(routinen, handle)
      norm_tls_vec = sqrt(dot_product(tls_vec, tls_vec))
      IF (norm_tls_vec /= 0.0_dp) THEN
         ALLOCATE (tls_norm(SIZE(tls_vec)))
 
         tls_norm = tls_vec/norm_tls_vec
         dimer_env%tsl%tls_vec => tls_norm
 
         dx = norm_tls_vec
         dx_thrs = gopt_param%cg_ls%max_step
         ! If curvature is positive let's make the largest step allowed
         IF (dimer_env%rot%curvature > 0) dx = dx_thrs
         x0 = x0 + dx*tls_norm
         CALL cp_eval_at(gopt_env, x0, opt_energy2, master=gopt_env%force_env%para_env%mepos, &
                         final_evaluation=.false., para_env=gopt_env%force_env%para_env)
         IF (dimer_env%rot%curvature > 0) THEN
            dx_min = 0.0_dp
            dx_min_save = dx
            dx_min_acc = dx
         ELSE
            ! First let's try to interpolate the minimum
            dx_min = -opt_energy/(opt_energy2 - opt_energy)*dx
            ! In case the solution is larger than the maximum threshold let's assume the maximum allowed
            ! step length
            dx_min_save = dx_min
            IF (abs(dx_min) > dx_thrs) dx_min = sign(1.0_dp, dx_min)*dx_thrs
            dx_min_acc = dx_min
            dx_min = dx_min - dx
         END IF
         x0 = x0 + dx_min*tls_norm
 
         ! Print out LS info
         IF (output_unit > 0) THEN
            WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
            WRITE (unit=output_unit, fmt="(T2,A,T24,A,T78,A)") &
               "***", "2PNT TRANSLATIONAL LINE SEARCH INFO", "***"
            WRITE (unit=output_unit, fmt="(T2,A,T78,A)") "***", "***"
            WRITE (unit=output_unit, fmt="(T2,A,3X,A,F12.6,T78,A)") &
               "***", "LOCAL CURVATURE =", dimer_env%rot%curvature, "***"
            WRITE (unit=output_unit, fmt="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
               "***", "DX (EVALUATED)=", dx, "DX (THRESHOLD)=", dx_thrs, "***"
            WRITE (unit=output_unit, fmt="(T2,A,3X,A,F12.6,T45,A,F12.6,T78,A)") &
               "***", "DX (FITTED   )=", dx_min_save, "DX (ACCEPTED )=", dx_min_acc, "***"
            WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
         END IF
 
         ! Here we compute the value of the energy in point zero..
         CALL cp_eval_at(gopt_env, x0, opt_energy, master=gopt_env%force_env%para_env%mepos, &
                         final_evaluation=.false., para_env=gopt_env%force_env%para_env)
 
         DEALLOCATE (tls_norm)
      ELSE
         ! Do Nothing, since.. if the effective force is 0 means that we are already
         ! in the saddle point..
      END IF
      CALL timestop(handle)
 
   END SUBROUTINE tslmin_2pnt
 
! **************************************************************************************************
!> \brief Rotational minimization for the Dimer Method - 2 pnt LS
!> \param gopt_env ...
!> \param dimer_env ...
!> \param x0 ...
!> \param theta ...
!> \param opt_energy ...
!> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
! **************************************************************************************************
   SUBROUTINE rotmin_2pnt(gopt_env, dimer_env, x0, theta, opt_energy)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      TYPE(dimer_env_type), POINTER                      :: dimer_env
      REAL(kind=dp), DIMENSION(:), POINTER               :: x0, theta
      REAL(kind=dp), INTENT(INOUT)                       :: opt_energy
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'rotmin_2pnt'
 
      INTEGER                                            :: handle
      REAL(kind=dp)                                      :: a0, a1, angle, b1, curvature0, &
                                                            curvature1, curvature2, dcdp, f
      REAL(kind=dp), DIMENSION(:), POINTER               :: work
 
      CALL timeset(routinen, handle)
      curvature0 = dimer_env%rot%curvature
      dcdp = dimer_env%rot%dCdp
      b1 = 0.5_dp*dcdp
      angle = -0.5_dp*atan(dcdp/(2.0_dp*abs(curvature0)))
      dimer_env%rot%angle1 = angle
      dimer_env%cg_rot%nvec_old = dimer_env%nvec
      IF (angle > dimer_env%rot%angle_tol) THEN
         ! Rotating the dimer of dtheta degrees
         CALL rotate_dimer(dimer_env%nvec, theta, angle)
         ! Re-compute energy, gradients and rotation vector for new R1
         CALL cp_eval_at(gopt_env, x0, f, master=gopt_env%force_env%para_env%mepos, &
                         final_evaluation=.false., para_env=gopt_env%force_env%para_env)
 
         curvature1 = dimer_env%rot%curvature
         a1 = (curvature0 - curvature1 + b1*sin(2.0_dp*angle))/(1.0_dp - cos(2.0_dp*angle))
         a0 = 2.0_dp*(curvature0 - a1)
         angle = 0.5_dp*atan(b1/a1)
         curvature2 = a0/2.0_dp + a1*cos(2.0_dp*angle) + b1*sin(2.0_dp*angle)
         IF (curvature2 > curvature0) THEN
            angle = angle + pi/2.0_dp
            curvature2 = a0/2.0_dp + a1*cos(2.0_dp*angle) + b1*sin(2.0_dp*angle)
         END IF
         dimer_env%rot%angle2 = angle
         dimer_env%rot%curvature = curvature2
         ! Rotating the dimer the optimized (in plane) vector position
         dimer_env%nvec = dimer_env%cg_rot%nvec_old
         CALL rotate_dimer(dimer_env%nvec, theta, angle)
 
         ! Evaluate (by interpolation) the norm of the rotational force in the
         ! minimum of the rotational search (this is for print-out only)
         ALLOCATE (work(SIZE(dimer_env%nvec)))
         work = dimer_env%rot%g1
         work = sin(dimer_env%rot%angle1 - dimer_env%rot%angle2)/sin(dimer_env%rot%angle1)*dimer_env%rot%g1 + &
                sin(dimer_env%rot%angle2)/sin(dimer_env%rot%angle1)*dimer_env%rot%g1p + &
                (1.0_dp - cos(dimer_env%rot%angle2) - sin(dimer_env%rot%angle2)*tan(dimer_env%rot%angle1/2.0_dp))* &
                dimer_env%rot%g0
         work = -2.0_dp*(work - dimer_env%rot%g0)
         work = work - dot_product(work, dimer_env%nvec)*dimer_env%nvec
         opt_energy = sqrt(dot_product(work, work))
         DEALLOCATE (work)
      END IF
      dimer_env%rot%angle2 = angle
      CALL timestop(handle)
 
   END SUBROUTINE rotmin_2pnt
 
! **************************************************************************************************
!> \brief Line Minimization - Fit
!> \param gopt_env ...
!> \param xvec ...
!> \param xi ...
!> \param opt_energy ...
!> \param brack_limit ...
!> \param step ...
!> \param output_unit ...
!> \param gopt_param ...
!> \param globenv ...
!> \par History
!>      10.2005 created [tlaino]
!> \author Teodoro Laino
!> \note
!>      Given as input the vector XVEC and XI, finds the scalar
!>      xmin that minimizes the energy XVEC+xmin*XI. Replace step
!>      with the optimal value. Enhanced Version
! **************************************************************************************************
   SUBROUTINE linmin_fit(gopt_env, xvec, xi, opt_energy, &
                         brack_limit, step, output_unit, gopt_param, globenv)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(kind=dp), DIMENSION(:), POINTER               :: xvec, xi
      REAL(kind=dp)                                      :: opt_energy, brack_limit, step
      INTEGER                                            :: output_unit
      TYPE(gopt_param_type), POINTER                     :: gopt_param
      TYPE(global_environment_type), POINTER             :: globenv
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'linmin_fit'
 
      INTEGER                                            :: handle, loc_iter, odim
      LOGICAL                                            :: should_stop
      REAL(kind=dp)                                      :: ax, bx, fprev, rms_dr, rms_force, scale, &
                                                            xmin, xx
      REAL(kind=dp), DIMENSION(:), POINTER               :: pcom, xicom
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: hist
 
      CALL timeset(routinen, handle)
 
      NULLIFY (pcom, xicom, hist)
      rms_dr = gopt_param%rms_dr
      rms_force = gopt_param%rms_force
      ALLOCATE (pcom(SIZE(xvec)))
      ALLOCATE (xicom(SIZE(xvec)))
 
      pcom = xvec
      xicom = xi
      xicom = xicom/sqrt(dot_product(xicom, xicom))
      step = step*0.8_dp ! target a little before the minimum for the first point
      ax = 0.0_dp
      xx = step
      CALL cg_mnbrak(gopt_env, ax, xx, bx, pcom, xicom, brack_limit, output_unit, &
                     histpoint=hist, globenv=globenv)
      !
      fprev = 0.0_dp
      opt_energy = minval(hist(:, 2))
      odim = SIZE(hist, 1)
      scale = 0.25_dp
      loc_iter = 0
      DO WHILE (abs(hist(odim, 3)) > rms_force*scale .OR. abs(hist(odim, 1) - hist(odim - 1, 1)) > scale*rms_dr)
         CALL external_control(should_stop, "LINFIT", globenv=globenv)
         IF (should_stop) EXIT
         !
         loc_iter = loc_iter + 1
         fprev = opt_energy
         xmin = findmin(hist(:, 1), hist(:, 2), hist(:, 3))
         CALL reallocate(hist, 1, odim + 1, 1, 3)
         hist(odim + 1, 1) = xmin
         hist(odim + 1, 3) = cg_deval1d(gopt_env, xmin, pcom, xicom, opt_energy)
         hist(odim + 1, 2) = opt_energy
         odim = SIZE(hist, 1)
      END DO
      !
      xicom = xmin*xicom
      step = xmin
      xvec = xvec + xicom
      DEALLOCATE (pcom)
      DEALLOCATE (xicom)
      DEALLOCATE (hist)
      IF (output_unit > 0) THEN
         WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
         WRITE (unit=output_unit, fmt="(T2,A,T22,A,I7,T78,A)") &
            "***", "FIT LS  - NUMBER OF ENERGY EVALUATIONS : ", loc_iter, "***"
         WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
      END IF
      CALL timestop(handle)
 
   END SUBROUTINE linmin_fit
 
! **************************************************************************************************
!> \brief Line Minimization routine - GOLD
!> \param gopt_env ...
!> \param xvec ...
!> \param xi ...
!> \param opt_energy ...
!> \param brent_tol ...
!> \param brent_max_iter ...
!> \param brack_limit ...
!> \param step ...
!> \param output_unit ...
!> \param globenv ...
!> \par History
!>      10.2005 created [tlaino]
!> \author Teodoro Laino
!> \note
!>      Given as input the vector XVEC and XI, finds the scalar
!>      xmin that minimizes the energy XVEC+xmin*XI. Replaces XMIN
!>      with the optimal value
! **************************************************************************************************
   SUBROUTINE linmin_gold(gopt_env, xvec, xi, opt_energy, brent_tol, brent_max_iter, &
                          brack_limit, step, output_unit, globenv)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(kind=dp), DIMENSION(:), POINTER               :: xvec, xi
      REAL(kind=dp)                                      :: opt_energy, brent_tol
      INTEGER                                            :: brent_max_iter
      REAL(kind=dp)                                      :: brack_limit, step
      INTEGER                                            :: output_unit
      TYPE(global_environment_type), POINTER             :: globenv
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'linmin_gold'
 
      INTEGER                                            :: handle
      REAL(kind=dp)                                      :: ax, bx, xmin, xx
      REAL(kind=dp), DIMENSION(:), POINTER               :: pcom, xicom
 
      CALL timeset(routinen, handle)
 
      NULLIFY (pcom, xicom)
      ALLOCATE (pcom(SIZE(xvec)))
      ALLOCATE (xicom(SIZE(xvec)))
 
      pcom = xvec
      xicom = xi
      xicom = xicom/sqrt(dot_product(xicom, xicom))
      step = step*0.8_dp ! target a little before the minimum for the first point
      ax = 0.0_dp
      xx = step
      CALL cg_mnbrak(gopt_env, ax, xx, bx, pcom, xicom, brack_limit, output_unit, &
                     globenv=globenv)
 
      opt_energy = cg_dbrent(gopt_env, ax, xx, bx, brent_tol, brent_max_iter, &
                             xmin, pcom, xicom, output_unit, globenv)
      xicom = xmin*xicom
      step = xmin
      xvec = xvec + xicom
      DEALLOCATE (pcom)
      DEALLOCATE (xicom)
      CALL timestop(handle)
   END SUBROUTINE linmin_gold
 
! **************************************************************************************************
!> \brief Routine for initially bracketing a minimum based on the golden search
!>      minimum
!> \param gopt_env ...
!> \param ax ...
!> \param bx ...
!> \param cx ...
!> \param pcom ...
!> \param xicom ...
!> \param brack_limit ...
!> \param output_unit ...
!> \param histpoint ...
!> \param globenv ...
!> \par History
!>      10.2005 created [tlaino]
!> \author Teodoro Laino
!> \note
!>      Given two distinct initial points ax and bx this routine searches
!>      in the downhill direction and returns new points ax, bx, cx that
!>      bracket the minimum of the function
! **************************************************************************************************
   SUBROUTINE cg_mnbrak(gopt_env, ax, bx, cx, pcom, xicom, brack_limit, output_unit, &
                        histpoint, globenv)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(kind=dp)                                      :: ax, bx, cx
      REAL(kind=dp), DIMENSION(:), POINTER               :: pcom, xicom
      REAL(kind=dp)                                      :: brack_limit
      INTEGER                                            :: output_unit
      REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: histpoint
      TYPE(global_environment_type), POINTER             :: globenv
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'cg_mnbrak'
 
      INTEGER                                            :: handle, loc_iter, odim
      LOGICAL                                            :: hist, should_stop
      REAL(kind=dp)                                      :: dum, fa, fb, fc, fu, gold, q, r, u, ulim
 
      CALL timeset(routinen, handle)
      hist = PRESENT(histpoint)
      IF (hist) THEN
         cpassert(.NOT. ASSOCIATED(histpoint))
         ALLOCATE (histpoint(3, 3))
      END IF
      gold = (1.0_dp + sqrt(5.0_dp))/2.0_dp
      IF (hist) THEN
         histpoint(1, 1) = ax
         histpoint(1, 3) = cg_deval1d(gopt_env, ax, pcom, xicom, fa)
         histpoint(1, 2) = fa
         histpoint(2, 1) = bx
         histpoint(2, 3) = cg_deval1d(gopt_env, bx, pcom, xicom, fb)
         histpoint(2, 2) = fb
      ELSE
         fa = cg_eval1d(gopt_env, ax, pcom, xicom)
         fb = cg_eval1d(gopt_env, bx, pcom, xicom)
      END IF
      IF (fb .GT. fa) THEN
         dum = ax
         ax = bx
         bx = dum
         dum = fb
         fb = fa
         fa = dum
      END IF
      cx = bx + gold*(bx - ax)
      IF (hist) THEN
         histpoint(3, 1) = cx
         histpoint(3, 3) = cg_deval1d(gopt_env, cx, pcom, xicom, fc)
         histpoint(3, 2) = fc
      ELSE
         fc = cg_eval1d(gopt_env, cx, pcom, xicom)
      END IF
      loc_iter = 3
      DO WHILE (fb .GE. fc)
         CALL external_control(should_stop, "MNBRACK", globenv=globenv)
         IF (should_stop) EXIT
         !
         r = (bx - ax)*(fb - fc)
         q = (bx - cx)*(fb - fa)
         u = bx - ((bx - cx)*q - (bx - ax)*r)/(2.0_dp*sign(max(abs(q - r), tiny(0.0_dp)), q - r))
         ulim = bx + brack_limit*(cx - bx)
         IF ((bx - u)*(u - cx) .GT. 0.0_dp) THEN
            IF (hist) THEN
               odim = SIZE(histpoint, 1)
               CALL reallocate(histpoint, 1, odim + 1, 1, 3)
               histpoint(odim + 1, 1) = u
               histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
               histpoint(odim + 1, 2) = fu
            ELSE
               fu = cg_eval1d(gopt_env, u, pcom, xicom)
            END IF
            loc_iter = loc_iter + 1
            IF (fu .LT. fc) THEN
               ax = bx
               fa = fb
               bx = u
               fb = fu
               EXIT
            ELSE IF (fu .GT. fb) THEN
               cx = u
               fc = fu
               EXIT
            END IF
            u = cx + gold*(cx - bx)
            IF (hist) THEN
               odim = SIZE(histpoint, 1)
               CALL reallocate(histpoint, 1, odim + 1, 1, 3)
               histpoint(odim + 1, 1) = u
               histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
               histpoint(odim + 1, 2) = fu
            ELSE
               fu = cg_eval1d(gopt_env, u, pcom, xicom)
            END IF
            loc_iter = loc_iter + 1
         ELSE IF ((cx - u)*(u - ulim) .GT. 0.) THEN
            IF (hist) THEN
               odim = SIZE(histpoint, 1)
               CALL reallocate(histpoint, 1, odim + 1, 1, 3)
               histpoint(odim + 1, 1) = u
               histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
               histpoint(odim + 1, 2) = fu
            ELSE
               fu = cg_eval1d(gopt_env, u, pcom, xicom)
            END IF
            loc_iter = loc_iter + 1
            IF (fu .LT. fc) THEN
               bx = cx
               cx = u
               u = cx + gold*(cx - bx)
               fb = fc
               fc = fu
               IF (hist) THEN
                  odim = SIZE(histpoint, 1)
                  CALL reallocate(histpoint, 1, odim + 1, 1, 3)
                  histpoint(odim + 1, 1) = u
                  histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
                  histpoint(odim + 1, 2) = fu
               ELSE
                  fu = cg_eval1d(gopt_env, u, pcom, xicom)
               END IF
               loc_iter = loc_iter + 1
            END IF
         ELSE IF ((u - ulim)*(ulim - cx) .GE. 0.) THEN
            u = ulim
            IF (hist) THEN
               odim = SIZE(histpoint, 1)
               CALL reallocate(histpoint, 1, odim + 1, 1, 3)
               histpoint(odim + 1, 1) = u
               histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
               histpoint(odim + 1, 2) = fu
            ELSE
               fu = cg_eval1d(gopt_env, u, pcom, xicom)
            END IF
            loc_iter = loc_iter + 1
         ELSE
            u = cx + gold*(cx - bx)
            IF (hist) THEN
               odim = SIZE(histpoint, 1)
               CALL reallocate(histpoint, 1, odim + 1, 1, 3)
               histpoint(odim + 1, 1) = u
               histpoint(odim + 1, 3) = cg_deval1d(gopt_env, u, pcom, xicom, fu)
               histpoint(odim + 1, 2) = fu
            ELSE
               fu = cg_eval1d(gopt_env, u, pcom, xicom)
            END IF
            loc_iter = loc_iter + 1
         END IF
         ax = bx
         bx = cx
         cx = u
         fa = fb
         fb = fc
         fc = fu
      END DO
      IF (output_unit > 0) THEN
         WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
         WRITE (unit=output_unit, fmt="(T2,A,T22,A,I7,T78,A)") &
            "***", "MNBRACK - NUMBER OF ENERGY EVALUATIONS : ", loc_iter, "***"
         WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
      END IF
      CALL timestop(handle)
   END SUBROUTINE cg_mnbrak
 
! **************************************************************************************************
!> \brief Routine implementing the Brent Method
!>      Brent,R.P. Algorithm for Minimization without Derivatives, Chapt.5
!>      1973
!>      Extension in the use of derivatives
!> \param gopt_env ...
!> \param ax ...
!> \param bx ...
!> \param cx ...
!> \param tol ...
!> \param itmax ...
!> \param xmin ...
!> \param pcom ...
!> \param xicom ...
!> \param output_unit ...
!> \param globenv ...
!> \return ...
!> \par History
!>      10.2005 created [tlaino]
!> \author Teodoro Laino
!> \note
!>      Given a bracketing  triplet of abscissas ax, bx, cx (such that bx
!>      is between ax and cx and energy of bx is less than energy of ax and cx),
!>      this routine isolates the minimum to a precision of about tol using
!>      Brent method. This routine implements the extension of the Brent Method
!>      using derivatives
! **************************************************************************************************
   FUNCTION cg_dbrent(gopt_env, ax, bx, cx, tol, itmax, xmin, pcom, xicom, output_unit, &
                      globenv) RESULT(dbrent)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(kind=dp)                                      :: ax, bx, cx, tol
      INTEGER                                            :: itmax
      REAL(kind=dp)                                      :: xmin
      REAL(kind=dp), DIMENSION(:), POINTER               :: pcom, xicom
      INTEGER                                            :: output_unit
      TYPE(global_environment_type), POINTER             :: globenv
      REAL(kind=dp)                                      :: dbrent
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'cg_dbrent'
      REAL(kind=dp), PARAMETER                           :: zeps = 1.0e-8_dp
 
      INTEGER                                            :: handle, iter, loc_iter
      LOGICAL                                            :: ok1, ok2, should_stop, skip0, skip1
      REAL(kind=dp)                                      :: a, b, d, d1, d2, du, dv, dw, dx, e, fu, &
                                                            fv, fw, fx, olde, tol1, tol2, u, u1, &
                                                            u2, v, w, x, xm
 
      CALL timeset(routinen, handle)
      a = min(ax, cx)
      b = max(ax, cx)
      v = bx; w = v; x = v
      e = 0.0_dp
      dx = cg_deval1d(gopt_env, x, pcom, xicom, fx)
      fv = fx
      fw = fx
      dv = dx
      dw = dx
      loc_iter = 1
      DO iter = 1, itmax
         CALL external_control(should_stop, "BRENT", globenv=globenv)
         IF (should_stop) EXIT
         !
         xm = 0.5_dp*(a + b)
         tol1 = tol*abs(x) + zeps
         tol2 = 2.0_dp*tol1
         skip0 = .false.
         skip1 = .false.
         IF (abs(x - xm) .LE. (tol2 - 0.5_dp*(b - a))) EXIT
         IF (abs(e) .GT. tol1) THEN
            d1 = 2.0_dp*(b - a)
            d2 = d1
            IF (dw .NE. dx) d1 = (w - x)*dx/(dx - dw)
            IF (dv .NE. dx) d2 = (v - x)*dx/(dx - dv)
            u1 = x + d1
            u2 = x + d2
            ok1 = ((a - u1)*(u1 - b) .GT. 0.0_dp) .AND. (dx*d1 .LE. 0.0_dp)
            ok2 = ((a - u2)*(u2 - b) .GT. 0.0_dp) .AND. (dx*d2 .LE. 0.0_dp)
            olde = e
            e = d
            IF (.NOT. (ok1 .OR. ok2)) THEN
               skip0 = .true.
            ELSE IF (ok1 .AND. ok2) THEN
               IF (abs(d1) .LT. abs(d2)) THEN
                  d = d1
               ELSE
                  d = d2
               END IF
            ELSE IF (ok1) THEN
               d = d1
            ELSE
               d = d2
            END IF
            IF (.NOT. skip0) THEN
               IF (abs(d) .GT. abs(0.5_dp*olde)) skip0 = .true.
               IF (.NOT. skip0) THEN
                  u = x + d
                  IF ((u - a) .LT. tol2 .OR. (b - u) .LT. tol2) d = sign(tol1, xm - x)
                  skip1 = .true.
               END IF
            END IF
         END IF
         IF (.NOT. skip1) THEN
            IF (dx .GE. 0.0_dp) THEN
               e = a - x
            ELSE
               e = b - x
            END IF
            d = 0.5_dp*e
         END IF
         IF (abs(d) .GE. tol1) THEN
            u = x + d
            du = cg_deval1d(gopt_env, u, pcom, xicom, fu)
            loc_iter = loc_iter + 1
         ELSE
            u = x + sign(tol1, d)
            du = cg_deval1d(gopt_env, u, pcom, xicom, fu)
            loc_iter = loc_iter + 1
            IF (fu .GT. fx) EXIT
         END IF
         IF (fu .LE. fx) THEN
            IF (u .GE. x) THEN
               a = x
            ELSE
               b = x
            END IF
            v = w; fv = fw; dv = dw; w = x
            fw = fx; dw = dx; x = u; fx = fu; dx = du
         ELSE
            IF (u .LT. x) THEN
               a = u
            ELSE
               b = u
            END IF
            IF (fu .LE. fw .OR. w .EQ. x) THEN
               v = w; fv = fw; dv = dw
               w = u; fw = fu; dw = du
            ELSE IF (fu .LE. fv .OR. v .EQ. x .OR. v .EQ. w) THEN
               v = u
               fv = fu
               dv = du
            END IF
         END IF
      END DO
      IF (output_unit > 0) THEN
         WRITE (unit=output_unit, fmt="(/,T2,A)") repeat("*", 79)
         WRITE (unit=output_unit, fmt="(T2,A,T22,A,I7,T78,A)") &
            "***", "BRENT   - NUMBER OF ENERGY EVALUATIONS : ", loc_iter, "***"
         IF (iter == itmax + 1) &
            WRITE (unit=output_unit, fmt="(T2,A,T22,A,T78,A)") &
            "***", "BRENT - NUMBER OF ITERATIONS EXCEEDED ", "***"
         WRITE (unit=output_unit, fmt="(T2,A)") repeat("*", 79)
      END IF
      cpassert(iter /= itmax + 1)
      xmin = x
      dbrent = fx
      CALL timestop(handle)
 
   END FUNCTION cg_dbrent
 
! **************************************************************************************************
!> \brief Evaluates energy in one dimensional space defined by the point
!>      pcom and with direction xicom, position x
!> \param gopt_env ...
!> \param x ...
!> \param pcom ...
!> \param xicom ...
!> \return ...
!> \par History
!>      10.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   FUNCTION cg_eval1d(gopt_env, x, pcom, xicom) RESULT(my_val)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(kind=dp)                                      :: x
      REAL(kind=dp), DIMENSION(:), POINTER               :: pcom, xicom
      REAL(kind=dp)                                      :: my_val
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'cg_eval1d'
 
      INTEGER                                            :: handle
      REAL(kind=dp), DIMENSION(:), POINTER               :: xvec
 
      CALL timeset(routinen, handle)
 
      ALLOCATE (xvec(SIZE(pcom)))
      xvec = pcom + x*xicom
      CALL cp_eval_at(gopt_env, xvec, my_val, master=gopt_env%force_env%para_env%mepos, &
                      final_evaluation=.false., para_env=gopt_env%force_env%para_env)
      DEALLOCATE (xvec)
 
      CALL timestop(handle)
 
   END FUNCTION cg_eval1d
 
! **************************************************************************************************
!> \brief Evaluates derivatives in one dimensional space defined by the point
!>      pcom and with direction xicom, position x
!> \param gopt_env ...
!> \param x ...
!> \param pcom ...
!> \param xicom ...
!> \param fval ...
!> \return ...
!> \par History
!>      10.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   FUNCTION cg_deval1d(gopt_env, x, pcom, xicom, fval) RESULT(my_val)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(kind=dp)                                      :: x
      REAL(kind=dp), DIMENSION(:), POINTER               :: pcom, xicom
      REAL(kind=dp)                                      :: fval, my_val
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'cg_deval1d'
 
      INTEGER                                            :: handle
      REAL(kind=dp)                                      :: energy
      REAL(kind=dp), DIMENSION(:), POINTER               :: grad, xvec
 
      CALL timeset(routinen, handle)
 
      ALLOCATE (xvec(SIZE(pcom)))
      ALLOCATE (grad(SIZE(pcom)))
      xvec = pcom + x*xicom
      CALL cp_eval_at(gopt_env, xvec, energy, grad, master=gopt_env%force_env%para_env%mepos, &
                      final_evaluation=.false., para_env=gopt_env%force_env%para_env)
      my_val = dot_product(grad, xicom)
      fval = energy
      DEALLOCATE (xvec)
      DEALLOCATE (grad)
      CALL timestop(handle)
 
   END FUNCTION cg_deval1d
 
! **************************************************************************************************
!> \brief Find the minimum of a parabolic function obtained with a least square fit
!> \param x ...
!> \param y ...
!> \param dy ...
!> \return ...
!> \par History
!>      10.2005 created [fawzi]
!> \author Fawzi Mohamed
! **************************************************************************************************
   FUNCTION findmin(x, y, dy) RESULT(res)
      REAL(kind=dp), DIMENSION(:)                        :: x, y, dy
      REAL(kind=dp)                                      :: res
 
      INTEGER                                            :: i, info, iwork(8*3), lwork, min_pos, np
      REAL(kind=dp)                                      :: diag(3), res1(3), res2(3), res3(3), &
                                                            spread, sum_x, sum_xx, tmpw(1), &
                                                            vt(3, 3)
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: work
      REAL(kind=dp), DIMENSION(2*SIZE(x), 3)             :: f
      REAL(kind=dp), DIMENSION(2*SIZE(x))                :: b, w
      REAL(kind=dp)                                      :: u(2*SIZE(x), 3)
 
      np = SIZE(x)
      cpassert(np > 1)
      sum_x = 0._dp
      sum_xx = 0._dp
      min_pos = 1
      DO i = 1, np
         sum_xx = sum_xx + x(i)**2
         sum_x = sum_x + x(i)
         IF (y(min_pos) > y(i)) min_pos = i
      END DO
      spread = sqrt(sum_xx/real(np, dp) - (sum_x/real(np, dp))**2)
      DO i = 1, np
         w(i) = exp(-(real(np - i, dp))**2/(real(2*9, dp)))
         w(i + np) = 2._dp*w(i)
      END DO
      DO i = 1, np
         f(i, 1) = w(i)
         f(i, 2) = x(i)*w(i)
         f(i, 3) = x(i)**2*w(i)
         f(i + np, 1) = 0
         f(i + np, 2) = w(i + np)
         f(i + np, 3) = 2*x(i)*w(i + np)
      END DO
      DO i = 1, np
         b(i) = y(i)*w(i)
         b(i + np) = dy(i)*w(i + np)
      END DO
      lwork = -1
      CALL dgesdd('S', SIZE(f, 1), SIZE(f, 2), f, SIZE(f, 1), diag, u, SIZE(u, 1), vt, SIZE(vt, 1), tmpw, lwork, &
                  iwork, info)
      lwork = ceiling(tmpw(1))
      ALLOCATE (work(lwork))
      CALL dgesdd('S', SIZE(f, 1), SIZE(f, 2), f, SIZE(f, 1), diag, u, SIZE(u, 1), vt, SIZE(vt, 1), work, lwork, &
                  iwork, info)
      DEALLOCATE (work)
      CALL dgemv('T', SIZE(u, 1), SIZE(u, 2), 1._dp, u, SIZE(u, 1), b, 1, 0._dp, res1, 1)
      DO i = 1, 3
         res2(i) = res1(i)/diag(i)
      END DO
      CALL dgemv('T', 3, 3, 1._dp, vt, SIZE(vt, 1), res2, 1, 0._dp, res3, 1)
      res = -0.5*res3(2)/res3(3)
   END FUNCTION findmin
 
! **************************************************************************************************
!> \brief Computes the Conjugate direction for the next search
!> \param gopt_env ...
!> \param Fletcher_Reeves ...
!> \param g contains the theta  of the previous step.. (norm 1.0 vector)
!> \param xi contains the -theta of the present step.. (norm 1.0 vector)
!> \param h contains the search direction of the previous step (must be orthogonal
!>            to nvec of the previous step (nvec_old))
!> \par   Info for DIMER method
!> \par History
!>      10.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE get_conjugate_direction(gopt_env, Fletcher_Reeves, g, xi, h)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      LOGICAL, INTENT(IN)                                :: fletcher_reeves
      REAL(kind=dp), DIMENSION(:), POINTER               :: g, xi, h
 
      CHARACTER(len=*), PARAMETER :: routinen = 'get_conjugate_direction'
 
      INTEGER                                            :: handle
      LOGICAL                                            :: check
      REAL(kind=dp)                                      :: dgg, gam, gg, norm, norm_h
      TYPE(dimer_env_type), POINTER                      :: dimer_env
 
      CALL timeset(routinen, handle)
      NULLIFY (dimer_env)
      IF (.NOT. gopt_env%dimer_rotation) THEN
         gg = dot_product(g, g)
         IF (fletcher_reeves) THEN
            dgg = dot_product(xi, xi)
         ELSE
            dgg = dot_product((xi + g), xi)
         END IF
         gam = dgg/gg
         g = h
         h = -xi + gam*h
      ELSE
         dimer_env => gopt_env%dimer_env
         check = abs(dot_product(g, g) - 1.0_dp) < max(1.0e-9_dp, dimer_thrs)
         cpassert(check)
 
         check = abs(dot_product(xi, xi) - 1.0_dp) < max(1.0e-9_dp, dimer_thrs)
         cpassert(check)
 
         check = abs(dot_product(h, dimer_env%cg_rot%nvec_old)) < max(1.0e-9_dp, dimer_thrs)
         cpassert(check)
         gg = dimer_env%cg_rot%norm_theta_old**2
         IF (fletcher_reeves) THEN
            dgg = dimer_env%cg_rot%norm_theta**2
         ELSE
            norm = dimer_env%cg_rot%norm_theta*dimer_env%cg_rot%norm_theta_old
            dgg = dimer_env%cg_rot%norm_theta**2 + dot_product(g, xi)*norm
         END IF
         ! Compute Theta** and store it in nvec_old
         CALL rotate_dimer(dimer_env%cg_rot%nvec_old, g, dimer_env%rot%angle2 + pi/2.0_dp)
         gam = dgg/gg
         g = h
         h = -xi*dimer_env%cg_rot%norm_theta + gam*dimer_env%cg_rot%norm_h*dimer_env%cg_rot%nvec_old
         h = h - dot_product(h, dimer_env%nvec)*dimer_env%nvec
         norm_h = sqrt(dot_product(h, h))
         IF (norm_h < epsilon(0.0_dp)) THEN
            h = 0.0_dp
         ELSE
            h = h/norm_h
         END IF
         dimer_env%cg_rot%norm_h = norm_h
      END IF
      CALL timestop(handle)
 
   END SUBROUTINE get_conjugate_direction
 
END MODULE cg_utils