Merge branch 'cleaning_dft' of https://github.com/QuantumPackage/qp2 into cleaning_dft

src/cipsi/EZFIO.cfg

@@ -1,9 +1,3 @@
-[pert_2rdm]
-type: logical
-doc: If true, computes the one- and two-body rdms with perturbation theory
-interface: ezfio,provider,ocaml
-default: False
-
 [save_wf_after_selection]
 type: logical
 doc: If true, saves the wave function after the selection, before the diagonalization

src/cipsi/NEED

@@ -2,5 +2,4 @@ perturbation
 zmq
 mpi
 iterations
-two_body_rdm
 csf

src/cipsi/pert_rdm_providers.irp.f

@@ -1,183 +0,0 @@
-
-use bitmasks
-use omp_lib
-
-BEGIN_PROVIDER [ integer(omp_lock_kind), pert_2rdm_lock]
-  use f77_zmq
-  implicit none
-  call omp_init_lock(pert_2rdm_lock)
-END_PROVIDER
-
-BEGIN_PROVIDER [integer, n_orb_pert_rdm]
- implicit none
- n_orb_pert_rdm = n_act_orb
-END_PROVIDER 
-
-BEGIN_PROVIDER [integer, list_orb_reverse_pert_rdm, (mo_num)]
- implicit none
- list_orb_reverse_pert_rdm = list_act_reverse
-
-END_PROVIDER 
-
-BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)]
- implicit none
- list_orb_pert_rdm = list_act
-
-END_PROVIDER
-
-BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)]
- implicit none
- pert_2rdm_provider = 0.d0
-
-END_PROVIDER
-
-subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, psi_det_connection, psi_coef_connection_reverse, n_det_connection)
-  use bitmasks
-  use selection_types
-  implicit none
-  
-  integer, intent(in)           :: n_det_connection
-  double precision, intent(in)  :: psi_coef_connection_reverse(N_states,n_det_connection)
-  integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection)
-  integer, intent(in) :: i_generator, sp, h1, h2
-  double precision, intent(in) :: mat(N_states, mo_num, mo_num)
-  logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num)
-  double precision, intent(in)           :: fock_diag_tmp(mo_num)
-  double precision, intent(in)    :: E0(N_states)
-  type(pt2_type), intent(inout)   :: pt2_data
-  type(selection_buffer), intent(inout) :: buf
-  logical :: ok
-  integer :: s1, s2, p1, p2, ib, j, istate, jstate
-  integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
-  double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states)
-  double precision, external :: diag_H_mat_elem_fock
-  double precision :: E_shift
-
-  logical, external :: detEq
-  double precision, allocatable :: values(:)
-  integer, allocatable          :: keys(:,:)
-  integer                       :: nkeys
-  integer :: sze_buff
-  sze_buff = 5 * mo_num ** 2 
-  allocate(keys(4,sze_buff),values(sze_buff))
-  nkeys = 0
-  if(sp == 3) then
-    s1 = 1
-    s2 = 2
-  else
-    s1 = sp
-    s2 = sp
-  end if
-  call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int)
-  E_shift = 0.d0
-
-  if (h0_type == 'CFG') then
-    j = det_to_configuration(i_generator)
-    E_shift = psi_det_Hii(i_generator) - psi_configuration_Hii(j)
-  endif
-
-  do p1=1,mo_num
-    if(bannedOrb(p1, s1)) cycle
-    ib = 1
-    if(sp /= 3) ib = p1+1
-
-    do p2=ib,mo_num
-
-! -----
-! /!\ Generating only single excited determinants doesn't work because a
-! determinant generated by a single excitation may be doubly excited wrt
-! to a determinant of the future. In that case, the determinant will be
-! detected as already generated when generating in the future with a
-! double excitation.
-!     
-!      if (.not.do_singles) then
-!        if ((h1 == p1) .or. (h2 == p2)) then
-!          cycle
-!        endif
-!      endif
-!
-!      if (.not.do_doubles) then
-!        if ((h1 /= p1).and.(h2 /= p2)) then
-!          cycle
-!        endif
-!      endif
-! -----
-
-      if(bannedOrb(p2, s2)) cycle
-      if(banned(p1,p2)) cycle
-
-
-      if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle
-      call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
-
-      if (do_only_cas) then
-        integer, external :: number_of_holes, number_of_particles
-        if (number_of_particles(det)>0) then
-          cycle
-        endif
-        if (number_of_holes(det)>0) then
-          cycle
-        endif
-      endif
-
-      if (do_ddci) then
-        logical, external  :: is_a_two_holes_two_particles
-        if (is_a_two_holes_two_particles(det)) then
-          cycle
-        endif
-      endif
-
-      if (do_only_1h1p) then
-        logical, external :: is_a_1h1p
-        if (.not.is_a_1h1p(det)) cycle
-      endif
-
-
-      Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
-
-      sum_e_pert = 0d0
-      integer :: degree
-      call get_excitation_degree(det,HF_bitmask,degree,N_int)
-      if(degree == 2)cycle
-      do istate=1,N_states
-        delta_E = E0(istate) - Hii + E_shift
-        alpha_h_psi = mat(istate, p1, p2)
-        val = alpha_h_psi + alpha_h_psi
-        tmp = dsqrt(delta_E * delta_E + val * val)
-        if (delta_E < 0.d0) then
-            tmp = -tmp
-        endif
-        e_pert = 0.5d0 * (tmp - delta_E)
-        coef(istate) = e_pert / alpha_h_psi
-        print*,e_pert,coef,alpha_h_psi
-        pt2_data % pt2(istate) += e_pert
-        pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi
-      enddo
-
-      do istate=1,N_states
-        alpha_h_psi = mat(istate, p1, p2)
-        e_pert = coef(istate) * alpha_h_psi
-        do jstate=1,N_states
-          pt2_data % overlap(jstate,jstate) = coef(istate) * coef(jstate)
-        enddo
-
-        if (weight_selection /= 5) then
-          ! Energy selection
-          sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
-
-        else
-          ! Variance selection
-          sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
-        endif
-      end do
-      call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff)
-
-      if(sum_e_pert <= buf%mini) then
-        call add_to_selection_buffer(buf, det, sum_e_pert)
-      end if
-    end do
-  end do
-  call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
-end
-
-

src/cipsi/pt2_stoch_routines.irp.f

@@ -133,7 +133,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
   PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted
   PROVIDE psi_det_hii selection_weight pseudo_sym
   PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
-  PROVIDE pert_2rdm excitation_beta_max  excitation_alpha_max excitation_max
+  PROVIDE excitation_beta_max  excitation_alpha_max excitation_max
 
   if (h0_type == 'CFG') then
     PROVIDE psi_configuration_hii det_to_configuration

src/cipsi/selection.irp.f

@@ -464,14 +464,14 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
       allocate (fullminilist (N_int, 2, fullinteresting(0)), &
                     minilist (N_int, 2,     interesting(0)) )
