diff --git a/plugins/FOBOCI/dress_simple.irp.f b/plugins/FOBOCI/dress_simple.irp.f index 99566a8e..021aa422 100644 --- a/plugins/FOBOCI/dress_simple.irp.f +++ b/plugins/FOBOCI/dress_simple.irp.f @@ -279,9 +279,11 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener do k = 1, N_states accu = 0.d0 do j =1, Ndet_generators + print*,'',eigvectors(j,i) , psi_coef_ref(j,k) accu += eigvectors(j,i) * psi_coef_ref(j,k) enddo - if(dabs(accu).ge.0.8d0)then + print*,'accu = ',accu + if(dabs(accu).ge.0.72d0)then i_good_state(0) +=1 i_good_state(i_good_state(0)) = i endif diff --git a/plugins/FOBOCI/fobo_scf.irp.f b/plugins/FOBOCI/fobo_scf.irp.f index 8656b633..0b0902b0 100644 --- a/plugins/FOBOCI/fobo_scf.irp.f +++ b/plugins/FOBOCI/fobo_scf.irp.f @@ -1,6 +1,6 @@ program foboscf implicit none - call run_prepare +!call run_prepare no_oa_or_av_opt = .True. touch no_oa_or_av_opt call routine_fobo_scf diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f index acc78912..9ccbb981 100644 --- a/src/AO_Basis/aos.irp.f +++ b/src/AO_Basis/aos.irp.f @@ -17,7 +17,7 @@ END_PROVIDER call ezfio_get_ao_basis_ao_prim_num_max(ao_prim_num_max) integer :: align_double ao_prim_num_max_align = align_double(ao_prim_num_max) - END_PROVIDER +END_PROVIDER BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num_align,ao_prim_num_max) ] &BEGIN_PROVIDER [ double precision, ao_coef_normalization_factor, (ao_num) ] @@ -109,6 +109,7 @@ END_PROVIDER BEGIN_PROVIDER [ integer, ao_l, (ao_num) ] &BEGIN_PROVIDER [ integer, ao_l_max ] +&BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ] implicit none BEGIN_DOC ! ao_l = l value of the AO: a+b+c in x^a y^b z^c @@ -116,6 +117,7 @@ END_PROVIDER integer :: i do i=1,ao_num ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) + ao_l_char(i) = l_to_charater(ao_l(i)) enddo ao_l_max = maxval(ao_l) END_PROVIDER @@ -143,20 +145,6 @@ integer function ao_power_index(nx,ny,nz) ao_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1 end - BEGIN_PROVIDER [ integer, ao_l, (ao_num) ] -&BEGIN_PROVIDER [ integer, ao_l_max ] -&BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ] - implicit none - BEGIN_DOC -! ao_l = l value of the AO: a+b+c in x^a y^b z^c - END_DOC - integer :: i - do i=1,ao_num - ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) - ao_l_char(i) = l_to_charater(ao_l(i)) - enddo - ao_l_max = maxval(ao_l) -END_PROVIDER BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)] BEGIN_DOC diff --git a/src/Determinants/two_body_dm_map.irp.f b/src/Determinants/two_body_dm_map.irp.f new file mode 100644 index 00000000..eccb9741 --- /dev/null +++ b/src/Determinants/two_body_dm_map.irp.f @@ -0,0 +1,210 @@ + +use map_module + +BEGIN_PROVIDER [ type(map_type), two_body_dm_ab_map ] + implicit none + BEGIN_DOC + ! Map of the two body density matrix elements for the alpha/beta elements + END_DOC + integer(key_kind) :: key_max + integer(map_size_kind) :: sze + call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max) + sze = key_max + call map_init(two_body_dm_ab_map,sze) + print*, 'two_body_dm_ab_map initialized' +END_PROVIDER + +subroutine insert_into_two_body_dm_ab_map(n_product,buffer_i, buffer_values, thr) + use map_module + implicit none + + BEGIN_DOC + ! Create new entry into two_body_dm_ab_map, or accumulate in an existing entry + END_DOC + + integer, intent(in) :: n_product + integer(key_kind), intent(inout) :: buffer_i(n_product) + real(integral_kind), intent(inout) :: buffer_values(n_product) + real(integral_kind), intent(in) :: thr + call map_update(two_body_dm_ab_map, buffer_i, buffer_values, n_product, thr) +end + +double precision function get_two_body_dm_ab_map_element(i,j,k,l,map) + use map_module + implicit none + BEGIN_DOC + ! Returns one value of the wo body density matrix \rho_{ijkl}^{\alpha \beta} defined as : + ! \rho_{ijkl}^{\alpha \beta } = <\Psi|a^{\dagger}_{i\alpha} a^{\dagger}_{j\beta} a_{k\beta} a_{l\alpha}|\Psi> + END_DOC + PROVIDE two_body_dm_ab_map + + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx + type(map_type), intent(inout) :: map + real(integral_kind) :: tmp + PROVIDE two_body_dm_in_map + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(two_body_dm_ab_map,idx,tmp) + get_two_body_dm_ab_map_element = dble(tmp) +end + +subroutine get_get_two_body_dm_ab_map_elements(j,k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple elements of the \rho_{ijkl}^{\alpha \beta }, all + ! i for j,k,l fixed. + END_DOC + integer, intent(in) :: j,k,l, sze + double precision, intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer(key_kind) :: hash(sze) + real(integral_kind) :: tmp_val(sze) + PROVIDE two_body_dm_in_map + + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,hash(i)) + enddo + + if (key_kind == 8) then + call map_get_many(two_body_dm_ab_map, hash, out_val, sze) + else + call map_get_many(two_body_dm_ab_map, hash, tmp_val, sze) + ! Conversion to double precision + do i=1,sze + out_val(i) = dble(tmp_val(i)) + enddo + endif +end + +BEGIN_PROVIDER [ logical, two_body_dm_in_map ] + implicit none + + BEGIN_DOC + ! If True, the map of the two body density matrix alpha/beta is provided + END_DOC + + two_body_dm_in_map = .True. + call add_values_to_two_body_dm_map(full_ijkl_bitmask_4) +END_PROVIDER + +subroutine add_values_to_two_body_dm_map(mask_ijkl) + use bitmasks + use map_module + implicit none + + BEGIN_DOC + ! Adds values to the map of two_body_dm according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijkl(N_int,4) + integer :: degree + + PROVIDE mo_coef psi_coef psi_det + + integer :: exc(0:2,2,2) + integer :: h1,h2,p1,p2,s1,s2 + double precision :: phase + double precision :: contrib + integer(key_kind),allocatable :: buffer_i(:) + double precision ,allocatable :: buffer_value(:) + integer :: size_buffer + integer :: n_elements + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,k,l,m + + size_buffer = min(mo_tot_num*mo_tot_num*mo_tot_num,16000000) + + allocate(buffer_i(size_buffer),buffer_value(size_buffer)) + + n_elements = 0 + do i = 1, N_det ! i == |I> + call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) + do j = i+1, N_det ! j == 2)cycle + call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int) + contrib = psi_coef(i,1) * psi_coef(j,1) * phase + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************ + + if(s1==s2)cycle ! Only the alpha/beta two body density matrix + ! * c_I * c_J + n_elements += 1 + buffer_value(n_elements) = contrib + !DEC$ FORCEINLINE + call mo_bielec_integrals_index(h1,h2,p1,p2,buffer_i(n_elements)) + if (n_elements == size_buffer) then + call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_elements = 0 + endif + + else ! case of the SINGLE EXCITATIONS *************************************************** + + if(s1==1)then ! Mono alpha : + do k = 1, elec_beta_num + m = occ(k,2) + n_elements += 1 + buffer_value(n_elements) = contrib + ! * c_I * c_J + call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements)) + if (n_elements == size_buffer) then + call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_elements = 0 + endif + enddo + else ! Mono Beta : + do k = 1, elec_alpha_num + m = occ(k,1) + n_elements += 1 + buffer_value(n_elements) = contrib + ! * c_I * c_J + call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements)) + if (n_elements == size_buffer) then + call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_elements = 0 + endif + enddo + endif + + endif + enddo + enddo + print*,'n_elements = ',n_elements + call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + +end + +BEGIN_PROVIDER [double precision, two_body_dm_ab_diag, (mo_tot_num, mo_tot_num)] + implicit none + integer :: i,j,k,l,m + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + double precision :: contrib + BEGIN_DOC + ! two_body_dm_ab_diag(k,m) = <\Psi | n_(k\alpha) n_(m\beta) | \Psi> + + END_DOC + two_body_dm_ab_diag = 0.d0 + do i = 1, N_det ! i == |I> + call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) + contrib = psi_coef(i,1)**2 + do j = 1, elec_beta_num + k = occ(j,2) + do l = 1, elec_beta_num + m = occ(l,1) + two_body_dm_ab_diag(k,m) += contrib + enddo + enddo + enddo +END_PROVIDER +