Merge branch 'master' of github.com:LCPQ/quantum_package

src/BiInts/map_integrals.irp.f

@@ -98,11 +98,12 @@ subroutine get_ao_bielec_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_
   PROVIDE ao_bielec_integrals_in_map
   thresh = ao_integrals_threshold
   
+  non_zero_int = 0
   if (ao_overlap_abs(j,l) < thresh) then
     out_val = 0.d0
     return
   endif
-  
+ 
   non_zero_int = 0
   do i=1,sze
     !DIR$ FORCEINLINE

src/CISD_SC2_selected/H_apply.irp.f

@@ -6,5 +6,9 @@ from perturbation import perturbations
 s = H_apply("PT2",SingleRef=True)
 s.set_perturbation("epstein_nesbet_sc2_projected")
 print s
+
+s = H_apply("PT2_en_sc2",SingleRef=True)
+s.set_perturbation("epstein_nesbet_sc2")
+print s
 END_SHELL
 

src/CISD_SC2_selected/NEEDED_MODULES

@@ -1 +1 @@
-AOs BiInts Bitmask CISD SC2 CISD_selected Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils
+AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils

src/CISD_SC2_selected/cisd_sc2.ezfio_config

@@ -0,0 +1,4 @@
+cisd_sc2_selected
+  n_det_max_cisd_sc2  integer 
+  pt2_max             double precision
+  do_pt2_end          logical 

src/CISD_SC2_selected/cisd_sc2_selection.irp.f

@@ -13,9 +13,25 @@ program cisd_sc2_selected
   pt2 = 1.d0
   perturbation = "epstein_nesbet_sc2_projected"
   E_old(1) = HF_energy
-  davidson_threshold = 1.d-6 
+  davidson_threshold = 1.d-8
+  if (N_det > n_det_max_cisd_sc2) then
+    call diagonalize_CI_SC2
+    call save_wavefunction
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    N_det = n_det_max_cisd_sc2
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_SC2_energy
+    print *,  'E+PT2    = ', CI_SC2_energy+pt2
+    print *,  '-----'
+  endif
 
-  do while (maxval(abs(pt2(1:N_st))) > 1.d-4)
+  do while (N_det < n_det_max_cisd_sc2.and.maxval(abs(pt2(1:N_st))) > pt2_max)
     print*,'----'
     print*,''
     call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag,  N_st)
@@ -37,37 +53,43 @@ program cisd_sc2_selected
       exit
     endif
   enddo
-  pt2 = 0.d0
+  N_det = min(n_det_max_cisd_sc2,N_det)
-  call H_apply_PT2(pt2, norm_pert, H_pert_diag,  N_st)
+  touch N_det psi_det psi_coef
   davidson_threshold = 1.d-10
   touch davidson_threshold davidson_criterion
-  do i = 1, N_st
+  call diagonalize_CI_SC2
-   max = 0.d0
+  pt2 = 0.d0
-
+  if(do_pt2_end)then
-    print*,''
+   threshold_selectors = 1.d0
-    print*,'-------------'
+   call H_apply_PT2(pt2, norm_pert, H_pert_diag,  N_st)
-    print*,'for state ',i
+   do i = 1, N_st
-    print*,''
+    max = 0.d0
-   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)       = ', (CI_SC2_energy(i)+pt2(i))
-   if(i==1)then
-    print *,  'E(SC2)+PT2(projctd)SC2       = ', (CI_SC2_energy(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
+     print*,''
+     print*,'-------------'
+     print*,'for state ',i
+     print*,''
+     print*,'N_det = ',N_det
+    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)       = ', CI_SC2_energy(i)+pt2(i)
+    print *,  'E_before(SC2)+PT2(SC2)_new   = ', CI_SC2_energy(i)+pt2(i)*(1.d0+norm_pert)
+   
+    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
   print*,'coucou'
+  endif
+  call save_wavefunction
   deallocate(pt2,norm_pert,H_pert_diag)
 end

src/CISD_SC2_selected/options.irp.f

@@ -0,0 +1,51 @@
+  BEGIN_PROVIDER [ integer, n_det_max_cisd_sc2 ]
+   implicit none
+   BEGIN_DOC
+  ! Get n_det_max_cisd_sc2 from EZFIO file
+   END_DOC
+   logical :: has_n_det_max_cisd_sc2
+   PROVIDE ezfio_filename
+   call ezfio_has_cisd_sc2_selected_n_det_max_cisd_sc2(has_n_det_max_cisd_sc2)
+   if (has_n_det_max_cisd_sc2) then
+     call ezfio_get_cisd_sc2_selected_n_det_max_cisd_sc2(n_det_max_cisd_sc2)
+   else
+     n_det_max_cisd_sc2 = 1000
+     call ezfio_set_cisd_sc2_selected_n_det_max_cisd_sc2(n_det_max_cisd_sc2)
+   endif
+   print*,'n_det_max_cisd_sc2 = ',n_det_max_cisd_sc2
+  END_PROVIDER
+
+  BEGIN_PROVIDER [ double precision , pt2_max ]
+   implicit none
+   BEGIN_DOC
+  ! Get pt2_max from EZFIO file
+   END_DOC
+   logical :: has_pt2_max
+   PROVIDE ezfio_filename
+   call ezfio_has_cisd_sc2_selected_pt2_max(has_pt2_max)
+   if (has_pt2_max) then
+     call ezfio_get_cisd_sc2_selected_pt2_max(pt2_max)
+   else
+     pt2_max = 1.d-3
+     call ezfio_set_cisd_sc2_selected_pt2_max(pt2_max)
+   endif
+   print*,'pt2_max = ',pt2_max
+  END_PROVIDER
+
+  BEGIN_PROVIDER [ logical, do_pt2_end ]
+   implicit none
+   BEGIN_DOC
+  ! Get do_pt2_end from EZFIO file
+   END_DOC
+   logical :: has_do_pt2_end
+   PROVIDE ezfio_filename
+   call ezfio_has_cisd_sc2_selected_do_pt2_end(has_do_pt2_end)
+   if (has_do_pt2_end) then
+     call ezfio_get_cisd_sc2_selected_do_pt2_end(do_pt2_end)
+   else
+     do_pt2_end = .True.
+     call ezfio_set_cisd_sc2_selected_do_pt2_end(do_pt2_end)
+   endif
+   print*,'do_pt2_end = ',do_pt2_end
+  END_PROVIDER
+

src/Dets/diagonalize_CI_SC2.irp.f

@@ -33,6 +33,7 @@ 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) = CI_eigenvectors(i,j)
     enddo
     CI_SC2_electronic_energy(j) = CI_electronic_energy(j) 

