From 04bf05d8ad904ec8c07717468d4ee53c4ea74a55 Mon Sep 17 00:00:00 2001 From: Manu Date: Sun, 1 Jun 2014 22:03:26 +0200 Subject: [PATCH] CISD_SC2_selected works better, new pertrubation sc2, better selection --- src/CISD_SC2/README.rst | 35 +++++++ src/CISD_SC2/SC2.irp.f | 12 +-- src/CISD_SC2/diagonalize_CI_SC2.irp.f | 18 ---- .../cisd_sc2_selection.irp.f | 37 +++++++- src/Dets/README.rst | 21 +++-- src/Dets/determinants.ezfio_config | 1 + src/Dets/determinants.irp.f | 18 ++++ src/Dets/diagonalize_CI.irp.f | 2 +- src/Perturbation/README.rst | 93 ++----------------- src/Perturbation/epstein_nesbet.irp.f | 5 +- .../{pert_sc2.irp.f.BROKEN => pert_sc2.irp.f} | 19 +++- src/Perturbation/selection.irp.f | 2 +- src/Utils/LinearAlgebra.irp.f | 2 + 13 files changed, 136 insertions(+), 129 deletions(-) rename src/Perturbation/{pert_sc2.irp.f.BROKEN => pert_sc2.irp.f} (85%) diff --git a/src/CISD_SC2/README.rst b/src/CISD_SC2/README.rst index 938b01aa..96e489ed 100644 --- a/src/CISD_SC2/README.rst +++ b/src/CISD_SC2/README.rst @@ -8,6 +8,41 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +`cisd_sc2 `_ + CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + Initial guess vectors are not necessarily orthonormal + +`repeat_excitation `_ + Undocumented + +`cisd `_ + Undocumented + +`ci_sc2_eigenvectors `_ + Eigenvectors/values of the CI matrix + +`ci_sc2_electronic_energy `_ + Eigenvectors/values of the CI matrix + +`ci_sc2_energy `_ + N_states lowest eigenvalues of the CI matrix + +`diagonalize_ci_sc2 `_ + Replace the coefficients of the CI states by the coefficients of the + eigenstates of the CI matrix + Needed Modules diff --git a/src/CISD_SC2/SC2.irp.f b/src/CISD_SC2/SC2.irp.f index 64eb3934..16937cc5 100644 --- a/src/CISD_SC2/SC2.irp.f +++ b/src/CISD_SC2/SC2.irp.f @@ -39,10 +39,10 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit) integer :: degree_exc(sze) integer :: i_ok double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:) - if(sze<500)then - allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) - allocate (H_matrix_tmp(size(H_matrix_all_dets,1),N_det)) - allocate (eigenvalues(N_det)) + if(sze.le.1000)then + allocate (eigenvectors(size(H_matrix_all_dets,1),sze)) + allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze)) + allocate (eigenvalues(sze)) do i = 1, sze do j = 1, sze H_matrix_tmp(i,j) = H_matrix_all_dets(i,j) @@ -119,14 +119,14 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit) H_jj_dressed(i) += accu enddo - if(sze>500)then + if(sze>1000)then call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint,output_CISD_SC2) else do i = 1,sze H_matrix_tmp(i,i) = H_jj_dressed(i) enddo call lapack_diag(eigenvalues,eigenvectors, & - H_matrix_tmp,size(H_matrix_all_dets,1),N_det) + H_matrix_tmp,size(H_matrix_all_dets,1),sze) do j=1,min(N_states,sze) do i=1,sze u_in(i,j) = eigenvectors(i,j) diff --git a/src/CISD_SC2/diagonalize_CI_SC2.irp.f b/src/CISD_SC2/diagonalize_CI_SC2.irp.f index cd01296a..c4bfa376 100644 --- a/src/CISD_SC2/diagonalize_CI_SC2.irp.f +++ b/src/CISD_SC2/diagonalize_CI_SC2.irp.f @@ -1,20 +1,3 @@ -BEGIN_PROVIDER [ character*(64), diag_algorithm ] - implicit none - BEGIN_DOC - ! Diagonalization algorithm (Davidson or Lapack) - END_DOC - if (N_det > 500) then - diag_algorithm = "Davidson" - else - diag_algorithm = "Lapack" - endif - - if (N_det < N_states) then - diag_algorithm = "Lapack" - endif - -END_PROVIDER - BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states) ] implicit none BEGIN_DOC @@ -46,7 +29,6 @@ END_PROVIDER enddo CI_SC2_electronic_energy(j) = CI_electronic_energy(j) enddo - print*,'E(CISD) = ',CI_electronic_energy(1) call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, & diff --git a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f index 14f92fad..2eec29eb 100644 --- a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f +++ b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f @@ -1,10 +1,11 @@ program cisd_sc2_selected implicit none integer :: i,k + use bitmasks double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_old(:) - integer :: N_st, iter + integer :: N_st, iter,degree character*(64) :: perturbation N_st = N_states allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st),E_old(N_st)) @@ -12,25 +13,53 @@ program cisd_sc2_selected pt2 = 1.d0 perturbation = "epstein_nesbet_sc2_projected" E_old(1) = HF_energy - do while (maxval(abs(pt2(1:N_st))) > 1.d-9) + do while (maxval(abs(pt2(1:N_st))) > 1.d-10) print*,'----' print*,'' call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st) call diagonalize_CI_SC2 print *, 'N_det = ', N_det do i = 1, N_st + print*,'state ',i print *, 'PT2(SC2) = ', pt2(i) print *, 'E(SC2) = ', CI_SC2_energy(i) - print *, 'E_before(SC2)+PT2(SC2) = ', (E_old(i)+pt2(i)) + print *, 'E_before(SC2)+PT2(SC2) = ', (E_old(i)+pt2(i)) if(i==1)then print *, 'E(SC2)+PT2(projctd)SC2 = ', (E_old(i)+H_pert_diag(i)) endif + E_old(i) = CI_SC2_energy(i) enddo ! print *, 'E corr = ', (E_old(1)) - HF_energy - E_old = CI_SC2_energy if (abort_all) then exit endif enddo + 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) = ', (E_old(i)+pt2(i)) + if(i==1)then + print *, 'E(SC2)+PT2(projctd)SC2 = ', (E_old(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 + + enddo deallocate(pt2,norm_pert,H_pert_diag) end diff --git a/src/Dets/README.rst b/src/Dets/README.rst index 208042cf..349a545b 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -133,33 +133,36 @@ Documentation `n_det `_ Number of determinants in the wave function -`n_det_reference `_ +`n_det_max_jacobi `_ + Maximum number of determinants diagonalized my jacobi + +`n_det_reference `_ Number of determinants in the reference wave function `n_states `_ Number of states to consider -`psi_average_norm_contrib `_ +`psi_average_norm_contrib `_ Contribution of determinants to the state-averaged density -`psi_average_norm_contrib_sorted `_ +`psi_average_norm_contrib_sorted `_ Wave function sorted by determinants (state-averaged) -`psi_coef `_ +`psi_coef `_ The wave function. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_coef_sorted `_ +`psi_coef_sorted `_ Wave function sorted by determinants (state-averaged) -`psi_det `_ +`psi_det `_ The wave function. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_det_size `_ +`psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_det_sorted `_ +`psi_det_sorted `_ Wave function sorted by determinants (state-averaged) `double_exc_bitmask `_ @@ -298,7 +301,7 @@ Documentation Returns where i and j are determinants `i_h_psi `_ - for the various Nstate + for the various Nstates `i_h_psi_sc2 `_ for the various Nstate diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config index 0a6ead21..f4b8b908 100644 --- a/src/Dets/determinants.ezfio_config +++ b/src/Dets/determinants.ezfio_config @@ -6,4 +6,5 @@ determinants psi_coef double precision (determinants_n_det,determinants_n_states) psi_det integer*8 (determinants_N_int*determinants_bit_kind/8,2,determinants_n_det) H_apply_threshold double precision + n_det_max_jacobi integer diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index b33bf664..45f0c3f0 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -35,6 +35,24 @@ BEGIN_PROVIDER [ integer, N_det ] END_PROVIDER +BEGIN_PROVIDER [ integer, N_det_max_jacobi ] + implicit none + BEGIN_DOC + ! Maximum number of determinants diagonalized my jacobi + END_DOC + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_determinants_n_det_max_jacobi(exists) + if (exists) then + call ezfio_get_determinants_n_det_max_jacobi(N_det_max_jacobi) + else + N_det_max_jacobi = 1500 + call ezfio_set_determinants_n_det_max_jacobi(N_det_max_jacobi) + endif + ASSERT (N_det_max_jacobi > 0) +END_PROVIDER + + BEGIN_PROVIDER [ integer, psi_det_size ] implicit none BEGIN_DOC diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f index 0dc39448..8ac7aba7 100644 --- a/src/Dets/diagonalize_CI.irp.f +++ b/src/Dets/diagonalize_CI.irp.f @@ -3,7 +3,7 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ] BEGIN_DOC ! Diagonalization algorithm (Davidson or Lapack) END_DOC - if (N_det > 500) then + if (N_det > N_det_max_jacobi) then diag_algorithm = "Davidson" else diag_algorithm = "Lapack" diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst index 0e910cdd..806ad292 100644 --- a/src/Perturbation/README.rst +++ b/src/Perturbation/README.rst @@ -95,110 +95,29 @@ Documentation `pt2_epstein_nesbet `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br - for the various n_st states. + for the various N_st states. .br c_pert(i) = /( E(i) - ) .br 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. + for the various N_st states. .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_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 CI_electronic_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 CI_electronic_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 CI_electronic_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 - -`fill_h_apply_buffer_selection `_ +`fill_h_apply_buffer_selection `_ Fill the H_apply buffer with determiants for the selection -`n_det_selectors `_ - Undocumented - -`psi_selectors `_ - On what we apply for perturbation. If selection, it may be 0.9 of the norm. - -`psi_selectors_coef `_ - On what we apply for perturbation. If selection, it may be 0.9 of the norm. - -`psi_selectors_size `_ - Undocumented - `remove_small_contributions `_ - Remove determinants with small contributions + Remove determinants with small contributions. N_states is assumed to be + provided. `selection_criterion `_ Threshold to select determinants. Set by selection routines. diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f index 36d9ab1f..eba7530e 100644 --- a/src/Perturbation/epstein_nesbet.irp.f +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -24,7 +24,10 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_s 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 (dabs(CI_electronic_energy(i) - h) > 1.d-6) then + if(CI_electronic_energy(i)>h.and.CI_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_electronic_energy(i) - h) > 1.d-6) then c_pert(i) = i_H_psi_array(i) / (CI_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) diff --git a/src/Perturbation/pert_sc2.irp.f.BROKEN b/src/Perturbation/pert_sc2.irp.f similarity index 85% rename from src/Perturbation/pert_sc2.irp.f.BROKEN rename to src/Perturbation/pert_sc2.irp.f index d9c99cd2..15b81377 100644 --- a/src/Perturbation/pert_sc2.irp.f.BROKEN +++ b/src/Perturbation/pert_sc2.irp.f @@ -34,7 +34,7 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag integer :: i,j,degree double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E - double precision :: repeat_all_e_corr,accu_e_corr_tmp + 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) @@ -51,9 +51,14 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag c_pert(1) = i_H_psi_array(1) * delta_E e_2_pert(1) = i_H_psi_array(1) * c_pert(1) H_pert_diag(1) = c_pert(1) * h0j/coef_hf_selector + e_2_pert_fonda = H_pert_diag(1) + do i =2,N_st H_pert_diag(i) = h - if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then + 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 c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h) e_2_pert(i) = c_pert(i) * i_H_psi_array(i) else @@ -61,6 +66,16 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag e_2_pert(i) = -dabs(i_H_psi_array(i)) endif enddo + + call get_excitation_degree(ref_bitmask,det_pert,degree,Nint) + if(degree==2)then + ! + 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) + enddo + enddo + endif end diff --git a/src/Perturbation/selection.irp.f b/src/Perturbation/selection.irp.f index 203dccec..254ee8a6 100644 --- a/src/Perturbation/selection.irp.f +++ b/src/Perturbation/selection.irp.f @@ -72,7 +72,7 @@ end BEGIN_DOC ! Threshold to select determinants. Set by selection routines. END_DOC - selection_criterion = 1.d0 + selection_criterion = 10.d0 selection_criterion_factor = 0.01d0 selection_criterion_min = selection_criterion diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index a5e12b31..b8a02470 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -201,6 +201,8 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) integer :: lwork, info, i,j,l,k, liwork allocate(A(nmax,n),eigenvalues(n)) +! print*,'Diagonalization by jacobi' +! print*,'n = ',n A=H lwork = 2*n*n + 6*n+ 1