Formally finished to code the MRPT_Utils

plugins/MRPT_Utils/MRPT_Utils.main.irp.f

@@ -0,0 +1,60 @@
+program MRPT_Utils
+  implicit none
+  read_wf = .True.
+  touch read_wf
+! call routine
+! call routine_2
+  call routine_3
+end
+
+
+subroutine routine_3
+ implicit none
+ provide one_creation
+
+end
+
+subroutine routine_2
+ implicit none
+ integer :: i
+ do i = 1, n_core_inact_orb
+  print*,fock_core_inactive_total(i,1),fock_core_inactive(i)
+ enddo
+ double precision :: accu
+ accu = 0.d0
+ do i = 1, n_act_orb
+  integer :: j_act_orb
+  j_act_orb = list_act(i)
+  accu += one_body_dm_mo_alpha(j_act_orb,j_act_orb)
+  print*,one_body_dm_mo_alpha(j_act_orb,j_act_orb),one_body_dm_mo_beta(j_act_orb,j_act_orb) 
+ enddo
+ print*,'accu = ',accu
+
+end
+
+subroutine routine
+ implicit none
+ integer :: i,j
+ integer :: orb, spin_exc
+ integer :: hole_particle
+ double precision, allocatable :: norm_out(:)
+ allocate(norm_out(N_states_diag))
+
+ orb = list_virt(10)
+ hole_particle = -1
+ spin_exc = 1 
+
+ call apply_exc_to_psi(orb,hole_particle,spin_exc, & 
+           norm_out,psi_det,psi_coef, n_det,psi_det_size,psi_det_size,N_states_diag)
+ do i = 1, N_det
+  if(psi_coef(i,1).ne.0.d0)then
+   print*, ''
+   call debug_det(psi_det(1,1,i),N_int)
+  print*, 'coef = ',psi_coef(i,1)
+  endif
+ enddo
+ print*,'norm_out = ',norm_out
+ 
+ deallocate(norm_out)
+
+end

plugins/MRPT_Utils/energies_cas.irp.f

