noL tested for Ne and O

src/bi_ort_ints/no_dressing_energy.irp.f

@@ -0,0 +1,66 @@
+
+! ---
+
+BEGIN_PROVIDER [double precision, energy_1e_noL_HF]
+
+  implicit none
+  integer :: i
+
+  PROVIDE mo_bi_ortho_tc_one_e
+
+  energy_1e_noL_HF = 0.d0
+  do i = 1, elec_beta_num
+    energy_1e_noL_HF += mo_bi_ortho_tc_one_e(i,i)
+  enddo
+  do i = 1, elec_alpha_num
+    energy_1e_noL_HF += mo_bi_ortho_tc_one_e(i,i)
+  enddo
+
+  print*, "energy_1e_noL_HF = ", energy_1e_noL_HF
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [double precision, energy_2e_noL_HF]
+
+  implicit none
+  integer :: i, j
+
+  PROVIDE mo_bi_ortho_tc_two_e
+
+  energy_2e_noL_HF = 0.d0
+  ! down-down & down-down
+  do i = 1, elec_beta_num
+    do j = 1, elec_beta_num
+      energy_2e_noL_HF += (mo_bi_ortho_tc_two_e(i,j,i,j) - mo_bi_ortho_tc_two_e(j,i,i,j))
+    enddo
+  enddo
+  ! down-down & up-up
+  do i = 1, elec_beta_num
+    do j = 1, elec_alpha_num
+      energy_2e_noL_HF += mo_bi_ortho_tc_two_e(i,j,i,j)
+    enddo
+  enddo
+  ! up-up & down-down
+  do i = 1, elec_alpha_num
+    do j = 1, elec_beta_num
+      energy_2e_noL_HF += mo_bi_ortho_tc_two_e(i,j,i,j)
+    enddo
+  enddo
+  ! up-up & up-up
+  do i = 1, elec_alpha_num
+    do j = 1, elec_alpha_num
+      energy_2e_noL_HF += (mo_bi_ortho_tc_two_e(i,j,i,j) - mo_bi_ortho_tc_two_e(j,i,i,j))
+    enddo
+  enddo
+
+  ! 0.5 x is in the Slater-Condon rules and not in the integrals
+  energy_2e_noL_HF = 0.5d0 * energy_2e_noL_HF
+
+  print*, "energy_2e_noL_HF = ", energy_2e_noL_HF
+
+END_PROVIDER
+
+! ---
+

src/bi_ort_ints/no_dressing_naive.irp.f

@@ -89,7 +89,7 @@ BEGIN_PROVIDER [double precision, noL_0e_naive]
   !$OMP END DO
   !$OMP END PARALLEL
 
-  noL_0e_naive = -1.d0 * (-sum(tmp)) / 6.d0
+  noL_0e_naive = -1.d0 * (sum(tmp)) / 6.d0
 
   deallocate(tmp)
 
@@ -182,9 +182,8 @@ BEGIN_PROVIDER [double precision, noL_1e_naive, (mo_num, mo_num)]
                                                 , j, sigma_j, s, sigma_s, i, sigma_i &
                                                 , I_pij_jsi)
 
-          ! x (-1) because integrals are over -L
           ! x 0.5  because we consider 0.5 (up + down)
-          noL_1e_naive(p,s) = noL_1e_naive(p,s) + 0.25d0 * (I_pij_sji - I_pij_sij + I_pij_jis - I_pij_ijs + I_pij_isj - I_pij_jsi)
+          noL_1e_naive(p,s) = noL_1e_naive(p,s) - 0.25d0 * (I_pij_sji - I_pij_sij + I_pij_jis - I_pij_ijs + I_pij_isj - I_pij_jsi)
         enddo ! j
       enddo ! i
     enddo ! s
@@ -254,9 +253,8 @@ BEGIN_PROVIDER [double precision, noL_1e_naive, (mo_num, mo_num)]
                                                 , j, sigma_j, s, sigma_s, i, sigma_i &
                                                 , I_pij_jsi)
 
-          ! x (-1) because integrals are over -L
           ! x 0.5  because we consider 0.5 (up + down)
