Accelerated gaussian product

bin/to_ezfio.py

@@ -19,6 +19,7 @@ else:
 firstArg = sys.argv[1]
 
 file = getFile(firstArg)
+file.convert_to_cartesian()
 print firstArg, 'recognized as', str(file).split('.')[-1].split()[0]
 
 from ezfio import ezfio
@@ -52,7 +53,6 @@ def write_ezfioFile(res,filename):
 
 # AO Basis
   import string
-  is_cartesian = True
   at = []
   num_prim = []
   magnetic_number = []
@@ -62,29 +62,18 @@ def write_ezfioFile(res,filename):
   power_z = []
   coefficient = []
   exponent = []
-  for b in res.basis:
-    if '+' in b.sym or '-' in b.sym:
-      is_cartesian = False
-  names = ["s","p","d","f","g","h","i","j"]
   for b in res.basis:
     c = b.center
     for i,atom in enumerate(res.geometry):
       if atom.coord == c:
          at.append(i+1)
     num_prim.append(len(b.prim))
-    if is_cartesian:
-      s = b.sym
-      power_x.append( string.count(s,"x") ) 
-      power_y.append( string.count(s,"y") )
-      power_z.append( string.count(s,"z") )
-    else:
-      magnetic_number.append(names.index(b.sym[0]))
-      angular_number.append(int(b.sym[1:]))
+    s = b.sym
+    power_x.append( string.count(s,"x") ) 
+    power_y.append( string.count(s,"y") )
+    power_z.append( string.count(s,"z") )
     coefficient.append(  b.coef )
     exponent.append( [ p.expo for p in b.prim ] )
-  if not is_cartesian:
-     print 'Only cartesian basis functions work...'
-     sys.exit(0)
   ezfio.ao_basis_ao_num = len(res.basis)
   ezfio.ao_basis_ao_nucl = at
   ezfio.ao_basis_ao_prim_num = num_prim

src/eplf_function.irp.f

@@ -80,13 +80,13 @@ END_PROVIDER
 
  do j=1,elec_beta_num
   do i=1,elec_beta_num
-    eplf_up_up = eplf_up_up + 2.d0*mo_value_p(i)* (  &
+    eplf_up_up += 2.d0*mo_value_p(i)* (  &
       mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
       mo_value_p(j)*mo_eplf_integral_matrix(i,j) )
   enddo
 
   do i=elec_beta_num+1,elec_alpha_num
-    eplf_up_up = eplf_up_up + mo_value_p(i)* (  &
+    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(i,j) )
   enddo
@@ -94,7 +94,7 @@ END_PROVIDER
 
  do j=elec_beta_num+1,elec_alpha_num
   do i=1,elec_alpha_num
-    eplf_up_up = eplf_up_up + mo_value_p(i)* (  &
+    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(i,j) )
   enddo
@@ -102,7 +102,7 @@ END_PROVIDER
 
  do j=1,elec_beta_num
   do i=1,elec_alpha_num
-    eplf_up_dn = eplf_up_dn + mo_value_p(i)**2 * &
+    eplf_up_dn += mo_value_p(i)**2 * &
       mo_eplf_integral_matrix(j,j) 
   enddo
  enddo
@@ -162,7 +162,7 @@ double precision function ao_eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,
  x = xmin
  integer :: k
  do k=1,Npoints
-  dtemp = dtemp + &
+  dtemp += &
    (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2))
   x = x+dx
  enddo
@@ -294,6 +294,7 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
  ! Gaussian product
  call gaussian_product(a,xa,b,xb,c1,p1,xp1)
  call gaussian_product(p1,xp1,gmma,xr,c,p,xp)
+
  inv_p = 1.d0/p
  S00 = dsqrt(pi*inv_p)*c*c1
  xpa = xp-xa
@@ -314,21 +315,21 @@ recursive double precision function ObaraS(i,j,xpa,xpb,inv_p,S00) result(res)
   else ! (j>0)
    res = xpb*ObaraS(0,j-1,xpa,xpb,inv_p,S00)
    if (j>1) then
-    res = res + 0.5d0*dble(j-1)*inv_p*ObaraS(0,j-2,xpa,xpb,inv_p,S00)
+    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 = res + 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,0,xpa,xpb,inv_p,S00) 
+     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 = res + 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00) 
+    res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00) 
    endif
-   res = res + 0.5d0*dble(j)*inv_p*ObaraS(i-1,j-1,xpa,xpb,inv_p,S00) 
+   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)
 
@@ -378,7 +379,7 @@ double precision function ao_eplf_integral(i,j,gmma,center)
     ao_power(j,3),                         &
     gmma,                                  &
     center(3))
-   ao_eplf_integral = ao_eplf_integral + integral*ao_coef(i,p)*ao_coef(j,q)
+   ao_eplf_integral += integral*ao_coef(i,p)*ao_coef(j,q)
   enddo
  enddo
  
@@ -394,7 +395,7 @@ double precision function mo_eplf_integral(i,j)
    do k=1,ao_num
      if (mo_coef(k,i) /= 0.) then
        do l=1,ao_num
-        mo_eplf_integral = mo_eplf_integral + &
+        mo_eplf_integral +=  &
           mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(k,l)
        enddo
      endif

src/overlap.irp.f

@@ -122,13 +122,15 @@ subroutine gaussian_product(a,xa,b,xb,k,p,xp)
  ASSERT (a>0.)
  ASSERT (b>0.)
 
+ double precision :: t, xab, ab
+
+ p_inv = 1.d0/(a+b)
  p = a+b
- xp = (a*xa+b*xb)
-
- p_inv = 1.d0/p
-
- xp = xp*p_inv
- k = dexp(-a*b*p_inv*(xa-xb)**2)
+ ab = a*b
+ t = (a*xa+b*xb)
+ xab = xa-xb
+ xp = t*p_inv
+ k = dexp(-ab*p_inv*xab**2)
 
 end subroutine