CIS_DT cleaned, add Full_ci/parameters.irp.f

src/CIS_dressed/CIS_DT_lapack.irp.f

@@ -0,0 +1,68 @@
+program CIS_DT
+ implicit none
+ integer :: i
+ print*,'MP2_dresssing=',mp2_dressing
+ print*,'standard_doubles=',standard_doubles
+ print*,'n_state_CIS=',n_state_CIS
+ print*,'n_core_cis=',n_core_cis
+ print*,'n_act_cis=',n_act_cis
+ print*,''
+ print*,'nuc repulsion E=',nuclear_repulsion
+if (mp2_dressing==.true.) then
+ print*,'correlation E=',MP2_corr_energy
+else
+ print*,'correlation E=',EN2_corr_energy
+endif
+ do i = 1,n_state_CIS
+   print*,''
+   print*,'i = ',i
+   print*,'CIS     = ',eigenvalues_CIS(i)
+   print*,'CIS(DdT)= ',eigenvalues_CIS_dress_D_dt(i)
+   print*,'s2(DdT)    = ',s_2_CIS_dress_D_dt(i)
+   print*,'<x>   = ',CIS_states_properties(1,i)
+   print*,'<y>   = ',CIS_states_properties(2,i)
+   print*,'<z>   = ',CIS_states_properties(3,i)
+   print*,'<xx>  = ',CIS_states_properties(4,i)
+   print*,'<yy>  = ',CIS_states_properties(5,i)
+   print*,'<zz>  = ',CIS_states_properties(6,i)
+   print*,''
+ enddo
+ double precision :: delta_E_CIS,delta_E_CIS_DT,convert
+
+ convert = 1.d0
+ print*,'Excitation energies :      CIS           CIS_DT         (Hartree)' 
+ do i = 2, n_state_CIS
+  delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1)
+  delta_E_CIS_DT = eigenvalues_CIS_dress_D_dt(i) - eigenvalues_CIS_dress_D_dt(1)
+  write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.5,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert
+ enddo
+
+ convert = 27.2114d0
+ print*,'Excitation energies :      CIS           CIS_DT         (eV)' 
+ do i = 2, n_state_CIS
+  delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1)
+  delta_E_CIS_DT = eigenvalues_CIS_dress_D_dt(i) - eigenvalues_CIS_dress_D_dt(1)
+  write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.6,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert
+ enddo
+
+
+ convert = 219475d0
+ print*,'Excitation energies :      CIS           CIS_DT         (cm-1)' 
+ do i = 2, n_state_CIS
+  delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1)
+  delta_E_CIS_DT = eigenvalues_CIS_dress_D_dt(i) - eigenvalues_CIS_dress_D_dt(1)
+  write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.1,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert
+ enddo
+
+ convert = 627.51d0
+ print*,'Excitation energies :      CIS           CIS_DT         (Kcal/mol)' 
+ do i = 2, n_state_CIS
+  delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1)
+  delta_E_CIS_DT = eigenvalues_CIS_dress_D_dt(i) - eigenvalues_CIS_dress_D_dt(1)
+  write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.5,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert
+ enddo
+
+!if(save_all_dm_cis)then
+! call save_all_density_matrix
+!endif
+end

src/CIS_dressed/CIS_providers.irp.f

@@ -207,13 +207,13 @@
  implicit none
  double precision,allocatable  :: delta_H_matrix_doub(:,:)
  double precision,allocatable  :: eigvalues(:),eigvectors(:,:)
- double precision :: overlap,max_overlap,s2
+ double precision :: overlap,max_overlap,s2,e_corr
  integer :: i_overlap,i,j,k
  allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS))
  allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS))
   do i = 1,n_state_CIS
-!  call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles
+   call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles
-   delta_H_matrix_doub = 0.d0
+!  delta_H_matrix_doub = 0.d0
    
    do j = 1,size_psi_CIS
     do k = 1,size_psi_CIS
@@ -246,7 +246,6 @@
    enddo
    print*,i,i_overlap
    print*,'overlap = ',max_overlap
-   i_overlap = i
    overlap_Ddt=max_overlap
    do k = 1,size_psi_CIS
     eigenvectors_CIS_dress_D_dt(k,i) = eigvectors(k,i_overlap)

src/CIS_dressed/EN2.irp.f

@@ -1,8 +1,8 @@
 
   BEGIN_PROVIDER [double precision, EN2_corr_energy]
