improvement in casscf with CISD, CISDTQ and so on

src/casscf/EZFIO.cfg

@@ -10,4 +10,10 @@ doc: Calculated |FCI| energy + |PT2|
 interface: ezfio
 size: (determinants.n_states)
 
+[cisd_guess]
+type: logical
+doc: If true, the CASSCF starts with a CISD wave function
+interface: ezfio,provider,ocaml
+default: True 
+
 

src/casscf/NEED

@@ -1,4 +1,4 @@
 cipsi
 selectors_full
-generators_cas 
+generators_fluid
 two_body_rdm

src/casscf/casscf.irp.f

@@ -4,45 +4,86 @@ program casscf
 ! TODO : Put the documentation of the program here
   END_DOC
   no_vvvv_integrals = .True.
-  pt2_max = 0.02
-  SOFT_TOUCH no_vvvv_integrals pt2_max
-  call run
-end
-
-subroutine run
-  implicit none
-  double precision               :: energy_old, energy
+  SOFT_TOUCH no_vvvv_integrals 
+  threshold_davidson = 1.d-7
+  touch threshold_davidson 
+  if(cisd_guess)then
     logical :: converged 
     integer :: iteration
+    double precision :: energy
+    print*,'*******************************'
+    print*,'*******************************'
+    print*,'*******************************'
+    print*,'USING A CISD WAVE FUNCTION AS GUESS FOR THE MCSCF WF'
+    print*,'*******************************'
+    print*,'*******************************'
     converged = .False.
-
-  energy = 0.d0
-  mo_label = "MCSCF"
-  iteration = 1
+    iteration = 0
+    generators_type = "HF"
+    touch generators_type
+    read_wf = .False.
+    touch read_wf 
+    logical :: do_cisdtq
+    do_cisdtq  = .True.
+    double precision :: thr
+    thr = 5.d-3
     do while (.not.converged)
-    call run_stochastic_cipsi
-    energy_old = energy
-    energy = eone+etwo+ecore
-
-    call write_time(6)
-    call write_int(6,iteration,'CAS-SCF iteration')
-    call write_double(6,energy,'CAS-SCF energy')
-    call write_double(6,energy_improvement, 'Predicted energy improvement')
-
-    converged = dabs(energy_improvement) < thresh_scf
-    pt2_max = dabs(energy_improvement / pt2_relative_error)
-
-    mo_coef = NewOrbs
-    call save_mos
-    call map_deinit(mo_integrals_map)
-    iteration += 1
-    N_det = N_det/2
-    psi_det = psi_det_sorted
-    psi_coef = psi_coef_sorted
-    read_wf = .True.
-    FREE mo_integrals_map mo_two_e_integrals_in_map
-    SOFT_TOUCH mo_coef N_det pt2_max  psi_det psi_coef
-
+     call cisd_scf_iteration(converged,iteration,energy,thr)
+     if(HF_index.ne.1.and.iteration.gt.0)then
+      print*,'*******************************'
+      print*,'*******************************'
+      print*,'The HF determinant is not the dominant determinant in the CISD WF ...'
+      print*,'Therefore we skip the CISD WF ..'
+      print*,'*******************************'
+      print*,'*******************************'
+      do_cisdtq = .False.
+      exit
+     endif
+     if(iteration.gt.15.and..not.converged)then
+      print*,'It seems that the orbital optimization for the CISD WAVE FUNCTION CANNOT CONVERGE ...' 
+      print*,'Passing to CISDTQ WAVE FUNCTION'
+      exit
+     endif
     enddo
