mirror of https://gitlab.com/scemama/eplf
Introduced second order density matrix
This commit is contained in:
Anthony Scemama
parent
15082646ca
commit
cd0cc7b1a5
255
src/det.irp.f
255
src/det.irp.f
|
|
@ -255,3 +255,258 @@ BEGIN_PROVIDER [ real, ci_mo, (mo_num,mo_num,3) ]
|
|||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, two_e_density_num_max ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of factors containing the Slater rules
|
||||
END_DOC
|
||||
|
||||
two_e_density_num_max = 2*mo_num
|
||||
|
||||
integer :: k,l
|
||||
integer :: exc(3), nact, nact2, p, p2
|
||||
integer :: det_exc
|
||||
do k=1,det_num
|
||||
do l=k,det_num
|
||||
exc(1) = abs(det_exc(k,l,1))
|
||||
exc(2) = abs(det_exc(k,l,2))
|
||||
exc(3) = exc(1)+exc(2)
|
||||
|
||||
do p=1,2
|
||||
p2 = 1+mod(p,2)
|
||||
nact = elec_num_2(p) -mo_closed_num
|
||||
nact2 = elec_num_2(p2)-mo_closed_num
|
||||
if ( exc(3) == 0 ) then
|
||||
two_e_density_num_max += 2*nact*mo_num
|
||||
else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
|
||||
two_e_density_num_max += 2*mo_num
|
||||
else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
|
||||
two_e_density_num_max += 2
|
||||
else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
|
||||
two_e_density_num_max += 1
|
||||
endif
|
||||
enddo
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, two_e_density_indice, (4,two_e_density_num_max) ]
|
||||
&BEGIN_PROVIDER [ real, two_e_density_value, (2,two_e_density_num_max) ]
|
||||
&BEGIN_PROVIDER [ integer, two_e_density_num ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Compact representation of eplf factors
|
||||
END_DOC
|
||||
|
||||
integer :: i,j,k,l,m
|
||||
|
||||
integer :: n,p,p2,q
|
||||
integer :: ik,il,jk,jl, idx(4)
|
||||
real :: phase
|
||||
integer :: exc(4), nact, nact2
|
||||
real :: det_kl
|
||||
integer :: det_exc
|
||||
|
||||
two_e_density_num = 0
|
||||
|
||||
PROVIDE det
|
||||
|
||||
do k=1,det_num
|
||||
do l=k,det_num
|
||||
|
||||
exc(1) = det_exc(k,l,1)
|
||||
exc(2) = det_exc(k,l,2)
|
||||
exc(4) = exc(1)*exc(2)
|
||||
exc(1) = abs(exc(1))
|
||||
exc(2) = abs(exc(2))
|
||||
exc(3) = exc(1)+exc(2)
|
||||
if (exc(4) /= 0) then
|
||||
exc(4) = exc(4)/abs(exc(4))
|
||||
else
|
||||
exc(4) = 1
|
||||
endif
|
||||
phase = dble(exc(4))
|
||||
|
||||
det_kl = phase*det_coef(k)*det_coef(l)
|
||||
if (k /= l) then
|
||||
det_kl += det_kl
|
||||
endif
|
||||
|
||||
logical :: notfound
|
||||
BEGIN_SHELL [ /usr/bin/python ]
|
||||
code = """
|
||||
notfound = .True.
|
||||
idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /)
|
||||
do q=1,two_e_density_num
|
||||
if (sum(abs(two_e_density_indice(:,q)-idx))) then
|
||||
two_e_density_value(1,q) += det_kl
|
||||
two_e_density_value(2,q) += det_kl
|
||||
notfound = .False.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (notfound) then
|
||||
two_e_density_num += 1
|
||||
two_e_density_indice(:,two_e_density_num)=idx
|
||||
two_e_density_value(1,two_e_density_num) = det_kl
|
||||
two_e_density_value(2,two_e_density_num) = det_kl
|
||||
endif
|
||||
|
||||
notfound = .True.
|
||||
idx = (/ min(%(I)s,%(K)s), max(%(I)s,%(K)s), min(%(J)s,%(L)s), max(%(J)s,%(L)s) /)
|
||||
do q=1,two_e_density_num
|
||||
if (sum(abs(two_e_density_indice(:,q)-idx))) then
|
||||
two_e_density_value(1,q) -= det_kl
|
||||
notfound = .False.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (notfound) then
|
||||
two_e_density_num += 1
|
||||
two_e_density_indice(:,two_e_density_num)=idx
|
||||
two_e_density_value(1,two_e_density_num) = -det_kl
|
||||
two_e_density_value(2,two_e_density_num) = 0.
|
||||
endif
|
||||
"""
|
||||
|
||||
code1 = """
|
||||
idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /)
|
||||
notfound = .True.
|
||||
do q=1,two_e_density_num
|
||||
if (sum(abs(two_e_density_indice(:,q)-idx))) then
|
||||
two_e_density_value(1,q) += det_kl
|
||||
notfound = .False.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (notfound) then
|
||||
two_e_density_num += 1
|
||||
two_e_density_indice(:,two_e_density_num)=idx
|
||||
two_e_density_value(1,two_e_density_num) = det_kl
|
||||
two_e_density_value(2,two_e_density_num) = 0.
|
||||
endif
|
||||
|
||||
notfound = .True.
|
||||
idx = (/ min(%(I)s,%(K)s), max(%(I)s,%(K)s), min(%(J)s,%(L)s), max(%(J)s,%(L)s) /)
|
||||
do q=1,two_e_density_num
|
||||
if (sum(abs(two_e_density_indice(:,q)-idx))) then
|
||||
two_e_density_value(1,q) -= det_kl
|
||||
notfound = .False.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (notfound) then
|
||||
two_e_density_num += 1
|
||||
two_e_density_indice(:,two_e_density_num)=idx
|
||||
two_e_density_value(1,two_e_density_num) = -det_kl
|
||||
two_e_density_value(2,two_e_density_num) = 0.
|
||||
endif
|
||||
"""
|
||||
|
||||
code2 = """
|
||||
notfound = .True.
|
||||
idx = (/ min(%(I)s,%(J)s), max(%(I)s,%(J)s), min(%(K)s,%(L)s), max(%(K)s,%(L)s) /)
|
||||
do q=1,two_e_density_num
|
||||
if (sum(abs(two_e_density_indice(:,q)-idx))) then
|
||||
two_e_density_value(2,q) += det_kl
|
||||
notfound = .False.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (notfound) then
|
||||
two_e_density_num += 1
|
||||
two_e_density_indice(:,two_e_density_num)=idx
|
||||
two_e_density_value(1,two_e_density_num) = 0.
|
||||
two_e_density_value(2,two_e_density_num) = det_kl
|
||||
endif
|
||||
"""
|
||||
|
||||
rep = { \
|
||||
'CLOSED' : code%{ 'I':'ik', 'J':'il', 'K':'j', 'L':'j' },
|
||||
'OPEN_CLOSED' : code%{ 'I':'j', 'J':'j', 'K':'ik', 'L':'il' },
|
||||
'OPEN_OPEN_1' : code1%{ 'I':'ik', 'J':'il', 'K':'jk', 'L':'jl' },
|
||||
'OPEN_OPEN_2' : code2%{ 'I':'ik', 'J':'il', 'K':'jk', 'L':'jl' }
|
||||
}
|
||||
|
||||
|
||||
print """
|
||||
|
||||
do p=1,2
|
||||
p2 = 1+mod(p,2)
|
||||
nact = elec_num_2(p) -mo_closed_num
|
||||
nact2 = elec_num_2(p2)-mo_closed_num
|
||||
|
||||
if ( exc(3) == 0 ) then
|
||||
do n=1,nact
|
||||
ik = det(n,k,p)
|
||||
il = det(n,l,p)
|
||||
do j=1,mo_closed_num
|
||||
! Closed-open shell interactions
|
||||
%(CLOSED)s
|
||||
!- Open-closed shell interactions
|
||||
%(OPEN_CLOSED)s
|
||||
enddo
|
||||
|
||||
!- Open-open shell interactions
|
||||
do m=1,nact
|
||||
jk = det(m,k,p)
|
||||
jl = det(m,l,p)
|
||||
%(OPEN_OPEN_1)s
|
||||
enddo
|
||||
do m=1,nact2
|
||||
jk = det(m,k,p2)
|
||||
jl = det(m,l,p2)
|
||||
%(OPEN_OPEN_2)s
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
|
||||
|
||||
! Sum over only the sigma-sigma interactions involving the excitation
|
||||
call get_single_excitation(k,l,ik,il,p)
|
||||
|
||||
do j=1,mo_closed_num
|
||||
!- Open-closed shell interactions
|
||||
%(CLOSED)s
|
||||
!- Closed-open shell interactions
|
||||
%(OPEN_CLOSED)s
|
||||
enddo
|
||||
|
||||
!- Open-open shell interactions
|
||||
do m=1,nact
|
||||
jk = det(m,k,p)
|
||||
jl = det(m,l,p)
|
||||
%(OPEN_OPEN_1)s
|
||||
enddo
|
||||
do m=1,nact2
|
||||
jk = det(m,k,p2)
|
||||
jl = det(m,l,p2)
|
||||
%(OPEN_OPEN_2)s
|
||||
enddo
|
||||
|
||||
else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
|
||||
|
||||
! Consider only the double excitations of same-spin electrons
|
||||
call get_double_excitation(k,l,ik,il,jk,jl,p)
|
||||
%(OPEN_OPEN_1)s
|
||||
|
||||
else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
|
||||
|
||||
! Consider only the double excitations of opposite-spin electrons
|
||||
call get_single_excitation(k,l,ik,il,p)
|
||||
call get_single_excitation(k,l,jk,jl,p2)
|
||||
%(OPEN_OPEN_2)s
|
||||
|
||||
endif
|
||||
|
||||
enddo
|
||||
"""%(rep)
|
||||
END_SHELL
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
|
|
|||
|
|
@ -96,228 +96,18 @@ END_PROVIDER
|
|||
|
||||
integer :: k,l,m
|
||||
|
||||
do m=1,eplf_factor_num_max
|
||||
i=eplf_factor_indice(1,m)
|
||||
j=eplf_factor_indice(2,m)
|
||||
k=eplf_factor_indice(3,m)
|
||||
l=eplf_factor_indice(4,m)
|
||||
do m=1,two_e_density_num
|
||||
i=two_e_density_indice(1,m)
|
||||
j=two_e_density_indice(2,m)
|
||||
k=two_e_density_indice(3,m)
|
||||
l=two_e_density_indice(4,m)
|
||||
temp = mo_value_prod_p(i,j)*mo_eplf_integral_matrix(k,l)
|
||||
eplf_up_up += eplf_factor_value(1,m)*temp
|
||||
eplf_up_dn += eplf_factor_value(2,m)*temp
|
||||
eplf_up_up += two_e_density_value(1,m)*temp
|
||||
eplf_up_dn += two_e_density_value(2,m)*temp
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, eplf_factor_num_max ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of factors containing the Slater rules
|
||||
END_DOC
|
||||
|
||||
eplf_factor_num_max = 0
|
||||
|
||||
integer :: k,l
|
||||
integer :: exc(3), nact, nact2, p, p2
|
||||
integer :: det_exc
|
||||
do k=1,det_num
|
||||
do l=k,det_num
|
||||
exc(1) = det_exc(k,l,1)
|
||||
exc(2) = det_exc(k,l,2)
|
||||
exc(4) = exc(1)*exc(2)
|
||||
exc(1) = abs(exc(1))
|
||||
exc(2) = abs(exc(2))
|
||||
exc(3) = exc(1)+exc(2)
|
||||
|
||||
do p=1,2
|
||||
p2 = 1+mod(p,2)
|
||||
nact = elec_num_2(p) -mo_closed_num
|
||||
nact2 = elec_num_2(p2)-mo_closed_num
|
||||
if ( exc(3) == 0 ) then
|
||||
eplf_factor_num_max += 2*nact*mo_num
|
||||
else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
|
||||
eplf_factor_num_max += 2*mo_num
|
||||
else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
|
||||
eplf_factor_num_max += 2
|
||||
else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
|
||||
eplf_factor_num_max += 1
|
||||
endif
|
||||
enddo
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, eplf_factor_indice, (4,eplf_factor_num_max) ]
|
||||
&BEGIN_PROVIDER [ real, eplf_factor_value, (2,eplf_factor_num_max) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Compact representation of eplf factors
|
||||
END_DOC
|
||||
|
||||
integer :: i,j,k,l,m
|
||||
|
||||
m=1
|
||||
do i=1,mo_num
|
||||
do j=1,mo_num
|
||||
do k=1,mo_num
|
||||
do l=1,mo_num
|
||||
if ( (eplf_factor(1,l,k,j,i) /= 0.).or. &
|
||||
(eplf_factor(2,l,k,j,i) /= 0.) ) then
|
||||
eplf_factor_indice(1,m) = l
|
||||
eplf_factor_indice(2,m) = k
|
||||
eplf_factor_indice(3,m) = j
|
||||
eplf_factor_indice(4,m) = i
|
||||
eplf_factor_value(1,m) = eplf_factor(1,l,k,j,i)
|
||||
eplf_factor_value(2,m) = eplf_factor(2,l,k,j,i)
|
||||
m += 1
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
FREE eplf_factor
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, eplf_factor, (2,mo_num,mo_num,mo_num,mo_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Factors containing the Slater rules
|
||||
END_DOC
|
||||
|
||||
integer :: i, j
|
||||
|
||||
integer :: k,l,m,n,p,p2
|
||||
integer :: ik,il,jk,jl
|
||||
real :: phase
|
||||
integer :: exc(4), nact, nact2
|
||||
real :: det_kl
|
||||
integer :: det_exc
|
||||
|
||||
do m=1,2
|
||||
do i=1,mo_num
|
||||
do j=1,mo_num
|
||||
do k=1,mo_num
|
||||
do l=1,mo_num
|
||||
eplf_factor(m,l,k,j,i) = 0.
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
PROVIDE det
|
||||
|
||||
do k=1,det_num
|
||||
do l=k,det_num
|
||||
|
||||
exc(1) = det_exc(k,l,1)
|
||||
exc(2) = det_exc(k,l,2)
|
||||
exc(4) = exc(1)*exc(2)
|
||||
exc(1) = abs(exc(1))
|
||||
exc(2) = abs(exc(2))
|
||||
exc(3) = exc(1)+exc(2)
|
||||
if (exc(4) /= 0) then
|
||||
exc(4) = exc(4)/abs(exc(4))
|
||||
else
|
||||
exc(4) = 1
|
||||
endif
|
||||
phase = dble(exc(4))
|
||||
|
||||
det_kl = phase*det_coef(k)*det_coef(l)
|
||||
if (k /= l) then
|
||||
det_kl += det_kl
|
||||
endif
|
||||
|
||||
do p=1,2
|
||||
p2 = 1+mod(p,2)
|
||||
nact = elec_num_2(p) -mo_closed_num
|
||||
nact2 = elec_num_2(p2)-mo_closed_num
|
||||
if ( exc(3) == 0 ) then
|
||||
do n=1,nact
|
||||
jk = det(n,k,p)
|
||||
jl = det(n,l,p)
|
||||
|
||||
do i=1,mo_closed_num
|
||||
! Closed-open shell interactions
|
||||
eplf_factor(1,jk,jl,i,i) += det_kl
|
||||
eplf_factor(2,jk,jl,i,i) += det_kl
|
||||
eplf_factor(1,i,jl,jk,i) -= det_kl
|
||||
|
||||
!- Open-closed shell interactions
|
||||
eplf_factor(1,i,i,jk,jl) += det_kl
|
||||
eplf_factor(2,i,i,jk,jl) += det_kl
|
||||
eplf_factor(1,jk,i,i,jl) -= det_kl
|
||||
enddo
|
||||
|
||||
!- Open-open shell interactions
|
||||
do m=1,nact
|
||||
ik = det(m,k,p)
|
||||
il = det(m,l,p)
|
||||
eplf_factor(1,ik,il,jk,jl) += det_kl
|
||||
eplf_factor(1,jk,il,ik,jl) -= det_kl
|
||||
enddo
|
||||
do m=1,nact2
|
||||
ik = det(m,k,p2)
|
||||
il = det(m,l,p2)
|
||||
eplf_factor(2,ik,il,jk,jl) += det_kl
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
else if ( (exc(3) == 1).and.(exc(p) == 1) ) then
|
||||
|
||||
! Sum over only the sigma-sigma interactions involving the excitation
|
||||
call get_single_excitation(k,l,ik,il,p)
|
||||
|
||||
do i=1,mo_closed_num
|
||||
!- Open-closed shell interactions
|
||||
eplf_factor(1,ik,il,i,i) += det_kl
|
||||
eplf_factor(2,ik,il,i,i) += det_kl
|
||||
eplf_factor(1,i,il,ik,i) -= det_kl
|
||||
|
||||
!- Closed-open shell interactions
|
||||
eplf_factor(1,i,i,jk,jl) += det_kl
|
||||
eplf_factor(2,i,i,jk,jl) += det_kl
|
||||
eplf_factor(1,jk,i,i,jl) -= det_kl
|
||||
enddo
|
||||
|
||||
!- Open-open shell interactions
|
||||
do m=1,nact
|
||||
jk = det(m,k,p)
|
||||
jl = det(m,l,p)
|
||||
eplf_factor(1,ik,il,jk,jl) += det_kl
|
||||
eplf_factor(1,jk,il,ik,jl) -= det_kl
|
||||
enddo
|
||||
do m=1,nact2
|
||||
jk = det(m,k,p2)
|
||||
jl = det(m,l,p2)
|
||||
eplf_factor(2,ik,il,jk,jl) += det_kl
|
||||
enddo
|
||||
|
||||
else if ( (exc(3) == 2).and.(exc(p) == 2) ) then
|
||||
|
||||
! Consider only the double excitations of same-spin electrons
|
||||
call get_double_excitation(k,l,ik,il,jk,jl,p)
|
||||
eplf_factor(1,ik,il,jk,jl) += det_kl
|
||||
eplf_factor(1,jk,il,ik,jl) -= det_kl
|
||||
|
||||
else if ( (exc(3) == 2).and.(exc(p) == 1) ) then
|
||||
|
||||
! Consider only the double excitations of opposite-spin electrons
|
||||
call get_single_excitation(k,l,ik,il,p)
|
||||
call get_single_excitation(k,l,jk,jl,p2)
|
||||
eplf_factor(2,ik,il,jk,jl) += det_kl
|
||||
|
||||
endif
|
||||
enddo
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ real, eplf_value_p ]
|
||||
implicit none
|
||||
|
|
|
|||
Loading…
Reference in New Issue