Merge branch 'master' of github.com:LCPQ/quantum_package

Conflicts:
	src/Hartree_Fock/HF_density_matrix_ao.irp.f
	src/Hartree_Fock/mo_SCF_iterations.irp.f
	src/MOs/cholesky_mo.irp.f

src/CID/H_apply.irp.f

@@ -0,0 +1,9 @@
+! Generates subroutine H_apply_cisd
+! ----------------------------------
+
+BEGIN_SHELL [ /usr/bin/env python ]
+from generate_h_apply import H_apply
+H = H_apply("cisd")
+print H
+END_SHELL
+

src/CID/Makefile

@@ -0,0 +1,8 @@
+default: all
+
+# Define here all new external source files and objects.Don't forget to prefix the
+# object files with IRPF90_temp/
+SRC=
+OBJ=
+
+include $(QPACKAGE_ROOT)/src/Makefile.common

src/CID/NEEDED_MODULES

@@ -0,0 +1,2 @@
+AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output SingleRefMethod Utils Selectors_full
+

src/CID/README.rst

@@ -0,0 +1,43 @@
+CISD
+====
+
+This is a test directory which builds a CISD by setting the follwoing rules:
+
+* The only generator determinant is the Hartee-Fock (single-reference method)
+* All generated determinants are included in the wave function (no perturbative
+  selection)
+
+These rules are set in the ``H_apply.irp.f`` file.
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
+* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
+* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
+* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
+* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
+* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd_lapack.irp.f#L1>`_
+  Undocumented
+
+
+

src/CID/cid.irp.f

@@ -0,0 +1,20 @@
+program cisd
+  implicit none
+  integer :: i
+
+  print *,   'HF      = ', HF_energy
+  print *,  'N_states = ', N_states
+  N_det = 1
+  touch psi_det psi_coef N_det
+  call H_apply_cisd
+  print *,  'N_det = ', N_det
+  do i = 1,N_states
+   print *,  'energy  = ',CI_energy(i) 
+   print *,  'E_corr  = ',CI_electronic_energy(i) - ref_bitmask_energy
+  enddo
+
+! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int)
+! do i = 1, N_states
+!  print*,'eigvalues(i) = ',eigvalues(i)
+! enddo
+end

src/CID/cid_lapack.irp.f

@@ -0,0 +1,16 @@
+program cisd
+  implicit none
+  integer :: i
+
+  diag_algorithm = "Lapack"
+  touch diag_algorithm
+  print *,   'HF      = ', HF_energy
+  print *,  'N_states = ', N_states
+  call H_apply_cisd
+  print *,  'N_det = ', N_det
+  do i = 1,N_states
+   print *,  'energy  = ',CI_energy(i) 
+   print *,  'E_corr  = ',CI_electronic_energy(i) - ref_bitmask_energy
+  enddo
+
+end

src/CID/do_mono_double.irp.f

@@ -0,0 +1,19 @@
+BEGIN_PROVIDER [logical, do_double_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the double excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_double_excitations = .True.
+
+END_PROVIDER
+
+BEGIN_PROVIDER [logical, do_mono_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the mono excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_mono_excitations = .False.
+
+END_PROVIDER

src/CID_SC2_selected/H_apply.irp.f

@@ -0,0 +1,10 @@
+use bitmasks
+BEGIN_SHELL [ /usr/bin/env python ]
+from generate_h_apply import *
+from perturbation import perturbations
+
+s = H_apply("PT2",SingleRef=True)
+s.set_perturbation("epstein_nesbet_sc2_projected")
+print s
+END_SHELL
+

src/CID_SC2_selected/Makefile

@@ -0,0 +1,8 @@
+default: all
+
+# Define here all new external source files and objects.Don't forget to prefix the
+# object files with IRPF90_temp/
+SRC=
+OBJ=
+
+include $(QPACKAGE_ROOT)/src/Makefile.common

src/CID_SC2_selected/NEEDED_MODULES

@@ -0,0 +1 @@
+AOs BiInts Bitmask CISD SC2 CISD_selected Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils

src/CID_SC2_selected/README.rst

@@ -0,0 +1,40 @@
+========================
+CISD_SC2_selected Module
+========================
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+`cisd_sc2_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2_selected/cisd_sc2_selection.irp.f#L1>`_
+  Undocumented
+
+
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
+* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
+* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
+* `SC2 <http://github.com/LCPQ/quantum_package/tree/master/src/SC2>`_
+* `CISD_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected>`_
+* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
+* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
+* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
+* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+

src/CID_SC2_selected/cid_sc2_selection.irp.f

@@ -0,0 +1,73 @@
+program cisd_sc2_selected
+  implicit none
+  integer                        :: i,k
+  use bitmasks
+
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:),E_old(:)
+  integer                        :: N_st, iter,degree
+  character*(64)                 :: perturbation
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st),E_old(N_st))
+  
+  pt2 = 1.d0
+  perturbation = "epstein_nesbet_sc2_projected"
+  E_old(1) = HF_energy
+  davidson_threshold = 1.d-6 
+
+  do while (maxval(abs(pt2(1:N_st))) > 1.d-4)
+    print*,'----'
+    print*,''
+    call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag,  N_st)
+!    soft_touch det_connections
+    call diagonalize_CI_SC2
+    print *,  'N_det                        = ', N_det
+    do i = 1, N_st
+     print*,'state ',i
+     print *,  'PT2(SC2)                     = ', pt2(i)
+     print *,  'E(SC2)                       = ', CI_SC2_energy(i) 
+     print *,  'E_before(SC2)+PT2(SC2)       = ', (E_old(i)+pt2(i)) 
+     if(i==1)then
+      print *,  'E(SC2)+PT2(projctd)SC2       = ', (E_old(i)+H_pert_diag(i)) 
+     endif
+     E_old(i) = CI_SC2_energy(i)
+    enddo
+!   print *,  'E corr           = ', (E_old(1)) - HF_energy
+    if (abort_all) then
+      exit
+    endif
+  enddo
+  pt2 = 0.d0
+  call H_apply_PT2(pt2, norm_pert, H_pert_diag,  N_st)
+  davidson_threshold = 1.d-10
+  touch davidson_threshold davidson_criterion
+  do i = 1, N_st
+   max = 0.d0
+
+    print*,''
+    print*,'-------------'
+    print*,'for state ',i
+    print*,''
+   do k = 1, N_det
+    if(dabs(psi_coef(k,i)).gt.max)then
+     max = dabs(psi_coef(k,i))
+     imax = k
+    endif
+   enddo
+   double precision :: max
+   integer :: imax
+   print *,  'PT2(SC2)                     = ', pt2(i)
+   print *,  'E(SC2)                       = ', CI_SC2_energy(i)
+   print *,  'E_before(SC2)+PT2(SC2)       = ', (CI_SC2_energy(i)+pt2(i))
+   if(i==1)then
+    print *,  'E(SC2)+PT2(projctd)SC2       = ', (CI_SC2_energy(i)+H_pert_diag(i))
+   endif
+
+   print*,'greater coeficient of the state : ',dabs(psi_coef(imax,i))
+   call get_excitation_degree(ref_bitmask,psi_det(1,1,imax),degree,N_int)
+   print*,'degree of excitation of such determinant : ',degree
+   
+  enddo
+  print*,'coucou'
+  deallocate(pt2,norm_pert,H_pert_diag)
+end