-
+    if(do_cisdtq)then
+     print*,'*******************************'
+     print*,'*******************************'
+     print*,'*******************************'
+     print*,'SWITCHING WITH A CISDTQ WAVE FUNCTION AS GUESS FOR THE MCSCF WF'
+     print*,'*******************************'
+     print*,'*******************************'
+     converged = .False.
+     iteration = 0
+     read_wf = .False.
+     touch read_wf 
+     pt2_max = 0.01d0
+     touch pt2_max 
+     energy = 0.d0
+     do while (.not.converged)
+      call cisdtq_scf_iteration(converged,iteration,energy,thr)
+      if(HF_index.ne.1.and.iteration.gt.0)then
+       print*,'*******************************'
+       print*,'*******************************'
+       print*,'The HF determinant is not the dominant determinant in the CISDTQ WF ...'
+       print*,'Therefore we skip the CISDTQ WF ..'
+       print*,'*******************************'
+       print*,'*******************************'
+       exit
+      endif
+      if(iteration.gt.15.and..not.converged)then
+       print*,'It seems that the orbital optimization for the CISDTQ WAVE FUNCTION CANNOT CONVERGE ...' 
+       print*,'Passing to CISDTQ WAVE FUNCTION'
+       exit
+      endif
+     enddo
+    endif
+  endif
+  generators_type = "CAS"
+  touch generators_type
+  read_wf = .False.
+  touch read_wf 
+  pt2_max = 0.015d0
+  touch pt2_max 
+  call run_cipsi_scf
 end
+

src/casscf/change_bitmasks.irp.f

@@ -0,0 +1,14 @@
+subroutine only_act_bitmask
+ implicit none
+  integer :: i,j,k
+  do k = 1, N_generators_bitmask
+   do j = 1, 6
+    do i = 1, N_int
+    generators_bitmask(i,1,j,k) = act_bitmask(i,1)
+    generators_bitmask(i,2,j,k) = act_bitmask(i,2)
+    enddo
+   enddo
+  enddo
+  touch generators_bitmask
+end
+

src/casscf/cipsi_routines.irp.f

@@ -0,0 +1,72 @@
+subroutine run_cipsi_scf
+  implicit none
+  double precision               :: energy_old, energy, extrap,extrap_old,pt2_max_begin
+  logical                        :: converged
+  integer                        :: iteration
+  print*,'*********************************'
+  print*,'*********************************'
+  print*,'   DOING THE CIPSI-SCF '
+  print*,'*********************************'
+  converged = .False.
+  pt2_max_begin = pt2_max
+  energy = 0.d0
+  extrap = 0.d0
+  mo_label = "MCSCF"
+  iteration = 1
+  threshold_davidson = 1.d-09
+  touch threshold_davidson 
+  do while (.not.converged)
+    print*,''
+    call write_int(6,iteration,'CI STEP OF THE ITERATION = ')
+    call write_double(6,pt2_max,'PT2 MAX = ')
+    call run_stochastic_cipsi
+    call change_orb_cipsi(converged,iteration,energy)
+    if(iteration.gt.n_it_scf_max.and..not.converged)then
+     print*,'It seems that the orbital optimization for the CISDTQ WAVE FUNCTION CANNOT CONVERGE ...' 
+     print*,'The required delta E was :',thresh_scf
+     print*,'The obtained delta E was :',extrap - extrap_old
+     print*,'After ',iteration,'iterations ...'
+     print*,'Getting out of the SCF loop ...'
+     exit
+    endif
+    iteration += 1
+  enddo
+
+end
+
+subroutine change_orb_cipsi(converged,iteration,energy)
+  implicit none
+  double precision               :: energy_old, extrap,extrap_old,pt2_max_begin
+  double precision, intent(inout):: energy
+  logical, intent(out)           :: converged
+  integer, intent(in)            :: iteration
+    extrap_old = energy
+    energy = eone+etwo+ecore
+    extrap = extrapolated_energy(2,1)
+
+    call write_time(6)
+    call write_int(6,iteration,'CAS-SCF iteration')
+    call write_double(6,energy,'CAS-SCF variational energy')
+    call write_double(6,extrap,'CAS-SCF extrapolated energy')
+    call write_double(6,extrap - extrap_old,'Change in extrapolated energy')
+    energy = extrap 
+    call write_double(6,energy_improvement, 'Predicted energy improvement')
+
+    converged = dabs(extrap - extrap_old) < thresh_scf
+    pt2_max = dabs(extrap - extrap_old) * 10.d0
+    pt2_max = min(pt2_max,1.d-2)
+    pt2_max = max(pt2_max,1.d-10)
+    if(N_det.gt.10**6)then
+     pt2_max = max(pt2_max,1.d-2)
+    endif
+
+    mo_coef = NewOrbs
+    call save_mos
+    call map_deinit(mo_integrals_map)
+    N_det = N_det/2
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    read_wf = .True.
+    FREE mo_integrals_map mo_two_e_integrals_in_map
+    SOFT_TOUCH mo_coef N_det pt2_max  psi_det psi_coef
+end

