diff --git a/plugins/dress_zmq/alpha_factory.irp.f b/plugins/dress_zmq/alpha_factory.irp.f index acd2e3ca..cd9df4be 100644 --- a/plugins/dress_zmq/alpha_factory.irp.f +++ b/plugins/dress_zmq/alpha_factory.irp.f @@ -1,16 +1,6 @@ use bitmasks -subroutine assert(cond, msg) - character(*), intent(in) :: msg - logical, intent(in) :: cond - - if(.not. cond) then - print *, "assert failed: "//msg - stop - end if -end subroutine - subroutine get_mask_phase(det, phasemask) use bitmasks @@ -74,6 +64,8 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index integer,allocatable :: preinteresting(:), prefullinteresting(:), interesting(:), fullinteresting(:) integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) logical, allocatable :: banned(:,:,:), bannedOrb(:,:) + integer, allocatable :: counted(:,:), countedOrb(:,:) + integer :: countedGlob double precision, allocatable :: mat(:,:,:) @@ -85,7 +77,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index allocate(abuf(N_int, 2, mo_tot_num**2)) allocate(preinteresting_det(N_int,2,N_det)) - + PROVIDE fragment_count @@ -202,6 +194,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index allocate(minilist(N_int, 2, N_det_selectors), fullminilist(N_int, 2, N_det)) allocate(banned(mo_tot_num, mo_tot_num,2), bannedOrb(mo_tot_num, 2)) + allocate(counted(mo_tot_num, mo_tot_num), countedOrb(mo_tot_num, 2)) allocate (mat(N_states, mo_tot_num, mo_tot_num)) maskInd = -1 integer :: nb_count @@ -360,11 +353,11 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting) if(fullMatch) cycle - !call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) + call count_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, countedGlob, countedOrb, counted, interesting) + call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, countedGlob, countedOrb, counted, interesting) + !call create_alpha_buffer(i_generator, sp, h1, h2, bannedOrb, banned, abuf, n) - call create_alpha_buffer(i_generator, sp, h1, h2, bannedOrb, banned, abuf, n) - - call dress_with_alpha_buffer(delta_ij_loc, minilist, interesting(0), abuf, n) + !call dress_with_alpha_buffer(delta_ij_loc, minilist, interesting(0), abuf, n) end if enddo if(s1 /= s2) monoBdo = .false. @@ -374,220 +367,9 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index end subroutine -subroutine create_alpha_buffer(i_generator, sp, h1, h2, bannedOrb, banned, abuf, n) - use bitmasks - implicit none - - integer, intent(in) :: i_generator, sp, h1, h2 - logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num) - integer(bit_kind),intent(inout) :: abuf(N_int, 2, *) ! mo_tot_num**2 - integer, intent(out) :: n - logical :: ok - integer :: s1, s2, p1, p2, ib, j, istate - integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - - n = 0 - - if(sp == 3) then - s1 = 1 - s2 = 2 - else - s1 = sp - s2 = sp - end if - - call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) - - do p1=1,mo_tot_num - if(bannedOrb(p1, s1)) cycle - ib = 1 - if(sp /= 3) ib = p1+1 - do p2=ib,mo_tot_num - if(bannedOrb(p2, s2)) cycle - if(banned(p1,p2)) cycle - n += 1 - call apply_particles(mask, s1, p1, s2, p2, abuf(1,1,n), ok, N_int) - if(.not. ok) stop "error in create_alpha_buffer" - end do - end do -end -double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2) - use bitmasks - implicit none - - integer, intent(in) :: phasemask(2,*) - integer, intent(in) :: s1, s2, h1, h2, p1, p2 - logical :: change - integer :: np1 - integer :: np - double precision, save :: res(0:1) = (/1d0, -1d0/) - - np1 = phasemask(s1,h1) + phasemask(s1,p1) + phasemask(s2,h2) + phasemask(s2,p2) - np = np1 - if(p1 < h1) np = np + 1 - if(p2 < h2) np = np + 1 - - if(s1 == s2 .and. max(h1, p1) > min(h2, p2)) np = np + 1 - get_phase_bi = res(iand(np,1)) -end - - -subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) - use bitmasks - implicit none - - integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2) - integer, intent(in) :: phasemask(2,N_int*bit_kind_size) - logical, intent(in) :: bannedOrb(mo_tot_num) - double precision, intent(in) :: coefs(N_states) - double precision, intent(inout) :: vect(N_states, mo_tot_num) - integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2) - integer :: i, j, h1, h2, p1, p2, sfix, hfix, pfix, hmob, pmob, puti - double precision :: hij - double precision, external :: get_phase_bi, mo_bielec_integral - - integer, parameter :: turn3_2(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) - integer, parameter :: turn2(2) = (/2,1/) - - if(h(0,sp) == 2) then - h1 = h(1, sp) - h2 = h(2, sp) - do i=1,3 - puti = p(i, sp) - if(bannedOrb(puti)) cycle - p1 = p(turn3_2(1,i), sp) - p2 = p(turn3_2(2,i), sp) - hij = mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2, p1, h1, h2) - hij *= get_phase_bi(phasemask, sp, sp, h1, p1, h2, p2) - vect(:, puti) += hij * coefs - end do - else if(h(0,sp) == 1) then - sfix = turn2(sp) - hfix = h(1,sfix) - pfix = p(1,sfix) - hmob = h(1,sp) - do j=1,2 - puti = p(j, sp) - if(bannedOrb(puti)) cycle - pmob = p(turn2(j), sp) - hij = mo_bielec_integral(pfix, pmob, hfix, hmob) - hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix) - vect(:, puti) += hij * coefs - end do - else - puti = p(1,sp) - if(.not. bannedOrb(puti)) then - sfix = turn2(sp) - p1 = p(1,sfix) - p2 = p(2,sfix) - h1 = h(1,sfix) - h2 = h(2,sfix) - hij = (mo_bielec_integral(p1,p2,h1,h2) - mo_bielec_integral(p2,p1,h1,h2)) - hij *= get_phase_bi(phasemask, sfix, sfix, h1, p1, h2, p2) - vect(:, puti) += hij * coefs - end if - end if -end - - -subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) - use bitmasks - implicit none - - integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2) - integer, intent(in) :: phasemask(2,N_int*bit_kind_size) - logical, intent(in) :: bannedOrb(mo_tot_num) - double precision, intent(in) :: coefs(N_states) - double precision, intent(inout) :: vect(N_states, mo_tot_num) - integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2) - integer :: i, hole, p1, p2, sh - logical :: ok - - logical, allocatable :: lbanned(:) - integer(bit_kind) :: det(N_int, 2) - double precision :: hij - double precision, external :: get_phase_bi, mo_bielec_integral - - allocate (lbanned(mo_tot_num)) - lbanned = bannedOrb - sh = 1 - if(h(0,2) == 1) sh = 2 - hole = h(1, sh) - lbanned(p(1,sp)) = .true. - if(p(0,sp) == 2) lbanned(p(2,sp)) = .true. - !print *, "SPm1", sp, sh - - p1 = p(1, sp) - - if(sp == sh) then - p2 = p(2, sp) - lbanned(p2) = .true. - - do i=1,hole-1 - if(lbanned(i)) cycle - hij = (mo_bielec_integral(p1, p2, i, hole) - mo_bielec_integral(p2, p1, i, hole)) - hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2) - vect(1:N_states,i) += hij * coefs(1:N_states) - end do - do i=hole+1,mo_tot_num - if(lbanned(i)) cycle - hij = (mo_bielec_integral(p1, p2, hole, i) - mo_bielec_integral(p2, p1, hole, i)) - hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2) - vect(1:N_states,i) += hij * coefs(1:N_states) - end do - - call apply_particle(mask, sp, p2, det, ok, N_int) - call i_h_j(gen, det, N_int, hij) - vect(1:N_states, p2) += hij * coefs(1:N_states) - else - p2 = p(1, sh) - do i=1,mo_tot_num - if(lbanned(i)) cycle - hij = mo_bielec_integral(p1, p2, i, hole) - hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2) - vect(1:N_states,i) += hij * coefs(1:N_states) - end do - end if - - call apply_particle(mask, sp, p1, det, ok, N_int) - call i_h_j(gen, det, N_int, hij) - vect(1:N_states, p1) += hij * coefs(1:N_states) -end - - -subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) - use bitmasks - implicit none - - integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2) - integer, intent(in) :: phasemask(2,N_int*bit_kind_size) - logical, intent(in) :: bannedOrb(mo_tot_num) - double precision, intent(in) :: coefs(N_states) - double precision, intent(inout) :: vect(N_states, mo_tot_num) - integer, intent(in) :: sp, h(0:2, 2), p(0:3, 2) - integer :: i - logical :: ok - - logical, allocatable :: lbanned(:) - integer(bit_kind) :: det(N_int, 2) - double precision :: hij - - allocate(lbanned(mo_tot_num)) - lbanned = bannedOrb - lbanned(p(1,sp)) = .true. - do i=1,mo_tot_num - if(lbanned(i)) cycle - call apply_particle(mask, sp, i, det, ok, N_int) - call i_h_j(gen, det, N_int, hij) - vect(1:N_states, i) += hij * coefs(1:N_states) - end do -end - - - -subroutine count_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting) +subroutine count_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, countedGlob, countedOrb, counted, interesting) use bitmasks implicit none @@ -595,14 +377,16 @@ subroutine count_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interes 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_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2) - integer, intent(inout) :: mat(mo_tot_num, mo_tot_num) + integer, intent(inout) :: countedGlob, countedOrb(mo_tot_num, 2), counted(mo_tot_num, mo_tot_num) + - integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt + integer :: i, s, 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 :: phasemask(2,N_int*bit_kind_size) PROVIDE psi_selectors_coef_transp - mat = 0 + countedGlob = 0 + countedOrb = 0 + counted = 0 do i=1,N_int negMask(i,1) = not(mask(i,1)) @@ -631,19 +415,34 @@ subroutine count_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interes if(nt > 4) cycle if (interesting(i) == i_gen) then - if(sp == 3) then - do j=1,mo_tot_num - do k=1,mo_tot_num - banned(j,k,2) = banned(k,j,1) - enddo - enddo - else - do k=1,mo_tot_num - do l=k+1,mo_tot_num - banned(l,k,1) = banned(k,l,1) - end do - end do + do s=1,2 + do j=1,mo_tot_num + if(bannedOrb(j, s)) then + if(sp == 3 .and. s == 1) then + banned(j, :, 1) = .true. + else if(sp == 3 .and. s == 2) then + banned(:, j, 1) = .true. + else if(s == sp) then + banned(j,:,1) = .true. + banned(:,j,1) = .true. + end if end if + end do + end do + + if(sp == 3) then + do j=1,mo_tot_num + do k=1,mo_tot_num + banned(j,k,2) = banned(k,j,1) + enddo + enddo + else + do k=1,mo_tot_num + do l=k+1,mo_tot_num + banned(l,k,1) = banned(k,l,1) + end do + end do + end if end if call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) @@ -660,133 +459,157 @@ subroutine count_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interes call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int) if (interesting(i) >= i_gen) then - call get_mask_phase(psi_det_sorted(1,1,interesting(i)), phasemask) if(nt == 4) then - call count_d2(mat, p, sp) + call count_d2(counted, p, sp) else if(nt == 3) then - call count_d1(mat, p, sp) + call count_d1(countedOrb, p) else - mat(:,:) = mat(:,:) + 1 + countedGlob += 1 end if else if(nt == 4) call past_d2(banned, p, sp) if(nt == 3) call past_d1(bannedOrb, p) + if(nt < 3) stop "past_d0 ?" end if end do + do i=1,mo_tot_num + if(bannedOrb(i,1)) countedOrb(i,1) = 0 + if(bannedOrb(i,2)) countedOrb(i,2) = 0 + do j=1,mo_tot_num + if(banned(i,j,1)) counted(i,j) = 0 + end do + end do + + if(sp /= 3) then + countedOrb(:, mod(sp, 2)+1) = 0 + end if + + ! USELESS? + !do j=1,2 + ! do i=1,mo_tot_num + ! if(bannedOrb(i, j)) then + ! countedOrb(i, j) = 0 + ! if(sp == 3 .and. j == 1) then + ! counted(i, :) = 0 + ! else if(sp == 3 .and. j == 2) then + ! counted(:, i) = 0 + ! else if(j == sp) then + ! counted(i,:) = 0 + ! counted(:,i) = 0 + ! end if + ! end if + ! end do + !end do +end + + + +!subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting) +subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, countedGlob, countedOrb, counted, interesting) + 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_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2) + integer, intent(inout) :: countedGlob, countedOrb(mo_tot_num, 2), counted(mo_tot_num, mo_tot_num) + integer :: counted2(mo_tot_num, mo_tot_num), countedOrb2(mo_tot_num, 2) + integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt, s + integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2) + integer :: phasemask(2,N_int*bit_kind_size) + + counted2 = counted + countedOrb2 = countedOrb + PROVIDE psi_selectors_coef_transp + 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 ! interesting(0) + !i = interesting(ii) + if (interesting(i) < 0) then + stop 'prefetch interesting(i)' + endif + if(interesting(i) < i_gen) cycle + + + 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)) + end do + + if(nt > 4) cycle + + + 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 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) + + if (interesting(i) >= i_gen) then + !call get_mask_phase(psi_det_sorted(1,1,interesting(i)), phasemask) + if(nt == 4) then + call get_d2(det(1,1,i), phasemask, bannedOrb, banned, counted, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) + else if(nt == 3) then + call get_d1(det(1,1,i), phasemask, bannedOrb, banned, counted, countedOrb, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) + else + countedGlob -= 1 + !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) call past_d2(banned, p, sp) + ! if(nt == 3) call past_d1(bannedOrb, p) + end if + end do + + if(countedGlob /= 0) stop "nonul glob" + + do s=1,2 + do i=1,mo_tot_num + if(countedOrb(i, s) /= 0) then + print *, "COUNTEDORB", sp, s, bannedOrb(i,s), countedOrb2(i, s), countedOrb(i, s) + !stop "COUNERe" + end if + end do + end do + do i=1,mo_tot_num do j=1,mo_tot_num - if(banned(i,j,1)) mat(i,j) = 0 - end do - end do - - if(sp == 3) then - do i=1,mo_tot_num - if(bannedOrb(i, 1)) mat(i, :) = 0 - if(bannedOrb(i, 2)) mat(:, i) = 0 - end do - else - do i=1,mo_tot_num - if(bannedOrb(i, sp)) then - mat(:,i) = 0 - mat(i,:) = 0 + if(counted2(i,j) /= 0) then + if(counted(i,j) /= 0) then + print *, counted(i,j) + stop "nonul" + end if + else + if(counted(i,j) /= 0) then + print *, counted2(i,j), counted(i,j) + stop "UNCOUNTED" end if - end do - end if -end - - - -subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting) - 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_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2) - double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_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 :: phasemask(2,N_int*bit_kind_size) - - PROVIDE psi_selectors_coef_transp - mat = 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 ! interesting(0) - !i = interesting(ii) - if (interesting(i) < 0) then - stop 'prefetch interesting(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)) - end do - - if(nt > 4) cycle - - if (interesting(i) == i_gen) then - if(sp == 3) then - do j=1,mo_tot_num - do k=1,mo_tot_num - banned(j,k,2) = banned(k,j,1) - enddo - enddo - else - do k=1,mo_tot_num - do l=k+1,mo_tot_num - banned(l,k,1) = banned(k,l,1) - end do - end do - end if - end if - - 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 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) - - if (interesting(i) >= i_gen) then - call get_mask_phase(psi_det_sorted(1,1,interesting(i)), phasemask) - if(nt == 4) then - 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(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) - else - 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) call past_d2(banned, p, sp) - if(nt == 3) call past_d1(bannedOrb, p) end if end do -end + end do +end subroutine -subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) +subroutine get_d2(gen, phasemask, bannedOrb, banned, counted, mask, h, p, sp, coefs) use bitmasks implicit none @@ -794,7 +617,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer, intent(in) :: phasemask(2,N_int*bit_kind_size) logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2) double precision, intent(in) :: coefs(N_states) - double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) + integer, intent(inout) :: counted(mo_tot_num, mo_tot_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp double precision, external :: get_phase_bi, mo_bielec_integral @@ -834,9 +657,9 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) if(ma == 1) then - mat(:, putj, puti) += coefs(:) * hij + counted(putj, puti) -= 1!coefs(:) * hij else - mat(:, puti, putj) += coefs(:) * hij + counted(puti, putj) -= 1!coefs(:) * hij end if end do else @@ -852,7 +675,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) p1 = p(turn2(i), 1) hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2) - mat(:, puti, putj) += coefs(:) * hij + counted(puti, putj) -= 1!coefs(:) * hij end do end do end if @@ -872,7 +695,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) p1 = p(i1, ma) p2 = p(i2, ma) hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2) - mat(:, puti, putj) += coefs(:) * hij + counted(puti, putj) -= 1!coefs(:) * hij end do end do else if(tip == 3) then @@ -886,7 +709,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) p2 = p(i, ma) hij = mo_bielec_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2) - mat(:, min(puti, putj), max(puti, putj)) += coefs(:) * hij + counted(min(puti, putj), max(puti, putj)) -= 1!coefs(:) * hij end do else ! tip == 4 puti = p(1, sp) @@ -897,14 +720,14 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) h1 = h(1, mi) h2 = h(2, mi) hij = (mo_bielec_integral(p1, p2, h1, h2) - mo_bielec_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2) - mat(:, puti, putj) += coefs(:) * hij + counted(puti, putj) -= 1!coefs(:) * hij end if end if end if end -subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) +subroutine get_d1(gen, phasemask, bannedOrb, banned, counted, countedOrb, mask, h, p, sp, coefs) use bitmasks implicit none @@ -913,7 +736,8 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2) integer(bit_kind) :: det(N_int, 2) double precision, intent(in) :: coefs(N_states) - double precision, intent(inout) :: mat(N_states, mo_tot_num, mo_tot_num) + integer, intent(inout) :: counted(mo_tot_num, mo_tot_num) + integer, intent(inout) :: countedOrb(mo_tot_num, 2) integer, intent(in) :: h(0:2,2), p(0:4,2), sp double precision :: hij, tmp_row(N_states, mo_tot_num), tmp_row2(N_states, mo_tot_num) double precision, external :: get_phase_bi, mo_bielec_integral @@ -966,9 +790,11 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) end do if(ma == 1) then - mat(1:N_states,1:mo_tot_num,puti) += tmp_row(1:N_states,1:mo_tot_num) + !mat(1:N_states,1:mo_tot_num,puti) += tmp_row(1:N_states,1:mo_tot_num) + countedOrb(puti, 2) -= 1 else - mat(1:N_states,puti,1:mo_tot_num) += tmp_row(1:N_states,1:mo_tot_num) + !mat(1:N_states,puti,1:mo_tot_num) += tmp_row(1:N_states,1:mo_tot_num) + countedOrb(puti, 1) -= 1 end if end if @@ -993,11 +819,15 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) end do if(mi == 1) then - mat(:,:,p1) += tmp_row(:,:) - mat(:,:,p2) += tmp_row2(:,:) + !mat(:,:,p1) += tmp_row(:,:) + if(.not. bannedOrb(p1, 2)) countedOrb(p1, 2) = countedOrb(p1,2) - 1 + !mat(:,:,p2) += tmp_row2(:,:) + if(.not. bannedOrb(p2, 2)) countedOrb(p2, 2) = countedOrb(p2,2)-1 else - mat(:,p1,:) += tmp_row(:,:) - mat(:,p2,:) += tmp_row2(:,:) + !mat(:,p1,:) += tmp_row(:,:) + if(.not. bannedOrb(p1, 1)) countedOrb(p1, 1) = countedOrb(p1,1)-1 + !mat(:,p2,:) += tmp_row2(:,:) + if(.not. bannedOrb(p2, 1)) countedOrb(p2, 1) = countedOrb(p2,1)-1 end if else if(p(0,ma) == 3) then @@ -1018,8 +848,11 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) tmp_row(:,putj) += hij * coefs(:) end do - mat(:, :puti-1, puti) += tmp_row(:,:puti-1) - mat(:, puti, puti:) += tmp_row(:,puti:) + !mat(:, :puti-1, puti) += tmp_row(:,:puti-1) + !mat(:, puti, puti:) += tmp_row(:,puti:) + if(.not. bannedOrb(puti, sp)) then + countedOrb(puti, sp) -= 1 + end if end do else hfix = h(1,mi) @@ -1042,10 +875,16 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) tmp_row2(:,puti) += hij * coefs(:) end if end do - mat(:,:p2-1,p2) += tmp_row(:,:p2-1) - mat(:,p2,p2:) += tmp_row(:,p2:) - mat(:,:p1-1,p1) += tmp_row2(:,:p1-1) - mat(:,p1,p1:) += tmp_row2(:,p1:) + !mat(:,:p2-1,p2) += tmp_row(:,:p2-1) + !mat(:,p2,p2:) += tmp_row(:,p2:) + if(.not. bannedOrb(p2, sp)) then + countedOrb(p2, sp) -= 1 + end if + !mat(:,:p1-1,p1) += tmp_row2(:,:p1-1) + !mat(:,p1,p1:) += tmp_row2(:,p1:) + if(.not. bannedOrb(p1, sp)) then + countedOrb(p1, sp) -= 1 + end if end if end if @@ -1067,7 +906,8 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) call i_h_j(gen, det, N_int, hij) - mat(:, p1, p2) += coefs(:) * hij + !mat(:, p1, p2) += coefs(:) * hij + !!!!!!!! DUPLICTATE counted(p1, p2) !!!!!!!!!!!!!!!!!!!! end do end do end @@ -1149,65 +989,6 @@ subroutine past_d1(bannedOrb, p) end -subroutine count_d1(mat, p, sp) - use bitmasks - implicit none - - integer, intent(inout) :: mat(mo_tot_num, mo_tot_num) - integer, intent(in) :: p(0:4, 2), sp - integer :: i,s,j - - - if(sp == 3) then - do i=1,p(0,1) - mat(p(i,1), :) += 1 - end do - do i=1,p(0,2) - mat(:, p(i,2)) += 1 - end do - - do i=1,p(0,1) - do j=1,p(0,2) - mat(p(i,1), p(j,2)) -= 1 - end do - end do - else - if(sp == 1 .and. p(0,2) /= 0) stop "count_d1 bug" - if(sp == 2 .and. p(0,1) /= 0) stop "count_d1 bug" - do i=1,p(0,sp) - mat(:p(i,sp), p(i,sp)) += 1 - mat(p(i,sp), p(i,sp):) += 1 - mat(p(i,sp), p(i,sp)) -= 1 - end do - end if -end - - -subroutine count_d2(mat, p, sp) - use bitmasks - implicit none - - integer, intent(inout) :: mat(mo_tot_num, mo_tot_num) - integer, intent(in) :: p(0:4, 2), sp - integer :: i,j - - if(sp == 3) then - do i=1,p(0,1) - do j=1,p(0,2) - mat(p(i,1), p(j,2)) += 1 - end do - end do - else - do i=1,p(0, sp) - do j=1,i-1 - mat(p(j,sp), p(i,sp)) += 1 - mat(p(i,sp), p(j,sp)) += 1 - end do - end do - end if -end - - subroutine past_d2(banned, p, sp) use bitmasks implicit none @@ -1225,6 +1006,7 @@ subroutine past_d2(banned, p, sp) else do i=1,p(0, sp) do j=1,i-1 + if(p(j,sp) > p(i,sp)) stop "PPPPPPPP" banned(p(j,sp), p(i,sp)) = .true. banned(p(i,sp), p(j,sp)) = .true. end do @@ -1233,6 +1015,48 @@ subroutine past_d2(banned, p, sp) end +subroutine count_d1(countedOrb, p) + use bitmasks + implicit none + + integer, intent(inout) :: countedOrb(mo_tot_num, 2) + integer, intent(in) :: p(0:4, 2) + integer :: i,s + + do s = 1, 2 + do i = 1, p(0, s) + countedOrb(p(i, s), s) += 1 + end do + end do +end + + +subroutine count_d2(counted, p, sp) + use bitmasks + implicit none + + integer, intent(inout) :: counted(mo_tot_num, mo_tot_num) + integer, intent(in) :: p(0:4, 2), sp + integer :: i,j + + if(sp == 3) then + do i=1,p(0,1) + do j=1,p(0,2) + counted(p(i,1), p(j,2)) += 1 + end do + end do + else + do i=1,p(0, sp) + do j=1,i-1 + counted(p(j,sp), p(i,sp)) += 1 + !counted(p(i,sp), p(j,sp)) += 1 + end do + end do + end if +end + + + subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting) use bitmasks implicit none @@ -1304,3 +1128,65 @@ subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint) enddo end + + + +subroutine create_alpha_buffer(i_generator, sp, h1, h2, bannedOrb, banned, abuf, n) + use bitmasks + implicit none + + integer, intent(in) :: i_generator, sp, h1, h2 + logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num) + integer(bit_kind),intent(inout) :: abuf(N_int, 2, *) ! mo_tot_num**2 + integer, intent(out) :: n + logical :: ok + integer :: s1, s2, p1, p2, ib, j, istate + integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) + + n = 0 + + if(sp == 3) then + s1 = 1 + s2 = 2 + else + s1 = sp + s2 = sp + end if + + call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) + + do p1=1,mo_tot_num + if(bannedOrb(p1, s1)) cycle + ib = 1 + if(sp /= 3) ib = p1+1 + do p2=ib,mo_tot_num + if(bannedOrb(p2, s2)) cycle + if(banned(p1,p2)) cycle + n += 1 + call apply_particles(mask, s1, p1, s2, p2, abuf(1,1,n), ok, N_int) + if(.not. ok) stop "error in create_alpha_buffer" + end do + end do +end + + +double precision function get_phase_bi(phasemask, s1, s2, h1, p1, h2, p2) + use bitmasks + implicit none + + integer, intent(in) :: phasemask(2,*) + integer, intent(in) :: s1, s2, h1, h2, p1, p2 + logical :: change + integer :: np1 + integer :: np + double precision, save :: res(0:1) = (/1d0, -1d0/) + + np1 = phasemask(s1,h1) + phasemask(s1,p1) + phasemask(s2,h2) + phasemask(s2,p2) + np = np1 + if(p1 < h1) np = np + 1 + if(p2 < h2) np = np + 1 + + if(s1 == s2 .and. max(h1, p1) > min(h2, p2)) np = np + 1 + get_phase_bi = res(iand(np,1)) +end + diff --git a/plugins/mrcc_sto/mrcc_sto.irp.f b/plugins/mrcc_sto/mrcc_sto.irp.f index 205c480b..173d8d26 100644 --- a/plugins/mrcc_sto/mrcc_sto.irp.f +++ b/plugins/mrcc_sto/mrcc_sto.irp.f @@ -3,12 +3,12 @@ program mrcc_sto BEGIN_DOC ! TODO END_DOC - print *, "!!!!!!========================!!!!!!" - print *, "!!!!!!========================!!!!!!" - print *, "!!!!!!========================!!!!!!" + print *, "========================" + print *, "========================" + print *, "========================" print *, "MRCC_STO not implemented - acts as a unittest for dress_zmq" - print *, "!!!!!!========================!!!!!!" - print *, "!!!!!!========================!!!!!!" - print *, "!!!!!!========================!!!!!!" + print *, "========================" + print *, "========================" + print *, "========================" call dress_zmq() end