@@ -0,0 +1,457 @@
+BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)]
+ implicit none
+ integer :: i 
+ double precision :: energies(N_states_diag)
+ do i = 1, N_states
+  call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i)
+  energy_cas_dyall(i) = energies(i)
+ enddo
+END_PROVIDER 
+
+
+BEGIN_PROVIDER [ double precision, one_creation, (n_act_orb,2)]
+ implicit none
+ integer :: i,j
+ integer :: ispin
+ integer :: orb, hole_particle,spin_exc
+ double precision  :: norm_out(N_states_diag)
+ integer(bit_kind) :: psi_in_out(N_int,2,n_det)
+ double precision :: psi_in_out_coef(n_det,N_states_diag)
+ use bitmasks
+
+ integer :: iorb
+ do iorb = 1,n_act_orb
+  do ispin = 1,2
+   orb = list_act(iorb)
+   hole_particle = 1
+   spin_exc = ispin 
+   do i = 1, n_det
+    do j = 1, n_states_diag
+      psi_in_out_coef(i,j) = psi_coef(i,j)
+    enddo
+    do j = 1, N_int
+     psi_in_out(j,1,i) =  psi_active(j,1,i) 
+     psi_in_out(j,1,i) =  psi_active(j,2,i) 
+    enddo
+   enddo
+   integer :: state_target
+   state_target = 1
+   double precision :: energies(n_states_diag)
+   call apply_exc_to_psi(orb,hole_particle,spin_exc, & 
+           norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+   call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
+   one_creation(iorb,ispin) = energy_cas_dyall(state_target)  - energies(state_target)
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, one_anhilation, (n_act_orb,2)]
+ implicit none
+ integer :: i,j
+ integer :: ispin
+ integer :: orb, hole_particle,spin_exc
+ double precision  :: norm_out(N_states_diag)
+ integer(bit_kind) :: psi_in_out(N_int,2,n_det)
+ double precision :: psi_in_out_coef(n_det,N_states_diag)
+ use bitmasks
+
+ integer :: iorb
+ do iorb = 1,n_act_orb
+  do ispin = 1,2
+   orb = list_act(iorb)
+   hole_particle = -1
+   spin_exc = ispin 
+   do i = 1, n_det
+    do j = 1, n_states_diag
+      psi_in_out_coef(i,j) = psi_coef(i,j)
+    enddo
+    do j = 1, N_int
+     psi_in_out(j,1,i) =  psi_active(j,1,i) 
+     psi_in_out(j,1,i) =  psi_active(j,2,i) 
+    enddo
+   enddo
+   integer :: state_target
+   state_target = 1
+   double precision :: energies(n_states_diag)
+   call apply_exc_to_psi(orb,hole_particle,spin_exc, & 
+           norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+   call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
+   one_anhilation(iorb,ispin) = energy_cas_dyall(state_target)  -  energies(state_target)
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, two_creation, (n_act_orb,n_act_orb,2,2)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin
+ integer :: orb_i, hole_particle_i,spin_exc_i
+ integer :: orb_j, hole_particle_j,spin_exc_j
+ double precision  :: norm_out(N_states_diag)
+ integer(bit_kind) :: psi_in_out(N_int,2,n_det)
+ double precision :: psi_in_out_coef(n_det,N_states_diag)
+ use bitmasks
+
+ integer :: iorb,jorb
+ integer :: state_target
+ state_target = 1
+ double precision :: energies(n_states_diag)
+ do iorb = 1,n_act_orb
+  do ispin = 1,2
+   orb_i = list_act(iorb)
+   hole_particle_i = 1
+   spin_exc_i = ispin 
+   do jorb = 1, n_act_orb
+    do jspin = 1,2
+     orb_j = list_act(jorb)
+     hole_particle_j = 1
+     spin_exc_j = jspin 
+     do i = 1, n_det
+      do j = 1, n_states_diag
+        psi_in_out_coef(i,j) = psi_coef(i,j)
+      enddo
+      do j = 1, N_int
+       psi_in_out(j,1,i) =  psi_active(j,1,i) 
+       psi_in_out(j,1,i) =  psi_active(j,2,i) 
+      enddo
+     enddo
+     call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & 
+             norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+     call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & 
+             norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+     call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
+     two_creation(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target)  -   energies(state_target)
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, two_anhilation, (n_act_orb,n_act_orb,2,2)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin
+ integer :: orb_i, hole_particle_i,spin_exc_i
+ integer :: orb_j, hole_particle_j,spin_exc_j
+ double precision  :: norm_out(N_states_diag)
+ integer(bit_kind) :: psi_in_out(N_int,2,n_det)
+ double precision :: psi_in_out_coef(n_det,N_states_diag)
+ use bitmasks
+
+ integer :: iorb,jorb
+ integer :: state_target
+ state_target = 1
+ double precision :: energies(n_states_diag)
+ do iorb = 1,n_act_orb
+  do ispin = 1,2
+   orb_i = list_act(iorb)
+   hole_particle_i = -1
+   spin_exc_i = ispin 
+   do jorb = 1, n_act_orb
+    do jspin = 1,2
+     orb_j = list_act(jorb)
+     hole_particle_j = -1
+     spin_exc_j = jspin 
+     do i = 1, n_det
+      do j = 1, n_states_diag
+        psi_in_out_coef(i,j) = psi_coef(i,j)
+      enddo
+      do j = 1, N_int
+       psi_in_out(j,1,i) =  psi_active(j,1,i) 
+       psi_in_out(j,1,i) =  psi_active(j,2,i) 
+      enddo
+     enddo
+     call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & 
+             norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+     call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & 
+             norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+     call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
+     two_anhilation(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target)  -   energies(state_target)
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, one_anhilation_one_creation, (n_act_orb,n_act_orb,2,2)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin
+ integer :: orb_i, hole_particle_i,spin_exc_i
+ integer :: orb_j, hole_particle_j,spin_exc_j
+ double precision  :: norm_out(N_states_diag)
+ integer(bit_kind) :: psi_in_out(N_int,2,n_det)
+ double precision :: psi_in_out_coef(n_det,N_states_diag)
+ use bitmasks
+
+ integer :: iorb,jorb
+ integer :: state_target
+ state_target = 1
+ double precision :: energies(n_states_diag)
+ do iorb = 1,n_act_orb
+  do ispin = 1,2
+   orb_i = list_act(iorb)
+   hole_particle_i =  1
+   spin_exc_i = ispin 
+   do jorb = 1, n_act_orb
+    do jspin = 1,2
+     orb_j = list_act(jorb)
+     hole_particle_j = -1
+     spin_exc_j = jspin 
+     do i = 1, n_det
+      do j = 1, n_states_diag
+        psi_in_out_coef(i,j) = psi_coef(i,j)
+      enddo
+      do j = 1, N_int
+       psi_in_out(j,1,i) =  psi_active(j,1,i) 
+       psi_in_out(j,1,i) =  psi_active(j,2,i) 
+      enddo
+     enddo
+     call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & 
+             norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+     call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & 
+             norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+     call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
+     one_anhilation_one_creation(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target)  -   energies(state_target)
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, two_anhilation_one_creation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin,kspin
+ integer :: orb_i, hole_particle_i,spin_exc_i
+ integer :: orb_j, hole_particle_j,spin_exc_j
+ integer :: orb_k, hole_particle_k,spin_exc_k
+ double precision  :: norm_out(N_states_diag)
+ integer(bit_kind) :: psi_in_out(N_int,2,n_det)
+ double precision :: psi_in_out_coef(n_det,N_states_diag)
+ use bitmasks
+
+ integer :: iorb,jorb
+ integer :: korb
+ integer :: state_target
+ state_target = 1
+ double precision :: energies(n_states_diag)
+ do iorb = 1,n_act_orb
+  do ispin = 1,2
+   orb_i = list_act(iorb)
+   hole_particle_i =  1
+   spin_exc_i = ispin 
+   do jorb = 1, n_act_orb
+    do jspin = 1,2
+     orb_j = list_act(jorb)
+     hole_particle_j = -1
+     spin_exc_j = jspin 
+     do korb = 1, n_act_orb
+      do kspin = 1,2
+       orb_k = list_act(korb)
+       hole_particle_k = -1
+       spin_exc_k = kspin 
+       do i = 1, n_det
+        do j = 1, n_states_diag
+          psi_in_out_coef(i,j) = psi_coef(i,j)
+        enddo
+        do j = 1, N_int
+         psi_in_out(j,1,i) =  psi_active(j,1,i) 
+         psi_in_out(j,1,i) =  psi_active(j,2,i) 
+        enddo
+       enddo
+       call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
+       two_anhilation_one_creation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target)  -  energies(state_target)
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, two_creation_one_anhilation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin,kspin
+ integer :: orb_i, hole_particle_i,spin_exc_i
+ integer :: orb_j, hole_particle_j,spin_exc_j
+ integer :: orb_k, hole_particle_k,spin_exc_k
+ double precision  :: norm_out(N_states_diag)
+ integer(bit_kind) :: psi_in_out(N_int,2,n_det)
+ double precision :: psi_in_out_coef(n_det,N_states_diag)
+ use bitmasks
+
+ integer :: iorb,jorb
+ integer :: korb
+ integer :: state_target
+ state_target = 1
+ double precision :: energies(n_states_diag)
+ do iorb = 1,n_act_orb
+  do ispin = 1,2
+   orb_i = list_act(iorb)
+   hole_particle_i =  1
+   spin_exc_i = ispin 
+   do jorb = 1, n_act_orb
+    do jspin = 1,2
+     orb_j = list_act(jorb)
+     hole_particle_j =  1
+     spin_exc_j = jspin 
+     do korb = 1, n_act_orb
+      do kspin = 1,2
+       orb_k = list_act(korb)
+       hole_particle_k = -1
+       spin_exc_k = kspin 
+       do i = 1, n_det
+        do j = 1, n_states_diag
+          psi_in_out_coef(i,j) = psi_coef(i,j)
+        enddo
+        do j = 1, N_int
+         psi_in_out(j,1,i) =  psi_active(j,1,i) 
+         psi_in_out(j,1,i) =  psi_active(j,2,i) 
+        enddo
+       enddo
+       call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
+       two_creation_one_anhilation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target)  -  energies(state_target)
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, three_creation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin,kspin
+ integer :: orb_i, hole_particle_i,spin_exc_i
+ integer :: orb_j, hole_particle_j,spin_exc_j
+ integer :: orb_k, hole_particle_k,spin_exc_k
+ double precision  :: norm_out(N_states_diag)
+ integer(bit_kind) :: psi_in_out(N_int,2,n_det)
+ double precision :: psi_in_out_coef(n_det,N_states_diag)
+ use bitmasks
+
+ integer :: iorb,jorb
+ integer :: korb
+ integer :: state_target
+ state_target = 1
+ double precision :: energies(n_states_diag)
+ do iorb = 1,n_act_orb
+  do ispin = 1,2
+   orb_i = list_act(iorb)
+   hole_particle_i =  1
+   spin_exc_i = ispin 
+   do jorb = 1, n_act_orb
+    do jspin = 1,2
+     orb_j = list_act(jorb)
+     hole_particle_j =  1
+     spin_exc_j = jspin 
+     do korb = 1, n_act_orb
+      do kspin = 1,2
+       orb_k = list_act(korb)
+       hole_particle_k =  1
+       spin_exc_k = kspin 
+       do i = 1, n_det
+        do j = 1, n_states_diag
+          psi_in_out_coef(i,j) = psi_coef(i,j)
+        enddo
+        do j = 1, N_int
+         psi_in_out(j,1,i) =  psi_active(j,1,i) 
+         psi_in_out(j,1,i) =  psi_active(j,2,i) 
+        enddo
+       enddo
+       call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
+       three_creation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target)  -  energies(state_target)
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, three_anhilation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin,kspin
+ integer :: orb_i, hole_particle_i,spin_exc_i
+ integer :: orb_j, hole_particle_j,spin_exc_j
+ integer :: orb_k, hole_particle_k,spin_exc_k
+ double precision  :: norm_out(N_states_diag)
+ integer(bit_kind) :: psi_in_out(N_int,2,n_det)
+ double precision :: psi_in_out_coef(n_det,N_states_diag)
+ use bitmasks
+
+ integer :: iorb,jorb
+ integer :: korb
+ integer :: state_target
+ state_target = 1
+ double precision :: energies(n_states_diag)
+ do iorb = 1,n_act_orb
+  do ispin = 1,2
+   orb_i = list_act(iorb)
+   hole_particle_i = -1
+   spin_exc_i = ispin 
+   do jorb = 1, n_act_orb
+    do jspin = 1,2
+     orb_j = list_act(jorb)
+     hole_particle_j = -1
+     spin_exc_j = jspin 
+     do korb = 1, n_act_orb
+      do kspin = 1,2
+       orb_k = list_act(korb)
+       hole_particle_k = -1
+       spin_exc_k = kspin 
+       do i = 1, n_det
+        do j = 1, n_states_diag
+          psi_in_out_coef(i,j) = psi_coef(i,j)
+        enddo
+        do j = 1, N_int
+         psi_in_out(j,1,i) =  psi_active(j,1,i) 
+         psi_in_out(j,1,i) =  psi_active(j,2,i) 
+        enddo
+       enddo
+       call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & 
+               norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
+       call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
+       three_anhilation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target)  -  energies(state_target)
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER

plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f

@@ -0,0 +1,190 @@
+
+BEGIN_PROVIDER [ double precision, one_creation_spin_averaged, (n_act_orb)]
+ implicit none
+ integer :: i
+ do i = 1, n_act_orb
+  one_creation_spin_averaged(i) = one_creation(i,1) + one_creation(i,2)
+  one_creation_spin_averaged(i) = 0.5d0 * one_creation_spin_averaged(i)
+ enddo
+END_PROVIDER 
+
+
+BEGIN_PROVIDER [ double precision, one_anhilation_spin_averaged, (n_act_orb)]
+ implicit none
+ integer :: i
+ do i = 1, n_act_orb
+  one_anhilation_spin_averaged(i) = one_anhilation(i,1) + one_anhilation(i,2)
+  one_anhilation_spin_averaged(i) = 0.5d0 * one_anhilation_spin_averaged(i)
+ enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, two_creation_spin_averaged, (n_act_orb,n_act_orb)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin
+ double precision :: counting
+ do i = 1, n_act_orb
+  do j = 1, n_act_orb 
+   two_creation_spin_averaged(j,i) = 0.d0
+   counting = 0.d0
+   do ispin = 1, 2
+    do jspin = 1,2 
+     two_creation_spin_averaged(j,i) += two_creation(j,i,ispin,jspin)
+     counting += 1.d0 
+    enddo
+   enddo 
+   two_creation_spin_averaged(j,i) = two_creation_spin_averaged(j,i) / counting
+  enddo
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, two_anhilation_spin_averaged, (n_act_orb,n_act_orb)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin
+ double precision :: counting
+ do i = 1, n_act_orb
+  do j = 1, n_act_orb 
+   two_anhilation_spin_averaged(j,i) = 0.d0
+   counting = 0.d0
+   do ispin = 1, 2
+    do jspin = 1,2 
+     two_anhilation_spin_averaged(j,i) += two_anhilation(j,i,ispin,jspin)
+     counting += 1.d0 
+    enddo
+   enddo 
+   two_anhilation_spin_averaged(j,i) = two_anhilation_spin_averaged(j,i) / counting
+  enddo
+ enddo
+END_PROVIDER 
+
+
+BEGIN_PROVIDER [ double precision, one_anhilation_one_creation_spin_averaged, (n_act_orb,n_act_orb)]
+ implicit none
+ integer :: i,j
+ integer :: ispin,jspin
+ double precision :: counting
+ do i = 1, n_act_orb
+  do j = 1, n_act_orb 
+   one_anhilation_one_creation_spin_averaged(j,i) = 0.d0
+   counting = 0.d0
+   do ispin = 1, 2
+    do jspin = 1,2 
+     one_anhilation_one_creation_spin_averaged(j,i) += one_anhilation_one_creation(j,i,jspin,ispin)
+     counting += 1.d0 
+    enddo
+   enddo
+   one_anhilation_one_creation_spin_averaged(j,i) = one_anhilation_one_creation_spin_averaged(j,i) / counting
+  enddo
+ enddo
+  
+END_PROVIDER
+ 
+
+BEGIN_PROVIDER [ double precision, two_anhilation_one_creation_spin_averaged, (n_act_orb,n_act_orb,n_act_orb)]
+ implicit none
+ integer :: i,j,k
+ integer :: ispin,jspin,kspin
+ double precision :: counting
+ 
+ do i = 1, n_act_orb 
+  do j = 1, n_act_orb 
+   do k = 1, n_act_orb
+    two_anhilation_one_creation_spin_averaged(k,j,i) = 0.d0
+    counting = 0.d0
+    do ispin = 1, 2 
+     do jspin = 1,2
+      do kspin = 1,2
+       two_anhilation_one_creation_spin_averaged(k,j,i)  += two_anhilation_one_creation(k,j,i,kspin,jspin,ispin)
+       counting += 1.d0
+      enddo 
+     enddo
+     two_anhilation_one_creation_spin_averaged(k,j,i) = two_anhilation_one_creation_spin_averaged(k,j,i) / counting
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, two_creation_one_anhilation_spin_averaged, (n_act_orb,n_act_orb,n_act_orb)]
+ implicit none
+ integer :: i,j,k
+ integer :: ispin,jspin,kspin
+ double precision :: counting
+ 
+ do i = 1, n_act_orb 
+  do j = 1, n_act_orb 
+   do k = 1, n_act_orb
+    two_creation_one_anhilation_spin_averaged(k,j,i) = 0.d0
+    counting = 0.d0
+    do ispin = 1, 2 
+     do jspin = 1,2
+      do kspin = 1,2
+       two_creation_one_anhilation_spin_averaged(k,j,i)  += two_creation_one_anhilation(k,j,i,kspin,jspin,ispin)
+       counting += 1.d0
+      enddo 
+     enddo
+     two_creation_one_anhilation_spin_averaged(k,j,i) = two_creation_one_anhilation_spin_averaged(k,j,i) / counting
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
+
+BEGIN_PROVIDER [ double precision, three_creation_spin_averaged, (n_act_orb,n_act_orb,n_act_orb)]
+ implicit none
+ integer :: i,j,k
+ integer :: ispin,jspin,kspin
+ double precision :: counting
+ 
+ do i = 1, n_act_orb 
+  do j = 1, n_act_orb 
+   do k = 1, n_act_orb
+    three_creation_spin_averaged(k,j,i) = 0.d0
+    counting = 0.d0
+    do ispin = 1, 2 
+     do jspin = 1,2
+      do kspin = 1,2
+       three_creation_spin_averaged(k,j,i)  += three_creation(k,j,i,kspin,jspin,ispin)
+       counting += 1.d0
+      enddo 
+     enddo
+     three_creation_spin_averaged(k,j,i) = three_creation_spin_averaged(k,j,i) / counting
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
+
+BEGIN_PROVIDER [ double precision, three_anhilation_spin_averaged, (n_act_orb,n_act_orb,n_act_orb)]
+ implicit none
+ integer :: i,j,k
+ integer :: ispin,jspin,kspin
+ double precision :: counting
+ 
+ do i = 1, n_act_orb 
+  do j = 1, n_act_orb 
+   do k = 1, n_act_orb
+    three_anhilation_spin_averaged(k,j,i) = 0.d0
+    counting = 0.d0
+    do ispin = 1, 2 
+     do jspin = 1,2
+      do kspin = 1,2
+       three_anhilation_spin_averaged(k,j,i)  += three_anhilation(k,j,i,kspin,jspin,ispin)
+       counting += 1.d0
+      enddo 
+     enddo
+     three_anhilation_spin_averaged(k,j,i) = three_anhilation_spin_averaged(k,j,i) / counting
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+

