diff --git a/plugins/MRPT_Utils/EZFIO.cfg b/plugins/MRPT_Utils/EZFIO.cfg index 2fcc26ad..948aa735 100644 --- a/plugins/MRPT_Utils/EZFIO.cfg +++ b/plugins/MRPT_Utils/EZFIO.cfg @@ -5,3 +5,10 @@ interface: ezfio,provider,ocaml default: True +[pure_state_specific_mrpt2] +type: logical +doc: If true, diagonalize the dressed matrix for each state and do a state following of the initial states +interface: ezfio,provider,ocaml +default: True + + diff --git a/plugins/MRPT_Utils/MRPT_Utils.main.irp.f b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f index 13c8228a..72750f8c 100644 --- a/plugins/MRPT_Utils/MRPT_Utils.main.irp.f +++ b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f @@ -10,16 +10,20 @@ end subroutine routine_3 implicit none + integer :: i !provide fock_virt_total_spin_trace provide delta_ij print *, 'N_det = ', N_det print *, 'N_states = ', N_states - print *, 'PT2 = ', second_order_pt_new(1) - print *, 'E = ', CI_energy(1) - print *, 'E+PT2 = ', CI_energy(1)+second_order_pt_new(1) - print *,'****** DIAGONALIZATION OF DRESSED MATRIX ******' - print *, 'E dressed= ', CI_dressed_pt2_new_energy(1) + do i = 1, N_States + print*,'State',i + write(*,'(A12,X,I3,A3,XX,F16.10)') ' PT2 ', i,' = ', second_order_pt_new(i) + write(*,'(A12,X,I3,A3,XX,F16.09)') ' E ', i,' = ', CI_energy(i) + write(*,'(A12,X,I3,A3,XX,F16.09)') ' E+PT2 ', i,' = ', CI_energy(i)+second_order_pt_new(i) + write(*,'(A12,X,I3,A3,XX,F16.09)') ' E dressed ', i,' = ', CI_dressed_pt2_new_energy(i) + write(*,'(A12,X,I3,A3,XX,F16.09)') ' S^2 ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i) + enddo end diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index f09c30cb..8f29717c 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -241,13 +241,13 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) 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) - if(orb_i == orb_j .and. ispin .ne. jspin)then + !if(orb_i == orb_j .and. ispin .ne. jspin)then call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) - else - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) - one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) - endif + !else + ! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + ! one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + !endif enddo enddo enddo @@ -527,7 +527,7 @@ END_PROVIDER double precision :: thresh_norm - thresh_norm = 1.d-10 + thresh_norm = 1.d-20 @@ -552,14 +552,7 @@ END_PROVIDER call debug_det(psi_in_out,N_int) print*, 'pb, i_ok ne 0 !!!' endif -! call i_H_j_no_k_operators_from_act(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij_test) call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij) -! if(i==1.and.dabs(hij)>1.d-8)then -! if(dabs(hij)>1.d-8)then -! print*, ispin,vorb,iorb -! print*, i,hij,hij_test -! pause -! endif do j = 1, n_states double precision :: coef,contrib coef = psi_coef(i,j) !* psi_coef(i,j) @@ -635,7 +628,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta double precision :: thresh_norm - thresh_norm = 1.d-10 + thresh_norm = 1.d-20 do aorb = 1,n_act_orb orb_a = list_act(aorb) diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index 275af0e4..60bb2b69 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -84,7 +84,11 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip do i_state = 1, N_states coef_array(i_state) = psi_coef(index_i,i_state) enddo - call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) + if(dabs(hialpha).le.1.d-10)then + delta_e = 1.d+20 + else + call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) + endif hij_array(index_i) = hialpha call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int) ! phase_array(index_i) = phase diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 4e8bc7d0..8ac8e3e0 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -121,8 +121,8 @@ ! 1h2p delta_ij_tmp = 0.d0 -!call give_1h2p_contrib(delta_ij_tmp) - call H_apply_mrpt_1h2p(delta_ij_tmp,N_det) + call give_1h2p_contrib(delta_ij_tmp) +!!!!call H_apply_mrpt_1h2p(delta_ij_tmp,N_det) accu = 0.d0 do i_state = 1, N_states do i = 1, N_det @@ -137,8 +137,8 @@ ! 2h1p delta_ij_tmp = 0.d0 -!call give_2h1p_contrib(delta_ij_tmp) - call H_apply_mrpt_2h1p(delta_ij_tmp,N_det) + call give_2h1p_contrib(delta_ij_tmp) +!!!! call H_apply_mrpt_2h1p(delta_ij_tmp,N_det) accu = 0.d0 do i_state = 1, N_states do i = 1, N_det @@ -152,19 +152,19 @@ print*, '2h1p = ',accu ! 2h2p -!delta_ij_tmp = 0.d0 -!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det) -!accu = 0.d0 -!do i_state = 1, N_states -!do i = 1, N_det -! do j = 1, N_det -! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) -! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) -! enddo -!enddo -!second_order_pt_new_2h2p(i_state) = accu(i_state) -!enddo -!print*, '2h2p = ',accu + delta_ij_tmp = 0.d0 +!!!!!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det) + accu = 0.d0 + do i_state = 1, N_states + do i = 1, N_det + do j = 1, N_det + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + enddo + enddo + second_order_pt_new_2h2p(i_state) = accu(i_state) + enddo + print*, '2h2p = ',accu double precision :: contrib_2h2p(N_states) call give_2h2p(contrib_2h2p) @@ -236,13 +236,15 @@ END_PROVIDER logical, allocatable :: good_state_array(:) double precision, allocatable :: s2_values_tmp(:) integer :: i_other_state - double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:), hmatrix_tmp(:,:) integer :: i_state double precision :: s2,e_0 integer :: i,j,k double precision, allocatable :: s2_eigvalues(:) double precision, allocatable :: e_array(:) integer, allocatable :: iorder(:) + double precision :: overlap(N_det) + double precision, allocatable :: psi_tmp(:) ! Guess values for the "N_states_diag" states of the CI_dressed_pt2_new_eigenvectors do j=1,min(N_states,N_det) @@ -265,82 +267,131 @@ END_PROVIDER else if (diag_algorithm == "Lapack") then allocate (eigenvectors(N_det,N_det)) allocate (eigenvalues(N_det)) - call lapack_diag(eigenvalues,eigenvectors, & - Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det,N_det) - CI_electronic_dressed_pt2_new_energy(:) = 0.d0 - if (s2_eig) then - i_state = 0 - allocate (s2_eigvalues(N_det)) - allocate(index_good_state_array(N_det),good_state_array(N_det)) - good_state_array = .False. - call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,& - N_det,size(eigenvectors,1)) - do j=1,N_det - ! Select at least n_states states with S^2 values closed to "expected_s2" - if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then - i_state += 1 - index_good_state_array(i_state) = j - good_state_array(j) = .True. - endif - if (i_state==N_states) then - exit - endif + if(pure_state_specific_mrpt2)then + allocate (hmatrix_tmp(N_det,N_det)) + allocate (iorder(N_det)) + allocate (psi_tmp(N_det)) + print*,'' + print*,'***************************' + do i_state = 1, N_states !! Big loop over states + print*,'' + print*,'Diagonalizing with the dressing for state',i_state + do i = 1, N_det + do j = 1, N_det + hmatrix_tmp(j,i) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) + enddo enddo - if (i_state /= 0) then - ! Fill the first "i_state" states that have a correct S^2 value - do j = 1, i_state - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) - enddo - CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(index_good_state_array(j)) - CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) - enddo - i_other_state = 0 - do j = 1, N_det - if(good_state_array(j))cycle - i_other_state +=1 - if(i_state+i_other_state.gt.n_states_diag)then - exit - endif - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) - enddo - CI_electronic_dressed_pt2_new_energy(i_state+i_other_state) = eigenvalues(j) - CI_dressed_pt2_new_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) - enddo - - else - print*,'' - print*,'!!!!!!!! WARNING !!!!!!!!!' - print*,' Within the ',N_det,'determinants selected' - print*,' and the ',N_states_diag,'states requested' - print*,' We did not find any state with S^2 values close to ',expected_s2 - print*,' We will then set the first N_states eigenvectors of the H matrix' - print*,' as the CI_dressed_pt2_new_eigenvectors' - print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' - print*,'' - do j=1,min(N_states_diag,N_det) - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) - enddo - CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) - CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j) - enddo - endif - deallocate(index_good_state_array,good_state_array) - deallocate(s2_eigvalues) - else - call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,& - min(N_det,N_states_diag),size(eigenvectors,1)) - ! Select the "N_states_diag" states of lowest energy - do j=1,min(N_det,N_states_diag) - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) - enddo - CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) + call lapack_diag(eigenvalues,eigenvectors, & + Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det,N_det) + write(*,'(A86)')'Looking for the most overlapping state within all eigenvectors of the dressed matrix' + print*,'' + print*,'Calculating the overlap for ...' + do i = 1, N_det + overlap(i) = 0.d0 + iorder(i) = i + print*,'eigenvector',i + do j = 1, N_det + overlap(i)+= psi_coef(j,i_state) * eigenvectors(j,i) + enddo + overlap(i) = -dabs(overlap(i)) + print*,'energy = ',eigenvalues(i) + nuclear_repulsion + print*,'overlap = ',dabs(overlap(i)) enddo + print*,'' + print*,'Sorting the eigenvectors per overlap' + call dsort(overlap,iorder,n_states) + print*,'' + print*,'The most overlapping state is the ',iorder(1) + print*,'with the overlap of ',dabs(overlap(1)) + print*,'and an energy of ',eigenvalues(iorder(1)) + nuclear_repulsion + print*,'Calculating the S^2 value ...' + do i=1,N_det + CI_dressed_pt2_new_eigenvectors(i,i_state) = eigenvectors(i,iorder(1)) + psi_tmp(i) = eigenvectors(i,iorder(1)) + enddo + CI_electronic_dressed_pt2_new_energy(i_state) = eigenvalues(iorder(1)) + call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2(i_state),psi_tmp,N_det,psi_det,N_int,1,N_det) + print*,'S^2 = ', CI_dressed_pt2_new_eigenvectors_s2(i_state) + enddo + else + call lapack_diag(eigenvalues,eigenvectors, & + Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det,N_det) + CI_electronic_dressed_pt2_new_energy(:) = 0.d0 + if (s2_eig) then + i_state = 0 + allocate (s2_eigvalues(N_det)) + allocate(index_good_state_array(N_det),good_state_array(N_det)) + good_state_array = .False. + call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,& + N_det,size(eigenvectors,1)) + do j=1,N_det + ! Select at least n_states states with S^2 values closed to "expected_s2" + print*, eigenvalues(j)+nuclear_repulsion, s2_eigvalues(j) + if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then + i_state += 1 + index_good_state_array(i_state) = j + good_state_array(j) = .True. + endif + if (i_state==N_states) then + exit + endif + enddo + if (i_state /= 0) then + ! Fill the first "i_state" states that have a correct S^2 value + do j = 1, i_state + do i=1,N_det + CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) + enddo + CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(index_good_state_array(j)) + CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) + enddo + i_other_state = 0 + do j = 1, N_det + if(good_state_array(j))cycle + i_other_state +=1 + if(i_state+i_other_state.gt.n_states)then + exit + endif + do i=1,N_det + CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) + enddo + CI_electronic_dressed_pt2_new_energy(i_state+i_other_state) = eigenvalues(j) + CI_dressed_pt2_new_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) + enddo + + else + print*,'' + print*,'!!!!!!!! WARNING !!!!!!!!!' + print*,' Within the ',N_det,'determinants selected' + print*,' and the ',N_states_diag,'states requested' + print*,' We did not find any state with S^2 values close to ',expected_s2 + print*,' We will then set the first N_states eigenvectors of the H matrix' + print*,' as the CI_dressed_pt2_new_eigenvectors' + print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' + print*,'' + do j=1,min(N_states_diag,N_det) + do i=1,N_det + CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) + enddo + CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) + CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j) + enddo + endif + deallocate(index_good_state_array,good_state_array) + deallocate(s2_eigvalues) + else + call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,& + min(N_det,N_states_diag),size(eigenvectors,1)) + ! Select the "N_states_diag" states of lowest energy + do j=1,min(N_det,N_states) + do i=1,N_det + CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) + enddo + CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) + enddo + endif + deallocate(eigenvectors,eigenvalues) endif - deallocate(eigenvectors,eigenvalues) endif @@ -361,7 +412,7 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ] CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion write(st,'(I4)') j call write_double(output_determinants,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st)) - call write_double(output_determinants,CI_eigenvectors_s2(j),'S^2 of state '//trim(st)) + call write_double(output_determinants, CI_dressed_pt2_new_eigenvectors_s2(j) ,'S^2 of state '//trim(st)) enddo END_PROVIDER diff --git a/plugins/MRPT_Utils/print_1h2p.irp.f b/plugins/MRPT_Utils/print_1h2p.irp.f index 03851d8a..747e2817 100644 --- a/plugins/MRPT_Utils/print_1h2p.irp.f +++ b/plugins/MRPT_Utils/print_1h2p.irp.f @@ -8,8 +8,10 @@ end subroutine routine_2 implicit none integer :: i,j - do i =1, n_inact_orb - write(*,'(100(F16.10,X))')one_anhil_one_creat_inact_virt(i,:,1) + do i =1, n_act_orb +!do i =1, 2 + write(*,'(I3,x,100(F16.10,X))')i,one_anhil_one_creat(i,:,:,:,1) +! write(*,'(I3,x,100(F16.10,X))')i,one_anhil_one_creat(1,4,1,2,1) enddo diff --git a/plugins/MRPT_Utils/special_hij.irp.f b/plugins/MRPT_Utils/special_hij.irp.f deleted file mode 100644 index 597d8ee3..00000000 --- a/plugins/MRPT_Utils/special_hij.irp.f +++ /dev/null @@ -1,183 +0,0 @@ - - -subroutine i_H_j_no_k_operators_from_act(key_i,key_j,Nint,hij) - use bitmasks - implicit none - BEGIN_DOC - ! Returns where i and j are determinants - END_DOC - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) - double precision, intent(out) :: hij - - integer :: exc(0:2,2,2) - integer :: degree - double precision :: get_mo_bielec_integral, phase - integer :: m,n,p,q - integer :: i,j,k - integer :: occ(Nint*bit_kind_size,2) - double precision :: diag_H_mat_elem - integer :: n_occ_ab(2) - logical :: has_mipi(Nint*bit_kind_size) - double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size), miip_other(Nint*bit_kind_size) - PROVIDE mo_bielec_integrals_in_map mo_integrals_map - - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) - ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) - ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) - ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) - - hij = 0.d0 - !DIR$ FORCEINLINE - call get_excitation_degree(key_i,key_j,degree,Nint) - select case (degree) - case (2) - call get_double_excitation(key_i,key_j,exc,phase,Nint) - if (exc(0,1,1) == 1) then - ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral( & - exc(1,1,1), & - exc(1,1,2), & - exc(1,2,1), & - exc(1,2,2) ,mo_integrals_map) - else if (exc(0,1,1) == 2) then - ! Double alpha - hij = phase*(get_mo_bielec_integral( & - exc(1,1,1), & - exc(2,1,1), & - exc(1,2,1), & - exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & - exc(1,1,1), & - exc(2,1,1), & - exc(2,2,1), & - exc(1,2,1) ,mo_integrals_map) ) - else if (exc(0,1,2) == 2) then - ! Double beta - hij = phase*(get_mo_bielec_integral( & - exc(1,1,2), & - exc(2,1,2), & - exc(1,2,2), & - exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & - exc(1,1,2), & - exc(2,1,2), & - exc(2,2,2), & - exc(1,2,2) ,mo_integrals_map) ) - endif - case (1) - call get_mono_excitation(key_i,key_j,exc,phase,Nint) - !DIR$ FORCEINLINE - call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) - has_mipi = .False. - logical :: is_i_in_active - double precision :: accu_a, accu_b, accu_core - accu_a = 0.d0 - accu_b = 0.d0 - if (exc(0,1,1) == 1) then - ! Mono alpha - m = exc(1,1,1) - p = exc(1,2,1) - do k = 1, elec_alpha_num - i = occ(k,1) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif - if(.not.is_i_in_active(i))then - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - else -! print*, i,get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - miip(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - accu_a += miip(i) - endif - enddo - do k = 1, elec_beta_num - i = occ(k,2) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif - if(.not.is_i_in_active(i))then - miip_other(i) = 0.d0 - else -! print*, i,get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - miip_other(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - accu_b += miip(i) - endif - enddo -! print*, accu_a,accu_b,accu_a + accu_b - accu_a = 0.d0 - accu_b = 0.d0 - accu_core = mo_mono_elec_integral(m,p) - - do k = 1, elec_alpha_num - hij = hij + mipi(occ(k,1)) - miip(occ(k,1)) - accu_Core += mipi(occ(k,1)) - if(is_i_in_active(occ(k,1)))then - accu_a += miip(occ(k,1)) - else - accu_Core -= miip(occ(k,1)) - endif - enddo -! print*, hij,accu_core - do k = 1, elec_beta_num - hij = hij + mipi(occ(k,2)) - miip_other(occ(k,2)) - accu_Core += mipi(occ(k,2)) - if(is_i_in_active(occ(k,2)))then - accu_b += miip_other(occ(k,2)) - else - accu_Core -= miip_other(occ(k,2)) - endif - enddo -! print*, hij,accu_core,accu_core - accu_a - accu_b -! print*, accu_a,accu_b,accu_a + accu_b -! pause - - else - ! Mono beta - m = exc(1,1,2) - p = exc(1,2,2) - do k = 1, elec_beta_num - i = occ(k,2) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif - if(.not.is_i_in_active(i))then - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - else - miip(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - endif - enddo - do k = 1, elec_alpha_num - i = occ(k,1) - if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - has_mipi(i) = .True. - endif - if(.not.is_i_in_active(i))then - miip_other(i) = 0.d0 - else - miip_other(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map) - endif - enddo - do k = 1, elec_alpha_num - hij = hij + mipi(occ(k,1)) - miip_other(occ(k,1)) - enddo - do k = 1, elec_beta_num - hij = hij + mipi(occ(k,2)) - miip(occ(k,2)) - enddo - - endif - hij = phase*(hij + mo_mono_elec_integral(m,p)) - - case (0) - hij = diag_H_mat_elem(key_i,Nint) - end select -end - - - diff --git a/plugins/loc_cele/loc_exchange_int.irp.f b/plugins/loc_cele/loc_exchange_int.irp.f index d7cc5c65..fa3cf8bf 100644 --- a/plugins/loc_cele/loc_exchange_int.irp.f +++ b/plugins/loc_cele/loc_exchange_int.irp.f @@ -13,21 +13,21 @@ program loc_int iorb = list_core_inact(i) exchange_int = 0.d0 iorder = 0 - print*,'' if(list_core_inact_check(iorb) == .False.)cycle do j = i+1, n_core_inact_orb jorb = list_core_inact(j) iorder(jorb) = jorb - exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + if(list_core_inact_check(jorb) == .False.)then + exchange_int(jorb) = 0.d0 + else + exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + endif enddo n_rot += 1 call dsort(exchange_int,iorder,mo_tot_num) indices(n_rot,1) = iorb indices(n_rot,2) = iorder(1) list_core_inact_check(iorder(1)) = .False. - print*,indices(n_rot,1),indices(n_rot,2) - print*,'' - print*,'' enddo print*,'****************************' print*,'-+++++++++++++++++++++++++' @@ -45,21 +45,21 @@ program loc_int iorb = list_act(i) exchange_int = 0.d0 iorder = 0 - print*,'' if(list_core_inact_check(iorb) == .False.)cycle do j = i+1, n_act_orb jorb = list_act(j) iorder(jorb) = jorb - exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + if(list_core_inact_check(jorb) == .False.)then + exchange_int(jorb) = 0.d0 + else + exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + endif enddo n_rot += 1 call dsort(exchange_int,iorder,mo_tot_num) indices(n_rot,1) = iorb indices(n_rot,2) = iorder(1) list_core_inact_check(iorder(1)) = .False. - print*,indices(n_rot,1),indices(n_rot,2) - print*,'' - print*,'' enddo print*,'****************************' print*,'-+++++++++++++++++++++++++' @@ -82,16 +82,17 @@ program loc_int do j = i+1, n_virt_orb jorb = list_virt(j) iorder(jorb) = jorb - exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + if(list_core_inact_check(jorb) == .False.)then + exchange_int(jorb) = 0.d0 + else + exchange_int(jorb) = -mo_bielec_integral_jj_exchange(iorb,jorb) + endif enddo n_rot += 1 call dsort(exchange_int,iorder,mo_tot_num) indices(n_rot,1) = iorb indices(n_rot,2) = iorder(1) list_core_inact_check(iorder(1)) = .False. - print*,indices(n_rot,1),indices(n_rot,2) - print*,'' - print*,'' enddo print*,'****************************' print*,'-+++++++++++++++++++++++++' diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index ed299447..47d02758 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -2168,9 +2168,27 @@ subroutine H_u_0_stored(v_0,u_0,hmatrix,sze) double precision, intent(in) :: u_0(sze) v_0 = 0.d0 call matrix_vector_product(u_0,v_0,hmatrix,sze,sze) - end +subroutine H_s2_u_0_stored(v_0,u_0,hmatrix,s2matrix,sze) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> + ! + ! n : number of determinants + ! + ! uses the big_matrix_stored array + END_DOC + integer, intent(in) :: sze + double precision, intent(in) :: hmatrix(sze,sze),s2matrix(sze,sze) + double precision, intent(out) :: v_0(sze) + double precision, intent(in) :: u_0(sze) + v_0 = 0.d0 + call matrix_vector_product(u_0,v_0,hmatrix,sze,sze) +end + + subroutine u_0_H_u_0_stored(e_0,u_0,hmatrix,sze) use bitmasks implicit none