mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-06-25 22:12:05 +02:00
naive two rdm in tc works for He in cisd and bi ortho orbitals
This commit is contained in:
eginer
parent
a15055e964
commit
5dc4fb2928
60
src/tc_bi_ortho/test_tc_two_rdm.irp.f
Normal file
60
src/tc_bi_ortho/test_tc_two_rdm.irp.f
Normal file
|
|
@ -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
|
||||
132
src/tc_bi_ortho/two_rdm_naive.irp.f
Normal file
132
src/tc_bi_ortho/two_rdm_naive.irp.f
Normal file
|
|
@ -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) = <Phi| a^dager_p
|
||||
END_DOC
|
||||
integer :: i,j,istate,m,mm,nn
|
||||
integer :: exc(0:2,2,2)
|
||||
double precision :: phase
|
||||
double precision :: contrib
|
||||
integer :: h1,p1,s1,h2,p2,s2,degree
|
||||
integer, allocatable :: occ(:,:)
|
||||
integer :: n_occ_ab(2),other_spin(2)
|
||||
other_spin(1) = 2
|
||||
other_spin(2) = 1
|
||||
allocate(occ(N_int*bit_kind_size,2))
|
||||
tc_two_rdm = 0.d0
|
||||
|
||||
do i = 1, N_det ! psi_left
|
||||
do j = 1, N_det ! psi_right
|
||||
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
|
||||
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
|
||||
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
|
||||
Loading…
Reference in New Issue
Block a user