src/CID_SC2_selected/do_mono_double.irp.f

@@ -0,0 +1,19 @@
+BEGIN_PROVIDER [logical, do_double_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the double excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_double_excitations = .True.
+
+END_PROVIDER
+
+BEGIN_PROVIDER [logical, do_mono_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the mono excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_mono_excitations = .False.
+
+END_PROVIDER

src/CID_selected/H_apply.irp.f

@@ -0,0 +1,31 @@
+use bitmasks
+BEGIN_SHELL [ /usr/bin/env python ]
+from generate_h_apply import *
+from perturbation import perturbations
+
+for perturbation in perturbations:
+  s = H_apply("cisd_selection_"+perturbation)
+  s.set_selection_pt2(perturbation)
+  print s
+END_SHELL
+
+
+subroutine H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st)
+  implicit none
+  character*(64), intent(in) :: perturbation
+  integer, intent(in)            :: N_st
+  double precision, intent(inout):: pt2(N_st), norm_pert(N_st), H_pert_diag(N_st)
+
+BEGIN_SHELL [  /usr/bin/env python ]
+from perturbation import perturbations
+
+for perturbation in perturbations:
+  print """
+  if (perturbation == '%s') then
+    call H_apply_cisd_selection_%s(pt2, norm_pert, H_pert_diag,  N_st)
+  endif
+  """%(perturbation,perturbation)
+END_SHELL
+
+
+end

src/CID_selected/Makefile

@@ -0,0 +1,8 @@
+default: all
+
+# Define here all new external source files and objects.Don't forget to prefix the
+# object files with IRPF90_temp/
+SRC=
+OBJ=
+
+include $(QPACKAGE_ROOT)/src/Makefile.common

src/CID_selected/NEEDED_MODULES

@@ -0,0 +1 @@
+AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Perturbation SingleRefMethod Utils Selectors_full

src/CID_selected/README.rst

@@ -0,0 +1,41 @@
+====================
+CISD_selected Module
+====================
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+`h_apply_cisd_selection <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/H_apply.irp.f#L13>`_
+  Undocumented
+
+`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/cisd_selection.irp.f#L1>`_
+  Undocumented
+
+
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
+* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
+* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
+* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
+* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
+* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
+* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
+* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
+* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
+* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
+

src/CID_selected/cid_selection.irp.f

@@ -0,0 +1,34 @@
+program cisd
+  implicit none
+  integer                        :: i,k
+
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:),E_old(:)
+  integer                        :: N_st, iter
+  character*(64)                 :: perturbation
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st),E_old(N_st))
+  
+  pt2 = 1.d0
+  perturbation = "epstein_nesbet"
+  E_old(1) = HF_energy
+  do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
+    print*,'----'
+    print*,''
+    call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag,  N_st)
+    call diagonalize_CI
+    print*,'N_det = ',N_det
+    do i = 1, N_st
+     print*,'state ',i
+     print *,  'PT2                          = ', pt2(i)
+     print *,  'E                            = ', CI_energy(i) 
+     print *,  'E_before     +PT2            = ', (E_old(i)+pt2(i)) 
+!   print *,  'E+PT2_new= ', (E_old(1)+1.d0*pt2(1)+H_pert_diag(1))/(1.d0 +norm_pert(1))
+    enddo
+    E_old = CI_energy
+    if (abort_all) then
+      exit
+    endif
+  enddo
+  deallocate(pt2,norm_pert,H_pert_diag)
+end

src/CID_selected/do_mono_double.irp.f

