10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-26 07:02:14 +02:00

added OMP - excitations as integer2

This commit is contained in:
Yann Garniron 2016-05-20 11:27:39 +02:00
parent 33bd506328
commit f8ece7d40b
3 changed files with 123 additions and 32 deletions

View File

@ -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

View File

@ -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(:,:,:)

View File

@ -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