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

Added spin dependent two-rdm.

Browse Source
This commit is contained in:
eginer 2023-09-13 12:58:26 +02:00
parent 2ff59e6869
commit 8b14a2b7ab
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_d = 1.d-6
thr_nd = 1.d-6 thr_nd = 1.d-6
thr_deg = 1.d-3 thr_deg = 1.d-3
if(n_core_orb.ne.0)then ! if(n_core_orb.ne.0)then
! print*,'core orbitals' ! call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg &
! pause ! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg & ! else
, natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval) ! call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg &
else ! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & ! endif
, 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 )
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 accu = 0.d0
do i = 1, mo_num do i = 1, mo_num
print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i) print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)

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

@ -25,44 +25,42 @@ end
subroutine test subroutine test
implicit none implicit none
integer :: h1,p1,h2,p2,i,j,istate integer :: h1,p1,h2,p2,i,j,istate,s1,s2
double precision :: rdm, integral, accu,ref double precision :: rdm, integral, accu,ref, accu_new ,rdm_new
double precision :: hmono, htwoe, hthree, htot double precision :: hmono, htwoe, hthree, htot
accu = 0.d0 accu = 0.d0
accu_new = 0.d0
do h1 = 1, mo_num do h1 = 1, mo_num
do p1 = 1, mo_num do p1 = 1, mo_num
do h2 = 1, mo_num do h2 = 1, mo_num
do p2 = 1, mo_num do p2 = 1, mo_num
integral = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) integral = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
rdm = tc_two_rdm(p1,h1,p2,h2) 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 ! if(dabs(rdm).gt.1.d-10)then
! print*,h1,p1,h2,p2 ! print*,h1,p1,h2,p2
! print*,rdm,integral,rdm*integral ! print*,rdm,integral,rdm*integral
! endif ! endif
accu += integral * rdm
enddo enddo
enddo enddo
enddo enddo
enddo enddo
accu *= 0.5d0 accu *= 0.5d0
accu_new *= 0.5d0
print*,'accu = ',accu print*,'accu = ',accu
! print*,mo_bi_ortho_tc_two_e(2,15,2,1) print*,'accu_new = ',accu_new
! 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)
ref = 0.d0 ref = 0.d0
do i = 1, N_det do i = 1, N_det
do j = 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) 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 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 ref += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe
enddo enddo
enddo enddo

25
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 implicit none
BEGIN_DOC BEGIN_DOC
! tc_two_rdm_chemist(p,s,q,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = CHEMIST NOTATION ! 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 other_spin(2) = 1
allocate(occ(N_int*bit_kind_size,2)) allocate(occ(N_int*bit_kind_size,2))
tc_two_rdm_chemist = 0.d0 tc_two_rdm_chemist = 0.d0
tc_two_rdm_chemist_s1s2 = 0.d0
do i = 1, N_det ! psi_left do i = 1, N_det ! psi_left
do j = 1, N_det ! psi_right do j = 1, N_det ! psi_right
@ -21,6 +23,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
if(degree.gt.2)cycle if(degree.gt.2)cycle
if(degree.gt.0)then if(degree.gt.0)then
! get excitation operators: from psi_det(j) --> psi_det(i) ! get excitation operators: from psi_det(j) --> psi_det(i)
! 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 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) 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) contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * phase * state_average_weight(1)
@ -29,6 +32,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
enddo enddo
if(degree == 2)then 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,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 else if(degree==1)then
! occupation of the determinant psi_det(j) ! occupation of the determinant psi_det(j)
call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) 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 h2 = m
p2 = 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,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
! run over the electrons of same spin than the excitation ! run over the electrons of same spin than the excitation
s2 = s1 s2 = s1
@ -48,6 +53,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
h2 = m h2 = m
p2 = 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,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 endif
else if(degree == 0)then 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 h2 = m
p2 = 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,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
! run over the couple of alpha-alpha electrons ! run over the couple of alpha-alpha electrons
s2 = s1 s2 = s1
@ -78,6 +85,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
p2 = m p2 = m
if(h2.le.h1)cycle 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,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
enddo enddo
s1 = 2 s1 = 2
@ -92,6 +100,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
p2 = m p2 = m
if(h2.le.h1)cycle 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,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
enddo enddo
endif endif
@ -125,11 +134,12 @@ 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 implicit none
BEGIN_DOC BEGIN_DOC
! tc_two_rdm(p,q,s,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = PHYSICIST NOTATION ! tc_two_rdm(p,q,s,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = PHYSICIST NOTATION
END_DOC END_DOC
integer :: p,q,r,s integer :: p,q,r,s,s1,s2
do r = 1, mo_num do r = 1, mo_num
do q = 1, mo_num do q = 1, mo_num
do s = 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 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 END_PROVIDER