density based mrpt2

plugins/MRPT/print_1h2p.irp.f

@@ -2,7 +2,7 @@ program print_1h2p
  implicit none
  read_wf = .True.
  touch read_wf
- call routine_2
+ call routine
 end
 
 subroutine routine_2
@@ -35,7 +35,8 @@ subroutine routine
  integer :: i,j,istate
  accu = 0.d0
  matrix_1h2p = 0.d0
- call H_apply_mrpt_2p(matrix_1h2p,N_det_ref)
+!call H_apply_mrpt_1h2p(matrix_1h2p,N_det_ref)
+ call give_1h2p_contrib(matrix_1h2p)
  do istate = 1, N_states
   do i = 1, N_det
    do j = 1, N_det 
@@ -44,6 +45,8 @@ subroutine routine
   enddo
   print*,accu(istate)
  enddo
+ call contrib_1h2p_dm_based(accu)
+ print*,accu(:)
 
  deallocate (matrix_1h2p)
 end

plugins/MRPT_Utils/density_matrix_based.irp.f

@@ -0,0 +1,193 @@
+subroutine contrib_1h2p_dm_based(accu)
+ implicit none
+ integer :: i_i,i_r,i_v,i_a,i_b
+ integer :: i,r,v,a,b
+ integer :: ispin,jspin
+ integer :: istate
+ double precision, intent(out) :: accu(N_states)
+ double precision :: active_int(n_act_orb,2)
+ double precision :: delta_e(n_act_orb,2,N_states)
+ double precision :: get_mo_bielec_integral
+ accu = 0.d0
+!do i_i = 1, 1
+ do i_i = 1, n_inact_orb
+  i = list_inact(i_i)
+! do i_r = 1, 1
+  do i_r = 1, n_virt_orb
+   r = list_virt(i_r)
+!  do i_v = 1, 1
+   do i_v = 1, n_virt_orb
+    v = list_virt(i_v)
+    do i_a = 1, n_act_orb
+     a = list_act(i_a)
+     active_int(i_a,1) = get_mo_bielec_integral(i,a,r,v,mo_integrals_map) ! direct
+     active_int(i_a,2) = get_mo_bielec_integral(i,a,v,r,mo_integrals_map) ! exchange
+     do istate = 1, N_states
+      do jspin=1, 2
+       delta_e(i_a,jspin,istate) = one_anhil(i_a,jspin,istate)                        &
+                                 - fock_virt_total_spin_trace(r,istate)               & 
+                                 - fock_virt_total_spin_trace(v,istate)               & 
+                                 + fock_core_inactive_total_spin_trace(i,istate)        
+       delta_e(i_a,jspin,istate) = 1.d0/delta_e(i_a,jspin,istate)  
+      enddo
+     enddo
+    enddo
+    do i_a = 1, n_act_orb
+     a = list_act(i_a)
+     do i_b = 1, n_act_orb
+!    do i_b = i_a, i_a
+      b = list_act(i_b)
+      do ispin = 1, 2 ! spin of (i --> r)
+       do jspin = 1, 2 ! spin of (a --> v)
+        if(ispin == jspin .and. r.le.v)cycle ! condition not to double count 
+        do istate = 1, N_states
+         if(ispin == jspin)then
+          accu(istate) += (active_int(i_a,1) - active_int(i_a,2)) * one_body_dm_mo_spin_index(a,b,istate,ispin)   &
+                                                          * (active_int(i_b,1) - active_int(i_b,2)) & 
+                                                          * delta_e(i_a,jspin,istate)
+         else 
+          accu(istate) += active_int(i_a,1)  * one_body_dm_mo_spin_index(a,b,istate,ispin) * delta_e(i_a,ispin,istate) & 
+                        * active_int(i_b,1) 
+         endif
+        enddo
+       enddo
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+
+
+end
+
+subroutine contrib_2h1p_dm_based(accu)
+ implicit none
+ integer :: i_i,i_j,i_v,i_a,i_b
+ integer :: i,j,v,a,b
+ integer :: ispin,jspin
+ integer :: istate
+ double precision, intent(out) :: accu(N_states)
+ double precision :: active_int(n_act_orb,2)
+ double precision :: delta_e(n_act_orb,2,N_states)
+ double precision :: get_mo_bielec_integral
+ accu = 0.d0
+ do i_i = 1, n_inact_orb
+  i = list_inact(i_i)
+  do i_j = 1, n_inact_orb
+   j = list_inact(i_j)
+   do i_v = 1, n_virt_orb
+    v = list_virt(i_v)
+    do i_a = 1, n_act_orb
+     a = list_act(i_a)
+     active_int(i_a,1) = get_mo_bielec_integral(i,j,v,a,mo_integrals_map) ! direct
+     active_int(i_a,2) = get_mo_bielec_integral(i,j,a,v,mo_integrals_map) ! exchange
+     do istate = 1, N_states
+      do jspin=1, 2
+!      delta_e(i_a,jspin,istate) = 
+!                                
+       delta_e(i_a,jspin,istate) = one_creat(i_a,jspin,istate)  - fock_virt_total_spin_trace(v,istate)  &
+                                 + fock_core_inactive_total_spin_trace(i,istate)       &
+                                 + fock_core_inactive_total_spin_trace(j,istate)        
+       delta_e(i_a,jspin,istate) = 1.d0/delta_e(i_a,jspin,istate)  
+      enddo
+     enddo
+    enddo
+    do i_a = 1, n_act_orb
+     a = list_act(i_a)
+     do i_b = 1, n_act_orb
+!    do i_b = i_a, i_a
+      b = list_act(i_b)
+      do ispin = 1, 2 ! spin of (i --> v)
+       do jspin = 1, 2 ! spin of (j --> a)
+        if(ispin == jspin .and. i.le.j)cycle ! condition not to double count 
+        do istate = 1, N_states
+         if(ispin == jspin)then
+          accu(istate) += (active_int(i_a,1) - active_int(i_a,2)) * one_body_dm_dagger_mo_spin_index(a,b,istate,ispin)   &
+                                                          * (active_int(i_b,1) - active_int(i_b,2)) & 
+                                                          * delta_e(i_a,jspin,istate)
+         else 
+          accu(istate) += active_int(i_a,1)  * one_body_dm_dagger_mo_spin_index(a,b,istate,ispin) * delta_e(i_a,ispin,istate) & 
+                        * active_int(i_b,1) 
+         endif
+        enddo
+       enddo
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+
+
+end
+
+
+subroutine contrib_2p_dm_based(accu)
+ implicit none
+ integer :: i_r,i_v,i_a,i_b,i_c,i_d
+ integer :: r,v,a,b,c,d
+ integer :: ispin,jspin
+ integer :: istate
+ double precision, intent(out) :: accu(N_states)
+ double precision :: active_int(n_act_orb,n_act_orb,2)
+ double precision :: delta_e(n_act_orb,n_act_orb,2,2,N_states)
+ double precision :: get_mo_bielec_integral
+ accu = 0.d0
+ do i_r = 1, n_virt_orb
+  r = list_virt(i_r)
+  do i_v = 1, n_virt_orb
+    v = list_virt(i_v)
+    do i_a = 1, n_act_orb
+     a = list_act(i_a)
+     do i_b = 1, n_act_orb
+      b = list_act(i_b)
+      active_int(i_a,i_b,1) = get_mo_bielec_integral(a,b,r,v,mo_integrals_map) ! direct
+      active_int(i_a,i_b,2) = get_mo_bielec_integral(a,b,v,r,mo_integrals_map) ! direct
+      do istate = 1, N_states
+       do jspin=1, 2 ! spin of i_a
+        do ispin = 1, 2 ! spin of i_b
+         delta_e(i_a,i_b,jspin,ispin,istate) = two_anhil(i_a,i_b,jspin,ispin,istate)    &
+                                   - fock_virt_total_spin_trace(r,istate)               & 
+                                   - fock_virt_total_spin_trace(v,istate)                 
+         delta_e(i_a,i_b,jspin,ispin,istate) = 1.d0/delta_e(i_a,i_b,jspin,ispin,istate)  
+        enddo
+       enddo
+      enddo
+     enddo
+    enddo
+    ! diagonal terms 
+    do i_a = 1, n_act_orb
+     a = list_act(i_a)
+     do i_b = 1, n_act_orb
+      b = list_act(i_b)
+      do ispin = 1, 2 ! spin of (a --> r)
+       do jspin = 1, 2 ! spin of (b --> v)
+        if(ispin == jspin .and. r.le.v)cycle ! condition not to double count 
+        if(ispin == jspin .and. a.le.b)cycle ! condition not to double count 
+        do istate = 1, N_states
+         if(ispin == jspin)then
+          double precision :: contrib_spin
+          if(ispin == 1)then
+           contrib_spin = two_body_dm_aa_diag_act(i_a,i_b)
+          else
+           contrib_spin = two_body_dm_bb_diag_act(i_a,i_b)
+          endif
+          accu(istate) += (active_int(i_a,i_b,1) - active_int(i_a,i_b,2)) * contrib_spin  &
+                        * (active_int(i_a,i_b,1) - active_int(i_a,i_b,2)) & 
+                        * delta_e(i_a,i_b,ispin,jspin,istate)
+         else 
+          accu(istate) += 0.5d0 * active_int(i_a,i_b,1)  * two_body_dm_ab_diag_act(i_a,i_b) * delta_e(i_a,i_b,ispin,jspin,istate) & 
+                                * active_int(i_a,i_b,1) 
+         endif
+        enddo
+       enddo
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
+
+
+end
+

