mirror of https://gitlab.com/scemama/eplf
576 lines
15 KiB
Fortran
576 lines
15 KiB
Fortran
BEGIN_PROVIDER [ real, eplf_gamma ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Value of the gaussian for the EPLF
|
|
END_DOC
|
|
include 'constants.F'
|
|
real :: eps, N
|
|
N = 0.1
|
|
eps = -real(dlog(tiny(1.d0)))
|
|
eplf_gamma = (4./3.*pi*density_value_p/N)**(2./3.) * eps
|
|
!eplf_gamma = 1.e10
|
|
!eplf_gamma = 1.e5
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF
|
|
END_DOC
|
|
integer :: i, j
|
|
double precision :: ao_eplf_integral
|
|
do i=1,ao_num
|
|
do j=i,ao_num
|
|
ao_eplf_integral_matrix(j,i) = ao_eplf_integral(j,i,eplf_gamma,point)
|
|
ao_eplf_integral_matrix(i,j) = ao_eplf_integral_matrix(j,i)
|
|
enddo
|
|
enddo
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Array of all the <phi_i phi_j | exp(-gamma r^2)> for EPLF
|
|
END_DOC
|
|
integer :: i, j, k, l
|
|
double precision :: t
|
|
PROVIDE ao_eplf_integral_matrix
|
|
PROVIDE mo_coef
|
|
do i=1,mo_num
|
|
do j=i,mo_num
|
|
mo_eplf_integral_matrix(j,i) = 0.d0
|
|
enddo
|
|
|
|
|
|
do k=1,ao_num
|
|
if (mo_coef(k,i) /= 0.) then
|
|
do l=1,ao_num
|
|
if (abs(ao_eplf_integral_matrix(l,k))>1.d-16) then
|
|
t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k)
|
|
do j=i,mo_num
|
|
mo_eplf_integral_matrix(j,i) += t*mo_coef_transp(j,l)
|
|
enddo
|
|
endif
|
|
enddo
|
|
endif
|
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
do i=1,mo_num
|
|
do j=i+1,mo_num
|
|
mo_eplf_integral_matrix(i,j) = mo_eplf_integral_matrix(j,i)
|
|
enddo
|
|
enddo
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, eplf_up_up ]
|
|
&BEGIN_PROVIDER [ double precision, eplf_up_dn ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Value of the d_upup and d_updn quantities needed for EPLF
|
|
END_DOC
|
|
|
|
integer :: i, j
|
|
|
|
eplf_up_up = 0.d0
|
|
eplf_up_dn = 0.d0
|
|
|
|
PROVIDE mo_coef_transp
|
|
|
|
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)
|
|
enddo
|
|
enddo
|
|
eplf_up_up *= 2.d0
|
|
eplf_up_dn *= 2.d0
|
|
|
|
integer :: k,l,m,n,p,p2
|
|
integer :: ik,il,jk,jl
|
|
double precision :: phase,dtemp(2)
|
|
integer :: exc(4)
|
|
|
|
PROVIDE det
|
|
PROVIDE elec_num_2
|
|
|
|
do k=1,det_num
|
|
do l=1,det_num
|
|
|
|
exc(1) = abs(det_exc(k,l,1))
|
|
exc(2) = abs(det_exc(k,l,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
|
|
|
|
dtemp(1) = 0.d0
|
|
dtemp(2) = 0.d0
|
|
do p=1,2
|
|
p2 = 1+mod(p,2)
|
|
if ( exc(3) == 0 ) then
|
|
! Closed-open shell interactions
|
|
do j=1,mo_closed_num
|
|
do n=1,elec_num_2(p)-mo_closed_num
|
|
ik = det(n,k,p)
|
|
il = det(n,l,p)
|
|
dtemp(1) += mo_value_p(ik)* ( &
|
|
mo_value_p(il)*mo_eplf_integral_matrix(j,j) - &
|
|
mo_value_p(j)*mo_eplf_integral_matrix(j,il) )
|
|
dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(j,j)
|
|
enddo
|
|
enddo
|
|
|
|
!- Open-closed shell interactions
|
|
do m=1,elec_num_2(p)-mo_closed_num
|
|
jk = det(m,k,p)
|
|
jl = det(m,l,p)
|
|
do i=1,mo_closed_num
|
|
dtemp(1) += mo_value_p(i)* ( &
|
|
mo_value_p(i)*mo_eplf_integral_matrix(jk,jl) - &
|
|
mo_value_p(jl)*mo_eplf_integral_matrix(jk,i) )
|
|
dtemp(2) += mo_value_p(i)*mo_value_p(i)*mo_eplf_integral_matrix(jk,jl)
|
|
enddo
|
|
enddo
|
|
|
|
!- Open-open shell interactions
|
|
do m=1,elec_num_2(p)-mo_closed_num
|
|
jk = det(m,k,p)
|
|
jl = det(m,l,p)
|
|
do n=1,elec_num_2(p)-mo_closed_num
|
|
ik = det(n,k,p)
|
|
il = det(n,l,p)
|
|
dtemp(1) += mo_value_p(ik)* ( &
|
|
mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
|
|
mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
|
|
enddo
|
|
do n=1,elec_num_2(p2)-mo_closed_num
|
|
ik = det(n,k,p2)
|
|
il = det(n,l,p2)
|
|
dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(jk,jl)
|
|
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)
|
|
|
|
!- Open-closed shell interactions
|
|
do j=1,mo_closed_num
|
|
dtemp(1) += mo_value_p(ik)* ( &
|
|
mo_value_p(il)*mo_eplf_integral_matrix(j,j) - &
|
|
mo_value_p(j)*mo_eplf_integral_matrix(j,il) )
|
|
dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(j,j)
|
|
enddo
|
|
|
|
!- Closed-open shell interactions
|
|
do i=1,mo_closed_num
|
|
dtemp(1) += mo_value_p(i)* ( &
|
|
mo_value_p(i)*mo_eplf_integral_matrix(jk,jl) - &
|
|
mo_value_p(jl)*mo_eplf_integral_matrix(jk,i) )
|
|
dtemp(2) += mo_value_p(i)*mo_value_p(i)*mo_eplf_integral_matrix(jk,jl)
|
|
enddo
|
|
|
|
!- Open-open shell interactions
|
|
do m=1,elec_num_2(p)-mo_closed_num
|
|
jk = det(m,k,p)
|
|
jl = det(m,l,p)
|
|
dtemp(1) += mo_value_p(ik)* ( &
|
|
mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
|
|
mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
|
|
enddo
|
|
do m=1,elec_num_2(p2)-mo_closed_num
|
|
jk = det(m,k,p2)
|
|
jl = det(m,l,p2)
|
|
dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(jk,jl)
|
|
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_p(ik)* ( &
|
|
mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
|
|
mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
|
|
|
|
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_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(jk,jl)
|
|
|
|
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)
|
|
|
|
enddo
|
|
enddo
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ real, eplf_value_p ]
|
|
implicit none
|
|
BEGIN_DOC
|
|
! Value of the EPLF at the current point.
|
|
END_DOC
|
|
double precision :: aa, ab
|
|
double precision, parameter :: eps = tiny(1.d0)
|
|
|
|
aa = eplf_up_up
|
|
ab = eplf_up_dn
|
|
if ( (aa > 0.d0).and.(ab > 0.d0) ) then
|
|
aa = min(1.d0,aa)
|
|
ab = min(1.d0,ab)
|
|
aa = -dlog(aa)/eplf_gamma
|
|
ab = -dlog(ab)/eplf_gamma
|
|
aa = dsqrt(aa)
|
|
ab = dsqrt(ab)
|
|
eplf_value_p = (aa-ab)/(aa+ab+eps)
|
|
else
|
|
eplf_value_p = 0.d0
|
|
endif
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
double precision function ao_eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,gmma,xr)
|
|
implicit none
|
|
include 'constants.F'
|
|
|
|
real, intent(in) :: a,b,gmma ! Exponents
|
|
real, intent(in) :: xa,xb,xr ! Centers
|
|
integer, intent(in) :: i,j ! Powers of xa and xb
|
|
integer,parameter :: Npoints=10000
|
|
real :: x, xmin, xmax, dx
|
|
|
|
ASSERT (a>0.)
|
|
ASSERT (b>0.)
|
|
ASSERT (i>=0)
|
|
ASSERT (j>=0)
|
|
|
|
xmin = min(xa,xb)
|
|
xmax = max(xa,xb)
|
|
xmin = min(xmin,xr) - 10.
|
|
xmax = max(xmax,xr) + 10.
|
|
dx = (xmax-xmin)/real(Npoints)
|
|
|
|
real :: dtemp
|
|
dtemp = 0.
|
|
x = xmin
|
|
integer :: k
|
|
do k=1,Npoints
|
|
dtemp += &
|
|
(x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2))
|
|
x = x+dx
|
|
enddo
|
|
ao_eplf_integral_primitive_oneD_numeric = dtemp*dx
|
|
|
|
end function
|
|
|
|
double precision function ao_eplf_integral_numeric(i,j,gmma,center)
|
|
implicit none
|
|
integer, intent(in) :: i, j
|
|
integer :: p,q,k
|
|
double precision :: integral
|
|
double precision :: ao_eplf_integral_primitive_oneD_numeric
|
|
real :: gmma, center(3), c
|
|
|
|
|
|
ao_eplf_integral_numeric = 0.d0
|
|
do q=1,ao_prim_num(j)
|
|
do p=1,ao_prim_num(i)
|
|
c = ao_coef(i,p)*ao_coef(j,q)
|
|
integral = &
|
|
ao_eplf_integral_primitive_oneD_numeric( &
|
|
ao_expo(i,p), &
|
|
nucl_coord(ao_nucl(i),1), &
|
|
ao_power(i,1), &
|
|
ao_expo(j,q), &
|
|
nucl_coord(ao_nucl(j),1), &
|
|
ao_power(j,1), &
|
|
gmma, &
|
|
center(1)) * &
|
|
ao_eplf_integral_primitive_oneD_numeric( &
|
|
ao_expo(i,p), &
|
|
nucl_coord(ao_nucl(i),2), &
|
|
ao_power(i,2), &
|
|
ao_expo(j,q), &
|
|
nucl_coord(ao_nucl(j),2), &
|
|
ao_power(j,2), &
|
|
gmma, &
|
|
center(2)) * &
|
|
ao_eplf_integral_primitive_oneD_numeric( &
|
|
ao_expo(i,p), &
|
|
nucl_coord(ao_nucl(i),3), &
|
|
ao_power(i,3), &
|
|
ao_expo(j,q), &
|
|
nucl_coord(ao_nucl(j),3), &
|
|
ao_power(j,3), &
|
|
gmma, &
|
|
center(3))
|
|
ao_eplf_integral_numeric = ao_eplf_integral_numeric + c*integral
|
|
enddo
|
|
enddo
|
|
|
|
end function
|
|
|
|
|
|
double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
|
implicit none
|
|
include 'constants.F'
|
|
|
|
real, intent(in) :: a,b,gmma ! Exponents
|
|
real, intent(in) :: xa,xb,xr ! Centers
|
|
integer, intent(in) :: i,j ! Powers of xa and xb
|
|
integer :: ii, jj, kk, ll
|
|
real :: xp1,xp
|
|
real :: p1,p
|
|
double precision :: S(0:i+1,0:j+1)
|
|
double precision :: inv_p(2), di(max(i,j)), dj(j), c
|
|
|
|
ASSERT (a>0.)
|
|
ASSERT (b>0.)
|
|
ASSERT (i>=0)
|
|
ASSERT (j>=0)
|
|
|
|
! Gaussian product
|
|
! Inlined Gaussian products (same as call gaussian_product)
|
|
real :: t(2), xab(2), ab(2)
|
|
inv_p(1) = 1.d0/(a+b)
|
|
p1 = a+b
|
|
ab(1) = a*b
|
|
inv_p(2) = 1.d0/(p1+gmma)
|
|
t(1) = (a*xa+b*xb)
|
|
xab(1) = xa-xb
|
|
xp1 = t(1)*inv_p(1)
|
|
p = p1+gmma
|
|
ab(2) = p1*gmma
|
|
t(2) = (p1*xp1+gmma*xr)
|
|
xab(2) = xp1-xr
|
|
xp = t(2)*inv_p(2)
|
|
c = real(ab(1)*inv_p(1)*xab(1)**2 + &
|
|
ab(2)*inv_p(2)*xab(2)**2)
|
|
if ( c > 32.d0 ) then
|
|
ao_eplf_integral_primitive_oneD = 0.d0
|
|
return
|
|
endif
|
|
c = exp(-c)
|
|
!S(0,0) = dsqrt(pi*inv_p(2))*c
|
|
S(0,0) = 1.d0 ! Factor is applied at the end
|
|
|
|
! Obara-Saika recursion
|
|
|
|
do ii=1,max(i,j)
|
|
di(ii) = 0.5d0*inv_p(2)*dble(ii)
|
|
enddo
|
|
|
|
xab(1) = xp-xa
|
|
xab(2) = xp-xb
|
|
S(1,0) = xab(1) * S(0,0)
|
|
if (i>1) then
|
|
do ii=1,i-1
|
|
S(ii+1,0) = xab(1) * S(ii,0) + di(ii)*S(ii-1,0)
|
|
enddo
|
|
endif
|
|
|
|
S(0,1) = xab(2) * S(0,0)
|
|
if (j>1) then
|
|
do jj=1,j-1
|
|
S(0,jj+1) = xab(2) * S(0,jj) + di(jj)*S(0,jj-1)
|
|
enddo
|
|
endif
|
|
|
|
do jj=1,j
|
|
S(1,jj) = xab(1) * S(0,jj) + di(jj) * S(0,jj-1)
|
|
do ii=2,i
|
|
S(ii,jj) = xab(1) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
|
|
enddo
|
|
enddo
|
|
|
|
ao_eplf_integral_primitive_oneD = dsqrt(pi*inv_p(2))*S(i,j)*c ! Application of the factor of S(0,0)
|
|
|
|
end function
|
|
|
|
!double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
|
! implicit none
|
|
! include 'constants.F'
|
|
!!
|
|
! real, intent(in) :: a,b,gmma ! Exponents
|
|
! real, intent(in) :: xa,xb,xr ! Centers
|
|
! integer, intent(in) :: i,j ! Powers of xa and xb
|
|
! integer :: ii, jj, kk, ll
|
|
! real :: xp1,xp
|
|
! real :: p1,p
|
|
! double precision :: xpa, xpb
|
|
! double precision :: inv_p(2),S00, c
|
|
! double precision :: ObaraS
|
|
!!
|
|
! ASSERT (a>0.)
|
|
! ASSERT (b>0.)
|
|
! ASSERT (i>=0)
|
|
! ASSERT (j>=0)
|
|
!!
|
|
! ! Inlined Gaussian products (same as call gaussian_product)
|
|
! real :: t(2), xab(2), ab(2)
|
|
! inv_p(1) = 1.d0/(a+b)
|
|
! p1 = a+b
|
|
! ab(1) = a*b
|
|
! inv_p(2) = 1.d0/(p1+gmma)
|
|
! t(1) = (a*xa+b*xb)
|
|
! xab(1) = xa-xb
|
|
! xp1 = t(1)*inv_p(1)
|
|
! p = p1+gmma
|
|
! ab(2) = p1*gmma
|
|
! t(2) = (p1*xp1+gmma*xr)
|
|
! xab(2) = xp1-xr
|
|
! xp = t(2)*inv_p(2)
|
|
! c = exp(- real(ab(1)*inv_p(1)*xab(1)**2 + &
|
|
! ab(2)*inv_p(2)*xab(2)**2) )
|
|
!!
|
|
! xpa = xp-xa
|
|
! xpb = xp-xb
|
|
! S00 = sqrt(real(pi*inv_p(2)))*c
|
|
! ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p(2),S00)
|
|
!!
|
|
!end function
|
|
!!
|
|
!recursive double precision function ObaraS(i,j,xpa,xpb,inv_p,S00) result(res)
|
|
! implicit none
|
|
! integer, intent(in) :: i, j
|
|
! double precision, intent(in) :: xpa, xpb, inv_p
|
|
! double precision,intent(in) :: S00
|
|
!!
|
|
! if (i == 0) then
|
|
! if (j == 0) then
|
|
! res = S00
|
|
! else ! (j>0)
|
|
! res = xpb*ObaraS(0,j-1,xpa,xpb,inv_p,S00)
|
|
! if (j>1) then
|
|
! res += 0.5d0*dble(j-1)*inv_p*ObaraS(0,j-2,xpa,xpb,inv_p,S00)
|
|
! endif
|
|
! endif ! (i==0).and.(j>0)
|
|
! else ! (i>0)
|
|
! if (j==0) then
|
|
! res = xpa*ObaraS(i-1,0,xpa,xpb,inv_p,S00)
|
|
! if (i>1) then
|
|
! res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,0,xpa,xpb,inv_p,S00)
|
|
! endif
|
|
! else ! (i>0).and.(j>0)
|
|
! res = xpa * ObaraS(i-1,j,xpa,xpb,inv_p,S00)
|
|
! if (i>1) then
|
|
! res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00)
|
|
! endif
|
|
! res += 0.5d0*dble(j)*inv_p*ObaraS(i-1,j-1,xpa,xpb,inv_p,S00)
|
|
! endif ! (i>0).and.(j>0)
|
|
! endif ! (i>0)
|
|
!!
|
|
!end function
|
|
|
|
double precision function ao_eplf_integral(i,j,gmma,center)
|
|
implicit none
|
|
integer, intent(in) :: i, j
|
|
integer :: p,q,k
|
|
double precision :: integral
|
|
!DEC$ ATTRIBUTES FORCEINLINE
|
|
double precision :: ao_eplf_integral_primitive_oneD
|
|
real :: gmma, center(3)
|
|
double precision :: buffer(100)
|
|
|
|
|
|
ASSERT(i>0)
|
|
ASSERT(j>0)
|
|
ASSERT(i<=ao_num)
|
|
ASSERT(j<=ao_num)
|
|
|
|
ao_eplf_integral = 0.d0
|
|
do p=1,ao_prim_num_max
|
|
buffer(p) = 0.d0
|
|
enddo
|
|
|
|
do q=1,ao_prim_num(j)
|
|
do p=1,ao_prim_num(i)
|
|
integral = &
|
|
ao_eplf_integral_primitive_oneD( &
|
|
ao_expo(i,p), &
|
|
nucl_coord(ao_nucl(i),1), &
|
|
ao_power(i,1), &
|
|
ao_expo(j,q), &
|
|
nucl_coord(ao_nucl(j),1), &
|
|
ao_power(j,1), &
|
|
gmma, &
|
|
center(1))
|
|
if (integral /= 0.d0) then
|
|
integral *= &
|
|
ao_eplf_integral_primitive_oneD( &
|
|
ao_expo(i,p), &
|
|
nucl_coord(ao_nucl(i),2), &
|
|
ao_power(i,2), &
|
|
ao_expo(j,q), &
|
|
nucl_coord(ao_nucl(j),2), &
|
|
ao_power(j,2), &
|
|
gmma, &
|
|
center(2))
|
|
if (integral /= 0.d0) then
|
|
integral *= &
|
|
ao_eplf_integral_primitive_oneD( &
|
|
ao_expo(i,p), &
|
|
nucl_coord(ao_nucl(i),3), &
|
|
ao_power(i,3), &
|
|
ao_expo(j,q), &
|
|
nucl_coord(ao_nucl(j),3), &
|
|
ao_power(j,3), &
|
|
gmma, &
|
|
center(3))
|
|
buffer(p) += integral*ao_coef(i,p)*ao_coef(j,q)
|
|
endif
|
|
endif
|
|
enddo
|
|
enddo
|
|
do p=1,ao_prim_num_max
|
|
ao_eplf_integral += buffer(p)
|
|
enddo
|
|
|
|
end function
|
|
|
|
double precision function mo_eplf_integral(i,j)
|
|
implicit none
|
|
integer :: i, j, k, l
|
|
PROVIDE ao_eplf_integral_matrix
|
|
PROVIDE mo_coef
|
|
mo_eplf_integral = 0.d0
|
|
|
|
do k=1,ao_num
|
|
if (mo_coef(k,i) /= 0.) then
|
|
do l=1,ao_num
|
|
mo_eplf_integral += &
|
|
mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(k,l)
|
|
enddo
|
|
endif
|
|
enddo
|
|
|
|
end function
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|