mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-22 20:35:19 +01:00
106 lines
2.9 KiB
Fortran
106 lines
2.9 KiB
Fortran
program test_two_bod
|
|
implicit none
|
|
read_wf = .True.
|
|
touch read_wf
|
|
call routine
|
|
end
|
|
subroutine routine
|
|
implicit none
|
|
integer :: i,j,k,l
|
|
integer :: h1,p1,h2,p2,s1,s2
|
|
double precision :: accu,get_two_body_dm_ab_map_element,get_mo_bielec_integral
|
|
accu = 0.d0
|
|
|
|
! Diag part of the core two body dm
|
|
do i = 1, n_core_orb
|
|
h1 = list_core(i)
|
|
do j = 1, n_core_orb
|
|
h2 = list_core(j)
|
|
accu += two_body_dm_ab_diag_core(j,i) * mo_bielec_integral_jj(h1,h2)
|
|
enddo
|
|
enddo
|
|
|
|
! Diag part of the active two body dm
|
|
do i = 1, n_act_orb
|
|
h1 = list_act(i)
|
|
do j = 1, n_act_orb
|
|
h2 = list_act(j)
|
|
accu += two_body_dm_ab_diag_act(j,i) * mo_bielec_integral_jj(h1,h2)
|
|
enddo
|
|
enddo
|
|
|
|
! Diag part of the active <-> core two body dm
|
|
do i = 1, n_act_orb
|
|
h1 = list_act(i)
|
|
do j = 1, n_core_orb
|
|
h2 = list_core(j)
|
|
accu += two_body_dm_diag_core_act(j,i) * mo_bielec_integral_jj(h1,h2)
|
|
enddo
|
|
enddo
|
|
print*,'BI ELECTRONIC = ',accu
|
|
|
|
double precision :: accu_extra_diag
|
|
accu_extra_diag = 0.d0
|
|
! purely active part of the two body dm
|
|
do l = 1, n_act_orb ! p2
|
|
p2 = list_act(l)
|
|
do k = 1, n_act_orb ! h2
|
|
h2 = list_act(k)
|
|
do j = 1, n_act_orb ! p1
|
|
p1 = list_act(j)
|
|
do i = 1,n_act_orb ! h1
|
|
h1 = list_act(i)
|
|
accu_extra_diag += two_body_dm_ab_big_array_act(i,j,k,l) * get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
! core <-> active part of the two body dm
|
|
do l = 1, n_act_orb ! p1
|
|
p1 = list_act(l)
|
|
do k = 1, n_act_orb ! h1
|
|
h1 = list_act(k)
|
|
do i = 1,n_core_orb ! h2
|
|
h2 = list_core(i)
|
|
accu_extra_diag += two_body_dm_ab_big_array_core_act(i,k,l) * get_mo_bielec_integral(h1,h2,p1,h2,mo_integrals_map)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
print*,'extra_diag = ',accu_extra_diag
|
|
double precision :: average_mono
|
|
call get_average(mo_mono_elec_integral,one_body_dm_mo,average_mono)
|
|
print*,'BI ELECTRONIC = ',accu+accu_extra_diag
|
|
print*,'MONO ELECTRONIC = ',average_mono
|
|
print*,'Total elec = ',accu+average_mono + accu_extra_diag
|
|
print*,'Total = ',accu+average_mono+nuclear_repulsion +accu_extra_diag
|
|
double precision :: e_0,hij
|
|
call u_0_H_u_0(e_0,psi_coef,n_det,psi_det,N_int)
|
|
print*,'<Psi| H |Psi> = ',e_0 + nuclear_repulsion
|
|
integer :: degree
|
|
integer :: exc(0:2,2,2)
|
|
double precision :: phase
|
|
integer :: n_elements
|
|
n_elements = 0
|
|
accu = 0.d0
|
|
do i = 1, N_det
|
|
do j = i+1, N_det
|
|
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
|
|
if(degree.gt.2)cycle
|
|
! if(degree.ne.1)cycle
|
|
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
|
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
if(s1.eq.s2)cycle
|
|
n_elements += 1
|
|
call i_H_j(psi_det(1,1,i),psi_det(1,1,j),N_int,hij)
|
|
accu += 2.d0 * hij * psi_coef(i,1) * psi_coef(j,1)
|
|
enddo
|
|
enddo
|
|
print*,'n_elements = ',n_elements
|
|
print*,'<Psi| extra diag ',accu
|
|
print*,'dm ',accu_extra_diag
|
|
|
|
|
|
end
|