src/casscf/cisd_routine.irp.f

@@ -0,0 +1,64 @@
+subroutine cisd_scf_iteration(converged,iteration,energy,thr)
+ implicit none 
+ double precision, intent(in) :: thr
+ logical, intent(out) :: converged
+ integer, intent(inout) :: iteration
+ double precision, intent(out) :: energy
+ converged = .False.
+ call only_act_bitmask
+ call run_cisd
+ call change_orb_cisd(converged,iteration,energy,thr)
+end
+
+subroutine change_orb_cisd(converged,iteration,energy,thr)
+ implicit none 
+ double precision, intent(in) :: thr
+ logical, intent(inout) :: converged
+ integer, intent(inout) :: iteration
+ double precision, intent(inout) :: energy
+ double precision :: energy_old
+ energy_old = energy
+
+ energy = eone+etwo+ecore
+
+ call write_time(6)
+ call write_int(6,iteration,'CISD-SCF iteration')
+ call write_double(6,energy,'CISD-SCF energy')
+ call write_double(6,energy_improvement, 'Predicted energy improvement')
+ converged = dabs(energy_improvement) < thr
+
+ mo_coef = NewOrbs
+ call save_mos
+ call map_deinit(mo_integrals_map)
+ FREE mo_integrals_map mo_two_e_integrals_in_map
+ iteration += 1
+
+end
+
+subroutine run_cisd
+  implicit none
+  integer :: i
+
+  if(pseudo_sym)then
+   call H_apply_cisd_sym
+  else
+   call H_apply_cisd
+  endif
+  print *,  'N_det = ', N_det
+  print*,'******************************'
+  print *,  'Energies  of the states:'
+  do i = 1,N_states
+    print *,  i, CI_energy(i)
+  enddo
+  if (N_states > 1) then
+    print*,'******************************'
+    print*,'Excitation energies '
+    do i = 2, N_states
+      print*, i ,CI_energy(i) - CI_energy(1)
+    enddo
+  endif
+  psi_coef = ci_eigenvectors
+  SOFT_TOUCH psi_coef
+  call save_wavefunction
+
+end

src/casscf/cisdtq_routine.irp.f

@@ -0,0 +1,47 @@
+subroutine cisdtq_scf_iteration(converged,iteration,energy,thr)
+ implicit none 
+ double precision, intent(in) :: thr
+ logical, intent(out) :: converged
+ integer, intent(inout) :: iteration
+ double precision, intent(inout) :: energy
+ converged = .False.
+ call only_act_bitmask
+ generators_type = "HF_SD"
+ threshold_generators = 0.99d0
+ touch threshold_generators
+ touch generators_type
+ selection_factor = 5
+ touch selection_factor
+ call run_stochastic_cipsi
+ call change_orb_cisdtq(converged,iteration,energy,thr)
+end
+
+subroutine change_orb_cisdtq(converged,iteration,energy,thr)
+ implicit none 
+ double precision, intent(in) :: thr
+ logical, intent(inout) :: converged
+ integer, intent(inout) :: iteration
+ double precision, intent(inout) :: energy
+ double precision :: extrap,extrap_old,pt2_max_begin
+ extrap_old = energy 
+ extrap = extrapolated_energy(2,1)
+ energy = extrap
+
+ call write_time(6)
+ call write_int(6,iteration,'CISDTQ-SCF iteration')
+ call write_double(6,energy,'CISDTQ-SCF variational energy')
+ call write_double(6,extrap,'CISDTQ-SCF extrapolated energy')
+ call write_double(6,extrap - extrap_old,'Change in extrapolated energy')
+
+ converged = dabs(extrap - extrap_old) < thr
+ pt2_max = dabs(extrap - extrap_old) * 10.d0
+ pt2_max = max(pt2_max,1.d-10)
+
+ mo_coef = NewOrbs
+ call save_mos
+ call map_deinit(mo_integrals_map)
+ FREE mo_integrals_map mo_two_e_integrals_in_map
+ iteration += 1
+
+end
+

