LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 03:23:29 +01:00
Code Releases Activity

Removed += in selection

Browse Source
This commit is contained in:
Anthony Scemama 2021-01-21 10:00:53 +01:00
parent 362dd42ddc
commit 340f7076c1
2 changed files with 149 additions and 86 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
src/cipsi/selection.irp.f
View File

@ -737,12 +737,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
alpha_h_psi = mat(istate, p1, p2)
do jstate=1,N_states
pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate)
enddo
pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi
pt2_data % pt2(istate) += e_pert(istate)
pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate)
pt2_data % variance(istate) = pt2_data % variance(istate) + alpha_h_psi * alpha_h_psi
pt2_data % pt2(istate) = pt2_data % pt2(istate) + e_pert(istate)
!!!DEBUG
! delta_E = E0(istate) - Hii + E_shift
@ -1573,7 +1570,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = mo_two_e_integral(p1, p2, h1, h2) * phase
end if
mat(:, p1, p2) += coefs(:) * hij
mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij
end do
end do
else ! AA BB
@ -1590,7 +1587,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
else
hij = (mo_two_e_integral(p1, p2, puti, putj) - mo_two_e_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
end if
mat(:, puti, putj) += coefs(:) * hij
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do
end do
end if
@ -1649,18 +1646,18 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
do putj=1, hfix-1
if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
end do
do putj=hfix+1, mo_num
if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
end do
if(ma == 1) then
mat(1:N_states,1:mo_num,puti) += tmp_row(1:N_states,1:mo_num)
mat(1:N_states,1:mo_num,puti) = mat(1:N_states,1:mo_num,puti) + tmp_row(1:N_states,1:mo_num)
else
mat(1:N_states,puti,1:mo_num) += tmp_row(1:N_states,1:mo_num)
mat(1:N_states,puti,1:mo_num) = mat(1:N_states,puti,1:mo_num) + tmp_row(1:N_states,1:mo_num)
end if
end if
@ -1674,22 +1671,22 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
putj = p1
if(.not. banned(putj,puti,bant)) then
hij = mo_two_e_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
tmp_row(:,puti) += hij * coefs(:)
tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:)
end if
putj = p2
if(.not. banned(putj,puti,bant)) then
hij = mo_two_e_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
tmp_row2(:,puti) += hij * coefs(:)
tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:)
end if
end do
if(mi == 1) then
mat(:,:,p1) += tmp_row(:,:)
mat(:,:,p2) += tmp_row2(:,:)
mat(:,:,p1) = mat(:,:,p1) + tmp_row(:,:)
mat(:,:,p2) = mat(:,:,p2) + tmp_row2(:,:)
else
mat(:,p1,:) += tmp_row(:,:)
mat(:,p2,:) += tmp_row2(:,:)
mat(:,p1,:) = mat(:,p1,:) + tmp_row(:,:)
mat(:,p2,:) = mat(:,p2,:) + tmp_row2(:,:)
end if
else
if(p(0,ma) == 3) then
@ -1702,16 +1699,16 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
do putj=1,hfix-1
if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
tmp_row(:,putj) += hij * coefs(:)
tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
end do
do putj=hfix+1,mo_num
if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
tmp_row(:,putj) += hij * coefs(:)
tmp_row(:,putj) = 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) = mat(:, :puti-1, puti) + tmp_row(:,:puti-1)
mat(:, puti, puti:) = mat(:, puti, puti:) + tmp_row(:,puti:)
end do
else
hfix = h(1,mi)
@ -1725,19 +1722,19 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
putj = p2
if(.not. banned(puti,putj,1)) then
hij = mo_two_e_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
tmp_row(:,puti) += hij * coefs(:)
tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:)
end if
putj = p1
if(.not. banned(puti,putj,1)) then
hij = mo_two_e_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
tmp_row2(:,puti) += hij * coefs(:)
tmp_row2(:,puti) = 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) = mat(:,:p2-1,p2) + tmp_row(:,:p2-1)
mat(:,p2,p2:) = mat(:,p2,p2:) + tmp_row(:,p2:)
mat(:,:p1-1,p1) = mat(:,:p1-1,p1) + tmp_row2(:,:p1-1)
mat(:,p1,p1:) = mat(:,p1,p1:) + tmp_row2(:,p1:)
end if
end if
deallocate(lbanned)
@ -1760,7 +1757,7 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
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) = mat(:, p1, p2) + coefs(:) * hij
end do
end do
end
@ -1813,9 +1810,9 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
if(ma == 1) then
mat(:, putj, puti) += coefs(:) * hij
mat(:, putj, puti) = mat(:, putj, puti) + coefs(:) * hij
else
mat(:, puti, putj) += coefs(:) * hij
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end if
end do
else
@ -1831,7 +1828,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
p1 = p(turn2(i), 1)
hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int)
mat(:, puti, putj) += coefs(:) * hij
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do
end do
end if
@ -1851,7 +1848,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
p1 = p(i1, ma)
p2 = p(i2, ma)
hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2,N_int)
mat(:, puti, putj) += coefs(:) * hij
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do
end do
else if(tip == 3) then
@ -1865,7 +1862,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
p2 = p(i, ma)
hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int)
mat(:, min(puti, putj), max(puti, putj)) += coefs(:) * hij
mat(:, min(puti, putj), max(puti, putj)) = mat(:, min(puti, putj), max(puti, putj)) + coefs(:) * hij
end do
else ! tip == 4
puti = p(1, sp)
@ -1876,7 +1873,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
h1 = h(1, mi)
h2 = h(2, mi)
hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2,N_int)
mat(:, puti, putj) += coefs(:) * hij
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end if
end if
end if

