diff --git a/plugins/CAS_SD/H_apply.irp.f b/plugins/CAS_SD/H_apply.irp.f index 35c45fb6..f1d0c66b 100644 --- a/plugins/CAS_SD/H_apply.irp.f +++ b/plugins/CAS_SD/H_apply.irp.f @@ -3,6 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * s = H_apply("CAS_SD") +s.unset_skip() print s s = H_apply("CAS_SD_selected_no_skip") @@ -12,6 +13,7 @@ print s s = H_apply("CAS_SD_selected") s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() print s s = H_apply("CAS_SD_PT2") @@ -22,13 +24,9 @@ print s s = H_apply("CAS_S",do_double_exc=False) print s -s = H_apply("CAS_S_selected_no_skip",do_double_exc=False) -s.set_selection_pt2("epstein_nesbet_2x2") -s.unset_skip() -print s - s = H_apply("CAS_S_selected",do_double_exc=False) s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() print s s = H_apply("CAS_S_PT2",do_double_exc=False) diff --git a/plugins/CAS_SD/cas_s_selected.irp.f b/plugins/CAS_SD/cas_s_selected.irp.f index 802de171..7c77b529 100644 --- a/plugins/CAS_SD/cas_s_selected.irp.f +++ b/plugins/CAS_SD/cas_s_selected.irp.f @@ -12,6 +12,7 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction @@ -28,49 +29,84 @@ program full_ci print *, 'E+PT2 = ', CI_energy+pt2 print *, '-----' endif + double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states) + double precision :: E_CI_before(N_states) + if(read_wf)then + call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array) + h = diag_H_mat_elem(psi_det(1,1,N_det),N_int) + selection_criterion = dabs(psi_coef(N_det,1) * (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0 + soft_touch selection_criterion + endif + + integer :: n_det_before + print*,'Beginning the selection ...' + E_CI_before(1:N_states) = CI_energy(1:N_states) do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) - call H_apply_CAS_S_selected_no_skip(pt2, norm_pert, H_pert_diag, N_st) + n_det_before = N_det + call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag, N_st) PROVIDE psi_coef PROVIDE psi_det PROVIDE psi_det_sorted - if (N_det > N_det_max) then - psi_det = psi_det_sorted - psi_coef = psi_coef_sorted - N_det = N_det_max - soft_touch N_det psi_det psi_coef - endif call diagonalize_CI + + if (N_det > N_det_max) then + N_det = N_det_max + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted + endif + + call save_wavefunction - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + do k = 1, N_states + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', CI_energy(k) + print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k) + enddo print *, '-----' + if(N_states.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_states + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_states + print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1)) + enddo + endif + E_CI_before(1:N_states) = CI_energy(1:N_states) call ezfio_set_cas_sd_energy(CI_energy(1)) enddo - call diagonalize_CI - + N_det = min(N_det_max,N_det) + touch N_det psi_det psi_coef + call diagonalize_CI if(do_pt2_end)then print*,'Last iteration only to compute the PT2' threshold_selectors = 1.d0 threshold_generators = 0.999d0 - call H_apply_CAS_S_PT2(pt2, norm_pert, H_pert_diag, N_st) + call H_apply_CAS_SD_PT2(pt2, norm_pert, H_pert_diag, N_st) print *, 'Final step' print *, 'N_det = ', N_det print *, 'N_states = ', N_states print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 + print *, 'E = ', CI_energy(1:N_states) + print *, 'E+PT2 = ', CI_energy(1:N_states)+pt2(1:N_states) print *, '-----' call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1)) endif - integer :: exc_max, degree_min exc_max = 0 print *, 'CAS determinants : ', N_det_cas @@ -79,6 +115,7 @@ program full_ci call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int) exc_max = max(exc_max,degree) enddo + print *, psi_coef_cas_diagonalized(i,:) call debug_det(psi_cas(1,1,i),N_int) print *, '' enddo diff --git a/plugins/CAS_SD/cas_sd.irp.f b/plugins/CAS_SD/cas_sd.irp.f index a5fc39b2..e2e8cb1f 100644 --- a/plugins/CAS_SD/cas_sd.irp.f +++ b/plugins/CAS_SD/cas_sd.irp.f @@ -1,7 +1,6 @@ program full_ci implicit none integer :: i,k - integer :: N_det_old double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) @@ -11,9 +10,9 @@ program full_ci character*(64) :: perturbation PROVIDE N_det_cas - N_det_old = 0 pt2 = 1.d0 diag_algorithm = "Lapack" + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction @@ -30,36 +29,68 @@ program full_ci print *, 'E+PT2 = ', CI_energy+pt2 print *, '-----' endif + double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states) + double precision :: E_CI_before(N_states) + if(read_wf)then + call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array) + h = diag_H_mat_elem(psi_det(1,1,N_det),N_int) + selection_criterion = dabs(psi_coef(N_det,1) * (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0 + soft_touch selection_criterion + endif + + integer :: n_det_before + print*,'Beginning the selection ...' + E_CI_before(1:N_states) = CI_energy(1:N_states) do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) - N_det_old = N_det + n_det_before = N_det call H_apply_CAS_SD(pt2, norm_pert, H_pert_diag, N_st) PROVIDE psi_coef PROVIDE psi_det PROVIDE psi_det_sorted - if (N_det > N_det_max) then - psi_det = psi_det_sorted - psi_coef = psi_coef_sorted - N_det = N_det_max - soft_touch N_det psi_det psi_coef - endif call diagonalize_CI - call save_wavefunction - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 - print *, '-----' - call ezfio_set_cas_sd_energy(CI_energy(1)) - if (N_det == N_det_old) then - exit - endif - enddo - call diagonalize_CI + if (N_det > N_det_max) then + N_det = N_det_max + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted + endif + + + call save_wavefunction + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + do k = 1, N_states + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', CI_energy(k) + print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k) + enddo + print *, '-----' + if(N_states.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_states + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_states + print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1)) + enddo + endif + E_CI_before(1:N_states) = CI_energy(1:N_states) + call ezfio_set_cas_sd_energy(CI_energy(1)) + enddo + N_det = min(N_det_max,N_det) + touch N_det psi_det psi_coef + call diagonalize_CI if(do_pt2_end)then print*,'Last iteration only to compute the PT2' threshold_selectors = 1.d0 @@ -70,13 +101,12 @@ program full_ci print *, 'N_det = ', N_det print *, 'N_states = ', N_states print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 + print *, 'E = ', CI_energy(1:N_states) + print *, 'E+PT2 = ', CI_energy(1:N_states)+pt2(1:N_states) print *, '-----' call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1)) endif - integer :: exc_max, degree_min exc_max = 0 print *, 'CAS determinants : ', N_det_cas @@ -85,6 +115,7 @@ program full_ci call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int) exc_max = max(exc_max,degree) enddo + print *, psi_coef_cas_diagonalized(i,:) call debug_det(psi_cas(1,1,i),N_int) print *, '' enddo diff --git a/plugins/CAS_SD/cas_sd_selected.irp.f b/plugins/CAS_SD/cas_sd_selected.irp.f index be893da4..7c77b529 100644 --- a/plugins/CAS_SD/cas_sd_selected.irp.f +++ b/plugins/CAS_SD/cas_sd_selected.irp.f @@ -12,6 +12,7 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction @@ -28,32 +29,68 @@ program full_ci print *, 'E+PT2 = ', CI_energy+pt2 print *, '-----' endif + double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states) + double precision :: E_CI_before(N_states) + if(read_wf)then + call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array) + h = diag_H_mat_elem(psi_det(1,1,N_det),N_int) + selection_criterion = dabs(psi_coef(N_det,1) * (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0 + soft_touch selection_criterion + endif + + integer :: n_det_before + print*,'Beginning the selection ...' + E_CI_before(1:N_states) = CI_energy(1:N_states) do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + n_det_before = N_det call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag, N_st) PROVIDE psi_coef PROVIDE psi_det PROVIDE psi_det_sorted - if (N_det > N_det_max) then - psi_det = psi_det_sorted - psi_coef = psi_coef_sorted - N_det = N_det_max - soft_touch N_det psi_det psi_coef - endif call diagonalize_CI + + if (N_det > N_det_max) then + N_det = N_det_max + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted + endif + + call save_wavefunction - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy(1:N_states) - print *, 'E+PT2 = ', CI_energy(1:N_states)+pt2(1:N_states) + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + do k = 1, N_states + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', CI_energy(k) + print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k) + enddo print *, '-----' + if(N_states.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_states + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_states + print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1)) + enddo + endif + E_CI_before(1:N_states) = CI_energy(1:N_states) call ezfio_set_cas_sd_energy(CI_energy(1)) enddo - call diagonalize_CI - + N_det = min(N_det_max,N_det) + touch N_det psi_det psi_coef + call diagonalize_CI if(do_pt2_end)then print*,'Last iteration only to compute the PT2' threshold_selectors = 1.d0 @@ -70,7 +107,6 @@ program full_ci call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1)) endif - integer :: exc_max, degree_min exc_max = 0 print *, 'CAS determinants : ', N_det_cas @@ -79,6 +115,7 @@ program full_ci call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int) exc_max = max(exc_max,degree) enddo + print *, psi_coef_cas_diagonalized(i,:) call debug_det(psi_cas(1,1,i),N_int) print *, '' enddo diff --git a/plugins/Full_CI/full_ci.irp.f b/plugins/Full_CI/full_ci.irp.f index 70a6ec43..42e773eb 100644 --- a/plugins/Full_CI/full_ci.irp.f +++ b/plugins/Full_CI/full_ci.irp.f @@ -11,7 +11,7 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" - + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction diff --git a/src/Davidson/EZFIO.cfg b/src/Davidson/EZFIO.cfg index de2a6954..415e359e 100644 --- a/src/Davidson/EZFIO.cfg +++ b/src/Davidson/EZFIO.cfg @@ -7,6 +7,6 @@ default: 1.e-12 [n_states_diag] type: States_number doc: n_states_diag -default: 1 +default: 10 interface: ezfio,provider,ocaml diff --git a/src/Davidson/diagonalization_hs2.irp.f b/src/Davidson/diagonalization_hs2.irp.f index e250c7d3..6c82ce1c 100644 --- a/src/Davidson/diagonalization_hs2.irp.f +++ b/src/Davidson/diagonalization_hs2.irp.f @@ -71,7 +71,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s ! ! N_st : Number of eigenstates ! - ! N_st_diag : Number of states in which H is diagonalized + ! N_st_diag : Number of states in which H is diagonalized. Assumed > sze ! ! iunit : Unit for the I/O ! @@ -107,6 +107,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s include 'constants.include.F' !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, S, y, h, lambda + if (N_st_diag > sze) then + stop 'error in Davidson : N_st_diag > sze' + endif PROVIDE nuclear_repulsion @@ -174,16 +177,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s converged = .False. - do k=1,N_st_diag - - if (k > N_st) then + do k=N_st+1,N_st_diag do i=1,sze double precision :: r1, r2 call random_number(r1) call random_number(r2) u_in(i,k) = dsqrt(-2.d0*dlog(r1))*dcos(dtwo_pi*r2) enddo - endif ! Gram-Schmidt ! ------------ @@ -247,12 +247,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s S(i,k,iter+1) = 0.d0 enddo enddo + ! do k=1,N_st_diag ! do iter2=1,iter ! do l=1,N_st_diag ! do i=1,sze ! U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1) ! W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1) +! S(i,k,iter+1) = W(i,k,iter+1) + S(i,l,iter2)*y(l,iter2,k,1) ! enddo ! enddo ! enddo @@ -276,13 +278,16 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s do k=1,N_st_diag do i=1,sze R(i,k) = (lambda(k) * U(i,k,iter+1) - W(i,k,iter+1) ) & - * (1.d0 + s2(k) * U(i,k,iter+1) - S(i,k,iter+1) ) + * (1.d0 + s2(k) * U(i,k,iter+1) - S(i,k,iter+1) - S_z2_Sz) enddo if (k <= N_st) then residual_norm(k) = u_dot_u(R(1,k),sze) to_print(1,k) = lambda(k) + nuclear_repulsion to_print(2,k) = s2(k) to_print(3,k) = residual_norm(k) + if (residual_norm(k) > 1.e9) then + stop 'Davidson failed' + endif endif enddo diff --git a/src/Davidson/diagonalize_CI.irp.f b/src/Davidson/diagonalize_CI.irp.f index 8e9a0543..57bee09d 100644 --- a/src/Davidson/diagonalize_CI.irp.f +++ b/src/Davidson/diagonalize_CI.irp.f @@ -57,10 +57,10 @@ END_PROVIDER call davidson_diag_HS2(psi_det,CI_eigenvectors, & size(CI_eigenvectors,1),CI_electronic_energy, & - N_det,N_states,N_states_diag,N_int,output_determinants) + N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants) call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,& - N_states_diag,size(CI_eigenvectors,1)) + min(N_det,N_states_diag),size(CI_eigenvectors,1)) else if (diag_algorithm == "Lapack") then