From d6ec67e8408ab2408dafd2db1d167b4503647184 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 6 May 2015 17:01:45 +0200 Subject: [PATCH] Fixed Full CAS-SD --- src/CAS_SD/cas_sd.irp.f | 14 ----- src/Determinants/H_apply.irp.f | 95 ++++++++++++++++++++++---------- src/Determinants/README.rst | 40 +++++++------- src/Determinants/s2.irp.f | 3 +- src/MRCC/README.rst | 4 +- src/MonoInts/pot_ao_ints.irp.f | 2 +- src/Perturbation/selection.irp.f | 2 +- 7 files changed, 92 insertions(+), 68 deletions(-) diff --git a/src/CAS_SD/cas_sd.irp.f b/src/CAS_SD/cas_sd.irp.f index 144eec88..0b34e7bd 100644 --- a/src/CAS_SD/cas_sd.irp.f +++ b/src/CAS_SD/cas_sd.irp.f @@ -69,18 +69,4 @@ program full_ci print *, '' enddo print *, 'Max excitation degree in the CAS :', exc_max - do i=1,N_det - in_cas = .False. - degree_min = 1000 - do k=1,N_det_generators - call get_excitation_degree(psi_det_generators(1,1,k),psi_det(1,1,i),degree,N_int) - degree_min = min(degree_min,degree) - enddo - if (degree_min > 2) then - print *, 'Error : This is not a CAS-SD : ' - print *, 'Excited determinant:', degree_min - call debug_det(psi_det(1,1,k),N_int) - stop - endif - enddo end diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f index d230c765..dab4390a 100644 --- a/src/Determinants/H_apply.irp.f +++ b/src/Determinants/H_apply.irp.f @@ -53,13 +53,17 @@ subroutine resize_H_apply_buffer(new_size,iproc) double precision, pointer :: buffer_e2(:,:) integer :: i,j,k integer :: Ndet + + BEGIN_DOC +! Resizes the H_apply buffer of proc iproc. The buffer lock should +! be set before calling this function. + END_DOC PROVIDE H_apply_buffer_allocated ASSERT (new_size > 0) ASSERT (iproc >= 0) ASSERT (iproc < nproc) - call omp_set_lock(H_apply_buffer_lock(1,iproc)) allocate ( buffer_det(N_int,2,new_size), & buffer_coef(new_size,N_states), & buffer_e2(new_size,N_states) ) @@ -93,7 +97,6 @@ subroutine resize_H_apply_buffer(new_size,iproc) H_apply_buffer(iproc)%sze = new_size H_apply_buffer(iproc)%N_det = min(new_size,H_apply_buffer(iproc)%N_det) - call omp_unset_lock(H_apply_buffer_lock(1,iproc)) end @@ -101,8 +104,7 @@ subroutine copy_H_apply_buffer_to_wf use omp_lib implicit none BEGIN_DOC -! Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det -! after calling this function. +! Copies the H_apply buffer to psi_coef. ! After calling this subroutine, N_det, psi_det and psi_coef need to be touched END_DOC integer(bit_kind), allocatable :: buffer_det(:,:,:) @@ -181,42 +183,77 @@ subroutine copy_H_apply_buffer_to_wf call normalize(psi_coef,N_det) SOFT_TOUCH N_det psi_det psi_coef - call debug_unicity_of_determinants + logical :: found_duplicates + call remove_duplicates_in_psi_det(found_duplicates) end -subroutine debug_unicity_of_determinants +subroutine remove_duplicates_in_psi_det(found_duplicates) implicit none + logical, intent(out) :: found_duplicates BEGIN_DOC -! This subroutine checks that there are no repetitions in the wave function +! Removes duplicate determinants in the wave function. END_DOC - logical :: same, failed - integer :: i,k - print *, "======= DEBUG UNICITY =========" - failed = .False. - do i=2,N_det - same = .True. - do k=1,N_int - if ( psi_det_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,i-1) ) then - same = .False. - exit + integer :: i,j,k + integer(bit_kind), allocatable :: bit_tmp(:) + logical,allocatable :: duplicate(:) + + allocate (duplicate(N_det), bit_tmp(N_det)) + + do i=1,N_det + integer, external :: det_search_key + !$DIR FORCEINLINE + bit_tmp(i) = det_search_key(psi_det_sorted_bit(1,1,i),N_int) + duplicate(i) = .False. + enddo + + do i=1,N_det-1 + if (duplicate(i)) then + cycle + endif + j = i+1 + do while (bit_tmp(j)==bit_tmp(i)) + if (duplicate(j)) then + j += 1 + cycle endif - if ( psi_det_sorted_bit(k,2,i) /= psi_det_sorted_bit(k,2,i-1) ) then - same = .False. + duplicate(j) = .True. + do k=1,N_int + if ( (psi_det_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,j) ) & + .or. (psi_det_sorted_bit(k,2,i) /= psi_det_sorted_bit(k,2,j) ) ) then + duplicate(j) = .False. + exit + endif + enddo + j += 1 + if (j > N_det) then exit endif enddo - if (same) then - failed = .True. - call debug_det(psi_det_sorted_bit(1,1,i)) + enddo + + found_duplicates = .False. + do i=1,N_det + if (duplicate(i)) then + found_duplicates = .True. + exit endif enddo - if (failed) then - print *, '======= Determinants not unique : Failed ! =========' - stop - else - print *, '======= Determinants are unique : OK ! =========' + if (found_duplicates) then + call write_bool(output_determinants,found_duplicates,'Found duplicate determinants') + k=0 + do i=1,N_det + if (.not.duplicate(i)) then + k += 1 + psi_det(:,:,k) = psi_det_sorted_bit (:,:,i) + psi_coef(k,:) = psi_coef_sorted_bit(i,:) + endif + enddo + print *, N_det,k + N_det = k + TOUCH N_det psi_det psi_coef endif + deallocate (duplicate,bit_tmp) end subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) @@ -231,11 +268,11 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) integer :: i,j,k integer :: new_size PROVIDE H_apply_buffer_allocated + call omp_set_lock(H_apply_buffer_lock(1,iproc)) new_size = H_apply_buffer(iproc)%N_det + n_selected if (new_size > H_apply_buffer(iproc)%sze) then call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc) endif - call omp_set_lock(H_apply_buffer_lock(1,iproc)) do i=1,H_apply_buffer(iproc)%N_det ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) @@ -250,7 +287,7 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) enddo do j=1,N_states do i=1,N_selected - H_apply_buffer(iproc)%coef(i,j) = 0.d0 + H_apply_buffer(iproc)%coef(i+H_apply_buffer(iproc)%N_det,j) = 0.d0 enddo enddo H_apply_buffer(iproc)%N_det = new_size diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index bf1e8c13..4b8c0472 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -187,10 +187,10 @@ Documentation `det_svd `_ Computes the SVD of the Alpha x Beta determinant coefficient matrix -`filter_3_highest_electrons `_ +`filter_3_highest_electrons `_ Returns a determinant with only the 3 highest electrons -`int_of_3_highest_electrons `_ +`int_of_3_highest_electrons `_ Returns an integer*8 as : .br |_<--- 21 bits ---><--- 21 bits ---><--- 21 bits --->| @@ -207,26 +207,26 @@ Documentation `n_det `_ Number of determinants in the wave function -`psi_average_norm_contrib `_ +`psi_average_norm_contrib `_ Contribution of determinants to the state-averaged density -`psi_average_norm_contrib_sorted `_ +`psi_average_norm_contrib_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef `_ +`psi_coef `_ The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_coef_sorted `_ +`psi_coef_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef_sorted_ab `_ +`psi_coef_sorted_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`psi_coef_sorted_bit `_ +`psi_coef_sorted_bit `_ Determinants on which we apply for perturbation. They are sorted by determinants interpreted as integers. Useful to accelerate the search of a random determinant in the wave @@ -239,46 +239,46 @@ Documentation `psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_det_sorted `_ +`psi_det_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_det_sorted_ab `_ +`psi_det_sorted_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`psi_det_sorted_bit `_ +`psi_det_sorted_bit `_ Determinants on which we apply for perturbation. They are sorted by determinants interpreted as integers. Useful to accelerate the search of a random determinant in the wave function. -`psi_det_sorted_next_ab `_ +`psi_det_sorted_next_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`read_dets `_ +`read_dets `_ Reads the determinants from the EZFIO file -`save_wavefunction `_ +`save_wavefunction `_ Save the wave function into the EZFIO file -`save_wavefunction_general `_ +`save_wavefunction_general `_ Save the wave function into the EZFIO file -`save_wavefunction_unsorted `_ +`save_wavefunction_unsorted `_ Save the wave function into the EZFIO file -`sort_dets_by_3_highest_electrons `_ +`sort_dets_by_3_highest_electrons `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`sort_dets_by_det_search_key `_ +`sort_dets_by_det_search_key `_ Determinants are sorted are sorted according to their det_search_key. Useful to accelerate the search of a random determinant in the wave function. @@ -316,7 +316,7 @@ Documentation `diag_algorithm `_ Diagonalization algorithm (Davidson or Lapack) -`diagonalize_ci `_ +`diagonalize_ci `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix @@ -345,7 +345,7 @@ Documentation `ci_electronic_energy_mono `_ Eigenvectors/values of the CI matrix -`diagonalize_ci_mono `_ +`diagonalize_ci_mono `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index 376528d7..760893e0 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -88,7 +88,7 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) double precision, intent(out) :: s2 integer :: i,j,l double precision :: s2_tmp - s2 = S_z2_Sz + s2 = 0.d0 !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP PRIVATE(i,j,s2_tmp) SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int) & !$OMP REDUCTION(+:s2) SCHEDULE(dynamic) @@ -101,5 +101,6 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) enddo enddo !$OMP END PARALLEL DO + s2 += S_z2_Sz end diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst index 38137667..5bb2009c 100644 --- a/src/MRCC/README.rst +++ b/src/MRCC/README.rst @@ -51,7 +51,7 @@ Documentation `ci_electronic_energy_dressed `_ Eigenvectors/values of the CI matrix -`ci_energy_dressed `_ +`ci_energy_dressed `_ N_states lowest eigenvalues of the dressed CI matrix `delta_ij `_ @@ -60,7 +60,7 @@ Documentation `delta_ij_non_cas `_ Dressing matrix in SD basis -`diagonalize_ci_dressed `_ +`diagonalize_ci_dressed `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index c93301e5..a9e72e57 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -10,7 +10,7 @@ integer :: i,j,k,l,n_pt_in,m double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - if (do_pseudo.eqv..TRUE.) then + if (do_pseudo) then ao_nucl_elec_integral = ao_nucl_elec_integral_pseudo else ao_nucl_elec_integral = 0.d0 diff --git a/src/Perturbation/selection.irp.f b/src/Perturbation/selection.irp.f index 77313888..84cc59ae 100644 --- a/src/Perturbation/selection.irp.f +++ b/src/Perturbation/selection.irp.f @@ -19,6 +19,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c ASSERT (Nint > 0) ASSERT (N_int == N_int) ASSERT (N_selected >= 0) + call omp_set_lock(H_apply_buffer_lock(1,iproc)) smax = selection_criterion smin = selection_criterion_min new_size = H_apply_buffer(iproc)%N_det + n_selected @@ -26,7 +27,6 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c if (new_size > h_apply_buffer(iproc)%sze) then call resize_h_apply_buffer(max(h_apply_buffer(iproc)%sze*2,new_size),iproc) endif - call omp_set_lock(H_apply_buffer_lock(1,iproc)) do i=1,H_apply_buffer(iproc)%N_det ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)