-      if(pert_2rdm)then
+!      if(pert_2rdm)then
-        allocate(coef_fullminilist_rev(N_states,fullinteresting(0)))
+!        allocate(coef_fullminilist_rev(N_states,fullinteresting(0)))
-        do i=1,fullinteresting(0)
+!        do i=1,fullinteresting(0)
-          do j = 1, N_states
+!          do j = 1, N_states
-            coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j)
+!            coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j)
-          enddo
+!          enddo
-        enddo
+!        enddo
-      endif
+!      endif
 
       do i=1,fullinteresting(0)
         do k=1,N_int
@@ -531,19 +531,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
             call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
 
-            if(.not.pert_2rdm)then
+!            if(.not.pert_2rdm)then
              call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
-            else
+!            else
-             call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
+!             call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
-            endif
+!            endif
           end if
         enddo
         if(s1 /= s2) monoBdo = .false.
       enddo
       deallocate(fullminilist,minilist)
-      if(pert_2rdm)then
+!      if(pert_2rdm)then
-        deallocate(coef_fullminilist_rev)
+!        deallocate(coef_fullminilist_rev)
-      endif
+!      endif
     enddo
   enddo
   deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)

src/cipsi/update_2rdm.irp.f

@@ -1,223 +0,0 @@
-use bitmasks
-
-subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff)
- implicit none
- integer, intent(in) :: n_det_connection,sze_buff
- double precision, intent(in) :: coef(N_states)
- integer(bit_kind), intent(in) :: det(N_int,2)
- integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) 
- double precision, intent(in)  :: psi_coef_connection_reverse(N_states,n_det_connection)
- integer,           intent(inout) :: keys(4,sze_buff),nkeys
- double precision,  intent(inout) :: values(sze_buff)
- integer :: i,j
- integer                        :: exc(0:2,2,2)
- integer                        :: degree
- double precision               :: phase, contrib
- do i = 1, n_det_connection
-  call get_excitation(det,psi_det_connection(1,1,i),exc,degree,phase,N_int) 
-  if(degree.gt.2)cycle
-  contrib = 0.d0
-  do j = 1, N_states
-   contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i) * phase * coef(j)
-  enddo
-  ! case of single excitations 
-  if(degree == 1)then
-   if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then
-    call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
-    nkeys = 0
-   endif
-   call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff)
-  else 
- !! case of double excitations 
- ! if (nkeys + 4 .ge. sze_buff)then
- !  call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
- !  nkeys = 0
- ! endif
- ! call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
-  endif
- enddo
-!call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
-!nkeys = 0
-
-end
-
-subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff)
- implicit none
- integer, intent(in) :: sze_buff
- integer(bit_kind), intent(in) :: det1(N_int,2)
- integer(bit_kind), intent(in) :: det2(N_int,2)
- integer,intent(in)             :: exc(0:2,2,2)
- double precision,intent(in)    :: phase, contrib
- integer, intent(inout)         :: nkeys, keys(4,sze_buff)
- double precision, intent(inout):: values(sze_buff)
- 
- integer                        :: occ(N_int*bit_kind_size,2)
- integer                        :: n_occ_ab(2),ispin,other_spin
- integer :: h1,h2,p1,p2,i
- call bitstring_to_list_ab(det1, occ, n_occ_ab, N_int)
-  
- if (exc(0,1,1) == 1) then
-  ! Mono alpha
-  h1 = exc(1,1,1)
-  p1 = exc(1,2,1)
-  ispin = 1
-  other_spin = 2
- else
-  ! Mono beta
-  h1 = exc(1,1,2)
-  p1 = exc(1,2,2)
-  ispin = 2
-  other_spin = 1
- endif
- if(list_orb_reverse_pert_rdm(h1).lt.0)return
- h1 = list_orb_reverse_pert_rdm(h1)
- if(list_orb_reverse_pert_rdm(p1).lt.0)return
- p1 = list_orb_reverse_pert_rdm(p1)
- !update the alpha/beta part
- do i = 1, n_occ_ab(other_spin)
-  h2 = occ(i,other_spin)
-  if(list_orb_reverse_pert_rdm(h2).lt.0)return
-  h2 = list_orb_reverse_pert_rdm(h2)
-
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h1
-  keys(2,nkeys) = h2
-  keys(3,nkeys) = p1
-  keys(4,nkeys) = h2
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h2
-  keys(2,nkeys) = h1
-  keys(3,nkeys) = h2
-  keys(4,nkeys) = p1
- enddo 
- !update the same spin part 
-!do i = 1, n_occ_ab(ispin)
-! h2 = occ(i,ispin)
-! if(list_orb_reverse_pert_rdm(h2).lt.0)return
-! h2 = list_orb_reverse_pert_rdm(h2)
-
-! nkeys += 1
-! values(nkeys) = 0.5d0 * contrib * phase
-! keys(1,nkeys) = h1
-! keys(2,nkeys) = h2
-! keys(3,nkeys) = p1
-! keys(4,nkeys) = h2
-
-! nkeys += 1
-! values(nkeys) = - 0.5d0 * contrib * phase
-! keys(1,nkeys) = h1
-! keys(2,nkeys) = h2
-! keys(3,nkeys) = h2
-! keys(4,nkeys) = p1
-!
-! nkeys += 1
-! values(nkeys) = 0.5d0 * contrib * phase
-! keys(1,nkeys) = h2
-! keys(2,nkeys) = h1
-! keys(3,nkeys) = h2
-! keys(4,nkeys) = p1
-
-! nkeys += 1
-! values(nkeys) = - 0.5d0 * contrib * phase
-! keys(1,nkeys) = h2
-! keys(2,nkeys) = h1
-! keys(3,nkeys) = p1
-! keys(4,nkeys) = h2
-!enddo 
-
-end
-
-subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
- implicit none
- integer, intent(in) :: sze_buff
- integer,intent(in)             :: exc(0:2,2,2)
- double precision,intent(in)    :: phase, contrib
- integer, intent(inout)         :: nkeys, keys(4,sze_buff)
- double precision, intent(inout):: values(sze_buff)
- integer :: h1,h2,p1,p2
-
- if (exc(0,1,1) == 1) then
-  ! Double alpha/beta
-  h1 = exc(1,1,1) 
-  h2 = exc(1,1,2) 
-  p1 = exc(1,2,1)
-  p2 = exc(1,2,2)
-  ! check if the orbitals involved are within the orbital range
-  if(list_orb_reverse_pert_rdm(h1).lt.0)return
-  h1 = list_orb_reverse_pert_rdm(h1)
-  if(list_orb_reverse_pert_rdm(h2).lt.0)return
-  h2 = list_orb_reverse_pert_rdm(h2)
-  if(list_orb_reverse_pert_rdm(p1).lt.0)return
-  p1 = list_orb_reverse_pert_rdm(p1)
-  if(list_orb_reverse_pert_rdm(p2).lt.0)return
-  p2 = list_orb_reverse_pert_rdm(p2)
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h1
-  keys(2,nkeys) = h2
-  keys(3,nkeys) = p1
-  keys(4,nkeys) = p2
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = p1
-  keys(2,nkeys) = p2
-  keys(3,nkeys) = h1
-  keys(4,nkeys) = h2 
-
- else 
-  if (exc(0,1,1) == 2) then
-    ! Double alpha/alpha
-   h1 =  exc(1,1,1)
-   h2 =  exc(2,1,1)
-   p1 =  exc(1,2,1)
-   p2 =  exc(2,2,1)
-  else if (exc(0,1,2) == 2) then
-   ! Double beta
-   h1 = exc(1,1,2)
-   h2 = exc(2,1,2)
-   p1 = exc(1,2,2)
-   p2 = exc(2,2,2)
-  endif
-  ! check if the orbitals involved are within the orbital range
-  if(list_orb_reverse_pert_rdm(h1).lt.0)return
-  h1 = list_orb_reverse_pert_rdm(h1)
-  if(list_orb_reverse_pert_rdm(h2).lt.0)return
-  h2 = list_orb_reverse_pert_rdm(h2)
-  if(list_orb_reverse_pert_rdm(p1).lt.0)return
-  p1 = list_orb_reverse_pert_rdm(p1)
-  if(list_orb_reverse_pert_rdm(p2).lt.0)return
-  p2 = list_orb_reverse_pert_rdm(p2)
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h1
-  keys(2,nkeys) = h2
-  keys(3,nkeys) = p1
-  keys(4,nkeys) = p2
- 
-  nkeys += 1
-  values(nkeys) = - 0.5d0 * contrib * phase
-  keys(1,nkeys) = h1
-  keys(2,nkeys) = h2
-  keys(3,nkeys) = p2
-  keys(4,nkeys) = p1
-                                        
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h2
-  keys(2,nkeys) = h1
-  keys(3,nkeys) = p2
-  keys(4,nkeys) = p1
- 
-  nkeys += 1
-  values(nkeys) = - 0.5d0 * contrib * phase
-  keys(1,nkeys) = h2
-  keys(2,nkeys) = h1
-  keys(3,nkeys) = p1
-  keys(4,nkeys) = p2
- endif
- 
-end
-
-

