From 919662ee0b64abb6a20a3e61c1d1bb261fcb05c4 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 4 Jul 2019 16:16:57 +0200 Subject: [PATCH] beginning to rewrite two_rdm --- src/casscf/get_energy.irp.f | 19 + src/davidson/u0_wee_u0.irp.f | 2 +- src/two_body_rdm/ab_only_routines.irp.f | 18 +- src/two_body_rdm/all_2rdm_routines.irp.f | 18 +- src/two_body_rdm/all_states_2_rdm.irp.f | 8 +- src/two_body_rdm/all_states_routines.irp.f | 20 +- src/two_body_rdm/orb_range_2_rdm.irp.f | 8 +- src/two_body_rdm/orb_range_2_rdm_openmp.irp.f | 87 +++ src/two_body_rdm/orb_range_routines.irp.f | 20 +- .../orb_range_routines_openmp.irp.f | 544 ++++++++++++++ ...utines_compute_2rdm_orb_range_openmp.irp.f | 683 ++++++++++++++++++ src/two_body_rdm/two_rdm.irp.f | 2 +- 12 files changed, 1381 insertions(+), 48 deletions(-) create mode 100644 src/two_body_rdm/orb_range_2_rdm_openmp.irp.f create mode 100644 src/two_body_rdm/orb_range_routines_openmp.irp.f create mode 100644 src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f index 2a595fe7..384ff804 100644 --- a/src/casscf/get_energy.irp.f +++ b/src/casscf/get_energy.irp.f @@ -33,4 +33,23 @@ subroutine routine enddo enddo print*,'accu = ',accu(1) + + accu = 0.d0 + do ll = 1, n_act_orb + l = list_act(ll) + do kk = 1, n_act_orb + k = list_act(kk) + do jj = 1, n_act_orb + j = list_act(jj) + do ii = 1, n_act_orb + i = list_act(ii) + integral = get_two_e_integral(i,j,k,l,mo_integrals_map) + accu(1) += state_av_act_two_rdm_openmp_spin_trace_mo(ii,jj,kk,ll) * integral + enddo + enddo + enddo + enddo + print*,'accu = ',accu(1) + print*,'psi_energy_two_e = ',psi_energy_two_e + end diff --git a/src/davidson/u0_wee_u0.irp.f b/src/davidson/u0_wee_u0.irp.f index c1f163d4..0c543aca 100644 --- a/src/davidson/u0_wee_u0.irp.f +++ b/src/davidson/u0_wee_u0.irp.f @@ -6,7 +6,7 @@ BEGIN_PROVIDER [ double precision, psi_energy_two_e, (N_states) ] integer :: i,j call u_0_H_u_0_two_e(psi_energy_two_e,psi_coef,N_det,psi_det,N_int,N_states,psi_det_size) do i=N_det+1,N_states - psi_energy(i) = 0.d0 + psi_energy_two_e(i) = 0.d0 enddo END_PROVIDER diff --git a/src/two_body_rdm/ab_only_routines.irp.f b/src/two_body_rdm/ab_only_routines.irp.f index 9041c753..fb3c421c 100644 --- a/src/two_body_rdm/ab_only_routines.irp.f +++ b/src/two_body_rdm/ab_only_routines.irp.f @@ -1,5 +1,5 @@ - subroutine two_rdm_ab_nstates_openmp(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze) + subroutine two_rdm_ab_nstates(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze) use bitmasks implicit none BEGIN_DOC @@ -27,7 +27,7 @@ size(u_t, 1), & N_det, N_st) - call two_rdm_ab_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1) + call two_rdm_ab_nstates_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) do k=1,N_st @@ -37,7 +37,7 @@ end - subroutine two_rdm_ab_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + subroutine two_rdm_ab_nstates_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC @@ -55,20 +55,20 @@ select case (N_int) case (1) - call two_rdm_ab_nstates_openmp_work_1(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_work_1(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (2) - call two_rdm_ab_nstates_openmp_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (3) - call two_rdm_ab_nstates_openmp_work_3(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_work_3(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (4) - call two_rdm_ab_nstates_openmp_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case default - call two_rdm_ab_nstates_openmp_work_N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_work_N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) end select end BEGIN_TEMPLATE - subroutine two_rdm_ab_nstates_openmp_work_$N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + subroutine two_rdm_ab_nstates_work_$N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none integer, intent(in) :: N_st,sze,istart,iend,ishift,istep diff --git a/src/two_body_rdm/all_2rdm_routines.irp.f b/src/two_body_rdm/all_2rdm_routines.irp.f index 3f08b18f..5127e31f 100644 --- a/src/two_body_rdm/all_2rdm_routines.irp.f +++ b/src/two_body_rdm/all_2rdm_routines.irp.f @@ -1,4 +1,4 @@ -subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_0,N_st,sze) +subroutine all_two_rdm_dm_nstates(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_0,N_st,sze) use bitmasks implicit none BEGIN_DOC @@ -28,7 +28,7 @@ subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab, size(u_t, 1), & N_det, N_st) - call all_two_rdm_dm_nstates_openmp_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1) + call all_two_rdm_dm_nstates_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) do k=1,N_st @@ -38,7 +38,7 @@ subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab, end -subroutine all_two_rdm_dm_nstates_openmp_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine all_two_rdm_dm_nstates_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC @@ -58,21 +58,21 @@ subroutine all_two_rdm_dm_nstates_openmp_work(big_array_aa,big_array_bb,big_arra select case (N_int) case (1) - call all_two_rdm_dm_nstates_openmp_work_1(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call all_two_rdm_dm_nstates_work_1(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (2) - call all_two_rdm_dm_nstates_openmp_work_2(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call all_two_rdm_dm_nstates_work_2(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (3) - call all_two_rdm_dm_nstates_openmp_work_3(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call all_two_rdm_dm_nstates_work_3(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (4) - call all_two_rdm_dm_nstates_openmp_work_4(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call all_two_rdm_dm_nstates_work_4(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case default - call all_two_rdm_dm_nstates_openmp_work_N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call all_two_rdm_dm_nstates_work_N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) end select end BEGIN_TEMPLATE -subroutine all_two_rdm_dm_nstates_openmp_work_$N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine all_two_rdm_dm_nstates_work_$N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC diff --git a/src/two_body_rdm/all_states_2_rdm.irp.f b/src/two_body_rdm/all_states_2_rdm.irp.f index cd74758f..bc503223 100644 --- a/src/two_body_rdm/all_states_2_rdm.irp.f +++ b/src/two_body_rdm/all_states_2_rdm.irp.f @@ -14,7 +14,7 @@ ! condition for alpha/beta spin ispin = 1 all_states_act_two_rdm_alpha_alpha_mo = 0.D0 - call orb_range_all_states_two_rdm_openmp(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 orb_range_all_states_two_rdm(all_states_act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) END_PROVIDER @@ -31,7 +31,7 @@ ! condition for alpha/beta spin ispin = 2 all_states_act_two_rdm_beta_beta_mo = 0.d0 - call orb_range_all_states_two_rdm_openmp(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 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)) END_PROVIDER @@ -53,7 +53,7 @@ ispin = 3 print*,'ispin = ',ispin all_states_act_two_rdm_alpha_beta_mo = 0.d0 - call orb_range_all_states_two_rdm_openmp(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 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)) END_PROVIDER @@ -77,7 +77,7 @@ all_states_act_two_rdm_spin_trace_mo = 0.d0 integer :: i - call orb_range_all_states_two_rdm_openmp(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 orb_range_all_states_two_rdm(all_states_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) END_PROVIDER diff --git a/src/two_body_rdm/all_states_routines.irp.f b/src/two_body_rdm/all_states_routines.irp.f index 3084dd5b..af7cafc2 100644 --- a/src/two_body_rdm/all_states_routines.irp.f +++ b/src/two_body_rdm/all_states_routines.irp.f @@ -1,4 +1,4 @@ -subroutine orb_range_all_states_two_rdm_openmp(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_0,N_st,sze) +subroutine orb_range_all_states_two_rdm(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_0,N_st,sze) use bitmasks implicit none BEGIN_DOC @@ -31,7 +31,7 @@ subroutine orb_range_all_states_two_rdm_openmp(big_array,dim1,norb,list_orb,list size(u_t, 1), & N_det, N_st) - call orb_range_all_states_two_rdm_openmp_work(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,1,N_det,0,1) + call orb_range_all_states_two_rdm_work(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) do k=1,N_st @@ -40,7 +40,7 @@ subroutine orb_range_all_states_two_rdm_openmp(big_array,dim1,norb,list_orb,list end -subroutine orb_range_all_states_two_rdm_openmp_work(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine orb_range_all_states_two_rdm_work(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC @@ -60,15 +60,15 @@ subroutine orb_range_all_states_two_rdm_openmp_work(big_array,dim1,norb,list_orb select case (N_int) case (1) - call orb_range_all_states_two_rdm_openmp_work_1(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_all_states_two_rdm_work_1(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (2) - call orb_range_all_states_two_rdm_openmp_work_2(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_all_states_two_rdm_work_2(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (3) - call orb_range_all_states_two_rdm_openmp_work_3(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_all_states_two_rdm_work_3(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (4) - call orb_range_all_states_two_rdm_openmp_work_4(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_all_states_two_rdm_work_4(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case default - call orb_range_all_states_two_rdm_openmp_work_N_int(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_all_states_two_rdm_work_N_int(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) end select end @@ -76,7 +76,7 @@ end BEGIN_TEMPLATE -subroutine orb_range_all_states_two_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine orb_range_all_states_two_rdm_work_$N_int(big_array,dim1,norb,list_orb,list_orb_reverse,ispin,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC @@ -129,7 +129,7 @@ subroutine orb_range_all_states_two_rdm_openmp_work_$N_int(big_array,dim1,norb,l else if(ispin == 4)then spin_trace = .True. else - print*,'Wrong parameter for ispin in general_two_rdm_dm_nstates_openmp_work' + print*,'Wrong parameter for ispin in general_two_rdm_dm_nstates_work' print*,'ispin = ',ispin stop endif diff --git a/src/two_body_rdm/orb_range_2_rdm.irp.f b/src/two_body_rdm/orb_range_2_rdm.irp.f index 02df58d8..d441e1df 100644 --- a/src/two_body_rdm/orb_range_2_rdm.irp.f +++ b/src/two_body_rdm/orb_range_2_rdm.irp.f @@ -14,7 +14,7 @@ ! condition for alpha/beta spin ispin = 1 state_av_act_two_rdm_alpha_alpha_mo = 0.D0 - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_state_av(state_av_act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) END_PROVIDER @@ -31,7 +31,7 @@ ! condition for alpha/beta spin ispin = 2 state_av_act_two_rdm_beta_beta_mo = 0.d0 - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_beta_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_state_av(state_av_act_two_rdm_beta_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) END_PROVIDER @@ -53,7 +53,7 @@ ispin = 3 print*,'ispin = ',ispin state_av_act_two_rdm_alpha_beta_mo = 0.d0 - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_alpha_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_state_av(state_av_act_two_rdm_alpha_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) END_PROVIDER @@ -79,7 +79,7 @@ double precision :: wall_0,wall_1 call wall_time(wall_0) print*,'providing the state average TWO-RDM ...' - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_state_av(state_av_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_1) print*,'Time to provide the state average TWO-RDM',wall_1 - wall_0 diff --git a/src/two_body_rdm/orb_range_2_rdm_openmp.irp.f b/src/two_body_rdm/orb_range_2_rdm_openmp.irp.f new file mode 100644 index 00000000..70bf0201 --- /dev/null +++ b/src/two_body_rdm/orb_range_2_rdm_openmp.irp.f @@ -0,0 +1,87 @@ + + + + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + BEGIN_DOC +! state_av_act_two_rdm_openmp_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs +! = + END_DOC + allocate(state_weights(N_states)) + state_weights = 1.d0/dble(N_states) + integer :: ispin + ! condition for alpha/beta spin + ispin = 1 + state_av_act_two_rdm_openmp_alpha_alpha_mo = 0.D0 + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + BEGIN_DOC +! state_av_act_two_rdm_openmp_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs +! = + END_DOC + allocate(state_weights(N_states)) + state_weights = 1.d0/dble(N_states) + integer :: ispin + ! condition for alpha/beta spin + ispin = 2 + state_av_act_two_rdm_openmp_beta_beta_mo = 0.d0 + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_beta_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + BEGIN_DOC +! state_av_act_two_rdm_openmp_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs +! = + END_DOC + allocate(state_weights(N_states)) + state_weights = 1.d0/dble(N_states) + integer :: ispin + ! condition for alpha/beta spin + print*,'' + print*,'' + print*,'' + print*,'providint state_av_act_two_rdm_openmp_alpha_beta_mo ' + ispin = 3 + print*,'ispin = ',ispin + state_av_act_two_rdm_openmp_alpha_beta_mo = 0.d0 + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_alpha_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + BEGIN_DOC +! state_av_act_two_rdm_openmp_spin_trace_mo(i,j,k,l) = state average physicist spin trace two-body rdm restricted to the ACTIVE indices +! The active part of the two-electron energy can be computed as: +! +! \sum_{i,j,k,l = 1, n_act_orb} state_av_act_two_rdm_openmp_spin_trace_mo(i,j,k,l) * < ii jj | kk ll > +! +! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l) + END_DOC + double precision, allocatable :: state_weights(:) + allocate(state_weights(N_states)) + state_weights = 1.d0/dble(N_states) + integer :: ispin + ! condition for alpha/beta spin + ispin = 4 + state_av_act_two_rdm_openmp_spin_trace_mo = 0.d0 + integer :: i + double precision :: wall_0,wall_1 + call wall_time(wall_0) + print*,'providing the state average TWO-RDM ...' + call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + call wall_time(wall_1) + print*,'Time to provide the state average TWO-RDM',wall_1 - wall_0 + END_PROVIDER + diff --git a/src/two_body_rdm/orb_range_routines.irp.f b/src/two_body_rdm/orb_range_routines.irp.f index b82c4799..d5bd7d1c 100644 --- a/src/two_body_rdm/orb_range_routines.irp.f +++ b/src/two_body_rdm/orb_range_routines.irp.f @@ -1,4 +1,4 @@ -subroutine orb_range_two_rdm_state_av_openmp(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_0,N_st,sze) +subroutine orb_range_two_rdm_state_av(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_0,N_st,sze) use bitmasks implicit none BEGIN_DOC @@ -31,7 +31,7 @@ subroutine orb_range_two_rdm_state_av_openmp(big_array,dim1,norb,list_orb,list_o size(u_t, 1), & N_det, N_st) - call orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1) + call orb_range_two_rdm_state_av_work(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) do k=1,N_st @@ -40,7 +40,7 @@ subroutine orb_range_two_rdm_state_av_openmp(big_array,dim1,norb,list_orb,list_o end -subroutine orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine orb_range_two_rdm_state_av_work(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC @@ -60,15 +60,15 @@ subroutine orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,l select case (N_int) case (1) - call orb_range_two_rdm_state_av_openmp_work_1(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_state_av_work_1(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (2) - call orb_range_two_rdm_state_av_openmp_work_2(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_state_av_work_2(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (3) - call orb_range_two_rdm_state_av_openmp_work_3(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_state_av_work_3(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (4) - call orb_range_two_rdm_state_av_openmp_work_4(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_state_av_work_4(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case default - call orb_range_two_rdm_state_av_openmp_work_N_int(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_state_av_work_N_int(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) end select end @@ -76,7 +76,7 @@ end BEGIN_TEMPLATE -subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine orb_range_two_rdm_state_av_work_$N_int(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC @@ -130,7 +130,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis else if(ispin == 4)then spin_trace = .True. else - print*,'Wrong parameter for ispin in general_two_rdm_state_av_openmp_work' + print*,'Wrong parameter for ispin in general_two_rdm_state_av_work' print*,'ispin = ',ispin stop endif diff --git a/src/two_body_rdm/orb_range_routines_openmp.irp.f b/src/two_body_rdm/orb_range_routines_openmp.irp.f new file mode 100644 index 00000000..7d791f7c --- /dev/null +++ b/src/two_body_rdm/orb_range_routines_openmp.irp.f @@ -0,0 +1,544 @@ +subroutine orb_range_two_rdm_state_av_openmp(big_array,dim1,norb,list_orb,state_weights,ispin,u_0,N_st,sze) + use bitmasks + implicit none + BEGIN_DOC + ! if ispin == 1 :: alpha/alpha 2rdm + ! == 2 :: beta /beta 2rdm + ! == 3 :: alpha/beta 2rdm + ! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba)) + ! + ! Assumes that the determinants are in psi_det + ! + ! istart, iend, ishift, istep are used in ZMQ parallelization. + END_DOC + integer, intent(in) :: N_st,sze + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + double precision, intent(in) :: u_0(sze,N_st),state_weights(N_st) + + integer :: k + double precision, allocatable :: u_t(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t + allocate(u_t(N_st,N_det)) + do k=1,N_st + call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det) + enddo + call dtranspose( & + u_0, & + size(u_0, 1), & + u_t, & + size(u_t, 1), & + N_det, N_st) + + call orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1) + deallocate(u_t) + + do k=1,N_st + call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) + enddo + +end + +subroutine orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + use bitmasks + implicit none + BEGIN_DOC + ! Computes two-rdm + ! + ! Default should be 1,N_det,0,1 + END_DOC + integer, intent(in) :: N_st,sze,istart,iend,ishift,istep + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + double precision, intent(in) :: u_t(N_st,N_det),state_weights(N_st) + + integer :: k + + PROVIDE N_int + + select case (N_int) + case (1) + call orb_range_two_rdm_state_av_openmp_work_1(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case (2) + call orb_range_two_rdm_state_av_openmp_work_2(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case (3) + call orb_range_two_rdm_state_av_openmp_work_3(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case (4) + call orb_range_two_rdm_state_av_openmp_work_4(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + case default + call orb_range_two_rdm_state_av_openmp_work_N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + end select +end + + + + + BEGIN_TEMPLATE +subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + use bitmasks + implicit none + BEGIN_DOC + ! Computes the two rdm for the N_st vectors |u_t> + ! if ispin == 1 :: alpha/alpha 2rdm + ! == 2 :: beta /beta 2rdm + ! == 3 :: alpha/beta 2rdm + ! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba)) + ! The 2rdm will be computed only on the list of orbitals list_orb, which contains norb + ! In any cases, the state average weights will be used with an array state_weights + ! Default should be 1,N_det,0,1 for istart,iend,ishift,istep + END_DOC + integer, intent(in) :: N_st,sze,istart,iend,ishift,istep + double precision, intent(in) :: u_t(N_st,N_det),state_weights(N_st) + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + + integer :: i,j,k,l + integer :: k_a, k_b, l_a, l_b, m_a, m_b + integer :: istate + integer :: krow, kcol, krow_b, kcol_b + integer :: lrow, lcol + integer :: mrow, mcol + integer(bit_kind) :: spindet($N_int) + integer(bit_kind) :: tmp_det($N_int,2) + integer(bit_kind) :: tmp_det2($N_int,2) + integer(bit_kind) :: tmp_det3($N_int,2) + integer(bit_kind), allocatable :: buffer(:,:) + integer :: n_doubles + integer, allocatable :: doubles(:) + integer, allocatable :: singles_a(:) + integer, allocatable :: singles_b(:) + integer, allocatable :: idx(:), idx0(:) + integer :: maxab, n_singles_a, n_singles_b, kcol_prev + integer*8 :: k8 + double precision :: c_average + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + integer(bit_kind) :: orb_bitmask($N_int) + integer :: list_orb_reverse(mo_num) + integer, allocatable :: keys(:,:) + double precision, allocatable :: values(:) + integer :: nkeys,sze_buff + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + else + print*,'Wrong parameter for ispin in general_two_rdm_state_av_openmp_work' + print*,'ispin = ',ispin + stop + endif + + !do i = 1, N_int + ! det_1_act(i,1) = iand(det_1(i,1),orb_bitmask(i)) + ! det_1_act(i,2) = iand(det_1(i,2),orb_bitmask(i)) + !enddo + + + PROVIDE N_int + + call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) + sze_buff = norb ** 3 + list_orb_reverse = -1000 + do i = 1, norb + list_orb_reverse(list_orb(i)) = i + enddo + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + allocate(idx0(maxab)) + + do i=1,maxab + idx0(i) = i + enddo + + ! Prepare the array of all alpha single excitations + ! ------------------------------------------------- + + PROVIDE N_int nthreads_davidson + !!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) & + ! !$OMP SHARED(psi_bilinear_matrix_rows, N_det, & + ! !$OMP psi_bilinear_matrix_columns, & + ! !$OMP psi_det_alpha_unique, psi_det_beta_unique,& + ! !$OMP n_det_alpha_unique, n_det_beta_unique, N_int,& + ! !$OMP psi_bilinear_matrix_transp_rows, & + ! !$OMP psi_bilinear_matrix_transp_columns, & + ! !$OMP psi_bilinear_matrix_transp_order, N_st, & + ! !$OMP psi_bilinear_matrix_order_transp_reverse, & + ! !$OMP psi_bilinear_matrix_columns_loc, & + ! !$OMP psi_bilinear_matrix_transp_rows_loc, & + ! !$OMP istart, iend, istep, irp_here, v_t, s_t, & + ! !$OMP ishift, idx0, u_t, maxab) & + ! !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,& + ! !$OMP lcol, lrow, l_a, l_b, & + ! !$OMP buffer, doubles, n_doubles, & + ! !$OMP tmp_det2, idx, l, kcol_prev, & + ! !$OMP singles_a, n_singles_a, singles_b, & + ! !$OMP n_singles_b, k8) + + ! Alpha/Beta double excitations + ! ============================= + nkeys = 0 + allocate( keys(4,sze_buff), values(sze_buff)) + allocate( buffer($N_int,maxab), & + singles_a(maxab), & + singles_b(maxab), & + doubles(maxab), & + idx(maxab)) + + kcol_prev=-1 + + ASSERT (iend <= N_det) + ASSERT (istart > 0) + ASSERT (istep > 0) + + !!$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + if (kcol /= kcol_prev) then + call get_all_spin_singles_$N_int( & + psi_det_beta_unique, idx0, & + tmp_det(1,2), N_det_beta_unique, & + singles_b, n_singles_b) + endif + kcol_prev = kcol + + ! Loop over singly excited beta columns + ! ------------------------------------- + + do i=1,n_singles_b + lcol = singles_b(i) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) + + l_a = psi_bilinear_matrix_columns_loc(lcol) + ASSERT (l_a <= N_det) + + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) + + ASSERT (l_a <= N_det) + idx(j) = l_a + l_a = l_a+1 + enddo + j = j-1 + + call get_all_spin_singles_$N_int( & + buffer, idx, tmp_det(1,1), j, & + singles_a, n_singles_a ) + + ! Loop over alpha singles + ! ----------------------- + + if(alpha_beta.or.spin_trace)then + do k = 1,n_singles_a + l_a = singles_a(k) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + if(alpha_beta)then + ! only ONE contribution + if (nkeys+1 .ge. size(values)) then + call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) + nkeys = 0 + endif + else if (spin_trace)then + ! TWO contributions + if (nkeys+2 .ge. size(values)) then + call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) + nkeys = 0 + endif + endif + call orb_range_off_diag_double_to_two_rdm_ab_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + + enddo + endif + + enddo + + enddo + ! !$OMP END DO + + ! !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + + ! Single and double alpha exitations + ! =================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + ! Initial determinant is at k_b in beta-major representation + ! ---------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + spindet(1:$N_int) = tmp_det(1:$N_int,1) + + ! Loop inside the beta column to gather all the connected alphas + lcol = psi_bilinear_matrix_columns(k_a) + l_a = psi_bilinear_matrix_columns_loc(lcol) + do i=1,N_det_alpha_unique + if (l_a > N_det) exit + lcol = psi_bilinear_matrix_columns(l_a) + if (lcol /= kcol) exit + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) + idx(i) = l_a + l_a = l_a+1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles_$N_int( & + buffer, idx, spindet, i, & + singles_a, doubles, n_singles_a, n_doubles ) + + ! Compute Hij for all alpha singles + ! ---------------------------------- + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + do i=1,n_singles_a + l_a = singles_a(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + if(alpha_beta.or.spin_trace.or.alpha_alpha)then + ! increment the alpha/beta part for single excitations +!!!! call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + ! increment the alpha/alpha part for single excitations +!!!! call orb_range_off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + endif + + enddo + + + ! Compute Hij for all alpha doubles + ! ---------------------------------- + + if(alpha_alpha.or.spin_trace)then + do i=1,n_doubles + l_a = doubles(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo +!!!! call orb_range_off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + enddo + endif + + + ! Single and double beta excitations + ! ================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + kcol = psi_bilinear_matrix_columns(k_a) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + spindet(1:$N_int) = tmp_det(1:$N_int,2) + + ! Initial determinant is at k_b in beta-major representation + ! ----------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + ! Loop inside the alpha row to gather all the connected betas + lrow = psi_bilinear_matrix_transp_rows(k_b) + l_b = psi_bilinear_matrix_transp_rows_loc(lrow) + do i=1,N_det_beta_unique + if (l_b > N_det) exit + lrow = psi_bilinear_matrix_transp_rows(l_b) + if (lrow /= krow) exit + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol) + idx(i) = l_b + l_b = l_b+1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles_$N_int( & + buffer, idx, spindet, i, & + singles_b, doubles, n_singles_b, n_doubles ) + + ! Compute Hij for all beta singles + ! ---------------------------------- + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + do i=1,n_singles_b + l_b = singles_b(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) + l_a = psi_bilinear_matrix_transp_order(l_b) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + if(alpha_beta.or.spin_trace.or.beta_beta)then + ! increment the alpha/beta part for single excitations +!!!! call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + ! increment the beta /beta part for single excitations +!!!! call orb_range_off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + endif + enddo + + ! Compute Hij for all beta doubles + ! ---------------------------------- + + if(beta_beta.or.spin_trace)then + do i=1,n_doubles + l_b = doubles(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + l_a = psi_bilinear_matrix_transp_order(l_b) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo +!!!! call orb_range_off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + ASSERT (l_a <= N_det) + + enddo + endif + + + ! Diagonal contribution + ! ===================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + double precision, external :: diag_wee_mat_elem, diag_S_mat_elem + + double precision :: c_1(N_states),c_2(N_states) + c_average = 0.d0 + do l = 1, N_states + c_1(l) = u_t(l,k_a) + c_average += c_1(l) * c_1(l) * state_weights(l) + enddo + + call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) + nkeys = 0 + call orb_range_diag_to_all_two_rdm_dm_buffer(tmp_det,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) + nkeys = 0 + + end do + !!$OMP END DO + deallocate(buffer, singles_a, singles_b, doubles, idx) + !!$OMP END PARALLEL + +end + + SUBST [ N_int ] + + 1;; + 2;; + 3;; + 4;; + N_int;; + + END_TEMPLATE + + +subroutine update_keys_values(keys,values,size_buff,nkeys,dim1,big_array) + implicit none + integer, intent(in) :: size_buff,nkeys,dim1 + integer, intent(in) :: keys(4,size_buff) + double precision, intent(in) :: values(size_buff) + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + + integer :: i,h1,h2,p1,p2 + do i = 1, nkeys + h1 = keys(1,i) + h2 = keys(2,i) + p1 = keys(3,i) + p2 = keys(4,i) + big_array(h1,h2,p1,p2) += values(i) + enddo + +end diff --git a/src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f b/src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f new file mode 100644 index 00000000..ae1c3a54 --- /dev/null +++ b/src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f @@ -0,0 +1,683 @@ + subroutine orb_range_diag_to_all_two_rdm_dm_buffer(det_1,c_1,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the DIAGONAL PART of the two body rdms in a specific range of orbitals for a given determinant det_1 + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: det_1(N_int,2) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1 + double precision, intent(out) :: values(sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,istate + integer(bit_kind) :: det_1_act(N_int,2) + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + do i = 1, N_int + det_1_act(i,1) = iand(det_1(i,1),orb_bitmask(i)) + det_1_act(i,2) = iand(det_1(i,2),orb_bitmask(i)) + enddo + + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call bitstring_to_list_ab(det_1_act, occ, n_occ_ab, N_int) + logical :: is_integer_in_string + integer :: i1,i2 + if(alpha_beta)then + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + values(nkeys) = c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + enddo + enddo + else if (alpha_alpha)then + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(1) + i2 = occ(j,1) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = -0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + else if (beta_beta)then + do i = 1, n_occ_ab(2) + i1 = occ(i,2) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = -0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + else if(spin_trace)then + ! 0.5 * (alpha beta + beta alpha) + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(1) + i2 = occ(j,1) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = -0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + do i = 1, n_occ_ab(2) + i1 = occ(i,2) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = -0.5d0 * c_1 + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + endif + end + + + subroutine orb_range_off_diag_double_to_two_rdm_ab_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC +! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for +! +! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another +! +! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +! +! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +! +! ispin determines which spin-spin component of the two-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 3 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1 + double precision, intent(out) :: values(sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call get_double_excitation(det_1,det_2,exc,phase,N_int) + h1 = exc(1,1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 = exc(1,1,2) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 = exc(1,2,1) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 = exc(1,2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) + if(alpha_beta)then + nkeys += 1 + values(nkeys) = c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + else if(spin_trace)then + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = p1 + keys(2,nkeys) = p2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + endif + end + +! subroutine orb_range_off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) +! use bitmasks +! BEGIN_DOC +!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for +!! +!! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another +!! +!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +!! +!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation +!! +!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +!! +!! ispin determines which spin-spin component of the two-rdm you will update +!! +!! ispin == 1 :: alpha/ alpha +!! ispin == 2 :: beta / beta +!! ispin == 3 :: alpha/ beta +!! ispin == 4 :: spin traced <=> total two-rdm +!! +!! here, only ispin == 3 or 4 will do something +! END_DOC +! implicit none +! integer, intent(in) :: dim1,ispin +! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) +! integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) +! integer(bit_kind), intent(in) :: orb_bitmask(N_int) +! integer, intent(in) :: list_orb_reverse(mo_num) +! double precision, intent(in) :: c_1 +! +! integer :: occ(N_int*bit_kind_size,2) +! integer :: n_occ_ab(2) +! integer :: i,j,h1,h2,istate,p1 +! integer :: exc(0:2,2,2) +! double precision :: phase +! +! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace +! logical :: is_integer_in_string +! alpha_alpha = .False. +! beta_beta = .False. +! alpha_beta = .False. +! spin_trace = .False. +! if( ispin == 1)then +! alpha_alpha = .True. +! else if(ispin == 2)then +! beta_beta = .True. +! else if(ispin == 3)then +! alpha_beta = .True. +! else if(ispin == 4)then +! spin_trace = .True. +! endif +! +! call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) +! call get_single_excitation(det_1,det_2,exc,phase,N_int) +! if(alpha_beta)then +! if (exc(0,1,1) == 1) then +! ! Mono alpha +! h1 = exc(1,1,1) +! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return +! h1 = list_orb_reverse(h1) +! p1 = exc(1,2,1) +! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return +! p1 = list_orb_reverse(p1) +! do i = 1, n_occ_ab(2) +! h2 = occ(i,2) +! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle +! h2 = list_orb_reverse(h2) +! big_array(h1,h2,p1,h2) += c_1 * phase +! enddo +! else +! ! Mono beta +! h1 = exc(1,1,2) +! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return +! h1 = list_orb_reverse(h1) +! p1 = exc(1,2,2) +! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return +! p1 = list_orb_reverse(p1) +! do i = 1, n_occ_ab(1) +! h2 = occ(i,1) +! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle +! h2 = list_orb_reverse(h2) +! big_array(h2,h1,h2,p1) += c_1 * phase +! enddo +! endif +! else if(spin_trace)then +! if (exc(0,1,1) == 1) then +! ! Mono alpha +! h1 = exc(1,1,1) +! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return +! h1 = list_orb_reverse(h1) +! p1 = exc(1,2,1) +! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return +! p1 = list_orb_reverse(p1) +! do i = 1, n_occ_ab(2) +! h2 = occ(i,2) +! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle +! h2 = list_orb_reverse(h2) +! big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase +! big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase +! enddo +! else +! ! Mono beta +! h1 = exc(1,1,2) +! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return +! h1 = list_orb_reverse(h1) +! p1 = exc(1,2,2) +! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return +! p1 = list_orb_reverse(p1) +! do i = 1, n_occ_ab(1) +! h2 = occ(i,1) +! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle +! h2 = list_orb_reverse(h2) +! big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase +! big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase +! enddo +! endif +! endif +! end + +! subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) +! BEGIN_DOC +!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for +!! +!! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another +!! +!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +!! +!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation +!! +!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +!! +!! ispin determines which spin-spin component of the two-rdm you will update +!! +!! ispin == 1 :: alpha/ alpha +!! ispin == 2 :: beta / beta +!! ispin == 3 :: alpha/ beta +!! ispin == 4 :: spin traced <=> total two-rdm +!! +!! here, only ispin == 1 or 4 will do something +! END_DOC +! use bitmasks +! implicit none +! integer, intent(in) :: dim1,ispin +! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) +! integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) +! integer(bit_kind), intent(in) :: orb_bitmask(N_int) +! integer, intent(in) :: list_orb_reverse(mo_num) +! double precision, intent(in) :: c_1 +! +! integer :: occ(N_int*bit_kind_size,2) +! integer :: n_occ_ab(2) +! integer :: i,j,h1,h2,istate,p1 +! integer :: exc(0:2,2,2) +! double precision :: phase +! +! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace +! logical :: is_integer_in_string +! alpha_alpha = .False. +! beta_beta = .False. +! alpha_beta = .False. +! spin_trace = .False. +! if( ispin == 1)then +! alpha_alpha = .True. +! else if(ispin == 2)then +! beta_beta = .True. +! else if(ispin == 3)then +! alpha_beta = .True. +! else if(ispin == 4)then +! spin_trace = .True. +! endif +! +! call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) +! call get_single_excitation(det_1,det_2,exc,phase,N_int) +! if(alpha_alpha.or.spin_trace)then +! if (exc(0,1,1) == 1) then +! ! Mono alpha +! h1 = exc(1,1,1) +! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return +! h1 = list_orb_reverse(h1) +! p1 = exc(1,2,1) +! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return +! p1 = list_orb_reverse(p1) +! do i = 1, n_occ_ab(1) +! h2 = occ(i,1) +! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle +! h2 = list_orb_reverse(h2) +! big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase +! big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase +! +! big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase +! big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase +! enddo +! else +! return +! endif +! endif +! end + +! subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) +! use bitmasks +! BEGIN_DOC +!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for +!! +!! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another +!! +!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +!! +!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation +!! +!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +!! +!! ispin determines which spin-spin component of the two-rdm you will update +!! +!! ispin == 1 :: alpha/ alpha +!! ispin == 2 :: beta / beta +!! ispin == 3 :: alpha/ beta +!! ispin == 4 :: spin traced <=> total two-rdm +!! +!! here, only ispin == 2 or 4 will do something +! END_DOC +! implicit none +! integer, intent(in) :: dim1,ispin +! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) +! integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) +! integer(bit_kind), intent(in) :: orb_bitmask(N_int) +! integer, intent(in) :: list_orb_reverse(mo_num) +! double precision, intent(in) :: c_1 +! +! +! integer :: occ(N_int*bit_kind_size,2) +! integer :: n_occ_ab(2) +! integer :: i,j,h1,h2,istate,p1 +! integer :: exc(0:2,2,2) +! double precision :: phase +! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace +! logical :: is_integer_in_string +! alpha_alpha = .False. +! beta_beta = .False. +! alpha_beta = .False. +! spin_trace = .False. +! if( ispin == 1)then +! alpha_alpha = .True. +! else if(ispin == 2)then +! beta_beta = .True. +! else if(ispin == 3)then +! alpha_beta = .True. +! else if(ispin == 4)then +! spin_trace = .True. +! endif +! +! +! call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) +! call get_single_excitation(det_1,det_2,exc,phase,N_int) +! if(beta_beta.or.spin_trace)then +! if (exc(0,1,1) == 1) then +! return +! else +! ! Mono beta +! h1 = exc(1,1,2) +! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return +! h1 = list_orb_reverse(h1) +! p1 = exc(1,2,2) +! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return +! p1 = list_orb_reverse(p1) +! do istate = 1, N_states +! do i = 1, n_occ_ab(2) +! h2 = occ(i,2) +! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle +! h2 = list_orb_reverse(h2) +! big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase +! big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase +! +! big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase +! big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase +! enddo +! enddo +! endif +! endif +! end + + +! subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) +! use bitmasks +! BEGIN_DOC +!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for +!! +!! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another +!! +!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +!! +!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation +!! +!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +!! +!! ispin determines which spin-spin component of the two-rdm you will update +!! +!! ispin == 1 :: alpha/ alpha +!! ispin == 2 :: beta / beta +!! ispin == 3 :: alpha/ beta +!! ispin == 4 :: spin traced <=> total two-rdm +!! +!! here, only ispin == 1 or 4 will do something +! END_DOC +! implicit none +! integer, intent(in) :: dim1,ispin +! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) +! integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) +! integer(bit_kind), intent(in) :: orb_bitmask(N_int) +! integer, intent(in) :: list_orb_reverse(mo_num) +! double precision, intent(in) :: c_1 +! +! integer :: i,j,h1,h2,p1,p2,istate +! integer :: exc(0:2,2) +! double precision :: phase +! +! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace +! logical :: is_integer_in_string +! alpha_alpha = .False. +! beta_beta = .False. +! alpha_beta = .False. +! spin_trace = .False. +! if( ispin == 1)then +! alpha_alpha = .True. +! else if(ispin == 2)then +! beta_beta = .True. +! else if(ispin == 3)then +! alpha_beta = .True. +! else if(ispin == 4)then +! spin_trace = .True. +! endif +! call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) +! h1 =exc(1,1) +! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return +! h1 = list_orb_reverse(h1) +! h2 =exc(2,1) +! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return +! h2 = list_orb_reverse(h2) +! p1 =exc(1,2) +! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return +! p1 = list_orb_reverse(p1) +! p2 =exc(2,2) +! if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return +! p2 = list_orb_reverse(p2) +! if(alpha_alpha.or.spin_trace)then +! do istate = 1, N_states +! big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase +! big_array(h1,h2,p2,p1) -= 0.5d0 * c_1 * phase +! +! big_array(h2,h1,p2,p1) += 0.5d0 * c_1 * phase +! big_array(h2,h1,p1,p2) -= 0.5d0 * c_1 * phase +! enddo +! endif +! end + +! subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) +! use bitmasks +! BEGIN_DOC +!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for +!! +!! a given couple of determinant det_1, det_2 being a BETA /BETA DOUBLE excitation with respect to one another +!! +!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +!! +!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation +!! +!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +!! +!! ispin determines which spin-spin component of the two-rdm you will update +!! +!! ispin == 1 :: alpha/ alpha +!! ispin == 2 :: beta / beta +!! ispin == 3 :: alpha/ beta +!! ispin == 4 :: spin traced <=> total two-rdm +!! +!! here, only ispin == 2 or 4 will do something +! END_DOC +! implicit none +! +! integer, intent(in) :: dim1,ispin +! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) +! integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) +! integer(bit_kind), intent(in) :: orb_bitmask(N_int) +! integer, intent(in) :: list_orb_reverse(mo_num) +! double precision, intent(in) :: c_1 +! +! integer :: i,j,h1,h2,p1,p2,istate +! integer :: exc(0:2,2) +! double precision :: phase +! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace +! logical :: is_integer_in_string +! alpha_alpha = .False. +! beta_beta = .False. +! alpha_beta = .False. +! spin_trace = .False. +! if( ispin == 1)then +! alpha_alpha = .True. +! else if(ispin == 2)then +! beta_beta = .True. +! else if(ispin == 3)then +! alpha_beta = .True. +! else if(ispin == 4)then +! spin_trace = .True. +! endif +! +! call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) +! h1 =exc(1,1) +! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return +! h1 = list_orb_reverse(h1) +! h2 =exc(2,1) +! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return +! h2 = list_orb_reverse(h2) +! p1 =exc(1,2) +! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return +! p1 = list_orb_reverse(p1) +! p2 =exc(2,2) +! if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return +! p2 = list_orb_reverse(p2) +! if(beta_beta.or.spin_trace)then +! big_array(h1,h2,p1,p2) += 0.5d0 * c_1* phase +! big_array(h1,h2,p2,p1) -= 0.5d0 * c_1* phase +! +! big_array(h2,h1,p2,p1) += 0.5d0 * c_1* phase +! big_array(h2,h1,p1,p2) -= 0.5d0 * c_1* phase +! endif +! end + diff --git a/src/two_body_rdm/two_rdm.irp.f b/src/two_body_rdm/two_rdm.irp.f index 06a8e1e6..c162f365 100644 --- a/src/two_body_rdm/two_rdm.irp.f +++ b/src/two_body_rdm/two_rdm.irp.f @@ -19,7 +19,7 @@ two_rdm_beta_beta_mo = 0.d0 print*,'providing two_rdm_alpha_beta ...' call wall_time(cpu_0) - call all_two_rdm_dm_nstates_openmp(two_rdm_alpha_alpha_mo,two_rdm_beta_beta_mo,two_rdm_alpha_beta_mo,dim1,dim2,dim3,dim4,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call all_two_rdm_dm_nstates(two_rdm_alpha_alpha_mo,two_rdm_beta_beta_mo,two_rdm_alpha_beta_mo,dim1,dim2,dim3,dim4,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(cpu_1) print*,'two_rdm_alpha_beta provided in',dabs(cpu_1-cpu_0)