From 419ed79c49a89382f3e473df647c9624a4f3e759 Mon Sep 17 00:00:00 2001 From: eginer Date: Sat, 10 Feb 2024 12:48:29 +0100 Subject: [PATCH] added transition two rdm --- src/davidson/u0_wee_u0.irp.f | 22 + src/two_body_rdm/act_2_transition_rdm.irp.f | 39 + src/two_body_rdm/example.irp.f | 88 ++ src/two_body_rdm/io_two_rdm.irp.f | 34 + src/two_body_rdm/test_2_rdm.irp.f | 1 + .../davidson_like_trans_2rdm.irp.f | 585 ++++++++++ src/two_rdm_routines/update_trans_rdm.irp.f | 1002 +++++++++++++++++ 7 files changed, 1771 insertions(+) create mode 100644 src/two_body_rdm/act_2_transition_rdm.irp.f create mode 100644 src/two_rdm_routines/davidson_like_trans_2rdm.irp.f create mode 100644 src/two_rdm_routines/update_trans_rdm.irp.f diff --git a/src/davidson/u0_wee_u0.irp.f b/src/davidson/u0_wee_u0.irp.f index 0c543aca..bd3525e1 100644 --- a/src/davidson/u0_wee_u0.irp.f +++ b/src/davidson/u0_wee_u0.irp.f @@ -492,3 +492,25 @@ subroutine u_0_H_u_0_two_e(e_0,u_0,n,keys_tmp,Nint,N_st,sze) deallocate (s_0, v_0) end +BEGIN_PROVIDER [double precision, psi_energy_two_e_trans, (N_states, N_states)] + implicit none + BEGIN_DOC +! psi_energy_two_e_trans(istate,jstate) = + END_dOC + integer :: i,j,istate,jstate + double precision :: hij, coef_i, coef_j + psi_energy_two_e_trans = 0.d0 + do i = 1, N_det + do j = 1, N_det + call i_H_j_two_e(psi_det(1,1,i),psi_det(1,1,j),N_int,hij) + do istate = 1, N_states + coef_i = psi_coef(i,istate) + do jstate = 1, N_states + coef_j = psi_coef(j,jstate) + psi_energy_two_e_trans(jstate,istate) += coef_i * coef_j * hij + enddo + enddo + enddo + enddo + +END_PROVIDER diff --git a/src/two_body_rdm/act_2_transition_rdm.irp.f b/src/two_body_rdm/act_2_transition_rdm.irp.f new file mode 100644 index 00000000..3d08b084 --- /dev/null +++ b/src/two_body_rdm/act_2_transition_rdm.irp.f @@ -0,0 +1,39 @@ + BEGIN_PROVIDER [double precision, act_2_rdm_trans_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states,N_states)] + implicit none + BEGIN_DOC +! act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2rdm_trans +! +! \sum_{\sigma,\sigma'} +! +! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" +! +! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec}^{act} * (N_{elec}^{act} - 1) +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" + END_DOC + integer :: ispin + double precision :: wall_1, wall_2 + ! condition for beta/beta spin + print*,'' + print*,'Providing act_2_rdm_trans_spin_trace_mo ' + character*(128) :: name_file + name_file = 'act_2_rdm_trans_spin_trace_mo' + ispin = 4 + act_2_rdm_trans_spin_trace_mo = 0.d0 + call wall_time(wall_1) +! if(read_two_body_rdm_trans_spin_trace)then +! print*,'Reading act_2_rdm_trans_spin_trace_mo from disk ...' +! call read_array_two_rdm_trans(n_act_orb,N_states,act_2_rdm_trans_spin_trace_mo,name_file) +! else + call orb_range_2_trans_rdm_openmp(act_2_rdm_trans_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) +! endif +! if(write_two_body_rdm_trans_spin_trace)then +! print*,'Writing act_2_rdm_trans_spin_trace_mo on disk ...' +! call write_array_two_rdm_trans(n_act_orb,n_states,act_2_rdm_trans_spin_trace_mo,name_file) +! call ezfio_set_two_body_rdm_trans_io_two_body_rdm_trans_spin_trace("Read") +! endif + + act_2_rdm_trans_spin_trace_mo *= 2.d0 + call wall_time(wall_2) + print*,'Wall time to provide act_2_rdm_trans_spin_trace_mo',wall_2 - wall_1 + END_PROVIDER diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 30e2685a..38510fe9 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -365,3 +365,91 @@ subroutine routine_full_mos end + +subroutine routine_active_only_trans + implicit none + integer :: i,j,k,l,iorb,jorb,korb,lorb,istate,jstate + 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_tot(N_states,N_states),rdm_transtot + double precision :: spin_trace + double precision :: accu_tot + + wee_tot = 0.d0 + + + iorb = 1 + jorb = 1 + korb = 1 + lorb = 1 + vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + provide act_2_rdm_trans_spin_trace_mo + i = 1 + j = 2 + + print*,'**************************' + print*,'**************************' + do jstate = 1, N_states + do istate = 1, N_states + !! PURE ACTIVE PART + !! + accu_tot = 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) + ! 1 2 1 2 2 1 2 1 +! if(dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(j,i,l,k,istate,jstate)).gt.1.d-10)then +! print*,'Error in act_2_rdm_trans_spin_trace_mo' +! print*,"dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l) - act_2_rdm_trans_spin_trace_mo(j,i,l,k)).gt.1.d-10" +! print*,i,j,k,l +! print*,act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate),act_2_rdm_trans_spin_trace_mo(j,i,l,k,istate,jstate),dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(j,i,l,k,istate,jstate)) +! endif + + ! 1 2 1 2 1 2 1 2 +! if(dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate)).gt.1.d-10)then +! print*,'Error in act_2_rdm_trans_spin_trace_mo' +! print*,"dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate)).gt.1.d-10" +! print*,i,j,k,l +! print*,act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate),act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate),dabs(act_2_rdm_trans_spin_trace_mo(i,j,k,l,istate,jstate) - act_2_rdm_trans_spin_trace_mo(k,l,i,j,istate,jstate)) +! endif + + vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + + + rdm_transtot = act_2_rdm_trans_spin_trace_mo(l,k,j,i,istate,jstate) + + wee_tot(istate,jstate) += 0.5d0 * vijkl * rdm_transtot + + enddo + enddo + enddo + enddo + print*,'' + print*,'' + print*,'Active space only energy for state ',istate,jstate + print*,'wee_tot = ',wee_tot(istate,jstate) + print*,'Full energy ' + print*,'psi_energy_two_e(istate,jstate)= ',psi_energy_two_e_trans(istate,jstate) + print*,'--------------------------' + enddo + enddo + print*,'Wee from DM ' + do istate = 1,N_states + write(*,'(100(F16.10,X))')wee_tot(1:N_states,istate) + enddo + print*,'Wee from Psi det' + do istate = 1,N_states + write(*,'(100(F16.10,X))')psi_energy_two_e_trans(1:N_states,istate) + enddo + +end + diff --git a/src/two_body_rdm/io_two_rdm.irp.f b/src/two_body_rdm/io_two_rdm.irp.f index bdd8a4f9..0b30d76f 100644 --- a/src/two_body_rdm/io_two_rdm.irp.f +++ b/src/two_body_rdm/io_two_rdm.irp.f @@ -31,3 +31,37 @@ subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file) close(unit=i_unit_output) end + +subroutine write_array_two_trans_rdm(n_orb,nstates,array_tmp,name_file) + implicit none + integer, intent(in) :: n_orb,nstates + character*(128), intent(in) :: name_file + double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,nstates,nstates) + + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + PROVIDE ezfio_filename + output=trim(ezfio_filename)//'/work/'//trim(name_file) + i_unit_output = getUnitAndOpen(output,'W') + call lock_io() + write(i_unit_output)array_tmp + call unlock_io() + close(unit=i_unit_output) +end + +subroutine read_array_two_trans_rdm(n_orb,nstates,array_tmp,name_file) + implicit none + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + integer, intent(in) :: n_orb,nstates + character*(128), intent(in) :: name_file + double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,N_states,nstates) + PROVIDE ezfio_filename + output=trim(ezfio_filename)//'/work/'//trim(name_file) + i_unit_output = getUnitAndOpen(output,'R') + call lock_io() + read(i_unit_output)array_tmp + call unlock_io() + close(unit=i_unit_output) +end + diff --git a/src/two_body_rdm/test_2_rdm.irp.f b/src/two_body_rdm/test_2_rdm.irp.f index 123261d8..de2606a7 100644 --- a/src/two_body_rdm/test_2_rdm.irp.f +++ b/src/two_body_rdm/test_2_rdm.irp.f @@ -4,5 +4,6 @@ program test_2_rdm touch read_wf call routine_active_only call routine_full_mos + call routine_active_only_trans end diff --git a/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f b/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f new file mode 100644 index 00000000..9e68a0e1 --- /dev/null +++ b/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f @@ -0,0 +1,585 @@ +subroutine orb_range_2_trans_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sze) + use bitmasks + implicit none + BEGIN_DOC + ! if ispin == 1 :: alpha/alpha 2_rdm + ! == 2 :: beta /beta 2_rdm + ! == 3 :: alpha/beta + beta/alpha 2trans_rdm + ! == 4 :: spin traced 2_rdm :: aa + bb + ab + ba + ! + ! notice that here it is the TRANSITION RDM THAT IS COMPUTED + ! + ! THE DIAGONAL PART IS THE USUAL ONE FOR A GIVEN STATE + ! 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,N_st,N_st) + double precision, intent(in) :: u_0(sze,N_st) + + integer :: k + double precision, allocatable :: u_t(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t + PROVIDE mo_two_e_integrals_in_map + 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_2_trans_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 + call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) + enddo + +end + +subroutine orb_range_2_trans_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 + ! Computes two-trans_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,N_st,N_st) + double precision, intent(in) :: u_t(N_st,N_det) + + integer :: k + + PROVIDE N_int + + select case (N_int) + case (1) + call orb_range_2_trans_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_2_trans_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_2_trans_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_2_trans_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_2_trans_rdm_openmp_work_N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + end select +end + + + BEGIN_TEMPLATE +subroutine orb_range_2_trans_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 + BEGIN_DOC + ! Computes the two trans_rdm for the N_st vectors |u_t> + ! if ispin == 1 :: alpha/alpha 2trans_rdm + ! == 2 :: beta /beta 2trans_rdm + ! == 3 :: alpha/beta 2trans_rdm + ! == 4 :: spin traced 2trans_rdm :: aa + bb + 0.5 (ab + ba)) + ! The 2trans_rdm will be computed only on the list of orbitals list_orb, which contains norb + ! 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) + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st) + + integer(omp_lock_kind) :: lock_2trans_rdm + integer :: i,j,k,l + integer :: k_a, k_b, l_a, l_b + integer :: krow, kcol + integer :: lrow, lcol + integer(bit_kind) :: spindet($N_int) + integer(bit_kind) :: tmp_det($N_int,2) + integer(bit_kind) :: tmp_det2($N_int,2) + integer(bit_kind) :: tmp_det3($N_int,2) + integer(bit_kind), allocatable :: buffer(:,:) + integer :: n_doubles + integer, allocatable :: doubles(:) + integer, allocatable :: singles_a(:) + integer, allocatable :: singles_b(:) + integer, allocatable :: idx(:), idx0(:) + integer :: maxab, n_singles_a, n_singles_b, kcol_prev + + 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 + integer :: ll + 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_2_trans_rdm_state_av_openmp_work' + print*,'ispin = ',ispin + stop + endif + + + PROVIDE N_int + + call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) + sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 + list_orb_reverse = -1000 + do i = 1, norb + list_orb_reverse(list_orb(i)) = i + enddo + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + allocate(idx0(maxab)) + + do i=1,maxab + idx0(i) = i + enddo + call omp_init_lock(lock_2trans_rdm) + + ! Prepare the array of all alpha single excitations + ! ------------------------------------------------- + + PROVIDE N_int nthreads_davidson elec_alpha_num + !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) & + !$OMP SHARED(psi_bilinear_matrix_rows, N_det,lock_2trans_rdm,& + !$OMP psi_bilinear_matrix_columns, & + !$OMP psi_det_alpha_unique, psi_det_beta_unique,& + !$OMP n_det_alpha_unique, n_det_beta_unique, N_int,& + !$OMP psi_bilinear_matrix_transp_rows, & + !$OMP psi_bilinear_matrix_transp_columns, & + !$OMP psi_bilinear_matrix_transp_order, N_st, & + !$OMP psi_bilinear_matrix_order_transp_reverse, & + !$OMP psi_bilinear_matrix_columns_loc, & + !$OMP psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, & + !$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, & + !$OMP ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask) & + !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, & + !$OMP lcol, lrow, l_a, l_b, & + !$OMP buffer, doubles, n_doubles, & + !$OMP tmp_det2, idx, l, kcol_prev, & + !$OMP singles_a, n_singles_a, singles_b, & + !$OMP n_singles_b, nkeys, keys, values) + + ! Alpha/Beta double excitations + ! ============================= + nkeys = 0 + allocate( keys(4,sze_buff), values(n_st,n_st,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) +! print*,'nkeys before = ',nkeys + do ll = 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if(alpha_beta)then + ! only ONE contribution + if (nkeys+1 .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + else if (spin_trace)then + ! TWO contributions + if (nkeys+2 .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + endif + call orb_range_off_diag_double_to_all_states_ab_trans_rdm_buffer(tmp_det,tmp_det2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + + enddo + endif + + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + 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) + do ll= 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if(alpha_beta.or.spin_trace.or.alpha_alpha)then + ! increment the alpha/beta part for single excitations + if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_ab_trans_rdm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + ! increment the alpha/alpha part for single excitations + if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_aa_trans_rdm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + endif + + enddo + + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + ! 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) + + do ll= 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if (nkeys+4 .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_double_to_all_states_aa_trans_rdm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + enddo + endif + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + + ! 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) + do ll= 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if(alpha_beta.or.spin_trace.or.beta_beta)then + ! increment the alpha/beta part for single excitations + if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_ab_trans_rdm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + ! increment the beta /beta part for single excitations + if (nkeys+4 * elec_alpha_num .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_bb_trans_rdm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + endif + enddo + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + ! 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) + do ll= 1, N_states + do l= 1, N_states + c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) + enddo + enddo + if (nkeys+4 .ge. sze_buff) then + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + endif + call orb_range_off_diag_double_to_all_states_trans_rdm_bb_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) +! print*,'to do orb_range_off_diag_double_to_2_trans_rdm_bb_dm_buffer' + ASSERT (l_a <= N_det) + + enddo + endif + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + + ! 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,N_states) + do ll = 1, N_states + do l = 1, N_states + c_1(l,ll) = u_t(ll,k_a) * u_t(l,k_a) + enddo + enddo + + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + call orb_range_diag_to_all_states_2_rdm_trans_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) + nkeys = 0 + + end do + !$OMP END DO + deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values) + !$OMP END PARALLEL + +end + + SUBST [ N_int ] + + 1;; + 2;; + 3;; + 4;; + N_int;; + + END_TEMPLATE + +subroutine update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + use omp_lib + implicit none + integer, intent(in) :: n_st,nkeys,dim1 + integer, intent(in) :: keys(4,nkeys) + double precision, intent(in) :: values(n_st,n_st,nkeys) + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,n_st,n_st) + + integer(omp_lock_kind),intent(inout):: lock_2rdm + + integer :: i,h1,h2,p1,p2,istate,jstate + call omp_set_lock(lock_2rdm) + +! print*,'*************' +! print*,'updating' +! print*,'nkeys',nkeys + do i = 1, nkeys + h1 = keys(1,i) + h2 = keys(2,i) + p1 = keys(3,i) + p2 = keys(4,i) + do jstate = 1, N_st + do istate = 1, N_st +!! print*,h1,h2,p1,p2,values(istate,i) + big_array(h1,h2,p1,p2,istate,jstate) += values(istate,jstate,i) + enddo + enddo + enddo + call omp_unset_lock(lock_2rdm) + +end + diff --git a/src/two_rdm_routines/update_trans_rdm.irp.f b/src/two_rdm_routines/update_trans_rdm.irp.f new file mode 100644 index 00000000..9f7077a2 --- /dev/null +++ b/src/two_rdm_routines/update_trans_rdm.irp.f @@ -0,0 +1,1002 @@ + subroutine orb_range_diag_to_all_states_2_rdm_trans_buffer(det_1,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for a given determinant det_1 + ! + ! c_1 is the array of the contributions to the trans_rdm for all states + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: det_1(N_int,2) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2 + integer(bit_kind) :: det_1_act(N_int,2) + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + do i = 1, N_int + det_1_act(i,1) = iand(det_1(i,1),orb_bitmask(i)) + det_1_act(i,2) = iand(det_1(i,2),orb_bitmask(i)) + enddo + + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call bitstring_to_list_ab(det_1_act, occ, n_occ_ab, N_int) + logical :: is_integer_in_string + integer :: i1,i2,istate + integer :: jstate + if(alpha_beta)then + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + ! If alpha/beta, electron 1 is alpha, electron 2 is beta + ! Therefore you don't necessayr have symmetry between electron 1 and 2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + + else if (alpha_alpha)then + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(1) + i2 = occ(j,1) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + else if (beta_beta)then + do i = 1, n_occ_ab(2) + i1 = occ(i,2) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + else if(spin_trace)then + ! 0.5 * (alpha beta + beta alpha) + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + do i = 1, n_occ_ab(1) + i1 = occ(i,1) + do j = 1, n_occ_ab(1) + i2 = occ(j,1) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + do i = 1, n_occ_ab(2) + i1 = occ(i,2) + do j = 1, n_occ_ab(2) + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate) + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 + enddo + enddo + endif + end + + + subroutine orb_range_off_diag_double_to_all_states_ab_trans_rdm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC +! routine that update the OFF DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for +! +! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another +! +! c_1 is the array of the contributions to the trans_rdm for all states +! +! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +! +! ispin determines which spin-spin component of the two-trans_rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-trans_rdm +! +! here, only ispin == 3 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + integer :: jstate + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call get_double_excitation(det_1,det_2,exc,phase,N_int) + h1 = exc(1,1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 = exc(1,1,2) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 = exc(1,2,1) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 = exc(1,2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) + if(alpha_beta)then + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + else if(spin_trace)then + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + endif + end + + subroutine orb_range_off_diag_single_to_all_states_ab_trans_rdm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_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 trans_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-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 3 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: jstate + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,p1,istate + integer :: exc(0:2,2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + 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) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + else + ! Mono beta + h1 = exc(1,1,2) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + 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) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + else + ! Mono beta + h1 = exc(1,1,2) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + enddo + endif + endif + end + + subroutine orb_range_off_diag_single_to_all_states_aa_trans_rdm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_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 trans_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-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 1 or 4 will do something + END_DOC + use bitmasks + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: jstate + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,p1,istate + integer :: exc(0:2,2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + 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) + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + else + return + endif + endif + end + + subroutine orb_range_off_diag_single_to_all_states_bb_trans_rdm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_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 trans_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-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 2 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: jstate + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,p1,istate + integer :: exc(0:2,2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + 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 i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + endif + endif + end + + + subroutine orb_range_off_diag_double_to_all_states_aa_trans_rdm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_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 trans_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-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 1 or 4 will do something + END_DOC + implicit none + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2) + double precision :: phase + + integer :: jstate + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) + h1 =exc(1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 =exc(2,1) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 =exc(1,2) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 =exc(2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) + if(alpha_alpha.or.spin_trace)then + nkeys += 1 + + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + endif + end + + subroutine orb_range_off_diag_double_to_all_states_trans_rdm_bb_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body trans_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 trans_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-trans_rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-trans_rdm + ! + ! here, only ispin == 2 or 4 will do something + END_DOC + implicit none + + integer, intent(in) :: ispin,sze_buff,N_st + integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) + double precision, intent(in) :: c_1(N_st,N_st) + double precision, intent(out) :: values(N_st,N_st,sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: jstate + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + + call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) + h1 =exc(1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 =exc(2,1) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 =exc(1,2) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 =exc(2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) + if(beta_beta.or.spin_trace)then + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + do jstate = 1, N_st + do istate = 1, N_st + values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase + enddo + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + endif + end +