pert 2 rdm works for a single determinant

src/casscf/50.casscf.bats

@@ -0,0 +1,49 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
+
+
+function run_stoch() {
+  thresh=$2
+  test_exe casscf || skip
+  qp set perturbation do_pt2 True
+  qp set determinants n_det_max $3
+  qp set davidson threshold_davidson 1.e-10
+  qp set davidson n_states_diag 4
+  qp run casscf | tee casscf.out 
+  energy1="$(ezfio get casscf energy_pt2 | tr '[]' ' ' | cut -d ',' -f 1)"
+  eq $energy1 $1 $thresh
+}
+
+@test "F2" { # 18.0198s
+  rm -rf f2_casscf
+  qp_create_ezfio -b aug-cc-pvdz ../input/f2.zmt  -o f2_casscf
+  qp set_file f2_casscf
+  qp run scf 
+  qp set_mo_class --core="[1-6,8-9]" --act="[7,10]" --virt="[11-46]"
+  run_stoch  -198.773366970 1.e-4  100000
+}
+
+@test "N2" { # 18.0198s
+  rm -rf n2_casscf
+  qp_create_ezfio -b aug-cc-pvdz ../input/n2.xyz  -o n2_casscf
+  qp set_file n2_casscf
+  qp run scf 
+  qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]"
+  run_stoch  -109.0961643162 1.e-4  100000
+}
+
+@test "N2_stretched" {
+  rm -rf n2_stretched_casscf
+  qp_create_ezfio -b aug-cc-pvdz -m 7 ../input/n2_stretched.xyz  -o n2_stretched_casscf
+  qp set_file n2_stretched_casscf 
+  qp run scf  | tee scf.out 
+  qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]"
+  qp set electrons elec_alpha_num 7 
+  qp set electrons elec_beta_num 7 
+  run_stoch  -108.7860471300 1.e-4  100000
+#   
+
+}
+

src/casscf/test_pert_2rdm.irp.f

@@ -2,7 +2,13 @@ program test_pert_2rdm
  implicit none
  read_wf = .True.
  touch read_wf
-!call get_pert_2rdm
+
+!provide is_pert_2rdm_provided
+ call test
+ end
+
+ subroutine test
+ implicit none
  integer :: i,j,k,l,ii,jj,kk,ll
  double precision :: accu , get_two_e_integral, integral
  accu = 0.d0
@@ -25,5 +31,12 @@ program test_pert_2rdm
    enddo
   enddo
  enddo
+ do ii = 1, n_orb_pert_rdm
+  i = list_orb_pert_rdm(ii) 
+  do jj = 1, n_orb_pert_rdm
+   j = list_orb_pert_rdm(jj) 
+   accu +=  pert_1rdm_provider(j,i) * mo_one_e_integrals(j,i)
+  enddo
+ enddo
  print*,'accu = ',accu
 end

src/cipsi/pert_rdm_providers.irp.f

@@ -10,7 +10,23 @@ END_PROVIDER
 
 BEGIN_PROVIDER [logical , pert_2rdm ]
  implicit none
- pert_2rdm = .False.
+ pert_2rdm = .True.
+END_PROVIDER 
+
+BEGIN_PROVIDER [logical, is_pert_2rdm_provided ]
+ implicit none
+ is_pert_2rdm_provided = .True.
+ provide pert_2rdm_provider
+ if(.True.)then
+
+  double precision :: pt2(N_states),relative_error, error(N_states),variance(N_states),norm(N_states)
+  relative_error = 0.d0
+  pert_2rdm_provider = 0.d0
+  call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, &                                         
+          norm,0) ! Stochastic PT2
+  print*,'is_pert_2rdm_provided = ',is_pert_2rdm_provided
+  print*,'pt2 = ',pt2
+ endif
 END_PROVIDER 
 
 BEGIN_PROVIDER [integer, n_orb_pert_rdm]
@@ -30,9 +46,24 @@ BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)]
 
 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)]
+BEGIN_PROVIDER [double precision, pert_1rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm)]
+ implicit none
+ integer :: i,j,k,l
+ if(is_pert_2rdm_provided)then
+  pert_1rdm_provider = 0.d0
+  do i = 1, n_orb_pert_rdm
+   do j = 1, n_orb_pert_rdm
+    do k = 1, n_orb_pert_rdm
+     pert_1rdm_provider(j,i) += 2.d0 * pert_2rdm_provider(i,k,j,k)
+    enddo
+   enddo
+  enddo
+ endif
+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)]
+&BEGIN_PROVIDER [double precision, pert_1rdm_provider_bis, (n_orb_pert_rdm,n_orb_pert_rdm)]
  implicit none
- pert_2rdm_provider = 0.d0
 
 END_PROVIDER
 

src/cipsi/update_2rdm.irp.f

@@ -28,12 +28,12 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio
    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)
+  ! 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)
@@ -69,6 +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
  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