Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable

bin/qp_test

@@ -46,7 +46,7 @@ def main(arguments):
                 append_bats(dirname, filenames)
     else:
         for (dirname, _, filenames) in os.walk(os.getcwd(), followlinks=False):
-            if "IRPF90_temp" not in dirname:
+            if "IRPF90_temp" not in dirname and "external" not in dirname:
                 append_bats(dirname, filenames)
     l_bats = [y for _, y in sorted(l_bats)]
 
@@ -67,6 +67,7 @@ def main(arguments):
             os.system(test+" python3 bats_to_sh.py "+bats_file+
                                     "| bash")
         else:
+#           print(" ".join(["bats", "--verbose-run", "--trace", bats_file]))
             subprocess.check_call(["bats", "--verbose-run", "--trace", bats_file], env=os.environ)
 
 

etc/qp.rc

@@ -110,6 +110,11 @@ function qp()
       unset COMMAND
       ;;
 
+    "test")
+      shift
+      qp_test $@
+      ;;
+
     *)
       which "qp_$1" &> /dev/null
       if [[ $? -eq 0 ]] ; then
@@ -183,7 +188,7 @@ _qp_Complete()
             ;;
         esac;;
       set_file)
-        COMPREPLY=( $(compgen -W "$(for i in * ; do [[ -f ${i}/ezfio/.version ]] && echo $i ; done)" -- ${cur} ) )
+        COMPREPLY=( $(compgen -W "$(for i in $(find . -name ezfio | sed 's/ezfio$/.version/') ; do [[ -f $i ]] && echo ${i%/.version} ; done)" -- ${cur} ) )
         return 0
         ;;
       plugins)
@@ -215,10 +220,15 @@ _qp_Complete()
             return 0
             ;;
         esac;;
+      test)
+        COMPREPLY=( $(compgen -W "-v -a " -- $cur ) )
+        return 0
+        ;;
       *)
         COMPREPLY=( $(compgen -W 'plugins set_file \
                                    unset_file man \
                                    create_ezfio \
+                                   test \
                                    convert_output_to_ezfio \
                                    -h update' -- $cur ) )
 

src/mo_localization/84.mo_localization.bats

@@ -0,0 +1,97 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
+
+zero () {
+    if [ -z "$1" ]; then echo 0.0; else echo $1; fi
+}
+
+function run() {
+  thresh1=1e-10
+  thresh2=1e-12
+  thresh3=1e-4
+  test_exe scf || skip
+  qp set_file $1
+  qp edit --check
+  qp reset -d
+  qp set_frozen_core
+  qp set localization localization_method boys
+  file="$(echo $1 | sed 's/.ezfio//g')"
+  energy="$(cat $1/hartree_fock/energy)"
+  fb_err1="$(qp run debug_gradient_loc | grep 'Max error' | tail -n 1 | awk '{print $3}')"
+  fb_err2="$(qp run debug_hessian_loc | grep 'Max error' | tail -n 1 | awk '{print $3}')"
+  qp run localization > $file.loc.out
+  fb_energy="$(qp run print_energy | grep -A 1 'Nuclear repulsion energy' | tail -n 1 )"
+  fb_c="$(cat $file.loc.out | grep 'Criterion:Core' | tail -n 1 | awk '{print $3}')i"
+  fb_i="$(cat $file.loc.out | grep 'Criterion:Inactive' | tail -n 1 | awk '{print $3}')"
+  fb_a="$(cat $file.loc.out | grep 'Criterion:Active' | tail -n 1 | awk '{print $3}')"
+  fb_v="$(cat $file.loc.out | grep 'Criterion:Virtual' | tail -n 1 | awk '{print $3}')"
+  qp reset -a
+  qp run scf
+  qp set_frozen_core
+  qp set localization localization_method pipek 
+  pm_err1="$(qp run debug_gradient_loc | grep 'Max error' | tail -n 1 | awk '{print $3}')"
+  pm_err2="$(qp run debug_hessian_loc | grep 'Max error' | tail -n 1 | awk '{print $3}')"
+  qp run localization > $file.loc.out
+  pm_c="$(cat $file.loc.out | grep 'Criterion:Core' | tail -n 1 | awk '{print $3}')i"
+  pm_i="$(cat $file.loc.out | grep 'Criterion:Inactive' | tail -n 1 | awk '{print $3}')"
+  pm_a="$(cat $file.loc.out | grep 'Criterion:Active' | tail -n 1 | awk '{print $3}')"
+  pm_v="$(cat $file.loc.out | grep 'Criterion:Virtual' | tail -n 1 | awk '{print $3}')"
+  pm_energy="$(qp run print_energy | grep -A 1 'Nuclear repulsion energy' | tail -n 1 )"
+  qp set localization localization_method boys
+  qp reset -a
+  qp run scf
+  qp set_frozen_core
+  eq $energy $fb_energy $thresh1 
+  eq $fb_err1 0.0 $thresh2
+  eq $fb_err2 0.0 $thresh2
+  eq $energy $pm_energy $thresh1
+  eq $pm_err1 0.0 $thresh2
+  eq $pm_err2 0.0 $thresh2
+  fb_c=$(zero $fb_c)
+  fb_i=$(zero $fb_i)
+  fb_a=$(zero $fb_a)
+  fb_v=$(zero $fb_v)
+  pm_c=$(zero $pm_c)
+  pm_i=$(zero $pm_i)
+  pm_a=$(zero $pm_a)
+  pm_v=$(zero $pm_v)
+  eq $fb_c $2 $thresh3
+  eq $fb_i $3 $thresh3
+  eq $fb_a $4 $thresh3
+  eq $fb_v $5 $thresh3
+  eq $pm_c $6 $thresh3
+  eq $pm_i $7 $thresh3
+  eq $pm_a $8 $thresh3
+  eq $pm_v $9 $thresh3
+}
+
+@test "b2_stretched" {
+run  b2_stretched.ezfio -32.1357551678876 -47.0041982094667 0.0 -223.470015856259 -1.99990778964451 -2.51376723927071 0.0 -12.8490602539275
+}
+
+@test "clo" {
+run  clo.ezfio -44.1624001765291 -32.4386660941387 0.0 -103.666309287187 -5.99985418946811 -5.46871580225222 0.0 -20.2480064922275
+}
+
+@test "clf" {
+run  clf.ezfio -47.5143398826967 -35.7206886315104 0.0 -107.043029033468 -5.99994222062230 -6.63916513458470 0.0 -19.7035159913484
+}
+
+@test "h2o2" {
+run  h2o2.ezfio -7.76848143170524 -30.9694344369829 0.0 -175.898343829453 -1.99990497554575 -5.62980322957485 0.0 -33.5699813186666
+}
+
+@test "h2o" {
+run  h2o.ezfio 0.0 -2.52317434969591 0.0 -45.3136377925359 0.0 -3.01248365356981 0.0 -22.4470831240924
+}
+
+@test "h3coh" {
+run  h3coh.ezfio -3.66763692804590 -24.0463089480870 0.0 -111.485948435075 -1.99714061342078 -4.89242181322988 0.0 -23.6405412057679
+}
+
+@test "n2h4" {
+run  n2h4.ezfio -7.46608163002070 -35.7632174051822 0.0 -305.913449004632 -1.99989326143356 -4.62496615892268 0.0 -51.5171904685553
+}
+

