add mo optimization

src/mo_optimization/83.mo_optimization.bats

@@ -0,0 +1,62 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
+
+
+function run() {
+  thresh=2e-3
+  test_exe scf || skip
+  qp set_file $1
+  qp edit --check
+  qp reset -a
+  qp run scf
+  qp set_frozen_core
+  qp set determinants n_states 2
+  qp set determinants read_wf true
+  qp set mo_two_e_ints io_mo_two_e_integrals None
+  file="$(echo $1 | sed 's/.ezfio//g')"
+  qp run cis
+  qp run debug_gradient_list_opt > $file.debug_g.out
+  err3="$(grep 'Max error:' $file.debug_g.out | awk '{print $3}')"
+  qp run debug_hessian_list_opt > $file.debug_h1.out
+  err1="$(grep 'Max error:' $file.debug_h1.out | awk '{print $3}')"
+  qp run orb_opt > $file.opt1.out
+  energy1="$(grep 'State average energy:' $file.opt1.out | tail -n 1 | awk '{print $4}')"
+  qp set orbital_optimization optimization_method diag
+  qp reset -d 
+  qp run scf
+  qp run cis
+  qp run debug_hessian_list_opt > $file.debug_h2.out
+  err2="$(grep 'Max error_H:' $file.debug_h2.out | awk '{print $3}')"
+  qp run orb_opt > $file.opt2.out
+  energy2="$(grep 'State average energy:' $file.opt2.out | tail -n 1 | awk '{print $4}')"
+  qp set orbital_optimization optimization_method full
+  qp reset -d
+  qp run scf
+  eq $energy1 $2 $thresh
+  eq $energy2 $3 $thresh
+  eq $err1 0.0 1e-12
+  eq $err2 0.0 1e-12
+  eq $err3 0.0 1e-12
+}
+
+@test "b2_stretched" {
+run  b2_stretched.ezfio -48.9852901484277 -48.9852937541510
+}
+
+@test "h2o" {
+run  h2o.ezfio -75.9025622449206 -75.8691844585879
+}
+
+@test "h2s" {
+run  h2s.ezfio -398.576255809878 -398.574145943928
+}
+
+@test "hbo" {
+run  hbo.ezfio -99.9234823022109 -99.9234763597840
+}
+
+@test "hco" {
+run  hco.ezfio -113.204915552241 -113.204905207050
+}

src/mo_optimization/EZFIO.cfg

@@ -0,0 +1,29 @@
+[optimization_method]
+type: character*(32)
+doc: Define the kind of hessian for the orbital optimization full : full hessian, diag : diagonal hessian, none : no hessian
+interface: ezfio,provider,ocaml
+default: full
+
+[n_det_start]
+type: integer
+doc: Number of determinants after which the orbital optimization will start, n_det_start must be greater than 1. The algorithm does a cipsi until n_det > n_det_start and the optimization starts after
+interface: ezfio,provider,ocaml
+default: 5
+
+[n_det_max_opt]
+type: integer
+doc: Maximal number of the determinants in the wf for the orbital optimization (to stop the optimization if n_det > n_det_max_opt)
+interface: ezfio,provider,ocaml
+default: 200000
+
+[optimization_max_nb_iter]
+type: integer
+doc: Maximal number of iterations for the orbital optimization
+interface: ezfio,provider,ocaml
+default: 20
+
+[thresh_opt_max_elem_grad]
+type: double precision
+doc: Threshold for the convergence, the optimization exits when the biggest element in the gradient is smaller than thresh_optimization_max_elem_grad
+interface: ezfio,provider,ocaml
+default: 1.e-5

src/mo_optimization/NEED

@@ -0,0 +1,7 @@
+two_body_rdm
+hartree_fock
+cipsi
+davidson_undressed
+selectors_full
+generators_full
+utils_trust_region

src/mo_optimization/README.md

@@ -0,0 +1,74 @@
+# Orbital optimization
+
+## Methods  
+Different methods are available:  
+- full hessian  
+``` 
+qp set orbital_optimization optimization_method full  
+```  
+- diagonal hessian  
+``` 
+qp set orbital_optimization optimization_method diag  
+``` 
+- identity matrix  
+``` 
+qp set orbital_optimization optimization_method none  
+``` 
+
+After the optimization the ezfio contains the optimized orbitals
+ 
+## For a fixed number of determinants
+To optimize the MOs for the actual determinants:  
+``` 
+qp run orb_opt
+``` 
+ 
+## For a complete optimization, i.e, with a larger and larger wave function
+To optimize the MOs with a larger and larger wave function:  
+``` 
+qp run optimization  
+``` 
+
+The results are stored in the EZFIO in "mo_optimization/result_opt",
+with the following format:  
+(1) (2) (3) (4)  
+1: Number of determinants in the wf,  
+2: Cispi energy before the optimization,   
+3: Cipsi energy after the optimization,  
+4: Energy difference between (2) and (3).  
+ 
+The optimization process if the following: 
+- we do a first cipsi step to obtain a small number of determinants in the wf 
+- we run an orbital optimization for this wf 
+- we do a new cipsi step to double the number of determinants in the wf 
+- we run an orbital optimization for this wf 
+- ... 
+- we do that until the energy difference between (2) and (3) is  
+  smaller than the targeted accuracy for the cispi (targeted_accuracy_cipsi in qp edit) 
+  or the wf is larger than a given size (n_det_max_opt in qp_edit) 
+- after that you can reset your determinants (qp reset -d) and run a clean Cispi calculation  
+  
+### End of the optimization
+You can choos the number of determinants after what the 
+optimization will stop:
+```
+qp set orbital_optimization n_det_max_opt 1e5 # or any number
+```
+## Weight of the states
+You can change the weights of the differents states directly in qp edit.  
+It will affect ths weights used in the orbital optimization.
+
+# Tests
+To run the tests:  
+``` 
+qp test
+``` 
+
+# Org files
+The org files are stored in the directory org in order to avoid overwriting on user changes.
+The org files can be modified, to export the change to the source code, run
+```
+./TANGLE_org_mode.sh
+mv *.irp.f ../.
+```
+

src/mo_optimization/class.irp.f

@@ -0,0 +1,12 @@
+ BEGIN_PROVIDER [ logical, do_only_1h1p ]
+&BEGIN_PROVIDER [ logical, do_only_cas  ]
+&BEGIN_PROVIDER [ logical, do_ddci ]
+ implicit none
+ BEGIN_DOC
+ ! In the FCI case, all those are always false
+ END_DOC
+ do_only_1h1p = .False.
+ do_only_cas  = .False.
+ do_ddci = .False.
+END_PROVIDER
+

src/mo_optimization/constants.h

@@ -0,0 +1 @@
+  logical, parameter :: debug=.False.

src/mo_optimization/debug_gradient_list_opt.irp.f

@@ -0,0 +1,78 @@
+! Debug the gradient
+
+! *Program to check the gradient*
+
+! The program compares the result of the first and last code for the
+! gradient.
+
+! Provided:
+! | mo_num | integer | number of MOs |
+
+! Internal:
+! | n          | integer          | number of orbitals pairs (p,q) p<q |
+! | v_grad(n)  | double precision | Original gradient                  |
+! | v_grad2(n) | double precision | Gradient                           |
+! | i          | integer          | index                              |
+! | threshold  | double precision | threshold for the errors           |
+! | max_error  | double precision | maximal error in the gradient      |
+! | nb_error   | integer          | number of error in the gradient    |
+
+
+program debug_gradient_list
+  
+  implicit none
+
+  ! Variables
+
+  double precision, allocatable :: v_grad(:), v_grad2(:)
+  integer                       :: n,m
+  integer                       :: i
+  double precision              :: threshold
+  double precision              :: max_error, max_elem, norm
+  integer                       :: nb_error
+  
+  m = dim_list_act_orb
+  ! Definition of n  
+  n = m*(m-1)/2
+
+  PROVIDE mo_two_e_integrals_in_map ! Vérifier pour suppression
+
+  ! Allocation
+  allocate(v_grad(n), v_grad2(n))
+
+  ! Calculation
+
+  call diagonalize_ci ! Vérifier pour suppression
+
+  ! Gradient  
+  call gradient_list_opt(n,m,list_act,v_grad,max_elem,norm)
+  call first_gradient_list_opt(n,m,list_act,v_grad2)
+  
+  
+  v_grad = v_grad - v_grad2
+  nb_error = 0
+  max_error = 0d0 
+  threshold = 1d-12 
+
+  do i = 1, n
+    if (ABS(v_grad(i)) > threshold) then
+       print*,i,v_grad(i)
+       nb_error = nb_error + 1
+
+       if (ABS(v_grad(i)) > max_error) then
+         max_error = v_grad(i)
+       endif
+
+    endif
+  enddo
+ 
+  print*,''
+  print*,'Check the gradient' 
+  print*,'Threshold:', threshold
+  print*,'Nb error:', nb_error
+  print*,'Max error:', max_error
+
+  ! Deallocation
+  deallocate(v_grad,v_grad2)
+
+end program

src/mo_optimization/debug_gradient_opt.irp.f

@@ -0,0 +1,76 @@
+! Debug the gradient
+
+! *Program to check the gradient*
+
+! The program compares the result of the first and last code for the
+! gradient.
+
+! Provided:
+! | mo_num | integer | number of MOs |
+
+! Internal:
+! | n          | integer          | number of orbitals pairs (p,q) p<q |
+! | v_grad(n)  | double precision | Original gradient                  |
+! | v_grad2(n) | double precision | Gradient                           |
+! | i          | integer          | index                              |
+! | threshold  | double precision | threshold for the errors           |
+! | max_error  | double precision | maximal error in the gradient      |
+! | nb_error   | integer          | number of error in the gradient    |
+
+
+program debug_gradient
+  
+  implicit none
+
+  ! Variables
+
+  double precision, allocatable :: v_grad(:), v_grad2(:)
+  integer                       :: n
+  integer                       :: i
+  double precision              :: threshold
+  double precision              :: max_error, max_elem
+  integer                       :: nb_error
+  
+  ! Definition of n  
+  n = mo_num*(mo_num-1)/2
+
+  PROVIDE mo_two_e_integrals_in_map ! Vérifier pour suppression
+
+  ! Allocation
+  allocate(v_grad(n), v_grad2(n))
+
+  ! Calculation
+
+  call diagonalize_ci ! Vérifier pour suppression
+
+  ! Gradient  
+  call first_gradient_opt(n,v_grad)
+  call gradient_opt(n,v_grad2,max_elem)
+  
+  v_grad = v_grad - v_grad2
+  nb_error = 0
+  max_error = 0d0 
+  threshold = 1d-12 
+
+  do i = 1, n
+    if (ABS(v_grad(i)) > threshold) then
+       print*,v_grad(i)
+       nb_error = nb_error + 1
+
+       if (ABS(v_grad(i)) > max_error) then
+         max_error = v_grad(i)
+       endif
+
+    endif
+  enddo
+ 
+  print*,''
+  print*,'Check the gradient' 
+  print*,'Threshold :', threshold
+  print*,'Nb error :', nb_error
+  print*,'Max error :', max_error
+
+  ! Deallocation
+  deallocate(v_grad,v_grad2)
+
+end program

src/mo_optimization/debug_hessian_list_opt.irp.f

@@ -0,0 +1,147 @@
+! Debug the hessian
+
+! *Program to check the hessian matrix*
+
+! The program compares the result of the first and last code for the
+! hessian. First of all the 4D hessian and after the 2D hessian.
+
+! Provided:
+! | mo_num                     | integer | number of MOs                        |
+! | optimization_method        | string  | Method for the orbital optimization: |
+! |                            |         | - 'full' -> full hessian             |
+! |                            |         | - 'diag' -> diagonal hessian         |
+! | dim_list_act_orb           | integer | number of active MOs                 |
+! | list_act(dim_list_act_orb) | integer | list of the actives MOs              |
+! |                            |         |                                      |
+
+! Internal:
+! | m                                 | integer          | number of MOs in the list          |
+! |                                   |                  | (active MOs)                       |
+! | n                                 | integer          | number of orbitals pairs (p,q) p<q |
+! |                                   |                  | n = m*(m-1)/2                      |
+! | H(n,n)                            | double precision | Original hessian matrix (2D)       |
+! | H2(n,n)                           | double precision | Hessian matrix (2D)                |
+! | h_f(mo_num,mo_num,mo_num,mo_num)  | double precision | Original hessian matrix (4D)       |
+! | h_f2(mo_num,mo_num,mo_num,mo_num) | double precision | Hessian matrix (4D)                |
+! | i,j,p,q,k                         | integer          | indexes                            |
+! | threshold                         | double precision | threshold for the errors           |
+! | max_error                         | double precision | maximal error in the 4D hessian    |
+! | max_error_H                       | double precision | maximal error in the 2D hessian    |
+! | nb_error                          | integer          | number of errors in the 4D hessian |
+! | nb_error_H                        | integer          | number of errors in the 2D hessian |
+
+
+program debug_hessian_list_opt
+
+  implicit none
+
+  ! Variables
+
+  double precision, allocatable :: H(:,:),H2(:,:), h_f(:,:,:,:), h_f2(:,:,:,:)
+  integer                       :: n,m
+  integer                       :: i,j,k,l
+  double precision              :: max_error, max_error_H
+  integer                       :: nb_error, nb_error_H
+  double precision              :: threshold
+  
+  m = dim_list_act_orb !mo_num
+
+  ! Definition of n  
+  n = m*(m-1)/2
+
+  PROVIDE mo_two_e_integrals_in_map ! Vérifier pour suppression
+
+  ! Hessian 
+  if (optimization_method == 'full') then
+    print*,'Use the full hessian matrix'
+    allocate(H(n,n),H2(n,n))  
+    allocate(h_f(m,m,m,m),h_f2(m,m,m,m))
+
+    call hessian_list_opt(n,m,list_act,H,h_f)
+    call first_hessian_list_opt(n,m,list_act,H2,h_f2)
+    !call hessian_opt(n,H2,h_f2)
+
+    ! Difference
+    h_f = h_f - h_f2
+    H = H - H2
+    max_error = 0d0
+    nb_error = 0    
+    threshold = 1d-12
+
+    do l = 1, m
+      do k= 1, m
+        do j = 1, m
+          do i = 1, m
+            if (ABS(h_f(i,j,k,l)) > threshold) then
+              print*,h_f(i,j,k,l)
+              nb_error = nb_error + 1
+              if (ABS(h_f(i,j,k,l)) > ABS(max_error)) then
+                max_error = h_f(i,j,k,l)
+              endif
+            endif
+          enddo
+        enddo
+      enddo
+    enddo
+
+    max_error_H = 0d0
+    nb_error_H = 0
+
+    do j = 1, n
+      do i = 1, n
+        if (ABS(H(i,j)) > threshold) then
+          print*, H(i,j)
+          nb_error_H = nb_error_H + 1
+
+          if (ABS(H(i,j)) > ABS(max_error_H)) then
+            max_error_H = H(i,j)
+          endif
+
+        endif
+      enddo
+    enddo 
+
+    ! Deallocation
+    deallocate(H, H2, h_f, h_f2)
+
+  else
+
+    print*, 'Use the diagonal hessian matrix'
+    allocate(H(n,1),H2(n,1))
+    call diag_hessian_list_opt(n,m,list_act,H)
+    call first_diag_hessian_list_opt(n,m,list_act,H2)
+    
+    H = H - H2
+  
+    max_error_H = 0d0
+    nb_error_H = 0
+ 
+    do i = 1, n
+      if (ABS(H(i,1)) > threshold) then
+        print*, H(i,1)
+        nb_error_H = nb_error_H + 1
+ 
+        if (ABS(H(i,1)) > ABS(max_error_H)) then
+          max_error_H = H(i,1)
+        endif
+ 
+      endif
+    enddo
+   
+ endif
+  
+  print*,''
+  if (optimization_method == 'full') then
+    print*,'Check of the full hessian'
+    print*,'Threshold:', threshold
+    print*,'Nb error:', nb_error
+    print*,'Max error:', max_error
+    print*,''
+  else
+    print*,'Check of the diagonal hessian'
+  endif
+   
+  print*,'Nb error_H:', nb_error_H
+  print*,'Max error_H:', max_error_H
+ 
+end program

src/mo_optimization/debug_hessian_opt.irp.f

@@ -0,0 +1,171 @@
+! Debug the hessian
+
+! *Program to check the hessian matrix*
+
+! The program compares the result of the first and last code for the
+! hessian. First of all the 4D hessian and after the 2D hessian.
+
+! Provided:
+! | mo_num | integer | number of MOs |
+
+! Internal:
+! | n                                 | integer          | number of orbitals pairs (p,q) p<q |
+! | H(n,n)                            | double precision | Original hessian matrix (2D)       |
+! | H2(n,n)                           | double precision | Hessian matrix (2D)                |
+! | h_f(mo_num,mo_num,mo_num,mo_num)  | double precision | Original hessian matrix (4D)       |
+! | h_f2(mo_num,mo_num,mo_num,mo_num) | double precision | Hessian matrix (4D)                |
+! | method                            | integer          | - 1: full hessian                  |
+! |                                   |                  | - 2: diagonal hessian              |
+! | i,j,p,q,k                         | integer          | indexes                            |
+! | threshold                         | double precision | threshold for the errors           |
+! | max_error                         | double precision | maximal error in the 4D hessian    |
+! | max_error_H                       | double precision | maximal error in the 2D hessian    |
+! | nb_error                          | integer          | number of errors in the 4D hessian |
+! | nb_error_H                        | integer          | number of errors in the 2D hessian |
+
+
+program debug_hessian
+
+  implicit none
+
+  ! Variables
+
+  double precision, allocatable :: H(:,:),H2(:,:), h_f(:,:,:,:), h_f2(:,:,:,:)
+  integer                       :: n
+  integer                       :: i,j,k,l
+  double precision              :: max_error, max_error_H
+  integer                       :: nb_error, nb_error_H
+  double precision              :: threshold
+  
+  ! Definition of n  
+  n = mo_num*(mo_num-1)/2
+
+  PROVIDE mo_two_e_integrals_in_map ! Vérifier pour suppression
+
+  ! Allocation
+  allocate(H(n,n),H2(n,n))  
+  allocate(h_f(mo_num,mo_num,mo_num,mo_num),h_f2(mo_num,mo_num,mo_num,mo_num))
+
+  ! Calculation
+  
+  ! Hessian 
+  if (optimization_method == 'full') then 
+
+    print*,'Use the full hessian matrix'
+    call hessian_opt(n,H,h_f)
+    call first_hessian_opt(n,H2,h_f2)
+
+    ! Difference
+    h_f = h_f - h_f2
+    H = H - H2
+    max_error = 0d0
+    nb_error = 0    
+    threshold = 1d-12
+
+    do l = 1, mo_num
+      do k= 1, mo_num
+        do j = 1, mo_num
+          do i = 1, mo_num
+            if (ABS(h_f(i,j,k,l)) > threshold) then
+              print*,h_f(i,j,k,l)
+              nb_error = nb_error + 1
+              if (ABS(h_f(i,j,k,l)) > ABS(max_error)) then
+                max_error = h_f(i,j,k,l)
+              endif
+            endif
+          enddo
+        enddo
+      enddo
+    enddo
+
+   max_error_H = 0d0
+   nb_error_H = 0
+
+   do j = 1, n
+     do i = 1, n
+       if (ABS(H(i,j)) > threshold) then
+         print*, H(i,j)
+         nb_error_H = nb_error_H + 1
+
+         if (ABS(H(i,j)) > ABS(max_error_H)) then
+           max_error_H = H(i,j)
+         endif
+
+       endif
+     enddo
+   enddo 
+
+  elseif (optimization_method == 'diag') then
+
+    print*, 'Use the diagonal hessian matrix'
+    call diag_hessian_opt(n,H,h_f)
+    call first_diag_hessian_opt(n,H2,h_f2)
+
+    h_f = h_f - h_f2
+    max_error = 0d0
+    nb_error = 0
+    threshold = 1d-12
+
+    do l = 1, mo_num
+      do k = 1, mo_num
+        do j = 1, mo_num
+          do i = 1, mo_num
+
+            if (ABS(h_f(i,j,k,l)) > threshold) then
+
+              print*,h_f(i,j,k,l)
+              nb_error = nb_error + 1
+
+              if (ABS(h_f(i,j,k,l)) > ABS(max_error)) then
+                max_error = h_f(i,j,k,l)
+              endif
+
+            endif
+
+          enddo
+        enddo
+      enddo
+    enddo
+
+    h=H-H2
+  
+    max_error_H = 0d0
+    nb_error_H = 0
+ 
+    do j = 1, n
+      do i = 1, n
+        if (ABS(H(i,j)) > threshold) then
+          print*, H(i,j)
+          nb_error_H = nb_error_H + 1
+ 
+          if (ABS(H(i,j)) > ABS(max_error_H)) then
+            max_error_H = H(i,j)
+          endif
+ 
+        endif
+      enddo
+    enddo
+  
+  else
+    print*,'Unknown optimization_method, please select full, diag'
+    call abort
+  endif
+  
+  print*,''
+  if (optimization_method == 'full') then
+    print*,'Check the full hessian'
+  else
+    print*,'Check the diagonal hessian'
+  endif
+   
+  print*,'Threshold :', threshold
+  print*,'Nb error :', nb_error
+  print*,'Max error :', max_error
+  print*,''
+  print*,'Nb error_H :', nb_error_H
+  print*,'Max error_H :', max_error_H
+ 
+  ! Deallocation
+  deallocate(H,H2,h_f,h_f2)
+
+end program

src/mo_optimization/diagonalization_hessian.irp.f

@@ -0,0 +1,136 @@
+! Diagonalization of the hessian
+
+! Just a matrix diagonalization using Lapack
+
+! Input:
+! | n      | integer          | mo_num*(mo_num-1)/2 |
+! | H(n,n) | double precision | hessian             |
+
+! Output:
+! | e_val(n) | double precision | eigenvalues of the hessian  |
+! | w(n,n)   | double precision | eigenvectors of the hessian |
+
+! Internal:
+! | nb_negative_nv | integer          | number of negative eigenvalues                  |
+! | lwork          | integer          | for Lapack                                      |
+! | work(lwork,n)  | double precision | temporary array for Lapack                      |
+! | info           | integer          | if 0 -> ok, else problem in the diagonalization |
+! | i,j            | integer          | dummy indexes                                   |
+
+
+subroutine diagonalization_hessian(n,H,e_val,w)
+
+  include 'constants.h'
+
+  implicit none
+
+  ! Variables
+
+  ! in
+  integer, intent(in) :: n
+  double precision, intent(in) :: H(n,n)
+
+  ! out
+  double precision, intent(out) :: e_val(n), w(n,n)
+
+  ! internal
+  double precision, allocatable :: work(:,:)
+  integer, allocatable          :: key(:)
+  integer                       :: info,lwork
+  integer                       :: i,j
+  integer                       :: nb_negative_vp
+  double precision              :: t1,t2,t3,max_elem
+
+  print*,''
+  print*,'---Diagonalization_hessian---'
+
+  call wall_time(t1)
+
+  if (optimization_method == 'full') then
+    ! Allocation
+    ! For Lapack
+    lwork=3*n-1
+
+    allocate(work(lwork,n))
+
+    ! Calculation
+
+    ! Copy the hessian matrix, the eigenvectors will be store in W
+    W=H
+
+    ! Diagonalization of the hessian
+    call dsyev('V','U',n,W,size(W,1),e_val,work,lwork,info)
+
+    if (info /= 0) then
+        print*, 'Error diagonalization : diagonalization_hessian'
+        print*, 'info = ', info
+        call ABORT
+    endif
+
+    if (debug) then
+      print *, 'vp Hess:'
+      write(*,'(100(F10.5))')  real(e_val(:))
+    endif
+
+    ! Number of negative eigenvalues
+    max_elem = 0d0
+    nb_negative_vp = 0
+    do i = 1, n
+      if (e_val(i) < 0d0) then
+        nb_negative_vp = nb_negative_vp + 1
+        if (e_val(i) < max_elem) then
+          max_elem = e_val(i)
+        endif
+        !print*,'e_val < 0 :', e_val(i)
+      endif
+    enddo
+    print*,'Number of negative eigenvalues:', nb_negative_vp
+    print*,'Lowest eigenvalue:',max_elem
+
+    !nb_negative_vp = 0
+    !do i = 1, n
+    !  if (e_val(i) < -thresh_eig) then
+    !    nb_negative_vp = nb_negative_vp + 1
+    !  endif
+    !enddo
+    !print*,'Number of negative eigenvalues <', -thresh_eig,':', nb_negative_vp  
+
+    ! Deallocation
+    deallocate(work)
+
+  elseif (optimization_method == 'diag') then
+    ! Diagonalization of the diagonal hessian by hands
+    allocate(key(n))
+    
+    do i = 1, n
+      e_val(i) = H(i,i)
+    enddo
+
+    ! Key list for dsort
+    do i = 1, n 
+      key(i) = i
+    enddo
+
+    ! Sort of the eigenvalues
+    call dsort(e_val, key, n)
+
+    ! Eigenvectors
+    W = 0d0
+    do i = 1, n
+      j = key(i)
+      W(j,i) = 1d0
+    enddo
+
+    deallocate(key)
+  else
+    print*,'Diagonalization_hessian, abort'
+    call abort
+  endif
+
+  call wall_time(t2)
+  t3 = t2 - t1
+  print*,'Time in diagonalization_hessian:', t3
+
+  print*,'---End diagonalization_hessian---'
+
+end subroutine

src/mo_optimization/first_diagonal_hessian_list_opt.irp.f