src/Dets/slater_rules.irp.f

@@ -571,6 +571,61 @@ subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx
 end
 
 
+subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat)
+  use bitmasks
+  BEGIN_DOC
+  ! <key|H|psi> for the various Nstate
+  !
+  ! returns in addition
+  !
+  ! the array of the index of the non connected determinants to key1
+  !
+  ! in order to know what double excitation can be repeated on key1
+  !
+  ! idx_repeat(0) is the number of determinants that can be used
+  !
+  ! to repeat the excitations
+  END_DOC
+  implicit none
+  integer, intent(in)            :: Nint, Ndet,Ndet_max,Nstate
+  integer(bit_kind), intent(in)  :: keys(Nint,2,Ndet)
+  integer(bit_kind), intent(in)  :: key(Nint,2)
+  double precision, intent(in)   :: coef(Ndet_max,Nstate)
+  double precision, intent(out)  :: i_H_psi_array(Nstate)
+  integer         , intent(out)  :: idx_repeat(0:Ndet)
+  
+  integer                        :: i, ii,j
+  double precision               :: phase
+  integer                        :: exc(0:2,2,2)
+  double precision               :: hij
+  integer                        :: idx(0:Ndet)
+  
+  ASSERT (Nint > 0)
+  ASSERT (N_int == Nint)
+  ASSERT (Nstate > 0)
+  ASSERT (Ndet > 0)
+  ASSERT (Ndet_max >= Ndet)
+  i_H_psi_array = 0.d0
+  call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat)
+  print*,'--------'
+  do ii=1,idx(0)
+    print*,'--'
+    i = idx(ii)
+    !DEC$ FORCEINLINE
+    call i_H_j(keys(1,1,i),key,Nint,hij)
+    if (i==1)then
+     print*,'i==1 !!'
+    endif
+    print*,coef(i,1) * hij,coef(i,1),hij
+    do j = 1, Nstate
+      i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
+    enddo
+    print*,i_H_psi_array(1)
+  enddo
+  print*,'------'
+end
+
+
 
 subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)
   use bitmasks

src/Perturbation/pert_sc2.irp.f

@@ -35,28 +35,25 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
   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
+  h =  diag_H_mat_elem(det_pert,Nint) + accu_e_corr
+  delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
 
-  delta_E = (CI_SC2_electronic_energy(1) - h)
-  delta_E = 1.d0/delta_E
 
-  c_pert(1) = i_H_psi_array(1) * delta_E
+  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(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
     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))
@@ -67,17 +64,16 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
   enddo
 
   degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) +                      &
-      popcnt(xor( ref_bitmask(1,2), det_pert(1,2)))
+           popcnt(xor( ref_bitmask(1,2), det_pert(1,2)))
   !DEC$ NOUNROLL
   do l=2,Nint
     degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) +            &
-        popcnt(xor( ref_bitmask(l,2), det_pert(l,2)))
+                     popcnt(xor( ref_bitmask(l,2), det_pert(l,2)))
   enddo
   if(degree==4)then
   ! <psi|delta_H|psi>