plugins/MRPT_Utils/fock_like_operators.irp.f

@@ -0,0 +1,170 @@
+ BEGIN_PROVIDER [double precision, fock_core_inactive, (mo_tot_num)]
+ BEGIN_DOC
+! inactive part of the fock operator with contributions only from the inactive
+ END_DOC
+ implicit none
+ integer :: i,j
+ double precision :: accu
+
+ integer :: j_inact_core_orb,i_inact_core_orb
+ do i = 1, n_core_inact_orb
+  accu = 0.d0
+  i_inact_core_orb = list_core_inact(i)
+  do j = 1, n_core_inact_orb
+! do j = 1, elec_alpha_num
+!  j_inact_core_orb = j
+   j_inact_core_orb = list_core_inact(j)
+   accu += 2.d0 * mo_bielec_integral_jj(i_inact_core_orb,j_inact_core_orb)  & 
+                - mo_bielec_integral_jj_exchange(i_inact_core_orb,j_inact_core_orb)
+  enddo
+  fock_core_inactive(i_inact_core_orb) = accu + mo_mono_elec_integral(i_inact_core_orb,i_inact_core_orb)
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, fock_virt_from_core_inact, (mo_tot_num)]
+ BEGIN_DOC
+! fock operator for the virtuals that comes from the doubly occupied orbitals
+ END_DOC
+ implicit none
+ integer :: i,j
+ double precision :: accu
+
+ integer :: j_inact_core_orb,i_virt_orb
+ do i = 1, n_virt_orb
+  accu = 0.d0
+  i_virt_orb = list_virt(i)
+  do j = 1, n_core_inact_orb
+! do j = 1, elec_alpha_num
+!  j_inact_core_orb = j
+   j_inact_core_orb = list_core_inact(j)
+   accu += 2.d0 * mo_bielec_integral_jj(i_virt_orb,j_inact_core_orb)  & 
+                - mo_bielec_integral_jj_exchange(i_virt_orb,j_inact_core_orb)
+  enddo
+  fock_virt_from_core_inact(i_virt_orb) = accu + mo_mono_elec_integral(i_virt_orb,i_virt_orb)
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, fock_core_inactive_from_act, (mo_tot_num,2)]
+ BEGIN_DOC
+! inactive part of the fock operator with contributions only from the active
+ END_DOC
+ implicit none
+ integer :: i,j
+ double precision :: accu_coulomb,accu_exchange(2)
+ double precision :: na,nb,ntot
+ double precision :: coulomb, exchange
+ integer :: j_act_orb,i_inact_core_orb
+
+ do i = 1, n_core_inact_orb
+  accu_coulomb = 0.d0
+  accu_exchange = 0.d0
+  i_inact_core_orb = list_core_inact(i)
+  do j = 1, n_act_orb
+   j_act_orb = list_act(j)
+   na = one_body_dm_mo_alpha(j_act_orb,j_act_orb)
+   nb = one_body_dm_mo_beta(j_act_orb,j_act_orb)
+   ntot = na + nb
+   coulomb = mo_bielec_integral_jj(i_inact_core_orb,j_act_orb)
+   exchange = mo_bielec_integral_jj_exchange(i_inact_core_orb,j_act_orb)
+   accu_coulomb += ntot * coulomb
+   accu_exchange(1) += na * exchange
+   accu_exchange(2) += nb * exchange
+  enddo
+  fock_core_inactive_from_act(i_inact_core_orb,1) = accu_coulomb + accu_exchange(1) 
+  fock_core_inactive_from_act(i_inact_core_orb,2) = accu_coulomb + accu_exchange(2) 
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, fock_virt_from_act, (mo_tot_num,2)]
+ BEGIN_DOC
+! virtual part of the fock operator with contributions only from the active
+ END_DOC
+ implicit none
+ integer :: i,j
+ double precision :: accu_coulomb,accu_exchange(2)
+ double precision :: na,nb,ntot
+ double precision :: coulomb, exchange
+ integer :: j_act_orb,i_virt_orb
+
+ do i = 1, n_virt_orb
+  accu_coulomb = 0.d0
+  accu_exchange = 0.d0
+  i_virt_orb = list_virt(i)
+  do j = 1, n_act_orb
+   j_act_orb = list_act(j)
+   na = one_body_dm_mo_alpha(j_act_orb,j_act_orb)
+   nb = one_body_dm_mo_beta(j_act_orb,j_act_orb)
+   ntot = na + nb
+   coulomb = mo_bielec_integral_jj(i_virt_orb,j_act_orb)
+   exchange = mo_bielec_integral_jj_exchange(i_virt_orb,j_act_orb)
+   accu_coulomb += ntot * coulomb
+   accu_exchange(1) += na * exchange
+   accu_exchange(2) += nb * exchange
+  enddo
+  fock_virt_from_act(i_virt_orb,1) = accu_coulomb + accu_exchange(1) 
+  fock_virt_from_act(i_virt_orb,2) = accu_coulomb + accu_exchange(2) 
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, fock_core_inactive_total, (mo_tot_num,2)]
+&BEGIN_PROVIDER [double precision, fock_core_inactive_total_spin_averaged, (mo_tot_num)]
+ BEGIN_DOC
+! inactive part of the fock operator 
+ END_DOC
+ implicit none
+ integer :: i
+ integer :: i_inact_core_orb
+ do i = 1, n_core_inact_orb
+  i_inact_core_orb = list_core_inact(i)
+  fock_core_inactive_total(i_inact_core_orb,1) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,1)
+  fock_core_inactive_total(i_inact_core_orb,2) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,2)
+  fock_core_inactive_total_spin_averaged(i_inact_core_orb) = 0.5d0 * (fock_core_inactive_total(i_inact_core_orb,1) + fock_core_inactive_total(i_inact_core_orb,2))
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, fock_virt_total, (mo_tot_num,2)]
+&BEGIN_PROVIDER [double precision, fock_virt_total_spin_averaged, (mo_tot_num)]
+ BEGIN_DOC
+! inactive part of the fock operator 
+ END_DOC
+ implicit none
+ integer :: i
+ integer :: i_virt_orb
+ do i = 1, n_virt_orb
+  i_virt_orb= list_virt(i)
+  fock_virt_total(i_virt_orb,1) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,1)
+  fock_virt_total(i_virt_orb,2) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,2)
+  fock_virt_total_spin_averaged(i_virt_orb) = 0.5d0 * ( fock_virt_total(i_virt_orb,1) + fock_virt_total(i_virt_orb,2) )
+ enddo
+ END_PROVIDER
+
+
+
+
+
+ BEGIN_PROVIDER [double precision, fock_operator_active_from_core_inact, (n_act_orb,n_act_orb)]
+ BEGIN_DOC
+! active part of the fock operator with contributions only from the inactive
+ END_DOC
+ implicit none
+ integer :: i,j,k
+ double precision :: accu
+ double precision               :: get_mo_bielec_integral,coulomb, exchange
+ PROVIDE mo_bielec_integrals_in_map
+ do i = 1, n_act_orb
+  do j = 1, n_act_orb
+   accu = 0.d0
+   do k = 1, n_core_inact_orb 
+    coulomb  = get_mo_bielec_integral(k,i,k,j,mo_integrals_map)
+    exchange = get_mo_bielec_integral(k,i,i,k,mo_integrals_map)
+    accu += 2.d0 * coulomb - exchange
+   enddo
+   fock_operator_active_from_core_inact(i,j) = accu
+  enddo
+ enddo
+
+ END_PROVIDER 
+
+ 
+
+