@@ -0,0 +1,372 @@
+subroutine first_diag_hessian_list_opt(tmp_n,m,list,H)!, h_tmpr)
+
+  include 'constants.h' 
+
+  implicit none
+
+  !===========================================================================
+  ! Compute the diagonal hessian of energy with respects to orbital rotations
+  !===========================================================================
+
+  !===========
+  ! Variables 
+  !===========
+ 
+  ! in
+  integer, intent(in)           :: tmp_n, m, list(m) 
+  ! tmp_n        : integer, tmp_n = m*(m-1)/2
+ 
+  ! out
+  double precision, intent(out) :: H(tmp_n)!, h_tmpr(m,m,m,m)
+  ! H        : n by n double precision matrix containing the 2D hessian
+  
+  ! internal
+  double precision, allocatable :: hessian(:,:,:,:), tmp(:,:),h_tmpr(:,:,:,:)
+  integer                       :: p,q, tmp_p,tmp_q
+  integer                       :: r,s,t,u,v,tmp_r,tmp_s,tmp_t,tmp_u,tmp_v
+  integer                       :: pq,rs,tmp_pq,tmp_rs
+  double precision              :: t1,t2,t3
+  ! hessian  : mo_num 4D double precision matrix containing the hessian before the permutations
+  ! h_tmpr   : mo_num 4D double precision matrix containing the hessian after the permutations
+  ! p,q,r,s  : integer, indexes of the 4D hessian matrix
+  ! t,u,v    : integer, indexes to compute hessian elements
+  ! pq,rs    : integer, indexes for the conversion from 4D to 2D hessian matrix
+  ! t1,t2,t3 : double precision, t3 = t2 - t1, time to compute the hessian 
+  
+  ! Function
+  double precision :: get_two_e_integral
+  ! get_two_e_integral : double precision function, two e integrals
+ 
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo_alpha, one_e_dm_mo_beta : one body density matrix
+  ! two_e_dm_mo : two body density matrix
+
+  print*,'---first_diag_hess_list---'
+
+  !============
+  ! Allocation
+  !============
+  
+  allocate(hessian(m,m,m,m),tmp(tmp_n,tmp_n),h_tmpr(mo_num,mo_num,mo_num,mo_num))
+
+  !=============
+  ! Calculation
+  !=============
+
+  ! From Anderson et. al. (2014) 
+  ! The Journal of Chemical Physics 141, 244104 (2014); doi: 10.1063/1.4904384
+
+  ! LaTeX formula :
+
+  !\begin{align*}
+  !H_{pq,rs} &= \dfrac{\partial^2 E(x)}{\partial x_{pq}^2} \\
+  !&= \mathcal{P}_{pq} \mathcal{P}_{rs} [ \frac{1}{2} \sum_u [\delta_{qr}(h_p^u \gamma_u^s + h_u^s \gamma_p^u) 
+  !+ \delta_{ps}(h_r^u \gamma_u^q + h_u^q \gamma_u^r)]
+  !-(h_p^s \gamma_r^q + h_r^q \gamma_p^s) \\
+  !&+ \frac{1}{2} \sum_{tuv} [\delta_{qr}(v_{pt}^{uv} \Gamma_{uv}^{st} +v_{uv}^{st} \Gamma_{pt}^{uv}) 
+  !+ \delta_{ps}(v_{uv}^{qt} \Gamma_{rt}^{uv} + v_{rt}^{uv}\Gamma_{uv}^{qt})] \\
+  !&+ \sum_{uv} (v_{pr}^{uv} \Gamma_{uv}^{qs} + v_{uv}^{qs}  \Gamma_{ps}^{uv}) \\
+  !&- \sum_{tu} (v_{pu}^{st} \Gamma_{rt}^{qu}+v_{pu}^{tr} \Gamma_{tr}^{qu}+v_{rt}^{qu}\Gamma_{pu}^{st} + v_{tr}^{qu}\Gamma_{pu}^{ts}) 
+  !\end{align*} 
+
+  !================
+  ! Initialization
+  !================
+  hessian = 0d0
+
+  CALL wall_time(t1)
+
+  !========================
+  ! First line, first term
+  !========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+           
+            if (q==r) then
+              do u = 1, mo_num
+
+                hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                    mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
+                  + mo_one_e_integrals(s,u) * one_e_dm_mo(p,u))
+
+              enddo
+            endif
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! First line, second term
+  !=========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+            .or. ((p==s) .and. (q==r))) then
+
+             if (p==s) then
+               do u = 1, mo_num
+
+                    hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                      mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
+                    + mo_one_e_integrals(q,u) * one_e_dm_mo(r,u))
+               enddo
+             endif
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! First line, third term
+  !========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+         
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+
+            hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+            - mo_one_e_integrals(s,p) * one_e_dm_mo(r,q) &
+            - mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
+
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! Second line, first term
+  !=========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+
+              if (q==r) then
+                do t = 1, mo_num
+                  do u = 1, mo_num
+                    do v = 1, mo_num
+
+                      hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * (  &
+                        get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
+                      + get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v))
+
+                    enddo
+                  enddo
+                enddo
+              endif
+            endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !==========================
+  ! Second line, second term
+  !==========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+              .or. ((p==s) .and. (q==r))) then
+
+             if (p==s) then
+               do t = 1, mo_num
+                 do u = 1, mo_num
+                   do v = 1, mo_num
+
+                     hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                       get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
+                     + get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t))
+
+                   enddo
+                 enddo
+               enddo
+             endif
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! Third line, first term
+  !========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+                .or. ((p==s) .and. (q==r))) then
+
+            do u = 1, mo_num
+              do v = 1, mo_num
+
+                hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+                 + get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
+                 + get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v)
+
+              enddo
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! Third line, second term
+  !=========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+           .or. ((p==s) .and. (q==r))) then
+
+            do t = 1, mo_num
+              do u = 1, mo_num
+
+                hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+                 - get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
+                 - get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
+                 - get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
+                 - get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
+
+              enddo
+            enddo
+
+          endif     
+  
+        enddo
+      enddo
+    enddo
+  enddo
+
+    CALL wall_time(t2)
+    t2 = t2 - t1
+    print*, 'Time to compute the hessian :', t2
+
+  !==============
+  ! Permutations
+  !==============
+ 
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  do tmp_r = 1, m
+    do tmp_s = 1, m
+      do tmp_q = 1, m
+        do tmp_p = 1, m
+
+          h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) - hessian(tmp_q,tmp_p,tmp_r,tmp_s) &
+                                          - hessian(tmp_p,tmp_q,tmp_s,tmp_r) + hessian(tmp_q,tmp_p,tmp_s,tmp_r)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! 4D matrix -> 2D matrix
+  !========================
+  
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  ! 4D mo_num matrix to 2D n matrix
+  do tmp_rs = 1, tmp_n
+    call vec_to_mat_index(tmp_rs,tmp_r,tmp_s)
+    do tmp_pq = 1, tmp_n
+      call vec_to_mat_index(tmp_pq,tmp_p,tmp_q)
+      tmp(tmp_pq,tmp_rs) = h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s)   
+    enddo
+  enddo
+
+  do p = 1, tmp_n
+    H(p) = tmp(p,p)
+  enddo
+
+  ! Display
+  if (debug) then 
+    print*,'2D diag Hessian matrix'
+    do tmp_pq = 1, tmp_n
+      write(*,'(100(F10.5))') tmp(tmp_pq,:)
+    enddo 
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(hessian,h_tmpr,tmp)
+
+  print*,'---End first_diag_hess_list---'
+
+end subroutine

src/mo_optimization/first_diagonal_hessian_opt.irp.f

@@ -0,0 +1,344 @@
+subroutine first_diag_hessian_opt(n,H, h_tmpr)
+
+  include 'constants.h' 
+
+  implicit none
+
+  !===========================================================================
+  ! Compute the diagonal hessian of energy with respects to orbital rotations
+  !===========================================================================
+
+  !===========
+  ! Variables 
+  !===========
+ 
+  ! in
+  integer, intent(in)           :: n 
+  ! n        : integer, n = mo_num*(mo_num-1)/2
+ 
+  ! out
+  double precision, intent(out) :: H(n,n), h_tmpr(mo_num,mo_num,mo_num,mo_num)
+  ! H        : n by n double precision matrix containing the 2D hessian
+  
+  ! internal
+  double precision, allocatable :: hessian(:,:,:,:)
+  integer                       :: p,q
+  integer                       :: r,s,t,u,v
+  integer                       :: pq,rs
+  double precision              :: t1,t2,t3
+  ! hessian  : mo_num 4D double precision matrix containing the hessian before the permutations
+  ! h_tmpr   : mo_num 4D double precision matrix containing the hessian after the permutations
+  ! p,q,r,s  : integer, indexes of the 4D hessian matrix
+  ! t,u,v    : integer, indexes to compute hessian elements
+  ! pq,rs    : integer, indexes for the conversion from 4D to 2D hessian matrix
+  ! t1,t2,t3 : double precision, t3 = t2 - t1, time to compute the hessian 
+  
+  ! Function
+  double precision :: get_two_e_integral
+  ! get_two_e_integral : double precision function, two e integrals
+ 
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo_alpha, one_e_dm_mo_beta : one body density matrix
+  ! two_e_dm_mo : two body density matrix
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(hessian(mo_num,mo_num,mo_num,mo_num))!,h_tmpr(mo_num,mo_num,mo_num,mo_num))
+
+  !=============
+  ! Calculation
+  !=============
+
+  if (debug) then
+          print*,'Enter in first_diag_hessien'
+  endif
+
+  ! From Anderson et. al. (2014) 
+  ! The Journal of Chemical Physics 141, 244104 (2014); doi: 10.1063/1.4904384
+
+  ! LaTeX formula :
+
+  !\begin{align*}
+  !H_{pq,rs} &= \dfrac{\partial^2 E(x)}{\partial x_{pq}^2} \\
+  !&= \mathcal{P}_{pq} \mathcal{P}_{rs} [ \frac{1}{2} \sum_u [\delta_{qr}(h_p^u \gamma_u^s + h_u^s \gamma_p^u) 
+  !+ \delta_{ps}(h_r^u \gamma_u^q + h_u^q \gamma_u^r)]
+  !-(h_p^s \gamma_r^q + h_r^q \gamma_p^s) \\
+  !&+ \frac{1}{2} \sum_{tuv} [\delta_{qr}(v_{pt}^{uv} \Gamma_{uv}^{st} +v_{uv}^{st} \Gamma_{pt}^{uv}) 
+  !+ \delta_{ps}(v_{uv}^{qt} \Gamma_{rt}^{uv} + v_{rt}^{uv}\Gamma_{uv}^{qt})] \\
+  !&+ \sum_{uv} (v_{pr}^{uv} \Gamma_{uv}^{qs} + v_{uv}^{qs}  \Gamma_{ps}^{uv}) \\
+  !&- \sum_{tu} (v_{pu}^{st} \Gamma_{rt}^{qu}+v_{pu}^{tr} \Gamma_{tr}^{qu}+v_{rt}^{qu}\Gamma_{pu}^{st} + v_{tr}^{qu}\Gamma_{pu}^{ts}) 
+  !\end{align*} 
+
+  !================
+  ! Initialization
+  !================
+  hessian = 0d0
+
+  CALL wall_time(t1)
+
+  !========================
+  ! First line, first term
+  !========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+           
+            if (q==r) then
+              do u = 1, mo_num
+
+                hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                    mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
+                  + mo_one_e_integrals(s,u) * one_e_dm_mo(p,u))
+
+              enddo
+            endif
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! First line, second term
+  !=========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+            .or. ((p==s) .and. (q==r))) then
+
+             if (p==s) then
+               do u = 1, mo_num
+
+                    hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                      mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
+                    + mo_one_e_integrals(q,u) * one_e_dm_mo(r,u))
+               enddo
+             endif
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! First line, third term
+  !========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+         
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+
+            hessian(p,q,r,s) = hessian(p,q,r,s) &
+            - mo_one_e_integrals(s,p) * one_e_dm_mo(r,q) &
+            - mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
+
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! Second line, first term
+  !=========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+
+              if (q==r) then
+                do t = 1, mo_num
+                  do u = 1, mo_num
+                    do v = 1, mo_num
+
+                      hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * (  &
+                        get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
+                      + get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v))
+
+                    enddo
+                  enddo
+                enddo
+              endif
+            endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !==========================
+  ! Second line, second term
+  !==========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+              .or. ((p==s) .and. (q==r))) then
+
+             if (p==s) then
+               do t = 1, mo_num
+                 do u = 1, mo_num
+                   do v = 1, mo_num
+
+                     hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                       get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
+                     + get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t))
+
+                   enddo
+                 enddo
+               enddo
+             endif
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! Third line, first term
+  !========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+                .or. ((p==s) .and. (q==r))) then
+
+            do u = 1, mo_num
+              do v = 1, mo_num
+
+                hessian(p,q,r,s) = hessian(p,q,r,s) &
+                 + get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
+                 + get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v)
+
+              enddo
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! Third line, second term
+  !=========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+           .or. ((p==s) .and. (q==r))) then
+
+            do t = 1, mo_num
+              do u = 1, mo_num
+
+                hessian(p,q,r,s) = hessian(p,q,r,s) &
+                 - get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
+                 - get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
+                 - get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
+                 - get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
+
+              enddo
+            enddo
+
+          endif     
+  
+        enddo
+      enddo
+    enddo
+  enddo
+
+    CALL wall_time(t2)
+    t2 = t2 - t1
+    print*, 'Time to compute the hessian :', t2
+
+  !==============
+  ! Permutations
+  !==============
+ 
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  do r = 1, mo_num
+    do s = 1, mo_num
+      do q = 1, mo_num
+        do p = 1, mo_num
+
+          h_tmpr(p,q,r,s) = (hessian(p,q,r,s) - hessian(q,p,r,s) - hessian(p,q,s,r) + hessian(q,p,s,r))
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! 4D matrix -> 2D matrix
+  !========================
+  
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  ! 4D mo_num matrix to 2D n matrix
+  do rs = 1, n
+    call vec_to_mat_index(rs,r,s)
+    do pq = 1, n
+      call vec_to_mat_index(pq,p,q)
+      H(pq,rs) = h_tmpr(p,q,r,s)   
+    enddo
+  enddo
+
+  ! Display
+  if (debug) then 
+    print*,'2D diag Hessian matrix'
+    do pq = 1, n
+      write(*,'(100(F10.5))') H(pq,:)
+    enddo 
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(hessian)
+
+  if (debug) then
+    print*,'Leave first_diag_hessien'
+  endif
+
+end subroutine

src/mo_optimization/first_gradient_list_opt.irp.f

@@ -0,0 +1,125 @@
+! First gradient
+
+subroutine first_gradient_list_opt(tmp_n,m,list,v_grad)
+
+  include 'constants.h'
+
+  implicit none
+
+  !===================================================================
+  ! Compute the gradient of energy with respects to orbital rotations
+  !===================================================================
+
+  ! Check if read_wf = true, else :
+  ! qp set determinant read_wf true
+
+  ! in
+  integer, intent(in) :: tmp_n,m,list(m)
+  ! n : integer, n = m*(m-1)/2
+  ! m = list_size
+  
+  ! out
+  double precision, intent(out) :: v_grad(tmp_n)
+  ! v_grad : double precision vector of length n containeing the gradient
+
+  ! internal
+  double precision, allocatable :: grad(:,:),A(:,:)
+  double precision :: norm
+  integer :: i,p,q,r,s,t,tmp_i,tmp_p,tmp_q,tmp_r,tmp_s,tmp_t
+  ! grad : double precision matrix containing the gradient before the permutation
+  ! A : double precision matrix containing the gradient after the permutation
+  ! norm : double precision number, the norm of the vector gradient
+  ! i,p,q,r,s,t : integer, indexes 
+  ! istate : integer, the electronic state
+
+  ! Function
+  double precision :: get_two_e_integral, norm2
+  ! get_two_e_integral :  double precision function that gives the two e integrals
+  ! norm2 : double precision function that gives the norm of a vector
+ 
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo : one body density matrix (state average)
+  ! two_e_dm_mo : two body density matrix (state average)
+
+   print*,'---first_gradient_list---'
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(grad(m,m),A(m,m))
+
+  !=============
+  ! Calculation
+  !=============
+
+  v_grad = 0d0
+  grad = 0d0
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      !grad(tmp_p,tmp_q) = 0d0
+      do r = 1, mo_num
+        grad(tmp_p,tmp_q) = grad(tmp_p,tmp_q) + mo_one_e_integrals(p,r) * one_e_dm_mo(r,q) &
+                               - mo_one_e_integrals(r,q) * one_e_dm_mo(p,r)
+
+      enddo
+
+      do r = 1, mo_num
+        do s = 1, mo_num
+          do t = 1, mo_num
+
+            grad(tmp_p,tmp_q) = grad(tmp_p,tmp_q) &
+                + get_two_e_integral(p,t,r,s,mo_integrals_map) * two_e_dm_mo(r,s,q,t) &
+                - get_two_e_integral(r,s,q,t,mo_integrals_map) * two_e_dm_mo(p,t,r,s)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ! Conversion mo_num*mo_num matrix to mo_num(mo_num-1)/2 vector
+  do tmp_i = 1, tmp_n
+    call vec_to_mat_index(tmp_i,tmp_p,tmp_q)
+    v_grad(tmp_i)=(grad(tmp_p,tmp_q) - grad(tmp_q,tmp_p))
+  enddo  
+
+  ! Display, vector containing the gradient elements 
+  if (debug) then  
+    print*,'Vector containing the gradient :'
+    write(*,'(100(F10.5))') v_grad(1:tmp_n)
+  endif  
+
+  ! Norm of the vector
+  norm = norm2(v_grad)
+  print*, 'Norm : ', norm
+
+  ! Matrix gradient
+  A = 0d0
+  do tmp_q = 1, m
+    do tmp_p = 1, m
+      A(tmp_p,tmp_q) = grad(tmp_p,tmp_q) - grad(tmp_q,tmp_p)
+    enddo
+  enddo
+
+  ! Display, matrix containting the gradient elements
+  if (debug) then
+    print*,'Matrix containing the gradient :'
+    do tmp_i = 1, m
+      write(*,'(100(E12.5))') A(tmp_i,1:m)
+    enddo
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(grad,A)
+
+  print*,'---End first_gradient_list---'
+
+end subroutine

src/mo_optimization/first_gradient_opt.irp.f

@@ -0,0 +1,128 @@
+! First gradient
+
+subroutine first_gradient_opt(n,v_grad)
+
+  include 'constants.h'
+
+  implicit none
+
+  !===================================================================
+  ! Compute the gradient of energy with respects to orbital rotations
+  !===================================================================
+
+  ! Check if read_wf = true, else :
+  ! qp set determinant read_wf true
+
+  END_DOC
+
+  ! in
+  integer, intent(in) :: n
+  ! n : integer, n = mo_num*(mo_num-1)/2
+  
+  ! out
+  double precision, intent(out) :: v_grad(n)
+  ! v_grad : double precision vector of length n containeing the gradient
+
+  ! internal
+  double precision, allocatable :: grad(:,:),A(:,:)
+  double precision :: norm
+  integer :: i,p,q,r,s,t
+  integer :: istate
+  ! grad : double precision matrix containing the gradient before the permutation
+  ! A : double precision matrix containing the gradient after the permutation
+  ! norm : double precision number, the norm of the vector gradient
+  ! i,p,q,r,s,t : integer, indexes 
+  ! istate : integer, the electronic state
+
+  ! Function
+  double precision :: get_two_e_integral, norm2
+  ! get_two_e_integral :  double precision function that gives the two e integrals
+  ! norm2 : double precision function that gives the norm of a vector
+ 
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo : one body density matrix (state average)
+  ! two_e_dm_mo : two body density matrix (state average)
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(grad(mo_num,mo_num),A(mo_num,mo_num))
+
+  !=============
+  ! Calculation
+  !=============
+
+  if (debug) then
+    print*,'---first_gradient---'
+  endif
+
+  v_grad = 0d0
+
+  do p = 1, mo_num
+    do q = 1, mo_num
+      grad(p,q) = 0d0
+      do r = 1, mo_num
+        grad(p,q) = grad(p,q) + mo_one_e_integrals(p,r) * one_e_dm_mo(r,q) &
+                               - mo_one_e_integrals(r,q) * one_e_dm_mo(p,r)
+
+      enddo
+
+      do r = 1, mo_num
+        do s = 1, mo_num
+          do t= 1, mo_num
+
+            grad(p,q) = grad(p,q) &
+                + get_two_e_integral(p,t,r,s,mo_integrals_map) * two_e_dm_mo(r,s,q,t) &
+                - get_two_e_integral(r,s,q,t,mo_integrals_map) * two_e_dm_mo(p,t,r,s)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ! Conversion mo_num*mo_num matrix to mo_num(mo_num-1)/2 vector
+  do i=1,n
+    call vec_to_mat_index(i,p,q)
+    v_grad(i)=(grad(p,q) - grad(q,p))
+  enddo  
+
+  ! Display, vector containing the gradient elements 
+  if (debug) then  
+    print*,'Vector containing the gradient :'
+    write(*,'(100(F10.5))') v_grad(1:n)
+  endif  
+
+  ! Norm of the vector
+  norm = norm2(v_grad)
+  print*, 'Norm : ', norm
+
+  ! Matrix gradient
+  A = 0d0
+  do q=1,mo_num
+    do p=1,mo_num
+      A(p,q) = grad(p,q) - grad(q,p)
+    enddo
+  enddo
+
+  ! Display, matrix containting the gradient elements
+  if (debug) then    
+    print*,'Matrix containing the gradient :'
+    do i = 1, mo_num
+      write(*,'(100(E12.5))') A(i,1:mo_num)
+    enddo
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(grad,A)
+
+  if (debug) then
+    print*,'---End first_gradient---'
+  endif
+
+end subroutine

src/mo_optimization/first_hessian_list_opt.irp.f

@@ -0,0 +1,365 @@
+subroutine first_hessian_list_opt(tmp_n,m,list,H,h_tmpr)
+
+  include 'constants.h' 
+
+  implicit none
+
+  !==================================================================
+  ! Compute the hessian of energy with respects to orbital rotations
+  !==================================================================
+
+  !===========
+  ! Variables  
+  !===========
+
+  ! in
+  integer, intent(in)           :: tmp_n, m, list(m)
+  !tmp_n         : integer, tmp_n = m*(m-1)/2
+  
+  ! out
+  double precision, intent(out) :: H(tmp_n,tmp_n),h_tmpr(m,m,m,m)
+  ! H        : n by n double precision matrix containing the 2D hessian
+ 
+  ! internal
+  double precision, allocatable :: hessian(:,:,:,:)
+  integer                       :: p,q, tmp_p,tmp_q
+  integer                       :: r,s,t,u,v,tmp_r,tmp_s,tmp_t,tmp_u,tmp_v
+  integer                       :: pq,rs,tmp_pq,tmp_rs
+  double precision              :: t1,t2,t3,t4,t5,t6
+  ! hessian  : mo_num 4D double precision matrix containing the hessian before the permutations
+  ! h_tmpr   : mo_num 4D double precision matrix containing the hessian after the permutations
+  ! p,q,r,s  : integer, indexes of the 4D hessian matrix
+  ! t,u,v    : integer, indexes to compute hessian elements
+  ! pq,rs    : integer, indexes for the conversion from 4D to 2D hessian matrix
+  ! t1,t2,t3 : double precision, t3 = t2 - t1, time to compute the hessian 
+
+  ! Funtion 
+  double precision              :: get_two_e_integral
+  ! get_two_e_integral :  double precision function, two e integrals 
+
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo_alpha, one_e_dm_mo_beta : one body density matrix
+  ! two_e_dm_mo : two body density matrix
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(hessian(m,m,m,m))
+
+  !=============
+  ! Calculation
+  !=============
+
+  print*,'---first_hess_list---'
+
+  ! From Anderson et. al. (2014) 
+  ! The Journal of Chemical Physics 141, 244104 (2014); doi: 10.1063/1.4904384
+
+  CALL wall_time(t1)
+
+  ! Initialization
+  hessian = 0d0
+
+  !========================
+  ! First line, first term
+  !========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+            if (q==r) then
+              do u = 1, mo_num
+
+                hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * (  &
+                  mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
+                + mo_one_e_integrals(s,u) * one_e_dm_mo(p,u))
+
+              enddo
+            endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 1 :', t6
+
+  !=========================
+  ! First line, second term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          if (p==s) then
+            do u = 1, mo_num
+
+                  hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                    mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
+                  + mo_one_e_integrals(q,u) * one_e_dm_mo(r,u))
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+  
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 2 :', t6
+
+  !========================
+  ! First line, third term
+  !========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+          - mo_one_e_integrals(s,p) * one_e_dm_mo(r,q)&
+          - mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 3 :', t6
+
+
+  !=========================
+  ! Second line, first term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+           if (q==r) then
+             do t = 1, mo_num
+               do u = 1, mo_num
+                 do v = 1, mo_num
+
+                   hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * (  &
+                     get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
+                   + get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v))
+
+                 enddo
+               enddo
+             enddo
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l2 1 :', t6
+
+  !==========================
+  ! Second line, second term
+  !==========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          if (p==s) then
+            do t = 1, mo_num
+              do u = 1, mo_num
+                do v = 1, mo_num
+
+                  hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                    get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
+                  + get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t))
+
+                enddo
+              enddo
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l2 2 :', t6
+
+  !========================
+  ! Third line, first term
+  !========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          do u = 1, mo_num
+            do v = 1, mo_num
+
+              hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+               + get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
+               + get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v)
+
+            enddo
+          enddo
+
+        enddo
+      enddo
+    enddo
+  enddo
+ 
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l3 1 :', t6
+
+  !=========================
+  ! Third line, second term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          do t = 1, mo_num
+            do u = 1, mo_num
+
+              hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+               - get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
+               - get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
+               - get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
+               - get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
+
+            enddo
+          enddo
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l3 2 :', t6
+
+  CALL wall_time(t2)
+  t3 = t2 -t1
+  print*,'Time to compute the hessian : ', t3
+
+  !==============
+  ! Permutations 
+  !==============
+
+  ! Hessian(p,q,r,s) = P_pq P_rs [ ...]
+  ! => Hessian(p,q,r,s) = (p,q,r,s) - (q,p,r,s) - (p,q,s,r) + (q,p,s,r) 
+
+  do tmp_s = 1, m
+    do tmp_r = 1, m
+      do tmp_q = 1, m
+        do tmp_p = 1, m
+
+          h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s) = (hessian(tmp_p,tmp_q,tmp_r,tmp_s) - hessian(tmp_q,tmp_p,tmp_r,tmp_s) &
+                                           - hessian(tmp_p,tmp_q,tmp_s,tmp_r) + hessian(tmp_q,tmp_p,tmp_s,tmp_r))
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! 4D matrix to 2D matrix
+  !========================
+
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  ! 4D mo_num matrix to 2D n matrix
+  do tmp_pq = 1, tmp_n
+    call vec_to_mat_index(tmp_pq,tmp_p,tmp_q)
+    do tmp_rs = 1, tmp_n
+      call vec_to_mat_index(tmp_rs,tmp_r,tmp_s)
+      H(tmp_pq,tmp_rs) = h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s)   
+    enddo
+  enddo
+
+  ! Display
+  if (debug) then 
+    print*,'2D Hessian matrix'
+    do tmp_pq = 1, tmp_n
+      write(*,'(100(F10.5))') H(tmp_pq,:)
+    enddo 
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(hessian)
+
+  print*,'---End first_hess_list---'
+
+end subroutine

src/mo_optimization/first_hessian_opt.irp.f

