diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f index 1b5a66f3..a72d356a 100644 --- a/src/tc_bi_ortho/tc_natorb.irp.f +++ b/src/tc_bi_ortho/tc_natorb.irp.f @@ -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) diff --git a/src/tc_bi_ortho/test_tc_two_rdm.irp.f b/src/tc_bi_ortho/test_tc_two_rdm.irp.f index 3e556312..044c31e0 100644 --- a/src/tc_bi_ortho/test_tc_two_rdm.irp.f +++ b/src/tc_bi_ortho/test_tc_two_rdm.irp.f @@ -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 diff --git a/src/tc_bi_ortho/two_rdm_naive.irp.f b/src/tc_bi_ortho/two_rdm_naive.irp.f index 3963d09e..d21d6a87 100644 --- a/src/tc_bi_ortho/two_rdm_naive.irp.f +++ b/src/tc_bi_ortho/two_rdm_naive.irp.f @@ -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) = = 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) = = 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