plugins/MRPT_Utils/energies_cas.irp.f

@@ -648,10 +648,11 @@ END_PROVIDER
         double precision ::  coef,contrib
         coef = psi_ref_coef(i,j) !* psi_ref_coef(i,j)
         psi_in_out_coef(i,j) = coef * hij 
-!       if(vorb == 1.and. iorb == 1)then
+        if(orb_i == 5 .and. orb_v == 20)then
-!       if(vorb == 1.and. iorb == 3)then
+!       if(orb_i == 2 .and. orb_v == 6 )then
-!       print*, i,hij,coef
+         print*, i, ispin
-!       endif
+         print*, coef * hij,coef,hij
+        endif
         norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) 
       enddo
      enddo
@@ -664,6 +665,10 @@ END_PROVIDER
        norm_no_inv(j,ispin) = norm(j,ispin)
        one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin)
        norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
+       if(orb_i == 5 .and. orb_v == 20)then
+!      if(orb_i == 2 .and. orb_v == 6 )then
+        print*,ispin ,norm(j,ispin)
+       endif
       endif
      enddo
      do i = 1, N_det_ref
@@ -679,21 +684,32 @@ END_PROVIDER
      do state_target = 1, N_states
       energies_alpha_beta(state_target, ispin) = 0.d0
       if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