src/cipsi/zmq_selection.irp.f

@@ -22,7 +22,7 @@ subroutine ZMQ_selection(N_in, pt2_data)
     PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
     PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym
     PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
-    PROVIDE pert_2rdm excitation_beta_max  excitation_alpha_max excitation_max
+    PROVIDE excitation_beta_max  excitation_alpha_max excitation_max
 
     call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')
 

src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f

@@ -0,0 +1,481 @@
+
+subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc)
+  use mmap_module
+  implicit none
+  BEGIN_DOC
+  ! Generic Davidson diagonalization with ONE DIAGONAL DRESSING OPERATOR 
+  !
+  ! Dress_jj : DIAGONAL DRESSING of the Hamiltonian 
+  !
+  ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
+  !
+  ! u_in : guess coefficients on the various states. Overwritten on exit
+  !
+  ! sze : leftmost dimension of u_in
+  !
+  ! sze : Number of determinants
+  !
+  ! N_st : Number of eigenstates
+  !
+  ! N_st_diag_in : Number of states in which H is diagonalized. Assumed > sze
+  !
+  ! Initial guess vectors are not necessarily orthonormal
+  !
+  ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output)
+  END_DOC
+  integer, intent(in)             :: sze, N_st, N_st_diag_in
+  double precision,  intent(in)   :: H_jj(sze),Dress_jj(sze)
+  double precision, intent(inout) :: u_in(sze,N_st_diag_in)
+  double precision, intent(out)   :: energies(N_st)
+  external hcalc
+
+  integer                        :: iter, N_st_diag
+  integer                        :: i,j,k,l,m
+  logical, intent(inout)         :: converged
+
+  double precision, external     :: u_dot_v, u_dot_u
+
+  integer                        :: k_pairs, kl
+
+  integer                        :: iter2, itertot
+  double precision, allocatable  :: y(:,:), h(:,:), lambda(:)
+  double precision, allocatable  :: residual_norm(:)
+  character*(16384)              :: write_buffer
+  double precision               :: to_print(2,N_st)
+  double precision               :: cpu, wall
+  integer                        :: shift, shift2, itermax, istate
+  double precision               :: r1, r2, alpha
+  integer                        :: nproc_target
+  integer                        :: order(N_st_diag_in)
+  double precision               :: cmax
+  double precision, allocatable  :: U(:,:), overlap(:,:)!, S_d(:,:)
+  double precision, pointer      :: W(:,:)
+  logical                        :: disk_based
+  double precision               :: energy_shift(N_st_diag_in*davidson_sze_max)
+
+  include 'constants.include.F'
+
+  N_st_diag = N_st_diag_in
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda
+  if (N_st_diag*3 > sze) then
+    print *,  'error in Davidson :'
+    print *,  'Increase n_det_max_full to ', N_st_diag*3
+    stop -1
+  endif
+
+  itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1
+  itertot = 0
+
+  if (state_following) then
+    allocate(overlap(N_st_diag*itermax, N_st_diag*itermax))
+  else
+    allocate(overlap(1,1))  ! avoid 'if' for deallocate
+  endif
+  overlap = 0.d0
+
+  provide threshold_davidson !nthreads_davidson
+  call write_time(6)
+  write(6,'(A)') ''
+  write(6,'(A)') 'Davidson Diagonalization'
+  write(6,'(A)') '------------------------'
+  write(6,'(A)') ''
+
+  ! Find max number of cores to fit in memory
+  ! -----------------------------------------
+
+  nproc_target = nproc
+  double precision :: rss
+  integer :: maxab
+  maxab = sze 
+
+  m=1
+  disk_based = .False.
+  call resident_memory(rss)
+  do
+    r1 = 8.d0 *                                   &! bytes
+         ( dble(sze)*(N_st_diag*itermax)          &! U
+         + 1.d0*dble(sze*m)*(N_st_diag*itermax)  &! W
+         + 2.0d0*(N_st_diag*itermax)**2           &! h,y
+         + 2.d0*(N_st_diag*itermax)               &! s2,lambda
+         + 1.d0*(N_st_diag)                       &! residual_norm
+                                                   ! In H_S2_u_0_nstates_zmq
+         + 3.d0*(N_st_diag*N_det)                 &! u_t, v_t, s_t on collector
+         + 3.d0*(N_st_diag*N_det)                 &! u_t, v_t, s_t on slave
+         + 0.5d0*maxab                            &! idx0 in H_S2_u_0_nstates_openmp_work_*
+         + nproc_target *                         &! In OMP section
+           ( 1.d0*(N_int*maxab)                   &! buffer
+           + 3.5d0*(maxab) )                      &! singles_a, singles_b, doubles, idx
+         ) / 1024.d0**3
+
+    if (nproc_target == 0) then
+      call check_mem(r1,irp_here)
+      nproc_target = 1
+      exit
+    endif
+
+    if (r1+rss < qp_max_mem) then
+      exit
+    endif
+
+    if (itermax > 4) then
+      itermax = itermax - 1
+    else if (m==1.and.disk_based_davidson) then
+      m=0
+      disk_based = .True.
+      itermax = 6
+    else
+      nproc_target = nproc_target - 1
+    endif
+
+  enddo
+  nthreads_davidson = nproc_target
+  TOUCH nthreads_davidson
+  call write_int(6,N_st,'Number of states')
+  call write_int(6,N_st_diag,'Number of states in diagonalization')
+  call write_int(6,sze,'Number of basis functions')
+  call write_int(6,nproc_target,'Number of threads for diagonalization')
+  call write_double(6, r1, 'Memory(Gb)')
+  if (disk_based) then
+    print *, 'Using swap space to reduce RAM'
+  endif
+
+  double precision, allocatable   :: H_jj_tmp(:)
+  ASSERT (N_st > 0)
+  ASSERT (sze > 0)
+  allocate(H_jj_tmp(sze))
+
+  do i=1,sze
+    H_jj_tmp(i)  = H_jj(i) + Dress_jj(i)
+  enddo
+
+  !---------------
+
+  write(6,'(A)') ''
+  write_buffer = '====='
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//' ================  ==========='
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+  write_buffer = 'Iter'
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//'       Energy         Residual '
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+  write_buffer = '====='
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//' ================  ==========='
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+
+
+  allocate(W(sze,N_st_diag*itermax))
+
+  allocate(                                                          &
+      ! Large
+      U(sze,N_st_diag*itermax),                                      &
+      ! Small
+      h(N_st_diag*itermax,N_st_diag*itermax),                        &
+      y(N_st_diag*itermax,N_st_diag*itermax),                        &
+      residual_norm(N_st_diag),                                      &
+      lambda(N_st_diag*itermax))
+
+  h = 0.d0
+  U = 0.d0
+  y = 0.d0
+
+
+  ASSERT (N_st > 0)
+  ASSERT (N_st_diag >= N_st)
+  ASSERT (sze > 0)
+
+  ! Davidson iterations
+  ! ===================
+
+  converged = .False.
+
+  ! Initialize from N_st to N_st_diat with gaussian random numbers
+  ! to be sure to have overlap with any eigenvectors
+  do k=N_st+1,N_st_diag
+    u_in(k,k) = 10.d0
+    do i=1,sze
+      call random_number(r1)
+      call random_number(r2)
+      r1 = dsqrt(-2.d0*dlog(r1))
+      r2 = dtwo_pi*r2
+      u_in(i,k) = r1*dcos(r2)
+    enddo
+  enddo
+  ! Normalize all states 
+  do k=1,N_st_diag
+    call normalize(u_in(1,k),sze)
+  enddo
+
+  ! Copy from the guess input "u_in" to the working vectors "U"
+  do k=1,N_st_diag
+    do i=1,sze
+      U(i,k) = u_in(i,k)
+    enddo
+  enddo
+
+
+  do while (.not.converged)
+    itertot = itertot+1
+    if (itertot == 8) then
+      exit
+    endif
+
+    do iter=1,itermax-1
+
+      shift  = N_st_diag*(iter-1)
+      shift2 = N_st_diag*iter
+
+      if ((iter > 1).or.(itertot == 1)) then
+        ! Compute |W_k> = \sum_i |i><i|H|u_k>
+        ! -----------------------------------
+
+         ! Gram-Schmidt to orthogonalize all new guess with the previous vectors 
+          call ortho_qr(U,size(U,1),sze,shift2)
+          call ortho_qr(U,size(U,1),sze,shift2)
+          !    it does W = H U with W(sze,N_st_diag),U(sze,N_st_diag)
+          !    where sze is the size of the vector, N_st_diag is the number of states 
+          call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag,sze)
+          ! Compute then the DIAGONAL PART OF THE DRESSING 
+          ! <i|W_k> += Dress_jj(i) * <i|U>
+          call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze)
+      else
+         ! Already computed in update below
+         continue
+      endif
+
+      ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
+      ! -------------------------------------------
+
+      call dgemm('T','N', shift2, shift2, sze,                       &
+          1.d0, U, size(U,1), W, size(W,1),                          &
+          0.d0, h, size(h,1))
+
+      ! Diagonalize h y = lambda y
+      ! ---------------
+
+      call lapack_diag(lambda,y,h,size(h,1),shift2)
+
+      if (state_following) then
+
+        overlap = -1.d0
+        do k=1,shift2
+          do i=1,shift2
+            overlap(k,i) = dabs(y(k,i))
+          enddo
+        enddo
+        do k=1,N_st
+          cmax = -1.d0
+          do i=1,N_st
+            if (overlap(i,k) > cmax) then
+              cmax = overlap(i,k)
+              order(k) = i
+            endif
+          enddo
+          do i=1,N_st_diag
+            overlap(order(k),i) = -1.d0
+          enddo
+        enddo
+        overlap = y
+        do k=1,N_st
+          l = order(k)
+          if (k /= l) then
+            y(1:shift2,k) = overlap(1:shift2,l)
+          endif
+        enddo
+        do k=1,N_st
+          overlap(k,1) = lambda(k)
+        enddo
+        do k=1,N_st
+          l = order(k)
+          if (k /= l) then
+            lambda(k) = overlap(l,1)
+          endif
+        enddo
+
+      endif
+
+
+      ! Express eigenvectors of h in the determinant basis
+      ! --------------------------------------------------
+
+      call dgemm('N','N', sze, N_st_diag, shift2,                    &
+          1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1))
+      call dgemm('N','N', sze, N_st_diag, shift2,                    &
+          1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1))
+
+      ! Compute residual vector and davidson step
+      ! -----------------------------------------
+
+      !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
+      do k=1,N_st_diag
+        do i=1,sze
+          U(i,shift2+k) =  &
+            (lambda(k) * U(i,shift2+k) - W(i,shift2+k) )      &
+              /max(H_jj_tmp(i) - lambda (k),1.d-2)
+        enddo
+
+        if (k <= N_st) then
+          residual_norm(k) = u_dot_u(U(1,shift2+k),sze)
+          to_print(1,k) = lambda(k) 
+          to_print(2,k) = residual_norm(k)
+        endif
+      enddo
+      !$OMP END PARALLEL DO
+
+
+      if ((itertot>1).and.(iter == 1)) then
+        !don't print 
+        continue
+      else
+        write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter-1, to_print(1:2,1:N_st)
+      endif
+
+      ! Check convergence
+      if (iter > 1) then
+          converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson
+      endif   
+      
+
+      do k=1,N_st
+        if (residual_norm(k) > 1.e8) then
+          print *, 'Davidson failed'
+          stop -1
+        endif
+      enddo
+      if (converged) then
+        exit
+      endif
+
+      logical, external :: qp_stop
+      if (qp_stop()) then
+        converged = .True.
+        exit
+      endif
+
+
+    enddo
+
+    call dgemm('N','N', sze, N_st_diag, shift2, 1.d0,      &
+        W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
+    do k=1,N_st_diag
+      do i=1,sze
+        W(i,k) = u_in(i,k)
+      enddo
+    enddo
+
+    call dgemm('N','N', sze, N_st_diag, shift2, 1.d0,      &
+        U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
+    do k=1,N_st_diag
+      do i=1,sze
+        U(i,k) = u_in(i,k)
+      enddo
+    enddo
+      call ortho_qr(U,size(U,1),sze,N_st_diag)
+      call ortho_qr(U,size(U,1),sze,N_st_diag)
+    do j=1,N_st_diag
+      k=1
+      do while ((k<sze).and.(U(k,j) == 0.d0))
+        k = k+1
+      enddo
+      if (U(k,j) * u_in(k,j) < 0.d0) then
+        do i=1,sze
+          W(i,j) = -W(i,j)
+        enddo
+      endif
+    enddo
+  enddo
+
+  do k=1,N_st
+    energies(k) = lambda(k)
+  enddo
+  write_buffer = '====='
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//' ================  ==========='
+  enddo
+  write(6,'(A)') trim(write_buffer)
+  write(6,'(A)') ''
+  call write_time(6)
+
+  deallocate(W)
+
+  deallocate (                                                       &
+      residual_norm,                                                 &
+      U, h,                                                 &
+      y,                                                             &
+      lambda                                                         &
+      )
+  deallocate(overlap)
+  FREE nthreads_davidson
+end
+
+subroutine hcalc_template(v,u,N_st,sze)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Template of routine for the application of H
+  !
+  ! Here, it is done with the Hamiltonian matrix 
+  !
+  ! on the set of determinants of psi_det 
+  !
+  ! Computes $v = H | u \rangle$ 
+  !
+  END_DOC
+  integer, intent(in)              :: N_st,sze
+  double precision, intent(in)     :: u(sze,N_st)
+  double precision, intent(inout)  :: v(sze,N_st)
+  integer :: i,j,istate
+  v = 0.d0
+  do istate = 1, N_st
+   do i = 1, sze
+    do j = 1, sze
+      v(i,istate) += H_matrix_all_dets(j,i) * u(j,istate)
+    enddo
+   enddo
+   do i = 1, sze
+    v(i,istate) += u(i,istate) * nuclear_repulsion
+   enddo
+  enddo
+end
+
+subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze)
+  implicit none
+  BEGIN_DOC
+  ! Routine that computes the diagonal part of the dressing 
+  !
+  ! v(i) += u(i) * dress_diag(i)
+  !
+  ! !!!!!!!! WARNING !!!!!!!! the vector v is not initialized 
+  !
+  ! !!!!!!!! SO MAKE SURE THERE ARE SOME MEANINGFUL VALUES IN THERE
+  END_DOC
+  integer, intent(in)              :: N_st,sze
+  double precision, intent(in)     :: u(sze,N_st),dress_diag(sze)
+  double precision, intent(inout)  :: v(sze,N_st)
+  integer :: i,istate
+  do istate = 1, N_st
+   do i = 1, sze
+     v(i,istate) += dress_diag(i) * u(i,istate)
+   enddo
+  enddo
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

src/dav_general_mat/dav_double_dress_ext_rout.irp.f

@@ -0,0 +1,518 @@
+subroutine dav_double_dressed(u_in,H_jj,Dress_jj,Dressing_vec,idx_dress,energies,sze,N_st,N_st_diag,converged,hcalc)
+  use mmap_module
+  BEGIN_DOC
+  ! Generic Davidson diagonalization with TWO DRESSING VECTORS 
+  !
+  ! Dress_jj : DIAGONAL DRESSING of the Hamiltonian 
+  !
+  ! Dressing_vec : COLUMN / LINE DRESSING VECTOR 
+  !
+  ! idx_dress : position of the basis function used to use the Dressing_vec (usually the largest coeff)
+  !
+  ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
+  !
+  ! u_in : guess coefficients on the various states. Overwritten on exit
+  !
+  ! sze : leftmost dimension of u_in
+  !
+  ! sze : Number of determinants
+  !
+  ! N_st : Number of eigenstates
+  !
+  ! N_st_diag : Number of states in which H is diagonalized. Assumed > sze
+  !
+  ! Initial guess vectors are not necessarily orthonormal
+  !
+  ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output)
+  END_DOC
+  implicit none
+  integer, intent(in)             :: sze, N_st, N_st_diag, idx_dress
+  double precision,  intent(in)   :: H_jj(sze),Dress_jj(sze),Dressing_vec(sze,N_st)
+  double precision, intent(inout) :: u_in(sze,N_st_diag)
+  double precision, intent(out)   :: energies(N_st_diag)
+  logical, intent(out)            :: converged
+  external hcalc
+
+  double precision, allocatable   :: H_jj_tmp(:)
+  ASSERT (N_st > 0)
+  ASSERT (sze > 0)
+  allocate(H_jj_tmp(sze))
+
+  do i=1,sze
+    H_jj_tmp(i)  = H_jj(i) + Dress_jj(i)
+  enddo
+  do k=1,N_st
+    do i=1,sze
+      H_jj_tmp(i)  += u_in(i,k) * Dressing_vec(i,k)
+    enddo
+  enddo
+
+  integer                        :: iter
+  integer                        :: i,j,k,l,m
+
+  double precision, external     :: u_dot_v, u_dot_u
+
+  integer                        :: k_pairs, kl
+
+  integer                        :: iter2, itertot
+  double precision, allocatable  :: y(:,:), h(:,:), lambda(:)
+  double precision, allocatable  :: s_tmp(:,:)
+  double precision, allocatable  :: residual_norm(:),inv_c_idx_dress_vec(:)
+  character*(16384)              :: write_buffer
+  double precision               :: to_print(2,N_st),inv_c_idx_dress
+  double precision               :: cpu, wall
+  integer                        :: shift, shift2, itermax, istate
+  double precision               :: r1, r2, alpha
+  logical                        :: state_ok(N_st_diag*davidson_sze_max)
+  integer                        :: nproc_target
+  integer                        :: order(N_st_diag)
+  double precision               :: cmax
+  double precision, allocatable  :: U(:,:), overlap(:,:)
+  double precision, pointer      :: W(:,:)
+  logical                        :: disk_based
+  double precision               :: energy_shift(N_st_diag*davidson_sze_max)
+
+
+  allocate(inv_c_idx_dress_vec(N_st))
+  inv_c_idx_dress = 1.d0/u_in(idx_dress,1)
+  do i = 1, N_st
+   inv_c_idx_dress_vec(i) = 1.d0/u_in(idx_dress,i)
+  enddo
+  include 'constants.include.F'
+
+  integer :: N_st_diag_in 
+  N_st_diag_in = N_st_diag
+
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda
+  if (N_st_diag_in*3 > sze) then
+    print *,  'error in Davidson :'
+    print *,  'Increase n_det_max_full to ', N_st_diag_in*3
+    stop -1
+  endif
+
+  itermax = max(2,min(davidson_sze_max, sze/N_st_diag_in))+1
+  itertot = 0
+
+  if (state_following) then
+    allocate(overlap(N_st_diag_in*itermax, N_st_diag_in*itermax))
+  else
+    allocate(overlap(1,1))  ! avoid 'if' for deallocate
+  endif
+  overlap = 0.d0
+
+  call write_time(6)
+  write(6,'(A)') ''
+  write(6,'(A)') 'Davidson Diagonalization'
+  write(6,'(A)') '------------------------'
+  write(6,'(A)') ''
+
+  ! Find max number of cores to fit in memory
+  ! -----------------------------------------
+
+  nproc_target = nproc
+  double precision :: rss
+  integer :: maxab
+  maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
+
+  m=1
+  disk_based = .False.
+  call resident_memory(rss)
+  do
+    r1 = 8.d0 *                                   &! bytes
+         ( dble(sze)*(N_st_diag_in*itermax)          &! U
+         + 1.0d0*dble(sze*m)*(N_st_diag_in*itermax)  &! W
+         + 3.0d0*(N_st_diag_in*itermax)**2           &! h,y,s_tmp
+         + 1.d0*(N_st_diag_in*itermax)               &! lambda
+         + 1.d0*(N_st_diag_in)                       &! residual_norm
+                                                   ! In H_u_0_nstates_zmq
+         + 2.d0*(N_st_diag_in*N_det)                 &! u_t, v_t, on collector
+         + 2.d0*(N_st_diag_in*N_det)                 &! u_t, v_t, on slave
+         + 0.5d0*maxab                            &! idx0 in H_u_0_nstates_openmp_work_*
+         + nproc_target *                         &! In OMP section
+           ( 1.d0*(N_int*maxab)                   &! buffer
+           + 3.5d0*(maxab) )                      &! singles_a, singles_b, doubles, idx
+         ) / 1024.d0**3
+
+    if (nproc_target == 0) then
+      call check_mem(r1,irp_here)
+      nproc_target = 1
+      exit
+    endif
+
+    if (r1+rss < qp_max_mem) then
+      exit
+    endif
+
+    if (itermax > 4) then
+      itermax = itermax - 1
+    else if (m==1.and.disk_based_davidson) then
+      m=0
+      disk_based = .True.
+      itermax = 6
+    else
+      nproc_target = nproc_target - 1
+    endif
+
+  enddo
+  nthreads_davidson = nproc_target
+  TOUCH nthreads_davidson
+  call write_int(6,N_st,'Number of states')
+  call write_int(6,N_st_diag_in,'Number of states in diagonalization')
+  call write_int(6,sze,'Number of basis functions ')
+  call write_int(6,nproc_target,'Number of threads for diagonalization')
+  call write_double(6, r1, 'Memory(Gb)')
+  if (disk_based) then
+    print *, 'Using swap space to reduce RAM'
+  endif
+
+  !---------------
+
+  write(6,'(A)') ''
+  write_buffer = '====='
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//' ================ ==========='
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+  write_buffer = 'Iter'
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//'       Energy        Residual '
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+  write_buffer = '====='
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//' ================ ==========='
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+
+
+  allocate(W(sze,N_st_diag_in*itermax))
+
+  allocate(                                                          &
+      ! Large
+      U(sze,N_st_diag_in*itermax),                                      &
+
+      ! Small
+      h(N_st_diag_in*itermax,N_st_diag_in*itermax),                        &
+      y(N_st_diag_in*itermax,N_st_diag_in*itermax),                        &
+      s_tmp(N_st_diag_in*itermax,N_st_diag_in*itermax),                    &
+      residual_norm(N_st_diag_in),                                      &
+      lambda(N_st_diag_in*itermax))
+
+  h = 0.d0
+  U = 0.d0
+  y = 0.d0
+  s_tmp = 0.d0
+
+
+  ASSERT (N_st > 0)
+  ASSERT (N_st_diag_in >= N_st)
+  ASSERT (sze > 0)
+
+  ! Davidson iterations
+  ! ===================
+
+  converged = .False.
+
+  do k=N_st+1,N_st_diag_in
+    do i=1,sze
+        call random_number(r1)
+        call random_number(r2)
+        r1 = dsqrt(-2.d0*dlog(r1))
+        r2 = dtwo_pi*r2
+        u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st)
+    enddo
+    u_in(k,k) = u_in(k,k) + 10.d0
+  enddo
+  do k=1,N_st_diag_in
+    call normalize(u_in(1,k),sze)
+  enddo
+
+  do k=1,N_st_diag_in
+    do i=1,sze
+      U(i,k) = u_in(i,k)
+    enddo
+  enddo
+
+
+  do while (.not.converged)
+    itertot = itertot+1
+    if (itertot == 2) then
+      exit
+    endif
+
+    do iter=1,itermax-1
+
+      shift  = N_st_diag_in*(iter-1)
+      shift2 = N_st_diag_in*iter
+
+      if ((iter > 1).or.(itertot == 1)) then
+        ! Compute |W_k> = \sum_i |i><i|H|u_k>
+        ! -----------------------------------
+        call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag_in,sze)
+        ! Compute then the DIAGONAL PART OF THE DRESSING 
+        ! <i|W_k> += Dress_jj(i) * <i|U>
+        call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze)
+      else
+         ! Already computed in update below
+         continue
+      endif
+
+
+        if (N_st == 1) then
+
+          l = idx_dress
+          double precision :: f
+          f = inv_c_idx_dress
+          do istate=1,N_st_diag_in
+            do i=1,sze
+              W(i,shift+istate) += Dressing_vec(i,1) *f * U(l,shift+istate)
+              W(l,shift+istate) += Dressing_vec(i,1) *f * U(i,shift+istate)
+            enddo
+          enddo
+
+        else
+         print*,'dav_double_dressed routine not yet implemented for N_st > 1'
+!
+!          call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, &
+!            psi_coef, size(psi_coef,1), &
+!            U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+!
+!          call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, &
+!            Dressing_vec, size(Dressing_vec,1), s_tmp, size(s_tmp,1), &
+!            1.d0, W(1,shift+1), size(W,1))
+!
+!
+!          call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, &
+!            Dressing_vec, size(Dressing_vec,1), &
+!            U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+!
+!          call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, &
+!            psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
+!            1.d0, W(1,shift+1), size(W,1))
+!
+        endif
+
+      ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
+      ! -------------------------------------------
+
+      call dgemm('T','N', shift2, shift2, sze,                       &
+          1.d0, U, size(U,1), W, size(W,1),                          &
+          0.d0, h, size(h,1))
+      call dgemm('T','N', shift2, shift2, sze,                       &
+          1.d0, U, size(U,1), U, size(U,1),                          &
+          0.d0, s_tmp, size(s_tmp,1))
+
+      ! Diagonalize h
+      ! ---------------
+
+       integer :: lwork, info
+       double precision, allocatable :: work(:)
+
+       y = h
+       lwork = -1
+       allocate(work(1))
+       call dsygv(1,'V','U',shift2,y,size(y,1), &
+          s_tmp,size(s_tmp,1), lambda, work,lwork,info)
+       lwork = int(work(1))
+       deallocate(work)
+       allocate(work(lwork))
+       call dsygv(1,'V','U',shift2,y,size(y,1), &
+          s_tmp,size(s_tmp,1), lambda, work,lwork,info)
+       deallocate(work)
+       if (info /= 0) then
+         stop 'DSYGV Diagonalization failed'
+       endif
+
+      ! Compute Energy for each eigenvector
+      ! -----------------------------------
+
+      call dgemm('N','N',shift2,shift2,shift2,                       &
+          1.d0, h, size(h,1), y, size(y,1),                          &
+          0.d0, s_tmp, size(s_tmp,1))
+
+      call dgemm('T','N',shift2,shift2,shift2,                       &
+          1.d0, y, size(y,1), s_tmp, size(s_tmp,1),                  &
+          0.d0, h, size(h,1))
+
+      do k=1,shift2
+        lambda(k) = h(k,k)
+      enddo
+
+      if (state_following) then
+
+        overlap = -1.d0
+        do k=1,shift2
+          do i=1,shift2
+            overlap(k,i) = dabs(y(k,i))
+          enddo
+        enddo
+        do k=1,N_st
+          cmax = -1.d0
+          do i=1,N_st
+            if (overlap(i,k) > cmax) then
+              cmax = overlap(i,k)
+              order(k) = i
+            endif
+          enddo
+          do i=1,N_st_diag_in
+            overlap(order(k),i) = -1.d0
+          enddo
+        enddo
+        overlap = y
+        do k=1,N_st
+          l = order(k)
+          if (k /= l) then
+            y(1:shift2,k) = overlap(1:shift2,l)
+          endif
+        enddo
+        do k=1,N_st
+          overlap(k,1) = lambda(k)
+        enddo
+
+      endif
+
+
+      ! Express eigenvectors of h in the determinant basis
+      ! --------------------------------------------------
+
+      call dgemm('N','N', sze, N_st_diag_in, shift2,                    &
+          1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1))
+      call dgemm('N','N', sze, N_st_diag_in, shift2,                    &
+          1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1))
+
+      ! Compute residual vector and davidson step
+      ! -----------------------------------------
+
+      !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
+      do k=1,N_st_diag_in
+        do i=1,sze
+          U(i,shift2+k) =  &
+            (lambda(k) * U(i,shift2+k) - W(i,shift2+k) )      &
+              /max(H_jj_tmp(i) - lambda (k),1.d-2)
+        enddo
+
+        if (k <= N_st) then
+          residual_norm(k) = u_dot_u(U(1,shift2+k),sze)
+          to_print(1,k) = lambda(k) 
+          to_print(2,k) = residual_norm(k)
+        endif
+      enddo
+      !$OMP END PARALLEL DO
+
+
+      if ((itertot>1).and.(iter == 1)) then
+        !don't print
+        continue
+      else
+        write(*,'(1X,I3,1X,100(1X,F16.10,1X,E11.3))') iter-1, to_print(1:2,1:N_st)
+      endif
+
+      ! Check convergence
+      if (iter > 1) then
+        converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson
+      endif
+
+      do k=1,N_st
+        if (residual_norm(k) > 1.d8) then
+          print *, 'Davidson failed'
+          stop -1
+        endif
+      enddo
+      if (converged) then
+        exit
+      endif
+
+      logical, external :: qp_stop
+      if (qp_stop()) then
+        converged = .True.
+        exit
+      endif
+
+
+    enddo
+
+    ! Re-contract U and update W
+    ! --------------------------------
+
+    call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0,      &
+        W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
+    do k=1,N_st_diag_in
+      do i=1,sze
+        W(i,k) = u_in(i,k)
+      enddo
+    enddo
+
+    call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0,      &
+        U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
+
+    do k=1,N_st_diag_in
+      do i=1,sze
+        U(i,k) = u_in(i,k)
+      enddo
+    enddo
+
+  enddo
+
+
+  call nullify_small_elements(sze,N_st_diag_in,U,size(U,1),threshold_davidson_pt2)
+  do k=1,N_st_diag_in
+    do i=1,sze
+      u_in(i,k) = U(i,k)
+    enddo
+  enddo
+
+  do k=1,N_st_diag_in
+    energies(k) = lambda(k)
+  enddo
+  write_buffer = '======'
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//' ================ ==========='
+  enddo
+  write(6,'(A)') trim(write_buffer)
+  write(6,'(A)') ''
+  call write_time(6)
+
+  deallocate(W)
+
+  deallocate (                                                       &
+      residual_norm,                                                 &
+      U, overlap,                                                    &
+      h, y, s_tmp,                                                   &
+      lambda                                                         &
+      )
+  FREE nthreads_davidson
+end
+
+
+subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze)
+  implicit none
+  BEGIN_DOC
+  ! Routine that computes the diagonal part of the dressing 
+  !
+  ! v(i) += u(i) * dress_diag(i)
+  !
+  ! !!!!!!!! WARNING !!!!!!!! the vector v is not initialized 
+  !
+  ! !!!!!!!! SO MAKE SURE THERE ARE SOME MEANINGFUL VALUES IN THERE
+  END_DOC
+  integer, intent(in)              :: N_st,sze
+  double precision, intent(in)     :: u(sze,N_st),dress_diag(sze)
+  double precision, intent(inout)  :: v(sze,N_st)
+  integer :: i,istate
+  do istate = 1, N_st
+   do i = 1, sze
+     v(i,istate) += dress_diag(i) * u(i,istate)
+   enddo
+  enddo
+end
+
+
+
+
+
+
+
+
+
+
+