src/mo_localization/EZFIO.cfg

@@ -0,0 +1,54 @@
+[localization_method]
+type: character*(32)
+doc: Method for the orbital localization. boys: Foster-Boys, pipek: Pipek-Mezey.
+interface: ezfio,provider,ocaml
+default: boys
+
+[localization_max_nb_iter]
+type: integer
+doc: Maximal number of iterations for the orbital localization.
+interface: ezfio,provider,ocaml
+default: 1000
+
+[localization_use_hessian]
+type: logical
+doc: If true, it uses the trust region algorithm with the gradient and the diagonal of the hessian. Else it computes the rotation between each pair of MOs that should be applied to maximize/minimize the localization criterion. The last option is not easy to converge. 
+interface: ezfio,provider,ocaml
+default: true
+
+[auto_mo_class]
+type: logical
+doc: If true, set automatically the classes.
+interface: ezfio,provider,ocaml
+default: true
+
+[thresh_loc_max_elem_grad]
+type: double precision
+doc: Threshold for the convergence, the localization exits when the largest element in the gradient is smaller than thresh_localization_max_elem_grad.
+interface: ezfio,provider,ocaml
+default: 1.e-6
+
+[kick_in_mos]
+type: logical
+doc: If True, it applies a rotation of an angle angle_pre_rot between the MOs of a same mo_class before the localization.
+interface: ezfio,provider,ocaml
+default: true
+
+[angle_pre_rot]
+type: double precision
+doc: To define the angle for the rotation of the MOs before the localization (in rad).
+interface: ezfio,provider,ocaml
+default: 0.1
+
+[sort_mos_by_e]
+type: logical
+doc: If True, the MOs are sorted using the diagonal elements of the Fock matrix.
+interface: ezfio,provider,ocaml
+default: false
+
+[debug_hf]
+type: logical
+doc: If True, prints the HF energy before/after the different steps of the localization. Only for debugging.
+interface: ezfio,provider,ocaml
+default: false
+

src/mo_localization/NEED

@@ -0,0 +1,3 @@
+hartree_fock
+utils_trust_region
+determinants

src/mo_localization/README.md

