diff --git a/plugins/MRPT_Utils/new_way_second_order_coef.irp.f b/plugins/MRPT_Utils/new_way_second_order_coef.irp.f index 4c12dbe1..9b25da95 100644 --- a/plugins/MRPT_Utils/new_way_second_order_coef.irp.f +++ b/plugins/MRPT_Utils/new_way_second_order_coef.irp.f @@ -210,10 +210,6 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p) ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb} hab = (fock_operator_local(aorb,borb,kspin) ) * phase - if(isnan(hab))then - print*, '1' - stop - endif ! < jdet | H | det_tmp_bis > = phase * (ir|cv) call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int) if(ispin == jspin)then diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index 173ceb57..59d69453 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -609,3 +609,332 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8) deallocate (shortcut, sort_idx, sorted, version, ut) end + + + + +subroutine H_S2_u_0_nstates_new(v_0,s_0,H_jj,S2_jj,N_st,sze_8) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + ! + ! n : number of determinants + ! + ! H_jj : array of + ! + ! S2_jj : array of + END_DOC + integer, intent(in) :: N_st,sze_8 + double precision, intent(out) :: v_0(sze_8,N_st), s_0(sze_8,N_st) + double precision, intent(in) :: H_jj(*), S2_jj(*) + + + PROVIDE ref_bitmask_energy + + double precision :: hij, s2 + integer :: i,j + integer :: k_a, k_b, l_a, l_b, m_a, m_b + integer :: degree, istate + integer :: krow, kcol, krow_b, kcol_b + integer :: lrow, lcol + integer :: mrow, mcol + integer(bit_kind) :: spindet(N_int) + integer(bit_kind) :: tmp_det(N_int,2) + integer(bit_kind) :: tmp_det2(N_int,2) + integer(bit_kind) :: tmp_det3(N_int,2) + integer(bit_kind), allocatable :: buffer(:,:) + double precision :: ck(N_st), cl(N_st), cm(N_st) + integer :: n_singles, n_doubles + integer, allocatable :: singles(:), doubles(:) + integer, allocatable :: idx(:), idx0(:) + logical, allocatable :: is_single_a(:) + + allocate( buffer(N_int,N_det_alpha_unique), & + singles(N_det_alpha_unique), doubles(N_det_alpha_unique), & + is_single_a(N_det_alpha_unique), & + idx(N_det_alpha_unique), idx0(N_det_alpha_unique) ) + + do k_a=1,N_det-1 + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + kcol = psi_bilinear_matrix_columns(k_a) + + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol) + + ! Initial determinant is at k_b in beta-major representation + ! ---------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_reverse(k_a) + + ! Get all single and double alpha excitations + ! =========================================== + + l_a = k_a+1 + spindet(1:N_int) = tmp_det(1:N_int,1) + + ! Loop inside the beta column to gather all the connected alphas + i=1 + lcol = kcol + do while ( (lcol == kcol).and.(l_a <= N_det) ) + lrow = psi_bilinear_matrix_rows(l_a) + lcol = psi_bilinear_matrix_columns(l_a) + buffer(1:N_int,i) = psi_det_alpha_unique(1:N_int, lrow) + idx(i) = lrow + i=i+1 + l_a = l_a + 1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles( & + buffer, idx, spindet, N_int, i, & + singles, doubles, n_singles, n_doubles ) + + ! Compute Hij for all alpha singles + ! ---------------------------------- + + l_a = k_a+1 + lrow = psi_bilinear_matrix_rows(l_a) + do i=1,n_singles + call i_H_j_mono_spin( tmp_det(1,1), psi_det_alpha_unique(1, singles(i)), N_int, 1, hij) + do while ( lrow < singles(i) ) + l_a = l_a+1 + lrow = psi_bilinear_matrix_rows(l_a) + enddo + v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st) + v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st) + + enddo + + ! Compute Hij for all alpha doubles + ! ---------------------------------- + + l_a = k_a+1 + lrow = psi_bilinear_matrix_rows(l_a) + do i=1,n_doubles + call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, doubles(i)), N_int, hij) + do while (lrow < doubles(i)) + l_a = l_a+1 + lrow = psi_bilinear_matrix_rows(l_a) + enddo + v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st) + v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st) + enddo + + + + ! Get all single and double beta excitations + ! =========================================== + + l_b = k_b+1 + spindet(1:N_int) = tmp_det(1:N_int,2) + + ! Loop inside the alpha row to gather all the connected betas + i=1 + lrow=krow + do while ( (lrow == krow).and.(l_b <= N_det) ) + lrow = psi_bilinear_matrix_transp_rows(l_b) + lcol = psi_bilinear_matrix_transp_columns(l_b) + buffer(1:N_int,i) = psi_det_beta_unique(1:N_int, lcol) + idx(i) = lcol + i=i+1 + l_b = l_b + 1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles( & + buffer, idx, spindet, N_int, i, & + singles, doubles, n_singles, n_doubles ) + + ! Compute Hij for all beta singles + ! ---------------------------------- + + l_b = k_b + lcol = psi_bilinear_matrix_transp_columns(l_b) + do i=1,n_singles + call i_H_j_mono_spin( tmp_det(1,2), psi_det_beta_unique(1, singles(i)), N_int, 2, hij) + do while ( lcol < singles(i) ) + l_b = l_b+1 + lcol = psi_bilinear_matrix_transp_columns(l_b) + enddo + l_a = psi_bilinear_matrix_order_reverse(l_b) + v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st) + v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st) + enddo + + ! Compute Hij for all beta doubles + ! ---------------------------------- + + l_b = k_b + lcol = psi_bilinear_matrix_transp_columns(l_b) + do i=1,n_doubles + call i_H_j_double_spin( tmp_det(1,2), psi_det_beta_unique(1, doubles(i)), N_int, hij) + do while (lcol < doubles(i)) + l_b = l_b+1 + lcol = psi_bilinear_matrix_transp_columns(l_b) + enddo + l_a = psi_bilinear_matrix_order_reverse(l_b) + v_0(l_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st) + v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(l_a,1:N_st) + enddo + + end do + + ! Alpha/Beta double excitations + ! ============================= + + do i=1,N_det_beta_unique + idx0(i) = i + enddo + is_single_a(:) = .False. + + k_a=1 + do i=1,N_det_beta_unique + + ! Select a beta determinant + ! ------------------------- + + spindet(1:N_int) = psi_det_beta_unique(1:N_int, i) + tmp_det(1:N_int,2) = spindet(1:N_int) + + call get_all_spin_singles( & + psi_det_beta_unique, idx0, spindet, N_int, N_det_beta_unique, & + singles, n_singles ) + + do j=1,n_singles + is_single_a( singles(j) ) = .True. + enddo + + ! For all alpha.beta pairs with the selected beta + ! ----------------------------------------------- + + kcol = psi_bilinear_matrix_columns(k_a) + do while (kcol < i) + k_a = k_a+1 + if (k_a > N_det) exit + kcol = psi_bilinear_matrix_columns(k_a) + enddo + + do while (kcol == i) + + krow = psi_bilinear_matrix_rows(k_a) + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) + + ! Loop over all alpha.beta pairs with a single exc alpha + ! ------------------------------------------------------ + + l_a = k_a+1 + if (l_a > N_det) exit + lrow = psi_bilinear_matrix_rows(l_a) + lcol = psi_bilinear_matrix_columns(l_a) + + do while (lrow == krow) + + ! Loop over all alpha.beta pairs with a single exc alpha + ! ------------------------------------------------------ + if (is_single_a(lrow)) then + + tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int,lrow) + + ! Build list of singly excited beta + ! --------------------------------- + + m_b = psi_bilinear_matrix_order_reverse(l_a) + m_b = m_b+1 + j=1 + do while ( (mrow == lrow) ) + mcol = psi_bilinear_matrix_transp_columns(m_b) + buffer(1:N_int,j) = psi_det_beta_unique(1:N_int,mcol) + idx(j) = mcol + j = j+1 + m_b = m_b+1 + if (m_b <= N_det) exit + mrow = psi_bilinear_matrix_transp_rows(m_b) + enddo + j=j-1 + + call get_all_spin_singles( & + buffer, idx, tmp_det(1,2), N_int, j, & + doubles, n_doubles) + + ! Compute Hij for all doubles + ! --------------------------- + + m_b = psi_bilinear_matrix_order(l_a)+1 + mcol = psi_bilinear_matrix_transp_columns(m_b) + do j=1,n_doubles + tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, doubles(j) ) + call i_H_j_double_alpha_beta(tmp_det,tmp_det2,N_int,hij) + do while (mcol /= doubles(j)) + m_b = m_b+1 + if (m_b > N_det) exit + mcol = psi_bilinear_matrix_transp_columns(m_b) + enddo + m_a = psi_bilinear_matrix_order_reverse(m_b) + v_0(m_a, 1:N_st) += hij * psi_bilinear_matrix_values(k_a,1:N_st) + v_0(k_a, 1:N_st) += hij * psi_bilinear_matrix_values(m_a,1:N_st) + enddo + + endif + l_a = l_a+1 + if (l_a > N_det) exit + lrow = psi_bilinear_matrix_rows(l_a) + lcol = psi_bilinear_matrix_columns(l_a) + enddo + + k_b = k_b+1 + if (k_b > N_det) exit + kcol = psi_bilinear_matrix_transp_columns(k_b) + enddo + + do j=1,n_singles + is_single_a( singles(j) ) = .False. + enddo + + enddo + + +end + + +subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8) + use bitmasks + implicit none + integer(bit_kind), intent(in) :: keys_tmp(Nint,2,n) + integer, intent(in) :: N_st,n,Nint, sze_8 + double precision, intent(out) :: v_0(sze_8,N_st), s_0(sze_8,N_st) + double precision, intent(in) :: u_0(sze_8,N_st) + double precision, intent(in) :: H_jj(n), S2_jj(n) + + PROVIDE ref_bitmask_energy + + double precision, allocatable :: vt(:,:) + integer, allocatable :: idx(:) + integer :: i,j, jj + double precision :: hij + + do i=1,n + v_0(i,:) = H_jj(i) * u_0(i,:) + enddo + + allocate(idx(0:n), vt(N_st,n)) + Vt = 0.d0 + do i=1,n + idx(0) = i + call filter_connected(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx) + do jj=1,idx(0) + j = idx(jj) + call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) + vt (:,i) = vt (:,i) + hij*u_0(j,:) + vt (:,j) = vt (:,j) + hij*u_0(i,:) + enddo + enddo + do i=1,n + v_0(i,:) = v_0(i,:) + vt(:,i) + enddo +end + diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 45458544..126e5c8f 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -2556,3 +2556,89 @@ subroutine get_mono_excitation_spin(det1,det2,exc,phase,Nint) enddo end +subroutine i_H_j_mono_spin(key_i,key_j,Nint,spin,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants differing by a single excitation + END_DOC + integer, intent(in) :: Nint, spin + integer(bit_kind), intent(in) :: key_i(Nint), key_j(Nint) + double precision, intent(out) :: hij + + integer :: exc(0:2,2) + double precision :: phase + + PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map + + call get_mono_excitation_spin(key_i,key_j,exc,phase,Nint) +if (exc(1,1) == 0) then + call debug_spindet(key_i, N_int) + call debug_spindet(key_j, N_int) + print *, exc(0:2,1:2) + stop +endif + call get_mono_excitation_from_fock(key_i,key_j,exc(1,2),exc(1,1),spin,phase,hij) +end + +subroutine i_H_j_double_spin(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants differing by a same-spin double excitation + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint), key_j(Nint) + double precision, intent(out) :: hij + + integer :: exc(0:2,2) + double precision :: phase + double precision, external :: get_mo_bielec_integral + + PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map + + call get_double_excitation_spin(key_i,key_j,exc,phase,Nint) + hij = phase*(get_mo_bielec_integral( & + exc(1,1), & + exc(2,1), & + exc(1,2), & + exc(2,2), mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1), & + exc(2,1), & + exc(2,2), & + exc(1,2), mo_integrals_map) ) +end + +subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants differing by an opposite-spin double excitation + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + double precision :: phase + double precision, external :: get_mo_bielec_integral + + PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map + + call get_mono_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint) + call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase,Nint) + if (exc(1,1,1) == exc(1,2,2)) then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) == exc(1,1,2)) then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif +end + + diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index da56388d..72218bc6 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -401,7 +401,6 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) 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 @@ -437,6 +436,7 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_ &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) ] +&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_reverse , (N_det) ] use bitmasks implicit none BEGIN_DOC @@ -474,6 +474,9 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_ do l=1,N_states call dset_order(psi_bilinear_matrix_transp_values(1,l),psi_bilinear_matrix_transp_order,N_det) enddo + do k=1,N_det + psi_bilinear_matrix_order_reverse(psi_bilinear_matrix_transp_order(k)) = k + enddo deallocate(to_sort) END_PROVIDER @@ -559,7 +562,7 @@ subroutine generate_all_alpha_beta_det_products ! Create a wave function from all possible alpha x beta determinants END_DOC integer :: i,j,k,l - integer :: idx, iproc + integer :: iproc integer, external :: get_index_in_psi_det_sorted_bit integer(bit_kind), allocatable :: tmp_det(:,:,:) logical, external :: is_in_wavefunction @@ -568,7 +571,7 @@ subroutine generate_all_alpha_beta_det_products !$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) + !$OMP PRIVATE(i,j,k,l,tmp_det,iproc) !$ iproc = omp_get_thread_num() allocate (tmp_det(N_int,2,N_det_alpha_unique)) !$OMP DO @@ -595,7 +598,7 @@ end -subroutine get_all_spin_singles_and_doubles(buffer, spindet, Nint, size_buffer, singles, doubles, n_singles, n_doubles) +subroutine get_all_spin_singles_and_doubles(buffer, idx, spindet, Nint, size_buffer, singles, doubles, n_singles, n_doubles) use bitmasks implicit none BEGIN_DOC @@ -606,7 +609,7 @@ subroutine get_all_spin_singles_and_doubles(buffer, spindet, Nint, size_buffer, ! /!\ : The buffer is transposed ! ! END_DOC - integer, intent(in) :: Nint, size_buffer + integer, intent(in) :: Nint, size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(Nint,size_buffer) integer(bit_kind), intent(in) :: spindet(Nint) integer, intent(out) :: singles(size_buffer) @@ -625,13 +628,13 @@ subroutine get_all_spin_singles_and_doubles(buffer, spindet, Nint, size_buffer, select case (Nint) case (1) - call get_all_spin_singles_and_doubles_1(buffer, spindet, size_buffer, singles, doubles, n_singles, n_doubles) + call get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) return case (2) - call get_all_spin_singles_and_doubles_2(buffer, spindet, size_buffer, singles, doubles, n_singles, n_doubles) + call get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) return case (3) - call get_all_spin_singles_and_doubles_3(buffer, spindet, size_buffer, singles, doubles, n_singles, n_doubles) + call get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) return end select @@ -667,11 +670,11 @@ subroutine get_all_spin_singles_and_doubles(buffer, spindet, Nint, size_buffer, n_doubles = 1 do i=1,size_buffer if ( degree(i) == 4 ) then - doubles(n_doubles) = i + doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 endif if ( degree(i) == 2 ) then - singles(n_singles) = i + singles(n_singles) = idx(i) n_singles = n_singles+1 endif enddo @@ -682,7 +685,7 @@ subroutine get_all_spin_singles_and_doubles(buffer, spindet, Nint, size_buffer, end -subroutine get_all_spin_singles(buffer, spindet, Nint, size_buffer, singles, n_singles) +subroutine get_all_spin_singles(buffer, idx, spindet, Nint, size_buffer, singles, n_singles) use bitmasks implicit none BEGIN_DOC @@ -691,7 +694,7 @@ subroutine get_all_spin_singles(buffer, spindet, Nint, size_buffer, singles, n_s ! unique alpha determinants. ! END_DOC - integer, intent(in) :: Nint, size_buffer + integer, intent(in) :: Nint, size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(Nint,size_buffer) integer(bit_kind), intent(in) :: spindet(Nint) integer, intent(out) :: singles(size_buffer) @@ -708,13 +711,13 @@ subroutine get_all_spin_singles(buffer, spindet, Nint, size_buffer, singles, n_s select case (Nint) case (1) - call get_all_spin_singles_1(buffer, spindet, size_buffer, singles, n_singles) + call get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_singles) return case (2) - call get_all_spin_singles_2(buffer, spindet, size_buffer, singles, n_singles) + call get_all_spin_singles_2(buffer, idx, spindet, size_buffer, singles, n_singles) return case (3) - call get_all_spin_singles_3(buffer, spindet, size_buffer, singles, n_singles) + call get_all_spin_singles_3(buffer, idx, spindet, size_buffer, singles, n_singles) return end select @@ -748,7 +751,7 @@ subroutine get_all_spin_singles(buffer, spindet, Nint, size_buffer, singles, n_s n_singles = 1 do i=1,size_buffer if ( degree(i) == 2 ) then - singles(n_singles) = i + singles(n_singles) = idx(i) n_singles = n_singles+1 endif enddo @@ -758,7 +761,7 @@ subroutine get_all_spin_singles(buffer, spindet, Nint, size_buffer, singles, n_s end -subroutine get_all_spin_doubles(buffer, spindet, Nint, size_buffer, doubles, n_doubles) +subroutine get_all_spin_doubles(buffer, idx, spindet, Nint, size_buffer, doubles, n_doubles) use bitmasks implicit none BEGIN_DOC @@ -767,7 +770,7 @@ subroutine get_all_spin_doubles(buffer, spindet, Nint, size_buffer, doubles, n_d ! unique alpha determinants. ! END_DOC - integer, intent(in) :: Nint, size_buffer + integer, intent(in) :: Nint, size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(Nint,size_buffer) integer(bit_kind), intent(in) :: spindet(Nint) integer, intent(out) :: doubles(size_buffer) @@ -784,13 +787,13 @@ subroutine get_all_spin_doubles(buffer, spindet, Nint, size_buffer, doubles, n_d select case (Nint) case (1) - call get_all_spin_doubles_1(buffer, spindet, size_buffer, doubles, n_doubles) + call get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_doubles) return case (2) - call get_all_spin_doubles_2(buffer, spindet, size_buffer, doubles, n_doubles) + call get_all_spin_doubles_2(buffer, idx, spindet, size_buffer, doubles, n_doubles) return case (3) - call get_all_spin_doubles_3(buffer, spindet, size_buffer, doubles, n_doubles) + call get_all_spin_doubles_3(buffer, idx, spindet, size_buffer, doubles, n_doubles) return end select @@ -824,7 +827,7 @@ subroutine get_all_spin_doubles(buffer, spindet, Nint, size_buffer, doubles, n_d n_doubles = 1 do i=1,size_buffer if ( degree(i) == 4 ) then - doubles(n_doubles) = i + doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 endif enddo @@ -833,7 +836,7 @@ subroutine get_all_spin_doubles(buffer, spindet, Nint, size_buffer, doubles, n_d end -subroutine get_all_spin_singles_and_doubles_1(buffer, spindet, size_buffer, singles, doubles, n_singles, n_doubles) +subroutine get_all_spin_singles_and_doubles_1(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) use bitmasks implicit none BEGIN_DOC @@ -845,6 +848,7 @@ subroutine get_all_spin_singles_and_doubles_1(buffer, spindet, size_buffer, sing ! END_DOC integer, intent(in) :: size_buffer + integer, intent(in) :: idx(size_buffer) integer(bit_kind), intent(in) :: buffer(size_buffer) integer(bit_kind), intent(in) :: spindet integer, intent(out) :: singles(size_buffer) @@ -872,29 +876,24 @@ subroutine get_all_spin_singles_and_doubles_1(buffer, spindet, size_buffer, sing n_doubles = 1 do i=1,size_buffer - if (xorvec(i) /= 0_8) then - degree = popcnt(xorvec(i)) - else - degree = 0 - endif - + degree = popcnt(xorvec(i)) if ( degree == 4 ) then - doubles(n_doubles) = i + doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 endif if ( degree == 2 ) then - singles(n_singles) = i + singles(n_singles) = idx(i) n_singles = n_singles+1 endif enddo n_singles = n_singles-1 n_doubles = n_doubles-1 + deallocate(xorvec) - end -subroutine get_all_spin_singles_1(buffer, spindet, size_buffer, singles, n_singles) +subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_singles) use bitmasks implicit none BEGIN_DOC @@ -903,7 +902,7 @@ subroutine get_all_spin_singles_1(buffer, spindet, size_buffer, singles, n_singl ! unique alpha determinants. ! END_DOC - integer, intent(in) :: size_buffer + integer, intent(in) :: size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(size_buffer) integer(bit_kind), intent(in) :: spindet integer, intent(out) :: singles(size_buffer) @@ -921,7 +920,7 @@ subroutine get_all_spin_singles_1(buffer, spindet, size_buffer, singles, n_singl n_singles = 1 do i=1,size_buffer if ( popcnt(xorvec(i)) == 2 ) then - singles(n_singles) = i + singles(n_singles) = idx(i) n_singles = n_singles+1 endif enddo @@ -931,7 +930,7 @@ subroutine get_all_spin_singles_1(buffer, spindet, size_buffer, singles, n_singl end -subroutine get_all_spin_doubles_1(buffer, spindet, size_buffer, doubles, n_doubles) +subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_doubles) use bitmasks implicit none BEGIN_DOC @@ -940,7 +939,7 @@ subroutine get_all_spin_doubles_1(buffer, spindet, size_buffer, doubles, n_doubl ! unique alpha determinants. ! END_DOC - integer, intent(in) :: size_buffer + integer, intent(in) :: size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(size_buffer) integer(bit_kind), intent(in) :: spindet integer, intent(out) :: doubles(size_buffer) @@ -961,7 +960,7 @@ subroutine get_all_spin_doubles_1(buffer, spindet, size_buffer, doubles, n_doubl do i=1,size_buffer if ( popcnt(xorvec(i)) == 4 ) then - doubles(n_doubles) = i + doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 endif enddo @@ -971,7 +970,7 @@ subroutine get_all_spin_doubles_1(buffer, spindet, size_buffer, doubles, n_doubl end -subroutine get_all_spin_singles_and_doubles_2(buffer, spindet, size_buffer, singles, doubles, n_singles, n_doubles) +subroutine get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) use bitmasks implicit none BEGIN_DOC @@ -982,7 +981,7 @@ subroutine get_all_spin_singles_and_doubles_2(buffer, spindet, size_buffer, sing ! /!\ : The buffer is transposed ! ! END_DOC - integer, intent(in) :: size_buffer + integer, intent(in) :: size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(2,size_buffer) integer(bit_kind), intent(in) :: spindet(2) integer, intent(out) :: singles(size_buffer) @@ -1027,11 +1026,11 @@ subroutine get_all_spin_singles_and_doubles_2(buffer, spindet, size_buffer, sing n_doubles = 1 do i=1,size_buffer if ( degree(i) == 4 ) then - doubles(n_doubles) = i + doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 endif if ( degree(i) == 2 ) then - singles(n_singles) = i + singles(n_singles) = idx(i) n_singles = n_singles+1 endif enddo @@ -1042,7 +1041,7 @@ subroutine get_all_spin_singles_and_doubles_2(buffer, spindet, size_buffer, sing end -subroutine get_all_spin_singles_2(buffer, spindet, size_buffer, singles, n_singles) +subroutine get_all_spin_singles_2(buffer, idx, spindet, size_buffer, singles, n_singles) use bitmasks implicit none BEGIN_DOC @@ -1051,7 +1050,7 @@ subroutine get_all_spin_singles_2(buffer, spindet, size_buffer, singles, n_singl ! unique alpha determinants. ! END_DOC - integer, intent(in) :: size_buffer + integer, intent(in) :: size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(2,size_buffer) integer(bit_kind), intent(in) :: spindet(2) integer, intent(out) :: singles(size_buffer) @@ -1093,7 +1092,7 @@ subroutine get_all_spin_singles_2(buffer, spindet, size_buffer, singles, n_singl n_singles = 1 do i=1,size_buffer if ( degree(i) == 2 ) then - singles(n_singles) = i + singles(n_singles) = idx(i) n_singles = n_singles+1 endif enddo @@ -1103,7 +1102,7 @@ subroutine get_all_spin_singles_2(buffer, spindet, size_buffer, singles, n_singl end -subroutine get_all_spin_doubles_2(buffer, spindet, size_buffer, doubles, n_doubles) +subroutine get_all_spin_doubles_2(buffer, idx, spindet, size_buffer, doubles, n_doubles) use bitmasks implicit none BEGIN_DOC @@ -1113,7 +1112,7 @@ subroutine get_all_spin_doubles_2(buffer, spindet, size_buffer, doubles, n_doubl ! END_DOC integer, intent(in) :: size_buffer - integer(bit_kind), intent(in) :: buffer(2,size_buffer) + integer(bit_kind), intent(in) :: buffer(2,size_buffer), idx(size_buffer) integer(bit_kind), intent(in) :: spindet(2) integer, intent(out) :: doubles(size_buffer) integer, intent(out) :: n_doubles @@ -1154,7 +1153,7 @@ subroutine get_all_spin_doubles_2(buffer, spindet, size_buffer, doubles, n_doubl n_doubles = 1 do i=1,size_buffer if ( degree(i) == 4 ) then - doubles(n_doubles) = i + doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 endif enddo @@ -1163,7 +1162,7 @@ subroutine get_all_spin_doubles_2(buffer, spindet, size_buffer, doubles, n_doubl end -subroutine get_all_spin_singles_and_doubles_3(buffer, spindet, size_buffer, singles, doubles, n_singles, n_doubles) +subroutine get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles) use bitmasks implicit none BEGIN_DOC @@ -1174,7 +1173,7 @@ subroutine get_all_spin_singles_and_doubles_3(buffer, spindet, size_buffer, sing ! /!\ : The buffer is transposed ! ! END_DOC - integer, intent(in) :: size_buffer + integer, intent(in) :: size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(3,size_buffer) integer(bit_kind), intent(in) :: spindet(3) integer, intent(out) :: singles(size_buffer) @@ -1226,11 +1225,11 @@ subroutine get_all_spin_singles_and_doubles_3(buffer, spindet, size_buffer, sing n_doubles = 1 do i=1,size_buffer if ( degree(i) == 4 ) then - doubles(n_doubles) = i + doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 endif if ( degree(i) == 2 ) then - singles(n_singles) = i + singles(n_singles) = idx(i) n_singles = n_singles+1 endif enddo @@ -1241,7 +1240,7 @@ subroutine get_all_spin_singles_and_doubles_3(buffer, spindet, size_buffer, sing end -subroutine get_all_spin_singles_3(buffer, spindet, size_buffer, singles, n_singles) +subroutine get_all_spin_singles_3(buffer, idx, spindet, size_buffer, singles, n_singles) use bitmasks implicit none BEGIN_DOC @@ -1250,7 +1249,7 @@ subroutine get_all_spin_singles_3(buffer, spindet, size_buffer, singles, n_singl ! unique alpha determinants. ! END_DOC - integer, intent(in) :: size_buffer + integer, intent(in) :: size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(3,size_buffer) integer(bit_kind), intent(in) :: spindet(3) integer, intent(out) :: singles(size_buffer) @@ -1299,7 +1298,7 @@ subroutine get_all_spin_singles_3(buffer, spindet, size_buffer, singles, n_singl n_singles = 1 do i=1,size_buffer if ( degree(i) == 2 ) then - singles(n_singles) = i + singles(n_singles) = idx(i) n_singles = n_singles+1 endif enddo @@ -1309,7 +1308,7 @@ subroutine get_all_spin_singles_3(buffer, spindet, size_buffer, singles, n_singl end -subroutine get_all_spin_doubles_3(buffer, spindet, size_buffer, doubles, n_doubles) +subroutine get_all_spin_doubles_3(buffer, idx, spindet, size_buffer, doubles, n_doubles) use bitmasks implicit none BEGIN_DOC @@ -1318,7 +1317,7 @@ subroutine get_all_spin_doubles_3(buffer, spindet, size_buffer, doubles, n_doubl ! unique alpha determinants. ! END_DOC - integer, intent(in) :: size_buffer + integer, intent(in) :: size_buffer, idx(size_buffer) integer(bit_kind), intent(in) :: buffer(3,size_buffer) integer(bit_kind), intent(in) :: spindet(3) integer, intent(out) :: doubles(size_buffer) @@ -1367,7 +1366,7 @@ subroutine get_all_spin_doubles_3(buffer, spindet, size_buffer, doubles, n_doubl n_doubles = 1 do i=1,size_buffer if ( degree(i) == 4 ) then - doubles(n_doubles) = i + doubles(n_doubles) = idx(i) n_doubles = n_doubles+1 endif enddo