LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-29 16:34:50 +02:00
Code Issues Releases Wiki Activity

Fixed PT2 stoch: fragments broken

Browse Source
This commit is contained in:
Anthony Scemama 2017-02-10 00:50:37 +01:00
parent 67f8fc8e3e
commit 92d7bbd57e
4 changed files with 58 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
View File

@ -192,7 +192,9 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
actually_computed(:) = computed(:) actually_computed(:) = computed(:)
parts_to_get(:) = 1 parts_to_get(:) = 1
if(fragment_first > 0) parts_to_get(1:fragment_first) = fragment_count if(fragment_first > 0) then
parts_to_get(1:fragment_first) = fragment_count
endif
do i=1,tbc(0) do i=1,tbc(0)
actually_computed(tbc(i)) = .false. actually_computed(tbc(i)) = .false.
@ -356,6 +358,10 @@ end subroutine
BEGIN_PROVIDER [ integer, size_tbc ] BEGIN_PROVIDER [ integer, size_tbc ]
implicit none
BEGIN_DOC
! Size of the tbc array
END_DOC
size_tbc = N_det_generators + fragment_count*fragment_first size_tbc = N_det_generators + fragment_count*fragment_first
END_PROVIDER END_PROVIDER
@ -522,7 +528,7 @@ end subroutine
iloc = N_det_generators iloc = N_det_generators
do i=comb_teeth, 1, -1 do i=comb_teeth, 1, -1
integer :: iloc integer :: iloc
iloc = pt2_find(stato, pt2_cweight, N_det_generators, 0, iloc) iloc = pt2_find(stato, pt2_cweight, N_det_generators, 1, iloc)
first_det_of_teeth(i) = iloc first_det_of_teeth(i) = iloc
stato -= comb_step stato -= comb_step
end do end do

55
plugins/Full_CI_ZMQ/selection.irp.f
View File

@ -2,7 +2,8 @@ use bitmasks
BEGIN_PROVIDER [ integer, fragment_count ] BEGIN_PROVIDER [ integer, fragment_count ]
implicit none implicit none
fragment_count = (elec_alpha_num-n_core_orb)**2 ! fragment_count = (elec_alpha_num-n_core_orb)**2
fragment_count = 1
END_PROVIDER END_PROVIDER
@ -44,7 +45,7 @@ subroutine assert(cond, msg)
logical, intent(in) :: cond logical, intent(in) :: cond
if(.not. cond) then if(.not. cond) then
print *, "assert fail: "//msg print *, "assert failed: "//msg
stop stop
end if end if
end end
@ -286,7 +287,7 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
logical :: fullMatch, ok logical :: fullMatch, ok
integer(bit_kind) :: mobMask(N_int, 2), negMask(N_int, 2) integer(bit_kind) :: mobMask(N_int, 2), negMask(N_int, 2)
integer,allocatable :: preinteresting(:), prefullinteresting(:), interesting(:), fullinteresting(:) integer,allocatable :: preinteresting(:), prefullinteresting(:), prefullinteresting_det(:,:,:), interesting(:), fullinteresting(:)
integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :)
logical :: monoAdo, monoBdo; logical :: monoAdo, monoBdo;
@ -296,7 +297,7 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
monoBdo = .true. monoBdo = .true.
allocate(minilist(N_int, 2, N_det_selectors), fullminilist(N_int, 2, N_det)) allocate(minilist(N_int, 2, N_det_selectors), fullminilist(N_int, 2, N_det))
allocate(preinteresting(0:N_det_selectors), prefullinteresting(0:N_det), interesting(0:N_det_selectors), fullinteresting(0:N_det)) allocate(preinteresting(0:N_det_selectors), prefullinteresting(0:N_det), interesting(0:N_det_selectors), fullinteresting(0:N_det), prefullinteresting_det(N_int,2,N_det))
do k=1,N_int do k=1,N_int
hole (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1)) hole (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1))
@ -312,19 +313,19 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
call bitstring_to_list_ab(hole , hole_list , N_holes , N_int) call bitstring_to_list_ab(hole , hole_list , N_holes , N_int)
call bitstring_to_list_ab(particle, particle_list, N_particles, N_int) call bitstring_to_list_ab(particle, particle_list, N_particles, N_int)
! ====== ! ! ======
! If the subset doesn't exist, return ! ! If the subset doesn't exist, return
logical :: will_compute ! logical :: will_compute
will_compute = subset == 0 ! will_compute = subset == 0
!
if (.not.will_compute) then ! if (.not.will_compute) then
maskInd = N_holes(1)*N_holes(2) + N_holes(2)*((N_holes(2)-1)/2) + N_holes(1)*((N_holes(1)-1)/2) ! maskInd = N_holes(1)*N_holes(2) + N_holes(2)*((N_holes(2)-1)/2) + N_holes(1)*((N_holes(1)-1)/2)
will_compute = (maskInd >= subset) ! will_compute = (maskInd >= subset)
if (.not.will_compute) then ! if (.not.will_compute) then
return ! return
endif ! endif
endif ! endif
! ====== ! ! ======
integer(bit_kind), allocatable:: preinteresting_det(:,:,:) integer(bit_kind), allocatable:: preinteresting_det(:,:,:)
@ -359,6 +360,10 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
else if(nt <= 2) then else if(nt <= 2) then
prefullinteresting(0) += 1 prefullinteresting(0) += 1
prefullinteresting(prefullinteresting(0)) = i prefullinteresting(prefullinteresting(0)) = i
do j=1,N_int
prefullinteresting_det(j,1,prefullinteresting(0)) = psi_det_sorted(j,1,i)
prefullinteresting_det(j,2,prefullinteresting(0)) = psi_det_sorted(j,2,i)
enddo
end if end if
end if end if
end do end do
@ -413,23 +418,23 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
do ii=1,prefullinteresting(0) do ii=1,prefullinteresting(0)
i = prefullinteresting(ii) i = prefullinteresting(ii)
nt = 0 nt = 0
mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i)) mobMask(1,1) = iand(negMask(1,1), prefullinteresting_det(1,1,ii))
mobMask(1,2) = iand(negMask(1,2), psi_det_sorted(1,2,i)) mobMask(1,2) = iand(negMask(1,2), prefullinteresting_det(1,1,ii))
nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2))
do j=2,N_int do j=2,N_int
mobMask(j,1) = iand(negMask(j,1), psi_det_sorted(j,1,i)) mobMask(j,1) = iand(negMask(j,1), prefullinteresting_det(j,1,ii))
mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i)) mobMask(j,2) = iand(negMask(j,2), prefullinteresting_det(j,2,ii))
nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
end do end do
if(nt <= 2) then if(nt <= 2) then
fullinteresting(0) += 1 fullinteresting(0) += 1
fullinteresting(fullinteresting(0)) = i fullinteresting(fullinteresting(0)) = i
fullminilist(1,1,fullinteresting(0)) = psi_det_sorted(1,1,i) fullminilist(1,1,fullinteresting(0)) = prefullinteresting_det(1,1,ii)
fullminilist(1,2,fullinteresting(0)) = psi_det_sorted(1,2,i) fullminilist(1,2,fullinteresting(0)) = prefullinteresting_det(1,2,ii)
do j=2,N_int do j=2,N_int
fullminilist(j,1,fullinteresting(0)) = psi_det_sorted(j,1,i) fullminilist(j,1,fullinteresting(0)) = prefullinteresting_det(j,1,ii)
fullminilist(j,2,fullinteresting(0)) = psi_det_sorted(j,2,i) fullminilist(j,2,fullinteresting(0)) = prefullinteresting_det(j,2,ii)
enddo enddo
end if end if
end do end do