plugins/MRPT_Utils/psi_active_prov.irp.f

@@ -0,0 +1,279 @@
+
+  use bitmasks
+BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)]
+ BEGIN_DOC
+! active part of psi
+ END_DOC
+ implicit none
+  use bitmasks
+ integer :: i,j,k,l
+ do i = 1, N_det
+  do j = 1, N_int
+   psi_active(j,1,i) = iand(psi_det(j,1,i),cas_bitmask(j,1,1))
+   psi_active(j,2,i) = iand(psi_det(j,2,i),cas_bitmask(j,1,1))
+  enddo
+ enddo
+END_PROVIDER
+
+
+subroutine give_holes_and_particles_in_active_space(det_1,det_2,n_holes_spin,n_particles_spin,n_holes,n_particles,& 
+                                                    holes_active_list,particles_active_list)
+ implicit none
+ use bitmasks
+ integer(bit_kind),intent(in)  :: det_1(N_int,2)
+ integer(bit_kind),intent(in ) :: det_2(N_int,2)
+ integer, intent(out) :: n_holes_spin(2),n_particles_spin(2)
+ integer, intent(out) :: n_holes,n_particles
+ integer, intent(out) :: holes_active_list(2 * n_act_orb,2)
+ integer, intent(out) :: particles_active_list(2 * n_act_orb,2)
+ integer :: i
+ integer(bit_kind) :: holes(N_int,2)
+ integer(bit_kind) :: particles(N_int,2)
+ integer(bit_kind) :: det_tmp_2(N_int,2),det_tmp_1(N_int,2)
+ BEGIN_DOC
+! returns the holes and particles operators WITHIN THE ACTIVE SPACE 
+! that connect det_1 and det_2. By definition, the holes/particles 
+! are such that one starts from det_1 and goes to det_2
+!        
+! n_holes is the total number of holes
+! n_particles is the total number of particles
+! n_holes_spin is the number of number of holes per spin (1=alpha, 2=beta)
+! n_particles_spin is the number of number of particles per spin (1=alpha, 2=beta)
+! holes_active_list is the index of the holes per spin, that ranges from 1 to n_act_orb
+! particles_active_list is the index of the particles per spin, that ranges from 1 to n_act_orb
+ END_DOC
+ 
+ call give_active_part_determinant(det_1,det_tmp_1)
+ call give_active_part_determinant(det_2,det_tmp_2)
+ do i = 1, N_int
+  holes(i,1) = iand(det_tmp_1(i,1),xor(det_tmp_1(i,1),det_tmp_2(i,1)))
+  holes(i,2) = iand(det_tmp_1(i,2),xor(det_tmp_1(i,2),det_tmp_2(i,2)))
+  particles(i,1) = iand(det_tmp_2(i,1),xor(det_tmp_1(i,1),det_tmp_2(i,1)))
+  particles(i,2) = iand(det_tmp_2(i,2),xor(det_tmp_1(i,2),det_tmp_2(i,2)))
+ enddo
+
+ integer :: holes_list(N_int*bit_kind_size,2)
+ holes_list = 0
+ call bitstring_to_list(holes(1,1), holes_list(1,1), n_holes_spin(1), N_int)
+ call bitstring_to_list(holes(1,2), holes_list(1,2), n_holes_spin(2), N_int)
+
+ n_holes = 0
+ do i = 1, n_holes_spin(1)
+  n_holes +=1
+  holes_active_list(i,1) = list_act_reverse(holes_list(i,1))
+ enddo
+ do i = 1, n_holes_spin(2)
+  n_holes +=1
+  holes_active_list(i,2) = list_act_reverse(holes_list(i,2))
+ enddo
+
+
+ integer :: particles_list(N_int*bit_kind_size,2)
+ particles_list = 0
+ call bitstring_to_list(particles(1,1), particles_list(1,1), n_particles_spin(1), N_int)
+ call bitstring_to_list(particles(1,2), particles_list(1,2), n_particles_spin(2), N_int)
+ n_particles = 0
+ do i = 1, n_particles_spin(1)
+  n_particles += 1
+  particles_active_list(i,1) = list_act_reverse(particles_list(i,1))
+ enddo
+ do i = 1, n_particles_spin(2)
+  n_particles += 1
+  particles_active_list(i,2) = list_act_reverse(particles_list(i,2))
+ enddo
+
+end
+
+subroutine give_holes_in_inactive_space(det_1,n_holes_spin,n_holes,holes_list)
+ BEGIN_DOC
+! returns the holes operators WITHIN THE INACTIVE SPACE 
+! that has lead to det_1.  
+!        
+! n_holes is the total number of holes
+! n_holes_spin is the number of number of holes per spin (1=alpha, 2=beta)
+! holes_inactive_list is the index of the holes per spin, that ranges from 1 to mo_tot_num
+ END_DOC
+ implicit none
+ use bitmasks
+ integer(bit_kind),intent(in)  :: det_1(N_int,2)
+ integer, intent(out) :: n_holes_spin(2)
+ integer, intent(out) :: n_holes
+ integer, intent(out) :: holes_list(N_int*bit_kind_size,2)
+ integer :: i
+ integer(bit_kind) :: holes(N_int,2)
+ integer(bit_kind) :: det_tmp_1(N_int,2)
+ 
+ call give_core_inactive_part_determinant(det_1,det_tmp_1)
+ 
+ do i = 1, N_int
+  holes(i,1) = iand(reunion_of_core_inact_bitmask(i,1),xor(det_tmp_1(i,1),reunion_of_core_inact_bitmask(i,1)))
+  holes(i,2) = iand(reunion_of_core_inact_bitmask(i,2),xor(det_tmp_1(i,2),reunion_of_core_inact_bitmask(i,2)))
+ enddo
+ holes_list = 0
+ call bitstring_to_list(holes(1,1), holes_list(1,1), n_holes_spin(1), N_int)
+ call bitstring_to_list(holes(1,2), holes_list(1,2), n_holes_spin(2), N_int)
+ n_holes =  n_holes_spin(1) +  n_holes_spin(2)
+
+end
+
+subroutine give_particles_in_virt_space(det_1,n_particles_spin,n_particles,particles_list)
+ BEGIN_DOC
+! returns the holes operators WITHIN THE VIRTUAL SPACE 
+! that has lead to det_1.  
+!        
+! n_particles is the total number of particles
+! n_particles_spin is the number of number of particles per spin (1=alpha, 2=beta)
+! particles_inactive_list is the index of the particles per spin, that ranges from 1 to mo_tot_num
+ END_DOC
+ implicit none
+ use bitmasks
+ integer(bit_kind),intent(in)  :: det_1(N_int,2)
+ integer, intent(out) :: n_particles_spin(2)
+ integer, intent(out) :: n_particles
+ integer, intent(out) :: particles_list(N_int*bit_kind_size,2)
+ integer :: i
+ integer(bit_kind) :: det_tmp_1(N_int,2)
+ integer(bit_kind) :: particles(N_int,2)
+
+ call give_virt_part_determinant(det_1,det_tmp_1)
+
+ do i = 1, N_int
+  particles(i,1) = iand(virt_bitmask(i,1),xor(det_tmp_1(i,1),virt_bitmask(i,1)))
+  particles(i,2) = iand(virt_bitmask(i,2),xor(det_tmp_1(i,2),virt_bitmask(i,2)))
+ enddo
+
+ particles_list = 0
+ call bitstring_to_list(particles(1,1), particles_list(1,1), n_particles_spin(1), N_int)
+ call bitstring_to_list(particles(1,2), particles_list(1,2), n_particles_spin(2), N_int)
+ n_particles = n_particles_spin(1) + n_particles_spin(2)
+ 
+
+end
+
+subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
+ implicit none
+  use bitmasks
+ double precision, intent(out) :: delta_e_final
+ integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
+ integer :: i,j,k,l
+ 
+ integer :: n_holes_spin(2)
+ integer :: n_holes
+ integer :: holes_list(N_int*bit_kind_size,2)
+
+
+ double precision :: delta_e_inactive
+ integer :: i_hole_inact
+
+ call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list)
+ delta_e_inactive = 0.d0
+ do i = 1, n_holes_spin(1)
+  i_hole_inact = holes_list(i,1)
+  delta_e_inactive += fock_core_inactive_total_spin_averaged(i_hole_inact)
+ enddo
+
+ do i = 1, n_holes_spin(2)
+  i_hole_inact = holes_list(i,2)
+  delta_e_inactive += fock_core_inactive_total_spin_averaged(i_hole_inact)
+ enddo
+
+ double precision :: delta_e_virt
+ integer :: i_part_virt
+ integer :: n_particles_spin(2)
+ integer :: n_particles
+ integer :: particles_list(N_int*bit_kind_size,2)
+ 
+ call give_particles_in_virt_space(det_2,n_particles_spin,n_particles,particles_list)
+ delta_e_virt = 0.d0
+ do i = 1, n_particles_spin(1)
+  i_part_virt = particles_list(i,1)
+  delta_e_virt += fock_virt_total_spin_averaged(i_part_virt)
+ enddo
+
+
+ integer :: n_holes_spin_act(2),n_particles_spin_act(2)
+ integer :: n_holes_act,n_particles_act
+ integer :: holes_active_list(2*n_act_orb,2)
+ integer :: holes_active_list_spin_traced(4*n_act_orb)
+ integer :: particles_active_list(2*n_act_orb,2)
+ integer :: particles_active_list_spin_traced(4*n_act_orb)
+ double precision :: delta_e_act
+ delta_e_act = 0.d0
+ call give_holes_and_particles_in_active_space(det_1,det_2,n_holes_spin_act,n_particles_spin_act, &
+                                               n_holes_act,n_particles_act,holes_active_list,particles_active_list)
+ integer :: icount
+ icount = 0
+ do i = 1, n_holes_spin_act(1)
+  icount += 1
+  holes_active_list_spin_traced(icount) = holes_active_list(i,1)
+ enddo
+ do i = 1, n_holes_spin_act(2)
+  icount += 1
+  holes_active_list_spin_traced(icount) = holes_active_list(i,2)
+ enddo
+ if(icount .ne. n_holes) then
+  print * , 'pb in holes_active_list_spin_traced !!' 
+  stop
+ endif
+
+ icount = 0
+ do i = 1, n_particles_spin_act(1)
+  icount += 1
+  particles_active_list_spin_traced(icount) = particles_active_list(i,1)
+ enddo
+ do i = 1, n_particles_spin_act(2)
+  icount += 1
+  particles_active_list_spin_traced(icount) = particles_active_list(i,2)
+ enddo
+ if(icount .ne. n_particles) then
+  print * , 'pb in particles_active_list_spin_traced !!' 
+  stop
+ endif
+
+
+ integer :: i_hole_act, j_hole_act, k_hole_act
+ integer :: i_particle_act, j_particle_act, k_particle_act
+ 
+
+ if      (n_holes_act == 1 .and. n_particles_act == 0) then
+  i_hole_act =  holes_active_list_spin_traced(1)
+  delta_e_act += one_creation_spin_averaged(i_hole_act)
+
+ else if (n_holes_act == 0 .and. n_particles_act == 1) then
+  i_particle_act =  particles_active_list_spin_traced(1)
+  delta_e_act += one_anhilation_spin_averaged(i_particle_act)
+
+ else if (n_holes_act == 1 .and. n_particles_act == 1) then
+  i_hole_act =  holes_active_list_spin_traced(1)
+  i_particle_act =  particles_active_list_spin_traced(1)
+  delta_e_act += one_anhilation_one_creation_spin_averaged(i_hole_act,i_particle_act)
+
+ else if (n_holes_act == 2 .and. n_particles_act == 1) then
+  i_hole_act =  holes_active_list_spin_traced(1)
+  j_hole_act =  holes_active_list_spin_traced(2)
+  i_particle_act =  particles_active_list_spin_traced(1)
+  delta_e_act += two_anhilation_one_creation_spin_averaged(i_hole_act,j_hole_act,i_particle_act)
+
+ else if (n_holes_act == 1 .and. n_particles_act == 2) then
+  i_hole_act =  holes_active_list_spin_traced(1)
+  i_particle_act =  particles_active_list_spin_traced(1)
+  j_particle_act =  particles_active_list_spin_traced(2)
+  delta_e_act += two_creation_one_anhilation_spin_averaged(i_hole_act,i_particle_act,j_particle_act)
+
+ else if (n_holes_act == 3 .and. n_particles_act == 0) then
+  i_hole_act =  holes_active_list_spin_traced(1)
+  j_hole_act =  holes_active_list_spin_traced(2)
+  k_hole_act =  holes_active_list_spin_traced(3)
+  delta_e_act += three_anhilation_spin_averaged(i_hole_act,j_hole_act,k_hole_act)
+
+ else if (n_holes_act == 0 .and. n_particles_act == 3) then
+  i_particle_act =  particles_active_list_spin_traced(1)
+  j_particle_act =  particles_active_list_spin_traced(2)
+  k_particle_act =  particles_active_list_spin_traced(3)
+  delta_e_act += three_creation_spin_averaged(i_particle_act,j_particle_act,k_particle_act)
+
+ endif
+
+ delta_e_final = delta_e_act + delta_e_virt + delta_e_inactive 
+
+end