-          noL_1e_naive(p,s) = noL_1e_naive(p,s) + 0.25d0 * (I_pij_sji - I_pij_sij + I_pij_jis - I_pij_ijs + I_pij_isj - I_pij_jsi)
+          noL_1e_naive(p,s) = noL_1e_naive(p,s) - 0.25d0 * (I_pij_sji - I_pij_sij + I_pij_jis - I_pij_ijs + I_pij_isj - I_pij_jsi)
         enddo ! j
       enddo ! i
     enddo ! s
@@ -335,9 +333,8 @@ BEGIN_PROVIDER [double precision, noL_2e_naive, (mo_num, mo_num, mo_num, mo_num)
                                                   , t, sigma_t, s, sigma_s, i, sigma_i &
                                                   , I_ipq_tsi)
 
-            ! x (-1) because integrals are over -L
             ! x 0.25  because we consider 0.25 (up-up + up-down + down-up + down-down)
-            noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) + 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
+            noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) - 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
           enddo ! i
         enddo ! p
       enddo ! q
@@ -389,9 +386,8 @@ BEGIN_PROVIDER [double precision, noL_2e_naive, (mo_num, mo_num, mo_num, mo_num)
                                                   , t, sigma_t, s, sigma_s, i, sigma_i &
                                                   , I_ipq_tsi)
 
-            ! x (-1) because integrals are over -L
             ! x 0.25  because we consider 0.25 (up-up + up-down + down-up + down-down)
-            noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) + 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
+            noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) - 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
           enddo ! i
         enddo ! p
       enddo ! q
@@ -443,9 +439,8 @@ BEGIN_PROVIDER [double precision, noL_2e_naive, (mo_num, mo_num, mo_num, mo_num)
                                                   , t, sigma_t, s, sigma_s, i, sigma_i &
                                                   , I_ipq_tsi)
 
-            ! x (-1) because integrals are over -L
             ! x 0.25  because we consider 0.25 (up-up + up-down + down-up + down-down)
-            noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) + 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
+            noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) - 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
           enddo ! i
         enddo ! p
       enddo ! q
@@ -497,9 +492,8 @@ BEGIN_PROVIDER [double precision, noL_2e_naive, (mo_num, mo_num, mo_num, mo_num)
                                                   , t, sigma_t, s, sigma_s, i, sigma_i &
                                                   , I_ipq_tsi)
 
-            ! x (-1) because integrals are over -L
             ! x 0.25  because we consider 0.25 (up-up + up-down + down-up + down-down)
-            noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) + 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
+            noL_2e_naive(p,q,s,t) = noL_2e_naive(p,q,s,t) - 0.125d0 * (I_ipq_ist - I_ipq_sit - I_ipq_tsi)
           enddo ! i
         enddo ! p
       enddo ! q

src/bi_ort_ints/no_dressing_v0.irp.f

@@ -40,7 +40,7 @@ BEGIN_PROVIDER [double precision, noL_0e]
     !$OMP END DO
     !$OMP END PARALLEL
 
-    noL_0e = -1.d0 * (-sum(tmp)) / 6.d0
+    noL_0e = -1.d0 * (sum(tmp)) / 6.d0
 
     deallocate(tmp)
 
@@ -114,7 +114,7 @@ BEGIN_PROVIDER [double precision, noL_0e]
     !$OMP END DO
     !$OMP END PARALLEL
 
-    noL_0e = -1.d0 * (-sum(tmp)) / 6.d0
+    noL_0e = -1.d0 * (sum(tmp)) / 6.d0
 
     deallocate(tmp)
 
@@ -131,12 +131,6 @@ END_PROVIDER
 
 BEGIN_PROVIDER [double precision, noL_1e, (mo_num, mo_num)]
 
-  BEGIN_DOC
-  !
-  ! x (-1) because integrals are over -L
-  !
-  END_DOC
-
   implicit none
   integer          :: p, s, i, j
   double precision :: I_pij_sij, I_pij_isj, I_pij_ijs, I_pij_sji, I_pij_jsi, I_pij_jis
@@ -167,7 +161,7 @@ BEGIN_PROVIDER [double precision, noL_1e, (mo_num, mo_num)]
             call give_integrals_3_body_bi_ort(p, i, j, i, j, s, I_pij_ijs)
             call give_integrals_3_body_bi_ort(p, i, j, s, j, i, I_pij_sji)
           
-            noL_1e(p,s) = noL_1e(p,s) - (2.d0*I_pij_sij - 2.d0*I_pij_isj + I_pij_ijs - I_pij_sji)
+            noL_1e(p,s) = noL_1e(p,s) + (2.d0*I_pij_sij - 2.d0*I_pij_isj + I_pij_ijs - I_pij_sji)
           enddo
         enddo
       enddo
