From 518520d619bb07828e75bb685513ea6c2db80f50 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 19 Feb 2016 17:32:35 +0100 Subject: [PATCH] Changes in FOBOCI --- plugins/FOBOCI/all_singles.irp.f | 4 +- plugins/FOBOCI/all_singles_split.irp.f | 44 ++++-- plugins/FOBOCI/new_approach.irp.f | 17 ++- plugins/FOBOCI/routines_dressing.irp.f | 189 ++++++++++++++++++++++--- scripts/generate_h_apply.py | 2 +- 5 files changed, 216 insertions(+), 40 deletions(-) diff --git a/plugins/FOBOCI/all_singles.irp.f b/plugins/FOBOCI/all_singles.irp.f index e2c4c01e..924e1c72 100644 --- a/plugins/FOBOCI/all_singles.irp.f +++ b/plugins/FOBOCI/all_singles.irp.f @@ -6,7 +6,7 @@ subroutine all_single double precision,allocatable :: E_before(:) N_st = N_states allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) - selection_criterion = 1.d-8 + selection_criterion = 0.d0 soft_touch selection_criterion threshold_davidson = 1.d-5 soft_touch threshold_davidson davidson_criterion @@ -79,6 +79,8 @@ subroutine all_single_no_1h_or_1p double precision,allocatable :: E_before(:) N_st = N_states allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + selection_criterion = 0.d0 + soft_touch selection_criterion threshold_davidson = 1.d-5 soft_touch threshold_davidson davidson_criterion i = 0 diff --git a/plugins/FOBOCI/all_singles_split.irp.f b/plugins/FOBOCI/all_singles_split.irp.f index e7b0943f..57f6ebde 100644 --- a/plugins/FOBOCI/all_singles_split.irp.f +++ b/plugins/FOBOCI/all_singles_split.irp.f @@ -32,46 +32,57 @@ subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,N end -subroutine all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) +subroutine all_single_for_1h(i_hole,dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_matrix_extra_1h_or_1p) implicit none use bitmasks + integer, intent(in) :: i_hole double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) double precision, intent(inout) :: dressing_matrix_2h1p(N_det_generators,N_det_generators) - integer :: i,i_hole + double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators) + integer :: i n_det_max_jacobi = 50 soft_touch n_det_max_jacobi - integer :: n_det_1h1p,n_det_2h1p + integer :: n_det_1h1p,n_det_2h1p,n_det_extra_1h_or_1p integer(bit_kind), allocatable :: psi_ref_out(:,:,:) integer(bit_kind), allocatable :: psi_1h1p(:,:,:) integer(bit_kind), allocatable :: psi_2h1p(:,:,:) + integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:) double precision, allocatable :: psi_ref_coef_out(:,:) double precision, allocatable :: psi_coef_1h1p(:,:) double precision, allocatable :: psi_coef_2h1p(:,:) - call all_single_no_1h_or_1p + double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:) +!call all_single_no_1h_or_1p + call all_single threshold_davidson = 1.d-12 soft_touch threshold_davidson davidson_criterion call diagonalize_CI - call give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p) + call give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p) allocate(psi_ref_out(N_int,2,N_det_generators)) allocate(psi_1h1p(N_int,2,n_det_1h1p)) allocate(psi_2h1p(N_int,2,n_det_2h1p)) + allocate(psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)) allocate(psi_ref_coef_out(N_det_generators,N_states)) allocate(psi_coef_1h1p(n_det_1h1p,N_states)) allocate(psi_coef_2h1p(n_det_2h1p,N_states)) - call split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p) + allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states)) + call split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & psi_1h1p,psi_coef_1h1p,n_det_1h1p) call provide_matrix_dressing_general(dressing_matrix_2h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & psi_2h1p,psi_coef_2h1p,n_det_2h1p) + call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p) deallocate(psi_ref_out) deallocate(psi_1h1p) deallocate(psi_2h1p) + deallocate(psi_extra_1h_or_1p) deallocate(psi_ref_coef_out) deallocate(psi_coef_1h1p) deallocate(psi_coef_2h1p) + deallocate(psi_coef_extra_1h_or_1p) end @@ -197,39 +208,48 @@ subroutine all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) soft_touch n_det_max_jacobi end -subroutine all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) +subroutine all_single_for_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p) implicit none use bitmasks + integer, intent(in ) :: i_particl double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators) double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators) - integer :: i,i_hole + double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators) + integer :: i n_det_max_jacobi = 50 soft_touch n_det_max_jacobi - integer :: n_det_1h1p,n_det_1h2p + integer :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p integer(bit_kind), allocatable :: psi_ref_out(:,:,:) integer(bit_kind), allocatable :: psi_1h1p(:,:,:) integer(bit_kind), allocatable :: psi_1h2p(:,:,:) + integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:) double precision, allocatable :: psi_ref_coef_out(:,:) double precision, allocatable :: psi_coef_1h1p(:,:) double precision, allocatable :: psi_coef_1h2p(:,:) - call all_single_no_1h_or_1p_or_2p + double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:) +!call all_single_no_1h_or_1p_or_2p + call all_single threshold_davidson = 1.d-12 soft_touch threshold_davidson davidson_criterion call diagonalize_CI - call give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p) + call give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p) allocate(psi_ref_out(N_int,2,N_det_generators)) allocate(psi_1h1p(N_int,2,n_det_1h1p)) allocate(psi_1h2p(N_int,2,n_det_1h2p)) + allocate(psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p)) allocate(psi_ref_coef_out(N_det_generators,N_states)) allocate(psi_coef_1h1p(n_det_1h1p,N_states)) allocate(psi_coef_1h2p(n_det_1h2p,N_states)) - call split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p) + allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states)) + call split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & psi_1h1p,psi_coef_1h1p,n_det_1h1p) call provide_matrix_dressing_general(dressing_matrix_1h2p,psi_ref_out,psi_ref_coef_out,N_det_generators, & psi_1h2p,psi_coef_1h2p,n_det_1h2p) + call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, & + psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p) deallocate(psi_ref_out) deallocate(psi_1h1p) diff --git a/plugins/FOBOCI/new_approach.irp.f b/plugins/FOBOCI/new_approach.irp.f index 49dcafc3..90b29faa 100644 --- a/plugins/FOBOCI/new_approach.irp.f +++ b/plugins/FOBOCI/new_approach.irp.f @@ -24,6 +24,7 @@ subroutine new_approach double precision, allocatable :: dressing_matrix_1h1p(:,:) double precision, allocatable :: dressing_matrix_2h1p(:,:) double precision, allocatable :: dressing_matrix_1h2p(:,:) + double precision, allocatable :: dressing_matrix_extra_1h_or_1p(:,:) double precision, allocatable :: H_matrix_tmp(:,:) logical :: verbose,is_ok @@ -81,12 +82,14 @@ subroutine new_approach ! so all the mono excitation on the new generators allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators)) allocate(dressing_matrix_2h1p(N_det_generators,N_det_generators)) + allocate(dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)) dressing_matrix_1h1p = 0.d0 dressing_matrix_2h1p = 0.d0 + dressing_matrix_extra_1h_or_1p = 0.d0 if(.not.do_it_perturbative)then n_good_hole +=1 ! call all_single_split_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) - call all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p) + call all_single_for_1h(i_hole_foboci,dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_matrix_extra_1h_or_1p) allocate(H_matrix_tmp(N_det_generators,N_det_generators)) do j = 1,N_det_generators do k = 1, N_det_generators @@ -96,7 +99,7 @@ subroutine new_approach enddo do j = 1, N_det_generators do k = 1, N_det_generators - H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_2h1p(j,k) + H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_2h1p(j,k) + dressing_matrix_extra_1h_or_1p(j,k) enddo enddo hjk = H_matrix_tmp(1,1) @@ -130,6 +133,7 @@ subroutine new_approach endif deallocate(dressing_matrix_1h1p) deallocate(dressing_matrix_2h1p) + deallocate(dressing_matrix_extra_1h_or_1p) enddo print*,'' @@ -155,12 +159,14 @@ subroutine new_approach ! so all the mono excitation on the new generators allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators)) allocate(dressing_matrix_1h2p(N_det_generators,N_det_generators)) + allocate(dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)) dressing_matrix_1h1p = 0.d0 dressing_matrix_1h2p = 0.d0 + dressing_matrix_extra_1h_or_1p = 0.d0 if(.not.do_it_perturbative)then n_good_hole +=1 ! call all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) - call all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p) + call all_single_for_1p(i_particl_osoci,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p) allocate(H_matrix_tmp(N_det_generators,N_det_generators)) do j = 1,N_det_generators do k = 1, N_det_generators @@ -170,7 +176,7 @@ subroutine new_approach enddo do j = 1, N_det_generators do k = 1, N_det_generators - H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_1h2p(j,k) + H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_1h2p(j,k) + dressing_matrix_extra_1h_or_1p(j,k) enddo enddo hjk = H_matrix_tmp(1,1) @@ -205,6 +211,7 @@ subroutine new_approach endif deallocate(dressing_matrix_1h1p) deallocate(dressing_matrix_1h2p) + deallocate(dressing_matrix_extra_1h_or_1p) enddo double precision, allocatable :: H_matrix_total(:,:) integer :: n_det_total @@ -221,7 +228,7 @@ subroutine new_approach !!! Adding the averaged dressing coming from the 1h1p that are redundant for each of the "n_good_hole" 1h H_matrix_total(i,j) += dressing_matrix_restart_1h1p(i,j)/dble(n_good_hole+n_good_particl) !!! Adding the dressing coming from the 2h1p that are not redundant for the any of CI calculations - H_matrix_total(i,j) += dressing_matrix_restart_2h1p(i,j) + H_matrix_total(i,j) += dressing_matrix_restart_2h1p(i,j) + dressing_matrix_restart_1h2p(i,j) enddo enddo do i = 1, n_good_det diff --git a/plugins/FOBOCI/routines_dressing.irp.f b/plugins/FOBOCI/routines_dressing.irp.f index 910f1109..b0edd949 100644 --- a/plugins/FOBOCI/routines_dressing.irp.f +++ b/plugins/FOBOCI/routines_dressing.irp.f @@ -81,6 +81,72 @@ subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det enddo !print*,'i_pert_count = ',i_pert_count end +subroutine provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix,psi_det_ref_input,psi_coef_ref_input,n_det_ref_input, & + psi_det_outer_input,psi_coef_outer_input,n_det_outer_input) + use bitmasks + implicit none + integer, intent(in) :: n_det_ref_input + integer(bit_kind), intent(in) :: psi_det_ref_input(N_int,2,n_det_ref_input) + double precision, intent(in) :: psi_coef_ref_input(n_det_ref_input,N_states) + integer, intent(in) :: n_det_outer_input + integer(bit_kind), intent(in) :: psi_det_outer_input(N_int,2,n_det_outer_input) + double precision, intent(in) :: psi_coef_outer_input(n_det_outer_input,N_states) + + double precision, intent(inout) :: dressing_matrix(n_det_ref_input,n_det_ref_input) + + + integer :: i_pert, i_pert_count,i,j,k + double precision :: f,href,hka,lambda_i + double precision :: H_array(n_det_ref_input),accu + integer :: n_h_out,n_p_out,n_p_in,n_h_in,number_of_holes,number_of_particles + call i_h_j(psi_det_ref_input(1,1,1),psi_det_ref_input(1,1,1),N_int,href) + i_pert_count = 0 + do i = 1, n_det_outer_input + call i_h_j(psi_det_outer_input(1,1,i),psi_det_outer_input(1,1,i),N_int,hka) + f = 1.d0/(href - hka) + H_array = 0.d0 + accu = 0.d0 +! n_h_out = number_of_holes(psi_det_outer_input(1,1,i)) +! n_p_out = number_of_particles(psi_det_outer_input(1,1,i)) + do j=1,n_det_ref_input + n_h_in = number_of_holes(psi_det_ref_input(1,1,j)) + n_p_in = number_of_particles(psi_det_ref_input(1,1,j)) +! if(n_h_in == 0 .and. n_h_in == 0)then + call i_h_j(psi_det_outer_input(1,1,i),psi_det_ref_input(1,1,j),N_int,hka) +! else +! hka = 0.d0 +! endif + H_array(j) = hka + accu += psi_coef_ref_input(j,1) * hka + enddo + lambda_i = psi_coef_outer_input(i,1)/accu + i_pert = 1 + if(accu * f / psi_coef_outer_input(i,1) .gt. 0.5d0 .and. accu * f/psi_coef_outer_input(i,1).gt.0.d0)then + i_pert = 0 + endif + do j = 1, n_det_ref_input + if(dabs(H_array(j)*lambda_i).gt.0.3d0)then + i_pert = 1 + exit + endif + enddo + if(i_pert==1)then + lambda_i = f + i_pert_count +=1 + endif + do k=1,n_det_ref_input + double precision :: contrib + contrib = H_array(k) * H_array(k) * lambda_i + dressing_matrix(k, k) += contrib + do j=k+1,n_det_ref_input + contrib = H_array(k) * H_array(j) * lambda_i + dressing_matrix(k, j) += contrib + dressing_matrix(j, k) += contrib + enddo + enddo + enddo +end + subroutine provide_matrix_dressing_general(dressing_matrix,psi_det_ref_input,psi_coef_ref_input,n_det_ref_input, & psi_det_outer_input,psi_coef_outer_input,n_det_outer_input) @@ -170,31 +236,41 @@ subroutine diag_dressed_matrix_and_set_to_psi_det(psi_det_generators_input,Ndet_ end -subroutine give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p) +subroutine give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p) use bitmasks implicit none - integer, intent(out) :: n_det_1h1p, n_det_2h1p + integer, intent(in) :: i_hole + integer, intent(out) :: n_det_1h1p, n_det_2h1p,n_det_extra_1h_or_1p integer :: i integer :: n_det_ref_restart_tmp,n_det_1h integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_hole_in_det n_det_ref_restart_tmp = 0 n_det_1h = 0 n_det_1h1p = 0 n_det_2h1p = 0 + n_det_extra_1h_or_1p = 0 do i = 1, N_det n_h = number_of_holes(psi_det(1,1,i)) n_p = number_of_particles(psi_det(1,1,i)) if(n_h == 0 .and. n_p == 0)then n_det_ref_restart_tmp +=1 else if (n_h ==1 .and. n_p==0)then - n_det_1h +=1 + if(is_the_hole_in_det(psi_det(1,1,i),1,i_hole).or.is_the_hole_in_det(psi_det(1,1,i),2,i_hole))then + n_det_1h +=1 + else + n_det_extra_1h_or_1p +=1 + endif + else if (n_h ==0 .and. n_p==1)then + n_det_extra_1h_or_1p +=1 else if (n_h ==1 .and. n_p==1)then n_det_1h1p +=1 else if (n_h ==2 .and. n_p==1)then n_det_2h1p +=1 else print*,'PB !!!!' - print*,'You have something else than a 1h, 1h1p or 2h1p' + print*,'You have something else than a 1h, 1p, 1h1p or 2h1p' + print*,'n_h,n_p = ',n_h,n_p call debug_det(psi_det(1,1,i),N_int) stop endif @@ -212,72 +288,89 @@ subroutine give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p) end -subroutine give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p) +subroutine give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p) use bitmasks implicit none - integer, intent(out) :: n_det_1h1p, n_det_1h2p + integer, intent(in) ::i_particl + integer, intent(out) :: n_det_1h1p, n_det_1h2p,n_det_extra_1h_or_1p integer :: i - integer :: n_det_ref_restart_tmp,n_det_1h + integer :: n_det_ref_restart_tmp,n_det_1p integer :: number_of_holes,n_h, number_of_particles,n_p + logical :: is_the_particl_in_det n_det_ref_restart_tmp = 0 - n_det_1h = 0 + n_det_1p = 0 n_det_1h1p = 0 n_det_1h2p = 0 + n_det_extra_1h_or_1p = 0 do i = 1, N_det n_h = number_of_holes(psi_det(1,1,i)) n_p = number_of_particles(psi_det(1,1,i)) if(n_h == 0 .and. n_p == 0)then n_det_ref_restart_tmp +=1 else if (n_h ==0 .and. n_p==1)then - n_det_1h +=1 + if(is_the_particl_in_det(psi_det(1,1,i),1,i_particl).or.is_the_particl_in_det(psi_det(1,1,i),2,i_particl))then + n_det_1p +=1 + else + n_det_extra_1h_or_1p +=1 + endif + else if (n_h ==1 .and. n_p==0)then + n_det_extra_1h_or_1p +=1 else if (n_h ==1 .and. n_p==1)then n_det_1h1p +=1 else if (n_h ==1 .and. n_p==2)then n_det_1h2p +=1 else print*,'PB !!!!' - print*,'You have something else than a 1p, 1h1p or 1h2p' + print*,'You have something else than a 1h, 1p, 1h1p or 1h2p' call debug_det(psi_det(1,1,i),N_int) stop endif enddo - if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then - print*,'PB !!!!' - print*,'You have forgotten something in your generators ... ' - stop - endif +!if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then +! print*,'PB !!!!' +! print*,'You have forgotten something in your generators ... ' +! stop +!endif end -subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p) +subroutine split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) use bitmasks implicit none - integer, intent(in) :: n_det_1h1p,n_det_2h1p + integer, intent(in) :: n_det_1h1p,n_det_2h1p,n_det_extra_1h_or_1p,i_hole integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators) integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p) integer(bit_kind), intent(out) :: psi_2h1p(N_int,2,n_det_2h1p) + integer(bit_kind), intent(out) :: psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p) double precision, intent(out) :: psi_ref_coef_out(N_det_generators,N_states) double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states) double precision, intent(out) :: psi_coef_2h1p(n_det_2h1p, N_states) + double precision, intent(out) :: psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p, N_states) integer :: i,j integer :: degree integer :: number_of_holes,n_h, number_of_particles,n_p - integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_2h1p_tmp + integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_2h1p_tmp,n_det_extra_1h_or_1p_tmp + integer :: n_det_extra_1h_tmp integer, allocatable :: index_generator(:) integer, allocatable :: index_1h1p(:) integer, allocatable :: index_2h1p(:) + integer, allocatable :: index_extra_1h_or_1p(:) + logical :: is_the_hole_in_det allocate(index_1h1p(n_det)) allocate(index_2h1p(n_det)) + allocate(index_extra_1h_or_1p(n_det)) allocate(index_generator(N_det)) n_det_generators_tmp = 0 n_det_1h1p_tmp = 0 n_det_2h1p_tmp = 0 + n_det_extra_1h_or_1p_tmp = 0 + n_det_extra_1h_tmp = 0 do i = 1, n_det n_h = number_of_holes(psi_det(1,1,i)) n_p = number_of_particles(psi_det(1,1,i)) @@ -287,6 +380,16 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_o else if (n_h ==2 .and. n_p==1)then n_det_2h1p_tmp +=1 index_2h1p(n_det_2h1p_tmp) = i + else if (n_h ==0 .and. n_p==1)then + n_det_extra_1h_or_1p_tmp +=1 + index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i + else if (n_h ==1 .and. n_p==0)then + if(is_the_hole_in_det(psi_det(1,1,i),1,i_hole).or.is_the_hole_in_det(psi_det(1,1,i),2,i_hole))then + n_det_extra_1h_tmp +=1 + else + n_det_extra_1h_or_1p_tmp +=1 + index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i + endif endif do j = 1, N_det_generators call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int) @@ -315,6 +418,12 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_o stop endif + if(n_det_extra_1h_or_1p.ne.n_det_extra_1h_or_1p_tmp)then + print*,'PB !!!' + print*,'n_det_extra_1h_or_1p.ne.n_det_extra_1h_or_1p_tmp' + stop + endif + do i = 1,N_det_generators do j = 1, N_int psi_ref_out(j,1,i) = psi_det(j,1,index_generator(i)) @@ -345,41 +454,59 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_o enddo enddo + do i = 1, n_det_extra_1h_or_1p + do j = 1, N_int + psi_extra_1h_or_1p(j,1,i) = psi_det(j,1,index_extra_1h_or_1p(i)) + psi_extra_1h_or_1p(j,2,i) = psi_det(j,2,index_extra_1h_or_1p(i)) + enddo + do j = 1, N_states + psi_coef_extra_1h_or_1p(i,j) = psi_coef(index_extra_1h_or_1p(i),j) + enddo + enddo deallocate(index_generator) deallocate(index_1h1p) deallocate(index_2h1p) + deallocate(index_extra_1h_or_1p) end -subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p) +subroutine split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p) use bitmasks implicit none - integer, intent(in) :: n_det_1h1p,n_det_1h2p + integer, intent(in) :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p,i_particl integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators) integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p) integer(bit_kind), intent(out) :: psi_1h2p(N_int,2,n_det_1h2p) + integer(bit_kind), intent(out) :: psi_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p) double precision, intent(out) :: psi_ref_coef_out(N_det_generators,N_states) double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states) double precision, intent(out) :: psi_coef_1h2p(n_det_1h2p, N_states) + double precision, intent(out) :: psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p, N_states) integer :: i,j integer :: degree integer :: number_of_holes,n_h, number_of_particles,n_p - integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_1h2p_tmp + integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_1h2p_tmp,n_det_extra_1h_or_1p_tmp integer, allocatable :: index_generator(:) integer, allocatable :: index_1h1p(:) integer, allocatable :: index_1h2p(:) + integer, allocatable :: index_extra_1h_or_1p(:) + logical :: is_the_particl_in_det + integer :: n_det_1p_tmp allocate(index_1h1p(n_det)) allocate(index_1h2p(n_det)) + allocate(index_extra_1h_or_1p(n_det)) allocate(index_generator(N_det)) n_det_generators_tmp = 0 n_det_1h1p_tmp = 0 n_det_1h2p_tmp = 0 + n_det_extra_1h_or_1p_tmp = 0 + n_det_1p_tmp = 0 do i = 1, n_det n_h = number_of_holes(psi_det(1,1,i)) n_p = number_of_particles(psi_det(1,1,i)) @@ -389,6 +516,15 @@ subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_o else if (n_h ==1 .and. n_p==2)then n_det_1h2p_tmp +=1 index_1h2p(n_det_1h2p_tmp) = i + else if (n_h ==1 .and. n_p==0)then + n_det_extra_1h_or_1p_tmp +=1 + index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i + else if (n_h ==0 .and. n_p==1)then + if(is_the_particl_in_det(psi_det(1,1,i),1,i_particl).or.is_the_particl_in_det(psi_det(1,1,i),2,i_particl))then + n_det_1p_tmp +=1 + else + n_det_extra_1h_or_1p_tmp +=1 + endif endif do j = 1, N_det_generators call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int) @@ -448,9 +584,20 @@ subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_o enddo + do i = 1, n_det_extra_1h_or_1p + do j = 1, N_int + psi_extra_1h_or_1p(j,1,i) = psi_det(j,1,index_extra_1h_or_1p(i)) + psi_extra_1h_or_1p(j,2,i) = psi_det(j,2,index_extra_1h_or_1p(i)) + enddo + do j = 1, N_states + psi_coef_extra_1h_or_1p(i,j) = psi_coef(index_extra_1h_or_1p(i),j) + enddo + enddo + deallocate(index_generator) deallocate(index_1h1p) deallocate(index_1h2p) + deallocate(index_extra_1h_or_1p) end diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index ca2be5d6..f9be4fa6 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -168,7 +168,7 @@ class H_apply(object): if (is_a_2p(hole)) cycle """ def filter_1p(self): - self["filter0p"] = """ + self["filter1p"] = """ ! ! DIR$ FORCEINLINE if (is_a_1p(hole)) cycle """