diff --git a/src/BiInts/ao_bi_integrals.irp.f b/src/BiInts/ao_bi_integrals.irp.f index 01bab63d..faac8a42 100644 --- a/src/BiInts/ao_bi_integrals.irp.f +++ b/src/BiInts/ao_bi_integrals.irp.f @@ -26,7 +26,6 @@ double precision function ao_bielec_integral(i,j,k,l) ao_bielec_integral = 0.d0 double precision :: thresh thresh = ao_integrals_threshold - thresh = 0.d0 if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then do p = 1, 3 diff --git a/src/BiInts/mo_bi_integrals.irp.f b/src/BiInts/mo_bi_integrals.irp.f index e7ffde60..23565b66 100644 --- a/src/BiInts/mo_bi_integrals.irp.f +++ b/src/BiInts/mo_bi_integrals.irp.f @@ -88,7 +88,6 @@ subroutine add_integrals_to_map(mask_ijkl) call wall_time(wall_1) call cpu_time(cpu_1) - mo_integrals_threshold = 0.d0 !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & !$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,& diff --git a/src/Bitmask/README.rst b/src/Bitmask/README.rst index 4169d164..057bd7eb 100644 --- a/src/Bitmask/README.rst +++ b/src/Bitmask/README.rst @@ -57,15 +57,40 @@ Documentation `full_ijkl_bitmask `_ Bitmask to include all possible MOs +`generators_bitmask `_ + Bitmasks for generator determinants. (N_int, alpha/beta, hole/particle, generator). + 3rd index is : + * 1 : hole for single exc + * 1 : particle for single exc + * 3 : hole for 1st exc of double + * 4 : particle for 1st exc of double + * 5 : hole for 2dn exc of double + * 6 : particle for 2dn exc of double + `hf_bitmask `_ Hartree Fock bit mask +`i_bitmask_gen `_ + Current bitmask for the generators + +`i_bitmask_ref `_ + Current bitmask for the reference + +`n_generators_bitmask `_ + Number of bitmasks for generators + `n_int `_ Number of 64-bit integers needed to represent determinants as binary strings +`n_reference_bitmask `_ + Number of bitmasks for reference + `ref_bitmask `_ Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask +`reference_bitmask `_ + Bitmasks for reference determinants. (N_int, alpha/beta, hole/particle, reference) + `bitstring_to_hexa `_ Transform a bit string to a string in hexadecimal format for printing diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 71833017..17d33c54 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -107,12 +107,17 @@ BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6,N_generators_ if (exists) then call ezfio_get_bitmasks_generators(generators_bitmask) else - generators_bitmask(:,:,s_hole ,1) = HF_bitmask - generators_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) - generators_bitmask(:,:,d_hole1,1) = HF_bitmask - generators_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) - generators_bitmask(:,:,d_hole2,1) = HF_bitmask - generators_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) + integer :: k, ispin + do k=1,N_generators_bitmask + do ispin=1,2 + generators_bitmask(:,ispin,s_hole ,k) = full_ijkl_bitmask(:,d_hole1) + generators_bitmask(:,ispin,s_part ,k) = full_ijkl_bitmask(:,d_part1) + generators_bitmask(:,ispin,d_hole1,k) = full_ijkl_bitmask(:,d_hole1) + generators_bitmask(:,ispin,d_part1,k) = full_ijkl_bitmask(:,d_part1) + generators_bitmask(:,ispin,d_hole2,k) = full_ijkl_bitmask(:,d_hole2) + generators_bitmask(:,ispin,d_part2,k) = full_ijkl_bitmask(:,d_part2) + enddo + enddo call ezfio_set_bitmasks_generators(generators_bitmask) endif diff --git a/src/Bitmask/bitmasks_module.f90 b/src/Bitmask/bitmasks_module.f90 index 5e5e2f02..c542fac6 100644 --- a/src/Bitmask/bitmasks_module.f90 +++ b/src/Bitmask/bitmasks_module.f90 @@ -2,10 +2,10 @@ module bitmasks integer, parameter :: bit_kind_shift = 6 ! 5: 32 bits, 6: 64 bits integer, parameter :: bit_kind_size = 64 integer, parameter :: bit_kind = 64/8 - integer, parameter :: s_hole = 1 - integer, parameter :: s_part = 2 - integer, parameter :: d_hole1 = 3 - integer, parameter :: d_part1 = 4 - integer, parameter :: d_hole2 = 5 - integer, parameter :: d_part2 = 6 + integer, parameter :: d_hole1 = 1 + integer, parameter :: d_part1 = 2 + integer, parameter :: d_hole2 = 3 + integer, parameter :: d_part2 = 4 + integer, parameter :: s_hole = 5 + integer, parameter :: s_part = 6 end module diff --git a/src/CISD_SC2_selected/README.rst b/src/CISD_SC2_selected/README.rst index e8bb9563..b6206850 100644 --- a/src/CISD_SC2_selected/README.rst +++ b/src/CISD_SC2_selected/README.rst @@ -23,7 +23,7 @@ Needed Modules * `BiInts `_ * `Bitmask `_ * `CISD `_ -* `CISD_SC2 `_ +* `SC2 `_ * `CISD_selected `_ * `Dets `_ * `Electrons `_ diff --git a/src/Dets/H_apply.irp.f b/src/Dets/H_apply.irp.f index c2302432..6e240a4b 100644 --- a/src/Dets/H_apply.irp.f +++ b/src/Dets/H_apply.irp.f @@ -19,7 +19,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ] ! Uninitialized. Filled by H_apply subroutines. END_DOC integer :: iproc, sze - sze = 100 + sze = 10000 if (.not.associated(H_apply_buffer)) then allocate(H_apply_buffer(0:nproc-1)) iproc = 0 diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 1e1500e2..c1dd13dc 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -47,14 +47,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene occ_hole_tmp(N_int*bit_kind_size,2)) $init_thread - !print*,'key_in !!' - !call print_key(key_in) - !print*,'hole_1, particl_1' - !call print_key(hole_1) - !call print_key(particl_1) - !print*,'hole_2, particl_2' - !call print_key(hole_2) - !call print_key(particl_2) !!!! First couple hole particle @@ -352,9 +344,11 @@ subroutine $subroutine($params_main) $decls_main - integer :: i_generator, k, nmax - double precision :: wall_0, wall_1, wall_2, d + integer :: i_generator, nmax + double precision :: wall_0, wall_1, wall_2 integer(omp_lock_kind) :: lck + integer(bit_kind), allocatable :: mask(:,:,:) + integer :: ispin, k PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_selectors psi_generators PROVIDE psi_det_sorted_bit coef_hf_selector @@ -363,22 +357,45 @@ subroutine $subroutine($params_main) call wall_time(wall_1) !$ call omp_init_lock(lck) !$OMP PARALLEL DEFAULT(SHARED) & - !$OMP PRIVATE(i_generator,wall_2) + !$OMP PRIVATE(i_generator,wall_2,ispin,k,mask) + allocate( mask(N_int,2,6) ) !$OMP DO SCHEDULE(guided) do i_generator=1,nmax if (abort_here) then cycle endif $skip + + ! Create bit masks for holes and particles + do ispin=1,2 + do k=1,N_int + mask(k,ispin,s_hole) = & + iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), & + psi_generators(k,ispin,i_generator) ) + mask(k,ispin,s_part) = & + iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), & + not(psi_generators(k,ispin,i_generator)) ) + mask(k,ispin,d_hole1) = & + iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), & + psi_generators(k,ispin,i_generator) ) + mask(k,ispin,d_part1) = & + iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), & + not(psi_generators(k,ispin,i_generator)) ) + mask(k,ispin,d_hole2) = & + iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), & + psi_generators(k,ispin,i_generator) ) + mask(k,ispin,d_part2) = & + iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), & + not (psi_generators(k,ispin,i_generator)) ) + enddo + enddo + call $subroutine_diexc(psi_generators(1,1,i_generator), & - generators_bitmask(1,1,d_hole1,i_bitmask_gen), & - generators_bitmask(1,1,d_part1,i_bitmask_gen), & - generators_bitmask(1,1,d_hole2,i_bitmask_gen), & - generators_bitmask(1,1,d_part2,i_bitmask_gen), & + mask(1,1,d_hole1), mask(1,1,d_part1), & + mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator $params_post) call $subroutine_monoexc(psi_generators(1,1,i_generator), & - generators_bitmask(1,1,s_hole ,i_bitmask_gen), & - generators_bitmask(1,1,s_part ,i_bitmask_gen), & + mask(1,1,s_hole ), mask(1,1,s_part ), & i_generator $params_post) !$ call omp_set_lock(lck) call wall_time(wall_2) @@ -390,23 +407,46 @@ subroutine $subroutine($params_main) !$ call omp_unset_lock(lck) enddo !$OMP END DO + deallocate( mask ) !$OMP END PARALLEL !$ call omp_destroy_lock(lck) + allocate( mask(N_int,2,6) ) do i_generator=nmax+1,N_det_generators if (abort_here) then exit endif $skip + + ! Create bit masks for holes and particles + do ispin=1,2 + do k=1,N_int + mask(k,ispin,s_hole) = & + iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), & + psi_generators(k,ispin,i_generator) ) + mask(k,ispin,s_part) = & + iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), & + not(psi_generators(k,ispin,i_generator)) ) + mask(k,ispin,d_hole1) = & + iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), & + psi_generators(k,ispin,i_generator) ) + mask(k,ispin,d_part1) = & + iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), & + not(psi_generators(k,ispin,i_generator)) ) + mask(k,ispin,d_hole2) = & + iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), & + psi_generators(k,ispin,i_generator) ) + mask(k,ispin,d_part2) = & + iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), & + not(psi_generators(k,ispin,i_generator)) ) + enddo + enddo call $subroutine_diexc(psi_generators(1,1,i_generator), & - generators_bitmask(1,1,d_hole1,i_bitmask_gen), & - generators_bitmask(1,1,d_part1,i_bitmask_gen), & - generators_bitmask(1,1,d_hole2,i_bitmask_gen), & - generators_bitmask(1,1,d_part2,i_bitmask_gen), & + mask(1,1,d_hole1), mask(1,1,d_part1), & + mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator $params_post) call $subroutine_monoexc(psi_generators(1,1,i_generator), & - generators_bitmask(1,1,s_hole ,i_bitmask_gen), & - generators_bitmask(1,1,s_part ,i_bitmask_gen), & + mask(1,1,s_hole ), mask(1,1,s_part ), & i_generator $params_post) call wall_time(wall_2) $printout_always @@ -419,6 +459,7 @@ subroutine $subroutine($params_main) $copy_buffer $generate_psi_guess abort_here = abort_all + deallocate( mask ) end diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f index 2b7c028c..6d7dcab1 100644 --- a/src/Dets/connected_to_ref.irp.f +++ b/src/Dets/connected_to_ref.irp.f @@ -56,16 +56,40 @@ logical function is_in_wavefunction(key,Nint,Ndet) endif enddo if (is_in_wavefunction) then - return + exit endif endif i += 1 if (i > N_det) then - exit + return + ! exit endif enddo + +! DEBUG is_in_wf +! if (is_in_wavefunction) then +! degree = 1 +! do i=1,N_det +! integer :: degree +! call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) +! if (degree == 0) then +! exit +! endif +! enddo +! if (degree /=0) then +! stop 'pouet 1' +! endif +! else +! do i=1,N_det +! call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) +! if (degree == 0) then +! stop 'pouet 2' +! endif +! enddo +! endif +! END DEBUG is_in_wf end integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) @@ -92,14 +116,10 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) N_past = max(1,N_past_in) if (Nint == 1) then - do i=N_past-1,1,-1 + do i=1,N_past-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 + if (degree_x2 > 4) then cycle else connected_to_ref = i @@ -112,16 +132,12 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) else if (Nint==2) then - do i=N_past-1,1,-1 + do i=1,N_past-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 + if (degree_x2 > 4) then cycle else connected_to_ref = i @@ -133,18 +149,14 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) else if (Nint==3) then - do i=N_past-1,1,-1 + do i=1,N_past-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 + if (degree_x2 > 4) then cycle else connected_to_ref = i @@ -156,7 +168,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) else - do i=N_past-1,1,-1 + do i=1,N_past-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) @@ -164,11 +176,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) 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 + if (degree_x2 > 4) then cycle else connected_to_ref = i diff --git a/src/Full_CI/H_apply.irp.f b/src/Full_CI/H_apply.irp.f index c9b44cbd..63ea8252 100644 --- a/src/Full_CI/H_apply.irp.f +++ b/src/Full_CI/H_apply.irp.f @@ -6,5 +6,9 @@ s = H_apply("FCI",openmp=True) s.set_selection_pt2("epstein_nesbet_2x2") print s +s = H_apply("FCI_PT2",openmp=True) +s.set_perturbation("epstein_nesbet_2x2") +print s + END_SHELL diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f index e5a60f25..8b30b577 100644 --- a/src/Full_CI/full_ci.irp.f +++ b/src/Full_CI/full_ci.irp.f @@ -5,9 +5,9 @@ program cisd double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) integer :: N_st, degree - character*(64) :: perturbation N_st = N_states allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + character*(64) :: perturbation pt2 = 1.d0 diag_algorithm = "Lapack" diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f index 68236da2..d2f6307a 100644 --- a/src/Generators_full/generators.irp.f +++ b/src/Generators_full/generators.irp.f @@ -33,8 +33,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ] N_det_generators = N_det do i=1,N_det norm = norm + psi_average_norm_contrib_sorted(i) - if (norm > threshold_generators) then - N_det_generators = i-1 + if (norm >= threshold_generators) then + N_det_generators = i exit endif enddo diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f index eba7530e..60d050f5 100644 --- a/src/Perturbation/epstein_nesbet.irp.f +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -66,14 +66,24 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) h = diag_H_mat_elem(det_pert,Nint) do i =1,N_st - delta_e = h - CI_electronic_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))) - if (dabs(i_H_psi_array(i)) > 1.d-6) then - c_pert(i) = e_2_pert(i)/i_H_psi_array(i) + if (i_H_psi_array(i) /= 0.d0) then + delta_e = h - CI_electronic_energy(i) + if (delta_e > 0.d0) then + 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))) + else + 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))) + endif + if (dabs(i_H_psi_array(i)) > 1.d-6) then + c_pert(i) = e_2_pert(i)/i_H_psi_array(i) + else + c_pert(i) = 0.d0 + endif + H_pert_diag(i) = h*c_pert(i)*c_pert(i) else - c_pert(i) = -1.d0 + e_2_pert(i) = 0.d0 + c_pert(i) = 0.d0 + H_pert_diag(i) = 0.d0 endif - H_pert_diag(i) = h*c_pert(i)*c_pert(i) enddo end diff --git a/src/SC2/README.rst b/src/SC2/README.rst index 96e489ed..30269d2d 100644 --- a/src/SC2/README.rst +++ b/src/SC2/README.rst @@ -8,7 +8,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`cisd_sc2 `_ +`cisd_sc2 `_ CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) .br dets_in : bitmasks corresponding to determinants @@ -24,25 +24,51 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`repeat_excitation `_ +`repeat_excitation `_ Undocumented -`cisd `_ +`cisd `_ Undocumented -`ci_sc2_eigenvectors `_ +`ci_sc2_eigenvectors `_ Eigenvectors/values of the CI matrix -`ci_sc2_electronic_energy `_ +`ci_sc2_electronic_energy `_ Eigenvectors/values of the CI matrix -`ci_sc2_energy `_ +`ci_sc2_energy `_ N_states lowest eigenvalues of the CI matrix -`diagonalize_ci_sc2 `_ +`diagonalize_ci_sc2 `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix +`pt2_epstein_nesbet_sc2_projected `_ + compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + .br + for the various N_st states, + .br + but with the correction in the denominator + .br + comming from the interaction of that determinant with all the others determinants + .br + that can be repeated by repeating all the double excitations + .br + : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + .br + that could be repeated to this determinant. + .br + In addition, for the perturbative energetic contribution you have the standard second order + .br + e_2_pert = ^2/(Delta_E) + .br + and also the purely projected contribution + .br + H_pert_diag = c_pert + +`repeat_all_e_corr `_ + Undocumented + Needed Modules diff --git a/src/Selectors_full/e_corr_selectors.irp.f b/src/Selectors_full/e_corr_selectors.irp.f index 580c478b..3ea78540 100644 --- a/src/Selectors_full/e_corr_selectors.irp.f +++ b/src/Selectors_full/e_corr_selectors.irp.f @@ -51,7 +51,11 @@ END_PROVIDER E_corr_per_selectors(i) = -1000.d0 endif enddo - inv_selectors_coef_hf = 1.d0/coef_hf_selector + if (dabs(coef_hf_selector) > 1.d-8) then + inv_selectors_coef_hf = 1.d0/coef_hf_selector + else + inv_selectors_coef_hf = 0.d0 + endif do i = 1,n_double_selectors E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf enddo