Slater rules seem to be OK. testing...

src/density.irp.f

@@ -40,21 +40,23 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ]
 ! TODO vectorization
  integer :: k,j,l, ik, il
  real    :: buffer
+ real    :: phase
  PROVIDE det
  PROVIDE elec_alpha_num
  do k=1,det_num
   do l=1,det_num
 
+   phase = dble(det_exc(k,l,4))
    if (det_exc(k,l,3) == 0) then
      buffer = 0.
      do i=1,elec_alpha_num-mo_closed_num
       buffer += mo_value_p(det(i,k,1))*mo_value_p(det(i,l,1))
      enddo
-     density_alpha_value_p += det_coef(k)*det_coef(l)*buffer
+     density_alpha_value_p += phase*det_coef(k)*det_coef(l)*buffer
    else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,1) == 1) ) then
      call get_single_excitation(k,l,ik,il,1)
      buffer = mo_value_p(ik)*mo_value_p(il)
-     density_alpha_value_p += det_coef(k)*det_coef(l)*buffer
+     density_alpha_value_p += phase*det_coef(k)*det_coef(l)*buffer
    endif
 
   enddo
@@ -77,20 +79,22 @@ BEGIN_PROVIDER [ real, density_beta_value_p ]
 ! TODO vectorization
  integer :: k,j,l, ik, il
  real    :: buffer
+ real    :: phase
  PROVIDE det
  PROVIDE elec_beta_num
  do k=1,det_num
   do l=1,det_num
+   phase = dble(det_exc(k,l,4))
    if (det_exc(k,l,3) == 0) then
      buffer = 0.
      do i=1,elec_beta_num-mo_closed_num
       buffer += mo_value_p(det(i,k,2))*mo_value_p(det(i,l,2))
      enddo
-     density_beta_value_p += det_coef(k)*det_coef(l)*buffer
+     density_beta_value_p += phase*det_coef(k)*det_coef(l)*buffer
    else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,2) == 1) ) then
      call get_single_excitation(k,l,ik,il,2)
      buffer = mo_value_p(ik)*mo_value_p(il)
-     density_beta_value_p += det_coef(k)*det_coef(l)*buffer
+     density_beta_value_p += phase*det_coef(k)*det_coef(l)*buffer
    endif
   enddo
  enddo

src/det.irp.f

@@ -33,10 +33,10 @@ BEGIN_PROVIDER  [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ]
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 3) ]
+BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 4) ]
  implicit none
  BEGIN_DOC  
-! Degree of excitation between two determinants. The sign is the phase.
+! Degree of excitation between two determinants. Indices are alpha, beta, alpha+beta, phase
  END_DOC
 
  integer :: p
@@ -66,42 +66,57 @@ BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 3) ]
 
     enddo
 
-    det_exc(l,k,p) = det_exc(k,l,p)
    enddo
   enddo
 
  enddo
 
  do l=1,det_num
+  det_exc(l,l,3) = 0
   do k=l+1,det_num
    det_exc(k,l,3) = det_exc(k,l,1) + det_exc(k,l,2)
   enddo
  enddo
 
  ! Phase
- do p=1,2
+
-  do i=mo_closed_num,mo_num
+ do l=1,det_num
-   integer :: det_pos(det_num)
+  det_exc(l,l,4) = 1
-   do k=1,det_num
+  do k=l+1,det_num
-    det_pos(k) = 0
+   integer :: nperm
-    do j=1,elec_num_2(p)-mo_closed_num
+   nperm = 0
-     if (det(j,k,p) == i) then
+   do p=1,2
-      det_pos(k) = j
+    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
    enddo
-   do k=1,det_num
+   det_exc(k,l,4) = 1-2*mod( nperm, 2 )
-    do l=k+1,det_num
-     det_exc(k,l,3) *= -2*mod( (det_pos(k)+det_pos(l)), 2 )+1
-    enddo
-   enddo
   enddo
-
  enddo
 
- do l=1,det_num
+
-  do k=l+1,det_num
+ do p=1,4
-   det_exc(l,k,3) = det_exc(k,l,3)
+  do l=1,det_num
+   do k=1,l-1
+    det_exc(k,l,p) = det_exc(l,k,p)
+   enddo
   enddo
  enddo
 

