From 78fe995f55360da3130ed1a1cba9ac0793e0276a Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 28 Jun 2019 20:45:07 +0200 Subject: [PATCH] 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