@@ -0,0 +1,19 @@
+BEGIN_PROVIDER [logical, do_double_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the double excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_double_excitations = .True.
+
+END_PROVIDER
+
+BEGIN_PROVIDER [logical, do_mono_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the mono excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_mono_excitations = .False.
+
+END_PROVIDER

src/CISD/do_mono_double.irp.f

@@ -0,0 +1,19 @@
+BEGIN_PROVIDER [logical, do_double_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the double excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_double_excitations = .True.
+
+END_PROVIDER
+
+BEGIN_PROVIDER [logical, do_mono_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the mono excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_mono_excitations = .True.
+
+END_PROVIDER

src/CISD_SC2_selected/cisd_sc2_selection.irp.f

@@ -13,9 +13,9 @@ program cisd_sc2_selected
   pt2 = 1.d0
   perturbation = "epstein_nesbet_sc2_projected"
   E_old(1) = HF_energy
-  davidson_threshold = 1.d-4
+  davidson_threshold = 1.d-6 
 
-  do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
+  do while (maxval(abs(pt2(1:N_st))) > 1.d-4)
     print*,'----'
     print*,''
     call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag,  N_st)
@@ -39,7 +39,7 @@ program cisd_sc2_selected
   enddo
   pt2 = 0.d0
   call H_apply_PT2(pt2, norm_pert, H_pert_diag,  N_st)
-  davidson_threshold = 1.d-8
+  davidson_threshold = 1.d-10
   touch davidson_threshold davidson_criterion
   do i = 1, N_st
    max = 0.d0

src/CISD_SC2_selected/do_mono_double.irp.f

@@ -0,0 +1,19 @@
+BEGIN_PROVIDER [logical, do_double_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the double excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_double_excitations = .True.
+
+END_PROVIDER
+
+BEGIN_PROVIDER [logical, do_mono_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the mono excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_mono_excitations = .True.
+
+END_PROVIDER

src/CISD_selected/do_mono_double.irp.f

@@ -0,0 +1,19 @@
+BEGIN_PROVIDER [logical, do_double_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the double excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_double_excitations = .True.
+
+END_PROVIDER
+
+BEGIN_PROVIDER [logical, do_mono_excitations]
+ implicit none
+ BEGIN_DOC
+ ! if True then the mono excitations are performed in the calculation
+ ! always true in the CISD
+ END_DOC
+ do_mono_excitations = .True.
+
+END_PROVIDER

src/Dets/H_apply_template.f

@@ -420,13 +420,17 @@ subroutine $subroutine($params_main)
       enddo
     enddo
 
+    if(do_double_excitations)then
     call $subroutine_diexc(psi_generators(1,1,i_generator),          &
         mask(1,1,d_hole1), mask(1,1,d_part1),                        &
         mask(1,1,d_hole2), mask(1,1,d_part2),                        &
         i_generator $params_post)
+    endif
+    if(do_mono_excitations)then
     call $subroutine_monoexc(psi_generators(1,1,i_generator),        &
         mask(1,1,s_hole ), mask(1,1,s_part ),                        &
         i_generator $params_post)
+    endif
     !$ call omp_set_lock(lck)
     call wall_time(wall_2)
     $printout_always
@@ -471,13 +475,17 @@ subroutine $subroutine($params_main)
             not(psi_generators(k,ispin,i_generator)) )
       enddo
     enddo
-    call $subroutine_diexc(psi_generators(1,1,i_generator),          &
-        mask(1,1,d_hole1), mask(1,1,d_part1),                        &
-        mask(1,1,d_hole2), mask(1,1,d_part2),                        &
-        i_generator $params_post)
-    call $subroutine_monoexc(psi_generators(1,1,i_generator),        &
-        mask(1,1,s_hole ), mask(1,1,s_part ),                        &
-        i_generator $params_post)
+    if(do_double_excitations)then
+     call $subroutine_diexc(psi_generators(1,1,i_generator),          &
+         mask(1,1,d_hole1), mask(1,1,d_part1),                        &
+         mask(1,1,d_hole2), mask(1,1,d_part2),                        &
+         i_generator $params_post)
+    endif
+    if(do_mono_excitations)then
+     call $subroutine_monoexc(psi_generators(1,1,i_generator),        &
+         mask(1,1,s_hole ), mask(1,1,s_part ),                        &
+         i_generator $params_post)
+    endif
     call wall_time(wall_2)
     $printout_always
     if (wall_2 - wall_0 > 2.d0) then

src/Dets/README.rst

@@ -50,24 +50,24 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L113>`_
+`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L/subroutine copy_H_apply_buffer_to_wf/;">`_
   Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det
   after calling this function.
 
-`fill_h_apply_buffer_no_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L199>`_
+`fill_h_apply_buffer_no_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L/subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)/;">`_
   Fill the H_apply buffer with determiants for CISD
 
-`h_apply_buffer_allocated <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L14>`_
+`h_apply_buffer_allocated <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L/BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ]/;">`_
   Buffer of determinants/coefficients/perturbative energy for H_apply.
   Uninitialized. Filled by H_apply subroutines.
 
-`h_apply_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L44>`_
+`h_apply_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L/BEGIN_PROVIDER [ double precision, H_apply_threshold ]/;">`_
   Theshold on | <Di|H|Dj> |
 
-`resize_h_apply_buffer <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L63>`_
+`resize_h_apply_buffer <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L/subroutine resize_H_apply_buffer(new_size,iproc)/;">`_
   Undocumented
 
-`cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/SC2.irp.f#L1>`_
+`cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/SC2.irp.f#L/subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence)/;">`_
   CISD+SC2 method              :: take off all the disconnected terms of a CISD (selected or not)
   .br
   dets_in : bitmasks corresponding to determinants
@@ -83,29 +83,29 @@ Documentation
   .br
   Initial guess vectors are not necessarily orthonormal
 
-`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/SC2.irp.f#L220>`_
+`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/SC2.irp.f#L/subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint)/;">`_
   Undocumented
 
-`connected_to_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L95>`_
+`connected_to_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L/integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)/;">`_
   Undocumented
 
-`det_is_not_or_may_be_in_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L191>`_
+`det_is_not_or_may_be_in_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L/logical function det_is_not_or_may_be_in_ref(key,Nint)/;">`_
   If true, det is not in ref
   If false, det may be in ref
 
-`is_in_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L1>`_
+`is_in_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L/logical function is_in_wavefunction(key,Nint,Ndet)/;">`_
   Undocumented
 
-`key_pattern_not_in_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L225>`_
+`key_pattern_not_in_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/connected_to_ref.irp.f#L/BEGIN_PROVIDER [ logical, key_pattern_not_in_ref, (-128:127,N_int,2) ]/;">`_
   Min and max values of the integers of the keys of the reference
 
-`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L383>`_
+`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L/subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged)/;">`_
   True if the Davidson algorithm is converged
 
-`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L373>`_
+`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L/BEGIN_PROVIDER [ character(64), davidson_criterion ]/;">`_
   Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
 
-`davidson_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L18>`_
+`davidson_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L/subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit)/;">`_
   Davidson diagonalization.
   .br
   dets_in : bitmasks corresponding to determinants
@@ -123,7 +123,7 @@ Documentation
   .br
   Initial guess vectors are not necessarily orthonormal
 
-`davidson_diag_hjj <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L68>`_
+`davidson_diag_hjj <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L/subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit)/;">`_
   Davidson diagonalization with specific diagonal elements of the H matrix
   .br
   H_jj : specific diagonal H matrix elements to diagonalize de Davidson
@@ -143,114 +143,114 @@ Documentation
   .br
   Initial guess vectors are not necessarily orthonormal
 
-`davidson_iter_max <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L1>`_
+`davidson_iter_max <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L/BEGIN_PROVIDER [ integer, davidson_iter_max ]/;">`_
   Max number of Davidson iterations
 
-`davidson_sze_max <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L9>`_
+`davidson_sze_max <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L/BEGIN_PROVIDER [ integer, davidson_sze_max ]/;">`_
   Max number of Davidson sizes
 
-`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L374>`_
+`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L/&BEGIN_PROVIDER [ double precision, davidson_threshold ]/;">`_
   Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
 
-`det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L273>`_
+`det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/integer*8 function det_search_key(det,Nint)/;">`_
   Return an integer*8 corresponding to a determinant index for searching
 
-`n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L20>`_
+`n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/BEGIN_PROVIDER [ integer, N_det ]/;">`_
   Number of determinants in the wave function
 
-`n_det_max_jacobi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L38>`_
+`n_det_max_jacobi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/BEGIN_PROVIDER [ integer, N_det_max_jacobi ]/;">`_
   Maximum number of determinants diagonalized my jacobi
 
-`n_states <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L3>`_
+`n_states <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/BEGIN_PROVIDER [ integer, N_states ]/;">`_
   Number of states to consider
 
-`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L186>`_
+`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (N_det) ]/;">`_
   Contribution of determinants to the state-averaged density
 
-`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L207>`_
+`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/&BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted, (N_det) ]/;">`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L150>`_
+`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]/;">`_
   The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file
   is empty
 
-`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L206>`_
+`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/&BEGIN_PROVIDER [ double precision, psi_coef_sorted, (N_det,N_states) ]/;">`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L237>`_
+`psi_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/&BEGIN_PROVIDER [ double precision, psi_coef_sorted_bit, (N_det,N_states) ]/;">`_
   Determinants on which we apply <i|H|psi> for perturbation.
   o They are sorted by determinants interpreted as integers. Useful
   to accelerate the search of a determinant
 
-`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L64>`_
+`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ]/;">`_
   The wave function determinants. Initialized with Hartree-Fock if the EZFIO file
   is empty
 
-`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L56>`_
+`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/BEGIN_PROVIDER [ integer, psi_det_size ]/;">`_
   Size of the psi_det/psi_coef arrays
 
-`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L205>`_
+`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted, (N_int,2,N_det) ]/;">`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L236>`_
+`psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_bit, (N_int,2,N_det) ]/;">`_
   Determinants on which we apply <i|H|psi> for perturbation.
   o They are sorted by determinants interpreted as integers. Useful
   to accelerate the search of a determinant
 
-`read_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L107>`_
+`read_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/subroutine read_dets(det,Nint,Ndet)/;">`_
   Reads the determinants from the EZFIO file
 
-`save_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L288>`_
+`save_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L/subroutine save_wavefunction/;">`_
   Save the wave function into the EZFIO file
 
-`double_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L40>`_
+`double_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L/BEGIN_PROVIDER [ integer(bit_kind), double_exc_bitmask, (N_int, 4, N_double_exc_bitmasks) ]/;">`_
   double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1
   double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1
   double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2
   double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2
   for a given couple of hole/particle excitations i.
 
-`n_double_exc_bitmasks <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L31>`_
+`n_double_exc_bitmasks <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L/BEGIN_PROVIDER [ integer, N_double_exc_bitmasks ]/;">`_
   Number of double excitation bitmasks
 
-`n_single_exc_bitmasks <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L8>`_
+`n_single_exc_bitmasks <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L/BEGIN_PROVIDER [ integer, N_single_exc_bitmasks ]/;">`_
   Number of single excitation bitmasks
 
-`single_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L17>`_
+`single_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L/BEGIN_PROVIDER [ integer(bit_kind), single_exc_bitmask, (N_int, 2, N_single_exc_bitmasks) ]/;">`_
   single_exc_bitmask(:,1,i) is the bitmask for holes
   single_exc_bitmask(:,2,i) is the bitmask for particles
   for a given couple of hole/particle excitations i.
 
-`ci_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L36>`_
+`ci_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L/&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states) ]/;">`_
   Eigenvectors/values of the CI matrix
 
-`ci_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L35>`_
+`ci_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L/BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states) ]/;">`_
   Eigenvectors/values of the CI matrix
 
-`ci_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L18>`_
+`ci_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L/BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ]/;">`_
   N_states lowest eigenvalues of the CI matrix
 
-`diag_algorithm <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L1>`_
+`diag_algorithm <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L/BEGIN_PROVIDER [ character*(64), diag_algorithm ]/;">`_
   Diagonalization algorithm (Davidson or Lapack)
 
-`diagonalize_ci <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L73>`_
+`diagonalize_ci <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI.irp.f#L/subroutine diagonalize_CI/;">`_
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 
-`ci_sc2_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L19>`_
+`ci_sc2_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L/&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states) ]/;">`_
   Eigenvectors/values of the CI matrix
 
-`ci_sc2_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L18>`_
+`ci_sc2_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L/BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states) ]/;">`_
   Eigenvectors/values of the CI matrix
 
-`ci_sc2_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L1>`_
+`ci_sc2_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L/BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states) ]/;">`_
   N_states lowest eigenvalues of the CI matrix
 
-`diagonalize_ci_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L38>`_
+`diagonalize_ci_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L/subroutine diagonalize_CI_SC2/;">`_
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 
-`filter_connected <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L2>`_
+`filter_connected <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L/subroutine filter_connected(key1,key2,Nint,sze,idx)/;">`_
   Filters out the determinants that are not connected by H
   .br
   returns the array idx which contains the index of the
@@ -261,7 +261,7 @@ Documentation
   .br
   idx(0) is the number of determinants that interact with key1
 
-`filter_connected_davidson <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L101>`_
+`filter_connected_davidson <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L/subroutine filter_connected_davidson(key1,key2,Nint,sze,idx)/;">`_
   Filters out the determinants that are not connected by H
   .br
   returns the array idx which contains the index of the
@@ -272,7 +272,7 @@ Documentation
   .br
   idx(0) is the number of determinants that interact with key1
 
-`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L233>`_
+`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L/subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)/;">`_
   returns the array idx which contains the index of the
   .br
   determinants in the array key1 that interact
@@ -281,7 +281,7 @@ Documentation
   .br
   idx(0) is the number of determinants that interact with key1
 
-`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L332>`_
+`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L/subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)/;">`_
   standard filter_connected_i_H_psi but returns in addition
   .br
   the array of the index of the non connected determinants to key1
@@ -292,69 +292,69 @@ Documentation
   .br
   to repeat the excitations
 
-`get_s2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L1>`_
+`get_s2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L/subroutine get_s2(key_i,key_j,phase,Nint)/;">`_
   Returns <S^2>
 
-`get_s2_u0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L46>`_
+`get_s2_u0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L/subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)/;">`_
   Undocumented
 
-`s_z <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L36>`_
+`s_z <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L/BEGIN_PROVIDER [ double precision, S_z ]/;">`_
   Undocumented
 
-`s_z2_sz <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L37>`_
+`s_z2_sz <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L/&BEGIN_PROVIDER [ double precision, S_z2_Sz ]/;">`_
   Undocumented
 
-`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L721>`_
+`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine a_operator(iorb,ispin,key,hjj,Nint,na,nb)/;">`_
   Needed for diag_H_mat_elem
 
-`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L766>`_
+`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine ac_operator(iorb,ispin,key,hjj,Nint,na,nb)/;">`_
   Needed for diag_H_mat_elem
 
-`decode_exc <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L76>`_
+`decode_exc <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)/;">`_
   Decodes the exc arrays returned by get_excitation.
   h1,h2 : Holes
   p1,p2 : Particles
   s1,s2 : Spins (1:alpha, 2:beta)
   degree : Degree of excitation
 
-`det_connections <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L898>`_
+`det_connections <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/BEGIN_PROVIDER [ integer*8, det_connections, (N_con_int,N_det) ]/;">`_
   .br
 
-`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L659>`_
+`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/double precision function diag_H_mat_elem(det_in,Nint)/;">`_
   Computes <i|H|i>
 
-`get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L141>`_
+`get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine get_double_excitation(det1,det2,exc,phase,Nint)/;">`_
   Returns the two excitation operators between two doubly excited determinants and the phase
 
-`get_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L30>`_
+`get_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine get_excitation(det1,det2,exc,degree,phase,Nint)/;">`_
   Returns the excitation operators between two determinants and the phase
 
-`get_excitation_degree <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L1>`_
+`get_excitation_degree <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine get_excitation_degree(key1,key2,degree,Nint)/;">`_
   Returns the excitation degree between two determinants
 
-`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L575>`_
+`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)/;">`_
   Applies get_excitation_degree to an array of determinants
 
-`get_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L274>`_
+`get_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine get_mono_excitation(det1,det2,exc,phase,Nint)/;">`_
   Returns the excitation operator between two singly excited determinants and the phase
 
-`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L814>`_
+`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine get_occ_from_key(key,occ,Nint)/;">`_
   Returns a list of occupation numbers from a bitstring
 
-`h_u_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L830>`_
+`h_u_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)/;">`_
   Computes v_0 = H|u_0>
   .br
   n : number of determinants
   .br
   H_jj : array of <j|H|j>
 
-`i_h_j <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L355>`_
+`i_h_j <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine i_H_j(key_i,key_j,Nint,hij)/;">`_
   Returns <i|H|j> where i and j are determinants
 
-`i_h_psi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L491>`_
+`i_h_psi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)/;">`_
   <key|H|psi> for the various Nstates
 
-`i_h_psi_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L527>`_
+`i_h_psi_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat)/;">`_
   <key|H|psi> for the various Nstate
   .br
   returns in addition
@@ -367,10 +367,10 @@ Documentation
   .br
   to repeat the excitations
 
-`n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L890>`_
+`n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L/BEGIN_PROVIDER [ integer, N_con_int ]/;">`_
   Number of integers to represent the connections between determinants
 
-`h_matrix_all_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/utils.irp.f#L1>`_
+`h_matrix_all_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/utils.irp.f#L/BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]/;">`_
   H matrix on the basis of the slater deter;inants defined by psi_det
 
 

src/Dets/diagonalize_CI_SC2.irp.f

@@ -15,6 +15,14 @@ BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states) ]
 
 END_PROVIDER
 
+ BEGIN_PROVIDER [ double precision, threshold_convergence_SC2]
+ implicit none
+ BEGIN_DOC
+ ! convergence of the correlation energy of SC2 iterations
+ END_DOC
+ threshold_convergence_SC2 = 1.d-8
+
+ END_PROVIDER
  BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states) ]
 &BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states) ]
   implicit none
