CISD_SC2_selected works better, new pertrubation sc2, better selection

src/CISD_SC2/README.rst

@@ -8,6 +8,41 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
+`cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/SC2.irp.f#L1>`_
+  CISD+SC2 method              :: take off all the disconnected terms of a CISD (selected or not)
+  .br
+  dets_in : bitmasks corresponding to determinants
+  .br
+  u_in : guess coefficients on the various states. Overwritten
+  on exit
+  .br
+  dim_in : leftmost dimension of u_in
+  .br
+  sze : Number of determinants
+  .br
+  N_st : Number of eigenstates
+  .br
+  Initial guess vectors are not necessarily orthonormal
+
+`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/SC2.irp.f#L169>`_
+  Undocumented
+
+`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/cisd_SC2.irp.f#L1>`_
+  Undocumented
+
+`ci_sc2_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/diagonalize_CI_SC2.irp.f#L19>`_
+  Eigenvectors/values of the CI matrix
+
+`ci_sc2_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/diagonalize_CI_SC2.irp.f#L18>`_
+  Eigenvectors/values of the CI matrix
+
+`ci_sc2_energy <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/diagonalize_CI_SC2.irp.f#L1>`_
+  N_states lowest eigenvalues of the CI matrix
+
+`diagonalize_ci_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/diagonalize_CI_SC2.irp.f#L38>`_
+  Replace the coefficients of the CI states by the coefficients of the
+  eigenstates of the CI matrix
+
 
 
 Needed Modules

src/CISD_SC2/SC2.irp.f

@@ -39,10 +39,10 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit)
   integer                        :: degree_exc(sze)
   integer                        :: i_ok
   double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:)
-  if(sze<500)then
+  if(sze.le.1000)then
-   allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
+   allocate (eigenvectors(size(H_matrix_all_dets,1),sze))
-   allocate (H_matrix_tmp(size(H_matrix_all_dets,1),N_det))
+   allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze))
-   allocate (eigenvalues(N_det))
+   allocate (eigenvalues(sze))
    do i = 1, sze
     do j = 1, sze
      H_matrix_tmp(i,j) = H_matrix_all_dets(i,j)
@@ -119,14 +119,14 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit)
       H_jj_dressed(i) += accu
     enddo
     
-    if(sze>500)then
+    if(sze>1000)then
      call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint,output_CISD_SC2)
     else
      do i = 1,sze
       H_matrix_tmp(i,i) = H_jj_dressed(i)
      enddo
      call lapack_diag(eigenvalues,eigenvectors,                       &
-         H_matrix_tmp,size(H_matrix_all_dets,1),N_det)
+         H_matrix_tmp,size(H_matrix_all_dets,1),sze)
      do j=1,min(N_states,sze)
        do i=1,sze
          u_in(i,j) = eigenvectors(i,j)

src/CISD_SC2/diagonalize_CI_SC2.irp.f

@@ -1,20 +1,3 @@
-BEGIN_PROVIDER [ character*(64), diag_algorithm ]
-  implicit none
-  BEGIN_DOC
-  ! Diagonalization algorithm (Davidson or Lapack)
-  END_DOC
-  if (N_det > 500) then
-    diag_algorithm = "Davidson"
-  else
-    diag_algorithm = "Lapack"
-  endif
-
-  if (N_det < N_states) then
-    diag_algorithm = "Lapack"
-  endif
-  
-END_PROVIDER
-
 BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states) ]
   implicit none
   BEGIN_DOC
@@ -46,7 +29,6 @@ END_PROVIDER
     enddo
     CI_SC2_electronic_energy(j) = CI_electronic_energy(j) 
   enddo
-  print*,'E(CISD) = ',CI_electronic_energy(1)
   
     
   call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &

src/CISD_SC2_selected/cisd_sc2_selection.irp.f

@@ -1,10 +1,11 @@
 program cisd_sc2_selected
   implicit none
   integer                        :: i,k
+  use bitmasks
 
   
   double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:),E_old(:)
-  integer                        :: N_st, iter
+  integer                        :: N_st, iter,degree
   character*(64)                 :: perturbation
   N_st = N_states
   allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st),E_old(N_st))
@@ -12,25 +13,53 @@ program cisd_sc2_selected
   pt2 = 1.d0
   perturbation = "epstein_nesbet_sc2_projected"
   E_old(1) = HF_energy
-  do while (maxval(abs(pt2(1:N_st))) > 1.d-9)
+  do while (maxval(abs(pt2(1:N_st))) > 1.d-10)
     print*,'----'
     print*,''
     call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag,  N_st)
     call diagonalize_CI_SC2
     print *,  'N_det                        = ', N_det
     do i = 1, N_st
+     print*,'state ',i
      print *,  'PT2(SC2)                     = ', pt2(i)
      print *,  'E(SC2)                       = ', CI_SC2_energy(i) 
-     print *,  'E_before(SC2)+PT2(SC2)       = ', (E_old(i)+pt2(i))
+     print *,  'E_before(SC2)+PT2(SC2)       = ', (E_old(i)+pt2(i)) 
      if(i==1)then
       print *,  'E(SC2)+PT2(projctd)SC2       = ', (E_old(i)+H_pert_diag(i)) 
      endif
