From 6b282c042ce14126ab822ad7209626ba7b54e85c Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sun, 22 Mar 2020 17:15:39 +0100 Subject: [PATCH] all 2rdm clean and work with openmp --- src/two_body_rdm/act_2_rdm.irp.f | 130 +++++ src/two_body_rdm/all_states_prov.irp.f | 106 ---- src/two_body_rdm/example.irp.f | 261 +++++++++ ...ll_orb_prov.irp.f => full_orb_2_rdm.irp.f} | 176 +++--- ..._av_2rdm.irp.f => state_av_act_2rdm.irp.f} | 32 +- .../state_av_full_orb_2_rdm.irp.f | 528 ++++++++++++++++++ src/two_body_rdm/test_2_rdm.irp.f | 133 +---- .../all_states_david_openmp.irp.f | 18 +- 8 files changed, 1033 insertions(+), 351 deletions(-) create mode 100644 src/two_body_rdm/act_2_rdm.irp.f delete mode 100644 src/two_body_rdm/all_states_prov.irp.f create mode 100644 src/two_body_rdm/example.irp.f rename src/two_body_rdm/{all_states_full_orb_prov.irp.f => full_orb_2_rdm.irp.f} (57%) rename src/two_body_rdm/{state_av_2rdm.irp.f => state_av_act_2rdm.irp.f} (59%) create mode 100644 src/two_body_rdm/state_av_full_orb_2_rdm.irp.f diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f new file mode 100644 index 00000000..af22946f --- /dev/null +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -0,0 +1,130 @@ + + BEGIN_PROVIDER [double precision, act_two_rdm_ab_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] + implicit none + BEGIN_DOC +! act_two_rdm_ab_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 +! +! act_two_rdm_ab_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 + integer :: ispin + double precision :: wall_1, wall_2 + ! condition for alpha/beta spin + print*,'' + print*,'' + print*,'' + print*,'Providing act_two_rdm_ab_mo ' + ispin = 3 + print*,'ispin = ',ispin + act_two_rdm_ab_mo = 0.d0 + call wall_time(wall_1) + call orb_range_two_rdm_openmp(act_two_rdm_ab_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + call wall_time(wall_2) + print*,'Wall time to provide act_two_rdm_ab_mo',wall_2 - wall_1 + END_PROVIDER + + + BEGIN_PROVIDER [double precision, act_two_rdm_aa_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] + implicit none + BEGIN_DOC +! act_two_rdm_aa_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 +! +! act_two_rdm_aa_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 + integer :: ispin + double precision :: wall_1, wall_2 + ! condition for alpha/beta spin + print*,'' + print*,'' + print*,'' + print*,'Providing act_two_rdm_aa_mo ' + ispin = 1 + print*,'ispin = ',ispin + act_two_rdm_aa_mo = 0.d0 + call wall_time(wall_1) + call orb_range_two_rdm_openmp(act_two_rdm_aa_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + call wall_time(wall_2) + print*,'Wall time to provide act_two_rdm_aa_mo',wall_2 - wall_1 + END_PROVIDER + + + BEGIN_PROVIDER [double precision, act_two_rdm_bb_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] + implicit none + BEGIN_DOC +! act_two_rdm_bb_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" +! +! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is beta, electron 2 is beta +! +! act_two_rdm_bb_mo(i,j,k,l,istate) = i:beta, j:beta, j:beta, l:beta +! +! Therefore you don't necessayr have symmetry between electron 1 and 2 + END_DOC + integer :: ispin + double precision :: wall_1, wall_2 + ! condition for beta/beta spin + print*,'' + print*,'' + print*,'' + print*,'Providing act_two_rdm_bb_mo ' + ispin = 2 + print*,'ispin = ',ispin + act_two_rdm_bb_mo = 0.d0 + call wall_time(wall_1) + call orb_range_two_rdm_openmp(act_two_rdm_bb_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + call wall_time(wall_2) + print*,'Wall time to provide act_two_rdm_bb_mo',wall_2 - wall_1 + END_PROVIDER + + BEGIN_PROVIDER [double precision, act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] + implicit none + BEGIN_DOC +! act_two_rdm_spin_trace_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" +! +! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is beta, electron 2 is beta +! +! act_two_rdm_spin_trace_mo(i,j,k,l,istate) = i:beta, j:beta, j:beta, l:beta +! +! Therefore you don't necessayr have symmetry between electron 1 and 2 + END_DOC + integer :: ispin + double precision :: wall_1, wall_2 + ! condition for beta/beta spin + print*,'' + print*,'' + print*,'' + print*,'Providing act_two_rdm_spin_trace_mo ' + ispin = 4 + print*,'ispin = ',ispin + act_two_rdm_spin_trace_mo = 0.d0 + call wall_time(wall_1) + call orb_range_two_rdm_openmp(act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + call wall_time(wall_2) + print*,'Wall time to provide act_two_rdm_spin_trace_mo',wall_2 - wall_1 + END_PROVIDER diff --git a/src/two_body_rdm/all_states_prov.irp.f b/src/two_body_rdm/all_states_prov.irp.f deleted file mode 100644 index 5523faa2..00000000 --- a/src/two_body_rdm/all_states_prov.irp.f +++ /dev/null @@ -1,106 +0,0 @@ - - - - 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 - 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 - integer :: ispin - ! condition for alpha/beta spin - ispin = 1 - all_states_act_two_rdm_alpha_alpha_mo = 0.D0 - double precision :: wall_1,wall_2 - call wall_time(wall_1) - print*,'providing all_states_act_two_rdm_alpha_alpha_mo ...' - 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)) - call wall_time(wall_2) - print*,'time to provide all_states_act_two_rdm_alpha_alpha_mo',wall_2 - wall_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)] - implicit none - BEGIN_DOC -! all_states_act_two_rdm_beta_beta_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta electrons -! -! -! -! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" - END_DOC - integer :: ispin - ! condition for alpha/beta spin - ispin = 2 - all_states_act_two_rdm_beta_beta_mo = 0.d0 - double precision :: wall_1,wall_2 - call wall_time(wall_1) - print*,'providing all_states_act_two_rdm_beta_beta_mo ...' - call orb_range_all_states_two_rdm(all_states_act_two_rdm_beta_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) - call wall_time(wall_2) - print*,'time to provide all_states_act_two_rdm_beta_beta_mo',wall_2 - wall_1 - - END_PROVIDER - - 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 - 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 - integer :: ispin - double precision :: wall_1, wall_2 - ! condition for alpha/beta spin - call wall_time(wall_1) - print*,'providing all_states_act_two_rdm_alpha_beta_mo ...' - ispin = 3 - print*,'ispin = ',ispin - all_states_act_two_rdm_alpha_beta_mo = 0.d0 - call wall_time(wall_1) - call orb_range_all_states_two_rdm(all_states_act_two_rdm_alpha_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) - - call wall_time(wall_2) - print*,'time to provide all_states_act_two_rdm_alpha_beta_mo',wall_2 - wall_1 - END_PROVIDER - - - BEGIN_PROVIDER [double precision, all_states_act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] - implicit none - BEGIN_DOC -! all_states_act_two_rdm_spin_trace_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM -! -! \sum_{\sigma, \sigma'} -! -! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" -! -! The active part of the two-electron energy for the state istate can be computed as: -! -! \sum_{i,j,k,l = 1, n_act_orb} all_states_act_two_rdm_spin_trace_mo(i,j,k,l,istate) * < ii jj | kk ll > -! -! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l) - END_DOC - integer :: ispin,i,j,k,l,istate - ! condition for alpha/beta spin - ispin = 4 - all_states_act_two_rdm_spin_trace_mo = 0.d0 - - double precision :: wall_1,wall_2 - call wall_time(wall_1) - print*,'providing all_states_act_two_rdm_spin_trace_mo ...' - 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)) - call wall_time(wall_2) - print*,'time to provide all_states_act_two_rdm_spin_trace_mo',wall_2 - wall_1 - END_PROVIDER - diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f new file mode 100644 index 00000000..92e24af2 --- /dev/null +++ b/src/two_body_rdm/example.irp.f @@ -0,0 +1,261 @@ + +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,get_two_e_integral + double precision :: wee_ab(N_states),rdmab + double precision :: wee_bb(N_states),rdmbb + double precision :: wee_aa(N_states),rdmaa + double precision :: wee_tot(N_states),rdmtot + double precision :: wee_aa_st_av, rdm_aa_st_av + double precision :: wee_bb_st_av, rdm_bb_st_av + double precision :: wee_ab_st_av, rdm_ab_st_av + double precision :: wee_tot_st_av, rdm_tot_st_av + double precision :: wee_aa_st_av_2,wee_ab_st_av_2,wee_bb_st_av_2,wee_tot_st_av_2,wee_tot_st_av_3 + + wee_ab = 0.d0 + wee_bb = 0.d0 + wee_aa = 0.d0 + wee_tot = 0.d0 + + wee_aa_st_av_2 = 0.d0 + wee_bb_st_av_2 = 0.d0 + wee_ab_st_av_2 = 0.d0 + wee_tot_st_av_2 = 0.d0 + wee_tot_st_av_3 = 0.d0 + + + iorb = 1 + jorb = 1 + korb = 1 + lorb = 1 + vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + provide act_two_rdm_ab_mo act_two_rdm_aa_mo act_two_rdm_bb_mo act_two_rdm_spin_trace_mo + provide state_av_act_two_rdm_ab_mo state_av_act_two_rdm_aa_mo + provide state_av_act_two_rdm_bb_mo state_av_act_two_rdm_spin_trace_mo + print*,'**************************' + print*,'**************************' + 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) + + + rdmab = act_two_rdm_ab_mo(l,k,j,i,istate) + rdmaa = act_two_rdm_aa_mo(l,k,j,i,istate) + rdmbb = act_two_rdm_bb_mo(l,k,j,i,istate) + rdmtot = act_two_rdm_spin_trace_mo(l,k,j,i,istate) + + + wee_ab(istate) += vijkl * rdmab + wee_aa(istate) += vijkl * rdmaa + wee_bb(istate) += vijkl * rdmbb + wee_tot(istate) += vijkl * rdmtot + + enddo + enddo + enddo + enddo + wee_aa_st_av_2 += wee_aa(istate) * state_average_weight(istate) + wee_bb_st_av_2 += wee_aa(istate) * state_average_weight(istate) + wee_ab_st_av_2 += wee_aa(istate) * state_average_weight(istate) + wee_tot_st_av_2 += wee_tot(istate) * state_average_weight(istate) + wee_tot_st_av_3 += psi_energy_two_e(istate) * state_average_weight(istate) + print*,'' + print*,'' + print*,'Active space only energy for state ',istate + print*,'wee_aa(istate) = ',wee_aa(istate) + print*,'wee_bb(istate) = ',wee_bb(istate) + print*,'wee_ab(istate) = ',wee_ab(istate) + print*,'' + print*,'sum (istate) = ',wee_aa(istate) + wee_bb(istate) + wee_ab(istate) + print*,'wee_tot = ',wee_tot(istate) + print*,'Full energy ' + print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate) + enddo + + wee_aa_st_av = 0.d0 + wee_bb_st_av = 0.d0 + wee_ab_st_av = 0.d0 + wee_tot_st_av = 0.d0 + 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) + + rdm_aa_st_av = state_av_act_two_rdm_aa_mo(l,k,j,i) + rdm_bb_st_av = state_av_act_two_rdm_bb_mo(l,k,j,i) + rdm_ab_st_av = state_av_act_two_rdm_ab_mo(l,k,j,i) + rdm_tot_st_av = state_av_act_two_rdm_spin_trace_mo(l,k,j,i) + + wee_aa_st_av += vijkl * rdm_aa_st_av + wee_bb_st_av += vijkl * rdm_bb_st_av + wee_ab_st_av += vijkl * rdm_ab_st_av + wee_tot_st_av += vijkl * rdm_tot_st_av + enddo + enddo + enddo + enddo + print*,'' + print*,'' + print*,'' + print*,'STATE AVERAGE ENERGY ' + print*,'Active space only energy for state ',istate + print*,'wee_aa_st_av = ',wee_aa_st_av + print*,'wee_aa_st_av_2 = ',wee_aa_st_av_2 + print*,'wee_bb_st_av = ',wee_bb_st_av + print*,'wee_bb_st_av_2 = ',wee_bb_st_av_2 + print*,'wee_ab_st_av = ',wee_ab_st_av + print*,'wee_ab_st_av_2 = ',wee_ab_st_av_2 + print*,'Sum of components = ',wee_aa_st_av+wee_bb_st_av+wee_ab_st_av + print*,'Sum of components_2 = ',wee_aa_st_av_2+wee_bb_st_av_2+wee_ab_st_av_2 + print*,'' + print*,'Full energy ' + print*,'wee_tot_st_av = ',wee_tot_st_av + print*,'wee_tot_st_av_2 = ',wee_tot_st_av_2 + print*,'wee_tot_st_av_3 = ',wee_tot_st_av_3 + +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 :: wee_aa(N_states),wee_bb(N_states),wee_ab(N_states),wee_tot(N_states) + double precision :: wee_aa_st_av, rdm_aa_st_av + double precision :: wee_bb_st_av, rdm_bb_st_av + double precision :: wee_ab_st_av, rdm_ab_st_av + double precision :: wee_tot_st_av, rdm_tot_st_av + double precision :: wee_aa_st_av_2,wee_ab_st_av_2,wee_bb_st_av_2,wee_tot_st_av_2,wee_tot_st_av_3 + wee_aa = 0.d0 + wee_ab = 0.d0 + wee_bb = 0.d0 + wee_tot = 0.d0 + + wee_aa_st_av_2 = 0.d0 + wee_bb_st_av_2 = 0.d0 + wee_ab_st_av_2 = 0.d0 + wee_tot_st_av_2 = 0.d0 + wee_tot_st_av_3 = 0.d0 + + + iorb = 1 + jorb = 1 + korb = 1 + lorb = 1 + vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + provide full_occ_two_rdm_ab_mo full_occ_two_rdm_aa_mo full_occ_two_rdm_bb_mo full_occ_two_rdm_spin_trace_mo + print*,'**************************' + print*,'**************************' + 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 = full_occ_two_rdm_aa_mo(l,k,j,i,istate) + rdmab = full_occ_two_rdm_ab_mo(l,k,j,i,istate) + rdmbb = full_occ_two_rdm_bb_mo(l,k,j,i,istate) + rdmtot = full_occ_two_rdm_spin_trace_mo(l,k,j,i,istate) + + wee_ab(istate) += vijkl * rdmab + wee_aa(istate) += vijkl * rdmaa + wee_bb(istate) += vijkl * rdmbb + wee_tot(istate)+= vijkl * rdmtot + enddo + enddo + enddo + enddo + wee_aa_st_av_2 += wee_aa(istate) * state_average_weight(istate) + wee_bb_st_av_2 += wee_bb(istate) * state_average_weight(istate) + wee_ab_st_av_2 += wee_ab(istate) * state_average_weight(istate) + wee_tot_st_av_2 += wee_tot(istate) * state_average_weight(istate) + wee_tot_st_av_3 += psi_energy_two_e(istate) * state_average_weight(istate) + print*,'' + print*,'' + print*,'Full energy for state ',istate + print*,'wee_aa(istate) = ',wee_aa(istate) + print*,'wee_bb(istate) = ',wee_bb(istate) + print*,'wee_ab(istate) = ',wee_ab(istate) + print*,'' + print*,'sum (istate) = ',wee_aa(istate) + wee_bb(istate) + wee_ab(istate) + print*,'wee_tot(istate) = ',wee_tot(istate) + print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate) + enddo + + wee_aa_st_av = 0.d0 + wee_bb_st_av = 0.d0 + wee_ab_st_av = 0.d0 + wee_tot_st_av = 0.d0 + 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) + + rdm_aa_st_av = state_av_full_occ_two_rdm_aa_mo(l,k,j,i) + rdm_bb_st_av = state_av_full_occ_two_rdm_bb_mo(l,k,j,i) + rdm_ab_st_av = state_av_full_occ_two_rdm_ab_mo(l,k,j,i) + rdm_tot_st_av = state_av_full_occ_two_rdm_spin_trace_mo(l,k,j,i) + + wee_aa_st_av += vijkl * rdm_aa_st_av + wee_bb_st_av += vijkl * rdm_bb_st_av + wee_ab_st_av += vijkl * rdm_ab_st_av + wee_tot_st_av += vijkl * rdm_tot_st_av + enddo + enddo + enddo + enddo + print*,'' + print*,'' + print*,'' + print*,'STATE AVERAGE ENERGY ' + print*,'wee_aa_st_av = ',wee_aa_st_av + print*,'wee_aa_st_av_2 = ',wee_aa_st_av_2 + print*,'wee_bb_st_av = ',wee_bb_st_av + print*,'wee_bb_st_av_2 = ',wee_bb_st_av_2 + print*,'wee_ab_st_av = ',wee_ab_st_av + print*,'wee_ab_st_av_2 = ',wee_ab_st_av_2 + print*,'Sum of components = ',wee_aa_st_av + wee_bb_st_av + wee_ab_st_av + print*,'Sum of components_2 = ',wee_aa_st_av_2 + wee_bb_st_av_2 + wee_ab_st_av_2 + print*,'' + print*,'Full energy ' + print*,'wee_tot_st_av = ',wee_tot_st_av + print*,'wee_tot_st_av_2 = ',wee_tot_st_av_2 + print*,'wee_tot_st_av_3 = ',wee_tot_st_av_3 + +end diff --git a/src/two_body_rdm/all_states_full_orb_prov.irp.f b/src/two_body_rdm/full_orb_2_rdm.irp.f similarity index 57% rename from src/two_body_rdm/all_states_full_orb_prov.irp.f rename to src/two_body_rdm/full_orb_2_rdm.irp.f index 55fa78ca..1843cd3c 100644 --- a/src/two_body_rdm/all_states_full_orb_prov.irp.f +++ b/src/two_body_rdm/full_orb_2_rdm.irp.f @@ -1,10 +1,10 @@ - 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)] + BEGIN_PROVIDER [double precision, full_occ_two_rdm_ab_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 + full_occ_two_rdm_ab_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 +! full_occ_two_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons ! ! ! @@ -14,13 +14,13 @@ ! ! !!!!! 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 +! act_two_rdm_ab_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 + full_occ_two_rdm_ab_mo = 0.d0 do istate = 1, N_states !! PURE ACTIVE PART ALPHA-BETA !! @@ -33,8 +33,8 @@ 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) + full_occ_two_rdm_ab_mo(lorb,korb,jorb,iorb,istate) = & + act_two_rdm_ab_mo(l,k,j,i,istate) enddo enddo enddo @@ -48,7 +48,7 @@ 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) + full_occ_two_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -62,7 +62,7 @@ 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) + full_occ_two_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -74,7 +74,7 @@ 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 + full_occ_two_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0 enddo enddo @@ -91,7 +91,7 @@ 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) + full_occ_two_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -105,7 +105,7 @@ 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) + full_occ_two_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -117,7 +117,7 @@ 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 + full_occ_two_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0 enddo enddo endif @@ -126,12 +126,12 @@ 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)] + BEGIN_PROVIDER [double precision, full_occ_two_rdm_aa_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 + full_occ_two_rdm_aa_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 +! full_occ_two_rdm_aa_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/alpha electrons ! ! ! @@ -153,8 +153,8 @@ 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) + full_occ_two_rdm_aa_mo(lorb,korb,jorb,iorb,istate) = & + act_two_rdm_aa_mo(l,k,j,i,istate) enddo enddo enddo @@ -168,11 +168,11 @@ 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) + full_occ_two_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_two_rdm_aa_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) + full_occ_two_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_two_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -182,8 +182,8 @@ 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 + full_occ_two_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + full_occ_two_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 enddo enddo @@ -199,11 +199,11 @@ 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) + full_occ_two_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_two_rdm_aa_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) + full_occ_two_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_two_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -213,8 +213,8 @@ 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 + full_occ_two_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + full_occ_two_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 enddo enddo endif @@ -222,12 +222,12 @@ 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)] + BEGIN_PROVIDER [double precision, full_occ_two_rdm_bb_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 + full_occ_two_rdm_bb_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 +! full_occ_two_rdm_bb_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta/beta electrons ! ! ! @@ -249,8 +249,8 @@ 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) + full_occ_two_rdm_bb_mo(lorb,korb,jorb,iorb,istate) = & + act_two_rdm_bb_mo(l,k,j,i,istate) enddo enddo enddo @@ -264,11 +264,11 @@ 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) + full_occ_two_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_two_rdm_bb_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) + full_occ_two_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_two_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -278,8 +278,8 @@ 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 + full_occ_two_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + full_occ_two_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 enddo enddo @@ -295,11 +295,11 @@ 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) + full_occ_two_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_two_rdm_bb_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) + full_occ_two_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_two_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -309,8 +309,8 @@ 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 + full_occ_two_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + full_occ_two_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 enddo enddo endif @@ -318,12 +318,12 @@ 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)] + BEGIN_PROVIDER [double precision, full_occ_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 + full_occ_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 +! full_occ_two_rdm_bb_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta/beta electrons ! ! ! @@ -346,8 +346,8 @@ 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) + full_occ_two_rdm_spin_trace_mo(lorb,korb,jorb,iorb,istate) += & + act_two_rdm_spin_trace_mo(l,k,j,i,istate) enddo enddo enddo @@ -364,11 +364,11 @@ 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) + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_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) + full_occ_two_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_two_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -377,8 +377,8 @@ 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 + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + full_occ_two_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 enddo enddo if (.not.no_core_density)then @@ -390,11 +390,11 @@ 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) + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_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) + full_occ_two_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) + full_occ_two_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate) enddo enddo enddo @@ -403,8 +403,8 @@ 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 + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + full_occ_two_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 enddo enddo endif @@ -420,11 +420,11 @@ 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) + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_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) + full_occ_two_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_two_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -433,8 +433,8 @@ 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 + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + full_occ_two_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 enddo enddo if (.not.no_core_density)then @@ -446,11 +446,11 @@ 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) + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_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) + full_occ_two_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) + full_occ_two_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -459,8 +459,8 @@ 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 + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0 + full_occ_two_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0 enddo enddo endif @@ -476,14 +476,14 @@ 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) + full_occ_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) + full_occ_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) + full_occ_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) + full_occ_two_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -493,8 +493,8 @@ 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 + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0 + full_occ_two_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0 enddo enddo @@ -510,14 +510,14 @@ 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) + full_occ_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) + full_occ_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) + full_occ_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) + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate) enddo enddo enddo @@ -527,8 +527,8 @@ 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 + full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0 + full_occ_two_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0 enddo enddo diff --git a/src/two_body_rdm/state_av_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f similarity index 59% rename from src/two_body_rdm/state_av_2rdm.irp.f rename to src/two_body_rdm/state_av_act_2rdm.irp.f index c8c625bc..640137b5 100644 --- a/src/two_body_rdm/state_av_2rdm.irp.f +++ b/src/two_body_rdm/state_av_act_2rdm.irp.f @@ -1,9 +1,9 @@ - BEGIN_PROVIDER [double precision, state_av_act_two_rdm_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_aa_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_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs +! state_av_act_two_rdm_aa_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)) @@ -11,20 +11,20 @@ integer :: ispin ! condition for alpha/beta spin ispin = 1 - state_av_act_two_rdm_alpha_alpha_mo = 0.D0 + state_av_act_two_rdm_aa_mo = 0.D0 call wall_time(wall_1) double precision :: wall_1, wall_2 - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_aa_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_2) - print*,'Wall time to provide state_av_act_two_rdm_alpha_alpha_mo',wall_2 - wall_1 + print*,'Wall time to provide state_av_act_two_rdm_aa_mo',wall_2 - wall_1 END_PROVIDER - BEGIN_PROVIDER [double precision, state_av_act_two_rdm_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_bb_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_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs +! state_av_act_two_rdm_bb_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)) @@ -32,20 +32,20 @@ integer :: ispin ! condition for alpha/beta spin ispin = 2 - state_av_act_two_rdm_beta_beta_mo = 0.d0 + state_av_act_two_rdm_bb_mo = 0.d0 call wall_time(wall_1) double precision :: wall_1, wall_2 - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_beta_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_bb_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_2) - print*,'Wall time to provide state_av_act_two_rdm_beta_beta_mo',wall_2 - wall_1 + print*,'Wall time to provide state_av_act_two_rdm_bb_mo',wall_2 - wall_1 END_PROVIDER - BEGIN_PROVIDER [double precision, state_av_act_two_rdm_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_ab_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_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs +! state_av_act_two_rdm_ab_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)) @@ -55,15 +55,15 @@ print*,'' print*,'' print*,'' - print*,'providint state_av_act_two_rdm_alpha_beta_mo ' + print*,'providint state_av_act_two_rdm_ab_mo ' ispin = 3 print*,'ispin = ',ispin - state_av_act_two_rdm_alpha_beta_mo = 0.d0 + state_av_act_two_rdm_ab_mo = 0.d0 call wall_time(wall_1) double precision :: wall_1, wall_2 - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_alpha_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_ab_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_2) - print*,'Wall time to provide state_av_act_two_rdm_alpha_beta_mo',wall_2 - wall_1 + print*,'Wall time to provide state_av_act_two_rdm_ab_mo',wall_2 - wall_1 END_PROVIDER diff --git a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f new file mode 100644 index 00000000..f5b3e18c --- /dev/null +++ b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f @@ -0,0 +1,528 @@ + + BEGIN_PROVIDER [double precision, state_av_full_occ_two_rdm_ab_mo, (n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb)] + implicit none + state_av_full_occ_two_rdm_ab_mo = 0.d0 + integer :: i,j,k,l,iorb,jorb,korb,lorb + BEGIN_DOC +! state_av_full_occ_two_rdm_ab_mo(i,j,k,l) = STATE AVERAGE 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 +! +! state_av_full_occ_two_rdm_ab_mo(i,j,k,l) = 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 + state_av_full_occ_two_rdm_ab_mo = 0.d0 + 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 + state_av_full_occ_two_rdm_ab_mo(lorb,korb,jorb,iorb) = & + state_av_act_two_rdm_ab_mo(l,k,j,i) + 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 + state_av_full_occ_two_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb) + 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 + state_av_full_occ_two_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb) + 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 + state_av_full_occ_two_rdm_ab_mo(korb,jorb,korb,jorb) = 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 + state_av_full_occ_two_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb) + 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 + state_av_full_occ_two_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb) + 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 + state_av_full_occ_two_rdm_ab_mo(korb,jorb,korb,jorb) = 1.D0 + enddo + enddo + endif + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, state_av_full_occ_two_rdm_aa_mo, (n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb)] + implicit none + state_av_full_occ_two_rdm_aa_mo = 0.d0 + integer :: i,j,k,l,iorb,jorb,korb,lorb + BEGIN_DOC +! state_av_full_occ_two_rdm_aa_mo(i,j,k,l) = STATE AVERAGE 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 + + !! 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) + state_av_full_occ_two_rdm_aa_mo(lorb,korb,jorb,iorb) = & + state_av_act_two_rdm_aa_mo(l,k,j,i) + 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 + state_av_full_occ_two_rdm_aa_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_two_rdm_aa_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + ! 1 2 1 2 : EXCHANGE TERM + state_av_full_occ_two_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_two_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + 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) + state_av_full_occ_two_rdm_aa_mo(korb,jorb,korb,jorb) += 0.5d0 + state_av_full_occ_two_rdm_aa_mo(korb,jorb,jorb,korb) -= 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 + state_av_full_occ_two_rdm_aa_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_two_rdm_aa_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + ! 1 2 1 2 : EXCHANGE TERM + state_av_full_occ_two_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_two_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + 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) + state_av_full_occ_two_rdm_aa_mo(korb,jorb,korb,jorb) += 0.5d0 + state_av_full_occ_two_rdm_aa_mo(korb,jorb,jorb,korb) -= 0.5d0 + enddo + enddo + endif + + END_PROVIDER + + BEGIN_PROVIDER [double precision, state_av_full_occ_two_rdm_bb_mo, (n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb)] + implicit none + state_av_full_occ_two_rdm_bb_mo = 0.d0 + integer :: i,j,k,l,iorb,jorb,korb,lorb + BEGIN_DOC +! state_av_full_occ_two_rdm_bb_mo(i,j,k,l) = STATE AVERAGE 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 + + !! 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) + state_av_full_occ_two_rdm_bb_mo(lorb,korb,jorb,iorb) = & + state_av_act_two_rdm_bb_mo(l,k,j,i) + 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 + state_av_full_occ_two_rdm_bb_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_two_rdm_bb_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + ! 1 2 1 2 : EXCHANGE TERM + state_av_full_occ_two_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_two_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + 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) + state_av_full_occ_two_rdm_bb_mo(korb,jorb,korb,jorb) += 0.5d0 + state_av_full_occ_two_rdm_bb_mo(korb,jorb,jorb,korb) -= 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 + state_av_full_occ_two_rdm_bb_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_two_rdm_bb_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + ! 1 2 1 2 : EXCHANGE TERM + state_av_full_occ_two_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_two_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + 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) + state_av_full_occ_two_rdm_bb_mo(korb,jorb,korb,jorb) += 0.5d0 + state_av_full_occ_two_rdm_bb_mo(korb,jorb,jorb,korb) -= 0.5d0 + enddo + enddo + endif + + END_PROVIDER + + BEGIN_PROVIDER [double precision, state_av_full_occ_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)] + implicit none + state_av_full_occ_two_rdm_spin_trace_mo = 0.d0 + integer :: i,j,k,l,iorb,jorb,korb,lorb + BEGIN_DOC +! state_av_full_occ_two_rdm_bb_mo(i,j,k,l) = STATE AVERAGE 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 + + !!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!! + !! 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) + state_av_full_occ_two_rdm_spin_trace_mo(lorb,korb,jorb,iorb) += & + state_av_act_two_rdm_spin_trace_mo(l,k,j,i) + 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 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + ! 1 2 1 2 : EXCHANGE TERM + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + 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) + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 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 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + ! 1 2 1 2 : EXCHANGE TERM + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + 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) + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 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 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + ! 1 2 1 2 : EXCHANGE TERM + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + 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) + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 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 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + ! 1 2 1 2 : EXCHANGE TERM + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + 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) + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 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 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + ! beta alph beta alph + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb) + ! BETA INACTIVE - ALPHA ACTIVE + ! beta alph beta alpha + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + ! alph beta alph beta + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb) + 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 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0 + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 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 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb) + ! beta alph beta alph + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb) + !! ALPHA ACTIVE - BETA CORE + ! alph beta alph beta + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb) + ! beta alph beta alph + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb) + 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 + state_av_full_occ_two_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0 + state_av_full_occ_two_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 0.5D0 + enddo + enddo + + endif + + END_PROVIDER diff --git a/src/two_body_rdm/test_2_rdm.irp.f b/src/two_body_rdm/test_2_rdm.irp.f index 9ba57e37..123261d8 100644 --- a/src/two_body_rdm/test_2_rdm.irp.f +++ b/src/two_body_rdm/test_2_rdm.irp.f @@ -2,138 +2,7 @@ program test_2_rdm implicit none read_wf = .True. touch read_wf -! call routine_full_mos call routine_active_only + call routine_full_mos 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) - double precision :: accu_ab_omp(N_states),rdmab_omp - double precision :: accu_bb_omp(N_states),rdmbb_omp - double precision :: accu_aa_omp(N_states),rdmaa_omp - double precision :: accu_tot_omp(N_states),rdmtot_omp - accu_ab_omp = 0.d0 - accu_bb_omp = 0.d0 - accu_aa_omp = 0.d0 - accu_tot_omp = 0.d0 - 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) - -! rdmab_omp = state_av_act_two_rdm_alpha_beta_mo(l,k,j,i) -! rdmbb_omp = state_av_act_two_rdm_beta_beta_mo(l,k,j,i) -! rdmaa_omp = state_av_act_two_rdm_alpha_alpha_mo(l,k,j,i) -! rdmtot_omp = state_av_act_two_rdm_spin_trace_mo(l,k,j,i) - - -! rdmab_omp = all_states_openmp_act_two_rdm_alpha_beta_mo(l,k,j,i,istate) -! rdmaa_omp = all_states_openmp_act_two_rdm_alpha_alpha_mo(l,k,j,i,istate) -! rdmbb_omp = all_states_openmp_act_two_rdm_beta_beta_mo(l,k,j,i,istate) - rdmtot_omp = all_states_openmp_act_two_rdm_spin_trace_mo(l,k,j,i,istate) - -! rdmab = all_states_act_two_rdm_alpha_beta_mo(l,k,j,i,istate) -! 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) - rdmtot = all_states_act_two_rdm_spin_trace_mo(l,k,j,i,istate) - -! accu_ab_omp(istate) += vijkl * rdmab_omp -! accu_aa_omp(istate) += vijkl * rdmaa_omp -! accu_bb_omp(istate) += vijkl * rdmbb_omp - accu_tot_omp(istate) += vijkl * rdmtot_omp - -! 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_aa_omp = ',accu_aa_omp(istate) -! print*,'accu_bb(istate) = ',accu_bb(istate) -! print*,'accu_bb_omp = ',accu_bb_omp(istate) -! print*,'accu_ab(istate) = ',accu_ab(istate) -! print*,'accu_ab_omp = ',accu_ab_omp(istate) - print*,'' -! print*,'sum (istate) = ',accu_aa(istate) + accu_bb(istate) + accu_ab(istate) - print*,'accu_tot(istate) = ',accu_tot(istate) - print*,'accu_tot_omp = ',accu_tot_omp(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 diff --git a/src/two_rdm_routines/all_states_david_openmp.irp.f b/src/two_rdm_routines/all_states_david_openmp.irp.f index 0da7e205..412c2f12 100644 --- a/src/two_rdm_routines/all_states_david_openmp.irp.f +++ b/src/two_rdm_routines/all_states_david_openmp.irp.f @@ -1,4 +1,4 @@ -subroutine orb_range_two_rdm_all_states_openmp_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sze) +subroutine orb_range_two_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sze) use bitmasks implicit none BEGIN_DOC @@ -30,7 +30,7 @@ subroutine orb_range_two_rdm_all_states_openmp_openmp(big_array,dim1,norb,list_o size(u_t, 1), & N_det, N_st) - call orb_range_two_rdm_all_states_openmp_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,1,N_det,0,1) + call orb_range_two_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) do k=1,N_st @@ -39,7 +39,7 @@ subroutine orb_range_two_rdm_all_states_openmp_openmp(big_array,dim1,norb,list_o end -subroutine orb_range_two_rdm_all_states_openmp_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine orb_range_two_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC @@ -58,15 +58,15 @@ subroutine orb_range_two_rdm_all_states_openmp_openmp_work(big_array,dim1,norb,l select case (N_int) case (1) - call orb_range_two_rdm_all_states_openmp_openmp_work_1(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_openmp_work_1(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (2) - call orb_range_two_rdm_all_states_openmp_openmp_work_2(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_openmp_work_2(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (3) - call orb_range_two_rdm_all_states_openmp_openmp_work_3(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_openmp_work_3(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (4) - call orb_range_two_rdm_all_states_openmp_openmp_work_4(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_openmp_work_4(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case default - call orb_range_two_rdm_all_states_openmp_openmp_work_N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_openmp_work_N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) end select end @@ -74,7 +74,7 @@ end BEGIN_TEMPLATE -subroutine orb_range_two_rdm_all_states_openmp_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine orb_range_two_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks use omp_lib implicit none