10
src/Bitmask/bitmasks.irp.f
View File

@ -6,6 +6,7 @@ BEGIN_PROVIDER [ integer, N_int ]
! Number of 64-bit integers needed to represent determinants as binary strings ! Number of 64-bit integers needed to represent determinants as binary strings
END_DOC END_DOC
N_int = (mo_tot_num-1)/bit_kind_size + 1 N_int = (mo_tot_num-1)/bit_kind_size + 1
call write_int(6,N_int, 'N_int')
END_PROVIDER END_PROVIDER
@ -386,6 +387,8 @@ END_PROVIDER
n_virt_orb += popcnt(virt_bitmask(i,1)) n_virt_orb += popcnt(virt_bitmask(i,1))
enddo enddo
endif endif
call write_int(6,n_inact_orb, 'Number of inactive MOs')
call write_int(6,n_virt_orb, 'Number of virtual MOs')
END_PROVIDER END_PROVIDER
@ -559,10 +562,11 @@ END_PROVIDER
integer :: i,j integer :: i,j
n_core_orb = 0 n_core_orb = 0
do i = 1, N_int do i = 1, N_int
core_bitmask(i,1) = xor(HF_bitmask(i,1),ior(reunion_of_cas_inact_bitmask(i,1),virt_bitmask(i,1))) core_bitmask(i,1) = xor(full_ijkl_bitmask(i),ior(reunion_of_cas_inact_bitmask(i,1),virt_bitmask(i,1)))
core_bitmask(i,2) = xor(HF_bitmask(i,2),ior(reunion_of_cas_inact_bitmask(i,2),virt_bitmask(i,1))) core_bitmask(i,2) = xor(full_ijkl_bitmask(i),ior(reunion_of_cas_inact_bitmask(i,2),virt_bitmask(i,1)))
n_core_orb += popcnt(core_bitmask(i,1)) n_core_orb += popcnt(core_bitmask(i,1))
enddo enddo
call write_int(6,n_core_orb,'Number of core MOs')
END_PROVIDER END_PROVIDER
@ -597,7 +601,7 @@ BEGIN_PROVIDER [ integer, n_act_orb]
do i = 1, N_int do i = 1, N_int
n_act_orb += popcnt(cas_bitmask(i,1,1)) n_act_orb += popcnt(cas_bitmask(i,1,1))
enddo enddo
print*,'n_act_orb = ',n_act_orb call write_int(6,n_act_orb, 'Number of active MOs')
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [integer, list_act, (n_act_orb)] BEGIN_PROVIDER [integer, list_act, (n_act_orb)]