@@ -0,0 +1,360 @@
+subroutine first_hessian_opt(n,H,h_tmpr)
+
+  include 'constants.h' 
+
+  implicit none
+
+  !==================================================================
+  ! Compute the hessian of energy with respects to orbital rotations
+  !==================================================================
+
+  !===========
+  ! Variables  
+  !===========
+
+  ! in
+  integer, intent(in)           :: n 
+  !n         : integer, n = mo_num*(mo_num-1)/2
+  
+  ! out
+  double precision, intent(out) :: H(n,n),h_tmpr(mo_num,mo_num,mo_num,mo_num)
+  ! H        : n by n double precision matrix containing the 2D hessian
+ 
+  ! internal
+  double precision, allocatable :: hessian(:,:,:,:)
+  integer                       :: p,q
+  integer                       :: r,s,t,u,v
+  integer                       :: pq,rs
+  double precision              :: t1,t2,t3,t4,t5,t6
+  ! hessian  : mo_num 4D double precision matrix containing the hessian before the permutations
+  ! h_tmpr   : mo_num 4D double precision matrix containing the hessian after the permutations
+  ! p,q,r,s  : integer, indexes of the 4D hessian matrix
+  ! t,u,v    : integer, indexes to compute hessian elements
+  ! pq,rs    : integer, indexes for the conversion from 4D to 2D hessian matrix
+  ! t1,t2,t3 : double precision, t3 = t2 - t1, time to compute the hessian 
+
+  ! Funtion 
+  double precision              :: get_two_e_integral
+  ! get_two_e_integral :  double precision function, two e integrals 
+
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo_alpha, one_e_dm_mo_beta : one body density matrix
+  ! two_e_dm_mo : two body density matrix
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(hessian(mo_num,mo_num,mo_num,mo_num))
+
+  !=============
+  ! Calculation
+  !=============
+
+  if (debug) then
+    print*,'Enter in first_hess'
+  endif
+
+  ! From Anderson et. al. (2014) 
+  ! The Journal of Chemical Physics 141, 244104 (2014); doi: 10.1063/1.4904384
+
+  CALL wall_time(t1)
+
+  ! Initialization
+  hessian = 0d0
+
+  !========================
+  ! First line, first term
+  !========================
+
+  CALL wall_time(t4)
+
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+            if (q==r) then
+              do u = 1, mo_num
+
+                hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * (  &
+                  mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
+                + mo_one_e_integrals(s,u) * one_e_dm_mo(p,u))
+
+              enddo
+            endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 1 :', t6
+
+  !=========================
+  ! First line, second term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          if (p==s) then
+            do u = 1, mo_num
+
+                  hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                    mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
+                  + mo_one_e_integrals(q,u) * one_e_dm_mo(r,u))
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+  
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 2 :', t6
+
+  !========================
+  ! First line, third term
+  !========================
+
+  CALL wall_time(t4)
+
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          hessian(p,q,r,s) = hessian(p,q,r,s) &
+          - mo_one_e_integrals(s,p) * one_e_dm_mo(r,q)&
+          - mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 3 :', t6
+
+
+  !=========================
+  ! Second line, first term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do s = 1, mo_num
+    do r = 1, mo_num
+      do q = 1, mo_num
+        do p = 1, mo_num
+
+!  do p = 1, mo_num
+!    do q = 1, mo_num
+!      do r = 1, mo_num
+!        do s = 1, mo_num
+
+           if (q==r) then
+             do t = 1, mo_num
+               do u = 1, mo_num
+                 do v = 1, mo_num
+
+                   hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * (  &
+                     get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
+                   + get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v))
+
+                 enddo
+               enddo
+             enddo
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l2 1 :', t6
+
+  !==========================
+  ! Second line, second term
+  !==========================
+
+  CALL wall_time(t4)
+
+  do s = 1, mo_num
+    do r = 1, mo_num
+      do q = 1, mo_num
+        do p = 1, mo_num
+
+!  do p = 1, mo_num
+!    do q = 1, mo_num
+!      do r = 1, mo_num
+!        do s = 1, mo_num
+
+          if (p==s) then
+            do t = 1, mo_num
+              do u = 1, mo_num
+                do v = 1, mo_num
+
+                  hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                    get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
+                  + get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t))
+
+                enddo
+              enddo
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l2 2 :', t6
+
+  !========================
+  ! Third line, first term
+  !========================
+
+  CALL wall_time(t4)
+
+      do s = 1, mo_num
+        do r = 1, mo_num
+          do q = 1, mo_num
+            do p = 1, mo_num
+
+!  do p = 1, mo_num
+!    do q = 1, mo_num
+!      do r = 1, mo_num
+!        do s = 1, mo_num
+
+          do u = 1, mo_num
+            do v = 1, mo_num
+
+              hessian(p,q,r,s) = hessian(p,q,r,s) &
+               + get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
+               + get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v)
+
+            enddo
+          enddo
+
+        enddo
+      enddo
+    enddo
+  enddo
+ 
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l3 1 :', t6
+
+  !=========================
+  ! Third line, second term
+  !=========================
+
+  CALL wall_time(t4)
+
+      do s = 1, mo_num
+        do r = 1, mo_num
+          do q = 1, mo_num
+            do p = 1, mo_num
+
+!  do p = 1, mo_num
+!    do q = 1, mo_num
+!      do r = 1, mo_num
+!        do s = 1, mo_num
+
+          do t = 1, mo_num
+            do u = 1, mo_num
+
+              hessian(p,q,r,s) = hessian(p,q,r,s) &
+               - get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
+               - get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
+               - get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
+               - get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
+
+            enddo
+          enddo
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l3 2 :', t6
+
+  CALL wall_time(t2)
+  t3 = t2 -t1
+  print*,'Time to compute the hessian : ', t3
+
+  !==============
+  ! Permutations 
+  !==============
+
+  ! Hessian(p,q,r,s) = P_pq P_rs [ ...]
+  ! => Hessian(p,q,r,s) = (p,q,r,s) - (q,p,r,s) - (p,q,s,r) + (q,p,s,r) 
+
+  do s = 1, mo_num
+    do r = 1, mo_num
+      do q = 1, mo_num
+        do p = 1, mo_num
+
+          h_tmpr(p,q,r,s) = (hessian(p,q,r,s) - hessian(q,p,r,s) - hessian(p,q,s,r) + hessian(q,p,s,r))
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! 4D matrix to 2D matrix
+  !========================
+
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  ! 4D mo_num matrix to 2D n matrix
+  do pq = 1, n
+    call vec_to_mat_index(pq,p,q)
+    do rs = 1, n
+      call vec_to_mat_index(rs,r,s)
+      H(pq,rs) = h_tmpr(p,q,r,s)   
+    enddo
+  enddo
+
+  ! Display
+  if (debug) then 
+    print*,'2D Hessian matrix'
+    do pq = 1, n
+      write(*,'(100(F10.5))') H(pq,:)
+    enddo 
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(hessian)
+
+  if (debug) then
+    print*,'Leave first_hess'
+  endif
+
+end subroutine

src/mo_optimization/gradient_list_opt.irp.f

@@ -0,0 +1,381 @@
+! Gradient
+
+! The gradient of the CI energy with respects to the orbital rotation
+! is:
+! (C-c C-x C-l)
+! $$
+! G(p,q) = \mathcal{P}_{pq} \left[ \sum_r (h_p^r \gamma_r^q - h_r^q \gamma_p^r) +
+! \sum_{rst}(v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs})
+! \right]
+! $$
+
+
+! $$
+! \mathcal{P}_{pq}= 1 - (p \leftrightarrow q)
+! $$
+
+! $$
+! G(p,q) = \left[
+! \sum_r (h_p^r \gamma_r^q - h_r^q \gamma_p^r) +
+! \sum_{rst}(v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs})
+! \right] - 
+! \left[
+! \sum_r (h_q^r \gamma_r^p - h_r^p \gamma_q^r) +
+! \sum_{rst}(v_{qt}^{rs} \Gamma_{rs}^{pt} - v_{rs}^{pt}
+! \Gamma_{qt}^{rs})
+! \right]
+! $$
+
+! Where p,q,r,s,t are general spatial orbitals
+! mo_num : the number of molecular orbitals
+! $$h$$ : One electron integrals
+! $$\gamma$$ : One body density matrix (state average in our case)
+! $$v$$ : Two electron integrals
+! $$\Gamma$$ : Two body density matrice (state average in our case)
+
+! The gradient is a mo_num by mo_num matrix, p,q,r,s,t take all the
+! values between 1 and mo_num (1 and mo_num include).
+
+! To do that we compute $$G(p,q)$$ for all the pairs (p,q).
+
+! Source :
+! Seniority-based coupled cluster theory
+! J. Chem. Phys. 141, 244104 (2014); https://doi.org/10.1063/1.4904384
+! Thomas M. Henderson, Ireneusz W. Bulik, Tamar Stein, and Gustavo
+! E. Scuseria
+
+! *Compute the gradient of energy with respects to orbital rotations*
+
+! Provided:
+! | mo_num                                   | integer          | number of MOs             |
+! | mo_one_e_integrals(mo_num,mo_num)        | double precision | mono_electronic integrals |
+! | one_e_dm_mo(mo_num,mo_num)               | double precision | one e- density matrix     |
+! | two_e_dm_mo(mo_num,mo_num,mo_num,mo_num) | double precision | two e- density matrix     |
+
+! Input:
+! | n | integer | mo_num*(mo_num-1)/2 |
+
+! Output:
+! | v_grad(n) | double precision | the gradient                    |
+! | max_elem  | double precision | maximum element of the gradient |
+
+! Internal:
+! | grad(mo_num,mo_num)                | double precison   | gradient before the tranformation in a vector             |
+! | A((mo_num,mo_num)                  | doubre precision  | gradient after the permutations                           |
+! | norm                               | double precision  | norm of the gradient                                      |
+! | p, q                               | integer           | indexes of the element in the matrix grad                 |
+! | i                                  | integer           | index for the tranformation in a vector                   |
+! | r, s, t                            | integer           | indexes dor the sums                                      |
+! | t1, t2, t3                         | double precision  | t3 = t2 - t1, time to compute the gradient                |
+! | t4, t5, t6                         | double precission | t6 = t5 - t4, time to compute each element                |
+! | tmp_bi_int_3(mo_num,mo_num,mo_num) | double precision  | 3 indexes temporary array for the bi-electronic integrals |
+! | tmp_2rdm_3(mo_num,mo_num,mo_num)   | double precision  | 3 indexes temporary array for the two e- density matrix   |
+! | tmp_accu(mo_num,mo_num)            | double precision  | temporary array                                           |
+
+! Function:
+! | get_two_e_integral | double precision | bi-electronic integrals |
+! | dnrm2              | double precision | (Lapack) norm           |
+
+
+subroutine gradient_list_opt(n,m,list,v_grad,max_elem,norm)
+  use omp_lib
+  include 'constants.h'
+
+  implicit none
+  
+  ! Variables
+  
+  ! in
+  integer, intent(in) :: n,m,list(m)
+  
+  ! out
+  double precision, intent(out) :: v_grad(n), max_elem, norm
+
+  ! internal
+  double precision, allocatable :: grad(:,:),A(:,:)
+  integer                       :: i,p,q,r,s,t, tmp_p, tmp_q, tmp_i
+  double precision              :: t1,t2,t3,t4,t5,t6
+
+  double precision, allocatable :: tmp_accu(:,:), tmp_mo_one_e_integrals(:,:),tmp_one_e_dm_mo(:,:)
+  double precision, allocatable :: tmp_bi_int_3(:,:,:), tmp_2rdm_3(:,:,:)
+
+  ! Functions
+  double precision :: get_two_e_integral, dnrm2
+
+  
+  print*,''
+  print*,'---gradient---'
+
+  ! Allocation of shared arrays
+  allocate(grad(m,m),A(m,m)) 
+  allocate(tmp_mo_one_e_integrals(m,mo_num),tmp_one_e_dm_mo(mo_num,m))
+
+
+  ! Initialization omp
+  call omp_set_max_active_levels(1)
+
+  !$OMP PARALLEL                                                 &
+      !$OMP PRIVATE(                                             &
+      !$OMP   p,q,r,s,t,tmp_p,tmp_q,                             &
+      !$OMP   tmp_accu,tmp_bi_int_3, tmp_2rdm_3)                &
+      !$OMP SHARED(grad, one_e_dm_mo,m,list,mo_num,mo_one_e_integrals, &
+      !$OMP mo_integrals_map,tmp_one_e_dm_mo, tmp_mo_one_e_integrals,t4,t5,t6)                           &
+      !$OMP DEFAULT(SHARED)
+ 
+  ! Allocation of private arrays
+  allocate(tmp_accu(m,m))
+  allocate(tmp_bi_int_3(mo_num,mo_num,m))
+  allocate(tmp_2rdm_3(mo_num,mo_num,m))
+
+! Initialization
+
+!$OMP DO
+do tmp_q = 1, m
+  do tmp_p = 1, m
+    grad(tmp_p,tmp_q) = 0d0
+  enddo
+enddo
+!$OMP END DO
+
+! Term 1
+  
+! Without optimization the term 1 is :
+
+! do p = 1, mo_num
+!   do q = 1, mo_num
+!      do r = 1, mo_num
+!        grad(p,q) = grad(p,q) &
+!                + mo_one_e_integrals(p,r) * one_e_dm_mo(r,q) &
+!                - mo_one_e_integrals(r,q) * one_e_dm_mo(p,r)
+!     enddo
+!   enddo
+! enddo
+   
+! Since the matrix multiplication A.B is defined like :
+! \begin{equation}
+! c_{ij} = \sum_k a_{ik}.b_{kj}
+! \end{equation}
+! The previous equation can be rewritten as a matrix multplication  
+  
+
+!****************
+  ! Opt first term
+  !****************
+
+  !$OMP DO
+  do r = 1, mo_num
+    do tmp_p = 1, m
+      p = list(tmp_p)
+      tmp_mo_one_e_integrals(tmp_p,r) = mo_one_e_integrals(p,r)
+    enddo
+  enddo
+  !$OMP END DO
+
+  !$OMP DO
+  do tmp_q = 1, m
+    q = list(tmp_q) 
+    do r = 1, mo_num
+      tmp_one_e_dm_mo(r,tmp_q) = one_e_dm_mo(r,q)
+    enddo
+  enddo
+  !$OMP END DO
+
+  call dgemm('N','N',m,m,mo_num,1d0,&
+  tmp_mo_one_e_integrals, size(tmp_mo_one_e_integrals,1),&
+  tmp_one_e_dm_mo,size(tmp_one_e_dm_mo,1),0d0,tmp_accu,size(tmp_accu,1))
+
+  !$OMP DO
+  do tmp_q = 1, m
+    do tmp_p = 1, m
+
+      grad(tmp_p,tmp_q) = grad(tmp_p,tmp_q) + (tmp_accu(tmp_p,tmp_q) - tmp_accu(tmp_q,tmp_p))
+
+    enddo
+  enddo 
+  !$OMP END DO
+
+  !$OMP MASTER
+  CALL wall_TIME(t4)
+  !$OMP END MASTER
+
+!  call dgemm('N','N',mo_num,mo_num,mo_num,1d0,mo_one_e_integrals,&
+!  mo_num,one_e_dm_mo,mo_num,0d0,tmp_accu,mo_num)
+!  
+!  !$OMP DO
+!  do q = 1, mo_num
+!    do p = 1, mo_num
+!
+!      grad(p,q) = grad(p,q) + (tmp_accu(p,q) - tmp_accu(q,p))
+!
+!    enddo
+!  enddo 
+!  !$OMP END DO
+  
+  !$OMP MASTER
+  CALL wall_TIME(t5)
+  t6 = t5-t4
+  print*,'Gradient, first term (s) :', t6 
+  !$OMP END MASTER
+
+! Term 2
+ 
+! Without optimization the second term is : 
+
+! do p = 1, mo_num
+!   do q = 1, mo_num 
+!     do r = 1, mo_num
+!       do s = 1, mo_num
+!         do t= 1, mo_num
+
+!         grad(p,q) = grad(p,q) &
+!                 + get_two_e_integral(p,t,r,s,mo_integrals_map) * two_e_dm_mo(r,s,q,t) &
+!                 - get_two_e_integral(r,s,q,t,mo_integrals_map) * two_e_dm_mo(p,t,r,s)
+!        enddo
+!       enddo
+!     enddo
+!   enddo
+! enddo
+
+! Using the bielectronic integral properties :
+! get_two_e_integral(p,t,r,s,mo_integrals_map) = get_two_e_integral(r,s,p,t,mo_integrals_map)
+
+! Using the two body matrix properties :
+! two_e_dm_mo(p,t,r,s) = two_e_dm_mo(r,s,p,t)
+
+! t is one the right, we can put it on the external loop and create 3
+! indexes temporary array 
+! r,s can be seen as one index
+
+! By doing so, a matrix multiplication appears 
+
+
+!*****************
+! Opt second term  
+!*****************
+
+!$OMP MASTER
+CALL wall_TIME(t4)
+!$OMP END MASTER 
+
+!$OMP DO
+do t = 1, mo_num
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do s = 1, mo_num
+      do r = 1, mo_num
+
+        tmp_bi_int_3(r,s,tmp_p) = get_two_e_integral(r,s,p,t,mo_integrals_map)
+
+      enddo
+    enddo
+  enddo
+
+  do tmp_q = 1, m
+    q = list(tmp_q)
+    do s = 1, mo_num
+      do r = 1, mo_num
+
+         tmp_2rdm_3(r,s,tmp_q) = two_e_dm_mo(r,s,q,t)
+
+      enddo
+    enddo
+  enddo
+
+  call dgemm('T','N',m,m,mo_num*mo_num,1d0,tmp_bi_int_3,&
+    mo_num*mo_num,tmp_2rdm_3,mo_num*mo_num,0d0,tmp_accu,size(tmp_accu,1))
+
+  !$OMP CRITICAL   
+  do tmp_q = 1, m
+    do tmp_p = 1, m
+
+      grad(tmp_p,tmp_q) = grad(tmp_p,tmp_q) + tmp_accu(tmp_p,tmp_q) - tmp_accu(tmp_q,tmp_p)
+
+    enddo
+  enddo
+  !$OMP END CRITICAL
+
+enddo
+!$OMP END DO
+
+!$OMP MASTER
+CALL wall_TIME(t5)
+t6 = t5-t4
+print*,'Gradient second term (s) : ', t6
+!$OMP END MASTER
+
+! Deallocation of private arrays
+
+deallocate(tmp_bi_int_3,tmp_2rdm_3,tmp_accu)
+
+!$OMP END PARALLEL
+
+call omp_set_max_active_levels(4)
+
+! Permutation, 2D matrix -> vector, transformation
+! In addition there is a permutation in the gradient formula :
+! \begin{equation}
+! P_{pq} = 1 - (p <-> q) 
+! \end{equation}
+
+! We need a vector to use the gradient. Here the gradient is a
+! antisymetric matrix so we can transform it in a vector of length
+! mo_num*(mo_num-1)/2.
+
+! Here we do these two things at the same time.
+
+
+do i=1,n
+  call vec_to_mat_index(i,p,q)
+  v_grad(i)=(grad(p,q) - grad(q,p))
+enddo  
+
+! Debug, diplay the vector containing the gradient elements 
+if (debug) then  
+  print*,'Vector containing the gradient :'
+  write(*,'(100(F10.5))') v_grad(1:n)
+endif
+
+! Norm of the gradient
+! The norm can be useful.
+
+norm = dnrm2(n,v_grad,1)
+print*, 'Gradient norm : ', norm
+
+! Maximum element in the gradient
+! The maximum element in the gradient is very important for the
+! convergence criterion of the Newton method.
+
+
+! Max element of the gradient
+max_elem = 0d0
+do i = 1, n
+  if (DABS(v_grad(i)) > DABS(max_elem)) then
+    max_elem = v_grad(i)
+  endif
+enddo
+
+print*,'Max element in the gradient :', max_elem  
+
+! Debug, display the matrix containting the gradient elements
+if (debug) then
+  ! Matrix gradient
+  A = 0d0
+  do q=1,m
+    do p=1,m
+      A(p,q) = grad(p,q) - grad(q,p)
+    enddo
+  enddo
+  print*,'Matrix containing the gradient :'
+  do i = 1, m
+    write(*,'(100(F10.5))') A(i,1:m)
+  enddo
+endif
+
+! Deallocation of shared arrays and end
+
+deallocate(grad,A, tmp_mo_one_e_integrals,tmp_one_e_dm_mo)
+
+print*,'---End gradient---'
+
+end subroutine

src/mo_optimization/gradient_opt.irp.f

@@ -0,0 +1,346 @@
+! Gradient
+
+! The gradient of the CI energy with respects to the orbital rotation
+! is:
+! (C-c C-x C-l)
+! $$
+! G(p,q) = \mathcal{P}_{pq} \left[ \sum_r (h_p^r \gamma_r^q - h_r^q \gamma_p^r) +
+! \sum_{rst}(v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs})
+! \right]
+! $$
+
+
+! $$
+! \mathcal{P}_{pq}= 1 - (p \leftrightarrow q)
+! $$
+
+! $$
+! G(p,q) = \left[
+! \sum_r (h_p^r \gamma_r^q - h_r^q \gamma_p^r) +
+! \sum_{rst}(v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs})
+! \right] - 
+! \left[
+! \sum_r (h_q^r \gamma_r^p - h_r^p \gamma_q^r) +
+! \sum_{rst}(v_{qt}^{rs} \Gamma_{rs}^{pt} - v_{rs}^{pt}
+! \Gamma_{qt}^{rs})
+! \right]
+! $$
+
+! Where p,q,r,s,t are general spatial orbitals
+! mo_num : the number of molecular orbitals
+! $$h$$ : One electron integrals
+! $$\gamma$$ : One body density matrix (state average in our case)
+! $$v$$ : Two electron integrals
+! $$\Gamma$$ : Two body density matrice (state average in our case)
+
+! The gradient is a mo_num by mo_num matrix, p,q,r,s,t take all the
+! values between 1 and mo_num (1 and mo_num include).
+
+! To do that we compute $$G(p,q)$$ for all the pairs (p,q).
+
+! Source :
+! Seniority-based coupled cluster theory
+! J. Chem. Phys. 141, 244104 (2014); https://doi.org/10.1063/1.4904384
+! Thomas M. Henderson, Ireneusz W. Bulik, Tamar Stein, and Gustavo
+! E. Scuseria
+
+! *Compute the gradient of energy with respects to orbital rotations*
+
+! Provided:
+! | mo_num                                   | integer          | number of MOs             |
+! | mo_one_e_integrals(mo_num,mo_num)        | double precision | mono_electronic integrals |
+! | one_e_dm_mo(mo_num,mo_num)               | double precision | one e- density matrix     |
+! | two_e_dm_mo(mo_num,mo_num,mo_num,mo_num) | double precision | two e- density matrix     |
+
+! Input:
+! | n | integer | mo_num*(mo_num-1)/2 |
+
+! Output:
+! | v_grad(n) | double precision | the gradient                    |
+! | max_elem  | double precision | maximum element of the gradient |
+
+! Internal:
+! | grad(mo_num,mo_num)                | double precison   | gradient before the tranformation in a vector             |
+! | A((mo_num,mo_num)                  | doubre precision  | gradient after the permutations                           |
+! | norm                               | double precision  | norm of the gradient                                      |
+! | p, q                               | integer           | indexes of the element in the matrix grad                 |
+! | i                                  | integer           | index for the tranformation in a vector                   |
+! | r, s, t                            | integer           | indexes dor the sums                                      |
+! | t1, t2, t3                         | double precision  | t3 = t2 - t1, time to compute the gradient                |
+! | t4, t5, t6                         | double precission | t6 = t5 - t4, time to compute each element                |
+! | tmp_bi_int_3(mo_num,mo_num,mo_num) | double precision  | 3 indexes temporary array for the bi-electronic integrals |
+! | tmp_2rdm_3(mo_num,mo_num,mo_num)   | double precision  | 3 indexes temporary array for the two e- density matrix   |
+! | tmp_accu(mo_num,mo_num)            | double precision  | temporary array                                           |
+
+! Function:
+! | get_two_e_integral | double precision | bi-electronic integrals |
+! | dnrm2              | double precision | (Lapack) norm           |
+
+
+subroutine gradient_opt(n,v_grad,max_elem)
+  use omp_lib
+  include 'constants.h'
+
+  implicit none
+  
+  ! Variables
+  
+  ! in
+  integer, intent(in) :: n
+  
+  ! out
+  double precision, intent(out) :: v_grad(n), max_elem
+
+  ! internal
+  double precision, allocatable :: grad(:,:),A(:,:)
+  double precision              :: norm
+  integer                       :: i,p,q,r,s,t
+  double precision              :: t1,t2,t3,t4,t5,t6
+
+  double precision, allocatable :: tmp_accu(:,:)
+  double precision, allocatable :: tmp_bi_int_3(:,:,:), tmp_2rdm_3(:,:,:)
+
+  ! Functions
+  double precision :: get_two_e_integral, dnrm2
+
+  
+  print*,''
+  print*,'---gradient---'
+
+  ! Allocation of shared arrays
+  allocate(grad(mo_num,mo_num),A(mo_num,mo_num)) 
+
+  ! Initialization omp
+  call omp_set_max_active_levels(1)
+
+  !$OMP PARALLEL                                                 &
+      !$OMP PRIVATE(                                             &
+      !$OMP   p,q,r,s,t,                                         &
+      !$OMP   tmp_accu, tmp_bi_int_3, tmp_2rdm_3)                &
+      !$OMP SHARED(grad, one_e_dm_mo, mo_num,mo_one_e_integrals, &
+      !$OMP mo_integrals_map,t4,t5,t6)                           &
+      !$OMP DEFAULT(SHARED)
+ 
+  ! Allocation of private arrays
+  allocate(tmp_accu(mo_num,mo_num))
+  allocate(tmp_bi_int_3(mo_num,mo_num,mo_num))
+  allocate(tmp_2rdm_3(mo_num,mo_num,mo_num))
+
+! Initialization
+
+!$OMP DO
+do q = 1, mo_num
+  do p = 1,mo_num
+    grad(p,q) = 0d0
+  enddo
+enddo
+!$OMP END DO
+
+! Term 1
+  
+! Without optimization the term 1 is :
+
+! do p = 1, mo_num
+!   do q = 1, mo_num
+!      do r = 1, mo_num
+!        grad(p,q) = grad(p,q) &
+!                + mo_one_e_integrals(p,r) * one_e_dm_mo(r,q) &
+!                - mo_one_e_integrals(r,q) * one_e_dm_mo(p,r)
+!     enddo
+!   enddo
+! enddo
+   
+! Since the matrix multiplication A.B is defined like :
+! \begin{equation}
+! c_{ij} = \sum_k a_{ik}.b_{kj}
+! \end{equation}
+! The previous equation can be rewritten as a matrix multplication  
+  
+
+!****************
+! Opt first term
+!****************
+
+!$OMP MASTER
+CALL wall_TIME(t4)
+!$OMP END MASTER
+
+call dgemm('N','N',mo_num,mo_num,mo_num,1d0,mo_one_e_integrals,&
+mo_num,one_e_dm_mo,mo_num,0d0,tmp_accu,mo_num)
+
+!$OMP DO
+do q = 1, mo_num
+  do p = 1, mo_num
+
+    grad(p,q) = grad(p,q) + (tmp_accu(p,q) - tmp_accu(q,p))
+
+  enddo
+enddo 
+!$OMP END DO
+
+!$OMP MASTER
+CALL wall_TIME(t5)
+t6 = t5-t4
+print*,'Gradient, first term (s) :', t6 
+!$OMP END MASTER
+
+! Term 2
+ 
+! Without optimization the second term is : 
+
+! do p = 1, mo_num
+!   do q = 1, mo_num 
+!     do r = 1, mo_num
+!       do s = 1, mo_num
+!         do t= 1, mo_num
+
+!         grad(p,q) = grad(p,q) &
+!                 + get_two_e_integral(p,t,r,s,mo_integrals_map) * two_e_dm_mo(r,s,q,t) &
+!                 - get_two_e_integral(r,s,q,t,mo_integrals_map) * two_e_dm_mo(p,t,r,s)
+!        enddo
+!       enddo
+!     enddo
+!   enddo
+! enddo
+
+! Using the bielectronic integral properties :
+! get_two_e_integral(p,t,r,s,mo_integrals_map) = get_two_e_integral(r,s,p,t,mo_integrals_map)
+
+! Using the two body matrix properties :
+! two_e_dm_mo(p,t,r,s) = two_e_dm_mo(r,s,p,t)
+
+! t is one the right, we can put it on the external loop and create 3
+! indexes temporary array 
+! r,s can be seen as one index
+
+! By doing so, a matrix multiplication appears 
+
+
+!*****************
+! Opt second term  
+!*****************
+
+!$OMP MASTER
+CALL wall_TIME(t4)
+!$OMP END MASTER 
+
+!$OMP DO
+do t = 1, mo_num
+
+  do p = 1, mo_num
+    do s = 1, mo_num
+      do r = 1, mo_num
+
+        tmp_bi_int_3(r,s,p) = get_two_e_integral(r,s,p,t,mo_integrals_map)
+
+      enddo
+    enddo
+  enddo
+
+  do q = 1, mo_num
+    do s = 1, mo_num
+      do r = 1, mo_num
+
+         tmp_2rdm_3(r,s,q) = two_e_dm_mo(r,s,q,t)
+
+      enddo
+    enddo
+  enddo
+
+  call dgemm('T','N',mo_num,mo_num,mo_num*mo_num,1d0,tmp_bi_int_3,&
+    mo_num*mo_num,tmp_2rdm_3,mo_num*mo_num,0d0,tmp_accu,mo_num)
+
+  !$OMP CRITICAL   
+  do q = 1, mo_num
+    do p = 1, mo_num
+
+      grad(p,q) = grad(p,q) + tmp_accu(p,q) - tmp_accu(q,p)
+
+    enddo
+  enddo
+  !$OMP END CRITICAL
+
+enddo
+!$OMP END DO
+
+!$OMP MASTER
+CALL wall_TIME(t5)
+t6 = t5-t4
+print*,'Gradient second term (s) : ', t6
+!$OMP END MASTER
+
+! Deallocation of private arrays
+
+deallocate(tmp_bi_int_3,tmp_2rdm_3,tmp_accu)
+
+!$OMP END PARALLEL
+
+call omp_set_max_active_levels(4)
+
+! Permutation, 2D matrix -> vector, transformation
+! In addition there is a permutation in the gradient formula :
+! \begin{equation}
+! P_{pq} = 1 - (p <-> q) 
+! \end{equation}
+
+! We need a vector to use the gradient. Here the gradient is a
+! antisymetric matrix so we can transform it in a vector of length
+! mo_num*(mo_num-1)/2.
+
+! Here we do these two things at the same time.
+
+
+do i=1,n
+  call vec_to_mat_index(i,p,q)
+  v_grad(i)=(grad(p,q) - grad(q,p))
+enddo  
+
+! Debug, diplay the vector containing the gradient elements 
+if (debug) then  
+  print*,'Vector containing the gradient :'
+  write(*,'(100(F10.5))') v_grad(1:n)
+endif
+
+! Norm of the gradient
+! The norm can be useful.
+
+norm = dnrm2(n,v_grad,1)
+print*, 'Gradient norm : ', norm
+
+! Maximum element in the gradient
+! The maximum element in the gradient is very important for the
+! convergence criterion of the Newton method.
+
+
+! Max element of the gradient
+max_elem = 0d0
+do i = 1, n
+  if (ABS(v_grad(i)) > ABS(max_elem)) then
+    max_elem = v_grad(i)
+  endif
+enddo
+
+print*,'Max element in the gradient :', max_elem  
+
+! Debug, display the matrix containting the gradient elements
+if (debug) then
+  ! Matrix gradient
+  A = 0d0
+  do q=1,mo_num
+    do p=1,mo_num
+      A(p,q) = grad(p,q) - grad(q,p)
+    enddo
+  enddo
+  print*,'Matrix containing the gradient :'
+  do i = 1, mo_num
+    write(*,'(100(F10.5))') A(i,1:mo_num)
+  enddo
+endif
+
+! Deallocation of shared arrays and end
+
+deallocate(grad,A)
+
+print*,'---End gradient---'
+
+end subroutine

src/mo_optimization/my_providers.irp.f

@@ -0,0 +1,141 @@
+! Dimensions of MOs
+
+
+BEGIN_PROVIDER [ integer, n_mo_dim ]
+  implicit none
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim = mo_num*(mo_num-1)/2
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_mo_dim_core ]
+  implicit none 
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of core MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim_core = dim_list_core_orb*(dim_list_core_orb-1)/2
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_mo_dim_act ]
+  implicit none
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of active MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim_act = dim_list_act_orb*(dim_list_act_orb-1)/2
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_mo_dim_inact ]
+  implicit none 
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of inactive MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim_inact = dim_list_inact_orb*(dim_list_inact_orb-1)/2
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, n_mo_dim_virt ]
+  implicit none 
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of virtual MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim_virt = dim_list_virt_orb*(dim_list_virt_orb-1)/2
+
+END_PROVIDER
+
+! Energies/criterions
+
+BEGIN_PROVIDER [ double precision, my_st_av_energy ]
+  implicit none
+  BEGIN_DOC
+  ! State average CI energy
+  END_DOC
+
+  !call update_st_av_ci_energy(my_st_av_energy)
+  call state_average_energy(my_st_av_energy)
+
+END_PROVIDER
+
+! With all the MOs
+
+BEGIN_PROVIDER [ double precision, my_gradient_opt, (n_mo_dim) ]
+&BEGIN_PROVIDER [ double precision, my_CC1_opt ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the energy with respect to the MO rotations, for all the MOs.
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision :: norm_grad
+
+  PROVIDE mo_two_e_integrals_in_map
+
+  call gradient_opt(n_mo_dim, my_gradient_opt, my_CC1_opt, norm_grad)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, my_hessian_opt, (n_mo_dim, n_mo_dim) ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the energy with respect to the MO rotations, for all the MOs.
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision, allocatable :: h_f(:,:,:,:)
+
+  PROVIDE mo_two_e_integrals_in_map
+
+  allocate(h_f(mo_num, mo_num, mo_num, mo_num))
+
+  call hessian_list_opt(n_mo_dim, my_hessian_opt, h_f)
+
+END_PROVIDER
+
+! With the list of active MOs
+! Can be generalized to any mo_class by changing the list/dimension
+
+BEGIN_PROVIDER [ double precision, my_gradient_list_opt, (n_mo_dim_act) ]
+&BEGIN_PROVIDER [ double precision, my_CC2_opt ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the energy with respect to the MO rotations, only for the active MOs !
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision :: norm_grad
+
+  PROVIDE mo_two_e_integrals_in_map !one_e_dm_mo two_e_dm_mo mo_one_e_integrals 
+
+  call gradient_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_gradient_list_opt, my_CC2_opt, norm_grad)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, my_hessian_list_opt, (n_mo_dim_act, n_mo_dim_act) ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the energy with respect to the MO rotations, only for the active MOs !
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision, allocatable :: h_f(:,:,:,:)
+
+  PROVIDE mo_two_e_integrals_in_map
+
+  allocate(h_f(dim_list_act_orb, dim_list_act_orb, dim_list_act_orb, dim_list_act_orb))
+
+  call hessian_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_hessian_list_opt, h_f)
+
+END_PROVIDER