-   call i_H_j(ref_bitmask,det_pert,Nint,h0j)
+   H_pert_diag(1) = e_2_pert(1)
-   H_pert_diag(1) = c_pert(1) * h0j/coef_hf_selector
+   e_2_pert_fonda = H_pert_diag(1) 
-   e_2_pert_fonda = H_pert_diag(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)
@@ -115,3 +111,46 @@ double precision function repeat_all_e_corr(key_in)
  repeat_all_e_corr = accu
 
 end
+
+
+subroutine pt2_epstein_nesbet_sc2(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)
+  
+  BEGIN_DOC
+  ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various N_st states, but with the CISD_SC2 energies and coefficients
+  !
+  ! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
+  !
+  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem, h
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  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(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) = selection_criterion*selection_criterion_factor*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)
+        H_pert_diag(i) = h*c_pert(i)*c_pert(i)
+        e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
+    else
+      c_pert(i) = -1.d0
+      e_2_pert(i) = -dabs(i_H_psi_array(i))
+      H_pert_diag(i) = h
+    endif
+  enddo
+  
+end

src/SC2/NEEDED_MODULES

@@ -1 +0,0 @@
-AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output SingleRefMethod Utils Selectors_full

src/SC2/README.rst

@@ -1,37 +0,0 @@
-===============
-CISD_SC2 Module
-===============
-
-Documentation
-=============
-
-.. Do not edit this section. It was auto-generated from the
-.. NEEDED_MODULES file.
-
-`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/cisd_SC2.irp.f#L1>`_
-  Undocumented
-
-
-
-Needed Modules
-==============
-
-.. Do not edit this section. It was auto-generated from the
-.. NEEDED_MODULES file.
-
-* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
-* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
-* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
-* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
-* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
-* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
-* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
-* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
-* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
-* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
-* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
-* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
-* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
-* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
-* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
-

src/SC2/cisd_SC2.irp.f

@@ -1,14 +0,0 @@
-program cisd
-  implicit none
-  integer :: i
-
-  print *,  ' HF      = ', HF_energy
-  print *,  'N_states = ', N_states
-  call H_apply_cisd
-  print *,  'N_det = ', N_det
-  do i = 1,N_states
-   print *,  'energy  = ',CI_SC2_energy(i) 
-   print *,  'E_corr  = ',CI_SC2_electronic_energy(i) - ref_bitmask_energy
-  enddo
-
-end

src/Selectors_full/e_corr_selectors.irp.f

@@ -25,7 +25,13 @@ use bitmasks
  enddo
 END_PROVIDER
   BEGIN_PROVIDER[double precision, coef_hf_selector]
+ &BEGIN_PROVIDER[double precision, inv_selectors_coef_hf]
+ &BEGIN_PROVIDER[double precision, inv_selectors_coef_hf_squared]
  &BEGIN_PROVIDER[double precision, E_corr_per_selectors, (N_det_selectors)]
+ &BEGIN_PROVIDER[double precision, i_H_HF_per_selectors, (N_det_selectors)]
+ &BEGIN_PROVIDER[double precision, Delta_E_per_selector, (N_det_selectors)]
+ &BEGIN_PROVIDER[double precision, E_corr_double_only ]
+ &BEGIN_PROVIDER[double precision, E_corr_second_order ]
  implicit none
  BEGIN_DOC
  ! energy of correlation per determinant respect to the Hartree Fock determinant
@@ -39,25 +45,33 @@ END_PROVIDER
  ! coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants
  END_DOC
  integer :: i,degree
- double precision :: hij,inv_selectors_coef_hf
+ double precision :: hij,diag_H_mat_elem
+ E_corr_double_only = 0.d0
+ E_corr_second_order = 0.d0
  do i = 1, N_det_selectors
   if(exc_degree_per_selectors(i)==2)then
    call i_H_j(ref_bitmask,psi_selectors(1,1,i),N_int,hij)
+   i_H_HF_per_selectors(i) = hij
    E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij
+   E_corr_double_only += E_corr_per_selectors(i)
+   E_corr_second_order += hij * hij /(ref_bitmask_energy - diag_H_mat_elem(psi_selectors(1,1,i),N_int))
   elseif(exc_degree_per_selectors(i) == 0)then
    coef_hf_selector = psi_selectors_coef(i,1)
    E_corr_per_selectors(i) = -1000.d0
+   Delta_E_per_selector(i) = 0.d0
   else
    E_corr_per_selectors(i) = -1000.d0
   endif
  enddo
  if (dabs(coef_hf_selector) > 1.d-8) then
    inv_selectors_coef_hf = 1.d0/coef_hf_selector
+   inv_selectors_coef_hf_squared = inv_selectors_coef_hf * inv_selectors_coef_hf
  else
    inv_selectors_coef_hf = 0.d0
+   inv_selectors_coef_hf_squared = 0.d0
  endif
  do i = 1,n_double_selectors
   E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf
  enddo
-
+ E_corr_double_only = E_corr_double_only * inv_selectors_coef_hf
  END_PROVIDER