src/eplf_function.irp.f

@@ -6,10 +6,8 @@ BEGIN_PROVIDER  [ real, eplf_gamma ]
  include 'constants.F'
  real            :: eps
  eps = -real(dlog(tiny(1.d0)))
-!real :: N
+ real :: N
-!N = 0.1
+ eplf_gamma = (4./3.*pi*density_value_p)**(2./3.) * eps
-!eplf_gamma = (4./(3.*N)*pi*density_value_p)**(2./3.) * eps
- eplf_gamma = density_value_p * eps
 !eplf_gamma = 1.e10
 !eplf_gamma = 1.e5
 END_PROVIDER
@@ -32,7 +30,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
  implicit none
  BEGIN_DOC  
-! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF
+! Array of all the <phi_i phi_j | exp(-gamma r^2)> for EPLF
  END_DOC
  integer :: i, j, k, l
  double precision :: t
@@ -47,12 +45,15 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
   do k=1,ao_num
    if (mo_coef(k,i) /= 0.) then
     do l=1,ao_num
-     t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k)
+     if (abs(ao_eplf_integral_matrix(l,k))>1.d-16) then
-     do j=i,mo_num
+      t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k)
-      mo_eplf_integral_matrix(j,i) += t*mo_coef_transp(j,l)
+      do j=i,mo_num
-     enddo
+       mo_eplf_integral_matrix(j,i) += t*mo_coef_transp(j,l)
+      enddo
+     endif
     enddo
    endif
+
   enddo
 
  enddo
@@ -78,20 +79,21 @@ END_PROVIDER
 
  PROVIDE mo_coef_transp
 
- do j=1,mo_closed_num
+ do i=1,mo_closed_num
-  do i=1,mo_closed_num
+  do j=1,mo_closed_num
-    eplf_up_up += 2.d0* 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) )
+      mo_value_p(j)*mo_eplf_integral_matrix(j,i) )
-    eplf_up_dn += 2.d0* mo_value_p(i)*mo_value_p(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
+ integer :: k,l,m,n,p,p2
  integer :: ik,il,jk,jl
- double precision :: ckl
+ double precision :: phase,dtemp(2)
- double precision :: phase
  integer :: exc
 
  PROVIDE det
@@ -100,101 +102,117 @@ END_PROVIDER
  do k=1,det_num
   do l=1,det_num
 
-   ckl = det_coef(k)*det_coef(l)
-
    exc = det_exc(k,l,3)
 
-   if ( exc < 0 ) then
+   dtemp(1) = 0.d0
-    phase = -1.0d0
+   dtemp(2) = 0.d0
-    exc = -exc
+   do p=1,2
-   else
+     p2 = 1+mod(p,2)
-    phase = 1.0d0
+     if ( exc == 0 ) then
-   endif
+      ! 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
 
-   if ( exc == 0 ) then
+      !- Open-closed shell interactions
-     ! Sum over all alpha-alpha and beta-beta interactions
+      do m=1,elec_num_2(p)-mo_closed_num
-     do p=1,2
+       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)
-         eplf_up_up += phase*ckl*mo_value_p(ik)* (  &
+         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
-     enddo
+ 
+     else if ( (exc == 1).and.(det_exc(k,l,p) == 1) ) then
 
-     ! Sum over all alpha-beta interactions
+     ! Sum over only the sigma-sigma interactions involving the excitation
-     do m=1,elec_beta_num-mo_closed_num
+      call get_single_excitation(k,l,ik,il,p)
-      jk = det(m,k,2)
+
-      jl = det(m,l,2)
+      !- Open-closed shell interactions
-      do n=1,elec_alpha_num-mo_closed_num
+      do j=1,mo_closed_num
-        ik = det(n,k,1)
+       dtemp(1) += mo_value_p(ik)* (  &
-        il = det(n,l,1)
+           mo_value_p(il)*mo_eplf_integral_matrix(j,j) - &
-        eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) &
+           mo_value_p(j)*mo_eplf_integral_matrix(j,il) )
-                            + mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) )
+       dtemp(2) += mo_value_p(ik)*mo_value_p(il)*mo_eplf_integral_matrix(j,j)
       enddo
