From a865a842d20e403e212c5b955501b661b1dd0796 Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Mon, 26 Mar 2018 13:08:26 +0200 Subject: [PATCH] undressing with s2_eig --- plugins/dress_zmq/NEEDED_CHILDREN_MODULES | 2 +- plugins/dress_zmq/alpha_factory.irp.f | 3 +- plugins/dress_zmq/dressing.irp.f | 2 +- plugins/shiftedbk/NEEDED_CHILDREN_MODULES | 2 +- plugins/shiftedbk/shifted_bk_routines.irp.f | 88 +++++++++++++-------- 5 files changed, 57 insertions(+), 40 deletions(-) diff --git a/plugins/dress_zmq/NEEDED_CHILDREN_MODULES b/plugins/dress_zmq/NEEDED_CHILDREN_MODULES index 55f8f738..96b2cfdc 100644 --- a/plugins/dress_zmq/NEEDED_CHILDREN_MODULES +++ b/plugins/dress_zmq/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Selectors_full Generators_full ZMQ +ZMQ diff --git a/plugins/dress_zmq/alpha_factory.irp.f b/plugins/dress_zmq/alpha_factory.irp.f index 190a94ad..ccbf177a 100644 --- a/plugins/dress_zmq/alpha_factory.irp.f +++ b/plugins/dress_zmq/alpha_factory.irp.f @@ -84,8 +84,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index !hole (k,1) = iand(psi_det_generators(k,1,i_generator), full_ijkl_bitmask(k)) !hole (k,2) = iand(psi_det_generators(k,2,i_generator), full_ijkl_bitmask(k)) !particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), full_ijkl_bitmask(k)) - !particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), full_ijkl_bitmask(k)) - + !particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), full_ijkl_bitmask(k)) enddo integer :: N_holes(2), N_particles(2) diff --git a/plugins/dress_zmq/dressing.irp.f b/plugins/dress_zmq/dressing.irp.f index ce89415d..9f4ede26 100644 --- a/plugins/dress_zmq/dressing.irp.f +++ b/plugins/dress_zmq/dressing.irp.f @@ -100,7 +100,7 @@ BEGIN_PROVIDER [ double precision, delta_ij_tmp, (N_states,N_det_delta_ij,2) ] ! else ! errr = 1d-4 ! end if - relative_error = 1.d-4 + relative_error = 1.d-5 call write_double(6,relative_error,"Convergence of the stochastic algorithm") call ZMQ_dress(E_CI_before, dress, del, del_s2, abs(relative_error)) diff --git a/plugins/shiftedbk/NEEDED_CHILDREN_MODULES b/plugins/shiftedbk/NEEDED_CHILDREN_MODULES index c3290687..4f09bfc8 100644 --- a/plugins/shiftedbk/NEEDED_CHILDREN_MODULES +++ b/plugins/shiftedbk/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -dress_zmq DavidsonDressed +dress_zmq DavidsonDressed Selectors_full Generators_full diff --git a/plugins/shiftedbk/shifted_bk_routines.irp.f b/plugins/shiftedbk/shifted_bk_routines.irp.f index 3bffcc0a..7e32568c 100644 --- a/plugins/shiftedbk/shifted_bk_routines.irp.f +++ b/plugins/shiftedbk/shifted_bk_routines.irp.f @@ -9,12 +9,12 @@ use selection_types implicit none integer :: i integer :: n_det_add - + N_det_increase_factor = 1d0 current_generator_(:) = 0 + n_det_add = max(1, int(float(N_det) * N_det_increase_factor)) do i=1,Nproc - n_det_add = max(1, int(float(N_det) * N_det_increase_factor)) call create_selection_buffer(n_det_add, n_det_add*2, sb(i)) end do a_h_i = 0d0 @@ -26,11 +26,13 @@ use selection_types subroutine delta_ij_done() use bitmasks implicit none - integer :: i, n_det_add + integer :: i, n_det_add, old_det_gen + integer(bit_kind), allocatable :: old_generators(:,:,:) - if(N_det /= N_det_delta_ij) stop "N_det /= N_det_delta_ij" - - + allocate(old_generators(N_int, 2, N_det_generators)) + old_generators(:,:,:) = psi_det_generators(:,:,:N_det_generators) + old_det_gen = N_det_generators + call sort_selection_buffer(sb(1)) do i=2,Nproc @@ -39,28 +41,19 @@ subroutine delta_ij_done() end do 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 - print *, "***" - print *, "*** WARNING - SHIFTED_BK currently does not support s2_eig ***" - print *, "***" - !call make_s2_eigenfunction + call make_s2_eigenfunction endif - !call save_wavefunction - n_det_add = max(1, int(float(N_det) * N_det_increase_factor)) - do i=1,Nproc - call delete_selection_buffer(sb(i)) - call create_selection_buffer(n_det_add, n_det_add*2, sb(i)) - end do + 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 -subroutine undress_with_alpha(alpha, n_alpha) +subroutine undress_with_alpha(old_generators, old_det_gen, alpha, n_alpha) use bitmasks implicit none @@ -69,20 +62,46 @@ subroutine undress_with_alpha(alpha, 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 :: exc(0:2,2,2), h1, h2, p1, p2, s1, s2 + integer :: i, j, k, ex, n_minilist, iproc, degree + double precision :: haa, contrib, phase + logical :: ok 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)) + + integer,intent(in) :: old_det_gen + integer(bit_kind), intent(in) :: old_generators(N_int, 2, old_det_gen) + + allocate(minilist(N_det_delta_ij), det_minilist(N_int, 2, N_det_delta_ij), delta_ij_loc(N_states, N_det_delta_ij, 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 + !$OMP PRIVATE(det_minilist, minilist, haa, contrib) & + !$OMP PRIVATE(exc, h1, h2, p1, p2, s1, s2, phase, degree, ok) + do i=n_alpha,1,-1 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 + ok = .false. + + do j=1, old_det_gen + call get_excitation_degree(alpha(1,1,i), old_generators(1,1,j), ex, N_int) + if(ex <= 2) then + call get_excitation(old_generators(1,1,j), alpha(1,1,i), exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + ok = (mo_class(h1)(1:1) == 'A' .or. mo_class(h1)(1:1) == 'I') .and. & + (mo_class(p1)(1:1) == 'A' .or. mo_class(p1)(1:1) == 'V') + if(ok .and. degree == 2) then + ok = (mo_class(h2)(1:1) == 'A' .or. mo_class(h2)(1:1) == 'I') .and. & + (mo_class(p2)(1:1) == 'A' .or. mo_class(p2)(1:1) == 'V') + end if + if(ok) exit + end if + end do + + if(.not. ok) cycle + + do j=1, N_det_delta_ij call get_excitation_degree(alpha(1,1,i), psi_det(1,1,j), ex, N_int) if(ex <= 2) then n_minilist += 1 @@ -90,16 +109,15 @@ subroutine undress_with_alpha(alpha, n_alpha) 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 + call i_h_j(alpha(1,1,i), alpha(1,1,i), N_int, haa) + call dress_with_alpha_(N_states, N_det_delta_ij, N_int, delta_ij_loc(1,1,1,iproc), & + minilist, det_minilist, n_minilist, alpha(1,1,i), haa, contrib, iproc) end do !$OMP END PARALLEL DO - do i=Nproc,1,-1 + do i=1,Nproc delta_ij_tmp(:,:,:) -= delta_ij_loc(:,:,:,i) + !print *, "DELTA_IJ_LOC", delta_ij_loc(:,2:5,2,i) end do end subroutine @@ -118,7 +136,7 @@ subroutine dress_with_alpha_(Nstates,Ndet,Nint,delta_ij_loc,minilist, det_minili 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, intent(inout) :: delta_ij_loc(Nstates,Ndet,2) double precision, intent(out) :: contrib double precision, intent(in) :: haa double precision :: hij, sij