diff --git a/src/Dets/README.rst b/src/Dets/README.rst index ab562b17..bf16b740 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -50,28 +50,45 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`copy_h_apply_buffer_to_wf `_ +`copy_h_apply_buffer_to_wf `_ Undocumented -`h_apply_buffer_coef `_ - Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines. +`h_apply_buffer_coef `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. -`h_apply_buffer_det `_ - Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines. +`h_apply_buffer_det `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. -`h_apply_buffer_n_det `_ - Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines. +`h_apply_buffer_e2 `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. -`h_apply_buffer_size `_ +`h_apply_buffer_n_det `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. + +`h_apply_buffer_size `_ Size of the H_apply buffer. -`h_apply_threshold `_ +`h_apply_threshold `_ Theshold on | | -`resize_h_apply_buffer_det `_ +`resize_h_apply_buffer_det `_ Undocumented -`davidson_diag `_ +`connected_to_ref `_ + Undocumented + +`det_is_not_or_may_be_in_ref `_ + If true, det is not in ref + If false, det may be in ref + +`key_pattern_not_in_ref `_ + Min and max values of the integers of the keys of the reference + +`davidson_diag `_ Davidson diagonalization. .br dets_in : bitmasks corresponding to determinants @@ -87,102 +104,164 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`davidson_iter_max `_ +`davidson_iter_max `_ Max number of Davidson iterations -`davidson_sze_max `_ +`davidson_sze_max `_ Max number of Davidson sizes -`n_det `_ +`n_det `_ Number of determinants in the wave function -`n_det_generators `_ +`n_det_generators `_ Number of generator determinants in the wave function -`n_states `_ +`n_states `_ Number of states to consider -`psi_coef `_ +`psi_coef `_ The wave function. Initialized with Hartree-Fock -`psi_det `_ +`psi_det `_ The wave function. Initialized with Hartree-Fock -`psi_det_size `_ +`psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_generators `_ +`psi_generators `_ Determinants on which H is applied -`double_exc_bitmask `_ +`psi_ref `_ + Determinants on which H is applied + +`psi_ref_coef `_ + Determinants on which H is applied + +`psi_ref_size `_ + Number of generator determinants in the wave function + +`double_exc_bitmask `_ double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1 double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1 double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2 double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2 for a given couple of hole/particle excitations i. -`n_double_exc_bitmasks `_ +`n_double_exc_bitmasks `_ Number of double excitation bitmasks -`n_single_exc_bitmasks `_ +`n_single_exc_bitmasks `_ Number of single excitation bitmasks -`single_exc_bitmask `_ +`single_exc_bitmask `_ single_exc_bitmask(:,1,i) is the bitmask for holes single_exc_bitmask(:,2,i) is the bitmask for particles for a given couple of hole/particle excitations i. -`get_s2 `_ +`filter_connected `_ + Filters out the determinants that are not connected by H + .br + returns the array idx which contains the index of the + .br + determinants in the array key1 that interact + .br + via the H operator with key2. + .br + idx(0) is the number of determinants that interact with key1 + +`filter_connected_i_h_psi0 `_ + returns the array idx which contains the index of the + .br + determinants in the array key1 that interact + .br + via the H operator with key2. + .br + idx(0) is the number of determinants that interact with key1 + +`filter_connected_i_h_psi0_sc2 `_ + standard filter_connected_i_H_psi but returns in addition + .br + the array of the index of the non connected determinants to key1 + .br + in order to know what double excitation can be repeated on key1 + .br + idx_repeat(0) is the number of determinants that can be used + .br + to repeat the excitations + +`get_s2 `_ Returns -`a_operator `_ +`get_s2_u0 `_ + Undocumented + +`s_z `_ + Undocumented + +`s_z2_sz `_ + Undocumented + +`a_operator `_ Needed for diag_H_mat_elem -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem -`decode_exc `_ +`decode_exc `_ Decodes the exc arrays returned by get_excitation. h1,h2 : Holes p1,p2 : Particles s1,s2 : Spins (1:alpha, 2:beta) degree : Degree of excitation -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes -`get_double_excitation `_ +`get_double_excitation `_ Returns the two excitation operators between two doubly excited determinants and the phase -`get_excitation `_ +`get_excitation `_ Returns the excitation operators between two determinants and the phase -`get_excitation_degree `_ +`get_excitation_degree `_ Returns the excitation degree between two determinants -`get_excitation_degree_vector `_ +`get_excitation_degree_vector `_ Applies get_excitation_degree to an array of determinants -`get_mono_excitation `_ +`get_mono_excitation `_ Returns the excitation operator between two singly excited determinants and the phase -`get_occ_from_key `_ +`get_occ_from_key `_ Returns a list of occupation numbers from a bitstring -`h_u_0 `_ +`h_u_0 `_ Computes v_0 = H|u_0> .br n : number of determinants .br H_jj : array of -`i_h_j `_ +`i_h_j `_ Returns where i and j are determinants -`i_h_psi `_ - Undocumented +`i_h_psi `_ + for the various Nstate -`h_matrix_all_dets `_ +`i_h_psi_sc2 `_ + for the various Nstate + .br + returns in addition + .br + the array of the index of the non connected determinants to key1 + .br + in order to know what double excitation can be repeated on key1 + .br + idx_repeat(0) is the number of determinants that can be used + .br + to repeat the excitations + +`h_matrix_all_dets `_ H matrix on the basis of the slater deter;inants defined by psi_det diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index 1cb89644..b46fdaaf 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -75,3 +75,31 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] END_PROVIDER +BEGIN_PROVIDER [ integer, psi_ref_size] + implicit none + BEGIN_DOC + ! Number of generator determinants in the wave function + END_DOC + psi_det_size = N_det +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_ref_size) ] +&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_ref_size,N_states) ] + implicit none + BEGIN_DOC + ! Determinants on which H is applied + END_DOC + integer :: i,k + + do k = 1, psi_ref_size + do i=1,N_int + psi_ref(i,1,k) = psi_det(i,1,k) + psi_ref(i,2,k) = psi_det(i,1,k) + enddo + do i = 1, N_states + psi_ref_coef(k,i) = psi_coef(k,i) + enddo + enddo + +END_PROVIDER + diff --git a/src/Dets/filter_connected.irp.f b/src/Dets/filter_connected.irp.f index 0353c9c6..d6caa977 100644 --- a/src/Dets/filter_connected.irp.f +++ b/src/Dets/filter_connected.irp.f @@ -3,7 +3,16 @@ subroutine filter_connected(key1,key2,Nint,sze,idx) use bitmasks implicit none BEGIN_DOC -! Filters out the determinants that are not connected by H + ! Filters out the determinants that are not connected by H + ! + ! returns the array idx which contains the index of the + ! + ! determinants in the array key1 that interact + ! + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + END_DOC END_DOC integer, intent(in) :: Nint, sze integer(bit_kind), intent(in) :: key1(Nint,2,sze) @@ -85,6 +94,15 @@ end subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) use bitmasks + BEGIN_DOC + ! returns the array idx which contains the index of the + ! + ! determinants in the array key1 that interact + ! + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + END_DOC implicit none integer, intent(in) :: Nint, sze integer(bit_kind), intent(in) :: key1(Nint,2,sze) @@ -173,3 +191,215 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) idx(0) = l-1 end +subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat) + use bitmasks + BEGIN_DOC + ! standard filter_connected_i_H_psi but returns in addition + ! + ! the array of the index of the non connected determinants to key1 + ! + ! in order to know what double excitation can be repeated on key1 + ! + ! idx_repeat(0) is the number of determinants that can be used + ! + ! to repeat the excitations + END_DOC + implicit none + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + integer, intent(out) :: idx_repeat(0:sze) + + integer :: i,l,l_repeat + integer :: degree_x2 + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sze > 0) + + l=1 + l_repeat=1 + call get_excitation_degree(ref_bitmask,key2,degree,Nint) + integer :: degree + + if(degree == 2)then + if (Nint==1) then + + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do l=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +& + popcnt(xor( key1(l,2,i), key2(l,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + endif + elseif(degree==1)then + if (Nint==1) then + + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do l=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +& + popcnt(xor( key1(l,2,i), key2(l,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + endif + + else + print*,'more than a double excitation, can not apply the ' + print*,'SC2 dressing of the diagonal element .....' + print*,'stop !!' + print*,'degree = ',degree + stop + endif + idx(0) = l-1 + idx_repeat(0) = l_repeat-1 +end + diff --git a/src/Dets/slater_rules.irp.f b/src/Dets/slater_rules.irp.f index 3f56eba3..3599b9a0 100644 --- a/src/Dets/slater_rules.irp.f +++ b/src/Dets/slater_rules.irp.f @@ -502,6 +502,9 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) integer :: exc(0:2,2,2) double precision :: hij integer :: idx(0:Ndet) + BEGIN_DOC + ! for the various Nstate + END_DOC ASSERT (Nint > 0) ASSERT (N_int == Nint) @@ -517,7 +520,55 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) do j = 1, Nstate i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij enddo - print *, 'x', coef(i,1), hij, i_H_psi_array(1) +! print *, 'x', coef(i,1), hij, i_H_psi_array(1) + enddo +end + + +subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat) + use bitmasks + BEGIN_DOC + ! for the various Nstate + ! + ! returns in addition + ! + ! the array of the index of the non connected determinants to key1 + ! + ! in order to know what double excitation can be repeated on key1 + ! + ! idx_repeat(0) is the number of determinants that can be used + ! + ! to repeat the excitations + END_DOC + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + integer , intent(out) :: idx_repeat(0:Ndet) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer :: idx(0:Ndet) + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + i_H_psi_array = 0.d0 + call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat) + do ii=1,idx(0) + i = idx(ii) + !DEC$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo +! print *, 'x', coef(i,1), hij, i_H_psi_array(1) enddo end diff --git a/src/Makefile.config.example b/src/Makefile.config.example new file mode 100644 index 00000000..f4cbb568 --- /dev/null +++ b/src/Makefile.config.example @@ -0,0 +1,29 @@ +OPENMP =1 +PROFILE =0 +DEBUG = 0 + +IRPF90_FLAGS+= --align=32 +FC = ifort -g +FCFLAGS= +FCFLAGS+= -xHost +#FCFLAGS+= -xAVX +FCFLAGS+= -O2 +FCFLAGS+= -ip +FCFLAGS+= -opt-prefetch +FCFLAGS+= -ftz +MKL=-mkl=parallel + +ifeq ($(PROFILE),1) +FC += -p -g +CXX += -pg +endif + +ifeq ($(OPENMP),1) +FC += -openmp +CXX += -fopenmp +endif + +ifeq ($(DEBUG),1) +FC += -C -traceback -fpe0 +#FCFLAGS =-O0 +endif diff --git a/src/Perturbation/Moller_plesset.irp.f b/src/Perturbation/Moller_plesset.irp.f new file mode 100644 index 00000000..4f99e9d2 --- /dev/null +++ b/src/Perturbation/Moller_plesset.irp.f @@ -0,0 +1,41 @@ +subroutine pt2_moller_plesset(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 Moller-Plesset perturbative first order coefficient and second order energetic contribution + ! + ! for the various n_st states. + ! + ! c_pert(i) = /(difference of orbital energies) + ! + ! e_2_pert(i) = ^2/(difference of orbital energies) + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase,delta_e + integer :: h1,h2,p1,p2,s1,s2 + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call get_excitation(det_pert,ref_bitmask,exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - & + Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2) + delta_e = 1.d0/delta_e + + 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) *delta_e + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + enddo + +end diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst index 83392b83..3acd8ef6 100644 --- a/src/Perturbation/README.rst +++ b/src/Perturbation/README.rst @@ -82,7 +82,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`pt2_epstein_nesbet `_ +`pt2_epstein_nesbet `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various n_st states. @@ -92,7 +92,7 @@ Documentation e_2_pert(i) = ^2/( E(i) - ) .br -`pt2_epstein_nesbet_2x2 `_ +`pt2_epstein_nesbet_2x2 `_ compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution .br for the various n_st states. @@ -102,6 +102,76 @@ Documentation c_pert(i) = e_2_pert(i)/ .br +`pt2_epstein_nesbet_2x2_sc2 `_ + compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution + .br + for the various n_st states. + .br + but with the correction in the denominator + .br + comming from the interaction of that determinant with all the others determinants + .br + that can be repeated by repeating all the double excitations + .br + : you repeat all the correlation energy already taken into account in reference_energy(1) + .br + that could be repeated to this determinant. + .br + ---> + delta_e_corr + .br + e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) + .br + c_pert(i) = e_2_pert(i)/ + .br + +`pt2_epstein_nesbet_sc2 `_ + compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + .br + for the various n_st states, + .br + but with the correction in the denominator + .br + comming from the interaction of that determinant with all the others determinants + .br + that can be repeated by repeating all the double excitations + .br + : you repeat all the correlation energy already taken into account in reference_energy(1) + .br + that could be repeated to this determinant. + .br + ---> + delta_e_corr + .br + c_pert(i) = /( E(i) - ( ) ) + .br + e_2_pert(i) = ^2/( E(i) - ( ) ) + .br + +`pt2_epstein_nesbet_sc2_projected `_ + compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + .br + for the various n_st states, + .br + but with the correction in the denominator + .br + comming from the interaction of that determinant with all the others determinants + .br + that can be repeated by repeating all the double excitations + .br + : you repeat all the correlation energy already taken into account in reference_energy(1) + .br + that could be repeated to this determinant. + .br + BUT on the contrary with ""pt2_epstein_nesbet_SC2"", you compute the energy by projection + .br + ---> + delta_e_corr + .br + c_pert(1) = 1/c_HF /( E(i) - ( ) ) + .br + e_2_pert(1) = c_pert(1) + .br + NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !! + .br + Needed Modules diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f index 4c1e7c45..37572ae5 100644 --- a/src/Perturbation/epstein_nesbet.irp.f +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -66,3 +66,164 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet print *, e_2_pert, delta_e , i_H_psi_array end + +subroutine pt2_epstein_nesbet_SC2(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) + integer :: idx_repeat(ndet) + + BEGIN_DOC + ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various n_st states, + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in reference_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! ---> + delta_e_corr + ! + ! c_pert(i) = /( E(i) - ( ) ) + ! + ! e_2_pert(i) = ^2/( E(i) - ( ) ) + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem,accu_e_corr,hij + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call i_H_psi_SC2(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + do i = 1, idx_repeat(0) + call i_H_j(psi_ref(1,1,idx_repeat(i)),det_pert,Nint,hij) + accu_e_corr = accu_e_corr + hij * psi_ref_coef(idx_repeat(i)) + enddo + accu_e_corr = accu_e_corr / psi_ref_coef(1) + H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + do i =1,n_st + c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag) + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + enddo + +end +subroutine pt2_epstein_nesbet_2x2_SC2(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) + integer :: idx_repeat(ndet) + + BEGIN_DOC + ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution + ! + ! for the various n_st states. + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in reference_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! ---> + delta_e_corr + ! + ! 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,accu_e_corr,hij,delta_e + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call i_H_psi_SC2(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + do i = 1, idx_repeat(0) + call i_H_j(psi_ref(1,1,idx_repeat(i)),det_pert,Nint,hij) + accu_e_corr = accu_e_corr + hij * psi_ref_coef(idx_repeat(i)) + enddo + accu_e_corr = accu_e_corr / psi_ref_coef(1) + H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + do i =1,n_st + delta_e = H_pert_diag - reference_energy(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 + +subroutine pt2_epstein_nesbet_SC2_projected(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) + integer :: idx_repeat(ndet) + + BEGIN_DOC + ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various n_st states, + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in reference_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! BUT on the contrary with ""pt2_epstein_nesbet_SC2"", you compute the energy by projection + ! + ! ---> + delta_e_corr + ! + ! c_pert(1) = 1/c_HF /( E(i) - ( ) ) + ! + ! e_2_pert(1) = c_pert(1) + ! + ! NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !! + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call i_H_psi_SC2(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + call i_H_j(ref_bitmask,det_pert,Nint,h0j) + do i = 1, idx_repeat(0) + call i_H_j(psi_ref(1,1,idx_repeat(i)),det_pert,Nint,hij) + accu_e_corr = accu_e_corr + hij * psi_ref_coef(idx_repeat(i)) + enddo + accu_e_corr = accu_e_corr / psi_ref_coef(1) + H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + + c_pert(1) = 1.d0/psi_ref_coef(1) * i_H_psi_array(1) / (reference_energy(i) - H_pert_diag) + e_2_pert(1) = c_pert(i) * h0j + do i =2,n_st + c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag) + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + enddo + +end