1
0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-11-07 22:53:42 +01:00
qp_plugins_scemama/devel/density_matrix/print_2rdm.irp.f

44 lines
1.3 KiB
Fortran

program print_2rdm
implicit none
integer :: i,j,k,l
double precision, external :: get_two_e_integral
read_wf = .True.
TOUCH read_wf
double precision :: e(10)
double precision, parameter :: thr = 1.d-15
e = 0.d0
print *, '1RDM'
do i=1,mo_num
do j=1,mo_num
if (dabs(one_e_dm_mo_alpha(i,j,1) + one_e_dm_mo_beta(i,j,1)) > thr) then
print *, i, j, one_e_dm_mo_alpha(i,j,1) + one_e_dm_mo_beta(i,j,1)
endif
e(4) += one_e_dm_mo_alpha(i,j,1) * mo_one_e_integrals(i,j)
e(4) += one_e_dm_mo_beta(i,j,1) * mo_one_e_integrals(i,j)
enddo
enddo
print *, '2RDM'
do i=1,mo_num
do j=1,mo_num
do k=1,mo_num
do l=1,mo_num
if (dabs(two_e_dm_aa(i,j,k,l) + two_e_dm_bb(i,j,k,l) + two_e_dm_ab(i,j,k,l)) > thr) then
print *, i, j, k, l, two_e_dm_aa(i,j,k,l) + two_e_dm_bb(i,j,k,l) + two_e_dm_ab(i,j,k,l)
endif
e(1) += two_e_dm_aa(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
e(2) += two_e_dm_bb(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
e(3) += two_e_dm_ab(i,j,k,l) * get_two_e_integral(i,j,k,l, mo_integrals_map)
enddo
enddo
enddo
enddo
print *, ''
print *, 'Energy ', sum(e(1:4)) + nuclear_repulsion
end