diff --git a/src/determinants/cas_one_e_rdm.irp.f b/src/determinants/cas_one_e_rdm.irp.f new file mode 100644 index 00000000..0471bde6 --- /dev/null +++ b/src/determinants/cas_one_e_rdm.irp.f @@ -0,0 +1,37 @@ + + BEGIN_PROVIDER [double precision, one_e_act_dm_beta_mo_for_dft, (n_act_orb,n_act_orb,N_states)] + implicit none + BEGIN_DOC + ! one_e_act_dm_beta_mo_for_dft = pure ACTIVE part of the ONE ELECTRON REDUCED DENSITY MATRIX for the BETA ELECTRONS + END_DOC + integer :: i,j,ii,jj,istate + do istate = 1, N_states + do ii = 1, n_act_orb + i = list_act(ii) + do jj = 1, n_act_orb + j = list_act(jj) + one_e_act_dm_beta_mo_for_dft(jj,ii,istate) = one_e_dm_mo_beta(j,i,istate) + enddo + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [double precision, one_e_act_dm_alpha_mo_for_dft, (n_act_orb,n_act_orb,N_states)] + implicit none + BEGIN_DOC + ! one_e_act_dm_alpha_mo_for_dft = pure ACTIVE part of the ONE ELECTRON REDUCED DENSITY MATRIX for the ALPHA ELECTRONS + END_DOC + integer :: i,j,ii,jj,istate + do istate = 1, N_states + do ii = 1, n_act_orb + i = list_act(ii) + do jj = 1, n_act_orb + j = list_act(jj) + one_e_act_dm_alpha_mo_for_dft(jj,ii,istate) = one_e_dm_mo_alpha(j,i,istate) + enddo + enddo + enddo + +END_PROVIDER + diff --git a/src/two_body_rdm/NEED b/src/two_body_rdm/NEED index 711fbf96..ca42c679 100644 --- a/src/two_body_rdm/NEED +++ b/src/two_body_rdm/NEED @@ -1 +1,2 @@ davidson_undressed +density_for_dft diff --git a/src/two_body_rdm/all_states_act_2_rdm_dav_routines.irp.f b/src/two_body_rdm/all_states_act_2_rdm_dav_routines.irp.f index 8f40f32a..9d29332e 100644 --- a/src/two_body_rdm/all_states_act_2_rdm_dav_routines.irp.f +++ b/src/two_body_rdm/all_states_act_2_rdm_dav_routines.irp.f @@ -474,6 +474,7 @@ subroutine orb_range_all_states_two_rdm_work_$N_int(big_array,dim1,norb,list_orb c_contrib(l) = c_1(l) * c_1(l) enddo + call orb_range_diagonal_contrib_to_all_two_rdm_dm_all_states(tmp_det,c_contrib,N_st,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) end do diff --git a/src/two_body_rdm/all_states_act_2_rdm_prov.irp.f b/src/two_body_rdm/all_states_act_2_rdm_prov.irp.f index fc6e4224..37a7d3fb 100644 --- a/src/two_body_rdm/all_states_act_2_rdm_prov.irp.f +++ b/src/two_body_rdm/all_states_act_2_rdm_prov.irp.f @@ -3,22 +3,20 @@ BEGIN_PROVIDER [double precision, all_states_act_two_rdm_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none - double precision, allocatable :: state_weights(:) BEGIN_DOC ! all_states_act_two_rdm_alpha_alpha_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha electrons ! -! +! 1/2 * ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC - allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) integer :: ispin ! condition for alpha/beta spin ispin = 1 all_states_act_two_rdm_alpha_alpha_mo = 0.D0 call orb_range_all_states_two_rdm(all_states_act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + END_PROVIDER BEGIN_PROVIDER [double precision, all_states_act_two_rdm_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] @@ -30,9 +28,6 @@ ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC - double precision, allocatable :: state_weights(:) - allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) integer :: ispin ! condition for alpha/beta spin ispin = 2 @@ -43,16 +38,19 @@ BEGIN_PROVIDER [double precision, all_states_act_two_rdm_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none - double precision, allocatable :: state_weights(:) BEGIN_DOC ! all_states_act_two_rdm_alpha_beta_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons ! ! ! ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" +! +! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is alpha, electron 2 is beta +! +! all_states_act_two_rdm_alpha_beta_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta +! +! Therefore you don't necessayr have symmetry between electron 1 and 2 END_DOC - allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) integer :: ispin ! condition for alpha/beta spin print*,'' @@ -82,16 +80,11 @@ ! ! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l) END_DOC - double precision, allocatable :: state_weights(:) - allocate(state_weights(N_states)) - state_weights = 1.d0/dble(N_states) - integer :: ispin + integer :: ispin,i,j,k,l,istate ! condition for alpha/beta spin ispin = 4 all_states_act_two_rdm_spin_trace_mo = 0.d0 - integer :: i call orb_range_all_states_two_rdm(all_states_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) - END_PROVIDER diff --git a/src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f b/src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f index 7606e353..a42f2d79 100644 --- a/src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f +++ b/src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f @@ -59,7 +59,7 @@ 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. @@ -73,6 +73,7 @@ else if(ispin == 4)then spin_trace = .True. endif +! call debug_det(det_1_act,N_int) call bitstring_to_list_ab(det_1_act, occ, n_occ_ab, N_int) logical :: is_integer_in_string integer :: i1,i2 @@ -84,7 +85,9 @@ i2 = occ(j,2) h1 = list_orb_reverse(i1) h2 = list_orb_reverse(i2) - big_array(h1,h2,h1,h2,istate) += c_1(istate) + ! If alpha/beta, electron 1 is alpha, electron 2 is beta + ! Therefore you don't necessayr have symmetry between electron 1 and 2 + big_array(h1,h2,h1,h2,istate) += 1.0d0 * c_1(istate) enddo enddo enddo @@ -101,6 +104,7 @@ enddo enddo enddo +! pause else if (beta_beta)then do istate = 1, N_st do i = 1, n_occ_ab(2) diff --git a/src/two_body_rdm/all_states_full_2_rdm_prov.irp.f b/src/two_body_rdm/all_states_full_2_rdm_prov.irp.f new file mode 100644 index 00000000..55fa78ca --- /dev/null +++ b/src/two_body_rdm/all_states_full_2_rdm_prov.irp.f @@ -0,0 +1,538 @@ + + BEGIN_PROVIDER [double precision, all_states_full_two_rdm_alpha_beta_mo, (n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,N_states)] + implicit none + all_states_full_two_rdm_alpha_beta_mo = 0.d0 + integer :: i,j,k,l,iorb,jorb,korb,lorb,istate + BEGIN_DOC +! all_states_full_two_rdm_alpha_beta_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" +! +! BUT THE STRUCTURE OF THE TWO-RDM ON THE RANGE OF OCCUPIED MOS (CORE+INACT+ACT) BECAUSE IT CAN BE CONVENIENT FOR SOME APPLICATIONS +! +! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA +! +! all_states_act_two_rdm_alpha_beta_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta +! +! Therefore you don't necessary have symmetry between electron 1 and 2 +! +! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero + END_DOC + all_states_full_two_rdm_alpha_beta_mo = 0.d0 + do istate = 1, N_states + !! PURE ACTIVE PART ALPHA-BETA + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_act_orb + korb = list_act(k) + do l = 1, n_act_orb + lorb = list_act(l) + ! alph beta alph beta + all_states_full_two_rdm_alpha_beta_mo(lorb,korb,jorb,iorb,istate) = & + all_states_act_two_rdm_alpha_beta_mo(l,k,j,i,istate) + enddo + enddo + enddo + enddo + !! BETA ACTIVE - ALPHA inactive + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_inact_orb + korb = list_inact(k) + ! alph beta alph beta + all_states_full_two_rdm_alpha_beta_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate) + enddo + enddo + enddo + + !! ALPHA ACTIVE - BETA inactive + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_inact_orb + korb = list_inact(k) + ! alph beta alph beta + all_states_full_two_rdm_alpha_beta_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate) + enddo + enddo + enddo + + !! ALPHA INACTIVE - BETA INACTIVE + !! + do j = 1, n_inact_orb + jorb = list_inact(j) + do k = 1, n_inact_orb + korb = list_inact(k) + ! alph beta alph beta + all_states_full_two_rdm_alpha_beta_mo(korb,jorb,korb,jorb,istate) = 1.D0 + enddo + enddo + +!!!!!!!!!!!! +!!!!!!!!!!!! if "no_core_density" then you don't put the core part +!!!!!!!!!!!! CAN BE USED + if (.not.no_core_density)then + !! BETA ACTIVE - ALPHA CORE + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_core_orb + korb = list_core(k) + ! alph beta alph beta + all_states_full_two_rdm_alpha_beta_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate) + enddo + enddo + enddo + + !! ALPHA ACTIVE - BETA CORE + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_core_orb + korb = list_core(k) + ! alph beta alph beta + all_states_full_two_rdm_alpha_beta_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate) + enddo + enddo + enddo + + !! ALPHA CORE - BETA CORE + !! + do j = 1, n_core_orb + jorb = list_core(j) + do k = 1, n_core_orb + korb = list_core(k) + ! alph beta alph beta + all_states_full_two_rdm_alpha_beta_mo(korb,jorb,korb,jorb,istate) = 1.D0 + enddo + enddo + endif + + enddo + END_PROVIDER + + + BEGIN_PROVIDER [double precision, all_states_full_two_rdm_alpha_alpha_mo, (n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,N_states)] + implicit none + all_states_full_two_rdm_alpha_alpha_mo = 0.d0 + integer :: i,j,k,l,iorb,jorb,korb,lorb,istate + BEGIN_DOC +! all_states_full_two_rdm_alpha_alpha_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/alpha electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" +! +! BUT THE STRUCTURE OF THE TWO-RDM ON THE FULL RANGE OF MOs IS IMPLEMENTED BECAUSE IT CAN BE CONVENIENT FOR SOME APPLICATIONS +! +! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero + END_DOC + + do istate = 1, N_states + !! PURE ACTIVE PART ALPHA-ALPHA + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_act_orb + korb = list_act(k) + do l = 1, n_act_orb + lorb = list_act(l) + all_states_full_two_rdm_alpha_alpha_mo(lorb,korb,jorb,iorb,istate) = & + all_states_act_two_rdm_alpha_alpha_mo(l,k,j,i,istate) + enddo + enddo + enddo + enddo + !! ALPHA ACTIVE - ALPHA inactive + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_inact_orb + korb = list_inact(k) + ! 1 2 1 2 : DIRECT TERM + all_states_full_two_rdm_alpha_alpha_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + all_states_full_two_rdm_alpha_alpha_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + ! 1 2 1 2 : EXCHANGE TERM + all_states_full_two_rdm_alpha_alpha_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + all_states_full_two_rdm_alpha_alpha_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + enddo + enddo + enddo + + !! ALPHA INACTIVE - ALPHA INACTIVE + do j = 1, n_inact_orb + jorb = list_inact(j) + do k = 1, n_inact_orb + korb = list_inact(k) + all_states_full_two_rdm_alpha_alpha_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + all_states_full_two_rdm_alpha_alpha_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + enddo + enddo + +!!!!!!!!!! +!!!!!!!!!! if "no_core_density" then you don't put the core part +!!!!!!!!!! CAN BE USED + if (.not.no_core_density)then + !! ALPHA ACTIVE - ALPHA CORE + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_core_orb + korb = list_core(k) + ! 1 2 1 2 : DIRECT TERM + all_states_full_two_rdm_alpha_alpha_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + all_states_full_two_rdm_alpha_alpha_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + ! 1 2 1 2 : EXCHANGE TERM + all_states_full_two_rdm_alpha_alpha_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + all_states_full_two_rdm_alpha_alpha_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + enddo + enddo + enddo + !! ALPHA CORE - ALPHA CORE + + do j = 1, n_core_orb + jorb = list_core(j) + do k = 1, n_core_orb + korb = list_core(k) + all_states_full_two_rdm_alpha_alpha_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + all_states_full_two_rdm_alpha_alpha_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + enddo + enddo + endif + enddo + + END_PROVIDER + + BEGIN_PROVIDER [double precision, all_states_full_two_rdm_beta_beta_mo, (n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,N_states)] + implicit none + all_states_full_two_rdm_beta_beta_mo = 0.d0 + integer :: i,j,k,l,iorb,jorb,korb,lorb,istate + BEGIN_DOC +! all_states_full_two_rdm_beta_beta_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta/beta electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" +! +! BUT THE STRUCTURE OF THE TWO-RDM ON THE FULL RANGE OF MOs IS IMPLEMENTED BECAUSE IT CAN BE CONVENIENT FOR SOME APPLICATIONS +! +! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero + END_DOC + + do istate = 1, N_states + !! PURE ACTIVE PART beta-beta + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_act_orb + korb = list_act(k) + do l = 1, n_act_orb + lorb = list_act(l) + all_states_full_two_rdm_beta_beta_mo(lorb,korb,jorb,iorb,istate) = & + all_states_act_two_rdm_beta_beta_mo(l,k,j,i,istate) + enddo + enddo + enddo + enddo + !! beta ACTIVE - beta inactive + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_inact_orb + korb = list_inact(k) + ! 1 2 1 2 : DIRECT TERM + all_states_full_two_rdm_beta_beta_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + all_states_full_two_rdm_beta_beta_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + ! 1 2 1 2 : EXCHANGE TERM + all_states_full_two_rdm_beta_beta_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + all_states_full_two_rdm_beta_beta_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + enddo + enddo + enddo + + !! beta INACTIVE - beta INACTIVE + do j = 1, n_inact_orb + jorb = list_inact(j) + do k = 1, n_inact_orb + korb = list_inact(k) + all_states_full_two_rdm_beta_beta_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + all_states_full_two_rdm_beta_beta_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + enddo + enddo + +!!!!!!!!!!!! +!!!!!!!!!!!! if "no_core_density" then you don't put the core part +!!!!!!!!!!!! CAN BE USED + if (.not.no_core_density)then + !! beta ACTIVE - beta CORE + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_core_orb + korb = list_core(k) + ! 1 2 1 2 : DIRECT TERM + all_states_full_two_rdm_beta_beta_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + all_states_full_two_rdm_beta_beta_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + ! 1 2 1 2 : EXCHANGE TERM + all_states_full_two_rdm_beta_beta_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + all_states_full_two_rdm_beta_beta_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + enddo + enddo + enddo + !! beta CORE - beta CORE + + do j = 1, n_core_orb + jorb = list_core(j) + do k = 1, n_core_orb + korb = list_core(k) + all_states_full_two_rdm_beta_beta_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + all_states_full_two_rdm_beta_beta_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + enddo + enddo + endif + enddo + + END_PROVIDER + + BEGIN_PROVIDER [double precision, all_states_full_two_rdm_spin_trace_mo, (n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,N_states)] + implicit none + all_states_full_two_rdm_spin_trace_mo = 0.d0 + integer :: i,j,k,l,iorb,jorb,korb,lorb,istate + BEGIN_DOC +! all_states_full_two_rdm_beta_beta_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta/beta electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" +! +! BUT THE STRUCTURE OF THE TWO-RDM ON THE FULL RANGE OF MOs IS IMPLEMENTED BECAUSE IT CAN BE CONVENIENT FOR SOME APPLICATIONS +! +! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero + END_DOC + + do istate = 1, N_states + !!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!! + !! PURE ACTIVE PART SPIN-TRACE + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_act_orb + korb = list_act(k) + do l = 1, n_act_orb + lorb = list_act(l) + all_states_full_two_rdm_spin_trace_mo(lorb,korb,jorb,iorb,istate) += & + all_states_act_two_rdm_spin_trace_mo(l,k,j,i,istate) + enddo + enddo + enddo + enddo + + !!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!! + !!!!! BETA-BETA !!!!! + !! beta ACTIVE - beta inactive + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_inact_orb + korb = list_inact(k) + ! 1 2 1 2 : DIRECT TERM + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + all_states_full_two_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + ! 1 2 1 2 : EXCHANGE TERM + all_states_full_two_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + all_states_full_two_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + enddo + enddo + enddo + !! beta INACTIVE - beta INACTIVE + do j = 1, n_inact_orb + jorb = list_inact(j) + do k = 1, n_inact_orb + korb = list_inact(k) + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + all_states_full_two_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + enddo + enddo + if (.not.no_core_density)then + !! beta ACTIVE - beta CORE + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_core_orb + korb = list_core(k) + ! 1 2 1 2 : DIRECT TERM + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + all_states_full_two_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + ! 1 2 1 2 : EXCHANGE TERM + all_states_full_two_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + all_states_full_two_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + enddo + enddo + enddo + !! beta CORE - beta CORE + do j = 1, n_core_orb + jorb = list_core(j) + do k = 1, n_core_orb + korb = list_core(k) + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + all_states_full_two_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + enddo + enddo + endif + + !!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!! + !!!!! ALPHA-ALPHA !!!!! + !! ALPHA ACTIVE - ALPHA inactive + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_inact_orb + korb = list_inact(k) + ! 1 2 1 2 : DIRECT TERM + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + all_states_full_two_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + ! 1 2 1 2 : EXCHANGE TERM + all_states_full_two_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + all_states_full_two_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + enddo + enddo + enddo + !! ALPHA INACTIVE - ALPHA INACTIVE + do j = 1, n_inact_orb + jorb = list_inact(j) + do k = 1, n_inact_orb + korb = list_inact(k) + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + all_states_full_two_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + enddo + enddo + if (.not.no_core_density)then + !! ALPHA ACTIVE - ALPHA CORE + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_core_orb + korb = list_core(k) + ! 1 2 1 2 : DIRECT TERM + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + all_states_full_two_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + ! 1 2 1 2 : EXCHANGE TERM + all_states_full_two_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + all_states_full_two_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + enddo + enddo + enddo + !! ALPHA CORE - ALPHA CORE + do j = 1, n_core_orb + jorb = list_core(j) + do k = 1, n_core_orb + korb = list_core(k) + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + all_states_full_two_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 + enddo + enddo + endif + + !!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!! + !!!!! ALPHA-BETA + BETA-ALPHA !!!!! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_inact_orb + korb = list_inact(k) + ! ALPHA INACTIVE - BETA ACTIVE + ! alph beta alph beta + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + ! beta alph beta alph + all_states_full_two_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + ! BETA INACTIVE - ALPHA ACTIVE + ! beta alph beta alpha + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + ! alph beta alph beta + all_states_full_two_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + enddo + enddo + enddo + !! ALPHA INACTIVE - BETA INACTIVE + do j = 1, n_inact_orb + jorb = list_inact(j) + do k = 1, n_inact_orb + korb = list_inact(k) + ! alph beta alph beta + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0 + all_states_full_two_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0 + enddo + enddo + +!!!!!!!!!!!! +!!!!!!!!!!!! if "no_core_density" then you don't put the core part +!!!!!!!!!!!! CAN BE USED + if (.not.no_core_density)then + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_core_orb + korb = list_core(k) + !! BETA ACTIVE - ALPHA CORE + ! alph beta alph beta + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate) + ! beta alph beta alph + all_states_full_two_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate) + !! ALPHA ACTIVE - BETA CORE + ! alph beta alph beta + all_states_full_two_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate) + ! beta alph beta alph + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate) + enddo + enddo + enddo + !! ALPHA CORE - BETA CORE + do j = 1, n_core_orb + jorb = list_core(j) + do k = 1, n_core_orb + korb = list_core(k) + ! alph beta alph beta + all_states_full_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0 + all_states_full_two_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0 + enddo + enddo + + endif + enddo + + END_PROVIDER diff --git a/src/two_body_rdm/compute_orb_range_omp.irp.f b/src/two_body_rdm/compute_orb_range_omp.irp.f new file mode 100644 index 00000000..0ba934d7 --- /dev/null +++ b/src/two_body_rdm/compute_orb_range_omp.irp.f @@ -0,0 +1,807 @@ + subroutine orb_range_diag_to_all_two_rdm_dm_buffer(det_1,c_1,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 supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! 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 + 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 + double precision, intent(out) :: values(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 + 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) + nkeys += 1 + values(nkeys) = c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + 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 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = -0.5d0 * c_1 + 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 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = -0.5d0 * c_1 + 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 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 + 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 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = -0.5d0 * c_1 + 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 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = -0.5d0 * c_1 + 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_two_rdm_ab_dm_buffer(det_1,det_2,c_1,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 supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +! +! 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 + 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 + double precision, intent(out) :: values(sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + integer :: i,j,h1,h2,p1,p2 + 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 + values(nkeys) = c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + else if(spin_trace)then + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = p1 + keys(2,nkeys) = p2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + endif + end + + subroutine orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(det_1,det_2,c_1,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 + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff + 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 + double precision, intent(out) :: values(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 + 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 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(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,1) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + nkeys += 1 + values(nkeys) = c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + else + ! Mono beta + h1 = exc(1,1,2) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,2) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + nkeys += 1 + values(nkeys) = c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + endif + else if(spin_trace)then + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,1) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + 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(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,2) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + !print*,'****************' + !print*,'****************' + !print*,'h1,p1',h1,p1 + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + ! print*,'h2 = ',h2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + 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_two_rdm_aa_dm_buffer(det_1,det_2,c_1,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 + END_DOC + use bitmasks + implicit none + integer, intent(in) :: ispin,sze_buff + 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 + double precision, intent(out) :: values(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 + 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 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(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,1) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + else + return + endif + endif + end + + subroutine orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(det_1,det_2,c_1,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 + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff + 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 + double precision, intent(out) :: values(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 + 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 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 + return + else + ! Mono beta + h1 = exc(1,1,2) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,2) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + 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_two_rdm_aa_dm_buffer(det_1,det_2,c_1,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 + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff + 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 + double precision, intent(out) :: values(sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + + integer :: i,j,h1,h2,p1,p2 + 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(alpha_alpha.or.spin_trace)then + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + 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_two_rdm_bb_dm_buffer(det_1,det_2,c_1,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 + 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 + double precision, intent(out) :: values(sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: i,j,h1,h2,p1,p2 + 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 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + endif + end + diff --git a/src/two_body_rdm/orb_range_omp.irp.f b/src/two_body_rdm/orb_range_omp.irp.f new file mode 100644 index 00000000..baa26ced --- /dev/null +++ b/src/two_body_rdm/orb_range_omp.irp.f @@ -0,0 +1,85 @@ + + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + BEGIN_DOC +! state_av_act_two_rdm_openmp_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs +! = + END_DOC + allocate(state_weights(N_states)) + state_weights = state_average_weight + integer :: ispin + ! condition for alpha/beta spin + ispin = 1 + state_av_act_two_rdm_openmp_alpha_alpha_mo = 0.D0 + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + BEGIN_DOC +! state_av_act_two_rdm_openmp_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs +! = + END_DOC + allocate(state_weights(N_states)) + state_weights = state_average_weight + integer :: ispin + ! condition for alpha/beta spin + ispin = 2 + state_av_act_two_rdm_openmp_beta_beta_mo = 0.d0 + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_beta_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + BEGIN_DOC +! state_av_act_two_rdm_openmp_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs +! = + END_DOC + allocate(state_weights(N_states)) + state_weights = state_average_weight + integer :: ispin + ! condition for alpha/beta spin + print*,'' + print*,'' + print*,'' + print*,'providint state_av_act_two_rdm_openmp_alpha_beta_mo ' + ispin = 3 + print*,'ispin = ',ispin + state_av_act_two_rdm_openmp_alpha_beta_mo = 0.d0 + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_alpha_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + BEGIN_DOC +! state_av_act_two_rdm_openmp_spin_trace_mo(i,j,k,l) = state average physicist spin trace two-body rdm restricted to the ACTIVE indices +! The active part of the two-electron energy can be computed as: +! +! \sum_{i,j,k,l = 1, n_act_orb} state_av_act_two_rdm_openmp_spin_trace_mo(i,j,k,l) * < ii jj | kk ll > +! +! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l) + END_DOC + double precision, allocatable :: state_weights(:) + allocate(state_weights(N_states)) + state_weights = state_average_weight + integer :: ispin + ! condition for alpha/beta spin + ispin = 4 + state_av_act_two_rdm_openmp_spin_trace_mo = 0.d0 + integer :: i + double precision :: wall_0,wall_1 + call wall_time(wall_0) + print*,'providing the state average TWO-RDM ...' + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + call wall_time(wall_1) + print*,'Time to provide the state average TWO-RDM',wall_1 - wall_0 + END_PROVIDER + diff --git a/src/two_body_rdm/orb_range_routines_omp.irp.f b/src/two_body_rdm/orb_range_routines_omp.irp.f new file mode 100644 index 00000000..b6e59540 --- /dev/null +++ b/src/two_body_rdm/orb_range_routines_omp.irp.f @@ -0,0 +1,568 @@ +subroutine orb_range_two_rdm_state_av_openmp(big_array,dim1,norb,list_orb,state_weights,ispin,u_0,N_st,sze) + use bitmasks + implicit none + BEGIN_DOC + ! if ispin == 1 :: alpha/alpha 2rdm + ! == 2 :: beta /beta 2rdm + ! == 3 :: alpha/beta 2rdm + ! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba)) + ! + ! Assumes that the determinants are in psi_det + ! + ! istart, iend, ishift, istep are used in ZMQ parallelization. + 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(in) :: u_0(sze,N_st),state_weights(N_st) + + integer :: k + double precision, allocatable :: u_t(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t + allocate(u_t(N_st,N_det)) + do k=1,N_st + call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det) + enddo + call dtranspose( & + u_0, & + size(u_0, 1), & + u_t, & + size(u_t, 1), & + N_det, N_st) + + call orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1) + deallocate(u_t) + + do k=1,N_st + call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) + enddo + +end + +subroutine orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + use bitmasks + implicit none + BEGIN_DOC + ! Computes two-rdm + ! + ! Default should be 1,N_det,0,1 + 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(in) :: u_t(N_st,N_det),state_weights(N_st) + + integer :: k + + PROVIDE N_int + + select case (N_int) + case (1) + call orb_range_two_rdm_state_av_openmp_work_1(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case (2) + call orb_range_two_rdm_state_av_openmp_work_2(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case (3) + call orb_range_two_rdm_state_av_openmp_work_3(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case (4) + call orb_range_two_rdm_state_av_openmp_work_4(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case default + call orb_range_two_rdm_state_av_openmp_work_N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + end select +end + + + + + BEGIN_TEMPLATE +subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + use bitmasks + use omp_lib + implicit none + BEGIN_DOC + ! Computes the two rdm for the N_st vectors |u_t> + ! if ispin == 1 :: alpha/alpha 2rdm + ! == 2 :: beta /beta 2rdm + ! == 3 :: alpha/beta 2rdm + ! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba)) + ! The 2rdm will be computed only on the list of orbitals list_orb, which contains norb + ! In any cases, the state average weights will be used with an array state_weights + ! Default should be 1,N_det,0,1 for istart,iend,ishift,istep + END_DOC + integer, intent(in) :: N_st,sze,istart,iend,ishift,istep + double precision, intent(in) :: u_t(N_st,N_det),state_weights(N_st) + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + + integer(omp_lock_kind) :: lock_2rdm + integer :: i,j,k,l + integer :: k_a, k_b, l_a, l_b + integer :: krow, kcol + integer :: lrow, lcol + integer(bit_kind) :: spindet($N_int) + integer(bit_kind) :: tmp_det($N_int,2) + integer(bit_kind) :: tmp_det2($N_int,2) + integer(bit_kind) :: tmp_det3($N_int,2) + integer(bit_kind), allocatable :: buffer(:,:) + integer :: n_doubles + integer, allocatable :: doubles(:) + integer, allocatable :: singles_a(:) + integer, allocatable :: singles_b(:) + integer, allocatable :: idx(:), idx0(:) + integer :: maxab, n_singles_a, n_singles_b, kcol_prev + double precision :: c_average + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + integer(bit_kind) :: orb_bitmask($N_int) + integer :: list_orb_reverse(mo_num) + integer, allocatable :: keys(:,:) + double precision, allocatable :: values(:) + integer :: nkeys,sze_buff + 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. + else + print*,'Wrong parameter for ispin in general_two_rdm_state_av_openmp_work' + print*,'ispin = ',ispin + stop + endif + + + PROVIDE N_int + + call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) + sze_buff = norb ** 3 + 6 * norb + list_orb_reverse = -1000 + do i = 1, norb + list_orb_reverse(list_orb(i)) = i + enddo + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + allocate(idx0(maxab)) + + do i=1,maxab + idx0(i) = i + enddo + call omp_init_lock(lock_2rdm) + + ! Prepare the array of all alpha single excitations + ! ------------------------------------------------- + + 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,& + !$OMP psi_bilinear_matrix_columns, & + !$OMP psi_det_alpha_unique, psi_det_beta_unique,& + !$OMP n_det_alpha_unique, n_det_beta_unique, N_int,& + !$OMP psi_bilinear_matrix_transp_rows, & + !$OMP psi_bilinear_matrix_transp_columns, & + !$OMP psi_bilinear_matrix_transp_order, N_st, & + !$OMP psi_bilinear_matrix_order_transp_reverse, & + !$OMP psi_bilinear_matrix_columns_loc, & + !$OMP psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, & + !$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, state_weights, dim1, & + !$OMP ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask) & + !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, c_2, & + !$OMP lcol, lrow, l_a, l_b, & + !$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, c_average) + + ! Alpha/Beta double excitations + ! ============================= + nkeys = 0 + allocate( keys(4,sze_buff), values(sze_buff)) + allocate( buffer($N_int,maxab), & + singles_a(maxab), & + singles_b(maxab), & + doubles(maxab), & + idx(maxab)) + + kcol_prev=-1 + + ASSERT (iend <= N_det) + ASSERT (istart > 0) + ASSERT (istep > 0) + + !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + if (kcol /= kcol_prev) then + call get_all_spin_singles_$N_int( & + psi_det_beta_unique, idx0, & + tmp_det(1,2), N_det_beta_unique, & + singles_b, n_singles_b) + endif + kcol_prev = kcol + + ! Loop over singly excited beta columns + ! ------------------------------------- + + do i=1,n_singles_b + lcol = singles_b(i) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) + + l_a = psi_bilinear_matrix_columns_loc(lcol) + ASSERT (l_a <= N_det) + + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) + + ASSERT (l_a <= N_det) + idx(j) = l_a + l_a = l_a+1 + enddo + j = j-1 + + call get_all_spin_singles_$N_int( & + buffer, idx, tmp_det(1,1), j, & + singles_a, n_singles_a ) + + ! Loop over alpha singles + ! ----------------------- + + if(alpha_beta.or.spin_trace)then + do k = 1,n_singles_a + l_a = singles_a(k) + ASSERT (l_a <= N_det) + + 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) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + if(alpha_beta)then + ! only ONE contribution + if (nkeys+1 .ge. size(values)) then + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + endif + else if (spin_trace)then + ! TWO contributions + if (nkeys+2 .ge. size(values)) then + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + endif + endif + call orb_range_off_diag_double_to_two_rdm_ab_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + + enddo + endif + + enddo + + enddo + !$OMP END DO + + !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + + ! Single and double alpha exitations + ! =================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + ! Initial determinant is at k_b in beta-major representation + ! ---------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + spindet(1:$N_int) = tmp_det(1:$N_int,1) + + ! Loop inside the beta column to gather all the connected alphas + lcol = psi_bilinear_matrix_columns(k_a) + l_a = psi_bilinear_matrix_columns_loc(lcol) + do i=1,N_det_alpha_unique + if (l_a > N_det) exit + lcol = psi_bilinear_matrix_columns(l_a) + if (lcol /= kcol) exit + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) + idx(i) = l_a + l_a = l_a+1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles_$N_int( & + buffer, idx, spindet, i, & + singles_a, doubles, n_singles_a, n_doubles ) + + ! Compute Hij for all alpha singles + ! ---------------------------------- + + 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) + + 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) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + 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(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,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(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + endif + + enddo + + + ! Compute Hij for all alpha doubles + ! ---------------------------------- + + if(alpha_alpha.or.spin_trace)then + do i=1,n_doubles + l_a = doubles(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + if (nkeys+4 .ge. sze_buff) then + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + endif + call orb_range_off_diag_double_to_two_rdm_aa_dm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + enddo + endif + + + ! Single and double beta excitations + ! ================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + kcol = psi_bilinear_matrix_columns(k_a) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + spindet(1:$N_int) = tmp_det(1:$N_int,2) + + ! Initial determinant is at k_b in beta-major representation + ! ----------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + ! Loop inside the alpha row to gather all the connected betas + lrow = psi_bilinear_matrix_transp_rows(k_b) + l_b = psi_bilinear_matrix_transp_rows_loc(lrow) + do i=1,N_det_beta_unique + if (l_b > N_det) exit + lrow = psi_bilinear_matrix_transp_rows(l_b) + if (lrow /= krow) exit + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol) + idx(i) = l_b + l_b = l_b+1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles_$N_int( & + buffer, idx, spindet, i, & + singles_b, doubles, n_singles_b, n_doubles ) + + ! Compute Hij for all beta singles + ! ---------------------------------- + + 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) + + 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) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + 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(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,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(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + endif + enddo + + ! Compute Hij for all beta doubles + ! ---------------------------------- + + if(beta_beta.or.spin_trace)then + do i=1,n_doubles + l_b = doubles(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + l_a = psi_bilinear_matrix_transp_order(l_b) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + if (nkeys+4 .ge. sze_buff) then + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + endif + call orb_range_off_diag_double_to_two_rdm_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + ASSERT (l_a <= N_det) + + enddo + endif + + + ! Diagonal contribution + ! ===================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + double precision, external :: diag_wee_mat_elem, diag_S_mat_elem + + double precision :: c_1(N_states),c_2(N_states) + c_average = 0.d0 + do l = 1, N_states + c_1(l) = u_t(l,k_a) + c_average += c_1(l) * c_1(l) * state_weights(l) + enddo + + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + call orb_range_diag_to_all_two_rdm_dm_buffer(tmp_det,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + nkeys = 0 + + end do + !$OMP END DO + deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values) + !$OMP END PARALLEL + +end + + SUBST [ N_int ] + + 1;; + 2;; + 3;; + 4;; + N_int;; + + END_TEMPLATE + + +subroutine update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + use omp_lib + implicit none + integer, intent(in) :: nkeys,dim1 + integer, intent(in) :: keys(4,nkeys) + double precision, intent(in) :: values(nkeys) + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + + integer(omp_lock_kind),intent(inout):: lock_2rdm + integer :: i,h1,h2,p1,p2 + call omp_set_lock(lock_2rdm) + do i = 1, nkeys + h1 = keys(1,i) + h2 = keys(2,i) + p1 = keys(3,i) + p2 = keys(4,i) + big_array(h1,h2,p1,p2) += values(i) + enddo + call omp_unset_lock(lock_2rdm) + +end + diff --git a/src/two_body_rdm/test_2_rdm.irp.f b/src/two_body_rdm/test_2_rdm.irp.f new file mode 100644 index 00000000..e993da24 --- /dev/null +++ b/src/two_body_rdm/test_2_rdm.irp.f @@ -0,0 +1,111 @@ +program test_2_rdm + implicit none + read_wf = .True. + touch read_wf + call routine_full_mos + call routine_active_only +end + +subroutine routine_active_only + implicit none + integer :: i,j,k,l,iorb,jorb,korb,lorb,istate + BEGIN_DOC +! This routine computes the two electron repulsion within the active space using various providers +! + END_DOC + + double precision :: vijkl,rdmaa,get_two_e_integral,rdmab,rdmbb,rdmtot + double precision :: accu_aa(N_states),accu_bb(N_states),accu_ab(N_states),accu_tot(N_states) + accu_aa = 0.d0 + accu_ab = 0.d0 + accu_bb = 0.d0 + accu_tot = 0.d0 + do istate = 1, N_states + !! PURE ACTIVE PART + !! + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_act_orb + jorb = list_act(j) + do k = 1, n_act_orb + korb = list_act(k) + do l = 1, n_act_orb + lorb = list_act(l) + + vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + + rdmaa = all_states_act_two_rdm_alpha_alpha_mo(l,k,j,i,istate) + rdmbb = all_states_act_two_rdm_beta_beta_mo(l,k,j,i,istate) + rdmab = all_states_act_two_rdm_alpha_beta_mo(l,k,j,i,istate) + rdmtot = all_states_act_two_rdm_spin_trace_mo(l,k,j,i,istate) + + accu_ab(istate) += vijkl * rdmab + accu_aa(istate) += vijkl * rdmaa + accu_bb(istate) += vijkl * rdmbb + accu_tot(istate)+= vijkl * rdmtot + enddo + enddo + enddo + enddo + print*,'' + print*,'Active space only energy ' + print*,'accu_aa(istate) = ',accu_aa(istate) + print*,'accu_bb(istate) = ',accu_bb(istate) + print*,'accu_ab(istate) = ',accu_ab(istate) + print*,'' + print*,'sum (istate) = ',accu_aa(istate) + accu_bb(istate) + accu_ab(istate) + print*,'accu_tot(istate) = ',accu_tot(istate) + print*,'psi_energy_two_e(istate) = ',psi_energy_two_e(istate) + enddo + +end + +subroutine routine_full_mos + implicit none + integer :: i,j,k,l,iorb,jorb,korb,lorb,istate + BEGIN_DOC +! This routine computes the two electron repulsion using various providers +! + END_DOC + + double precision :: vijkl,rdmaa,get_two_e_integral,rdmab,rdmbb,rdmtot + double precision :: accu_aa(N_states),accu_bb(N_states),accu_ab(N_states),accu_tot(N_states) + accu_aa = 0.d0 + accu_ab = 0.d0 + accu_bb = 0.d0 + accu_tot = 0.d0 + do istate = 1, N_states + do i = 1, n_core_inact_act_orb + iorb = list_core_inact_act(i) + do j = 1, n_core_inact_act_orb + jorb = list_core_inact_act(j) + do k = 1, n_core_inact_act_orb + korb = list_core_inact_act(k) + do l = 1, n_core_inact_act_orb + lorb = list_core_inact_act(l) + vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + + rdmaa = all_states_full_two_rdm_alpha_alpha_mo(l,k,j,i,istate) + rdmab = all_states_full_two_rdm_alpha_beta_mo(l,k,j,i,istate) + rdmbb = all_states_full_two_rdm_beta_beta_mo(l,k,j,i,istate) + rdmtot = all_states_full_two_rdm_spin_trace_mo(l,k,j,i,istate) + + accu_ab(istate) += vijkl * rdmab + accu_aa(istate) += vijkl * rdmaa + accu_bb(istate) += vijkl * rdmbb + accu_tot(istate)+= vijkl * rdmtot + enddo + enddo + enddo + enddo + print*,'Full energy ' + print*,'accu_aa(istate) = ',accu_aa(istate) + print*,'accu_bb(istate) = ',accu_bb(istate) + print*,'accu_ab(istate) = ',accu_ab(istate) + print*,'' + print*,'sum (istate) = ',accu_aa(istate) + accu_bb(istate) + accu_ab(istate) + print*,'accu_tot(istate) = ',accu_tot(istate) + print*,'psi_energy_two_e(istate) = ',psi_energy_two_e(istate) + enddo + +end