From 4b1a77c5c0fbd34825fd3617d6e69ce608731e9f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 23 May 2018 18:16:33 +0200 Subject: [PATCH] Iterative sbk, broken --- plugins/Psiref_Utils/psi_ref_utils.irp.f | 4 +- plugins/dress_zmq/dress_general.irp.f | 4 +- plugins/dress_zmq/dress_stoch_routines.irp.f | 12 +- plugins/dress_zmq/dressing.irp.f | 3 +- plugins/shiftedbk/EZFIO.cfg | 39 +++++ plugins/shiftedbk/selection_buffer.irp.f | 86 ++++++++++ plugins/shiftedbk/shifted_bk_iter.irp.f | 159 +++++++++++++++++++ plugins/shiftedbk/shifted_bk_routines.irp.f | 21 ++- src/Determinants/H_apply.irp.f | 4 +- src/Determinants/determinants.irp.f | 15 +- src/Determinants/psi_cas.irp.f | 4 +- 11 files changed, 314 insertions(+), 37 deletions(-) create mode 100644 plugins/shiftedbk/shifted_bk_iter.irp.f diff --git a/plugins/Psiref_Utils/psi_ref_utils.irp.f b/plugins/Psiref_Utils/psi_ref_utils.irp.f index c59bbd9f..a63b0ded 100644 --- a/plugins/Psiref_Utils/psi_ref_utils.irp.f +++ b/plugins/Psiref_Utils/psi_ref_utils.irp.f @@ -9,7 +9,7 @@ use bitmasks ! function. END_DOC call sort_dets_by_det_search_key(N_det_ref, psi_ref, psi_ref_coef, & - psi_ref_sorted_bit, psi_ref_coef_sorted_bit) + psi_ref_sorted_bit, psi_ref_coef_sorted_bit, N_states) END_PROVIDER @@ -152,7 +152,7 @@ END_PROVIDER ! function. END_DOC call sort_dets_by_det_search_key(N_det_ref, psi_non_ref, psi_non_ref_coef, & - psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit) + psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit, N_states) END_PROVIDER diff --git a/plugins/dress_zmq/dress_general.irp.f b/plugins/dress_zmq/dress_general.irp.f index 048605b8..b99eb1d7 100644 --- a/plugins/dress_zmq/dress_general.irp.f +++ b/plugins/dress_zmq/dress_general.irp.f @@ -38,10 +38,10 @@ subroutine run_dressing(N_st,energy) E_old = sum(psi_energy(:)) print *, 'Variational energy ' do i=1,N_st - print *, i, psi_energy(i) + print *, i, psi_energy(i)+nuclear_repulsion enddo !print *, "DELTA IJ", delta_ij(1,1,1) - if(.true.) dummy = delta_ij_tmp(1,1,1) + PROVIDE delta_ij_tmp if(.true.) call delta_ij_done() print *, 'Dressed energy ' do i=1,N_st diff --git a/plugins/dress_zmq/dress_stoch_routines.irp.f b/plugins/dress_zmq/dress_stoch_routines.irp.f index 17fe0c95..6b7bf396 100644 --- a/plugins/dress_zmq/dress_stoch_routines.irp.f +++ b/plugins/dress_zmq/dress_stoch_routines.irp.f @@ -229,10 +229,8 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, pullLoop : do while (loop) call pull_dress_results(zmq_socket_pull, ind, cur_cp, delta_loc, int_buf, double_buf, det_buf, N_buf, task_id, dress_mwen) - !print *, cur_cp, ind if(floop) then call wall_time(time) - print *, "first_pull", time-time0 time0 = time floop = .false. end if @@ -240,8 +238,6 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, call wall_time(time) end if - print *, cur_cp, ind - if(cur_cp == -1) then call dress_pulled(ind, int_buf, double_buf, det_buf, N_buf) if (zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,1,more) == -1) then @@ -295,7 +291,7 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, !print '(I2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', avg+E(istate)+E0, eqt, time-time0, '' print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', cps_N(cur_cp), avg+E0+E(istate), eqt, time-time0, '' - if ((dabs(eqt) < relative_error .and. cps_N(cur_cp) >= 30)) then + if ((dabs(eqt/(avg+E0+E(istate))) < relative_error .and. cps_N(cur_cp) >= 10)) then ! Termination print *, "TERMINATE" if (zmq_abort(zmq_to_qp_run_socket) == -1) then @@ -511,12 +507,6 @@ END_PROVIDER end do -print *, 'needed_by_cp' -do i=1,cur_cp - print *, i, needed_by_cp(i) -enddo - - under_det = 0 cp_first_tooth = 0 do i=1,N_dress_jobs diff --git a/plugins/dress_zmq/dressing.irp.f b/plugins/dress_zmq/dressing.irp.f index 12fef07a..bf2ab207 100644 --- a/plugins/dress_zmq/dressing.irp.f +++ b/plugins/dress_zmq/dressing.irp.f @@ -96,12 +96,13 @@ BEGIN_PROVIDER [ double precision, delta_ij_tmp, (N_states,N_det_delta_ij,2) ] E_CI_before(:) = psi_energy(:) + nuclear_repulsion threshold_selectors = 1.d0 threshold_generators = 1.d0 + SOFT_TOUCH threshold_selectors threshold_generators ! if(errr /= 0d0) then ! errr = errr / 2d0 ! else ! errr = 1d-4 ! end if - relative_error = 5.d-5 + relative_error = 1.d-2 call write_double(6,relative_error,"Convergence of the stochastic algorithm") diff --git a/plugins/shiftedbk/EZFIO.cfg b/plugins/shiftedbk/EZFIO.cfg index 001535b9..be4998dd 100644 --- a/plugins/shiftedbk/EZFIO.cfg +++ b/plugins/shiftedbk/EZFIO.cfg @@ -3,3 +3,42 @@ type: Perturbation doc: Type of zeroth-order Hamiltonian [ EN | Barycentric ] interface: ezfio,provider,ocaml default: EN +[energy] +type: double precision +doc: Calculated Selected FCI energy +interface: ezfio + +[energy_pt2] +type: double precision +doc: Calculated FCI energy + PT2 +interface: ezfio + +[iterative_save] +type: integer +doc: Save data at each iteration : 1(Append) | 2(Overwrite) | 3(NoSave) +interface: ezfio,ocaml +default: 2 + +[n_iter] +interface: ezfio +doc: number of iterations +type:integer + +[n_det_iter] +interface: ezfio +doc: number of determinants at iteration +type: integer +size: (full_ci_zmq.n_iter) + +[energy_iter] +interface: ezfio +doc: The energy without a pt2 correction for n_det +type: double precision +size: (determinants.n_states,full_ci_zmq.n_iter) + +[pt2_iter] +interface: ezfio +doc: The pt2 correction for n_det +type: double precision +size: (determinants.n_states,full_ci_zmq.n_iter) + diff --git a/plugins/shiftedbk/selection_buffer.irp.f b/plugins/shiftedbk/selection_buffer.irp.f index 8c3bee91..8b140666 100644 --- a/plugins/shiftedbk/selection_buffer.irp.f +++ b/plugins/shiftedbk/selection_buffer.irp.f @@ -156,3 +156,89 @@ end subroutine +subroutine unique_selection_buffer(b) + use selection_types + implicit none + BEGIN_DOC +! Removes duplicate determinants in the selection buffer + END_DOC + type(selection_buffer), intent(inout) :: b + integer, allocatable :: iorder(:) + integer(bit_kind), pointer :: detmp(:,:,:) + double precision, pointer :: val(:) + integer :: i,j,k + integer(bit_kind), allocatable :: bit_tmp(:) + logical,allocatable :: duplicate(:) + + logical :: found_duplicates + integer*8, external :: det_search_key + + if (b%N == 0 .or. b%cur == 0) return + allocate (duplicate(b%cur), val(size(b%val)), detmp(N_int, 2, size(b%val)), bit_tmp(b%cur)) + call sort_dets_by_det_search_key(b%cur, b%det, b%val, detmp, val, 1) + + deallocate(b%det, b%val) + do i=b%cur+1,b%N + val(i) = 0.d0 + detmp(1:N_int,1:2,i) = 0_bit_kind + enddo + b%det => detmp + b%val => val + + do i=1,b%cur + bit_tmp(i) = det_search_key(b%det(1,1,i),N_int) + duplicate(i) = .False. + enddo + + do i=1,b%cur-1 + if (duplicate(i)) then + cycle + endif + j = i+1 + do while (bit_tmp(j)==bit_tmp(i)) + if (duplicate(j)) then + j += 1 + if (j > b%cur) then + exit + else + cycle + endif + endif + duplicate(j) = .True. + do k=1,N_int + if ( (b%det(k,1,i) /= b%det(k,1,j) ) & + .or. (b%det(k,2,i) /= b%det(k,2,j) ) ) then + duplicate(j) = .False. + exit + endif + enddo + j += 1 + if (j > b%cur) then + exit + endif + enddo + enddo + + found_duplicates = .False. + do i=1,b%cur + if (duplicate(i)) then + found_duplicates = .True. + exit + endif + enddo + + if (found_duplicates) then + k=0 + do i=1,N_det + if (.not.duplicate(i)) then + k += 1 + b%det(:,:,k) = b%det(:,:,i) + b%val(k) = b%val(i) + endif + enddo + b%cur = k + endif + deallocate (duplicate,bit_tmp) +end + + diff --git a/plugins/shiftedbk/shifted_bk_iter.irp.f b/plugins/shiftedbk/shifted_bk_iter.irp.f new file mode 100644 index 00000000..429efa4b --- /dev/null +++ b/plugins/shiftedbk/shifted_bk_iter.irp.f @@ -0,0 +1,159 @@ +program shifted_bk + implicit none + integer :: i,j,k + double precision, allocatable :: pt2(:) + integer :: degree + integer :: n_det_before + double precision :: threshold_davidson_in + + allocate (pt2(N_states)) + + double precision :: hf_energy_ref + logical :: has + double precision :: relative_error, absolute_error + integer :: N_states_p + character*(512) :: fmt + + PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique + PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order + PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns + PROVIDE psi_bilinear_matrix_transp_order + + + pt2 = -huge(1.e0) + threshold_davidson_in = threshold_davidson + threshold_davidson = threshold_davidson_in * 100.d0 + SOFT_TOUCH threshold_davidson + + call diagonalize_CI_dressed + call save_wavefunction + + call ezfio_has_hartree_fock_energy(has) + if (has) then + call ezfio_get_hartree_fock_energy(hf_energy_ref) + else + hf_energy_ref = ref_bitmask_energy + endif + + 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 + call diagonalize_CI_dressed + call save_wavefunction + N_states_p = min(N_det,N_states) + endif + + n_det_before = 0 + + character*(8) :: pt2_string + double precision :: threshold_selectors_save, threshold_generators_save + threshold_selectors_save = threshold_selectors + threshold_generators_save = threshold_generators + double precision :: error(N_states), energy(N_states) + error = 0.d0 + + threshold_selectors = 1.d0 + threshold_generators = 1d0 + + if (.True.) then + pt2_string = '(sh-Bk) ' + do while ( (N_det < N_det_max) ) + write(*,'(A)') '--------------------------------------------------------------------------------' + + N_det_delta_ij = N_det + + do i=1,N_states + energy(i) = psi_energy(i)+nuclear_repulsion + enddo + + PROVIDE delta_ij_tmp + call delta_ij_done() + + call diagonalize_ci_dressed + do i=1,N_states + pt2(i) = ci_energy_dressed(i) - energy(i) + enddo + + N_states_p = min(N_det,N_states) + + print *, '' + print '(A,I12)', 'Summary at N_det = ', N_det + print '(A)', '-----------------------------------' + print *, '' + print *, '' + + write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' + write(*,fmt) + write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))' + write(*,fmt) ('State',k, k=1,N_states_p) + write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' + write(*,fmt) + write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))' + write(*,fmt) '# E ', energy(1:N_states_p) + if (N_states_p > 1) then + write(*,fmt) '# Excit. (au)', energy(1:N_states_p)-energy(1) + write(*,fmt) '# Excit. (eV)', (energy(1:N_states_p)-energy(1))*27.211396641308d0 + endif + write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))' + write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p) + write(*,'(A)') '#' + write(*,fmt) '# E+PT2 ', (energy(k)+pt2(k),error(k), k=1,N_states_p) + if (N_states_p > 1) then + write(*,fmt) '# Excit. (au)', ( (energy(k)+pt2(k)-energy(1)-pt2(1)), & + dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p) + write(*,fmt) '# Excit. (eV)', ( (energy(k)+pt2(k)-energy(1)-pt2(1))*27.211396641308d0, & + dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p) + endif + write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' + write(*,fmt) + print *, '' + + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + + do k=1, N_states_p + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', energy(k) + print *, 'E+PT2'//pt2_string//' = ', energy(k)+pt2(k) + enddo + + print *, '-----' + if(N_states.gt.1)then + print *, 'Variational Energy difference (au | eV)' + do i=2, N_states_p + print*,'Delta E = ', (energy(i) - energy(1)), & + (energy(i) - energy(1)) * 27.211396641308d0 + enddo + print *, '-----' + print*, 'Variational + perturbative Energy difference (au | eV)' + do i=2, N_states_p + print*,'Delta E = ', (energy(i)+ pt2(i) - (energy(1) + pt2(1))), & + (energy(i)+ pt2(i) - (energy(1) + pt2(1))) * 27.211396641308d0 + enddo + endif + call ezfio_set_shiftedbk_energy_pt2(energy(1)+pt2(1)) +! call dump_fci_iterations_value(N_det,energy,pt2) + + n_det_before = N_det + + PROVIDE psi_coef + PROVIDE psi_det + PROVIDE psi_det_sorted + + if (N_det >= N_det_max) then + threshold_davidson = threshold_davidson_in + end if + call save_wavefunction + call ezfio_set_shiftedbk_energy(energy(1)) + call ezfio_set_shiftedbk_energy_pt2(ci_energy_dressed(1)) + enddo + endif + + + + +end + diff --git a/plugins/shiftedbk/shifted_bk_routines.irp.f b/plugins/shiftedbk/shifted_bk_routines.irp.f index 3d4540bd..0663f75e 100644 --- a/plugins/shiftedbk/shifted_bk_routines.irp.f +++ b/plugins/shiftedbk/shifted_bk_routines.irp.f @@ -20,7 +20,7 @@ END_PROVIDER fock_diag_tmp_(:,:,:) = 0.d0 integer :: i - N_det_increase_factor = 1d0 + N_det_increase_factor = dble(N_states) n_det_add = max(1, int(float(N_det) * N_det_increase_factor)) @@ -120,17 +120,16 @@ subroutine delta_ij_done() old_det_gen = N_det_generators - if (dress_stoch_istate == N_states) then - ! Add buffer only when the last state is computed - call sort_selection_buffer(global_sb) - call fill_H_apply_buffer_no_selection(global_sb%cur,global_sb%det,N_int,0) - call copy_H_apply_buffer_to_wf() - if (s2_eig.or.(N_states > 1) ) then - call make_s2_eigenfunction - endif - call undress_with_alpha(old_generators, old_det_gen, psi_det(1,1,N_det_delta_ij+1), N_det-N_det_delta_ij) - call save_wavefunction + ! Add buffer only when the last state is computed + call unique_selection_buffer(global_sb) + call sort_selection_buffer(global_sb) + call fill_H_apply_buffer_no_selection(global_sb%cur,global_sb%det,N_int,0) + call copy_H_apply_buffer_to_wf() + if (s2_eig.or.(N_states > 1) ) then + call make_s2_eigenfunction endif + call undress_with_alpha(old_generators, old_det_gen, psi_det(1,1,N_det_delta_ij+1), N_det-N_det_delta_ij) + call save_wavefunction end subroutine diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f index cd1baa8f..3ba674f1 100644 --- a/src/Determinants/H_apply.irp.f +++ b/src/Determinants/H_apply.irp.f @@ -192,8 +192,8 @@ subroutine copy_H_apply_buffer_to_wf call normalize(psi_coef,N_det) SOFT_TOUCH N_det psi_det psi_coef - logical :: found_duplicates - call remove_duplicates_in_psi_det(found_duplicates) +! logical :: found_duplicates +! call remove_duplicates_in_psi_det(found_duplicates) end diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 8530fa64..d01c80ff 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -321,21 +321,24 @@ end subroutine END_DOC call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, & - psi_det_sorted_bit, psi_coef_sorted_bit) + psi_det_sorted_bit, psi_coef_sorted_bit, N_states) END_PROVIDER -subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out) +subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out, N_st) use bitmasks implicit none - integer, intent(in) :: Ndet + integer, intent(in) :: Ndet, N_st integer(bit_kind), intent(in) :: det_in (N_int,2,psi_det_size) - double precision , intent(in) :: coef_in(psi_det_size,N_states) + double precision , intent(in) :: coef_in(psi_det_size,N_st) integer(bit_kind), intent(out) :: det_out (N_int,2,psi_det_size) - double precision , intent(out) :: coef_out(psi_det_size,N_states) + double precision , intent(out) :: coef_out(psi_det_size,N_st) BEGIN_DOC ! Determinants are sorted are sorted according to their det_search_key. ! Useful to accelerate the search of a random determinant in the wave ! function. + ! + ! /!\ The first dimension of coef_out and coef_in need to be psi_det_size + ! END_DOC integer :: i,j,k integer, allocatable :: iorder(:) @@ -356,7 +359,7 @@ subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out) det_out(j,1,i) = det_in(j,1,iorder(i)) det_out(j,2,i) = det_in(j,2,iorder(i)) enddo - do k=1,N_states + do k=1,N_st coef_out(i,k) = coef_in(iorder(i),k) enddo enddo diff --git a/src/Determinants/psi_cas.irp.f b/src/Determinants/psi_cas.irp.f index 591843f7..58a71fbc 100644 --- a/src/Determinants/psi_cas.irp.f +++ b/src/Determinants/psi_cas.irp.f @@ -54,7 +54,7 @@ END_PROVIDER ! function. END_DOC call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, & - psi_cas_sorted_bit, psi_cas_coef_sorted_bit) + psi_cas_sorted_bit, psi_cas_coef_sorted_bit, N_states) END_PROVIDER @@ -107,7 +107,7 @@ END_PROVIDER ! function. END_DOC call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, & - psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit) + psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit, N_states) END_PROVIDER