plugins/MRPT_Utils/utils_bitmask.irp.f

@@ -0,0 +1,36 @@
+
+subroutine give_active_part_determinant(det_in,det_out)
+ implicit none
+ use bitmasks
+ integer(bit_kind),intent(in)  :: det_in(N_int,2)
+ integer(bit_kind),intent(out) :: det_out(N_int,2)
+ integer :: i
+ do i = 1,N_int
+  det_out(i,1) = iand(det_in(i,1),cas_bitmask(i,1,1))
+  det_out(i,2) = iand(det_in(i,2),cas_bitmask(i,1,1))
+ enddo
+end
+
+subroutine give_core_inactive_part_determinant(det_in,det_out)
+ implicit none
+ use bitmasks
+ integer(bit_kind),intent(in)  :: det_in(N_int,2)
+ integer(bit_kind),intent(out) :: det_out(N_int,2)
+ integer :: i
+ do i = 1,N_int
+  det_out(i,1) = iand(det_in(i,1),reunion_of_core_inact_bitmask(i,1))
+  det_out(i,2) = iand(det_in(i,2),reunion_of_core_inact_bitmask(i,1))
+ enddo
+end
+
+subroutine give_virt_part_determinant(det_in,det_out)
+ implicit none
+ use bitmasks
+ integer(bit_kind),intent(in)  :: det_in(N_int,2)
+ integer(bit_kind),intent(out) :: det_out(N_int,2)
+ integer :: i
+ do i = 1,N_int
+  det_out(i,1) = iand(det_in(i,1),virt_bitmask(i,1))
+  det_out(i,2) = iand(det_in(i,2),virt_bitmask(i,1))
+ enddo
+end