src/casscf/gradient.irp.f

@@ -171,11 +171,11 @@ BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
     norm_grad+=gradvec2(indx)*gradvec2(indx)
   end do
   norm_grad=sqrt(norm_grad)
-  if (bavard) then
+! if (bavard) then
     write(6,*)
     write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
     write(6,*)
-  endif
+! endif
   
 END_PROVIDER
 

src/casscf/h_apply.irp.f

@@ -0,0 +1,18 @@
+! Generates subroutine H_apply_cisd
+! ----------------------------------
+
+BEGIN_SHELL [ /usr/bin/env python2 ]
+from generate_h_apply import H_apply
+H = H_apply("cisd",do_double_exc=True)
+print H
+
+from generate_h_apply import H_apply
+H = H_apply("cisdtq",do_double_exc=True)
+H.set_selection_pt2("epstein_nesbet_2x2")
+print H
+
+H = H_apply("cisd_sym",do_double_exc=True)
+H.filter_only_connected_to_hf()
+print H
+END_SHELL
+

src/casscf/neworbs.irp.f

@@ -66,6 +66,7 @@ END_PROVIDER
   integer                        :: best_vector
   real*8                         :: best_overlap
   best_overlap=0.D0
+  best_vector = -1000
   do i=1,nMonoEx+1
     if (SXeigenval(i).lt.0.D0) then
       if (abs(SXeigenvec(1,i)).gt.best_overlap) then
@@ -74,7 +75,9 @@ END_PROVIDER
       end if
     end if
   end do
-  
+  if(best_vector.lt.0)then 
+   best_vector = minloc(SXeigenval,nMonoEx+1) 
+  endif
   energy_improvement = SXeigenval(best_vector)
   
   c0=SXeigenvec(1,best_vector)

src/cisd/cisd.irp.f

@@ -46,34 +46,6 @@ program cisd
   END_DOC
   read_wf = .False.
   SOFT_TOUCH read_wf
-  call run
-end
-
-subroutine run
-  implicit none
-  integer :: i
-
-  if(pseudo_sym)then
-   call H_apply_cisd_sym
-  else
-   call H_apply_cisd
-  endif
-  print *,  'N_det = ', N_det
-  print*,'******************************'
-  print *,  'Energies  of the states:'
-  do i = 1,N_states
-    print *,  i, CI_energy(i)
-  enddo
-  if (N_states > 1) then
-    print*,'******************************'
-    print*,'Excitation energies '
-    do i = 2, N_states
-      print*, i ,CI_energy(i) - CI_energy(1)
-    enddo
-  endif
-  psi_coef = ci_eigenvectors
-  SOFT_TOUCH psi_coef
-  call save_wavefunction
-  call ezfio_set_cisd_energy(CI_energy)
-
+  call only_act_bitmask
+  call run_cisd
 end

src/cisd/cisd_routine.irp.f

@@ -0,0 +1,42 @@
+subroutine only_act_bitmask
+ implicit none
+  integer :: i,j,k
+  do k = 1, N_generators_bitmask
+   do j = 1, 6
+    do i = 1, N_int
+    generators_bitmask(i,1,j,k) = act_bitmask(i,1)
+    generators_bitmask(i,2,j,k) = act_bitmask(i,2)
+    enddo
+   enddo
+  enddo
+  touch generators_bitmask
+end
+
+subroutine run_cisd
+  implicit none
+  integer :: i
+
+  if(pseudo_sym)then
+   call H_apply_cisd_sym
+  else
+   call H_apply_cisd
+  endif
+  print *,  'N_det = ', N_det
+  print*,'******************************'
+  print *,  'Energies  of the states:'
+  do i = 1,N_states
+    print *,  i, CI_energy(i)
+  enddo
+  if (N_states > 1) then
+    print*,'******************************'
+    print*,'Excitation energies '
+    do i = 2, N_states
+      print*, i ,CI_energy(i) - CI_energy(1)
+    enddo
+  endif
+  psi_coef = ci_eigenvectors
+  SOFT_TOUCH psi_coef
+  call save_wavefunction
+  call ezfio_set_cisd_energy(CI_energy)
+
+end

