From 8718f5fd35ded4a420a34fce9d18108d0ba962d3 Mon Sep 17 00:00:00 2001 From: Manu Date: Mon, 6 Oct 2014 15:49:16 +0200 Subject: [PATCH] minor changes --- src/CISD_SC2_selected/H_apply.irp.f | 2 +- .../cisd_sc2_selection.irp.f | 14 +++- src/Dets/determinants.ezfio_config | 1 + src/Dets/diagonalize_CI.irp.f | 21 ++++-- src/Dets/diagonalize_CI_SC2.irp.f | 6 +- src/Dets/options.irp.f | 1 + src/Dets/s2.irp.f | 15 ++++ src/Perturbation/pert_sc2.irp.f | 74 ++++++++++++++++++- src/Utils/LinearAlgebra.irp.f | 66 +++++++++++++++++ src/Utils/one_e_integration.irp.f | 8 +- src/Utils/util.irp.f | 8 +- 11 files changed, 194 insertions(+), 22 deletions(-) diff --git a/src/CISD_SC2_selected/H_apply.irp.f b/src/CISD_SC2_selected/H_apply.irp.f index c40d6d4a..85ecae95 100644 --- a/src/CISD_SC2_selected/H_apply.irp.f +++ b/src/CISD_SC2_selected/H_apply.irp.f @@ -4,7 +4,7 @@ from generate_h_apply import * from perturbation import perturbations s = H_apply("PT2",SingleRef=True) -s.set_perturbation("epstein_nesbet_sc2_projected") +s.set_perturbation("epstein_nesbet_sc2_no_projected") print s s = H_apply("PT2_en_sc2",SingleRef=True) diff --git a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f index acce82f8..b8770155 100644 --- a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f +++ b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f @@ -12,8 +12,9 @@ program cisd_sc2_selected pt2 = 1.d0 perturbation = "epstein_nesbet_sc2_projected" + E_old(1) = HF_energy - davidson_threshold = 1.d-8 + davidson_threshold = 1.d-10 if (N_det > n_det_max_cisd_sc2) then call diagonalize_CI_SC2 call save_wavefunction @@ -31,6 +32,8 @@ program cisd_sc2_selected print *, '-----' endif + integer :: i_count + i_count = 0 do while (N_det < n_det_max_cisd_sc2.and.maxval(abs(pt2(1:N_st))) > pt2_max) print*,'----' print*,'' @@ -49,6 +52,13 @@ program cisd_sc2_selected E_old(i) = CI_SC2_energy(i) enddo ! print *, 'E corr = ', (E_old(1)) - HF_energy + if(dabs(E_old(i) - CI_SC2_energy(i) ).le.1.d-12)then + i_count += 1 + selection_criterion_factor = selection_criterion_factor * 0.5d0 + if(i_count > 5)then + exit + endif + endif if (abort_all) then exit endif @@ -81,7 +91,7 @@ program cisd_sc2_selected 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 *, 'E_before(SC2)+PT2(SC2)_new = ', CI_SC2_energy(i)+pt2(i)* (1.d0 + norm_pert) - H_pert_diag(i) 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) diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config index 47fe30e8..ff97ca36 100644 --- a/src/Dets/determinants.ezfio_config +++ b/src/Dets/determinants.ezfio_config @@ -13,4 +13,5 @@ determinants det_occ integer (electrons_elec_alpha_num,determinants_det_num,2) det_coef double precision (determinants_det_num) read_wf logical + expected_s2 double precision diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f index f9cc3ef0..1329b581 100644 --- a/src/Dets/diagonalize_CI.irp.f +++ b/src/Dets/diagonalize_CI.irp.f @@ -59,12 +59,21 @@ END_PROVIDER call lapack_diag(eigenvalues,eigenvectors, & H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) CI_electronic_energy(:) = 0.d0 - do j=1,min(N_states,N_det) - do i=1,N_det - CI_eigenvectors(i,j) = eigenvectors(i,j) - enddo - CI_electronic_energy(j) = eigenvalues(j) - enddo + integer :: i_state + double precision :: s2 + j=0 + i_state = 0 +! do while(i_state.lt.N_states) + j+=1 +! call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) +! if(dabs(s2-expected_s2).le.0.1d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy(i_state) = eigenvalues(j) +! endif +! enddo deallocate(eigenvectors,eigenvalues) endif diff --git a/src/Dets/diagonalize_CI_SC2.irp.f b/src/Dets/diagonalize_CI_SC2.irp.f index 987586b2..412dc2ae 100644 --- a/src/Dets/diagonalize_CI_SC2.irp.f +++ b/src/Dets/diagonalize_CI_SC2.irp.f @@ -20,7 +20,7 @@ END_PROVIDER BEGIN_DOC ! convergence of the correlation energy of SC2 iterations END_DOC - threshold_convergence_SC2 = 1.d-8 + threshold_convergence_SC2 = 1.d-10 END_PROVIDER BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states) ] @@ -33,8 +33,8 @@ 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) + 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) enddo diff --git a/src/Dets/options.irp.f b/src/Dets/options.irp.f index 2a8ffe47..df2ba1b5 100644 --- a/src/Dets/options.irp.f +++ b/src/Dets/options.irp.f @@ -20,6 +20,7 @@ T.set_doc ( "If true, read the wave function from the EZFIO file" ) T.set_ezfio_name( "read_wf" ) T.set_output ( "output_dets" ) print T + END_SHELL diff --git a/src/Dets/s2.irp.f b/src/Dets/s2.irp.f index b218e7cb..6993675c 100644 --- a/src/Dets/s2.irp.f +++ b/src/Dets/s2.irp.f @@ -42,6 +42,21 @@ BEGIN_PROVIDER [ double precision, S_z ] END_PROVIDER +BEGIN_PROVIDER [ double precision, expected_s2] + implicit none + PROVIDE ezfio_filename + logical :: has_expected_s2 + + call ezfio_has_determinants_expected_s2(has_expected_s2) + if (has_expected_s2) then + call ezfio_get_determinants_expected_s2(expected_s2) + else + expected_s2 = elec_alpha_num - elec_beta_num + 0.5d0 * ((elec_alpha_num - elec_beta_num)**2*0.5d0 - (elec_alpha_num-elec_beta_num)) + call ezfio_set_determinants_expected_s2(expected_s2) + endif + +END_PROVIDER + subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) implicit none diff --git a/src/Perturbation/pert_sc2.irp.f b/src/Perturbation/pert_sc2.irp.f index e5094a2a..f0fc8a0f 100644 --- a/src/Perturbation/pert_sc2.irp.f +++ b/src/Perturbation/pert_sc2.irp.f @@ -72,8 +72,8 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag enddo if(degree==4)then ! - H_pert_diag(1) = e_2_pert(1) - e_2_pert_fonda = H_pert_diag(1) + e_2_pert_fonda = e_2_pert(1) + H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(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) @@ -83,6 +83,76 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag end + +subroutine pt2_epstein_nesbet_SC2_no_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(N_st) + double precision :: i_H_psi_array(N_st) + integer :: idx_repeat(0: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 CI_electronic_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! In addition, for the perturbative energetic contribution you have the standard second order + ! + ! e_2_pert = ^2/(Delta_E) + ! + ! and also the purely projected contribution + ! + ! H_pert_diag = c_pert + END_DOC + + 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 + delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + + + 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 (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)) + else + c_pert(i) = i_H_psi_array(i) + e_2_pert(i) = -dabs(i_H_psi_array(i)) + endif + enddo +end + + + + + double precision function repeat_all_e_corr(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 9f2336db..b29654b1 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -307,6 +307,72 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) deallocate(A,eigenvalues) end +subroutine lapack_diag_s2(eigvalues,eigvectors,H,nmax,n) + implicit none + BEGIN_DOC + ! Diagonalize matrix H + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! + END_DOC + integer, intent(in) :: n,nmax + double precision, intent(out) :: eigvectors(nmax,n) + double precision, intent(out) :: eigvalues(n) + double precision, intent(in) :: H(nmax,n) + double precision,allocatable :: eigenvalues(:) + double precision,allocatable :: work(:) + double precision,allocatable :: A(:,:) + 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 + allocate (work(lwork)) + + lwork = -1 + call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + info ) + if (info < 0) then + print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = int( work( 1 ) ) + deallocate (work) + + allocate (work(lwork)) + call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + info ) + deallocate(work) + + if (info < 0) then + print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value' + stop 2 + else if( info > 0 ) then + write(*,*)'DSYEV Failed' + stop 1 + end if + + eigvectors = 0.d0 + eigvalues = 0.d0 + do j = 1, n + eigvalues(j) = eigenvalues(j) + do i = 1, n + eigvectors(i,j) = A(i,j) + enddo + enddo + deallocate(A,eigenvalues) +end + + + + subroutine lapack_partial_diag(eigvalues,eigvectors,H,nmax,n,n_st) implicit none BEGIN_DOC diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f index d7e52d8e..273b4f48 100644 --- a/src/Utils/one_e_integration.irp.f +++ b/src/Utils/one_e_integration.irp.f @@ -15,10 +15,10 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_ call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,& beta,power_A,power_B,A_center,B_center,dim) - if(fact_p.lt.0.000001d0)then - overlap_gaussian_x = 0.d0 - return - endif +! if(fact_p.lt.0.000001d0)then +! overlap_gaussian_x = 0.d0 +! return +! endif overlap_gaussian_x = 0.d0 integer :: i diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index 7d81d49b..ab7ce9db 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -49,17 +49,17 @@ double precision function binom_func(i,j) end - BEGIN_PROVIDER [ double precision, binom, (0:20,0:20) ] -&BEGIN_PROVIDER [ double precision, binom_transp, (0:20,0:20) ] + BEGIN_PROVIDER [ double precision, binom, (0:40,0:40) ] +&BEGIN_PROVIDER [ double precision, binom_transp, (0:40,0:40) ] implicit none BEGIN_DOC ! Binomial coefficients END_DOC integer :: k,l double precision :: fact, f - do k=0,20 + do k=0,40 f = fact(k) - do l=0,20 + do l=0,40 binom(k,l) = f/(fact(l)*fact(k-l)) binom_transp(l,k) = binom(k,l) enddo