@@ -32,7 +40,7 @@ END_PROVIDER
   
   double precision :: convergence
   call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &
-        size(CI_SC2_eigenvectors,1),N_det,N_states,N_int,davidson_threshold)
+        size(CI_SC2_eigenvectors,1),N_det,N_states,N_int,threshold_convergence_SC2)
 END_PROVIDER
 
 subroutine diagonalize_CI_SC2

src/Dets/filter_connected.irp.f

@@ -415,14 +415,14 @@ subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)
            popcnt(xor( key1(2,2,i), key2(2,2))) +                     &
            popcnt(xor( key1(3,1,i), key2(3,1))) +                     &
            popcnt(xor( key1(3,2,i), key2(3,2)))
-       if (degree_x2 < 5) then
+       if(degree_x2>6)then
+         idx_repeat(l_repeat) = i
+         l_repeat = l_repeat + 1
+       else if (degree_x2 < 5) then
          if(degree_x2 .ne. 0)then
            idx(l) = i
            l = l+1
          endif
-       elseif(degree_x2>6)then
-         idx_repeat(l_repeat) = i
-         l_repeat = l_repeat + 1
        endif
      enddo
      

src/Hartree_Fock/Fock_matrix.irp.f

@@ -217,3 +217,46 @@ BEGIN_PROVIDER [ double precision, HF_energy ]
 END_PROVIDER
 
 
+BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Fock matrix in AO basis set
+ END_DOC
+ 
+ if (elec_alpha_num == elec_beta_num) then
+ integer :: i,j
+   do j=1,ao_num
+    !DIR$ VECTOR ALIGNED
+    do i=1,ao_num_align
+     Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j)
+    enddo
+   enddo
+ else
+   double precision, allocatable :: T(:,:), M(:,:)
+   ! F_ao = S C F_mo C^t S
+   allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num))
+   call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                   &
+       ao_overlap, size(ao_overlap,1),                               &
+       mo_coef, size(mo_coef,1),                                     &
+       0.d0,                                                         &
+       M, size(M,1))
+   call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,           &
+       M, size(M,1),                                                 &
+       Fock_matrix_mo, size(Fock_matrix_mo,1),                       &
+       0.d0,                                                         &
+       T, size(T,1))
+   call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0,           &
+       T, size(T,1),                                                 &
+       mo_coef, size(mo_coef,1),                                     &
+       0.d0,                                                         &
+       M, size(M,1))
+   call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                   &
+       M, size(M,1),                                                 &
+       ao_overlap, size(ao_overlap,1),                               &
+       0.d0,                                                         &
+       Fock_matrix_ao, size(Fock_matrix_ao,1))
+
+   deallocate(T)
+ endif
+END_PROVIDER
+

src/Hartree_Fock/HF_density_matrix_ao.irp.f

@@ -1,46 +1,27 @@
- BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ]
-&BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num_align,ao_num) ]
+BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! Alpha and Beta density matrix in the AO basis
+   ! Alpha density matrix in the AO basis
    END_DOC
-   integer                        :: i,j,k,l1,l2
-   integer, allocatable           :: occ(:,:)
    
-   allocate ( occ(elec_alpha_num,2) )
-   call bitstring_to_list( HF_bitmask(1,1), occ(1,1), j, N_int)
-   ASSERT ( j==elec_alpha_num )
+   call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, &
+        mo_coef, size(mo_coef,1), &
+        mo_coef, size(mo_coef,1), 0.d0, &
+        HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1))
 
-   call bitstring_to_list( HF_bitmask(1,2), occ(1,2), j, N_int)
-   ASSERT ( j==elec_beta_num )
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num_align,ao_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Beta density matrix in the AO basis
+   END_DOC
+   
+   call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, &
+        mo_coef, size(mo_coef,1), &
+        mo_coef, size(mo_coef,1), 0.d0, &
+        HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1))
 