src/generators_fluid/NEED

@@ -0,0 +1 @@
+determinants

src/generators_fluid/extract_cas.irp.f

@@ -0,0 +1,23 @@
+subroutine extract_cas
+  implicit none
+  BEGIN_DOC
+  ! Replaces the total wave function by the normalized projection on the CAS.
+  END_DOC
+
+  integer                        :: i,j,k
+  do k=1,N_states
+    do j=1,N_det_generators
+      psi_coef(j,k) = psi_coef_generators(j,k)
+    enddo
+  enddo
+
+  do j=1,N_det_generators
+    do k=1,N_int
+      psi_det(k,1,j) = psi_det_generators(k,1,j)
+      psi_det(k,2,j) = psi_det_generators(k,2,j)
+    enddo
+  enddo
+  N_det = N_det_generators
+
+  SOFT_TOUCH N_det psi_det psi_coef
+end

src/generators_fluid/generators.irp.f

@@ -0,0 +1,101 @@
+use bitmasks
+
+BEGIN_PROVIDER [ character*(32), generators_type]
+ implicit none
+ generators_type = trim("CAS")
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ integer, N_det_generators ]
+  implicit none
+  BEGIN_DOC
+  ! Number of generator detetrminants
+  END_DOC
+  if(generators_type == "CAS")then
+   N_det_generators = N_det_generators_CAS 
+  else if (generators_type == "HF")then
+   N_det_generators = N_det_generators_HF  
+  else if (generators_type == "HF_SD")then
+   N_det_generators = N_det_generators_HF_SD
+  endif
+  N_det_generators = max(N_det_generators,1)
+  call write_int(6,N_det_generators,'Number of generators')
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ 
+ if(generators_type == "CAS")then
+  psi_det_generators(1:N_int,1:2,1:N_det_generators_CAS) = psi_det_generators_CAS(1:N_int,1:2,1:N_det_generators_CAS)
+  psi_coef_generators(1:N_det_generators_CAS,1:N_states) = psi_coef_generators_CAS(1:N_det_generators_CAS,1:N_states)
+ else if (generators_type == "HF")then
+  psi_det_generators(1:N_int,1:2,1:N_det_generators_HF) = psi_det_generators_HF(1:N_int,1:2,1:N_det_generators_HF)
+  psi_coef_generators(1:N_det_generators_HF,1:N_states) = psi_coef_generators_HF(1:N_det_generators_HF,1:N_states)
+ else if (generators_type == "HF_SD")then
+  psi_det_generators(1:N_int,1:2,1:N_det_generators_HF_SD) = psi_det_generators_HF_SD(1:N_int,1:2,1:N_det_generators_HF_SD)
+  psi_coef_generators(1:N_det_generators_HF_SD,1:N_states) = psi_coef_generators_HF_SD(1:N_det_generators_HF_SD,1:N_states)
+ endif
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_order,     (psi_det_size)  ]
+
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+  if(generators_type == "CAS")then
+   psi_det_sorted_gen = psi_det_sorted_gen_CAS
+   psi_coef_sorted_gen = psi_coef_sorted_gen_CAS
+   psi_det_sorted_gen_order = psi_det_sorted_gen_CAS_order
+  else if(generators_type == "HF")then
+   psi_det_sorted_gen = 0_bit_kind
+   psi_coef_sorted_gen = 0.d0
+   psi_det_sorted_gen_order = 0
+  else if(generators_type == "HF_SD")then
+   psi_det_sorted_gen = psi_det_sorted_gen_HF_SD
+   psi_coef_sorted_gen = psi_coef_sorted_gen_HF_SD
+   psi_det_sorted_gen_order = psi_det_sorted_gen_HF_SD_order
+  endif
+END_PROVIDER
+
+
+BEGIN_PROVIDER [integer, degree_max_generators]
+ implicit none
+ BEGIN_DOC
+! Max degree of excitation (respect to HF) of the generators
+ END_DOC
+ integer :: i,degree
+  degree_max_generators = 0
+  do i = 1, N_det_generators
+   call get_excitation_degree(HF_bitmask,psi_det_generators(1,1,i),degree,N_int)
+   if(degree .gt. degree_max_generators)then
+    degree_max_generators = degree
+   endif
+  enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, size_select_max]
+ implicit none
+ BEGIN_DOC
+ ! Size of the select_max array
+ END_DOC
+ size_select_max = 10000
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ]
+ implicit none
+ BEGIN_DOC
+ ! Memo to skip useless selectors
+ END_DOC
+ select_max = huge(1.d0)
+END_PROVIDER
+

