EPLF MCSCF tested => OK !

bin/to_ezfio.py

@@ -140,7 +140,7 @@ def write_ezfioFile(res,filename):
   ezfio.mo_basis_mo_coef = MoMatrix
 
 # Determinants
-  det_thr = 1.e-16
+  det_thr = 1.e-6
   closed_mos = res.closed_mos
   nactive = ezfio.get_mo_basis_mo_active_num()
   dets_a = []
@@ -166,7 +166,7 @@ def write_ezfioFile(res,filename):
   if len(dets_a[0]) > 0:
    for i in xrange(len(coef),0,-1):
      i -= 1
-     if coef[i] == 0.:
+     if abs(coef[i]) < det_thr :
        dets_a.pop(i)
        dets_b.pop(i)
        coef.pop(i)

src/density.irp.f

@@ -41,19 +41,28 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ]
  integer :: k,j,l, ik, il
  real    :: buffer
  real    :: phase
+ integer :: exc(4)
  PROVIDE det
  PROVIDE elec_alpha_num
  do k=1,det_num
   do l=1,det_num
 
-   phase = dble(det_exc(k,l,4))
+   exc(1) = abs(det_exc(k,l,1))
-   if (det_exc(k,l,3) == 0) then
+   exc(2) = abs(det_exc(k,l,2))
+   exc(3) = exc(1)+exc(2)
+   exc(4) = exc(1)*exc(2)
+   if (exc(4) /= 0) then
+    exc(4) = exc(4)/abs(exc(4))
+   endif
+
+   phase = dble(exc(4))
+   if (exc(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 += phase*det_coef(k)*det_coef(l)*buffer
-   else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,1) == 1) ) then
+   else if ( (exc(3) == 1).and.(exc(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 += phase*det_coef(k)*det_coef(l)*buffer
@@ -80,18 +89,27 @@ BEGIN_PROVIDER [ real, density_beta_value_p ]
  integer :: k,j,l, ik, il
  real    :: buffer
  real    :: phase
+ integer :: exc(4)
  PROVIDE det
  PROVIDE elec_beta_num
  do k=1,det_num
   do l=1,det_num
-   phase = dble(det_exc(k,l,4))
+   exc(1) = abs(det_exc(k,l,1))
-   if (det_exc(k,l,3) == 0) then
+   exc(2) = abs(det_exc(k,l,2))
+   exc(3) = exc(1)+exc(2)
+   exc(4) = exc(1)*exc(2)
+   if (exc(4) /= 0) then
+    exc(4) = exc(4)/abs(exc(4))
+   endif
+
+   phase = dble(exc(4))
+   if (exc(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 += phase*det_coef(k)*det_coef(l)*buffer
-   else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,2) == 1) ) then
+   else if ( (exc(3) == 1).and.(exc(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 += phase*det_coef(k)*det_coef(l)*buffer

src/det.irp.f

@@ -33,10 +33,11 @@ BEGIN_PROVIDER  [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ]
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 4) ]
+BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 2) ]
  implicit none
  BEGIN_DOC  
-! Degree of excitation between two determinants. Indices are alpha, beta, alpha+beta, phase
+! Degree of excitation between two determinants. Indices are alpha, beta
+! The sign is the phase factor
  END_DOC
 
  integer :: p
@@ -71,17 +72,9 @@ BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 4) ]
 
  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 l=1,det_num
-  det_exc(l,l,4) = 1
   do k=l+1,det_num
    integer :: nperm
    nperm = 0
@@ -106,13 +99,12 @@ BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 4) ]
       nperm += 1
      endif
     enddo
+    det_exc(k,l,p) *= (1-2*mod( nperm, 2 ))
    enddo
-   det_exc(k,l,4) = 1-2*mod( nperm, 2 )
   enddo
  enddo
 
-
+ do p=1,2
- do p=1,4
   do l=1,det_num
    do k=1,l-1
     det_exc(k,l,p) = det_exc(l,k,p)

src/eplf_function.irp.f

@@ -4,10 +4,10 @@ BEGIN_PROVIDER  [ real, eplf_gamma ]
 ! Value of the gaussian for the EPLF
  END_DOC
  include 'constants.F'
- real            :: eps
+ real            :: eps, N
+ N = 0.1
  eps = -real(dlog(tiny(1.d0)))
- real :: N
+ eplf_gamma = (4./3.*pi*density_value_p/N)**(2./3.) * eps
- eplf_gamma = (4./3.*pi*density_value_p)**(2./3.) * eps
 !eplf_gamma = 1.e10
 !eplf_gamma = 1.e5
 END_PROVIDER
@@ -94,7 +94,7 @@ END_PROVIDER
  integer :: k,l,m,n,p,p2
  integer :: ik,il,jk,jl
  double precision :: phase,dtemp(2)
- integer :: exc
+ integer :: exc(4)
 
  PROVIDE det
  PROVIDE elec_num_2
@@ -102,13 +102,21 @@ END_PROVIDER
  do k=1,det_num
   do l=1,det_num
 
-   exc = det_exc(k,l,3)
+   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 == 0 ) then
+     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
@@ -151,7 +159,7 @@ END_PROVIDER
        enddo
       enddo
  
-     else if ( (exc == 1).and.(det_exc(k,l,p) == 1) ) then
+     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)
@@ -186,7 +194,7 @@ END_PROVIDER
        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
+     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)
@@ -195,7 +203,7 @@ END_PROVIDER
              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
+     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)
@@ -206,7 +214,7 @@ END_PROVIDER
      endif
    enddo
 
-   phase = dble(det_exc(k,l,4))
+   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)