-!      call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
        call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
        energies_alpha_beta(state_target, ispin) +=  energies(state_target) 
+       if(orb_i == 5 .and. orb_v == 20)then
+!      if(orb_i == 2 .and. orb_v == 6 )then
+        print*, ispin, energy_cas_dyall_no_exchange(1) , energies_alpha_beta(state_target, ispin)
+        print*, ispin, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target, ispin)
+       endif
       endif
      enddo
     enddo ! ispin 
    do state_target = 1, N_states
     if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
-!    one_anhil_one_creat_inact_virt(iorb,vorb,state_target) =  energy_cas_dyall(state_target) - &
      one_anhil_one_creat_inact_virt(iorb,vorb,state_target) =  energy_cas_dyall_no_exchange(state_target) - &
       ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) & 
      /( norm_bis(state_target,1) +  norm_bis(state_target,2) )
     else  
      one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0
     endif
+    if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.30584271462d0) < 1.d-11)then
+     print*, orb_i,orb_v
+     print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,1) / norm_bis(state_target,1)
+     print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,2) / norm_bis(state_target,2)
+     print*, fock_core_inactive_total_spin_trace(orb_i,1)
+     print*, fock_virt_total_spin_trace(orb_v,1)
+     print*, one_anhil_one_creat_inact_virt(iorb,vorb,state_target)
+    endif
    enddo
   enddo
  enddo

plugins/MRPT_Utils/excitations_cas.irp.f

@@ -774,6 +774,7 @@ subroutine i_H_j_dyall_no_exchange(key_i,key_j,Nint,hij)
       
     case (0)
       hij = diag_H_mat_elem_no_elec_check_no_exchange(key_i,Nint)
+!     hij = 0.d0
   end select
 end
 
