PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-15 10:33:50 +01:00
Code Releases Activity

tc-two-rdm broken ...

Browse Source
This commit is contained in:
eginer 2023-09-13 12:58:26 +02:00
parent 2ff59e6869
commit 3f95bf40ed
3 changed files with 57 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/tc_bi_ortho/tc_natorb.irp.f
View File

@ -32,19 +32,17 @@
thr_d = 1.d-6
thr_nd = 1.d-6
thr_deg = 1.d-3
if(n_core_orb.ne.0)then
! print*,'core orbitals'
! pause
call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg &
, natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
else
call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg &
, natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
endif
! call non_hrmt_bieig( mo_num, dm_tmp&
! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo&
! , mo_num, natorb_tc_eigval )
! if(n_core_orb.ne.0)then
! call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg &
! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
! else
! call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg &
! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
! endif
call non_hrmt_bieig(mo_num, dm_tmp, thresh_biorthog_diag, thresh_biorthog_nondiag &
, natorb_tc_leigvec_mo, natorb_tc_reigvec_mo &
, mo_num, natorb_tc_eigval )
accu = 0.d0
do i = 1, mo_num
print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)

32
src/tc_bi_ortho/test_tc_two_rdm.irp.f
View File

@ -25,49 +25,47 @@ end
subroutine test
implicit none
integer :: h1,p1,h2,p2,i,j,istate
double precision :: rdm, integral, accu,ref
integer :: h1,p1,h2,p2,i,j,istate,s1,s2
double precision :: rdm, integral, accu,ref, accu_new ,rdm_new
double precision :: hmono, htwoe, hthree, htot
accu = 0.d0
accu_new = 0.d0
do h1 = 1, mo_num
do p1 = 1, mo_num
do h2 = 1, mo_num
do p2 = 1, mo_num
integral = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
rdm = tc_two_rdm(p1,h1,p2,h2)
accu += integral * rdm
rdm_new = 0.d0
do s2 = 1, 2
do s1 = 1, 2
rdm_new += tc_two_rdm_chemist_s1s2(p1,h1,p2,h2,s1,s2)
enddo
enddo
accu_new += integral * rdm_new
! if(dabs(rdm).gt.1.d-10)then
! print*,h1,p1,h2,p2
! print*,rdm,integral,rdm*integral
! endif
accu += integral * rdm
enddo
enddo
enddo
enddo
accu *= 0.5d0
print*,'accu = ',accu
! print*,mo_bi_ortho_tc_two_e(2,15,2,1)
! print*,mo_bi_ortho_tc_two_e(15,2,2,1)
! print*,mo_bi_ortho_tc_two_e(2,1,2,15)
! print*,mo_bi_ortho_tc_two_e(2,1,15,2)
accu_new *= 0.5d0
print*,'accu = ',accu
print*,'accu_new = ',accu_new
ref = 0.d0
do i = 1, N_det
do j = 1, N_det
! if(i.eq.j)cycle
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
do istate = 1,N_states
! print*,'i,j',i,j
! print*,psi_l_coef_bi_ortho(i,istate) , psi_r_coef_bi_ortho(j,istate) , htwoe
! print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe
! if(i.ne.j)then
! print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) , htwoe
! print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe
! endif
ref += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe
enddo
enddo
enddo
print*,' ref = ',ref
print*,' ref = ',ref
print*,'delta= ',ref-accu
end

41
src/tc_bi_ortho/two_rdm_naive.irp.f
View File

@ -1,4 +1,5 @@
BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist_s1s2, (mo_num, mo_num, mo_num, mo_num, 2,2)]
implicit none
BEGIN_DOC
! tc_two_rdm_chemist(p,s,q,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = CHEMIST NOTATION
@ -14,6 +15,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
other_spin(2) = 1
allocate(occ(N_int*bit_kind_size,2))
tc_two_rdm_chemist = 0.d0
tc_two_rdm_chemist_s1s2 = 0.d0
do i = 1, N_det ! psi_left
do j = 1, N_det ! psi_right
@ -21,14 +23,16 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
if(degree.gt.2)cycle
if(degree.gt.0)then
! get excitation operators: from psi_det(j) --> psi_det(i)
call get_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * phase * state_average_weight(1)
do istate = 2, N_states
contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * phase * state_average_weight(istate)
enddo
! T_{j-->i} = a^p1_s1 a_h1_s1
call get_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * phase * state_average_weight(1)
do istate = 2, N_states
contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * phase * state_average_weight(istate)
enddo
if(degree == 2)then
call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
! call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
else if(degree==1)then
! occupation of the determinant psi_det(j)
call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int)
@ -40,6 +44,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
h2 = m
p2 = m
call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
! call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
enddo
! run over the electrons of same spin than the excitation
s2 = s1
@ -48,6 +53,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
h2 = m
p2 = m
call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
! call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
enddo
endif
else if(degree == 0)then
@ -69,6 +75,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
h2 = m
p2 = m
call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
! call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
enddo
! run over the couple of alpha-alpha electrons
s2 = s1
@ -78,6 +85,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
p2 = m
if(h2.le.h1)cycle
call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
! call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
enddo
enddo
s1 = 2
@ -92,6 +100,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
p2 = m
if(h2.le.h1)cycle
call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
! call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
enddo
enddo
endif
@ -124,12 +133,13 @@ subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,contrib)
end
BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, tc_two_rdm_s1s2, (mo_num, mo_num, mo_num, mo_num,2,2)]
implicit none
BEGIN_DOC
! tc_two_rdm(p,q,s,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = PHYSICIST NOTATION
END_DOC
integer :: p,q,r,s
integer :: p,q,r,s,s1,s2
do r = 1, mo_num
do q = 1, mo_num
do s = 1, mo_num
@ -139,5 +149,18 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
enddo
enddo
enddo
do s2 = 1, 2
do s1 = 1, 2
do r = 1, mo_num
do q = 1, mo_num
do s = 1, mo_num
do p = 1, mo_num
tc_two_rdm_s1s2(p,q,s,r,s1,s2) = tc_two_rdm_chemist_s1s2(p,s,q,r,s1,s2)
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER