diff --git a/src/Util.irp.f b/src/Util.irp.f index 3907f3c..3093056 100644 --- a/src/Util.irp.f +++ b/src/Util.irp.f @@ -17,37 +17,54 @@ double precision function binom(k,n) -implicit double precision(a-h,o-z) - binom=fact(n)/(fact(k)*fact(n-k)) + implicit none + integer, intent(in) :: k,n + double precision, save :: memo(0:20,0:20) + integer, save :: ifirst + double precision :: fact + if (ifirst == 0) then + ifirst = 1 + integer :: i,j + do j=0,20 + do i=0,i + memo(i,j) = fact(j)/(fact(i)*fact(j-i)) + enddo + enddo + endif + + if (n<20) then + binom = memo(k,n) + else + binom=fact(n)/(fact(k)*fact(n-k)) + endif end double precision function fact2(n) implicit none integer :: n - double precision, save :: memo(1:100) + double precision, save :: memo(100) integer, save :: memomax = 1 ASSERT (mod(n,2) /= 0) - if (n<=memomax) then - if (n<3) then - fact2 = 1.d0 - else - fact2 = memo(n) - endif - return - endif - - integer :: i memo(1) = 1.d0 - do i=memomax+2,min(n,99),2 - memo(i) = memo(i-2)* float(i) - enddo - memomax = min(n,99) - fact2 = memo(memomax) + if (n<2) then + fact2 = 1.d0 + else if (n<=memomax) then + fact2 = memo(n) + else - do i=101,n,2 - fact2 = fact2*float(i) - enddo + integer :: i + do i=memomax+2,min(n,99),2 + memo(i) = memo(i-2)* float(i) + enddo + memomax = min(n,99) + fact2 = memo(memomax) + + do i=101,n,2 + fact2 *= float(i) + enddo + + endif end function @@ -55,28 +72,25 @@ end function double precision function fact(n) implicit none integer :: n - double precision, save :: memo(1:100) + double precision, save :: memo(0:100) integer, save :: memomax = 1 - if (n<=memomax) then - if (n<2) then - fact = 1.d0 - else - fact = memo(n) - endif - return - endif - - integer :: i + memo(0) = 1.d0 memo(1) = 1.d0 - do i=memomax+1,min(n,100) - memo(i) = memo(i-1)*float(i) - enddo - memomax = min(n,100) - fact = memo(memomax) - do i=101,n - fact = fact*float(i) - enddo + + if (n<=memomax) then + fact = memo(n) + else + integer :: i + do i=memomax+1,min(n,100) + memo(i) = memo(i-1)*float(i) + enddo + memomax = min(n,100) + fact = memo(memomax) + do i=101,n + fact *= float(i) + enddo + endif end function diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f index 242faa7..99399e7 100644 --- a/src/eplf_function.irp.f +++ b/src/eplf_function.irp.f @@ -30,33 +30,49 @@ BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ] enddo END_PROVIDER +BEGIN_PROVIDER [ integer, ao_eplf_integral_matrix_non_zero_idx, (0:ao_num,ao_num) ] + implicit none + BEGIN_DOC +! Non zero indices of the ao_eplf_integral_matrix + END_DOC + integer :: i, j + integer :: idx + do j=1,ao_num + idx = 0 + do i=1,ao_num + if (abs(ao_eplf_integral_matrix(i,j)) > 1.e-16) then + idx +=1 + ao_eplf_integral_matrix_non_zero_idx(idx,j) = i + endif + enddo + ao_eplf_integral_matrix_non_zero_idx(0,j) = idx + enddo +END_PROVIDER + BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ] implicit none BEGIN_DOC ! Array of all the for EPLF END_DOC - integer :: i, j, k, l - double precision :: t - PROVIDE ao_eplf_integral_matrix - PROVIDE mo_coef mo_coef_transp + integer :: i, j, k, l, kmo, kao + double precision :: t, moki do i=1,mo_num do j=i,mo_num mo_eplf_integral_matrix(j,i) = 0.d0 enddo + 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 + do k=1,ao_num + do kmo=1,mo_coef_transp_non_zero_idx(0,k) + i = mo_coef_transp_non_zero_idx(kmo,k) + moki = mo_coef_transp(i,k) + do kao = 1,ao_eplf_integral_matrix_non_zero_idx(0,k) + l = ao_eplf_integral_matrix_non_zero_idx(kao,k) + t = moki*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 enddo enddo @@ -66,6 +82,8 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ] mo_eplf_integral_matrix(i,j) = mo_eplf_integral_matrix(j,i) enddo enddo + + END_PROVIDER @@ -231,32 +249,37 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) ASSERT (i>=0) ASSERT (j>=0) - double precision :: alpha, beta, delta, c - alpha = a*xa*xa + b*xb*xb + gmma*xr*xr - delta = a*xa + b*xb + gmma*xr - beta = a + b + gmma - c = alpha-delta**2/beta + double precision :: alpha, beta, delta, c, beta_inv + double precision :: db_inv + + beta = a + b + gmma + beta_inv = 1.d0/beta + delta = a*xa + b*xb + gmma*xr + db_inv = delta*beta_inv + alpha = a*xa*xa + b*xb*xb + gmma*xr*xr + c = alpha-delta*db_inv - if ( c > 32.d0 ) then ! Cut-off on exp(-32) + if ( c > 20.d0 ) then ! Cut-off on exp(-32) ao_eplf_integral_primitive_oneD = 0.d0 return endif double precision :: u,v,w c = exp(-c) - w = 1.d0/sqrt(beta) - u=delta/beta-xa - v=delta/beta-xb + w = sqrt(beta_inv) + u=db_inv-xa + v=db_inv-xb double precision :: fact2, binom, f, ac ac = 0.d0 integer :: istart(0:1) - istart = (/ 0, 1 /) + istart(0) = 0 + istart(1) = 1 do ii=0,i do jj=istart(mod(ii,2)),j,2 - kk=(ii+jj)/2 + kk= -(ii+jj)/2 f=binom(ii,i)*binom(jj,j)*fact2(ii+jj-1) - ac += u**(i-ii)*v**(j-jj)*w**(ii+jj)*f/dble(2**kk) + ac += u**(i-ii) * v**(j-jj) * w**(ii+jj) * f * 2.d0**kk enddo enddo ao_eplf_integral_primitive_oneD = dsqrt_pi*w*c*ac @@ -283,50 +306,50 @@ double precision function ao_eplf_integral(i,j,gmma,center) ASSERT(j<=ao_num) ao_eplf_integral = 0.d0 - do p=1,ao_prim_num_max + do p=1,ao_prim_num(i) 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 + 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 + do p=1,ao_prim_num(i) ao_eplf_integral += buffer(p) enddo diff --git a/src/mo.irp.f b/src/mo.irp.f index 4b6fac9..b5d4753 100644 --- a/src/mo.irp.f +++ b/src/mo.irp.f @@ -78,6 +78,28 @@ BEGIN_PROVIDER [ real, mo_coef_transp, (mo_num,ao_num) ] END_PROVIDER + +BEGIN_PROVIDER [ integer, mo_coef_transp_non_zero_idx, (0:mo_num,ao_num) ] + implicit none + BEGIN_DOC +! Indices of the non-zero elements of the transpose of the Molecular orbital coefficients + END_DOC + integer :: i, j + + integer :: idx + do j=1,ao_num + idx = 0 + do i=1,mo_num + if (mo_coef_transp(i,j) /= 0.) then + idx += 1 + mo_coef_transp_non_zero_idx(idx,j) = i + endif + enddo + mo_coef_transp_non_zero_idx(0,j) = idx + enddo + +END_PROVIDER + BEGIN_PROVIDER [ logical, mo_is_closed, (mo_num) ] &BEGIN_PROVIDER [ logical, mo_is_active, (mo_num) ] implicit none