@@ -0,0 +1,113 @@
+# Orbital localisation
+To localize the MOs:  
+```
+qp run localization  
+```
+By default, the different otbital classes are automatically set by splitting  
+the orbitales in the following classes:  
+- Core -> Core  
+- Active, doubly occupied -> Inactive  
+- Active, singly occupied -> Active  
+- Active, empty -> Virtual  
+- Deleted -> Deleted  
+The orbitals will be localized among each class, excpect the deleted ones.
+If you want to choose another splitting, you can set
+```
+qp set mo_localization auto_mo_class false
+```
+and define the classes with
+```
+qp set_mo_class -c [] -a [] -v [] -i [] -d [] 
+```
+for more information
+```
+qp set_mo_class -q
+```
+We don't care about the name of the   
+mo classes. The algorithm just localizes all the MOs of  
+a given class between them, for all the classes, except the deleted MOs.  
+If you are using the last option don't forget to reset the initial mo classes  
+after the localization.
+
+Before the localization, a kick is done for each mo class  
+(except the deleted ones) to break the MOs. This is done by   
+doing a given rotation between the MOs.
+This feature can be removed by setting:
+```
+qp set localization kick_in_mos false
+```
+and the default angle for the rotation can be changed with:
+```
+qp set localization angle_pre_rot 1e-3 # or something else
+```
+
+After the localization, the MOs of each class (except the deleted ones)  
+can be sorted between them using the diagonal elements of  
+the fock matrix with:
+```
+qp set localization sort_mos_by_e true
+```
+
+You can check the Hartree-Fock energy before/during/after the localization  
+by putting (only for debugging):
+```
+qp set localization debug_hf true 
+```
+
+## Foster-Boys & Pipek-Mezey
+Foster-Boys:  
+``` 
+qp set localization localization_method boys 
+``` 
+ 
+Pipek-Mezey:  
+``` 
+qp set localization localization_method pipek 
+``` 
+
+# Break the spatial symmetry of the MOs
+This program work exactly as the localization.  
+To break the spatial symmetry of the MOs:   
+```
+qp run break_spatial_sym
+```
+The default angle for the rotations is too big for this kind of
+application, a value between 1e-3 and 1e-6 should break the spatial
+symmetry with just a small change in the energy:
+```
+qp set localization angle_pre_rot 1e-3
+``` 
+
+# With or without hessian + trust region
+With hessian +  trust region
+```
+qp set localization localisation_use_hessian true
+```
+It uses the trust region algorithm with the diagonal of the hessian of the 
+localization criterion with respect to the MO rotations.   
+
+Without the hessian and the trust region
+```
+qp set localization localisation_use_hessian false
+```
+By doing so it does not require to store the hessian but the
+convergence is not easy, in particular for virtual MOs.
+It seems that it not possible to converge with Pipek-Mezey
+localization with this approach.
+
+# Parameters
+Some other parameters are available for the localization (qp edit for more details).
+
+# 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_localization/break_spatial_sym.irp.f

@@ -0,0 +1,27 @@
+! ! A small program to break the spatial symmetry of the MOs.
+
+! ! You have to defined your MO classes or set security_mo_class to false
+! ! with:
+! ! qp set orbital_optimization security_mo_class false
+
+! ! The default angle for the rotations is too big for this kind of
+! ! application, a value between 1e-3 and 1e-6 should break the spatial
+! ! symmetry with just a small change in the energy. 
+
+
+program break_spatial_sym
+
+  !BEGIN_DOC
+  ! Break the symmetry of the MOs with a rotation
+  !END_DOC
+
+  implicit none
+
+  kick_in_mos = .True.
+  TOUCH kick_in_mos
+
+  call set_classes_loc 
+  call apply_pre_rotation
+  call unset_classes_loc 
+  
+end

src/mo_localization/debug_gradient_loc.irp.f

@@ -0,0 +1,65 @@
+program debug_gradient_loc
+
+  !BEGIN_DOC
+  ! Check if the gradient is correct
+  !END_DOC
+
+  implicit none
+
+  integer :: list_size, n
+  integer, allocatable :: list(:)
+  double precision, allocatable :: v_grad(:), v_grad2(:)
+  double precision :: norm, max_elem, threshold, max_error
+  integer :: i, nb_error
+
+  threshold = 1d-12
+
+  list_size = dim_list_act_orb
+
+  allocate(list(list_size))
+  
+  list = list_act
+
+  n = list_size*(list_size-1)/2
+  
+  allocate(v_grad(n),v_grad2(n))
+
+  if (localization_method == 'boys') then
+    print*,'Foster-Boys'
+    call gradient_FB(n,list_size,list,v_grad,max_elem,norm)
+    call gradient_FB_omp(n,list_size,list,v_grad2,max_elem,norm)
+  elseif (localization_method == 'pipek') then
+    print*,'Pipek-Mezey'
+    call gradient_PM(n,list_size,list,v_grad,max_elem,norm)
+    call gradient_PM(n,list_size,list,v_grad2,max_elem,norm) 
+  else
+    print*,'Unknown localization_method, please select boys or pipek'
+    call abort
+  endif
+ 
+  do i = 1, n
+    print*,i,v_grad(i)
+  enddo
+
+  v_grad = v_grad - v_grad2
+
+  nb_error = 0
+  max_elem = 0d0
+
+  do i = 1, n
+    if (dabs(v_grad(i)) > threshold) then
+      print*,v_grad(i)
+      nb_error = nb_error + 1
+      if (dabs(v_grad(i)) > max_elem) then
+        max_elem = v_grad(i)
+      endif
+    endif
+  enddo
+
+  print*,'Threshold error', threshold
+  print*, 'Nb error', nb_error
+  print*,'Max error', max_elem
+
+  deallocate(v_grad,v_grad2)
+ 
+end

src/mo_localization/debug_hessian_loc.irp.f

