diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 8873a940..42058145 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -288,7 +288,8 @@ logical function is_generable(det1, det2, Nint) implicit none integer, intent(in) :: Nint integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2) - integer :: degree, f, exc(0:2, 2, 2), h1, h2, p1, p2, s1, s2, t + integer :: degree, f, exc(0:2, 2, 2), t + integer*2 :: h1, h2, p1, p2, s1, s2 integer, external :: searchExc logical, external :: excEq double precision :: phase @@ -312,7 +313,7 @@ logical function is_generable(det1, det2, Nint) ! ! print *, f ! return - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + call decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2) if(degree == 1) then h2 = h1 @@ -323,14 +324,14 @@ logical function is_generable(det1, det2, Nint) s1 = 0 end if - if(h1 + s1*mo_tot_num < h2 + s2*mo_tot_num) then + if(h1 + (s1-1)*mo_tot_num < h2 + (s2-1)*mo_tot_num) then f = searchExc(hh_exists(1,1), (/s1, h1, s2, h2/), hh_shortcut(0)) else f = searchExc(hh_exists(1,1), (/s2, h2, s1, h1/), hh_shortcut(0)) end if if(f == -1) return - if(p1 + s1*mo_tot_num < p2 + s2*mo_tot_num) then + if(p1 + (s1-1)*mo_tot_num < p2 + (s2-1)*mo_tot_num) then f = searchExc(pp_exists(1,hh_shortcut(f)), (/s1, p1, s2, p2/), hh_shortcut(f+1)-hh_shortcut(f)) else f = searchExc(pp_exists(1,hh_shortcut(f)), (/s2, p2, s1, p1/), hh_shortcut(f+1)-hh_shortcut(f)) @@ -376,7 +377,7 @@ integer function searchExc(excs, exc, n) use bitmasks integer, intent(in) :: n - integer,intent(in) :: excs(4,n), exc(4) + integer*2,intent(in) :: excs(4,n), exc(4) integer :: l, h, c integer, external :: excCmp logical, external :: excEq @@ -441,8 +442,8 @@ subroutine sort_exc(key, N_key) integer, intent(in) :: N_key - integer,intent(inout) :: key(4,N_key) - integer :: tmp(4) + integer*2,intent(inout) :: key(4,N_key) + integer*2 :: tmp(4) integer :: i,ni @@ -464,7 +465,7 @@ end subroutine logical function exc_inf(exc1, exc2) implicit none - integer,intent(in) :: exc1(4), exc2(4) + integer*2,intent(in) :: exc1(4), exc2(4) integer :: i exc_inf = .false. do i=1,4 @@ -486,9 +487,9 @@ subroutine tamise_exc(key, no, n, N_key) ! Uncodumented : TODO END_DOC integer,intent(in) :: no, n, N_key - integer,intent(inout) :: key(4, N_key) + integer*2,intent(inout) :: key(4, N_key) integer :: k,j - integer :: tmp(4) + integer*2 :: tmp(4) logical :: exc_inf integer :: ni @@ -518,7 +519,7 @@ end subroutine subroutine dec_exc(exc, h1, h2, p1, p2) implicit none integer :: exc(0:2,2,2), s1, s2, degree - integer, intent(out) :: h1, h2, p1, p2 + integer*2, intent(out) :: h1, h2, p1, p2 degree = exc(0,1,1) + exc(0,1,2) @@ -529,7 +530,7 @@ subroutine dec_exc(exc, h1, h2, p1, p2) if(degree == 0) return - call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2) + call decode_exc_int2(exc, degree, h1, p1, h2, p2, s1, s2) h1 += mo_tot_num * (s1-1) p1 += mo_tot_num * (s1-1) @@ -556,12 +557,13 @@ subroutine dec_exc(exc, h1, h2, p1, p2) end subroutine - BEGIN_PROVIDER [ integer, hh_exists, (4, N_det_ref * N_det_non_ref) ] + BEGIN_PROVIDER [ integer*2, hh_exists, (4, N_det_ref * N_det_non_ref) ] &BEGIN_PROVIDER [ integer, hh_shortcut, (0:N_det_ref * N_det_non_ref + 1) ] -&BEGIN_PROVIDER [ integer, pp_exists, (4, N_det_ref * N_det_non_ref) ] +&BEGIN_PROVIDER [ integer*2, pp_exists, (4, N_det_ref * N_det_non_ref) ] implicit none - integer,allocatable :: num(:,:) - integer :: exc(0:2, 2, 2), degree, n, on, s, h1, h2, p1, p2, l, i + integer*2,allocatable :: num(:,:) + integer :: exc(0:2, 2, 2), degree, n, on, s, l, i + integer*2 :: h1, h2, p1, p2 double precision :: phase logical, external :: excEq @@ -587,19 +589,19 @@ end subroutine hh_shortcut(0) = 1 hh_shortcut(1) = 1 - hh_exists(:,1) = (/1, num(1,1), 1, num(2,1)/) - pp_exists(:,1) = (/1, num(3,1), 1, num(4,1)/) + hh_exists(:,1) = (/1_2, num(1,1), 1_2, num(2,1)/) + pp_exists(:,1) = (/1_2, num(3,1), 1_2, num(4,1)/) s = 1 do i=2,n if(.not. excEq(num(1,i), num(1,s))) then s += 1 num(:, s) = num(:, i) - pp_exists(:,s) = (/1, num(3,s), 1, num(4,s)/) + pp_exists(:,s) = (/1_2, num(3,s), 1_2, num(4,s)/) if(hh_exists(2, hh_shortcut(0)) /= num(1,s) .or. & hh_exists(4, hh_shortcut(0)) /= num(2,s)) then hh_shortcut(0) += 1 hh_shortcut(hh_shortcut(0)) = s - hh_exists(:,hh_shortcut(0)) = (/1, num(1,s), 1, num(2,s)/) + hh_exists(:,hh_shortcut(0)) = (/1_2, num(1,s), 1_2, num(2,s)/) end if end if end do @@ -629,7 +631,7 @@ END_PROVIDER logical function excEq(exc1, exc2) implicit none - integer, intent(in) :: exc1(4), exc2(4) + integer*2, intent(in) :: exc1(4), exc2(4) integer :: i excEq = .false. do i=1, 4 @@ -641,7 +643,7 @@ end function integer function excCmp(exc1, exc2) implicit none - integer, intent(in) :: exc1(4), exc2(4) + integer*2, intent(in) :: exc1(4), exc2(4) integer :: i excCmp = 0 do i=1, 4 @@ -660,8 +662,8 @@ subroutine apply_hole(det, exc, res, ok, Nint) use bitmasks implicit none integer, intent(in) :: Nint - integer, intent(in) :: exc(4) - integer :: s1, s2, h1, h2 + integer*2, intent(in) :: exc(4) + integer*2 :: s1, s2, h1, h2 integer(bit_kind),intent(in) :: det(Nint, 2) integer(bit_kind),intent(out) :: res(Nint, 2) logical, intent(out) :: ok @@ -695,8 +697,8 @@ subroutine apply_particle(det, exc, res, ok, Nint) use bitmasks implicit none integer, intent(in) :: Nint - integer, intent(in) :: exc(4) - integer :: s1, s2, p1, p2 + integer*2, intent(in) :: exc(4) + integer*2 :: s1, s2, p1, p2 integer(bit_kind),intent(in) :: det(Nint, 2) integer(bit_kind),intent(out) :: res(Nint, 2) logical, intent(out) :: ok diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 5d8d4850..4f99f6e1 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -6,16 +6,28 @@ use bitmasks &BEGIN_PROVIDER [ double precision, delta_ii_mrcc, (N_states, N_det_ref) ] use bitmasks implicit none - integer :: gen, h, p, i_state, n, t, i, h1, h2, p1, p2, s1, s2 - integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2), buf(N_int, 2, N_det_non_ref) + integer :: gen, h, p, i_state, n, t, i, h1, h2, p1, p2, s1, s2, iproc + integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2) + integer(bit_kind), allocatable :: buf(:,:,:) + double precision, allocatable :: delta_ij_mwen(:,:,:,:), delta_ii_mwen(:,:,:) logical :: ok logical, external :: detEq delta_ij_mrcc = 0d0 delta_ii_mrcc = 0d0 i_state = 1 + provide hh_shortcut psi_det_size + allocate(delta_ij_mwen(N_states,N_det_non_ref,N_det_ref,nproc), delta_ii_mwen(N_states,N_det_ref,nproc)) + allocate(buf(N_int, 2, N_det_non_ref)) + delta_ij_mwen = 0d0 + delta_ii_mwen = 0d0 + !$OMP PARALLEL DO default(none) schedule(dynamic) & + !$OMP shared(psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) & + !$OMP shared(N_states, N_det_non_ref, N_det_ref, delta_ii_mwen, delta_ij_mwen) & + !$OMP private(h, n, mask, omask, buf, ok, iproc) do gen= 1, N_det_generators + iproc = omp_get_thread_num() + 1 print *, gen, "/", N_det_generators do h=1, hh_shortcut(0) call apply_hole(psi_det_generators(1,1,gen), hh_exists(1, h), mask, ok, N_int) @@ -28,12 +40,23 @@ use bitmasks if(ok) n = n + 1 end do n = n - 1 - if(n /= 0) call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc,gen,n,buf,N_int,omask) + if(n /= 0) call mrcc_part_dress(delta_ij_mwen(1,1,1,iproc), delta_ii_mwen(1,1,iproc),gen,n,buf,N_int,omask) end do end do + !$OMP END PARALLEL DO + do iproc=1, nproc + delta_ij_mrcc = delta_ij_mrcc + delta_ij_mwen(:,:,:,iproc) + delta_ii_mrcc = delta_ii_mrcc + delta_ii_mwen(:,:,iproc) + end do END_PROVIDER +! subroutine blit(b1, b2) +! double precision :: b1(N_states,N_det_non_ref,N_det_ref), b2(N_states,N_det_non_ref,N_det_ref) +! b1 = b1 + b2 +! end subroutine + + subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffer,Nint,key_mask) use bitmasks implicit none @@ -463,7 +486,7 @@ END_PROVIDER double precision :: Hjk, Hki, Hij, pre(N_det_ref), wall integer :: i_state, degree, npre, ipre(N_det_ref), npres(N_det_ref) - provide lambda_mrcc +! provide lambda_mrcc npres = 0 delta_cas = 0d0 call wall_time(wall) @@ -605,8 +628,8 @@ end function call wall_time(wall) print *, "cepa0", wall - provide det_cepa0_active delta_cas lambda_mrcc - provide mo_bielec_integrals_in_map +! provide det_cepa0_active delta_cas lambda_mrcc +! provide mo_bielec_integrals_in_map allocate(idx_sorted_bit(N_det)) sortRef(:,:,:) = det_ref_active(:,:,:) diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 2b5ae4f1..3374dfb2 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -139,6 +139,72 @@ subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) end +subroutine decode_exc_int2(exc,degree,h1,p1,h2,p2,s1,s2) + use bitmasks + implicit none + BEGIN_DOC + ! Decodes the exc arrays returned by get_excitation. + ! h1,h2 : Holes + ! p1,p2 : Particles + ! s1,s2 : Spins (1:alpha, 2:beta) + ! degree : Degree of excitation + END_DOC + integer, intent(in) :: exc(0:2,2,2),degree + integer*2, intent(out) :: h1,h2,p1,p2,s1,s2 + ASSERT (degree > 0) + ASSERT (degree < 3) + + select case(degree) + case(2) + if (exc(0,1,1) == 2) then + h1 = exc(1,1,1) + h2 = exc(2,1,1) + p1 = exc(1,2,1) + p2 = exc(2,2,1) + s1 = 1 + s2 = 1 + else if (exc(0,1,2) == 2) then + h1 = exc(1,1,2) + h2 = exc(2,1,2) + p1 = exc(1,2,2) + p2 = exc(2,2,2) + s1 = 2 + s2 = 2 + else + h1 = exc(1,1,1) + h2 = exc(1,1,2) + p1 = exc(1,2,1) + p2 = exc(1,2,2) + s1 = 1 + s2 = 2 + endif + case(1) + if (exc(0,1,1) == 1) then + h1 = exc(1,1,1) + h2 = 0 + p1 = exc(1,2,1) + p2 = 0 + s1 = 1 + s2 = 0 + else + h1 = exc(1,1,2) + h2 = 0 + p1 = exc(1,2,2) + p2 = 0 + s1 = 2 + s2 = 0 + endif + case(0) + h1 = 0 + p1 = 0 + h2 = 0 + p2 = 0 + s1 = 0 + s2 = 0 + end select +end + + subroutine get_double_excitation(det1,det2,exc,phase,Nint) use bitmasks implicit none