@@ -197,7 +191,7 @@ BEGIN_PROVIDER [double precision, noL_1e, (mo_num, mo_num)]
             call give_integrals_3_body_bi_ort(p, i, j, i, j, s, I_pij_ijs)
             call give_integrals_3_body_bi_ort(p, i, j, s, j, i, I_pij_sji)
           
-            noL_1e(p,s) = noL_1e(p,s) - (2.d0*I_pij_sij - 2.d0*I_pij_isj + I_pij_ijs - I_pij_sji)
+            noL_1e(p,s) = noL_1e(p,s) + (2.d0*I_pij_sij - 2.d0*I_pij_isj + I_pij_ijs - I_pij_sji)
           enddo ! j
         enddo ! i
 
@@ -211,7 +205,7 @@ BEGIN_PROVIDER [double precision, noL_1e, (mo_num, mo_num)]
             call give_integrals_3_body_bi_ort(p, i, j, i, s, j, I_pij_isj)
             call give_integrals_3_body_bi_ort(p, i, j, i, j, s, I_pij_ijs)
           
-            noL_1e(p,s) = noL_1e(p,s) + 0.5d0 * (2.d0*I_pij_sji - I_pij_jsi + 2.d0*I_pij_jis - 4.d0*I_pij_sij + 2.d0*I_pij_isj - I_pij_ijs)
+            noL_1e(p,s) = noL_1e(p,s) - 0.5d0 * (2.d0*I_pij_sji - I_pij_jsi + 2.d0*I_pij_jis - 4.d0*I_pij_sij + 2.d0*I_pij_isj - I_pij_ijs)
           enddo ! j
 
           do j = elec_beta_num+1, elec_alpha_num
@@ -221,7 +215,7 @@ BEGIN_PROVIDER [double precision, noL_1e, (mo_num, mo_num)]
             call give_integrals_3_body_bi_ort(p, i, j, i, j, s, I_pij_ijs)
             call give_integrals_3_body_bi_ort(p, i, j, s, j, i, I_pij_sji)
           
-            noL_1e(p,s) = noL_1e(p,s) - 0.5d0 * (I_pij_sij - I_pij_isj + I_pij_ijs - I_pij_sji)
+            noL_1e(p,s) = noL_1e(p,s) + 0.5d0 * (I_pij_sij - I_pij_isj + I_pij_ijs - I_pij_sji)
           enddo ! j
         enddo ! i
 
@@ -241,12 +235,6 @@ END_PROVIDER
 
 BEGIN_PROVIDER [double precision, noL_2e, (mo_num, mo_num, mo_num, mo_num)]
 
-  BEGIN_DOC
-  !
-  ! x (-1) because integrals are over -L
-  !
-  END_DOC
-
   implicit none
   integer          :: p, q, s, t, i
   double precision :: I_ipq_sit, I_ipq_tsi, I_ipq_ist
@@ -276,7 +264,7 @@ BEGIN_PROVIDER [double precision, noL_2e, (mo_num, mo_num, mo_num, mo_num)]
               call give_integrals_3_body_bi_ort(i, p, q, t, s, i, I_ipq_tsi)
               call give_integrals_3_body_bi_ort(i, p, q, i, s, t, I_ipq_ist)
           
-              noL_2e(p,q,s,t) = noL_2e(p,q,s,t) - 0.5d0 * (I_ipq_sit + I_ipq_tsi - 2.d0*I_ipq_ist)
+              noL_2e(p,q,s,t) = noL_2e(p,q,s,t) + 0.5d0 * (I_ipq_sit + I_ipq_tsi - 2.d0*I_ipq_ist)
             enddo
           enddo
         enddo
@@ -306,7 +294,7 @@ BEGIN_PROVIDER [double precision, noL_2e, (mo_num, mo_num, mo_num, mo_num)]
               call give_integrals_3_body_bi_ort(i, p, q, t, s, i, I_ipq_tsi)
               call give_integrals_3_body_bi_ort(i, p, q, i, s, t, I_ipq_ist)
             
