LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-02 11:25:26 +02:00
Code Issues Releases Wiki Activity

2RDM

Browse Source
This commit is contained in:
Anthony Scemama 2018-10-16 16:40:04 +02:00
parent e2b4b1f58d
commit ee13470949
5 changed files with 184 additions and 268 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
plugins/DensityMatrix/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Bitmask
Bitmask Determinants

211
plugins/DensityMatrix/density_matrix.irp.f
View File

@ -1,211 +0,0 @@
use bitmasks
BEGIN_PROVIDER [ integer, iunit_two_body_dm_aa ]
&BEGIN_PROVIDER [ integer, iunit_two_body_dm_ab ]
&BEGIN_PROVIDER [ integer, iunit_two_body_dm_bb ]
implicit none
use bitmasks
BEGIN_DOC
! Temporary files for 2-body dm calculation
END_DOC
integer :: getUnitAndOpen
iunit_two_body_dm_aa = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_aa.tmp','w')
iunit_two_body_dm_ab = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_ab.tmp','w')
iunit_two_body_dm_bb = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_bb.tmp','w')
! Compute two body DM in file
integer :: k,l,degree, idx,i
integer :: exc(0:2,2,2),n_occ_alpha
double precision :: phase, coef
integer :: h1,h2,p1,p2,s1,s2
double precision :: ck, cl
character*(128), parameter :: f = '(i8,4(x,i5),x,d16.8)'
do k=1,det_num
ck = (det_coef_provider(k)+det_coef_provider(k))
do l=1,k-1
cl = det_coef_provider(l)
call get_excitation_degree(det_provider(1,1,k),det_provider(1,1,l),degree,N_int)
if (degree == 2) then
call get_double_excitation(det_provider(1,1,k),det_provider(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
call bielec_integrals_index(h1,h2,p1,p2,idx)
ckl = phase*ck*cl
select case (s1+s2)
case(2) ! alpha alpha
write(iunit_two_body_dm_aa,f) idx, h1,h2,p1,p2, ckl
call bielec_integrals_index(h1,h2,p2,p1,idx)
write(iunit_two_body_dm_aa,f) idx, h1,h2,p2,p1, -ckl
case(3) ! alpha beta
write(iunit_two_body_dm_ab,f) idx, h1,h2,p1,p2, ckl
case(4) ! beta beta
write(iunit_two_body_dm_bb,f) idx, h1,h2,p1,p2, ckl
call bielec_integrals_index(h1,h2,p2,p1,idx)
write(iunit_two_body_dm_bb,f) idx, h1,h2,p2,p1, -ckl
end select
else if (degree == 1) then
call get_mono_excitation(det_provider(1,1,k),det_provider(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
double precision :: ckl
ckl = phase*ck*cl
call bitstring_to_list(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int)
call bitstring_to_list(det_provider(1,2,k), occ(1,2), n_occ_alpha, N_int)
select case (s1)
case (1) ! Alpha single excitation
integer :: occ(N_int*bit_kind_size,2)
do i = 1, elec_alpha_num
p2=occ(i,1)
h2=p2
call bielec_integrals_index(h1,h2,p1,p2,idx)
write(iunit_two_body_dm_aa,f) idx, h1,h2,p1,p2, ckl
call bielec_integrals_index(h1,h2,p2,p1,idx)
write(iunit_two_body_dm_aa,f) idx, h1,h2,p2,p1, -ckl
enddo
do i = 1, elec_beta_num
p2=occ(i,2)
h2=p2
call bielec_integrals_index(h1,h2,p1,p2,idx)
write(iunit_two_body_dm_ab,f) idx, h1,h2,p1,p2, ckl
enddo
case (2) ! Beta single excitation
do i = 1, elec_alpha_num
p2=occ(i,1)
h2=p2
call bielec_integrals_index(h1,h2,p1,p2,idx)
write(iunit_two_body_dm_ab,f) idx, h1,h2,p1,p2, ckl
enddo
do i = 1, elec_beta_num
p2=occ(i,2)
h2=p2
call bielec_integrals_index(h1,h2,p1,p2,idx)
write(iunit_two_body_dm_bb,f) idx, h1,h2,p1,p2, ckl
call bielec_integrals_index(h1,h2,p2,p1,idx)
write(iunit_two_body_dm_bb,f) idx, h1,h2,p2,p1, -ckl
enddo
end select
endif
enddo
enddo
! Sort file
! Merge coefs
close(iunit_two_body_dm_aa)
close(iunit_two_body_dm_ab)
close(iunit_two_body_dm_bb)
character*(128) :: filename
filename = trim(ezfio_filename)//'/work/two_body_aa.tmp'
call system('sort -n '//trim(filename)//' > '//trim(filename)//'2 ; cp '//trim(filename)//'2 '//trim(filename))
filename = trim(ezfio_filename)//'/work/two_body_ab.tmp'
call system('sort -n '//trim(filename)//' > '//trim(filename)//'2 ; cp '//trim(filename)//'2 '//trim(filename))
filename = trim(ezfio_filename)//'/work/two_body_bb.tmp'
call system('sort -n '//trim(filename)//' > '//trim(filename)//'2 ; cp '//trim(filename)//'2 '//trim(filename))
iunit_two_body_dm_aa = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_aa.tmp','r')
iunit_two_body_dm_ab = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_ab.tmp','r')
iunit_two_body_dm_bb = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_bb.tmp','r')
END_PROVIDER
BEGIN_TEMPLATE
BEGIN_PROVIDER [ integer, size_two_body_dm_$AA ]
implicit none
use bitmasks
BEGIN_DOC
! Size of the two body $ALPHA density matrix
END_DOC
integer *8 :: key, key_old
rewind(iunit_two_body_dm_$AA)
size_two_body_dm_$AA = 0
key = 0_8
key_old = key
do while (.True.)
read(iunit_two_body_dm_$AA,*,END=99) key
if (key /= key_old) then
size_two_body_dm_$AA += 1
key_old = key
endif
end do
99 continue
END_PROVIDER
BEGIN_PROVIDER [ integer, two_body_dm_index_$AA, (4,size_two_body_dm_$AA) ]
&BEGIN_PROVIDER [ double precision, two_body_dm_value_$AA, (size_two_body_dm_$AA) ]
implicit none
use bitmasks
BEGIN_DOC
! Two body $ALPHA density matrix
END_DOC
rewind(iunit_two_body_dm_$AA)
integer *8 :: key, key_old
integer :: ii, i,j,k,l
double precision :: c
key = 0_8
key_old = key
ii = 0
do while (.True.)
read(iunit_two_body_dm_$AA,*,END=99) key, i,j,k,l, c
if (key /= key_old) then
ii += 1
two_body_dm_index_$AA(1,ii) = i
two_body_dm_index_$AA(2,ii) = j
two_body_dm_index_$AA(3,ii) = k
two_body_dm_index_$AA(4,ii) = l
two_body_dm_value_$AA(ii) = 0.d0
key_old = key
endif
two_body_dm_value_$AA(ii) += c
enddo
99 continue
close(iunit_two_body_dm_$AA, status='DELETE')
END_PROVIDER
SUBST [ AA, ALPHA ]
aa ; alpha-alpha ;;
ab ; alpha-beta ;;
bb ; beta-beta ;;
END_TEMPLATE
BEGIN_PROVIDER [ double precision, two_body_dm_diag_aa, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [ double precision, two_body_dm_diag_bb, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [ double precision, two_body_dm_diag_ab, (mo_tot_num,mo_tot_num)]
implicit none
use bitmasks
BEGIN_DOC
! diagonal part of the two body density matrix
END_DOC
integer :: i,j,k,e1,e2
integer :: occ(N_int*bit_kind_size,2)
double precision :: ck
integer :: n_occ_alpha
two_body_dm_diag_aa=0.d0
two_body_dm_diag_ab=0.d0
two_body_dm_diag_bb=0.d0
do k = 1, det_num
call bitstring_to_list(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int)
call bitstring_to_list(det_provider(1,2,k), occ(1,2), n_occ_alpha, N_int)
ck = det_coef_provider(k) * det_coef_provider(k)
do i = 1,elec_alpha_num
e1=occ(i,1)
do j = 1,elec_alpha_num
e2=occ(j,1)
! alpha-alpha
two_body_dm_diag_aa(e1,e2) = two_body_dm_diag_aa(e1,e2) + ck
enddo
do j = 1,elec_beta_num
e2=occ(j,2)
! alpha-beta
two_body_dm_diag_ab(e1,e2) = two_body_dm_diag_ab(e1,e2) + ck
enddo
enddo
do i = 1,elec_beta_num
e1=occ(i,2)
do j = 1,elec_beta_num
e2=occ(j,2)
! beta-beta
two_body_dm_diag_bb(e1,e2) = two_body_dm_diag_bb(e1,e2) + ck
enddo
enddo
enddo
END_PROVIDER

116
plugins/DensityMatrix/density_matrix_array.irp.f Normal file
View File

@ -0,0 +1,116 @@
use bitmasks
BEGIN_PROVIDER [ double precision, two_body_dm_aa, (mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, two_body_dm_bb, (mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, two_body_dm_ab, (mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num) ]
implicit none
use bitmasks
BEGIN_DOC
! Temporary files for 2-body dm calculation
END_DOC
integer :: getUnitAndOpen
! Compute two body DM in file
integer :: k,l,degree, idx,i,j
integer :: exc(0:2,2,2),n_occ_alpha
double precision :: phase, coef
integer :: h1,h2,p1,p2,s1,s2, e1, e2
double precision :: ck, cl
character*(128), parameter :: f = '(i8,4(x,i5),x,d16.8)'
integer :: istate
two_body_dm_aa = 0.d0
two_body_dm_ab = 0.d0
two_body_dm_bb = 0.d0
istate = 1
! OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,ck,ckl,i,j,e1,e2,cl,phase,h1,p1,h2,p2,s1,s2,occ)
! OMP DO SCHEDULE(dynamic,64)
do k=1,N_det
ck = psi_coef(k,istate)
call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int)
call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int)
ckl = psi_coef(k,istate) * psi_coef(k,istate)
do i = 1,elec_alpha_num
e1=occ(i,1)
do j = 1,elec_alpha_num
e2=occ(j,1)
! alpha-alpha
two_body_dm_aa(e1,e2,e1,e2) += 0.5d0*ckl
two_body_dm_aa(e1,e2,e2,e1) -= 0.5d0*ckl
enddo
do j = 1,elec_beta_num
e2=occ(j,2)
! alpha-beta
two_body_dm_ab(e1,e2,e1,e2) += ckl
enddo
enddo
do i = 1,elec_beta_num
e1=occ(i,2)
do j = 1,elec_beta_num
e2=occ(j,2)
! beta-beta
two_body_dm_bb(e1,e2,e1,e2) += 0.5d0*ckl
two_body_dm_bb(e1,e2,e2,e1) -= 0.5d0*ckl
enddo
enddo
do l=1,k-1
cl = 2.d0*psi_coef(l,istate)
call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int)
if (degree == 2) then
call get_double_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
ckl = phase*ck*cl
select case (s1+s2)
case(2) ! alpha alpha
two_body_dm_aa(h1,h2,p1,p2) += ckl
two_body_dm_aa(h1,h2,p2,p1) -= ckl
case(3) ! alpha beta
two_body_dm_ab(h1,h2,p1,p2) += ckl
case(4) ! beta beta
two_body_dm_bb(h1,h2,p1,p2) += ckl
two_body_dm_bb(h1,h2,p2,p1) -= ckl
end select
else if (degree == 1) then
call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
double precision :: ckl
ckl = phase*ck*cl
call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int)
call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int)
select case (s1)
case (1) ! Alpha single excitation
integer :: occ(N_int*bit_kind_size,2)
do i = 1, elec_alpha_num
p2=occ(i,1)
h2=p2
two_body_dm_aa(h1,h2,p1,p2) += ckl
two_body_dm_aa(h1,h2,p2,p1) -= ckl
enddo
do i = 1, elec_beta_num
p2=occ(i,2)
h2=p2
two_body_dm_ab(h1,h2,p1,p2) += ckl
enddo
case (2) ! Beta single excitation
do i = 1, elec_alpha_num
p2=occ(i,1)
h2=p2
two_body_dm_ab(h1,h2,p1,p2) += ckl
enddo
do i = 1, elec_beta_num
p2=occ(i,2)
h2=p2
two_body_dm_bb(h1,h2,p1,p2) += ckl
two_body_dm_bb(h1,h2,p2,p1) -= ckl
enddo
end select
endif
enddo
enddo
! OMP END DO
! OMP END PARALLEL
END_PROVIDER

56
plugins/DensityMatrix/det_num.irp.f
View File

@ -1,56 +0,0 @@
use bitmasks
BEGIN_PROVIDER [integer, det_num]
det_num = 10
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), det_provider, (N_int,2,det_num)]
&BEGIN_PROVIDER [ double precision , det_coef_provider, (det_num) ]
use bitmasks
implicit none
integer :: i
det_provider = 0
det_provider(1,1,1 ) = #001f ! 0000 0000 0001 1111
det_provider(1,1,2 ) = #003b ! 0000 0000 0011 1011
det_provider(1,1,3 ) = #008f ! 0000 0000 1000 1111
det_provider(1,1,4 ) = #0057 ! 0000 0000 0101 0111
det_provider(1,1,5 ) = #100f ! 0001 0000 0000 1111
det_provider(1,1,6 ) = #001f ! 0000 0000 0001 1111
det_provider(1,1,7 ) = #003b ! 0000 0000 0011 1011
det_provider(1,1,8 ) = #00c7 ! 0000 0000 1100 0111
det_provider(1,1,9 ) = #00ab ! 0000 0000 1010 1011
det_provider(1,1,10) = #0073 ! 0000 0000 0111 0011
det_provider(1,2,1 ) = #0007 ! 0000 0000 0001 0111
det_provider(1,2,2 ) = #0023 ! 0000 0000 0010 0011
det_provider(1,2,3 ) = #0023 ! 0000 0000 0010 0011
det_provider(1,2,4 ) = #0023 ! 0000 0000 0010 0011
det_provider(1,2,5 ) = #0015 ! 0000 0000 0001 0101
det_provider(1,2,6 ) = #000d ! 0000 0000 0000 1101
det_provider(1,2,7 ) = #0007 ! 0000 0000 0000 0111
det_provider(1,2,8 ) = #0007 ! 0000 0000 0000 0111
det_provider(1,2,9 ) = #0007 ! 0000 0000 0000 0111
det_provider(1,2,10) = #0007 ! 0000 0000 0000 0111
det_coef_provider = (/ &
0.993536117982429D+00, &
-0.556089064313864D-01, &
0.403074722590178D-01, &
0.403074717461626D-01, &
-0.340290975461932D-01, &
-0.340290958781670D-01, &
-0.333949939765448D-01, &
0.333418373363987D-01, &
-0.316337211787351D-01, &
-0.316337207748718D-01 &
/)
!do i=1,10
! call write_bitstring( 6, det_provider(1,1,i), N_int )
!enddo
!print *, ''
!do i=1,10
! call write_bitstring( 6, det_provider(1,2,i), N_int )
!enddo
!print *, ''
END_PROVIDER

67
plugins/DensityMatrix/print_2rdm.irp.f Normal file
View File

@ -0,0 +1,67 @@
program pouet
implicit none
integer :: i,j,k,l
double precision, external :: get_mo_bielec_integral
read_wf = .True.
TOUCH read_wf
double precision :: e(10)
e = 0.d0
print *, '1RDM ALPHA'
do i=1,mo_tot_num
do j=1,mo_tot_num
print *, i, j, one_body_dm_mo_alpha(i,j,1)
e(4) += one_body_dm_mo_alpha(i,j,1) * mo_mono_elec_integral(i,j)
enddo
enddo
print *, '1RDM BETA'
do i=1,mo_tot_num
do j=1,mo_tot_num
print *, i, j, one_body_dm_mo_beta(i,j,1)
e(4) += one_body_dm_mo_beta(i,j,1) * mo_mono_elec_integral(i,j)
enddo
enddo
print *, '2RDM ALPHA ALPHA'
do i=1,mo_tot_num
do j=1,mo_tot_num
do k=1,mo_tot_num
do l=1,mo_tot_num
print *, i, j, k, l, two_body_dm_aa(i,j,k,l)
e(1) += two_body_dm_aa(i,j,k,l) * get_mo_bielec_integral(i,j,k,l, mo_integrals_map)
enddo
enddo
enddo
enddo
print *, '2RDM BETA BETA'
do i=1,mo_tot_num
do j=1,mo_tot_num
do k=1,mo_tot_num
do l=1,mo_tot_num
print *, i, j, k, l, two_body_dm_bb(i,j,k,l)
e(2) += two_body_dm_bb(i,j,k,l) * get_mo_bielec_integral(i,j,k,l, mo_integrals_map)
enddo
enddo
enddo
enddo
print *, '2RDM ALPHA BETA'
do i=1,mo_tot_num
do j=1,mo_tot_num
do k=1,mo_tot_num
do l=1,mo_tot_num
print *, i, j, k, l, two_body_dm_ab(i,j,k,l)
e(3) += two_body_dm_ab(i,j,k,l) * get_mo_bielec_integral(i,j,k,l, mo_integrals_map)
enddo
enddo
enddo
enddo
print *, ''
print *, 'Energy ', sum(e(1:4)) + nuclear_repulsion
end