-&BEGIN_PROVIDER [double precision, p_imp_EN,(elec_alpha_num+1:n_act_cis)]
+  &BEGIN_PROVIDER [double precision, p_imp_EN,(elec_alpha_num+1:n_act_cis)]
-&BEGIN_PROVIDER [double precision, h_imp_EN,(n_core_cis+1:elec_alpha_num)]
+  &BEGIN_PROVIDER [double precision, h_imp_EN,(n_core_cis+1:elec_alpha_num)]
-&BEGIN_PROVIDER [double precision, hp_imp_EN,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)]
+  &BEGIN_PROVIDER [double precision, hp_imp_EN,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)]
  
   BEGIN_DOC
   !Calculation of the EN2 correlation energy (EN2_corr_energy)
@@ -18,7 +18,11 @@
   double precision :: delta_e_ik,delta_e_ikj
   double precision :: delta_e !delta epsilons
   double precision :: delta_e_tmp,H_jj_total
+  integer, allocatable :: key_in(:,:)
+  integer, allocatable :: key_out(:,:)
   integer :: ispin1,ispin2,i_ok
+  allocate(key_in(N_int,2))
+  allocate(key_out(N_int,2))
  
   print*,'EN2_corr_energy'
  
@@ -31,29 +35,25 @@
    hp_imp_EN(i,k)=0.d0
   enddo
  enddo
- print*,'HF_energy = ',HF_energy
- print*,'1'
  
   if(EN_2_2)then
    do i=n_core_cis+1,elec_alpha_num
-  h_imp_EN(i)=0.d0
    
-  e_i=diagonal_Fock_matrix_mo(i)
+    e_i=Fock_matrix_diag_mo(i)
    
     do k=elec_alpha_num+1,n_act_cis
-   hp_imp_EN(i,k)=0.d0
     
-   e_k=diagonal_Fock_matrix_mo(k)
+     e_k=Fock_matrix_diag_mo(k)
      delta_e_ik=e_i-e_k
   
      !same spin contribution for EN2_corr_energy 
      do j=i+1,elec_alpha_num
-    e_j=diagonal_Fock_matrix_mo(j)
+      e_j=Fock_matrix_diag_mo(j)
       delta_e_ikj=delta_e_ik+e_j
   
       !same spin contribution for EN2_corr_energy and a part of the contribution to p_imp_EN,h_imp_EN
       do l=k+1,n_act_cis
-     e_l=diagonal_Fock_matrix_mo(l)
+       e_l=Fock_matrix_diag_mo(l)
        delta_e=delta_e_ikj-e_l
        delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - &
                        mo_bielec_integral_jj_anti(k,l)
@@ -82,12 +82,12 @@
   
      do j=n_core_cis+1,elec_alpha_num
       if(j==i)cycle
-    e_j=diagonal_Fock_matrix_mo(j)
+      e_j=Fock_matrix_diag_mo(j)
       delta_e_ikj=delta_e_ik+e_j
    
       do l=elec_alpha_num+1,n_act_cis
        if(l==k)cycle
-     e_l=diagonal_Fock_matrix_mo(l)
+       e_l=Fock_matrix_diag_mo(l)
        delta_e=delta_e_ikj-e_l
        delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - &
                        mo_bielec_integral_jj_anti(k,l)
@@ -109,11 +109,11 @@
   
      !different spin contribution
      do j=n_core_cis+1,elec_beta_num
-    e_j=diagonal_Fock_matrix_mo(j)
+      e_j=Fock_matrix_diag_mo(j)
       delta_e_ikj=delta_e_ik+e_j
   
       do l=elec_beta_num+1,n_act_cis
-     e_l=diagonal_Fock_matrix_mo(l)
+       e_l=Fock_matrix_diag_mo(l)
        delta_e=delta_e_ikj-e_l 
        delta_e=delta_e-mo_bielec_integral_jj(i,j) - &
                        mo_bielec_integral_jj(k,l)
@@ -121,7 +121,6 @@
                        mo_bielec_integral_jj_anti(j,l) + &
                        mo_bielec_integral_jj(i,l) + &
                        mo_bielec_integral_jj(j,k)
-
        direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
   
        hij=direct*direct
@@ -136,28 +135,25 @@
      enddo
     enddo
    enddo
- print*,'2'
 
   else
    do i=n_core_cis+1,elec_alpha_num
-  h_imp_EN(i)=0.d0
    
-  e_i=diagonal_Fock_matrix_mo(i)
+    e_i=Fock_matrix_diag_mo(i)
    
     do k=elec_alpha_num+1,n_act_cis
-   hp_imp_EN(i,k)=0.d0
     
-   e_k=diagonal_Fock_matrix_mo(k)
+     e_k=Fock_matrix_diag_mo(k)
      delta_e_ik=e_i-e_k
   
      !same spin contribution for EN2_corr_energy 
      do j=i+1,elec_alpha_num