src/dav_general_mat/dav_dressed_ext_rout.irp.f

@@ -147,7 +147,7 @@ subroutine davidson_general_ext_rout_dressed(u_in,H_jj,energies,sze,N_st,N_st_di
   TOUCH nthreads_davidson
   call write_int(6,N_st,'Number of states')
   call write_int(6,N_st_diag,'Number of states in diagonalization')
-  call write_int(6,sze,'Number of determinants')
+  call write_int(6,sze,'Number of basis function')
   call write_int(6,nproc_target,'Number of threads for diagonalization')
   call write_double(6, r1, 'Memory(Gb)')
   if (disk_based) then
@@ -387,7 +387,7 @@ subroutine davidson_general_ext_rout_dressed(u_in,H_jj,energies,sze,N_st,N_st_di
 
         if (k <= N_st) then
           residual_norm(k) = u_dot_u(U(1,shift2+k),sze)
-          to_print(1,k) = lambda(k) + nuclear_repulsion
+          to_print(1,k) = lambda(k) 
           to_print(2,k) = residual_norm(k)
         endif
       enddo

src/dav_general_mat/dav_ext_rout.irp.f

@@ -3,7 +3,7 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co
   use mmap_module
   implicit none
   BEGIN_DOC
-  ! Davidson diagonalization with specific diagonal elements of the H matrix
+  ! Generic Davidson diagonalization 
   !
   ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
   !
@@ -221,7 +221,6 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co
       if ((iter > 1).or.(itertot == 1)) then
         ! Compute |W_k> = \sum_i |i><i|H|u_k>
         ! -----------------------------------
-
          ! Gram-Schmidt to orthogonalize all new guess with the previous vectors 
           call ortho_qr(U,size(U,1),sze,shift2)
           call ortho_qr(U,size(U,1),sze,shift2)
@@ -345,6 +344,9 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co
 
     enddo
 
+    ! Re-contract U and update W
+    ! --------------------------------
+
     call dgemm('N','N', sze, N_st_diag, shift2, 1.d0,      &
         W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
     do k=1,N_st_diag
@@ -360,8 +362,8 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,co
         U(i,k) = u_in(i,k)
       enddo
     enddo
-      call ortho_qr(U,size(U,1),sze,N_st_diag)
+    call ortho_qr(U,size(U,1),sze,N_st_diag)
-      call ortho_qr(U,size(U,1),sze,N_st_diag)
+    call ortho_qr(U,size(U,1),sze,N_st_diag)
     do j=1,N_st_diag
       k=1
       do while ((k<sze).and.(U(k,j) == 0.d0))

src/davidson_dressed/diagonalize_ci.irp.f

@@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
   enddo
   do j=1,min(N_det,N_states)
     write(st,'(I4)') j
-    call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st))
+    call write_double(6,CI_energy_dressed(j),'Energy dressed of state '//trim(st))
     call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st))
   enddo
 