src/mo_optimization/orb_opt.irp.f

@@ -0,0 +1,22 @@
+! Orbital optimization program
+
+! This is an optimization program for molecular orbitals. It produces
+! orbital rotations in order to lower the energy of a truncated wave
+! function.  
+! This program just optimize the orbitals for a fixed number of
+! determinants. This optimization process must be repeated for different
+! number of determinants.
+
+
+
+
+! Main program : orb_opt_trust
+
+
+program orb_opt
+  read_wf = .true. ! must be True for the orbital optimization !!!
+  TOUCH read_wf
+  io_mo_two_e_integrals = 'None'
+  TOUCH io_mo_two_e_integrals
+  call run_orb_opt_trust_v2
+end

src/mo_optimization/org/TANGLE_org_mode.sh

@@ -0,0 +1,7 @@
+#!/bin/sh   
+
+list='ls *.org'
+for element in $list    
+do   
+		emacs --batch $element -f org-babel-tangle
+done

src/mo_optimization/org/TODO.org

@@ -0,0 +1,17 @@
+TODO:
+** TODO Keep under surveillance the performance of rotation matrix
+- is the fix ok ?
+** DONE Provider state_average_weight
+** DONE Diagonal hessian for orbital optimization with a list of MOs
+** DONE Something to force the step cancellation if R.R^T > treshold
+** TODO Iterative method to compute the rotation matrix
+- doesn't work actually
+** DONE Test trust region with polynomial functions
+** DONE Optimization/Localization program using the template
+** DONE Correction OMP hessian shared/private arrays
+** DONE State average energy
+** DONE Correction of Rho
+** TODO Check the PROVIDE/FREE/TOUCH
+** TODO research of lambda without the power 2
+** DONE Clean the OMP sections
+

src/mo_optimization/org/debug_gradient_list_opt.org

@@ -0,0 +1,79 @@
+* Debug the gradient
+
+*Program to check the gradient*
+
+The program compares the result of the first and last code for the
+gradient.
+
+Provided:
+| mo_num | integer | number of MOs |
+
+Internal:
+| n          | integer          | number of orbitals pairs (p,q) p<q |
+| v_grad(n)  | double precision | Original gradient                  |
+| v_grad2(n) | double precision | Gradient                           |
+| i          | integer          | index                              |
+| threshold  | double precision | threshold for the errors           |
+| max_error  | double precision | maximal error in the gradient      |
+| nb_error   | integer          | number of error in the gradient    |
+
+#+BEGIN_SRC f90 :comments org :tangle debug_gradient_list_opt.irp.f
+program debug_gradient_list
+  
+  implicit none
+
+  ! Variables
+
+  double precision, allocatable :: v_grad(:), v_grad2(:)
+  integer                       :: n,m
+  integer                       :: i
+  double precision              :: threshold
+  double precision              :: max_error, max_elem, norm
+  integer                       :: nb_error
+  
+  m = dim_list_act_orb
+  ! Definition of n  
+  n = m*(m-1)/2
+
+  PROVIDE mo_two_e_integrals_in_map ! Vérifier pour suppression
+
+  ! Allocation
+  allocate(v_grad(n), v_grad2(n))
+
+  ! Calculation
+
+  call diagonalize_ci ! Vérifier pour suppression
+
+  ! Gradient  
+  call gradient_list_opt(n,m,list_act,v_grad,max_elem,norm)
+  call first_gradient_list_opt(n,m,list_act,v_grad2)
+  
+  
+  v_grad = v_grad - v_grad2
+  nb_error = 0
+  max_error = 0d0 
+  threshold = 1d-12 
+
+  do i = 1, n
+    if (ABS(v_grad(i)) > threshold) then
+       print*,i,v_grad(i)
+       nb_error = nb_error + 1
+
+       if (ABS(v_grad(i)) > max_error) then
+         max_error = v_grad(i)
+       endif
+
+    endif
+  enddo
+ 
+  print*,''
+  print*,'Check the gradient' 
+  print*,'Threshold:', threshold
+  print*,'Nb error:', nb_error
+  print*,'Max error:', max_error
+
+  ! Deallocation
+  deallocate(v_grad,v_grad2)
+
+end program
+#+END_SRC

src/mo_optimization/org/debug_gradient_opt.org

@@ -0,0 +1,77 @@
+* Debug the gradient
+
+*Program to check the gradient*
+
+The program compares the result of the first and last code for the
+gradient.
+
+Provided:
+| mo_num | integer | number of MOs |
+
+Internal:
+| n          | integer          | number of orbitals pairs (p,q) p<q |
+| v_grad(n)  | double precision | Original gradient                  |
+| v_grad2(n) | double precision | Gradient                           |
+| i          | integer          | index                              |
+| threshold  | double precision | threshold for the errors           |
+| max_error  | double precision | maximal error in the gradient      |
+| nb_error   | integer          | number of error in the gradient    |
+
+#+BEGIN_SRC f90 :comments org :tangle debug_gradient_opt.irp.f
+program debug_gradient
+  
+  implicit none
+
+  ! Variables
+
+  double precision, allocatable :: v_grad(:), v_grad2(:)
+  integer                       :: n
+  integer                       :: i
+  double precision              :: threshold
+  double precision              :: max_error, max_elem
+  integer                       :: nb_error
+  
+  ! Definition of n  
+  n = mo_num*(mo_num-1)/2
+
+  PROVIDE mo_two_e_integrals_in_map ! Vérifier pour suppression
+
+  ! Allocation
+  allocate(v_grad(n), v_grad2(n))
+
+  ! Calculation
+
+  call diagonalize_ci ! Vérifier pour suppression
+
+  ! Gradient  
+  call first_gradient_opt(n,v_grad)
+  call gradient_opt(n,v_grad2,max_elem)
+  
+  v_grad = v_grad - v_grad2
+  nb_error = 0
+  max_error = 0d0 
+  threshold = 1d-12 
+
+  do i = 1, n
+    if (ABS(v_grad(i)) > threshold) then
+       print*,v_grad(i)
+       nb_error = nb_error + 1
+
+       if (ABS(v_grad(i)) > max_error) then
+         max_error = v_grad(i)
+       endif
+
+    endif
+  enddo
+ 
+  print*,''
+  print*,'Check the gradient' 
+  print*,'Threshold :', threshold
+  print*,'Nb error :', nb_error
+  print*,'Max error :', max_error
+
+  ! Deallocation
+  deallocate(v_grad,v_grad2)
+
+end program
+#+END_SRC

src/mo_optimization/org/debug_hessian_list_opt.org

@@ -0,0 +1,148 @@
+* Debug the hessian
+
+*Program to check the hessian matrix*
+
+The program compares the result of the first and last code for the
+hessian. First of all the 4D hessian and after the 2D hessian.
+
+Provided:
+| mo_num                     | integer | number of MOs                        |
+| optimization_method        | string  | Method for the orbital optimization: |
+|                            |         | - 'full' -> full hessian             |
+|                            |         | - 'diag' -> diagonal hessian         |
+| dim_list_act_orb           | integer | number of active MOs                 |
+| list_act(dim_list_act_orb) | integer | list of the actives MOs              |
+|                            |         |                                      |
+
+Internal:
+| m                                 | integer          | number of MOs in the list          |
+|                                   |                  | (active MOs)                       |
+| n                                 | integer          | number of orbitals pairs (p,q) p<q |
+|                                   |                  | n = m*(m-1)/2                      |
+| H(n,n)                            | double precision | Original hessian matrix (2D)       |
+| H2(n,n)                           | double precision | Hessian matrix (2D)                |
+| h_f(mo_num,mo_num,mo_num,mo_num)  | double precision | Original hessian matrix (4D)       |
+| h_f2(mo_num,mo_num,mo_num,mo_num) | double precision | Hessian matrix (4D)                |
+| i,j,p,q,k                         | integer          | indexes                            |
+| threshold                         | double precision | threshold for the errors           |
+| max_error                         | double precision | maximal error in the 4D hessian    |
+| max_error_H                       | double precision | maximal error in the 2D hessian    |
+| nb_error                          | integer          | number of errors in the 4D hessian |
+| nb_error_H                        | integer          | number of errors in the 2D hessian |
+
+#+BEGIN_SRC f90 :comments org :tangle debug_hessian_list_opt.irp.f
+program debug_hessian_list_opt
+
+  implicit none
+
+  ! Variables
+
+  double precision, allocatable :: H(:,:),H2(:,:), h_f(:,:,:,:), h_f2(:,:,:,:)
+  integer                       :: n,m
+  integer                       :: i,j,k,l
+  double precision              :: max_error, max_error_H
+  integer                       :: nb_error, nb_error_H
+  double precision              :: threshold
+  
+  m = dim_list_act_orb !mo_num
+
+  ! Definition of n  
+  n = m*(m-1)/2
+
+  PROVIDE mo_two_e_integrals_in_map ! Vérifier pour suppression
+
+  ! Hessian 
+  if (optimization_method == 'full') then
+    print*,'Use the full hessian matrix'
+    allocate(H(n,n),H2(n,n))  
+    allocate(h_f(m,m,m,m),h_f2(m,m,m,m))
+
+    call hessian_list_opt(n,m,list_act,H,h_f)
+    call first_hessian_list_opt(n,m,list_act,H2,h_f2)
+    !call hessian_opt(n,H2,h_f2)
+
+    ! Difference
+    h_f = h_f - h_f2
+    H = H - H2
+    max_error = 0d0
+    nb_error = 0    
+    threshold = 1d-12
+
+    do l = 1, m
+      do k= 1, m
+        do j = 1, m
+          do i = 1, m
+            if (ABS(h_f(i,j,k,l)) > threshold) then
+              print*,h_f(i,j,k,l)
+              nb_error = nb_error + 1
+              if (ABS(h_f(i,j,k,l)) > ABS(max_error)) then
+                max_error = h_f(i,j,k,l)
+              endif
+            endif
+          enddo
+        enddo
+      enddo
+    enddo
+
+    max_error_H = 0d0
+    nb_error_H = 0
+
+    do j = 1, n
+      do i = 1, n
+        if (ABS(H(i,j)) > threshold) then
+          print*, H(i,j)
+          nb_error_H = nb_error_H + 1
+
+          if (ABS(H(i,j)) > ABS(max_error_H)) then
+            max_error_H = H(i,j)
+          endif
+
+        endif
+      enddo
+    enddo 
+
+    ! Deallocation
+    deallocate(H, H2, h_f, h_f2)
+
+  else
+
+    print*, 'Use the diagonal hessian matrix'
+    allocate(H(n,1),H2(n,1))
+    call diag_hessian_list_opt(n,m,list_act,H)
+    call first_diag_hessian_list_opt(n,m,list_act,H2)
+    
+    H = H - H2
+  
+    max_error_H = 0d0
+    nb_error_H = 0
+ 
+    do i = 1, n
+      if (ABS(H(i,1)) > threshold) then
+        print*, H(i,1)
+        nb_error_H = nb_error_H + 1
+ 
+        if (ABS(H(i,1)) > ABS(max_error_H)) then
+          max_error_H = H(i,1)
+        endif
+ 
+      endif
+    enddo
+   
+ endif
+  
+  print*,''
+  if (optimization_method == 'full') then
+    print*,'Check of the full hessian'
+    print*,'Threshold:', threshold
+    print*,'Nb error:', nb_error
+    print*,'Max error:', max_error
+    print*,''
+  else
+    print*,'Check of the diagonal hessian'
+  endif
+   
+  print*,'Nb error_H:', nb_error_H
+  print*,'Max error_H:', max_error_H
+ 
+end program
+#+END_SRC

src/mo_optimization/org/debug_hessian_opt.org

@@ -0,0 +1,172 @@
+* Debug the hessian
+
+*Program to check the hessian matrix*
+
+The program compares the result of the first and last code for the
+hessian. First of all the 4D hessian and after the 2D hessian.
+
+Provided:
+| mo_num | integer | number of MOs |
+
+Internal:
+| n                                 | integer          | number of orbitals pairs (p,q) p<q |
+| H(n,n)                            | double precision | Original hessian matrix (2D)       |
+| H2(n,n)                           | double precision | Hessian matrix (2D)                |
+| h_f(mo_num,mo_num,mo_num,mo_num)  | double precision | Original hessian matrix (4D)       |
+| h_f2(mo_num,mo_num,mo_num,mo_num) | double precision | Hessian matrix (4D)                |
+| method                            | integer          | - 1: full hessian                  |
+|                                   |                  | - 2: diagonal hessian              |
+| i,j,p,q,k                         | integer          | indexes                            |
+| threshold                         | double precision | threshold for the errors           |
+| max_error                         | double precision | maximal error in the 4D hessian    |
+| max_error_H                       | double precision | maximal error in the 2D hessian    |
+| nb_error                          | integer          | number of errors in the 4D hessian |
+| nb_error_H                        | integer          | number of errors in the 2D hessian |
+
+#+BEGIN_SRC f90 :comments org :tangle debug_hessian_opt.irp.f
+program debug_hessian
+
+  implicit none
+
+  ! Variables
+
+  double precision, allocatable :: H(:,:),H2(:,:), h_f(:,:,:,:), h_f2(:,:,:,:)
+  integer                       :: n
+  integer                       :: i,j,k,l
+  double precision              :: max_error, max_error_H
+  integer                       :: nb_error, nb_error_H
+  double precision              :: threshold
+  
+  ! Definition of n  
+  n = mo_num*(mo_num-1)/2
+
+  PROVIDE mo_two_e_integrals_in_map ! Vérifier pour suppression
+
+  ! Allocation
+  allocate(H(n,n),H2(n,n))  
+  allocate(h_f(mo_num,mo_num,mo_num,mo_num),h_f2(mo_num,mo_num,mo_num,mo_num))
+
+  ! Calculation
+  
+  ! Hessian 
+  if (optimization_method == 'full') then 
+
+    print*,'Use the full hessian matrix'
+    call hessian_opt(n,H,h_f)
+    call first_hessian_opt(n,H2,h_f2)
+
+    ! Difference
+    h_f = h_f - h_f2
+    H = H - H2
+    max_error = 0d0
+    nb_error = 0    
+    threshold = 1d-12
+
+    do l = 1, mo_num
+      do k= 1, mo_num
+        do j = 1, mo_num
+          do i = 1, mo_num
+            if (ABS(h_f(i,j,k,l)) > threshold) then
+              print*,h_f(i,j,k,l)
+              nb_error = nb_error + 1
+              if (ABS(h_f(i,j,k,l)) > ABS(max_error)) then
+                max_error = h_f(i,j,k,l)
+              endif
+            endif
+          enddo
+        enddo
+      enddo
+    enddo
+
+   max_error_H = 0d0
+   nb_error_H = 0
+
+   do j = 1, n
+     do i = 1, n
+       if (ABS(H(i,j)) > threshold) then
+         print*, H(i,j)
+         nb_error_H = nb_error_H + 1
+
+         if (ABS(H(i,j)) > ABS(max_error_H)) then
+           max_error_H = H(i,j)
+         endif
+
+       endif
+     enddo
+   enddo 
+
+  elseif (optimization_method == 'diag') then
+
+    print*, 'Use the diagonal hessian matrix'
+    call diag_hessian_opt(n,H,h_f)
+    call first_diag_hessian_opt(n,H2,h_f2)
+
+    h_f = h_f - h_f2
+    max_error = 0d0
+    nb_error = 0
+    threshold = 1d-12
+
+    do l = 1, mo_num
+      do k = 1, mo_num
+        do j = 1, mo_num
+          do i = 1, mo_num
+
+            if (ABS(h_f(i,j,k,l)) > threshold) then
+
+              print*,h_f(i,j,k,l)
+              nb_error = nb_error + 1
+
+              if (ABS(h_f(i,j,k,l)) > ABS(max_error)) then
+                max_error = h_f(i,j,k,l)
+              endif
+
+            endif
+
+          enddo
+        enddo
+      enddo
+    enddo
+
+    h=H-H2
+  
+    max_error_H = 0d0
+    nb_error_H = 0
+ 
+    do j = 1, n
+      do i = 1, n
+        if (ABS(H(i,j)) > threshold) then
+          print*, H(i,j)
+          nb_error_H = nb_error_H + 1
+ 
+          if (ABS(H(i,j)) > ABS(max_error_H)) then
+            max_error_H = H(i,j)
+          endif
+ 
+        endif
+      enddo
+    enddo
+  
+  else
+    print*,'Unknown optimization_method, please select full, diag'
+    call abort
+  endif
+  
+  print*,''
+  if (optimization_method == 'full') then
+    print*,'Check the full hessian'
+  else
+    print*,'Check the diagonal hessian'
+  endif
+   
+  print*,'Threshold :', threshold
+  print*,'Nb error :', nb_error
+  print*,'Max error :', max_error
+  print*,''
+  print*,'Nb error_H :', nb_error_H
+  print*,'Max error_H :', max_error_H
+ 
+  ! Deallocation
+  deallocate(H,H2,h_f,h_f2)
+
+end program
+#+END_SRC

src/mo_optimization/org/diagonalization_hessian.org

@@ -0,0 +1,138 @@
+* Diagonalization of the hessian
+
+Just a matrix diagonalization using Lapack
+
+Input:
+| n      | integer          | mo_num*(mo_num-1)/2 |
+| H(n,n) | double precision | hessian             |
+
+Output:
+| e_val(n) | double precision | eigenvalues of the hessian  |
+| w(n,n)   | double precision | eigenvectors of the hessian |
+
+Internal:
+| nb_negative_nv | integer          | number of negative eigenvalues                  |
+| lwork          | integer          | for Lapack                                      |
+| work(lwork,n)  | double precision | temporary array for Lapack                      |
+| info           | integer          | if 0 -> ok, else problem in the diagonalization |
+| i,j            | integer          | dummy indexes                                   |
+
+#+BEGIN_SRC f90 :comments org :tangle diagonalization_hessian.irp.f
+subroutine diagonalization_hessian(n,H,e_val,w)
+
+  include 'constants.h'
+
+  implicit none
+
+  ! Variables
+
+  ! in
+  integer, intent(in) :: n
+  double precision, intent(in) :: H(n,n)
+
+  ! out
+  double precision, intent(out) :: e_val(n), w(n,n)
+
+  ! internal
+  double precision, allocatable :: work(:,:)
+  integer, allocatable          :: key(:)
+  integer                       :: info,lwork
+  integer                       :: i,j
+  integer                       :: nb_negative_vp
+  double precision              :: t1,t2,t3,max_elem
+
+  print*,''
+  print*,'---Diagonalization_hessian---'
+
+  call wall_time(t1)
+
+  if (optimization_method == 'full') then
+    ! Allocation
+    ! For Lapack
+    lwork=3*n-1
+
+    allocate(work(lwork,n))
+
+    ! Calculation
+
+    ! Copy the hessian matrix, the eigenvectors will be store in W
+    W=H
+
+    ! Diagonalization of the hessian
+    call dsyev('V','U',n,W,size(W,1),e_val,work,lwork,info)
+
+    if (info /= 0) then
+        print*, 'Error diagonalization : diagonalization_hessian'
+        print*, 'info = ', info
+        call ABORT
+    endif
+
+    if (debug) then
+      print *, 'vp Hess:'
+      write(*,'(100(F10.5))')  real(e_val(:))
+    endif
+
+    ! Number of negative eigenvalues
+    max_elem = 0d0
+    nb_negative_vp = 0
+    do i = 1, n
+      if (e_val(i) < 0d0) then
+        nb_negative_vp = nb_negative_vp + 1
+        if (e_val(i) < max_elem) then
+          max_elem = e_val(i)
+        endif
+        !print*,'e_val < 0 :', e_val(i)
+      endif
+    enddo
+    print*,'Number of negative eigenvalues:', nb_negative_vp
+    print*,'Lowest eigenvalue:',max_elem
+
+    !nb_negative_vp = 0
+    !do i = 1, n
+    !  if (e_val(i) < -thresh_eig) then
+    !    nb_negative_vp = nb_negative_vp + 1
+    !  endif
+    !enddo
+    !print*,'Number of negative eigenvalues <', -thresh_eig,':', nb_negative_vp  
+
+    ! Deallocation
+    deallocate(work)
+
+  elseif (optimization_method == 'diag') then
+    ! Diagonalization of the diagonal hessian by hands
+    allocate(key(n))
+    
+    do i = 1, n
+      e_val(i) = H(i,i)
+    enddo
+
+    ! Key list for dsort
+    do i = 1, n 
+      key(i) = i
+    enddo
+
+    ! Sort of the eigenvalues
+    call dsort(e_val, key, n)
+
+    ! Eigenvectors
+    W = 0d0
+    do i = 1, n
+      j = key(i)
+      W(j,i) = 1d0
+    enddo
+
+    deallocate(key)
+  else
+    print*,'Diagonalization_hessian, abort'
+    call abort
+  endif
+
+  call wall_time(t2)
+  t3 = t2 - t1
+  print*,'Time in diagonalization_hessian:', t3
+
+  print*,'---End diagonalization_hessian---'
+
+end subroutine
+#+END_SRC
+