-   do j=1,ao_num
-     !DIR$ VECTOR ALIGNED
-     do i=1,ao_num_align
-       HF_density_matrix_ao_alpha(i,j) = 0.d0
-       HF_density_matrix_ao_beta (i,j) = 0.d0
-     enddo
-     do k=1,elec_beta_num
-       l1 = occ(k,1)
-       l2 = occ(k,2)
-       !DIR$ VECTOR ALIGNED
-       do i=1,ao_num
-         HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
-             mo_coef(i,l1) * mo_coef(j,l1)
-         HF_density_matrix_ao_beta (i,j) = HF_density_matrix_ao_beta (i,j) +&
-             mo_coef(i,l2) * mo_coef(j,l2)
-       enddo
-     enddo
-     do k=elec_beta_num+1,elec_alpha_num
-       l1 = occ(k,1)
-       !DIR$ VECTOR ALIGNED
-       do i=1,ao_num
-         HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
-             mo_coef(i,l1) * mo_coef(j,l1)
-       enddo
-     enddo
-   enddo
-   deallocate(occ)
 END_PROVIDER
  
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
@@ -48,40 +29,13 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
    BEGIN_DOC
    ! Density matrix in the AO basis
    END_DOC
-   integer                        :: i,j,k,l1,l2
-   integer, allocatable           :: occ(:,:)
+   ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1))
+   if (elec_alpha_num== elec_beta_num) then
+     HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_alpha
+   else
+     ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_beta ,1))
+     HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_beta
+   endif
    
-   allocate ( occ(elec_alpha_num,2) )
-   call bitstring_to_list( HF_bitmask(1,1), occ(1,1), j, N_int)
-   ASSERT ( j==elec_alpha_num )
-   
-   call bitstring_to_list( HF_bitmask(1,2), occ(1,2), j, N_int)
-   ASSERT ( j==elec_beta_num )
-   
-   do j=1,ao_num
-     !DIR$ VECTOR ALIGNED
-     do i=1,ao_num_align
-       HF_density_matrix_ao(i,j) = 0.d0
-     enddo
-     do k=1,elec_beta_num
-       l1 = occ(k,1)
-       l2 = occ(k,2)
-       !DIR$ VECTOR ALIGNED
-       do i=1,ao_num
-         HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) +           &
-             mo_coef(i,l1) * mo_coef(j,l1) +                         &
-             mo_coef(i,l2) * mo_coef(j,l2)
-       enddo
-     enddo
-     do k=elec_beta_num+1,elec_alpha_num
-       l1 = occ(k,1)
-       !DIR$ VECTOR ALIGNED
-       do i=1,ao_num
-         HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) +           &
-             mo_coef(i,l1) * mo_coef(j,l1)
-       enddo
-     enddo
-   enddo
-   deallocate(occ)
 END_PROVIDER
  

src/Hartree_Fock/README.rst

@@ -26,19 +26,22 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`fock_matrix_alpha_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L83>`_
+`fock_matrix_alpha_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_ao, (ao_num_align, ao_num) ]/;">`_
   Alpha Fock matrix in AO basis set
 
-`fock_matrix_alpha_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L172>`_
+`fock_matrix_alpha_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_mo, (mo_tot_num_align,mo_tot_num) ]/;">`_
   Fock matrix on the MO basis
 
-`fock_matrix_beta_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L84>`_
+`fock_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ]/;">`_
+  Fock matrix in AO basis set
+
+`fock_matrix_beta_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_ao,  (ao_num_align, ao_num) ]/;">`_
   Alpha Fock matrix in AO basis set
 
-`fock_matrix_beta_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L192>`_
+`fock_matrix_beta_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, Fock_matrix_beta_mo, (mo_tot_num_align,mo_tot_num) ]/;">`_
   Fock matrix on the MO basis
 
-`fock_matrix_diag_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L2>`_
+`fock_matrix_diag_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, Fock_matrix_diag_mo, (mo_tot_num)]/;">`_
   Fock matrix on the MO basis.
   For open shells, the ROHF Fock Matrix is
   .br
@@ -53,7 +56,7 @@ Documentation
   K = Fb - Fa
   .br
 
-`fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L1>`_
+`fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, Fock_matrix_mo, (mo_tot_num_align,mo_tot_num) ]/;">`_
   Fock matrix on the MO basis.
   For open shells, the ROHF Fock Matrix is
   .br
@@ -68,49 +71,70 @@ Documentation
   K = Fb - Fa
   .br
 
-`hf_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L211>`_
+`hf_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, HF_energy ]/;">`_
   Hartree-Fock energy
 
-`hf_density_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L46>`_
+`hf_density_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L/BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]/;">`_
   Density matrix in the AO basis
 
-`hf_density_matrix_ao_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L1>`_
-  Alpha and Beta density matrix in the AO basis
+`hf_density_matrix_ao_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L/BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ]/;">`_
+  Alpha density matrix in the AO basis
 
-`hf_density_matrix_ao_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L2>`_
-  Alpha and Beta density matrix in the AO basis
+`hf_density_matrix_ao_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L/BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num_align,ao_num) ]/;">`_
+  Beta density matrix in the AO basis
 
-`diagonal_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L1>`_
-  Diagonal Fock matrix in the MO basis
-
-`eigenvectors_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L2>`_
-  Diagonal Fock matrix in the MO basis
-
-`scf_iteration <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L1>`_
+`fock_mo_to_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)/;">`_
   Undocumented
 
-`do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L41>`_
+`insert_new_scf_density_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/subroutine insert_new_SCF_density_matrix/;">`_
+  Undocumented
+
+`it_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/BEGIN_PROVIDER [ integer, it_scf ]/;">`_
+  Number of the current SCF iteration
+
+`scf_density_matrices <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,n_it_scf_max) ]/;">`_
+  Density matrices at every SCF iteration
+
+`scf_energies <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ]/;">`_
+  Density matrices at every SCF iteration
+
+`scf_interpolation_step <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/subroutine SCF_interpolation_step/;">`_
+  Undocumented
+
+`scf_iterations <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L/subroutine scf_iterations/;">`_
+  Undocumented
+
+`diagonal_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L/BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (mo_tot_num) ]/;">`_
+  Diagonal Fock matrix in the MO basis
+
+`eigenvectors_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L/&BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num_align,mo_tot_num) ]/;">`_
+  Diagonal Fock matrix in the MO basis
+
+`xcf_iteration <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L/subroutine xcf_iteration/;">`_
+  Undocumented
+
+`do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L/BEGIN_PROVIDER [ logical, do_DIIS ]/;">`_
   If True, compute integrals on the fly
 
-`n_it_scf_max <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L22>`_
+`n_it_scf_max <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L/BEGIN_PROVIDER [ integer, n_it_scf_max]/;">`_
   Maximum number of SCF iterations
 
-`thresh_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L1>`_
+`thresh_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L/BEGIN_PROVIDER [ double precision,thresh_SCF ]/;">`_
   Threshold on the convergence of the Hartree Fock energy
 
-`bi_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L5>`_
+`bi_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/&BEGIN_PROVIDER [ double precision, bi_elec_ref_bitmask_energy ]/;">`_
   Energy of the reference bitmask used in Slater rules
 
