diff --git a/plugins/MRCC_Utils/README.rst b/plugins/MRCC_Utils/README.rst index 62eba2d1..56b519f2 100644 --- a/plugins/MRCC_Utils/README.rst +++ b/plugins/MRCC_Utils/README.rst @@ -258,7 +258,11 @@ Documentation N_states lowest eigenvalues of the dressed CI matrix -`davidson_diag_hjj_mrcc `_ +`create_minilist `_ + Undocumented + + +`davidson_diag_hjj_mrcc `_ Davidson diagonalization with specific diagonal elements of the H matrix .br H_jj : specific diagonal H matrix elements to diagonalize de Davidson @@ -279,7 +283,7 @@ Documentation Initial guess vectors are not necessarily orthonormal -`davidson_diag_mrcc `_ +`davidson_diag_mrcc `_ Davidson diagonalization. .br dets_in : bitmasks corresponding to determinants @@ -370,7 +374,7 @@ Documentation Find A.C = B -`find_triples_and_quadruples `_ +`find_triples_and_quadruples `_ Undocumented @@ -431,18 +435,26 @@ Documentation Undocumented -`h_apply_mrcc `_ +`h_apply_mrcc `_ Calls H_apply on the HF determinant and selects all connected single and double excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. -`h_apply_mrcc_diexc `_ +`h_apply_mrcc_diexc `_ + Undocumented + + +`h_apply_mrcc_diexcorg `_ Generate all double excitations of key_in using the bit masks of holes and particles. Assume N_int is already provided. -`h_apply_mrcc_monoexc `_ +`h_apply_mrcc_diexcp `_ + Undocumented + + +`h_apply_mrcc_monoexc `_ Generate all single excitations of key_in using the bit masks of holes and particles. Assume N_int is already provided. @@ -452,7 +464,15 @@ Documentation Dressed H with Delta_ij -`h_u_0_mrcc `_ +`h_u_0_mrcc `_ + Computes v_0 = H|u_0> + .br + n : number of determinants + .br + H_jj : array of + + +`h_u_0_mrcc_org `_ Computes v_0 = H|u_0> .br n : number of determinants @@ -757,11 +777,11 @@ Documentation n! -`mrcc_dress `_ +`mrcc_dress `_ Undocumented -`mrcc_dress_simple `_ +`mrcc_dress_simple `_ Undocumented @@ -855,7 +875,7 @@ Documentation Current status for displaying progress bars. Global variable. -`psi_ref_lock `_ +`psi_ref_lock `_ Locks on ref determinants to fill delta_ij diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index d9e697a1..2fa09575 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -1,3 +1,6 @@ + + + subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,istate) use bitmasks implicit none @@ -102,6 +105,9 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin double precision :: to_print(2,N_st) double precision :: cpu, wall + integer(bit_kind) :: dets_in_sorted(Nint,2,sze) + integer :: idx(sze), shortcut(0:sze+1),sh,ii,tmp + PROVIDE det_connections call write_time(iunit) @@ -147,6 +153,10 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin ! Initialization ! ============== + + dets_in_sorted(:,:,:) = dets_in(:,:,:) + call sort_dets_ab(dets_in_sorted, idx, shortcut, sze, Nint) + k_pairs=0 do l=1,N_st do k=1,l @@ -205,7 +215,8 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin ! ---------------------- do k=1,N_st - call H_u_0_mrcc(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint,istate) + call H_u_0_mrcc(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in_sorted,shortcut,idx,Nint,istate) + !call H_u_0_mrcc_org(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint,istate) enddo ! Compute h_kl = = @@ -357,7 +368,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin abort_here = abort_all end -subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) +subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,shortcut,sort_idx,Nint,istate) use bitmasks implicit none BEGIN_DOC @@ -377,6 +388,115 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) double precision, allocatable :: vt(:) integer :: i,j,k,l, jj,ii integer :: i0, j0 + + integer,intent(in) :: shortcut(0:n+1), sort_idx(n) + integer :: tmp, warp(2,0:n+1), sh, ni +! + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (n>0) + PROVIDE ref_bitmask_energy delta_ij + integer, parameter :: block_size = 157 + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,warp,tmp,sh) & + !$OMP SHARED(n_det_ref,n_det_non_ref,idx_ref,idx_non_ref,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij,shortcut,sort_idx) + + !$OMP DO SCHEDULE(static) + do i=1,n + v_0(i) = H_jj(i) * u_0(i) + enddo + !$OMP END DO + + allocate(idx(0:n), vt(n)) + Vt = 0.d0 + !$OMP DO SCHEDULE(dynamic) + do sh=1,shortcut(0) + warp(1,0) = 0 + do ii=1,sh!shortcut(0) + tmp = 0 + do ni=1,Nint + tmp = popcnt(xor(keys_tmp(ni,1, shortcut(ii)), keys_tmp(ni,1,shortcut(sh)))) + end do + if(tmp <= 4) then + tmp = warp(1,0) + 1 + warp(1,0) = tmp + warp(1,tmp) = shortcut(ii) + warp(2,tmp) = shortcut(ii+1)-1 + end if + end do + + do ii=shortcut(sh),shortcut(sh+1)-1 + idx(0) = ii + + + !call filter_connected_davidson_mwen(keys_tmp,shortcut,keys_tmp(1,1,ii),Nint,ii-1,idx) + call filter_connected_davidson_warp(keys_tmp,warp,keys_tmp(1,1,ii),Nint,ii-1,idx) + i = sort_idx(ii) + + do jj=1,idx(0) + j = sort_idx(idx(jj)) + !j = idx(jj) + if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then + call i_H_j(keys_tmp(1,1,idx(jj)),keys_tmp(1,1,ii),Nint,hij) + vt (i) = vt (i) + hij*u_0(j) + vt (j) = vt (j) + hij*u_0(i) + endif + enddo + enddo + enddo + !$OMP END DO + + + + !$OMP DO SCHEDULE(guided) + do ii=1,n_det_ref + i = idx_ref(ii) + do jj = 1, n_det_non_ref + j = idx_non_ref(jj) + vt (i) = vt (i) + delta_ij(ii,jj,istate)*u_0(j) + vt (j) = vt (j) + delta_ij(ii,jj,istate)*u_0(i) + enddo + enddo + !$OMP END DO + !$OMP CRITICAL + do i=1,n + v_0(i) = v_0(i) + vt(i) + enddo + !$OMP END CRITICAL + deallocate(idx,vt) + !$OMP END PARALLEL +end + + + +subroutine H_u_0_mrcc_org(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> + ! + ! n : number of determinants + ! + ! H_jj : array of + END_DOC + integer, intent(in) :: n,Nint,istate + double precision, intent(out) :: v_0(n) + double precision, intent(in) :: u_0(n) + double precision, intent(in) :: H_jj(n) + integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) + integer, allocatable :: idx(:) + double precision :: hij + double precision, allocatable :: vt(:) + integer :: i,j,k,l, jj,ii + integer :: i0, j0 + + + + + ASSERT (Nint > 0) ASSERT (Nint == N_int) ASSERT (n>0) @@ -427,4 +547,3 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate) !$OMP END PARALLEL end - diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index 3a4717ed..f427c81f 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -1,4 +1,5 @@ use omp_lib +use bitmasks BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ] implicit none @@ -94,7 +95,6 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n integer :: idx_miniList(N_det_non_ref), N_miniList - call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) allocate (dIa_hla(N_states,Ndet_non_ref)) @@ -285,6 +285,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq integer, intent(out) :: N_tq integer :: c_ref integer :: connected_to_ref + N_tq = 0 do i=1,N_selected diff --git a/plugins/Perturbation/README.rst b/plugins/Perturbation/README.rst index 507e843c..aa6ebf54 100644 --- a/plugins/Perturbation/README.rst +++ b/plugins/Perturbation/README.rst @@ -112,42 +112,42 @@ Documentation routine. -`perturb_buffer_by_mono_dipole_moment_z `_ +`perturb_buffer_by_mono_dipole_moment_z `_ Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_by_mono_epstein_nesbet `_ +`perturb_buffer_by_mono_epstein_nesbet `_ Applly pertubration ``epstein_nesbet`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_by_mono_epstein_nesbet_2x2 `_ +`perturb_buffer_by_mono_epstein_nesbet_2x2 `_ Applly pertubration ``epstein_nesbet_2x2`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_by_mono_epstein_nesbet_sc2 `_ +`perturb_buffer_by_mono_epstein_nesbet_sc2 `_ Applly pertubration ``epstein_nesbet_sc2`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_by_mono_epstein_nesbet_sc2_no_projected `_ +`perturb_buffer_by_mono_epstein_nesbet_sc2_no_projected `_ Applly pertubration ``epstein_nesbet_sc2_no_projected`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_by_mono_epstein_nesbet_sc2_projected `_ +`perturb_buffer_by_mono_epstein_nesbet_sc2_projected `_ Applly pertubration ``epstein_nesbet_sc2_projected`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_by_mono_h_core `_ +`perturb_buffer_by_mono_h_core `_ Applly pertubration ``h_core`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_by_mono_moller_plesset `_ +`perturb_buffer_by_mono_moller_plesset `_ Applly pertubration ``moller_plesset`` to the buffer of determinants generated in the H_apply routine. @@ -157,42 +157,42 @@ Documentation routine. -`perturb_buffer_dipole_moment_z `_ +`perturb_buffer_dipole_moment_z `_ Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_epstein_nesbet `_ +`perturb_buffer_epstein_nesbet `_ Applly pertubration ``epstein_nesbet`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_epstein_nesbet_2x2 `_ +`perturb_buffer_epstein_nesbet_2x2 `_ Applly pertubration ``epstein_nesbet_2x2`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_epstein_nesbet_sc2 `_ +`perturb_buffer_epstein_nesbet_sc2 `_ Applly pertubration ``epstein_nesbet_sc2`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_epstein_nesbet_sc2_no_projected `_ +`perturb_buffer_epstein_nesbet_sc2_no_projected `_ Applly pertubration ``epstein_nesbet_sc2_no_projected`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_epstein_nesbet_sc2_projected `_ +`perturb_buffer_epstein_nesbet_sc2_projected `_ Applly pertubration ``epstein_nesbet_sc2_projected`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_h_core `_ +`perturb_buffer_h_core `_ Applly pertubration ``h_core`` to the buffer of determinants generated in the H_apply routine. -`perturb_buffer_moller_plesset `_ +`perturb_buffer_moller_plesset `_ Applly pertubration ``moller_plesset`` to the buffer of determinants generated in the H_apply routine. diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 3d1de6b5..fe360a96 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -1,13 +1,15 @@ -subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters ) +subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters ) integer(bit_kind), intent(in) :: key_in(N_int, 2), hole_1(N_int, 2), hole_2(N_int, 2) integer(bit_kind), intent(in) :: particl_1(N_int, 2), particl_2(N_int, 2) - integer(bit_kind) :: p1_mask(N_int, 2), p2_mask(N_int, 2) + integer(bit_kind) :: p1_mask(N_int, 2), p2_mask(N_int, 2), tmp integer,intent(in) :: i_generator,iproc_in integer(bit_kind) :: status(N_int*bit_kind_size, 2) - integer :: highest, p1,p2,sp,ni,i,mi + integer :: highest, p1,p2,sp,ni,i,mi,nt,ns + + integer(bit_kind), intent(in) :: key_prev(N_int, 2, *) $declarations @@ -29,6 +31,32 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene end do end do +! GEL D'ELECTRONS +! nt = 0 +! do i=1,i_generator-1 +! if(key_in(1,1) == key_prev(1,1,i)) then +! tmp = xor(key_in(1,2), key_prev(1,2,i)) +! if(popcnt(tmp) == 2) then +! ns = 1+trailz(iand(tmp, key_in(1,2))) +! if(status(ns, 2) /= 0) then +! nt += 1 +! end if +! status(ns, 2) = 0 +! end if +! else if(key_in(1,2) == key_prev(1,2,i)) then +! tmp = xor(key_in(1,1), key_prev(1,1,i)) +! if(popcnt(tmp) == 2) then +! ns = 1+trailz(iand(tmp, key_in(1,1))) +! if(status(ns, 1) /= 0) then +! nt += 1 +! end if +! status(ns, 1) = 0 +! end if +! end if +! end do +! print *, "nt", nt, i_generator + + do sp=1,2 do p1=1,highest if(status(p1, sp) == 0) then @@ -76,35 +104,13 @@ subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, p1_mask(:,:) = 0_bit_kind p2_mask(:,:) = 0_bit_kind - p1_mask(fh1/bit_kind_size + 1, fs1) = 2**(mod(fh1-1,bit_kind_size)) - p2_mask(fh2/bit_kind_size + 1, fs2) = 2**(mod(fh2-1,bit_kind_size)) + p1_mask(ishft(fh1,-bit_kind_shift) + 1, fs1) = ishft(1,iand(fh1-1,bit_kind_size-1)) + p2_mask(ishft(fh2,-bit_kind_shift) + 1, fs2) = ishft(1,iand(fh2-1,bit_kind_size-1)) -! n_alpha = 0 -! n_beta = 0 key_mask(:,:) = key_in(:,:) - key_mask(fh1/bit_kind_size + 1, fs1) -= 2**(mod(fh1-1,bit_kind_size)) - key_mask(fh2/bit_kind_size + 1, fs2) -= 2**(mod(fh2-1,bit_kind_size)) -! do i=1,N_int -! n_alpha = n_alpha + popcnt(key_mask(i, 1)) -! n_beta = n_beta + popcnt(key_mask(i, 2)) -! end do -! -! do i=1, N_det -! deg(1) = n_alpha -! deg(2) = n_beta -! -! do ni = 1, N_int -! ! deg(1) = deg(1) - popcnt(iand(key_mask(ni, 1), psi_non_ref(ni, 1, i))) -! ! deg(2) = deg(2) - popcnt(iand(key_mask(ni, 2), psi_non_ref(ni, 2, i))) -! end do -! -! -! if(deg(1) + deg(2) <= 2) then -! ! ndet_out = ndet_out + 1 -! ! idx(ndet_out) = i -! end if -! end do + key_mask(ishft(fh1,-bit_kind_shift) + 1, fs1) -= ishft(1,iand(fh1-1,bit_kind_size-1)) + key_mask(ishft(fh2,-bit_kind_shift) + 1, fs2) -= ishft(1,iand(fh2-1,bit_kind_size-1)) call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, i_generator, iproc_in $parameters ) end subroutine @@ -529,7 +535,6 @@ subroutine $subroutine($params_main) nmax = mod( N_det_generators,nproc ) - !$ call omp_init_lock(lck) call start_progress(N_det_generators,'Selection (norm)',0.d0) @@ -572,6 +577,7 @@ subroutine $subroutine($params_main) enddo if($do_double_excitations)then call $subroutine_diexc(psi_det_generators(1,1,i_generator), & + psi_det_generators(1,1,1), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator, iproc $params_post) @@ -632,6 +638,7 @@ subroutine $subroutine($params_main) if($do_double_excitations)then call $subroutine_diexc(psi_det_generators(1,1,i_generator), & + psi_det_generators(1,1,1), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator, iproc $params_post) diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index 62b035fc..d12d8426 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -54,7 +54,7 @@ Documentation .. by the `update_README.py` script. -`a_operator `_ +`a_operator `_ Needed for diag_H_mat_elem @@ -66,7 +66,7 @@ Documentation Max and min values of the coefficients -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem @@ -157,11 +157,11 @@ Documentation of alpha and beta determinants -`davidson_converged `_ +`davidson_converged `_ True if the Davidson algorithm is converged -`davidson_criterion `_ +`davidson_criterion `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] @@ -184,7 +184,7 @@ Documentation Initial guess vectors are not necessarily orthonormal -`davidson_diag_hjj `_ +`davidson_diag_hjj `_ Davidson diagonalization with specific diagonal elements of the H matrix .br H_jj : specific diagonal H matrix elements to diagonalize de Davidson @@ -213,7 +213,7 @@ Documentation Max number of Davidson sizes -`davidson_threshold `_ +`davidson_threshold `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] @@ -229,10 +229,14 @@ Documentation det_coef -`det_connections `_ +`det_connections `_ Build connection proxy between determinants +`det_inf `_ + Undocumented + + `det_num `_ det_num @@ -253,7 +257,7 @@ Documentation Diagonalization algorithm (Davidson or Lapack) -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes @@ -313,7 +317,7 @@ Documentation idx(0) is the number of determinants that interact with key1 -`filter_connected_davidson `_ +`filter_connected_davidson `_ Filters out the determinants that are not connected by H returns the array idx which contains the index of the determinants in the array key1 that interact @@ -323,7 +327,27 @@ Documentation key1 should come from psi_det_sorted_ab. -`filter_connected_i_h_psi0 `_ +`filter_connected_davidson_shortcut `_ + Filters out the determinants that are not connected by H + returns the array idx which contains the index of the + determinants in the array key1 that interact + via the H operator with key2. + .br + idx(0) is the number of determinants that interact with key1 + key1 should come from psi_det_sorted_ab. + + +`filter_connected_davidson_warp `_ + Filters out the determinants that are not connected by H + returns the array idx which contains the index of the + determinants in the array key1 that interact + via the H operator with key2. + .br + idx(0) is the number of determinants that interact with key1 + key1 should come from psi_det_sorted_ab. + + +`filter_connected_i_h_psi0 `_ returns the array idx which contains the index of the .br determinants in the array key1 that interact @@ -333,7 +357,7 @@ Documentation idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0_sc2 `_ +`filter_connected_i_h_psi0_sc2 `_ standard filter_connected_i_H_psi but returns in addition .br the array of the index of the non connected determinants to key1 @@ -359,7 +383,7 @@ Documentation Create a wave function from all possible alpha x beta determinants -`get_double_excitation `_ +`get_double_excitation `_ Returns the two excitation operators between two doubly excited determinants and the phase @@ -371,7 +395,7 @@ Documentation Returns the excitation degree between two determinants -`get_excitation_degree_vector `_ +`get_excitation_degree_vector `_ Applies get_excitation_degree to an array of determinants @@ -387,11 +411,11 @@ Documentation Returns the index of the determinant in the ``psi_det_sorted_bit`` array -`get_mono_excitation `_ +`get_mono_excitation `_ Returns the excitation operator between two singly excited determinants and the phase -`get_occ_from_key `_ +`get_occ_from_key `_ Returns a list of occupation numbers from a bitstring @@ -425,7 +449,7 @@ Documentation Undocumented -`h_u_0 `_ +`h_u_0 `_ Computes v_0 = H|u_0> .br n : number of determinants @@ -433,23 +457,31 @@ Documentation H_jj : array of -`i_h_j `_ +`h_u_0_org `_ + Computes v_0 = H|u_0> + .br + n : number of determinants + .br + H_jj : array of + + +`i_h_j `_ Returns where i and j are determinants -`i_h_j_phase_out `_ +`i_h_j_phase_out `_ Returns where i and j are determinants -`i_h_j_verbose `_ +`i_h_j_verbose `_ Returns where i and j are determinants -`i_h_psi `_ +`i_h_psi `_ for the various Nstates -`i_h_psi_sc2 `_ +`i_h_psi_sc2 `_ for the various Nstate .br returns in addition @@ -463,7 +495,7 @@ Documentation to repeat the excitations -`i_h_psi_sc2_verbose `_ +`i_h_psi_sc2_verbose `_ for the various Nstate .br returns in addition @@ -477,7 +509,7 @@ Documentation to repeat the excitations -`i_h_psi_sec_ord `_ +`i_h_psi_sec_ord `_ for the various Nstates @@ -524,7 +556,7 @@ Documentation Energy of the reference bitmask used in Slater rules -`n_con_int `_ +`n_con_int `_ Number of integers to represent the connections between determinants @@ -889,6 +921,10 @@ Documentation for a given couple of hole/particle excitations i. +`sort_dets_ab `_ + Undocumented + + `sort_dets_by_3_highest_electrons `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, @@ -910,6 +946,10 @@ Documentation Weights in the state-average calculation of the density matrix +`tamiser `_ + Undocumented + + `threshold_convergence_sc2 `_ convergence of the correlation energy of SC2 iterations diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f index 0c5610ae..626ecec3 100644 --- a/src/Determinants/davidson.irp.f +++ b/src/Determinants/davidson.irp.f @@ -65,6 +65,106 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit) deallocate (H_jj) end + +logical function det_inf(key1, key2, Nint) + use bitmasks + implicit none + integer(bit_kind),intent(in) :: key1(Nint, 2), key2(Nint, 2) + integer,intent(in) :: Nint + integer :: i,j + + det_inf = .false. + + do i=1,2 + do j=Nint,1,-1 + if(key1(j,i) < key2(j,i)) then + det_inf = .true. + return + else if(key1(j,i) > key2(j,i)) then + return + end if + end do + end do +end function + + +subroutine tamiser(key, idx, no, n, Nint, N_key) + use bitmasks + + implicit none + integer(bit_kind),intent(inout) :: key(Nint, 2, N_key) + integer,intent(in) :: no, n, Nint, N_key + integer,intent(inout) :: idx(N_key) + integer :: k,j,tmpidx + integer(bit_kind) :: tmp(Nint, 2) + logical :: det_inf + + k = no + j = 2*k + do while(j <= n) + if(j < n .and. det_inf(key(:,:,j), key(:,:,j+1), Nint)) then + j = j+1 + end if + if(det_inf(key(:,:,k), key(:,:,j), Nint)) then + tmp(:,:) = key(:,:,k) + key(:,:,k) = key(:,:,j) + key(:,:,j) = tmp(:,:) + tmpidx = idx(k) + idx(k) = idx(j) + idx(j) = tmpidx + k = j + j = 2*k + else + return + end if + end do +end subroutine + + +subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint) + use bitmasks + implicit none + + integer(bit_kind),intent(inout) :: key(Nint,2,N_key) + integer,intent(out) :: idx(N_key) + integer,intent(out) :: shortcut(0:N_key+1) + integer, intent(in) :: Nint, N_key + integer(bit_kind) :: tmp(Nint, 2) + integer :: tmpidx,i,ni + + do i=1,N_key + idx(i) = i + end do + + do i=N_key/2,1,-1 + call tamiser(key, idx, i, N_key, Nint, N_key) + end do + + do i=N_key,2,-1 + tmp(:,:) = key(:,:,i) + key(:,:,i) = key(:,:,1) + key(:,:,1) = tmp(:,:) + tmpidx = idx(i) + idx(i) = idx(1) + idx(1) = tmpidx + call tamiser(key, idx, 1, i-1, Nint, N_key) + end do + + shortcut(0) = 1 + shortcut(1) = 1 + do i=2,N_key + do ni=1,nint + if(key(ni,1,i) /= key(ni,1,i-1)) then + shortcut(0) = shortcut(0) + 1 + shortcut(shortcut(0)) = i + exit + end if + end do + end do + shortcut(shortcut(0)+1) = N_key+1 +end subroutine + + subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit) use bitmasks implicit none @@ -106,7 +206,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun integer :: k_pairs, kl integer :: iter2 - double precision, allocatable :: W(:,:,:), U(:,:,:), R(:,:) + double precision, allocatable :: W(:,:,:), U(:,:,:), R(:,:), Wt(:) double precision, allocatable :: y(:,:,:,:), h(:,:,:,:), lambda(:) double precision :: diag_h_mat_elem double precision :: residual_norm(N_st) @@ -114,6 +214,9 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun double precision :: to_print(2,N_st) double precision :: cpu, wall + integer(bit_kind) :: dets_in_sorted(Nint, 2, sze) + integer :: idx(sze), shortcut(0:sze+1) + PROVIDE det_connections call write_time(iunit) @@ -145,6 +248,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun allocate( & kl_pairs(2,N_st*(N_st+1)/2), & W(sze,N_st,davidson_sze_max), & + Wt(sze), & U(sze,N_st,davidson_sze_max), & R(sze,N_st), & h(N_st,davidson_sze_max,N_st,davidson_sze_max), & @@ -159,6 +263,9 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! Initialization ! ============== + dets_in_sorted(:,:,:) = dets_in(:,:,:) + call sort_dets_ab(dets_in_sorted, idx, shortcut, sze, Nint) + k_pairs=0 do l=1,N_st do k=1,l @@ -170,7 +277,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, & - !$OMP Nint,dets_in,u_in) & + !$OMP Nint,dets_in_sorted,dets_in,u_in) & !$OMP PRIVATE(k,l,kl,i) @@ -217,9 +324,11 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! ---------------------- do k=1,N_st - call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint) +! call H_u_0_org(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint) + call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in_sorted,shortcut,idx,Nint) enddo + ! Compute h_kl = = ! ------------------------------------------- @@ -274,7 +383,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun to_print(1,k) = lambda(k) + nuclear_repulsion to_print(2,k) = residual_norm(k) enddo - + write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))'), iter, to_print(:,1:N_st) call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) if (converged) then @@ -360,6 +469,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun deallocate ( & kl_pairs, & W, & + Wt, & U, & R, & h, & diff --git a/src/Determinants/filter_connected.irp.f b/src/Determinants/filter_connected.irp.f index c3d333ad..f55643bd 100644 --- a/src/Determinants/filter_connected.irp.f +++ b/src/Determinants/filter_connected.irp.f @@ -158,7 +158,160 @@ subroutine filter_connected_sorted_ab(key1,key2,next,Nint,sze,idx) end +subroutine filter_connected_davidson_warp(key1,warp,key2,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Filters out the determinants that are not connected by H + ! returns the array idx which contains the index of the + ! determinants in the array key1 that interact + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + ! key1 should come from psi_det_sorted_ab. + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer,intent(in) :: warp(2,0:sze+1) + + integer :: i,j,k,l + integer :: degree_x2 + integer :: i_alpha, i_beta, exc_a, exc_b, endloop, ni + integer(bit_kind) :: tmp1, tmp2 + + ASSERT (Nint > 0) + ASSERT (sze >= 0) + l=1 + i_alpha = 0 + + + if (Nint /= 1) then + do while(i_alpha < warp(1,0) .and. warp(1,i_alpha+1) <= sze) + i_alpha = i_alpha + 1 + exc_a = 0 + do ni=1,Nint + exc_a += popcnt(xor(key1(ni,1,warp(1,i_alpha)), key2(ni,1))) + end do + endloop = min(warp(2,i_alpha), sze) + if(exc_a == 4) then + beta_loop : do i_beta=warp(1,i_alpha),endloop + do ni=1,Nint + if(key1(ni,2,i_beta) /= key2(ni,2)) then + cycle beta_loop + end if + end do + idx(l) = i_beta + l = l + 1 + end do beta_loop + else + do i_beta=warp(1,i_alpha),endloop + exc_b = 0 + do ni=1,Nint + exc_b += popcnt(xor(key1(ni,2,i_beta), key2(ni,2))) + end do + if(exc_b + exc_a <= 4) then + idx(l) = i_beta + l = l + 1 + end if + end do + end if + end do + else + do while(i_alpha < warp(1,0) .and. warp(1,i_alpha+1) <= sze) + i_alpha = i_alpha + 1 + exc_a = popcnt(xor(key1(1,1,warp(1,i_alpha)), key2(1,1))) + endloop = min(warp(2,i_alpha), sze) + if(exc_a == 4) then + do i_beta=warp(1,i_alpha),endloop + if(key1(1,2,i_beta) == key2(1,2)) then + idx(l) = i_beta + l = l + 1 + exit + end if + end do + else + do i_beta=warp(1,i_alpha),endloop + exc_b = popcnt(xor(key1(1,2,i_beta), key2(1,2))) + if(exc_b + exc_a <= 4) then + idx(l) = i_beta + l = l + 1 + end if + end do + end if + end do + end if + + idx(0) = l-1 +end + + +subroutine filter_connected_davidson_shortcut(key1,shortcut,key2,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Filters out the determinants that are not connected by H + ! returns the array idx which contains the index of the + ! determinants in the array key1 that interact + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + ! key1 should come from psi_det_sorted_ab. + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer,intent(in) :: shortcut(0:sze+1) + + integer :: i,j,k,l + integer :: degree_x2 + integer :: i_alpha, i_beta, exc_a, exc_b, endloop + integer(bit_kind) :: tmp1, tmp2 + + ASSERT (Nint > 0) + ASSERT (sze >= 0) + + l=1 + i_alpha = 0 + + if (Nint==1) then + do while(shortcut(i_alpha+1) < sze) + i_alpha = i_alpha + 1 + exc_a = popcnt(xor(key1(1,1,shortcut(i_alpha)), key2(1,1))) + if(exc_a > 4) then + cycle + end if + endloop = min(shortcut(i_alpha+1)-1, sze) + if(exc_a == 4) then + do i_beta = shortcut(i_alpha), endloop + if(key1(1,2,i_beta) == key2(1,2)) then + idx(l) = i_beta + l = l + 1 + exit + end if + end do + else + do i_beta = shortcut(i_alpha), endloop + exc_b = popcnt(xor(key1(1,2,i_beta), key2(1,2))) + if(exc_b + exc_a <= 4) then + idx(l) = i_beta + l = l + 1 + end if + end do + end if + end do + else + print *, "TBD : filter_connected_davidson_shortcut Nint>1" + stop + end if + + idx(0) = l-1 +end subroutine filter_connected_davidson(key1,key2,Nint,sze,idx) use bitmasks @@ -183,6 +336,7 @@ subroutine filter_connected_davidson(key1,key2,Nint,sze,idx) integer*8 :: itmp PROVIDE N_con_int det_connections + ASSERT (Nint > 0) ASSERT (sze >= 0) @@ -190,7 +344,7 @@ subroutine filter_connected_davidson(key1,key2,Nint,sze,idx) if (Nint==1) then - i = idx(0) + i = idx(0) ! lecture dans un intent(out) ? do j_int=1,N_con_int itmp = det_connections(j_int,i) do while (itmp /= 0_8) diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index edda4362..72615089 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -138,6 +138,7 @@ subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) end select end + subroutine get_double_excitation(det1,det2,exc,phase,Nint) use bitmasks implicit none @@ -1215,7 +1216,7 @@ subroutine get_occ_from_key(key,occ,Nint) end -subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) +subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,shortcut,sort_idx,Nint) use bitmasks implicit none BEGIN_DOC @@ -1233,29 +1234,55 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) integer, allocatable :: idx(:) double precision :: hij double precision, allocatable :: vt(:) - integer :: i,j,k,l, jj + integer :: i,j,k,l, jj,ii,sh integer :: i0, j0 + + integer,intent(in) :: shortcut(0:n+1), sort_idx(n) + integer :: tmp, warp(2,0:n+1), ni + + ASSERT (Nint > 0) ASSERT (Nint == N_int) ASSERT (n>0) PROVIDE ref_bitmask_energy !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,idx,jj,vt) & - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold) + !$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,warp,tmp,sh) & + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold,shortcut,sort_idx) allocate(idx(0:n), vt(n)) Vt = 0.d0 v_0 = 0.d0 - !$OMP DO SCHEDULE(guided) - do i=1,n - idx(0) = i - call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx) - do jj=1,idx(0) - j = idx(jj) - if ( dabs(u_0(j)) + dabs(u_0(i)) > davidson_threshold ) then - 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) - endif + !$OMP DO SCHEDULE(dynamic) + + + do sh=1,shortcut(0) + warp(1,0) = 0 + do ii=1,sh!shortcut(0) + tmp = 0 + do ni=1,Nint + tmp = popcnt(xor(keys_tmp(ni,1, shortcut(ii)), keys_tmp(ni,1,shortcut(sh)))) + end do + if(tmp <= 4) then + warp(1,0) = warp(1,0) + 1 + warp(1,warp(1,0)) = shortcut(ii) + warp(2,warp(1,0)) = shortcut(ii+1)-1 + end if + end do + + + do ii=shortcut(sh),shortcut(sh+1)-1 + idx(0) = ii + + call filter_connected_davidson_warp(keys_tmp,warp,keys_tmp(1,1,ii),Nint,ii-1,idx) + i = sort_idx(ii) + + do jj=1,idx(0) + j = sort_idx(idx(jj)) + if ( dabs(u_0(j)) + dabs(u_0(i)) > davidson_threshold ) then + call i_H_j(keys_tmp(1,1,idx(jj)),keys_tmp(1,1,ii),Nint,hij) + vt (i) = vt (i) + hij*u_0(j) + vt (j) = vt (j) + hij*u_0(i) + endif + enddo enddo enddo !$OMP END DO @@ -1267,7 +1294,75 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) deallocate(idx,vt) !$OMP END PARALLEL do i=1,n - v_0(i) += H_jj(i) * u_0(i) + v_0(i) += H_jj(i) * u_0(i) + enddo +end + + +subroutine H_u_0_org(v_0,u_0,H_jj,n,keys_tmp,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> + ! + ! n : number of determinants + ! + ! H_jj : array of + END_DOC + integer, intent(in) :: n,Nint + double precision, intent(out) :: v_0(n) + double precision, intent(in) :: u_0(n) + double precision, intent(in) :: H_jj(n) + integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) + integer, allocatable :: idx(:) + double precision :: hij + double precision, allocatable :: vt(:) + integer :: i,j,k,l, jj,ii,sh + integer :: i0, j0 + + + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (n>0) + PROVIDE ref_bitmask_energy + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii) & + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold) + allocate(idx(0:n), vt(n)) + Vt = 0.d0 + v_0 = 0.d0 + !$OMP DO SCHEDULE(guided) + + + + + + do ii=1,n + idx(0) = ii + i = ii + call filter_connected_davidson(keys_tmp,keys_tmp(1,1,ii),Nint,ii-1,idx) + + do jj=1,idx(0) + j = idx(jj) + if ( dabs(u_0(j)) + dabs(u_0(i)) > davidson_threshold ) then + call i_H_j(keys_tmp(1,1,idx(jj)),keys_tmp(1,1,ii),Nint,hij) + vt (i) = vt (i) + hij*u_0(j) + vt (j) = vt (j) + hij*u_0(i) + endif + enddo + enddo + + !$OMP END DO + !$OMP CRITICAL + do i=1,n + v_0(i) = v_0(i) + vt(i) + enddo + !$OMP END CRITICAL + deallocate(idx,vt) + !$OMP END PARALLEL + do i=1,n + v_0(i) += H_jj(i) * u_0(i) enddo end diff --git a/src/Ezfio_files/README.rst b/src/Ezfio_files/README.rst index f361c168..986d7216 100644 --- a/src/Ezfio_files/README.rst +++ b/src/Ezfio_files/README.rst @@ -203,99 +203,83 @@ Documentation Output file for Bitmask -`output_cisd `_ - Output file for CISD - - `output_cpu_time_0 `_ Initial CPU and wall times when printing in the output files -`output_determinants `_ +`output_determinants `_ Output file for Determinants -`output_electrons `_ +`output_electrons `_ Output file for Electrons -`output_ezfio_files `_ +`output_ezfio_files `_ Output file for Ezfio_files -`output_full_ci `_ - Output file for Full_CI - - -`output_generators_full `_ +`output_generators_full `_ Output file for Generators_full -`output_hartree_fock `_ +`output_hartree_fock `_ Output file for Hartree_Fock -`output_integrals_bielec `_ +`output_integrals_bielec `_ Output file for Integrals_Bielec -`output_integrals_monoelec `_ +`output_integrals_monoelec `_ Output file for Integrals_Monoelec -`output_mo_basis `_ +`output_mo_basis `_ Output file for MO_Basis -`output_moguess `_ +`output_moguess `_ Output file for MOGuess -`output_mrcc_cassd `_ +`output_mrcc_cassd `_ Output file for MRCC_CASSD -`output_mrcc_utils `_ +`output_mrcc_utils `_ Output file for MRCC_Utils -`output_myhartreefock `_ - Output file for MyHartreeFock - - -`output_nuclei `_ +`output_nuclei `_ Output file for Nuclei -`output_perturbation `_ +`output_perturbation `_ Output file for Perturbation -`output_properties `_ +`output_properties `_ Output file for Properties -`output_pseudo `_ +`output_pseudo `_ Output file for Pseudo -`output_psiref_cas `_ +`output_psiref_cas `_ Output file for Psiref_CAS -`output_psiref_utils `_ +`output_psiref_utils `_ Output file for Psiref_Utils -`output_selectors_full `_ +`output_selectors_full `_ Output file for Selectors_full -`output_singlerefmethod `_ - Output file for SingleRefMethod - - -`output_utils `_ +`output_utils `_ Output file for Utils