LCPQ
/
eplf
mirror of https://gitlab.com/scemama/eplf
10
0
Fork 0
Code Issues Releases Wiki Activity

Acceleration of multi-det

This commit is contained in:
Anthony Scemama 2011-02-11 09:11:15 +01:00
parent 8d290ed264
commit 15082646ca
6 changed files with 269 additions and 127 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
bin/to_ezfio.exe
View File

Binary file not shown.

173
src/det.irp.f
View File

@ -32,80 +32,70 @@ BEGIN_PROVIDER [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ]
END_PROVIDER
BEGIN_PROVIDER [ integer*1, det_exc, (det_num, det_num, 2) ]
BEGIN_PROVIDER [ real, mo_occ, (mo_tot_num) ]
implicit none
BEGIN_DOC
! Degree of excitation between two determinants. Indices are alpha, beta
! The sign is the phase factor
BEGIN_DOC
! Occupation numbers of molecular orbitals
END_DOC
integer :: p
do p=1,2
call get_mo_basis_mo_occ(mo_occ)
integer :: k, l
do l=1,det_num
det_exc(l,l,p) = 0
do k=l+1,det_num
det_exc(k,l,p) = 0
END_PROVIDER
! Excitation degree
integer :: i, j
do i=1,elec_num_2(p)-mo_closed_num
logical :: found
found = .False.
do j=1,elec_num_2(p)-mo_closed_num
if (det(j,l,p) == det(i,k,p)) then
found = .True.
exit
endif
enddo
if (.not.found) then
det_exc(k,l,p) += 1
endif
enddo
integer function det_exc(k,l,p)
implicit none
! Degree of excitation between two determinants. Indices are alpha, beta
! The sign is the phase factor
enddo
integer :: k,l,p
integer :: i, j
det_exc = 0
do i=1,elec_num_2(p)-mo_closed_num
logical :: found
found = .False.
do j=1,elec_num_2(p)-mo_closed_num
if (det(j,l,p) == det(i,k,p)) then
found = .True.
! exit
endif
enddo
if (.not.found) then
det_exc += 1
endif
enddo
! Phase
do l=1,det_num
do k=l+1,det_num
integer :: nperm
nperm = 0
do p=1,2
integer :: buffer(0:mo_num-mo_closed_num)
do i=1,elec_num_2(p)-mo_closed_num
buffer(i) = det(i,k,p)
enddo
do i=1,elec_num_2(p)-mo_closed_num
if (buffer(i) /= det(i,l,p)) then
integer :: m
m=elec_num_2(p)-mo_closed_num
do j=i+1,elec_num_2(p)-mo_closed_num
if (buffer(i) == det(j,l,p)) then ! found
m=j
exit
endif
enddo
buffer(0) = buffer(i)
buffer(i) = det(m,l,p)
buffer(m) = buffer(0)
nperm += 1
endif
enddo
det_exc(k,l,p) *= (1-2*mod( nperm, 2 ))
det_exc(l,k,p) = det_exc(k,l,p)
enddo
enddo
integer :: nperm
nperm = 0
integer :: buffer(0:mo_num-mo_closed_num)
do i=1,elec_num_2(p)-mo_closed_num
buffer(i) = det(i,k,p)
enddo
do i=1,elec_num_2(p)-mo_closed_num
if (buffer(i) /= det(i,l,p)) then
integer :: m
m=elec_num_2(p)-mo_closed_num
do j=i+1,elec_num_2(p)-mo_closed_num
if (buffer(i) == det(j,l,p)) then ! found
m=j
exit
endif
enddo
buffer(0) = buffer(i)
buffer(i) = det(m,l,p)
buffer(m) = buffer(0)
nperm += 1
endif
enddo
det_exc *= (1-2*mod( nperm, 2 ))
END_PROVIDER
end
subroutine get_single_excitation(k,l,m,n,p)
implicit none
@ -198,3 +188,70 @@ subroutine get_double_excitation(k,l,m,n,r,s,p)
end
BEGIN_PROVIDER [ real, ci_mo, (mo_num,mo_num,3) ]
implicit none
BEGIN_DOC
! Spin Density matrix in the AO basis
END_DOC
integer :: i,j,k,l,m,ispin, ik,il
do ispin=1,3
do j=1,mo_num
do i=1,mo_num
ci_mo(i,j,ispin) = 0.
enddo
enddo
do l=1,det_num
do m=1,det_num
real :: factor
factor = 2.*det_coef(l)*det_coef(m)
do il=1,mo_closed_num
do ik=1,mo_closed_num
ci_mo(ik,il,ispin) += factor
enddo
enddo
enddo
enddo
enddo
do l=1,det_num
do m=1,det_num
factor = det_coef(l)*det_coef(m)
do ispin=1,2
do j=mo_closed_num+1,elec_num_2(ispin)
ik = det(j-mo_closed_num,l,ispin)
do il=1,mo_closed_num
ci_mo(ik,il,ispin) += factor
ci_mo(il,ik,ispin) += factor
enddo
do i=mo_closed_num+1,elec_num_2(ispin)
il = det(i-mo_closed_num,m,ispin)
ci_mo(ik,il,ispin) += factor
ci_mo(il,ik,ispin) += factor
enddo
enddo
enddo
ispin=3
do j=mo_closed_num+1,elec_num_2(1)
ik = det(j-mo_closed_num,l,1)
do il=1,mo_closed_num
ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m)
ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m)
enddo
do i=mo_closed_num+1,elec_num_2(2)
il = det(i-mo_closed_num,m,2)
ci_mo(ik,il,ispin) += det_coef(l)*det_coef(m)
ci_mo(il,ik,ispin) += det_coef(l)*det_coef(m)
enddo
enddo
enddo
enddo
END_PROVIDER