+     E_old(i) = CI_SC2_energy(i)
     enddo
 !   print *,  'E corr           = ', (E_old(1)) - HF_energy
-    E_old = CI_SC2_energy
     if (abort_all) then
       exit
     endif
   enddo
+  do i = 1, N_st
+   max = 0.d0
+
+    print*,''
+    print*,'-------------'
+    print*,'for state ',i
+    print*,''
+   do k = 1, N_det
+    if(dabs(psi_coef(k,i)).gt.max)then
+     max = dabs(psi_coef(k,i))
+     imax = k
+    endif
+   enddo
+   double precision :: max
+   integer :: imax
+   print *,  'PT2(SC2)                     = ', pt2(i)
+   print *,  'E(SC2)                       = ', CI_SC2_energy(i)
+   print *,  'E_before(SC2)+PT2(SC2)       = ', (E_old(i)+pt2(i))
+   if(i==1)then
+    print *,  'E(SC2)+PT2(projctd)SC2       = ', (E_old(i)+H_pert_diag(i))
+   endif
+
+   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)
+   print*,'degree of excitation of such determinant : ',degree
+   
+  enddo
   deallocate(pt2,norm_pert,H_pert_diag)
 end

src/Dets/README.rst

@@ -133,33 +133,36 @@ Documentation
 `n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L20>`_
   Number of determinants in the wave function
 
-`n_det_reference <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L75>`_
+`n_det_max_jacobi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L38>`_
+  Maximum number of determinants diagonalized my jacobi
+
+`n_det_reference <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L93>`_
   Number of determinants in the reference wave function
 
 `n_states <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L3>`_
   Number of states to consider
 
-`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L84>`_
+`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L102>`_
   Contribution of determinants to the state-averaged density
 
-`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L105>`_
+`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L123>`_
   Wave function sorted by determinants (state-averaged)
 
-`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L47>`_
+`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L65>`_
   The wave function. Initialized with Hartree-Fock if the EZFIO file
   is empty
 
-`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L104>`_
+`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L122>`_
   Wave function sorted by determinants (state-averaged)
 
-`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L46>`_
+`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L64>`_
   The wave function. Initialized with Hartree-Fock if the EZFIO file
   is empty
 
-`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L38>`_
+`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L56>`_
   Size of the psi_det/psi_coef arrays
 
-`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L103>`_
+`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L121>`_
   Wave function sorted by determinants (state-averaged)
 
 `double_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L40>`_
@@ -298,7 +301,7 @@ Documentation
   Returns <i|H|j> where i and j are determinants
 
 `i_h_psi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L491>`_
-  <key|H|psi> for the various Nstate
+  <key|H|psi> for the various Nstates
 
 `i_h_psi_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L527>`_
   <key|H|psi> for the various Nstate

src/Dets/determinants.ezfio_config

@@ -6,4 +6,5 @@ determinants
   psi_coef      double precision (determinants_n_det,determinants_n_states)
   psi_det       integer*8 (determinants_N_int*determinants_bit_kind/8,2,determinants_n_det)
   H_apply_threshold  double precision
+  n_det_max_jacobi integer
 

src/Dets/determinants.irp.f

@@ -35,6 +35,24 @@ BEGIN_PROVIDER [ integer, N_det ]
 END_PROVIDER
 
 
+BEGIN_PROVIDER [ integer, N_det_max_jacobi ]
+ implicit none
+ BEGIN_DOC
+ ! Maximum number of determinants diagonalized my jacobi
+ END_DOC
+ logical                        :: exists
+ PROVIDE ezfio_filename
+ call ezfio_has_determinants_n_det_max_jacobi(exists)
+ if (exists) then
+   call ezfio_get_determinants_n_det_max_jacobi(N_det_max_jacobi)
+ else
+   N_det_max_jacobi = 1500
+   call ezfio_set_determinants_n_det_max_jacobi(N_det_max_jacobi)
+ endif
+ ASSERT (N_det_max_jacobi > 0)
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ integer, psi_det_size ]
  implicit none
  BEGIN_DOC

src/Dets/diagonalize_CI.irp.f

@@ -3,7 +3,7 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ]
   BEGIN_DOC
   ! Diagonalization algorithm (Davidson or Lapack)
   END_DOC
-  if (N_det > 500) then
+  if (N_det > N_det_max_jacobi) then
     diag_algorithm = "Davidson"
   else
     diag_algorithm = "Lapack"

src/Perturbation/README.rst

@@ -95,110 +95,29 @@ Documentation
 `pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L1>`_
   compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
   .br
-  for the various n_st states.
+  for the various N_st states.
   .br
   c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
   .br
   e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
   .br
 
-`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L38>`_
+`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L40>`_
   compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
   .br
-  for the various n_st states.
+  for the various N_st states.
   .br
   e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> -  E(i) )  - sqrt( ( <det_pert|H|det_pert> -  E(i)) ^2 + 4 <psi(i)|H|det_pert>^2  )
   .br
   c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
   .br
 
-`pt2_epstein_nesbet_2x2_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L129>`_
+`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L1>`_
-  compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
-  .br
-  for the various n_st states.
-  .br
-  but  with the correction in the denominator
-  .br
-  comming from the interaction of that determinant with all the others determinants
-  .br
-  that can be repeated by repeating all the double excitations
-  .br
-  : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
-  .br
-  that could be repeated to this determinant.
-  .br
-  <det_pert|H|det_pert> --->  <det_pert|H|det_pert> + delta_e_corr
-  .br
-  e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> -  E(i) )  - sqrt( ( <det_pert|H|det_pert> -  E(i)) ^2 + 4 <psi(i)|H|det_pert>^2  )
-  .br
-  c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
-  .br
-
-`pt2_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L75>`_
-  compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
-  .br
-  for the various n_st states,
-  .br
-  but  with the correction in the denominator
-  .br
-  comming from the interaction of that determinant with all the others determinants
-  .br
-  that can be repeated by repeating all the double excitations
-  .br
-  : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
-  .br
-  that could be repeated to this determinant.
-  .br
-  <det_pert|H|det_pert> --->  <det_pert|H|det_pert> + delta_e_corr
-  .br
-  c_pert(i) = <psi(i)|H|det_pert>/( E(i) - (<det_pert|H|det_pert> ) )
-  .br
-  e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - (<det_pert|H|det_pert> ) )
-  .br
-
-`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L185>`_
-  compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
-  .br
-  for the various n_st states,
-  .br
-  but  with the correction in the denominator
-  .br
-  comming from the interaction of that determinant with all the others determinants
-  .br
-  that can be repeated by repeating all the double excitations
-  .br
-  : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
-  .br
-  that could be repeated to this determinant.
-  .br
-  BUT on the contrary with ""pt2_epstein_nesbet_SC2"", you compute the energy by projection
-  .br
-  <det_pert|H|det_pert> --->  <det_pert|H|det_pert> + delta_e_corr
-  .br
-  c_pert(1) = 1/c_HF <psi(i)|H|det_pert>/( E(i) - (<det_pert|H|det_pert> ) )
-  .br
-  e_2_pert(1) = <HF|H|det_pert> c_pert(1)
-  .br
-  NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !!
-  .br
-
-`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L2>`_
   Fill the H_apply buffer with determiants for the selection
 
-`n_det_selectors <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L120>`_
-  Undocumented
-
-`psi_selectors <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L125>`_
-  On what we apply <i|H|psi> for perturbation. If selection, it may be 0.9 of the norm.
-
-`psi_selectors_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L126>`_
-  On what we apply <i|H|psi> for perturbation. If selection, it may be 0.9 of the norm.
-
-`psi_selectors_size <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L116>`_
-  Undocumented
-
 `remove_small_contributions <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L81>`_
-  Remove determinants with small contributions
+  Remove determinants with small contributions. N_states is assumed to be
+  provided.
 
 `selection_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L68>`_
   Threshold to select determinants. Set by selection routines.

src/Perturbation/epstein_nesbet.irp.f

@@ -24,7 +24,10 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_s
   call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
   h = diag_H_mat_elem(det_pert,Nint)
   do i =1,N_st
-    if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
+    if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
+      c_pert(i) = -1.d0
+      e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
+    else if  (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
         c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
         H_pert_diag(i) = h*c_pert(i)*c_pert(i)
         e_2_pert(i) = c_pert(i) * i_H_psi_array(i)

src/Perturbation/pert_sc2.irp.f

@@ -34,7 +34,7 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
   
   integer                        :: i,j,degree
   double precision               :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E
-  double precision               :: repeat_all_e_corr,accu_e_corr_tmp
+  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)
@@ -51,9 +51,14 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
   c_pert(1) = i_H_psi_array(1) * delta_E
   e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
   H_pert_diag(1) = c_pert(1) * h0j/coef_hf_selector
+  e_2_pert_fonda = H_pert_diag(1)
+
   do i =2,N_st
     H_pert_diag(i) = h
-    if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
+    if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then
+      c_pert(i) = -1.d0
+      e_2_pert(i) = -2.d0
+    else if  (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
       c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h)
       e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
     else
@@ -61,6 +66,16 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
       e_2_pert(i) = -dabs(i_H_psi_array(i))
     endif
   enddo
+
+  call get_excitation_degree(ref_bitmask,det_pert,degree,Nint)
+  if(degree==2)then
+  ! <psi|delta_H|psi>
+   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)
+    enddo
+   enddo
+  endif
   
 end
 

src/Perturbation/selection.irp.f

@@ -72,7 +72,7 @@ end
  BEGIN_DOC
  ! Threshold to select determinants. Set by selection routines.
  END_DOC
- selection_criterion =  1.d0 
+ selection_criterion =  10.d0 
  selection_criterion_factor = 0.01d0
  selection_criterion_min = selection_criterion
 

src/Utils/LinearAlgebra.irp.f

@@ -201,6 +201,8 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
   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