-    e_j=diagonal_Fock_matrix_mo(j)
+      e_j=Fock_matrix_diag_mo(j)
       delta_e_ikj=delta_e_ik+e_j
   
       !same spin contribution for EN2_corr_energy and a part of the contribution to p_imp_EN,h_imp_EN
       do l=k+1,n_act_cis
-     e_l=diagonal_Fock_matrix_mo(l)
+       e_l=Fock_matrix_diag_mo(l)
        delta_e=delta_e_ikj-e_l
        delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - &
                        mo_bielec_integral_jj_anti(k,l)
@@ -170,7 +166,7 @@
        exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map)
   
        hij=(direct-exchg)*(direct-exchg)
-     hij = hij*hij/delta_e
+       hij=hij/delta_e
   
        EN2_corr_energy=EN2_corr_energy+2*hij
        p_imp_EN(k)=p_imp_EN(k)+hij
@@ -186,12 +182,12 @@
   
      do j=n_core_cis+1,elec_alpha_num
       if(j==i)cycle
-    e_j=diagonal_Fock_matrix_mo(j)
+      e_j=Fock_matrix_diag_mo(j)
       delta_e_ikj=delta_e_ik+e_j
    
       do l=elec_alpha_num+1,n_act_cis
        if(l==k)cycle
-     e_l=diagonal_Fock_matrix_mo(l)
+       e_l=Fock_matrix_diag_mo(l)
        delta_e=delta_e_ikj-e_l
        delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - &
                        mo_bielec_integral_jj_anti(k,l)
@@ -204,7 +200,7 @@
        exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map)
   
        hij=(direct-exchg)*(direct-exchg)
-     hij = hij*hij/delta_e
+       hij=hij/delta_e
   
        hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij
   
@@ -213,11 +209,11 @@
   
      !different spin contribution
      do j=n_core_cis+1,elec_beta_num
-    e_j=diagonal_Fock_matrix_mo(j)
+      e_j=Fock_matrix_diag_mo(j)
       delta_e_ikj=delta_e_ik+e_j
   
       do l=elec_beta_num+1,n_act_cis
-     e_l=diagonal_Fock_matrix_mo(l)
+       e_l=Fock_matrix_diag_mo(l)
        delta_e=delta_e_ikj-e_l 
        delta_e=delta_e-mo_bielec_integral_jj(i,j) - &
                        mo_bielec_integral_jj(k,l)
@@ -225,11 +221,10 @@
                        mo_bielec_integral_jj_anti(j,l) + &
                        mo_bielec_integral_jj(i,l) + &
                        mo_bielec_integral_jj(j,k)
-
        direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
   
        hij=direct*direct
-     hij = hij*hij/delta_e
+       hij=hij/delta_e
   
        EN2_corr_energy=EN2_corr_energy+hij
        p_imp_EN(k)=p_imp_EN(k)+hij
@@ -240,13 +235,9 @@
      enddo
     enddo
    enddo
-
   endif
  
-
   print*,'EN correlation energy=',EN2_corr_energy
   print*,'EN correlation energy=',EN2_corr_energy
  
   END_PROVIDER
-
-

src/CIS_dressed/MP2.irp.f

@@ -29,22 +29,22 @@
  do i=n_core_cis+1,elec_alpha_num
   h_imp(i)=0.d0
  
-  e_i=diagonal_Fock_matrix_mo(i)
+  e_i=Fock_matrix_diag_mo(i)
  
   do k=elec_alpha_num+1,n_act_cis
    hp_imp(i,k)=0.d0
   
-   e_k=diagonal_Fock_matrix_mo(k)
+   e_k=Fock_matrix_diag_mo(k)
    delta_e_ik=e_i-e_k
 
    !same spin contribution for MP2_corr_energy 
    do j=i+1,elec_alpha_num
-    e_j=diagonal_Fock_matrix_mo(j)
+    e_j=Fock_matrix_diag_mo(j)
     delta_e_ikj=delta_e_ik+e_j
 
     !same spin contribution for MP2_corr_energy and a part of the contribution to p_imp and h_imp
     do l=k+1,n_act_cis
-     e_l=diagonal_Fock_matrix_mo(l)
+     e_l=Fock_matrix_diag_mo(l)
      delta_e=delta_e_ikj-e_l
 
      direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
@@ -60,7 +60,7 @@
 
     !remaining same spin contribution for p_imp    
     do l=elec_alpha_num+1,k-1
-     e_l=diagonal_Fock_matrix_mo(l)
+     e_l=Fock_matrix_diag_mo(l)
      delta_e=delta_e_ikj-e_l 
 
      direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