168
src/davidson/u0_h_u0.irp.f
View File

@ -211,6 +211,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
double precision :: rss, mem, ratio
double precision, allocatable :: utl(:,:)
integer, parameter :: block_size=128
logical :: u_is_sparse
! call resident_memory(rss)
! mem = dble(singles_beta_csc_size) / 1024.d0**3
@ -222,6 +223,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
! endif
compute_singles=.True.
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
@ -249,7 +251,7 @@ compute_singles=.True.
!$OMP singles_alpha_csc,singles_alpha_csc_idx, &
!$OMP singles_beta_csc,singles_beta_csc_idx) &
!$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, &
!$OMP lcol, lrow, l_a, l_b, utl, kk, &
!$OMP lcol, lrow, l_a, l_b, utl, kk, u_is_sparse, &
!$OMP buffer, doubles, n_doubles, umax, &
!$OMP tmp_det2, hij, sij, idx, l, kcol_prev, &
!$OMP singles_a, n_singles_a, singles_b, ratio, &
@ -266,6 +268,19 @@ compute_singles=.True.
kcol_prev=-1
! Check if u has multiple zeros
kk=0
do k=1,N_det
umax = 0.d0
do l=1,N_st
umax = max(umax, dabs(u_t(l,k)))
enddo
if (umax < 1.d-20) then
kk = kk+1
endif
enddo
u_is_sparse = N_det / kk < 10
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
@ -405,16 +420,25 @@ compute_singles=.True.
do k = 1,n_singles_a,block_size
umax = 0.d0
! Prefetch u_t(:,l_a)
do kk=0,block_size-1
if (k+kk > n_singles_a) exit
l_a = singles_a(k+kk)
ASSERT (l_a <= N_det)
if (u_is_sparse) then
do kk=0,block_size-1
if (k+kk > n_singles_a) exit
l_a = singles_a(k+kk)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
enddo
else
do kk=0,block_size-1
if (k+kk > n_singles_a) exit
l_a = singles_a(k+kk)
ASSERT (l_a <= N_det)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
@ -497,16 +521,25 @@ compute_singles=.True.
do i=1,n_singles_a,block_size
umax = 0.d0
! Prefetch u_t(:,l_a)
do kk=0,block_size-1
if (i+kk > n_singles_a) exit
l_a = singles_a(i+kk)
ASSERT (l_a <= N_det)
if (u_is_sparse) then
do kk=0,block_size-1
if (i+kk > n_singles_a) exit
l_a = singles_a(i+kk)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
enddo
else
do kk=0,block_size-1
if (i+kk > n_singles_a) exit
l_a = singles_a(i+kk)
ASSERT (l_a <= N_det)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
@ -534,16 +567,25 @@ compute_singles=.True.
do i=1,n_doubles,block_size
umax = 0.d0
! Prefetch u_t(:,l_a)
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_a = doubles(i+kk)
ASSERT (l_a <= N_det)
if (u_is_sparse) then
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_a = doubles(i+kk)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
enddo
else
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_a = doubles(i+kk)
ASSERT (l_a <= N_det)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
@ -611,19 +653,29 @@ compute_singles=.True.
!DIR$ LOOP COUNT avg(1000)
do i=1,n_singles_b,block_size
umax = 0.d0
do kk=0,block_size-1
if (i+kk > n_singles_b) exit
l_b = singles_b(i+kk)
ASSERT (l_b <= N_det)
if (u_is_sparse) then
do kk=0,block_size-1
if (i+kk > n_singles_b) exit
l_b = singles_b(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
enddo
else
do kk=0,block_size-1
if (i+kk > n_singles_b) exit
l_b = singles_b(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
@ -649,18 +701,28 @@ compute_singles=.True.
!DIR$ LOOP COUNT avg(50000)
do i=1,n_doubles,block_size
umax = 0.d0
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_b = doubles(i+kk)
ASSERT (l_b <= N_det)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
if (u_is_sparse) then
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_b = doubles(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
do l=1,N_st
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo
enddo
else
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_b = doubles(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle
do kk=0,block_size-1
@ -688,10 +750,14 @@ compute_singles=.True.
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
umax = 0.d0
do l=1,N_st
umax = max(umax, dabs(u_t(l,k_a)))
enddo
if (u_is_sparse) then
umax = 0.d0
do l=1,N_st
umax = max(umax, dabs(u_t(l,k_a)))
enddo
else
umax = 1.d0
endif
if (umax < 1.d-20) cycle
krow = psi_bilinear_matrix_rows(k_a)