@@ -798,7 +799,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
   ! alpha - alpha 
   do i = 1, elec_num_tab_local(1)
    iorb =  occ(i,1)
-   diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb)
+   diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb)
    do j = i+1, elec_num_tab_local(1)
     jorb = occ(j,1)
     diag_H_mat_elem_no_elec_check_no_exchange +=  mo_bielec_integral_jj(jorb,iorb)
@@ -808,7 +809,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
   ! beta - beta   
   do i = 1, elec_num_tab_local(2)
    iorb =  occ(i,2)
-   diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb)
+   diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb)
    do j = i+1, elec_num_tab_local(2)
     jorb = occ(j,2)
     diag_H_mat_elem_no_elec_check_no_exchange +=  mo_bielec_integral_jj(jorb,iorb)
@@ -826,6 +827,8 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
   enddo 
   
 
+! return
+
   ! alpha - core-act
   do i = 1, elec_num_tab_local(1)
    iorb =  occ(i,1)
@@ -833,6 +836,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
     jorb = list_core_inact(j)
     diag_H_mat_elem_no_elec_check_no_exchange +=  2.d0 * mo_bielec_integral_jj(jorb,iorb) 
     core_act_exchange(1) += - mo_bielec_integral_jj_exchange(jorb,iorb)
+!   diag_H_mat_elem_no_elec_check_no_exchange += core_act_exchange(1) 
    enddo
   enddo 
 
@@ -843,6 +847,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
     jorb = list_core_inact(j)
     diag_H_mat_elem_no_elec_check_no_exchange +=  2.d0 * mo_bielec_integral_jj(jorb,iorb) 
     core_act_exchange(2) += - mo_bielec_integral_jj_exchange(jorb,iorb)
+!   diag_H_mat_elem_no_elec_check_no_exchange += core_act_exchange(2) 
    enddo
   enddo 
   

plugins/MRPT_Utils/mrpt_dress.irp.f

@@ -113,19 +113,20 @@ subroutine mrpt_dress(delta_ij_,  Ndet,i_generator,n_selected,det_buffer,Nint,ip
       integer :: degree_scalar
 
       call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,index_i),degree_scalar,N_int)
+!     if(degree_scalar == 2)then
+!      hialpha = 0.d0
+!     endif
       if(dabs(hialpha).le.1.d-20)then
        do i_state = 1, N_states
         delta_e(i_state) = 1.d+20
        enddo
       else 
        call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
-      
 !      !!!!!!!!!!!!! SHIFTED BK 
 !      double precision :: hjj
 !      call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj)
 !      delta_e(1) = CI_electronic_energy(1) - hjj
 !      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
       endif
       hij_array(index_i) = hialpha
       do i_state = 1,N_states
@@ -139,6 +140,8 @@ subroutine mrpt_dress(delta_ij_,  Ndet,i_generator,n_selected,det_buffer,Nint,ip
       call omp_set_lock( psi_ref_bis_lock(index_i) )
       do j = 1, idx_alpha(0)
        index_j = idx_alpha(j)
+       !!!!!!!!!!!!!!!!!! WARNING TEST
+       if(index_j .ne. index_i)cycle
        do i_state=1,N_states
         ! standard dressing first order
          delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_inv_array(index_j,i_state)

plugins/MRPT_Utils/mrpt_utils.irp.f

@@ -104,23 +104,23 @@
   !print*, 'ionic   = ',ionic
    print*, '1h1p = ',accu
   
-   ! 1h1p third order
+  !! 1h1p third order
-   if(do_third_order_1h1p)then
+  !if(do_third_order_1h1p)then
-    delta_ij_tmp = 0.d0
+  ! delta_ij_tmp = 0.d0
-    call give_1h1p_sec_order_singles_contrib(delta_ij_tmp)
+  ! call give_1h1p_sec_order_singles_contrib(delta_ij_tmp)
-    accu = 0.d0
+  ! accu = 0.d0
-    do i_state = 1, N_states
+  ! do i_state = 1, N_states
-    do i = 1, N_det_ref
+  ! do i = 1, N_det_ref
-     write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,i_state)
+  !  write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,i_state)
-     do j = 1, N_det_ref
+  !  do j = 1, N_det_ref
-      accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
+  !   accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
-      delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
+  !   delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
-     enddo
+  !  enddo
-    enddo
+  ! enddo
-    second_order_pt_new_1h1p(i_state) = accu(i_state) 
+  ! second_order_pt_new_1h1p(i_state) = accu(i_state) 
-    enddo
+  ! enddo
-    print*, '1h1p(3)',accu
+  ! print*, '1h1p(3)',accu
-   endif
+  !endif
   
    ! 2h   
    delta_ij_tmp = 0.d0

