mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-06-30 00:44:37 +02:00
594 lines
19 KiB
Org Mode
594 lines
19 KiB
Org Mode
* Algorithm for the trust region
|
||
|
||
step_in_trust_region:
|
||
Computes the step in the trust region (delta)
|
||
(automatically sets at the iteration 0 and which evolves during the
|
||
process in function of the evolution of rho). The step is computing by
|
||
constraining its norm with a lagrange multiplier.
|
||
Since the calculation of the step is based on the Newton method, an
|
||
estimation of the gain in energy is given using the Taylors series
|
||
truncated at the second order (criterion_model).
|
||
If (DABS(criterion-criterion_model) < 1d-12) then
|
||
must_exit = .True.
|
||
else
|
||
must_exit = .False.
|
||
|
||
This estimation of the gain in energy is used by
|
||
is_step_cancel_trust_region to say if the step is accepted or cancelled.
|
||
|
||
If the step must be cancelled, the calculation restart from the same
|
||
hessian and gradient and recomputes the step but in a smaller trust
|
||
region and so on until the step is accepted. If the step is accepted
|
||
the hessian and the gradient are recomputed to produce a new step.
|
||
|
||
Example:
|
||
|
||
#+BEGIN_SRC f90 :comments org :tangle algo_trust.irp.f
|
||
! !### Initialization ###
|
||
! delta = 0d0
|
||
! nb_iter = 0 ! Must start at 0 !!!
|
||
! rho = 0.5d0
|
||
! not_converged = .True.
|
||
!
|
||
! ! ### TODO ###
|
||
! ! Compute the criterion before the loop
|
||
! call #your_criterion(prev_criterion)
|
||
!
|
||
! do while (not_converged)
|
||
! ! ### TODO ##
|
||
! ! Call your gradient
|
||
! ! Call you hessian
|
||
! call #your_gradient(v_grad) (1D array)
|
||
! call #your_hessian(H) (2D array)
|
||
!
|
||
! ! ### TODO ###
|
||
! ! Diagonalization of the hessian
|
||
! call diagonalization_hessian(n,H,e_val,w)
|
||
!
|
||
! cancel_step = .True. ! To enter in the loop just after
|
||
! ! Loop to Reduce the trust radius until the criterion decreases and rho >= thresh_rho
|
||
! do while (cancel_step)
|
||
!
|
||
! ! Hessian,gradient,Criterion -> x
|
||
! call trust_region_step_w_expected_e(tmp_n,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,tmp_x,must_exit)
|
||
!
|
||
! if (must_exit) then
|
||
! ! ### Message ###
|
||
! ! if step_in_trust_region sets must_exit on true for numerical reasons
|
||
! print*,'algo_trust1 sends the message : Exit'
|
||
! !### exit ###
|
||
! endif
|
||
!
|
||
! !### TODO ###
|
||
! ! Compute x -> m_x
|
||
! ! Compute m_x -> R
|
||
! ! Apply R and keep the previous MOs...
|
||
! ! Update/touch
|
||
! ! Compute the new criterion/energy -> criterion
|
||
!
|
||
! call #your_routine_1D_to_2D_antisymmetric_array(x,m_x)
|
||
! call #your_routine_2D_antisymmetric_array_to_rotation_matrix(m_x,R)
|
||
! call #your_routine_to_apply_the_rotation_matrix(R,prev_mos)
|
||
!
|
||
! TOUCH #your_variables
|
||
!
|
||
! call #your_criterion(criterion)
|
||
!
|
||
! ! Criterion -> step accepted or rejected
|
||
! call trust_region_is_step_cancelled(nb_iter,prev_criterion, criterion, criterion_model,rho,cancel_step)
|
||
!
|
||
! ! ### TODO ###
|
||
! !if (cancel_step) then
|
||
! ! Cancel the previous step (mo_coef = prev_mos if you keep them...)
|
||
! !endif
|
||
! #if (cancel_step) then
|
||
! #mo_coef = prev_mos
|
||
! #endif
|
||
!
|
||
! enddo
|
||
!
|
||
! !call save_mos() !### depend of the time for 1 iteration
|
||
!
|
||
! ! To exit the external loop if must_exit = .True.
|
||
! if (must_exit) then
|
||
! !### exit ###
|
||
! endif
|
||
!
|
||
! ! Step accepted, nb iteration + 1
|
||
! nb_iter = nb_iter + 1
|
||
!
|
||
! ! ### TODO ###
|
||
! !if (###Conditions###) then
|
||
! ! no_converged = .False.
|
||
! !endif
|
||
! #if (#your_conditions) then
|
||
! # not_converged = .False.
|
||
! #endif
|
||
!
|
||
! enddo
|
||
#+END_SRC
|
||
|
||
Variables:
|
||
|
||
Input:
|
||
| n | integer | m*(m-1)/2 |
|
||
| m | integer | number of mo in the mo_class |
|
||
| H(n,n) | double precision | Hessian |
|
||
| v_grad(n) | double precision | Gradient |
|
||
| W(n,n) | double precision | Eigenvectors of the hessian |
|
||
| e_val(n) | double precision | Eigenvalues of the hessian |
|
||
| criterion | double precision | Actual criterion |
|
||
| prev_criterion | double precision | Value of the criterion before the first iteration/after the previous iteration |
|
||
| rho | double precision | Given by is_step_cancel_trus_region |
|
||
| | | Agreement between the real function and the Taylor series (2nd order) |
|
||
| nb_iter | integer | Actual number of iterations |
|
||
|
||
Input/output:
|
||
| delta | double precision | Radius of the trust region |
|
||
|
||
Output:
|
||
| criterion_model | double precision | Predicted criterion after the rotation |
|
||
| x(n) | double precision | Step |
|
||
| must_exit | logical | If the program must exit the loop |
|
||
|
||
#+BEGIN_SRC f90 :comments org :tangle algo_trust.irp.f
|
||
subroutine trust_region_step_w_expected_e(n,H,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,x,must_exit)
|
||
|
||
include 'pi.h'
|
||
|
||
BEGIN_DOC
|
||
! Compute the step and the expected criterion/energy after the step
|
||
END_DOC
|
||
|
||
implicit none
|
||
|
||
! in
|
||
integer, intent(in) :: n, nb_iter
|
||
double precision, intent(in) :: H(n,n), W(n,n), v_grad(n)
|
||
double precision, intent(in) :: rho, prev_criterion
|
||
|
||
! inout
|
||
double precision, intent(inout) :: delta, e_val(n)
|
||
|
||
! out
|
||
double precision, intent(out) :: criterion_model, x(n)
|
||
logical, intent(out) :: must_exit
|
||
|
||
! internal
|
||
integer :: info
|
||
|
||
must_exit = .False.
|
||
|
||
call trust_region_step(n,nb_iter,v_grad,rho,e_val,W,x,delta)
|
||
|
||
call trust_region_expected_e(n,v_grad,H,x,prev_criterion,criterion_model)
|
||
|
||
! exit if DABS(prev_criterion - criterion_model) < 1d-12
|
||
if (DABS(prev_criterion - criterion_model) < thresh_model) then
|
||
print*,''
|
||
print*,'###############################################################################'
|
||
print*,'DABS(prev_criterion - criterion_model) <', thresh_model, 'stop the trust region'
|
||
print*,'###############################################################################'
|
||
print*,''
|
||
must_exit = .True.
|
||
endif
|
||
|
||
if (delta < thresh_delta) then
|
||
print*,''
|
||
print*,'##############################################'
|
||
print*,'Delta <', thresh_delta, 'stop the trust region'
|
||
print*,'##############################################'
|
||
print*,''
|
||
must_exit = .True.
|
||
endif
|
||
|
||
! Add after the call to this subroutine, a statement:
|
||
! "if (must_exit) then
|
||
! exit
|
||
! endif"
|
||
! in order to exit the optimization loop
|
||
|
||
end subroutine
|
||
#+END_SRC
|
||
|
||
Variables:
|
||
|
||
Input:
|
||
| nb_iter | integer | actual number of iterations |
|
||
| prev_criterion | double precision | criterion before the application of the step x |
|
||
| criterion | double precision | criterion after the application of the step x |
|
||
| criterion_model | double precision | predicted criterion after the application of x |
|
||
|
||
Output:
|
||
| rho | double precision | Agreement between the predicted criterion and the real new criterion |
|
||
| cancel_step | logical | If the step must be cancelled |
|
||
|
||
#+BEGIN_SRC f90 :comments org :tangle algo_trust.irp.f
|
||
subroutine trust_region_is_step_cancelled(nb_iter,prev_criterion, criterion, criterion_model,rho,cancel_step)
|
||
|
||
include 'pi.h'
|
||
|
||
BEGIN_DOC
|
||
! Compute if the step should be cancelled
|
||
END_DOC
|
||
|
||
implicit none
|
||
|
||
! in
|
||
double precision, intent(in) :: prev_criterion, criterion, criterion_model
|
||
|
||
! inout
|
||
integer, intent(inout) :: nb_iter
|
||
|
||
! out
|
||
logical, intent(out) :: cancel_step
|
||
double precision, intent(out) :: rho
|
||
|
||
! Computes rho
|
||
call trust_region_rho(prev_criterion,criterion,criterion_model,rho)
|
||
|
||
if (nb_iter == 0) then
|
||
nb_iter = 1 ! in order to enable the change of delta if the first iteration is cancelled
|
||
endif
|
||
|
||
! If rho < thresh_rho -> give something in output to cancel the step
|
||
if (rho >= thresh_rho) then !0.1d0) then
|
||
! The step is accepted
|
||
cancel_step = .False.
|
||
else
|
||
! The step is rejected
|
||
cancel_step = .True.
|
||
print*, '***********************'
|
||
print*, 'Step cancel : rho <', thresh_rho
|
||
print*, '***********************'
|
||
endif
|
||
|
||
end subroutine
|
||
#+END_SRC
|
||
|
||
** Template for MOs
|
||
#+BEGIN_SRC f90 :comments org :tangle trust_region_template_mos.txt
|
||
subroutine algo_trust_template(tmp_n, tmp_list_size, tmp_list)
|
||
|
||
implicit none
|
||
|
||
! Variables
|
||
|
||
! In
|
||
integer, intent(in) :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
|
||
|
||
! Out
|
||
! Rien ou un truc pour savoir si ça c'est bien passé
|
||
|
||
! Internal
|
||
double precision, allocatable :: e_val(:), W(:,:), tmp_R(:,:), R(:,:), tmp_x(:), tmp_m_x(:,:)
|
||
double precision, allocatable :: prev_mos(:,:)
|
||
double precision :: criterion, prev_criterion, criterion_model
|
||
double precision :: delta, rho
|
||
logical :: not_converged, cancel_step, must_exit, enforce_step_cancellation
|
||
integer :: nb_iter, info, nb_sub_iter
|
||
integer :: i,j,tmp_i,tmp_j
|
||
|
||
allocate(W(tmp_n, tmp_n),e_val(tmp_n),tmp_x(tmp_n),tmp_m_x(tmp_list_size, tmp_list_size))
|
||
allocate(tmp_R(tmp_list_size, tmp_list_size), R(mo_num, mo_num))
|
||
allocate(prev_mos(ao_num, mo_num))
|
||
|
||
! Provide the criterion, but unnecessary because it's done
|
||
! automatically
|
||
PROVIDE C_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER
|
||
|
||
! Initialization
|
||
delta = 0d0
|
||
nb_iter = 0 ! Must start at 0 !!!
|
||
rho = 0.5d0 ! Must start at 0.5
|
||
not_converged = .True. ! Must be true
|
||
|
||
! Compute the criterion before the loop
|
||
prev_criterion = C_PROVIDER
|
||
|
||
do while (not_converged)
|
||
|
||
print*,''
|
||
print*,'******************'
|
||
print*,'Iteration', nb_iter
|
||
print*,'******************'
|
||
print*,''
|
||
|
||
! The new hessian and gradient are computed at the end of the previous iteration
|
||
! Diagonalization of the hessian
|
||
call diagonalization_hessian(tmp_n, H_PROVIDER, e_val, W)
|
||
|
||
cancel_step = .True. ! To enter in the loop just after
|
||
nb_sub_iter = 0
|
||
|
||
! Loop to Reduce the trust radius until the criterion decreases and rho >= thresh_rho
|
||
do while (cancel_step)
|
||
|
||
print*,'-----------------------------'
|
||
print*,'Iteration:', nb_iter
|
||
print*,'Sub iteration:', nb_sub_iter
|
||
print*,'-----------------------------'
|
||
|
||
! Hessian,gradient,Criterion -> x
|
||
call trust_region_step_w_expected_e(tmp_n, H_PROVIDER, W, e_val, g_PROVIDER, &
|
||
prev_criterion, rho, nb_iter, delta, criterion_model, tmp_x, must_exit)
|
||
|
||
if (must_exit) then
|
||
! if step_in_trust_region sets must_exit on true for numerical reasons
|
||
print*,'trust_region_step_w_expected_e sent the message : Exit'
|
||
exit
|
||
endif
|
||
|
||
! 1D tmp -> 2D tmp
|
||
call vec_to_mat_v2(tmp_n, tmp_list_size, tmp_x, tmp_m_x)
|
||
|
||
! Rotation submatrix (square matrix tmp_list_size by tmp_list_size)
|
||
call rotation_matrix(tmp_m_x, tmp_list_size, tmp_R, tmp_list_size, tmp_list_size, info, enforce_step_cancellation)
|
||
|
||
if (enforce_step_cancellation) then
|
||
print*, 'Forces the step cancellation, too large error in the rotation matrix'
|
||
rho = 0d0
|
||
cycle
|
||
endif
|
||
|
||
! tmp_R to R, subspace to full space
|
||
call sub_to_full_rotation_matrix(tmp_list_size, tmp_list, tmp_R, R)
|
||
|
||
! Rotation of the MOs
|
||
call apply_mo_rotation(R, prev_mos)
|
||
|
||
! touch mo_coef
|
||
call clear_mo_map ! Only if you are using the bi-electronic integrals
|
||
! mo_coef becomes valid
|
||
! And avoid the recomputation of the providers which depend of mo_coef
|
||
TOUCH mo_coef C_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER
|
||
|
||
! To update the other parameters if needed
|
||
call #update_parameters()
|
||
|
||
! To enforce the program to provide new criterion after the update
|
||
! of the parameters
|
||
FREE C_PROVIDER
|
||
PROVIDE C_PROVIDER
|
||
criterion = C_PROVIDER
|
||
|
||
! Criterion -> step accepted or rejected
|
||
call trust_region_is_step_cancelled(nb_iter, prev_criterion, criterion, criterion_model, rho, cancel_step)
|
||
|
||
! Cancellation of the step ?
|
||
if (cancel_step) then
|
||
! Replacement by the previous MOs
|
||
mo_coef = prev_mos
|
||
! call save_mos() ! depends of the time for 1 iteration
|
||
|
||
! No need to clear_mo_map since we don't recompute the gradient and the hessian
|
||
! mo_coef becomes valid
|
||
! Avoid the recomputation of the providers which depend of mo_coef
|
||
TOUCH mo_coef H_PROVIDER g_PROVIDER C_PROVIDER cc_PROVIDER
|
||
else
|
||
! The step is accepted:
|
||
! criterion -> prev criterion
|
||
|
||
! The replacement "criterion -> prev criterion" is already done
|
||
! in trust_region_rho, so if the criterion does not have a reason
|
||
! to change, it will change nothing for the criterion and will
|
||
! force the program to provide the new hessian, gradient and
|
||
! convergence criterion for the next iteration.
|
||
! But in the case of orbital optimization we diagonalize the CI
|
||
! matrix after the "FREE" statement, so the criterion will change
|
||
|
||
FREE C_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER
|
||
PROVIDE C_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER
|
||
prev_criterion = C_PROVIDER
|
||
|
||
endif
|
||
|
||
nb_sub_iter = nb_sub_iter + 1
|
||
enddo
|
||
|
||
! call save_mos() ! depends of the time for 1 iteration
|
||
|
||
! To exit the external loop if must_exit = .True.
|
||
if (must_exit) then
|
||
exit
|
||
endif
|
||
|
||
! Step accepted, nb iteration + 1
|
||
nb_iter = nb_iter + 1
|
||
|
||
! Provide the convergence criterion
|
||
! Provide the gradient and the hessian for the next iteration
|
||
PROVIDE cc_PROVIDER
|
||
|
||
! To exit
|
||
if (dabs(cc_PROVIDER) < thresh_opt_max_elem_grad) then
|
||
not_converged = .False.
|
||
endif
|
||
|
||
if (nb_iter > optimization_max_nb_iter) then
|
||
not_converged = .False.
|
||
endif
|
||
|
||
if (delta < thresh_delta) then
|
||
not_converged = .False.
|
||
endif
|
||
|
||
enddo
|
||
|
||
! Save the final MOs
|
||
call save_mos()
|
||
|
||
! Diagonalization of the hessian
|
||
! (To see the eigenvalues at the end of the optimization)
|
||
call diagonalization_hessian(tmp_n, H_PROVIDER, e_val, W)
|
||
|
||
deallocate(e_val, W, tmp_R, R, tmp_x, prev_mos)
|
||
|
||
end
|
||
#+END_SRC
|
||
|
||
** Cartesian version
|
||
#+BEGIN_SRC f90 :comments org :tangle trust_region_template_xyz.txt
|
||
subroutine algo_trust_cartesian_template(tmp_n)
|
||
|
||
implicit none
|
||
|
||
! Variables
|
||
|
||
! In
|
||
integer, intent(in) :: tmp_n
|
||
|
||
! Out
|
||
! Rien ou un truc pour savoir si ça c'est bien passé
|
||
|
||
! Internal
|
||
double precision, allocatable :: e_val(:), W(:,:), tmp_x(:)
|
||
double precision :: criterion, prev_criterion, criterion_model
|
||
double precision :: delta, rho
|
||
logical :: not_converged, cancel_step, must_exit
|
||
integer :: nb_iter, nb_sub_iter
|
||
integer :: i,j
|
||
|
||
allocate(W(tmp_n, tmp_n),e_val(tmp_n),tmp_x(tmp_n))
|
||
|
||
PROVIDE C_PROVIDER X_PROVIDER H_PROVIDER g_PROVIDER
|
||
|
||
! Initialization
|
||
delta = 0d0
|
||
nb_iter = 0 ! Must start at 0 !!!
|
||
rho = 0.5d0 ! Must start at 0.5
|
||
not_converged = .True. ! Must be true
|
||
|
||
! Compute the criterion before the loop
|
||
prev_criterion = C_PROVIDER
|
||
|
||
do while (not_converged)
|
||
|
||
print*,''
|
||
print*,'******************'
|
||
print*,'Iteration', nb_iter
|
||
print*,'******************'
|
||
print*,''
|
||
|
||
if (nb_iter > 0) then
|
||
PROVIDE H_PROVIDER g_PROVIDER
|
||
endif
|
||
|
||
! Diagonalization of the hessian
|
||
call diagonalization_hessian(tmp_n, H_PROVIDER, e_val, W)
|
||
|
||
cancel_step = .True. ! To enter in the loop just after
|
||
nb_sub_iter = 0
|
||
|
||
! Loop to Reduce the trust radius until the criterion decreases and rho >= thresh_rho
|
||
do while (cancel_step)
|
||
|
||
print*,'-----------------------------'
|
||
print*,'Iteration:', nb_iter
|
||
print*,'Sub iteration:', nb_sub_iter
|
||
print*,'-----------------------------'
|
||
|
||
! Hessian,gradient,Criterion -> x
|
||
call trust_region_step_w_expected_e(tmp_n, H_PROVIDER, W, e_val, g_PROVIDER, &
|
||
prev_criterion, rho, nb_iter, delta, criterion_model, tmp_x, must_exit)
|
||
|
||
if (must_exit) then
|
||
! if step_in_trust_region sets must_exit on true for numerical reasons
|
||
print*,'trust_region_step_w_expected_e sent the message : Exit'
|
||
exit
|
||
endif
|
||
|
||
! New coordinates, check the sign
|
||
X_PROVIDER = X_PROVIDER - tmp_x
|
||
|
||
! touch X_PROVIDER
|
||
TOUCH X_PROVIDER H_PROVIDER g_PROVIDER cc_PROVIDER
|
||
|
||
! To update the other parameters if needed
|
||
call #update_parameters()
|
||
|
||
! New criterion
|
||
PROVIDE C_PROVIDER ! Unnecessary
|
||
criterion = C_PROVIDER
|
||
|
||
! Criterion -> step accepted or rejected
|
||
call trust_region_is_step_cancelled(nb_iter, prev_criterion, criterion, criterion_model, rho, cancel_step)
|
||
|
||
! Cancel the previous step
|
||
if (cancel_step) then
|
||
! Replacement by the previous coordinates, check the sign
|
||
X_PROVIDER = X_PROVIDER + tmp_x
|
||
|
||
! Avoid the recomputation of the hessian and the gradient
|
||
TOUCH X_PROVIDER H_PROVIDER g_PROVIDER C_PROVIDER cc_PROVIDER
|
||
endif
|
||
|
||
nb_sub_iter = nb_sub_iter + 1
|
||
enddo
|
||
|
||
! To exit the external loop if must_exit = .True.
|
||
if (must_exit) then
|
||
exit
|
||
endif
|
||
|
||
! Step accepted, nb iteration + 1
|
||
nb_iter = nb_iter + 1
|
||
|
||
PROVIDE cc_PROVIDER
|
||
|
||
! To exit
|
||
if (dabs(cc_PROVIDER) < thresh_opt_max_elem_grad) then
|
||
not_converged = .False.
|
||
endif
|
||
|
||
if (nb_iter > optimization_max_nb_iter) then
|
||
not_converged = .False.
|
||
endif
|
||
|
||
if (delta < thresh_delta) then
|
||
not_converged = .False.
|
||
endif
|
||
|
||
enddo
|
||
|
||
deallocate(e_val, W, tmp_x)
|
||
|
||
end
|
||
#+END_SRC
|
||
|
||
** Script template
|
||
#+BEGIN_SRC bash :tangle script_template_mos.sh
|
||
#!/bin/bash
|
||
|
||
your_file=
|
||
|
||
your_C_PROVIDER=
|
||
your_H_PROVIDER=
|
||
your_g_PROVIDER=
|
||
your_cc_PROVIDER=
|
||
|
||
sed "s/C_PROVIDER/$your_C_PROVIDER/g" trust_region_template_mos.txt > $your_file
|
||
sed -i "s/H_PROVIDER/$your_H_PROVIDER/g" $your_file
|
||
sed -i "s/g_PROVIDER/$your_g_PROVIDER/g" $your_file
|
||
sed -i "s/cc_PROVIDER/$your_cc_PROVIDER/g" $your_file
|
||
#+END_SRC
|
||
|
||
#+BEGIN_SRC bash :tangle script_template_xyz.sh
|
||
#!/bin/bash
|
||
|
||
your_file=
|
||
|
||
your_C_PROVIDER=
|
||
your_X_PROVIDER=
|
||
your_H_PROVIDER=
|
||
your_g_PROVIDER=
|
||
your_cc_PROVIDER=
|
||
|
||
sed "s/C_PROVIDER/$your_C_PROVIDER/g" trust_region_template_xyz.txt > $your_file
|
||
sed -i "s/X_PROVIDER/$your_X_PROVIDER/g" $your_file
|
||
sed -i "s/H_PROVIDER/$your_H_PROVIDER/g" $your_file
|
||
sed -i "s/g_PROVIDER/$your_g_PROVIDER/g" $your_file
|
||
sed -i "s/cc_PROVIDER/$your_cc_PROVIDER/g" $your_file
|
||
#+END_SRC
|
||
|