-              noL_2e(p,q,s,t) = noL_2e(p,q,s,t) - 0.5d0 * (I_ipq_sit + I_ipq_tsi - 2.d0*I_ipq_ist)
+              noL_2e(p,q,s,t) = noL_2e(p,q,s,t) + 0.5d0 * (I_ipq_sit + I_ipq_tsi - 2.d0*I_ipq_ist)
             enddo ! i
 
             do i = elec_beta_num+1, elec_alpha_num
@@ -315,7 +303,7 @@ BEGIN_PROVIDER [double precision, noL_2e, (mo_num, mo_num, mo_num, mo_num)]
               call give_integrals_3_body_bi_ort(i, p, q, t, s, i, I_ipq_tsi)
               call give_integrals_3_body_bi_ort(i, p, q, i, s, t, I_ipq_ist)
             
-              noL_2e(p,q,s,t) = noL_2e(p,q,s,t) - 0.25d0 * (I_ipq_sit + I_ipq_tsi - 2.d0*I_ipq_ist)
+              noL_2e(p,q,s,t) = noL_2e(p,q,s,t) + 0.25d0 * (I_ipq_sit + I_ipq_tsi - 2.d0*I_ipq_ist)
             enddo ! i
 
           enddo ! p

src/bi_ort_ints/one_e_bi_ort.irp.f

@@ -53,11 +53,13 @@ END_PROVIDER
  BEGIN_PROVIDER [double precision, mo_bi_orth_bipole_x , (mo_num,mo_num)]
 &BEGIN_PROVIDER [double precision, mo_bi_orth_bipole_y , (mo_num,mo_num)]
 &BEGIN_PROVIDER [double precision, mo_bi_orth_bipole_z , (mo_num,mo_num)]
+
   BEGIN_DOC
   ! array of the integrals of Left MO_i * x Right MO_j
   ! array of the integrals of Left MO_i * y Right MO_j
   ! array of the integrals of Left MO_i * z Right MO_j
   END_DOC
+ 
   implicit none
 
   call ao_to_mo_bi_ortho(                                                     &

src/bi_ort_ints/three_body_ints_bi_ort.irp.f

@@ -126,7 +126,7 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
 
   BEGIN_DOC
   !
-  ! < n l k | -L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS
+  ! < n l k | L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS
   !
   END_DOC
 

src/bi_ort_ints/total_twoe_pot.irp.f

@@ -258,7 +258,8 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e, (mo_num, mo_num, mo_num,
 
   if(noL_standard) then
     PROVIDE noL_2e
-    mo_bi_ortho_tc_two_e = mo_bi_ortho_tc_two_e + noL_2e
+    ! x 2 because of the Slater-Condon rules convention
+    mo_bi_ortho_tc_two_e = mo_bi_ortho_tc_two_e + 2.d0 * noL_2e
     FREE noL_2e
   endif
 
@@ -272,9 +273,11 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti,     (mo_num,mo_num)]
 
   BEGIN_DOC
+  !
   ! mo_bi_ortho_tc_two_e_jj         (i,j) = J_ij = <ji|W-K|ji>
   ! mo_bi_ortho_tc_two_e_jj_exchange(i,j) = K_ij = <ij|W-K|ji>
   ! mo_bi_ortho_tc_two_e_jj_anti    (i,j) = J_ij - K_ij
+  !
   END_DOC
 
   implicit none

src/tc_bi_ortho/slater_tc_opt.irp.f

@@ -90,8 +90,11 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
   hthree = 0.d0
 
   call get_excitation_degree(key_i, key_j, degree, Nint)
+
   if(.not.pure_three_body_h_tc) then
+
     if(degree .gt. 2) return
+    
     if(degree == 0) then
       call diag_htilde_mu_mat_fock_bi_ortho  (Nint, key_i, hmono, htwoe, hthree, htot)
     else if (degree == 1) then
@@ -99,8 +102,11 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
     else if(degree == 2) then
       call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
     endif
+
   else 
+
     if(degree .gt. 3) return
+
     if(degree == 0) then
       call diag_htilde_mu_mat_fock_bi_ortho  (Nint, key_i, hmono, htwoe, hthree, htot)
     else if (degree == 1) then
@@ -110,6 +116,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
     else
       call triple_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
     endif
+
   endif
 
   if(degree==0) then
@@ -161,3 +168,4 @@ subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot)
 end
 
 ! ---
+

src/tc_bi_ortho/slater_tc_opt_diag.irp.f

