LCPQ
/
quantum_package
mirror of https://github.com/LCPQ/quantum_package
10
0
Fork 0
Code Issues Releases Wiki Activity
quantum_package/src/Determinants/spindeterminants.irp.f

543 lines
14 KiB
Fortran
Raw Blame History

!==============================================================================!
! !
! Independent alpha/beta parts !
! !
!==============================================================================!
use bitmasks
integer*8 function spin_det_search_key(det,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Return an integer(8) corresponding to a determinant index for searching
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det(Nint)
integer(bit_kind), parameter :: unsigned_shift = -huge(1_bit_kind) ! 100...00
integer :: i
spin_det_search_key = det(1)
do i=2,Nint
spin_det_search_key = ieor(spin_det_search_key,det(i))
enddo
spin_det_search_key = spin_det_search_key+1_bit_kind-unsigned_shift
end
BEGIN_PROVIDER [ integer(bit_kind), psi_det_alpha, (N_int,psi_det_size) ]
implicit none
BEGIN_DOC
! List of alpha determinants of psi_det
END_DOC
integer :: i,k
do i=1,N_det
do k=1,N_int
psi_det_alpha(k,i) = psi_det(k,1,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_beta, (N_int,psi_det_size) ]
implicit none
BEGIN_DOC
! List of beta determinants of psi_det
END_DOC
integer :: i,k
do i=1,N_det
do k=1,N_int
psi_det_beta(k,i) = psi_det(k,2,i)
enddo
enddo
END_PROVIDER
BEGIN_TEMPLATE
BEGIN_PROVIDER [ integer(bit_kind), psi_det_$alpha_unique, (N_int,psi_det_size) ]
&BEGIN_PROVIDER [ integer, N_det_$alpha_unique ]
implicit none
BEGIN_DOC
! Unique $alpha determinants
END_DOC
integer :: i,j,k
integer, allocatable :: iorder(:)
integer(8), allocatable :: bit_tmp(:)
integer(8) :: last_key
integer(8), external :: spin_det_search_key
logical,allocatable :: duplicate(:)
allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) )
do i=1,N_det
iorder(i) = i
bit_tmp(i) = spin_det_search_key(psi_det_$alpha(1,i),N_int)
enddo
call i8sort(bit_tmp,iorder,N_det)
N_det_$alpha_unique = 0
last_key = 0_8
do i=1,N_det
last_key = bit_tmp(i)
N_det_$alpha_unique += 1
do k=1,N_int
psi_det_$alpha_unique(k,N_det_$alpha_unique) = psi_det_$alpha(k,iorder(i))
enddo
duplicate(i) = .False.
enddo
j=1
do i=1,N_det_$alpha_unique-1
if (duplicate(i)) then
cycle
endif
j = i+1
do while (bit_tmp(j)==bit_tmp(i))
if (duplicate(j)) then
j += 1
cycle
endif
duplicate(j) = .True.
do k=1,N_int
if (psi_det_$alpha_unique(k,i) /= psi_det_$alpha_unique(k,j)) then
duplicate(j) = .False.
exit
endif
enddo
j+=1
if (j > N_det_$alpha_unique) then
exit
endif
enddo
enddo
j=1
do i=2,N_det_$alpha_unique
if (duplicate(i)) then
cycle
else
j += 1
psi_det_$alpha_unique(:,j) = psi_det_$alpha_unique(:,i)
endif
enddo
N_det_$alpha_unique = j
deallocate (iorder, bit_tmp, duplicate)
END_PROVIDER
SUBST [ alpha ]
alpha ;;
beta ;;
END_TEMPLATE
integer function get_index_in_psi_det_alpha_unique(key,Nint)
use bitmasks
BEGIN_DOC
! Returns the index of the determinant in the ``psi_det_alpha_unique`` array
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key(Nint)
integer :: i, ibegin, iend, istep, l
integer(8) :: det_ref, det_search
integer(8), external :: spin_det_search_key
logical :: in_wavefunction
in_wavefunction = .False.
get_index_in_psi_det_alpha_unique = 0
ibegin = 1
iend = N_det_alpha_unique + 1
!DIR$ FORCEINLINE
det_ref = spin_det_search_key(key,Nint)
!DIR$ FORCEINLINE
det_search = spin_det_search_key(psi_det_alpha_unique(1,1),Nint)
istep = ishft(iend-ibegin,-1)
i=ibegin+istep
do while (istep > 0)
!DIR$ FORCEINLINE
det_search = spin_det_search_key(psi_det_alpha_unique(1,i),Nint)
if ( det_search > det_ref ) then
iend = i
else if ( det_search == det_ref ) then
exit
else
ibegin = i
endif
istep = ishft(iend-ibegin,-1)
i = ibegin + istep
end do
!DIR$ FORCEINLINE
do while (spin_det_search_key(psi_det_alpha_unique(1,i),Nint) == det_ref)
i = i-1
if (i == 0) then
exit
endif
enddo
i += 1
if (i > N_det_alpha_unique) then
return
endif
!DIR$ FORCEINLINE
do while (spin_det_search_key(psi_det_alpha_unique(1,i),Nint) == det_ref)
if (key(1) /= psi_det_alpha_unique(1,i)) then
continue
else
in_wavefunction = .True.
!DIR$ IVDEP
!DIR$ LOOP COUNT MIN(3)
do l=2,Nint
if (key(l) /= psi_det_alpha_unique(l,i)) then
in_wavefunction = .False.
endif
enddo
if (in_wavefunction) then
get_index_in_psi_det_alpha_unique = i
return
endif
endif
i += 1
if (i > N_det_alpha_unique) then
return
endif
enddo
end
integer function get_index_in_psi_det_beta_unique(key,Nint)
use bitmasks
BEGIN_DOC
! Returns the index of the determinant in the ``psi_det_beta_unique`` array
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key(Nint)
integer :: i, ibegin, iend, istep, l
integer(8) :: det_ref, det_search
integer(8), external :: spin_det_search_key
logical :: in_wavefunction
in_wavefunction = .False.
get_index_in_psi_det_beta_unique = 0
ibegin = 1
iend = N_det_beta_unique + 1
!DIR$ FORCEINLINE
det_ref = spin_det_search_key(key,Nint)
!DIR$ FORCEINLINE
det_search = spin_det_search_key(psi_det_beta_unique(1,1),Nint)
istep = ishft(iend-ibegin,-1)
i=ibegin+istep
do while (istep > 0)
!DIR$ FORCEINLINE
det_search = spin_det_search_key(psi_det_beta_unique(1,i),Nint)
if ( det_search > det_ref ) then
iend = i
else if ( det_search == det_ref ) then
exit
else
ibegin = i
endif
istep = ishft(iend-ibegin,-1)
i = ibegin + istep
end do
!DIR$ FORCEINLINE
do while (spin_det_search_key(psi_det_beta_unique(1,i),Nint) == det_ref)
i = i-1
if (i == 0) then
exit
endif
enddo
i += 1
if (i > N_det_beta_unique) then
return
endif
!DIR$ FORCEINLINE
do while (spin_det_search_key(psi_det_beta_unique(1,i),Nint) == det_ref)
if (key(1) /= psi_det_beta_unique(1,i)) then
continue
else
in_wavefunction = .True.
!DIR$ IVDEP
!DIR$ LOOP COUNT MIN(3)
do l=2,Nint
if (key(l) /= psi_det_beta_unique(l,i)) then
in_wavefunction = .False.
endif
enddo
if (in_wavefunction) then
get_index_in_psi_det_beta_unique = i
return
endif
endif
i += 1
if (i > N_det_beta_unique) then
return
endif
enddo
end
subroutine write_spindeterminants
use bitmasks
implicit none
integer(8), allocatable :: tmpdet(:,:)
integer :: N_int2
integer :: i,j,k
integer(8) :: det_8(100)
integer(bit_kind) :: det_bk((100*8)/bit_kind)
equivalence (det_8, det_bk)
N_int2 = (N_int*bit_kind)/8
call ezfio_set_spindeterminants_n_det_alpha(N_det_alpha_unique)
call ezfio_set_spindeterminants_n_det_beta(N_det_beta_unique)
call ezfio_set_spindeterminants_n_det(N_det)
call ezfio_set_spindeterminants_n_int(N_int)
call ezfio_set_spindeterminants_bit_kind(bit_kind)
call ezfio_set_spindeterminants_n_states(N_states)
allocate(tmpdet(N_int2,N_det_alpha_unique))
do i=1,N_det_alpha_unique
do k=1,N_int
det_bk(k) = psi_det_alpha_unique(k,i)
enddo
do k=1,N_int2
tmpdet(k,i) = det_8(k)
enddo
enddo
call ezfio_set_spindeterminants_psi_det_alpha(psi_det_alpha_unique)
deallocate(tmpdet)
allocate(tmpdet(N_int2,N_det_beta_unique))
do i=1,N_det_beta_unique
do k=1,N_int
det_bk(k) = psi_det_beta_unique(k,i)
enddo
do k=1,N_int2
tmpdet(k,i) = det_8(k)
enddo
enddo
call ezfio_set_spindeterminants_psi_det_beta(psi_det_beta_unique)
deallocate(tmpdet)
call ezfio_set_spindeterminants_psi_coef_matrix_values(psi_bilinear_matrix_values)
call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_bilinear_matrix_rows)
call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_bilinear_matrix_columns)
end
BEGIN_PROVIDER [ double precision, det_alpha_norm, (N_det_alpha_unique) ]
&BEGIN_PROVIDER [ double precision, det_beta_norm, (N_det_beta_unique) ]
implicit none
BEGIN_DOC
! Norm of the alpha and beta spin determinants in the wave function:
!
! ||Da||_i \sum_j C_{ij}**2
END_DOC
integer :: i,j,k,l
double precision :: f
det_alpha_norm = 0.d0
det_beta_norm = 0.d0
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
do l=1,N_states
f = psi_bilinear_matrix_values(k,l)*psi_bilinear_matrix_values(k,l)
enddo
det_alpha_norm(i) += f
det_beta_norm(j) += f
enddo
det_alpha_norm = det_alpha_norm / dble(N_states)
det_beta_norm = det_beta_norm / dble(N_states)
END_PROVIDER
!==============================================================================!
! !
! Alpha x Beta Matrix !
! !
!==============================================================================!
BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) ]
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows, (N_det) ]
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns, (N_det) ]
use bitmasks
implicit none
BEGIN_DOC
! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
! D_a^t C D_b
END_DOC
integer :: i,j,k, l
integer(bit_kind) :: tmp_det(N_int,2)
integer :: idx
integer, external :: get_index_in_psi_det_sorted_bit
PROVIDE psi_coef_sorted_bit
integer, allocatable :: iorder(:), to_sort(:)
integer, external :: get_index_in_psi_det_alpha_unique
integer, external :: get_index_in_psi_det_beta_unique
allocate(iorder(N_det), to_sort(N_det))
do k=1,N_det
i = get_index_in_psi_det_alpha_unique(psi_det(1,1,k),N_int)
j = get_index_in_psi_det_beta_unique (psi_det(1,2,k),N_int)
do l=1,N_states
psi_bilinear_matrix_values(k,l) = psi_coef(k,l)
enddo
psi_bilinear_matrix_rows(k) = i
psi_bilinear_matrix_columns(k) = j
to_sort(k) = N_det_alpha_unique * (j-1) + i
iorder(k) = k
enddo
call isort(to_sort, iorder, N_det)
call iset_order(psi_bilinear_matrix_rows,iorder,N_det)
call iset_order(psi_bilinear_matrix_columns,iorder,N_det)
call dset_order(psi_bilinear_matrix_values,iorder,N_det)
deallocate(iorder,to_sort)
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_bilinear_matrix, (N_det_alpha_unique,N_det_beta_unique,N_states) ]
implicit none
BEGIN_DOC
! Coefficient matrix if the wave function is expressed in a bilinear form :
! D_a^t C D_b
END_DOC
integer :: i,j,k,istate
psi_bilinear_matrix = 0.d0
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
do istate=1,N_states
psi_bilinear_matrix(i,j,istate) = psi_bilinear_matrix_values(k,istate)
enddo
enddo
END_PROVIDER
subroutine create_wf_of_psi_bilinear_matrix(truncate)
use bitmasks
implicit none
BEGIN_DOC
! Generate a wave function containing all possible products
! of alpha and beta determinants
END_DOC
logical, intent(in) :: truncate
integer :: i,j,k
integer(bit_kind) :: tmp_det(N_int,2)
integer :: idx
integer, external :: get_index_in_psi_det_sorted_bit
double precision :: norm(N_states)
call generate_all_alpha_beta_det_products
norm = 0.d0
do j=1,N_det_beta_unique
do k=1,N_int
tmp_det(k,2) = psi_det_beta_unique(k,j)
enddo
do i=1,N_det_alpha_unique
do k=1,N_int
tmp_det(k,1) = psi_det_alpha_unique(k,i)
enddo
idx = get_index_in_psi_det_sorted_bit(tmp_det,N_int)
if (idx > 0) then
do k=1,N_states
psi_coef_sorted_bit(idx,k) = psi_bilinear_matrix(i,j,k)
norm(k) += psi_bilinear_matrix(i,j,k)
enddo
endif
enddo
enddo
do k=1,N_states
norm(k) = 1.d0/dsqrt(norm(k))
do i=1,N_det
psi_coef_sorted_bit(i,k) = psi_coef_sorted_bit(i,k)*norm(k)
enddo
enddo
psi_det = psi_det_sorted_bit
psi_coef = psi_coef_sorted_bit
TOUCH psi_det psi_coef
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
norm(1) = 0.d0
do i=1,N_det
norm(1) += psi_average_norm_contrib_sorted(i)
if (truncate) then
if (norm(1) >= 0.999999d0) then
exit
endif
endif
enddo
N_det = min(i,N_det)
SOFT_TOUCH psi_det psi_coef N_det
end
subroutine generate_all_alpha_beta_det_products
implicit none
BEGIN_DOC
! Create a wave function from all possible alpha x beta determinants
END_DOC
integer :: i,j,k,l
integer :: idx, iproc
integer, external :: get_index_in_psi_det_sorted_bit
integer(bit_kind), allocatable :: tmp_det(:,:,:)
logical, external :: is_in_wavefunction
integer, external :: omp_get_thread_num
!$OMP PARALLEL DEFAULT(NONE) SHARED(psi_coef_sorted_bit,N_det_beta_unique,&
!$OMP N_det_alpha_unique, N_int, psi_det_alpha_unique, psi_det_beta_unique,&
!$OMP N_det) &
!$OMP PRIVATE(i,j,k,l,tmp_det,idx,iproc)
!$ iproc = omp_get_thread_num()
allocate (tmp_det(N_int,2,N_det_alpha_unique))
!$OMP DO
do j=1,N_det_beta_unique
l = 1
do i=1,N_det_alpha_unique
do k=1,N_int
tmp_det(k,1,l) = psi_det_alpha_unique(k,i)
tmp_det(k,2,l) = psi_det_beta_unique (k,j)
enddo
if (.not.is_in_wavefunction(tmp_det(1,1,l),N_int)) then
l = l+1
endif
enddo
call fill_H_apply_buffer_no_selection(l-1, tmp_det, N_int, iproc)
enddo
!$OMP END DO NOWAIT
deallocate(tmp_det)
!$OMP END PARALLEL
call copy_H_apply_buffer_to_wf
SOFT_TOUCH psi_det psi_coef N_det
end
Reference in New Issue View Git Blame Copy Permalink