diff --git a/scripts/check_dependencies.sh b/scripts/check_dependencies.sh index 943c5211..4e7e4bf1 100755 --- a/scripts/check_dependencies.sh +++ b/scripts/check_dependencies.sh @@ -51,7 +51,7 @@ DEPS=$(unique_list $DEPS_LONG) if [[ ! "$COMMAND_LINE" == "$DEPS" ]] then - DEPS=$(check_dependencies.sh $DEPS) + DEPS=$(${QPACKAGE_ROOT}/scripts/check_dependencies.sh ${DEPS}) fi echo $DEPS diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 74c35116..f40833d7 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -10,7 +10,8 @@ subroutine parameters initialization declarations -keys_work +keys_work_locked +keys_work_unlocked finalization """.split() @@ -24,7 +25,10 @@ class H_apply(object): self.openmp = openmp if openmp: s["subroutine"] += "_OpenMP" + + self.selection_pt2 = None self.perturbation = None + #s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) & s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, & @@ -34,7 +38,7 @@ class H_apply(object): !$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,& !$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, & !$OMP N_elec_in_key_hole_2,ia_ja_pairs) & - !$OMP SHARED(key_in,N_int,elec_num_tab, & + !$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, & !$OMP hole_1, particl_1, hole_2, particl_2, & !$OMP lck,thresh,elec_alpha_num)""" s["omp_init_lock"] = "call omp_init_lock(lck)" @@ -72,6 +76,8 @@ class H_apply(object): return buffer def set_perturbation(self,pert): + if self.perturbation is not None: + raise self.perturbation = pert if pert is not None: self.data["parameters"] = ",sum_e_2_pert_in,sum_norm_pert_in,sum_H_pert_diag_in,N_st,Nint" @@ -83,16 +89,17 @@ class H_apply(object): double precision :: sum_e_2_pert(N_st) double precision :: sum_norm_pert(N_st) double precision :: sum_H_pert_diag + double precision :: e_2_pert_buffer(N_st,size_max) + double precision :: coef_pert_buffer(N_st,size_max) """ self.data["size_max"] = "256" self.data["initialization"] = """ - E_ref = diag_H_mat_elem(key_in,N_int) sum_e_2_pert = sum_e_2_pert_in sum_norm_pert = sum_norm_pert_in sum_H_pert_diag = sum_H_pert_diag_in """ - self.data["keys_work"] += """ - call perturb_buffer_%s(keys_out,key_idx,sum_e_2_pert, & + self.data["keys_work_unlocked"] += """ + call perturb_buffer_%s(keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & sum_norm_pert,sum_H_pert_diag,N_st,Nint) """%(pert,) self.data["finalization"] = """ @@ -101,10 +108,29 @@ class H_apply(object): sum_H_pert_diag_in = sum_H_pert_diag """ if self.openmp: + self.data["omp_set_lock"] = "" + self.data["omp_unset_lock"] = "" self.data["omp_test_lock"] = ".False." self.data["omp_parallel"] += """& - !$OMP SHARED(N_st) & + !$OMP SHARED(N_st,Nint) PRIVATE(e_2_pert_buffer,coef_pert_buffer) & !$OMP REDUCTION(+:sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)""" + def set_selection_pt2(self,pert): + if self.selection_pt2 is not None: + raise + self.set_perturbation(pert) + self.selection_pt2 = pert + if pert is not None: + self.data["size_max"] = str(1024*128) + self.data["keys_work_unlocked"] = """ + e_2_pert_buffer = 0.d0 + coef_pert_buffer = 0.d0 + """ + self.data["keys_work_unlocked"] + self.data["keys_work_locked"] = """ + call fill_H_apply_buffer_selection(key_idx,keys_out,e_2_pert_buffer,coef_pert_buffer,N_st,N_int) + """ + self.data["omp_test_lock"] = "omp_test_lock(lck)" + self.data["omp_set_lock"] = "call omp_set_lock(lck)" + self.data["omp_unset_lock"] = "call omp_unset_lock(lck)" diff --git a/src/Dets/H_apply.irp.f b/src/Dets/H_apply.irp.f index 40fe4b1b..974e123a 100644 --- a/src/Dets/H_apply.irp.f +++ b/src/Dets/H_apply.irp.f @@ -33,11 +33,12 @@ subroutine resize_H_apply_buffer_det(new_size) integer, intent(in) :: new_size integer(bit_kind), allocatable :: buffer_det(:,:,:) double precision, allocatable :: buffer_coef(:,:) + double precision, allocatable :: buffer_e2(:,:) integer :: i,j,k integer :: Ndet ASSERT (new_size > 0) - allocate ( buffer_det(N_int,2,new_size), buffer_coef(new_size,N_states) ) + allocate ( buffer_det(N_int,2,new_size), buffer_coef(new_size,N_states), buffer_e2(new_size,N_states) ) do i=1,min(new_size,H_apply_buffer_N_det) do k=1,N_int @@ -48,9 +49,10 @@ subroutine resize_H_apply_buffer_det(new_size) ASSERT (sum(popcnt(H_apply_buffer_det(:,2,i))) == elec_beta_num ) enddo do k=1,N_states - do i=1,min(new_size,H_apply_buffer_N_det) - buffer_coef(i,k) = H_apply_buffer_coef(i,k) - enddo + do i=1,min(new_size,H_apply_buffer_N_det) + buffer_coef(i,k) = H_apply_buffer_coef(i,k) + buffer_e2(i,k) = H_apply_buffer_e2(i,k) + enddo enddo H_apply_buffer_size = new_size @@ -70,20 +72,23 @@ subroutine resize_H_apply_buffer_det(new_size) do k=1,N_states do i=1,H_apply_buffer_N_det H_apply_buffer_coef(i,k) = buffer_coef(i,k) + H_apply_buffer_e2(i,k) = buffer_e2(i,k) enddo enddo - deallocate (buffer_det, buffer_coef) - SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det + deallocate (buffer_det, buffer_coef, buffer_e2) + SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det H_apply_buffer_e2 end BEGIN_PROVIDER [ integer(bit_kind), H_apply_buffer_det,(N_int,2,H_apply_buffer_size) ] &BEGIN_PROVIDER [ double precision, H_apply_buffer_coef,(H_apply_buffer_size,N_states) ] +&BEGIN_PROVIDER [ double precision, H_apply_buffer_e2,(H_apply_buffer_size,N_states) ] &BEGIN_PROVIDER [ integer, H_apply_buffer_N_det ] implicit none BEGIN_DOC - ! Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines. + ! Buffer of determinants/coefficients/perturbative energy for H_apply. + ! Uninitialized. Filled by H_apply subroutines. END_DOC H_apply_buffer_N_det = 0 @@ -148,8 +153,9 @@ subroutine copy_H_apply_buffer_to_wf psi_coef(i+N_det_old,k) = H_apply_buffer_coef(i,k) enddo enddo + H_apply_buffer_N_det = 0 - SOFT_TOUCH psi_det psi_coef + SOFT_TOUCH psi_det psi_coef H_apply_buffer_N_det H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_e2 end diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 6a7efa8d..14fea836 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -7,6 +7,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame ! particles. ! Assume N_int is already provided. END_DOC + integer,parameter :: size_max = $size_max $declarations integer(omp_lock_kind) :: lck integer(bit_kind),intent(in) :: key_in(N_int,2) @@ -26,13 +27,12 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) - integer,parameter :: size_max = $size_max double precision :: hij_elec, mo_bielec_integral, thresh integer, allocatable :: ia_ja_pairs(:,:,:) double precision :: diag_H_mat_elem - PROVIDE mo_integrals_map ref_bitmask_energy - PROVIDE mo_bielec_integrals_in_map + PROVIDE mo_integrals_map ref_bitmask_energy key_pattern_not_in_ref + PROVIDE mo_bielec_integrals_in_map reference_energy psi_ref_coef psi_ref $set_i_H_j_threshold @@ -156,8 +156,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame hij_tab(key_idx) = hij_elec ASSERT (key_idx <= size_max) if (key_idx == size_max) then + $keys_work_unlocked $omp_set_lock - $keys_work + $keys_work_locked $omp_unset_lock key_idx = 0 endif @@ -165,7 +166,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame enddo if (key_idx > ishft(size_max,-5)) then if ($omp_test_lock) then - $keys_work + $keys_work_unlocked + $keys_work_locked $omp_unset_lock key_idx = 0 endif @@ -204,8 +206,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame hij_tab(key_idx) = hij_elec ASSERT (key_idx <= size_max) if (key_idx == size_max) then + $keys_work_unlocked $omp_set_lock - $keys_work + $keys_work_locked $omp_unset_lock key_idx = 0 endif @@ -213,7 +216,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame enddo if (key_idx > ishft(size_max,-5)) then if ($omp_test_lock) then - $keys_work + $keys_work_locked + $keys_work_unlocked $omp_unset_lock key_idx = 0 endif @@ -222,8 +226,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame enddo ! ii $omp_enddo enddo ! ispin + $keys_work_unlocked $omp_set_lock - $keys_work + $keys_work_locked $omp_unset_lock deallocate (keys_out,hij_tab,ia_ja_pairs) $omp_end_parallel @@ -241,6 +246,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) ! particles. ! Assume N_int is already provided. END_DOC + integer,parameter :: size_max = $size_max $declarations integer(omp_lock_kind) :: lck integer(bit_kind),intent(in) :: key_in(N_int,2) @@ -260,13 +266,12 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) - integer,parameter :: size_max = $size_max double precision :: hij_elec, thresh integer, allocatable :: ia_ja_pairs(:,:,:) double precision :: diag_H_mat_elem - PROVIDE mo_integrals_map ref_bitmask_energy - PROVIDE mo_bielec_integrals_in_map + PROVIDE mo_integrals_map ref_bitmask_energy key_pattern_not_in_ref + PROVIDE mo_bielec_integrals_in_map reference_energy psi_ref_coef psi_ref $set_i_H_j_threshold @@ -311,7 +316,6 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) integer(bit_kind) :: test(N_int,2) double precision :: accu accu = 0.d0 - hij_elec = 0.d0 do ispin=1,2 other_spin = iand(ispin,1)+1 $omp_do @@ -325,33 +329,33 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) k_a = ishft(j_a-1,-bit_kind_shift)+1 l_a = j_a-ishft(k_a-1,bit_kind_shift)-1 hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) - call i_H_j(hole,key_in,N_int,hij_elec) - if(dabs(hij_elec) .ge. thresh)then - key_idx += 1 - do k=1,N_int - keys_out(k,1,key_idx) = hole(k,1) - keys_out(k,2,key_idx) = hole(k,2) - enddo - hij_tab(key_idx) = hij_elec - if (key_idx > ishft(size_max,-5)) then - if ($omp_test_lock) then - $keys_work - $omp_unset_lock - key_idx = 0 - endif - endif - if (key_idx == size_max) then - $omp_set_lock - $keys_work + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = hole(k,1) + keys_out(k,2,key_idx) = hole(k,2) + enddo + hij_tab(key_idx) = hij_elec + if (key_idx > ishft(size_max,-5)) then + if ($omp_test_lock) then + $keys_work_unlocked + $keys_work_locked $omp_unset_lock key_idx = 0 endif endif + if (key_idx == size_max) then + $keys_work_unlocked + $omp_set_lock + $keys_work_locked + $omp_unset_lock + key_idx = 0 + endif enddo ! ii $omp_enddo enddo ! ispin + $keys_work_unlocked $omp_set_lock - $keys_work + $keys_work_locked $omp_unset_lock deallocate (keys_out,hij_tab,ia_ja_pairs) $omp_end_parallel diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f index 81e4e12a..1fdaacf5 100644 --- a/src/Dets/connected_to_ref.irp.f +++ b/src/Dets/connected_to_ref.irp.f @@ -1,254 +1,256 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh) - use bitmasks - implicit none - integer, intent(in) :: Nint, N_past_in, Ndet - integer(bit_kind), intent(in) :: keys(ishft(Nint,-1),2,Ndet) - integer(bit_kind), intent(in) :: key(ishft(Nint,-1),2) - double precision, intent(in) :: thresh - - integer :: N_past - integer :: i, l - integer :: degree_x2 - logical :: det_is_not_or_may_be_in_ref, t - double precision :: hij_elec - - ! output : 0 : not connected - ! i : connected to determinant i of the past - ! -i : is the ith determinant of the refernce wf keys - - ASSERT (Nint == N_int) - - connected_to_ref = 0 - N_past = max(1,N_past_in) - if (Nint == 1) then - - do i=N_past-1,1,-1 - degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & - popcnt(xor( key(1,2), keys(1,2,i))) - if(degree_x2 == 0)then + use bitmasks + implicit none + integer, intent(in) :: Nint, N_past_in, Ndet + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: thresh + + integer :: N_past + integer :: i, l + integer :: degree_x2 + logical :: det_is_not_or_may_be_in_ref, t + double precision :: hij_elec + + ! output : 0 : not connected + ! i : connected to determinant i of the past + ! -i : is the ith determinant of the refernce wf keys + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + connected_to_ref = 0 + N_past = max(1,N_past_in) + if (Nint == 1) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + if(degree_x2 == 0)then + connected_to_ref = -i + return + endif + if (degree_x2 > 5) then + cycle + else + call i_H_j(keys(1,1,i),key,N_int,hij_elec) + if(dabs(hij_elec).lt.thresh)cycle + connected_to_ref = i + return + endif + enddo + + !DIR$ FORCEINLINE + t = det_is_not_or_may_be_in_ref(key,N_int) + if ( t ) then + return + endif + + do i=N_past,Ndet + if ( (key(1,1) /= keys(1,1,i)).or. & + (key(1,2) /= keys(1,2,i)) ) then + cycle + endif connected_to_ref = -i return - endif - if (degree_x2 > 5) then - cycle - else - call i_H_j(keys(1,1,i),key,N_int,hij_elec) - if(dabs(hij_elec).lt.thresh)cycle - connected_to_ref = i - return - endif - enddo - - !DIR$ FORCEINLINE - t = det_is_not_or_may_be_in_ref(key,N_int) - if ( t ) then - return - endif - - do i=N_past,Ndet - if ( (key(1,1) /= keys(1,1,i)).or. & - (key(1,2) /= keys(1,2,i)) ) then - cycle - endif - connected_to_ref = -i - return - enddo - return - - else if (Nint==2) then - - do i=N_past-1,1,-1 - degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & - popcnt(xor( key(1,2), keys(1,2,i))) + & - popcnt(xor( key(2,1), keys(2,1,i))) + & - popcnt(xor( key(2,2), keys(2,2,i))) - if(degree_x2 == 0)then + enddo + return + + else if (Nint==2) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + if(degree_x2 == 0)then + connected_to_ref = -i + return + endif + if (degree_x2 > 5) then + cycle + else + call i_H_j(keys(1,1,i),key,N_int,hij_elec) + if(dabs(hij_elec).lt.thresh)cycle + connected_to_ref = i + return + endif + enddo + + !DIR$ FORCEINLINE + t = det_is_not_or_may_be_in_ref(key,N_int) + if ( t ) then + return + endif + + !DIR$ LOOP COUNT (1000) + do i=N_past+1,Ndet + if ( (key(1,1) /= keys(1,1,i)).or. & + (key(1,2) /= keys(1,2,i)).or. & + (key(2,1) /= keys(2,1,i)).or. & + (key(2,2) /= keys(2,2,i)) ) then + cycle + endif connected_to_ref = -i return - endif - if (degree_x2 > 5) then - cycle - else - call i_H_j(keys(1,1,i),key,N_int,hij_elec) - if(dabs(hij_elec).lt.thresh)cycle - connected_to_ref = i - return - endif - enddo - - !DIR$ FORCEINLINE - t = det_is_not_or_may_be_in_ref(key,N_int) - if ( t ) then - return - endif - - !DIR$ LOOP COUNT (1000) - do i=N_past+1,Ndet - if ( (key(1,1) /= keys(1,1,i)).or. & - (key(1,2) /= keys(1,2,i)).or. & - (key(2,1) /= keys(2,1,i)).or. & - (key(2,2) /= keys(2,2,i)) ) then - cycle - endif - connected_to_ref = -i - return - enddo - return - - else if (Nint==3) then - - do i=N_past-1,1,-1 - degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & - popcnt(xor( key(1,2), keys(1,2,i))) + & - popcnt(xor( key(2,1), keys(2,1,i))) + & - popcnt(xor( key(2,2), keys(2,2,i))) + & - popcnt(xor( key(3,1), keys(3,1,i))) + & - popcnt(xor( key(3,2), keys(3,2,i))) - if(degree_x2 == 0)then + enddo + return + + else if (Nint==3) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + & + popcnt(xor( key(3,1), keys(3,1,i))) + & + popcnt(xor( key(3,2), keys(3,2,i))) + if(degree_x2 == 0)then + connected_to_ref = -i + return + endif + if (degree_x2 > 5) then + cycle + else + call i_H_j(keys(1,1,i),key,N_int,hij_elec) + if(dabs(hij_elec).lt.thresh)cycle + connected_to_ref = i + return + endif + enddo + + !DIR$ FORCEINLINE + t = det_is_not_or_may_be_in_ref(key,N_int) + if ( t ) then + return + endif + + do i=N_past+1,Ndet + if ( (key(1,1) /= keys(1,1,i)).or. & + (key(1,2) /= keys(1,2,i)).or. & + (key(2,1) /= keys(2,1,i)).or. & + (key(2,2) /= keys(2,2,i)).or. & + (key(3,1) /= keys(3,1,i)).or. & + (key(3,2) /= keys(3,2,i)) ) then + cycle + endif connected_to_ref = -i return - endif - if (degree_x2 > 5) then - cycle - else - call i_H_j(keys(1,1,i),key,N_int,hij_elec) - if(dabs(hij_elec).lt.thresh)cycle - connected_to_ref = i - return - endif - enddo - - !DIR$ FORCEINLINE - t = det_is_not_or_may_be_in_ref(key,N_int) - if ( t ) then - return - endif - - do i=N_past+1,Ndet - if ( (key(1,1) /= keys(1,1,i)).or. & - (key(1,2) /= keys(1,2,i)).or. & - (key(2,1) /= keys(2,1,i)).or. & - (key(2,2) /= keys(2,2,i)).or. & - (key(3,1) /= keys(3,1,i)).or. & - (key(3,2) /= keys(3,2,i)) ) then - cycle - endif - connected_to_ref = -i - return - enddo - return - - else - - do i=N_past-1,1,-1 - degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & - popcnt(xor( key(1,2), keys(1,2,i))) - !DEC$ LOOP COUNT MIN(3) - do l=2,Nint - degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) + & - popcnt(xor( key(l,2), keys(l,2,i))) - enddo - if(degree_x2 == 0)then - connected_to_ref = -i + enddo + return + + else + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + !DEC$ LOOP COUNT MIN(3) + do l=2,Nint + degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +& + popcnt(xor( key(l,2), keys(l,2,i))) + enddo + if(degree_x2 == 0)then + connected_to_ref = -i + return + endif + if (degree_x2 > 5) then + cycle + else + call i_H_j(keys(1,1,i),key,N_int,hij_elec) + if(dabs(hij_elec).lt.thresh)cycle + connected_to_ref = i + return + endif + enddo + + !DIR$ FORCEINLINE + t = det_is_not_or_may_be_in_ref(key,N_int) + if ( t ) then return - endif - if (degree_x2 > 5) then - cycle - else - call i_H_j(keys(1,1,i),key,N_int,hij_elec) - if(dabs(hij_elec).lt.thresh)cycle - connected_to_ref = i - return - endif - enddo - - !DIR$ FORCEINLINE - t = det_is_not_or_may_be_in_ref(key,N_int) - if ( t ) then - return - endif - - do i=N_past+1,Ndet - if ( (key(1,1) /= keys(1,1,i)).or. & - (key(1,2) /= keys(1,2,i)) ) then - cycle - else - connected_to_ref = -1 - !DEC$ LOOP COUNT MIN(3) - do l=2,Nint - if ( (key(l,1) /= keys(l,1,i)).or. & - (key(l,2) /= keys(l,2,i)) ) then - connected_to_ref = 0 - exit - endif - enddo - if (connected_to_ref /= 0) then - return - endif - endif - enddo - - endif - + endif + + do i=N_past+1,Ndet + if ( (key(1,1) /= keys(1,1,i)).or. & + (key(1,2) /= keys(1,2,i)) ) then + cycle + else + connected_to_ref = -1 + !DEC$ LOOP COUNT MIN(3) + do l=2,Nint + if ( (key(l,1) /= keys(l,1,i)).or. & + (key(l,2) /= keys(l,2,i)) ) then + connected_to_ref = 0 + exit + endif + enddo + if (connected_to_ref /= 0) then + return + endif + endif + enddo + + endif + end logical function det_is_not_or_may_be_in_ref(key,Nint) - implicit none - ! If true, det is not in ref - ! If false, det may be in ref - - integer, intent(in) :: key(Nint,2), Nint - integer :: key_int - integer*1 :: key_short(4) - !DIR$ ATTRIBUTES ALIGN : 32 :: key_short - equivalence (key_int,key_short) - - integer :: i, ispin - - det_is_not_or_may_be_in_ref = .False. - do ispin=1,2 - do i=1,Nint - key_int = key(i,ispin) - if ( & - key_pattern_not_in_ref(key_short(1), i, ispin) & - .or.key_pattern_not_in_ref(key_short(2), i, ispin) & - .or.key_pattern_not_in_ref(key_short(3), i, ispin) & - .or.key_pattern_not_in_ref(key_short(4), i, ispin) & - ) then - det_is_not_or_may_be_in_ref = .True. + use bitmasks + implicit none + BEGIN_DOC + ! If true, det is not in ref + ! If false, det may be in ref + END_DOC + integer(bit_kind), intent(in) :: key(Nint,2), Nint + integer(bit_kind) :: key_int + integer*1 :: key_short(bit_kind) + !DIR$ ATTRIBUTES ALIGN : 32 :: key_short + equivalence (key_int,key_short) + + integer :: i, ispin, k + + det_is_not_or_may_be_in_ref = .False. + do ispin=1,2 + do i=1,Nint + key_int = key(i,ispin) + do k=1,bit_kind + det_is_not_or_may_be_in_ref = & + det_is_not_or_may_be_in_ref .or. & + key_pattern_not_in_ref(key_short(k), i, ispin) + enddo + if(det_is_not_or_may_be_in_ref) then return - endif + endif + enddo enddo - enddo - + end BEGIN_PROVIDER [ logical, key_pattern_not_in_ref, (-128:127,N_int,2) ] - implicit none - BEGIN_DOC -! Min and max values of the integers of the keys of the reference - END_DOC - - integer :: i, j, ispin - integer :: key - integer*1 :: key_short(4) - equivalence (key,key_short) - integer :: idx - - key_pattern_not_in_ref = .True. - - do j=1,N_det - do ispin=1,2 - do i=1,N_int - key = psi_det(i,ispin,j) - key_pattern_not_in_ref( key_short(1), i, ispin ) = .False. - key_pattern_not_in_ref( key_short(2), i, ispin ) = .False. - key_pattern_not_in_ref( key_short(3), i, ispin ) = .False. - key_pattern_not_in_ref( key_short(4), i, ispin ) = .False. - enddo + use bitmasks + implicit none + BEGIN_DOC + ! Min and max values of the integers of the keys of the reference + END_DOC + + integer :: i, j, ispin + integer(bit_kind) :: key + integer*1 :: key_short(bit_kind) + equivalence (key,key_short) + integer :: idx, k + + key_pattern_not_in_ref = .True. + + do j=1,N_det + do ispin=1,2 + do i=1,N_int + key = psi_det(i,ispin,j) + do k=1,bit_kind + key_pattern_not_in_ref( key_short(k), i, ispin ) = .False. + enddo + enddo + enddo enddo - enddo - + END_PROVIDER diff --git a/src/Dets/filter_connected.irp.f b/src/Dets/filter_connected.irp.f index 18acd748..0353c9c6 100644 --- a/src/Dets/filter_connected.irp.f +++ b/src/Dets/filter_connected.irp.f @@ -95,6 +95,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) integer :: degree_x2 ASSERT (Nint > 0) + ASSERT (Nint == N_int) ASSERT (sze > 0) l=1 @@ -118,8 +119,8 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) !DIR$ LOOP COUNT (1000) do i=1,sze degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & - popcnt(xor( key1(1,2,i), key2(1,2))) + & popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & popcnt(xor( key1(2,2,i), key2(2,2))) if (degree_x2 < 5) then if(degree_x2 .ne. 0)then diff --git a/src/Dets/slater_rules.irp.f b/src/Dets/slater_rules.irp.f index c160b0c2..3f56eba3 100644 --- a/src/Dets/slater_rules.irp.f +++ b/src/Dets/slater_rules.irp.f @@ -492,9 +492,9 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) use bitmasks implicit none integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate - integer, intent(in) :: keys(Nint,2,Ndet_max) + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) double precision, intent(in) :: coef(Ndet_max,Nstate) - integer, intent(in) :: key(Nint,2) double precision, intent(out) :: i_H_psi_array(Nstate) integer :: i, ii,j @@ -503,6 +503,11 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) double precision :: hij integer :: idx(0:Ndet) + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) i_H_psi_array = 0.d0 call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) do ii=1,idx(0) @@ -512,6 +517,7 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) do j = 1, Nstate i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij enddo + print *, 'x', coef(i,1), hij, i_H_psi_array(1) enddo end diff --git a/src/Hartree_Fock_AOs/README.rst b/src/Hartree_Fock_AOs/README.rst deleted file mode 100644 index 0418aca5..00000000 --- a/src/Hartree_Fock_AOs/README.rst +++ /dev/null @@ -1,30 +0,0 @@ -=================== -Hartree-Fock Module -=================== - - -Needed Modules -============== - -.. Do not edit this section. It was auto-generated from the -.. NEEDED_MODULES file. - -* `AOs `_ -* `BiInts `_ -* `Bitmask `_ -* `Electrons `_ -* `Ezfio_files `_ -* `MonoInts `_ -* `MOs `_ -* `Nuclei `_ -* `Output `_ -* `Utils `_ - -Documentation -============= - -.. Do not edit this section. It was auto-generated from the -.. NEEDED_MODULES file. - - - diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index d85c28b1..d09721ef 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation +AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation Selection diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst index 41a8d6a9..83392b83 100644 --- a/src/Perturbation/README.rst +++ b/src/Perturbation/README.rst @@ -82,6 +82,26 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +`pt2_epstein_nesbet `_ + compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + .br + for the various n_st states. + .br + c_pert(i) = /( E(i) - ) + .br + e_2_pert(i) = ^2/( E(i) - ) + .br + +`pt2_epstein_nesbet_2x2 `_ + compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution + .br + for the various n_st states. + .br + e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) + .br + c_pert(i) = e_2_pert(i)/ + .br + Needed Modules diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f index 7f8474b0..4c1e7c45 100644 --- a/src/Perturbation/epstein_nesbet.irp.f +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -24,7 +24,7 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array) H_pert_diag = diag_H_mat_elem(det_pert,Nint) do i =1,n_st - c_pert(i) = i_H_psi_array(i) / (E_ref(i) - H_pert_diag) + c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag) e_2_pert(i) = c_pert(i) * i_H_psi_array(i) enddo @@ -53,12 +53,16 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet double precision :: diag_H_mat_elem,delta_e ASSERT (Nint == N_int) ASSERT (Nint > 0) - call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array) + print *, 'coefs' + print *, psi_ref_coef(1:N_det_ref,1) + print *, '-----' + call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,N_det_ref,psi_ref_size,n_st,i_H_psi_array) H_pert_diag = diag_H_mat_elem(det_pert,Nint) do i =1,n_st - delta_e = H_pert_diag - E_ref(i) + delta_e = H_pert_diag - reference_energy(i) e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) c_pert(i) = e_2_pert(i)/i_H_psi_array(i) enddo - + print *, e_2_pert, delta_e , i_H_psi_array + end diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f index 5e672529..3365e67c 100644 --- a/src/Perturbation/perturbation_template.f +++ b/src/Perturbation/perturbation_template.f @@ -2,7 +2,7 @@ BEGIN_SHELL [ /usr/bin/env python ] import perturbation END_SHELL -subroutine perturb_buffer_$PERT(buffer,buffer_size,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint) +subroutine perturb_buffer_$PERT(buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint) implicit none BEGIN_DOC ! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply @@ -11,7 +11,8 @@ subroutine perturb_buffer_$PERT(buffer,buffer_size,sum_e_2_pert,sum_norm_pert,su integer, intent(in) :: Nint, N_st, buffer_size integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size) - double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st),sum_H_pert_diag(N_st) + double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st) + double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st) double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag integer :: i,k, c_ref integer :: connected_to_ref @@ -33,6 +34,8 @@ subroutine perturb_buffer_$PERT(buffer,buffer_size,sum_e_2_pert,sum_norm_pert,su c_pert,e_2_pert,H_pert_diag,Nint,n_det_ref,n_st) do k = 1,N_st + e_2_pert_buffer(k,i) = e_2_pert(k) + coef_pert_buffer(k,i) = c_pert(k) sum_norm_pert(k) += c_pert(k) * c_pert(k) sum_e_2_pert(k) += e_2_pert(k) sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag