minor changes

src/CISD_SC2_selected/H_apply.irp.f

@@ -4,7 +4,7 @@ from generate_h_apply import *
 from perturbation import perturbations
 
 s = H_apply("PT2",SingleRef=True)
-s.set_perturbation("epstein_nesbet_sc2_projected")
+s.set_perturbation("epstein_nesbet_sc2_no_projected")
 print s
 
 s = H_apply("PT2_en_sc2",SingleRef=True)

src/CISD_SC2_selected/cisd_sc2_selection.irp.f

@@ -12,8 +12,9 @@ program cisd_sc2_selected
   
   pt2 = 1.d0
   perturbation = "epstein_nesbet_sc2_projected"
+                 
   E_old(1) = HF_energy
-  davidson_threshold = 1.d-8
+  davidson_threshold = 1.d-10
   if (N_det > n_det_max_cisd_sc2) then
     call diagonalize_CI_SC2
     call save_wavefunction
@@ -31,6 +32,8 @@ program cisd_sc2_selected
     print *,  '-----'
   endif
 
+  integer :: i_count
+  i_count = 0
   do while (N_det < n_det_max_cisd_sc2.and.maxval(abs(pt2(1:N_st))) > pt2_max)
     print*,'----'
     print*,''
@@ -49,6 +52,13 @@ program cisd_sc2_selected
      E_old(i) = CI_SC2_energy(i)
     enddo
 !   print *,  'E corr           = ', (E_old(1)) - HF_energy
+     if(dabs(E_old(i) - CI_SC2_energy(i) ).le.1.d-12)then
+      i_count += 1
+      selection_criterion_factor = selection_criterion_factor * 0.5d0
+      if(i_count > 5)then
+       exit
+      endif
+     endif
     if (abort_all) then
       exit
     endif
@@ -81,7 +91,7 @@ program cisd_sc2_selected
     print *,  'PT2(SC2)                     = ', pt2(i)
     print *,  'E(SC2)                       = ', CI_SC2_energy(i)
     print *,  'E_before(SC2)+PT2(SC2)       = ', CI_SC2_energy(i)+pt2(i)
-    print *,  'E_before(SC2)+PT2(SC2)_new   = ', CI_SC2_energy(i)+pt2(i)*(1.d0+norm_pert)
+    print *,  'E_before(SC2)+PT2(SC2)_new   = ', CI_SC2_energy(i)+pt2(i)* (1.d0 + norm_pert) - H_pert_diag(i)
    
     print*,'greater coeficient of the state : ',dabs(psi_coef(imax,i))
     call get_excitation_degree(ref_bitmask,psi_det(1,1,imax),degree,N_int)

src/Dets/determinants.ezfio_config

@@ -13,4 +13,5 @@ determinants
   det_occ                integer           (electrons_elec_alpha_num,determinants_det_num,2)
   det_coef               double precision  (determinants_det_num)
   read_wf                logical 
+  expected_s2            double precision 
 

src/Dets/diagonalize_CI.irp.f

@@ -59,12 +59,21 @@ END_PROVIDER
     call lapack_diag(eigenvalues,eigenvectors,                       &
         H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
     CI_electronic_energy(:) = 0.d0
-    do j=1,min(N_states,N_det)
+    integer :: i_state
+    double precision :: s2
+    j=0
+    i_state = 0
+!   do while(i_state.lt.N_states)
+      j+=1
+!     call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
+!     if(dabs(s2-expected_s2).le.0.1d0)then
+       i_state += 1
        do i=1,N_det
-        CI_eigenvectors(i,j) = eigenvectors(i,j)
+         CI_eigenvectors(i,i_state) = eigenvectors(i,j)
-      enddo
-      CI_electronic_energy(j) = eigenvalues(j)
        enddo
+       CI_electronic_energy(i_state) = eigenvalues(j)
+!     endif
+!   enddo
     deallocate(eigenvectors,eigenvalues)
   endif
   

src/Dets/diagonalize_CI_SC2.irp.f

@@ -20,7 +20,7 @@ END_PROVIDER
  BEGIN_DOC
  ! convergence of the correlation energy of SC2 iterations
  END_DOC
- threshold_convergence_SC2 = 1.d-8
+ threshold_convergence_SC2 = 1.d-10
 
  END_PROVIDER
  BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states) ]
@@ -33,8 +33,8 @@ END_PROVIDER
   
   do j=1,N_states
     do i=1,N_det
-!     CI_SC2_eigenvectors(i,j) = psi_coef(i,j)
+      CI_SC2_eigenvectors(i,j) = psi_coef(i,j)
-      CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j)
+!     CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j)
     enddo
     CI_SC2_electronic_energy(j) = CI_electronic_energy(j) 
   enddo

src/Dets/options.irp.f

@@ -20,6 +20,7 @@ T.set_doc       ( "If true, read the wave function from the EZFIO file" )
 T.set_ezfio_name( "read_wf" )
 T.set_output    ( "output_dets" )
 print T
+
 END_SHELL
 
 

src/Dets/s2.irp.f

@@ -42,6 +42,21 @@ BEGIN_PROVIDER [ double precision, S_z ]
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, expected_s2]
+ implicit none
+   PROVIDE ezfio_filename
+   logical :: has_expected_s2
+
+   call ezfio_has_determinants_expected_s2(has_expected_s2)
+   if (has_expected_s2) then
+     call ezfio_get_determinants_expected_s2(expected_s2)
+   else
+     expected_s2 = elec_alpha_num - elec_beta_num + 0.5d0 * ((elec_alpha_num - elec_beta_num)**2*0.5d0 - (elec_alpha_num-elec_beta_num))
+     call ezfio_set_determinants_expected_s2(expected_s2)
+   endif
+
+END_PROVIDER 
+
 
 subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
  implicit none

src/Perturbation/pert_sc2.irp.f

@@ -72,8 +72,8 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
   enddo
   if(degree==4)then
   ! <psi|delta_H|psi>
-   H_pert_diag(1) = e_2_pert(1)
+   e_2_pert_fonda = e_2_pert(1) 
-   e_2_pert_fonda = H_pert_diag(1) 
+   H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1)
    do i = 1, N_st
     do j = 1, idx_repeat(0)
      e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i)
@@ -83,6 +83,76 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
   
 end
 
+
+subroutine pt2_epstein_nesbet_SC2_no_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: Nint,ndet,N_st
+  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
+  double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
+  double precision               :: i_H_psi_array(N_st)
+  integer                        :: idx_repeat(0:ndet)
+  
+  BEGIN_DOC
+  ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various N_st states, 
+  ! 
+  ! but  with the correction in the denominator 
+  ! 
+  ! comming from the interaction of that determinant with all the others determinants 
+  ! 
+  ! that can be repeated by repeating all the double excitations
+  !
+  ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
+  ! 
+  ! that could be repeated to this determinant.
+  !
+  ! In addition, for the perturbative energetic contribution you have the standard second order
+  !
+  ! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
+  !
+  ! and also the purely projected contribution 
+  !
+  ! H_pert_diag = <HF|H|det_pert> c_pert
+  END_DOC
+  
+  integer                        :: i,j,degree,l
+  double precision               :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E
+  double precision               :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
+
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+
+  call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
+  accu_e_corr = 0.d0
+  !$IVDEP
+  do i = 1, idx_repeat(0)
+   accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
+  enddo
+  h =  diag_H_mat_elem(det_pert,Nint) + accu_e_corr
+  delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
+
+
+  c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
+  e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
+
+  do i =2,N_st
+    H_pert_diag(i) = h
+    if  (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
+      c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
+      e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
+    else
+      c_pert(i) = i_H_psi_array(i)
+      e_2_pert(i) = -dabs(i_H_psi_array(i))
+    endif
+  enddo
+end
+
+
+
+
+
 double precision function repeat_all_e_corr(key_in)
  implicit none
  integer(bit_kind), intent(in)  :: key_in(N_int,2)

src/Utils/LinearAlgebra.irp.f

@@ -307,6 +307,72 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
   deallocate(A,eigenvalues)
 end
 
+subroutine lapack_diag_s2(eigvalues,eigvectors,H,nmax,n)
+  implicit none
+  BEGIN_DOC
+  ! Diagonalize matrix H
+  !
+  ! H is untouched between input and ouptut
+  !
+  ! eigevalues(i) = ith lowest eigenvalue of the H matrix
+  !
+  ! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
+  !
+  END_DOC
+  integer, intent(in)            :: n,nmax
+  double precision, intent(out)  :: eigvectors(nmax,n)
+  double precision, intent(out)  :: eigvalues(n)
+  double precision, intent(in)   :: H(nmax,n)
+  double precision,allocatable   :: eigenvalues(:)
+  double precision,allocatable   :: work(:)
+  double precision,allocatable   :: A(:,:)
+  integer                        :: lwork, info, i,j,l,k, liwork
+
+  allocate(A(nmax,n),eigenvalues(n))
+! print*,'Diagonalization by jacobi'
+! print*,'n = ',n
+
+  A=H
+  lwork = 2*n*n + 6*n+ 1
+  allocate (work(lwork))
+
+  lwork = -1
+  call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    info )
+  if (info < 0) then
+    print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value'
+    stop 2
+  endif
+  lwork  = int( work( 1 ) )
+  deallocate (work)
+
+  allocate (work(lwork))
+  call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    info )
+  deallocate(work)
+
+  if (info < 0) then
+    print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value'
+    stop 2
+  else if( info > 0  ) then
+     write(*,*)'DSYEV Failed'
+     stop 1
+  end if
+
+  eigvectors = 0.d0
+  eigvalues = 0.d0
+  do j = 1, n
+    eigvalues(j) = eigenvalues(j)
+    do i = 1, n
+      eigvectors(i,j) = A(i,j)
+    enddo
+  enddo
+  deallocate(A,eigenvalues)
+end
+
+
+
+
 subroutine lapack_partial_diag(eigvalues,eigvectors,H,nmax,n,n_st)
   implicit none
   BEGIN_DOC

src/Utils/one_e_integration.irp.f

@@ -15,10 +15,10 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
   call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,&
       beta,power_A,power_B,A_center,B_center,dim)
   
-  if(fact_p.lt.0.000001d0)then
+! if(fact_p.lt.0.000001d0)then
-    overlap_gaussian_x = 0.d0
+!   overlap_gaussian_x = 0.d0
-    return
+!   return
-  endif
+! endif
   
   overlap_gaussian_x = 0.d0
   integer                        :: i

src/Utils/util.irp.f

@@ -49,17 +49,17 @@ double precision function binom_func(i,j)
 end
 
 
- BEGIN_PROVIDER [ double precision, binom, (0:20,0:20) ]
+ BEGIN_PROVIDER [ double precision, binom, (0:40,0:40) ]
-&BEGIN_PROVIDER [ double precision, binom_transp, (0:20,0:20) ]
+&BEGIN_PROVIDER [ double precision, binom_transp, (0:40,0:40) ]
   implicit none
   BEGIN_DOC
   ! Binomial coefficients
   END_DOC
   integer                        :: k,l
   double precision               :: fact, f
-  do k=0,20
+  do k=0,40
     f = fact(k)
-    do l=0,20
+    do l=0,40
       binom(k,l) = f/(fact(l)*fact(k-l))
       binom_transp(l,k) = binom(k,l)
     enddo