-`kinetic_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L3>`_
+`kinetic_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/&BEGIN_PROVIDER [ double precision, kinetic_ref_bitmask_energy ]/;">`_
   Energy of the reference bitmask used in Slater rules
 
-`mono_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L2>`_
+`mono_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/&BEGIN_PROVIDER [ double precision, mono_elec_ref_bitmask_energy ]/;">`_
   Energy of the reference bitmask used in Slater rules
 
-`nucl_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L4>`_
+`nucl_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/&BEGIN_PROVIDER [ double precision, nucl_elec_ref_bitmask_energy ]/;">`_
   Energy of the reference bitmask used in Slater rules
 
-`ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L1>`_
+`ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/BEGIN_PROVIDER [ double precision, ref_bitmask_energy ]/;">`_
   Energy of the reference bitmask used in Slater rules
 
 

src/Hartree_Fock/SCF.irp.f

@@ -0,0 +1,107 @@
+BEGIN_PROVIDER [ integer, it_scf ]
+  implicit none
+  BEGIN_DOC
+  ! Number of the current SCF iteration
+  END_DOC
+  it_scf = 0
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,n_it_scf_max) ]
+&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ]
+  implicit none
+  BEGIN_DOC
+  ! Density matrices at every SCF iteration
+  END_DOC
+  SCF_density_matrices = 0.d0
+  SCF_energies = 0.d0
+END_PROVIDER
+
+subroutine insert_new_SCF_density_matrix
+  implicit none
+  integer                        :: i,j
+  do j=1,ao_num
+    do i=1,ao_num
+      SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j)
+      SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j)
+    enddo
+  enddo
+  SCF_energies(it_scf) = HF_energy
+end
+
+subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
+  implicit none
+  integer, intent(in)            :: LDFMO ! size(FMO,1)
+  integer, intent(in)            :: LDFAO ! size(FAO,1)
+  double precision, intent(in)   :: FMO(LDFMO,*)
+  double precision, intent(out)  :: FAO(LDFAO,*)
+  
+  double precision, allocatable  :: T(:,:), M(:,:)
+  ! F_ao = S C F_mo C^t S
+  allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num))
+  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
+      ao_overlap, size(ao_overlap,1),                                &
+      mo_coef, size(mo_coef,1),                                      &
+      0.d0,                                                          &
+      M, size(M,1))
+  call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
+      M, size(M,1),                                                  &
+      FMO, size(FMO,1),                                              &
+      0.d0,                                                          &
+      T, size(T,1))
+  call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0,            &
+      T, size(T,1),                                                  &
+      mo_coef, size(mo_coef,1),                                      &
+      0.d0,                                                          &
+      M, size(M,1))
+  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
+      M, size(M,1),                                                  &
+      ao_overlap, size(ao_overlap,1),                                &
+      0.d0,                                                          &
+      FAO, size(FAO,1))
+  deallocate(T,M)
+end
+
+subroutine SCF_interpolation_step
+  implicit none
+  integer                        :: i,j
+  double precision               :: c
+  
+  if (it_scf == 1) then
+    return
+  endif
+  call random_number(c)
+  do j=1,ao_num
+    do i=1,ao_num
+      HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf)+SCF_density_matrices(i,j,1,it_scf-1) * (1.d0 - c)
+      HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf)+SCF_density_matrices(i,j,2,it_scf-1) * (1.d0 - c)
+    enddo
+  enddo
+  TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta
+  
+  ! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1),&
+      !      Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) )
+  ! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1),&
+      !      Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) )
+  ! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta
+end
+
+subroutine scf_iterations
+  implicit none
+  integer                        :: i,j
+  do i=1,n_it_scf_max
+    it_scf += 1
+    SOFT_TOUCH it_scf
+    mo_coef = eigenvectors_Fock_matrix_mo
+    TOUCH mo_coef
+    call insert_new_SCF_density_matrix
+    print *,  HF_energy
+    if (SCF_energies(it_scf)>SCF_energies(it_scf-1)) then
+      call SCF_interpolation_step
+    endif
+    if (it_scf>1 ) then
+      if (dabs(SCF_energies(it_scf)-SCF_energies(it_scf-1)) < thresh_SCF) then
+        exit
+      endif
+    endif
+  enddo
+end

src/Hartree_Fock/diagonalize_fock.irp.f

@@ -5,16 +5,52 @@
    ! Diagonal Fock matrix in the MO basis
    END_DOC
    
-   double precision, allocatable  :: R(:,:)
-   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: R
+   integer                        :: i,j
+   integer                        :: liwork, lwork, n, info
+   integer, allocatable           :: iwork(:)
+   double precision, allocatable  :: work(:), F(:,:), S(:,:)
    
-   allocate(R(mo_tot_num_align,mo_tot_num))
+   allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
+   do j=1,ao_num
+     do i=1,ao_num
+       S(i,j) = ao_overlap(i,j)
+       F(i,j) = Fock_matrix_ao(i,j)
+     enddo
+   enddo
 
-   call lapack_diag(diagonal_Fock_matrix_mo,R,Fock_matrix_mo,size(Fock_matrix_mo,1),&
-       mo_tot_num)
+   n = ao_num
+   lwork = 1+6*n + 2*n*n
+   liwork = 3 + 5*n
    
-   call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num,1.d0,mo_coef,size(mo_coef,1),&
-       R,size(R,1),0.d0,eigenvectors_Fock_matrix_mo,size(eigenvectors_Fock_matrix_mo,1))
-   deallocate(R)
+   allocate(work(lwork), iwork(liwork) )
+
+   lwork = -1
+   liwork = -1
+   call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),&
+     diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
+
+   if (info /= 0) then
+     print *,  irp_here//' failed'
+     stop 1
+   endif
+   lwork = int(work(1))
+   liwork = iwork(1)
+   deallocate(work,iwork)
+   allocate(work(lwork), iwork(liwork) )
+
+   call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),&
+     diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
+
+   if (info /= 0) then
+     print *,  irp_here//' failed'
+     stop 1
+   endif
+   do j=1,ao_num
+     do i=1,ao_num
+       eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
+     enddo
+   enddo
+
+   deallocate(work, iwork, F, S)
 END_PROVIDER
  

src/Hartree_Fock/mo_SCF_iterations.irp.f

@@ -1,71 +1,15 @@
-program scf_iteration
+program xcf_iteration
   use bitmasks
   implicit none
   double precision               :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral
   double precision               :: E0
   integer                        :: i_it, i, j, k
-  integer, allocatable           :: iorder(:)
-    double precision, allocatable :: DM_occ(:,:), E_new(:), R(:,:)
    
   E0 = HF_energy 