@@ -75,11 +75,11 @@
 
    !remaining same spin contribution for h_imp
    do j=n_core_cis+1,i-1
-    e_j=diagonal_Fock_matrix_mo(j)
+    e_j=Fock_matrix_diag_mo(j)
     delta_e_ikj=delta_e_ik+e_j
 
     do l=k+1,n_act_cis
-     e_l=diagonal_Fock_matrix_mo(l)
+     e_l=Fock_matrix_diag_mo(l)
      delta_e=delta_e_ikj-e_l
 
      direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
@@ -94,11 +94,11 @@
 
    !same spin contribution for hp_imp
    do j=n_core_cis+1,elec_alpha_num
-    e_j=diagonal_Fock_matrix_mo(j)
+    e_j=Fock_matrix_diag_mo(j)
     delta_e_ikj=delta_e_ik+e_j
  
     do l=elec_alpha_num+1,n_act_cis
-     e_l=diagonal_Fock_matrix_mo(l)
+     e_l=Fock_matrix_diag_mo(l)
      delta_e=delta_e_ikj-e_l
 
      direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
@@ -113,11 +113,11 @@
 
    !different spin contribution
    do j=n_core_cis+1,elec_beta_num
-    e_j=diagonal_Fock_matrix_mo(j)
+    e_j=Fock_matrix_diag_mo(j)
     delta_e_ikj=delta_e_ik+e_j
 
     do l=elec_beta_num+1,n_act_cis
-     e_l=diagonal_Fock_matrix_mo(l)
+     e_l=Fock_matrix_diag_mo(l)
      delta_e=delta_e_ikj-e_l 
 
      direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
@@ -385,25 +385,20 @@
   enddo
 
  else !EN Dressing
-  print*,'coucou !'
   dress_T_discon_array_CIS(1)=EN2_corr_energy
-  print*,'coucou !'
   
-! do i=n_core_cis+1,elec_alpha_num
+  do i=n_core_cis+1,elec_alpha_num
-!  print*,'i',i,n_core_cis
+   do k=elec_alpha_num+1,n_act_cis
-!  do k=elec_alpha_num+1,n_act_cis
+    key=psi_CIS_adress(i,k)
-!   print*,'k',k,n_act_cis
-!   key=psi_CIS_adress(i,k)
-!   
-!   dress_T_discon(i,k)=EN2_corr_energy-p_imp_EN(k)-h_imp_EN(i)+hp_imp_EN(i,k)
     
-!   dress_T_discon_array_CIS(key) = dress_T_discon(i,k)
+    dress_T_discon(i,k)=EN2_corr_energy-p_imp_EN(k)-h_imp_EN(i)+hp_imp_EN(i,k)
-!   dress_T_discon_array_CIS(key+1) = dress_T_discon(i,k)
 
-!  enddo
+    dress_T_discon_array_CIS(key) = dress_T_discon(i,k)
-! enddo
+    dress_T_discon_array_CIS(key+1) = dress_T_discon(i,k)
+
+   enddo
+  enddo
  end if
- print*,'end'
 
  END_PROVIDER
 
@@ -809,24 +804,24 @@
  delta_H_trip=0.d0
 !do i=5,6
  do i=n_core_cis+1,elec_alpha_num
-  e_i=diagonal_Fock_matrix_mo(i)
+  e_i=Fock_matrix_diag_mo(i)
  !r=7
   do r=elec_alpha_num+1,n_act_cis
-   e_r=diagonal_Fock_matrix_mo(r)
+   e_r=Fock_matrix_diag_mo(r)
    delta_e_ir=e_i-e_r
    key_ir=psi_CIS_adress(i,r)
    do j=i+1,elec_alpha_num
-    e_j=diagonal_Fock_matrix_mo(j)
+    e_j=Fock_matrix_diag_mo(j)
     delta_e_irj=delta_e_ir+e_j
     do s=r+1,n_act_cis
-     e_s=diagonal_Fock_matrix_mo(s)
+     e_s=Fock_matrix_diag_mo(s)
      delta_e_irjs=delta_e_irj-e_s
      !alpha-alpha-alpha
      do k=j+1,elec_alpha_num
-      e_k=diagonal_Fock_matrix_mo(k)
+      e_k=Fock_matrix_diag_mo(k)
       delta_e_irjsk=delta_e_irjs+e_k
       do t=s+1,n_act_cis
-       e_t=diagonal_Fock_matrix_mo(t)
+       e_t=Fock_matrix_diag_mo(t)
        delta_e=delta_e_irjsk-e_t
        energy=1.d0/(ref_energy-delta_e)
        occ=i