@@ -7,7 +7,9 @@
 &BEGIN_PROVIDER [ double precision, ref_tc_energy_3e]
 
   BEGIN_DOC
+  !
   ! Various component of the TC energy for the reference "HF" Slater determinant
+  !
   END_DOC 
 
   implicit none
@@ -41,7 +43,9 @@ END_PROVIDER
 subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree, htot)
 
   BEGIN_DOC
+  !
   ! Computes $\langle i|H|i \rangle$.
+  !
   END_DOC
 
   implicit none
@@ -124,6 +128,7 @@ end
 subroutine ac_tc_operator(iorb, ispin, key, hmono, htwoe, hthree, Nint, na, nb)
 
   BEGIN_DOC
+  !
   ! Routine that computes one- and two-body energy corresponding 
   ! 
   ! to the ADDITION of an electron in an orbital 'iorb' of spin 'ispin' 
@@ -133,6 +138,7 @@ subroutine ac_tc_operator(iorb, ispin, key, hmono, htwoe, hthree, Nint, na, nb)
   ! in output, the determinant key is changed by the ADDITION of that electron 
   !
   ! and the quantities hmono,htwoe,hthree are INCREMENTED 
+  !
   END_DOC
 
   use bitmasks
@@ -188,8 +194,8 @@ subroutine ac_tc_operator(iorb, ispin, key, hmono, htwoe, hthree, Nint, na, nb)
   enddo
 
   if(three_body_h_tc .and. (elec_num.gt.2) .and. three_e_3_idx_term) then
-
     !!!!! 3-e part 
+
     !! same-spin/same-spin
     do j = 1, na
       jj = occ(j,ispin)
@@ -227,9 +233,12 @@ end
 ! ---
 
 subroutine a_tc_operator(iorb, ispin, key, hmono, htwoe, hthree, Nint, na, nb)
+
   use bitmasks
   implicit none
+
   BEGIN_DOC
+  !
   ! Routine that computes one- and two-body energy corresponding 
   ! 
   ! to the REMOVAL of an electron in an orbital 'iorb' of spin 'ispin' 
@@ -239,7 +248,9 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
   ! in output, the determinant key is changed by the REMOVAL of that electron 
   !
   ! and the quantities hmono,htwoe,hthree are INCREMENTED 
+  !
   END_DOC
+
   integer,           intent(in)    :: iorb, ispin, Nint
   integer,           intent(inout) :: na, nb
   integer(bit_kind), intent(inout) :: key(Nint,2)
@@ -280,6 +291,7 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
 
   if(three_body_h_tc .and. elec_num.gt.2 .and. three_e_3_idx_term) then
     !!!!! 3-e part 
+
     !! same-spin/same-spin
     do j = 1, na
       jj = occ(j,ispin)
@@ -312,17 +324,19 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
 
 end
 
+! ---
 
 subroutine diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, det_in,htot)
-  implicit none
+
   BEGIN_DOC
   ! Computes $\langle i|H|i \rangle$. WITHOUT ANY CONTRIBUTIONS FROM 3E TERMS
   END_DOC
+
+  implicit none
   integer,           intent(in)  :: Nint
   integer(bit_kind), intent(in)  :: det_in(Nint,2)
   double precision,  intent(out) :: htot
   double precision               :: hmono, htwoe
-
   integer(bit_kind)              :: hole(Nint,2)
   integer(bit_kind)              :: particle(Nint,2)
   integer                        :: i, nexc(2), ispin
@@ -367,8 +381,8 @@ subroutine diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, det_in,htot)
   ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha
   ASSERT (tmp(2) == nexc(2)) ! Number of holes beta 
 
-  
   det_tmp = ref_bitmask
+  
   hmono = ref_tc_energy_1e
   htwoe = ref_tc_energy_2e 
   do ispin=1,2

src/tc_bi_ortho/slater_tc_opt_double.irp.f

@@ -1,4 +1,6 @@
 
+! ---
+
 subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
 
   BEGIN_DOC
@@ -32,6 +34,7 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
   if(degree .ne. 2) then
     return
   endif
+
   integer :: degree_i, degree_j
   call get_excitation_degree(ref_bitmask, key_i, degree_i, N_int)
   call get_excitation_degree(ref_bitmask, key_j, degree_j, N_int)
@@ -40,44 +43,65 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
 
   if(s1 .ne. s2) then
     ! opposite spin two-body
+
     htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
