diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index 1b8c511a..dfcca27c 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/src/BiInts/map_integrals.irp.f @@ -98,11 +98,12 @@ subroutine get_ao_bielec_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_ PROVIDE ao_bielec_integrals_in_map thresh = ao_integrals_threshold + non_zero_int = 0 if (ao_overlap_abs(j,l) < thresh) then out_val = 0.d0 return endif - + non_zero_int = 0 do i=1,sze !DIR$ FORCEINLINE diff --git a/src/CISD_SC2_selected/H_apply.irp.f b/src/CISD_SC2_selected/H_apply.irp.f index 79668af7..c40d6d4a 100644 --- a/src/CISD_SC2_selected/H_apply.irp.f +++ b/src/CISD_SC2_selected/H_apply.irp.f @@ -6,5 +6,9 @@ from perturbation import perturbations s = H_apply("PT2",SingleRef=True) s.set_perturbation("epstein_nesbet_sc2_projected") print s + +s = H_apply("PT2_en_sc2",SingleRef=True) +s.set_perturbation("epstein_nesbet_sc2") +print s END_SHELL diff --git a/src/CISD_SC2_selected/NEEDED_MODULES b/src/CISD_SC2_selected/NEEDED_MODULES index bb98e702..c0f534d2 100644 --- a/src/CISD_SC2_selected/NEEDED_MODULES +++ b/src/CISD_SC2_selected/NEEDED_MODULES @@ -1 +1 @@ -AOs BiInts Bitmask CISD SC2 CISD_selected Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils +AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils diff --git a/src/CISD_SC2_selected/cisd_sc2.ezfio_config b/src/CISD_SC2_selected/cisd_sc2.ezfio_config new file mode 100644 index 00000000..35008968 --- /dev/null +++ b/src/CISD_SC2_selected/cisd_sc2.ezfio_config @@ -0,0 +1,4 @@ +cisd_sc2_selected + n_det_max_cisd_sc2 integer + pt2_max double precision + do_pt2_end logical diff --git a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f index 3610a0fe..acce82f8 100644 --- a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f +++ b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f @@ -13,9 +13,25 @@ program cisd_sc2_selected pt2 = 1.d0 perturbation = "epstein_nesbet_sc2_projected" E_old(1) = HF_energy - davidson_threshold = 1.d-6 + davidson_threshold = 1.d-8 + if (N_det > n_det_max_cisd_sc2) then + call diagonalize_CI_SC2 + call save_wavefunction + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = n_det_max_cisd_sc2 + soft_touch N_det psi_det psi_coef + call diagonalize_CI + call save_wavefunction + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_SC2_energy + print *, 'E+PT2 = ', CI_SC2_energy+pt2 + print *, '-----' + endif - do while (maxval(abs(pt2(1:N_st))) > 1.d-4) + do while (N_det < n_det_max_cisd_sc2.and.maxval(abs(pt2(1:N_st))) > pt2_max) print*,'----' print*,'' call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st) @@ -37,37 +53,43 @@ program cisd_sc2_selected exit endif enddo - pt2 = 0.d0 - call H_apply_PT2(pt2, norm_pert, H_pert_diag, N_st) + N_det = min(n_det_max_cisd_sc2,N_det) + touch N_det psi_det psi_coef davidson_threshold = 1.d-10 touch davidson_threshold davidson_criterion - do i = 1, N_st - max = 0.d0 - - print*,'' - print*,'-------------' - print*,'for state ',i - print*,'' - do k = 1, N_det - if(dabs(psi_coef(k,i)).gt.max)then - max = dabs(psi_coef(k,i)) - imax = k - endif - enddo - double precision :: max - integer :: imax - print *, 'PT2(SC2) = ', pt2(i) - print *, 'E(SC2) = ', CI_SC2_energy(i) - print *, 'E_before(SC2)+PT2(SC2) = ', (CI_SC2_energy(i)+pt2(i)) - if(i==1)then - print *, 'E(SC2)+PT2(projctd)SC2 = ', (CI_SC2_energy(i)+H_pert_diag(i)) - endif - - print*,'greater coeficient of the state : ',dabs(psi_coef(imax,i)) - call get_excitation_degree(ref_bitmask,psi_det(1,1,imax),degree,N_int) - print*,'degree of excitation of such determinant : ',degree + call diagonalize_CI_SC2 + pt2 = 0.d0 + if(do_pt2_end)then + threshold_selectors = 1.d0 + call H_apply_PT2(pt2, norm_pert, H_pert_diag, N_st) + do i = 1, N_st + max = 0.d0 - enddo + print*,'' + print*,'-------------' + print*,'for state ',i + print*,'' + print*,'N_det = ',N_det + do k = 1, N_det + if(dabs(psi_coef(k,i)).gt.max)then + max = dabs(psi_coef(k,i)) + imax = k + endif + enddo + double precision :: max + integer :: imax + print *, 'PT2(SC2) = ', pt2(i) + print *, 'E(SC2) = ', CI_SC2_energy(i) + print *, 'E_before(SC2)+PT2(SC2) = ', CI_SC2_energy(i)+pt2(i) + print *, 'E_before(SC2)+PT2(SC2)_new = ', CI_SC2_energy(i)+pt2(i)*(1.d0+norm_pert) + + print*,'greater coeficient of the state : ',dabs(psi_coef(imax,i)) + call get_excitation_degree(ref_bitmask,psi_det(1,1,imax),degree,N_int) + print*,'degree of excitation of such determinant : ',degree + + enddo print*,'coucou' + endif + call save_wavefunction deallocate(pt2,norm_pert,H_pert_diag) end diff --git a/src/CISD_SC2_selected/options.irp.f b/src/CISD_SC2_selected/options.irp.f new file mode 100644 index 00000000..1fe3e158 --- /dev/null +++ b/src/CISD_SC2_selected/options.irp.f @@ -0,0 +1,51 @@ + BEGIN_PROVIDER [ integer, n_det_max_cisd_sc2 ] + implicit none + BEGIN_DOC + ! Get n_det_max_cisd_sc2 from EZFIO file + END_DOC + logical :: has_n_det_max_cisd_sc2 + PROVIDE ezfio_filename + call ezfio_has_cisd_sc2_selected_n_det_max_cisd_sc2(has_n_det_max_cisd_sc2) + if (has_n_det_max_cisd_sc2) then + call ezfio_get_cisd_sc2_selected_n_det_max_cisd_sc2(n_det_max_cisd_sc2) + else + n_det_max_cisd_sc2 = 1000 + call ezfio_set_cisd_sc2_selected_n_det_max_cisd_sc2(n_det_max_cisd_sc2) + endif + print*,'n_det_max_cisd_sc2 = ',n_det_max_cisd_sc2 + END_PROVIDER + + BEGIN_PROVIDER [ double precision , pt2_max ] + implicit none + BEGIN_DOC + ! Get pt2_max from EZFIO file + END_DOC + logical :: has_pt2_max + PROVIDE ezfio_filename + call ezfio_has_cisd_sc2_selected_pt2_max(has_pt2_max) + if (has_pt2_max) then + call ezfio_get_cisd_sc2_selected_pt2_max(pt2_max) + else + pt2_max = 1.d-3 + call ezfio_set_cisd_sc2_selected_pt2_max(pt2_max) + endif + print*,'pt2_max = ',pt2_max + END_PROVIDER + + BEGIN_PROVIDER [ logical, do_pt2_end ] + implicit none + BEGIN_DOC + ! Get do_pt2_end from EZFIO file + END_DOC + logical :: has_do_pt2_end + PROVIDE ezfio_filename + call ezfio_has_cisd_sc2_selected_do_pt2_end(has_do_pt2_end) + if (has_do_pt2_end) then + call ezfio_get_cisd_sc2_selected_do_pt2_end(do_pt2_end) + else + do_pt2_end = .True. + call ezfio_set_cisd_sc2_selected_do_pt2_end(do_pt2_end) + endif + print*,'do_pt2_end = ',do_pt2_end + END_PROVIDER + diff --git a/src/Dets/diagonalize_CI_SC2.irp.f b/src/Dets/diagonalize_CI_SC2.irp.f index c1b67e54..987586b2 100644 --- a/src/Dets/diagonalize_CI_SC2.irp.f +++ b/src/Dets/diagonalize_CI_SC2.irp.f @@ -33,6 +33,7 @@ END_PROVIDER do j=1,N_states do i=1,N_det +! CI_SC2_eigenvectors(i,j) = psi_coef(i,j) CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j) enddo CI_SC2_electronic_energy(j) = CI_electronic_energy(j) diff --git a/src/Dets/slater_rules.irp.f b/src/Dets/slater_rules.irp.f index 741a1e01..d14c30d9 100644 --- a/src/Dets/slater_rules.irp.f +++ b/src/Dets/slater_rules.irp.f @@ -571,6 +571,61 @@ subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx end +subroutine i_H_psi_SC2_verbose(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) + print*,'--------' + do ii=1,idx(0) + print*,'--' + i = idx(ii) + !DEC$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + if (i==1)then + print*,'i==1 !!' + endif + print*,coef(i,1) * hij,coef(i,1),hij + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo + print*,i_H_psi_array(1) + enddo + print*,'------' +end + + subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) use bitmasks diff --git a/src/Perturbation/pert_sc2.irp.f b/src/Perturbation/pert_sc2.irp.f index b0bc1c93..e5094a2a 100644 --- a/src/Perturbation/pert_sc2.irp.f +++ b/src/Perturbation/pert_sc2.irp.f @@ -35,28 +35,25 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag integer :: i,j,degree,l double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda + ASSERT (Nint == N_int) ASSERT (Nint > 0) + call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) accu_e_corr = 0.d0 !$IVDEP do i = 1, idx_repeat(0) accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i)) enddo - h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) - delta_E = (CI_SC2_electronic_energy(1) - h) - delta_E = 1.d0/delta_E - c_pert(1) = i_H_psi_array(1) * delta_E + c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) e_2_pert(1) = i_H_psi_array(1) * c_pert(1) do i =2,N_st H_pert_diag(i) = h -! if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then -! c_pert(i) = -1.d0 -! e_2_pert(i) = -2.d0 -! else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) @@ -67,17 +64,16 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag enddo degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) + & - popcnt(xor( ref_bitmask(1,2), det_pert(1,2))) + popcnt(xor( ref_bitmask(1,2), det_pert(1,2))) !DEC$ NOUNROLL do l=2,Nint degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) + & - popcnt(xor( ref_bitmask(l,2), det_pert(l,2))) + popcnt(xor( ref_bitmask(l,2), det_pert(l,2))) enddo if(degree==4)then ! - call i_H_j(ref_bitmask,det_pert,Nint,h0j) - H_pert_diag(1) = c_pert(1) * h0j/coef_hf_selector - e_2_pert_fonda = H_pert_diag(1) + H_pert_diag(1) = e_2_pert(1) + e_2_pert_fonda = H_pert_diag(1) do i = 1, N_st do j = 1, idx_repeat(0) e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i) @@ -115,3 +111,46 @@ double precision function repeat_all_e_corr(key_in) repeat_all_e_corr = accu 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(N_st) + 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, but with the CISD_SC2 energies and coefficients + ! + ! c_pert(i) = /( E(i) - ) + ! + ! e_2_pert(i) = ^2/( E(i) - ) + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem, h + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) + h = diag_H_mat_elem(det_pert,Nint) + do i =1,N_st + if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then + c_pert(i) = -1.d0 + e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 + else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h) + H_pert_diag(i) = h*c_pert(i)*c_pert(i) + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + else + c_pert(i) = -1.d0 + e_2_pert(i) = -dabs(i_H_psi_array(i)) + H_pert_diag(i) = h + endif + enddo + +end diff --git a/src/SC2/ASSUMPTIONS.rst b/src/SC2/ASSUMPTIONS.rst deleted file mode 100644 index e69de29b..00000000 diff --git a/src/SC2/NEEDED_MODULES b/src/SC2/NEEDED_MODULES deleted file mode 100644 index 9c1880bb..00000000 --- a/src/SC2/NEEDED_MODULES +++ /dev/null @@ -1 +0,0 @@ -AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output SingleRefMethod Utils Selectors_full diff --git a/src/SC2/README.rst b/src/SC2/README.rst deleted file mode 100644 index 32a5c02f..00000000 --- a/src/SC2/README.rst +++ /dev/null @@ -1,37 +0,0 @@ -=============== -CISD_SC2 Module -=============== - -Documentation -============= - -.. Do not edit this section. It was auto-generated from the -.. NEEDED_MODULES file. - -`cisd `_ - Undocumented - - - -Needed Modules -============== - -.. Do not edit this section. It was auto-generated from the -.. NEEDED_MODULES file. - -* `AOs `_ -* `BiInts `_ -* `Bitmask `_ -* `CISD `_ -* `Dets `_ -* `Electrons `_ -* `Ezfio_files `_ -* `Hartree_Fock `_ -* `MonoInts `_ -* `MOs `_ -* `Nuclei `_ -* `Output `_ -* `SingleRefMethod `_ -* `Utils `_ -* `Selectors_full `_ - diff --git a/src/SC2/cisd_SC2.irp.f b/src/SC2/cisd_SC2.irp.f deleted file mode 100644 index 6cae2e9c..00000000 --- a/src/SC2/cisd_SC2.irp.f +++ /dev/null @@ -1,14 +0,0 @@ -program cisd - implicit none - integer :: i - - print *, ' HF = ', HF_energy - print *, 'N_states = ', N_states - call H_apply_cisd - print *, 'N_det = ', N_det - do i = 1,N_states - print *, 'energy = ',CI_SC2_energy(i) - print *, 'E_corr = ',CI_SC2_electronic_energy(i) - ref_bitmask_energy - enddo - -end diff --git a/src/Selectors_full/e_corr_selectors.irp.f b/src/Selectors_full/e_corr_selectors.irp.f index 3ea78540..e210c482 100644 --- a/src/Selectors_full/e_corr_selectors.irp.f +++ b/src/Selectors_full/e_corr_selectors.irp.f @@ -25,7 +25,13 @@ use bitmasks enddo END_PROVIDER BEGIN_PROVIDER[double precision, coef_hf_selector] + &BEGIN_PROVIDER[double precision, inv_selectors_coef_hf] + &BEGIN_PROVIDER[double precision, inv_selectors_coef_hf_squared] &BEGIN_PROVIDER[double precision, E_corr_per_selectors, (N_det_selectors)] + &BEGIN_PROVIDER[double precision, i_H_HF_per_selectors, (N_det_selectors)] + &BEGIN_PROVIDER[double precision, Delta_E_per_selector, (N_det_selectors)] + &BEGIN_PROVIDER[double precision, E_corr_double_only ] + &BEGIN_PROVIDER[double precision, E_corr_second_order ] implicit none BEGIN_DOC ! energy of correlation per determinant respect to the Hartree Fock determinant @@ -39,25 +45,33 @@ END_PROVIDER ! coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants END_DOC integer :: i,degree - double precision :: hij,inv_selectors_coef_hf + double precision :: hij,diag_H_mat_elem + E_corr_double_only = 0.d0 + E_corr_second_order = 0.d0 do i = 1, N_det_selectors if(exc_degree_per_selectors(i)==2)then call i_H_j(ref_bitmask,psi_selectors(1,1,i),N_int,hij) + i_H_HF_per_selectors(i) = hij E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij + E_corr_double_only += E_corr_per_selectors(i) + E_corr_second_order += hij * hij /(ref_bitmask_energy - diag_H_mat_elem(psi_selectors(1,1,i),N_int)) elseif(exc_degree_per_selectors(i) == 0)then coef_hf_selector = psi_selectors_coef(i,1) E_corr_per_selectors(i) = -1000.d0 + Delta_E_per_selector(i) = 0.d0 else E_corr_per_selectors(i) = -1000.d0 endif enddo if (dabs(coef_hf_selector) > 1.d-8) then inv_selectors_coef_hf = 1.d0/coef_hf_selector + inv_selectors_coef_hf_squared = inv_selectors_coef_hf * inv_selectors_coef_hf else inv_selectors_coef_hf = 0.d0 + inv_selectors_coef_hf_squared = 0.d0 endif do i = 1,n_double_selectors E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf enddo - + E_corr_double_only = E_corr_double_only * inv_selectors_coef_hf END_PROVIDER