@@ -0,0 +1,65 @@
+program debug_hessian_loc
+
+  !BEGIN_DOC
+  ! Check if the hessian is correct
+  !END_DOC
+  
+  implicit none
+
+  integer :: list_size, n
+  integer, allocatable :: list(:)
+  double precision, allocatable :: H(:), H2(:)
+  double precision :: threshold, max_error, max_elem
+  integer :: i, nb_error
+
+  threshold = 1d-12
+
+  list_size = dim_list_act_orb
+
+  allocate(list(list_size))
+  
+  list = list_act
+
+  n = list_size*(list_size-1)/2
+  
+  allocate(H(n),H2(n))
+
+  if (localization_method == 'boys') then
+    print*,'Foster-Boys'
+    call hessian_FB(n,list_size,list,H)
+    call hessian_FB_omp(n,list_size,list,H2)
+  elseif(localization_method == 'pipek') then
+    print*,'Pipek-Mezey'
+    call hessian_PM(n,list_size,list,H)
+    call hessian_PM(n,list_size,list,H2)
+  else
+    print*,'Unknown localization_method, please select boys or pipek'
+    call abort
+  endif
+ 
+  do i = 1, n
+    print*,i,H(i)
+  enddo
+
+  H = H - H2
+
+  nb_error = 0
+  max_elem = 0d0
+
+  do i = 1, n
+    if (dabs(H(i)) > threshold) then
+      print*,H(i)
+      nb_error = nb_error + 1
+      if (dabs(H(i)) > max_elem) then
+        max_elem = H(i)
+      endif
+    endif
+  enddo
+
+  print*,'Threshold error', threshold
+  print*, 'Nb error', nb_error
+  print*,'Max error', max_elem
+
+  deallocate(H,H2)
+  
+end

src/mo_localization/kick_the_mos.irp.f

@@ -0,0 +1,16 @@
+program kick_the_mos
+
+  !BEGIN_DOC
+  ! To do a small rotation of the MOs
+  !END_DOC
+
+  implicit none
+
+  kick_in_mos = .True.
+  TOUCH kick_in_mos
+
+  call set_classes_loc
+  call apply_pre_rotation
+  call unset_classes_loc
+  
+end

src/mo_localization/localization.irp.f