-    if(three_body_h_tc.and.elec_num.gt.2)then
+
+    if(three_body_h_tc .and. (elec_num .gt. 2)) then
+      ! add 3-e term
+
       if(.not.double_normal_ord .and. three_e_5_idx_term) then
+        ! 5-idx approx
+
         if(degree_i > degree_j) then
           call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree)
         else
           call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
         endif
+
       elseif(double_normal_ord) then
+        ! noL a la Manu
+
         htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)
       endif
     endif
+
   else
     ! same spin two-body
+
     ! direct terms
     htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
+
     ! exchange terms
     htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1)
-   if(three_body_h_tc.and.elec_num.gt.2)then
+
+    if(three_body_h_tc .and. (elec_num .gt. 2)) then
+      ! add 3-e term
+
       if(.not.double_normal_ord.and.three_e_5_idx_term)then
+        ! 5-idx approx
+
         if(degree_i > degree_j) then
           call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree)
         else
           call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
         endif
+
       elseif(double_normal_ord) then
+        ! noL a la Manu
+
         htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)
         htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)
       endif
     endif
   endif
+
   hthree *= phase
   htwoe  *= phase
   htot    = htwoe + hthree
 
 end
 
-
+! ---
 
 subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
  implicit none

src/tc_bi_ortho/slater_tc_opt_single.irp.f

@@ -1,12 +1,16 @@
 
+! ---
 
 subroutine  single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
+
   BEGIN_DOC
+  !
   ! <key_j |H_tilde | key_i> for single excitation ONLY FOR ONE- AND TWO-BODY TERMS 
   !!
   !! WARNING !!
   ! 
   ! Non hermitian !!
+  !
   END_DOC
 
   use bitmasks
@@ -31,25 +35,31 @@ subroutine  single_htilde_mu_mat_fock_bi_ortho (Nint, key_j, key_i, hmono, htwoe
   htwoe  = 0.d0
   hthree = 0.d0
   htot   = 0.d0
+
   call get_excitation_degree(key_i, key_j, degree, Nint)
   if(degree .ne. 1) then
     return
   endif
-  call bitstring_to_list_ab(key_i, occ, Ne, Nint)
 
+  call bitstring_to_list_ab(key_i, occ, Ne, Nint)
   call get_single_excitation(key_i, key_j, exc, phase, Nint)
   call decode_exc(exc, 1, h1, p1, h2, p2, s1, s2)
   call get_single_excitation_from_fock_tc(key_i, key_j, h1, p1, s1, phase, hmono, htwoe, hthree, htot)
+
 end
 
+! ---
 
 subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, hmono, htwoe, hthree, htot)
+
   use bitmasks
+
   implicit none
   integer,           intent(in)  :: h, p, spin
   double precision,  intent(in)  :: phase
   integer(bit_kind), intent(in)  :: key_i(N_int,2), key_j(N_int,2)
   double precision,  intent(out) :: hmono, htwoe, hthree, htot
+
   integer(bit_kind)              :: differences(N_int,2)
   integer(bit_kind)              :: hole(N_int,2)
   integer(bit_kind)              :: partcl(N_int,2)
@@ -58,10 +68,12 @@ subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,h
   integer                        :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
   integer                        :: i0,i
   double precision               :: buffer_c(mo_num),buffer_x(mo_num)
+
   do i = 1, mo_num
     buffer_c(i) = tc_2e_3idx_coulomb_integrals (i,p,h)
     buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
   enddo
+
   do i = 1, N_int
     differences(i,1) = xor(key_i(i,1), ref_closed_shell_bitmask(i,1))
     differences(i,2) = xor(key_i(i,2), ref_closed_shell_bitmask(i,2))
@@ -70,10 +82,12 @@ subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,h
     partcl     (i,1) = iand(differences(i,1), key_i(i,1))
     partcl     (i,2) = iand(differences(i,2), key_i(i,2))
   enddo
+
   call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int)
   call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int)
   hmono = mo_bi_ortho_tc_one_e(p,h)
   htwoe = fock_op_2_e_tc_closed_shell(p,h)
+
   ! holes :: direct terms
   do i0 = 1, n_occ_ab_hole(1)
     i = occ_hole(i0,1)
@@ -105,12 +119,12 @@ subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,h
     i = occ_partcl(i0,spin)
     htwoe -= buffer_x(i)
   enddo
