mirror of
https://github.com/LCPQ/quantum_package
synced 2025-04-10 02:31:45 +02:00
Accelerated selection
This commit is contained in:
parent
b6ea2a8a45
commit
60164de0c0
@ -671,10 +671,10 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
|
||||
if(mat(1, p1, p2) == 0d0) cycle
|
||||
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
|
||||
logical, external :: is_in_wavefunction
|
||||
if (is_in_wavefunction(det,N_int)) then
|
||||
stop 'is_in_wf'
|
||||
cycle
|
||||
endif
|
||||
! if (is_in_wavefunction(det,N_int)) then
|
||||
! stop 'is_in_wf'
|
||||
! cycle
|
||||
! endif
|
||||
|
||||
if (do_ddci) then
|
||||
integer, external :: is_a_two_holes_two_particles
|
||||
|
||||
@ -41,31 +41,30 @@ subroutine sort_selection_buffer(b)
|
||||
implicit none
|
||||
|
||||
type(selection_buffer), intent(inout) :: b
|
||||
double precision, allocatable :: vals(:), absval(:)
|
||||
double precision, allocatable:: absval(:)
|
||||
integer, allocatable :: iorder(:)
|
||||
integer(bit_kind), allocatable :: detmp(:,:,:)
|
||||
double precision, pointer :: vals(:)
|
||||
integer(bit_kind), pointer :: detmp(:,:,:)
|
||||
integer :: i, nmwen
|
||||
logical, external :: detEq
|
||||
nmwen = min(b%N, b%cur)
|
||||
|
||||
|
||||
allocate(iorder(b%cur), detmp(N_int, 2, nmwen), absval(b%cur), vals(nmwen))
|
||||
allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)), absval(b%cur), vals(size(b%val)))
|
||||
absval = -dabs(b%val(:b%cur))
|
||||
do i=1,b%cur
|
||||
iorder(i) = i
|
||||
end do
|
||||
call dsort(absval, iorder, b%cur)
|
||||
|
||||
! Optimal for almost sorted data
|
||||
call insertion_dsort(absval, iorder, b%cur)
|
||||
do i=1, nmwen
|
||||
detmp(1:N_int,1,i) = b%det(1:N_int,1,iorder(i))
|
||||
detmp(1:N_int,2,i) = b%det(1:N_int,2,iorder(i))
|
||||
vals(i) = b%val(iorder(i))
|
||||
end do
|
||||
b%det = 0_bit_kind
|
||||
b%val = 0d0
|
||||
b%det(1:N_int,1,1:nmwen) = detmp(1:N_int,1,1:nmwen)
|
||||
b%det(1:N_int,2,1:nmwen) = detmp(1:N_int,2,1:nmwen)
|
||||
b%val(1:nmwen) = vals(1:nmwen)
|
||||
deallocate(b%det, b%val)
|
||||
b%det => detmp
|
||||
b%val => vals
|
||||
b%mini = max(b%mini,dabs(b%val(b%N)))
|
||||
b%cur = nmwen
|
||||
end subroutine
|
||||
|
||||
@ -1,9 +1,9 @@
|
||||
module selection_types
|
||||
type selection_buffer
|
||||
integer :: N, cur
|
||||
integer(8), allocatable :: det(:,:,:)
|
||||
double precision, allocatable :: val(:)
|
||||
double precision :: mini
|
||||
integer(8) , pointer :: det(:,:,:)
|
||||
double precision, pointer :: val(:)
|
||||
double precision :: mini
|
||||
endtype
|
||||
end module
|
||||
|
||||
|
||||
@ -18,7 +18,7 @@ interface: ezfio
|
||||
type: logical
|
||||
doc: Compute perturbative contribution of the Triples
|
||||
interface: ezfio,provider,ocaml
|
||||
default: true
|
||||
default: false
|
||||
|
||||
[energy]
|
||||
type: double precision
|
||||
|
||||
@ -2,11 +2,16 @@ use bitmasks
|
||||
|
||||
BEGIN_PROVIDER [ integer, N_int ]
|
||||
implicit none
|
||||
include 'Utils/constants.include.F'
|
||||
BEGIN_DOC
|
||||
! Number of 64-bit integers needed to represent determinants as binary strings
|
||||
END_DOC
|
||||
N_int = (mo_tot_num-1)/bit_kind_size + 1
|
||||
call write_int(6,N_int, 'N_int')
|
||||
call write_int(6,N_int, 'N_int')
|
||||
if (N_int > N_int_max) then
|
||||
stop 'N_int > N_int_max'
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
@ -444,7 +444,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
!$OMP DO SCHEDULE(static,4)
|
||||
do sh=1,shortcut(0,2)
|
||||
do i=shortcut(sh,2),shortcut(sh+1,2)-1
|
||||
org_i = sort_idx(i,2)
|
||||
@ -477,7 +477,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
!$OMP DO SCHEDULE(static,4)
|
||||
do sh=1,shortcut(0,1)
|
||||
do sh2=1,shortcut(0,1)
|
||||
if (sh==sh2) cycle
|
||||
|
||||
@ -32,102 +32,105 @@ END_PROVIDER
|
||||
double precision, allocatable :: tmp_a(:,:,:), tmp_b(:,:,:)
|
||||
integer :: krow, kcol, lrow, lcol
|
||||
|
||||
one_body_dm_mo_alpha = 0.d0
|
||||
one_body_dm_mo_beta = 0.d0
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
|
||||
!$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
|
||||
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
|
||||
!$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
|
||||
!$OMP mo_tot_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns, &
|
||||
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, &
|
||||
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values)
|
||||
allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
|
||||
tmp_a = 0.d0
|
||||
tmp_b = 0.d0
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
do k=1,N_det
|
||||
krow = psi_bilinear_matrix_rows(k)
|
||||
kcol = psi_bilinear_matrix_columns(k)
|
||||
tmp_det(:,1) = psi_det(:,1, krow)
|
||||
tmp_det(:,2) = psi_det(:,2, kcol)
|
||||
call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
|
||||
do m=1,N_states
|
||||
ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m)
|
||||
do l=1,elec_alpha_num
|
||||
j = occ(l,1)
|
||||
tmp_a(j,j,m) += ck
|
||||
enddo
|
||||
do l=1,elec_beta_num
|
||||
j = occ(l,2)
|
||||
tmp_b(j,j,m) += ck
|
||||
enddo
|
||||
enddo
|
||||
PROVIDE psi_det
|
||||
|
||||
l = k+1
|
||||
lrow = psi_bilinear_matrix_rows(l)
|
||||
lcol = psi_bilinear_matrix_columns(l)
|
||||
do while ( lcol == kcol )
|
||||
tmp_det2(:,1) = psi_det(:,1, lrow)
|
||||
tmp_det2(:,2) = psi_det(:,2, lcol)
|
||||
call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
|
||||
if (degree == 1) then
|
||||
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
do m=1,N_states
|
||||
ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase
|
||||
if (s1==1) then
|
||||
tmp_a(h1,p1,m) += ckl
|
||||
tmp_a(p1,h1,m) += ckl
|
||||
else
|
||||
tmp_b(h1,p1,m) += ckl
|
||||
tmp_b(p1,h1,m) += ckl
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
l = l+1
|
||||
if (l>N_det) exit
|
||||
lrow = psi_bilinear_matrix_rows(l)
|
||||
lcol = psi_bilinear_matrix_columns(l)
|
||||
enddo
|
||||
one_body_dm_mo_alpha = 0.d0
|
||||
one_body_dm_mo_beta = 0.d0
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
|
||||
!$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
|
||||
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,&
|
||||
!$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
|
||||
!$OMP mo_tot_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns, &
|
||||
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, &
|
||||
!$OMP psi_bilinear_matrix_transp_order, psi_det_alpha_unique, psi_det_beta_unique, &
|
||||
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values)
|
||||
allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) )
|
||||
tmp_a = 0.d0
|
||||
tmp_b = 0.d0
|
||||
!$OMP DO SCHEDULE(guided)
|
||||
do k=1,N_det
|
||||
krow = psi_bilinear_matrix_rows(k)
|
||||
kcol = psi_bilinear_matrix_columns(k)
|
||||
tmp_det(:,1) = psi_det_alpha_unique(:,krow)
|
||||
tmp_det(:,2) = psi_det_beta_unique (:,kcol)
|
||||
call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int)
|
||||
do m=1,N_states
|
||||
ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m)
|
||||
do l=1,elec_alpha_num
|
||||
j = occ(l,1)
|
||||
tmp_a(j,j,m) += ck
|
||||
enddo
|
||||
do l=1,elec_beta_num
|
||||
j = occ(l,2)
|
||||
tmp_b(j,j,m) += ck
|
||||
enddo
|
||||
enddo
|
||||
|
||||
l = k+1
|
||||
lrow = psi_bilinear_matrix_transp_rows(l)
|
||||
lcol = psi_bilinear_matrix_transp_columns(l)
|
||||
do while ( lrow == krow )
|
||||
tmp_det2(:,1) = psi_det(:,1, lrow)
|
||||
tmp_det2(:,2) = psi_det(:,2, lcol)
|
||||
call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
|
||||
if (degree == 1) then
|
||||
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
do m=1,N_states
|
||||
ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_transp_values(l,m) * phase
|
||||
if (s1==1) then
|
||||
tmp_a(h1,p1,m) += ckl
|
||||
tmp_a(p1,h1,m) += ckl
|
||||
else
|
||||
tmp_b(h1,p1,m) += ckl
|
||||
tmp_b(p1,h1,m) += ckl
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
l = l+1
|
||||
if (l>N_det) exit
|
||||
lrow = psi_bilinear_matrix_transp_rows(l)
|
||||
lcol = psi_bilinear_matrix_transp_columns(l)
|
||||
enddo
|
||||
l = k+1
|
||||
lrow = psi_bilinear_matrix_rows(l)
|
||||
lcol = psi_bilinear_matrix_columns(l)
|
||||
do while ( lcol == kcol )
|
||||
tmp_det2(:,1) = psi_det_alpha_unique(:, lrow)
|
||||
tmp_det2(:,2) = psi_det_beta_unique (:, lcol)
|
||||
call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
|
||||
if (degree == 1) then
|
||||
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
do m=1,N_states
|
||||
ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase
|
||||
if (s1==1) then
|
||||
tmp_a(h1,p1,m) += ckl
|
||||
tmp_a(p1,h1,m) += ckl
|
||||
else
|
||||
tmp_b(h1,p1,m) += ckl
|
||||
tmp_b(p1,h1,m) += ckl
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
l = l+1
|
||||
if (l>N_det) exit
|
||||
lrow = psi_bilinear_matrix_rows(l)
|
||||
lcol = psi_bilinear_matrix_columns(l)
|
||||
enddo
|
||||
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
!$OMP CRITICAL
|
||||
one_body_dm_mo_alpha(:,:,:) = one_body_dm_mo_alpha(:,:,:) + tmp_a(:,:,:)
|
||||
!$OMP END CRITICAL
|
||||
!$OMP CRITICAL
|
||||
one_body_dm_mo_beta(:,:,:) = one_body_dm_mo_beta(:,:,:) + tmp_b(:,:,:)
|
||||
!$OMP END CRITICAL
|
||||
deallocate(tmp_a,tmp_b)
|
||||
!$OMP END PARALLEL
|
||||
l = psi_bilinear_matrix_transp_order(k)+1
|
||||
lrow = psi_bilinear_matrix_transp_rows(l)
|
||||
lcol = psi_bilinear_matrix_transp_columns(l)
|
||||
do while ( lrow == krow )
|
||||
tmp_det2(:,1) = psi_det_alpha_unique(:, lrow)
|
||||
tmp_det2(:,2) = psi_det_beta_unique (:, lcol)
|
||||
call get_excitation_degree(tmp_det,tmp_det2,degree,N_int)
|
||||
if (degree == 1) then
|
||||
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
do m=1,N_states
|
||||
ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_transp_values(l,m) * phase
|
||||
if (s1==1) then
|
||||
tmp_a(h1,p1,m) += ckl
|
||||
tmp_a(p1,h1,m) += ckl
|
||||
else
|
||||
tmp_b(h1,p1,m) += ckl
|
||||
tmp_b(p1,h1,m) += ckl
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
l = l+1
|
||||
if (l>N_det) exit
|
||||
lrow = psi_bilinear_matrix_transp_rows(l)
|
||||
lcol = psi_bilinear_matrix_transp_columns(l)
|
||||
enddo
|
||||
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
!$OMP CRITICAL
|
||||
one_body_dm_mo_alpha(:,:,:) = one_body_dm_mo_alpha(:,:,:) + tmp_a(:,:,:)
|
||||
!$OMP END CRITICAL
|
||||
!$OMP CRITICAL
|
||||
one_body_dm_mo_beta(:,:,:) = one_body_dm_mo_beta(:,:,:) + tmp_b(:,:,:)
|
||||
!$OMP END CRITICAL
|
||||
deallocate(tmp_a,tmp_b)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
@ -1,32 +1,50 @@
|
||||
subroutine get_excitation_degree(key1,key2,degree,Nint)
|
||||
use bitmasks
|
||||
include 'Utils/constants.include.F'
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the excitation degree between two determinants
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key1(Nint,2)
|
||||
integer(bit_kind), intent(in) :: key2(Nint,2)
|
||||
integer(bit_kind), intent(in) :: key1(Nint*2)
|
||||
integer(bit_kind), intent(in) :: key2(Nint*2)
|
||||
integer, intent(out) :: degree
|
||||
|
||||
integer(bit_kind) :: xorvec(2*N_int_max)
|
||||
integer :: l
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
|
||||
degree = popcnt(xor( key1(1,1), key2(1,1))) + &
|
||||
popcnt(xor( key1(1,2), key2(1,2)))
|
||||
!DIR$ NOUNROLL
|
||||
do l=2,Nint
|
||||
degree = degree+ popcnt(xor( key1(l,1), key2(l,1))) + &
|
||||
popcnt(xor( key1(l,2), key2(l,2)))
|
||||
enddo
|
||||
ASSERT (degree >= 0)
|
||||
select case (Nint)
|
||||
|
||||
case (1)
|
||||
xorvec(1:2) = xor( key1(1:2), key2(1:2))
|
||||
degree = sum(popcnt(xorvec(1:2)))
|
||||
|
||||
case (2)
|
||||
xorvec(1:4) = xor( key1(1:4), key2(1:4))
|
||||
degree = sum(popcnt(xorvec(1:4)))
|
||||
|
||||
case (3)
|
||||
xorvec(1:6) = xor( key1(1:6), key2(1:6))
|
||||
degree = sum(popcnt(xorvec(1:6)))
|
||||
|
||||
case (4)
|
||||
xorvec(1:8) = xor( key1(1:8), key2(1:8))
|
||||
degree = sum(popcnt(xorvec(1:8)))
|
||||
|
||||
case default
|
||||
l = ishft(Nint,1)
|
||||
xorvec(1:l) = xor( key1(1:l), key2(1:l))
|
||||
degree = sum(popcnt(xorvec(1:l)))
|
||||
|
||||
end select
|
||||
|
||||
degree = ishft(degree,-1)
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine get_excitation(det1,det2,exc,degree,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
@ -2206,3 +2224,335 @@ subroutine u_0_H_u_0_stored(e_0,u_0,hmatrix,sze)
|
||||
call matrix_vector_product(u_0,v_0,hmatrix,sze,sze)
|
||||
e_0 = u_dot_v(v_0,u_0,sze)
|
||||
end
|
||||
|
||||
|
||||
|
||||
! Spin-determinant routines
|
||||
! -------------------------
|
||||
|
||||
subroutine get_excitation_degree_spin(key1,key2,degree,Nint)
|
||||
use bitmasks
|
||||
include 'Utils/constants.include.F'
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the excitation degree between two determinants
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key1(Nint)
|
||||
integer(bit_kind), intent(in) :: key2(Nint)
|
||||
integer, intent(out) :: degree
|
||||
|
||||
integer(bit_kind) :: xorvec(N_int_max)
|
||||
integer :: l
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
|
||||
select case (Nint)
|
||||
|
||||
case (1)
|
||||
xorvec(1) = xor( key1(1), key2(1))
|
||||
degree = popcnt(xorvec(1))
|
||||
|
||||
case (2)
|
||||
xorvec(1:2) = xor( key1(1:2), key2(1:2))
|
||||
degree = sum(popcnt(xorvec(1:2)))
|
||||
|
||||
case (3)
|
||||
xorvec(1:3) = xor( key1(1:3), key2(1:3))
|
||||
degree = sum(popcnt(xorvec(1:3)))
|
||||
|
||||
case (4)
|
||||
xorvec(1:4) = xor( key1(1:4), key2(1:4))
|
||||
degree = sum(popcnt(xorvec(1:4)))
|
||||
|
||||
case default
|
||||
xorvec(1:Nint) = xor( key1(1:Nint), key2(1:Nint))
|
||||
degree = sum(popcnt(xorvec(1:Nint)))
|
||||
|
||||
end select
|
||||
|
||||
degree = ishft(degree,-1)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine get_excitation_spin(det1,det2,exc,degree,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the excitation operators between two determinants and the phase
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: det1(Nint)
|
||||
integer(bit_kind), intent(in) :: det2(Nint)
|
||||
integer, intent(out) :: exc(0:2,2)
|
||||
integer, intent(out) :: degree
|
||||
double precision, intent(out) :: phase
|
||||
! exc(number,hole/particle)
|
||||
! ex :
|
||||
! exc(0,1) = number of holes
|
||||
! exc(0,2) = number of particles
|
||||
! exc(1,2) = first particle
|
||||
! exc(1,1) = first hole
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call get_excitation_degree_spin(det1,det2,degree,Nint)
|
||||
select case (degree)
|
||||
|
||||
case (3:)
|
||||
degree = -1
|
||||
return
|
||||
|
||||
case (2)
|
||||
call get_double_excitation_spin(det1,det2,exc,phase,Nint)
|
||||
return
|
||||
|
||||
case (1)
|
||||
call get_mono_excitation_spin(det1,det2,exc,phase,Nint)
|
||||
return
|
||||
|
||||
case(0)
|
||||
return
|
||||
|
||||
end select
|
||||
end
|
||||
|
||||
subroutine decode_exc_spin(exc,h1,p1,h2,p2)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Decodes the exc arrays returned by get_excitation.
|
||||
! h1,h2 : Holes
|
||||
! p1,p2 : Particles
|
||||
END_DOC
|
||||
integer, intent(in) :: exc(0:2,2)
|
||||
integer, intent(out) :: h1,h2,p1,p2
|
||||
|
||||
select case (exc(0,1))
|
||||
case(2)
|
||||
h1 = exc(1,1)
|
||||
h2 = exc(2,1)
|
||||
p1 = exc(1,2)
|
||||
p2 = exc(2,2)
|
||||
case(1)
|
||||
h1 = exc(1,1)
|
||||
h2 = 0
|
||||
p1 = exc(1,2)
|
||||
p2 = 0
|
||||
case default
|
||||
h1 = 0
|
||||
p1 = 0
|
||||
h2 = 0
|
||||
p2 = 0
|
||||
end select
|
||||
end
|
||||
|
||||
|
||||
subroutine get_double_excitation_spin(det1,det2,exc,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the two excitation operators between two doubly excited spin-determinants
|
||||
! and the phase
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: det1(Nint)
|
||||
integer(bit_kind), intent(in) :: det2(Nint)
|
||||
integer, intent(out) :: exc(0:2,2)
|
||||
double precision, intent(out) :: phase
|
||||
integer :: tz
|
||||
integer :: l, idx_hole, idx_particle, ishift
|
||||
integer :: nperm
|
||||
integer :: i,j,k,m,n
|
||||
integer :: high, low
|
||||
integer :: a,b,c,d
|
||||
integer(bit_kind) :: hole, particle, tmp
|
||||
double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
nperm = 0
|
||||
exc(0,1) = 0
|
||||
exc(0,2) = 0
|
||||
|
||||
idx_particle = 0
|
||||
idx_hole = 0
|
||||
ishift = 1-bit_kind_size
|
||||
do l=1,Nint
|
||||
ishift = ishift + bit_kind_size
|
||||
if (det1(l) == det2(l)) then
|
||||
cycle
|
||||
endif
|
||||
tmp = xor( det1(l), det2(l) )
|
||||
particle = iand(tmp, det2(l))
|
||||
hole = iand(tmp, det1(l))
|
||||
do while (particle /= 0_bit_kind)
|
||||
tz = trailz(particle)
|
||||
idx_particle = idx_particle + 1
|
||||
exc(0,2) = exc(0,2) + 1
|
||||
exc(idx_particle,2) = tz+ishift
|
||||
particle = iand(particle,particle-1_bit_kind)
|
||||
enddo
|
||||
if (iand(exc(0,1),exc(0,2))==2) then ! exc(0,1)==2 or exc(0,2)==2
|
||||
exit
|
||||
endif
|
||||
do while (hole /= 0_bit_kind)
|
||||
tz = trailz(hole)
|
||||
idx_hole = idx_hole + 1
|
||||
exc(0,1) = exc(0,1) + 1
|
||||
exc(idx_hole,1) = tz+ishift
|
||||
hole = iand(hole,hole-1_bit_kind)
|
||||
enddo
|
||||
if (iand(exc(0,1),exc(0,2))==2) then ! exc(0,1)==2 or exc(0,2)==2
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
|
||||
select case (exc(0,1))
|
||||
|
||||
case(1)
|
||||
low = min(exc(1,1), exc(1,2))
|
||||
high = max(exc(1,1), exc(1,2))
|
||||
|
||||
ASSERT (low > 0)
|
||||
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
|
||||
n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size)
|
||||
ASSERT (high > 0)
|
||||
k = ishft(high-1,-bit_kind_shift)+1
|
||||
m = iand(high-1,bit_kind_size-1)+1
|
||||
|
||||
if (j==k) then
|
||||
nperm = nperm + popcnt(iand(det1(j), &
|
||||
iand( ibset(0_bit_kind,m-1)-1_bit_kind, &
|
||||
ibclr(-1_bit_kind,n)+1_bit_kind ) ))
|
||||
else
|
||||
nperm = nperm + popcnt(iand(det1(k), &
|
||||
ibset(0_bit_kind,m-1)-1_bit_kind))
|
||||
if (n < bit_kind_size) then
|
||||
nperm = nperm + popcnt(iand(det1(j), ibclr(-1_bit_kind,n) +1_bit_kind))
|
||||
endif
|
||||
do i=j+1,k-1
|
||||
nperm = nperm + popcnt(det1(i))
|
||||
end do
|
||||
endif
|
||||
|
||||
case (2)
|
||||
|
||||
do i=1,2
|
||||
low = min(exc(i,1), exc(i,2))
|
||||
high = max(exc(i,1), exc(i,2))
|
||||
|
||||
ASSERT (low > 0)
|
||||
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
|
||||
n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size)
|
||||
ASSERT (high > 0)
|
||||
k = ishft(high-1,-bit_kind_shift)+1
|
||||
m = iand(high-1,bit_kind_size-1)+1
|
||||
|
||||
if (j==k) then
|
||||
nperm = nperm + popcnt(iand(det1(j), &
|
||||
iand( ibset(0_bit_kind,m-1)-1_bit_kind, &
|
||||
ibclr(-1_bit_kind,n)+1_bit_kind ) ))
|
||||
else
|
||||
nperm = nperm + popcnt(iand(det1(k), &
|
||||
ibset(0_bit_kind,m-1)-1_bit_kind))
|
||||
if (n < bit_kind_size) then
|
||||
nperm = nperm + popcnt(iand(det1(j), ibclr(-1_bit_kind,n) +1_bit_kind))
|
||||
endif
|
||||
do l=j+1,k-1
|
||||
nperm = nperm + popcnt(det1(l))
|
||||
end do
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
a = min(exc(1,1), exc(1,2))
|
||||
b = max(exc(1,1), exc(1,2))
|
||||
c = min(exc(2,1), exc(2,2))
|
||||
d = max(exc(2,1), exc(2,2))
|
||||
if (c>a .and. c<b .and. d>b) then
|
||||
nperm = nperm + 1
|
||||
endif
|
||||
end select
|
||||
|
||||
phase = phase_dble(iand(nperm,1))
|
||||
|
||||
end
|
||||
|
||||
subroutine get_mono_excitation_spin(det1,det2,exc,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the excitation operator between two singly excited determinants and the phase
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: det1(Nint)
|
||||
integer(bit_kind), intent(in) :: det2(Nint)
|
||||
integer, intent(out) :: exc(0:2,2)
|
||||
double precision, intent(out) :: phase
|
||||
integer :: tz
|
||||
integer :: l, idx_hole, idx_particle, ishift
|
||||
integer :: nperm
|
||||
integer :: i,j,k,m,n
|
||||
integer :: high, low
|
||||
integer :: a,b,c,d
|
||||
integer(bit_kind) :: hole, particle, tmp
|
||||
double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
nperm = 0
|
||||
exc(0,1) = 0
|
||||
exc(0,2) = 0
|
||||
|
||||
ishift = 1-bit_kind_size
|
||||
do l=1,Nint
|
||||
ishift = ishift + bit_kind_size
|
||||
if (det1(l) == det2(l)) then
|
||||
cycle
|
||||
endif
|
||||
tmp = xor( det1(l), det2(l) )
|
||||
particle = iand(tmp, det2(l))
|
||||
hole = iand(tmp, det1(l))
|
||||
if (particle /= 0_bit_kind) then
|
||||
tz = trailz(particle)
|
||||
exc(0,2) = 1
|
||||
exc(1,2) = tz+ishift
|
||||
endif
|
||||
if (hole /= 0_bit_kind) then
|
||||
tz = trailz(hole)
|
||||
exc(0,1) = 1
|
||||
exc(1,1) = tz+ishift
|
||||
endif
|
||||
|
||||
if ( iand(exc(0,1),exc(0,2)) /= 1) then ! exc(0,1)/=1 and exc(0,2) /= 1
|
||||
cycle
|
||||
endif
|
||||
|
||||
low = min(exc(1,1),exc(1,2))
|
||||
high = max(exc(1,1),exc(1,2))
|
||||
|
||||
ASSERT (low > 0)
|
||||
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
|
||||
n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size)
|
||||
ASSERT (high > 0)
|
||||
k = ishft(high-1,-bit_kind_shift)+1
|
||||
m = iand(high-1,bit_kind_size-1)+1
|
||||
if (j==k) then
|
||||
nperm = popcnt(iand(det1(j), &
|
||||
iand(ibset(0_bit_kind,m-1)-1_bit_kind,ibclr(-1_bit_kind,n)+1_bit_kind)))
|
||||
else
|
||||
nperm = nperm + popcnt(iand(det1(k),ibset(0_bit_kind,m-1)-1_bit_kind))
|
||||
if (n < bit_kind_size) then
|
||||
nperm = nperm + popcnt(iand(det1(j),ibclr(-1_bit_kind,n)+1_bit_kind))
|
||||
endif
|
||||
do i=j+1,k-1
|
||||
nperm = nperm + popcnt(det1(i))
|
||||
end do
|
||||
endif
|
||||
phase = phase_dble(iand(nperm,1))
|
||||
return
|
||||
|
||||
enddo
|
||||
end
|
||||
|
||||
|
||||
@ -386,8 +386,9 @@ END_PROVIDER
|
||||
!==============================================================================!
|
||||
|
||||
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_rows , (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns, (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order , (N_det) ]
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -395,6 +396,8 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
|
||||
! D_a^t C D_b
|
||||
!
|
||||
! Rows are alpha determinants and columns are beta.
|
||||
!
|
||||
! Order refers to psi_det
|
||||
END_DOC
|
||||
integer :: i,j,k, l
|
||||
integer(bit_kind) :: tmp_det(N_int,2)
|
||||
@ -404,10 +407,10 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
|
||||
|
||||
PROVIDE psi_coef_sorted_bit
|
||||
|
||||
integer, allocatable :: iorder(:), to_sort(:)
|
||||
integer, allocatable :: 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))
|
||||
allocate(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)
|
||||
@ -418,36 +421,40 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
|
||||
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
|
||||
psi_bilinear_matrix_order(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 isort(to_sort, psi_bilinear_matrix_order, N_det)
|
||||
call iset_order(psi_bilinear_matrix_rows,psi_bilinear_matrix_order,N_det)
|
||||
call iset_order(psi_bilinear_matrix_columns,psi_bilinear_matrix_order,N_det)
|
||||
do l=1,N_states
|
||||
call dset_order(psi_bilinear_matrix_values(1,l),iorder,N_det)
|
||||
call dset_order(psi_bilinear_matrix_values(1,l),psi_bilinear_matrix_order,N_det)
|
||||
enddo
|
||||
deallocate(iorder,to_sort)
|
||||
deallocate(to_sort)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows, (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows , (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_columns, (N_det) ]
|
||||
&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_order , (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
|
||||
!
|
||||
! Rows are Beta determinants and columns are alpha
|
||||
! Rows are Alpha determinants and columns are beta, but the matrix is stored in row major
|
||||
! format
|
||||
!
|
||||
! Order refers to psi_bilinear_matrix
|
||||
END_DOC
|
||||
integer :: i,j,k,l
|
||||
|
||||
|
||||
PROVIDE psi_coef_sorted_bit
|
||||
|
||||
integer, allocatable :: iorder(:), to_sort(:)
|
||||
allocate(iorder(N_det), to_sort(N_det))
|
||||
integer, allocatable :: to_sort(:)
|
||||
allocate(to_sort(N_det))
|
||||
do l=1,N_states
|
||||
do k=1,N_det
|
||||
psi_bilinear_matrix_transp_values (k,l) = psi_bilinear_matrix_values (k,l)
|
||||
@ -459,15 +466,15 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_
|
||||
i = psi_bilinear_matrix_transp_columns(k)
|
||||
j = psi_bilinear_matrix_transp_rows (k)
|
||||
to_sort(k) = N_det_beta_unique * (j-1) + i
|
||||
iorder(k) = k
|
||||
psi_bilinear_matrix_transp_order(k) = k
|
||||
enddo
|
||||
call isort(to_sort, iorder, N_det)
|
||||
call iset_order(psi_bilinear_matrix_transp_rows,iorder,N_det)
|
||||
call iset_order(psi_bilinear_matrix_transp_columns,iorder,N_det)
|
||||
call isort(to_sort, psi_bilinear_matrix_transp_order, N_det)
|
||||
call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det)
|
||||
call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det)
|
||||
do l=1,N_states
|
||||
call dset_order(psi_bilinear_matrix_transp_values(1,l),iorder,N_det)
|
||||
call dset_order(psi_bilinear_matrix_transp_values(1,l),psi_bilinear_matrix_transp_order,N_det)
|
||||
enddo
|
||||
deallocate(iorder,to_sort)
|
||||
deallocate(to_sort)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
@ -1,5 +1,6 @@
|
||||
integer, parameter :: max_dim = 511
|
||||
integer, parameter :: SIMD_vector = 32
|
||||
integer, parameter :: N_int_max = 16
|
||||
|
||||
double precision, parameter :: pi = dacos(-1.d0)
|
||||
double precision, parameter :: sqpi = dsqrt(dacos(-1.d0))
|
||||
|
||||
@ -13,7 +13,7 @@ source $QP_ROOT/tests/bats/common.bats.sh
|
||||
qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]"
|
||||
qp_run cassd_zmq $INPUT
|
||||
energy="$(ezfio get cas_sd_zmq energy_pt2)"
|
||||
eq $energy -76.231084536315 5.E-5
|
||||
eq $energy -76.231248286858 5.E-5
|
||||
|
||||
ezfio set determinants n_det_max 1024
|
||||
ezfio set determinants read_wf True
|
||||
@ -21,6 +21,6 @@ source $QP_ROOT/tests/bats/common.bats.sh
|
||||
qp_run cassd_zmq $INPUT
|
||||
ezfio set determinants read_wf False
|
||||
energy="$(ezfio get cas_sd_zmq energy)"
|
||||
eq $energy -76.2225863580749 2.E-5
|
||||
eq $energy -76.2225678834779 2.E-5
|
||||
}
|
||||
|
||||
|
||||
@ -42,11 +42,13 @@ function run_FCI_ZMQ() {
|
||||
qp_set_mo_class h2o.ezfio -core "[1]" -act "[2-12]" -del "[13-24]"
|
||||
}
|
||||
@test "FCI H2O cc-pVDZ" {
|
||||
run_FCI h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02
|
||||
run_FCI h2o.ezfio 2000 -76.1253758241716 -76.1258130146102
|
||||
}
|
||||
|
||||
|
||||
|
||||
@test "FCI-ZMQ H2O cc-pVDZ" {
|
||||
run_FCI_ZMQ h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02
|
||||
run_FCI_ZMQ h2o.ezfio 2000 -76.1250552686394 -76.1258817228809
|
||||
}
|
||||
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user