LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-18 11:23:38 +01:00
Code Releases Activity

Merge branch 'dev' of github.com:QuantumPackage/qp2 into dev

Browse Source
This commit is contained in:
Anthony Scemama 2020-05-15 15:58:10 +02:00
commit 764aed423e
1 changed files with 14 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -29,7 +29,7 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
double precision, intent(in) :: variance(N_st) double precision, intent(in) :: variance(N_st)
double precision, intent(in) :: norm(N_st) double precision, intent(in) :: norm(N_st)
double precision :: avg, pt2_rpt2(N_st), element, dt, x double precision :: avg, rpt2(N_st), element, dt, x
integer :: k integer :: k
integer, save :: i_iter=0 integer, save :: i_iter=0
integer, parameter :: i_itermax = 1 integer, parameter :: i_itermax = 1
@ -46,16 +46,16 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
i_iter = 1 i_iter = 1
endif endif
dt = 1.0d0 dt = 2.0d0
do k=1,N_st do k=1,N_st
! PT2 + rPT2 ! rPT2
pt2_rpt2(k) = pt2(k)* (1.d0 + 1.d0/(1.d0 + norm(k))) rpt2(k) = pt2(k)/(1.d0 + norm(k))
enddo enddo
avg = sum(pt2_rpt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero avg = sum(pt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero
do k=1,N_st do k=1,N_st
element = exp(dt*(pt2_rpt2(k)/avg -1.d0)) element = exp(dt*(pt2(k)/avg -1.d0))
element = min(2.0d0 , element) element = min(2.0d0 , element)
element = max(0.5d0 , element) element = max(0.5d0 , element)
memo_pt2(k,i_iter) = element memo_pt2(k,i_iter) = element
@ -71,6 +71,12 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
variance_match_weight(k) *= product(memo_variance(k,:)) variance_match_weight(k) *= product(memo_variance(k,:))
enddo enddo
if (N_det < 100) then
! For tiny wave functions, weights are 1.d0
pt2_match_weight(:) = 1.d0
variance_match_weight(:) = 1.d0
endif
threshold_davidson_pt2 = min(1.d-6, & threshold_davidson_pt2 = min(1.d-6, &
max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) ) max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) )
@ -326,7 +332,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
i = psi_bilinear_matrix_rows(l_a) i = psi_bilinear_matrix_rows(l_a)
if (nt + exc_degree(i) <= 4) then if (nt + exc_degree(i) <= 4) then
idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a)) idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a))
if (psi_average_norm_contrib_sorted(idx) > 0.d0) then if (psi_average_norm_contrib_sorted(idx) > 1.d-20) then
indices(k) = idx indices(k) = idx
k=k+1 k=k+1
endif endif
@ -350,7 +356,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
idx = psi_det_sorted_order( & idx = psi_det_sorted_order( &
psi_bilinear_matrix_order( & psi_bilinear_matrix_order( &
psi_bilinear_matrix_transp_order(l_a))) psi_bilinear_matrix_transp_order(l_a)))
if (psi_average_norm_contrib_sorted(idx) > 0.d0) then if (psi_average_norm_contrib_sorted(idx) > 1.d-20) then
indices(k) = idx indices(k) = idx
k=k+1 k=k+1
endif endif