@@ -898,10 +893,10 @@
      enddo
      !alpha-alpha-beta
      do k=n_core_cis+1,elec_alpha_num
-      e_k=diagonal_Fock_matrix_mo(k)
+      e_k=Fock_matrix_diag_mo(k)
       delta_e_irjsk=delta_e_irjs+e_k
       do t=elec_alpha_num+1,n_act_cis
-       e_t=diagonal_Fock_matrix_mo(t)
+       e_t=Fock_matrix_diag_mo(t)
        delta_e=delta_e_irjsk-e_t
        energy=1.d0/(ref_energy-delta_e)
        occ=i
@@ -992,23 +987,23 @@
 !!generating the Singles included in the Triples
 !! do m=1,n_state_CIS
 !do i=n_core_cis+1,elec_alpha_num
-! e_i=diagonal_Fock_matrix_mo(i)
+! e_i=Fock_matrix_diag_mo(i)
 
 
 
 !!do r=elec_alpha_num+1,n_state_CIS
 ! do r=elec_alpha_num+1,n_act_cis
-!  e_r=diagonal_Fock_matrix_mo(r)
+!  e_r=Fock_matrix_diag_mo(r)
 !  delta_e_ir=e_i-e_r
 !key_ir=psi_CIS_adress(i,r)
 
 !  !alpha-alpha-x (=beta-beta-x)
 !  do j=i+1,elec_alpha_num
-!   e_j=diagonal_Fock_matrix_mo(j)
+!   e_j=Fock_matrix_diag_mo(j)
 !   delta_e_irj=delta_e_ir+e_j
 !key_jr=psi_CIS_adress(j,r)
 !   do s=r+1,n_act_cis
-!    e_s=diagonal_Fock_matrix_mo(s)
+!    e_s=Fock_matrix_diag_mo(s)
 !    delta_e_irjs=delta_e_irj-e_s
 
 ! key_is=psi_CIS_adress(i,s)
@@ -1017,13 +1012,13 @@
 
 !    !alpha-alpha-alpha (=beta-beta-beta)
 !    do k=j+1,elec_alpha_num
-!     e_k=diagonal_Fock_matrix_mo(k)
+!     e_k=Fock_matrix_diag_mo(k)
 !     delta_e_irjsk=delta_e_irjs+e_k
 !   key_kr=psi_CIS_adress(k,r)
 !   key_ks=psi_CIS_adress(k,s)
 
 !     do t=s+1,n_act_cis
-!      e_t=diagonal_Fock_matrix_mo(t)
+!      e_t=Fock_matrix_diag_mo(t)
 !      delta_e=delta_e_irjsk-e_t
 
 !     key_it=psi_CIS_adress(i,t)
@@ -1157,7 +1152,7 @@
 
 !    !alpha-alpha-beta (=beta-beta-alpha)
 !    do k=n_core_cis+1,elec_beta_num
-!     e_k=diagonal_Fock_matrix_mo(k)
+!     e_k=Fock_matrix_diag_mo(k)
 !     delta_e_irjsk=delta_e_irjs+e_k
 
 !  key_kr=psi_CIS_adress(k,r)
@@ -1166,7 +1161,7 @@
 
 
 !     do t=elec_beta_num+1,n_act_cis
-!      e_t=diagonal_Fock_matrix_mo(t)
+!      e_t=Fock_matrix_diag_mo(t)
 !      delta_e=delta_e_irjsk-e_t
 
 !   key_it=psi_CIS_adress(i,t)
@@ -1298,53 +1293,3 @@
 !enddo
 
  END