src/mo_optimization/org/first_diagonal_hessian_list_opt.org

@@ -0,0 +1,376 @@
+* First diagonal hessian
+
+#+BEGIN_SRC f90 :comments :tangle first_diagonal_hessian_list_opt.irp.f
+subroutine first_diag_hessian_list_opt(tmp_n,m,list,H)!, h_tmpr)
+
+  include 'constants.h' 
+
+  implicit none
+
+  !===========================================================================
+  ! Compute the diagonal hessian of energy with respects to orbital rotations
+  !===========================================================================
+
+  !===========
+  ! Variables 
+  !===========
+ 
+  ! in
+  integer, intent(in)           :: tmp_n, m, list(m) 
+  ! tmp_n        : integer, tmp_n = m*(m-1)/2
+ 
+  ! out
+  double precision, intent(out) :: H(tmp_n)!, h_tmpr(m,m,m,m)
+  ! H        : n by n double precision matrix containing the 2D hessian
+  
+  ! internal
+  double precision, allocatable :: hessian(:,:,:,:), tmp(:,:),h_tmpr(:,:,:,:)
+  integer                       :: p,q, tmp_p,tmp_q
+  integer                       :: r,s,t,u,v,tmp_r,tmp_s,tmp_t,tmp_u,tmp_v
+  integer                       :: pq,rs,tmp_pq,tmp_rs
+  double precision              :: t1,t2,t3
+  ! hessian  : mo_num 4D double precision matrix containing the hessian before the permutations
+  ! h_tmpr   : mo_num 4D double precision matrix containing the hessian after the permutations
+  ! p,q,r,s  : integer, indexes of the 4D hessian matrix
+  ! t,u,v    : integer, indexes to compute hessian elements
+  ! pq,rs    : integer, indexes for the conversion from 4D to 2D hessian matrix
+  ! t1,t2,t3 : double precision, t3 = t2 - t1, time to compute the hessian 
+  
+  ! Function
+  double precision :: get_two_e_integral
+  ! get_two_e_integral : double precision function, two e integrals
+ 
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo_alpha, one_e_dm_mo_beta : one body density matrix
+  ! two_e_dm_mo : two body density matrix
+
+  print*,'---first_diag_hess_list---'
+
+  !============
+  ! Allocation
+  !============
+  
+  allocate(hessian(m,m,m,m),tmp(tmp_n,tmp_n),h_tmpr(mo_num,mo_num,mo_num,mo_num))
+
+  !=============
+  ! Calculation
+  !=============
+
+  ! From Anderson et. al. (2014) 
+  ! The Journal of Chemical Physics 141, 244104 (2014); doi: 10.1063/1.4904384
+
+  ! LaTeX formula :
+
+  !\begin{align*}
+  !H_{pq,rs} &= \dfrac{\partial^2 E(x)}{\partial x_{pq}^2} \\
+  !&= \mathcal{P}_{pq} \mathcal{P}_{rs} [ \frac{1}{2} \sum_u [\delta_{qr}(h_p^u \gamma_u^s + h_u^s \gamma_p^u) 
+  !+ \delta_{ps}(h_r^u \gamma_u^q + h_u^q \gamma_u^r)]
+  !-(h_p^s \gamma_r^q + h_r^q \gamma_p^s) \\
+  !&+ \frac{1}{2} \sum_{tuv} [\delta_{qr}(v_{pt}^{uv} \Gamma_{uv}^{st} +v_{uv}^{st} \Gamma_{pt}^{uv}) 
+  !+ \delta_{ps}(v_{uv}^{qt} \Gamma_{rt}^{uv} + v_{rt}^{uv}\Gamma_{uv}^{qt})] \\
+  !&+ \sum_{uv} (v_{pr}^{uv} \Gamma_{uv}^{qs} + v_{uv}^{qs}  \Gamma_{ps}^{uv}) \\
+  !&- \sum_{tu} (v_{pu}^{st} \Gamma_{rt}^{qu}+v_{pu}^{tr} \Gamma_{tr}^{qu}+v_{rt}^{qu}\Gamma_{pu}^{st} + v_{tr}^{qu}\Gamma_{pu}^{ts}) 
+  !\end{align*} 
+
+  !================
+  ! Initialization
+  !================
+  hessian = 0d0
+
+  CALL wall_time(t1)
+
+  !========================
+  ! First line, first term
+  !========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+           
+            if (q==r) then
+              do u = 1, mo_num
+
+                hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                    mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
+                  + mo_one_e_integrals(s,u) * one_e_dm_mo(p,u))
+
+              enddo
+            endif
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! First line, second term
+  !=========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+            .or. ((p==s) .and. (q==r))) then
+
+             if (p==s) then
+               do u = 1, mo_num
+
+                    hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                      mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
+                    + mo_one_e_integrals(q,u) * one_e_dm_mo(r,u))
+               enddo
+             endif
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! First line, third term
+  !========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+         
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+
+            hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+            - mo_one_e_integrals(s,p) * one_e_dm_mo(r,q) &
+            - mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
+
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! Second line, first term
+  !=========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+
+              if (q==r) then
+                do t = 1, mo_num
+                  do u = 1, mo_num
+                    do v = 1, mo_num
+
+                      hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * (  &
+                        get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
+                      + get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v))
+
+                    enddo
+                  enddo
+                enddo
+              endif
+            endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !==========================
+  ! Second line, second term
+  !==========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+              .or. ((p==s) .and. (q==r))) then
+
+             if (p==s) then
+               do t = 1, mo_num
+                 do u = 1, mo_num
+                   do v = 1, mo_num
+
+                     hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                       get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
+                     + get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t))
+
+                   enddo
+                 enddo
+               enddo
+             endif
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! Third line, first term
+  !========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+                .or. ((p==s) .and. (q==r))) then
+
+            do u = 1, mo_num
+              do v = 1, mo_num
+
+                hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+                 + get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
+                 + get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v)
+
+              enddo
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! Third line, second term
+  !=========================
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+           .or. ((p==s) .and. (q==r))) then
+
+            do t = 1, mo_num
+              do u = 1, mo_num
+
+                hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+                 - get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
+                 - get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
+                 - get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
+                 - get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
+
+              enddo
+            enddo
+
+          endif     
+  
+        enddo
+      enddo
+    enddo
+  enddo
+
+    CALL wall_time(t2)
+    t2 = t2 - t1
+    print*, 'Time to compute the hessian :', t2
+
+  !==============
+  ! Permutations
+  !==============
+ 
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  do tmp_r = 1, m
+    do tmp_s = 1, m
+      do tmp_q = 1, m
+        do tmp_p = 1, m
+
+          h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) - hessian(tmp_q,tmp_p,tmp_r,tmp_s) &
+                                          - hessian(tmp_p,tmp_q,tmp_s,tmp_r) + hessian(tmp_q,tmp_p,tmp_s,tmp_r)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! 4D matrix -> 2D matrix
+  !========================
+  
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  ! 4D mo_num matrix to 2D n matrix
+  do tmp_rs = 1, tmp_n
+    call vec_to_mat_index(tmp_rs,tmp_r,tmp_s)
+    do tmp_pq = 1, tmp_n
+      call vec_to_mat_index(tmp_pq,tmp_p,tmp_q)
+      tmp(tmp_pq,tmp_rs) = h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s)   
+    enddo
+  enddo
+
+  do p = 1, tmp_n
+    H(p) = tmp(p,p)
+  enddo
+
+  ! Display
+  if (debug) then 
+    print*,'2D diag Hessian matrix'
+    do tmp_pq = 1, tmp_n
+      write(*,'(100(F10.5))') tmp(tmp_pq,:)
+    enddo 
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(hessian,h_tmpr,tmp)
+
+  print*,'---End first_diag_hess_list---'
+
+end subroutine
+#+END_SRC

src/mo_optimization/org/first_diagonal_hessian_opt.org

@@ -0,0 +1,348 @@
+* First diagonal hessian
+
+#+BEGIN_SRC f90 :comments :tangle first_diagonal_hessian_opt.irp.f
+subroutine first_diag_hessian_opt(n,H, h_tmpr)
+
+  include 'constants.h' 
+
+  implicit none
+
+  !===========================================================================
+  ! Compute the diagonal hessian of energy with respects to orbital rotations
+  !===========================================================================
+
+  !===========
+  ! Variables 
+  !===========
+ 
+  ! in
+  integer, intent(in)           :: n 
+  ! n        : integer, n = mo_num*(mo_num-1)/2
+ 
+  ! out
+  double precision, intent(out) :: H(n,n), h_tmpr(mo_num,mo_num,mo_num,mo_num)
+  ! H        : n by n double precision matrix containing the 2D hessian
+  
+  ! internal
+  double precision, allocatable :: hessian(:,:,:,:)
+  integer                       :: p,q
+  integer                       :: r,s,t,u,v
+  integer                       :: pq,rs
+  double precision              :: t1,t2,t3
+  ! hessian  : mo_num 4D double precision matrix containing the hessian before the permutations
+  ! h_tmpr   : mo_num 4D double precision matrix containing the hessian after the permutations
+  ! p,q,r,s  : integer, indexes of the 4D hessian matrix
+  ! t,u,v    : integer, indexes to compute hessian elements
+  ! pq,rs    : integer, indexes for the conversion from 4D to 2D hessian matrix
+  ! t1,t2,t3 : double precision, t3 = t2 - t1, time to compute the hessian 
+  
+  ! Function
+  double precision :: get_two_e_integral
+  ! get_two_e_integral : double precision function, two e integrals
+ 
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo_alpha, one_e_dm_mo_beta : one body density matrix
+  ! two_e_dm_mo : two body density matrix
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(hessian(mo_num,mo_num,mo_num,mo_num))!,h_tmpr(mo_num,mo_num,mo_num,mo_num))
+
+  !=============
+  ! Calculation
+  !=============
+
+  if (debug) then
+          print*,'Enter in first_diag_hessien'
+  endif
+
+  ! From Anderson et. al. (2014) 
+  ! The Journal of Chemical Physics 141, 244104 (2014); doi: 10.1063/1.4904384
+
+  ! LaTeX formula :
+
+  !\begin{align*}
+  !H_{pq,rs} &= \dfrac{\partial^2 E(x)}{\partial x_{pq}^2} \\
+  !&= \mathcal{P}_{pq} \mathcal{P}_{rs} [ \frac{1}{2} \sum_u [\delta_{qr}(h_p^u \gamma_u^s + h_u^s \gamma_p^u) 
+  !+ \delta_{ps}(h_r^u \gamma_u^q + h_u^q \gamma_u^r)]
+  !-(h_p^s \gamma_r^q + h_r^q \gamma_p^s) \\
+  !&+ \frac{1}{2} \sum_{tuv} [\delta_{qr}(v_{pt}^{uv} \Gamma_{uv}^{st} +v_{uv}^{st} \Gamma_{pt}^{uv}) 
+  !+ \delta_{ps}(v_{uv}^{qt} \Gamma_{rt}^{uv} + v_{rt}^{uv}\Gamma_{uv}^{qt})] \\
+  !&+ \sum_{uv} (v_{pr}^{uv} \Gamma_{uv}^{qs} + v_{uv}^{qs}  \Gamma_{ps}^{uv}) \\
+  !&- \sum_{tu} (v_{pu}^{st} \Gamma_{rt}^{qu}+v_{pu}^{tr} \Gamma_{tr}^{qu}+v_{rt}^{qu}\Gamma_{pu}^{st} + v_{tr}^{qu}\Gamma_{pu}^{ts}) 
+  !\end{align*} 
+
+  !================
+  ! Initialization
+  !================
+  hessian = 0d0
+
+  CALL wall_time(t1)
+
+  !========================
+  ! First line, first term
+  !========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+           
+            if (q==r) then
+              do u = 1, mo_num
+
+                hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                    mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
+                  + mo_one_e_integrals(s,u) * one_e_dm_mo(p,u))
+
+              enddo
+            endif
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! First line, second term
+  !=========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+            .or. ((p==s) .and. (q==r))) then
+
+             if (p==s) then
+               do u = 1, mo_num
+
+                    hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                      mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
+                    + mo_one_e_integrals(q,u) * one_e_dm_mo(r,u))
+               enddo
+             endif
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! First line, third term
+  !========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+         
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+
+            hessian(p,q,r,s) = hessian(p,q,r,s) &
+            - mo_one_e_integrals(s,p) * one_e_dm_mo(r,q) &
+            - mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
+
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! Second line, first term
+  !=========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+             .or. ((p==s) .and. (q==r))) then
+
+              if (q==r) then
+                do t = 1, mo_num
+                  do u = 1, mo_num
+                    do v = 1, mo_num
+
+                      hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * (  &
+                        get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
+                      + get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v))
+
+                    enddo
+                  enddo
+                enddo
+              endif
+            endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !==========================
+  ! Second line, second term
+  !==========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+              .or. ((p==s) .and. (q==r))) then
+
+             if (p==s) then
+               do t = 1, mo_num
+                 do u = 1, mo_num
+                   do v = 1, mo_num
+
+                     hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                       get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
+                     + get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t))
+
+                   enddo
+                 enddo
+               enddo
+             endif
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! Third line, first term
+  !========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+
+           ! Permutations 
+           if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+                .or. ((p==s) .and. (q==r))) then
+
+            do u = 1, mo_num
+              do v = 1, mo_num
+
+                hessian(p,q,r,s) = hessian(p,q,r,s) &
+                 + get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
+                 + get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v)
+
+              enddo
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !=========================
+  ! Third line, second term
+  !=========================
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          ! Permutations 
+          if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
+           .or. ((p==s) .and. (q==r))) then
+
+            do t = 1, mo_num
+              do u = 1, mo_num
+
+                hessian(p,q,r,s) = hessian(p,q,r,s) &
+                 - get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
+                 - get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
+                 - get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
+                 - get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
+
+              enddo
+            enddo
+
+          endif     
+  
+        enddo
+      enddo
+    enddo
+  enddo
+
+    CALL wall_time(t2)
+    t2 = t2 - t1
+    print*, 'Time to compute the hessian :', t2
+
+  !==============
+  ! Permutations
+  !==============
+ 
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  do r = 1, mo_num
+    do s = 1, mo_num
+      do q = 1, mo_num
+        do p = 1, mo_num
+
+          h_tmpr(p,q,r,s) = (hessian(p,q,r,s) - hessian(q,p,r,s) - hessian(p,q,s,r) + hessian(q,p,s,r))
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! 4D matrix -> 2D matrix
+  !========================
+  
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  ! 4D mo_num matrix to 2D n matrix
+  do rs = 1, n
+    call vec_to_mat_index(rs,r,s)
+    do pq = 1, n
+      call vec_to_mat_index(pq,p,q)
+      H(pq,rs) = h_tmpr(p,q,r,s)   
+    enddo
+  enddo
+
+  ! Display
+  if (debug) then 
+    print*,'2D diag Hessian matrix'
+    do pq = 1, n
+      write(*,'(100(F10.5))') H(pq,:)
+    enddo 
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(hessian)
+
+  if (debug) then
+    print*,'Leave first_diag_hessien'
+  endif
+
+end subroutine
+#+END_SRC

src/mo_optimization/org/first_gradient_list_opt.org

@@ -0,0 +1,127 @@
+* First gradient
+#+BEGIN_SRC f90 :comments org :tangle first_gradient_list_opt.irp.f
+subroutine first_gradient_list_opt(tmp_n,m,list,v_grad)
+
+  include 'constants.h'
+
+  implicit none
+
+  !===================================================================
+  ! Compute the gradient of energy with respects to orbital rotations
+  !===================================================================
+
+  ! Check if read_wf = true, else :
+  ! qp set determinant read_wf true
+
+  ! in
+  integer, intent(in) :: tmp_n,m,list(m)
+  ! n : integer, n = m*(m-1)/2
+  ! m = list_size
+  
+  ! out
+  double precision, intent(out) :: v_grad(tmp_n)
+  ! v_grad : double precision vector of length n containeing the gradient
+
+  ! internal
+  double precision, allocatable :: grad(:,:),A(:,:)
+  double precision :: norm
+  integer :: i,p,q,r,s,t,tmp_i,tmp_p,tmp_q,tmp_r,tmp_s,tmp_t
+  ! grad : double precision matrix containing the gradient before the permutation
+  ! A : double precision matrix containing the gradient after the permutation
+  ! norm : double precision number, the norm of the vector gradient
+  ! i,p,q,r,s,t : integer, indexes 
+  ! istate : integer, the electronic state
+
+  ! Function
+  double precision :: get_two_e_integral, norm2
+  ! get_two_e_integral :  double precision function that gives the two e integrals
+  ! norm2 : double precision function that gives the norm of a vector
+ 
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo : one body density matrix (state average)
+  ! two_e_dm_mo : two body density matrix (state average)
+
+   print*,'---first_gradient_list---'
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(grad(m,m),A(m,m))
+
+  !=============
+  ! Calculation
+  !=============
+
+  v_grad = 0d0
+  grad = 0d0
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      !grad(tmp_p,tmp_q) = 0d0
+      do r = 1, mo_num
+        grad(tmp_p,tmp_q) = grad(tmp_p,tmp_q) + mo_one_e_integrals(p,r) * one_e_dm_mo(r,q) &
+                               - mo_one_e_integrals(r,q) * one_e_dm_mo(p,r)
+
+      enddo
+
+      do r = 1, mo_num
+        do s = 1, mo_num
+          do t = 1, mo_num
+
+            grad(tmp_p,tmp_q) = grad(tmp_p,tmp_q) &
+                + get_two_e_integral(p,t,r,s,mo_integrals_map) * two_e_dm_mo(r,s,q,t) &
+                - get_two_e_integral(r,s,q,t,mo_integrals_map) * two_e_dm_mo(p,t,r,s)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ! Conversion mo_num*mo_num matrix to mo_num(mo_num-1)/2 vector
+  do tmp_i = 1, tmp_n
+    call vec_to_mat_index(tmp_i,tmp_p,tmp_q)
+    v_grad(tmp_i)=(grad(tmp_p,tmp_q) - grad(tmp_q,tmp_p))
+  enddo  
+
+  ! Display, vector containing the gradient elements 
+  if (debug) then  
+    print*,'Vector containing the gradient :'
+    write(*,'(100(F10.5))') v_grad(1:tmp_n)
+  endif  
+
+  ! Norm of the vector
+  norm = norm2(v_grad)
+  print*, 'Norm : ', norm
+
+  ! Matrix gradient
+  A = 0d0
+  do tmp_q = 1, m
+    do tmp_p = 1, m
+      A(tmp_p,tmp_q) = grad(tmp_p,tmp_q) - grad(tmp_q,tmp_p)
+    enddo
+  enddo
+
+  ! Display, matrix containting the gradient elements
+  if (debug) then
+    print*,'Matrix containing the gradient :'
+    do tmp_i = 1, m
+      write(*,'(100(E12.5))') A(tmp_i,1:m)
+    enddo
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(grad,A)
+
+  print*,'---End first_gradient_list---'
+
+end subroutine
+
+#+END_SRC

src/mo_optimization/org/first_gradient_opt.org

@@ -0,0 +1,130 @@
+* First gradient
+#+BEGIN_SRC f90 :comments org :tangle first_gradient_opt.irp.f
+subroutine first_gradient_opt(n,v_grad)
+
+  include 'constants.h'
+
+  implicit none
+
+  !===================================================================
+  ! Compute the gradient of energy with respects to orbital rotations
+  !===================================================================
+
+  ! Check if read_wf = true, else :
+  ! qp set determinant read_wf true
+
+  END_DOC
+
+  ! in
+  integer, intent(in) :: n
+  ! n : integer, n = mo_num*(mo_num-1)/2
+  
+  ! out
+  double precision, intent(out) :: v_grad(n)
+  ! v_grad : double precision vector of length n containeing the gradient
+
+  ! internal
+  double precision, allocatable :: grad(:,:),A(:,:)
+  double precision :: norm
+  integer :: i,p,q,r,s,t
+  integer :: istate
+  ! grad : double precision matrix containing the gradient before the permutation
+  ! A : double precision matrix containing the gradient after the permutation
+  ! norm : double precision number, the norm of the vector gradient
+  ! i,p,q,r,s,t : integer, indexes 
+  ! istate : integer, the electronic state
+
+  ! Function
+  double precision :: get_two_e_integral, norm2
+  ! get_two_e_integral :  double precision function that gives the two e integrals
+  ! norm2 : double precision function that gives the norm of a vector
+ 
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo : one body density matrix (state average)
+  ! two_e_dm_mo : two body density matrix (state average)
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(grad(mo_num,mo_num),A(mo_num,mo_num))
+
+  !=============
+  ! Calculation
+  !=============
+
+  if (debug) then
+    print*,'---first_gradient---'
+  endif
+
+  v_grad = 0d0
+
+  do p = 1, mo_num
+    do q = 1, mo_num
+      grad(p,q) = 0d0
+      do r = 1, mo_num
+        grad(p,q) = grad(p,q) + mo_one_e_integrals(p,r) * one_e_dm_mo(r,q) &
+                               - mo_one_e_integrals(r,q) * one_e_dm_mo(p,r)
+
+      enddo
+
+      do r = 1, mo_num
+        do s = 1, mo_num
+          do t= 1, mo_num
+
+            grad(p,q) = grad(p,q) &
+                + get_two_e_integral(p,t,r,s,mo_integrals_map) * two_e_dm_mo(r,s,q,t) &
+                - get_two_e_integral(r,s,q,t,mo_integrals_map) * two_e_dm_mo(p,t,r,s)
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ! Conversion mo_num*mo_num matrix to mo_num(mo_num-1)/2 vector
+  do i=1,n
+    call vec_to_mat_index(i,p,q)
+    v_grad(i)=(grad(p,q) - grad(q,p))
+  enddo  
+
+  ! Display, vector containing the gradient elements 
+  if (debug) then  
+    print*,'Vector containing the gradient :'
+    write(*,'(100(F10.5))') v_grad(1:n)
+  endif  
+
+  ! Norm of the vector
+  norm = norm2(v_grad)
+  print*, 'Norm : ', norm
+
+  ! Matrix gradient
+  A = 0d0
+  do q=1,mo_num
+    do p=1,mo_num
+      A(p,q) = grad(p,q) - grad(q,p)
+    enddo
+  enddo
+
+  ! Display, matrix containting the gradient elements
+  if (debug) then    
+    print*,'Matrix containing the gradient :'
+    do i = 1, mo_num
+      write(*,'(100(E12.5))') A(i,1:mo_num)
+    enddo
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(grad,A)
+
+  if (debug) then
+    print*,'---End first_gradient---'
+  endif
+
+end subroutine
+
+#+END_SRC

src/mo_optimization/org/first_hessian_list_opt.org

@@ -0,0 +1,370 @@
+* First hessian
+
+#+BEGIN_SRC f90 :comments :tangle first_hessian_list_opt.irp.f
+subroutine first_hessian_list_opt(tmp_n,m,list,H,h_tmpr)
+
+  include 'constants.h' 
+
+  implicit none
+
+  !==================================================================
+  ! Compute the hessian of energy with respects to orbital rotations
+  !==================================================================
+
+  !===========
+  ! Variables  
+  !===========
+
+  ! in
+  integer, intent(in)           :: tmp_n, m, list(m)
+  !tmp_n         : integer, tmp_n = m*(m-1)/2
+  
+  ! out
+  double precision, intent(out) :: H(tmp_n,tmp_n),h_tmpr(m,m,m,m)
+  ! H        : n by n double precision matrix containing the 2D hessian
+ 
+  ! internal
+  double precision, allocatable :: hessian(:,:,:,:)
+  integer                       :: p,q, tmp_p,tmp_q
+  integer                       :: r,s,t,u,v,tmp_r,tmp_s,tmp_t,tmp_u,tmp_v
+  integer                       :: pq,rs,tmp_pq,tmp_rs
+  double precision              :: t1,t2,t3,t4,t5,t6
+  ! hessian  : mo_num 4D double precision matrix containing the hessian before the permutations
+  ! h_tmpr   : mo_num 4D double precision matrix containing the hessian after the permutations
+  ! p,q,r,s  : integer, indexes of the 4D hessian matrix
+  ! t,u,v    : integer, indexes to compute hessian elements
+  ! pq,rs    : integer, indexes for the conversion from 4D to 2D hessian matrix
+  ! t1,t2,t3 : double precision, t3 = t2 - t1, time to compute the hessian 
+
+  ! Funtion 
+  double precision              :: get_two_e_integral
+  ! get_two_e_integral :  double precision function, two e integrals 
+
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo_alpha, one_e_dm_mo_beta : one body density matrix
+  ! two_e_dm_mo : two body density matrix
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(hessian(m,m,m,m))
+
+  !=============
+  ! Calculation
+  !=============
+
+  print*,'---first_hess_list---'
+
+  ! From Anderson et. al. (2014) 
+  ! The Journal of Chemical Physics 141, 244104 (2014); doi: 10.1063/1.4904384
+
+  CALL wall_time(t1)
+
+  ! Initialization
+  hessian = 0d0
+
+  !========================
+  ! First line, first term
+  !========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+            if (q==r) then
+              do u = 1, mo_num
+
+                hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * (  &
+                  mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
+                + mo_one_e_integrals(s,u) * one_e_dm_mo(p,u))
+
+              enddo
+            endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 1 :', t6
+
+  !=========================
+  ! First line, second term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          if (p==s) then
+            do u = 1, mo_num
+
+                  hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                    mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
+                  + mo_one_e_integrals(q,u) * one_e_dm_mo(r,u))
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+  
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 2 :', t6
+
+  !========================
+  ! First line, third term
+  !========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+          - mo_one_e_integrals(s,p) * one_e_dm_mo(r,q)&
+          - mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 3 :', t6
+
+
+  !=========================
+  ! Second line, first term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+           if (q==r) then
+             do t = 1, mo_num
+               do u = 1, mo_num
+                 do v = 1, mo_num
+
+                   hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * (  &
+                     get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
+                   + get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v))
+
+                 enddo
+               enddo
+             enddo
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l2 1 :', t6
+
+  !==========================
+  ! Second line, second term
+  !==========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          if (p==s) then
+            do t = 1, mo_num
+              do u = 1, mo_num
+                do v = 1, mo_num
+
+                  hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) + 0.5d0 * ( &
+                    get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
+                  + get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t))
+
+                enddo
+              enddo
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l2 2 :', t6
+
+  !========================
+  ! Third line, first term
+  !========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          do u = 1, mo_num
+            do v = 1, mo_num
+
+              hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+               + get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
+               + get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v)
+
+            enddo
+          enddo
+
+        enddo
+      enddo
+    enddo
+  enddo
+ 
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l3 1 :', t6
+
+  !=========================
+  ! Third line, second term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do tmp_p = 1, m
+    p = list(tmp_p)
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do tmp_r = 1, m
+        r = list(tmp_r)
+        do tmp_s = 1, m
+          s = list(tmp_s)
+
+          do t = 1, mo_num
+            do u = 1, mo_num
+
+              hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
+               - get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
+               - get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
+               - get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
+               - get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
+
+            enddo
+          enddo
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l3 2 :', t6
+
+  CALL wall_time(t2)
+  t3 = t2 -t1
+  print*,'Time to compute the hessian : ', t3
+
+  !==============
+  ! Permutations 
+  !==============
+
+  ! Hessian(p,q,r,s) = P_pq P_rs [ ...]
+  ! => Hessian(p,q,r,s) = (p,q,r,s) - (q,p,r,s) - (p,q,s,r) + (q,p,s,r) 
+
+  do tmp_s = 1, m
+    do tmp_r = 1, m
+      do tmp_q = 1, m
+        do tmp_p = 1, m
+
+          h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s) = (hessian(tmp_p,tmp_q,tmp_r,tmp_s) - hessian(tmp_q,tmp_p,tmp_r,tmp_s) &
+                                           - hessian(tmp_p,tmp_q,tmp_s,tmp_r) + hessian(tmp_q,tmp_p,tmp_s,tmp_r))
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! 4D matrix to 2D matrix
+  !========================
+
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  ! 4D mo_num matrix to 2D n matrix
+  do tmp_pq = 1, tmp_n
+    call vec_to_mat_index(tmp_pq,tmp_p,tmp_q)
+    do tmp_rs = 1, tmp_n
+      call vec_to_mat_index(tmp_rs,tmp_r,tmp_s)
+      H(tmp_pq,tmp_rs) = h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s)   
+    enddo
+  enddo
+
+  ! Display
+  if (debug) then 
+    print*,'2D Hessian matrix'
+    do tmp_pq = 1, tmp_n
+      write(*,'(100(F10.5))') H(tmp_pq,:)
+    enddo 
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(hessian)
+
+  print*,'---End first_hess_list---'
+
+end subroutine
+
+#+END_SRC

