LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 19:13:29 +01:00
Code Releases Activity

all 2rdm clean and work with openmp

Browse Source
This commit is contained in:
Emmanuel Giner 2020-03-22 17:15:39 +01:00
parent d96108f772
commit 6b282c042c
8 changed files with 1033 additions and 351 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

130
src/two_body_rdm/act_2_rdm.irp.f Normal file
View File

@ -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
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
!
! !!!!! 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
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
!
! !!!!! 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
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
!
! !!!!! 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
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
!
! !!!!! 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

106
src/two_body_rdm/all_states_prov.irp.f
View File

@ -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 * <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \alpha} |Psi>
!
! !!!!! 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
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
!
! !!!!! 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
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
!
! !!!!! 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'} <Psi| a^{\dagger}_{i \sigma} a^{\dagger}_{j \sigma'} a_{l \sigma'} a_{k \sigma} |Psi>
!
! !!!!! 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

261
src/two_body_rdm/example.irp.f Normal file
View File

@ -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

176
src/two_body_rdm/all_states_full_orb_prov.irp.f → src/two_body_rdm/full_orb_2_rdm.irp.f
View File

@ -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
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
!
@ -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
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \alpha} |Psi>
!
@ -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
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
!
@ -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
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
!
@ -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

32
src/two_body_rdm/state_av_2rdm.irp.f → src/two_body_rdm/state_av_act_2rdm.irp.f
View File

@ -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
! = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
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
! = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
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
! = <Psi| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi>
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

528
src/two_body_rdm/state_av_full_orb_2_rdm.irp.f Normal file
View File

@ -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
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
!
! !!!!! 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
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \alpha} |Psi>
!
! !!!!! 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
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
!
! !!!!! 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
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
!
! !!!!! 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

133
src/two_body_rdm/test_2_rdm.irp.f
View File

@ -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

18
src/two_rdm_routines/all_states_david_openmp.irp.f
View File

@ -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