src/determinants/utils.irp.f

@@ -8,6 +8,7 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
  double precision :: hij
  integer :: degree(N_det),idx(0:N_det)
  call  i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij)
+ print*,'Providing the H_matrix_all_dets ...'
  !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hij,degree,idx,k) &
  !$OMP SHARED (N_det, psi_det, N_int,H_matrix_all_dets)
  do i =1,N_det
@@ -18,6 +19,26 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
   enddo
  enddo
  !$OMP END PARALLEL DO
+ print*,'H_matrix_all_dets done '
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, H_matrix_diag_all_dets,(N_det) ]
+  use bitmasks
+ implicit none
+ BEGIN_DOC
+ ! |H| matrix on the basis of the Slater determinants defined by psi_det
+ END_DOC
+ integer :: i
+ double precision :: hij
+ integer :: degree(N_det)
+ call  i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij)
+ !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,hij,degree) &
+ !$OMP SHARED (N_det, psi_det, N_int,H_matrix_diag_all_dets)
+ do i =1,N_det
+   call  i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij)
+   H_matrix_diag_all_dets(i) = hij
+ enddo
+ !$OMP END PARALLEL DO
 END_PROVIDER
 
 

src/dft_one_e/mu_erf_dft.irp.f

@@ -72,6 +72,7 @@ BEGIN_PROVIDER [double precision, mu_of_r_dft_average]
   r(3) = final_grid_points(3,i)
   call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
   rho = dm_a + dm_b