src/Bitmask/bitmasks.irp.f

@@ -431,12 +431,40 @@ END_PROVIDER
 
  END_PROVIDER 
 
+
+ BEGIN_PROVIDER [ integer, list_core_inact, (n_core_inact_orb)]
+&BEGIN_PROVIDER [ integer, list_core_inact_reverse, (mo_tot_num)]
+
+ implicit none
+ integer :: occ_inact(N_int*bit_kind_size)
+ integer :: itest,i
+ occ_inact = 0
+
+ call bitstring_to_list(reunion_of_core_inact_bitmask(1,1), occ_inact(1), itest, N_int)
+
+ list_core_inact_reverse = 0
+ do i = 1, n_core_inact_orb
+  list_core_inact(i) = occ_inact(i)
+  list_core_inact_reverse(occ_inact(i)) = i
+ enddo
+
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [ integer, n_core_inact_orb ]
+ implicit none
+ integer :: i
+ n_core_inact_orb = 0
+ do i = 1, N_int
+  n_core_inact_orb += popcnt(reunion_of_core_inact_bitmask(i,1))
+ enddo
+ ENd_PROVIDER
+
  BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)]
  implicit none
  BEGIN_DOC
  ! Reunion of the core and inactive and virtual bitmasks
  END_DOC
- integer :: i,j
+ integer :: i
  do i = 1, N_int
   reunion_of_core_inact_bitmask(i,1) = ior(core_bitmask(i,1),inact_bitmask(i,1))
   reunion_of_core_inact_bitmask(i,2) = ior(core_bitmask(i,2),inact_bitmask(i,2))