src/mo_optimization/org/first_hessian_opt.org

@@ -0,0 +1,365 @@
+* First hessian
+
+#+BEGIN_SRC f90 :comments :tangle first_hessian_opt.irp.f
+subroutine first_hessian_opt(n,H,h_tmpr)
+
+  include 'constants.h' 
+
+  implicit none
+
+  !==================================================================
+  ! Compute the hessian of energy with respects to orbital rotations
+  !==================================================================
+
+  !===========
+  ! Variables  
+  !===========
+
+  ! in
+  integer, intent(in)           :: n 
+  !n         : integer, n = mo_num*(mo_num-1)/2
+  
+  ! out
+  double precision, intent(out) :: H(n,n),h_tmpr(mo_num,mo_num,mo_num,mo_num)
+  ! H        : n by n double precision matrix containing the 2D hessian
+ 
+  ! internal
+  double precision, allocatable :: hessian(:,:,:,:)
+  integer                       :: p,q
+  integer                       :: r,s,t,u,v
+  integer                       :: pq,rs
+  double precision              :: t1,t2,t3,t4,t5,t6
+  ! hessian  : mo_num 4D double precision matrix containing the hessian before the permutations
+  ! h_tmpr   : mo_num 4D double precision matrix containing the hessian after the permutations
+  ! p,q,r,s  : integer, indexes of the 4D hessian matrix
+  ! t,u,v    : integer, indexes to compute hessian elements
+  ! pq,rs    : integer, indexes for the conversion from 4D to 2D hessian matrix
+  ! t1,t2,t3 : double precision, t3 = t2 - t1, time to compute the hessian 
+
+  ! Funtion 
+  double precision              :: get_two_e_integral
+  ! get_two_e_integral :  double precision function, two e integrals 
+
+  ! Provided :
+  ! mo_one_e_integrals : mono e- integrals
+  ! get_two_e_integral : two e- integrals
+  ! one_e_dm_mo_alpha, one_e_dm_mo_beta : one body density matrix
+  ! two_e_dm_mo : two body density matrix
+
+  !============
+  ! Allocation
+  !============
+
+  allocate(hessian(mo_num,mo_num,mo_num,mo_num))
+
+  !=============
+  ! Calculation
+  !=============
+
+  if (debug) then
+    print*,'Enter in first_hess'
+  endif
+
+  ! From Anderson et. al. (2014) 
+  ! The Journal of Chemical Physics 141, 244104 (2014); doi: 10.1063/1.4904384
+
+  CALL wall_time(t1)
+
+  ! Initialization
+  hessian = 0d0
+
+  !========================
+  ! First line, first term
+  !========================
+
+  CALL wall_time(t4)
+
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+            if (q==r) then
+              do u = 1, mo_num
+
+                hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * (  &
+                  mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
+                + mo_one_e_integrals(s,u) * one_e_dm_mo(p,u))
+
+              enddo
+            endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 1 :', t6
+
+  !=========================
+  ! First line, second term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          if (p==s) then
+            do u = 1, mo_num
+
+                  hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                    mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
+                  + mo_one_e_integrals(q,u) * one_e_dm_mo(r,u))
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+  
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 2 :', t6
+
+  !========================
+  ! First line, third term
+  !========================
+
+  CALL wall_time(t4)
+
+  do p = 1, mo_num
+    do q = 1, mo_num
+      do r = 1, mo_num
+        do s = 1, mo_num
+
+          hessian(p,q,r,s) = hessian(p,q,r,s) &
+          - mo_one_e_integrals(s,p) * one_e_dm_mo(r,q)&
+          - mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l1 3 :', t6
+
+
+  !=========================
+  ! Second line, first term
+  !=========================
+
+  CALL wall_time(t4)
+
+  do s = 1, mo_num
+    do r = 1, mo_num
+      do q = 1, mo_num
+        do p = 1, mo_num
+
+!  do p = 1, mo_num
+!    do q = 1, mo_num
+!      do r = 1, mo_num
+!        do s = 1, mo_num
+
+           if (q==r) then
+             do t = 1, mo_num
+               do u = 1, mo_num
+                 do v = 1, mo_num
+
+                   hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * (  &
+                     get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
+                   + get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v))
+
+                 enddo
+               enddo
+             enddo
+           endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l2 1 :', t6
+
+  !==========================
+  ! Second line, second term
+  !==========================
+
+  CALL wall_time(t4)
+
+  do s = 1, mo_num
+    do r = 1, mo_num
+      do q = 1, mo_num
+        do p = 1, mo_num
+
+!  do p = 1, mo_num
+!    do q = 1, mo_num
+!      do r = 1, mo_num
+!        do s = 1, mo_num
+
+          if (p==s) then
+            do t = 1, mo_num
+              do u = 1, mo_num
+                do v = 1, mo_num
+
+                  hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( &
+                    get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
+                  + get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t))
+
+                enddo
+              enddo
+            enddo
+          endif
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l2 2 :', t6
+
+  !========================
+  ! Third line, first term
+  !========================
+
+  CALL wall_time(t4)
+
+      do s = 1, mo_num
+        do r = 1, mo_num
+          do q = 1, mo_num
+            do p = 1, mo_num
+
+!  do p = 1, mo_num
+!    do q = 1, mo_num
+!      do r = 1, mo_num
+!        do s = 1, mo_num
+
+          do u = 1, mo_num
+            do v = 1, mo_num
+
+              hessian(p,q,r,s) = hessian(p,q,r,s) &
+               + get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
+               + get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v)
+
+            enddo
+          enddo
+
+        enddo
+      enddo
+    enddo
+  enddo
+ 
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l3 1 :', t6
+
+  !=========================
+  ! Third line, second term
+  !=========================
+
+  CALL wall_time(t4)
+
+      do s = 1, mo_num
+        do r = 1, mo_num
+          do q = 1, mo_num
+            do p = 1, mo_num
+
+!  do p = 1, mo_num
+!    do q = 1, mo_num
+!      do r = 1, mo_num
+!        do s = 1, mo_num
+
+          do t = 1, mo_num
+            do u = 1, mo_num
+
+              hessian(p,q,r,s) = hessian(p,q,r,s) &
+               - get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
+               - get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
+               - get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
+               - get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
+
+            enddo
+          enddo
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  CALL wall_time(t5)
+  t6 = t5-t4
+  print*,'l3 2 :', t6
+
+  CALL wall_time(t2)
+  t3 = t2 -t1
+  print*,'Time to compute the hessian : ', t3
+
+  !==============
+  ! Permutations 
+  !==============
+
+  ! Hessian(p,q,r,s) = P_pq P_rs [ ...]
+  ! => Hessian(p,q,r,s) = (p,q,r,s) - (q,p,r,s) - (p,q,s,r) + (q,p,s,r) 
+
+  do s = 1, mo_num
+    do r = 1, mo_num
+      do q = 1, mo_num
+        do p = 1, mo_num
+
+          h_tmpr(p,q,r,s) = (hessian(p,q,r,s) - hessian(q,p,r,s) - hessian(p,q,s,r) + hessian(q,p,s,r))
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !========================
+  ! 4D matrix to 2D matrix
+  !========================
+
+  ! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
+  ! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
+  ! H(pq,rs) : p<q and r<s
+
+  ! 4D mo_num matrix to 2D n matrix
+  do pq = 1, n
+    call vec_to_mat_index(pq,p,q)
+    do rs = 1, n
+      call vec_to_mat_index(rs,r,s)
+      H(pq,rs) = h_tmpr(p,q,r,s)   
+    enddo
+  enddo
+
+  ! Display
+  if (debug) then 
+    print*,'2D Hessian matrix'
+    do pq = 1, n
+      write(*,'(100(F10.5))') H(pq,:)
+    enddo 
+  endif
+
+  !==============
+  ! Deallocation
+  !==============
+
+  deallocate(hessian)
+
+  if (debug) then
+    print*,'Leave first_hess'
+  endif
+
+end subroutine
+
+#+END_SRC

src/mo_optimization/org/gradient_list_opt.org

@@ -0,0 +1,393 @@
+* Gradient
+
+The gradient of the CI energy with respects to the orbital rotation
+is:
+(C-c C-x C-l)
+$$
+G(p,q) = \mathcal{P}_{pq} \left[ \sum_r (h_p^r \gamma_r^q - h_r^q \gamma_p^r) +
+\sum_{rst}(v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs})
+\right]
+$$
+
+
+$$
+\mathcal{P}_{pq}= 1 - (p \leftrightarrow q)
+$$
+
+$$
+G(p,q) = \left[
+\sum_r (h_p^r \gamma_r^q - h_r^q \gamma_p^r) +
+\sum_{rst}(v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs})
+\right] - 
+\left[
+\sum_r (h_q^r \gamma_r^p - h_r^p \gamma_q^r) +
+\sum_{rst}(v_{qt}^{rs} \Gamma_{rs}^{pt} - v_{rs}^{pt}
+\Gamma_{qt}^{rs})
+\right]
+$$
+
+Where p,q,r,s,t are general spatial orbitals
+mo_num : the number of molecular orbitals
+$$h$$ : One electron integrals
+$$\gamma$$ : One body density matrix (state average in our case)
+$$v$$ : Two electron integrals
+$$\Gamma$$ : Two body density matrice (state average in our case)
+
+The gradient is a mo_num by mo_num matrix, p,q,r,s,t take all the
+values between 1 and mo_num (1 and mo_num include).
+
+To do that we compute $$G(p,q)$$ for all the pairs (p,q).
+
+Source :
+Seniority-based coupled cluster theory
+J. Chem. Phys. 141, 244104 (2014); https://doi.org/10.1063/1.4904384
+Thomas M. Henderson, Ireneusz W. Bulik, Tamar Stein, and Gustavo
+E. Scuseria
+
+*Compute the gradient of energy with respects to orbital rotations*
+
+Provided:
+| mo_num                                   | integer          | number of MOs             |
+| mo_one_e_integrals(mo_num,mo_num)        | double precision | mono_electronic integrals |
+| one_e_dm_mo(mo_num,mo_num)               | double precision | one e- density matrix     |
+| two_e_dm_mo(mo_num,mo_num,mo_num,mo_num) | double precision | two e- density matrix     |
+
+Input:
+| n | integer | mo_num*(mo_num-1)/2 |
+
+Output:
+| v_grad(n) | double precision | the gradient                    |
+| max_elem  | double precision | maximum element of the gradient |
+
+Internal:
+| grad(mo_num,mo_num)                | double precison   | gradient before the tranformation in a vector             |
+| A((mo_num,mo_num)                  | doubre precision  | gradient after the permutations                           |
+| norm                               | double precision  | norm of the gradient                                      |
+| p, q                               | integer           | indexes of the element in the matrix grad                 |
+| i                                  | integer           | index for the tranformation in a vector                   |
+| r, s, t                            | integer           | indexes dor the sums                                      |
+| t1, t2, t3                         | double precision  | t3 = t2 - t1, time to compute the gradient                |
+| t4, t5, t6                         | double precission | t6 = t5 - t4, time to compute each element                |
+| tmp_bi_int_3(mo_num,mo_num,mo_num) | double precision  | 3 indexes temporary array for the bi-electronic integrals |
+| tmp_2rdm_3(mo_num,mo_num,mo_num)   | double precision  | 3 indexes temporary array for the two e- density matrix   |
+| tmp_accu(mo_num,mo_num)            | double precision  | temporary array                                           |
+
+Function:
+| get_two_e_integral | double precision | bi-electronic integrals |
+| dnrm2              | double precision | (Lapack) norm           |
+
+#+BEGIN_SRC f90 :comments org :tangle gradient_list_opt.irp.f
+subroutine gradient_list_opt(n,m,list,v_grad,max_elem,norm)
+  use omp_lib
+  include 'constants.h'
+
+  implicit none
+  
+  ! Variables
+  
+  ! in
+  integer, intent(in) :: n,m,list(m)
+  
+  ! out
+  double precision, intent(out) :: v_grad(n), max_elem, norm
+
+  ! internal
+  double precision, allocatable :: grad(:,:),A(:,:)
+  integer                       :: i,p,q,r,s,t, tmp_p, tmp_q, tmp_i
+  double precision              :: t1,t2,t3,t4,t5,t6
+
+  double precision, allocatable :: tmp_accu(:,:), tmp_mo_one_e_integrals(:,:),tmp_one_e_dm_mo(:,:)
+  double precision, allocatable :: tmp_bi_int_3(:,:,:), tmp_2rdm_3(:,:,:)
+
+  ! Functions
+  double precision :: get_two_e_integral, dnrm2
+
+  
+  print*,''
+  print*,'---gradient---'
+
+  ! Allocation of shared arrays
+  allocate(grad(m,m),A(m,m)) 
+  allocate(tmp_mo_one_e_integrals(m,mo_num),tmp_one_e_dm_mo(mo_num,m))
+
+
+  ! Initialization omp
+  call omp_set_max_active_levels(1)
+
+  !$OMP PARALLEL                                                 &
+      !$OMP PRIVATE(                                             &
+      !$OMP   p,q,r,s,t,tmp_p,tmp_q,                             &
+      !$OMP   tmp_accu,tmp_bi_int_3, tmp_2rdm_3)                &
+      !$OMP SHARED(grad, one_e_dm_mo,m,list,mo_num,mo_one_e_integrals, &
+      !$OMP mo_integrals_map,tmp_one_e_dm_mo, tmp_mo_one_e_integrals,t4,t5,t6)                           &
+      !$OMP DEFAULT(SHARED)
+ 
+  ! Allocation of private arrays
+  allocate(tmp_accu(m,m))
+  allocate(tmp_bi_int_3(mo_num,mo_num,m))
+  allocate(tmp_2rdm_3(mo_num,mo_num,m))
+#+END_SRC
+  
+** Calculation
+*** Initialization
+#+BEGIN_SRC f90 :comments org :tangle gradient_list_opt.irp.f
+  !$OMP DO
+  do tmp_q = 1, m
+    do tmp_p = 1, m
+      grad(tmp_p,tmp_q) = 0d0
+    enddo
+  enddo
+  !$OMP END DO
+#+END_SRC
+
+*** Term 1
+  
+Without optimization the term 1 is :
+
+do p = 1, mo_num
+  do q = 1, mo_num
+     do r = 1, mo_num
+       grad(p,q) = grad(p,q) &
+               + mo_one_e_integrals(p,r) * one_e_dm_mo(r,q) &
+               - mo_one_e_integrals(r,q) * one_e_dm_mo(p,r)
+    enddo
+  enddo
+enddo
+   
+Since the matrix multiplication A.B is defined like :
+\begin{equation}
+c_{ij} = \sum_k a_{ik}.b_{kj}
+\end{equation}
+The previous equation can be rewritten as a matrix multplication  
+  
+#+BEGIN_SRC f90 :comments org :tangle gradient_list_opt.irp.f
+  !****************
+  ! Opt first term
+  !****************
+
+  !$OMP DO
+  do r = 1, mo_num
+    do tmp_p = 1, m
+      p = list(tmp_p)
+      tmp_mo_one_e_integrals(tmp_p,r) = mo_one_e_integrals(p,r)
+    enddo
+  enddo
+  !$OMP END DO
+
+  !$OMP DO
+  do tmp_q = 1, m
+    q = list(tmp_q) 
+    do r = 1, mo_num
+      tmp_one_e_dm_mo(r,tmp_q) = one_e_dm_mo(r,q)
+    enddo
+  enddo
+  !$OMP END DO
+
+  call dgemm('N','N',m,m,mo_num,1d0,&
+  tmp_mo_one_e_integrals, size(tmp_mo_one_e_integrals,1),&
+  tmp_one_e_dm_mo,size(tmp_one_e_dm_mo,1),0d0,tmp_accu,size(tmp_accu,1))
+
+  !$OMP DO
+  do tmp_q = 1, m
+    do tmp_p = 1, m
+
+      grad(tmp_p,tmp_q) = grad(tmp_p,tmp_q) + (tmp_accu(tmp_p,tmp_q) - tmp_accu(tmp_q,tmp_p))
+
+    enddo
+  enddo 
+  !$OMP END DO
+
+  !$OMP MASTER
+  CALL wall_TIME(t4)
+  !$OMP END MASTER
+
+!  call dgemm('N','N',mo_num,mo_num,mo_num,1d0,mo_one_e_integrals,&
+!  mo_num,one_e_dm_mo,mo_num,0d0,tmp_accu,mo_num)
+!  
+!  !$OMP DO
+!  do q = 1, mo_num
+!    do p = 1, mo_num
+!
+!      grad(p,q) = grad(p,q) + (tmp_accu(p,q) - tmp_accu(q,p))
+!
+!    enddo
+!  enddo 
+!  !$OMP END DO
+  
+  !$OMP MASTER
+  CALL wall_TIME(t5)
+  t6 = t5-t4
+  print*,'Gradient, first term (s) :', t6 
+  !$OMP END MASTER
+#+END_SRC
+
+*** Term 2
+ 
+Without optimization the second term is : 
+
+do p = 1, mo_num
+  do q = 1, mo_num 
+    do r = 1, mo_num
+      do s = 1, mo_num
+        do t= 1, mo_num
+
+        grad(p,q) = grad(p,q) &
+                + get_two_e_integral(p,t,r,s,mo_integrals_map) * two_e_dm_mo(r,s,q,t) &
+                - get_two_e_integral(r,s,q,t,mo_integrals_map) * two_e_dm_mo(p,t,r,s)
+       enddo
+      enddo
+    enddo
+  enddo
+enddo
+
+Using the bielectronic integral properties :
+get_two_e_integral(p,t,r,s,mo_integrals_map) = get_two_e_integral(r,s,p,t,mo_integrals_map)
+
+Using the two body matrix properties :
+two_e_dm_mo(p,t,r,s) = two_e_dm_mo(r,s,p,t)
+
+t is one the right, we can put it on the external loop and create 3
+indexes temporary array 
+r,s can be seen as one index
+
+By doing so, a matrix multiplication appears 
+
+#+BEGIN_SRC f90 :comments org :tangle gradient_list_opt.irp.f
+  !*****************
+  ! Opt second term  
+  !*****************
+
+  !$OMP MASTER
+  CALL wall_TIME(t4)
+  !$OMP END MASTER 
+ 
+  !$OMP DO
+  do t = 1, mo_num
+    
+    do tmp_p = 1, m
+      p = list(tmp_p)
+      do s = 1, mo_num
+        do r = 1, mo_num
+            
+          tmp_bi_int_3(r,s,tmp_p) = get_two_e_integral(r,s,p,t,mo_integrals_map)
+         
+        enddo
+      enddo
+    enddo
+
+    do tmp_q = 1, m
+      q = list(tmp_q)
+      do s = 1, mo_num
+        do r = 1, mo_num
+             
+           tmp_2rdm_3(r,s,tmp_q) = two_e_dm_mo(r,s,q,t)
+  
+        enddo
+      enddo
+    enddo
+
+    call dgemm('T','N',m,m,mo_num*mo_num,1d0,tmp_bi_int_3,&
+      mo_num*mo_num,tmp_2rdm_3,mo_num*mo_num,0d0,tmp_accu,size(tmp_accu,1))
+
+    !$OMP CRITICAL   
+    do tmp_q = 1, m
+      do tmp_p = 1, m
+
+        grad(tmp_p,tmp_q) = grad(tmp_p,tmp_q) + tmp_accu(tmp_p,tmp_q) - tmp_accu(tmp_q,tmp_p)
+
+      enddo
+    enddo
+    !$OMP END CRITICAL
+
+  enddo
+  !$OMP END DO
+
+  !$OMP MASTER
+  CALL wall_TIME(t5)
+  t6 = t5-t4
+  print*,'Gradient second term (s) : ', t6
+  !$OMP END MASTER  
+#+END_SRC
+
+*** Deallocation of private arrays
+#+BEGIN_SRC f90 :comments org :tangle gradient_list_opt.irp.f
+  deallocate(tmp_bi_int_3,tmp_2rdm_3,tmp_accu)
+
+  !$OMP END PARALLEL
+
+  call omp_set_max_active_levels(4)
+#+END_SRC
+
+*** Permutation, 2D matrix -> vector, transformation
+In addition there is a permutation in the gradient formula :
+\begin{equation}
+P_{pq} = 1 - (p <-> q) 
+\end{equation}
+
+We need a vector to use the gradient. Here the gradient is a
+antisymetric matrix so we can transform it in a vector of length
+mo_num*(mo_num-1)/2.
+
+Here we do these two things at the same time.
+
+#+BEGIN_SRC f90 :comments org :tangle gradient_list_opt.irp.f
+  do i=1,n
+    call vec_to_mat_index(i,p,q)
+    v_grad(i)=(grad(p,q) - grad(q,p))
+  enddo  
+
+  ! Debug, diplay the vector containing the gradient elements 
+  if (debug) then  
+    print*,'Vector containing the gradient :'
+    write(*,'(100(F10.5))') v_grad(1:n)
+  endif  
+#+END_SRC
+
+*** Norm of the gradient
+The norm can be useful.
+#+BEGIN_SRC f90 :comments org :tangle gradient_list_opt.irp.f
+  norm = dnrm2(n,v_grad,1)
+  print*, 'Gradient norm : ', norm
+#+END_SRC
+
+*** Maximum element in the gradient
+The maximum element in the gradient is very important for the
+convergence criterion of the Newton method.
+
+#+BEGIN_SRC f90 :comments org :tangle gradient_list_opt.irp.f
+  ! Max element of the gradient
+  max_elem = 0d0
+  do i = 1, n
+    if (DABS(v_grad(i)) > DABS(max_elem)) then
+      max_elem = v_grad(i)
+    endif
+  enddo
+
+  print*,'Max element in the gradient :', max_elem  
+
+  ! Debug, display the matrix containting the gradient elements
+  if (debug) then
+    ! Matrix gradient
+    A = 0d0
+    do q=1,m
+      do p=1,m
+        A(p,q) = grad(p,q) - grad(q,p)
+      enddo
+    enddo
+    print*,'Matrix containing the gradient :'
+    do i = 1, m
+      write(*,'(100(F10.5))') A(i,1:m)
+    enddo
+  endif
+#+END_SRC
+
+*** Deallocation of shared arrays and end
+#+BEGIN_SRC f90 :comments org :tangle gradient_list_opt.irp.f
+  deallocate(grad,A, tmp_mo_one_e_integrals,tmp_one_e_dm_mo)
+
+  print*,'---End gradient---'
+
+  end subroutine
+
+#+END_SRC
+  

src/mo_optimization/org/gradient_opt.org

@@ -0,0 +1,358 @@
+* Gradient
+
+The gradient of the CI energy with respects to the orbital rotation
+is:
+(C-c C-x C-l)
+$$
+G(p,q) = \mathcal{P}_{pq} \left[ \sum_r (h_p^r \gamma_r^q - h_r^q \gamma_p^r) +
+\sum_{rst}(v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs})
+\right]
+$$
+
+
+$$
+\mathcal{P}_{pq}= 1 - (p \leftrightarrow q)
+$$
+
+$$
+G(p,q) = \left[
+\sum_r (h_p^r \gamma_r^q - h_r^q \gamma_p^r) +
+\sum_{rst}(v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs})
+\right] - 
+\left[
+\sum_r (h_q^r \gamma_r^p - h_r^p \gamma_q^r) +
+\sum_{rst}(v_{qt}^{rs} \Gamma_{rs}^{pt} - v_{rs}^{pt}
+\Gamma_{qt}^{rs})
+\right]
+$$
+
+Where p,q,r,s,t are general spatial orbitals
+mo_num : the number of molecular orbitals
+$$h$$ : One electron integrals
+$$\gamma$$ : One body density matrix (state average in our case)
+$$v$$ : Two electron integrals
+$$\Gamma$$ : Two body density matrice (state average in our case)
+
+The gradient is a mo_num by mo_num matrix, p,q,r,s,t take all the
+values between 1 and mo_num (1 and mo_num include).
+
+To do that we compute $$G(p,q)$$ for all the pairs (p,q).
+
+Source :
+Seniority-based coupled cluster theory
+J. Chem. Phys. 141, 244104 (2014); https://doi.org/10.1063/1.4904384
+Thomas M. Henderson, Ireneusz W. Bulik, Tamar Stein, and Gustavo
+E. Scuseria
+
+*Compute the gradient of energy with respects to orbital rotations*
+
+Provided:
+| mo_num                                   | integer          | number of MOs             |
+| mo_one_e_integrals(mo_num,mo_num)        | double precision | mono_electronic integrals |
+| one_e_dm_mo(mo_num,mo_num)               | double precision | one e- density matrix     |
+| two_e_dm_mo(mo_num,mo_num,mo_num,mo_num) | double precision | two e- density matrix     |
+
+Input:
+| n | integer | mo_num*(mo_num-1)/2 |
+
+Output:
+| v_grad(n) | double precision | the gradient                    |
+| max_elem  | double precision | maximum element of the gradient |
+
+Internal:
+| grad(mo_num,mo_num)                | double precison   | gradient before the tranformation in a vector             |
+| A((mo_num,mo_num)                  | doubre precision  | gradient after the permutations                           |
+| norm                               | double precision  | norm of the gradient                                      |
+| p, q                               | integer           | indexes of the element in the matrix grad                 |
+| i                                  | integer           | index for the tranformation in a vector                   |
+| r, s, t                            | integer           | indexes dor the sums                                      |
+| t1, t2, t3                         | double precision  | t3 = t2 - t1, time to compute the gradient                |
+| t4, t5, t6                         | double precission | t6 = t5 - t4, time to compute each element                |
+| tmp_bi_int_3(mo_num,mo_num,mo_num) | double precision  | 3 indexes temporary array for the bi-electronic integrals |
+| tmp_2rdm_3(mo_num,mo_num,mo_num)   | double precision  | 3 indexes temporary array for the two e- density matrix   |
+| tmp_accu(mo_num,mo_num)            | double precision  | temporary array                                           |
+
+Function:
+| get_two_e_integral | double precision | bi-electronic integrals |
+| dnrm2              | double precision | (Lapack) norm           |
+
+#+BEGIN_SRC f90 :comments org :tangle gradient_opt.irp.f
+subroutine gradient_opt(n,v_grad,max_elem)
+  use omp_lib
+  include 'constants.h'
+
+  implicit none
+  
+  ! Variables
+  
+  ! in
+  integer, intent(in) :: n
+  
+  ! out
+  double precision, intent(out) :: v_grad(n), max_elem
+
+  ! internal
+  double precision, allocatable :: grad(:,:),A(:,:)
+  double precision              :: norm
+  integer                       :: i,p,q,r,s,t
+  double precision              :: t1,t2,t3,t4,t5,t6
+
+  double precision, allocatable :: tmp_accu(:,:)
+  double precision, allocatable :: tmp_bi_int_3(:,:,:), tmp_2rdm_3(:,:,:)
+
+  ! Functions
+  double precision :: get_two_e_integral, dnrm2
+
+  
+  print*,''
+  print*,'---gradient---'
+
+  ! Allocation of shared arrays
+  allocate(grad(mo_num,mo_num),A(mo_num,mo_num)) 
+
+  ! Initialization omp
+  call omp_set_max_active_levels(1)
+
+  !$OMP PARALLEL                                                 &
+      !$OMP PRIVATE(                                             &
+      !$OMP   p,q,r,s,t,                                         &
+      !$OMP   tmp_accu, tmp_bi_int_3, tmp_2rdm_3)                &
+      !$OMP SHARED(grad, one_e_dm_mo, mo_num,mo_one_e_integrals, &
+      !$OMP mo_integrals_map,t4,t5,t6)                           &
+      !$OMP DEFAULT(SHARED)
+ 
+  ! Allocation of private arrays
+  allocate(tmp_accu(mo_num,mo_num))
+  allocate(tmp_bi_int_3(mo_num,mo_num,mo_num))
+  allocate(tmp_2rdm_3(mo_num,mo_num,mo_num))
+#+END_SRC
+  
+** Calculation
+*** Initialization
+#+BEGIN_SRC f90 :comments org :tangle gradient_opt.irp.f
+  !$OMP DO
+  do q = 1, mo_num
+    do p = 1,mo_num
+      grad(p,q) = 0d0
+    enddo
+  enddo
+  !$OMP END DO
+#+END_SRC
+
+*** Term 1
+  
+Without optimization the term 1 is :
+
+do p = 1, mo_num
+  do q = 1, mo_num
+     do r = 1, mo_num
+       grad(p,q) = grad(p,q) &
+               + mo_one_e_integrals(p,r) * one_e_dm_mo(r,q) &
+               - mo_one_e_integrals(r,q) * one_e_dm_mo(p,r)
+    enddo
+  enddo
+enddo
+   
+Since the matrix multiplication A.B is defined like :
+\begin{equation}
+c_{ij} = \sum_k a_{ik}.b_{kj}
+\end{equation}
+The previous equation can be rewritten as a matrix multplication  
+  
+#+BEGIN_SRC f90 :comments org :tangle gradient_opt.irp.f
+  !****************
+  ! Opt first term
+  !****************
+
+  !$OMP MASTER
+  CALL wall_TIME(t4)
+  !$OMP END MASTER
+
+  call dgemm('N','N',mo_num,mo_num,mo_num,1d0,mo_one_e_integrals,&
+  mo_num,one_e_dm_mo,mo_num,0d0,tmp_accu,mo_num)
+  
+  !$OMP DO
+  do q = 1, mo_num
+    do p = 1, mo_num
+
+      grad(p,q) = grad(p,q) + (tmp_accu(p,q) - tmp_accu(q,p))
+
+    enddo
+  enddo 
+  !$OMP END DO
+  
+  !$OMP MASTER
+  CALL wall_TIME(t5)
+  t6 = t5-t4
+  print*,'Gradient, first term (s) :', t6 
+  !$OMP END MASTER
+#+END_SRC
+
+*** Term 2
+ 
+Without optimization the second term is : 
+
+do p = 1, mo_num
+  do q = 1, mo_num 
+    do r = 1, mo_num
+      do s = 1, mo_num
+        do t= 1, mo_num
+
+        grad(p,q) = grad(p,q) &
+                + get_two_e_integral(p,t,r,s,mo_integrals_map) * two_e_dm_mo(r,s,q,t) &
+                - get_two_e_integral(r,s,q,t,mo_integrals_map) * two_e_dm_mo(p,t,r,s)
+       enddo
+      enddo
+    enddo
+  enddo
+enddo
+
+Using the bielectronic integral properties :
+get_two_e_integral(p,t,r,s,mo_integrals_map) = get_two_e_integral(r,s,p,t,mo_integrals_map)
+
+Using the two body matrix properties :
+two_e_dm_mo(p,t,r,s) = two_e_dm_mo(r,s,p,t)
+
+t is one the right, we can put it on the external loop and create 3
+indexes temporary array 
+r,s can be seen as one index
+
+By doing so, a matrix multiplication appears 
+
+#+BEGIN_SRC f90 :comments org :tangle gradient_opt.irp.f
+  !*****************
+  ! Opt second term  
+  !*****************
+
+  !$OMP MASTER
+  CALL wall_TIME(t4)
+  !$OMP END MASTER 
+ 
+  !$OMP DO
+  do t = 1, mo_num
+    
+    do p = 1, mo_num
+      do s = 1, mo_num
+        do r = 1, mo_num
+            
+          tmp_bi_int_3(r,s,p) = get_two_e_integral(r,s,p,t,mo_integrals_map)
+         
+        enddo
+      enddo
+    enddo
+
+    do q = 1, mo_num
+      do s = 1, mo_num
+        do r = 1, mo_num
+             
+           tmp_2rdm_3(r,s,q) = two_e_dm_mo(r,s,q,t)
+  
+        enddo
+      enddo
+    enddo
+
+    call dgemm('T','N',mo_num,mo_num,mo_num*mo_num,1d0,tmp_bi_int_3,&
+      mo_num*mo_num,tmp_2rdm_3,mo_num*mo_num,0d0,tmp_accu,mo_num)
+
+    !$OMP CRITICAL   
+    do q = 1, mo_num
+      do p = 1, mo_num
+
+        grad(p,q) = grad(p,q) + tmp_accu(p,q) - tmp_accu(q,p)
+
+      enddo
+    enddo
+    !$OMP END CRITICAL
+
+  enddo
+  !$OMP END DO
+
+  !$OMP MASTER
+  CALL wall_TIME(t5)
+  t6 = t5-t4
+  print*,'Gradient second term (s) : ', t6
+  !$OMP END MASTER  
+#+END_SRC
+
+*** Deallocation of private arrays
+#+BEGIN_SRC f90 :comments org :tangle gradient_opt.irp.f
+  deallocate(tmp_bi_int_3,tmp_2rdm_3,tmp_accu)
+
+  !$OMP END PARALLEL
+
+  call omp_set_max_active_levels(4)
+#+END_SRC
+
+*** Permutation, 2D matrix -> vector, transformation
+In addition there is a permutation in the gradient formula :
+\begin{equation}
+P_{pq} = 1 - (p <-> q) 
+\end{equation}
+
+We need a vector to use the gradient. Here the gradient is a
+antisymetric matrix so we can transform it in a vector of length
+mo_num*(mo_num-1)/2.
+
+Here we do these two things at the same time.
+
+#+BEGIN_SRC f90 :comments org :tangle gradient_opt.irp.f
+  do i=1,n
+    call vec_to_mat_index(i,p,q)
+    v_grad(i)=(grad(p,q) - grad(q,p))
+  enddo  
+
+  ! Debug, diplay the vector containing the gradient elements 
+  if (debug) then  
+    print*,'Vector containing the gradient :'
+    write(*,'(100(F10.5))') v_grad(1:n)
+  endif  
+#+END_SRC
+
+*** Norm of the gradient
+The norm can be useful.
+#+BEGIN_SRC f90 :comments org :tangle gradient_opt.irp.f
+  norm = dnrm2(n,v_grad,1)
+  print*, 'Gradient norm : ', norm
+#+END_SRC
+
+*** Maximum element in the gradient
+The maximum element in the gradient is very important for the
+convergence criterion of the Newton method.
+
+#+BEGIN_SRC f90 :comments org :tangle gradient_opt.irp.f
+  ! Max element of the gradient
+  max_elem = 0d0
+  do i = 1, n
+    if (ABS(v_grad(i)) > ABS(max_elem)) then
+      max_elem = v_grad(i)
+    endif
+  enddo
+
+  print*,'Max element in the gradient :', max_elem  
+
+  ! Debug, display the matrix containting the gradient elements
+  if (debug) then
+    ! Matrix gradient
+    A = 0d0
+    do q=1,mo_num
+      do p=1,mo_num
+        A(p,q) = grad(p,q) - grad(q,p)
+      enddo
+    enddo
+    print*,'Matrix containing the gradient :'
+    do i = 1, mo_num
+      write(*,'(100(F10.5))') A(i,1:mo_num)
+    enddo
+  endif
+#+END_SRC
+
+*** Deallocation of shared arrays and end
+#+BEGIN_SRC f90 :comments org :tangle gradient_opt.irp.f
+  deallocate(grad,A)
+
+  print*,'---End gradient---'
+
+  end subroutine
+
+#+END_SRC
+  

src/mo_optimization/org/my_providers.org

@@ -0,0 +1,308 @@
+* Providers
+** Dimensions of MOs
+
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ integer, n_mo_dim ]
+  implicit none
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim = mo_num*(mo_num-1)/2
+
+END_PROVIDER
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ integer, n_mo_dim_core ]
+  implicit none 
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of core MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim_core = dim_list_core_orb*(dim_list_core_orb-1)/2
+
+END_PROVIDER
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ integer, n_mo_dim_act ]
+  implicit none
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of active MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim_act = dim_list_act_orb*(dim_list_act_orb-1)/2
+
+END_PROVIDER
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ integer, n_mo_dim_inact ]
+  implicit none 
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of inactive MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim_inact = dim_list_inact_orb*(dim_list_inact_orb-1)/2
+
+END_PROVIDER
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ integer, n_mo_dim_virt ]
+  implicit none 
+  BEGIN_DOC
+  ! Number of different pairs (i,j) of virtual MOs we can build,
+  ! with i>j
+  END_DOC
+
+  n_mo_dim_virt = dim_list_virt_orb*(dim_list_virt_orb-1)/2
+
+END_PROVIDER
+#+END_SRC
+
+** Energies/criterions
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_st_av_energy ]
+  implicit none
+  BEGIN_DOC
+  ! State average CI energy
+  END_DOC
+
+  !call update_st_av_ci_energy(my_st_av_energy)
+  call state_average_energy(my_st_av_energy)
+
+END_PROVIDER
+#+END_SRC
+
+** Gradient/hessian
+*** Orbital optimization
+**** With all the MOs
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_gradient_opt, (n_mo_dim) ]
+&BEGIN_PROVIDER [ double precision, my_CC1_opt ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the energy with respect to the MO rotations, for all the MOs.
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision :: norm_grad
+
+  PROVIDE mo_two_e_integrals_in_map
+
+  call gradient_opt(n_mo_dim, my_gradient_opt, my_CC1_opt, norm_grad)
+
+END_PROVIDER
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_hessian_opt, (n_mo_dim, n_mo_dim) ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the energy with respect to the MO rotations, for all the MOs.
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision, allocatable :: h_f(:,:,:,:)
+
+  PROVIDE mo_two_e_integrals_in_map
+
+  allocate(h_f(mo_num, mo_num, mo_num, mo_num))
+
+  call hessian_list_opt(n_mo_dim, my_hessian_opt, h_f)
+
+END_PROVIDER
+#+END_SRC
+
+**** With the list of active MOs
+Can be generalized to any mo_class by changing the list/dimension
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_gradient_list_opt, (n_mo_dim_act) ]
+&BEGIN_PROVIDER [ double precision, my_CC2_opt ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the energy with respect to the MO rotations, only for the active MOs !
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision :: norm_grad
+
+  PROVIDE mo_two_e_integrals_in_map !one_e_dm_mo two_e_dm_mo mo_one_e_integrals 
+
+  call gradient_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_gradient_list_opt, my_CC2_opt, norm_grad)
+
+END_PROVIDER
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_hessian_list_opt, (n_mo_dim_act, n_mo_dim_act) ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the energy with respect to the MO rotations, only for the active MOs !
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision, allocatable :: h_f(:,:,:,:)
+
+  PROVIDE mo_two_e_integrals_in_map
+
+  allocate(h_f(dim_list_act_orb, dim_list_act_orb, dim_list_act_orb, dim_list_act_orb))
+
+  call hessian_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_hessian_list_opt, h_f)
+
+END_PROVIDER
+#+END_SRC
+
+*** Orbital localization
+**** Gradient
+***** Core MOs
+#+BEGIN_SRC f90 :comments org 
+!:tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_gradient_loc_core, (n_mo_dim_core) ]
+&BEGIN_PROVIDER [ double precision, my_CC_loc_core ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the MO localization with respect to the MO rotations for the core MOs
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision :: norm_grad
+
+  !PROVIDE something ?
+
+  call gradient_localization(n_mo_dim_core, dim_list_core_orb, list_core, my_gradient_loc_core, my_CC_loc_core , norm_grad)
+
+END_PROVIDER
+#+END_SRC
+
+***** Active MOs 
+#+BEGIN_SRC f90 :comments org 
+!:tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_gradient_loc_act, (n_mo_dim_act) ]
+&BEGIN_PROVIDER [ double precision, my_CC_loc_act ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the MO localization with respect to the MO rotations for the active MOs
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision :: norm_grad
+
+  !PROVIDE something ?
+
+  call gradient_localization(n_mo_dim_act, dim_list_act_orb, list_act, my_gradient_loc_act, my_CC_loc_act , norm_grad)
+
+END_PROVIDER
+#+END_SRC
+
+***** Inactive MOs
+#+BEGIN_SRC f90 :comments org !
+:tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_gradient_loc_inact, (n_mo_dim_inact) ]
+&BEGIN_PROVIDER [ double precision, my_CC_loc_inact ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the MO localization with respect to the MO rotations for the inactive MOs
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision :: norm_grad
+
+  !PROVIDE something ?
+
+  call gradient_localization(n_mo_dim_inact, dim_list_inact_orb, list_inact, my_gradient_loc_inact, my_CC_loc_inact , norm_grad)
+
+END_PROVIDER
+#+END_SRC
+
+***** Virtual MOs
+#+BEGIN_SRC f90 :comments org 
+!:tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_gradient_loc_virt, (n_mo_dim_virt) ]
+&BEGIN_PROVIDER [ double precision, my_CC_loc_virt ]
+  implicit none
+  BEGIN_DOC
+  ! - Gradient of the MO localization with respect to the MO rotations for the virtual MOs
+  ! - Maximal element of the gradient in absolute value 
+  END_DOC
+
+  double precision :: norm_grad
+
+  !PROVIDE something ?
+
+  call gradient_localization(n_mo_dim_virt, dim_list_virt_orb, list_virt, my_gradient_loc_virt, my_CC_loc_virt , norm_grad)
+
+END_PROVIDER
+#+END_SRC
+
+**** Hessian
+***** Core MOs
+#+BEGIN_SRC f90 :comments org 
+!:tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_hessian_loc_core, (n_mo_dim_core) ]
+  implicit none
+  BEGIN_DOC
+  ! - Hessian of the MO localization with respect to the MO rotations for the core MOs
+  END_DOC
+
+  !PROVIDE something ?
+
+  call hessian_localization(n_mo_dim_core, dim_list_core_orb, list_core, my_hessian_loc_core)
+
+END_PROVIDER
+#+END_SRC
+
+***** Active MOs
+#+BEGIN_SRC f90 :comments org 
+!:tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_hessian_loc_act, (n_mo_dim_act) ]
+  implicit none
+  BEGIN_DOC
+  ! - Hessian of the MO localization with respect to the MO rotations for the active MOs
+  END_DOC
+
+  !PROVIDE something ?
+
+  call hessian_localization(n_mo_dim_act, dim_list_act_orb, list_act, my_hessian_loc_act)
+
+END_PROVIDER
+#+END_SRC
+
+***** Inactive MOs
+#+BEGIN_SRC f90 :comments org 
+!:tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_hessian_loc_inact, (n_mo_dim_inact) ]
+  implicit none
+  BEGIN_DOC
+  ! - Hessian of the MO localization with respect to the MO rotations for the inactive MOs
+  END_DOC
+
+  !PROVIDE something ?
+
+  call hessian_localization(n_mo_dim_inact, dim_list_inact_orb, list_inact, my_hessian_loc_inact)
+
+END_PROVIDER
+#+END_SRC
+
+***** Virtual MOs
+#+BEGIN_SRC f90 :comments org 
+!:tangle my_providers.irp.f
+BEGIN_PROVIDER [ double precision, my_hessian_loc_virt, (n_mo_dim_virt) ]
+  implicit none
+  BEGIN_DOC
+  ! - Hessian of the MO localization with respect to the MO rotations for the virtual MOs
+  END_DOC
+
+  !PROVIDE something ?
+
+  call hessian_localization(n_mo_dim_virt, dim_list_virt_orb, list_virt, my_hessian_loc_virt)
+
+END_PROVIDER
+#+END_SRC
+

