Perturbation OK. Need to add tests

src/Perturbation/ASSUMPTIONS.rst

@@ -0,0 +1,6 @@
+* This is not allowed:
+
+  subroutine &
+    pt2_....
+
+

src/Perturbation/Makefile

@@ -0,0 +1,8 @@
+default: all
+
+# Define here all new external source files and objects.Don't forget to prefix the
+# object files with IRPF90_temp/
+SRC=perturbation_template.f
+OBJ=
+
+include $(QPACKAGE_ROOT)/src/Makefile.common

src/Perturbation/NEEDED_MODULES

@@ -0,0 +1 @@
+AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils

src/Perturbation/epstein_nesbet.irp.f

@@ -0,0 +1,64 @@
+subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: Nint,ndet,n_st
+  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
+  double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
+  double precision               :: i_H_psi_array(N_st)
+  
+  BEGIN_DOC
+  ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various n_st states.
+  !
+  ! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
+  !
+  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array)
+  H_pert_diag = diag_H_mat_elem(det_pert,Nint)
+  do i =1,n_st
+    c_pert(i) = i_H_psi_array(i) / (E_ref(i) - H_pert_diag)
+    e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
+  enddo
+  
+end
+
+subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: Nint,ndet,n_st
+  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
+  double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
+  double precision               :: i_H_psi_array(N_st)
+  
+  BEGIN_DOC
+  ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
+  !
+  ! for the various n_st states.
+  !
+  ! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> -  E(i) )  - sqrt( ( <det_pert|H|det_pert> -  E(i)) ^2 + 4 <psi(i)|H|det_pert>^2  )
+  !
+  ! c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem,delta_e
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array)
+  H_pert_diag = diag_H_mat_elem(det_pert,Nint)
+  do i =1,n_st
+    delta_e = H_pert_diag - E_ref(i)
+    e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
+    c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
+  enddo
+  
+end

src/Perturbation/perturbation_template.f

@@ -0,0 +1,44 @@
+BEGIN_SHELL [ /usr/bin/env python ]
+import perturbation
+END_SHELL
+
+subroutine perturb_buffer_$PERT(buffer,buffer_size,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint)
+  implicit none
+  BEGIN_DOC
+  !  Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
+  ! routine.
+  END_DOC
+  
+  integer, intent(in)            :: Nint, N_st, buffer_size
+  integer(bit_kind), intent(in)  :: buffer(Nint,2,buffer_size)
+  double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st),sum_H_pert_diag(N_st)
+  double precision               :: c_pert(N_st), e_2_pert(N_st), H_pert_diag
+  integer                        :: i,k, c_ref
+  integer                        :: connected_to_ref
+  
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+  ASSERT (buffer_size >= 0)
+  ASSERT (minval(sum_norm_pert) >= 0.d0)
+  ASSERT (N_st > 0)
+  do i = 1,buffer_size
+
+    c_ref = connected_to_ref(buffer(1,1,i),psi_det,Nint,N_det_ref,N_det,h_apply_threshold)
+
+    if (c_ref /= 0) then
+      cycle
+    endif
+    
+    call pt2_$PERT(buffer(1,1,i),         &
+        c_pert,e_2_pert,H_pert_diag,Nint,n_det_ref,n_st)
+
+    do k = 1,N_st
+      sum_norm_pert(k) += c_pert(k) * c_pert(k)
+      sum_e_2_pert(k) += e_2_pert(k)
+      sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag
+    enddo
+    
+  enddo
+  
+end
+