-
-subroutine diexcitation(i,j,k,l,ispin1,ispin2,key_in,key_out,i_ok,Nint)
-  implicit none
-  use bitmasks
-  ! realize the double excitation i-->k (ispin1) +  j-->l (ispin2)  on key_in
-  ! returns key_out and i_ok (i_ok = 0 means not possible, i_ok = 1 means the excitation was possible)
-  integer, intent(in) :: ispin1,ispin2,i,j,k,l,Nint
-  integer(bit_kind), intent(in) :: key_in(Nint,2)
-  integer, intent(out):: i_ok
-  integer(bit_kind), intent(out):: key_out(Nint,2)
-  integer :: k_hole,j_hole,k_particl,j_particl,i_nint,Nelec_alpha,Nelec_beta
-  integer :: i_test_hole,i_test_particl
-  key_out = key_in
-
-  k_hole = ishft(i-1,-bit_kind_shift)+1
-  j_hole = i-ishft(k_hole-1,bit_kind_shift)-1
-  i_test_hole = ibset(0,j_hole)
-  if(iand(key_in(k_hole,ispin1),i_test_hole).ne.i_test_hole)then
-   i_ok = 0
-   return
-  endif
-  key_out(k_hole,ispin1) = ibclr(key_out(k_hole,ispin1),j_hole)
-  k_particl = ishft(k-1,-bit_kind_shift)+1
-  j_particl = k-ishft(k_particl-1,bit_kind_shift)-1
-  i_test_particl= ibset(0,j_particl)
-  if(iand(key_in(k_particl,ispin1),i_test_particl).ne.0)then
-   i_ok = 0
-   return
-  endif
-  key_out(k_particl,ispin1) = ibset(key_out(k_particl,ispin1),j_particl)
-
-  k_hole = ishft(j-1,-bit_kind_shift)+1
-  j_hole = j-ishft(k_hole-1,bit_kind_shift)-1
-  i_test_hole = ibset(0,j_hole)
-  if(iand(key_in(k_hole,ispin2),i_test_hole).ne.i_test_hole)then
-   i_ok = 0
-   return
-  endif
-  key_out(k_hole,ispin2) = ibclr(key_out(k_hole,ispin2),j_hole)
-  k_particl = ishft(l-1,-bit_kind_shift)+1
-  j_particl = l-ishft(k_particl-1,bit_kind_shift)-1
-  i_test_particl = ibset(0,j_particl)
-  if(iand(key_in(k_particl,ispin2),i_test_particl).ne.0)then
-   i_ok = 0
-   return
-  endif
-  key_out(k_particl,ispin2) = ibset(key_out(k_particl,ispin2),j_particl)
-  i_ok = 1
-  end
-

src/CIS_dressed/natural_particl_orbitals.irp.f

@@ -56,7 +56,7 @@
   enddo
   do k = 1, mo_tot_num
    do l = 1, mo_tot_num
-    c_k = eigvectors(k,j) * eigvectors(l,j)
+    c_k = eigvectors(k,i) * eigvectors(l,i)
     particle_natural_orb_CIS_properties(1,i) += c_k * mo_dipole_x(k,l)
     particle_natural_orb_CIS_properties(2,i) += c_k * mo_dipole_y(k,l)
     particle_natural_orb_CIS_properties(3,i) += c_k * mo_dipole_z(k,l)

src/CIS_dressed/options.irp.f

@@ -10,7 +10,8 @@ BEGIN_PROVIDER [ integer , n_state_cis ]
   if (has) then
     call ezfio_get_cis_dressed_n_state_cis(n_state_cis)
   else
-    n_state_cis = 5
+    n_state_cis = 10
+    call ezfio_set_cis_dressed_n_state_cis(n_state_cis)
   endif
 
 END_PROVIDER

src/CIS_dressed/repeat_all_doubles.irp.f