src/generators_fluid/generators_cas.irp.f

@@ -0,0 +1,69 @@
+use bitmasks
+
+BEGIN_PROVIDER [ integer, N_det_generators_CAS ]
+  implicit none
+  BEGIN_DOC
+  ! Number of generator detetrminants
+  END_DOC
+  integer                        :: i,k,l
+  logical                        :: good
+  integer, external :: number_of_holes,number_of_particles
+  call write_time(6)
+  N_det_generators_CAS = 0
+  do i=1,N_det
+    good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 )
+    if (good) then
+      N_det_generators_CAS += 1
+    endif
+  enddo
+  N_det_generators_CAS = max(N_det_generators_CAS,1)
+  call write_int(6,N_det_generators_CAS,'Number of generators_CAS')
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_CAS, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators_CAS, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen_CAS, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen_CAS, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_CAS_order, (psi_det_size) ]
+  implicit none
+  BEGIN_DOC
+  ! For Single reference wave functions, the gen_CASerator is the
+  ! Hartree-Fock determinant
+  END_DOC
+  integer                        :: i, k, l, m
+  logical                        :: good
+  integer, external :: number_of_holes,number_of_particles
+  integer, allocatable :: nongen_CAS(:)
+  integer :: inongen_CAS
+
+  allocate(nongen_CAS(N_det))
+
+  inongen_CAS = 0
+  m=0
+  do i=1,N_det
+    good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 )
+    if (good) then
+      m = m+1
+      psi_det_sorted_gen_CAS_order(i) = m
+      do k=1,N_int
+        psi_det_generators_CAS(k,1,m) = psi_det_sorted(k,1,i)
+        psi_det_generators_CAS(k,2,m) = psi_det_sorted(k,2,i)
+      enddo
+      psi_coef_generators_CAS(m,:) = psi_coef_sorted(i,:)
+    else
+      inongen_CAS += 1
+      nongen_CAS(inongen_CAS) = i
+    endif
+  enddo
+  ASSERT (m == N_det_generators_CAS)
+
+  psi_det_sorted_gen_CAS(:,:,:N_det_generators_CAS) = psi_det_generators_CAS(:,:,:N_det_generators_CAS)
+  psi_coef_sorted_gen_CAS(:N_det_generators_CAS, :) = psi_coef_generators_CAS(:N_det_generators_CAS, :)
+  do i=1,inongen_CAS
+    psi_det_sorted_gen_CAS_order(nongen_CAS(i)) = N_det_generators_CAS+i
+    psi_det_sorted_gen_CAS(:,:,N_det_generators_CAS+i) = psi_det_sorted(:,:,nongen_CAS(i))
+    psi_coef_sorted_gen_CAS(N_det_generators_CAS+i, :) = psi_coef_sorted(nongen_CAS(i),:)
+  end do
+
+END_PROVIDER
+

src/generators_fluid/generators_hf.irp.f

@@ -0,0 +1,51 @@
+
+use bitmasks
+
+BEGIN_PROVIDER [ integer, N_det_generators_HF ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the number of generators is 1 : the
+ ! Hartree-Fock determinant
+ END_DOC
+ N_det_generators_HF = 1
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_HF, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators_HF, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_det_generators_HF = 0_bit_kind
+ integer :: i,j
+ integer :: degree
+
+ do i=1,N_int
+   psi_det_generators_HF(i,1,1) = HF_bitmask(i,1)
+   psi_det_generators_HF(i,2,1) = HF_bitmask(i,2)
+ enddo
+
+ do j=1,N_det
+   call get_excitation_degree(HF_bitmask,psi_det(1,1,j),degree,N_int)
+   if (degree == 0) then
+     exit
+   endif
+ end do
+
+ psi_det_generators_HF(:,:,1) = psi_det(:,:,j)
+ psi_coef_generators_HF(1,:) = psi_coef_generators_HF(j,:)
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer          , HF_index ]
+ implicit none
+ integer :: j,degree
+ do j=1,N_det
+   call get_excitation_degree(HF_bitmask,psi_det_sorted(1,1,j),degree,N_int)
+   if (degree == 0) then
+     HF_index = j 
+     exit
+   endif
+ end do
+END_PROVIDER 

