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

now new version is default, but added the possibility to compute the pt2 the old way and to check the new version

Browse Source
This commit is contained in:
eginer 2023-01-30 18:40:49 +01:00
parent 57a2ea1f67
commit 33a6b33a76
2 changed files with 46 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

81
src/cipsi_tc_bi_ortho/selection.irp.f
View File

@ -772,56 +772,57 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
! -------------------------------------------
istate = 1
call htilde_mu_mat_bi_ortho_tot(det, det, N_int, Hii)
! call htilde_mu_mat_bi_ortho_tot(det, det, N_int, Hii)
double precision :: hmono, htwoe, hthree
call diag_htilde_mu_mat_fock_bi_ortho(N_int, det, hmono, htwoe, hthree, hii)
delta_E = E0(istate) - Hii + E_shift
!delta_E = 1.d0
call get_excitation_degree( HF_bitmask, det, degree, N_int)
double precision :: alpha_h_psi_tmp, psi_h_alpha_tmp
psi_h_alpha_tmp = mat_l(istate, p1, p2)
alpha_h_psi_tmp = mat_r(istate, p1, p2)
!
psi_h_alpha = 0.d0
alpha_h_psi = 0.d0
if(debug_tc_pt2 == 2)then !! Using the old version
psi_h_alpha = 0.d0
alpha_h_psi = 0.d0
do iii = 1, N_det
! call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha)
! call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i)
call htilde_mu_mat_bi_ortho_tot(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
call htilde_mu_mat_bi_ortho_tot(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function
alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function
! psi_h_alpha += i_h_alpha * 1.d0
! alpha_h_psi += alpha_h_i * 1.d0
enddo
!!! print*,'---',p1,p2
!!! call debug_det(det,N_int)
!!! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi
!!! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp
!!! if(dabs(psi_h_alpha_tmp*alpha_h_psi_tmp).gt.1.d+10)then
else if(debug_tc_pt2 == 1)then !! debugging the new version
psi_h_alpha_tmp = mat_l(istate, p1, p2) ! new version
alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version
psi_h_alpha = 0.d0
alpha_h_psi = 0.d0
do iii = 1, N_det ! old version
call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function
alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function
enddo
if(dabs(psi_h_alpha).gt.1.d-10.or.dabs(alpha_h_psi).gt.1.d-10)then
if(dabs(psi_h_alpha - psi_h_alpha_tmp).gt.1.d-10 .or. dabs(alpha_h_psi - alpha_h_psi_tmp).gt.1.d-10)then
! if(dabs(psi_h_alpha_tmp).gt.1.d-10.or.dabs(alpha_h_psi_tmp).gt.1.d-10)then
! if(degree==2)then
call debug_det(det,N_int)
print*,'psi_h_alpha,alpha_h_psi'
print*,psi_h_alpha,alpha_h_psi
print*,psi_h_alpha_tmp,alpha_h_psi_tmp
print*,dabs(psi_h_alpha - psi_h_alpha_tmp),dabs(alpha_h_psi - alpha_h_psi_tmp)
do iii = 1, N_det
call get_excitation_degree( psi_det(1,1,iii), det, degree, N_int)
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha)
call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i)
if(dabs(i_h_alpha).gt.1.d-10.or.dabs(alpha_h_i).gt.1.d-10)then
call debug_det(psi_det(1,1,iii),N_int)
print*,iii,degree,i_h_alpha,alpha_h_i
print*,leigvec_tc_bi_orth(iii,1),reigvec_tc_bi_orth(iii,1)
endif
enddo
stop
endif
endif
if(dabs(psi_h_alpha - psi_h_alpha_tmp).gt.1.d-10 .or. dabs(alpha_h_psi - alpha_h_psi_tmp).gt.1.d-10)then
call debug_det(det,N_int)
print*,'psi_h_alpha,alpha_h_psi'
print*,psi_h_alpha,alpha_h_psi
print*,psi_h_alpha_tmp,alpha_h_psi_tmp
print*,dabs(psi_h_alpha - psi_h_alpha_tmp),dabs(alpha_h_psi - alpha_h_psi_tmp)
do iii = 1, N_det
call get_excitation_degree( psi_det(1,1,iii), det, degree, N_int)
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha)
call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i)
if(dabs(i_h_alpha).gt.1.d-10.or.dabs(alpha_h_i).gt.1.d-10)then
call debug_det(psi_det(1,1,iii),N_int)
print*,iii,degree,i_h_alpha,alpha_h_i
print*,leigvec_tc_bi_orth(iii,1),reigvec_tc_bi_orth(iii,1)
endif
enddo
stop
endif
endif
else
psi_h_alpha = mat_l(istate, p1, p2)
alpha_h_psi = mat_r(istate, p1, p2)
endif
!if(alpha_h_psi*psi_h_alpha/delta_E.gt.1.d-10)then
! print*, 'E0,Hii,E_shift'

5
src/tc_keywords/EZFIO.cfg
View File

@ -184,3 +184,8 @@ doc: Thresholds on the Imag part of energy
interface: ezfio,provider,ocaml
default: 1.e-7
[debug_tc_pt2]
type: integer
doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body
interface: ezfio,provider,ocaml
default: -1