@@ -0,0 +1,178 @@
+subroutine repeat_all_doubles(key_in, e_corr)
+ implicit none
+ integer(bit_kind), intent(in) :: key_in(N_int,2)
+ double precision, intent(out) :: e_corr
+ integer :: i,j,k,l
+ integer :: s1,s2
+ integer :: i_ok
+ double precision :: hij,get_mo_bielec_integral,diag_H_mat_elem,delta_e
+ integer(bit_kind) :: key_out(N_int,2)
+ ! same spin (alpha)
+ if(mp2_dressing)then
+ s1 = 1
+ e_corr = 0.d0
+ do i = n_core_cis + 1, elec_alpha_num
+  do j = i + 1, elec_alpha_num
+   do k = elec_alpha_num + 1, n_act_cis
+    do l = k+1, n_act_cis
+     ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in>
+      call diexcitation(i,j,k,l,s1,s1,key_in,key_out,i_ok,N_int)
+      if(i_ok.ne.0)then
+       hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) & 
+            -get_mo_bielec_integral(i,j,l,k,mo_integrals_map)
+       e_corr += hij*hij/(Fock_matrix_diag_mo(i)+Fock_matrix_diag_mo(j) & 
+                         -Fock_matrix_diag_mo(l)-Fock_matrix_diag_mo(k))
+      endif
+    enddo
+   enddo
+  enddo
+ enddo
+
+ s1 = 2
+ do i = n_core_cis + 1, elec_alpha_num
+  do j = i + 1, elec_alpha_num
+   do k = elec_alpha_num + 1, n_act_cis
+    do l = k+1, n_act_cis
+     ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in>
+      call diexcitation(i,j,k,l,s1,s1,key_in,key_out,i_ok,N_int)
+      if(i_ok.ne.0)then
+       hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) & 
+            -get_mo_bielec_integral(i,j,l,k,mo_integrals_map)
+       e_corr += hij*hij/(Fock_matrix_diag_mo(i)+Fock_matrix_diag_mo(j) & 
+                         -Fock_matrix_diag_mo(l)-Fock_matrix_diag_mo(k))
+      endif
+    enddo
+   enddo
+  enddo
+ enddo
+
+ s1 = 1
+ s2 = 2
+ do i = n_core_cis + 1, elec_alpha_num
+  do j = n_core_cis + 1, elec_alpha_num
+   do k = elec_alpha_num + 1, n_act_cis
+    do l = elec_alpha_num + 1, n_act_cis
+     ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in>
+      call diexcitation(i,j,k,l,s1,s2,key_in,key_out,i_ok,N_int)
+      if(i_ok.ne.0)then
+       hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) 
+       e_corr += hij*hij/(Fock_matrix_diag_mo(i)+Fock_matrix_diag_mo(j) & 
+                         -Fock_matrix_diag_mo(l)-Fock_matrix_diag_mo(k))
+      endif
+    enddo
+   enddo
+  enddo
+ enddo
+
+ else
+  ! same spin (alpha)
+  s1 = 1
+  e_corr = 0.d0
+  do i = n_core_cis + 1, elec_alpha_num
+   do j = i + 1, elec_alpha_num
+    do k = elec_alpha_num + 1, n_act_cis
+     do l = k+1, n_act_cis
+      ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in>
+       call diexcitation(i,j,k,l,s1,s1,key_in,key_out,i_ok,N_int)
+       if(i_ok.ne.0)then
+        hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) & 
+             -get_mo_bielec_integral(i,j,l,k,mo_integrals_map)
+        call diexcitation(i,j,k,l,s1,s1,ref_bitmask,key_out,i_ok,N_int)
+        delta_e = HF_energy -  diag_H_mat_elem(key_out,N_int)
+        e_corr += hij*hij/delta_e
+       endif
+     enddo
+    enddo
+   enddo
+  enddo
+
+  s1 = 2
+  do i = n_core_cis + 1, elec_alpha_num
+   do j = i + 1, elec_alpha_num
+    do k = elec_alpha_num + 1, n_act_cis
+     do l = k+1, n_act_cis
+      ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in>
+       call diexcitation(i,j,k,l,s1,s1,key_in,key_out,i_ok,N_int)
+       if(i_ok.ne.0)then
+        hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) & 
+             -get_mo_bielec_integral(i,j,l,k,mo_integrals_map)
+        call diexcitation(i,j,k,l,s1,s1,ref_bitmask,key_out,i_ok,N_int)
+        delta_e = HF_energy -  diag_H_mat_elem(key_out,N_int)
+        e_corr += hij*hij/delta_e
+       endif
+     enddo
+    enddo
+   enddo
+  enddo
+
+  s1 = 1
+  s2 = 2
+  do i = n_core_cis + 1, elec_alpha_num
+   do j = n_core_cis + 1, elec_alpha_num
+    do k = elec_alpha_num + 1, n_act_cis
+     do l = elec_alpha_num + 1, n_act_cis
+      ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in>
+       call diexcitation(i,j,k,l,s1,s2,key_in,key_out,i_ok,N_int)
+       if(i_ok.ne.0)then
+        hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) 
+        call diexcitation(i,j,k,l,s1,s2,ref_bitmask,key_out,i_ok,N_int)
+        delta_e = HF_energy -  diag_H_mat_elem(key_out,N_int)
+        e_corr += hij*hij/delta_e
+       endif
+     enddo
+    enddo
+   enddo
+  enddo
+ endif
+end
+
+subroutine diexcitation(i,j,k,l,ispin1,ispin2,key_in,key_out,i_ok,Nint)
+  implicit none
+  use bitmasks
+  ! realize the double excitation i-->k (ispin1) +  j-->l (ispin2)  on key_in
+  ! returns key_out and i_ok (i_ok = 0 means not possible, i_ok = 1 means the excitation was possible)
+  integer, intent(in) :: ispin1,ispin2,i,j,k,l,Nint
+  integer(bit_kind), intent(in) :: key_in(Nint,2)
+  integer, intent(out):: i_ok
+  integer(bit_kind), intent(out):: key_out(Nint,2)
+  integer :: k_hole,j_hole,k_particl,j_particl,i_nint,Nelec_alpha,Nelec_beta
+  integer(bit_kind) :: i_test_hole,i_test_particl
+  character*(512) :: output(2)
+  key_out = key_in
+
+  k_hole = ishft(i-1,-bit_kind_shift)+1
+  j_hole = i-ishft(k_hole-1,bit_kind_shift)-1
+  i_test_hole = ibset(0_bit_kind,j_hole)
+  if(iand(key_in(k_hole,ispin1),i_test_hole).ne.i_test_hole)then
+   i_ok = 0
+   return
+  endif
+  key_out(k_hole,ispin1) = ibclr(key_out(k_hole,ispin1),j_hole)
+  k_particl = ishft(k-1,-bit_kind_shift)+1
+  j_particl = k-ishft(k_particl-1,bit_kind_shift)-1
+  i_test_particl= ibset(0_bit_kind,j_particl)
+  if(iand(key_in(k_particl,ispin1),i_test_particl).ne.0_bit_kind)then
+   i_ok = 0
+   return
+  endif
+  key_out(k_particl,ispin1) = ibset(key_out(k_particl,ispin1),j_particl)
+
+  k_hole = ishft(j-1,-bit_kind_shift)+1
+  j_hole = j-ishft(k_hole-1,bit_kind_shift)-1
+  i_test_hole = ibset(0_bit_kind,j_hole)
+  if(iand(key_in(k_hole,ispin2),i_test_hole).ne.i_test_hole)then
+   i_ok = 0
+   return
+  endif
+  key_out(k_hole,ispin2) = ibclr(key_out(k_hole,ispin2),j_hole)
+  k_particl = ishft(l-1,-bit_kind_shift)+1
+  j_particl = l-ishft(k_particl-1,bit_kind_shift)-1
+  i_test_particl = ibset(0_bit_kind,j_particl)
+  if(iand(key_in(k_particl,ispin2),i_test_particl).ne.0_bit_kind)then
+   i_ok = 0
+   return
+  endif
+  key_out(k_particl,ispin2) = ibset(key_out(k_particl,ispin2),j_particl)
+  i_ok = 1
+  end
+

