diff --git a/plugins/CAS_SD_ZMQ/selection.irp.f b/plugins/CAS_SD_ZMQ/selection.irp.f index 76ba04d6..bdce5df5 100644 --- a/plugins/CAS_SD_ZMQ/selection.irp.f +++ b/plugins/CAS_SD_ZMQ/selection.irp.f @@ -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 diff --git a/plugins/Full_CI_ZMQ/selection_buffer.irp.f b/plugins/Full_CI_ZMQ/selection_buffer.irp.f index d296b399..4ca3830e 100644 --- a/plugins/Full_CI_ZMQ/selection_buffer.irp.f +++ b/plugins/Full_CI_ZMQ/selection_buffer.irp.f @@ -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 diff --git a/plugins/Full_CI_ZMQ/selection_types.f90 b/plugins/Full_CI_ZMQ/selection_types.f90 index 9506629c..29e48524 100644 --- a/plugins/Full_CI_ZMQ/selection_types.f90 +++ b/plugins/Full_CI_ZMQ/selection_types.f90 @@ -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 diff --git a/plugins/mrcepa0/EZFIO.cfg b/plugins/mrcepa0/EZFIO.cfg index a2dc1bb3..53519ec7 100644 --- a/plugins/mrcepa0/EZFIO.cfg +++ b/plugins/mrcepa0/EZFIO.cfg @@ -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 diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 10ab6f67..e50cf25a 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -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 diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index 8ef574d4..173ceb57 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -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 diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index e0ec6a6b..dbd934d7 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -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 diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index eebf9611..45458544 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -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. cb) 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 + diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index acc49d50..d013513e 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -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 diff --git a/src/Utils/constants.include.F b/src/Utils/constants.include.F index 991ef80a..4974fd8e 100644 --- a/src/Utils/constants.include.F +++ b/src/Utils/constants.include.F @@ -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)) diff --git a/tests/bats/cassd.bats b/tests/bats/cassd.bats index 1b845e91..67c35235 100644 --- a/tests/bats/cassd.bats +++ b/tests/bats/cassd.bats @@ -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 } diff --git a/tests/bats/fci.bats b/tests/bats/fci.bats index 79ff91ab..6cded581 100644 --- a/tests/bats/fci.bats +++ b/tests/bats/fci.bats @@ -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 }