LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-18 12:03:57 +01:00
Code Issues Releases Wiki Activity

Davidson optimization

Browse Source
This commit is contained in:
Anthony Scemama 2015-11-13 17:36:29 +01:00
parent b7e20b3493
commit 6a640567e8
3 changed files with 184 additions and 135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

85
src/Determinants/davidson.irp.f
View File

@ -90,27 +90,36 @@ end function
subroutine tamiser(key, idx, no, n, Nint, N_key) subroutine tamiser(key, idx, no, n, Nint, N_key)
use bitmasks use bitmasks
implicit none implicit none
integer(bit_kind),intent(inout) :: key(Nint, 2, N_key)
BEGIN_DOC
! Uncodumented : TODO
END_DOC
integer,intent(in) :: no, n, Nint, N_key integer,intent(in) :: no, n, Nint, N_key
integer(bit_kind),intent(inout) :: key(Nint, 2, N_key)
integer,intent(inout) :: idx(N_key) integer,intent(inout) :: idx(N_key)
integer :: k,j,tmpidx integer :: k,j,tmpidx
integer(bit_kind) :: tmp(Nint, 2) integer(bit_kind) :: tmp(Nint, 2)
logical :: det_inf logical :: det_inf
integer :: ni
k = no k = no
j = 2*k j = 2*k
do while(j <= n) do while(j <= n)
if(j < n) then if(j < n) then
if (det_inf(key(:,:,j), key(:,:,j+1), Nint)) then if (det_inf(key(1,1,j), key(1,1,j+1), Nint)) then
j = j+1 j = j+1
endif endif
endif endif
if(det_inf(key(:,:,k), key(:,:,j), Nint)) then if(det_inf(key(1,1,k), key(1,1,j), Nint)) then
tmp(:,:) = key(:,:,k) do ni=1,Nint
key(:,:,k) = key(:,:,j) tmp(ni,1) = key(ni,1,k)
key(:,:,j) = tmp(:,:) tmp(ni,2) = key(ni,2,k)
key(ni,1,k) = key(ni,1,j)
key(ni,2,k) = key(ni,2,j)
key(ni,1,j) = tmp(ni,1)
key(ni,2,j) = tmp(ni,2)
enddo
tmpidx = idx(k) tmpidx = idx(k)
idx(k) = idx(j) idx(k) = idx(j)
idx(j) = tmpidx idx(j) = tmpidx
@ -126,17 +135,25 @@ end subroutine
subroutine sort_dets_ba_v(key_in, key_out, idx, shortcut, version, N_key, Nint) subroutine sort_dets_ba_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
use bitmasks use bitmasks
implicit none implicit none
integer, intent(in) :: Nint, N_key
integer(bit_kind),intent(in) :: key_in(Nint,2,N_key) integer(bit_kind),intent(in) :: key_in(Nint,2,N_key)
integer(bit_kind) :: key(Nint,2,N_key) integer(bit_kind) :: key(Nint,2,N_key)
integer(bit_kind),intent(out) :: key_out(Nint,N_key) integer(bit_kind),intent(out) :: key_out(Nint,N_key)
integer,intent(out) :: idx(N_key) integer,intent(out) :: idx(N_key)
integer,intent(out) :: shortcut(0:N_key+1) integer,intent(out) :: shortcut(0:N_key+1)
integer(bit_kind),intent(out) :: version(Nint,N_key+1) integer(bit_kind),intent(out) :: version(Nint,N_key+1)
integer, intent(in) :: Nint, N_key
integer(bit_kind) :: tmp(Nint, 2,N_key) integer(bit_kind) :: tmp(Nint, 2,N_key)
integer :: i,ni
key(:,1,:N_key) = key_in(:,2,:N_key) BEGIN_DOC
key(:,2,:N_key) = key_in(:,1,:N_key) ! Uncodumented : TODO
END_DOC
do i=1,N_key
do ni=1,Nint
key(ni,1,i) = key_in(ni,2,i)
key(ni,2,i) = key_in(ni,1,i)
enddo
enddo
call sort_dets_ab_v(key, key_out, idx, shortcut, version, N_key, Nint) call sort_dets_ab_v(key, key_out, idx, shortcut, version, N_key, Nint)
@ -148,18 +165,24 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
use bitmasks use bitmasks
implicit none implicit none
BEGIN_DOC
! Uncodumented : TODO
END_DOC
integer, intent(in) :: Nint, N_key
integer(bit_kind),intent(in) :: key_in(Nint,2,N_key) integer(bit_kind),intent(in) :: key_in(Nint,2,N_key)
integer(bit_kind) :: key(Nint,2,N_key) integer(bit_kind) :: key(Nint,2,N_key)
integer(bit_kind),intent(out) :: key_out(Nint,N_key) integer(bit_kind),intent(out) :: key_out(Nint,N_key)
integer,intent(out) :: idx(N_key) integer,intent(out) :: idx(N_key)
integer,intent(out) :: shortcut(0:N_key+1) integer,intent(out) :: shortcut(0:N_key+1)
integer(bit_kind),intent(out) :: version(Nint,N_key+1) integer(bit_kind),intent(out) :: version(Nint,N_key+1)
integer, intent(in) :: Nint, N_key
integer(bit_kind) :: tmp(Nint, 2) integer(bit_kind) :: tmp(Nint, 2)
integer :: tmpidx,i,ni integer :: tmpidx,i,ni
key(:,:,:) = key_in(:,:,:)
do i=1,N_key do i=1,N_key
do ni=1,Nint
key(ni,1,i) = key_in(ni,1,i)
key(ni,2,i) = key_in(ni,2,i)
enddo
idx(i) = i idx(i) = i
end do end do
@ -168,9 +191,14 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
end do end do
do i=N_key,2,-1 do i=N_key,2,-1
tmp(:,:) = key(:,:,i) do ni=1,Nint
key(:,:,i) = key(:,:,1) tmp(ni,1) = key(ni,1,i)
key(:,:,1) = tmp(:,:) tmp(ni,2) = key(ni,2,i)
key(ni,1,i) = key(ni,1,1)
key(ni,2,i) = key(ni,2,1)
key(ni,1,1) = tmp(ni,1)
key(ni,2,1) = tmp(ni,2)
enddo
tmpidx = idx(i) tmpidx = idx(i)
idx(i) = idx(1) idx(i) = idx(1)
idx(1) = tmpidx idx(1) = tmpidx
@ -179,7 +207,9 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
shortcut(0) = 1 shortcut(0) = 1
shortcut(1) = 1 shortcut(1) = 1
version(:,1) = key(:,1,1) do ni=1,Nint
version(ni,1) = key(ni,1,1)
enddo
do i=2,N_key do i=2,N_key
do ni=1,nint do ni=1,nint
if(key(ni,1,i) /= key(ni,1,i-1)) then if(key(ni,1,i) /= key(ni,1,i-1)) then
@ -191,15 +221,22 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
end do end do
end do end do
shortcut(shortcut(0)+1) = N_key+1 shortcut(shortcut(0)+1) = N_key+1
key_out(:,:) = key(:,2,:) do i=1,N_key
do ni=1,Nint
key_out(ni,i) = key(ni,2,i)
enddo
enddo
end subroutine end subroutine
c
subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint) subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint)
use bitmasks use bitmasks
implicit none implicit none
BEGIN_DOC
! Uncodumented : TODO
END_DOC
integer(bit_kind),intent(inout) :: key(Nint,2,N_key) integer(bit_kind),intent(inout) :: key(Nint,2,N_key)
integer,intent(out) :: idx(N_key) integer,intent(out) :: idx(N_key)
integer,intent(out) :: shortcut(0:N_key+1) integer,intent(out) :: shortcut(0:N_key+1)
@ -216,9 +253,15 @@ subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint)
end do end do
do i=N_key,2,-1 do i=N_key,2,-1
tmp(:,:) = key(:,:,i) do ni=1,Nint
key(:,:,i) = key(:,:,1) tmp(ni,1) = key(ni,1,i)
key(:,:,1) = tmp(:,:) tmp(ni,2) = key(ni,2,i)
key(ni,1,i) = key(ni,1,1)
key(ni,2,i) = key(ni,2,1)
key(ni,1,1) = tmp(ni,1)
key(ni,2,1) = tmp(ni,2)
enddo
tmpidx = idx(i) tmpidx = idx(i)
idx(i) = idx(1) idx(i) = idx(1)
idx(1) = tmpidx idx(1) = tmpidx