203
src/eplf_function.irp.f
View File

@ -68,6 +68,7 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, eplf_up_up ]
&BEGIN_PROVIDER [ double precision, eplf_up_dn ]
implicit none
@ -84,40 +85,151 @@ END_PROVIDER
do i=1,mo_closed_num
do j=1,mo_closed_num
eplf_up_up += mo_value_p(i)* ( &
mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
mo_value_p(j)*mo_eplf_integral_matrix(j,i) )
eplf_up_dn += mo_value_p(i)*mo_value_p(i)* &
mo_eplf_integral_matrix(j,j)
double precision :: temp
temp = mo_value_prod_p(i,i)*mo_eplf_integral_matrix(j,j)
eplf_up_up += temp - mo_value_prod_p(j,i)*mo_eplf_integral_matrix(j,i)
eplf_up_dn += temp
enddo
enddo
eplf_up_up *= 2.d0
eplf_up_dn *= 2.d0
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)
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
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
double precision :: phase,dtemp(2)
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
PROVIDE elec_num_2
PROVIDE mo_value_prod_p
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(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)
exc(4) = det_exc(k,l,1)*det_exc(k,l,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
dtemp(1) = 0.d0
dtemp(2) = 0.d0
do p=1,2
p2 = 1+mod(p,2)
nact = elec_num_2(p) -mo_closed_num
@ -129,30 +241,27 @@ END_PROVIDER
do i=1,mo_closed_num
! Closed-open shell interactions
dtemp(1) += ( &
mo_value_prod_p(jl,jk)*mo_eplf_integral_matrix(i,i) - &
mo_value_prod_p(i,jk)*mo_eplf_integral_matrix(i,jl) )
dtemp(2) += mo_value_prod_p(jl,jk)*mo_eplf_integral_matrix(i,i)
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
dtemp(1) += ( &
mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk) - &
mo_value_prod_p(i,jl)*mo_eplf_integral_matrix(i,jk) )
dtemp(2) += mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk)
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)
dtemp(1) += ( &
mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(jl,jk) - &
mo_value_prod_p(jl,ik)*mo_eplf_integral_matrix(il,jk) )
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)
dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(jl,jk)
eplf_factor(2,ik,il,jk,jl) += det_kl
enddo
enddo
@ -164,60 +273,46 @@ END_PROVIDER
do i=1,mo_closed_num
!- Open-closed shell interactions
dtemp(1) += ( &
mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(i,i) - &
mo_value_prod_p(i,ik)*mo_eplf_integral_matrix(i,il) )
dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(i,i)
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
dtemp(1) += ( &
mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk) - &
mo_value_prod_p(i,jl)*mo_eplf_integral_matrix(i,jk) )
dtemp(2) += mo_value_prod_p(i,i)*mo_eplf_integral_matrix(jl,jk)
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)
dtemp(1) += ( &
mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(jl,jk) - &
mo_value_prod_p(jl,ik)*mo_eplf_integral_matrix(il,jk) )
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)
dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(jl,jk)
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)
dtemp(1) += ( &
mo_value_prod_p(il,ik)*mo_eplf_integral_matrix(jl,jk) - &
mo_value_prod_p(jl,ik)*mo_eplf_integral_matrix(il,jk) )
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)
dtemp(2) += mo_value_prod_p(ik,il)*mo_eplf_integral_matrix(jl,jk)
! 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
phase = dble(exc(4))
eplf_up_up += phase * det_coef(k)*det_coef(l) * dtemp(1)
eplf_up_dn += phase * det_coef(k)*det_coef(l) * dtemp(2)
if (k /= l) then
eplf_up_up += phase * det_coef(k)*det_coef(l) * dtemp(1)
eplf_up_dn += phase * det_coef(k)*det_coef(l) * dtemp(2)
endif
enddo
enddo

10
src/mo.irp.f
View File

@ -43,16 +43,6 @@ BEGIN_PROVIDER [ integer, mo_num ]
END_PROVIDER
BEGIN_PROVIDER [ real, mo_occ, (mo_tot_num) ]
implicit none
BEGIN_DOC
! Occupation numbers of molecular orbitals
END_DOC
call get_mo_basis_mo_occ(mo_occ)
END_PROVIDER
BEGIN_PROVIDER [ real, mo_coef, (ao_num,mo_num) ]
implicit none
BEGIN_DOC

8
src/overlap.irp.f
View File

@ -155,10 +155,10 @@ double precision function primitive_overlap_oneD(a,xa,i,b,xb,j)
xp = xp*inv_p
c = a*b*inv_p*(xa-xb)**2
!if ( c > 32.d0 ) then ! Cut-off on exp(-32)
! primitive_overlap_oneD = 0.d0
! return
!endif
if ( c > 32.d0 ) then ! Cut-off on exp(-32)
primitive_overlap_oneD = 0.d0
return
endif
c = dexp(-c)

2
src/test_1d.irp.f
View File

@ -10,7 +10,7 @@ subroutine run
point(1) = 0.
point(2) = 0.
integer :: i
do i=-40,40
do i=-40,60
point(3) = real(i)/20.
TOUCH point
print *, point(3), eplf_value_p, eplf_up_up, eplf_up_dn