15
src/Davidson/u0Hu0.irp.f
View File

@ -405,6 +405,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
double precision, allocatable :: vt(:,:), ut(:,:), st(:,:) double precision, allocatable :: vt(:,:), ut(:,:), st(:,:)
integer :: i,j,k,l, jj,ii integer :: i,j,k,l, jj,ii
integer :: i0, j0 integer :: i0, j0
logical, allocatable :: utloop(:)
integer, allocatable :: shortcut(:,:), sort_idx(:,:) integer, allocatable :: shortcut(:,:), sort_idx(:,:)
integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:) integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:)
@ -427,7 +428,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
PROVIDE ref_bitmask_energy PROVIDE ref_bitmask_energy
allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
allocate( ut(N_st_8,n)) allocate( ut(N_st_8,n), utloop(n) )
v_0 = 0.d0 v_0 = 0.d0
s_0 = 0.d0 s_0 = 0.d0
@ -437,16 +438,19 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,s2,j,k,jj,vt,st,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& !$OMP PRIVATE(i,hij,s2,j,k,jj,vt,st,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)&
!$OMP SHARED(n,keys_tmp,ut,Nint,u_0,v_0,s_0,sorted,shortcut,sort_idx,version,N_st,N_st_8) !$OMP SHARED(n,keys_tmp,ut,Nint,u_0,v_0,s_0,sorted,shortcut,sort_idx,version,N_st,N_st_8,utloop)
allocate(vt(N_st_8,n),st(N_st_8,n)) allocate(vt(N_st_8,n),st(N_st_8,n))
Vt = 0.d0 Vt = 0.d0
St = 0.d0 St = 0.d0
!$OMP DO !$OMP DO
do i=1,n do i=1,n
utloop(i) = .False.
do istate=1,N_st do istate=1,N_st
ut(istate,i) = u_0(sort_idx(i,2),istate) ut(istate,i) = u_0(sort_idx(i,2),istate)
utloop(i) = utloop(i) .or. (dabs(u_0(sort_idx(i,2),istate)) > 1.d-20)
enddo enddo
utloop(i) = .not.utloop(i)
enddo enddo
!$OMP END DO !$OMP END DO
@ -455,6 +459,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
do i=shortcut(sh,2),shortcut(sh+1,2)-1 do i=shortcut(sh,2),shortcut(sh+1,2)-1
org_i = sort_idx(i,2) org_i = sort_idx(i,2)
do j=shortcut(sh,2),shortcut(sh+1,2)-1 do j=shortcut(sh,2),shortcut(sh+1,2)-1
if (utloop(j)) cycle
org_j = sort_idx(j,2) org_j = sort_idx(j,2)
ext = popcnt(xor(sorted(1,i,2), sorted(1,j,2))) ext = popcnt(xor(sorted(1,i,2), sorted(1,j,2)))
if (ext > 4) cycle if (ext > 4) cycle
@ -477,9 +482,12 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
!$OMP DO !$OMP DO
do i=1,n do i=1,n
utloop(i) = .False.
do istate=1,N_st do istate=1,N_st
ut(istate,i) = u_0(sort_idx(i,1),istate) ut(istate,i) = u_0(sort_idx(i,1),istate)
utloop(i) = utloop(i) .or. (dabs(u_0(sort_idx(i,2),istate)) > 1.d-20)
enddo enddo
utloop(i) = .not.utloop(i)
enddo enddo
!$OMP END DO !$OMP END DO
@ -503,6 +511,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
enddo enddo
do j=shortcut(sh2,1),shortcut(sh2+1,1)-1 do j=shortcut(sh2,1),shortcut(sh2+1,1)-1
if (utloop(j)) cycle
ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1)))
if (ext > 4) cycle if (ext > 4) cycle
do ni=2,Nint do ni=2,Nint
@ -540,6 +549,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
enddo enddo
do j=shortcut(sh,1),i-1 do j=shortcut(sh,1),i-1
if (utloop(j)) cycle
ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1)))
if (ext > 4) cycle if (ext > 4) cycle
do ni=2,Nint do ni=2,Nint
@ -566,6 +576,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
enddo enddo
do j=i+1,shortcut(sh+1,1)-1 do j=i+1,shortcut(sh+1,1)-1
if (utloop(j)) cycle
ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1)))
if (ext > 4) cycle if (ext > 4) cycle
do ni=2,Nint do ni=2,Nint