working on pert rdms

src/cipsi/pert_rdm_providers.irp.f

@@ -0,0 +1,166 @@
+
+use bitmasks
+
+BEGIN_PROVIDER [logical , pert_2rdm ]
+ implicit none
+ pert_2rdm = .False.
+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
+
+subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf, psi_det_connection, psi_coef_connection, n_det_connection)
+  use bitmasks
+  use selection_types
+  implicit none
+  
+  integer, intent(in)           :: n_det_connection
+  double precision, intent(in)  :: psi_coef_connection(n_det_connection,N_states)
+  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)
+  double precision, intent(inout) :: pt2(N_states)
+  double precision, intent(inout) :: variance(N_states)
+  double precision, intent(inout) :: norm(N_states)
+  type(selection_buffer), intent(inout) :: buf
+  logical :: ok
+  integer :: s1, s2, p1, p2, ib, j, istate
+  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 == 'SOP') then
+    j = det_to_occ_pattern(i_generator)
+    E_shift = psi_det_Hii(i_generator) - psi_occ_pattern_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
+
+      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
+        pt2(istate) = pt2(istate) + e_pert
+        variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi
+        norm(istate) = norm(istate) + coef(istate) * coef(istate)
+
+        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,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
+end
+
+

src/cipsi/selection.irp.f

@@ -1,9 +1,5 @@
-use bitmasks
 
-BEGIN_PROVIDER [logical , pert_2rdm ]
+use bitmasks
- implicit none
- pert_2rdm = .False.
-END_PROVIDER 
 
 BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ]
  implicit none
@@ -768,148 +764,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   end do
 end
 
-
-subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf, psi_det_connection, psi_coef_connection, n_det_connection)
-  use bitmasks
-  use selection_types
-  implicit none
-  
-  integer, intent(in)           :: n_det_connection
-  double precision, intent(in)  :: psi_coef_connection(n_det_connection,N_states)
-  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)
-  double precision, intent(inout) :: pt2(N_states)
-  double precision, intent(inout) :: variance(N_states)
-  double precision, intent(inout) :: norm(N_states)
-  type(selection_buffer), intent(inout) :: buf
-  logical :: ok
-  integer :: s1, s2, p1, p2, ib, j, istate
-  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 == 'SOP') then
-    j = det_to_occ_pattern(i_generator)
-    E_shift = psi_det_Hii(i_generator) - psi_occ_pattern_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
-
-      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
-        pt2(istate) = pt2(istate) + e_pert
-        variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi
-        norm(istate) = norm(istate) + coef(istate) * coef(istate)
-
-        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,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
-end
-
-
 subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting)
   use bitmasks
   implicit none

src/cipsi/update_2rdm.irp.f

@@ -9,5 +9,57 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio
  double precision, intent(in)  :: psi_coef_connection(n_det_connection, N_states)
  integer,           intent(inout) :: keys(4,sze_buff),sze_buff
  double precision,  intent(inout) :: values(sze_buff)
+ integer :: i
+ 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(i,j) * phase * coef(j)
+  enddo
+  ! case of single excitations 
+  if(degree == 1)then
+   if (nkeys+ 2 * elec_alpha_num .ge. sze_buff)then
+    call update_rdms(nkeys,keys,values,sze_buff)
+    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_rdms(nkeys,keys,values,sze_buff)
+    nkeys = 0
+   endif
+   call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
+  endif
+ enddo
+
+end
+
+subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff)
+ implicit none
+ integer, intent(in) :: nkeys,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)
+ 
+
+
+end
+
+subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
+ implicit none
+ integer, intent(in) :: nkeys,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)
+
 
 end

src/two_body_rdm/orb_range_routines_openmp.irp.f

@@ -166,7 +166,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
             !$OMP          psi_bilinear_matrix_transp_order, N_st,   &
             !$OMP          psi_bilinear_matrix_order_transp_reverse, &
             !$OMP          psi_bilinear_matrix_columns_loc,          &
-            !$OMP          psi_bilinear_matrix_transp_rows_loc,norb, &
+            !$OMP          psi_bilinear_matrix_transp_rows_loc, &
             !$OMP          istart, iend, istep, irp_here,list_orb_reverse, n_states, state_weights, dim1, & 
             !$OMP          ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask)     &
             !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, c_2, &
@@ -348,13 +348,13 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
        enddo
        if(alpha_beta.or.spin_trace.or.alpha_alpha)then
        ! increment the alpha/beta part for single excitations
-       if (nkeys+ 2 * norb .ge. size(values)) then
+       if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then
          call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
          nkeys = 0
        endif
        call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
        ! increment the alpha/alpha part for single excitations
-       if (nkeys+4 * norb .ge. size(values)) then
+       if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then
          call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
          nkeys = 0
        endif
@@ -381,7 +381,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
           c_2(l) = u_t(l,k_a)
           c_average += c_1(l) * c_2(l) * state_weights(l)
         enddo
-        if (nkeys+4 .ge. size(values)) then
+        if (nkeys+4 .ge. sze_buff) then
           call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
           nkeys = 0
         endif
@@ -452,13 +452,13 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
        enddo
        if(alpha_beta.or.spin_trace.or.beta_beta)then
         ! increment the alpha/beta  part for single excitations
-        if (nkeys+2 * norb .ge. size(values)) then
+        if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then
           call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
           nkeys = 0
         endif
         call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
         ! increment the beta /beta  part for single excitations
-        if (nkeys+4 * norb .ge. size(values)) then
+        if (nkeys+4 * elec_alpha_num .ge. sze_buff) then
           call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
           nkeys = 0
         endif
@@ -484,7 +484,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
           c_2(l) = u_t(l,k_a)
           c_average += c_1(l) * c_2(l) * state_weights(l)
         enddo
-        if (nkeys+4 .ge. size(values)) then
+        if (nkeys+4 .ge. sze_buff) then
           call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
           nkeys = 0
         endif