-     enddo
 
-   else if ( exc == 1 ) then
+      !- 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 == 2).and.(det_exc(k,l,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 == 2).and.(det_exc(k,l,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)
     
-     do p=1,2
-      if ( det_exc(k,l,p) == 1 ) then
-      ! Sum over only the sigma-sigma interactions involving the excitation
-       call get_single_excitation(k,l,ik,il,p)
-       do m=1,elec_num_2(p)-mo_closed_num
-        jk = det(m,k,p)
-        jl = det(m,l,p)
-        eplf_up_up += phase*ckl*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
-
-       ! Sum over only the sigma-(sigma_bar) interactions involving the excitation
-       integer :: p2
-       p2 = 1+mod(p,2)
-       do m=1,elec_num_2(p2)-mo_closed_num
-        jk = det(m,k,p2)
-        jl = det(m,l,p2)
-        eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) &
-                            + mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) )
-       enddo
-      endif
-     enddo
-
-   else if (exc == 2) then
-
-     if ( ( det_exc(k,l,1) == 2 ).or.( det_exc(k,l,2) == 2 ) ) then
-     
-      ! Consider only the double excitations of same-spin electrons
-      if ( det_exc(k,l,1) == 2 ) then
-       call get_double_excitation(k,l,ik,jk,il,jl,1)
-      else if ( det_exc(k,l,2) == 2 ) then
-       call get_double_excitation(k,l,ik,jk,il,jl,2)
-      endif
-
-      eplf_up_up += phase*ckl*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 ( det_exc(k,l,1) == 1 ) then
-      ! Consider only the double excitations of opposite-spin electrons
-      call get_single_excitation(k,l,ik,jk,1)
-      call get_single_excitation(k,l,il,jl,2)
-
-      eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) &
-                  + mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) )
      endif
+   enddo
 
-   endif
+   phase = dble(det_exc(k,l,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
 
 
@@ -304,7 +322,7 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center)
   
 end function
  
- double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
+double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
   implicit none
   include 'constants.F'
  
@@ -348,24 +366,26 @@ end function
    di(ii) = 0.5d0*inv_p(2)*dble(ii)
   enddo
  
-  S(1,0) = (xp-xa) * S(0,0)
+  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) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0)
+    S(ii+1,0) = xab(1) * S(ii,0) + di(ii)*S(ii-1,0)
    enddo
   endif
  
-  S(0,1) = (xp-xb) * S(0,0)
+  S(0,1) = xab(2) * S(0,0)
   if (j>1) then
    do jj=1,j-1
-    S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-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) = (xp-xa) * S(0,jj) + di(jj) * S(0,jj-1)
+   S(1,jj) = xab(1) * S(0,jj) + di(jj) * S(0,jj-1)
    do ii=2,i
-    S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
+    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
  

src/overlap.irp.f

@@ -8,13 +8,9 @@ BEGIN_PROVIDER [ double precision, ao_overlap_matrix, (ao_num,ao_num) ]
  integer :: i, j
  double precision :: ao_overlap
  do j=1,ao_num
-  do i=j,ao_num
+  do i=1,j
    ao_overlap_matrix(i,j) = ao_overlap(i,j)
-  enddo
+   ao_overlap_matrix(j,i) = ao_overlap_matrix(i,j)
- enddo
- do j=1,ao_num
-  do i=1,j-1
-   ao_overlap_matrix(i,j) = ao_overlap(j,i)
   enddo
  enddo
 END_PROVIDER
@@ -237,14 +233,14 @@ double precision function ao_overlap(i,j)
 
  do q=1,ao_prim_num(j)
   do p=1,ao_prim_num(i)
-   integral(p,q) = integral(p,q)*ao_coef(p,i)*ao_coef(q,j)
+   integral(p,q) *= ao_coef(p,i)*ao_coef(q,j)
   enddo
  enddo
 
- ao_overlap = 0.
+ ao_overlap = 0.d0
  do q=1,ao_prim_num(j)
   do p=1,ao_prim_num(i)
-   ao_overlap = ao_overlap + integral(p,q)
+   ao_overlap += integral(p,q)
   enddo
  enddo
  

src/test_1d.irp.f

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