Accelerations

src/Util.irp.f

@@ -17,28 +17,43 @@
 
 
 double precision function binom(k,n)
-implicit double precision(a-h,o-z)
+ 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
+  memo(1) = 1.d0
-    if (n<3) then
+  if (n<2) then
       fact2 = 1.d0
-    else
+  else if (n<=memomax) then
       fact2 = memo(n)
-    endif
+  else
-    return
-  endif
 
     integer :: i
-  memo(1) = 1.d0
     do i=memomax+2,min(n,99),2
       memo(i) = memo(i-2)* float(i)
     enddo
@@ -46,37 +61,36 @@ double precision function fact2(n)
     fact2 = memo(memomax)
 
     do i=101,n,2
-    fact2 = fact2*float(i)
+      fact2 *= float(i)
     enddo
 
+  endif
+
 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
+  memo(0) = 1.d0
-    if (n<2) then
-      fact = 1.d0
-    else
-      fact = memo(n)
-    endif
-    return
-  endif
-
-  integer :: i
   memo(1) = 1.d0
+
+  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 = fact*float(i)
+      fact *= float(i)
     enddo
+  endif
 end function
 
 

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 <phi_i phi_j | exp(-gamma r^2)> for EPLF
  END_DOC
- integer :: i, j, k, l
+ integer :: i, j, k, l, kmo, kao
- double precision :: t
+ double precision :: t, moki
- PROVIDE ao_eplf_integral_matrix
- PROVIDE mo_coef mo_coef_transp
  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 kmo=1,mo_coef_transp_non_zero_idx(0,k)
-    do l=1,ao_num
+   i = mo_coef_transp_non_zero_idx(kmo,k)
-     if (abs(ao_eplf_integral_matrix(l,k))>1.d-16) then
+   moki = mo_coef_transp(i,k)
-      t = mo_coef(k,i)*ao_eplf_integral_matrix(l,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
@@ -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
+  double precision :: alpha, beta, delta, c, beta_inv
-  alpha = a*xa*xa + b*xb*xb + gmma*xr*xr
+  double precision :: db_inv
-  delta = a*xa    + b*xb    + gmma*xr
-  beta  = a       + b       + gmma
-  c = alpha-delta**2/beta
 
-  if ( c > 32.d0 ) then   ! Cut-off on exp(-32)
+  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 > 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)
+  w = sqrt(beta_inv)
-  u=delta/beta-xa
+  u=db_inv-xa
-  v=delta/beta-xb
+  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,7 +306,7 @@ 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
 
@@ -326,7 +349,7 @@ double precision function ao_eplf_integral(i,j,gmma,center)
        endif
     enddo
  enddo
- do p=1,ao_prim_num_max
+ do p=1,ao_prim_num(i)
    ao_eplf_integral += buffer(p)
  enddo
  

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