From 3e0ada95387d4c606087e9494349d33528314c83 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 18 Mar 2020 15:13:49 +0100 Subject: [PATCH] beginning the cleaning of two_body_rdm --- src/density_for_dft/EZFIO.cfg | 6 +- src/density_for_dft/density_for_dft.irp.f | 4 +- src/two_body_rdm/README.rst | 1 + src/two_body_rdm/ab_only_routines.irp.f | 402 --------- src/two_body_rdm/all_2rdm_routines.irp.f | 442 ---------- ...> all_states_act_2_rdm_dav_routines.irp.f} | 0 ....irp.f => all_states_act_2_rdm_prov.irp.f} | 34 +- ...ll_states_act_2_rdm_update_routines.irp.f} | 0 src/two_body_rdm/compute.irp.f | 269 ------ src/two_body_rdm/compute_orb_range_omp.irp.f | 807 ------------------ src/two_body_rdm/orb_range_omp.irp.f | 85 -- src/two_body_rdm/orb_range_routines_omp.irp.f | 568 ------------ ... => state_av_act_2_rdm_dav_routines.irp.f} | 0 ...ge.irp.f => state_av_act_2_rdm_prov.irp.f} | 32 +- ... state_av_act_2_rdm_update_routines.irp.f} | 0 src/two_body_rdm/two_rdm.irp.f | 62 -- 16 files changed, 50 insertions(+), 2662 deletions(-) delete mode 100644 src/two_body_rdm/ab_only_routines.irp.f delete mode 100644 src/two_body_rdm/all_2rdm_routines.irp.f rename src/two_body_rdm/{all_states_routines.irp.f => all_states_act_2_rdm_dav_routines.irp.f} (100%) rename src/two_body_rdm/{all_states_2_rdm.irp.f => all_states_act_2_rdm_prov.irp.f} (64%) rename src/two_body_rdm/{compute_all_states.irp.f => all_states_act_2_rdm_update_routines.irp.f} (100%) delete mode 100644 src/two_body_rdm/compute.irp.f delete mode 100644 src/two_body_rdm/compute_orb_range_omp.irp.f delete mode 100644 src/two_body_rdm/orb_range_omp.irp.f delete mode 100644 src/two_body_rdm/orb_range_routines_omp.irp.f rename src/two_body_rdm/{orb_range_routines.irp.f => state_av_act_2_rdm_dav_routines.irp.f} (100%) rename src/two_body_rdm/{orb_range.irp.f => state_av_act_2_rdm_prov.irp.f} (68%) rename src/two_body_rdm/{compute_orb_range.irp.f => state_av_act_2_rdm_update_routines.irp.f} (100%) delete mode 100644 src/two_body_rdm/two_rdm.irp.f diff --git a/src/density_for_dft/EZFIO.cfg b/src/density_for_dft/EZFIO.cfg index 42b8eab4..63d6bc08 100644 --- a/src/density_for_dft/EZFIO.cfg +++ b/src/density_for_dft/EZFIO.cfg @@ -11,10 +11,10 @@ interface: ezfio,provider,ocaml default: 0.5 [no_core_density] -type: character*(32) -doc: Type of density. If [no_core_dm] then all elements of the density matrix involving at least one orbital set as core are set to zero +type: logical +doc: If [no_core_density] then all elements of the density matrix involving at least one orbital set as core are set to zero. The default is False in order to take all the density. interface: ezfio, provider, ocaml -default: full_density +default: False [normalize_dm] type: logical diff --git a/src/density_for_dft/density_for_dft.irp.f b/src/density_for_dft/density_for_dft.irp.f index c925bdf8..ee70cd38 100644 --- a/src/density_for_dft/density_for_dft.irp.f +++ b/src/density_for_dft/density_for_dft.irp.f @@ -22,7 +22,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_alpha_for_dft, (mo_num,mo_num, N_s one_e_dm_mo_alpha_for_dft(:,:,1) = one_e_dm_mo_alpha_average(:,:) endif - if(no_core_density .EQ. "no_core_dm")then + if(no_core_density)then integer :: ii,i,j do ii = 1, n_core_orb i = list_core(ii) @@ -73,7 +73,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st one_e_dm_mo_beta_for_dft(:,:,1) = one_e_dm_mo_beta_average(:,:) endif - if(no_core_density .EQ. "no_core_dm")then + if(no_core_density)then integer :: ii,i,j do ii = 1, n_core_orb i = list_core(ii) diff --git a/src/two_body_rdm/README.rst b/src/two_body_rdm/README.rst index 978240c9..69f1f2a4 100644 --- a/src/two_body_rdm/README.rst +++ b/src/two_body_rdm/README.rst @@ -6,3 +6,4 @@ Contains the two rdms $\alpha\alpha$, $\beta\beta$ and $\alpha\beta$ stored as arrays, with pysicists notation, consistent with the two-electron integrals in the MO basis. + diff --git a/src/two_body_rdm/ab_only_routines.irp.f b/src/two_body_rdm/ab_only_routines.irp.f deleted file mode 100644 index fb3c421c..00000000 --- a/src/two_body_rdm/ab_only_routines.irp.f +++ /dev/null @@ -1,402 +0,0 @@ - - subroutine two_rdm_ab_nstates(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze) - use bitmasks - implicit none - BEGIN_DOC - ! Computes the alpha/beta part of the two-body density matrix IN CHEMIST NOTATIONS - ! - ! 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,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: u_0(sze,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 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 - call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) - enddo - - end - - - 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 - ! Computes the alpha/beta part of the two-body density matrix - ! - ! Default should be 1,N_det,0,1 - END_DOC - integer, intent(in) :: N_st,sze,istart,iend,ishift,istep - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - double precision, intent(in) :: u_t(N_st,N_det) - - - PROVIDE N_int - - select case (N_int) - case (1) - 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_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) - case (3) - 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_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) - case default - 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_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 - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - double precision, intent(in) :: u_t(N_st,N_det) - - double precision :: hij, sij - 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, nmax - integer*8 :: k8 - - 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 - - ! Alpha/Beta double excitations - ! ============================= - - 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) - - 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 - ! ----------------------- - - 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) - !!!!!!!!!!!!!!!!!! ALPHA BETA - do l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - enddo - call off_diagonal_double_to_two_rdm_ab_dm(tmp_det,tmp_det2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) - enddo - - enddo - - enddo - - - do k_a=istart+ishift,iend,istep - - - ! Single and double alpha excitations - ! =================================== - - - ! 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) - - 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) - !!!! MONO SPIN - do l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - enddo - call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) - - enddo - - - !! Compute Hij for all alpha doubles - !! ---------------------------------- - ! - !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) - - ! call i_H_j_double_spin_erf( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij) - ! do l=1,N_st - ! v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a) - ! ! same spin => sij = 0 - ! enddo - !enddo - - - - ! 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) - - ! 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 l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - enddo - call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) - ASSERT (l_a <= N_det) - enddo - ! - !! Compute Hij for all beta doubles - !! ---------------------------------- - ! - !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) - - ! call i_H_j_double_spin_erf( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij) - ! l_a = psi_bilinear_matrix_transp_order(l_b) - ! ASSERT (l_a <= N_det) - - ! do l=1,N_st - ! v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a) - ! ! same spin => sij = 0 - ! enddo - !enddo - - - ! 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_H_mat_elem_erf, diag_S_mat_elem - double precision :: c_1(N_states),c_2(N_states) - do l = 1, N_states - c_1(l) = u_t(l,k_a) - enddo - - call diagonal_contrib_to_two_rdm_ab_dm(tmp_det,c_1,big_array,dim1,dim2,dim3,dim4) - - end do - deallocate(buffer, singles_a, singles_b, doubles, idx) - - end - - SUBST [ N_int ] - - 1;; - 2;; - 3;; - 4;; - N_int;; - - END_TEMPLATE diff --git a/src/two_body_rdm/all_2rdm_routines.irp.f b/src/two_body_rdm/all_2rdm_routines.irp.f deleted file mode 100644 index fa036e6a..00000000 --- a/src/two_body_rdm/all_2rdm_routines.irp.f +++ /dev/null @@ -1,442 +0,0 @@ -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 - ! Computes the alpha/alpha, beta/beta and alpha/beta part of the two-body density matrix IN CHEMIST NOTATIONS - ! - ! 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,dim2,dim3,dim4 - double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: u_0(sze,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 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 - call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) - enddo - -end - - -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 - ! 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,dim2,dim3,dim4 - double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states) - double precision, intent(in) :: u_t(N_st,N_det) - - - PROVIDE N_int - - select case (N_int) - case (1) - 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_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_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_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_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_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 - ! Computes $v_t = H | u_t \\rangle$ and $s_t = S^2 | u_t \\rangle$ - ! - ! Default should be 1,N_det,0,1 - 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,dim2,dim3,dim4 - double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states) - - 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 - - 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 - ! ============================= - - 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 - ! ----------------------- - - 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) - !call i_H_j_double_alpha_beta(tmp_det,tmp_det2,$N_int,hij) - do l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - enddo - call off_diagonal_double_to_two_rdm_ab_dm(tmp_det,tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4) - enddo - - 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 l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - enddo - ! increment the alpha/beta part for single excitations - call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4) - ! increment the alpha/alpha part for single excitations - call off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_1,c_2,big_array_aa,dim1,dim2,dim3,dim4) - - enddo - - - ! Compute Hij for all alpha doubles - ! ---------------------------------- - - 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 l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - enddo - call off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,c_2,big_array_aa,dim1,dim2,dim3,dim4) - enddo - - - ! 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 l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - enddo - ! increment the alpha/beta part for single excitations - call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4) - ! increment the beta /beta part for single excitations - call off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4) - enddo - - ! Compute Hij for all beta doubles - ! ---------------------------------- - - 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 l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - enddo - call off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4) - ASSERT (l_a <= N_det) - - enddo - - - ! 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) - do l = 1, N_states - c_1(l) = u_t(l,k_a) - enddo - - call diagonal_contrib_to_all_two_rdm_dm(tmp_det,c_1,big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4) - - 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 - diff --git a/src/two_body_rdm/all_states_routines.irp.f b/src/two_body_rdm/all_states_act_2_rdm_dav_routines.irp.f similarity index 100% rename from src/two_body_rdm/all_states_routines.irp.f rename to src/two_body_rdm/all_states_act_2_rdm_dav_routines.irp.f diff --git a/src/two_body_rdm/all_states_2_rdm.irp.f b/src/two_body_rdm/all_states_act_2_rdm_prov.irp.f similarity index 64% rename from src/two_body_rdm/all_states_2_rdm.irp.f rename to src/two_body_rdm/all_states_act_2_rdm_prov.irp.f index bc503223..fc6e4224 100644 --- a/src/two_body_rdm/all_states_2_rdm.irp.f +++ b/src/two_body_rdm/all_states_act_2_rdm_prov.irp.f @@ -5,8 +5,11 @@ implicit none double precision, allocatable :: state_weights(:) BEGIN_DOC -! all_states_act_two_rdm_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs -! = +! all_states_act_two_rdm_alpha_alpha_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) @@ -20,11 +23,14 @@ BEGIN_PROVIDER [double precision, all_states_act_two_rdm_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none - double precision, allocatable :: state_weights(:) BEGIN_DOC -! all_states_act_two_rdm_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs -! = +! all_states_act_two_rdm_beta_beta_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC + double precision, allocatable :: state_weights(:) allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) integer :: ispin @@ -39,8 +45,11 @@ implicit none double precision, allocatable :: state_weights(:) BEGIN_DOC -! all_states_act_two_rdm_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs -! = +! all_states_act_two_rdm_alpha_beta_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) @@ -61,10 +70,15 @@ BEGIN_PROVIDER [double precision, all_states_act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none BEGIN_DOC -! all_states_act_two_rdm_spin_trace_mo(i,j,k,l) = 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: +! all_states_act_two_rdm_spin_trace_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM +! +! \sum_{\sigma, \sigma'} ! -! \sum_{i,j,k,l = 1, n_act_orb} all_states_act_two_rdm_spin_trace_mo(i,j,k,l) * < ii jj | kk ll > +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" +! +! The active part of the two-electron energy for the state istate can be computed as: +! +! \sum_{i,j,k,l = 1, n_act_orb} all_states_act_two_rdm_spin_trace_mo(i,j,k,l,istate) * < ii jj | kk ll > ! ! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l) END_DOC diff --git a/src/two_body_rdm/compute_all_states.irp.f b/src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f similarity index 100% rename from src/two_body_rdm/compute_all_states.irp.f rename to src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f diff --git a/src/two_body_rdm/compute.irp.f b/src/two_body_rdm/compute.irp.f deleted file mode 100644 index 112d2e36..00000000 --- a/src/two_body_rdm/compute.irp.f +++ /dev/null @@ -1,269 +0,0 @@ - - - subroutine diagonal_contrib_to_two_rdm_ab_dm(det_1,c_1,big_array,dim1,dim2,dim3,dim4) - use bitmasks - BEGIN_DOC -! routine that update the DIAGONAL PART of the alpha/beta two body rdm IN CHEMIST NOTATIONS - END_DOC - implicit none - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - integer(bit_kind), intent(in) :: det_1(N_int,2) - double precision, intent(in) :: c_1(N_states) - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2) - integer :: i,j,h1,h2,istate - double precision :: c_1_bis - call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) - do istate = 1, N_states - c_1_bis = c_1(istate) * c_1(istate) - do i = 1, n_occ_ab(1) - h1 = occ(i,1) - do j = 1, n_occ_ab(2) - h2 = occ(j,2) - big_array(h1,h1,h2,h2,istate) += c_1_bis - enddo - enddo - enddo - end - - - subroutine diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4) - use bitmasks - BEGIN_DOC -! routine that update the DIAGONAL PART of ALL THREE two body rdm IN CHEMIST NOTATIONS - END_DOC - implicit none - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states) - double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states) - integer(bit_kind), intent(in) :: det_1(N_int,2) - double precision, intent(in) :: c_1(N_states) - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2) - integer :: i,j,h1,h2,istate - double precision :: c_1_bis - BEGIN_DOC -! no factor 1/2 have to be taken into account as the permutations are already taken into account - END_DOC - call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) - do istate = 1, N_states - c_1_bis = c_1(istate) * c_1(istate) - do i = 1, n_occ_ab(1) - h1 = occ(i,1) - do j = 1, n_occ_ab(2) - h2 = occ(j,2) - big_array_ab(h1,h1,h2,h2,istate) += c_1_bis - enddo - do j = 1, n_occ_ab(1) - h2 = occ(j,1) - big_array_aa(h1,h1,h2,h2,istate) += 0.5d0 * c_1_bis - big_array_aa(h1,h2,h2,h1,istate) -= 0.5d0 * c_1_bis - enddo - enddo - do i = 1, n_occ_ab(2) - h1 = occ(i,2) - do j = 1, n_occ_ab(2) - h2 = occ(j,2) - big_array_bb(h1,h1,h2,h2,istate) += 0.5d0 * c_1_bis - big_array_bb(h1,h2,h2,h1,istate) -= 0.5d0 * c_1_bis - enddo - enddo - enddo - end - - - subroutine off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) - use bitmasks - BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS IN CHEMIST NOTATIONS - END_DOC - implicit none - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - double precision, intent(in) :: c_1(N_states),c_2(N_states) - integer :: i,j,h1,h2,p1,p2,istate - integer :: exc(0:2,2,2) - double precision :: phase - call get_double_excitation(det_1,det_2,exc,phase,N_int) - h1 = exc(1,1,1) - h2 = exc(1,1,2) - p1 = exc(1,2,1) - p2 = exc(1,2,2) - do istate = 1, N_states - big_array(h1,p1,h2,p2,istate) += c_1(istate) * phase * c_2(istate) -! big_array(p1,h1,p2,h2,istate) += c_1(istate) * phase * c_2(istate) - enddo - end - - subroutine off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) - use bitmasks - BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS - END_DOC - implicit none - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - double precision, intent(in) :: c_1(N_states),c_2(N_states) - 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 - 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 (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - p1 = exc(1,2,1) - do istate = 1, N_states - do i = 1, n_occ_ab(2) - h2 = occ(i,2) - big_array(h1,p1,h2,h2,istate) += 1.d0 * c_1(istate) * c_2(istate) * phase - enddo - enddo - else - ! Mono beta - h1 = exc(1,1,2) - p1 = exc(1,2,2) - do istate = 1, N_states - do i = 1, n_occ_ab(1) - h2 = occ(i,1) - big_array(h2,h2,h1,p1,istate) += 1.d0 * c_1(istate) * c_2(istate) * phase - enddo - enddo - endif - end - - subroutine off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) - BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS - END_DOC - use bitmasks - implicit none - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - double precision, intent(in) :: c_1(N_states),c_2(N_states) - 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 - 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 (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - p1 = exc(1,2,1) - do istate = 1, N_states - do i = 1, n_occ_ab(1) - h2 = occ(i,1) - big_array(h1,p1,h2,h2,istate) += 0.5d0 * c_1(istate) * c_2(istate) * phase - big_array(h1,h2,h2,p1,istate) -= 0.5d0 * c_1(istate) * c_2(istate) * phase - - big_array(h2,h2,h1,p1,istate) += 0.5d0 * c_1(istate) * c_2(istate) * phase - big_array(h2,p1,h1,h2,istate) -= 0.5d0 * c_1(istate) * c_2(istate) * phase - enddo - enddo - else - return - endif - end - - subroutine off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) - use bitmasks - BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS - END_DOC - implicit none - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - double precision, intent(in) :: c_1(N_states),c_2(N_states) - 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 - 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 (exc(0,1,1) == 1) then - return - else - ! Mono beta - h1 = exc(1,1,2) - p1 = exc(1,2,2) - do istate = 1, N_states - do i = 1, n_occ_ab(2) - h2 = occ(i,2) - big_array(h1,p1,h2,h2,istate) += 0.5d0 * c_1(istate) * c_2(istate) * phase - big_array(h1,h2,h2,p1,istate) -= 0.5d0 * c_1(istate) * c_2(istate) * phase - - big_array(h2,h2,h1,p1,istate) += 0.5d0 * c_1(istate) * c_2(istate) * phase - big_array(h2,p1,h1,h2,istate) -= 0.5d0 * c_1(istate) * c_2(istate) * phase - enddo - enddo - endif - end - - - subroutine off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) - use bitmasks - BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS IN CHEMIST NOTATIONS - END_DOC - implicit none - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) - double precision, intent(in) :: c_1(N_states),c_2(N_states) - integer :: i,j,h1,h2,p1,p2,istate - integer :: exc(0:2,2) - double precision :: phase - call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) - h1 =exc(1,1) - h2 =exc(2,1) - p1 =exc(1,2) - p2 =exc(2,2) -!print*,'h1,p1,h2,p2',h1,p1,h2,p2,c_1(istate) * phase * c_2(istate) - do istate = 1, N_states - big_array(h1,p1,h2,p2,istate) += 0.5d0 * c_1(istate) * phase * c_2(istate) - big_array(h1,p2,h2,p1,istate) -= 0.5d0 * c_1(istate) * phase * c_2(istate) - - big_array(h2,p2,h1,p1,istate) += 0.5d0 * c_1(istate) * phase * c_2(istate) - big_array(h2,p1,h1,p2,istate) -= 0.5d0 * c_1(istate) * phase * c_2(istate) - enddo - end - - subroutine off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4) - use bitmasks - BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for DOUBLE EXCITATIONS - END_DOC - implicit none - integer, intent(in) :: dim1,dim2,dim3,dim4 - double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states) - integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) - double precision, intent(in) :: c_1(N_states),c_2(N_states) - integer :: i,j,h1,h2,p1,p2,istate - integer :: exc(0:2,2) - double precision :: phase - call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) - h1 =exc(1,1) - h2 =exc(2,1) - p1 =exc(1,2) - p2 =exc(2,2) -!print*,'h1,p1,h2,p2',h1,p1,h2,p2,c_1(istate) * phase * c_2(istate) - do istate = 1, N_states - big_array(h1,p1,h2,p2,istate) += 0.5d0 * c_1(istate) * phase * c_2(istate) - big_array(h1,p2,h2,p1,istate) -= 0.5d0 * c_1(istate) * phase * c_2(istate) - - big_array(h2,p2,h1,p1,istate) += 0.5d0 * c_1(istate) * phase * c_2(istate) - big_array(h2,p1,h1,p2,istate) -= 0.5d0 * c_1(istate) * phase * c_2(istate) - enddo - end - diff --git a/src/two_body_rdm/compute_orb_range_omp.irp.f b/src/two_body_rdm/compute_orb_range_omp.irp.f deleted file mode 100644 index 0ba934d7..00000000 --- a/src/two_body_rdm/compute_orb_range_omp.irp.f +++ /dev/null @@ -1,807 +0,0 @@ - subroutine orb_range_diag_to_all_two_rdm_dm_buffer(det_1,c_1,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - use bitmasks - BEGIN_DOC - ! routine that update the DIAGONAL PART of the two body rdms in a specific range of orbitals for a given determinant det_1 - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - END_DOC - implicit none - integer, intent(in) :: ispin,sze_buff - integer, intent(in) :: list_orb_reverse(mo_num) - integer(bit_kind), intent(in) :: det_1(N_int,2) - integer(bit_kind), intent(in) :: orb_bitmask(N_int) - double precision, intent(in) :: c_1 - double precision, intent(out) :: values(sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2) - integer :: i,j,h1,h2 - integer(bit_kind) :: det_1_act(N_int,2) - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - do i = 1, N_int - det_1_act(i,1) = iand(det_1(i,1),orb_bitmask(i)) - det_1_act(i,2) = iand(det_1(i,2),orb_bitmask(i)) - enddo - - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - call bitstring_to_list_ab(det_1_act, occ, n_occ_ab, N_int) - logical :: is_integer_in_string - integer :: i1,i2 - if(alpha_beta)then - do i = 1, n_occ_ab(1) - i1 = occ(i,1) - do j = 1, n_occ_ab(2) - i2 = occ(j,2) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - values(nkeys) = c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - enddo - enddo - else if (alpha_alpha)then - do i = 1, n_occ_ab(1) - i1 = occ(i,1) - do j = 1, n_occ_ab(1) - i2 = occ(j,1) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = -0.5d0 * c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - else if (beta_beta)then - do i = 1, n_occ_ab(2) - i1 = occ(i,2) - do j = 1, n_occ_ab(2) - i2 = occ(j,2) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = -0.5d0 * c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - else if(spin_trace)then - ! 0.5 * (alpha beta + beta alpha) - do i = 1, n_occ_ab(1) - i1 = occ(i,1) - do j = 1, n_occ_ab(2) - i2 = occ(j,2) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - do i = 1, n_occ_ab(1) - i1 = occ(i,1) - do j = 1, n_occ_ab(1) - i2 = occ(j,1) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = -0.5d0 * c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - do i = 1, n_occ_ab(2) - i1 = occ(i,2) - do j = 1, n_occ_ab(2) - i2 = occ(j,2) - h1 = list_orb_reverse(i1) - h2 = list_orb_reverse(i2) - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = -0.5d0 * c_1 - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = h1 - enddo - enddo - endif - end - - - subroutine orb_range_off_diag_double_to_two_rdm_ab_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - use bitmasks - BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for -! -! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another -! -! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 -! -! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals -! -! ispin determines which spin-spin component of the two-rdm you will update -! -! ispin == 1 :: alpha/ alpha -! ispin == 2 :: beta / beta -! ispin == 3 :: alpha/ beta -! ispin == 4 :: spin traced <=> total two-rdm -! -! here, only ispin == 3 or 4 will do something - END_DOC - implicit none - integer, intent(in) :: ispin,sze_buff - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - integer, intent(in) :: list_orb_reverse(mo_num) - double precision, intent(in) :: c_1 - double precision, intent(out) :: values(sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - integer :: i,j,h1,h2,p1,p2 - integer :: exc(0:2,2,2) - double precision :: phase - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - call get_double_excitation(det_1,det_2,exc,phase,N_int) - h1 = exc(1,1,1) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - h2 = exc(1,1,2) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - p1 = exc(1,2,1) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - p2 = exc(1,2,2) - if(list_orb_reverse(p2).lt.0)return - p2 = list_orb_reverse(p2) - if(alpha_beta)then - nkeys += 1 - values(nkeys) = c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - else if(spin_trace)then - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = p1 - keys(2,nkeys) = p2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - endif - end - - subroutine orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - use bitmasks - BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 3 or 4 will do something - END_DOC - implicit none - integer, intent(in) :: ispin,sze_buff - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - integer, intent(in) :: list_orb_reverse(mo_num) - double precision, intent(in) :: c_1 - double precision, intent(out) :: values(sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2) - integer :: i,j,h1,h2,p1 - integer :: exc(0:2,2,2) - double precision :: phase - - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - - call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) - call get_single_excitation(det_1,det_2,exc,phase,N_int) - if(alpha_beta)then - if (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,1) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(2) - h2 = occ(i,2) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - nkeys += 1 - values(nkeys) = c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - enddo - else - ! Mono beta - h1 = exc(1,1,2) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,2) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(1) - h2 = occ(i,1) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - nkeys += 1 - values(nkeys) = c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - enddo - endif - else if(spin_trace)then - if (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,1) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(2) - h2 = occ(i,2) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - enddo - else - ! Mono beta - h1 = exc(1,1,2) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,2) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - !print*,'****************' - !print*,'****************' - !print*,'h1,p1',h1,p1 - do i = 1, n_occ_ab(1) - h2 = occ(i,1) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - ! print*,'h2 = ',h2 - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - enddo - endif - endif - end - - subroutine orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 1 or 4 will do something - END_DOC - use bitmasks - implicit none - integer, intent(in) :: ispin,sze_buff - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - integer, intent(in) :: list_orb_reverse(mo_num) - double precision, intent(in) :: c_1 - double precision, intent(out) :: values(sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2) - integer :: i,j,h1,h2,p1 - integer :: exc(0:2,2,2) - double precision :: phase - - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - - call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) - call get_single_excitation(det_1,det_2,exc,phase,N_int) - if(alpha_alpha.or.spin_trace)then - if (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,1) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(1) - h2 = occ(i,1) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - - nkeys += 1 - values(nkeys) = - 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = - 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - enddo - else - return - endif - endif - end - - subroutine orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - use bitmasks - BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 2 or 4 will do something - END_DOC - implicit none - integer, intent(in) :: ispin,sze_buff - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - integer, intent(in) :: list_orb_reverse(mo_num) - double precision, intent(in) :: c_1 - double precision, intent(out) :: values(sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2) - integer :: i,j,h1,h2,p1 - integer :: exc(0:2,2,2) - double precision :: phase - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - - - call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) - call get_single_excitation(det_1,det_2,exc,phase,N_int) - if(beta_beta.or.spin_trace)then - if (exc(0,1,1) == 1) then - return - else - ! Mono beta - h1 = exc(1,1,2) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - p1 = exc(1,2,2) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - do i = 1, n_occ_ab(2) - h2 = occ(i,2) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - - nkeys += 1 - values(nkeys) = - 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = - 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - enddo - endif - endif - end - - - subroutine orb_range_off_diag_double_to_two_rdm_aa_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - use bitmasks - BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 1 or 4 will do something - END_DOC - implicit none - integer, intent(in) :: ispin,sze_buff - integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) - integer, intent(in) :: list_orb_reverse(mo_num) - double precision, intent(in) :: c_1 - double precision, intent(out) :: values(sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - - - integer :: i,j,h1,h2,p1,p2 - integer :: exc(0:2,2) - double precision :: phase - - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) - h1 =exc(1,1) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - h2 =exc(2,1) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - p1 =exc(1,2) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - p2 =exc(2,2) - if(list_orb_reverse(p2).lt.0)return - p2 = list_orb_reverse(p2) - if(alpha_alpha.or.spin_trace)then - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - - nkeys += 1 - values(nkeys) = - 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = - 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - endif - end - - subroutine orb_range_off_diag_double_to_two_rdm_bb_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - use bitmasks - BEGIN_DOC - ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for - ! - ! a given couple of determinant det_1, det_2 being a BETA /BETA DOUBLE excitation with respect to one another - ! - ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 - ! - ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation - ! - ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals - ! - ! ispin determines which spin-spin component of the two-rdm you will update - ! - ! ispin == 1 :: alpha/ alpha - ! ispin == 2 :: beta / beta - ! ispin == 3 :: alpha/ beta - ! ispin == 4 :: spin traced <=> total two-rdm - ! - ! here, only ispin == 2 or 4 will do something - END_DOC - implicit none - - integer, intent(in) :: ispin,sze_buff - integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) - integer, intent(in) :: list_orb_reverse(mo_num) - double precision, intent(in) :: c_1 - double precision, intent(out) :: values(sze_buff) - integer , intent(out) :: keys(4,sze_buff) - integer , intent(inout):: nkeys - - integer :: i,j,h1,h2,p1,p2 - integer :: exc(0:2,2) - double precision :: phase - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - logical :: is_integer_in_string - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - endif - - call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) - h1 =exc(1,1) - if(list_orb_reverse(h1).lt.0)return - h1 = list_orb_reverse(h1) - h2 =exc(2,1) - if(list_orb_reverse(h2).lt.0)return - h2 = list_orb_reverse(h2) - p1 =exc(1,2) - if(list_orb_reverse(p1).lt.0)return - p1 = list_orb_reverse(p1) - p2 =exc(2,2) - if(list_orb_reverse(p2).lt.0)return - p2 = list_orb_reverse(p2) - if(beta_beta.or.spin_trace)then - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - - nkeys += 1 - values(nkeys) = - 0.5d0 * c_1 * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = - 0.5d0 * c_1 * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - endif - end - diff --git a/src/two_body_rdm/orb_range_omp.irp.f b/src/two_body_rdm/orb_range_omp.irp.f deleted file mode 100644 index baa26ced..00000000 --- a/src/two_body_rdm/orb_range_omp.irp.f +++ /dev/null @@ -1,85 +0,0 @@ - - BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] - implicit none - double precision, allocatable :: state_weights(:) - BEGIN_DOC -! state_av_act_two_rdm_openmp_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs -! = - END_DOC - allocate(state_weights(N_states)) - state_weights = state_average_weight - integer :: ispin - ! condition for alpha/beta spin - ispin = 1 - state_av_act_two_rdm_openmp_alpha_alpha_mo = 0.D0 - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) - - END_PROVIDER - - BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] - implicit none - double precision, allocatable :: state_weights(:) - BEGIN_DOC -! state_av_act_two_rdm_openmp_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs -! = - END_DOC - allocate(state_weights(N_states)) - state_weights = state_average_weight - integer :: ispin - ! condition for alpha/beta spin - ispin = 2 - state_av_act_two_rdm_openmp_beta_beta_mo = 0.d0 - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_beta_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) - - END_PROVIDER - - BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] - implicit none - double precision, allocatable :: state_weights(:) - BEGIN_DOC -! state_av_act_two_rdm_openmp_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs -! = - END_DOC - allocate(state_weights(N_states)) - state_weights = state_average_weight - integer :: ispin - ! condition for alpha/beta spin - print*,'' - print*,'' - print*,'' - print*,'providint state_av_act_two_rdm_openmp_alpha_beta_mo ' - ispin = 3 - print*,'ispin = ',ispin - state_av_act_two_rdm_openmp_alpha_beta_mo = 0.d0 - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_alpha_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) - - END_PROVIDER - - - BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] - implicit none - BEGIN_DOC -! state_av_act_two_rdm_openmp_spin_trace_mo(i,j,k,l) = state average physicist spin trace two-body rdm restricted to the ACTIVE indices -! The active part of the two-electron energy can be computed as: -! -! \sum_{i,j,k,l = 1, n_act_orb} state_av_act_two_rdm_openmp_spin_trace_mo(i,j,k,l) * < ii jj | kk ll > -! -! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l) - END_DOC - double precision, allocatable :: state_weights(:) - allocate(state_weights(N_states)) - state_weights = state_average_weight - integer :: ispin - ! condition for alpha/beta spin - ispin = 4 - state_av_act_two_rdm_openmp_spin_trace_mo = 0.d0 - integer :: i - double precision :: wall_0,wall_1 - call wall_time(wall_0) - print*,'providing the state average TWO-RDM ...' - call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) - - call wall_time(wall_1) - print*,'Time to provide the state average TWO-RDM',wall_1 - wall_0 - END_PROVIDER - diff --git a/src/two_body_rdm/orb_range_routines_omp.irp.f b/src/two_body_rdm/orb_range_routines_omp.irp.f deleted file mode 100644 index b6e59540..00000000 --- a/src/two_body_rdm/orb_range_routines_omp.irp.f +++ /dev/null @@ -1,568 +0,0 @@ -subroutine orb_range_two_rdm_state_av_openmp(big_array,dim1,norb,list_orb,state_weights,ispin,u_0,N_st,sze) - use bitmasks - implicit none - BEGIN_DOC - ! if ispin == 1 :: alpha/alpha 2rdm - ! == 2 :: beta /beta 2rdm - ! == 3 :: alpha/beta 2rdm - ! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba)) - ! - ! Assumes that the determinants are in psi_det - ! - ! istart, iend, ishift, istep are used in ZMQ parallelization. - END_DOC - integer, intent(in) :: N_st,sze - integer, intent(in) :: dim1,norb,list_orb(norb),ispin - double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) - double precision, intent(in) :: u_0(sze,N_st),state_weights(N_st) - - integer :: k - double precision, allocatable :: u_t(:,:) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t - allocate(u_t(N_st,N_det)) - do k=1,N_st - call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det) - enddo - call dtranspose( & - u_0, & - size(u_0, 1), & - u_t, & - size(u_t, 1), & - N_det, N_st) - - call orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1) - deallocate(u_t) - - do k=1,N_st - call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) - enddo - -end - -subroutine orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) - use bitmasks - implicit none - BEGIN_DOC - ! Computes two-rdm - ! - ! Default should be 1,N_det,0,1 - END_DOC - integer, intent(in) :: N_st,sze,istart,iend,ishift,istep - integer, intent(in) :: dim1,norb,list_orb(norb),ispin - double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) - double precision, intent(in) :: u_t(N_st,N_det),state_weights(N_st) - - integer :: k - - PROVIDE N_int - - select case (N_int) - case (1) - call orb_range_two_rdm_state_av_openmp_work_1(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) - case (2) - call orb_range_two_rdm_state_av_openmp_work_2(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) - case (3) - call orb_range_two_rdm_state_av_openmp_work_3(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) - case (4) - call orb_range_two_rdm_state_av_openmp_work_4(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) - case default - call orb_range_two_rdm_state_av_openmp_work_N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) - end select -end - - - - - BEGIN_TEMPLATE -subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) - use bitmasks - use omp_lib - implicit none - BEGIN_DOC - ! Computes the two rdm for the N_st vectors |u_t> - ! if ispin == 1 :: alpha/alpha 2rdm - ! == 2 :: beta /beta 2rdm - ! == 3 :: alpha/beta 2rdm - ! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba)) - ! The 2rdm will be computed only on the list of orbitals list_orb, which contains norb - ! In any cases, the state average weights will be used with an array state_weights - ! Default should be 1,N_det,0,1 for istart,iend,ishift,istep - END_DOC - integer, intent(in) :: N_st,sze,istart,iend,ishift,istep - double precision, intent(in) :: u_t(N_st,N_det),state_weights(N_st) - integer, intent(in) :: dim1,norb,list_orb(norb),ispin - double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) - - integer(omp_lock_kind) :: lock_2rdm - integer :: i,j,k,l - integer :: k_a, k_b, l_a, l_b - integer :: krow, kcol - integer :: lrow, lcol - integer(bit_kind) :: spindet($N_int) - integer(bit_kind) :: tmp_det($N_int,2) - integer(bit_kind) :: tmp_det2($N_int,2) - integer(bit_kind) :: tmp_det3($N_int,2) - integer(bit_kind), allocatable :: buffer(:,:) - integer :: n_doubles - integer, allocatable :: doubles(:) - integer, allocatable :: singles_a(:) - integer, allocatable :: singles_b(:) - integer, allocatable :: idx(:), idx0(:) - integer :: maxab, n_singles_a, n_singles_b, kcol_prev - double precision :: c_average - - logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace - integer(bit_kind) :: orb_bitmask($N_int) - integer :: list_orb_reverse(mo_num) - integer, allocatable :: keys(:,:) - double precision, allocatable :: values(:) - integer :: nkeys,sze_buff - alpha_alpha = .False. - beta_beta = .False. - alpha_beta = .False. - spin_trace = .False. - if( ispin == 1)then - alpha_alpha = .True. - else if(ispin == 2)then - beta_beta = .True. - else if(ispin == 3)then - alpha_beta = .True. - else if(ispin == 4)then - spin_trace = .True. - else - print*,'Wrong parameter for ispin in general_two_rdm_state_av_openmp_work' - print*,'ispin = ',ispin - stop - endif - - - PROVIDE N_int - - call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) - sze_buff = norb ** 3 + 6 * norb - list_orb_reverse = -1000 - do i = 1, norb - list_orb_reverse(list_orb(i)) = i - enddo - maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 - allocate(idx0(maxab)) - - do i=1,maxab - idx0(i) = i - enddo - call omp_init_lock(lock_2rdm) - - ! Prepare the array of all alpha single excitations - ! ------------------------------------------------- - - PROVIDE N_int nthreads_davidson elec_alpha_num - !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) & - !$OMP SHARED(psi_bilinear_matrix_rows, N_det,lock_2rdm,& - !$OMP psi_bilinear_matrix_columns, & - !$OMP psi_det_alpha_unique, psi_det_beta_unique,& - !$OMP n_det_alpha_unique, n_det_beta_unique, N_int,& - !$OMP psi_bilinear_matrix_transp_rows, & - !$OMP psi_bilinear_matrix_transp_columns, & - !$OMP psi_bilinear_matrix_transp_order, N_st, & - !$OMP psi_bilinear_matrix_order_transp_reverse, & - !$OMP psi_bilinear_matrix_columns_loc, & - !$OMP psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, & - !$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, state_weights, dim1, & - !$OMP ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask) & - !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, c_2, & - !$OMP lcol, lrow, l_a, l_b, & - !$OMP buffer, doubles, n_doubles, & - !$OMP tmp_det2, idx, l, kcol_prev, & - !$OMP singles_a, n_singles_a, singles_b, & - !$OMP n_singles_b, nkeys, keys, values, c_average) - - ! Alpha/Beta double excitations - ! ============================= - nkeys = 0 - allocate( keys(4,sze_buff), values(sze_buff)) - allocate( buffer($N_int,maxab), & - singles_a(maxab), & - singles_b(maxab), & - doubles(maxab), & - idx(maxab)) - - kcol_prev=-1 - - ASSERT (iend <= N_det) - ASSERT (istart > 0) - ASSERT (istep > 0) - - !$OMP DO SCHEDULE(dynamic,64) - do k_a=istart+ishift,iend,istep - - krow = psi_bilinear_matrix_rows(k_a) - ASSERT (krow <= N_det_alpha_unique) - - kcol = psi_bilinear_matrix_columns(k_a) - ASSERT (kcol <= N_det_beta_unique) - - tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) - tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - - if (kcol /= kcol_prev) then - call get_all_spin_singles_$N_int( & - psi_det_beta_unique, idx0, & - tmp_det(1,2), N_det_beta_unique, & - singles_b, n_singles_b) - endif - kcol_prev = kcol - - ! Loop over singly excited beta columns - ! ------------------------------------- - - do i=1,n_singles_b - lcol = singles_b(i) - - tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) - - l_a = psi_bilinear_matrix_columns_loc(lcol) - ASSERT (l_a <= N_det) - - do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a - lrow = psi_bilinear_matrix_rows(l_a) - ASSERT (lrow <= N_det_alpha_unique) - - buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) - - ASSERT (l_a <= N_det) - idx(j) = l_a - l_a = l_a+1 - enddo - j = j-1 - - call get_all_spin_singles_$N_int( & - buffer, idx, tmp_det(1,1), j, & - singles_a, n_singles_a ) - - ! Loop over alpha singles - ! ----------------------- - - if(alpha_beta.or.spin_trace)then - do k = 1,n_singles_a - l_a = singles_a(k) - ASSERT (l_a <= N_det) - - lrow = psi_bilinear_matrix_rows(l_a) - ASSERT (lrow <= N_det_alpha_unique) - - tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) - c_average = 0.d0 - do l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - c_average += c_1(l) * c_2(l) * state_weights(l) - enddo - if(alpha_beta)then - ! only ONE contribution - if (nkeys+1 .ge. size(values)) then - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - endif - else if (spin_trace)then - ! TWO contributions - if (nkeys+2 .ge. size(values)) then - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - endif - endif - call orb_range_off_diag_double_to_two_rdm_ab_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - - enddo - endif - - enddo - - enddo - !$OMP END DO - - !$OMP DO SCHEDULE(dynamic,64) - do k_a=istart+ishift,iend,istep - - - ! Single and double alpha exitations - ! =================================== - - - ! Initial determinant is at k_a in alpha-major representation - ! ----------------------------------------------------------------------- - - krow = psi_bilinear_matrix_rows(k_a) - ASSERT (krow <= N_det_alpha_unique) - - kcol = psi_bilinear_matrix_columns(k_a) - ASSERT (kcol <= N_det_beta_unique) - - tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) - tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - - ! Initial determinant is at k_b in beta-major representation - ! ---------------------------------------------------------------------- - - k_b = psi_bilinear_matrix_order_transp_reverse(k_a) - ASSERT (k_b <= N_det) - - spindet(1:$N_int) = tmp_det(1:$N_int,1) - - ! Loop inside the beta column to gather all the connected alphas - lcol = psi_bilinear_matrix_columns(k_a) - l_a = psi_bilinear_matrix_columns_loc(lcol) - do i=1,N_det_alpha_unique - if (l_a > N_det) exit - lcol = psi_bilinear_matrix_columns(l_a) - if (lcol /= kcol) exit - lrow = psi_bilinear_matrix_rows(l_a) - ASSERT (lrow <= N_det_alpha_unique) - - buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) - idx(i) = l_a - l_a = l_a+1 - enddo - i = i-1 - - call get_all_spin_singles_and_doubles_$N_int( & - buffer, idx, spindet, i, & - singles_a, doubles, n_singles_a, n_doubles ) - - ! Compute Hij for all alpha singles - ! ---------------------------------- - - tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - do i=1,n_singles_a - l_a = singles_a(i) - ASSERT (l_a <= N_det) - - lrow = psi_bilinear_matrix_rows(l_a) - ASSERT (lrow <= N_det_alpha_unique) - - tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) - c_average = 0.d0 - do l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - c_average += c_1(l) * c_2(l) * state_weights(l) - enddo - if(alpha_beta.or.spin_trace.or.alpha_alpha)then - ! increment the alpha/beta part for single excitations - if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - ! increment the alpha/alpha part for single excitations - if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - endif - - enddo - - - ! Compute Hij for all alpha doubles - ! ---------------------------------- - - if(alpha_alpha.or.spin_trace)then - do i=1,n_doubles - l_a = doubles(i) - ASSERT (l_a <= N_det) - - lrow = psi_bilinear_matrix_rows(l_a) - ASSERT (lrow <= N_det_alpha_unique) - - c_average = 0.d0 - do l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - c_average += c_1(l) * c_2(l) * state_weights(l) - enddo - if (nkeys+4 .ge. sze_buff) then - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_double_to_two_rdm_aa_dm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - enddo - endif - - - ! Single and double beta excitations - ! ================================== - - - ! Initial determinant is at k_a in alpha-major representation - ! ----------------------------------------------------------------------- - - krow = psi_bilinear_matrix_rows(k_a) - kcol = psi_bilinear_matrix_columns(k_a) - - tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) - tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - - spindet(1:$N_int) = tmp_det(1:$N_int,2) - - ! Initial determinant is at k_b in beta-major representation - ! ----------------------------------------------------------------------- - - k_b = psi_bilinear_matrix_order_transp_reverse(k_a) - ASSERT (k_b <= N_det) - - ! Loop inside the alpha row to gather all the connected betas - lrow = psi_bilinear_matrix_transp_rows(k_b) - l_b = psi_bilinear_matrix_transp_rows_loc(lrow) - do i=1,N_det_beta_unique - if (l_b > N_det) exit - lrow = psi_bilinear_matrix_transp_rows(l_b) - if (lrow /= krow) exit - lcol = psi_bilinear_matrix_transp_columns(l_b) - ASSERT (lcol <= N_det_beta_unique) - - buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol) - idx(i) = l_b - l_b = l_b+1 - enddo - i = i-1 - - call get_all_spin_singles_and_doubles_$N_int( & - buffer, idx, spindet, i, & - singles_b, doubles, n_singles_b, n_doubles ) - - ! Compute Hij for all beta singles - ! ---------------------------------- - - tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) - do i=1,n_singles_b - l_b = singles_b(i) - ASSERT (l_b <= N_det) - - lcol = psi_bilinear_matrix_transp_columns(l_b) - ASSERT (lcol <= N_det_beta_unique) - - tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) - l_a = psi_bilinear_matrix_transp_order(l_b) - c_average = 0.d0 - do l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - c_average += c_1(l) * c_2(l) * state_weights(l) - enddo - if(alpha_beta.or.spin_trace.or.beta_beta)then - ! increment the alpha/beta part for single excitations - if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - ! increment the beta /beta part for single excitations - if (nkeys+4 * elec_alpha_num .ge. sze_buff) then - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - endif - enddo - - ! Compute Hij for all beta doubles - ! ---------------------------------- - - if(beta_beta.or.spin_trace)then - do i=1,n_doubles - l_b = doubles(i) - ASSERT (l_b <= N_det) - - lcol = psi_bilinear_matrix_transp_columns(l_b) - ASSERT (lcol <= N_det_beta_unique) - - l_a = psi_bilinear_matrix_transp_order(l_b) - c_average = 0.d0 - do l= 1, N_states - c_1(l) = u_t(l,l_a) - c_2(l) = u_t(l,k_a) - c_average += c_1(l) * c_2(l) * state_weights(l) - enddo - if (nkeys+4 .ge. sze_buff) then - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_double_to_two_rdm_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - ASSERT (l_a <= N_det) - - enddo - endif - - - ! Diagonal contribution - ! ===================== - - - ! Initial determinant is at k_a in alpha-major representation - ! ----------------------------------------------------------------------- - - krow = psi_bilinear_matrix_rows(k_a) - ASSERT (krow <= N_det_alpha_unique) - - kcol = psi_bilinear_matrix_columns(k_a) - ASSERT (kcol <= N_det_beta_unique) - - tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) - tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - - double precision, external :: diag_wee_mat_elem, diag_S_mat_elem - - double precision :: c_1(N_states),c_2(N_states) - c_average = 0.d0 - do l = 1, N_states - c_1(l) = u_t(l,k_a) - c_average += c_1(l) * c_1(l) * state_weights(l) - enddo - - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - call orb_range_diag_to_all_two_rdm_dm_buffer(tmp_det,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - nkeys = 0 - - end do - !$OMP END DO - deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values) - !$OMP END PARALLEL - -end - - SUBST [ N_int ] - - 1;; - 2;; - 3;; - 4;; - N_int;; - - END_TEMPLATE - - -subroutine update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) - use omp_lib - implicit none - integer, intent(in) :: nkeys,dim1 - integer, intent(in) :: keys(4,nkeys) - double precision, intent(in) :: values(nkeys) - double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) - - integer(omp_lock_kind),intent(inout):: lock_2rdm - integer :: i,h1,h2,p1,p2 - call omp_set_lock(lock_2rdm) - do i = 1, nkeys - h1 = keys(1,i) - h2 = keys(2,i) - p1 = keys(3,i) - p2 = keys(4,i) - big_array(h1,h2,p1,p2) += values(i) - enddo - call omp_unset_lock(lock_2rdm) - -end - diff --git a/src/two_body_rdm/orb_range_routines.irp.f b/src/two_body_rdm/state_av_act_2_rdm_dav_routines.irp.f similarity index 100% rename from src/two_body_rdm/orb_range_routines.irp.f rename to src/two_body_rdm/state_av_act_2_rdm_dav_routines.irp.f diff --git a/src/two_body_rdm/orb_range.irp.f b/src/two_body_rdm/state_av_act_2_rdm_prov.irp.f similarity index 68% rename from src/two_body_rdm/orb_range.irp.f rename to src/two_body_rdm/state_av_act_2_rdm_prov.irp.f index 2bcd04dc..420a0264 100644 --- a/src/two_body_rdm/orb_range.irp.f +++ b/src/two_body_rdm/state_av_act_2_rdm_prov.irp.f @@ -5,8 +5,11 @@ implicit none double precision, allocatable :: state_weights(:) BEGIN_DOC -! state_av_act_two_rdm_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs -! = +! state_av_act_two_rdm_alpha_alpha_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC allocate(state_weights(N_states)) state_weights = state_average_weight @@ -22,8 +25,11 @@ implicit none double precision, allocatable :: state_weights(:) BEGIN_DOC -! state_av_act_two_rdm_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs -! = +! state_av_act_two_rdm_beta_beta_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of beta electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC allocate(state_weights(N_states)) state_weights = state_average_weight @@ -39,8 +45,11 @@ implicit none double precision, allocatable :: state_weights(:) BEGIN_DOC -! state_av_act_two_rdm_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs -! = +! state_av_act_two_rdm_alpha_beta_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha/beta electrons +! +! +! +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" END_DOC allocate(state_weights(N_states)) state_weights = state_average_weight @@ -61,13 +70,12 @@ BEGIN_PROVIDER [double precision, state_av_act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] implicit none BEGIN_DOC -! state_av_act_two_rdm_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: +! state_av_act_two_rdm_spin_trace_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM +! +! \sum_{\sigma, \sigma'} ! -! \sum_{i,j,k,l = 1, n_act_orb} state_av_act_two_rdm_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 +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" + END_DOC double precision, allocatable :: state_weights(:) allocate(state_weights(N_states)) state_weights = state_average_weight diff --git a/src/two_body_rdm/compute_orb_range.irp.f b/src/two_body_rdm/state_av_act_2_rdm_update_routines.irp.f similarity index 100% rename from src/two_body_rdm/compute_orb_range.irp.f rename to src/two_body_rdm/state_av_act_2_rdm_update_routines.irp.f diff --git a/src/two_body_rdm/two_rdm.irp.f b/src/two_body_rdm/two_rdm.irp.f deleted file mode 100644 index c162f365..00000000 --- a/src/two_body_rdm/two_rdm.irp.f +++ /dev/null @@ -1,62 +0,0 @@ - BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] -&BEGIN_PROVIDER [double precision, two_rdm_alpha_alpha_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] -&BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] - implicit none - BEGIN_DOC - ! two_rdm_alpha_beta(i,j,k,l) = - ! 1 1 2 2 = chemist notations - ! note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry - ! - END_DOC - integer :: dim1,dim2,dim3,dim4 - double precision :: cpu_0,cpu_1 - dim1 = mo_num - dim2 = mo_num - dim3 = mo_num - dim4 = mo_num - two_rdm_alpha_beta_mo = 0.d0 - two_rdm_alpha_alpha_mo= 0.d0 - two_rdm_beta_beta_mo = 0.d0 - print*,'providing two_rdm_alpha_beta ...' - call wall_time(cpu_0) - 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) - -END_PROVIDER - - - BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)] -&BEGIN_PROVIDER [double precision, two_rdm_alpha_alpha_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)] -&BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)] - implicit none - BEGIN_DOC - ! two_rdm_alpha_beta_mo_physicist,(i,j,k,l) = - ! 1 2 1 2 = physicist notations - ! note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry - ! - END_DOC - integer :: i,j,k,l,istate - double precision :: cpu_0,cpu_1 - two_rdm_alpha_beta_mo_physicist = 0.d0 - print*,'providing two_rdm_alpha_beta_mo_physicist ...' - call wall_time(cpu_0) - do istate = 1, N_states - do i = 1, mo_num - do j = 1, mo_num - do k = 1, mo_num - do l = 1, mo_num - ! 1 2 1 2 1 1 2 2 - two_rdm_alpha_beta_mo_physicist(l,k,i,j,istate) = two_rdm_alpha_beta_mo(i,l,j,k,istate) - two_rdm_alpha_alpha_mo_physicist(l,k,i,j,istate) = two_rdm_alpha_alpha_mo(i,l,j,k,istate) - two_rdm_beta_beta_mo_physicist(l,k,i,j,istate) = two_rdm_beta_beta_mo(i,l,j,k,istate) - enddo - enddo - enddo - enddo - enddo - call wall_time(cpu_1) - print*,'two_rdm_alpha_beta_mo_physicist provided in',dabs(cpu_1-cpu_0) - -END_PROVIDER -