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

plugins/DavidsonDressed/diagonalize_CI.irp.f

@@ -63,8 +63,8 @@ END_PROVIDER
     call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,&
         size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,&
         N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1)
-    call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
+!    call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
-        N_states_diag,size(CI_eigenvectors_dressed,1))
+!        N_states_diag,size(CI_eigenvectors_dressed,1))
     
     
   else if (diag_algorithm == "Lapack") then

plugins/dress_zmq/dress_general.irp.f

@@ -29,27 +29,26 @@ subroutine run_dressing(N_st,energy)
     E_new = 0.d0
     delta_E = 1.d0
     iteration = 0
-    do while (delta_E > thresh_dress)
+    do iteration=1,n_it_dress_max
-      iteration += 1
       print *,  '===============================================' 
       print *,  'Iteration', iteration, '/', n_it_dress_max
       print *,  '===============================================' 
       print *,  ''
-      E_old = dress_e0_denominator(1) !sum(ci_energy_dressed(1:N_states))
+      E_old = sum(psi_energy(:))
       do i=1,N_st
         call write_double(6,ci_energy_dressed(i),"Energy")
       enddo
       call diagonalize_ci_dressed
-      E_new = dress_e0_denominator(1) !sum(ci_energy_dressed(1:N_states))
+      E_new = sum(psi_energy(:))
 
       delta_E = (E_new - E_old)/dble(N_states)
       print *,  ''
       call write_double(6,thresh_dress,"thresh_dress")
-      call write_double(6,delta_E,"delta_E")
+      call write_double(6,delta_E,"delta_E (undressed)")
       delta_E = dabs(delta_E)
       call save_wavefunction
 !     call ezfio_set_dress_zmq_energy(ci_energy_dressed(1))
-      if (iteration >= n_it_dress_max) then
+      if (delta_E < thresh_dress) then
         exit
       endif
     enddo

plugins/dress_zmq/dressing.irp.f

@@ -78,7 +78,7 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det,2) ]
 
   delta_ij = 0d0
 
-  E_CI_before(:) = dress_E0_denominator(:) + nuclear_repulsion
+  E_CI_before(:) = psi_energy(:) + nuclear_repulsion
   threshold_selectors = 1.d0
   threshold_generators = 1d0 
 !  if(errr /= 0d0) then

plugins/dress_zmq/dressing_vector.irp.f

@@ -18,12 +18,6 @@
      dressing_column_h(j,k) = delta_ij(k,j,1) 
      dressing_column_s(j,k) = delta_ij(k,j,2) 
    enddo
-!   tmp = u_dot_v(dressing_column_h(1,k), psi_coef(1,k), N_det) &
-!     - dressing_column_h(l,k) * psi_coef(l,k)
-!   dressing_column_h(l,k) -= tmp * f
-!   tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) &
-!     - dressing_column_s(l,k) * psi_coef(l,k)
-!   dressing_column_s(l,k) -= tmp * f
  enddo
 
 END_PROVIDER

plugins/dress_zmq/energy.irp.f

@@ -13,13 +13,24 @@ BEGIN_PROVIDER [ double precision, dress_E0_denominator, (N_states) ]
  END_DOC
  integer :: i
  if (initialize_dress_E0_denominator) then
-  call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1))
+   if (h0_type == "EN") then
-  do i=N_det+1,N_states
+     dress_E0_denominator(1:N_states) = psi_energy(1:N_states)
-    dress_E0_denominator(i) = 0.d0
+   else if (h0_type == "Barycentric") then
-  enddo
+!     dress_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
+     dress_E0_denominator(1:N_states) = minval(diagonal_H_matrix_on_psi_det(1:N_det))
+   else
+     print *,  h0_type, ' not implemented'
+     stop
+   endif
+!  call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1))
+!  do i=N_det+1,N_states
+!    dress_E0_denominator(i) = 0.d0
+!  enddo
   call write_double(6,dress_E0_denominator(1)+nuclear_repulsion, 'dress Energy denominator')
  else
    dress_E0_denominator = -huge(1.d0)
  endif
 END_PROVIDER
 