+  if(mu_of_r_dft(i).gt.1.d+3)cycle
   mu_of_r_dft_average += rho * mu_of_r_dft(i) * final_weight_at_r_vector(i)
  enddo
  mu_of_r_dft_average = mu_of_r_dft_average / dble(elec_alpha_num + elec_beta_num)

src/dft_utils_in_r/ints_grad.irp.f

@@ -0,0 +1,39 @@
+ BEGIN_PROVIDER [ double precision, mo_grad_ints, (mo_num, mo_num,3)]
+ implicit none
+ BEGIN_DOC
+! mo_grad_ints(i,j,m) = <phi_i^MO | d/dx | phi_j^MO>
+ END_DOC
+ integer :: i,j,ipoint,m
+ double precision :: weight
+ mo_grad_ints = 0.d0
+ do m = 1, 3
+  do ipoint = 1, n_points_final_grid
+   weight = final_weight_at_r_vector(ipoint)
+   do j = 1, mo_num
+    do i = 1, mo_num
+     mo_grad_ints(i,j,m) += mos_grad_in_r_array(j,ipoint,m) * mos_in_r_array(i,ipoint) * weight
+    enddo
+   enddo
+  enddo
+ enddo
+
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ double precision, mo_grad_ints_transp, (3,mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! mo_grad_ints(i,j,m) = <phi_i^MO | d/dx | phi_j^MO>
+ END_DOC
+ integer :: i,j,ipoint,m
+ double precision :: weight
+ do m = 1, 3
+  do j = 1, mo_num
+   do i = 1, mo_num
+    mo_grad_ints_transp(m,i,j) = mo_grad_ints(i,j,m)
+   enddo
+  enddo
+ enddo
+
+
+END_PROVIDER 

src/iterations/print_summary.irp.f

@@ -98,7 +98,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
     enddo
   endif
 
-  call print_energy_components()
+!  call print_energy_components()
 
 end subroutine
 

src/mo_two_e_ints/map_integrals.irp.f

@@ -302,21 +302,21 @@ end
  integer(key_kind)              :: idx
  double precision :: tmp
 
- icount = 1 ! Avoid division by zero
+!icount = 1 ! Avoid division by zero
- do j=1,mo_num
+!do j=1,mo_num
-   do i=1,j-1
+!  do i=1,j-1
-    call two_e_integrals_index(i,j,j,i,idx)
+!   call two_e_integrals_index(i,j,j,i,idx)
-    !DIR$ FORCEINLINE
+!   !DIR$ FORCEINLINE
-    call map_get(mo_integrals_map,idx,tmp)
+!   call map_get(mo_integrals_map,idx,tmp)
-    banned_excitation(i,j) = dabs(tmp) < 1.d-14
+!   banned_excitation(i,j) = dabs(tmp) < 1.d-14
-    banned_excitation(j,i) = banned_excitation(i,j)
+!   banned_excitation(j,i) = banned_excitation(i,j)
-    if (banned_excitation(i,j)) icount = icount+2
+!   if (banned_excitation(i,j)) icount = icount+2
-  enddo
+! enddo
- enddo
+!enddo
- use_banned_excitation =  (mo_num*mo_num) / icount <= 100  !1%
+!use_banned_excitation =  (mo_num*mo_num) / icount <= 100  !1%
- if (use_banned_excitation) then
+!if (use_banned_excitation) then
-   print *, 'Using sparsity of exchange integrals'
+!  print *, 'Using sparsity of exchange integrals'
- endif
+!endif
 
 END_PROVIDER
 

src/tools/NEED

@@ -2,3 +2,4 @@ fci
 mo_two_e_erf_ints
 aux_quantities
 hartree_fock
+two_body_rdm