@@ -0,0 +1,520 @@
+program localization
+
+  implicit none
+
+  call set_classes_loc
+  call run_localization
+  call unset_classes_loc
+  
+end
+
+
+
+
+! Variables:
+! | pre_rot(mo_num, mo_num)   | double precision | Matrix for the pre rotation                       |
+! | R(mo_num,mo_num)          | double precision | Rotation matrix                                   |
+! | tmp_R(:,:)                | double precision | Rottation matrix in a subsapce                    |
+! | prev_mos(ao_num, mo_num)  | double precision | Previous mo_coef                                  |
+! | spatial_extent(mo_num)    | double precision | Spatial extent of the orbitals                    |
+! | criterion                 | double precision | Localization criterion                            |
+! | prev_criterion            | double precision | Previous criterion                                |
+! | criterion_model           | double precision | Estimated next criterion                          |
+! | rho                       | double precision | Ratio to measure the agreement between the model  |
+! |                           |                  | and the reality                                   |
+! | delta                     | double precision | Radisu of the trust region                        |
+! | norm_grad                 | double precision | Norm of the gradient                              |
+! | info                      | integer          | for dsyev from Lapack                             |
+! | max_elem                  | double precision | maximal element in the gradient                   |
+! | v_grad(:)                 | double precision | Gradient                                          |
+! | H(:,:)                    | double precision | Hessian (diagonal)                                |
+! | e_val(:)                  | double precision | Eigenvalues of the hessian                        |
+! | W(:,:)                    | double precision | Eigenvectors of the hessian                       |
+! | tmp_x(:)                  | double precision | Step in 1D (in a subaspace)                       |
+! | tmp_m_x(:,:)              | double precision | Step in 2D (in a subaspace)                       |
+! | tmp_list(:)               | double precision | List of MOs in a mo_class                         |
+! | i,j,k                     | integer          | Indexes in the full MO space                      |
+! | tmp_i, tmp_j, tmp_k       | integer          | Indexes in a subspace                             |
+! | l                         | integer          | Index for the mo_class                            |
+! | key(:)                    | integer          | Key to sort the eigenvalues of the hessian        |
+! | nb_iter                   | integer          | Number of iterations                              |
+! | must_exit                 | logical          | To exit the trust region loop                     |
+! | cancel_step               | logical          | To cancel a step                                  |
+! | not_*converged            | logical          | To localize the different mo classes              |
+! | t*                        | double precision | To measure the time                               |
+! | n                         | integer          | mo_num*(mo_num-1)/2, number of orbital parameters |
+! | tmp_n                     | integer          | dim_subspace*(dim_subspace-1)/2                   |
+! |                           |                  | Number of dimension in the subspace               |
+
+! Variables in qp_edit for the localization:
+! | localization_method      |
+! | localization_max_nb_iter |
+! | default_mo_class         |
+! | thresh_loc_max_elem_grad |
+! | kick_in_mos              |
+! | angle_pre_rot            |
+
+! + all the variables for the trust region
+
+! Cf. qp_edit orbital optimization
+
+
+subroutine run_localization
+
+  include 'pi.h'
+
+  BEGIN_DOC
+  ! Orbital localization
+  END_DOC
+
+  implicit none
+
+  ! Variables
+  double precision, allocatable :: pre_rot(:,:), R(:,:)
+  double precision, allocatable :: prev_mos(:,:), spatial_extent(:), tmp_R(:,:)
+  double precision              :: criterion, norm_grad
+  integer                       :: i,j,k,l,p, tmp_i, tmp_j, tmp_k
+  integer                       :: info
+  integer                       :: n, tmp_n, tmp_list_size
+  double precision, allocatable :: v_grad(:), H(:), tmp_m_x(:,:), tmp_x(:),W(:),e_val(:)
+  double precision              :: max_elem, t1, t2, t3, t4, t5, t6
+  integer, allocatable          :: tmp_list(:), key(:)
+  double precision              :: prev_criterion, rho, delta, criterion_model
+  integer                       :: nb_iter, nb_sub_iter
+  logical                       :: not_converged, not_core_converged
+  logical                       :: not_act_converged, not_inact_converged, not_virt_converged
+  logical                       :: use_trust_region, must_exit, cancel_step,enforce_step_cancellation
+
+  n = mo_num*(mo_num-1)/2
+
+  ! Allocation
+  allocate(spatial_extent(mo_num))
+  allocate(pre_rot(mo_num, mo_num), R(mo_num, mo_num))
+  allocate(prev_mos(ao_num, mo_num))
+
+  ! Locality before the localization
+  call compute_spatial_extent(spatial_extent)
+
+  ! Choice of the method
+  print*,''
+  print*,'Localization method:',localization_method
+  if (localization_method == 'boys') then
+    print*,'Foster-Boys localization'
+  elseif (localization_method == 'pipek') then
+    print*,'Pipek-Mezey localization'
+  else
+    print*,'Unknown localization_method, please select boys or pipek'
+    call abort
+  endif
+  print*,''
+
+  ! Localization criterion (FB, PM, ...) for each mo_class
+  print*,'### Before the pre rotation'
+  
+  ! Debug
+  if (debug_hf) then
+    print*,'HF energy:', HF_energy
+  endif
+  
+  do l = 1, 4
+    if (l==1) then ! core
+      tmp_list_size = dim_list_core_orb
+    elseif (l==2) then ! act
+      tmp_list_size = dim_list_act_orb
+    elseif (l==3) then ! inact
+      tmp_list_size = dim_list_inact_orb
+    else ! virt
+      tmp_list_size = dim_list_virt_orb
+    endif
+
+     ! Allocation tmp array  
+    allocate(tmp_list(tmp_list_size))
+
+    ! To give the list of MOs in a mo_class
+    if (l==1) then ! core
+      tmp_list = list_core
+    elseif (l==2) then
+      tmp_list = list_act
+    elseif (l==3) then
+      tmp_list = list_inact
+    else
+      tmp_list = list_virt
+    endif
+
+    if (tmp_list_size >= 2) then
+      call criterion_localization(tmp_list_size, tmp_list,criterion)
+      print*,'Criterion:', criterion, mo_class(tmp_list(1))
+    endif
+
+    deallocate(tmp_list)    
+
+  enddo
+
+  ! Debug
+  !print*,'HF', HF_energy
+
+! Loc
+
+! Pre rotation, to give a little kick in the MOs
+  call apply_pre_rotation()
+
+  ! Criterion after the pre rotation
+  ! Localization criterion (FB, PM, ...) for each mo_class
+  print*,'### After the pre rotation'
+  
+  ! Debug
+  if (debug_hf) then
+    touch mo_coef
+    print*,'HF energy:', HF_energy
+  endif
+  
+  do l = 1, 4
+    if (l==1) then ! core
+      tmp_list_size = dim_list_core_orb
+    elseif (l==2) then ! act
+      tmp_list_size = dim_list_act_orb
+    elseif (l==3) then ! inact
+      tmp_list_size = dim_list_inact_orb
+    else ! virt
+      tmp_list_size = dim_list_virt_orb
+    endif
+    
+    if (tmp_list_size >= 2) then
+      ! Allocation tmp array  
+      allocate(tmp_list(tmp_list_size))
+
+      ! To give the list of MOs in a mo_class
+      if (l==1) then ! core
+        tmp_list = list_core
+      elseif (l==2) then
+        tmp_list = list_act
+      elseif (l==3) then
+        tmp_list = list_inact
+      else
+        tmp_list = list_virt
+      endif
+
+      call criterion_localization(tmp_list_size, tmp_list,criterion)
+      print*,'Criterion:', criterion, trim(mo_class(tmp_list(1)))
+      
+      deallocate(tmp_list)    
+    endif
+    
+  enddo
+
+  ! Debug
+  !print*,'HF', HF_energy
+
+  print*,''
+  print*,'========================'
+  print*,'  Orbital localization'
+  print*,'========================'
+  print*,'' 
+
+  !Initialization
+  not_converged = .TRUE.
+
+  ! To do the localization only if there is at least 2 MOs
+  if (dim_list_core_orb >= 2) then
+    not_core_converged = .TRUE.
+  else
+    not_core_converged = .FALSE. 
+  endif
+  
+  if (dim_list_act_orb >= 2) then
+    not_act_converged = .TRUE.
+  else
+    not_act_converged = .FALSE.
+  endif
+
+  if (dim_list_inact_orb >= 2) then
+    not_inact_converged = .TRUE.
+  else
+    not_inact_converged = .FALSE.
+  endif
+
+  if (dim_list_virt_orb >= 2) then
+    not_virt_converged = .TRUE.
+  else
+    not_virt_converged = .FALSE.
+  endif
+ 
+  ! Loop over the mo_classes
+  do l = 1, 4
+
+    if (l==1) then ! core
+      not_converged = not_core_converged
+      tmp_list_size = dim_list_core_orb
+    elseif (l==2) then ! act
+      not_converged = not_act_converged
+      tmp_list_size = dim_list_act_orb
+    elseif (l==3) then ! inact
+      not_converged = not_inact_converged
+      tmp_list_size = dim_list_inact_orb
+    else ! virt
+      not_converged = not_virt_converged
+      tmp_list_size = dim_list_virt_orb
+    endif
+
+    ! Next iteration if converged = true 
+    if (.not. not_converged) then
+      cycle
+    endif
+ 
+    ! Allocation tmp array  
+    allocate(tmp_list(tmp_list_size))
+
+    ! To give the list of MOs in a mo_class
+    if (l==1) then ! core
+      tmp_list = list_core
+    elseif (l==2) then
+      tmp_list = list_act
+    elseif (l==3) then
+      tmp_list = list_inact
+    else
+      tmp_list = list_virt
+    endif
+
+    ! Display
+    if (not_converged) then
+      print*,''
+      print*,'###', trim(mo_class(tmp_list(1))), 'MOs ###'
+      print*,''
+    endif
+
+    ! Size for the 2D -> 1D transformation 
+    tmp_n = tmp_list_size * (tmp_list_size - 1)/2
+
+    ! Without hessian + trust region 
+    if (.not. localization_use_hessian) then
+       
+      ! Allocation of temporary arrays
+      allocate(v_grad(tmp_n), tmp_m_x(tmp_list_size, tmp_list_size))
+      allocate(tmp_R(tmp_list_size, tmp_list_size), tmp_x(tmp_n))
+
+      ! Criterion
+      call criterion_localization(tmp_list_size, tmp_list, prev_criterion)
+
+      ! Init
+      nb_iter = 0
+      delta = 1d0
+
+      !Loop
+      do while (not_converged)
+         
+        print*,''
+        print*,'***********************'
+        print*,'Iteration', nb_iter
+        print*,'***********************'
+        print*,''
+        
+        ! Angles of rotation
+        call theta_localization(tmp_list, tmp_list_size, tmp_m_x, max_elem)
+        tmp_m_x = - tmp_m_x * delta
+
+        ! Rotation submatrix
+        call rotation_matrix(tmp_m_x, tmp_list_size, tmp_R, tmp_list_size, tmp_list_size, &
+             info, enforce_step_cancellation)
+
+        ! To ensure that the rotation matrix is unitary
+        if (enforce_step_cancellation) then
+          print*, 'Step cancellation, too large error in the rotation matrix'
+          delta = delta * 0.5d0
+          cycle
+        else
+          delta = min(delta * 2d0, 1d0)
+        endif
+
+        ! Full rotation matrix and application of the rotation
+        call sub_to_full_rotation_matrix(tmp_list_size, tmp_list, tmp_R, R)
+        call apply_mo_rotation(R, prev_mos)
+
+        ! Update the needed data
+        call update_data_localization()
+
+        ! New criterion
+        call criterion_localization(tmp_list_size, tmp_list, criterion)
+        print*,'Criterion:', trim(mo_class(tmp_list(1))), nb_iter, criterion
+        print*,'Max elem :', max_elem
+        print*,'Delta    :', delta
+        
+        nb_iter = nb_iter + 1
+
+        ! Exit
+        if (nb_iter >= localization_max_nb_iter .or. dabs(max_elem) < thresh_loc_max_elem_grad) then
+           not_converged = .False.
+        endif
+      enddo
+
+      ! Save the changes
+      call update_data_localization()
+      call save_mos()
+      TOUCH mo_coef
+
+      ! Deallocate
+      deallocate(v_grad, tmp_m_x, tmp_list)
+      deallocate(tmp_R, tmp_x)
+
+    ! Trust region
+    else
+    
+      ! Allocation of temporary arrays
+      allocate(v_grad(tmp_n), H(tmp_n), tmp_m_x(tmp_list_size, tmp_list_size))
+      allocate(tmp_R(tmp_list_size, tmp_list_size))
+      allocate(tmp_x(tmp_n), W(tmp_n), e_val(tmp_n), key(tmp_n))
+
+      ! ### Initialization ###
+      delta = 0d0 ! can be deleted (normally)
+      nb_iter = 0 ! Must start at 0 !!!
+      rho = 0.5d0 ! Must be 0.5
+
+      ! Compute the criterion before the loop
+      call criterion_localization(tmp_list_size, tmp_list, prev_criterion)
+
+      ! Loop until the convergence
+      do while (not_converged)
+
+        print*,''
+        print*,'***********************'
+        print*,'Iteration', nb_iter
+        print*,'***********************'
+        print*,''
+  
+        ! Gradient
+        call gradient_localization(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad)
+        ! Diagonal hessian
+        call hessian_localization(tmp_n, tmp_list_size, tmp_list, H)
+        
+        ! Diagonalization of the diagonal hessian by hands
+        !call diagonalization_hessian(tmp_n,H,e_val,w)
+        do i = 1, tmp_n
+          e_val(i) = H(i)
+        enddo
+
+        ! Key list for dsort
+        do i = 1, tmp_n 
+          key(i) = i
+        enddo
+
+        ! Sort of the eigenvalues
+        call dsort(e_val, key, tmp_n)
+
+        ! Eigenvectors
+        W = 0d0
+        do i = 1, tmp_n
+          W(i) = dble(key(i))
+        enddo
+
+        ! To enter in the loop just after
+        cancel_step = .True.
+        nb_sub_iter = 0 
+
+        ! Loop to reduce the trust radius until the criterion decreases and rho >= thresh_rho
+        do while (cancel_step)
+          print*,'-----------------------------'
+          print*, mo_class(tmp_list(1))
+          print*,'Iteration:', nb_iter
+          print*,'Sub iteration:', nb_sub_iter
+          print*,'Max elem grad:', max_elem
+          print*,'-----------------------------'
+
+          ! Hessian,gradient,Criterion -> x 
+          call trust_region_step_w_expected_e(tmp_n,1, H, W, e_val, v_grad, prev_criterion, &
+               rho, nb_iter, delta, criterion_model, tmp_x, must_exit)
+
+          ! Internal loop exit condition
+          if (must_exit) then
+            print*,'trust_region_step_w_expected_e sent: Exit'
+            exit 
+          endif
+
+          ! 1D tmp -> 2D tmp 
+          call vec_to_mat_v2(tmp_n, tmp_list_size, tmp_x, tmp_m_x)
+
+          ! Rotation submatrix (square matrix tmp_list_size by tmp_list_size)
+          call rotation_matrix(tmp_m_x, tmp_list_size, tmp_R, tmp_list_size, tmp_list_size, &
+               info, enforce_step_cancellation)
+
+          if (enforce_step_cancellation) then
+            print*, 'Step cancellation, too large error in the rotation matrix'
+            rho = 0d0
+            cycle
+          endif
+
+          ! tmp_R to R, subspace to full space
+          call sub_to_full_rotation_matrix(tmp_list_size, tmp_list, tmp_R, R)
+        
+          ! Rotation of the MOs
+          call apply_mo_rotation(R, prev_mos)   
+
+          ! Update the things related to mo_coef
+          call update_data_localization()
+
+          ! Update the criterion
+          call criterion_localization(tmp_list_size, tmp_list, criterion)
+          print*,'Criterion:', trim(mo_class(tmp_list(1))), nb_iter, 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, previous MOs
+          if (cancel_step) then
+            mo_coef = prev_mos
+          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_exti = .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) < thresh_loc_max_elem_grad) then
+          not_converged = .False.
+        endif
+        if (nb_iter > localization_max_nb_iter) then
+          not_converged = .False.
+        endif
+      enddo
+
+      ! Deallocation of temporary arrays
+      deallocate(v_grad, H, tmp_m_x, tmp_R, tmp_list, tmp_x, W, e_val, key)
+      
+      ! Save the MOs
+      call save_mos()
+      TOUCH mo_coef
+ 
+      ! Debug
+      if (debug_hf) then
+        touch mo_coef
+        print*,'HF energy:', HF_energy
+      endif
+      
+    endif
+  enddo
+
+  ! Seems unecessary
+  TOUCH mo_coef 
+
+  ! To sort the MOs using the diagonal elements of the Fock matrix
+  if (sort_mos_by_e) then
+    call run_sort_by_fock_energies()
+  endif
+
+  ! Debug
+  if (debug_hf) then
+    touch mo_coef
+    print*,'HF energy:', HF_energy
+  endif
+
+  ! Locality after the localization
+  call compute_spatial_extent(spatial_extent)
+
+end

