diff --git a/plugins/shiftedbk/shifted_bk_routines.irp.f b/plugins/shiftedbk/shifted_bk_routines.irp.f index 98faf2ec..3bffcc0a 100644 --- a/plugins/shiftedbk/shifted_bk_routines.irp.f +++ b/plugins/shiftedbk/shifted_bk_routines.irp.f @@ -24,9 +24,13 @@ use selection_types subroutine delta_ij_done() + use bitmasks implicit none integer :: i, n_det_add + if(N_det /= N_det_delta_ij) stop "N_det /= N_det_delta_ij" + + call sort_selection_buffer(sb(1)) do i=2,Nproc @@ -36,10 +40,16 @@ subroutine delta_ij_done() call sort_selection_buffer(sb(1)) + call undress_with_alpha(sb(1)%det, sb(1)%cur) + call fill_H_apply_buffer_no_selection(sb(1)%cur,sb(1)%det,N_int,0) call copy_H_apply_buffer_to_wf() + if(N_det == N_det_delta_ij) stop "N_det == N_det_delta_ij" if (s2_eig.or.(N_states > 1) ) then - call make_s2_eigenfunction + print *, "***" + print *, "*** WARNING - SHIFTED_BK currently does not support s2_eig ***" + print *, "***" + !call make_s2_eigenfunction endif !call save_wavefunction n_det_add = max(1, int(float(N_det) * N_det_increase_factor)) @@ -47,13 +57,56 @@ subroutine delta_ij_done() call delete_selection_buffer(sb(i)) call create_selection_buffer(n_det_add, n_det_add*2, sb(i)) end do - !delta_ij = 0d0 end subroutine -subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minilist, det_minilist, n_minilist, alpha, iproc) - use bitmasks - implicit none +subroutine undress_with_alpha(alpha, n_alpha) + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: alpha(N_int,2,n_alpha) + integer, intent(in) :: n_alpha + integer, allocatable :: minilist(:) + integer(bit_kind), allocatable :: det_minilist(:,:,:) + double precision, allocatable :: delta_ij_loc(:,:,:,:) + integer :: i, j, k, ex, n_minilist, iproc + double precision :: haa, contrib + integer, external :: omp_get_thread_num + allocate(minilist(N_det), det_minilist(N_int, 2, N_det), delta_ij_loc(N_states, N_det, 2, Nproc)) + delta_ij_loc = 0d0 + print *, "UNDRESSING..." + + !$OMP PARALLEL DO DEFAULT(SHARED) SCHEDULE(STATIC,1) PRIVATE(i, j, iproc, n_minilist, ex) & + !$OMP PRIVATE(det_minilist, minilist, haa, contrib) + do i=1, n_alpha + iproc = omp_get_thread_num()+1 + if(mod(i,10000) == 0) print *, "UNDRESSING", i, "/", n_alpha, iproc + n_minilist = 0 + do j=1, N_det + call get_excitation_degree(alpha(1,1,i), psi_det(1,1,j), ex, N_int) + if(ex <= 2) then + n_minilist += 1 + det_minilist(:,:,n_minilist) = psi_det(:,:,j) + minilist(n_minilist) = j + end if + end do + if(n_minilist > 0) then + call i_h_j(alpha(1,1,i), alpha(1,1,i), N_int, haa) + call dress_with_alpha_(N_states, N_det, N_int, delta_ij_loc(1,1,1,iproc), & + minilist, det_minilist, n_minilist, alpha(1,1,i), haa, contrib, iproc) + end if + end do + !$OMP END PARALLEL DO + + do i=Nproc,1,-1 + delta_ij_tmp(:,:,:) -= delta_ij_loc(:,:,:,i) + end do +end subroutine + + +subroutine dress_with_alpha_(Nstates,Ndet,Nint,delta_ij_loc,minilist, det_minilist, n_minilist, alpha, haa, contrib, iproc) + use bitmasks + implicit none BEGIN_DOC !delta_ij_loc(:,:,1) : dressing column for H !delta_ij_loc(:,:,2) : dressing column for S2 @@ -62,25 +115,20 @@ subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minili !n_minilist : size of minilist !alpha : alpha determinant END_DOC - integer, intent(in) :: Nint, Ndet, Nstates, n_minilist, iproc, i_gen + integer, intent(in) :: Nint, Ndet, Nstates, n_minilist, iproc integer(bit_kind), intent(in) :: alpha(Nint,2), det_minilist(Nint, 2, n_minilist) integer,intent(in) :: minilist(n_minilist) double precision, intent(inout) :: delta_ij_loc(Nstates,N_det,2) - double precision :: haa, hij, sij - double precision, external :: diag_H_mat_elem_fock + double precision, intent(out) :: contrib + double precision, intent(in) :: haa + double precision :: hij, sij integer :: i,j,k,l,m, l_sd double precision :: hdress, sdress - double precision :: de, a_h_psi(Nstates), c_alpha, contrib - - + double precision :: de, a_h_psi(Nstates), c_alpha + + contrib = 0d0 a_h_psi = 0d0 - if(current_generator_(iproc) /= i_gen) then - current_generator_(iproc) = i_gen - call build_fock_tmp(fock_diag_tmp_(1,1,iproc),psi_det_generators(1,1,i_gen),N_int) - end if - - haa = diag_H_mat_elem_fock(psi_det_generators(1,1,i_gen),alpha,fock_diag_tmp_(1,1,iproc),N_int) do l_sd=1,n_minilist call i_h_j_s2(alpha,det_minilist(1,1,l_sd),N_int,hij, sij) @@ -107,7 +155,38 @@ subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minili delta_ij_loc(i, minilist(l_sd), 2) += sdress end do end do +end subroutine + +subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minilist, det_minilist, n_minilist, alpha, iproc) + use bitmasks + implicit none + BEGIN_DOC + !delta_ij_loc(:,:,1) : dressing column for H + !delta_ij_loc(:,:,2) : dressing column for S2 + !i_gen : generator index in psi_det_generators + !minilist : indices of determinants connected to alpha ( in psi_det_sorted ) + !n_minilist : size of minilist + !alpha : alpha determinant + END_DOC + integer, intent(in) :: Nint, Ndet, Nstates, n_minilist, iproc, i_gen + integer(bit_kind), intent(in) :: alpha(Nint,2), det_minilist(Nint, 2, n_minilist) + integer,intent(in) :: minilist(n_minilist) + double precision, intent(inout) :: delta_ij_loc(Nstates,N_det,2) + double precision, external :: diag_H_mat_elem_fock + double precision :: haa, contrib + + + + if(current_generator_(iproc) /= i_gen) then + current_generator_(iproc) = i_gen + call build_fock_tmp(fock_diag_tmp_(1,1,iproc),psi_det_generators(1,1,i_gen),N_int) + end if + + haa = diag_H_mat_elem_fock(psi_det_generators(1,1,i_gen),alpha,fock_diag_tmp_(1,1,iproc),N_int) + + call dress_with_alpha_(Nstates, Ndet, Nint, delta_ij_loc, minilist, det_minilist, n_minilist, alpha, haa, contrib, iproc) + call add_to_selection_buffer(sb(iproc), alpha, contrib) end subroutine