From 92d7bbd57ebc1a9f631b471e1042a952a397b81e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 10 Feb 2017 00:50:37 +0100 Subject: [PATCH] Fixed PT2 stoch: fragments broken --- plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f | 10 +++- plugins/Full_CI_ZMQ/selection.irp.f | 55 +++++++++++--------- src/Bitmask/bitmasks.irp.f | 10 ++-- src/Davidson/u0Hu0.irp.f | 15 +++++- 4 files changed, 58 insertions(+), 32 deletions(-) diff --git a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f index d7c98933..ac98bd9c 100644 --- a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f +++ b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f @@ -192,7 +192,9 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su actually_computed(:) = computed(:) 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) actually_computed(tbc(i)) = .false. @@ -356,6 +358,10 @@ end subroutine 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 END_PROVIDER @@ -522,7 +528,7 @@ end subroutine iloc = N_det_generators do i=comb_teeth, 1, -1 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 stato -= comb_step end do diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index de7c93f8..96b47e53 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/plugins/Full_CI_ZMQ/selection.irp.f @@ -2,7 +2,8 @@ use bitmasks BEGIN_PROVIDER [ integer, fragment_count ] 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 @@ -44,7 +45,7 @@ subroutine assert(cond, msg) logical, intent(in) :: cond if(.not. cond) then - print *, "assert fail: "//msg + print *, "assert failed: "//msg stop end if end @@ -286,7 +287,7 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p logical :: fullMatch, ok 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(:, :, :) logical :: monoAdo, monoBdo; @@ -296,7 +297,7 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p monoBdo = .true. 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 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(particle, particle_list, N_particles, N_int) - ! ====== - ! If the subset doesn't exist, return - logical :: will_compute - will_compute = subset == 0 - - 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) - will_compute = (maskInd >= subset) - if (.not.will_compute) then - return - endif - endif - ! ====== +! ! ====== +! ! If the subset doesn't exist, return +! logical :: will_compute +! will_compute = subset == 0 +! +! 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) +! will_compute = (maskInd >= subset) +! if (.not.will_compute) then +! return +! endif +! endif +! ! ====== 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 prefullinteresting(0) += 1 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 do @@ -413,23 +418,23 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p do ii=1,prefullinteresting(0) i = prefullinteresting(ii) nt = 0 - mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i)) - mobMask(1,2) = iand(negMask(1,2), psi_det_sorted(1,2,i)) + mobMask(1,1) = iand(negMask(1,1), prefullinteresting_det(1,1,ii)) + mobMask(1,2) = iand(negMask(1,2), prefullinteresting_det(1,1,ii)) nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) do j=2,N_int - mobMask(j,1) = iand(negMask(j,1), psi_det_sorted(j,1,i)) - mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i)) + mobMask(j,1) = iand(negMask(j,1), prefullinteresting_det(j,1,ii)) + mobMask(j,2) = iand(negMask(j,2), prefullinteresting_det(j,2,ii)) nt = nt+ popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) end do if(nt <= 2) then fullinteresting(0) += 1 fullinteresting(fullinteresting(0)) = i - fullminilist(1,1,fullinteresting(0)) = psi_det_sorted(1,1,i) - fullminilist(1,2,fullinteresting(0)) = psi_det_sorted(1,2,i) + fullminilist(1,1,fullinteresting(0)) = prefullinteresting_det(1,1,ii) + fullminilist(1,2,fullinteresting(0)) = prefullinteresting_det(1,2,ii) do j=2,N_int - fullminilist(j,1,fullinteresting(0)) = psi_det_sorted(j,1,i) - fullminilist(j,2,fullinteresting(0)) = psi_det_sorted(j,2,i) + fullminilist(j,1,fullinteresting(0)) = prefullinteresting_det(j,1,ii) + fullminilist(j,2,fullinteresting(0)) = prefullinteresting_det(j,2,ii) enddo end if end do diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 6fe201ff..10ab6f67 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -6,6 +6,7 @@ BEGIN_PROVIDER [ integer, N_int ] ! Number of 64-bit integers needed to represent determinants as binary strings END_DOC N_int = (mo_tot_num-1)/bit_kind_size + 1 + call write_int(6,N_int, 'N_int') END_PROVIDER @@ -386,6 +387,8 @@ END_PROVIDER n_virt_orb += popcnt(virt_bitmask(i,1)) enddo 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 @@ -559,10 +562,11 @@ END_PROVIDER integer :: i,j n_core_orb = 0 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,2) = xor(HF_bitmask(i,2),ior(reunion_of_cas_inact_bitmask(i,2),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(full_ijkl_bitmask(i),ior(reunion_of_cas_inact_bitmask(i,2),virt_bitmask(i,1))) n_core_orb += popcnt(core_bitmask(i,1)) enddo + call write_int(6,n_core_orb,'Number of core MOs') END_PROVIDER @@ -597,7 +601,7 @@ BEGIN_PROVIDER [ integer, n_act_orb] do i = 1, N_int n_act_orb += popcnt(cas_bitmask(i,1,1)) enddo - print*,'n_act_orb = ',n_act_orb + call write_int(6,n_act_orb, 'Number of active MOs') END_PROVIDER BEGIN_PROVIDER [integer, list_act, (n_act_orb)] diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index 026921d0..d9481886 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -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(:,:) integer :: i,j,k,l, jj,ii integer :: i0, j0 + logical, allocatable :: utloop(:) integer, allocatable :: shortcut(:,:), sort_idx(:,:) 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 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 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 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)) Vt = 0.d0 St = 0.d0 !$OMP DO do i=1,n + utloop(i) = .False. do istate=1,N_st 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 + utloop(i) = .not.utloop(i) enddo !$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 org_i = sort_idx(i,2) do j=shortcut(sh,2),shortcut(sh+1,2)-1 + if (utloop(j)) cycle org_j = sort_idx(j,2) ext = popcnt(xor(sorted(1,i,2), sorted(1,j,2))) 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 do i=1,n + utloop(i) = .False. do istate=1,N_st 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 + utloop(i) = .not.utloop(i) enddo !$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 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))) if (ext > 4) cycle 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 do j=shortcut(sh,1),i-1 + if (utloop(j)) cycle ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) if (ext > 4) cycle 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 do j=i+1,shortcut(sh+1,1)-1 + if (utloop(j)) cycle ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) if (ext > 4) cycle do ni=2,Nint