-  i_it = 1
-  n_it_scf_max = 100
-  SCF_energy_before = 0.d0
-  SCF_energy_after = E0
-  print *,  E0
-  mo_label = "Canonical"
-  thresh_SCF = 1.d-10
-  DM_occ = mo_density_matrix
-  allocate (DM_occ(size(mo_density_matrix,1),mo_tot_num), &
-            E_new(mo_tot_num), R(mo_tot_num,mo_tot_num), iorder(mo_tot_num))
-  do while (i_it < n_it_scf_max .and. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF)
-    if (SCF_energy_after <= SCF_energy_before+1.d-4) then
-      mo_coef = eigenvectors_Fock_matrix_mo
-      TOUCH mo_coef mo_label
-      DM_occ = mo_density_matrix
-    else
-      DM_occ = mo_density_matrix
-      mo_coef = eigenvectors_Fock_matrix_mo
-      TOUCH mo_coef mo_label mo_integrals_map
-      DM_occ = DM_occ + 0.0d0*mo_density_matrix
-      integer :: rank
-      call cholesky_mo(ao_num,mo_tot_num,DM_occ,size(DM_occ,1),mo_coef,size(mo_coef,1),-1.d0,rank)
-      print *,  rank
-      TOUCH mo_coef mo_label
-      call orthonormalize_mos
-      call find_rotation(eigenvectors_Fock_matrix_mo,mo_tot_num_align,mo_coef,ao_num,R, mo_tot_num)
-       do i=1,mo_tot_num
-         iorder(i) = i
-         E_new(i) = 0.d0
-         do k=1,mo_tot_num
-           E_new(i) += R(k,i)*R(k,i)*diagonal_fock_matrix_mo(k)
-         enddo
-       enddo
-       call dsort(E_new(1),iorder(1),mo_tot_num)
-       eigenvectors_Fock_matrix_mo = mo_coef
-       do j=1,mo_tot_num
-         do i=1,ao_num
-           mo_coef(i,j) = eigenvectors_Fock_matrix_mo(i,iorder(j))
-         enddo
-       enddo
-      TOUCH mo_coef mo_label mo_integrals_map
-    endif
-    call clear_mo_map
-    SCF_energy_before = SCF_energy_after
-    SCF_energy_after = HF_energy
-    print*,SCF_energy_after, dabs(SCF_energy_before - SCF_energy_after) 
-    i_it +=1
-    if(i_it > n_it_scf_max)exit
-  enddo
   
-  if (i_it >= n_it_scf_max) then
-    stop 'Failed'
-  endif
-  if (SCF_energy_after - E0 > thresh_SCF) then
-    stop 'Failed'
-  endif
+  thresh_SCF = 1.d-10
+  call scf_iterations
   mo_label = "Canonical"
-  deallocate (DM_occ)
   TOUCH mo_label mo_coef
   call save_mos
   

src/MOGuess/README.rst

@@ -22,19 +22,19 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`h_core_guess <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/H_CORE_guess.irp.f#L1>`_
+`h_core_guess <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/H_CORE_guess.irp.f#L/subroutine h_core_guess/;">`_
   Undocumented
 
-`ao_ortho_lowdin_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L2>`_
+`ao_ortho_lowdin_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L/BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)]/;">`_
   matrix of the coefficients of the mos generated by the
   orthonormalization by the S^{-1/2} canonical transformation of the aos
   ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital
 
-`ao_ortho_lowdin_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L26>`_
+`ao_ortho_lowdin_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L/BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)]/;">`_
   overlap matrix of the ao_ortho_lowdin
   supposed to be the Identity
 
-`ao_ortho_lowdin_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f#L1>`_
+`ao_ortho_lowdin_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f#L/BEGIN_PROVIDER [double precision, ao_ortho_lowdin_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)]/;">`_
   Undocumented
 
 

src/MOs/cholesky_mo.irp.f

@@ -1,7 +1,7 @@
 subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
  implicit none
  BEGIN_DOC
-! Cholesky decomposition of MO Density matrix to
+! Cholesky decomposition of AO Density matrix to
 ! generate MOs
  END_DOC
  integer, intent(in) :: n,m, LDC, LDP

src/MOs/mos.ezfio_config

@@ -1,4 +1,5 @@
 mo_basis
+  ao_num             integer
   mo_tot_num         integer
   mo_coef            double precision  (ao_basis_ao_num,mo_basis_mo_tot_num)
   mo_label           character*(64)

src/Perturbation/README.rst

@@ -82,7 +82,7 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`pt2_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/Moller_plesset.irp.f#L1>`_
+`pt2_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/Moller_plesset.irp.f#L/subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)/;">`_
   compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
   .br
   for the various n_st states.
@@ -92,7 +92,7 @@ Documentation
   e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies)
   .br
 
-`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L1>`_
+`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L/subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)/;">`_
   compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
   .br
   for the various N_st states.
@@ -102,7 +102,7 @@ Documentation
   e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
   .br
 
-`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L40>`_
+`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L/subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)/;">`_
   compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
   .br
   for the various N_st states.
@@ -112,20 +112,46 @@ Documentation
   c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
   .br
 
-`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L1>`_
+`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/pert_sc2.irp.f#L/subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)/;">`_
+  compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  .br
+  for the various N_st states,
+  .br
+  but  with the correction in the denominator
+  .br
+  comming from the interaction of that determinant with all the others determinants
+  .br
+  that can be repeated by repeating all the double excitations
+  .br
+  : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
+  .br
+  that could be repeated to this determinant.
+  .br
+  In addition, for the perturbative energetic contribution you have the standard second order
+  .br
+  e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
+  .br
+  and also the purely projected contribution
+  .br
+  H_pert_diag = <HF|H|det_pert> c_pert
+
+`repeat_all_e_corr <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/pert_sc2.irp.f#L/double precision function repeat_all_e_corr(key_in)/;">`_
+  Undocumented
+
+`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,coef_pert_buffer, &>`_
   Fill the H_apply buffer with determiants for the selection
 
-`remove_small_contributions <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L81>`_
+`remove_small_contributions <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/subroutine remove_small_contributions/;">`_
   Remove determinants with small contributions. N_states is assumed to be
   provided.
 
-`selection_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L68>`_
+`selection_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/BEGIN_PROVIDER [ double precision, selection_criterion ]/;">`_
   Threshold to select determinants. Set by selection routines.
 
-`selection_criterion_factor <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L70>`_
+`selection_criterion_factor <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/&BEGIN_PROVIDER [ double precision, selection_criterion_factor ]/;">`_
   Threshold to select determinants. Set by selection routines.
 
-`selection_criterion_min <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L69>`_
+`selection_criterion_min <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/&BEGIN_PROVIDER [ double precision, selection_criterion_min ]/;">`_
   Threshold to select determinants. Set by selection routines.