+
   hthree = 0.d0
   if (three_body_h_tc .and. elec_num.gt.2 .and. three_e_4_idx_term) then
     call three_comp_fock_elem(key_i, h, p, spin, hthree)
   endif
 
-
   htwoe  = htwoe * phase
   hmono  = hmono * phase
   hthree = hthree * phase
@@ -118,6 +132,8 @@ subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,h
 
 end
 
+! ---
+
 subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree)
  implicit none
  integer,intent(in) :: h_fock,p_fock,ispin_fock

src/tc_bi_ortho/slater_tc_slow.irp.f

@@ -81,8 +81,14 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree,
   endif
 
   htot = hmono + htwoe + hthree
+
   if(degree==0) then
     htot += nuclear_repulsion
+
+    if(noL_standard) then
+      PROVIDE noL_0e
+      htot += noL_0e
+    endif
   endif
  
 end
@@ -92,7 +98,9 @@ end
 subroutine diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
 
   BEGIN_DOC
+  !
   ! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS 
+  !
   END_DOC
 
   use bitmasks
@@ -108,45 +116,19 @@ subroutine diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
 
   PROVIDE mo_bi_ortho_tc_two_e
 
-!  PROVIDE mo_two_e_integrals_tc_int_in_map mo_bi_ortho_tc_two_e
-!
-!  PROVIDE mo_integrals_erf_map core_energy nuclear_repulsion core_bitmask
-!  PROVIDE core_fock_operator
-!
-!  PROVIDE j1b_gauss
-
-!  if(core_tc_op)then
-!   print*,'core_tc_op not already taken into account for bi ortho'
-!   print*,'stopping ...'
-!   stop
-!   do i = 1, Nint
-!    key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1))
-!    key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2))
-!   enddo
-!   call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
-!   hmono = core_energy - nuclear_repulsion
-!  else
-   call bitstring_to_list_ab(key_i, occ, Ne, Nint)
   hmono = 0.d0
-!  endif
   htwoe = 0.d0
   htot  = 0.d0
 
+  call bitstring_to_list_ab(key_i, occ, Ne, Nint)
+
   do ispin = 1, 2 
-   do i = 1, Ne(ispin) ! 
+    do i = 1, Ne(ispin)
       ii = occ(i,ispin) 
       hmono += mo_bi_ortho_tc_one_e(ii,ii)
-
-!    if(core_tc_op)then
-!   print*,'core_tc_op not already taken into account for bi ortho'
-!   print*,'stopping ...'
-!   stop
-!     hmono += core_fock_operator(ii,ii) ! add the usual Coulomb - Exchange from the core 
-!    endif
     enddo
   enddo
 
-
   ! alpha/beta two-body
   ispin = 1
   jspin = 2 
@@ -175,11 +157,12 @@ subroutine diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
       htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
     enddo
   enddo
+
   htot = hmono + htwoe 
 
 end
 
-
+! ---
 
 subroutine double_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
 

src/tc_bi_ortho/tc_hmat.irp.f

@@ -1,10 +1,14 @@
 
+! ---
+
 BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)]
 
   BEGIN_DOC
+  !
   ! htilde_matrix_elmt_bi_ortho(j,i) = <J| H^tilde |I> 
   !
   ! WARNING !!!!!!!!! IT IS NOT HERMITIAN !!!!!!!!!
+  !
   END_DOC
  
   implicit none
@@ -34,11 +38,16 @@ END_PROVIDER
 ! ---
 
 BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho_tranp, (N_det,N_det)]
+
   implicit none
   integer ::i,j
+
   do i = 1, N_det
     do j = 1, N_det
       htilde_matrix_elmt_bi_ortho_tranp(j,i) = htilde_matrix_elmt_bi_ortho(i,j)
     enddo
   enddo
 END_PROVIDER 
+
+! ---
+

src/tc_bi_ortho/test_tc_bi_ortho.irp.f

@@ -557,6 +557,8 @@ subroutine test_no_1()
 
   print*, ' accu (%) = ', 100.d0*accu/norm
 
+  PROVIDE energy_1e_noL_HF
+
   return
 end
 
@@ -572,6 +574,7 @@ subroutine test_no_2()
 
   PROVIDE noL_2e_naive
   PROVIDE noL_2e
+  PROVIDE energy_2e_noL_HF
 
   thr = 1d-8