plugins/MRPT_Utils/new_way.irp.f

@@ -129,6 +129,7 @@ subroutine give_2h1p_contrib(matrix_2h1p)
        integer :: kspin
        do jdet = 1, idx(0)
         if(idx(jdet).ne.idet)then
+!         cycle
         ! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator 
         ! are connected by the presence of the perturbers determinants |det_tmp>
          aorb  = index_orb_act_mono(idx(jdet),1)   !  a^{\dagger}_{aorb}
@@ -213,16 +214,18 @@ subroutine give_1h2p_contrib(matrix_1h2p)
  double precision :: active_int(n_act_orb,2)
  double precision :: hij,phase
 !matrix_1h2p = 0.d0
- 
  elec_num_tab_local = 0
  do inint = 1, N_int
   elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1))
   elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1))
  enddo
+!do i = 1, 1            ! First inactive 
  do i = 1, n_inact_orb  ! First inactive 
   iorb = list_inact(i)
+!  do v = 1, 1
    do v = 1, n_virt_orb  ! First virtual 
     vorb = list_virt(v)
+!   do r = 1, 1
     do r = 1, n_virt_orb    ! Second virtual
      rorb = list_virt(r) 
      ! take all the integral you will need for i,j,r fixed

plugins/MRPT_Utils/psi_active_prov.irp.f

@@ -14,6 +14,7 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)]
    psi_active(j,1,i) = iand(psi_ref(j,1,i),cas_bitmask(j,1,1))
    psi_active(j,2,i) = iand(psi_ref(j,2,i),cas_bitmask(j,1,1))
   enddo
+  call debug_det(psi_active(1,1,i),N_int)
  enddo
 END_PROVIDER
 

src/Determinants/density_matrix.irp.f

@@ -15,6 +15,57 @@
    enddo
 END_PROVIDER
 
+ BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ]
+ implicit none 
+ integer :: i,j,ispin,istate
+ ispin = 1
+  do istate = 1, N_states
+   do j = 1, mo_tot_num
+    do i = 1, mo_tot_num
+     one_body_dm_mo_spin_index(i,j,istate,ispin) = one_body_dm_mo_alpha(i,j,istate)
+    enddo
+   enddo
+  enddo
+
+ ispin = 2
+  do istate = 1, N_states
+   do j = 1, mo_tot_num
+    do i = 1, mo_tot_num
+     one_body_dm_mo_spin_index(i,j,istate,ispin) = one_body_dm_mo_beta(i,j,istate)
+    enddo
+   enddo
+  enddo
+
+ END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, one_body_dm_dagger_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ]
+ implicit none 
+ integer :: i,j,ispin,istate
+ ispin = 1
+  do istate = 1, N_states
+   do j = 1, mo_tot_num
+    one_body_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_body_dm_mo_alpha(j,j,istate)
+    do i = j+1, mo_tot_num
+     one_body_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_body_dm_mo_alpha(i,j,istate)
+     one_body_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_body_dm_mo_alpha(i,j,istate)
+    enddo
+   enddo
+  enddo
+
+ ispin = 2
+  do istate = 1, N_states
+   do j = 1, mo_tot_num
+    one_body_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_body_dm_mo_beta(j,j,istate)
+    do i = j+1, mo_tot_num
+     one_body_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_body_dm_mo_beta(i,j,istate)
+     one_body_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_body_dm_mo_beta(i,j,istate)
+    enddo
+   enddo
+  enddo
+
+ END_PROVIDER
+
  BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num,N_states) ]
 &BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num,N_states) ]
    implicit none

src/Determinants/two_body_dm_map.irp.f

@@ -194,6 +194,8 @@ subroutine add_values_to_two_body_dm_map(mask_ijkl)
 end
 
  BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_act, (n_act_orb, n_act_orb)]