src/mo_localization/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_localization/org/break_spatial_sym.org

@@ -0,0 +1,28 @@
+! A small program to break the spatial symmetry of the MOs.
+
+! You have to defined your MO classes or set security_mo_class to false
+! with:
+! qp set orbital_optimization security_mo_class false
+
+! The default angle for the rotations is too big for this kind of
+! application, a value between 1e-3 and 1e-6 should break the spatial
+! symmetry with just a small change in the energy. 
+
+#+BEGIN_SRC f90 :comments org :tangle break_spatial_sym.irp.f
+program break_spatial_sym
+
+  !BEGIN_DOC
+  ! Break the symmetry of the MOs with a rotation
+  !END_DOC
+
+  implicit none
+
+  kick_in_mos = .True.
+  TOUCH kick_in_mos
+
+  call set_classes_loc 
+  call apply_pre_rotation
+  call unset_classes_loc 
+  
+end
+#+END_SRC

src/mo_localization/org/debug_gradient_loc.org

@@ -0,0 +1,67 @@
+#+BEGIN_SRC f90 :comments org :tangle debug_gradient_loc.irp.f
+program debug_gradient_loc
+
+  !BEGIN_DOC
+  ! Check if the gradient is correct
+  !END_DOC
+
+  implicit none
+
+  integer :: list_size, n
+  integer, allocatable :: list(:)
+  double precision, allocatable :: v_grad(:), v_grad2(:)
+  double precision :: norm, max_elem, threshold, max_error
+  integer :: i, nb_error
+
+  threshold = 1d-12
+
+  list_size = dim_list_act_orb
+
+  allocate(list(list_size))
+  
+  list = list_act
+
+  n = list_size*(list_size-1)/2
+  
+  allocate(v_grad(n),v_grad2(n))
+
+  if (localization_method == 'boys') then
+    print*,'Foster-Boys'
+    call gradient_FB(n,list_size,list,v_grad,max_elem,norm)
+    call gradient_FB_omp(n,list_size,list,v_grad2,max_elem,norm)
+  elseif (localization_method == 'pipek') then
+    print*,'Pipek-Mezey'
+    call gradient_PM(n,list_size,list,v_grad,max_elem,norm)
+    call gradient_PM(n,list_size,list,v_grad2,max_elem,norm) 
+  else
+    print*,'Unknown localization_method, please select boys or pipek'
+    call abort
+  endif
+ 
+  do i = 1, n
+    print*,i,v_grad(i)
+  enddo
+
+  v_grad = v_grad - v_grad2
+
+  nb_error = 0
+  max_elem = 0d0
+
+  do i = 1, n
+    if (dabs(v_grad(i)) > threshold) then
+      print*,v_grad(i)
+      nb_error = nb_error + 1
+      if (dabs(v_grad(i)) > max_elem) then
+        max_elem = v_grad(i)
+      endif
+    endif
+  enddo
+
+  print*,'Threshold error', threshold
+  print*, 'Nb error', nb_error
+  print*,'Max error', max_elem
+
+  deallocate(v_grad,v_grad2)
+ 
+end
+#+END_SRC