+
+

plugins/shiftedbk/shifted_bk_routines.irp.f

@@ -4,6 +4,7 @@
 &BEGIN_PROVIDER [ double precision, a_s2_i, (N_det, Nproc) ]
   implicit none
   current_generator_(:) = 0
+  fock_diag_tmp_(:,:,:) = 0.d0
   a_h_i = 0d0
   a_s2_i = 0d0
  END_PROVIDER
@@ -52,7 +53,7 @@ subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minili
 
 
   do i=1,Nstates
-    de = E0_denominator(i) - haa
+    de = dress_E0_denominator(i) - haa
     if(DABS(de) < 1D-5) cycle
 
     c_alpha = a_h_psi(i) / de
@@ -76,23 +77,4 @@ BEGIN_PROVIDER [ logical, initialize_E0_denominator ]
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, E0_denominator, (N_states) ]
-  implicit none
-  BEGIN_DOC
-  ! E0 in the denominator of the PT2
-  END_DOC
-  if (initialize_E0_denominator) then
-   if (h0_type == "EN") then
-     E0_denominator(1:N_states) = psi_energy(1:N_states)
-   else if (h0_type == "Barycentric") then
-     E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
-   else
-     print *,  h0_type, ' not implemented'
-     stop
-   endif
-  else
-    E0_denominator = -huge(1.d0)
-  endif
-END_PROVIDER
-
 

src/Davidson/diagonalization_hs2_dressed.irp.f

@@ -35,7 +35,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
   double precision, intent(inout) :: u_in(dim_in,N_st_diag)
   double precision, intent(out)  :: energies(N_st_diag), s2_out(N_st_diag)
   integer, intent(in)            :: dressing_state
-  double precision, allocatable  :: H_jj(:), S2_jj(:)
+  double precision, allocatable  :: H_jj(:)
   
   double precision, external     :: diag_H_mat_elem, diag_S_mat_elem
   integer                        :: i,k
@@ -44,7 +44,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
   PROVIDE mo_bielec_integrals_in_map
-  allocate(H_jj(sze),S2_jj(sze))
+  allocate(H_jj(sze))
   
   H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)
   !$OMP PARALLEL DEFAULT(NONE)                                       &
@@ -60,13 +60,13 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
   if (dressing_state > 0) then
     do k=1,N_st
       do i=1,sze
-        H_jj(i) += u_in(i,k) * dressing_column_h(i,k)
+        H_jj(i)  += u_in(i,k) * dressing_column_h(i,k)
       enddo
     enddo
   endif
 
   call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state)
-  deallocate (H_jj,S2_jj)
+  deallocate (H_jj)
 end
 
 
@@ -254,9 +254,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
           dressing_column_h, size(dressing_column_h,1), s_tmp, size(s_tmp,1), &
           1.d0, W(1,shift+1), size(W,1))
 
-        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
+!        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
-          dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), &
+!          dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), &
-          1.d0, S(1,shift+1), size(S,1))
+!          1.d0, S(1,shift+1), size(S,1))
 
 
         call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
@@ -267,13 +267,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
           psi_coef, size(psi_coef,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, sze, 1.d0, &
+!        call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
-          dressing_column_s, size(dressing_column_s,1), &
+!          dressing_column_s, size(dressing_column_s,1), &
-          U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+!          U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
-
+!
-        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
+!        call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
-          psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
+!          psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
-          1.d0, S(1,shift+1), size(S,1))
+!          1.d0, S(1,shift+1), size(S,1))
 
       endif
 

src/Determinants/slater_rules.irp.f

@@ -509,11 +509,12 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
   select case (degree)
     case (2)
       call get_double_excitation(key_i,key_j,exc,phase,Nint)
-      
+      ! Mono alpha, mono beta
       if (exc(0,1,1) == 1) then
-        ! Mono alpha, mono beta
+        if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then
+          s2 =  -phase
+        endif
         if(exc(1,1,1) == exc(1,2,2) )then