20
src/Determinants/s2.irp.f
View File

@ -126,16 +126,12 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
s2 = 0.d0 s2 = 0.d0
davidson_threshold_bis = davidson_threshold davidson_threshold_bis = davidson_threshold
call sort_dets_ab_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int)
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,s2_tmp,idx,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass) & !$OMP PRIVATE(i,j,s2_tmp,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass)&
!$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold,shortcut,sorted,sort_idx,version)& !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold,shortcut,sorted,sort_idx,version)&
!$OMP REDUCTION(+:s2) !$OMP REDUCTION(+:s2)
allocate(idx(0:n))
!$OMP SINGLE
call sort_dets_ab_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int)
!$OMP END SINGLE
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0) do sh=1,shortcut(0)
@ -177,9 +173,14 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP SINGLE !$OMP END PARALLEL
call sort_dets_ba_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int) call sort_dets_ba_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int)
!$OMP END SINGLE
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,s2_tmp,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass)&
!$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold,shortcut,sorted,sort_idx,version)&
!$OMP REDUCTION(+:s2)
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0) do sh=1,shortcut(0)
@ -204,7 +205,6 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
enddo enddo
!$OMP END DO !$OMP END DO
deallocate(idx)
!$OMP END PARALLEL !$OMP END PARALLEL
s2 = s2+s2 s2 = s2+s2
do i=1,n do i=1,n