@@ -474,6 +502,7 @@ END_PROVIDER
 
 
  BEGIN_PROVIDER [ integer(bit_kind), inact_virt_bitmask, (N_int,2)]
+&BEGIN_PROVIDER [ integer(bit_kind), core_inact_virt_bitmask, (N_int,2)]
  implicit none
  BEGIN_DOC
  ! Reunion of the inactive and virtual bitmasks
@@ -482,6 +511,8 @@ END_PROVIDER
  do i = 1, N_int
   inact_virt_bitmask(i,1) = ior(inact_bitmask(i,1),virt_bitmask(i,1))
   inact_virt_bitmask(i,2) = ior(inact_bitmask(i,2),virt_bitmask(i,2)) 
+  core_inact_virt_bitmask(i,1) = ior(core_bitmask(i,1),inact_virt_bitmask(i,1))
+  core_inact_virt_bitmask(i,2) = ior(core_bitmask(i,2),inact_virt_bitmask(i,2)) 
  enddo
  END_PROVIDER
 

src/Determinants/create_excitations.irp.f

@@ -45,3 +45,16 @@ subroutine set_bit_to_integer(i_physical,key,Nint)
  j = i_physical-ishft(k-1,bit_kind_shift)-1
  key(k) = ibset(key(k),j)
 end
+
+
+subroutine clear_bit_to_integer(i_physical,key,Nint)
+ use bitmasks
+ implicit none
+ integer, intent(in) :: i_physical,Nint
+ integer(bit_kind), intent(inout) :: key(Nint)
+ integer :: k,j,i
+ k = ishft(i_physical-1,-bit_kind_shift)+1
+ j = i_physical-ishft(k-1,bit_kind_shift)-1
+ key(k) = ibclr(key(k),j)
+end
+