-         if(exc(1,1,2) == exc(1,2,1)) s2 = -phase !!!!!
          hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
         else if (exc(1,2,1) ==exc(1,1,2))then
          hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
@@ -524,8 +525,8 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
              exc(1,2,1),                                              &
              exc(1,2,2) ,mo_integrals_map)
         endif
+      ! Double alpha
       else if (exc(0,1,1) == 2) then
-        ! Double alpha
         hij = phase*(get_mo_bielec_integral(                         &
             exc(1,1,1),                                              &
             exc(2,1,1),                                              &
@@ -536,8 +537,8 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
             exc(2,1,1),                                              &
             exc(2,2,1),                                              &
             exc(1,2,1) ,mo_integrals_map) )
+      ! Double beta
       else if (exc(0,1,2) == 2) then
-        ! Double beta
         hij = phase*(get_mo_bielec_integral(                         &
             exc(1,1,2),                                              &
             exc(2,1,2),                                              &
@@ -553,13 +554,13 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
       call get_mono_excitation(key_i,key_j,exc,phase,Nint)
       !DIR$ FORCEINLINE
       call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
+      ! Mono alpha
       if (exc(0,1,1) == 1) then
-        ! Mono alpha
         m = exc(1,1,1)
         p = exc(1,2,1)
         spin = 1
+      ! Mono beta
       else
-        ! Mono beta
         m = exc(1,1,2)
         p = exc(1,2,2)
         spin = 2
@@ -567,7 +568,7 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
       call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij)
       
     case (0)
-      print *," ZERO"
+      print *,irp_here,": ZERO"
       double precision, external :: diag_S_mat_elem
       s2 = diag_S_mat_elem(key_i,Nint)
       hij = diag_H_mat_elem(key_i,Nint)

src/Utils/map_functions.irp.f

@@ -9,6 +9,9 @@ subroutine map_save_to_disk(filename,map)
   integer*8                      :: i,k
   integer                        :: j
 
+  integer*8 :: n_elements
+  n_elements = int(map % n_elements,8)
+
 
   if (map % consolidated) then
     stop 'map already consolidated'
@@ -17,11 +20,11 @@ subroutine map_save_to_disk(filename,map)
   call mmap(trim(filename)//'_consolidated_idx', (/ map % map_size + 2_8 /), 8, fd(1), .False., c_pointer(1))
   call c_f_pointer(c_pointer(1),map % consolidated_idx, (/ map % map_size +2_8/))
 
-  call mmap(trim(filename)//'_consolidated_key', (/ map % n_elements /), cache_key_kind, fd(2), .False., c_pointer(2))
+  call mmap(trim(filename)//'_consolidated_key', (/ n_elements /), cache_key_kind, fd(2), .False., c_pointer(2))
-  call c_f_pointer(c_pointer(2),map % consolidated_key, (/ map % n_elements /))
+  call c_f_pointer(c_pointer(2),map % consolidated_key, (/ n_elements /))
 
-  call mmap(trim(filename)//'_consolidated_value', (/ map % n_elements /), integral_kind, fd(3), .False., c_pointer(3))
+  call mmap(trim(filename)//'_consolidated_value', (/ n_elements /), integral_kind, fd(3), .False., c_pointer(3))
-  call c_f_pointer(c_pointer(3),map % consolidated_value, (/ map % n_elements /))
+  call c_f_pointer(c_pointer(3),map % consolidated_value, (/ n_elements /))
 
   if (.not.associated(map%consolidated_key)) then
     stop 'cannot consolidate map : consolidated_key not associated'
@@ -52,9 +55,6 @@ subroutine map_save_to_disk(filename,map)
   map % consolidated_idx (map % map_size + 2_8) = k
   map % consolidated = .True.
 
-  integer*8 :: n_elements
-  n_elements = int(map % n_elements,8)
-
   print *,  'Writing data to disk...'
   call msync ( (/ map % map_size + 2_8 /),    8, fd(1), c_pointer(1))
   call msync ( (/ n_elements /), cache_key_kind, fd(2), c_pointer(2))