src/mo_localization/org/debug_hessian_loc.org

@@ -0,0 +1,67 @@
+#+BEGIN_SRC f90 :comments org :tangle debug_hessian_loc.irp.f
+program debug_hessian_loc
+
+  !BEGIN_DOC
+  ! Check if the hessian is correct
+  !END_DOC
+  
+  implicit none
+
+  integer :: list_size, n
+  integer, allocatable :: list(:)
+  double precision, allocatable :: H(:), H2(:)
+  double precision :: threshold, max_error, max_elem
+  integer :: i, nb_error
+
+  threshold = 1d-12
+
+  list_size = dim_list_act_orb
+
+  allocate(list(list_size))
+  
+  list = list_act
+
+  n = list_size*(list_size-1)/2
+  
+  allocate(H(n),H2(n))
+
+  if (localization_method == 'boys') then
+    print*,'Foster-Boys'
+    call hessian_FB(n,list_size,list,H)
+    call hessian_FB_omp(n,list_size,list,H2)
+  elseif(localization_method == 'pipek') then
+    print*,'Pipek-Mezey'
+    call hessian_PM(n,list_size,list,H)
+    call hessian_PM(n,list_size,list,H2)
+  else
+    print*,'Unknown localization_method, please select boys or pipek'
+    call abort
+  endif
+ 
+  do i = 1, n
+    print*,i,H(i)
+  enddo
+
+  H = H - H2
+
+  nb_error = 0
+  max_elem = 0d0
+
+  do i = 1, n
+    if (dabs(H(i)) > threshold) then
+      print*,H(i)
+      nb_error = nb_error + 1
+      if (dabs(H(i)) > max_elem) then
+        max_elem = H(i)
+      endif
+    endif
+  enddo
+
+  print*,'Threshold error', threshold
+  print*, 'Nb error', nb_error
+  print*,'Max error', max_elem
+
+  deallocate(H,H2)
+  
+end
+#+END_SRC

src/mo_localization/org/kick_the_mos.org

@@ -0,0 +1,18 @@
+#+BEGIN_SRC f90 :comments org :tangle kick_the_mos.irp.f
+program kick_the_mos
+
+  !BEGIN_DOC
+  ! To do a small rotation of the MOs
+  !END_DOC
+
+  implicit none
+
+  kick_in_mos = .True.
+  TOUCH kick_in_mos
+
+  call set_classes_loc
+  call apply_pre_rotation
+  call unset_classes_loc
+  
+end
+#+END_SRC

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,23 @@
+[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_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/optimization.irp.f

@@ -0,0 +1,86 @@
+program optimization
+  
+  read_wf = .true. ! must be True for the orbital optimization !!!
+  TOUCH read_wf 
+  call run_optimization
+
+end
+
+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

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
+