diff --git a/src/Perturbation/ASSUMPTIONS.rst b/src/Perturbation/ASSUMPTIONS.rst new file mode 100644 index 00000000..fd22e8a3 --- /dev/null +++ b/src/Perturbation/ASSUMPTIONS.rst @@ -0,0 +1,6 @@ +* This is not allowed: + + subroutine & + pt2_.... + + diff --git a/src/Perturbation/Makefile b/src/Perturbation/Makefile new file mode 100644 index 00000000..ea95cbf1 --- /dev/null +++ b/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 diff --git a/src/Perturbation/NEEDED_MODULES b/src/Perturbation/NEEDED_MODULES new file mode 100644 index 00000000..f6551dd9 --- /dev/null +++ b/src/Perturbation/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f new file mode 100644 index 00000000..7f8474b0 --- /dev/null +++ b/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) = /( E(i) - ) + ! + ! e_2_pert(i) = ^2/( E(i) - ) + ! + 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 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) + ! + ! c_pert(i) = e_2_pert(i)/ + ! + 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 diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f new file mode 100644 index 00000000..5e672529 --- /dev/null +++ b/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 +