src/generators_fluid/generators_hf_sd.irp.f

@@ -0,0 +1,80 @@
+
+use bitmasks
+
+BEGIN_PROVIDER [ integer, N_det_generators_HF_SD ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the number of generators is 1 : the
+ ! Hartree-Fock determinant
+ END_DOC
+ N_det_generators_HF_SD = 0
+ integer :: i,degree
+ double precision :: thr
+ double precision :: accu
+ accu = 0.d0
+ thr = threshold_generators
+ do i = 1, N_det
+  call get_excitation_degree(HF_bitmask,psi_det_sorted(1,1,i),degree,N_int)
+  if(degree.le.2.and. accu .le. thr )then
+   accu += psi_coef_sorted(i,1)**2
+   N_det_generators_HF_SD += 1
+  endif
+ enddo
+!print*,''
+!print*,'N_det_generators_HF_SD = ',N_det_generators_HF_SD
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_HF_SD, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators_HF_SD, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen_HF_SD, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen_HF_SD, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_HF_SD_order, (psi_det_size) ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_det_generators_HF_SD = 0_bit_kind
+ integer :: i,j,k
+ integer :: degree
+ double precision :: thr
+ double precision :: accu
+ integer, allocatable :: nongen(:)
+ integer :: inongen
+
+ allocate(nongen(N_det))
+
+ thr = threshold_generators
+
+ accu = 0.d0
+ k = 0
+ inongen = 0
+ do j=1,N_det
+   call get_excitation_degree(HF_bitmask,psi_det_sorted(1,1,j),degree,N_int)
+   if(degree.le.2.and. accu.le.thr )then
+    accu += psi_coef_sorted(j,1)**2
+    k += 1
+    psi_det_sorted_gen_HF_SD_order(j) = k
+    do i = 1, N_int
+     psi_det_generators_HF_SD(i,1,k) = psi_det_sorted(i,1,j)
+     psi_det_generators_HF_SD(i,2,k) = psi_det_sorted(i,2,j)
+    enddo
+    do i = 1, N_states
+     psi_coef_generators_HF_SD(k,i) = psi_coef_sorted(j,i)
+    enddo
+   else
+    inongen += 1
+    nongen(inongen) = j
+   endif
+ end do
+
+ psi_det_sorted_gen_HF_SD(:,:,:N_det_generators_HF_SD) = psi_det_generators_HF_SD(:,:,:N_det_generators_HF_SD)
+ psi_coef_sorted_gen_HF_SD(:N_det_generators_HF_SD, :) = psi_coef_generators_HF_SD(:N_det_generators_HF_SD, :)
+ do i=1,inongen
+   psi_det_sorted_gen_HF_SD_order(nongen(i)) = N_det_generators_HF_SD+i
+   psi_det_sorted_gen_HF_SD(:,:,N_det_generators_HF_SD+i) = psi_det_sorted(:,:,nongen(i))
+   psi_coef_sorted_gen_HF_SD(N_det_generators_HF_SD+i, :) = psi_coef_sorted(nongen(i),:)
+ end do
+
+END_PROVIDER
+

src/two_body_rdm/orb_range_2_rdm.irp.f

@@ -76,8 +76,12 @@
  ispin = 4 
  state_av_act_two_rdm_spin_trace_mo = 0.d0
  integer :: i
-
+ double precision :: wall_0,wall_1
+ call wall_time(wall_0)
+ print*,'providing the  state average TWO-RDM ...'
  call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
 
+ call wall_time(wall_1)
+ print*,'Time to provide the state average TWO-RDM',wall_1 - wall_0
  END_PROVIDER