mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 04:43:50 +01:00
887 lines
28 KiB
Fortran
887 lines
28 KiB
Fortran
subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input)
|
|
use bitmasks
|
|
implicit none
|
|
integer, intent(in) :: ndet_generators_input
|
|
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,ndet_generators_input)
|
|
double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input)
|
|
double precision :: H_array(N_det),hka
|
|
logical :: is_a_ref_det(N_det)
|
|
integer :: i,j,n_det_ref_tmp
|
|
integer :: connected_to_ref_by_mono,degree
|
|
double precision :: coef_ref(Ndet_generators_input)
|
|
double precision :: accu,lambda_i
|
|
integer :: k
|
|
integer :: index_ref_tmp(N_det)
|
|
is_a_ref_det = .False.
|
|
n_det_ref_tmp = 0
|
|
do i = 1, N_det
|
|
do j = 1, Ndet_generators_input
|
|
call get_excitation_degree(psi_det(1,1,i),psi_det_generators_input(1,1,j),degree,N_int)
|
|
if(degree == 0)then
|
|
is_a_ref_det(i) = .True.
|
|
n_det_ref_tmp +=1
|
|
index_ref_tmp(n_det_ref_tmp) = i
|
|
coef_ref(n_det_ref_tmp) = psi_coef(i,1)
|
|
exit
|
|
endif
|
|
enddo
|
|
enddo
|
|
if( ndet_generators_input .ne. n_det_ref_tmp)then
|
|
print*,'Problem !!!! '
|
|
print*,' ndet_generators .ne. n_det_ref_tmp !!!'
|
|
print*,'ndet_generators,n_det_ref_tmp'
|
|
print*,ndet_generators_input,n_det_ref_tmp
|
|
stop
|
|
endif
|
|
|
|
call i_h_j(psi_det_generators_input(1,1,1),psi_det_generators_input(1,1,1),N_int,href)
|
|
integer :: i_pert, i_pert_count
|
|
i_pert_count = 0
|
|
do i = 1, N_det
|
|
if(is_a_ref_det(i))cycle
|
|
call i_h_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hka)
|
|
double precision :: f,href
|
|
f = 1.d0/(href - hka)
|
|
H_array = 0.d0
|
|
accu = 0.d0
|
|
do j=1,ndet_generators_input
|
|
call i_h_j(psi_det(1,1,i),psi_det_generators_input(1,1,j),N_int,hka)
|
|
H_array(j) = hka
|
|
accu += coef_ref(j) * hka
|
|
enddo
|
|
lambda_i = psi_coef(i,1)/accu
|
|
i_pert = 1
|
|
if(accu * f / psi_coef(i,1) .gt. 0.5d0 .and. accu * f/psi_coef(i,1).gt.0.d0)then
|
|
i_pert = 0
|
|
endif
|
|
do j = 1, ndet_generators_input
|
|
if(dabs(H_array(j)*lambda_i).gt.0.1d0)then
|
|
i_pert = 1
|
|
exit
|
|
endif
|
|
enddo
|
|
if(i_pert==1)then
|
|
lambda_i = f
|
|
i_pert_count +=1
|
|
endif
|
|
do k=1,ndet_generators_input
|
|
double precision :: contrib
|
|
contrib = H_array(k) * H_array(k) * lambda_i
|
|
dressing_matrix(k, k) += contrib
|
|
do j=k+1,ndet_generators_input
|
|
contrib = H_array(k) * H_array(j) * lambda_i
|
|
dressing_matrix(k, j) += contrib
|
|
dressing_matrix(j, k) += contrib
|
|
enddo
|
|
enddo
|
|
enddo
|
|
href = dressing_matrix(1,1)
|
|
print*,'Diagonal part of the dressing'
|
|
do i = 1, ndet_generators_input
|
|
print*,'delta e = ',dressing_matrix(i,i) - href
|
|
enddo
|
|
!print*,'i_pert_count = ',i_pert_count
|
|
end
|
|
|
|
subroutine update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,H_jj_in)
|
|
use bitmasks
|
|
implicit none
|
|
integer, intent(in) :: ndet_generators_input
|
|
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,ndet_generators_input)
|
|
double precision, intent(in) :: H_jj_in(N_det)
|
|
double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input)
|
|
integer :: i,j,n_det_ref_tmp,degree
|
|
double precision :: href
|
|
n_det_ref_tmp = 0
|
|
do i = 1, N_det
|
|
do j = 1, Ndet_generators_input
|
|
call get_excitation_degree(psi_det(1,1,i),psi_det_generators_input(1,1,j),degree,N_int)
|
|
if(degree == 0)then
|
|
dressing_matrix(j,j) += H_jj_in(i)
|
|
n_det_ref_tmp +=1
|
|
exit
|
|
endif
|
|
enddo
|
|
enddo
|
|
if( ndet_generators_input .ne. n_det_ref_tmp)then
|
|
print*,'Problem !!!! '
|
|
print*,' ndet_generators .ne. n_det_ref_tmp !!!'
|
|
print*,'ndet_generators,n_det_ref_tmp'
|
|
print*,ndet_generators_input,n_det_ref_tmp
|
|
stop
|
|
endif
|
|
|
|
href = dressing_matrix(1,1)
|
|
print*,''
|
|
print*,'Update with the SC2 dressing'
|
|
print*,''
|
|
print*,'Diagonal part of the dressing'
|
|
do i = 1, ndet_generators_input
|
|
print*,'delta e = ',dressing_matrix(i,i) - href
|
|
enddo
|
|
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.5d0)then
|
|
i_pert = 1
|
|
exit
|
|
endif
|
|
enddo
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
! i_pert = 0
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
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)
|
|
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
|
|
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
|
|
do j=1,n_det_ref_input
|
|
call i_h_j(psi_det_outer_input(1,1,i),psi_det_ref_input(1,1,j),N_int,hka)
|
|
H_array(j) = hka
|
|
accu += psi_coef_ref_input(j,1) * hka
|
|
enddo
|
|
lambda_i = psi_coef_outer_input(i,1)/accu
|
|
i_pert = 0
|
|
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.5d0)then
|
|
i_pert = 1
|
|
exit
|
|
endif
|
|
enddo
|
|
! i_pert = 0
|
|
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 diag_dressed_matrix_and_set_to_psi_det(psi_det_generators_input,Ndet_generators_input,dressing_matrix)
|
|
use bitmasks
|
|
implicit none
|
|
integer, intent(in) :: ndet_generators_input
|
|
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,ndet_generators_input)
|
|
double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input)
|
|
integer :: i,j
|
|
|
|
double precision :: eigenvectors(Ndet_generators_input,Ndet_generators_input), eigenvalues(Ndet_generators_input)
|
|
|
|
call lapack_diag(eigenvalues,eigenvectors,dressing_matrix,Ndet_generators_input,Ndet_generators_input)
|
|
print*,'Dressed eigenvalue, to be compared with the CI one'
|
|
print*,'E = ',eigenvalues(1)+nuclear_repulsion
|
|
print*,'Dressed matrix, to be compared to the intermediate Hamiltonian one'
|
|
do i = 1, Ndet_generators_input
|
|
write(*,'(100(F12.5,X))')dressing_matrix(i,:)
|
|
enddo
|
|
n_det = Ndet_generators_input
|
|
do i = 1, Ndet_generators_input
|
|
psi_coef(i,1) = eigenvectors(i,1)
|
|
do j = 1, N_int
|
|
psi_det(j,1,i) = psi_det_generators_input(j,1,i)
|
|
psi_det(j,2,i) = psi_det_generators_input(j,2,i)
|
|
enddo
|
|
enddo
|
|
|
|
touch N_det psi_coef psi_det
|
|
|
|
end
|
|
|
|
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(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
|
|
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, 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
|
|
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_2h1p + n_det_1h1p + n_det_extra_1h_or_1p + n_det_generators .ne. N_det)then
|
|
print*,'PB !!!!'
|
|
print*,'You have forgotten something in your generators ... '
|
|
stop
|
|
endif
|
|
|
|
end
|
|
|
|
subroutine give_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p)
|
|
use bitmasks
|
|
implicit none
|
|
integer, intent(out) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p
|
|
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_1h_1p = 0
|
|
n_det_2h1p = 0
|
|
n_det_1h1p = 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_1h_1p +=1
|
|
else if (n_h ==1 .and. n_p==0)then
|
|
n_det_ref_1h_1p +=1
|
|
else if (n_h ==0 .and. n_p==1)then
|
|
n_det_ref_1h_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, 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
|
|
enddo
|
|
|
|
end
|
|
|
|
subroutine give_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p)
|
|
use bitmasks
|
|
implicit none
|
|
integer, intent(out) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p
|
|
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_1h_1p = 0
|
|
n_det_1h2p = 0
|
|
n_det_1h1p = 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_1h_1p +=1
|
|
else if (n_h ==1 .and. n_p==0)then
|
|
n_det_ref_1h_1p +=1
|
|
else if (n_h ==0 .and. n_p==1)then
|
|
n_det_ref_1h_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 1h, 1p, 1h1p or 1h2p'
|
|
print*,'n_h,n_p = ',n_h,n_p
|
|
call debug_det(psi_det(1,1,i),N_int)
|
|
stop
|
|
endif
|
|
enddo
|
|
|
|
end
|
|
|
|
subroutine give_wf_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,&
|
|
psi_det_2h1p,psi_coef_2h1p,psi_det_1h1p,psi_coef_1h1p)
|
|
use bitmasks
|
|
implicit none
|
|
integer, intent(in) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p
|
|
integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p)
|
|
integer(bit_kind), intent(out) :: psi_det_2h1p(N_int,2,n_det_2h1p)
|
|
integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p)
|
|
double precision, intent(out) :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)
|
|
double precision, intent(out) :: psi_coef_2h1p(n_det_2h1p,N_states)
|
|
double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p,N_states)
|
|
integer :: n_det_ref_1h_1p_tmp,n_det_2h1p_tmp,n_det_1h1p_tmp
|
|
integer :: i,j
|
|
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
|
|
integer, allocatable :: index_ref_1h_1p(:)
|
|
integer, allocatable :: index_2h1p(:)
|
|
integer, allocatable :: index_1h1p(:)
|
|
allocate(index_ref_1h_1p(n_det))
|
|
allocate(index_2h1p(n_det))
|
|
allocate(index_1h1p(n_det))
|
|
n_det_ref_1h_1p_tmp = 0
|
|
n_det_2h1p_tmp = 0
|
|
n_det_1h1p_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))
|
|
if(n_h == 0 .and. n_p == 0)then
|
|
n_det_ref_1h_1p_tmp +=1
|
|
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
|
|
else if (n_h ==1 .and. n_p==0)then
|
|
n_det_ref_1h_1p_tmp +=1
|
|
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
|
|
else if (n_h ==0 .and. n_p==1)then
|
|
n_det_ref_1h_1p_tmp +=1
|
|
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
|
|
else if (n_h ==1 .and. n_p==1)then
|
|
n_det_1h1p_tmp +=1
|
|
index_1h1p(n_det_1h1p_tmp) = i
|
|
else if (n_h ==2 .and. n_p==1)then
|
|
n_det_2h1p_tmp +=1
|
|
index_2h1p(n_det_2h1p_tmp) = i
|
|
else
|
|
print*,'PB !!!!'
|
|
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
|
|
enddo
|
|
do i = 1, n_det_2h1p
|
|
do j = 1, N_int
|
|
psi_det_2h1p(j,1,i) = psi_det(j,1,index_2h1p(i))
|
|
psi_det_2h1p(j,2,i) = psi_det(j,2,index_2h1p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_2h1p(i,j) = psi_coef(index_2h1p(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
do i = 1, n_det_1h1p
|
|
do j = 1, N_int
|
|
psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
|
|
psi_det_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
do i = 1, n_det_ref_1h_1p
|
|
do j = 1, N_int
|
|
psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i))
|
|
psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
end
|
|
|
|
subroutine give_wf_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,&
|
|
psi_det_1h2p,psi_coef_1h2p,psi_det_1h1p,psi_coef_1h1p)
|
|
use bitmasks
|
|
implicit none
|
|
integer, intent(in) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p
|
|
integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p)
|
|
integer(bit_kind), intent(out) :: psi_det_1h2p(N_int,2,n_det_1h2p)
|
|
integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p)
|
|
double precision, intent(out) :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)
|
|
double precision, intent(out) :: psi_coef_1h2p(n_det_1h2p,N_states)
|
|
double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p,N_states)
|
|
integer :: n_det_ref_1h_1p_tmp,n_det_1h2p_tmp,n_det_1h1p_tmp
|
|
integer :: i,j
|
|
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
|
|
integer, allocatable :: index_ref_1h_1p(:)
|
|
integer, allocatable :: index_1h2p(:)
|
|
integer, allocatable :: index_1h1p(:)
|
|
allocate(index_ref_1h_1p(n_det))
|
|
allocate(index_1h2p(n_det))
|
|
allocate(index_1h1p(n_det))
|
|
n_det_ref_1h_1p_tmp = 0
|
|
n_det_1h2p_tmp = 0
|
|
n_det_1h1p_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))
|
|
if(n_h == 0 .and. n_p == 0)then
|
|
n_det_ref_1h_1p_tmp +=1
|
|
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
|
|
else if (n_h ==1 .and. n_p==0)then
|
|
n_det_ref_1h_1p_tmp +=1
|
|
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
|
|
else if (n_h ==0 .and. n_p==1)then
|
|
n_det_ref_1h_1p_tmp +=1
|
|
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
|
|
else if (n_h ==1 .and. n_p==1)then
|
|
n_det_1h1p_tmp +=1
|
|
index_1h1p(n_det_1h1p_tmp) = i
|
|
else if (n_h ==1 .and. n_p==2)then
|
|
n_det_1h2p_tmp +=1
|
|
index_1h2p(n_det_1h2p_tmp) = i
|
|
else
|
|
print*,'PB !!!!'
|
|
print*,'You have something else than a 1h, 1p, 1h1p or 1h2p'
|
|
print*,'n_h,n_p = ',n_h,n_p
|
|
call debug_det(psi_det(1,1,i),N_int)
|
|
stop
|
|
endif
|
|
enddo
|
|
do i = 1, n_det_1h2p
|
|
do j = 1, N_int
|
|
psi_det_1h2p(j,1,i) = psi_det(j,1,index_1h2p(i))
|
|
psi_det_1h2p(j,2,i) = psi_det(j,2,index_1h2p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_1h2p(i,j) = psi_coef(index_1h2p(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
do i = 1, n_det_1h1p
|
|
do j = 1, N_int
|
|
psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
|
|
psi_det_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
do i = 1, n_det_ref_1h_1p
|
|
do j = 1, N_int
|
|
psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i))
|
|
psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
end
|
|
|
|
|
|
|
|
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(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_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_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
|
|
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 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
|
|
|
|
|
|
end
|
|
|
|
|
|
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,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,n_det_extra_1h_or_1p_tmp
|
|
integer :: n_det_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_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))
|
|
if (n_h ==1 .and. n_p==1)then
|
|
n_det_1h1p_tmp +=1
|
|
index_1h1p(n_det_1h1p_tmp) = i
|
|
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_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)
|
|
if(degree == 0)then
|
|
n_det_generators_tmp +=1
|
|
index_generator(n_det_generators_tmp) = i
|
|
endif
|
|
enddo
|
|
enddo
|
|
if(n_det_1h1p_tmp.ne.n_det_1h1p)then
|
|
print*,'PB !!!'
|
|
print*,'n_det_1h1p_tmp.ne.n_det_1h1p)'
|
|
stop
|
|
endif
|
|
|
|
|
|
if(n_det_2h1p_tmp.ne.n_det_2h1p)then
|
|
print*,'PB !!!'
|
|
print*,'n_det_2h1p_tmp.ne.n_det_2h1p)'
|
|
stop
|
|
endif
|
|
|
|
if(N_det_generators_tmp.ne.n_det_generators)then
|
|
print*,'PB !!!'
|
|
print*,'N_det_generators_tmp.ne.n_det_generators'
|
|
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))
|
|
psi_ref_out(j,2,i) = psi_det(j,2,index_generator(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_ref_coef_out(i,j) = psi_coef(index_generator(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
do i = 1, n_det_1h1p
|
|
do j = 1, N_int
|
|
psi_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
|
|
psi_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
do i = 1, n_det_2h1p
|
|
do j = 1, N_int
|
|
psi_2h1p(j,1,i) = psi_det(j,1,index_2h1p(i))
|
|
psi_2h1p(j,2,i) = psi_det(j,2,index_2h1p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_2h1p(i,j) = psi_coef(index_2h1p(i),j)
|
|
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(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,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,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))
|
|
if (n_h ==1 .and. n_p==1)then
|
|
n_det_1h1p_tmp +=1
|
|
index_1h1p(n_det_1h1p_tmp) = i
|
|
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)
|
|
if(degree == 0)then
|
|
n_det_generators_tmp +=1
|
|
index_generator(n_det_generators_tmp) = i
|
|
endif
|
|
enddo
|
|
enddo
|
|
if(n_det_1h1p_tmp.ne.n_det_1h1p)then
|
|
print*,'PB !!!'
|
|
print*,'n_det_1h1p_tmp.ne.n_det_1h1p)'
|
|
stop
|
|
endif
|
|
|
|
|
|
if(n_det_1h2p_tmp.ne.n_det_1h2p)then
|
|
print*,'PB !!!'
|
|
print*,'n_det_1h2p_tmp.ne.n_det_1h2p)'
|
|
stop
|
|
endif
|
|
|
|
if(N_det_generators_tmp.ne.n_det_generators)then
|
|
print*,'PB !!!'
|
|
print*,'N_det_generators_tmp.ne.n_det_generators'
|
|
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))
|
|
psi_ref_out(j,2,i) = psi_det(j,2,index_generator(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_ref_coef_out(i,j) = psi_coef(index_generator(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
do i = 1, n_det_1h1p
|
|
do j = 1, N_int
|
|
psi_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
|
|
psi_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j)
|
|
enddo
|
|
enddo
|
|
|
|
do i = 1, n_det_1h2p
|
|
do j = 1, N_int
|
|
psi_1h2p(j,1,i) = psi_det(j,1,index_1h2p(i))
|
|
psi_1h2p(j,2,i) = psi_det(j,2,index_1h2p(i))
|
|
enddo
|
|
do j = 1, N_states
|
|
psi_coef_1h2p(i,j) = psi_coef(index_1h2p(i),j)
|
|
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_1h2p)
|
|
deallocate(index_extra_1h_or_1p)
|
|
|
|
end
|
|
|
|
|