10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-16 18:25:27 +02:00

Changes in FOBOCI

This commit is contained in:
Emmanuel Giner 2016-02-19 17:32:35 +01:00
parent 87ae288904
commit 518520d619
5 changed files with 216 additions and 40 deletions

View File

@ -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

View File

@ -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)

View File

@ -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

View File

@ -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

View File

@ -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
"""