42
src/Determinants/slater_rules.irp.f
View File

@ -1074,8 +1074,8 @@ double precision function diag_H_mat_elem(det_in,Nint)
particle(i,2) = iand(hole(i,2),det_in(i,2)) particle(i,2) = iand(hole(i,2),det_in(i,2))
hole(i,1) = iand(hole(i,1),ref_bitmask(i,1)) hole(i,1) = iand(hole(i,1),ref_bitmask(i,1))
hole(i,2) = iand(hole(i,2),ref_bitmask(i,2)) hole(i,2) = iand(hole(i,2),ref_bitmask(i,2))
nexc(1) += popcnt(hole(i,1)) nexc(1) = nexc(1) + popcnt(hole(i,1))
nexc(2) += popcnt(hole(i,2)) nexc(2) = nexc(2) + popcnt(hole(i,2))
enddo enddo
diag_H_mat_elem = ref_bitmask_energy diag_H_mat_elem = ref_bitmask_energy
@ -1239,10 +1239,11 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
integer :: shortcut(0:n+1), sort_idx(n) integer :: shortcut(0:n+1), sort_idx(n)
integer(bit_kind) :: sorted(Nint,n), version(Nint,n) integer(bit_kind) :: sorted(Nint,n), version(Nint,n)
integer(bit_kind) :: sorted_i(Nint)
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi
! double precision :: local_threshold
ASSERT (Nint > 0) ASSERT (Nint > 0)
@ -1254,9 +1255,9 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint) call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi) & !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold,sorted_i)&
!$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold,sorted,shortcut,sort_idx,version) !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold,sorted,shortcut,sort_idx,version)
allocate(idx(0:n), vt(n)) allocate(vt(n))
Vt = 0.d0 Vt = 0.d0
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
@ -1264,28 +1265,32 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
do sh2=1,sh do sh2=1,sh
exa = 0 exa = 0
do ni=1,Nint do ni=1,Nint
exa += popcnt(xor(version(ni,sh), version(ni,sh2))) exa = exa + popcnt(xor(version(ni,sh), version(ni,sh2)))
end do end do
if(exa > 2) then if(exa > 2) then
cycle cycle
end if end if
do i=shortcut(sh),shortcut(sh+1)-1 do i=shortcut(sh),shortcut(sh+1)-1
org_i = sort_idx(i)
local_threshold = davidson_threshold - dabs(u_0(org_i))
if(sh==sh2) then if(sh==sh2) then
endi = i-1 endi = i-1
else else
endi = shortcut(sh2+1)-1 endi = shortcut(sh2+1)-1
end if end if
do ni=1,Nint
sorted_i(ni) = sorted(ni,i)
enddo
do j=shortcut(sh2),endi do j=shortcut(sh2),endi
org_j = sort_idx(j)
if ( dabs(u_0(org_j)) > local_threshold ) then
ext = exa ext = exa
do ni=1,Nint do ni=1,Nint
ext += popcnt(xor(sorted(ni,i), sorted(ni,j))) ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j)))
end do end do
if(ext <= 4) then if(ext <= 4) then
org_i = sort_idx(i)
org_j = sort_idx(j)
if ( dabs(u_0(org_j)) + dabs(u_0(org_i)) > davidson_threshold ) then
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
vt (org_i) = vt (org_i) + hij*u_0(org_j) vt (org_i) = vt (org_i) + hij*u_0(org_j)
vt (org_j) = vt (org_j) + hij*u_0(org_i) vt (org_j) = vt (org_j) + hij*u_0(org_i)
@ -1303,29 +1308,30 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
enddo enddo
!$OMP END CRITICAL !$OMP END CRITICAL
deallocate(idx,vt) deallocate(vt)
!$OMP END PARALLEL !$OMP END PARALLEL
call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint) call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi) & !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold)&
!$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold,sorted,shortcut,sort_idx,version) !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold,sorted,shortcut,sort_idx,version)
allocate(idx(0:n), vt(n)) allocate(vt(n))
Vt = 0.d0 Vt = 0.d0
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0) do sh=1,shortcut(0)
do i=shortcut(sh),shortcut(sh+1)-1 do i=shortcut(sh),shortcut(sh+1)-1
local_threshold = davidson_threshold - dabs(u_0(org_i))
org_i = sort_idx(i)
do j=shortcut(sh),i-1 do j=shortcut(sh),i-1
org_j = sort_idx(j)
if ( dabs(u_0(org_j)) > local_threshold ) then
ext = 0 ext = 0
do ni=1,Nint do ni=1,Nint
ext += popcnt(xor(sorted(ni,i), sorted(ni,j))) ext = ext + popcnt(xor(sorted(ni,i), sorted(ni,j)))
end do end do
if(ext == 4) then if(ext == 4) then
org_i = sort_idx(i)
org_j = sort_idx(j)
if ( dabs(u_0(org_j)) + dabs(u_0(org_i)) > davidson_threshold ) then
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
vt (org_i) = vt (org_i) + hij*u_0(org_j) vt (org_i) = vt (org_i) + hij*u_0(org_j)
vt (org_j) = vt (org_j) + hij*u_0(org_i) vt (org_j) = vt (org_j) + hij*u_0(org_i)
@ -1341,7 +1347,7 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
v_0(i) = v_0(i) + vt(i) v_0(i) = v_0(i) + vt(i)
enddo enddo
!$OMP END CRITICAL !$OMP END CRITICAL
deallocate(idx,vt) deallocate(vt)
!$OMP END PARALLEL !$OMP END PARALLEL
do i=1,n do i=1,n