From e23dba89ba7c84f6b382b62fae16a5f62ce3e0fb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 29 Jul 2015 18:27:07 +0200 Subject: [PATCH] Accelerated MRCC --- plugins/MRCC_CASSD/.gitignore | 1 + plugins/MRCC_Utils_new/README.rst | 34 +++---- plugins/MRCC_Utils_new/mrcc_dress.irp.f | 28 ++--- plugins/Perturbation/perturbation.template.f | 4 +- plugins/Psiref_Utils/README.rst | 8 ++ plugins/Psiref_Utils/psi_ref_utils.irp.f | 101 +++++++++++++++++++ src/Determinants/README.rst | 26 ++--- src/Determinants/connected_to_ref.irp.f | 16 +-- src/Determinants/occ_pattern.irp.f | 2 +- src/Determinants/spindeterminants.irp.f | 24 ++--- src/Ezfio_files/README.rst | 4 +- 11 files changed, 179 insertions(+), 69 deletions(-) diff --git a/plugins/MRCC_CASSD/.gitignore b/plugins/MRCC_CASSD/.gitignore index 4f0460a1..de65845a 100644 --- a/plugins/MRCC_CASSD/.gitignore +++ b/plugins/MRCC_CASSD/.gitignore @@ -28,4 +28,5 @@ Utils ezfio_interface.irp.f irpf90.make irpf90_entities +mrcc_cassd tags \ No newline at end of file diff --git a/plugins/MRCC_Utils_new/README.rst b/plugins/MRCC_Utils_new/README.rst index c3c8a85f..5d2f243e 100644 --- a/plugins/MRCC_Utils_new/README.rst +++ b/plugins/MRCC_Utils_new/README.rst @@ -23,7 +23,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. by the `update_README.py` script. -`apply_excitation_operator `_ +`apply_excitation_operator `_ Undocumented @@ -108,27 +108,27 @@ Documentation .br excitation_operators(:,i) represents the holes and particles that link the ith connected determinant to det_ref - if :: + if :: excitation_operators(5,i) = 2 :: double excitation alpha excitation_operators(5,i) = -2 :: double excitation beta - !! excitation_operators(1,i) :: hole 1 - !! excitation_operators(2,i) :: particle 1 - !! excitation_operators(3,i) :: hole 2 - !! excitation_operators(4,i) :: particle 2 - else if :: + !! excitation_operators(1,i) :: hole 1 + !! excitation_operators(2,i) :: particle 1 + !! excitation_operators(3,i) :: hole 2 + !! excitation_operators(4,i) :: particle 2 + else if :: excitation_operators(5,i) = 1 :: single excitation alpha - !! excitation_operators(1,i) :: hole 1 - !! excitation_operators(2,i) :: particle 1 - else if :: + !! excitation_operators(1,i) :: hole 1 + !! excitation_operators(2,i) :: particle 1 + else if :: excitation_operators(5,i) = -1 :: single excitation beta - !! excitation_operators(3,i) :: hole 1 - !! excitation_operators(4,i) :: particle 1 - else if :: + !! excitation_operators(3,i) :: hole 1 + !! excitation_operators(4,i) :: particle 1 + else if :: !! excitation_operators(5,i) = 0 :: double excitation alpha/beta - !! excitation_operators(1,i) :: hole 1 alpha - !! excitation_operators(2,i) :: particle 1 alpha - !! excitation_operators(3,i) :: hole 2 beta - !! excitation_operators(4,i) :: particle 2 beta + !! excitation_operators(1,i) :: hole 1 alpha + !! excitation_operators(2,i) :: particle 1 alpha + !! excitation_operators(3,i) :: hole 2 beta + !! excitation_operators(4,i) :: particle 2 beta `h_matrix_dressed `_ diff --git a/plugins/MRCC_Utils_new/mrcc_dress.irp.f b/plugins/MRCC_Utils_new/mrcc_dress.irp.f index 5f287850..ee998995 100644 --- a/plugins/MRCC_Utils_new/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils_new/mrcc_dress.irp.f @@ -3,7 +3,7 @@ subroutine mrcc_dress(ndetref,ndetnonref,nstates,delta_ij_,delta_ii_) implicit none integer, intent(in) :: ndetref,nstates,ndetnonref double precision, intent(inout) :: delta_ii_(ndetref,nstates),delta_ij_(ndetref,ndetnonref,nstates) - integer :: i,j,k,l + integer :: i,j,k,l,m integer :: i_state integer :: N_connect_ref integer*2,allocatable :: excitation_operators(:,:) @@ -12,13 +12,14 @@ subroutine mrcc_dress(ndetref,ndetnonref,nstates,delta_ij_,delta_ii_) integer(bit_kind), allocatable :: key_test(:,:) integer, allocatable :: index_connected(:) integer :: i_hole,i_particle,ispin,i_ok,connected_to_ref,index_wf - integer, allocatable :: idx_vector(:), degree_vector(:) + integer, allocatable :: idx_vector(:) double precision :: phase_ij double precision :: dij,phase_la double precision :: hij,phase integer :: exc(0:2,2,2),degree logical :: is_in_wavefunction double precision, allocatable :: delta_ij_tmp(:,:,:), delta_ii_tmp(:,:) + logical, external :: is_in_psi_ref i_state = 1 allocate(excitation_operators(5,N_det_non_ref)) @@ -29,15 +30,14 @@ subroutine mrcc_dress(ndetref,ndetnonref,nstates,delta_ij_,delta_ii_) !$OMP SHARED(N_det_ref, N_det_non_ref, psi_ref, i_state, & !$OMP N_connect_ref,index_connected,psi_non_ref, & !$OMP excitation_operators,amplitudes_phase_less, & - !$OMP psi_non_ref_coef,N_int,lambda_mrcc,N_det, & + !$OMP psi_non_ref_coef,N_int,lambda_mrcc, & !$OMP delta_ii_,delta_ij_,psi_ref_coef,nstates, & - !$OMP mo_integrals_threshold) & + !$OMP mo_integrals_threshold,idx_non_ref_rev) & !$OMP PRIVATE(i,j,k,l,hil,phase_il,exc,degree,t_il, & - !$OMP key_test,i_ok,phase_la,hij,phase_ij, & - !$OMP dij,degree_vector,idx_vector,delta_ij_tmp, & + !$OMP key_test,i_ok,phase_la,hij,phase_ij,m, & + !$OMP dij,idx_vector,delta_ij_tmp, & !$OMP delta_ii_tmp,phase) allocate(idx_vector(0:N_det_non_ref)) - allocate(degree_vector(N_det_non_ref)) allocate(key_test(N_int,2)) allocate(delta_ij_tmp(size(delta_ij_,1),size(delta_ij_,2),nstates)) allocate(delta_ii_tmp(size(delta_ij_,1),nstates)) @@ -52,8 +52,9 @@ subroutine mrcc_dress(ndetref,ndetnonref,nstates,delta_ij_,delta_ii_) !$OMP END SINGLE !$OMP BARRIER - !$OMP DO SCHEDULE(guided) + !$OMP DO SCHEDULE(dynamic) do l = 1, N_det_non_ref +! print *, l, '/', N_det_non_ref double precision :: t_il,phase_il,hil call i_H_j_phase_out(psi_ref(1,1,i),psi_non_ref(1,1,l),N_int,hil,phase_il,exc,degree) t_il = hil * lambda_mrcc(i_state,l) @@ -75,7 +76,7 @@ subroutine mrcc_dress(ndetref,ndetnonref,nstates,delta_ij_,delta_ii_) if(i_ok.ne.1)cycle ! we check if such determinant is already in the wave function - if(is_in_wavefunction(key_test,N_int,N_det))cycle + if(is_in_wavefunction(key_test,N_int))cycle ! we get the phase for psi_non_ref(l) -> T_I->j |psi_non_ref(l)> call get_excitation(psi_non_ref(1,1,l),key_test,exc,degree,phase_la,N_int) @@ -90,12 +91,14 @@ subroutine mrcc_dress(ndetref,ndetnonref,nstates,delta_ij_,delta_ii_) endif ! we compute the interaction of such determinant with all the non_ref dets - call get_excitation_degree_vector(psi_non_ref,key_test,degree_vector,N_int,N_det_non_ref,idx_vector) + call filter_connected(psi_non_ref,key_test,N_int,N_det_non_ref,idx_vector) do k = 1, idx_vector(0) - call i_H_j_phase_out(key_test,psi_non_ref(1,1,idx_vector(k)),N_int,hij,phase,exc,degree) - delta_ij_tmp(i,idx_vector(k),i_state) += hij * dij + m = idx_vector(k) + call i_H_j_phase_out(key_test,psi_non_ref(1,1,m),N_int,hij,phase,exc,degree) + delta_ij_tmp(i,m,i_state) += hij * dij enddo + enddo @@ -117,7 +120,6 @@ subroutine mrcc_dress(ndetref,ndetnonref,nstates,delta_ij_,delta_ii_) deallocate(delta_ii_tmp,delta_ij_tmp) deallocate(idx_vector) deallocate(key_test) - deallocate(degree_vector) !$OMP END PARALLEL deallocate(excitation_operators) diff --git a/plugins/Perturbation/perturbation.template.f b/plugins/Perturbation/perturbation.template.f index 9f50b636..a5ab12e7 100644 --- a/plugins/Perturbation/perturbation.template.f +++ b/plugins/Perturbation/perturbation.template.f @@ -31,7 +31,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c cycle endif - if (is_in_wavefunction(buffer(1,1,i),Nint,N_det)) then + if (is_in_wavefunction(buffer(1,1,i),Nint)) then cycle endif @@ -82,7 +82,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ cycle endif - if (is_in_wavefunction(buffer(1,1,i),Nint,N_det)) then + if (is_in_wavefunction(buffer(1,1,i),Nint)) then cycle endif diff --git a/plugins/Psiref_Utils/README.rst b/plugins/Psiref_Utils/README.rst index 26c2e87d..d00ee877 100644 --- a/plugins/Psiref_Utils/README.rst +++ b/plugins/Psiref_Utils/README.rst @@ -13,6 +13,10 @@ Documentation .. Do not edit this section. It was auto-generated from the .. by the `update_README.py` script. +`get_index_in_psi_ref_sorted_bit `_ + Returns the index of the determinant in the ``psi_ref_sorted_bit`` array + + `h_matrix_ref `_ Undocumented @@ -30,6 +34,10 @@ Documentation idx_non_ref gives the indice of the determinant in psi_det. +`is_in_psi_ref `_ + True if the determinant ``det`` is in the wave function + + `n_det_non_ref `_ Set of determinants which are not part of the reference, defined from the application of the reference bitmask on the determinants. diff --git a/plugins/Psiref_Utils/psi_ref_utils.irp.f b/plugins/Psiref_Utils/psi_ref_utils.irp.f index c5e80199..21c58056 100644 --- a/plugins/Psiref_Utils/psi_ref_utils.irp.f +++ b/plugins/Psiref_Utils/psi_ref_utils.irp.f @@ -18,17 +18,20 @@ END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref, (N_int,2,psi_det_size) ] &BEGIN_PROVIDER [ double precision, psi_non_ref_coef, (psi_det_size,n_states) ] &BEGIN_PROVIDER [ integer, idx_non_ref, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, idx_non_ref_rev, (psi_det_size) ] &BEGIN_PROVIDER [ integer, N_det_non_ref ] implicit none BEGIN_DOC ! Set of determinants which are not part of the reference, defined from the application ! of the reference bitmask on the determinants. ! idx_non_ref gives the indice of the determinant in psi_det. + ! idx_non_ref_rev gives the reverse. END_DOC integer :: i_non_ref,j,k integer :: degree logical :: in_ref i_non_ref =0 + idx_non_ref_rev = 0 do k=1,N_det in_ref = .False. do j=1,N_det_ref @@ -49,6 +52,7 @@ END_PROVIDER psi_non_ref_coef(i_non_ref,j) = psi_coef(k,j) enddo idx_non_ref(i_non_ref) = k + idx_non_ref_rev(k) = i_non_ref endif enddo N_det_non_ref = i_non_ref @@ -119,5 +123,102 @@ END_PROVIDER END_PROVIDER +logical function is_in_psi_ref(key,Nint) + use bitmasks + implicit none + BEGIN_DOC +! True if the determinant ``det`` is in the wave function + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key(Nint,2) + integer, external :: get_index_in_psi_ref_sorted_bit + !DIR$ FORCEINLINE + is_in_psi_ref = get_index_in_psi_ref_sorted_bit(key,Nint) > 0 +end + +integer function get_index_in_psi_ref_sorted_bit(key,Nint) + use bitmasks + BEGIN_DOC +! Returns the index of the determinant in the ``psi_ref_sorted_bit`` array + END_DOC + implicit none + + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key(Nint,2) + + integer :: i, ibegin, iend, istep, l + integer*8 :: det_ref, det_search + integer*8, external :: det_search_key + logical :: in_wavefunction + + in_wavefunction = .False. + get_index_in_psi_ref_sorted_bit = 0 + ibegin = 1 + iend = N_det+1 + + !DIR$ FORCEINLINE + det_ref = det_search_key(key,Nint) + !DIR$ FORCEINLINE + det_search = det_search_key(psi_ref_sorted_bit(1,1,1),Nint) + + istep = ishft(iend-ibegin,-1) + i=ibegin+istep + do while (istep > 0) + !DIR$ FORCEINLINE + det_search = det_search_key(psi_ref_sorted_bit(1,1,i),Nint) + if ( det_search > det_ref ) then + iend = i + else if ( det_search == det_ref ) then + exit + else + ibegin = i + endif + istep = ishft(iend-ibegin,-1) + i = ibegin + istep + end do + + !DIR$ FORCEINLINE + do while (det_search_key(psi_ref_sorted_bit(1,1,i),Nint) == det_ref) + i = i-1 + if (i == 0) then + exit + endif + enddo + i += 1 + + if (i > N_det) then + return + endif + + !DIR$ FORCEINLINE + do while (det_search_key(psi_ref_sorted_bit(1,1,i),Nint) == det_ref) + if ( (key(1,1) /= psi_ref_sorted_bit(1,1,i)).or. & + (key(1,2) /= psi_ref_sorted_bit(1,2,i)) ) then + continue + else + in_wavefunction = .True. + !DIR$ IVDEP + !DIR$ LOOP COUNT MIN(3) + do l=2,Nint + if ( (key(l,1) /= psi_ref_sorted_bit(l,1,i)).or. & + (key(l,2) /= psi_ref_sorted_bit(l,2,i)) ) then + in_wavefunction = .False. + endif + enddo + if (in_wavefunction) then + get_index_in_psi_ref_sorted_bit = i +! exit + return + endif + endif + i += 1 + if (i > N_det) then +! exit + return + endif + + enddo + +end diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index 1e5751ef..c7e455bf 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -43,7 +43,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. by the `update_README.py` script. -`a_operator `_ +`a_operator `_ Needed for diag_H_mat_elem @@ -55,7 +55,7 @@ Documentation Max and min values of the coefficients -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem @@ -141,7 +141,7 @@ Documentation After calling this subroutine, N_det, psi_det and psi_coef need to be touched -`create_wf_of_psi_svd_matrix `_ +`create_wf_of_psi_svd_matrix `_ Matrix of wf coefficients. Outer product of alpha and beta determinants @@ -217,7 +217,7 @@ Documentation det_coef -`det_connections `_ +`det_connections `_ Build connection proxy between determinants @@ -245,7 +245,7 @@ Documentation Diagonalization algorithm (Davidson or Lapack) -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes @@ -347,7 +347,7 @@ Documentation Determinants are taken from the psi_det_sorted_ab array -`generate_all_alpha_beta_det_products `_ +`generate_all_alpha_beta_det_products `_ Create a wave function from all possible alpha x beta determinants @@ -383,7 +383,7 @@ Documentation 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 @@ -417,7 +417,7 @@ Documentation Undocumented -`h_u_0 `_ +`h_u_0 `_ Computes v_0 = H|u_0> .br n : number of determinants @@ -516,7 +516,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 @@ -774,19 +774,19 @@ Documentation psi_occ_pattern(:,2,j) = jth occ_pattern of the wave function : represent all the double occupation -`psi_svd_alpha `_ +`psi_svd_alpha `_ SVD wave function -`psi_svd_beta `_ +`psi_svd_beta `_ SVD wave function -`psi_svd_coefs `_ +`psi_svd_coefs `_ SVD wave function -`psi_svd_matrix `_ +`psi_svd_matrix `_ Matrix of wf coefficients. Outer product of alpha and beta determinants diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index 3f4b62d2..49a3604a 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -33,13 +33,13 @@ end -logical function is_in_wavefunction(key,Nint,Ndet) +logical function is_in_wavefunction(key,Nint) use bitmasks implicit none BEGIN_DOC ! True if the determinant ``det`` is in the wave function END_DOC - integer, intent(in) :: Nint, Ndet + integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key(Nint,2) integer, external :: get_index_in_psi_det_sorted_bit @@ -60,9 +60,9 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) integer :: i, ibegin, iend, istep, l integer*8 :: det_ref, det_search integer*8, external :: det_search_key - logical :: is_in_wavefunction + logical :: in_wavefunction - is_in_wavefunction = .False. + in_wavefunction = .False. get_index_in_psi_det_sorted_bit = 0 ibegin = 1 iend = N_det+1 @@ -107,16 +107,16 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) (key(1,2) /= psi_det_sorted_bit(1,2,i)) ) then continue else - is_in_wavefunction = .True. + in_wavefunction = .True. !DIR$ IVDEP !DIR$ LOOP COUNT MIN(3) do l=2,Nint if ( (key(l,1) /= psi_det_sorted_bit(l,1,i)).or. & (key(l,2) /= psi_det_sorted_bit(l,2,i)) ) then - is_in_wavefunction = .False. + in_wavefunction = .False. endif enddo - if (is_in_wavefunction) then + if (in_wavefunction) then get_index_in_psi_det_sorted_bit = i ! exit return @@ -131,7 +131,7 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) enddo ! DEBUG is_in_wf -! if (is_in_wavefunction) then +! if (in_wavefunction) then ! degree = 1 ! do i=1,N_det ! integer :: degree diff --git a/src/Determinants/occ_pattern.irp.f b/src/Determinants/occ_pattern.irp.f index 8027cf60..aa059870 100644 --- a/src/Determinants/occ_pattern.irp.f +++ b/src/Determinants/occ_pattern.irp.f @@ -292,7 +292,7 @@ subroutine make_s2_eigenfunction endif call occ_pattern_to_dets(psi_occ_pattern(1,1,i),d,s,elec_alpha_num,N_int) do j=1,s - if (.not. is_in_wavefunction( d(1,1,j), N_int, N_det)) then + if (.not. is_in_wavefunction(d(1,1,j), N_int) ) then N_det_new += 1 do k=1,N_int det_buffer(k,1,N_det_new) = d(k,1,j) diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index ad9d9960..5279dc2c 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -151,9 +151,9 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint) integer :: i, ibegin, iend, istep, l integer*8 :: det_ref, det_search integer*8, external :: spin_det_search_key - logical :: is_in_wavefunction + logical :: in_wavefunction - is_in_wavefunction = .False. + in_wavefunction = .False. get_index_in_psi_det_alpha_unique = 0 ibegin = 1 iend = N_det_alpha_unique + 1 @@ -198,15 +198,15 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint) if (key(1) /= psi_det_alpha_unique(1,i)) then continue else - is_in_wavefunction = .True. + in_wavefunction = .True. !DIR$ IVDEP !DIR$ LOOP COUNT MIN(3) do l=2,Nint if (key(l) /= psi_det_alpha_unique(l,i)) then - is_in_wavefunction = .False. + in_wavefunction = .False. endif enddo - if (is_in_wavefunction) then + if (in_wavefunction) then get_index_in_psi_det_alpha_unique = i return endif @@ -233,9 +233,9 @@ integer function get_index_in_psi_det_beta_unique(key,Nint) integer :: i, ibegin, iend, istep, l integer*8 :: det_ref, det_search integer*8, external :: spin_det_search_key - logical :: is_in_wavefunction + logical :: in_wavefunction - is_in_wavefunction = .False. + in_wavefunction = .False. get_index_in_psi_det_beta_unique = 0 ibegin = 1 iend = N_det_beta_unique + 1 @@ -279,15 +279,15 @@ integer function get_index_in_psi_det_beta_unique(key,Nint) if (key(1) /= psi_det_beta_unique(1,i)) then continue else - is_in_wavefunction = .True. + in_wavefunction = .True. !DIR$ IVDEP !DIR$ LOOP COUNT MIN(3) do l=2,Nint if (key(l) /= psi_det_beta_unique(l,i)) then - is_in_wavefunction = .False. + in_wavefunction = .False. endif enddo - if (is_in_wavefunction) then + if (in_wavefunction) then get_index_in_psi_det_beta_unique = i return endif @@ -369,7 +369,6 @@ BEGIN_PROVIDER [ double precision, psi_svd_matrix_values, (N_det,N_states) ] integer(bit_kind) :: tmp_det(N_int,2) integer :: idx integer, external :: get_index_in_psi_det_sorted_bit - logical, external :: is_in_wavefunction PROVIDE psi_coef_sorted_bit @@ -423,7 +422,6 @@ subroutine create_wf_of_psi_svd_matrix integer(bit_kind) :: tmp_det(N_int,2) integer :: idx integer, external :: get_index_in_psi_det_sorted_bit - logical, external :: is_in_wavefunction double precision :: norm(N_states) call generate_all_alpha_beta_det_products @@ -494,7 +492,7 @@ subroutine generate_all_alpha_beta_det_products tmp_det(k,1,l) = psi_det_alpha_unique(k,i) tmp_det(k,2,l) = psi_det_beta_unique (k,j) enddo - if (.not.is_in_wavefunction(tmp_det(1,1,l),N_int,N_det)) then + if (.not.is_in_wavefunction(tmp_det(1,1,l),N_int)) then l = l+1 endif enddo diff --git a/src/Ezfio_files/README.rst b/src/Ezfio_files/README.rst index 0d825518..c97e6268 100644 --- a/src/Ezfio_files/README.rst +++ b/src/Ezfio_files/README.rst @@ -106,8 +106,8 @@ Documentation Output file for MRCC_CASSD -`output_mrcc_utils `_ - Output file for MRCC_Utils +`output_mrcc_utils_new `_ + Output file for MRCC_Utils_new `output_nuclei `_