diff --git a/plugins/local/basis_correction/print_routine.irp.f b/plugins/local/basis_correction/print_routine.irp.f index b3b38673..8879fd5d 100644 --- a/plugins/local/basis_correction/print_routine.irp.f +++ b/plugins/local/basis_correction/print_routine.irp.f @@ -22,53 +22,58 @@ subroutine print_basis_correction print*, '****************************************' print*, '****************************************' print*, 'mu_of_r_potential = ',mu_of_r_potential - if(mu_of_r_potential.EQ."hf".or.mu_of_r_potential.EQ."hf_old".or.mu_of_r_potential.EQ."hf_sparse")then - print*, '' - print*,'Using a HF-like two-body density to define mu(r)' - print*,'This assumes that HF is a qualitative representation of the wave function ' - print*,'********************************************' - print*,'Functionals more suited for weak correlation' - print*,'********************************************' - print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) ' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate) - enddo - print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) - enddo + if(mu_of_r_potential.EQ."hf".or. & + mu_of_r_potential.EQ."hf_old".or.& + mu_of_r_potential.EQ."hf_sparse".or.& + mu_of_r_potential.EQ."proj")then + print*, '' + print*,'Using a HF-like two-body density to define mu(r)' + print*,'This assumes that HF is a qualitative representation of the wave function ' + print*,'********************************************' + print*,'Functionals more suited for weak correlation' + print*,'********************************************' + print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) ' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate) + enddo + print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) + enddo - else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then - print*, '' - print*,'Using a CAS-like two-body density to define mu(r)' - print*,'This assumes that the CAS is a qualitative representation of the wave function ' - print*,'********************************************' - print*,'Functionals more suited for weak correlation' - print*,'********************************************' - print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) ' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate) - enddo - print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) - enddo - print*,'' - print*,'********************************************' - print*,'********************************************' - print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) ' - print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate) - enddo - print*,'' - print*,'********************************************' - print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)' - print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate) - enddo - print*,'' + else if(mu_of_r_potential.EQ."cas_full".or. & + mu_of_r_potential.EQ."cas_truncated".or. & + mu_of_r_potential.EQ."pure_act") then + print*, '' + print*,'Using a CAS-like two-body density to define mu(r)' + print*,'This assumes that the CAS is a qualitative representation of the wave function ' + print*,'********************************************' + print*,'Functionals more suited for weak correlation' + print*,'********************************************' + print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) ' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate) + enddo + print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) + enddo + print*,'' + print*,'********************************************' + print*,'********************************************' + print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) ' + print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate) + enddo + print*,'' + print*,'********************************************' + print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)' + print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate) + enddo + print*,'' endif print*,'' diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index d15ebdc3..36061ef0 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -178,7 +178,7 @@ END_PROVIDER rank_max = np ! Avoid too large arrays when there are many electrons if (elec_num > 10) then - rank_max = min(np,20*elec_num*elec_num) + rank_max = min(np,25*elec_num*elec_num) endif call mmap_create_d('', (/ ndim8, rank_max /), .False., .True., map) diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index fb376ce1..1cb7617e 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -54,6 +54,7 @@ double precision function ao_two_e_integral(i, j, k, l) else if (use_only_lr) then ao_two_e_integral = ao_two_e_integral_erf(i, j, k, l) + return else if (do_schwartz_accel(i,j,k,l)) then diff --git a/src/bitmask/bitmasks_routines.irp.f b/src/bitmask/bitmasks_routines.irp.f index 9c6f4f0c..be0a2dd1 100644 --- a/src/bitmask/bitmasks_routines.irp.f +++ b/src/bitmask/bitmasks_routines.irp.f @@ -283,33 +283,16 @@ subroutine print_det_one_dimension(string,Nint) end -logical function is_integer_in_string(bite,string,Nint) - use bitmasks +logical function is_integer_in_string(orb,bitmask,Nint) + use bitmasks implicit none - integer, intent(in) :: bite,Nint - integer(bit_kind), intent(in) :: string(Nint) - integer(bit_kind) :: string_bite(Nint) - integer :: i,itot,itot_and - character*(2048) :: output(1) - string_bite = 0_bit_kind - call set_bit_to_integer(bite,string_bite,Nint) - itot = 0 - itot_and = 0 - is_integer_in_string = .False. -!print*,'' -!print*,'' -!print*,'bite = ',bite -!call bitstring_to_str( output(1), string_bite, Nint ) -! print *, trim(output(1)) -!call bitstring_to_str( output(1), string, Nint ) -! print *, trim(output(1)) - do i = 1, Nint - itot += popcnt(string(i)) - itot_and += popcnt(ior(string(i),string_bite(i))) - enddo -!print*,'itot,itot_and',itot,itot_and - if(itot == itot_and)then - is_integer_in_string = .True. - endif -!pause + BEGIN_DOC +! Checks is the orbital orb is set to 1 in the bit string + END_DOC + integer, intent(in) :: orb, Nint + integer(bit_kind), intent(in) :: bitmask(Nint) + integer :: j, k + k = ishft(orb-1,-bit_kind_shift)+1 + j = orb-ishft(k-1,bit_kind_shift)-1 + is_integer_in_string = iand(bitmask(k), ibset(0_bit_kind, j)) /= 0_bit_kind end diff --git a/src/cas_based_on_top/on_top_cas_prov.irp.f b/src/cas_based_on_top/on_top_cas_prov.irp.f index dd93ed40..9f9a1f06 100644 --- a/src/cas_based_on_top/on_top_cas_prov.irp.f +++ b/src/cas_based_on_top/on_top_cas_prov.irp.f @@ -15,14 +15,17 @@ pure_act_on_top_of_r = 0.d0 do l = 1, n_act_orb phi_l = act_mos_in_r_array(l,ipoint) + if (dabs(phi_l) < 1.d-12) cycle do k = 1, n_act_orb - phi_k = act_mos_in_r_array(k,ipoint) + phi_k = act_mos_in_r_array(k,ipoint) * phi_l + if (dabs(phi_k) < 1.d-12) cycle do j = 1, n_act_orb - phi_j = act_mos_in_r_array(j,ipoint) + phi_j = act_mos_in_r_array(j,ipoint) * phi_k + if (dabs(phi_j) < 1.d-12) cycle do i = 1, n_act_orb - phi_i = act_mos_in_r_array(i,ipoint) - ! 1 2 1 2 - pure_act_on_top_of_r += act_2_rdm_ab_mo(i,j,k,l,istate) * phi_i * phi_j * phi_k * phi_l + phi_i = act_mos_in_r_array(i,ipoint) * phi_j + ! 1 2 1 2 + pure_act_on_top_of_r = pure_act_on_top_of_r + act_2_rdm_ab_mo(i,j,k,l,istate) * phi_i !* phi_j * phi_k * phi_l enddo enddo enddo diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index c8822776..ffd8c9d5 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -8,21 +8,14 @@ BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)] END_DOC implicit none - integer :: i, j - double precision :: tmp_array(ao_num), r(3) + integer :: i !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,tmp_array,j) & - !$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points) + !$OMP PRIVATE (i) & + !$OMP SHARED(aos_in_r_array,n_points_final_grid,final_grid_points) do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_aos_at_r(r, tmp_array) - do j = 1, ao_num - aos_in_r_array(j,i) = tmp_array(j) - enddo + call give_all_aos_at_r(final_grid_points(1,i), aos_in_r_array(1,i)) enddo !$OMP END PARALLEL DO @@ -42,7 +35,7 @@ BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_ do i = 1, n_points_final_grid do j = 1, ao_num - aos_in_r_array_transp(i,j) = aos_in_r_array(j,i) + aos_in_r_array_transp(i,j) = aos_in_r_array(j,i) enddo enddo @@ -62,27 +55,29 @@ BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_gri implicit none integer :: i, j, m - double precision :: aos_array(ao_num), r(3) - double precision :: aos_grad_array(3,ao_num) + double precision :: r(3) + double precision, allocatable :: aos_grad_array(:,:), aos_array(:) - !$OMP PARALLEL DO & + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) & + !$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) & !$OMP SHARED(aos_grad_in_r_array,n_points_final_grid,ao_num,final_grid_points) + allocate(aos_grad_array(3,ao_num), aos_array(ao_num)) + + !$OMP DO do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) + call give_all_aos_and_grad_at_r(final_grid_points(1,i),aos_array,aos_grad_array) do m = 1, 3 do j = 1, ao_num aos_grad_in_r_array(j,i,m) = aos_grad_array(m,j) enddo enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO + deallocate(aos_grad_array,aos_array) + !$OMP END PARALLEL -END_PROVIDER +END_PROVIDER ! --- @@ -116,7 +111,7 @@ END_PROVIDER enddo enddo enddo - END_PROVIDER + END_PROVIDER BEGIN_PROVIDER [double precision, aos_lapl_in_r_array, (3,ao_num,n_points_final_grid)] implicit none @@ -126,32 +121,32 @@ END_PROVIDER ! k = 1 : x, k= 2, y, k 3, z END_DOC integer :: i,j,m - double precision :: aos_array(ao_num), r(3) - double precision :: aos_grad_array(3,ao_num) - double precision :: aos_lapl_array(3,ao_num) - !$OMP PARALLEL DO & + double precision, allocatable :: aos_lapl_array(:,:), aos_grad_array(:,:), aos_array(:) + + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,aos_array,aos_grad_array,aos_lapl_array,j,m) & + !$OMP PRIVATE (i,aos_array,aos_grad_array,aos_lapl_array,j,m) & !$OMP SHARED(aos_lapl_in_r_array,n_points_final_grid,ao_num,final_grid_points) + allocate( aos_array(ao_num), aos_grad_array(3,ao_num), aos_lapl_array(3,ao_num)) + !$OMP DO do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array) + call give_all_aos_and_grad_and_lapl_at_r(final_grid_points(1,i),aos_array,aos_grad_array,aos_lapl_array) do j = 1, ao_num do m = 1, 3 aos_lapl_in_r_array(m,j,i) = aos_lapl_array(m,j) enddo enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO + deallocate( aos_array, aos_grad_array, aos_lapl_array) + !$OMP END PARALLEL END_PROVIDER BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_bis, (n_points_final_grid,ao_num,3)] implicit none BEGIN_DOC -! Transposed gradients -! +! Transposed gradients +! END_DOC integer :: i,j,m double precision :: aos_array(ao_num), r(3) @@ -169,8 +164,8 @@ END_PROVIDER BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_3, (3,n_points_final_grid,ao_num)] implicit none BEGIN_DOC -! Transposed gradients -! +! Transposed gradients +! END_DOC integer :: i,j,m double precision :: aos_array(ao_num), r(3) @@ -187,22 +182,14 @@ END_PROVIDER BEGIN_PROVIDER[double precision, aos_in_r_array_extra, (ao_num,n_points_extra_final_grid)] implicit none BEGIN_DOC - ! aos_in_r_array_extra(i,j) = value of the ith ao on the jth grid point of the EXTRA grid + ! aos_in_r_array_extra(i,j) = value of the ith ao on the jth grid point of the EXTRA grid END_DOC - integer :: i,j - double precision :: aos_array(ao_num), r(3) + integer :: i !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,aos_array,j) & - !$OMP SHARED(aos_in_r_array_extra,n_points_extra_final_grid,ao_num,final_grid_points_extra) + !$OMP DEFAULT (NONE) PRIVATE (i) & + !$OMP SHARED(aos_in_r_array_extra,n_points_extra_final_grid,final_grid_points_extra) do i = 1, n_points_extra_final_grid - r(1) = final_grid_points_extra(1,i) - r(2) = final_grid_points_extra(2,i) - r(3) = final_grid_points_extra(3,i) - call give_all_aos_at_r(r,aos_array) - do j = 1, ao_num - aos_in_r_array_extra(j,i) = aos_array(j) - enddo + call give_all_aos_at_r(final_grid_points_extra(1,i),aos_in_r_array_extra(1,i)) enddo !$OMP END PARALLEL DO @@ -214,9 +201,9 @@ END_PROVIDER BEGIN_PROVIDER[double precision, aos_in_r_array_extra_transp, (n_points_extra_final_grid,ao_num)] BEGIN_DOC - ! aos_in_r_array_extra_transp(i,j) = value of the jth ao on the ith grid point of the EXTRA grid + ! aos_in_r_array_extra_transp(i,j) = value of the jth ao on the ith grid point of the EXTRA grid END_DOC - + implicit none integer :: i, j double precision :: aos_array(ao_num), r(3) @@ -235,27 +222,28 @@ BEGIN_PROVIDER[double precision, aos_grad_in_r_array_extra, (ao_num,n_points_ext implicit none integer :: i, j, m - double precision :: aos_array(ao_num), r(3) - double precision :: aos_grad_array(3,ao_num) + double precision, allocatable :: aos_array(:), aos_grad_array(:,:) - !$OMP PARALLEL DO & + + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) & + !$OMP PRIVATE (i,j,m,aos_array,aos_grad_array) & !$OMP SHARED(aos_grad_in_r_array_extra,n_points_extra_final_grid,ao_num,final_grid_points_extra) + allocate(aos_array(ao_num), aos_grad_array(3,ao_num)) + !$OMP DO do i = 1, n_points_extra_final_grid - r(1) = final_grid_points_extra(1,i) - r(2) = final_grid_points_extra(2,i) - r(3) = final_grid_points_extra(3,i) - call give_all_aos_and_grad_at_r(r, aos_array, aos_grad_array) + call give_all_aos_and_grad_at_r(final_grid_points_extra(1,i), aos_array, aos_grad_array) do m = 1, 3 do j = 1, ao_num aos_grad_in_r_array_extra(j,i,m) = aos_grad_array(m,j) enddo enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO + deallocate(aos_array,aos_grad_array) + !$OMP END PARALLEL -END_PROVIDER +END_PROVIDER ! --- diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index 623de4f8..f3cc5559 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -21,20 +21,11 @@ BEGIN_DOC ! mos_in_r_array(i,j) = value of the ith mo on the jth grid point END_DOC - integer :: i,j - double precision :: mos_array(mo_num), r(3) - !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,mos_array,j) & + integer :: i + !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE (i) & !$OMP SHARED(mos_in_r_array_omp,n_points_final_grid,mo_num,final_grid_points) do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_mos_at_r(r,mos_array) - do j = 1, mo_num - mos_in_r_array_omp(j,i) = mos_array(j) - enddo + call give_all_mos_at_r(final_grid_points(1,i),mos_in_r_array_omp(1,i)) enddo !$OMP END PARALLEL DO END_PROVIDER diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f index 88dad8c3..ba3675d4 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -22,22 +22,32 @@ endif do istate = 1, N_states - do ipoint = 1, n_points_final_grid if(mu_of_r_potential.EQ."hf")then - mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) + enddo else if(mu_of_r_potential.EQ."hf_old")then - mu_of_r_prov(ipoint,istate) = mu_of_r_hf_old(ipoint) + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_hf_old(ipoint) + enddo else if(mu_of_r_potential.EQ."hf_sparse")then - mu_of_r_prov(ipoint,istate) = mu_of_r_hf_sparse(ipoint) + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_hf_sparse(ipoint) + enddo else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then - mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) + enddo + else if(mu_of_r_potential.EQ."proj")then + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_projector_mo(ipoint) + enddo else print*,'you requested the following mu_of_r_potential' print*,mu_of_r_potential print*,'which does not correspond to any of the options for such keyword' stop endif - enddo enddo if (write_mu_of_r) then @@ -225,3 +235,66 @@ enddo END_PROVIDER + +BEGIN_PROVIDER [double precision, mu_of_r_projector_mo, (n_points_final_grid) ] + implicit none + BEGIN_DOC + ! mu(r) computed with the projector onto the atomic basis + ! P_B(\mathbf{r},\mathbf{r}') = \sum_{ij} | + ! \chi_{i} \rangle \left[S^{-1}\right]_{ij} \langle \chi_{j} | + ! \] where $i$ and $j$ denote all atomic orbitals. + END_DOC + + double precision, parameter :: factor = dsqrt(2.d0*dacos(-1.d0)) + double precision, allocatable :: tmp(:,:) + integer :: ipoint + + + do ipoint=1,n_points_final_grid + mu_of_r_projector_mo(ipoint) = 0.d0 + integer :: i,j + do j=1,n_inact_act_orb + i = list_inact_act(j) + mu_of_r_projector_mo(ipoint) = mu_of_r_projector_mo(ipoint) + & + mos_in_r_array_omp(i,ipoint) * mos_in_r_array_omp(i,ipoint) + enddo + do j=1,n_virt_orb + i = list_virt(j) + mu_of_r_projector_mo(ipoint) = mu_of_r_projector_mo(ipoint) + & + mos_in_r_array_omp(i,ipoint) * mos_in_r_array_omp(i,ipoint) + enddo + enddo + + do ipoint=1,n_points_final_grid + ! epsilon + mu_of_r_projector_mo(ipoint) = 1.d0/(2.d0*dacos(-1.d0) * mu_of_r_projector_mo(ipoint)**(2.d0/3.d0)) + ! mu + mu_of_r_projector_mo(ipoint) = 1.d0/dsqrt( 2.d0*mu_of_r_projector_mo(ipoint) ) + enddo +END_PROVIDER + + + +BEGIN_PROVIDER [double precision, mu_average_proj, (N_states)] + implicit none + BEGIN_DOC + ! average value of mu(r) weighted with the total one-e density and divided by the number of electrons + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the one- and two-body density matrix are excluded + END_DOC + integer :: ipoint,istate + double precision :: weight,density + do istate = 1, N_states + mu_average_proj(istate) = 0.d0 + do ipoint = 1, n_points_final_grid + weight =final_weight_at_r_vector(ipoint) + density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) & + + one_e_dm_and_grad_beta_in_r(4,ipoint,istate) + mu_average_proj(istate) += mu_of_r_projector_mo(ipoint) * weight * density + enddo + mu_average_proj(istate) = mu_average_proj(istate) / elec_num_grid_becke(istate) + enddo +END_PROVIDER + diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f index 9e2ea018..6cad72be 100644 --- a/src/two_body_rdm/act_2_rdm.irp.f +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -145,6 +145,7 @@ print*,'' print*,'Providing act_2_rdm_spin_trace_mo ' character*(128) :: name_file + PROVIDE all_mo_integrals name_file = 'act_2_rdm_spin_trace_mo' ispin = 4 act_2_rdm_spin_trace_mo = 0.d0 diff --git a/src/two_rdm_routines/davidson_like_2rdm.irp.f b/src/two_rdm_routines/davidson_like_2rdm.irp.f index 09436663..d92f1924 100644 --- a/src/two_rdm_routines/davidson_like_2rdm.irp.f +++ b/src/two_rdm_routines/davidson_like_2rdm.irp.f @@ -13,7 +13,7 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz END_DOC integer, intent(in) :: N_st,sze integer, intent(in) :: dim1,norb,list_orb(norb),ispin - double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st) double precision, intent(in) :: u_0(sze,N_st) integer :: k @@ -50,7 +50,7 @@ subroutine orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_ END_DOC integer, intent(in) :: N_st,sze,istart,iend,ishift,istep integer, intent(in) :: dim1,norb,list_orb(norb),ispin - double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st) double precision, intent(in) :: u_t(N_st,N_det) integer :: k @@ -91,7 +91,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin integer, intent(in) :: N_st,sze,istart,iend,ishift,istep double precision, intent(in) :: u_t(N_st,N_det) integer, intent(in) :: dim1,norb,list_orb(norb),ispin - double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st) integer(omp_lock_kind) :: lock_2rdm integer :: i,j,k,l @@ -139,6 +139,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 + sze_buff = sze_buff*100 list_orb_reverse = -1000 do i = 1, norb list_orb_reverse(list_orb(i)) = i @@ -154,6 +155,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin ! Prepare the array of all alpha single excitations ! ------------------------------------------------- + double precision, allocatable :: big_array_local(:,:,:,:,:) + PROVIDE N_int nthreads_davidson elec_alpha_num !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) & !$OMP SHARED(psi_bilinear_matrix_rows, N_det,lock_2rdm,& @@ -173,7 +176,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin !$OMP buffer, doubles, n_doubles, & !$OMP tmp_det2, idx, l, kcol_prev, & !$OMP singles_a, n_singles_a, singles_b, & - !$OMP n_singles_b, nkeys, keys, values) + !$OMP n_singles_b, nkeys, keys, values, big_array_local) ! Alpha/Beta double excitations ! ============================= @@ -184,6 +187,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin singles_b(maxab), & doubles(maxab), & idx(maxab)) + allocate( big_array_local(N_states,dim1, dim1, dim1, dim1) ) + big_array_local(:,:,:,:,:) = 0.d0 kcol_prev=-1 @@ -191,8 +196,9 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin ASSERT (istart > 0) ASSERT (istep > 0) - !$OMP DO SCHEDULE(dynamic,64) + !$OMP DO SCHEDULE(dynamic) do k_a=istart+ishift,iend,istep +!print *, 'aa', k_a, '/', iend krow = psi_bilinear_matrix_rows(k_a) ASSERT (krow <= N_det_alpha_unique) @@ -254,33 +260,36 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin do l= 1, N_states c_1(l) = u_t(l,l_a) * u_t(l,k_a) enddo - if(alpha_beta)then - ! only ONE contribution - if (nkeys+1 .ge. sze_buff) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 - endif - else if (spin_trace)then - ! TWO contributions +! if(alpha_beta)then +! ! only ONE contribution +! if (nkeys+1 .ge. sze_buff) then +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 +! endif +! else if (spin_trace)then +! ! TWO contributions if (nkeys+2 .ge. sze_buff) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) nkeys = 0 endif - endif +! endif call orb_range_off_diag_double_to_all_states_ab_dm_buffer(tmp_det,tmp_det2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) enddo endif - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 enddo enddo - !$OMP END DO + !$OMP END DO NOWAIT +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) + nkeys = 0 - !$OMP DO SCHEDULE(dynamic,64) + !$OMP DO SCHEDULE(dynamic) do k_a=istart+ishift,iend,istep +!print *, 'ab', k_a, '/', iend ! Single and double alpha exitations @@ -331,36 +340,39 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin ! ---------------------------------- tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - do i=1,n_singles_a - l_a = singles_a(i) - ASSERT (l_a <= N_det) + if(alpha_beta.or.spin_trace.or.alpha_alpha)then + do i=1,n_singles_a + l_a = singles_a(i) + ASSERT (l_a <= N_det) - lrow = psi_bilinear_matrix_rows(l_a) - ASSERT (lrow <= N_det_alpha_unique) + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) + do l= 1, N_states + c_1(l) = u_t(l,l_a) * u_t(l,k_a) + enddo + + ! increment the alpha/beta part for single excitations + if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + ! increment the alpha/alpha part for single excitations + if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_aa_dm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) - do l= 1, N_states - c_1(l) = u_t(l,l_a) * u_t(l,k_a) enddo - if(alpha_beta.or.spin_trace.or.alpha_alpha)then - ! increment the alpha/beta part for single excitations - if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - ! increment the alpha/alpha part for single excitations - if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_single_to_all_states_aa_dm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - endif + endif - enddo - - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 ! Compute Hij for all alpha doubles ! ---------------------------------- @@ -377,14 +389,15 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin c_1(l) = u_t(l,l_a) * u_t(l,k_a) enddo if (nkeys+4 .ge. sze_buff) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) nkeys = 0 endif call orb_range_off_diag_double_to_all_states_aa_dm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) enddo endif - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 ! Single and double beta excitations @@ -432,35 +445,39 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin ! ---------------------------------- tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) - do i=1,n_singles_b - l_b = singles_b(i) - ASSERT (l_b <= N_det) + if(alpha_beta.or.spin_trace.or.beta_beta)then + do i=1,n_singles_b + l_b = singles_b(i) + ASSERT (l_b <= N_det) - lcol = psi_bilinear_matrix_transp_columns(l_b) - ASSERT (lcol <= N_det_beta_unique) + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) + l_a = psi_bilinear_matrix_transp_order(l_b) + do l= 1, N_states + c_1(l) = u_t(l,l_a) * u_t(l,k_a) + enddo - tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) - l_a = psi_bilinear_matrix_transp_order(l_b) - do l= 1, N_states - c_1(l) = u_t(l,l_a) * u_t(l,k_a) - enddo - if(alpha_beta.or.spin_trace.or.beta_beta)then ! increment the alpha/beta part for single excitations if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) nkeys = 0 endif call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) ! increment the beta /beta part for single excitations if (nkeys+4 * elec_alpha_num .ge. sze_buff) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) nkeys = 0 endif call orb_range_off_diag_single_to_all_states_bb_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - endif - enddo - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 + enddo + endif + +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 ! Compute Hij for all beta doubles ! ---------------------------------- @@ -478,7 +495,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin c_1(l) = u_t(l,l_a) * u_t(l,k_a) enddo if (nkeys+4 .ge. sze_buff) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) nkeys = 0 endif call orb_range_off_diag_double_to_all_states_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) @@ -487,8 +505,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin enddo endif - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 ! Diagonal contribution @@ -514,16 +532,28 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin c_1(l) = u_t(l,k_a) * u_t(l,k_a) enddo - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 + if (nkeys+elec_alpha_num*elec_alpha_num .ge. sze_buff) then +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) + nkeys = 0 + endif + call orb_range_diag_to_all_states_2_rdm_dm_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) nkeys = 0 end do - !$OMP END DO + !$OMP END DO NOWAIT deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values) - !$OMP END PARALLEL + !$OMP CRITICAL + do i=1,N_states + big_array(:,:,:,:,i) = big_array(:,:,:,:,i) + big_array_local(i,:,:,:,:) + enddo + !$OMP END CRITICAL + deallocate(big_array_local) + !$OMP END PARALLEL end @@ -550,22 +580,66 @@ subroutine update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,loc integer :: istate integer :: i,h1,h2,p1,p2 - call omp_set_lock(lock_2rdm) + integer, allocatable :: iorder(:) + integer*8, allocatable :: to_sort(:) + + allocate(iorder(nkeys)) + do i=1,nkeys + iorder(i) = i + enddo + + ! If the lock is already taken, sort the keys while waiting for a faster access + if (.not.omp_test_lock(lock_2rdm)) then + allocate(to_sort(nkeys)) + do i=1,nkeys + h1 = keys(1,iorder(i)) + h2 = keys(2,iorder(i))-1 + p1 = keys(3,iorder(i))-1 + p2 = keys(4,iorder(i))-1 + to_sort(i) = int(h1,8) + int(dim1,8)*(int(h2,8) + int(dim1,8)*(int(p1,8) + int(dim1,8)*int(p2,8))) + enddo + call i8sort(to_sort, iorder, nkeys) + deallocate(to_sort) + call omp_set_lock(lock_2rdm) + endif ! print*,'*************' ! print*,'updating' ! print*,'nkeys',nkeys + do istate = 1, N_st + do i = 1, nkeys + h1 = keys(1,iorder(i)) + h2 = keys(2,iorder(i)) + p1 = keys(3,iorder(i)) + p2 = keys(4,iorder(i)) + big_array(h1,h2,p1,p2,istate) = big_array(h1,h2,p1,p2,istate) + values(istate,iorder(i)) + enddo + enddo + call omp_unset_lock(lock_2rdm) + deallocate(iorder) + +end + +subroutine update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local) + use omp_lib + implicit none + integer, intent(in) :: n_st,nkeys,dim1 + integer, intent(in) :: keys(4,nkeys) + double precision, intent(in) :: values(n_st,nkeys) + double precision, intent(inout) :: big_array_local(n_st,dim1,dim1,dim1,dim1) + + integer :: istate + integer :: i,h1,h2,p1,p2 + do i = 1, nkeys h1 = keys(1,i) h2 = keys(2,i) p1 = keys(3,i) p2 = keys(4,i) do istate = 1, N_st -! print*,h1,h2,p1,p2,values(istate,i) - big_array(h1,h2,p1,p2,istate) += values(istate,i) + big_array_local(istate,h1,h2,p1,p2) = big_array_local(istate,h1,h2,p1,p2) + values(istate,i) enddo enddo - call omp_unset_lock(lock_2rdm) end diff --git a/src/two_rdm_routines/update_rdm.irp.f b/src/two_rdm_routines/update_rdm.irp.f index 8aeb0699..f01706f0 100644 --- a/src/two_rdm_routines/update_rdm.irp.f +++ b/src/two_rdm_routines/update_rdm.irp.f @@ -1,762 +1,897 @@ - subroutine orb_range_diag_to_all_states_2_rdm_dm_buffer(det_1,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) +logical function is_integer_in_string_local(orb,bitmask,Nint) + use bitmasks + implicit none + integer, intent(in) :: orb, Nint + integer(bit_kind), intent(in) :: bitmask(Nint) + integer :: j, k + k = ishft(orb-1,-bit_kind_shift)+1 + j = orb-ishft(k-1,bit_kind_shift)-1 + is_integer_in_string_local = iand(bitmask(k), ibset(0_bit_kind, j)) /= 0_bit_kind +end + +subroutine orb_range_diag_to_all_states_2_rdm_dm_buffer(det_1,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) use bitmasks BEGIN_DOC - ! routine that update the DIAGONAL PART of the two body rdms in a specific range of orbitals for a given determinant det_1 - ! - ! c_1 is the array of the contributions to the rdm for all states - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - END_DOC - implicit none - integer, intent(in) :: ispin,sze_buff,N_st - integer, intent(in) :: list_orb_reverse(mo_num) - integer(bit_kind), intent(in) :: det_1(N_int,2) - integer(bit_kind), intent(in) :: orb_bitmask(N_int) - double precision, intent(in) :: c_1(N_st) - double precision, intent(out) :: values(N_st,sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2) - integer :: i,j,h1,h2 - integer(bit_kind) :: det_1_act(N_int,2) - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - do i = 1, N_int - det_1_act(i,1) = iand(det_1(i,1),orb_bitmask(i)) - det_1_act(i,2) = iand(det_1(i,2),orb_bitmask(i)) - enddo - - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - call bitstring_to_list_ab(det_1_act, occ, n_occ_ab, N_int) - logical :: is_integer_in_string - integer :: i1,i2,istate - if(alpha_beta)then - do i = 1, n_occ_ab(1) - i1 = occ(i,1) - do j = 1, n_occ_ab(2) - i2 = occ(j,2) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - ! If alpha/beta, electron 1 is alpha, electron 2 is beta - ! Therefore you don't necessayr have symmetry between electron 1 and 2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - - else if (alpha_alpha)then - do i = 1, n_occ_ab(1) - i1 = occ(i,1) - do j = 1, n_occ_ab(1) - i2 = occ(j,1) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = -0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - else if (beta_beta)then - do i = 1, n_occ_ab(2) - i1 = occ(i,2) - do j = 1, n_occ_ab(2) - i2 = occ(j,2) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = -0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - else if(spin_trace)then - ! 0.5 * (alpha beta + beta alpha) - do i = 1, n_occ_ab(1) - i1 = occ(i,1) - do j = 1, n_occ_ab(2) - i2 = occ(j,2) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - do i = 1, n_occ_ab(1) - i1 = occ(i,1) - do j = 1, n_occ_ab(1) - i2 = occ(j,1) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = -0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - do i = 1, n_occ_ab(2) - i1 = occ(i,2) - do j = 1, n_occ_ab(2) - i2 = occ(j,2) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = -0.5d0 * c_1(istate) - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - endif - end - - - subroutine orb_range_off_diag_double_to_all_states_ab_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - use bitmasks - BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for -! -! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another -! -! c_1 is the array of the contributions to the rdm for all states -! -! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals -! -! ispin determines which spin-spin component of the two-rdm you will update -! -! ispin == 1 :: alpha/ alpha -! ispin == 2 :: beta / beta -! ispin == 3 :: alpha/ beta -! ispin == 4 :: spin traced <=> total two-rdm -! -! here, only ispin == 3 or 4 will do something - END_DOC - implicit none - integer, intent(in) :: ispin,sze_buff,N_st - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - integer, intent(in) :: list_orb_reverse(mo_num) - double precision, intent(in) :: c_1(N_st) - double precision, intent(out) :: values(N_st,sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - integer :: i,j,h1,h2,p1,p2,istate - integer :: exc(0:2,2,2) - double precision :: phase - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - call get_double_excitation(det_1,det_2,exc,phase,N_int) - h1 = exc(1,1,1) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - h2 = exc(1,1,2) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - p1 = exc(1,2,1) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - p2 = exc(1,2,2) - if(list_orb_reverse(p2).lt.0)return - p2 = list_orb_reverse(p2) - if(alpha_beta)then - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - else if(spin_trace)then - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - endif - end - - subroutine orb_range_off_diag_single_to_all_states_ab_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - use bitmasks - BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 3 or 4 will do something + ! routine that update the DIAGONAL PART of the two body rdms in a specific range of orbitals for a given determinant det_1 + ! + ! c_1 is the array of the contributions to the rdm for all states + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm END_DOC implicit none - integer, intent(in) :: ispin,sze_buff,N_st - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - integer, intent(in) :: list_orb_reverse(mo_num) - integer(bit_kind), intent(in) :: orb_bitmask(N_int) - double precision, intent(in) :: c_1(N_st) - double precision, intent(out) :: values(N_st,sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - + integer, intent(in) :: ispin,sze_buff,N_st + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: det_1(N_int,2) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st) + double precision, intent(out) :: values(N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout) :: nkeys + integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,p1,istate - integer :: exc(0:2,2,2) - double precision :: phase - + integer :: i,j,h1,h2 + integer(bit_kind) :: det_1_act(N_int,2) logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. + do i = 1, N_int + det_1_act(i,1) = iand(det_1(i,1),orb_bitmask(i)) + det_1_act(i,2) = iand(det_1(i,2),orb_bitmask(i)) + enddo + + alpha_alpha = .False. beta_beta = .False. alpha_beta = .False. spin_trace = .False. if( ispin == 1)then - alpha_alpha = .True. + alpha_alpha = .True. else if(ispin == 2)then - beta_beta = .True. + beta_beta = .True. else if(ispin == 3)then - alpha_beta = .True. + alpha_beta = .True. else if(ispin == 4)then - spin_trace = .True. + spin_trace = .True. endif - + + call bitstring_to_list_ab(det_1_act, occ, n_occ_ab, N_int) + logical :: is_integer_in_string_local + integer :: i1,i2,istate + if(alpha_beta)then + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + h1 = list_orb_reverse(i1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h2 = list_orb_reverse(i2) + ! If alpha/beta, electron 1 is alpha, electron 2 is beta + ! Therefore you don't necessayr have symmetry between electron 1 and 2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + + else if (alpha_alpha)then + + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + h1 = list_orb_reverse(i1) + do j = 1, n_occ_ab(1) + i2 = occ(j,1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = -0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + + else if (beta_beta)then + + do i = 1, n_occ_ab(2) + i1 = occ(i,2) + h1 = list_orb_reverse(i1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h2 = list_orb_reverse(i2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = -0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + + else if(spin_trace)then + + ! 0.5 * (alpha beta + beta alpha) + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + h1 = list_orb_reverse(i1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h2 = list_orb_reverse(i2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + do j = 1, n_occ_ab(1) + i2 = occ(j,1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = -0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + do i = 1, n_occ_ab(2) + i1 = occ(i,2) + h1 = list_orb_reverse(i1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h2 = list_orb_reverse(i2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = -0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + endif +end + + +subroutine orb_range_off_diag_double_to_all_states_ab_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that updates the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another + ! + ! c_1 is the array of the contributions to the rdm for all states + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 3 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1(N_st) + double precision, intent(out) :: values(N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout) :: nkeys + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string_local + if (ispin <= 2) return + +! alpha_alpha = .False. +! beta_beta = .False. +! alpha_beta = .False. +! spin_trace = .False. +! if( ispin == 1)then +! alpha_alpha = .True. +! else if(ispin == 2)then +! beta_beta = .True. +! else if(ispin == 3)then +! alpha_beta = .True. +! else if(ispin == 4)then +! spin_trace = .True. +! endif + call get_double_excitation(det_1,det_2,exc,phase,N_int) + + h1 = exc(1,1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 = exc(1,1,2) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 = exc(1,2,1) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 = exc(1,2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) +! if(alpha_beta)then + nkeys += 1 + phase = phase * 0.5d0 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 +! else if(spin_trace)then +! nkeys += 1 +! do istate = 1, N_st +! values(istate,nkeys) = 0.5d0 * c_1(istate) * phase +! enddo +! keys(1,nkeys) = h1 +! keys(2,nkeys) = h2 +! keys(3,nkeys) = p1 +! keys(4,nkeys) = p2 +! nkeys += 1 +! do istate = 1, N_st +! values(istate,nkeys) = 0.5d0 * c_1(istate) * phase +! enddo +! keys(1,nkeys) = h2 +! keys(2,nkeys) = h1 +! keys(3,nkeys) = p2 +! keys(4,nkeys) = p1 +! endif +end + +subroutine orb_range_off_diag_single_to_all_states_ab_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that updates the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 3 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st) + double precision, intent(out) :: values(N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout) :: nkeys + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,p1,istate + integer :: exc(0:2,2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string_local + if (ispin <= 2) return + + alpha_beta = .False. + spin_trace = .False. + if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) call get_single_excitation(det_1,det_2,exc,phase,N_int) if(alpha_beta)then - if (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,1) - if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(2) - h2 = occ(i,2) - if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle - h2 = list_orb_reverse(h2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase + + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return + p1 = exc(1,2,1) + if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return + + h1 = list_orb_reverse(h1) + p1 = list_orb_reverse(p1) + + phase = 0.5d0 * phase + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase + else + ! Mono beta + h1 = exc(1,1,2) + if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return + p1 = exc(1,2,2) + if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return + + h1 = list_orb_reverse(h1) + p1 = list_orb_reverse(p1) + phase = 0.5d0 * phase + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - enddo - else - ! Mono beta - h1 = exc(1,1,2) - if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,2) - if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(1) - h2 = occ(i,1) - if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle - h2 = list_orb_reverse(h2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - enddo - endif + endif + else if(spin_trace)then - if (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,1) - if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(2) - h2 = occ(i,2) - if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle - h2 = list_orb_reverse(h2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - enddo - else - ! Mono beta - h1 = exc(1,1,2) - if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,2) - if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(1) - h2 = occ(i,1) - if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle - h2 = list_orb_reverse(h2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - enddo - endif - endif - end - subroutine orb_range_off_diag_single_to_all_states_aa_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return + p1 = exc(1,2,1) + if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return + + h1 = list_orb_reverse(h1) + p1 = list_orb_reverse(p1) + phase = 0.5d0 * phase + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + + else + + ! Mono beta + h1 = exc(1,1,2) + if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return + p1 = exc(1,2,2) + if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return + + h1 = list_orb_reverse(h1) + p1 = list_orb_reverse(p1) + + phase = 0.5d0 * phase + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + endif + + endif +end + +subroutine orb_range_off_diag_single_to_all_states_aa_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 1 or 4 will do something + ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 1 or 4 will do something END_DOC use bitmasks implicit none - integer, intent(in) :: ispin,sze_buff,N_st - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - integer, intent(in) :: list_orb_reverse(mo_num) - integer(bit_kind), intent(in) :: orb_bitmask(N_int) - double precision, intent(in) :: c_1(N_st) - double precision, intent(out) :: values(N_st,sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st) + double precision, intent(out) :: values(N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout) :: nkeys + integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,p1,istate + integer :: i,j,h1,h2,p1,istate integer :: exc(0:2,2,2) double precision :: phase - + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. + logical :: is_integer_in_string_local + alpha_alpha = .False. +! beta_beta = .False. +! alpha_beta = .False. spin_trace = .False. if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. + alpha_alpha = .True. else if(ispin == 4)then - spin_trace = .True. - endif - - call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) - call get_single_excitation(det_1,det_2,exc,phase,N_int) - if(alpha_alpha.or.spin_trace)then - if (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,1) - if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(1) - h2 = occ(i,1) - if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle - h2 = list_orb_reverse(h2) - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - enddo - else + spin_trace = .True. + else return - endif endif - end - subroutine orb_range_off_diag_single_to_all_states_bb_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) +! if(alpha_alpha.or.spin_trace)then + call get_single_excitation(det_1,det_2,exc,phase,N_int) + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return + p1 = exc(1,2,1) + if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return + + call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) + + h1 = list_orb_reverse(h1) + p1 = list_orb_reverse(p1) + + phase = 0.5d0*phase + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = - c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = - c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + else + return + endif +! endif +end + +subroutine orb_range_off_diag_single_to_all_states_bb_dm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) use bitmasks BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 2 or 4 will do something + ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 2 or 4 will do something END_DOC implicit none - integer, intent(in) :: ispin,sze_buff,N_st - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - integer, intent(in) :: list_orb_reverse(mo_num) - integer(bit_kind), intent(in) :: orb_bitmask(N_int) - double precision, intent(in) :: c_1(N_st) - double precision, intent(out) :: values(N_st,sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st) + double precision, intent(out) :: values(N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout) :: nkeys + integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,p1,istate + integer :: i,j,h1,h2,p1,istate integer :: exc(0:2,2,2) double precision :: phase logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. + logical :: is_integer_in_string_local +! alpha_alpha = .False. beta_beta = .False. - alpha_beta = .False. +! alpha_beta = .False. spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. + if(ispin == 2)then + beta_beta = .True. else if(ispin == 4)then - spin_trace = .True. - endif - - - call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) - call get_single_excitation(det_1,det_2,exc,phase,N_int) - if(beta_beta.or.spin_trace)then - if (exc(0,1,1) == 1) then + spin_trace = .True. + else return - else - ! Mono beta - h1 = exc(1,1,2) - if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,2) - if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(2) - h2 = occ(i,2) - if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle - h2 = list_orb_reverse(h2) - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - enddo - endif endif - end - subroutine orb_range_off_diag_double_to_all_states_aa_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) +! if(beta_beta.or.spin_trace)then + call get_single_excitation(det_1,det_2,exc,phase,N_int) + if (exc(0,1,1) == 1) then + return + else + ! Mono beta + h1 = exc(1,1,2) + if(.not.is_integer_in_string_local(h1,orb_bitmask,N_int))return + p1 = exc(1,2,2) + if(.not.is_integer_in_string_local(p1,orb_bitmask,N_int))return + + call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) + h1 = list_orb_reverse(h1) + p1 = list_orb_reverse(p1) + + phase = 0.5d0*phase + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(.not.is_integer_in_string_local(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = - c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = - c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + endif +! endif +end + + +subroutine orb_range_off_diag_double_to_all_states_aa_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) use bitmasks BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 1 or 4 will do something + ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 1 or 4 will do something END_DOC implicit none - integer, intent(in) :: ispin,sze_buff,N_st - integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) - integer, intent(in) :: list_orb_reverse(mo_num) - double precision, intent(in) :: c_1(N_st) - double precision, intent(out) :: values(N_st,sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - - - integer :: i,j,h1,h2,p1,p2,istate + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1(N_st) + double precision, intent(out) :: values(N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout) :: nkeys + + + integer :: i,j,h1,h2,p1,p2,istate integer :: exc(0:2,2) double precision :: phase - + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. + logical :: is_integer_in_string_local + alpha_alpha = .False. +! beta_beta = .False. +! alpha_beta = .False. spin_trace = .False. if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. + alpha_alpha = .True. else if(ispin == 4)then - spin_trace = .True. + spin_trace = .True. + else + return endif + call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) + h1 =exc(1,1) + if(list_orb_reverse(h1).lt.0)return + h2 =exc(2,1) + if(list_orb_reverse(h2).lt.0)return + p1 =exc(1,2) + if(list_orb_reverse(p1).lt.0)return + p2 =exc(2,2) + if(list_orb_reverse(p2).lt.0)return + + h1 = list_orb_reverse(h1) + h2 = list_orb_reverse(h2) + p1 = list_orb_reverse(p1) + p2 = list_orb_reverse(p2) + + phase = 0.5d0*phase +! if(alpha_alpha.or.spin_trace)then + nkeys += 1 + + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = - c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = - c_1(istate) * phase + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 +! endif +end + +subroutine orb_range_off_diag_double_to_all_states_bb_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a BETA /BETA DOUBLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 2 or 4 will do something + END_DOC + implicit none + + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1(N_st) + double precision, intent(out) :: values(N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout) :: nkeys + + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string_local +! alpha_alpha = .False. + beta_beta = .False. +! alpha_beta = .False. + spin_trace = .False. + if(ispin == 2)then + beta_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + else + return + endif + call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) h1 =exc(1,1) if(list_orb_reverse(h1).lt.0)return @@ -770,11 +905,12 @@ p2 =exc(2,2) if(list_orb_reverse(p2).lt.0)return p2 = list_orb_reverse(p2) - if(alpha_alpha.or.spin_trace)then + +! if(beta_beta.or.spin_trace)then + phase = 0.5d0*phase nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase + values(istate,nkeys) = c_1(istate) * phase enddo keys(1,nkeys) = h1 keys(2,nkeys) = h2 @@ -783,133 +919,30 @@ nkeys += 1 do istate = 1, N_st - values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase + values(istate,nkeys) = - c_1(istate) * phase enddo keys(1,nkeys) = h2 keys(2,nkeys) = h1 keys(3,nkeys) = p1 keys(4,nkeys) = p2 - endif - end - - subroutine orb_range_off_diag_double_to_all_states_bb_dm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - use bitmasks - BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a BETA /BETA DOUBLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 2 or 4 will do something - END_DOC - implicit none - - integer, intent(in) :: ispin,sze_buff,N_st - integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) - integer, intent(in) :: list_orb_reverse(mo_num) - double precision, intent(in) :: c_1(N_st) - double precision, intent(out) :: values(N_st,sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - - integer :: i,j,h1,h2,p1,p2,istate - integer :: exc(0:2,2) - double precision :: phase - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - - call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) - h1 =exc(1,1) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - h2 =exc(2,1) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - p1 =exc(1,2) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - p2 =exc(2,2) - if(list_orb_reverse(p2).lt.0)return - p2 = list_orb_reverse(p2) - if(beta_beta.or.spin_trace)then - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 nkeys += 1 do istate = 1, N_st - values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase + values(istate,nkeys) = - c_1(istate) * phase enddo keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = p2 keys(4,nkeys) = p1 - + nkeys += 1 do istate = 1, N_st - values(istate,nkeys) = 0.5d0 * c_1(istate) * phase + values(istate,nkeys) = c_1(istate) * phase enddo keys(1,nkeys) = h2 keys(2,nkeys) = h1 keys(3,nkeys) = p2 keys(4,nkeys) = p1 - - nkeys += 1 - do istate = 1, N_st - values(istate,nkeys) = - 0.5d0 * c_1(istate) * phase - enddo - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - endif - end +! endif +end