diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 62d7c52c..6f40a809 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -312,9 +312,6 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d end do deallocate(indices) -! !$OMP CRITICAL -! print *, 'Step1: ', i_generator, preinteresting(0) -! !$OMP END CRITICAL allocate(banned(mo_num, mo_num,2), bannedOrb(mo_num, 2)) allocate (mat(N_states, mo_num, mo_num)) @@ -466,17 +463,8 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d fullinteresting(sze+1) = i end if end do - allocate (fullminilist (N_int, 2, fullinteresting(0)), & minilist (N_int, 2, interesting(0)) ) -! if(pert_2rdm)then -! allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) -! do i=1,fullinteresting(0) -! do j = 1, N_states -! coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) -! enddo -! enddo -! endif do i=1,fullinteresting(0) fullminilist(:,:,i) = psi_det_sorted(:,:,fullinteresting(i)) @@ -524,33 +512,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting) if(fullMatch) cycle -! !$OMP CRITICAL -! print *, 'Step3: ', i_generator, h1, interesting(0) -! !$OMP END CRITICAL call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) - -! if(.not.pert_2rdm)then - call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) -! else -! call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) -! endif + call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) end if enddo if(s1 /= s2) monoBdo = .false. enddo deallocate(fullminilist,minilist) -! if(pert_2rdm)then -! deallocate(coef_fullminilist_rev) -! endif enddo enddo deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) deallocate(banned, bannedOrb,mat) end subroutine - - - subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) use bitmasks use selection_types @@ -606,18 +580,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d ! to a determinant of the future. In that case, the determinant will be ! detected as already generated when generating in the future with a ! double excitation. -! -! if (.not.do_singles) then -! if ((h1 == p1) .or. (h2 == p2)) then -! cycle -! endif -! endif -! -! if (.not.do_doubles) then -! if ((h1 /= p1).and.(h2 /= p2)) then -! cycle -! endif -! endif ! ----- if(bannedOrb(p2, s2)) cycle @@ -974,13 +936,10 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere call get_mask_phase(psi_det_sorted(1,1,interesting(i)), phasemask,N_int) if(nt == 4) then -! call get_d2_reference(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) call get_d2(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) else if(nt == 3) then -! call get_d1_reference(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) else -! call get_d0_reference(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) end if else if(nt == 4) then @@ -1540,8 +1499,6 @@ subroutine past_d2(banned, p, sp) end if end - - subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting) use bitmasks implicit none diff --git a/src/cipsi_tc_bi_ortho/cipsi.irp.f b/src/cipsi_tc_bi_ortho/cipsi.irp.f index b1941068..fb907cb3 100644 --- a/src/cipsi_tc_bi_ortho/cipsi.irp.f +++ b/src/cipsi_tc_bi_ortho/cipsi.irp.f @@ -64,7 +64,7 @@ subroutine run_cipsi endif if (N_det > N_det_max) then - psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det) + psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det) psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states) N_det = N_det_max soft_touch N_det psi_det psi_coef diff --git a/src/cipsi_tc_bi_ortho/pt2.irp.f b/src/cipsi_tc_bi_ortho/pt2.irp.f index e7dca456..13b4dff4 100644 --- a/src/cipsi_tc_bi_ortho/pt2.irp.f +++ b/src/cipsi_tc_bi_ortho/pt2.irp.f @@ -52,7 +52,7 @@ subroutine pt2_tc_bi_ortho ! call routine_save_right if (N_det > N_det_max) then - psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det) + psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det) psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states) N_det = N_det_max soft_touch N_det psi_det psi_coef diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f index 633ca815..393023f2 100644 --- a/src/cipsi_tc_bi_ortho/selection.irp.f +++ b/src/cipsi_tc_bi_ortho/selection.irp.f @@ -91,7 +91,6 @@ subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset) end subroutine select_connected - double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2, Nint) use bitmasks implicit none @@ -136,7 +135,7 @@ double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2, Nint) end -subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock_diag_tmp, E0, pt2_data, buf, subset, csubset) +subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,buf,subset,csubset) use bitmasks use selection_types implicit none @@ -266,7 +265,6 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock enddo do k = 1, nmax - i = indices(k) mobMask(1,1) = iand(negMask(1,1), psi_det_sorted_tc(1,1,i)) mobMask(1,2) = iand(negMask(1,2), psi_det_sorted_tc(1,2,i)) @@ -304,10 +302,10 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock prefullinteresting(sze+1) = i endif endif - enddo deallocate(indices) + allocate( banned(mo_num, mo_num,2), bannedOrb(mo_num, 2) ) allocate( mat(N_states, mo_num, mo_num) ) allocate( mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num) ) @@ -463,17 +461,11 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock minilist (N_int, 2, interesting(0)) ) do i = 1, fullinteresting(0) - do k = 1, N_int - fullminilist(k,1,i) = psi_det_sorted_tc(k,1,fullinteresting(i)) - fullminilist(k,2,i) = psi_det_sorted_tc(k,2,fullinteresting(i)) - enddo + fullminilist(:,:,i) = psi_det_sorted_tc(:,:,fullinteresting(i)) enddo do i = 1, interesting(0) - do k = 1, N_int - minilist(k,1,i) = psi_det_sorted_tc(k,1,interesting(i)) - minilist(k,2,i) = psi_det_sorted_tc(k,2,interesting(i)) - enddo + minilist(:,:,i) = psi_det_sorted_tc(:,:,interesting(i)) enddo do s2 = s1, 2 @@ -516,196 +508,19 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock if(fullMatch) cycle call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting, mat_l, mat_r) - call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, mat_l, mat_r) endif - enddo - if(s1 /= s2) monoBdo = .false. enddo - deallocate(fullminilist, minilist) - enddo enddo - deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) deallocate(banned, bannedOrb,mat) deallocate(mat_l, mat_r) - end subroutine select_singles_and_doubles - -! --- - -subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting) - - use bitmasks - implicit none - - BEGIN_DOC - ! Identify the determinants in det which are in the internal space. These are - ! the determinants that can be produced by creating two particles on the mask. - END_DOC - - integer, intent(in) :: i_gen, N - integer, intent(in) :: interesting(0:N) - integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N) - logical, intent(inout) :: banned(mo_num, mo_num) - logical, intent(out) :: fullMatch - - integer :: i, j, na, nb, list(3) - integer(bit_kind) :: myMask(N_int, 2), negMask(N_int, 2) - - fullMatch = .false. - - do i=1,N_int - negMask(i,1) = not(mask(i,1)) - negMask(i,2) = not(mask(i,2)) - end do - - genl : do i=1, N - ! If det(i) can't be generated by the mask, cycle - do j=1, N_int - if(iand(det(j,1,i), mask(j,1)) /= mask(j, 1)) cycle genl - if(iand(det(j,2,i), mask(j,2)) /= mask(j, 2)) cycle genl - end do - - ! If det(i) < det(i_gen), it hs already been considered - if(interesting(i) < i_gen) then - fullMatch = .true. - return - end if - - ! Identify the particles - do j=1, N_int - myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1)) - myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2)) - end do - - call bitstring_to_list_in_selection(myMask(1,1), list(1), na, N_int) - call bitstring_to_list_in_selection(myMask(1,2), list(na+1), nb, N_int) - banned(list(1), list(2)) = .true. - end do genl - -end subroutine spot_isinwf - -! --- - -subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, mat_l, mat_r) - - BEGIN_DOC - ! Computes the contributions A(r,s) by - ! comparing the external determinant to all the internal determinants det(i). - ! an applying two particles (r,s) to the mask. - END_DOC - - use bitmasks - implicit none - - integer, intent(in) :: sp, i_gen, N_sel - integer, intent(in) :: interesting(0:N_sel) - integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel) - logical, intent(inout) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2) - double precision, intent(inout) :: mat(N_states, mo_num, mo_num) - double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num) - - integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt - integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2) - integer(bit_kind) :: phasemask(N_int,2) - - - PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc - PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp - - - mat = 0d0 - mat_l = 0d0 - mat_r = 0d0 - - do i = 1, N_int - negMask(i,1) = not(mask(i,1)) - negMask(i,2) = not(mask(i,2)) - end do - - do i = 1, N_sel - if(interesting(i) < 0) then - stop 'prefetch interesting(i) and det(i)' - endif - - mobMask(1,1) = iand(negMask(1,1), det(1,1,i)) - mobMask(1,2) = iand(negMask(1,2), det(1,2,i)) - nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) - - if(nt > 4) cycle - - do j = 2, N_int - mobMask(j,1) = iand(negMask(j,1), det(j,1,i)) - mobMask(j,2) = iand(negMask(j,2), det(j,2,i)) - nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) - enddo - - if(nt > 4) cycle - - if (interesting(i) == i_gen) then - if(sp == 3) then - do k = 1, mo_num - do j = 1, mo_num - banned(j,k,2) = banned(k,j,1) - enddo - enddo - else - do k = 1, mo_num - do l = k+1, mo_num - banned(l,k,1) = banned(k,l,1) - enddo - enddo - endif - endif - - if (interesting(i) >= i_gen) then - - call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) - call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) - - perMask(1,1) = iand(mask(1,1), not(det(1,1,i))) - perMask(1,2) = iand(mask(1,2), not(det(1,2,i))) - do j=2,N_int - perMask(j,1) = iand(mask(j,1), not(det(j,1,i))) - perMask(j,2) = iand(mask(j,2), not(det(j,2,i))) - end do -! call get_d3_h ( det(1,1,i), bannedOrb, banned, mat , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) ) -! call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_r, mat_l, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) & -! , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) ) - - call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int) - call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int) - - call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int) - if(nt == 4) then - call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) - elseif(nt == 3) then - call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) - else - call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) - endif - elseif(nt == 4) then - call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) - call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) - call past_d2(banned, p, sp) - elseif(nt == 3) then - call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) - call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) - call past_d1(bannedOrb, p) - endif - enddo - -end subroutine splash_pq - -! --- - subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, mat_l, mat_r) - use bitmasks use selection_types implicit none @@ -740,7 +555,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs - do jstate = 1, N_states do istate = 1, N_states s_weight(istate,jstate) = dsqrt(selection_weight(istate)*selection_weight(jstate)) @@ -780,11 +594,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d ! endif !endif - ! MANU: ERREUR dans les calculs puisque < I | H | J > = 0 - ! n'implique pas < I | H_TC | J > = 0 ?? - !val = maxval(abs(mat(1:N_states, p1, p2))) - !if( val == 0d0) cycle - call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) if(do_only_cas) then @@ -811,7 +620,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if(excitation_max >= 0) then do_cycle = .True. if(excitation_ref == 1) then - call get_excitation_degree(HF_bitmask, det(1,1), degree, N_int) + call get_excitation_degree(HF_bitmask,det(1,1),degree,N_int) do_cycle = do_cycle .and. (degree > excitation_max) elseif(excitation_ref == 2) then do k = 1, N_dominant_dets_of_cfgs @@ -995,12 +804,118 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d enddo ! end do p2 enddo ! end do p1 - end subroutine fill_buffer_double - ! --- +subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, mat_l, mat_r) + + BEGIN_DOC + ! Computes the contributions A(r,s) by + ! comparing the external determinant to all the internal determinants det(i). + ! an applying two particles (r,s) to the mask. + END_DOC + + use bitmasks + implicit none + + integer, intent(in) :: sp, i_gen, N_sel + integer, intent(in) :: interesting(0:N_sel) + integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel) + logical, intent(inout) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2) + double precision, intent(inout) :: mat(N_states, mo_num, mo_num) + double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num) + + integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt + integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2) + integer(bit_kind) :: phasemask(N_int,2) + PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc + PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp + + + mat = 0d0 + mat_l = 0d0 + mat_r = 0d0 + + do i = 1, N_int + negMask(i,1) = not(mask(i,1)) + negMask(i,2) = not(mask(i,2)) + end do + + do i = 1, N_sel + if(interesting(i) < 0) then + stop 'prefetch interesting(i) and det(i)' + endif + + mobMask(1,1) = iand(negMask(1,1), det(1,1,i)) + mobMask(1,2) = iand(negMask(1,2), det(1,2,i)) + nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) + + if(nt > 4) cycle + + do j = 2, N_int + mobMask(j,1) = iand(negMask(j,1), det(j,1,i)) + mobMask(j,2) = iand(negMask(j,2), det(j,2,i)) + nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) + enddo + + if(nt > 4) cycle + + if (interesting(i) == i_gen) then + if(sp == 3) then + do k = 1, mo_num + do j = 1, mo_num + banned(j,k,2) = banned(k,j,1) + enddo + enddo + else + do k = 1, mo_num + do l = k+1, mo_num + banned(l,k,1) = banned(k,l,1) + enddo + enddo + endif + endif + + if (interesting(i) >= i_gen) then + + call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) + call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) + + perMask(1,1) = iand(mask(1,1), not(det(1,1,i))) + perMask(1,2) = iand(mask(1,2), not(det(1,2,i))) + do j=2,N_int + perMask(j,1) = iand(mask(j,1), not(det(j,1,i))) + perMask(j,2) = iand(mask(j,2), not(det(j,2,i))) + end do +! call get_d3_h ( det(1,1,i), bannedOrb, banned, mat , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) ) +! call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_r, mat_l, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) & +! , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) ) + + call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int) + call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int) + + call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int) + if(nt == 4) then + call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) + elseif(nt == 3) then + call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) + else + call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) + endif + elseif(nt == 4) then + call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) + call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) + call past_d2(banned, p, sp) + elseif(nt == 3) then + call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) + call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) + call past_d1(bannedOrb, p) + endif + enddo + +end subroutine splash_pq +! --- subroutine past_d1(bannedOrb, p) use bitmasks @@ -1043,9 +958,61 @@ subroutine past_d2(banned, p, sp) end do end do end if - end subroutine past_d2 +subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting) + + use bitmasks + implicit none + + BEGIN_DOC + ! Identify the determinants in det which are in the internal space. These are + ! the determinants that can be produced by creating two particles on the mask. + END_DOC + + integer, intent(in) :: i_gen, N + integer, intent(in) :: interesting(0:N) + integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N) + logical, intent(inout) :: banned(mo_num, mo_num) + logical, intent(out) :: fullMatch + + integer :: i, j, na, nb, list(3) + integer(bit_kind) :: myMask(N_int, 2), negMask(N_int, 2) + + fullMatch = .false. + + do i=1,N_int + negMask(i,1) = not(mask(i,1)) + negMask(i,2) = not(mask(i,2)) + end do + + genl : do i=1, N + ! If det(i) can't be generated by the mask, cycle + do j=1, N_int + if(iand(det(j,1,i), mask(j,1)) /= mask(j, 1)) cycle genl + if(iand(det(j,2,i), mask(j,2)) /= mask(j, 2)) cycle genl + end do + + ! If det(i) < det(i_gen), it hs already been considered + if(interesting(i) < i_gen) then + fullMatch = .true. + return + end if + + ! Identify the particles + do j=1, N_int + myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1)) + myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2)) + end do + +! call debug_det(myMask, N_int) + call bitstring_to_list_in_selection(myMask(1,1), list(1), na, N_int) + call bitstring_to_list_in_selection(myMask(1,2), list(na+1), nb, N_int) + banned(list(1), list(2)) = .true. + end do genl + +end subroutine spot_isinwf + ! --- subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint) diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index e67287d3..64e7e6ba 100644 --- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -54,7 +54,7 @@ subroutine run_stochastic_cipsi ! if (N_det > N_det_max) then -! psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted_tc_gen(1:N_int,1:2,1:N_det) +! psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det) ! psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states) ! N_det = N_det_max ! soft_touch N_det psi_det psi_coef diff --git a/src/fci_tc_bi/diagonalize_ci.irp.f b/src/fci_tc_bi/diagonalize_ci.irp.f index 56c561ac..c8369e93 100644 --- a/src/fci_tc_bi/diagonalize_ci.irp.f +++ b/src/fci_tc_bi/diagonalize_ci.irp.f @@ -49,9 +49,8 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) psi_coef(i,j) = dabs(psi_l_coef_bi_ortho(i,j) * psi_r_coef_bi_ortho(i,j)) enddo enddo - SOFT_TOUCH eigval_left_tc_bi_orth eigval_right_tc_bi_orth leigvec_tc_bi_orth reigvec_tc_bi_orth norm_ground_left_right_bi_orth psi_coef psi_l_coef_bi_ortho psi_r_coef_bi_ortho - - + SOFT_TOUCH eigval_left_tc_bi_orth eigval_right_tc_bi_orth leigvec_tc_bi_orth reigvec_tc_bi_orth norm_ground_left_right_bi_orth + SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho psi_coef call save_tc_bi_ortho_wavefunction end diff --git a/src/fci_tc_bi/generators.irp.f b/src/fci_tc_bi/generators.irp.f index 55c0cbb9..bf972423 100644 --- a/src/fci_tc_bi/generators.irp.f +++ b/src/fci_tc_bi/generators.irp.f @@ -43,9 +43,14 @@ END_PROVIDER ! For Single reference wave functions, the generator is the ! Hartree-Fock determinant END_DOC - psi_det_sorted_tc_gen = psi_det_sorted_tc + psi_det_sorted_tc_gen = psi_det_sorted_tc psi_coef_sorted_tc_gen = psi_coef_sorted_tc psi_det_sorted_tc_gen_order = psi_det_sorted_tc_order + integer :: i +! do i = 1,N_det +! print*,'i = ',i +! call debug_det(psi_det_sorted_tc(1,1,i),N_int) +! enddo END_PROVIDER diff --git a/src/fci_tc_bi/selectors.irp.f b/src/fci_tc_bi/selectors.irp.f index 3830927b..94aa4b01 100644 --- a/src/fci_tc_bi/selectors.irp.f +++ b/src/fci_tc_bi/selectors.irp.f @@ -47,13 +47,20 @@ END_PROVIDER enddo do k=1,N_states do i=1,N_det_selectors - psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k) + psi_selectors_coef(i,k) = dsqrt(dabs(psi_l_coef_sorted_bi_ortho(i,k) * psi_r_coef_sorted_bi_ortho(i,k))) psi_selectors_coef_tc(i,1,k) = psi_l_coef_sorted_bi_ortho(i,k) psi_selectors_coef_tc(i,2,k) = psi_r_coef_sorted_bi_ortho(i,k) +! psi_selectors_coef_tc(i,1,k) = psi_l_coef_bi_ortho(i,k) +! psi_selectors_coef_tc(i,2,k) = psi_r_coef_bi_ortho(i,k) ! psi_selectors_coef_tc(i,1,k) = 1.d0 ! psi_selectors_coef_tc(i,2,k) = 1.d0 enddo enddo + print*,'selectors ' + do i = 1, N_det_selectors + print*,i,dabs(psi_selectors_coef_tc(i,1,1)*psi_selectors_coef_tc(i,2,1)) + call debug_det(psi_selectors(1,1,i),N_int) + enddo END_PROVIDER diff --git a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f index 42617557..35c78468 100644 --- a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f +++ b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f @@ -10,7 +10,10 @@ BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_tc, (psi_det_size) ] psi_average_norm_contrib_tc(:) = 0.d0 do k=1,N_states + print*,'in psi_average_norm_contrib_tc' do i=1,N_det + print*,i,dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k)) + call debug_det(psi_det(1,1,i),N_int) psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i) + & dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k))*state_average_weight(k) enddo @@ -26,11 +29,18 @@ END_PROVIDER &BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc, (psi_det_size,N_states) ] &BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted_tc, (psi_det_size) ] &BEGIN_PROVIDER [ integer, psi_det_sorted_tc_order, (psi_det_size) ] +&BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho, (psi_det_size, N_states)] +&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho, (psi_det_size, N_states)] implicit none BEGIN_DOC ! Wave function sorted by determinants contribution to the norm (state-averaged) ! ! psi_det_sorted_tc_order(i) -> k : index in psi_det + ! + ! psi_r_coef_sorted_bi_ortho : right coefficients corresponding to psi_det_sorted_tc + ! + ! psi_l_coef_sorted_bi_ortho : left coefficients corresponding to psi_det_sorted_tc + END_DOC END_DOC integer :: i,j,k integer, allocatable :: iorder(:) @@ -39,7 +49,7 @@ END_PROVIDER psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_tc(i) iorder(i) = i enddo -! call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det) + call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det) do i=1,N_det do j=1,N_int psi_det_sorted_tc(j,1,i) = psi_det(j,1,iorder(i)) @@ -66,24 +76,13 @@ END_PROVIDER psi_average_norm_contrib_sorted_tc(N_det+1:psi_det_size) = 0.d0 psi_det_sorted_tc_order(N_det+1:psi_det_size) = 0 - deallocate(iorder) - -END_PROVIDER - - BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho, (psi_det_size, N_states)] -&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho, (psi_det_size, N_states)] - BEGIN_DOC - ! psi_r_coef_sorted_bi_ortho : right coefficients corresponding to psi_det_sorted_tc - ! psi_l_coef_sorted_bi_ortho : left coefficients corresponding to psi_det_sorted_tc - END_DOC - implicit none - integer :: i, j, k psi_r_coef_sorted_bi_ortho = 0.d0 psi_l_coef_sorted_bi_ortho = 0.d0 do i = 1, N_det - psi_r_coef_sorted_bi_ortho(i,1) = psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1) - psi_l_coef_sorted_bi_ortho(i,1) = psi_l_coef_bi_ortho(psi_det_sorted_tc_order(i),1) + psi_r_coef_sorted_bi_ortho(i,1:N_states) = psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1:N_states) + psi_l_coef_sorted_bi_ortho(i,1:N_states) = psi_l_coef_bi_ortho(psi_det_sorted_tc_order(i),1:N_states) enddo + deallocate(iorder) END_PROVIDER diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f index d39b7a29..d12bbb4e 100644 --- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f +++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f @@ -133,7 +133,10 @@ call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,size(reigvec_tc_bi_orth,1),N_det,N_states) print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) norm_ground_left_right_bi_orth = 0.d0 + print*,'In diago' do j = 1, N_det + print*,j,dabs(leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1)) + call debug_det(psi_det(1,1,j),N_int) norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1) enddo print*,'norm l/r = ',norm_ground_left_right_bi_orth