src/Full_CI/parameters.irp.f

@@ -0,0 +1,51 @@
+BEGIN_PROVIDER [ integer, N_det_max_fci ]
+ implicit none
+ BEGIN_DOC
+! Max number od determinants in the wave function
+ END_DOC
+ logical                        :: exists
+ PROVIDE ezfio_filename
+ call ezfio_has_full_ci_n_det_max_fci(exists)
+ if (exists) then
+   call ezfio_get_full_ci_n_det_max_fci(n_det_max_fci)
+ else
+   n_det_max_fci = 10000
+ endif
+ call write_int(output_full_ci,n_det_max_fci,'Max number of determinants ')
+ ASSERT (n_det_max_fci > 0)
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ logical , do_pt2_end ]
+ implicit none
+ BEGIN_DOC
+! if True then compute the PT2 when the selection process is finished
+ END_DOC
+ logical                        :: exists
+ PROVIDE ezfio_filename
+ call ezfio_has_full_ci_do_pt2_end(exists)
+ if (exists) then
+   call ezfio_get_full_ci_do_pt2_end(do_pt2_end)
+ else
+   do_pt2_end = .True.
+ endif
+!call write_i(output_full_ci,do_pt2_end,' computes the PT2 at the end of the selection ')
+ ASSERT (do_pt2_end > 0)
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision , pt2_max ]
+ implicit none
+ BEGIN_DOC
+! The selection process stops when the largest PT2 (for all the states) is lower than pt2_max
+! in absolute value
+ END_DOC
+ logical                        :: exists
+ PROVIDE ezfio_filename
+ call ezfio_has_full_ci_pt2_max(exists)
+ if (exists) then
+   call ezfio_get_full_ci_pt2_max(pt2_max)
+ else
+   pt2_max = 0.1d0
+ endif
+END_PROVIDER

src/Hartree_Fock/Fock_matrix.irp.f

@@ -115,6 +115,7 @@ END_PROVIDER
  if (do_direct_integrals) then
    PROVIDE all_utils ao_overlap_abs ao_integrals_threshold
    PROVIDE HF_density_matrix_ao_alpha  HF_density_matrix_ao_beta
+   PROVIDE ao_bi_elec_integral_alpha
    !$OMP PARALLEL DEFAULT(NONE) &
    !$OMP PRIVATE(i,j,l,k1,k,integral) &
    !$OMP SHARED(ao_num,ao_bi_elec_integral_alpha,ao_mono_elec_integral,&
@@ -147,6 +148,8 @@ END_PROVIDER
    !$OMP END PARALLEL
 
  else
+   PROVIDE ao_bielec_integrals_in_map
+           
    !$OMP PARALLEL DEFAULT(NONE) &
    !$OMP PRIVATE(i,j,l,k1,k,integral,ao_ints_val,ao_ints_idx,kmax) &
    !$OMP SHARED(ao_num,ao_bi_elec_integral_alpha,ao_mono_elec_integral,&