src/mo_optimization/org/optimization.org

@@ -0,0 +1,91 @@
+#+BEGIN_SRC f90 :comments org :tangle optimization.irp.f
+program optimization
+  
+  read_wf = .true. ! must be True for the orbital optimization !!!
+  TOUCH read_wf 
+  call run_optimization
+
+end
+#+END_SRC
+
+#+BEGIN_SRC f90 :comments org :tangle optimization.irp.f
+subroutine run_optimization
+
+  implicit none
+
+  double precision :: e_cipsi, e_opt, delta_e
+  integer :: nb_iter,i
+  logical :: not_converged
+  character (len=100) :: filename
+
+  PROVIDE psi_det psi_coef mo_two_e_integrals_in_map
+
+  not_converged = .True.
+  nb_iter = 0
+
+  ! To start from the wf
+  N_det_max = max(n_det,5)
+  TOUCH N_det_max
+
+  open(unit=10, file=trim(ezfio_filename)//'/mo_optimization/result_opt')
+  write(10,*) "   Ndet        E_cipsi         E_opt          Delta_e"
+  call state_average_energy(e_cipsi)
+  write(10,'(I10, 3F15.7)') n_det, e_cipsi, e_cipsi, 0d0
+  close(10)
+
+  do while (not_converged)
+      print*,''
+      print*,'======================'
+      print*,' Cipsi step:', nb_iter
+      print*,'======================'
+      print*,'' 
+      print*,'********** cipsi step **********'
+      ! cispi calculation
+      call run_stochastic_cipsi
+
+      ! State average energy after the cipsi step
+      call state_average_energy(e_cipsi)
+
+      print*,''
+      print*,'********** optimization step **********'
+      ! orbital optimization
+      call run_orb_opt_trust_v2
+
+      ! State average energy after the orbital optimization
+      call state_average_energy(e_opt)
+
+      print*,''
+      print*,'********** diff step **********'
+      ! Gain in energy
+      delta_e = e_opt - e_cipsi
+      print*, 'Gain in energy during the orbital optimization:', delta_e
+
+      open(unit=10, file=trim(ezfio_filename)//'/mo_optimization/result_opt', position='append')
+      write(10,'(I10, 3F15.7)') n_det, e_cipsi, e_opt, delta_e
+      close(10)
+
+      ! Exit
+      if (delta_e > 1d-12) then
+          print*, 'WARNING, something wrong happened'
+          print*, 'The gain (delta_e) in energy during the optimization process'
+          print*, 'is > 0, but it must be < 0'
+          print*, 'The program will exit'
+          exit
+      endif
+
+      if (n_det > n_det_max_opt) then
+          print*, 'The number of determinants in the wf > n_det_max_opt'
+          print*, 'The program will exit'
+          exit
+      endif
+      
+      ! To double the number of determinants in the wf
+      N_det_max = int(dble(n_det * 2)*0.9)
+      TOUCH N_det_max
+
+      nb_iter = nb_iter + 1
+  enddo
+
+end
+
+#+END_SRC

src/mo_optimization/org/orb_opt_trust_v2.org

@@ -0,0 +1,349 @@
+* Orbital optimization program
+
+This is an optimization program for molecular orbitals. It produces
+orbital rotations in order to lower the energy of a truncated wave
+function.  
+This program just optimize the orbitals for a fixed number of
+determinants. This optimization process must be repeated for different
+number of determinants.
+
+#+BEGIN_SRC f90 :comments org :tangle orb_opt.irp.f
+#+END_SRC
+
+* Main program : orb_opt_trust
+
+#+BEGIN_SRC f90 :comments org :tangle orb_opt.irp.f 
+program orb_opt
+  read_wf = .true. ! must be True for the orbital optimization !!!
+  TOUCH read_wf
+  io_mo_two_e_integrals = 'None'
+  TOUCH io_mo_two_e_integrals
+  call run_orb_opt_trust_v2
+end
+#+END_SRC
+
+* Subroutine : run_orb_opt_trust
+
+#+BEGIN_SRC f90 :comments org :tangle run_orb_opt_trust_v2.irp.f
+
+#+END_SRC
+
+Subroutine to optimize the MOs using a trust region algorithm:
+- choice of the method
+- initialization
+- optimization until convergence
+
+The optimization use the trust region algorithm, the different parts
+are explained in the corresponding subroutine files.
+
+qp_edit:
+| thresh_opt_max_elem_grad |
+| optimization_max_nb_iter |
+| optimization_method      |
+
+Provided:
+| mo_num                         | integer          | number of MOs                  |
+| ao_num                         | integer          | number of AOs                  |
+| N_states                       | integer          | number of states               |
+| ci_energy(N_states)            | double precision | CI energies                    |
+| state_average_weight(N_states) | double precision | Weight of the different states |
+
+Variables:
+| m                         | integer          | number of active MOs                               |
+| tmp_n                     | integer          | m*(m-1)/2, number of MO parameters                 |
+| tmp_n2                    | integer          | m*(m-1)/2 or 1 if the hessian is diagonal          |
+| v_grad(tmp_n)             | double precision | gradient                                           |
+| H(tmp_n,tmp_n)            | double precision | hessian (2D)                                       |
+| h_f(m,m,m,m)              | double precision | hessian (4D)                                       |
+| e_val(m)                  | double precision | eigenvalues of the hessian                         |
+| w(m,m)                    | double precision | eigenvectors of the hessian                        |
+| x(m)                      | double precision | step given by the trust region                     |
+| m_x(m,m)                  | double precision | step given by the trust region after               |
+| tmp_R(m,m)                | double precision | rotation matrix for active MOs                     |
+| R(mo_num,mo_num)          | double precision | full rotation matrix                               |
+| prev_mos(ao_num,mo_num)   | double precision | previous MOs (before the rotation)                 |
+| new_mos(ao_num,mo_num)    | double precision | new MOs (after the roration)                       |
+| delta                     | double precision | radius of the trust region                         |
+| rho                       | double precision | agreement between the model and the exact function |
+| max_elem                  | double precision | maximum element in the gradient                    |
+| i                         | integer          | index                                              |
+| tmp_i,tmp_j               | integer          | indexes in the subspace containing only            |
+|                           |                  | the active MOs                                     |
+| converged                 | logical          | convergence of the algorithm                       |
+| cancel_step               | logical          | if the step must be cancelled                      |
+| nb_iter                   | integer          | number of iterations (accepted)                    |
+| nb_diag                   | integer          | number of diagonalizations of the CI matrix        |
+| nb_cancel                 | integer          | number of cancelled steps for the actual iteration |
+| nb_cancel_tot             | integer          | total number of cancel steps                       |
+| info                      | integer          | if 0 ok, else problem in the diagonalization of    |
+|                           |                  | the hessian with the Lapack routine                |
+| criterion                 | double precision | energy at a given step                             |
+| prev_criterion            | double precision | energy before the rotation                         |
+| criterion_model           | double precision | estimated energy after the rotation using          |
+|                           |                  | a Taylor series                                    |
+| must_exit                 | logical          | To exit the trust region algorithm when            |
+|                           |                  | criterion - criterion_model is too small           |
+| enforce_step_cancellation | logical          | To force the cancellation of the step if the       |
+|                           |                  | error in the rotation matrix is too large          |
+
+#+BEGIN_SRC f90 :comments org :tangle run_orb_opt_trust_v2.irp.f
+subroutine run_orb_opt_trust_v2
+
+  include 'constants.h'
+
+  implicit none
+
+  BEGIN_DOC
+  ! Orbital optimization
+  END_DOC
+
+  ! Variables
+
+  double precision, allocatable :: R(:,:)
+  double precision, allocatable :: H(:,:),h_f(:,:,:,:)
+  double precision, allocatable :: v_grad(:)
+  double precision, allocatable :: prev_mos(:,:),new_mos(:,:)
+  integer                       :: info
+  integer                       :: n
+  integer                       :: i,j,p,q,k
+  double precision              :: max_elem_grad, delta, rho, norm_grad, normalization_factor
+  logical                       :: cancel_step
+  integer                       :: nb_iter, nb_diag, nb_cancel, nb_cancel_tot, nb_sub_iter
+  double precision              :: t1, t2, t3
+  double precision              :: prev_criterion, criterion, criterion_model
+  logical                       :: not_converged, must_exit, enforce_step_cancellation
+  integer                       :: m, tmp_n, tmp_i, tmp_j, tmp_k, tmp_n2
+  integer,allocatable           :: tmp_list(:), key(:)
+  double precision, allocatable :: tmp_m_x(:,:),tmp_R(:,:), tmp_x(:), W(:,:), e_val(:)
+
+  PROVIDE mo_two_e_integrals_in_map ci_energy psi_det psi_coef
+#+END_SRC
+   
+** Allocation
+
+#+BEGIN_SRC f90 :comments org :tangle run_orb_opt_trust_v2.irp.f
+  allocate(R(mo_num,mo_num))  ! rotation matrix
+  allocate(prev_mos(ao_num,mo_num), new_mos(ao_num,mo_num)) ! old and new MOs
+  
+  ! Definition of m and tmp_n
+  m = dim_list_act_orb
+  tmp_n = m*(m-1)/2
+  
+  allocate(tmp_list(m))
+  allocate(tmp_R(m,m), tmp_m_x(m,m), tmp_x(tmp_n))
+  allocate(e_val(tmp_n),key(tmp_n),v_grad(tmp_n))
+
+#+END_SRC
+
+** Method
+   There are three different methods : 
+   - the "full" hessian, which uses all the elements of the hessian
+     matrix"
+   - the "diagonal" hessian, which uses only the diagonal elements of the
+     hessian
+   - without the hessian (hessian = identity matrix) 
+
+#+BEGIN_SRC f90 :comments org :tangle run_orb_opt_trust_v2.irp.f
+  !Display the method
+   print*, 'Method :', optimization_method
+  if (optimization_method == 'full') then 
+    print*, 'Full hessian'
+    allocate(H(tmp_n,tmp_n), h_f(m,m,m,m),W(tmp_n,tmp_n))
+    tmp_n2 = tmp_n
+  elseif (optimization_method == 'diag') then
+    print*,'Diagonal hessian'
+    allocate(H(tmp_n,1),W(tmp_n,1))
+    tmp_n2 = 1
+  elseif (optimization_method == 'none') then
+    print*,'No hessian'
+    allocate(H(tmp_n,1),W(tmp_n,1))
+    tmp_n2 = 1
+  else
+    print*,'Unknown optimization_method, please select full, diag or none'
+    call abort
+  endif
+  print*, 'Absolute value of the hessian:', absolute_eig
+#+END_SRC
+
+** Calculations
+*** Algorithm
+
+Here is the main algorithm of the optimization:
+- First of all we initialize some parameters and we compute the
+  criterion (the ci energy) before doing any MO rotations
+- We compute the gradient and the hessian for the active MOs
+- We diagonalize the hessian
+- We compute a step and loop to reduce the radius of the
+  trust region (and the size of the step by the way) until the step is
+  accepted 
+- We repeat the process until the convergence 
+  NB: the convergence criterion can be changed
+
+#+BEGIN_SRC f90 :comments org :tangle run_orb_opt_trust_v2.irp.f
+  ! Loop until the convergence of the optimization
+  ! call diagonalize_ci 
+
+  !### Initialization ###
+  nb_iter = 0
+  rho = 0.5d0
+  not_converged = .True.
+  tmp_list = list_act ! Optimization of the active MOs
+  nb_cancel_tot = 0
+
+  ! Renormalization of the weights of the states
+  call state_weight_normalization
+
+  ! Compute the criterion before the loop
+  call state_average_energy(prev_criterion)
+
+  do while (not_converged)
+    print*,''
+    print*,'******************'
+    print*,'Iteration', nb_iter
+    print*,'******************'
+    print*,''
+
+    ! Gradient
+    call gradient_list_opt(tmp_n, m, tmp_list, v_grad, max_elem_grad, norm_grad)
+    
+    ! Hessian
+    if (optimization_method == 'full') then
+      ! Full hessian
+      call hessian_list_opt(tmp_n, m, tmp_list, H, h_f)
+
+      ! Diagonalization of the hessian 
+      call diagonalization_hessian(tmp_n, H, e_val, w)
+
+    elseif (optimization_method == 'diag') then
+      ! Diagonal hessian 
+      call diag_hessian_list_opt(tmp_n, m, tmp_list, H)
+    else
+      ! Identity matrix 
+      do tmp_i = 1, tmp_n
+        H(tmp_i,1) = 1d0
+      enddo
+    endif
+
+    if (optimization_method /= 'full') then
+      ! Sort
+      do tmp_i = 1, tmp_n
+        key(tmp_i) = tmp_i
+        e_val(tmp_i) = H(tmp_i,1)
+      enddo
+      call dsort(e_val,key,tmp_n)
+
+      ! Eigenvalues and eigenvectors
+      do tmp_i = 1, tmp_n
+        w(tmp_i,1) = dble(key(tmp_i))
+      enddo
+      
+    endif
+ 
+    ! Init before the internal loop
+    cancel_step = .True. ! To enter in the loop just after 
+    nb_cancel = 0
+    nb_sub_iter = 0
+
+    ! Loop to reduce the trust radius until the criterion decreases and rho >= thresh_rho
+    do while (cancel_step)
+      print*,''
+      print*,'-----------------------------'
+      print*,'Iteration:    ', nb_iter
+      print*,'Sub iteration:', nb_sub_iter
+      print*,'Max elem grad:', max_elem_grad
+      print*,'-----------------------------'
+
+      ! Hessian,gradient,Criterion -> x 
+      call trust_region_step_w_expected_e(tmp_n,tmp_n2,H,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,tmp_x,must_exit) 
+      
+      if (must_exit) then
+        print*,'step_in_trust_region sends: Exit'
+        exit
+      endif
+
+      ! 1D tmp -> 2D tmp 
+      call vec_to_mat_v2(tmp_n, m, tmp_x, tmp_m_x)
+
+      ! Rotation matrix for the active MOs
+      call rotation_matrix(tmp_m_x, m, tmp_R, m, m, info, enforce_step_cancellation)
+
+      ! Security to ensure an unitary transformation
+      if (enforce_step_cancellation) then
+        print*, '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(m, tmp_list, tmp_R, R)
+    
+      ! MO rotations
+      call apply_mo_rotation(R, prev_mos)   
+
+      ! Update of the energy before the diagonalization of the hamiltonian
+      call clear_mo_map
+      TOUCH mo_coef psi_det psi_coef ci_energy two_e_dm_mo 
+      call state_average_energy(criterion)
+
+      ! 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 necessary
+      if (cancel_step) then
+        mo_coef = prev_mos
+        call save_mos()
+        nb_cancel = nb_cancel + 1
+        nb_cancel_tot = nb_cancel_tot + 1
+      else
+        ! Diagonalization of the hamiltonian
+        FREE ci_energy! To enforce the recomputation
+        call diagonalize_ci
+        call save_wavefunction_unsorted
+
+        ! Energy obtained after the diagonalization of the CI matrix
+        call state_average_energy(prev_criterion)
+      endif
+
+      nb_sub_iter = nb_sub_iter + 1
+    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
+
+    ! External loop exit conditions
+    if (DABS(max_elem_grad) < thresh_opt_max_elem_grad) then
+      print*,'Converged: DABS(max_elem_grad) < thresh_opt_max_elem_grad'
+      not_converged = .False.
+    endif
+    if (nb_iter >= optimization_max_nb_iter) then
+      print*,'Not converged: nb_iter >= optimization_max_nb_iter'
+      not_converged = .False. 
+    endif
+
+    if (.not. not_converged) then
+      print*,'#############################'
+      print*,'   End of the optimization'
+      print*,'#############################'
+    endif
+  enddo
+
+#+END_SRC
+    
+** Deallocation, end
+
+#+BEGIN_SRC f90 :comments org :tangle run_orb_opt_trust_v2.irp.f
+  deallocate(v_grad,H,R,W,e_val)
+  deallocate(prev_mos,new_mos)
+  if (optimization_method == 'full') then
+    deallocate(h_f)
+  endif
+
+end
+#+END_SRC
+

src/mo_optimization/org/state_average_energy.org

@@ -0,0 +1,73 @@
+* State average energy
+
+Calculation of the state average energy from the integrals and the
+density matrices.
+
+\begin{align*}
+E = \sum_{ij} h_{ij} \gamma_{ij} + \frac{1}{2} v_{ij}^{kl} \Gamma_{ij}^{kl}
+\end{align*}
+$h_{ij}$: mono-electronic integral
+$\gamma_{ij}$: one electron density matrix
+$v_{ij}^{kl}$: bi-electronic integral
+$\Gamma_{ij}^{kl}$: two electrons density matrix
+
+TODO: OMP version
+
+PROVIDED:
+| mo_one_e_integrals | double precision | mono-electronic integrals    |
+| get_two_e_integral | double precision | bi-electronic integrals      |
+| one_e_dm_mo        | double precision | one electron density matrix  |
+| two_e_dm_mo        | double precision | two electrons density matrix |
+| nuclear_repulsion  | double precision | nuclear repulsion            |
+| mo_num             | integer          | number of MOs                |
+
+Output:
+| energy | double precision | state average energy |
+
+Internal:
+| mono_e  | double precision | mono-electronic energy       |
+| bi_e    | double precision | bi-electronic energy         |
+| i,j,k,l | integer          | indexes to loop over the MOs |
+
+#+BEGIN_SRC f90 :comments org :tangle state_average_energy.irp.f
+subroutine state_average_energy(energy)
+
+  implicit none
+
+  double precision, intent(out) :: energy
+
+  double precision :: get_two_e_integral
+  double precision :: mono_e, bi_e
+  integer :: i,j,k,l
+  
+  ! mono electronic part
+  mono_e = 0d0
+  do j = 1, mo_num
+    do i = 1, mo_num
+      mono_e =  mono_e + mo_one_e_integrals(i,j) * one_e_dm_mo(i,j)
+    enddo
+  enddo
+
+  ! bi electronic part
+  bi_e = 0d0
+  do l = 1, mo_num
+    do k = 1, mo_num
+      do j = 1, mo_num
+        do i = 1, mo_num
+          bi_e = bi_e + get_two_e_integral(i,j,k,l,mo_integrals_map) * two_e_dm_mo(i,j,k,l)
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ! State average energy
+  energy = mono_e + 0.5d0 * bi_e + nuclear_repulsion
+
+  ! Check
+  !call print_energy_components
+  
+  print*,'State average energy:', energy
+  !print*,ci_energy
+
+end
+#+END_SRC

src/mo_optimization/org/state_weight_normalization.org

@@ -0,0 +1,31 @@
+#+BEGIN_SRC f90 :comments org :tangle state_weight_normalization.irp.f
+subroutine state_weight_normalization
+
+  implicit none
+
+  BEGIN_DOC
+  ! Renormalization of the state weights or enforcing state average
+  ! weights for orbital optimization
+  END_DOC
+
+  integer          :: i
+  double precision :: normalization_factor
+  
+  ! To normalize the sum of the state weights
+  normalization_factor = 0d0
+  do i = 1, N_states
+    normalization_factor = normalization_factor + state_average_weight(i)
+  enddo
+  normalization_factor = 1d0 / normalization_factor
+  
+  do i = 1, N_states
+    state_average_weight(i) = state_average_weight(i) * normalization_factor
+  enddo
+  TOUCH state_average_weight
+  
+  print*, 'Number of states:', N_states
+  print*, 'State average weights:'
+  print*, state_average_weight(:)
+
+end
+#+END_SRC

src/mo_optimization/org/update_parameters.org

@@ -0,0 +1,16 @@
+Subroutine toupdate the parameters.
+Ex: TOUCH mo_coef ...
+
+#+BEGIN_SRC f90 :comments org :tangle update_parameters.irp.f
+subroutine update_parameters()
+
+  implicit none
+
+  !### TODO
+  ! Touch yours parameters
+  call clear_mo_map
+  TOUCH mo_coef psi_det psi_coef
+  call diagonalize_ci
+  call save_wavefunction_unsorted
+end
+#+END_SRC

src/mo_optimization/org/update_st_av_ci_energy.org

@@ -0,0 +1,26 @@
+* Update the CI state average energy
+
+Computes the state average energy 
+\begin{align*}
+E =\sum_{i=1}^{N_{states}} E_i . w_i
+\end{align*}
+
+$E_i$: energy of state i
+$w_i$: weight of state i
+
+#+BEGIN_SRC f90 :comments org :tangle update_st_av_ci_energy.irp.f
+subroutine update_st_av_ci_energy(energy)
+
+  implicit none
+  
+  double precision, intent(out) :: energy
+  integer :: i
+
+  energy = 0d0
+  do i = 1, N_states
+    energy = energy + ci_energy(i) * state_average_weight(i)
+  enddo
+  
+  print*, 'ci_energy :', energy
+end
+#+END_SRC

src/mo_optimization/run_orb_opt_trust_v2.irp.f

@@ -0,0 +1,317 @@
+! Subroutine : run_orb_opt_trust
+
+
+
+
+
+
+! Subroutine to optimize the MOs using a trust region algorithm:
+! - choice of the method
+! - initialization
+! - optimization until convergence
+
+! The optimization use the trust region algorithm, the different parts
+! are explained in the corresponding subroutine files.
+
+! qp_edit:
+! | thresh_opt_max_elem_grad |
+! | optimization_max_nb_iter |
+! | optimization_method      |
+
+! Provided:
+! | mo_num                         | integer          | number of MOs                  |
+! | ao_num                         | integer          | number of AOs                  |
+! | N_states                       | integer          | number of states               |
+! | ci_energy(N_states)            | double precision | CI energies                    |
+! | state_average_weight(N_states) | double precision | Weight of the different states |
+
+! Variables:
+! | m                         | integer          | number of active MOs                               |
+! | tmp_n                     | integer          | m*(m-1)/2, number of MO parameters                 |
+! | tmp_n2                    | integer          | m*(m-1)/2 or 1 if the hessian is diagonal          |
+! | v_grad(tmp_n)             | double precision | gradient                                           |
+! | H(tmp_n,tmp_n)            | double precision | hessian (2D)                                       |
+! | h_f(m,m,m,m)              | double precision | hessian (4D)                                       |
+! | e_val(m)                  | double precision | eigenvalues of the hessian                         |
+! | w(m,m)                    | double precision | eigenvectors of the hessian                        |
+! | x(m)                      | double precision | step given by the trust region                     |
+! | m_x(m,m)                  | double precision | step given by the trust region after               |
+! | tmp_R(m,m)                | double precision | rotation matrix for active MOs                     |
+! | R(mo_num,mo_num)          | double precision | full rotation matrix                               |
+! | prev_mos(ao_num,mo_num)   | double precision | previous MOs (before the rotation)                 |
+! | new_mos(ao_num,mo_num)    | double precision | new MOs (after the roration)                       |
+! | delta                     | double precision | radius of the trust region                         |
+! | rho                       | double precision | agreement between the model and the exact function |
+! | max_elem                  | double precision | maximum element in the gradient                    |
+! | i                         | integer          | index                                              |
+! | tmp_i,tmp_j               | integer          | indexes in the subspace containing only            |
+! |                           |                  | the active MOs                                     |
+! | converged                 | logical          | convergence of the algorithm                       |
+! | cancel_step               | logical          | if the step must be cancelled                      |
+! | nb_iter                   | integer          | number of iterations (accepted)                    |
+! | nb_diag                   | integer          | number of diagonalizations of the CI matrix        |
+! | nb_cancel                 | integer          | number of cancelled steps for the actual iteration |
+! | nb_cancel_tot             | integer          | total number of cancel steps                       |
+! | info                      | integer          | if 0 ok, else problem in the diagonalization of    |
+! |                           |                  | the hessian with the Lapack routine                |
+! | criterion                 | double precision | energy at a given step                             |
+! | prev_criterion            | double precision | energy before the rotation                         |
+! | criterion_model           | double precision | estimated energy after the rotation using          |
+! |                           |                  | a Taylor series                                    |
+! | must_exit                 | logical          | To exit the trust region algorithm when            |
+! |                           |                  | criterion - criterion_model is too small           |
+! | enforce_step_cancellation | logical          | To force the cancellation of the step if the       |
+! |                           |                  | error in the rotation matrix is too large          |
+
+
+subroutine run_orb_opt_trust_v2
+
+  include 'constants.h'
+
+  implicit none
+
+  BEGIN_DOC
+  ! Orbital optimization
+  END_DOC
+
+  ! Variables
+
+  double precision, allocatable :: R(:,:)
+  double precision, allocatable :: H(:,:),h_f(:,:,:,:)
+  double precision, allocatable :: v_grad(:)
+  double precision, allocatable :: prev_mos(:,:),new_mos(:,:)
+  integer                       :: info
+  integer                       :: n
+  integer                       :: i,j,p,q,k
+  double precision              :: max_elem_grad, delta, rho, norm_grad, normalization_factor
+  logical                       :: cancel_step
+  integer                       :: nb_iter, nb_diag, nb_cancel, nb_cancel_tot, nb_sub_iter
+  double precision              :: t1, t2, t3
+  double precision              :: prev_criterion, criterion, criterion_model
+  logical                       :: not_converged, must_exit, enforce_step_cancellation
+  integer                       :: m, tmp_n, tmp_i, tmp_j, tmp_k, tmp_n2
+  integer,allocatable           :: tmp_list(:), key(:)
+  double precision, allocatable :: tmp_m_x(:,:),tmp_R(:,:), tmp_x(:), W(:,:), e_val(:)
+
+  PROVIDE mo_two_e_integrals_in_map ci_energy psi_det psi_coef
+
+! Allocation
+
+
+allocate(R(mo_num,mo_num))  ! rotation matrix
+allocate(prev_mos(ao_num,mo_num), new_mos(ao_num,mo_num)) ! old and new MOs
+
+! Definition of m and tmp_n
+m = dim_list_act_orb
+tmp_n = m*(m-1)/2
+
+allocate(tmp_list(m))
+allocate(tmp_R(m,m), tmp_m_x(m,m), tmp_x(tmp_n))
+allocate(e_val(tmp_n),key(tmp_n),v_grad(tmp_n))
+
+! Method
+!    There are three different methods : 
+!    - the "full" hessian, which uses all the elements of the hessian
+!      matrix"
+!    - the "diagonal" hessian, which uses only the diagonal elements of the
+!      hessian
+!    - without the hessian (hessian = identity matrix) 
+
+
+!Display the method
+ print*, 'Method :', optimization_method
+if (optimization_method == 'full') then 
+  print*, 'Full hessian'
+  allocate(H(tmp_n,tmp_n), h_f(m,m,m,m),W(tmp_n,tmp_n))
+  tmp_n2 = tmp_n
+elseif (optimization_method == 'diag') then
+  print*,'Diagonal hessian'
+  allocate(H(tmp_n,1),W(tmp_n,1))
+  tmp_n2 = 1
+elseif (optimization_method == 'none') then
+  print*,'No hessian'
+  allocate(H(tmp_n,1),W(tmp_n,1))
+  tmp_n2 = 1
+else
+  print*,'Unknown optimization_method, please select full, diag or none'
+  call abort
+endif
+print*, 'Absolute value of the hessian:', absolute_eig
+
+! Algorithm
+
+! Here is the main algorithm of the optimization:
+! - First of all we initialize some parameters and we compute the
+!   criterion (the ci energy) before doing any MO rotations
+! - We compute the gradient and the hessian for the active MOs
+! - We diagonalize the hessian
+! - We compute a step and loop to reduce the radius of the
+!   trust region (and the size of the step by the way) until the step is
+!   accepted 
+! - We repeat the process until the convergence 
+!   NB: the convergence criterion can be changed
+
+
+! Loop until the convergence of the optimization
+! call diagonalize_ci 
+
+!### Initialization ###
+nb_iter = 0
+rho = 0.5d0
+not_converged = .True.
+tmp_list = list_act ! Optimization of the active MOs
+nb_cancel_tot = 0
+
+! Renormalization of the weights of the states
+call state_weight_normalization
+
+! Compute the criterion before the loop
+call state_average_energy(prev_criterion)
+
+do while (not_converged)
+  print*,''
+  print*,'******************'
+  print*,'Iteration', nb_iter
+  print*,'******************'
+  print*,''
+
+  ! Gradient
+  call gradient_list_opt(tmp_n, m, tmp_list, v_grad, max_elem_grad, norm_grad)
+
+  ! Hessian
+  if (optimization_method == 'full') then
+    ! Full hessian
+    call hessian_list_opt(tmp_n, m, tmp_list, H, h_f)
+
+    ! Diagonalization of the hessian 
+    call diagonalization_hessian(tmp_n, H, e_val, w)
+
+  elseif (optimization_method == 'diag') then
+    ! Diagonal hessian 
+    call diag_hessian_list_opt(tmp_n, m, tmp_list, H)
+  else
+    ! Identity matrix 
+    do tmp_i = 1, tmp_n
+      H(tmp_i,1) = 1d0
+    enddo
+  endif
+
+  if (optimization_method /= 'full') then
+    ! Sort
+    do tmp_i = 1, tmp_n
+      key(tmp_i) = tmp_i
+      e_val(tmp_i) = H(tmp_i,1)
+    enddo
+    call dsort(e_val,key,tmp_n)
+
+    ! Eigenvalues and eigenvectors
+    do tmp_i = 1, tmp_n
+      w(tmp_i,1) = dble(key(tmp_i))
+    enddo
+
+  endif
+
+  ! Init before the internal loop
+  cancel_step = .True. ! To enter in the loop just after 
+  nb_cancel = 0
+  nb_sub_iter = 0
+
+  ! Loop to reduce the trust radius until the criterion decreases and rho >= thresh_rho
+  do while (cancel_step)
+    print*,''
+    print*,'-----------------------------'
+    print*,'Iteration:    ', nb_iter
+    print*,'Sub iteration:', nb_sub_iter
+    print*,'Max elem grad:', max_elem_grad
+    print*,'-----------------------------'
+
+    ! Hessian,gradient,Criterion -> x 
+    call trust_region_step_w_expected_e(tmp_n,tmp_n2,H,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,tmp_x,must_exit) 
+
+    if (must_exit) then
+      print*,'step_in_trust_region sends: Exit'
+      exit
+    endif
+
+    ! 1D tmp -> 2D tmp 
+    call vec_to_mat_v2(tmp_n, m, tmp_x, tmp_m_x)
+
+    ! Rotation matrix for the active MOs
+    call rotation_matrix(tmp_m_x, m, tmp_R, m, m, info, enforce_step_cancellation)
+
+    ! Security to ensure an unitary transformation
+    if (enforce_step_cancellation) then
+      print*, '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(m, tmp_list, tmp_R, R)
+
+    ! MO rotations
+    call apply_mo_rotation(R, prev_mos)   
+
+    ! Update of the energy before the diagonalization of the hamiltonian
+    call clear_mo_map
+    TOUCH mo_coef psi_det psi_coef ci_energy two_e_dm_mo 
+    call state_average_energy(criterion)
+
+    ! 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 necessary
+    if (cancel_step) then
+      mo_coef = prev_mos
+      call save_mos()
+      nb_cancel = nb_cancel + 1
+      nb_cancel_tot = nb_cancel_tot + 1
+    else
+      ! Diagonalization of the hamiltonian
+      FREE ci_energy! To enforce the recomputation
+      call diagonalize_ci
+      call save_wavefunction_unsorted
+
+      ! Energy obtained after the diagonalization of the CI matrix
+      call state_average_energy(prev_criterion)
+    endif
+
+    nb_sub_iter = nb_sub_iter + 1
+  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
+
+  ! External loop exit conditions
+  if (DABS(max_elem_grad) < thresh_opt_max_elem_grad) then
+    print*,'Converged: DABS(max_elem_grad) < thresh_opt_max_elem_grad'
+    not_converged = .False.
+  endif
+  if (nb_iter >= optimization_max_nb_iter) then
+    print*,'Not converged: nb_iter >= optimization_max_nb_iter'
+    not_converged = .False. 
+  endif
+
+  if (.not. not_converged) then
+    print*,'#############################'
+    print*,'   End of the optimization'
+    print*,'#############################'
+  endif
+enddo
+
+! Deallocation, end
+
+
+deallocate(v_grad,H,R,W,e_val)
+  deallocate(prev_mos,new_mos)
+  if (optimization_method == 'full') then
+    deallocate(h_f)
+  endif
+
+end

src/mo_optimization/save_energy.irp.f

@@ -0,0 +1,9 @@
+subroutine save_energy(E,pt2)
+  implicit none
+  BEGIN_DOC
+! Saves the energy in |EZFIO|.
+  END_DOC
+  double precision, intent(in) :: E(N_states), pt2(N_states)
+  call ezfio_set_fci_energy(E(1:N_states))
+  call ezfio_set_fci_energy_pt2(E(1:N_states)+pt2(1:N_states))
+end

src/mo_optimization/state_average_energy.irp.f

@@ -0,0 +1,72 @@
+! State average energy
+
+! Calculation of the state average energy from the integrals and the
+! density matrices.
+
+! \begin{align*}
+! E = \sum_{ij} h_{ij} \gamma_{ij} + \frac{1}{2} v_{ij}^{kl} \Gamma_{ij}^{kl}
+! \end{align*}
+! $h_{ij}$: mono-electronic integral
+! $\gamma_{ij}$: one electron density matrix
+! $v_{ij}^{kl}$: bi-electronic integral
+! $\Gamma_{ij}^{kl}$: two electrons density matrix
+
+! TODO: OMP version
+
+! PROVIDED:
+! | mo_one_e_integrals | double precision | mono-electronic integrals    |
+! | get_two_e_integral | double precision | bi-electronic integrals      |
+! | one_e_dm_mo        | double precision | one electron density matrix  |
+! | two_e_dm_mo        | double precision | two electrons density matrix |
+! | nuclear_repulsion  | double precision | nuclear repulsion            |
+! | mo_num             | integer          | number of MOs                |
+
+! Output:
+! | energy | double precision | state average energy |
+
+! Internal:
+! | mono_e  | double precision | mono-electronic energy       |
+! | bi_e    | double precision | bi-electronic energy         |
+! | i,j,k,l | integer          | indexes to loop over the MOs |
+
+
+subroutine state_average_energy(energy)
+
+  implicit none
+
+  double precision, intent(out) :: energy
+
+  double precision :: get_two_e_integral
+  double precision :: mono_e, bi_e
+  integer :: i,j,k,l
+  
+  ! mono electronic part
+  mono_e = 0d0
+  do j = 1, mo_num
+    do i = 1, mo_num
+      mono_e =  mono_e + mo_one_e_integrals(i,j) * one_e_dm_mo(i,j)
+    enddo
+  enddo
+
+  ! bi electronic part
+  bi_e = 0d0
+  do l = 1, mo_num
+    do k = 1, mo_num
+      do j = 1, mo_num
+        do i = 1, mo_num
+          bi_e = bi_e + get_two_e_integral(i,j,k,l,mo_integrals_map) * two_e_dm_mo(i,j,k,l)
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ! State average energy
+  energy = mono_e + 0.5d0 * bi_e + nuclear_repulsion
+
+  ! Check
+  !call print_energy_components
+  
+  print*,'State average energy:', energy
+  !print*,ci_energy
+
+end

src/mo_optimization/state_weight_normalization.irp.f

@@ -0,0 +1,29 @@
+subroutine state_weight_normalization
+
+  implicit none
+
+  BEGIN_DOC
+  ! Renormalization of the state weights or enforcing state average
+  ! weights for orbital optimization
+  END_DOC
+
+  integer          :: i
+  double precision :: normalization_factor
+  
+  ! To normalize the sum of the state weights
+  normalization_factor = 0d0
+  do i = 1, N_states
+    normalization_factor = normalization_factor + state_average_weight(i)
+  enddo
+  normalization_factor = 1d0 / normalization_factor
+  
+  do i = 1, N_states
+    state_average_weight(i) = state_average_weight(i) * normalization_factor
+  enddo
+  TOUCH state_average_weight
+  
+  print*, 'Number of states:', N_states
+  print*, 'State average weights:'
+  print*, state_average_weight(:)
+
+end

src/mo_optimization/update_parameters.irp.f

@@ -0,0 +1,15 @@
+! Subroutine toupdate the parameters.
+! Ex: TOUCH mo_coef ...
+
+
+subroutine update_parameters()
+
+  implicit none
+
+  !### TODO
+  ! Touch yours parameters
+  call clear_mo_map
+  TOUCH mo_coef psi_det psi_coef
+  call diagonalize_ci
+  call save_wavefunction_unsorted
+end

src/mo_optimization/update_st_av_ci_energy.irp.f

@@ -0,0 +1,25 @@
+! Update the CI state average energy
+
+! Computes the state average energy 
+! \begin{align*}
+! E =\sum_{i=1}^{N_{states}} E_i . w_i
+! \end{align*}
+
+! $E_i$: energy of state i
+! $w_i$: weight of state i
+
+
+subroutine update_st_av_ci_energy(energy)
+
+  implicit none
+  
+  double precision, intent(out) :: energy
+  integer :: i
+
+  energy = 0d0
+  do i = 1, N_states
+    energy = energy + ci_energy(i) * state_average_weight(i)
+  enddo
+  
+  print*, 'ci_energy :', energy
+end