+&BEGIN_PROVIDER [double precision, two_body_dm_aa_diag_act, (n_act_orb, n_act_orb)]
+&BEGIN_PROVIDER [double precision, two_body_dm_bb_diag_act, (n_act_orb, n_act_orb)]
 &BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_inact, (n_inact_orb_allocate, n_inact_orb_allocate)]
 &BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_core, (n_core_orb_allocate, n_core_orb_allocate)]
 &BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_all, (mo_tot_num, mo_tot_num)]
@@ -234,6 +236,8 @@ end
 
  two_body_dm_ab_diag_all = 0.d0
  two_body_dm_ab_diag_act = 0.d0
+ two_body_dm_aa_diag_act = 0.d0
+ two_body_dm_bb_diag_act = 0.d0
  two_body_dm_ab_diag_core = 0.d0
  two_body_dm_ab_diag_inact = 0.d0
  two_body_dm_diag_core_a_act_b = 0.d0
@@ -269,8 +273,20 @@ end
      two_body_dm_ab_diag_act(k,m) += 0.5d0 * contrib
      two_body_dm_ab_diag_act(m,k) += 0.5d0 * contrib
     enddo
+    do l = 1, n_occ_ab_act(2)
+     m = list_act_reverse(occ_act(l,2))
+     two_body_dm_bb_diag_act(k,m) += 0.5d0 * contrib
+     two_body_dm_bb_diag_act(m,k) += 0.5d0 * contrib
+    enddo
+   enddo
+   do j = 1,n_occ_ab_act(1)
+    k = list_act_reverse(occ_act(j,1))
+    do l = 1, n_occ_ab_act(1)
+     m = list_act_reverse(occ_act(l,1))
+     two_body_dm_aa_diag_act(k,m) += 0.5d0 * contrib
+     two_body_dm_aa_diag_act(m,k) += 0.5d0 * contrib
+    enddo
    enddo
-
    ! CORE PART of the diagonal part of the two body dm
    do j = 1, N_int
     key_tmp_core(j,1) = psi_det(j,1,i)
@@ -325,6 +341,8 @@ end
 END_PROVIDER 
 
  BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_act, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
+&BEGIN_PROVIDER [double precision, two_body_dm_aa_big_array_act, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
+&BEGIN_PROVIDER [double precision, two_body_dm_bb_big_array_act, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
 &BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_core_act, (n_core_orb_allocate,n_act_orb,n_act_orb)]
  implicit none
  use bitmasks
@@ -394,14 +412,22 @@ END_PROVIDER
    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
    contrib = 0.5d0 * psi_coef(i,1) * psi_coef(j,1) * phase
    if(degree==2)then  ! case of the DOUBLE EXCITATIONS  ************************************
-    if(s1==s2)cycle  ! Only the alpha/beta two body density matrix
     ! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J
     h1 = list_act_reverse(h1)
     h2 = list_act_reverse(h2)
     p1 = list_act_reverse(p1)
     p2 = list_act_reverse(p2)
+    if(s1==s2)then 
+     if(s1==1)then
+      call insert_into_two_body_dm_big_array( two_body_dm_aa_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2)
+!     call insert_into_two_body_dm_big_array( two_body_dm_aa_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,-contrib,h1,p2,h2,p1)
+     else
+      call insert_into_two_body_dm_big_array( two_body_dm_bb_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2)
+!     call insert_into_two_body_dm_big_array( two_body_dm_bb_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,-contrib,h1,p2,h2,p1)
+     endif
+    else ! alpha/beta two body density matrix
      call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2)
-
+    endif
    else if(degree==1)then! case of the SINGLE EXCITATIONS  ***************************************************
     print*,'h1 = ',h1
     h1 = list_act_reverse(h1)
@@ -417,6 +443,12 @@ END_PROVIDER
       ! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
       call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
      enddo
