From e9724fa8c74662b4d39335e6c4b185f7d04bbaf4 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 28 Jun 2019 15:17:04 +0200 Subject: [PATCH 1/6] beginning to work on general routine for 2rdm --- src/two_body_rdm/general_2rdm_routines.irp.f | 488 ++++++++++++++++++ src/two_body_rdm/orb_range_2_rdm.irp.f | 61 +++ .../routines_compute_2rdm_orb_range.irp.f | 430 +++++++++++++++ 3 files changed, 979 insertions(+) create mode 100644 src/two_body_rdm/general_2rdm_routines.irp.f create mode 100644 src/two_body_rdm/orb_range_2_rdm.irp.f create mode 100644 src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f diff --git a/src/two_body_rdm/general_2rdm_routines.irp.f b/src/two_body_rdm/general_2rdm_routines.irp.f new file mode 100644 index 00000000..a9fcd61a --- /dev/null +++ b/src/two_body_rdm/general_2rdm_routines.irp.f @@ -0,0 +1,488 @@ +subroutine orb_range_two_rdm_dm_nstates_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_dm_nstates_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_dm_nstates_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) + + + PROVIDE N_int + + select case (N_int) + case (1) + call orb_range_two_rdm_dm_nstates_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_dm_nstates_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_dm_nstates_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_dm_nstates_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_dm_nstates_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_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + use bitmasks + implicit none + BEGIN_DOC + ! Computes the two rdm for the N_st vectors |u_t> + ! if ispin == 1 :: alpha/alpha 2rdm + ! == 2 :: beta /beta 2rdm + ! == 3 :: alpha/beta 2rdm + ! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba)) + ! The 2rdm will be computed only on the list of orbitals list_orb, which contains norb + ! In any cases, the state average weights will be used with an array state_weights + ! Default should be 1,N_det,0,1 for istart,iend,ishift,istep + END_DOC + integer, intent(in) :: N_st,sze,istart,iend,ishift,istep + double precision, intent(in) :: u_t(N_st,N_det),state_weights(N_st) + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + + integer :: i,j,k,l + integer :: k_a, k_b, l_a, l_b, m_a, m_b + integer :: istate + integer :: krow, kcol, krow_b, kcol_b + integer :: lrow, lcol + integer :: mrow, mcol + integer(bit_kind) :: spindet($N_int) + integer(bit_kind) :: tmp_det($N_int,2) + integer(bit_kind) :: tmp_det2($N_int,2) + integer(bit_kind) :: tmp_det3($N_int,2) + integer(bit_kind), allocatable :: buffer(:,:) + integer :: n_doubles + integer, allocatable :: doubles(:) + integer, allocatable :: singles_a(:) + integer, allocatable :: singles_b(:) + integer, allocatable :: idx(:), idx0(:) + integer :: maxab, n_singles_a, n_singles_b, kcol_prev + integer*8 :: k8 + double precision :: c_average + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + 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_dm_nstates_openmp_work' + print*,'ispin = ',ispin + stop + endif + + + 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 + ! ----------------------- + + 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 + call orb_range_off_diagonal_double_to_two_rdm_ab_dm(tmp_det,tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + enddo + endif + + enddo + + enddo + ! !$OMP END DO + + ! !$OMP DO SCHEDULE(dynamic,64) + do k_a=istart+ishift,iend,istep + + + ! Single and double alpha exitations + ! =================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + ! Initial determinant is at k_b in beta-major representation + ! ---------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + spindet(1:$N_int) = tmp_det(1:$N_int,1) + + ! Loop inside the beta column to gather all the connected alphas + lcol = psi_bilinear_matrix_columns(k_a) + l_a = psi_bilinear_matrix_columns_loc(lcol) + do i=1,N_det_alpha_unique + if (l_a > N_det) exit + lcol = psi_bilinear_matrix_columns(l_a) + if (lcol /= kcol) exit + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) + idx(i) = l_a + l_a = l_a+1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles_$N_int( & + buffer, idx, spindet, i, & + singles_a, doubles, n_singles_a, n_doubles ) + + ! Compute Hij for all alpha singles + ! ---------------------------------- + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + do i=1,n_singles_a + l_a = singles_a(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + if(alpha_beta.or.spin_trace.or.alpha_alpha)then + ! increment the alpha/beta part for single excitations + call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + ! increment the alpha/alpha part for single excitations + call orb_range_off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + endif + + enddo + + + ! Compute Hij for all alpha doubles + ! ---------------------------------- + + if(alpha_alpha.or.spin_trace)then + do i=1,n_doubles + l_a = doubles(i) + ASSERT (l_a <= N_det) + + lrow = psi_bilinear_matrix_rows(l_a) + ASSERT (lrow <= N_det_alpha_unique) + + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + call orb_range_off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,big_array,dim1,norb,list_orb,ispin) + enddo + endif + + + ! Single and double beta excitations + ! ================================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + kcol = psi_bilinear_matrix_columns(k_a) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + spindet(1:$N_int) = tmp_det(1:$N_int,2) + + ! Initial determinant is at k_b in beta-major representation + ! ----------------------------------------------------------------------- + + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) + ASSERT (k_b <= N_det) + + ! Loop inside the alpha row to gather all the connected betas + lrow = psi_bilinear_matrix_transp_rows(k_b) + l_b = psi_bilinear_matrix_transp_rows_loc(lrow) + do i=1,N_det_beta_unique + if (l_b > N_det) exit + lrow = psi_bilinear_matrix_transp_rows(l_b) + if (lrow /= krow) exit + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol) + idx(i) = l_b + l_b = l_b+1 + enddo + i = i-1 + + call get_all_spin_singles_and_doubles_$N_int( & + buffer, idx, spindet, i, & + singles_b, doubles, n_singles_b, n_doubles ) + + ! Compute Hij for all beta singles + ! ---------------------------------- + + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + do i=1,n_singles_b + l_b = singles_b(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) + l_a = psi_bilinear_matrix_transp_order(l_b) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + if(alpha_beta.or.spin_trace.or.beta_beta)then + ! increment the alpha/beta part for single excitations + call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + ! increment the beta /beta part for single excitations + call orb_range_off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + endif + enddo + + ! Compute Hij for all beta doubles + ! ---------------------------------- + + if(beta_beta.or.spin_trace)then + do i=1,n_doubles + l_b = doubles(i) + ASSERT (l_b <= N_det) + + lcol = psi_bilinear_matrix_transp_columns(l_b) + ASSERT (lcol <= N_det_beta_unique) + + l_a = psi_bilinear_matrix_transp_order(l_b) + c_average = 0.d0 + do l= 1, N_states + c_1(l) = u_t(l,l_a) + c_2(l) = u_t(l,k_a) + c_average += c_1(l) * c_2(l) * state_weights(l) + enddo + call orb_range_off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_average,big_array,dim1,norb,list_orb,ispin) + ASSERT (l_a <= N_det) + + enddo + endif + + + ! Diagonal contribution + ! ===================== + + + ! Initial determinant is at k_a in alpha-major representation + ! ----------------------------------------------------------------------- + + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) + tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) + + double precision, external :: diag_wee_mat_elem, diag_S_mat_elem + + double precision :: c_1(N_states),c_2(N_states) + c_average = 0.d0 + do l = 1, N_states + c_1(l) = u_t(l,k_a) + c_average += c_1(l) * c_1(l) * state_weights(l) + enddo + + call orb_range_diagonal_contrib_to_all_two_rdm_dm(tmp_det,c_average,big_array,dim1,norb,list_orb,ispin) + + 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/orb_range_2_rdm.irp.f b/src/two_body_rdm/orb_range_2_rdm.irp.f new file mode 100644 index 00000000..621f6c4b --- /dev/null +++ b/src/two_body_rdm/orb_range_2_rdm.irp.f @@ -0,0 +1,61 @@ + + + + BEGIN_PROVIDER [double precision, act_two_rdm_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + allocate(state_weights(N_states)) + state_weights = 1.d0/dble(N_states) + integer :: ispin + ! condition for alpha/beta spin + ispin = 1 + act_two_rdm_alpha_alpha_mo = 0.D0 + call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + + BEGIN_PROVIDER [double precision, act_two_rdm_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + allocate(state_weights(N_states)) + state_weights = 1.d0/dble(N_states) + integer :: ispin + ! condition for alpha/beta spin + ispin = 2 + act_two_rdm_beta_beta_mo = 0.d0 + call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_beta_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + + BEGIN_PROVIDER [double precision, act_two_rdm_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + allocate(state_weights(N_states)) + state_weights = 1.d0/dble(N_states) + integer :: ispin + ! condition for alpha/beta spin + print*,'' + print*,'' + print*,'' + print*,'providint act_two_rdm_alpha_beta_mo ' + ispin = 3 + print*,'ispin = ',ispin + act_two_rdm_alpha_beta_mo = 0.d0 + call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_alpha_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + double precision, allocatable :: state_weights(:) + allocate(state_weights(N_states)) + state_weights = 1.d0/dble(N_states) + integer :: ispin + ! condition for alpha/beta spin + ispin = 4 + act_two_rdm_spin_trace_mo = 0.d0 + call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + + END_PROVIDER + diff --git a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f new file mode 100644 index 00000000..d115f1bd --- /dev/null +++ b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f @@ -0,0 +1,430 @@ + + subroutine orb_range_diagonal_contrib_to_two_rdm_ab_dm(det_1,c_1,big_array,dim1,norb,list_orb) + use bitmasks + BEGIN_DOC +! routine that update the DIAGONAL PART of the alpha/beta two body rdm in a specific range of orbitals +! c_1 is supposed to be a scalar quantity, such as state averaged coef + END_DOC + implicit none + integer, intent(in) :: dim1,norb,list_orb(norb) + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + integer(bit_kind), intent(in) :: det_1(N_int,2) + double precision, intent(in) :: c_1 + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,istate + call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) + 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) += c_1 + enddo + enddo + end + + + subroutine orb_range_diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array,dim1,norb,list_orb,ispin) + use bitmasks + BEGIN_DOC +! routine that update the DIAGONAL PART of ALL THREE two body rdm + END_DOC + implicit none + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + integer(bit_kind), intent(in) :: det_1(N_int,2) + double precision, intent(in) :: c_1 + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,istate + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + 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 + 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) + if(alpha_beta)then + 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) += c_1 + enddo + enddo + else if (alpha_alpha)then + do i = 1, n_occ_ab(1) + h1 = occ(i,1) + do j = 1, n_occ_ab(1) + h2 = occ(j,1) + big_array(h1,h1,h2,h2) += 0.5d0 * c_1 + big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 + enddo + enddo + else if (beta_beta)then + do i = 1, n_occ_ab(2) + h1 = occ(i,2) + do j = 1, n_occ_ab(2) + h2 = occ(j,2) + big_array(h1,h1,h2,h2) += 0.5d0 * c_1 + big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 + enddo + enddo + else if(spin_trace)then + ! 0.5 * (alpha beta + beta alpha) + 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) += 0.5d0 * (c_1 ) + big_array(h2,h2,h1,h1) += 0.5d0 * (c_1 ) + enddo + enddo + ! alpha alpha + do i = 1, n_occ_ab(1) + h1 = occ(i,1) + do j = 1, n_occ_ab(1) + h2 = occ(j,1) + big_array(h1,h1,h2,h2) += 0.5d0 * c_1 + big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 + enddo + enddo + ! beta beta + do i = 1, n_occ_ab(2) + h1 = occ(i,2) + do j = 1, n_occ_ab(2) + h2 = occ(j,2) + big_array(h1,h1,h2,h2) += 0.5d0 * c_1 + big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 + enddo + enddo + endif + end + + + subroutine orb_range_off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + use bitmasks + BEGIN_DOC +! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS + END_DOC + implicit none + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + double precision, intent(in) :: c_1 + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + 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) + h2 = exc(1,1,2) + p1 = exc(1,2,1) + p2 = exc(1,2,2) + if(alpha_beta)then + big_array(h1,p1,h2,p2) += c_1 * phase + else if(spin_trace)then + big_array(h1,p1,h2,p2) += 0.5d0 * c_1 * phase + big_array(h2,p2,h1,p1) += 0.5d0 * c_1 * phase + endif + end + + subroutine orb_range_off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + use bitmasks + BEGIN_DOC +! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS + END_DOC + implicit none + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + double precision, intent(in) :: c_1 + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,istate,p1 + integer :: exc(0:2,2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + 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) + p1 = exc(1,2,1) + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + big_array(h1,p1,h2,h2) += c_1 * phase + enddo + else + ! Mono beta + h1 = exc(1,1,2) + p1 = exc(1,2,2) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + big_array(h2,h2,h1,p1) += c_1 * phase + enddo + endif + else if(spin_trace)then + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + p1 = exc(1,2,1) + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + big_array(h1,p1,h2,h2) += 0.5d0 * c_1 * phase + big_array(h2,h2,h1,p1) += 0.5d0 * c_1 * phase + enddo + else + ! Mono beta + h1 = exc(1,1,2) + p1 = exc(1,2,2) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + big_array(h1,p1,h2,h2) += 0.5d0 * c_1 * phase + big_array(h2,h2,h1,p1) += 0.5d0 * c_1 * phase + enddo + endif + endif + end + + subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + BEGIN_DOC +! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS + END_DOC + use bitmasks + implicit none + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + double precision, intent(in) :: c_1 + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,istate,p1 + integer :: exc(0:2,2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + 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) + p1 = exc(1,2,1) + do i = 1, n_occ_ab(1) + h2 = occ(i,1) + big_array(h1,p1,h2,h2) += 0.5d0 * c_1 * phase + big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase + + big_array(h2,h2,h1,p1) += 0.5d0 * c_1 * phase + big_array(h2,p1,h1,h2) -= 0.5d0 * c_1 * phase + enddo + else + return + endif + endif + end + + subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + use bitmasks + BEGIN_DOC +! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for SINGLE EXCITATIONS + END_DOC + implicit none + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + double precision, intent(in) :: c_1 + + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,istate,p1 + integer :: exc(0:2,2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + 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) + 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) += 0.5d0 * c_1 * phase + big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase + + big_array(h2,h2,h1,p1) += 0.5d0 * c_1 * phase + big_array(h2,p1,h1,h2) -= 0.5d0 * c_1 * phase + enddo + enddo + endif + endif + end + + + subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + use bitmasks + BEGIN_DOC +! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS + END_DOC + implicit none + integer, intent(in) :: dim1,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + double precision, intent(in) :: c_1 + + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + 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) + h2 =exc(2,1) + p1 =exc(1,2) + p2 =exc(2,2) + if(alpha_alpha.or.spin_trace)then + do istate = 1, N_states + big_array(h1,p1,h2,p2) += 0.5d0 * c_1 * phase + big_array(h1,p2,h2,p1) -= 0.5d0 * c_1 * phase + + big_array(h2,p2,h1,p1) += 0.5d0 * c_1 * phase + big_array(h2,p1,h1,p2) -= 0.5d0 * c_1 * phase + enddo + endif + end + + subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + 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,norb,list_orb(norb),ispin + double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) + integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + double precision, intent(in) :: c_1 + + integer :: i,j,h1,h2,p1,p2,istate + integer :: exc(0:2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + 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) + h2 =exc(2,1) + p1 =exc(1,2) + p2 =exc(2,2) + if(beta_beta.or.spin_trace)then + big_array(h1,p1,h2,p2) += 0.5d0 * c_1* phase + big_array(h1,p2,h2,p1) -= 0.5d0 * c_1* phase + + big_array(h2,p2,h1,p1) += 0.5d0 * c_1* phase + big_array(h2,p1,h1,p2) -= 0.5d0 * c_1* phase + endif + end + From c90c49b56c100b367067d90edbf226dd784a8cc8 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 28 Jun 2019 15:55:32 +0200 Subject: [PATCH 2/6] beginning to do it directly in physicist --- .../routines_compute_2rdm_orb_range.irp.f | 64 +++++++++---------- 1 file changed, 31 insertions(+), 33 deletions(-) diff --git a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f index d115f1bd..d918932a 100644 --- a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f +++ b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f @@ -18,7 +18,7 @@ h1 = occ(i,1) do j = 1, n_occ_ab(2) h2 = occ(j,2) - big_array(h1,h1,h2,h2) += c_1 + big_array(h1,h2,h1,h2) += c_1 enddo enddo end @@ -61,7 +61,7 @@ h1 = occ(i,1) do j = 1, n_occ_ab(2) h2 = occ(j,2) - big_array(h1,h1,h2,h2) += c_1 + big_array(h1,h2,h1,h2) += c_1 enddo enddo else if (alpha_alpha)then @@ -69,7 +69,7 @@ h1 = occ(i,1) do j = 1, n_occ_ab(1) h2 = occ(j,1) - big_array(h1,h1,h2,h2) += 0.5d0 * c_1 + big_array(h1,h2,h1,h2) += 0.5d0 * c_1 big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 enddo enddo @@ -78,7 +78,7 @@ h1 = occ(i,2) do j = 1, n_occ_ab(2) h2 = occ(j,2) - big_array(h1,h1,h2,h2) += 0.5d0 * c_1 + big_array(h1,h2,h1,h2) += 0.5d0 * c_1 big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 enddo enddo @@ -88,25 +88,23 @@ h1 = occ(i,1) do j = 1, n_occ_ab(2) h2 = occ(j,2) - big_array(h1,h1,h2,h2) += 0.5d0 * (c_1 ) - big_array(h2,h2,h1,h1) += 0.5d0 * (c_1 ) + big_array(h1,h2,h1,h2) += 0.5d0 * (c_1 ) + big_array(h2,h1,h2,h1) += 0.5d0 * (c_1 ) enddo enddo - ! alpha alpha do i = 1, n_occ_ab(1) h1 = occ(i,1) do j = 1, n_occ_ab(1) h2 = occ(j,1) - big_array(h1,h1,h2,h2) += 0.5d0 * c_1 + big_array(h1,h2,h1,h2) += 0.5d0 * c_1 big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 enddo enddo - ! beta beta do i = 1, n_occ_ab(2) h1 = occ(i,2) do j = 1, n_occ_ab(2) h2 = occ(j,2) - big_array(h1,h1,h2,h2) += 0.5d0 * c_1 + big_array(h1,h2,h1,h2) += 0.5d0 * c_1 big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 enddo enddo @@ -147,10 +145,10 @@ p1 = exc(1,2,1) p2 = exc(1,2,2) if(alpha_beta)then - big_array(h1,p1,h2,p2) += c_1 * phase + big_array(h1,h2,p1,p2) += c_1 * phase else if(spin_trace)then - big_array(h1,p1,h2,p2) += 0.5d0 * c_1 * phase - big_array(h2,p2,h1,p1) += 0.5d0 * c_1 * phase + big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase + big_array(p1,p2,h1,h2) += 0.5d0 * c_1 * phase endif end @@ -195,7 +193,7 @@ p1 = exc(1,2,1) do i = 1, n_occ_ab(2) h2 = occ(i,2) - big_array(h1,p1,h2,h2) += c_1 * phase + big_array(h1,h2,p1,h2) += c_1 * phase enddo else ! Mono beta @@ -203,7 +201,7 @@ p1 = exc(1,2,2) do i = 1, n_occ_ab(1) h2 = occ(i,1) - big_array(h2,h2,h1,p1) += c_1 * phase + big_array(h2,h1,h2,p1) += c_1 * phase enddo endif else if(spin_trace)then @@ -213,8 +211,8 @@ p1 = exc(1,2,1) do i = 1, n_occ_ab(2) h2 = occ(i,2) - big_array(h1,p1,h2,h2) += 0.5d0 * c_1 * phase - big_array(h2,h2,h1,p1) += 0.5d0 * c_1 * phase + big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase + big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase enddo else ! Mono beta @@ -222,8 +220,8 @@ p1 = exc(1,2,2) do i = 1, n_occ_ab(1) h2 = occ(i,1) - big_array(h1,p1,h2,h2) += 0.5d0 * c_1 * phase - big_array(h2,h2,h1,p1) += 0.5d0 * c_1 * phase + big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase + big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase enddo endif endif @@ -270,11 +268,11 @@ p1 = exc(1,2,1) do i = 1, n_occ_ab(1) h2 = occ(i,1) - big_array(h1,p1,h2,h2) += 0.5d0 * c_1 * phase + big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase - big_array(h2,h2,h1,p1) += 0.5d0 * c_1 * phase - big_array(h2,p1,h1,h2) -= 0.5d0 * c_1 * phase + big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase + big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase enddo else return @@ -327,11 +325,11 @@ do istate = 1, N_states do i = 1, n_occ_ab(2) h2 = occ(i,2) - big_array(h1,p1,h2,h2) += 0.5d0 * c_1 * phase + big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase - big_array(h2,h2,h1,p1) += 0.5d0 * c_1 * phase - big_array(h2,p1,h1,h2) -= 0.5d0 * c_1 * phase + big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase + big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase enddo enddo endif @@ -375,11 +373,11 @@ p2 =exc(2,2) if(alpha_alpha.or.spin_trace)then do istate = 1, N_states - big_array(h1,p1,h2,p2) += 0.5d0 * c_1 * phase - big_array(h1,p2,h2,p1) -= 0.5d0 * c_1 * phase + big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase + big_array(h1,h2,p2,p1) -= 0.5d0 * c_1 * phase - big_array(h2,p2,h1,p1) += 0.5d0 * c_1 * phase - big_array(h2,p1,h1,p2) -= 0.5d0 * c_1 * phase + big_array(h2,h1,p2,p1) += 0.5d0 * c_1 * phase + big_array(h2,h1,p1,p2) -= 0.5d0 * c_1 * phase enddo endif end @@ -420,11 +418,11 @@ p1 =exc(1,2) p2 =exc(2,2) if(beta_beta.or.spin_trace)then - big_array(h1,p1,h2,p2) += 0.5d0 * c_1* phase - big_array(h1,p2,h2,p1) -= 0.5d0 * c_1* phase + big_array(h1,h2,p1,p2) += 0.5d0 * c_1* phase + big_array(h1,h2,p2,p1) -= 0.5d0 * c_1* phase - big_array(h2,p2,h1,p1) += 0.5d0 * c_1* phase - big_array(h2,p1,h1,p2) -= 0.5d0 * c_1* phase + big_array(h2,h1,p2,p1) += 0.5d0 * c_1* phase + big_array(h2,h1,p1,p2) -= 0.5d0 * c_1* phase endif end From de7e1f70950bebaa8315aa2c12046ef4814a4478 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 28 Jun 2019 16:51:16 +0200 Subject: [PATCH 3/6] added test for energy --- src/casscf/get_energy.irp.f | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 src/casscf/get_energy.irp.f diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f new file mode 100644 index 00000000..a7b53f13 --- /dev/null +++ b/src/casscf/get_energy.irp.f @@ -0,0 +1,21 @@ +program print_2rdm + implicit none + read_wf = .True. + touch read_wf + integer :: i,j,k,l + double precision :: accu(4),twodm,thr,act_twodm2,integral,get_two_e_integral + thr = 1.d-10 + + accu = 0.d0 + do l = 1, mo_num + do k = 1, mo_num + do j = 1, mo_num + do i = 1, mo_num + integral = get_two_e_integral(i,j,k,l,mo_integrals_map) + accu(1) += act_two_rdm_spin_trace_mo(i,j,k,l) * integral + enddo + enddo + enddo + enddo + print*,'accu = ',accu(1) +end From 78fe995f55360da3130ed1a1cba9ac0793e0276a Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 28 Jun 2019 20:45:07 +0200 Subject: [PATCH 4/6] getting there with active orbitals --- src/bitmask/bitmasks_routines.irp.f | 33 ++- src/casscf/get_energy.irp.f | 31 ++- src/two_body_rdm/general_2rdm_routines.irp.f | 42 ++-- src/two_body_rdm/orb_range_2_rdm.irp.f | 10 +- .../routines_compute_2rdm_orb_range.irp.f | 189 +++++++++++++++--- 5 files changed, 248 insertions(+), 57 deletions(-) diff --git a/src/bitmask/bitmasks_routines.irp.f b/src/bitmask/bitmasks_routines.irp.f index 378a3dcd..5c4bf347 100644 --- a/src/bitmask/bitmasks_routines.irp.f +++ b/src/bitmask/bitmasks_routines.irp.f @@ -33,7 +33,7 @@ subroutine bitstring_to_list( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint) @@ -213,3 +213,34 @@ subroutine print_spindet(string,Nint) print *, trim(output(1)) end + +logical function is_integer_in_string(bite,string,Nint) + use bitmasks + implicit none + integer, intent(in) :: bite,Nint + integer(bit_kind), intent(in) :: string(Nint) + integer(bit_kind) :: string_bite(Nint) + integer :: i,itot,itot_and + character*(2048) :: output(1) + string_bite = 0_bit_kind + call set_bit_to_integer(bite,string_bite,Nint) + itot = 0 + itot_and = 0 + is_integer_in_string = .False. +!print*,'' +!print*,'' +!print*,'bite = ',bite +!call bitstring_to_str( output(1), string_bite, Nint ) +! print *, trim(output(1)) +!call bitstring_to_str( output(1), string, Nint ) +! print *, trim(output(1)) + do i = 1, Nint + itot += popcnt(string(i)) + itot_and += popcnt(ior(string(i),string_bite(i))) + enddo +!print*,'itot,itot_and',itot,itot_and + if(itot == itot_and)then + is_integer_in_string = .True. + endif +!pause +end diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f index a7b53f13..29a12cad 100644 --- a/src/casscf/get_energy.irp.f +++ b/src/casscf/get_energy.irp.f @@ -2,20 +2,39 @@ program print_2rdm implicit none read_wf = .True. touch read_wf + call routine +end + +subroutine routine integer :: i,j,k,l + integer :: ii,jj,kk,ll double precision :: accu(4),twodm,thr,act_twodm2,integral,get_two_e_integral thr = 1.d-10 + accu = 0.d0 - do l = 1, mo_num - do k = 1, mo_num - do j = 1, mo_num - do i = 1, mo_num + do ll = 1, n_act_orb + l = list_act(ll) + do kk = 1, n_act_orb + k = list_act(kk) + do jj = 1, n_act_orb + j = list_act(jj) + do ii = 1, n_act_orb + i = list_act(ii) integral = get_two_e_integral(i,j,k,l,mo_integrals_map) - accu(1) += act_two_rdm_spin_trace_mo(i,j,k,l) * integral + accu(1) += act_two_rdm_spin_trace_mo(ii,jj,kk,ll) * integral + !if(dabs(act_two_rdm_spin_trace_mo(ii,jj,kk,ll)).gt.thr)then + !print*,'',ii,jj,kk,ll,act_two_rdm_spin_trace_mo(ii,jj,kk,ll)*integral + !print*,'accu',accu(1) + !endif enddo enddo enddo enddo - print*,'accu = ',accu(1) + print*,'accu = ',accu(1) + print*,'psi_energy_two_e = ',psi_energy_two_e +!double precision :: hij +!call i_H_j_double_alpha_beta(psi_det(1,1,1),psi_det(1,1,2),N_int,hij) +!print*,'hij * 2',hij * psi_coef(1,1) * psi_coef(2,1) * 2.d0 +!print*,'psi diag = ',psi_energy_two_e - hij * psi_coef(1,1) * psi_coef(2,1) * 2.d0 end diff --git a/src/two_body_rdm/general_2rdm_routines.irp.f b/src/two_body_rdm/general_2rdm_routines.irp.f index a9fcd61a..0157c46b 100644 --- a/src/two_body_rdm/general_2rdm_routines.irp.f +++ b/src/two_body_rdm/general_2rdm_routines.irp.f @@ -1,4 +1,4 @@ -subroutine orb_range_two_rdm_dm_nstates_openmp(big_array,dim1,norb,list_orb,state_weights,ispin,u_0,N_st,sze) +subroutine orb_range_two_rdm_dm_nstates_openmp(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_0,N_st,sze) use bitmasks implicit none BEGIN_DOC @@ -13,6 +13,7 @@ subroutine orb_range_two_rdm_dm_nstates_openmp(big_array,dim1,norb,list_orb,stat END_DOC integer, intent(in) :: N_st,sze integer, intent(in) :: dim1,norb,list_orb(norb),ispin + integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) double precision, intent(in) :: u_0(sze,N_st),state_weights(N_st) @@ -30,7 +31,7 @@ subroutine orb_range_two_rdm_dm_nstates_openmp(big_array,dim1,norb,list_orb,stat size(u_t, 1), & N_det, N_st) - call orb_range_two_rdm_dm_nstates_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1) + call orb_range_two_rdm_dm_nstates_openmp_work(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) do k=1,N_st @@ -39,7 +40,7 @@ subroutine orb_range_two_rdm_dm_nstates_openmp(big_array,dim1,norb,list_orb,stat end -subroutine orb_range_two_rdm_dm_nstates_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine orb_range_two_rdm_dm_nstates_openmp_work(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC @@ -49,23 +50,25 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work(big_array,dim1,norb,list_orb END_DOC integer, intent(in) :: N_st,sze,istart,iend,ishift,istep integer, intent(in) :: dim1,norb,list_orb(norb),ispin + integer, intent(in) :: list_orb_reverse(mo_num) 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_dm_nstates_openmp_work_1(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_dm_nstates_openmp_work_1(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (2) - call orb_range_two_rdm_dm_nstates_openmp_work_2(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_dm_nstates_openmp_work_2(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (3) - call orb_range_two_rdm_dm_nstates_openmp_work_3(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_dm_nstates_openmp_work_3(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case (4) - call orb_range_two_rdm_dm_nstates_openmp_work_4(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_dm_nstates_openmp_work_4(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) case default - call orb_range_two_rdm_dm_nstates_openmp_work_N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) + call orb_range_two_rdm_dm_nstates_openmp_work_N_int(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) end select end @@ -73,7 +76,7 @@ end BEGIN_TEMPLATE -subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) +subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC @@ -89,6 +92,7 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,l 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 + integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) integer :: i,j,k,l @@ -112,6 +116,7 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,l double precision :: c_average logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + integer(bit_kind) :: orb_bitmask($N_int) alpha_alpha = .False. beta_beta = .False. alpha_beta = .False. @@ -129,7 +134,10 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,l print*,'ispin = ',ispin stop endif + + PROVIDE N_int + call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 allocate(idx0(maxab)) @@ -242,7 +250,7 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,l c_2(l) = u_t(l,k_a) c_average += c_1(l) * c_2(l) * state_weights(l) enddo - call orb_range_off_diagonal_double_to_two_rdm_ab_dm(tmp_det,tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + call orb_range_off_diagonal_double_to_two_rdm_ab_dm(tmp_det,tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) enddo endif @@ -319,9 +327,9 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,l enddo if(alpha_beta.or.spin_trace.or.alpha_alpha)then ! increment the alpha/beta part for single excitations - call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) ! increment the alpha/alpha part for single excitations - call orb_range_off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + call orb_range_off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) endif enddo @@ -344,7 +352,7 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,l c_2(l) = u_t(l,k_a) c_average += c_1(l) * c_2(l) * state_weights(l) enddo - call orb_range_off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,big_array,dim1,norb,list_orb,ispin) + call orb_range_off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) enddo endif @@ -411,9 +419,9 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,l enddo if(alpha_beta.or.spin_trace.or.beta_beta)then ! increment the alpha/beta part for single excitations - call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + call orb_range_off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) ! increment the beta /beta part for single excitations - call orb_range_off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_average,big_array,dim1,norb,list_orb,ispin) + call orb_range_off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) endif enddo @@ -435,7 +443,7 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,l c_2(l) = u_t(l,k_a) c_average += c_1(l) * c_2(l) * state_weights(l) enddo - call orb_range_off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_average,big_array,dim1,norb,list_orb,ispin) + call orb_range_off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) ASSERT (l_a <= N_det) enddo @@ -467,7 +475,7 @@ subroutine orb_range_two_rdm_dm_nstates_openmp_work_$N_int(big_array,dim1,norb,l c_average += c_1(l) * c_1(l) * state_weights(l) enddo - call orb_range_diagonal_contrib_to_all_two_rdm_dm(tmp_det,c_average,big_array,dim1,norb,list_orb,ispin) + call orb_range_diagonal_contrib_to_all_two_rdm_dm(tmp_det,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) end do !!$OMP END DO diff --git a/src/two_body_rdm/orb_range_2_rdm.irp.f b/src/two_body_rdm/orb_range_2_rdm.irp.f index 621f6c4b..e98612c5 100644 --- a/src/two_body_rdm/orb_range_2_rdm.irp.f +++ b/src/two_body_rdm/orb_range_2_rdm.irp.f @@ -10,7 +10,7 @@ ! condition for alpha/beta spin ispin = 1 act_two_rdm_alpha_alpha_mo = 0.D0 - call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) END_PROVIDER @@ -23,7 +23,7 @@ ! condition for alpha/beta spin ispin = 2 act_two_rdm_beta_beta_mo = 0.d0 - call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_beta_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_beta_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) END_PROVIDER @@ -41,7 +41,7 @@ ispin = 3 print*,'ispin = ',ispin act_two_rdm_alpha_beta_mo = 0.d0 - call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_alpha_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_alpha_beta_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) END_PROVIDER @@ -55,7 +55,9 @@ ! condition for alpha/beta spin ispin = 4 act_two_rdm_spin_trace_mo = 0.d0 - call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + integer :: i + + call orb_range_two_rdm_dm_nstates_openmp(act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) END_PROVIDER diff --git a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f index d918932a..c2283fb2 100644 --- a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f +++ b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f @@ -1,14 +1,15 @@ - subroutine orb_range_diagonal_contrib_to_two_rdm_ab_dm(det_1,c_1,big_array,dim1,norb,list_orb) + subroutine orb_range_diagonal_contrib_to_two_rdm_ab_dm(det_1,c_1,big_array,dim1,orb_bitmask) use bitmasks BEGIN_DOC ! routine that update the DIAGONAL PART of the alpha/beta two body rdm in a specific range of orbitals ! c_1 is supposed to be a scalar quantity, such as state averaged coef END_DOC implicit none - integer, intent(in) :: dim1,norb,list_orb(norb) + integer, intent(in) :: dim1 double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) 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 integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) @@ -24,21 +25,32 @@ end - subroutine orb_range_diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array,dim1,norb,list_orb,ispin) + subroutine orb_range_diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC ! routine that update the DIAGONAL PART of ALL THREE two body rdm END_DOC implicit none - integer, intent(in) :: dim1,norb,list_orb(norb),ispin + integer, intent(in) :: dim1,ispin + integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) 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 integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) integer :: i,j,h1,h2,istate + integer(bit_kind) :: det_1_act(N_int,2) logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + do i = 1, N_int + det_1_act(i,1) = iand(det_1(i,1),orb_bitmask(i)) + det_1_act(i,2) = iand(det_1(i,2),orb_bitmask(i)) + enddo + +!print*,'ahah' +!call debug_det(det_1_act,N_int) +!pause alpha_alpha = .False. beta_beta = .False. alpha_beta = .False. @@ -55,29 +67,43 @@ 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) + 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) - h1 = occ(i,1) + i1 = occ(i,1) +! if(.not.is_integer_in_string(i1,orb_bitmask,N_int))cycle do j = 1, n_occ_ab(2) - h2 = occ(j,2) +! if(.not.is_integer_in_string(i2,orb_bitmask,N_int))cycle + i2 = occ(j,2) + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) big_array(h1,h2,h1,h2) += c_1 enddo enddo else if (alpha_alpha)then do i = 1, n_occ_ab(1) - h1 = occ(i,1) + i1 = occ(i,1) +! if(.not.is_integer_in_string(i1,orb_bitmask,N_int))cycle do j = 1, n_occ_ab(1) - h2 = occ(j,1) + i2 = occ(j,1) +! if(.not.is_integer_in_string(i2,orb_bitmask,N_int))cycle + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) big_array(h1,h2,h1,h2) += 0.5d0 * c_1 big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 enddo enddo else if (beta_beta)then do i = 1, n_occ_ab(2) - h1 = occ(i,2) + i1 = occ(i,2) +! if(.not.is_integer_in_string(i1,orb_bitmask,N_int))cycle do j = 1, n_occ_ab(2) - h2 = occ(j,2) + i2 = occ(j,2) +! if(.not.is_integer_in_string(i2,orb_bitmask,N_int))cycle + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) big_array(h1,h2,h1,h2) += 0.5d0 * c_1 big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 enddo @@ -85,25 +111,38 @@ else if(spin_trace)then ! 0.5 * (alpha beta + beta alpha) do i = 1, n_occ_ab(1) - h1 = occ(i,1) + i1 = occ(i,1) +! if(.not.is_integer_in_string(i1,orb_bitmask,N_int))cycle do j = 1, n_occ_ab(2) - h2 = occ(j,2) + i2 = occ(j,2) +! if(.not.is_integer_in_string(i2,orb_bitmask,N_int))cycle + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) big_array(h1,h2,h1,h2) += 0.5d0 * (c_1 ) big_array(h2,h1,h2,h1) += 0.5d0 * (c_1 ) enddo enddo + !stop do i = 1, n_occ_ab(1) - h1 = occ(i,1) + i1 = occ(i,1) +! if(.not.is_integer_in_string(i1,orb_bitmask,N_int))cycle do j = 1, n_occ_ab(1) - h2 = occ(j,1) + i2 = occ(j,1) +! if(.not.is_integer_in_string(i2,orb_bitmask,N_int))cycle + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) big_array(h1,h2,h1,h2) += 0.5d0 * c_1 big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 enddo enddo do i = 1, n_occ_ab(2) - h1 = occ(i,2) + i1 = occ(i,2) +! if(.not.is_integer_in_string(i1,orb_bitmask,N_int))cycle do j = 1, n_occ_ab(2) - h2 = occ(j,2) + i2 = occ(j,2) +! if(.not.is_integer_in_string(i2,orb_bitmask,N_int))cycle + h1 = list_orb_reverse(i1) + h2 = list_orb_reverse(i2) big_array(h1,h2,h1,h2) += 0.5d0 * c_1 big_array(h1,h2,h2,h1) -= 0.5d0 * c_1 enddo @@ -112,20 +151,23 @@ end - subroutine orb_range_off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + subroutine orb_range_off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC ! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS END_DOC implicit none - integer, intent(in) :: dim1,norb,list_orb(norb),ispin + integer, intent(in) :: dim1,ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(in) :: c_1 integer :: i,j,h1,h2,p1,p2,istate integer :: exc(0:2,2,2) double precision :: phase logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string alpha_alpha = .False. beta_beta = .False. alpha_beta = .False. @@ -139,28 +181,52 @@ else if(ispin == 4)then spin_trace = .True. endif +!print*,'' +!do i = 1, mo_num +! print*,'list_orb',i,list_orb_reverse(i) +!enddo call get_double_excitation(det_1,det_2,exc,phase,N_int) h1 = exc(1,1,1) +!print*,'h1',h1 + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) +!print*,'passed h1 = ',h1 h2 = exc(1,1,2) +!print*,'h2',h2 + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return + h2 = list_orb_reverse(h2) +!print*,'passed h2 = ',h2 p1 = exc(1,2,1) +!print*,'p1',p1 + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) +!print*,'passed p1 = ',p1 p2 = exc(1,2,2) +!print*,'p2',p2 + if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return + p2 = list_orb_reverse(p2) +!print*,'passed p2 = ',p2 if(alpha_beta)then big_array(h1,h2,p1,p2) += c_1 * phase else if(spin_trace)then big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase big_array(p1,p2,h1,h2) += 0.5d0 * c_1 * phase + !print*,'h1,h2,p1,p2',h1,h2,p1,p2 + !print*,'',big_array(h1,h2,p1,p2) endif end - subroutine orb_range_off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + subroutine orb_range_off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC ! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS END_DOC implicit none - integer, intent(in) :: dim1,norb,list_orb(norb),ispin + integer, intent(in) :: dim1,ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(in) :: c_1 integer :: occ(N_int*bit_kind_size,2) @@ -170,6 +236,7 @@ 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. @@ -190,17 +257,29 @@ if (exc(0,1,1) == 1) then ! Mono alpha h1 = exc(1,1,1) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) p1 = exc(1,2,1) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(2) h2 = occ(i,2) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) big_array(h1,h2,p1,h2) += c_1 * phase enddo else ! Mono beta h1 = exc(1,1,2) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) p1 = exc(1,2,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(1) h2 = occ(i,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) big_array(h2,h1,h2,p1) += c_1 * phase enddo endif @@ -208,18 +287,30 @@ if (exc(0,1,1) == 1) then ! Mono alpha h1 = exc(1,1,1) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) p1 = exc(1,2,1) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(2) h2 = occ(i,2) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase enddo else ! Mono beta h1 = exc(1,1,2) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) p1 = exc(1,2,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(1) h2 = occ(i,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase enddo @@ -227,15 +318,17 @@ endif end - subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) BEGIN_DOC ! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS END_DOC use bitmasks implicit none - integer, intent(in) :: dim1,norb,list_orb(norb),ispin + integer, intent(in) :: dim1,ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(in) :: c_1 integer :: occ(N_int*bit_kind_size,2) @@ -245,6 +338,7 @@ 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. @@ -265,9 +359,15 @@ if (exc(0,1,1) == 1) then ! Mono alpha h1 = exc(1,1,1) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) p1 = exc(1,2,1) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(1) h2 = occ(i,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase @@ -280,15 +380,17 @@ endif end - subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC ! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for SINGLE EXCITATIONS END_DOC implicit none - integer, intent(in) :: dim1,norb,list_orb(norb),ispin + integer, intent(in) :: dim1,ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(in) :: c_1 @@ -298,6 +400,7 @@ 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. @@ -321,10 +424,16 @@ else ! Mono beta h1 = exc(1,1,2) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) p1 = exc(1,2,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) do istate = 1, N_states do i = 1, n_occ_ab(2) h2 = occ(i,2) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase @@ -337,15 +446,17 @@ end - subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC ! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS END_DOC implicit none - integer, intent(in) :: dim1,norb,list_orb(norb),ispin + integer, intent(in) :: dim1,ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(in) :: c_1 integer :: i,j,h1,h2,p1,p2,istate @@ -353,6 +464,7 @@ 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. @@ -368,9 +480,17 @@ endif call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) h1 =exc(1,1) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) h2 =exc(2,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return + h2 = list_orb_reverse(h2) p1 =exc(1,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) p2 =exc(2,2) + if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return + p2 = list_orb_reverse(p2) if(alpha_alpha.or.spin_trace)then do istate = 1, N_states big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase @@ -382,22 +502,25 @@ endif end - subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,norb,list_orb,ispin) + subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) 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,norb,list_orb(norb),ispin + integer, intent(in) :: dim1,ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) + integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(in) :: c_1 integer :: i,j,h1,h2,p1,p2,istate integer :: exc(0:2,2) double precision :: phase logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string alpha_alpha = .False. beta_beta = .False. alpha_beta = .False. @@ -414,9 +537,17 @@ call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) h1 =exc(1,1) + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return + h1 = list_orb_reverse(h1) h2 =exc(2,1) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return + h2 = list_orb_reverse(h2) p1 =exc(1,2) + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return + p1 = list_orb_reverse(p1) p2 =exc(2,2) + if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return + p2 = list_orb_reverse(p2) if(beta_beta.or.spin_trace)then big_array(h1,h2,p1,p2) += 0.5d0 * c_1* phase big_array(h1,h2,p2,p1) -= 0.5d0 * c_1* phase From 57eabff6758b254d9c6e92b04e356205509ecc1d Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 29 Jun 2019 17:29:32 +0200 Subject: [PATCH 5/6] added documentation for the two-rdm --- src/casscf/densities.irp.f | 8 +- src/casscf/get_energy.irp.f | 14 +- src/two_body_rdm/ab_only_routines.irp.f | 22 +-- src/two_body_rdm/all_2rdm_routines.irp.f | 2 +- src/two_body_rdm/orb_range_2_rdm.irp.f | 20 +++ ...outines.irp.f => orb_range_routines.irp.f} | 0 src/two_body_rdm/routines_compute_2rdm.irp.f | 14 +- .../routines_compute_2rdm_orb_range.irp.f | 129 +++++++++++++++++- 8 files changed, 168 insertions(+), 41 deletions(-) rename src/two_body_rdm/{general_2rdm_routines.irp.f => orb_range_routines.irp.f} (100%) diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f index 7b243bb4..30a914f1 100644 --- a/src/casscf/densities.irp.f +++ b/src/casscf/densities.irp.f @@ -42,7 +42,7 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22 if (bavard) then - write(6,*) ' providing density matrix P0' + write(6,*) ' providing the 2 body RDM on the active part' endif P0tuvx= 0.d0 @@ -55,11 +55,7 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] uu = list_act(u) do t = 1, n_act_orb tt = list_act(t) - P0tuvx(t,u,v,x) = & - state_average_weight(istate) * & - ( two_rdm_alpha_beta_mo (tt,uu,vv,xx,istate) + & - two_rdm_alpha_alpha_mo(tt,uu,vv,xx,istate) + & - two_rdm_beta_beta_mo (tt,uu,vv,xx,istate) ) + P0tuvx(t,u,v,x) = act_two_rdm_spin_trace_mo(t,v,u,x) enddo enddo enddo diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f index 29a12cad..0a5cfb49 100644 --- a/src/casscf/get_energy.irp.f +++ b/src/casscf/get_energy.irp.f @@ -1,5 +1,10 @@ program print_2rdm implicit none + BEGIN_DOC + ! get the active part of the bielectronic energy on a given wave function. + ! + ! useful to test the active part of the spin trace 2 rdms + END_DOC read_wf = .True. touch read_wf call routine @@ -23,18 +28,9 @@ subroutine routine i = list_act(ii) integral = get_two_e_integral(i,j,k,l,mo_integrals_map) accu(1) += act_two_rdm_spin_trace_mo(ii,jj,kk,ll) * integral - !if(dabs(act_two_rdm_spin_trace_mo(ii,jj,kk,ll)).gt.thr)then - !print*,'',ii,jj,kk,ll,act_two_rdm_spin_trace_mo(ii,jj,kk,ll)*integral - !print*,'accu',accu(1) - !endif enddo enddo enddo enddo print*,'accu = ',accu(1) - print*,'psi_energy_two_e = ',psi_energy_two_e -!double precision :: hij -!call i_H_j_double_alpha_beta(psi_det(1,1,1),psi_det(1,1,2),N_int,hij) -!print*,'hij * 2',hij * psi_coef(1,1) * psi_coef(2,1) * 2.d0 -!print*,'psi diag = ',psi_energy_two_e - hij * psi_coef(1,1) * psi_coef(2,1) * 2.d0 end diff --git a/src/two_body_rdm/ab_only_routines.irp.f b/src/two_body_rdm/ab_only_routines.irp.f index 195f439a..9041c753 100644 --- a/src/two_body_rdm/ab_only_routines.irp.f +++ b/src/two_body_rdm/ab_only_routines.irp.f @@ -1,9 +1,9 @@ - subroutine two_rdm_dm_nstates_openmp(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze) + subroutine two_rdm_ab_nstates_openmp(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze) use bitmasks implicit none BEGIN_DOC - ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + ! Computes the alpha/beta part of the two-body density matrix IN CHEMIST NOTATIONS ! ! Assumes that the determinants are in psi_det ! @@ -27,7 +27,7 @@ size(u_t, 1), & N_det, N_st) - call two_rdm_dm_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1) + call two_rdm_ab_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) do k=1,N_st @@ -37,11 +37,11 @@ end - subroutine two_rdm_dm_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + subroutine two_rdm_ab_nstates_openmp_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none BEGIN_DOC - ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + ! Computes the alpha/beta part of the two-body density matrix ! ! Default should be 1,N_det,0,1 END_DOC @@ -55,20 +55,20 @@ select case (N_int) case (1) - call two_rdm_dm_nstates_openmp_work_1(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_work_1(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (2) - call two_rdm_dm_nstates_openmp_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (3) - call two_rdm_dm_nstates_openmp_work_3(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_work_3(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case (4) - call two_rdm_dm_nstates_openmp_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) case default - call two_rdm_dm_nstates_openmp_work_N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + call two_rdm_ab_nstates_openmp_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_dm_nstates_openmp_work_$N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) + subroutine two_rdm_ab_nstates_openmp_work_$N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks implicit none integer, intent(in) :: N_st,sze,istart,iend,ishift,istep diff --git a/src/two_body_rdm/all_2rdm_routines.irp.f b/src/two_body_rdm/all_2rdm_routines.irp.f index 75d71ded..3f08b18f 100644 --- a/src/two_body_rdm/all_2rdm_routines.irp.f +++ b/src/two_body_rdm/all_2rdm_routines.irp.f @@ -2,7 +2,7 @@ subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab, use bitmasks implicit none BEGIN_DOC - ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + ! 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 ! diff --git a/src/two_body_rdm/orb_range_2_rdm.irp.f b/src/two_body_rdm/orb_range_2_rdm.irp.f index e98612c5..c40c46d2 100644 --- a/src/two_body_rdm/orb_range_2_rdm.irp.f +++ b/src/two_body_rdm/orb_range_2_rdm.irp.f @@ -4,6 +4,10 @@ BEGIN_PROVIDER [double precision, act_two_rdm_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] implicit none double precision, allocatable :: state_weights(:) + BEGIN_DOC +! 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 +! = + END_DOC allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) integer :: ispin @@ -17,6 +21,10 @@ BEGIN_PROVIDER [double precision, act_two_rdm_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] implicit none double precision, allocatable :: state_weights(:) + BEGIN_DOC +! 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 +! = + END_DOC allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) integer :: ispin @@ -30,6 +38,10 @@ BEGIN_PROVIDER [double precision, act_two_rdm_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] implicit none double precision, allocatable :: state_weights(:) + BEGIN_DOC +! 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 +! = + END_DOC allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) integer :: ispin @@ -48,6 +60,14 @@ BEGIN_PROVIDER [double precision, act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] implicit none + BEGIN_DOC +! 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: +! +! \sum_{i,j,k,l = 1, n_act_orb} 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 double precision, allocatable :: state_weights(:) allocate(state_weights(N_states)) state_weights = 1.d0/dble(N_states) diff --git a/src/two_body_rdm/general_2rdm_routines.irp.f b/src/two_body_rdm/orb_range_routines.irp.f similarity index 100% rename from src/two_body_rdm/general_2rdm_routines.irp.f rename to src/two_body_rdm/orb_range_routines.irp.f diff --git a/src/two_body_rdm/routines_compute_2rdm.irp.f b/src/two_body_rdm/routines_compute_2rdm.irp.f index 7165576f..112d2e36 100644 --- a/src/two_body_rdm/routines_compute_2rdm.irp.f +++ b/src/two_body_rdm/routines_compute_2rdm.irp.f @@ -3,7 +3,7 @@ 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 +! 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 @@ -31,7 +31,7 @@ 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 +! 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 @@ -77,7 +77,7 @@ 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 +! 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 @@ -101,7 +101,7 @@ 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 +! 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 @@ -140,7 +140,7 @@ 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 +! 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 @@ -177,7 +177,7 @@ 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 +! 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 @@ -214,7 +214,7 @@ 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 +! 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 diff --git a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f index c2283fb2..a3c7a76d 100644 --- a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f +++ b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f @@ -28,7 +28,20 @@ subroutine orb_range_diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the DIAGONAL PART of ALL THREE two body rdm +! 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 +! +! 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 END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -154,7 +167,24 @@ subroutine orb_range_off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS +! 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 +! +! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation +! +! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +! +! ispin determines which spin-spin component of the two-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 3 or 4 will do something END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -219,7 +249,24 @@ subroutine orb_range_off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for +! +! a given couple of determinant det_1, det_2 being a SINGLE excitation with respect to one another +! +! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +! +! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation +! +! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +! +! ispin determines which spin-spin component of the two-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 3 or 4 will do something END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -320,7 +367,24 @@ subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS +! 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 @@ -383,7 +447,24 @@ subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for SINGLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for +! +! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another +! +! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +! +! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation +! +! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +! +! ispin determines which spin-spin component of the two-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 2 or 4 will do something END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -449,7 +530,24 @@ subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS +! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for +! +! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another +! +! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 +! +! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation +! +! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals +! +! ispin determines which spin-spin component of the two-rdm you will update +! +! ispin == 1 :: alpha/ alpha +! ispin == 2 :: beta / beta +! ispin == 3 :: alpha/ beta +! ispin == 4 :: spin traced <=> total two-rdm +! +! here, only ispin == 1 or 4 will do something END_DOC implicit none integer, intent(in) :: dim1,ispin @@ -505,7 +603,24 @@ subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) use bitmasks BEGIN_DOC -! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for DOUBLE EXCITATIONS +! 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 From 3c9728be99bc54ec98c217936e03861e8652a7c5 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 29 Jun 2019 17:34:20 +0200 Subject: [PATCH 6/6] comments --- src/casscf/densities.irp.f | 7 ++++--- src/two_body_rdm/README.rst | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f index 30a914f1..3cfd7583 100644 --- a/src/casscf/densities.irp.f +++ b/src/casscf/densities.irp.f @@ -19,14 +19,15 @@ END_PROVIDER BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] BEGIN_DOC - ! the second-order density matrix in the basis of the starting MOs - ! matrices are state averaged + ! The second-order density matrix in the basis of the starting MOs ONLY IN THE RANGE OF ACTIVE MOS + ! The values are state averaged ! - ! we use the spin-free generators of mono-excitations + ! We use the spin-free generators of mono-excitations ! E_pq destroys q and creates p ! D_pq = <0|E_pq|0> = D_qp ! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0> ! + ! P0tuvx(p,q,r,s) = chemist notation : 1/2 <0|E_pq E_rs - delta_qr E_ps|0> END_DOC implicit none integer :: t,u,v,x diff --git a/src/two_body_rdm/README.rst b/src/two_body_rdm/README.rst index ea5839e8..978240c9 100644 --- a/src/two_body_rdm/README.rst +++ b/src/two_body_rdm/README.rst @@ -3,6 +3,6 @@ two_body_rdm ============ Contains the two rdms $\alpha\alpha$, $\beta\beta$ and $\alpha\beta$ stored as -maps, with pysicists notation, consistent with the two-electron integrals in the +arrays, with pysicists notation, consistent with the two-electron integrals in the MO basis.