From 5dc4fb29284c496cf3b6921e9c9422e94b972a97 Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 16 Aug 2023 14:06:29 +0200 Subject: [PATCH] naive two rdm in tc works for He in cisd and bi ortho orbitals --- src/tc_bi_ortho/test_tc_two_rdm.irp.f | 60 ++++++++++++ src/tc_bi_ortho/two_rdm_naive.irp.f | 132 ++++++++++++++++++++++++++ 2 files changed, 192 insertions(+) create mode 100644 src/tc_bi_ortho/test_tc_two_rdm.irp.f create mode 100644 src/tc_bi_ortho/two_rdm_naive.irp.f diff --git a/src/tc_bi_ortho/test_tc_two_rdm.irp.f b/src/tc_bi_ortho/test_tc_two_rdm.irp.f new file mode 100644 index 00000000..ecdeef43 --- /dev/null +++ b/src/tc_bi_ortho/test_tc_two_rdm.irp.f @@ -0,0 +1,60 @@ +program test_tc_rdm + + BEGIN_DOC + ! + ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together + ! with the energy. Saves the left-right wave functions at the end. + ! + END_DOC + + my_grid_becke = .True. + PROVIDE tc_grid1_a tc_grid1_r + my_n_pt_r_grid = tc_grid1_r + my_n_pt_a_grid = tc_grid1_a + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + read_wf = .True. + touch read_wf + + print*, ' nb of states = ', N_states + print*, ' nb of det = ', N_det + + call test() + +end + +subroutine test + implicit none + integer :: h1,p1,h2,p2,i,j,istate + double precision :: rdm, integral, accu,ref + double precision :: hmono, htwoe, hthree, htot + accu = 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 + enddo + enddo + enddo + enddo + accu *= 0.5d0 + print*,'accu = ',accu +! print*,tc_two_rdm(1,1,1,1),mo_bi_ortho_tc_two_e(1,1,1,1) + ref = 0.d0 + do i = 1, N_det + do j = 1, N_det +! if(i.ne.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,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe +! print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,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 + print*,' ref = ',ref + +end diff --git a/src/tc_bi_ortho/two_rdm_naive.irp.f b/src/tc_bi_ortho/two_rdm_naive.irp.f new file mode 100644 index 00000000..9694c653 --- /dev/null +++ b/src/tc_bi_ortho/two_rdm_naive.irp.f @@ -0,0 +1,132 @@ +BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)] + implicit none + BEGIN_DOC + ! tc_two_rdm(p,s,q,r) = 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 + if(degree == 2)then + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + else if(degree==1)then +! cycle + ! occupation of the determinant psi_det(j) + call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) + + ! run over the electrons of opposite spin than the excitation + s2 = other_spin(s1) + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + ! run over the electrons of same spin than the excitation + s2 = s1 + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + if(h2.le.h1)cycle + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + endif + else if(degree == 0)then + contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * state_average_weight(1) +! print*,'contrib',contrib + do istate = 2, N_states + contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) + enddo + ! occupation of the determinant psi_det(j) + call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) + s1 = 1 ! alpha electrons + do nn = 1, n_occ_ab(s1) + h1 = occ(nn,s1) + p1 = occ(nn,s1) + ! run over the couple of alpha-beta electrons + s2 = other_spin(s1) + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + ! run over the couple of alpha-alpha electrons + s2 = s1 + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + if(h2.le.h1)cycle + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + enddo + s1 = 2 + do nn = 1, n_occ_ab(s1) + h1 = occ(nn,s1) + p1 = occ(nn,s1) + ! run over the couple of beta-beta electrons + s2 = s1 + do mm = 1, n_occ_ab(s2) + m = occ(mm,s2) + h2 = m + p2 = m + if(h2.le.h1)cycle + call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib) + enddo + enddo + endif + enddo + enddo + +END_PROVIDER + +subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,contrib) + implicit none + integer, intent(in) :: h1,p1,h2,p2,s1,s2,sze + double precision, intent(in) :: contrib + double precision, intent(inout) :: array(sze, sze, sze, sze) + integer :: istate + if(s1.ne.s2)then + array(p1,h1,p2,h2) += contrib + ! permutation for particle symmetry + array(p2,h2,p1,h1) += contrib + else ! same spin double excitation + array(p1,h1,p2,h2) += contrib + ! exchange + ! exchanging the holes + array(p2,h1,p1,h2) -= contrib + ! exchanging the particles + array(p1,h2,p2,h1) -= contrib + + ! permutation for particle symmetry + array(p2,h2,p1,h1) += contrib + ! exchange + ! exchanging the holes + array(p1,h2,p2,h1) -= contrib + ! exchanging the particles + array(p2,h1,p1,h2) -= contrib + endif + +end