+     do k = 1, n_occ_ab(1)
+      m = list_act_reverse(occ(k,1))
+      ! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
+      call insert_into_two_body_dm_big_array( two_body_dm_aa_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
+!     call insert_into_two_body_dm_big_array( two_body_dm_aa_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,-contrib,h1,m,p1,m)
+     enddo
      
      ! core <-> active part of the extra diagonal two body dm 
      do k = 1, n_occ_ab_core(2)
@@ -432,6 +464,12 @@ END_PROVIDER
       ! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
       call insert_into_two_body_dm_big_array(two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
      enddo
+     do k = 1, n_occ_ab(2)
+      m = list_act_reverse(occ(k,2))
+      ! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
+      call insert_into_two_body_dm_big_array(two_body_dm_bb_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
+!     call insert_into_two_body_dm_big_array(two_body_dm_bb_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,-contrib,h1,m,p1,m)
+     enddo
      
      ! core <-> active part of the extra diagonal two body dm 
      do k = 1, n_occ_ab_core(1)
@@ -464,156 +502,3 @@ subroutine insert_into_two_body_dm_big_array(big_array,dim1,dim2,dim3,dim4,contr
 
  
 end
-
-double precision function compute_extra_diag_two_body_dm_ab(r1,r2) 
- implicit none
- BEGIN_DOC
-! compute the extra diagonal contribution to the alpha/bet two body density at r1, r2
- END_DOC 
- double precision :: r1(3), r2(3)
- double precision :: compute_extra_diag_two_body_dm_ab_act,compute_extra_diag_two_body_dm_ab_core_act
- compute_extra_diag_two_body_dm_ab = compute_extra_diag_two_body_dm_ab_act(r1,r2)+compute_extra_diag_two_body_dm_ab_core_act(r1,r2)
-end
-
-double precision function compute_extra_diag_two_body_dm_ab_act(r1,r2)
- implicit none
- BEGIN_DOC
-! compute the extra diagonal contribution to the two body density at r1, r2
-! involving ONLY THE ACTIVE PART, which means that the four index of the excitations 
-! involved in the two body density matrix are ACTIVE 
- END_DOC 
- PROVIDE n_act_orb
- double precision, intent(in) :: r1(3),r2(3)
- integer :: i,j,k,l
- double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb)
- double precision :: contrib
- double precision :: contrib_tmp
-!print*,'n_act_orb = ',n_act_orb
- compute_extra_diag_two_body_dm_ab_act = 0.d0
- call give_all_act_mos_at_r(r1,mos_array_r1)
- call give_all_act_mos_at_r(r2,mos_array_r2)
- do l = 1, n_act_orb  ! p2 
-  do k = 1, n_act_orb  ! h2 
-   do j = 1, n_act_orb  ! p1 
-    do i = 1,n_act_orb   ! h1 
-     contrib_tmp = mos_array_r1(i) * mos_array_r1(j) * mos_array_r2(k) * mos_array_r2(l)
-     compute_extra_diag_two_body_dm_ab_act += two_body_dm_ab_big_array_act(i,j,k,l) * contrib_tmp
-    enddo
-   enddo
-  enddo
- enddo
-
-end
-
-double precision function compute_extra_diag_two_body_dm_ab_core_act(r1,r2)
- implicit none
- BEGIN_DOC
-! compute the extra diagonal contribution to the two body density at r1, r2
-! involving ONLY THE ACTIVE PART, which means that the four index of the excitations 
-! involved in the two body density matrix are ACTIVE 
- END_DOC 
- double precision, intent(in) :: r1(3),r2(3)
- integer :: i,j,k,l
- double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb)
- double precision :: mos_array_core_r1(n_core_orb),mos_array_core_r2(n_core_orb)
- double precision :: contrib_core_1,contrib_core_2
- double precision :: contrib_act_1,contrib_act_2
- double precision :: contrib_tmp
- compute_extra_diag_two_body_dm_ab_core_act = 0.d0
- call give_all_act_mos_at_r(r1,mos_array_act_r1)
- call give_all_act_mos_at_r(r2,mos_array_act_r2)
- call give_all_core_mos_at_r(r1,mos_array_core_r1)
- call give_all_core_mos_at_r(r2,mos_array_core_r2)
-  do i = 1, n_act_orb  ! h1 
-   do j = 1, n_act_orb  ! p1 
-    contrib_act_1 = mos_array_act_r1(i) * mos_array_act_r1(j) 
-    contrib_act_2 = mos_array_act_r2(i) * mos_array_act_r2(j) 
-    do k = 1,n_core_orb  ! h2 
-     contrib_core_1 = mos_array_core_r1(k)  * mos_array_core_r1(k)
-     contrib_core_2 = mos_array_core_r2(k)  * mos_array_core_r2(k)
-     contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_core_2 + contrib_act_2 * contrib_core_1)
-     compute_extra_diag_two_body_dm_ab_core_act += two_body_dm_ab_big_array_core_act(k,i,j) * contrib_tmp
-    enddo
-   enddo
-  enddo
-
-end
-
-double precision function compute_diag_two_body_dm_ab_core(r1,r2)
- implicit none
- double precision :: r1(3),r2(3)
- integer :: i,j,k,l
- double precision :: mos_array_r1(n_core_orb_allocate),mos_array_r2(n_core_orb_allocate)
- double precision :: contrib,contrib_tmp
- compute_diag_two_body_dm_ab_core = 0.d0
- call give_all_core_mos_at_r(r1,mos_array_r1)
- call give_all_core_mos_at_r(r2,mos_array_r2)
- do l = 1, n_core_orb  ! 
-  contrib = mos_array_r2(l)*mos_array_r2(l)
-! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
-  do k = 1, n_core_orb  ! 
-    contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k)
-!  if(dabs(contrib).lt.threshld_two_bod_dm)cycle
-     compute_diag_two_body_dm_ab_core += two_body_dm_ab_diag_core(k,l) * contrib_tmp
-  enddo
- enddo
-
-end
-
-
-double precision function compute_diag_two_body_dm_ab_act(r1,r2)
- implicit none
- double precision :: r1(3),r2(3)
- integer :: i,j,k,l
- double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb)
- double precision :: contrib,contrib_tmp
- compute_diag_two_body_dm_ab_act = 0.d0
- call give_all_act_mos_at_r(r1,mos_array_r1)
- call give_all_act_mos_at_r(r2,mos_array_r2)
- do l = 1, n_act_orb  ! 
-  contrib = mos_array_r2(l)*mos_array_r2(l)
-! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
-  do k = 1, n_act_orb  ! 
-    contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k)
-!  if(dabs(contrib).lt.threshld_two_bod_dm)cycle
-     compute_diag_two_body_dm_ab_act += two_body_dm_ab_diag_act(k,l) * contrib_tmp
-  enddo
- enddo
-end
-
-double precision function compute_diag_two_body_dm_ab_core_act(r1,r2)
- implicit none
- double precision :: r1(3),r2(3)
- integer :: i,j,k,l
- double precision :: mos_array_core_r1(n_core_orb_allocate),mos_array_core_r2(n_core_orb_allocate)
- double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb)
- double precision :: contrib_core_1,contrib_core_2
- double precision :: contrib_act_1,contrib_act_2
- double precision :: contrib_tmp
- compute_diag_two_body_dm_ab_core_act = 0.d0
- call give_all_act_mos_at_r(r1,mos_array_act_r1)
- call give_all_act_mos_at_r(r2,mos_array_act_r2)
- call give_all_core_mos_at_r(r1,mos_array_core_r1)
- call give_all_core_mos_at_r(r2,mos_array_core_r2)
-! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
- do k = 1, n_act_orb  ! 
-    contrib_act_1 = mos_array_act_r1(k) * mos_array_act_r1(k)
-    contrib_act_2 = mos_array_act_r2(k) * mos_array_act_r2(k)
-    contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_act_2 + contrib_act_2 * contrib_act_1)
-!  if(dabs(contrib).lt.threshld_two_bod_dm)cycle
-    do l = 1, n_core_orb  ! 
-     contrib_core_1 = mos_array_core_r1(l) * mos_array_core_r1(l)
-     contrib_core_2 = mos_array_core_r2(l) * mos_array_core_r2(l)
-     compute_diag_two_body_dm_ab_core_act += two_body_dm_diag_core_act(l,k) * contrib_tmp
-    enddo
- enddo
-end
-
-double precision function compute_diag_two_body_dm_ab(r1,r2)
- implicit none
- double precision,intent(in) :: r1(3),r2(3)
- double precision :: compute_diag_two_body_dm_ab_act,compute_diag_two_body_dm_ab_core
- double precision :: compute_diag_two_body_dm_ab_core_act
- compute_diag_two_body_dm_ab  =   compute_diag_two_body_dm_ab_act(r1,r2)+compute_diag_two_body_dm_ab_core(r1,r2) & 
-                                + compute_diag_two_body_dm_ab_core_act(r1,r2)
-end