pert 2rdm work for the Li atom with HF

src/casscf/test_pert_2rdm.irp.f

@@ -24,7 +24,7 @@ program test_pert_2rdm
      integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
 !    if(dabs(pert_2rdm_provider(ii,jj,kk,ll) * integral).gt.1.d-12)then
 !     print*,i,j,k,l
-!     print*,pert_2rdm_provider(ii,jj,kk,ll) * integral,pert_2rdm_provider(ii,jj,kk,ll), pert_2rdm_provider(ii,jj,kk,ll), integral
+!     print*,pert_2rdm_provider(ii,jj,kk,ll) , integral, pert_2rdm_provider(ii,jj,kk,ll)* integral
 !    endif
      accu    += pert_2rdm_provider(ii,jj,kk,ll) * integral
     enddo

src/cipsi/pert_rdm_providers.irp.f

@@ -174,9 +174,13 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo
       Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
 
       sum_e_pert = 0d0
+      integer                        :: exc(0:2,2,2)
       integer                        :: degree
+      double precision               :: phase
+      call get_excitation(det,HF_bitmask,exc,degree,phase,N_int) 
       call get_excitation_degree(det,HF_bitmask,degree,N_int)
-      if(degree == 2)cycle
+      if(degree == 1)cycle
+!     if (exc(0,1,1) .ne. 1) cycle !only double alpha/beta
       do istate=1,N_states
         delta_E = E0(istate) - Hii + E_shift
         alpha_h_psi = mat(istate, p1, p2)
@@ -187,7 +191,7 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo
         endif
         e_pert = 0.5d0 * (tmp - delta_E)
         coef(istate) = e_pert / alpha_h_psi
-        print*,e_pert,coef,alpha_h_psi
+!       print*,coef,alpha_h_psi,e_pert
         pt2(istate) = pt2(istate) + e_pert
         variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi
         norm(istate) = norm(istate) + coef(istate) * coef(istate)

src/cipsi/update_2rdm.irp.f

@@ -18,7 +18,7 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio
   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)
+   contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i)  * coef(j)
   enddo
   ! case of single excitations 
   if(degree == 1)then
@@ -36,8 +36,8 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio
    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)
+ call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
-!nkeys = 0
+ nkeys = 0
 
 end
 
@@ -69,7 +69,7 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v
   ispin = 2
   other_spin = 1
  endif
- print*,'h1,p1,s1',h1,p1,ispin
+!print*,'h1,p1,s1',h1,p1,ispin
  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
@@ -94,39 +94,41 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v
   keys(4,nkeys) = p1
  enddo 
  !update the same spin part 
-!do i = 1, n_occ_ab(ispin)
+ do i = 1, n_occ_ab(ispin)
-! h2 = occ(i,ispin)
+  h2 = occ(i,ispin)
-! if(list_orb_reverse_pert_rdm(h2).lt.0)return
+  if(list_orb_reverse_pert_rdm(h2).lt.0)return
-! h2 = list_orb_reverse_pert_rdm(h2)
+  h2 = list_orb_reverse_pert_rdm(h2)
 
-! nkeys += 1
+  nkeys += 1
-! values(nkeys) = 0.5d0 * contrib * phase
+  values(nkeys) = 0.5d0 * contrib * phase
-! keys(1,nkeys) = h1
+  keys(1,nkeys) = h1
-! keys(2,nkeys) = h2
+  keys(2,nkeys) = h2
-! keys(3,nkeys) = p1
+  keys(3,nkeys) = p1
-! keys(4,nkeys) = h2
+  keys(4,nkeys) = h2
 
-! nkeys += 1
+  nkeys += 1
-! values(nkeys) = - 0.5d0 * contrib * phase
+  values(nkeys) = - 0.5d0 * contrib * phase
-! keys(1,nkeys) = h1
+  keys(1,nkeys) = h1
-! keys(2,nkeys) = h2
+  keys(2,nkeys) = h2
-! keys(3,nkeys) = h2
+  keys(3,nkeys) = h2
-! keys(4,nkeys) = p1
+  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
+  nkeys += 1
-! values(nkeys) = - 0.5d0 * contrib * phase
+  values(nkeys) = 0.5d0 * contrib * phase
-! keys(1,nkeys) = h2
+  keys(1,nkeys) = h2
-! keys(2,nkeys) = h1
+  keys(2,nkeys) = h1
-! keys(3,nkeys) = p1
+  keys(3,nkeys) = h2
-! keys(4,nkeys) = h2
+  keys(4,nkeys) = p1
-!enddo 
+
+  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 
+ call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
+ nkeys = 0
 
 end
 
@@ -160,12 +162,14 @@ subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_
   keys(2,nkeys) = h2
   keys(3,nkeys) = p1
   keys(4,nkeys) = p2
+! print*,contrib 
   nkeys += 1
   values(nkeys) = 0.5d0 * contrib * phase
   keys(1,nkeys) = p1
   keys(2,nkeys) = p2
   keys(3,nkeys) = h1
   keys(4,nkeys) = h2 
+! print*,contrib 
 
  else 
   if (exc(0,1,1) == 2) then
@@ -218,6 +222,8 @@ subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_
   keys(3,nkeys) = p1
   keys(4,nkeys) = p2
  endif
+ call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
+ nkeys = 0
  
 end