diff --git a/src/tc_bi_ortho/EZFIO.cfg b/src/tc_bi_ortho/EZFIO.cfg
new file mode 100644
index 00000000..a34d2134
--- /dev/null
+++ b/src/tc_bi_ortho/EZFIO.cfg
@@ -0,0 +1,11 @@
+[psi_l_coef_bi_ortho]
+interface: ezfio
+doc: Coefficients for the left wave function
+type: double precision
+size: (determinants.n_det,determinants.n_states)
+
+[psi_r_coef_bi_ortho]
+interface: ezfio
+doc: Coefficients for the right wave function
+type: double precision
+size: (determinants.n_det,determinants.n_states)
diff --git a/src/tc_bi_ortho/NEED b/src/tc_bi_ortho/NEED
new file mode 100644
index 00000000..9a0c20ef
--- /dev/null
+++ b/src/tc_bi_ortho/NEED
@@ -0,0 +1,6 @@
+bi_ort_ints
+bi_ortho_mos
+tc_keywords
+non_hermit_dav
+dav_general_mat
+tc_scf
diff --git a/src/tc_bi_ortho/e_corr_bi_ortho.irp.f b/src/tc_bi_ortho/e_corr_bi_ortho.irp.f
new file mode 100644
index 00000000..ec66a8b5
--- /dev/null
+++ b/src/tc_bi_ortho/e_corr_bi_ortho.irp.f
@@ -0,0 +1,104 @@
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+ BEGIN_PROVIDER [ double precision, e_tilde_00]
+ implicit none
+ double precision :: hmono,htwoe,hthree,htot
+ call htilde_mu_mat_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,htot)
+ e_tilde_00 = htot
+ END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, e_pt2_tc_bi_orth]
+&BEGIN_PROVIDER [ double precision, e_pt2_tc_bi_orth_single]
+&BEGIN_PROVIDER [ double precision, e_pt2_tc_bi_orth_double]
+ implicit none
+ integer :: i,degree
+ double precision :: hmono,htwoe,hthree,htilde_ij,coef_pt1,e_i0,delta_e
+ e_pt2_tc_bi_orth = 0.d0
+ e_pt2_tc_bi_orth_single = 0.d0
+ e_pt2_tc_bi_orth_double = 0.d0
+ do i = 1, N_det
+ call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
+ if(degree == 1 .or. degree == 2)then
+ call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+ call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
+ delta_e = e_tilde_00 - e_i0
+ coef_pt1 = htilde_ij / delta_e
+ call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
+ e_pt2_tc_bi_orth += coef_pt1 * htilde_ij
+ if(degree == 1)then
+ e_pt2_tc_bi_orth_single += coef_pt1 * htilde_ij
+ else
+! print*,'coef_pt1, e_pt2',coef_pt1,coef_pt1 * htilde_ij
+ e_pt2_tc_bi_orth_double += coef_pt1 * htilde_ij
+ endif
+ endif
+ enddo
+ END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, e_tilde_bi_orth_00]
+ implicit none
+ double precision :: hmono,htwoe,hthree,htilde_ij
+ call htilde_mu_mat_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,e_tilde_bi_orth_00)
+ e_tilde_bi_orth_00 += nuclear_repulsion
+ END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, e_corr_bi_orth ]
+&BEGIN_PROVIDER [ double precision, e_corr_bi_orth_proj ]
+&BEGIN_PROVIDER [ double precision, e_corr_single_bi_orth ]
+&BEGIN_PROVIDER [ double precision, e_corr_double_bi_orth ]
+ implicit none
+ integer :: i,degree
+ double precision :: hmono,htwoe,hthree,htilde_ij
+
+ e_corr_bi_orth = 0.d0
+ e_corr_single_bi_orth = 0.d0
+ e_corr_double_bi_orth = 0.d0
+ do i = 1, N_det
+ call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
+ call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
+ if(degree == 1)then
+ e_corr_single_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
+ else if(degree == 2)then
+ e_corr_double_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
+! print*,'coef_wf , e_cor',reigvec_tc_bi_orth(i,1)/reigvec_tc_bi_orth(1,1), reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
+ endif
+ enddo
+ e_corr_bi_orth_proj = e_corr_single_bi_orth + e_corr_double_bi_orth
+ e_corr_bi_orth = eigval_right_tc_bi_orth(1) - e_tilde_bi_orth_00
+ END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, e_tc_left_right ]
+ implicit none
+ integer :: i,j
+ double precision :: hmono,htwoe,hthree,htilde_ij,accu
+ e_tc_left_right = 0.d0
+ accu = 0.d0
+ do i = 1, N_det
+ accu += reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(i,1)
+ do j = 1, N_det
+ call htilde_mu_mat_bi_ortho(psi_det(1,1,j),psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
+ e_tc_left_right += htilde_ij * reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(j,1)
+ enddo
+ enddo
+ e_tc_left_right *= 1.d0/accu
+ e_tc_left_right += nuclear_repulsion
+
+ END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, coef_pt1_bi_ortho, (N_det)]
+ implicit none
+ integer :: i,degree
+ double precision :: hmono,htwoe,hthree,htilde_ij,coef_pt1,e_i0,delta_e
+ do i = 1, N_det
+ call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
+ if(degree==0)then
+ coef_pt1_bi_ortho(i) = 1.d0
+ else
+ call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+ call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
+ delta_e = e_tilde_00 - e_i0
+ coef_pt1 = htilde_ij / delta_e
+ coef_pt1_bi_ortho(i)= coef_pt1
+ endif
+ enddo
+END_PROVIDER
diff --git a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f
new file mode 100644
index 00000000..b7129d36
--- /dev/null
+++ b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f
@@ -0,0 +1,92 @@
+subroutine htc_bi_ortho_calc_tdav(v, u, N_st, sze)
+
+ use bitmasks
+
+ BEGIN_DOC
+ ! Application of H_TC on a vector
+ !
+ ! v(i,istate) = \sum_j u(j,istate) H_TC(i,j), with:
+ ! H_TC(i,j) = < Di | H_TC | Dj >
+ !
+ END_DOC
+
+ implicit none
+
+ integer, intent(in) :: N_st, sze
+ double precision, intent(in) :: u(sze,N_st)
+ double precision, intent(inout) :: v(sze,N_st)
+
+ integer :: i, j, istate
+ double precision :: htot
+
+ PROVIDE N_int
+ PROVIDE psi_det
+
+
+ ! TODO : transform it with the bi-linear representation in terms of alpha-beta.
+
+ i = 1
+ j = 1
+ call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
+
+ v = 0.d0
+ !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
+ !$OMP SHARED(N_st, sze, N_int, psi_det, u, v) &
+ !$OMP PRIVATE(istate, i, j, htot)
+ do istate = 1, N_st
+ do i = 1, sze
+ do j = 1, sze
+ call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
+ v(i,istate) = v(i,istate) + htot * u(j,istate)
+ enddo
+ enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+end
+
+subroutine htcdag_bi_ortho_calc_tdav(v, u, N_st, sze)
+
+ use bitmasks
+
+ BEGIN_DOC
+ ! Application of (H_TC)^dagger on a vector
+ !
+ ! v(i,istate) = \sum_j u(j,istate) H_TC(j,i), with:
+ ! H_TC(i,j) = < Di | H_TC | Dj >
+ !
+ END_DOC
+
+ implicit none
+
+ integer, intent(in) :: N_st, sze
+ double precision, intent(in) :: u(sze,N_st)
+ double precision, intent(inout) :: v(sze,N_st)
+
+ integer :: i, j, istate
+ double precision :: htot
+
+ PROVIDE N_int
+ PROVIDE psi_det
+
+ i = 1
+ j = 1
+ call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
+
+ v = 0.d0
+
+ !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
+ !$OMP SHARED(N_st, sze, N_int, psi_det, u, v) &
+ !$OMP PRIVATE(istate, i, j, htot)
+ do istate = 1, N_st
+ do i = 1, sze
+ do j = 1, sze
+ call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot)
+ v(i,istate) = v(i,istate) + htot * u(j,istate)
+ enddo
+ enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+end
+
diff --git a/src/tc_bi_ortho/normal_ordered.irp.f b/src/tc_bi_ortho/normal_ordered.irp.f
new file mode 100644
index 00000000..81f5fb2c
--- /dev/null
+++ b/src/tc_bi_ortho/normal_ordered.irp.f
@@ -0,0 +1,319 @@
+BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_num, mo_num)]
+
+ BEGIN_DOC
+ ! Normal ordering of the three body interaction on the HF density
+ END_DOC
+
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+
+ implicit none
+
+ integer :: i,h1,p1,h2,p2
+ integer :: hh1,hh2,pp1,pp2
+ integer :: Ne(2)
+ integer, allocatable :: occ(:,:)
+ integer(bit_kind), allocatable :: key_i_core(:,:)
+ double precision :: hthree_aba,hthree_aaa,hthree_aab
+ double precision :: wall0,wall1
+
+ PROVIDE N_int
+
+ allocate( occ(N_int*bit_kind_size,2) )
+ allocate( key_i_core(N_int,2) )
+
+ if(core_tc_op) then
+ do i = 1, N_int
+ key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
+ key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
+ enddo
+ call bitstring_to_list_ab(key_i_core,occ,Ne,N_int)
+ else
+ call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
+ endif
+
+ normal_two_body_bi_orth = 0.d0
+ print*,'Providing normal_two_body_bi_orth ...'
+ call wall_time(wall0)
+
+ !$OMP PARALLEL &
+ !$OMP DEFAULT (NONE) &
+ !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aba, hthree_aab, hthree_aaa) &
+ !$OMP SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth)
+ !$OMP DO SCHEDULE (static)
+ do hh1 = 1, n_act_orb
+ h1 = list_act(hh1)
+ do pp1 = 1, n_act_orb
+ p1 = list_act(pp1)
+ do hh2 = 1, n_act_orb
+ h2 = list_act(hh2)
+ do pp2 = 1, n_act_orb
+ p2 = list_act(pp2)
+ ! opposite spin double excitations
+ call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aba)
+ ! same spin double excitations with opposite spin contributions
+ if(h1
h2
+ ! same spin double excitations with same spin contributions
+ if(Ne(2).ge.3)then
+ call give_aaa_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aaa) ! exchange h1<->h2
+ else
+ hthree_aaa = 0.d0
+ endif
+ else
+ call give_aab_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aab)
+ if(Ne(2).ge.3)then
+ call give_aaa_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aaa)
+ else
+ hthree_aaa = 0.d0
+ endif
+ endif
+ normal_two_body_bi_orth(p2,h2,p1,h1) = 0.5d0*(hthree_aba + hthree_aab + hthree_aaa)
+ enddo
+ enddo
+ enddo
+ enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+ call wall_time(wall1)
+ print*,'Wall time for normal_two_body_bi_orth ',wall1-wall0
+
+ deallocate( occ )
+ deallocate( key_i_core )
+
+END_PROVIDER
+
+
+
+subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
+
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+
+ implicit none
+ integer, intent(in) :: Nint, h1, h2, p1, p2
+ integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
+ double precision, intent(out) :: hthree
+ integer :: ii, i
+ double precision :: int_direct, int_exc_12, int_exc_13, integral
+
+ !!!! double alpha/beta
+ hthree = 0.d0
+ do ii = 1, Ne(2) ! purely closed shell part
+ i = occ(ii,2)
+ call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
+ int_direct = -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
+ int_exc_13 = -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p2, i,p1,i,h2,h1,integral)
+ int_exc_12 = -1.d0 * integral
+ hthree += 2.d0 * int_direct - 1.d0 * ( int_exc_13 + int_exc_12)
+ enddo
+ do ii = Ne(2) + 1, Ne(1) ! purely open-shell part
+ i = occ(ii,1)
+ call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
+ int_direct = -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
+ int_exc_13 = -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p2, i,p1,i,h2,h1,integral)
+ int_exc_12 = -1.d0 * integral
+ hthree += 1.d0 * int_direct - 0.5d0* ( int_exc_13 + int_exc_12)
+ enddo
+
+end subroutine give_aba_contraction
+
+
+
+BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num, mo_num, mo_num)]
+
+ BEGIN_DOC
+ ! Normal ordered two-body sector of the three-body terms for opposite spin double excitations
+ END_DOC
+
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+
+ implicit none
+ integer :: h1, p1, h2, p2, i
+ integer :: hh1, hh2, pp1, pp2
+ integer :: Ne(2)
+ integer, allocatable :: occ(:,:)
+ integer(bit_kind), allocatable :: key_i_core(:,:)
+ double precision :: hthree
+
+ PROVIDE N_int
+
+ allocate( key_i_core(N_int,2) )
+ allocate( occ(N_int*bit_kind_size,2) )
+
+ if(core_tc_op)then
+ do i = 1, N_int
+ key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
+ key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
+ enddo
+ call bitstring_to_list_ab(key_i_core,occ,Ne,N_int)
+ else
+ call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
+ endif
+ normal_two_body_bi_orth_ab = 0.d0
+ do hh1 = 1, n_act_orb
+ h1 = list_act(hh1)
+ do pp1 = 1, n_act_orb
+ p1 = list_act(pp1)
+ do hh2 = 1, n_act_orb
+ h2 = list_act(hh2)
+ do pp2 = 1, n_act_orb
+ p2 = list_act(pp2)
+ call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree)
+ normal_two_body_bi_orth_ab(p2,h2,p1,h1) = hthree
+ enddo
+ enddo
+ enddo
+ enddo
+
+ deallocate( key_i_core )
+ deallocate( occ )
+
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_aa_bb, (n_act_orb, n_act_orb, n_act_orb, n_act_orb)]
+
+ BEGIN_DOC
+ ! Normal ordered two-body sector of the three-body terms for same spin double excitations
+ END_DOC
+
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+
+ implicit none
+ integer :: i,ii,j,h1,p1,h2,p2
+ integer :: hh1,hh2,pp1,pp2
+ integer :: Ne(2)
+ integer, allocatable :: occ(:,:)
+ integer(bit_kind), allocatable :: key_i_core(:,:)
+ double precision :: hthree_aab, hthree_aaa
+
+ PROVIDE N_int
+
+ allocate( key_i_core(N_int,2) )
+ allocate( occ(N_int*bit_kind_size,2) )
+
+ if(core_tc_op)then
+ do i = 1, N_int
+ key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
+ key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
+ enddo
+ call bitstring_to_list_ab(key_i_core, occ, Ne, N_int)
+ else
+ call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
+ endif
+
+ normal_two_body_bi_orth_aa_bb = 0.d0
+ do hh1 = 1, n_act_orb
+ h1 = list_act(hh1)
+ do pp1 = 1 , n_act_orb
+ p1 = list_act(pp1)
+ do hh2 = 1, n_act_orb
+ h2 = list_act(hh2)
+ do pp2 = 1 , n_act_orb
+ p2 = list_act(pp2)
+ if(h1h2
+ if(Ne(2).ge.3)then
+ call give_aaa_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aaa) ! exchange h1<->h2
+ else
+ hthree_aaa = 0.d0
+ endif
+ else
+ call give_aab_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aab)
+ if(Ne(2).ge.3)then
+ call give_aaa_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aaa)
+ else
+ hthree_aaa = 0.d0
+ endif
+ endif
+ normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) = hthree_aab + hthree_aaa
+ enddo
+ enddo
+ enddo
+ enddo
+
+ deallocate( key_i_core )
+ deallocate( occ )
+
+END_PROVIDER
+
+
+
+subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
+
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+
+ implicit none
+ integer, intent(in) :: Nint, h1, h2, p1, p2
+ integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
+ double precision, intent(out) :: hthree
+ integer :: ii,i
+ double precision :: int_direct,int_exc_12,int_exc_13,int_exc_23
+ double precision :: integral,int_exc_l,int_exc_ll
+
+ hthree = 0.d0
+ do ii = 1, Ne(2) ! purely closed shell part
+ i = occ(ii,2)
+ call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
+ int_direct = -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p2,p1,i ,i,h2,h1,integral)
+ int_exc_l = -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p1,i ,p2,i,h2,h1,integral)
+ int_exc_ll= -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p2,i ,p1,i,h2,h1,integral)
+ int_exc_12= -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
+ int_exc_13= -1.d0 * integral
+ call give_integrals_3_body_bi_ort(i ,p1,p2,i,h2,h1,integral)
+ int_exc_23= -1.d0 * integral
+
+ hthree += 1.d0 * int_direct + int_exc_l + int_exc_ll -( int_exc_12+ int_exc_13+ int_exc_23 )
+ enddo
+ do ii = Ne(2)+1,Ne(1) ! purely open-shell part
+ i = occ(ii,1)
+ call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
+ int_direct = -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p2,p1,i ,i,h2,h1,integral)
+ int_exc_l = -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p1,i ,p2,i,h2,h1,integral)
+ int_exc_ll= -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p2,i ,p1,i,h2,h1,integral)
+ int_exc_12= -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
+ int_exc_13= -1.d0 * integral
+ call give_integrals_3_body_bi_ort(i ,p1,p2,i,h2,h1,integral)
+ int_exc_23= -1.d0 * integral
+
+ hthree += 1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll -( int_exc_12+ int_exc_13+ int_exc_23 ))
+ enddo
+
+end subroutine give_aaa_contraction
+
+
+
+subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
+ implicit none
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+ integer, intent(in) :: Nint, h1, h2, p1, p2
+ integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
+ double precision, intent(out) :: hthree
+ integer :: ii, i
+ double precision :: int_direct, int_exc_12, int_exc_13, int_exc_23
+ double precision :: integral, int_exc_l, int_exc_ll
+
+ hthree = 0.d0
+ do ii = 1, Ne(2) ! purely closed shell part
+ i = occ(ii,2)
+ call give_integrals_3_body_bi_ort(p2,p1,i,h2,h1,i,integral)
+ int_direct = -1.d0 * integral
+ call give_integrals_3_body_bi_ort(p1,p2,i,h2,h1,i,integral)
+ int_exc_23= -1.d0 * integral
+ hthree += 1.d0 * int_direct - int_exc_23
+ enddo
+
+end subroutine give_aab_contraction
diff --git a/src/tc_bi_ortho/print_tc_wf.irp.f b/src/tc_bi_ortho/print_tc_wf.irp.f
new file mode 100644
index 00000000..58a733a7
--- /dev/null
+++ b/src/tc_bi_ortho/print_tc_wf.irp.f
@@ -0,0 +1,104 @@
+program print_tc_bi_ortho
+ implicit none
+ BEGIN_DOC
+! TODO : Put the documentation of the program here
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+! if(three_body_h_tc)then
+! call provide_all_three_ints_bi_ortho
+! endif
+! call routine
+ call write_l_r_wf
+end
+
+subroutine write_l_r_wf
+ implicit none
+ character*(128) :: output
+ integer :: i_unit_output,getUnitAndOpen
+ output=trim(ezfio_filename)//'.tc_wf'
+ i_unit_output = getUnitAndOpen(output,'w')
+ integer :: i
+ print*,'Writing the left-right wf'
+ do i = 1, N_det
+ write(i_unit_output,*)i,psi_l_coef_sorted_bi_ortho_left(i),psi_r_coef_sorted_bi_ortho_right(i)
+ enddo
+
+
+end
+
+subroutine routine
+ implicit none
+ integer :: i,degree
+ integer :: exc(0:2,2,2),h1,p1,s1,h2,p2,s2
+ double precision :: hmono,htwoe,hthree,htilde_ij,coef_pt1,e_i0,delta_e,e_pt2
+ double precision :: contrib_pt,e_corr,coef,contrib,phase
+ double precision :: accu_positive,accu_positive_pt, accu_positive_core,accu_positive_core_pt
+ e_pt2 = 0.d0
+ accu_positive = 0.D0
+ accu_positive_pt = 0.D0
+ accu_positive_core = 0.d0
+ accu_positive_core_pt = 0.d0
+
+ do i = 1, N_det
+ call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
+ if(degree == 1 .or. degree == 2)then
+ call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+ call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
+ delta_e = e_tilde_00 - e_i0
+ coef_pt1 = htilde_ij / delta_e
+
+ call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
+ contrib_pt = coef_pt1 * htilde_ij
+ e_pt2 += contrib_pt
+
+ coef = psi_r_coef_bi_ortho(i,1)/psi_r_coef_bi_ortho(1,1)
+ contrib = coef * htilde_ij
+ e_corr += contrib
+ call get_excitation(HF_bitmask,psi_det(1,1,i),exc,degree,phase,N_int)
+ call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+ print*,'*********'
+ if(degree==1)then
+ print*,'s1',s1
+ print*,'h1,p1 = ',h1,p1
+ else if(degree ==2)then
+ print*,'s1',s1
+ print*,'h1,p1 = ',h1,p1
+ print*,'s2',s2
+ print*,'h2,p2 = ',h2,p2
+ endif
+ print*,'coef_pt1 = ',coef_pt1
+ print*,'coef = ',coef
+ print*,'contrib_pt ',contrib_pt
+ print*,'contrib = ',contrib
+ if(contrib.gt.0.d0)then
+ accu_positive += contrib
+ if(h1==1.or.h2==1)then
+ accu_positive_core += contrib
+ endif
+ if(dabs(contrib).gt.1.d-5)then
+ print*,'Found a positive contribution to correlation energy !!'
+ endif
+ endif
+ if(contrib_pt.gt.0.d0)then
+ accu_positive_pt += contrib_pt
+ if(h2==1.or.h1==1)then
+ accu_positive_core_pt += contrib_pt
+ endif
+ endif
+ endif
+ enddo
+ print*,''
+ print*,''
+ print*,'Total correlation energy = ',e_corr
+ print*,'Total correlation energy PT = ',e_pt2
+ print*,'Positive contribution to ecorr = ',accu_positive
+ print*,'Positive contribution to ecorr PT = ',accu_positive_pt
+ print*,'Pure core contribution = ',accu_positive_core
+ print*,'Pure core contribution PT = ',accu_positive_core_pt
+end
diff --git a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
new file mode 100644
index 00000000..212c8588
--- /dev/null
+++ b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f
@@ -0,0 +1,157 @@
+use bitmasks
+
+BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_tc, (psi_det_size) ]
+ implicit none
+ BEGIN_DOC
+ ! Contribution of determinants to the state-averaged density.
+ END_DOC
+ integer :: i,j,k
+ double precision :: f
+
+ psi_average_norm_contrib_tc(:) = 0.d0
+ do k=1,N_states
+ do i=1,N_det
+ psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i) + &
+ dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k))*state_average_weight(k)
+ enddo
+ enddo
+ f = 1.d0/sum(psi_average_norm_contrib_tc(1:N_det))
+ do i=1,N_det
+ psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i)*f
+ enddo
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc, (psi_det_size,N_states) ]
+&BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted_tc, (psi_det_size) ]
+&BEGIN_PROVIDER [ integer, psi_det_sorted_tc_order, (psi_det_size) ]
+ implicit none
+ BEGIN_DOC
+ ! Wave function sorted by determinants contribution to the norm (state-averaged)
+ !
+ ! psi_det_sorted_tc_order(i) -> k : index in psi_det
+ END_DOC
+ integer :: i,j,k
+ integer, allocatable :: iorder(:)
+ allocate ( iorder(N_det) )
+ do i=1,N_det
+ psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_tc(i)
+ iorder(i) = i
+ enddo
+ call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det)
+ do i=1,N_det
+ do j=1,N_int
+ psi_det_sorted_tc(j,1,i) = psi_det(j,1,iorder(i))
+ psi_det_sorted_tc(j,2,i) = psi_det(j,2,iorder(i))
+ enddo
+ psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_sorted_tc(i)
+ psi_det_sorted_tc_order(iorder(i)) = i
+ enddo
+ double precision :: accu
+ do k=1,N_states
+ accu = 0.d0
+ do i=1,N_det
+ psi_coef_sorted_tc(i,k) = dsqrt(dabs(psi_l_coef_bi_ortho(iorder(i),k)*psi_r_coef_bi_ortho(iorder(i),k)))
+ accu += psi_coef_sorted_tc(i,k)**2
+ enddo
+ accu = 1.d0/dsqrt(accu)
+ do i=1,N_det
+ psi_coef_sorted_tc(i,k) *= accu
+ enddo
+ enddo
+
+ psi_det_sorted_tc(:,:,N_det+1:psi_det_size) = 0_bit_kind
+ psi_coef_sorted_tc(N_det+1:psi_det_size,:) = 0.d0
+ psi_average_norm_contrib_sorted_tc(N_det+1:psi_det_size) = 0.d0
+ psi_det_sorted_tc_order(N_det+1:psi_det_size) = 0
+
+ deallocate(iorder)
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho, (psi_det_size, N_states)]
+&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho, (psi_det_size, N_states)]
+ BEGIN_DOC
+ ! psi_r_coef_sorted_bi_ortho : right coefficients corresponding to psi_det_sorted_tc
+ ! psi_l_coef_sorted_bi_ortho : left coefficients corresponding to psi_det_sorted_tc
+ END_DOC
+ implicit none
+ integer :: i, j, k
+ psi_r_coef_sorted_bi_ortho = 0.d0
+ psi_l_coef_sorted_bi_ortho = 0.d0
+ do i = 1, N_det
+ psi_r_coef_sorted_bi_ortho(i,1) = psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
+ psi_l_coef_sorted_bi_ortho(i,1) = psi_l_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
+ enddo
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_bit, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc_bit, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Determinants on which we apply $\langle i|H|psi \rangle$ for perturbation.
+ ! They are sorted by determinants interpreted as integers. Useful
+ ! to accelerate the search of a random determinant in the wave
+ ! function.
+ END_DOC
+
+ call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, size(psi_coef,1), &
+ psi_det_sorted_tc_bit, psi_coef_sorted_tc_bit, N_states)
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_right, (N_int,2,N_det) ]
+&BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho_right, (N_det)]
+ implicit none
+ BEGIN_DOC
+ ! psi_det_sorted_tc_right : Slater determinants sorted by decreasing value of |right- coefficients|
+ !
+ ! psi_r_coef_sorted_bi_ortho_right : right wave function according to psi_det_sorted_tc_right
+ END_DOC
+ integer, allocatable :: iorder(:)
+ double precision, allocatable :: coef(:)
+ integer :: i,j
+ allocate ( iorder(N_det) , coef(N_det))
+ do i=1,N_det
+ coef(i) = -dabs(psi_r_coef_bi_ortho(i,1)/psi_r_coef_bi_ortho(1,1))
+ iorder(i) = i
+ enddo
+ call dsort(coef,iorder,N_det)
+ do i=1,N_det
+ do j=1,N_int
+ psi_det_sorted_tc_right(j,1,i) = psi_det(j,1,iorder(i))
+ psi_det_sorted_tc_right(j,2,i) = psi_det(j,2,iorder(i))
+ enddo
+ psi_r_coef_sorted_bi_ortho_right(i) = psi_r_coef_bi_ortho(iorder(i),1)/psi_r_coef_bi_ortho(iorder(1),1)
+ enddo
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_left, (N_int,2,N_det) ]
+&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho_left, (N_det)]
+ implicit none
+ BEGIN_DOC
+ ! psi_det_sorted_tc_left : Slater determinants sorted by decreasing value of |LEFTt- coefficients|
+ !
+ ! psi_r_coef_sorted_bi_ortho_left : LEFT wave function according to psi_det_sorted_tc_left
+ END_DOC
+ integer, allocatable :: iorder(:)
+ double precision, allocatable :: coef(:)
+ integer :: i,j
+ allocate ( iorder(N_det) , coef(N_det))
+ do i=1,N_det
+ coef(i) = -dabs(psi_l_coef_bi_ortho(i,1)/psi_r_coef_bi_ortho(1,1))
+ iorder(i) = i
+ enddo
+ call dsort(coef,iorder,N_det)
+ do i=1,N_det
+ do j=1,N_int
+ psi_det_sorted_tc_left(j,1,i) = psi_det(j,1,iorder(i))
+ psi_det_sorted_tc_left(j,2,i) = psi_det(j,2,iorder(i))
+ enddo
+ psi_l_coef_sorted_bi_ortho_left(i) = psi_l_coef_bi_ortho(iorder(i),1)/psi_l_coef_bi_ortho(iorder(1),1)
+ enddo
+END_PROVIDER
diff --git a/src/tc_bi_ortho/psi_left_qmc.irp.f b/src/tc_bi_ortho/psi_left_qmc.irp.f
new file mode 100644
index 00000000..25048f82
--- /dev/null
+++ b/src/tc_bi_ortho/psi_left_qmc.irp.f
@@ -0,0 +1,44 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, psi_bitcleft_bilinear_matrix_values, (N_det,N_states) ]
+
+ BEGIN_DOC
+ ! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
+ ! $D_\alpha^\dagger.C.D_\beta$
+ !
+ ! Rows are $\alpha$ determinants and columns are $\beta$.
+ !
+ ! Order refers to psi_det
+ END_DOC
+
+ use bitmasks
+
+ implicit none
+ integer :: k, l
+
+ if(N_det .eq. 1) then
+
+ do l = 1, N_states
+ psi_bitcleft_bilinear_matrix_values(1,l) = 1.d0
+ enddo
+
+ else
+
+ do l = 1, N_states
+ do k = 1, N_det
+ psi_bitcleft_bilinear_matrix_values(k,l) = psi_l_coef_bi_ortho(k,l)
+ enddo
+ enddo
+
+ PROVIDE psi_bilinear_matrix_order
+ do l = 1, N_states
+ call dset_order(psi_bitcleft_bilinear_matrix_values(1,l), psi_bilinear_matrix_order, N_det)
+ enddo
+
+ endif
+
+END_PROVIDER
+
+! ---
+
diff --git a/src/tc_bi_ortho/psi_r_l_prov.irp.f b/src/tc_bi_ortho/psi_r_l_prov.irp.f
new file mode 100644
index 00000000..551da858
--- /dev/null
+++ b/src/tc_bi_ortho/psi_r_l_prov.irp.f
@@ -0,0 +1,217 @@
+use bitmasks
+
+BEGIN_PROVIDER [ double precision, psi_l_coef_bi_ortho, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file
+ ! is empty.
+ END_DOC
+
+ integer :: i,k, N_int2
+ logical :: exists
+ character*(64) :: label
+
+ PROVIDE read_wf N_det mo_label ezfio_filename nproc
+ psi_l_coef_bi_ortho = 0.d0
+ do i=1,min(N_states,N_det)
+ psi_l_coef_bi_ortho(i,i) = 1.d0
+ enddo
+
+ if (mpi_master) then
+ if (read_wf) then
+ call ezfio_has_tc_bi_ortho_psi_l_coef_bi_ortho(exists)
+! if (exists) then
+! call ezfio_has_tc_bi_ortho_mo_label(exists)
+! if (exists) then
+! call ezfio_get_tc_bi_ortho_mo_label(label)
+! exists = (label == mo_label)
+! endif
+! endif
+
+ if (exists) then
+
+ double precision, allocatable :: psi_l_coef_bi_ortho_read(:,:)
+ allocate (psi_l_coef_bi_ortho_read(N_det,N_states))
+ print *, 'Read psi_l_coef_bi_ortho', N_det, N_states
+ call ezfio_get_tc_bi_ortho_psi_l_coef_bi_ortho(psi_l_coef_bi_ortho_read)
+ do k=1,N_states
+ do i=1,N_det
+ psi_l_coef_bi_ortho(i,k) = psi_l_coef_bi_ortho_read(i,k)
+ enddo
+ enddo
+ deallocate(psi_l_coef_bi_ortho_read)
+
+ endif
+ endif
+ endif
+ IRP_IF MPI_DEBUG
+ print *, irp_here, mpi_rank
+ call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+ IRP_ENDIF
+ IRP_IF MPI
+ include 'mpif.h'
+ integer :: ierr
+ call MPI_BCAST( psi_l_coef_bi_ortho, size(psi_l_coef_bi_ortho), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+ if (ierr /= MPI_SUCCESS) then
+ stop 'Unable to read psi_l_coef_bi_ortho with MPI'
+ endif
+ IRP_ENDIF
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file
+ ! is empty.
+ END_DOC
+
+ integer :: i,k, N_int2
+ logical :: exists
+ character*(64) :: label
+
+ PROVIDE read_wf N_det mo_label ezfio_filename nproc
+ psi_r_coef_bi_ortho = 0.d0
+ do i=1,min(N_states,N_det)
+ psi_r_coef_bi_ortho(i,i) = 1.d0
+ enddo
+
+ if (mpi_master) then
+ if (read_wf) then
+ call ezfio_has_tc_bi_ortho_psi_r_coef_bi_ortho(exists)
+! if (exists) then
+! call ezfio_has_tc_bi_ortho_mo_label(exists)
+! if (exists) then
+! call ezfio_get_tc_bi_ortho_mo_label(label)
+! exists = (label == mo_label)
+! endif
+! endif
+
+ if (exists) then
+
+ double precision, allocatable :: psi_r_coef_bi_ortho_read(:,:)
+ allocate (psi_r_coef_bi_ortho_read(N_det,N_states))
+ print *, 'Read psi_r_coef_bi_ortho', N_det, N_states
+ call ezfio_get_tc_bi_ortho_psi_r_coef_bi_ortho(psi_r_coef_bi_ortho_read)
+ do k=1,N_states
+ do i=1,N_det
+ psi_r_coef_bi_ortho(i,k) = psi_r_coef_bi_ortho_read(i,k)
+ enddo
+ enddo
+ deallocate(psi_r_coef_bi_ortho_read)
+
+ endif
+ endif
+ endif
+ IRP_IF MPI_DEBUG
+ print *, irp_here, mpi_rank
+ call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+ IRP_ENDIF
+ IRP_IF MPI
+ include 'mpif.h'
+ integer :: ierr
+ call MPI_BCAST( psi_r_coef_bi_ortho, size(psi_r_coef_bi_ortho), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+ if (ierr /= MPI_SUCCESS) then
+ stop 'Unable to read psi_r_coef_bi_ortho with MPI'
+ endif
+ IRP_ENDIF
+END_PROVIDER
+
+
+subroutine save_tc_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psilcoef,psircoef)
+ implicit none
+ BEGIN_DOC
+ ! Save the wave function into the |EZFIO| file
+ END_DOC
+ use bitmasks
+ include 'constants.include.F'
+ integer, intent(in) :: ndet,nstates,dim_psicoef
+ integer(bit_kind), intent(in) :: psidet(N_int,2,ndet)
+ double precision, intent(in) :: psilcoef(dim_psicoef,nstates)
+ double precision, intent(in) :: psircoef(dim_psicoef,nstates)
+ integer*8, allocatable :: psi_det_save(:,:,:)
+ double precision, allocatable :: psil_coef_save(:,:)
+ double precision, allocatable :: psir_coef_save(:,:)
+
+ double precision :: accu_norm
+ integer :: i,j,k, ndet_qp_edit
+
+ if (mpi_master) then
+ ndet_qp_edit = min(ndet,N_det_qp_edit)
+
+ call ezfio_set_determinants_N_int(N_int)
+ call ezfio_set_determinants_bit_kind(bit_kind)
+ call ezfio_set_determinants_N_det(ndet)
+ call ezfio_set_determinants_N_det_qp_edit(ndet_qp_edit)
+ call ezfio_set_determinants_n_states(nstates)
+ call ezfio_set_determinants_mo_label(mo_label)
+
+ allocate (psi_det_save(N_int,2,ndet))
+ do i=1,ndet
+ do j=1,2
+ do k=1,N_int
+ psi_det_save(k,j,i) = transfer(psidet(k,j,i),1_8)
+ enddo
+ enddo
+ enddo
+ call ezfio_set_determinants_psi_det(psi_det_save)
+ call ezfio_set_determinants_psi_det_qp_edit(psi_det_save)
+ deallocate (psi_det_save)
+
+ allocate (psil_coef_save(ndet,nstates),psir_coef_save(ndet,nstates))
+ do k=1,nstates
+ do i=1,ndet
+ psil_coef_save(i,k) = psilcoef(i,k)
+ psir_coef_save(i,k) = psircoef(i,k)
+ enddo
+ enddo
+
+ call ezfio_set_tc_bi_ortho_psi_l_coef_bi_ortho(psil_coef_save)
+ call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(psir_coef_save)
+ deallocate (psil_coef_save,psir_coef_save)
+
+! allocate (psi_coef_save(ndet_qp_edit,nstates))
+! do k=1,nstates
+! do i=1,ndet_qp_edit
+! psi_coef_save(i,k) = psicoef(i,k)
+! enddo
+! enddo
+!
+! call ezfio_set_determinants_psi_coef_qp_edit(psi_coef_save)
+! deallocate (psi_coef_save)
+
+ call write_int(6,ndet,'Saved determinantsi and psi_r/psi_l coef')
+ endif
+end
+
+subroutine save_tc_bi_ortho_wavefunction
+ implicit none
+ call save_tc_wavefunction_general(N_det,N_states,psi_det,size(psi_l_coef_bi_ortho, 1),psi_l_coef_bi_ortho,psi_r_coef_bi_ortho)
+ call routine_save_right_bi_ortho
+end
+
+subroutine routine_save_right_bi_ortho
+ implicit none
+ double precision, allocatable :: coef_tmp(:,:)
+ integer :: i
+ N_states = 1
+ allocate(coef_tmp(N_det, N_states))
+ do i = 1, N_det
+ coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1)
+ enddo
+ call save_wavefunction_general_unormalized(N_det,N_states,psi_det,size(coef_tmp,1),coef_tmp(1,1))
+end
+
+subroutine routine_save_left_right_bi_ortho
+ implicit none
+ double precision, allocatable :: coef_tmp(:,:)
+ integer :: i,n_states_tmp
+ n_states_tmp = 2
+ allocate(coef_tmp(N_det, n_states_tmp))
+ do i = 1, N_det
+ coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1)
+ coef_tmp(i,2) = psi_l_coef_bi_ortho(i,1)
+ enddo
+ call save_wavefunction_general_unormalized(N_det,n_states_tmp,psi_det,size(coef_tmp,1),coef_tmp(1,1))
+end
+
diff --git a/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f b/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f
new file mode 100644
index 00000000..60201f5f
--- /dev/null
+++ b/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f
@@ -0,0 +1,61 @@
+program save_bitcpsileft_for_qmcchem
+
+ integer :: iunit
+ logical :: exists
+ double precision :: e_ref
+
+ print *, ' '
+ print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+ print *, ' call save_for_qmcchem before '
+ print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+ print *, ' '
+
+ call write_lr_spindeterminants()
+
+ e_ref = 0.d0
+ iunit = 13
+ open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write')
+ call ezfio_has_fci_energy_pt2(exists)
+
+ if(.not.exists) then
+ call ezfio_has_fci_energy(exists)
+
+ if(.not.exists) then
+ call ezfio_has_tc_scf_bitc_energy(exists)
+ if(exists) then
+ call ezfio_get_tc_scf_bitc_energy(e_ref)
+ endif
+ endif
+
+ endif
+ write(iunit,*) e_ref
+ close(iunit)
+
+end
+
+! --
+
+subroutine write_lr_spindeterminants()
+
+ use bitmasks
+
+ implicit none
+
+ integer :: k, l
+ double precision, allocatable :: buffer(:,:)
+
+ PROVIDE psi_bitcleft_bilinear_matrix_values
+
+ allocate(buffer(N_det,N_states))
+ do l = 1, N_states
+ do k = 1, N_det
+ buffer(k,l) = psi_bitcleft_bilinear_matrix_values(k,l)
+ enddo
+ enddo
+ call ezfio_set_spindeterminants_psi_left_coef_matrix_values(buffer)
+ deallocate(buffer)
+
+end subroutine write_lr_spindeterminants
+
+! ---
+
diff --git a/src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f b/src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f
new file mode 100644
index 00000000..5eb3c069
--- /dev/null
+++ b/src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f
@@ -0,0 +1,15 @@
+program tc_bi_ortho
+ implicit none
+ BEGIN_DOC
+! TODO : Put the documentation of the program here
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+ call routine_save_left_right_bi_ortho
+! call test
+end
diff --git a/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f b/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f
new file mode 100644
index 00000000..9b7b9f5a
--- /dev/null
+++ b/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f
@@ -0,0 +1,33 @@
+ program tc_natorb_bi_ortho
+ implicit none
+ BEGIN_DOC
+ ! TODO : Put the documentation of the program here
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+ call save_tc_natorb
+ end
+
+ subroutine save_tc_natorb
+ implicit none
+ print*,'Saving the natorbs '
+ provide natorb_tc_leigvec_ao natorb_tc_reigvec_ao
+ call ezfio_set_bi_ortho_mos_mo_l_coef(natorb_tc_leigvec_ao)
+ call ezfio_set_bi_ortho_mos_mo_r_coef(natorb_tc_reigvec_ao)
+ call save_ref_determinant_nstates_1
+ call ezfio_set_determinants_read_wf(.False.)
+ end
+
+ subroutine save_ref_determinant_nstates_1
+ implicit none
+ use bitmasks
+ double precision :: buffer(1,N_states)
+ buffer = 0.d0
+ buffer(1,1) = 1.d0
+ call save_wavefunction_general(1,1,ref_bitmask,1,buffer)
+ end
diff --git a/src/tc_bi_ortho/select_dets_bi_ortho.irp.f b/src/tc_bi_ortho/select_dets_bi_ortho.irp.f
new file mode 100644
index 00000000..e6bf3d6e
--- /dev/null
+++ b/src/tc_bi_ortho/select_dets_bi_ortho.irp.f
@@ -0,0 +1,61 @@
+program tc_bi_ortho
+ implicit none
+ BEGIN_DOC
+! TODO : Put the documentation of the program here
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+ !!!!!!!!!!!!!!! WARNING NO 3-BODY
+ !!!!!!!!!!!!!!! WARNING NO 3-BODY
+ three_body_h_tc = .False.
+ touch three_body_h_tc
+ !!!!!!!!!!!!!!! WARNING NO 3-BODY
+ !!!!!!!!!!!!!!! WARNING NO 3-BODY
+
+ call routine_test
+! call test
+end
+
+subroutine routine_test
+ implicit none
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+ integer :: i,n_good,degree
+ integer(bit_kind), allocatable :: dets(:,:,:)
+ integer, allocatable :: iorder(:)
+ double precision, allocatable :: coef(:),coef_new(:,:)
+ double precision :: thr
+ allocate(coef(N_det), iorder(N_det))
+ do i = 1, N_det
+ iorder(i) = i
+ call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
+ if(degree==1)then
+ coef(i) = -0.5d0
+ else
+ coef(i) = -dabs(coef_pt1_bi_ortho(i))
+ endif
+ enddo
+ call dsort(coef,iorder,N_det)
+ !thr = save_threshold
+ thr = 1d-15
+ n_good = 0
+ do i = 1, N_det
+ if(dabs(coef(i)).gt.thr)then
+ n_good += 1
+ endif
+ enddo
+ print*,'n_good = ',n_good
+ allocate(dets(N_int,2,n_good),coef_new(n_good,n_states))
+ do i = 1, n_good
+ dets(:,:,i) = psi_det(:,:,iorder(i))
+ coef_new(i,:) = psi_coef(iorder(i),:)
+ enddo
+ call save_wavefunction_general(n_good,n_states,dets,n_good,coef_new)
+
+
+end
diff --git a/src/tc_bi_ortho/slater_tc.irp.f b/src/tc_bi_ortho/slater_tc.irp.f
new file mode 100644
index 00000000..45115a40
--- /dev/null
+++ b/src/tc_bi_ortho/slater_tc.irp.f
@@ -0,0 +1,359 @@
+!!!!!!
+subroutine htilde_mu_mat_bi_ortho_tot(key_j, key_i, Nint, htot)
+
+ BEGIN_DOC
+ ! where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
+ !!
+ !! WARNING !!
+ !
+ ! Non hermitian !!
+ END_DOC
+
+ use bitmasks
+
+ implicit none
+ integer, intent(in) :: Nint
+ integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
+ double precision, intent(out) :: htot
+ double precision :: hmono,htwoe,hthree
+ integer :: degree
+ call get_excitation_degree(key_j, key_i, degree, Nint)
+ if(degree.gt.2)then
+ htot = 0.d0
+ else
+ call htilde_mu_mat_bi_ortho(key_j,key_i, Nint, hmono,htwoe,hthree,htot)
+ endif
+
+end subroutine htilde_mu_mat_tot
+
+! --
+
+subroutine htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
+ implicit none
+ use bitmasks
+ BEGIN_DOC
+ ! where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
+ !!
+ ! Returns the detail of the matrix element in terms of single, two and three electron contribution.
+ !! WARNING !!
+ !
+ ! Non hermitian !!
+ END_DOC
+ integer, intent(in) :: Nint
+ integer(bit_kind), intent(in) :: key_i(Nint,2),key_j(Nint,2)
+ double precision, intent(out) :: hmono,htwoe,hthree,htot
+ integer :: degree
+
+ hmono = 0.d0
+ htwoe= 0.d0
+ htot = 0.d0
+ hthree = 0.D0
+ call get_excitation_degree(key_i, key_j, degree, Nint)
+ if(degree.gt.2)return
+ if(degree == 0)then
+ call diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot)
+ else if (degree == 1)then
+ call single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
+ else if(degree == 2)then
+ call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
+ endif
+ if(three_body_h_tc) then
+ if(degree == 2) then
+ if(.not.double_normal_ord) then
+ call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
+ endif
+ else if(degree == 1)then
+ call single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
+ else if(degree == 0)then
+ call diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
+ endif
+ endif
+ htot = hmono + htwoe + hthree
+ if(degree==0)then
+ htot += nuclear_repulsion
+ endif
+
+end
+
+subroutine diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot)
+
+ BEGIN_DOC
+ ! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS
+ END_DOC
+
+ use bitmasks
+
+ implicit none
+ integer, intent(in) :: Nint
+ integer(bit_kind), intent(in) :: key_i(Nint,2)
+ double precision, intent(out) :: hmono,htwoe,htot
+ integer :: occ(Nint*bit_kind_size,2)
+ integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
+ double precision :: get_mo_two_e_integral_tc_int
+ integer(bit_kind) :: key_i_core(Nint,2)
+
+! 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
+
+ do ispin = 1, 2
+ do i = 1, Ne(ispin) !
+ ii = occ(i,ispin)
+ hmono += mo_bi_ortho_tc_one_e(ii,ii)
+
+! if(j1b_gauss .eq. 1) then
+! print*,'j1b not implemented for bi ortho TC'
+! print*,'stopping ....'
+! stop
+! !hmono += mo_j1b_gauss_hermI (ii,ii) &
+! ! + mo_j1b_gauss_hermII (ii,ii) &
+! ! + mo_j1b_gauss_nonherm(ii,ii)
+! endif
+
+! 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
+ do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1))
+ ii = occ(i,ispin)
+ do j = 1, Ne(jspin) ! electron 2
+ jj = occ(j,jspin)
+ htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii)
+ enddo
+ enddo
+
+ ! alpha/alpha two-body
+ do i = 1, Ne(ispin)
+ ii = occ(i,ispin)
+ do j = i+1, Ne(ispin)
+ jj = occ(j,ispin)
+ htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
+ enddo
+ enddo
+
+ ! beta/beta two-body
+ do i = 1, Ne(jspin)
+ ii = occ(i,jspin)
+ do j = i+1, Ne(jspin)
+ jj = occ(j,jspin)
+ 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(Nint, key_j, key_i, hmono, htwoe, htot)
+
+ BEGIN_DOC
+ ! for double excitation ONLY FOR ONE- AND TWO-BODY TERMS
+ !!
+ !! WARNING !!
+ !
+ ! Non hermitian !!
+ END_DOC
+
+ use bitmasks
+
+ implicit none
+ integer, intent(in) :: Nint
+ integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
+ double precision, intent(out) :: hmono, htwoe, htot
+ integer :: occ(Nint*bit_kind_size,2)
+ integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
+ integer :: degree,exc(0:2,2,2)
+ integer :: h1, p1, h2, p2, s1, s2
+ integer :: other_spin(2)
+ integer(bit_kind) :: key_i_core(Nint,2)
+ double precision :: get_mo_two_e_integral_tc_int,phase
+
+! PROVIDE mo_two_e_integrals_tc_int_in_map mo_bi_ortho_tc_two_e
+
+ other_spin(1) = 2
+ other_spin(2) = 1
+
+ call get_excitation_degree(key_i, key_j, degree, Nint)
+
+ hmono = 0.d0
+ htwoe= 0.d0
+ htot = 0.d0
+
+ if(degree.ne.2)then
+ return
+ endif
+
+! 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)
+! else
+ call bitstring_to_list_ab(key_i, occ, Ne, Nint)
+! endif
+ call get_double_excitation(key_i, key_j, exc, phase, Nint)
+ call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
+
+ if(s1.ne.s2)then
+ ! opposite spin two-body
+! key_j, key_i
+ htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
+ if(double_normal_ord.and.+Ne(1).gt.2)then
+ htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ???
+ 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(double_normal_ord.and.+Ne(1).gt.2)then
+ htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ???
+ htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ???
+ endif
+ endif
+ htwoe *= phase
+ htot = htwoe
+
+end
+
+
+subroutine single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
+
+ BEGIN_DOC
+ ! for single excitation ONLY FOR ONE- AND TWO-BODY TERMS
+ !!
+ !! WARNING !!
+ !
+ ! Non hermitian !!
+ END_DOC
+
+ use bitmasks
+
+ implicit none
+ integer, intent(in) :: Nint
+ integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
+ double precision, intent(out) :: hmono, htwoe, htot
+ integer :: occ(Nint*bit_kind_size,2)
+ integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
+ integer :: degree,exc(0:2,2,2)
+ integer :: h1, p1, h2, p2, s1, s2
+ double precision :: get_mo_two_e_integral_tc_int, phase
+ double precision :: direct_int, exchange_int_12, exchange_int_23, exchange_int_13
+ integer :: other_spin(2)
+ integer(bit_kind) :: key_j_core(Nint,2), key_i_core(Nint,2)
+
+! PROVIDE mo_two_e_integrals_tc_int_in_map mo_bi_ortho_tc_two_e
+!
+! PROVIDE core_bitmask core_fock_operator mo_integrals_erf_map
+
+! PROVIDE j1b_gauss
+
+ other_spin(1) = 2
+ other_spin(2) = 1
+
+ hmono = 0.d0
+ htwoe= 0.d0
+ htot = 0.d0
+ call get_excitation_degree(key_i, key_j, degree, Nint)
+ if(degree.ne.1)then
+ return
+ endif
+! 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))
+! key_j_core(i,1) = xor(key_j(i,1),core_bitmask(i,1))
+! key_j_core(i,2) = xor(key_j(i,2),core_bitmask(i,2))
+! enddo
+! call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
+! else
+ call bitstring_to_list_ab(key_i, occ, Ne, Nint)
+! endif
+
+ call get_single_excitation(key_i, key_j, exc, phase, Nint)
+ call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
+
+ hmono = mo_bi_ortho_tc_one_e(p1,h1) * phase
+
+! if(j1b_gauss .eq. 1) then
+! print*,'j1b not implemented for bi ortho TC'
+! print*,'stopping ....'
+! stop
+! !hmono += ( mo_j1b_gauss_hermI (h1,p1) &
+! ! + mo_j1b_gauss_hermII (h1,p1) &
+! ! + mo_j1b_gauss_nonherm(h1,p1) ) * phase
+! endif
+
+! if(core_tc_op)then
+! print*,'core_tc_op not already taken into account for bi ortho'
+! print*,'stopping ...'
+! stop
+! hmono += phase * core_fock_operator(h1,p1)
+! endif
+
+ ! alpha/beta two-body
+ ispin = other_spin(s1)
+ if(s1==1)then
+ ! single alpha
+ do i = 1, Ne(ispin) ! electron 2
+ ii = occ(i,ispin)
+ htwoe += mo_bi_ortho_tc_two_e(ii,p1,ii,h1)
+ enddo
+ else
+ ! single beta
+ do i = 1, Ne(ispin) ! electron 1
+ ii = occ(i,ispin)
+ htwoe += mo_bi_ortho_tc_two_e(p1,ii,h1,ii)
+ enddo
+ endif
+! ! same spin two-body
+ do i = 1, Ne(s1)
+ ii = occ(i,s1)
+ ! (h1p1|ii ii) - (h1 ii| p1 ii)
+ htwoe += mo_bi_ortho_tc_two_e(ii,p1,ii,h1) - mo_bi_ortho_tc_two_e(p1,ii,ii,h1)
+ enddo
+
+ htwoe *= phase
+ htot = hmono + htwoe
+
+end
+
+
diff --git a/src/tc_bi_ortho/slater_tc_3e.irp.f b/src/tc_bi_ortho/slater_tc_3e.irp.f
new file mode 100644
index 00000000..7c2c9c86
--- /dev/null
+++ b/src/tc_bi_ortho/slater_tc_3e.irp.f
@@ -0,0 +1,293 @@
+subroutine provide_all_three_ints_bi_ortho
+ implicit none
+ BEGIN_DOC
+! routine that provides all necessary three-electron integrals
+ END_DOC
+ if(three_body_h_tc)then
+ PROVIDE three_e_3_idx_direct_bi_ort three_e_3_idx_cycle_1_bi_ort three_e_3_idx_cycle_2_bi_ort
+ PROVIDE three_e_3_idx_exch23_bi_ort three_e_3_idx_exch13_bi_ort three_e_3_idx_exch12_bi_ort
+ PROVIDE three_e_4_idx_direct_bi_ort three_e_4_idx_cycle_1_bi_ort three_e_4_idx_cycle_2_bi_ort
+ PROVIDE three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort three_e_4_idx_exch12_bi_ort
+ endif
+if(.not.double_normal_ord)then
+ PROVIDE three_e_5_idx_direct_bi_ort three_e_5_idx_cycle_1_bi_ort three_e_5_idx_cycle_2_bi_ort
+ PROVIDE three_e_5_idx_exch23_bi_ort three_e_5_idx_exch13_bi_ort three_e_5_idx_exch12_bi_ort
+else
+ PROVIDE normal_two_body_bi_orth
+endif
+end
+
+subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
+
+ BEGIN_DOC
+ ! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
+ END_DOC
+
+ use bitmasks
+
+ implicit none
+ integer, intent(in) :: Nint
+ integer(bit_kind), intent(in) :: key_i(Nint,2)
+ double precision, intent(out) :: hthree
+ integer :: occ(Nint*bit_kind_size,2)
+ integer :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm
+ integer(bit_kind) :: key_i_core(Nint,2)
+ double precision :: direct_int, exchange_int
+ double precision :: sym_3_e_int_from_6_idx_tensor
+ double precision :: three_e_diag_parrallel_spin
+
+ if(core_tc_op)then
+ 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)
+ else
+ call bitstring_to_list_ab(key_i,occ,Ne,Nint)
+ endif
+ hthree = 0.d0
+
+ if(Ne(1)+Ne(2).ge.3)then
+!! ! alpha/alpha/beta three-body
+ do i = 1, Ne(1)
+ ii = occ(i,1)
+ do j = i+1, Ne(1)
+ jj = occ(j,1)
+ do m = 1, Ne(2)
+ mm = occ(m,2)
+! direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) USES THE 6-IDX TENSOR
+! exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) USES THE 6-IDX TENSOR
+ direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR
+ exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR
+ hthree += direct_int - exchange_int
+ enddo
+ enddo
+ enddo
+
+ ! beta/beta/alpha three-body
+ do i = 1, Ne(2)
+ ii = occ(i,2)
+ do j = i+1, Ne(2)
+ jj = occ(j,2)
+ do m = 1, Ne(1)
+ mm = occ(m,1)
+ direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii)
+ exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii)
+ hthree += direct_int - exchange_int
+ enddo
+ enddo
+ enddo
+
+ ! alpha/alpha/alpha three-body
+ do i = 1, Ne(1)
+ ii = occ(i,1) ! 1
+ do j = i+1, Ne(1)
+ jj = occ(j,1) ! 2
+ do m = j+1, Ne(1)
+ mm = occ(m,1) ! 3
+! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR
+ hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS
+ enddo
+ enddo
+ enddo
+
+ ! beta/beta/beta three-body
+ do i = 1, Ne(2)
+ ii = occ(i,2) ! 1
+ do j = i+1, Ne(2)
+ jj = occ(j,2) ! 2
+ do m = j+1, Ne(2)
+ mm = occ(m,2) ! 3
+! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR
+ hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS
+ enddo
+ enddo
+ enddo
+ endif
+
+end
+
+
+subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
+
+ BEGIN_DOC
+ ! for single excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
+ !!
+ !! WARNING !!
+ !
+ ! Non hermitian !!
+ END_DOC
+
+ use bitmasks
+
+ implicit none
+ integer, intent(in) :: Nint
+ integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
+ double precision, intent(out) :: hthree
+ integer :: occ(Nint*bit_kind_size,2)
+ integer :: Ne(2),i,j,ii,jj,ispin,jspin,k,kk
+ integer :: degree,exc(0:2,2,2)
+ integer :: h1, p1, h2, p2, s1, s2
+ double precision :: direct_int,phase,exchange_int,three_e_single_parrallel_spin
+ double precision :: sym_3_e_int_from_6_idx_tensor
+ integer :: other_spin(2)
+ integer(bit_kind) :: key_j_core(Nint,2),key_i_core(Nint,2)
+
+ other_spin(1) = 2
+ other_spin(2) = 1
+
+
+ hthree = 0.d0
+ call get_excitation_degree(key_i,key_j,degree,Nint)
+ if(degree.ne.1)then
+ return
+ endif
+ if(core_tc_op)then
+ 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))
+ key_j_core(i,1) = xor(key_j(i,1),core_bitmask(i,1))
+ key_j_core(i,2) = xor(key_j(i,2),core_bitmask(i,2))
+ enddo
+ call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
+ else
+ call bitstring_to_list_ab(key_i, occ, Ne, Nint)
+ endif
+
+ call get_single_excitation(key_i, key_j, exc, phase, Nint)
+ call decode_exc(exc, 1, h1, p1, h2, p2, s1, s2)
+
+ ! alpha/alpha/beta three-body
+! print*,'IN SLAT RULES'
+ if(Ne(1)+Ne(2).ge.3)then
+ ! hole of spin s1 :: contribution from purely other spin
+ ispin = other_spin(s1) ! ispin is the other spin than s1
+ do i = 1, Ne(ispin) ! i is the orbitals of the other spin than s1
+ ii = occ(i,ispin)
+ do j = i+1, Ne(ispin) ! j has the same spin than s1
+ jj = occ(j,ispin)
+ ! is == ispin in ::: s1 is is s1 is is s1 is is s1 is is
+ ! < h1 j i | p1 j i > - < h1 j i | p1 i j >
+ !
+! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor
+! exchange_int = three_body_ints_bi_ort(jj,ii,p1,ii,jj,h1) ! USES THE 6-IDX tensor
+ direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1)
+ exchange_int = three_e_4_idx_exch23_bi_ort(jj,ii,p1,h1)
+ hthree += direct_int - exchange_int
+ enddo
+ enddo
+
+ ! hole of spin s1 :: contribution from mixed other spin / same spin
+ do i = 1, Ne(ispin) ! other spin
+ ii = occ(i,ispin) ! other spin
+ do j = 1, Ne(s1) ! same spin
+ jj = occ(j,s1) ! same spin
+! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor
+! exchange_int = three_body_ints_bi_ort(jj,ii,p1,h1,ii,jj) ! exchange the spin s1 :: 6-IDX tensor
+ direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1)
+ exchange_int = three_e_4_idx_exch13_bi_ort(jj,ii,p1,h1)
+ ! < h1 j i | p1 j i > - < h1 j i | j p1 i >
+ hthree += direct_int - exchange_int
+! print*,'h1,p1,ii,jj = ',h1,p1,ii,jj
+! print*,direct_int, exchange_int
+ enddo
+ enddo
+!
+ ! hole of spin s1 :: PURE SAME SPIN CONTRIBUTIONS !!!
+ do i = 1, Ne(s1)
+ ii = occ(i,s1)
+ do j = i+1, Ne(s1)
+ jj = occ(j,s1)
+! ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1)
+ hthree += three_e_single_parrallel_spin(jj,ii,p1,h1) ! USES THE 4-IDX TENSOR
+ enddo
+ enddo
+ endif
+ hthree *= phase
+
+end
+
+subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
+
+ BEGIN_DOC
+ ! for double excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
+ !!
+ !! WARNING !!
+ !
+ ! Non hermitian !!
+ END_DOC
+
+ use bitmasks
+
+ implicit none
+ integer, intent(in) :: Nint
+ integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
+ double precision, intent(out) :: hthree
+ integer :: occ(Nint*bit_kind_size,2)
+ integer :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm
+ integer :: degree,exc(0:2,2,2)
+ integer :: h1, p1, h2, p2, s1, s2
+ double precision :: phase
+ integer :: other_spin(2)
+ integer(bit_kind) :: key_i_core(Nint,2)
+ double precision :: direct_int,exchange_int,sym_3_e_int_from_6_idx_tensor
+ double precision :: three_e_double_parrallel_spin
+
+ other_spin(1) = 2
+ other_spin(2) = 1
+
+ call get_excitation_degree(key_i, key_j, degree, Nint)
+
+ hthree = 0.d0
+
+ if(degree.ne.2)then
+ return
+ endif
+
+ if(core_tc_op)then
+ 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)
+ else
+ call bitstring_to_list_ab(key_i, occ, Ne, Nint)
+ endif
+ call get_double_excitation(key_i, key_j, exc, phase, Nint)
+ call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
+
+
+ if(Ne(1)+Ne(2).ge.3)then
+ if(s1==s2)then ! same spin excitation
+ ispin = other_spin(s1)
+! print*,'htilde ij'
+ do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1)
+ mm = occ(m,ispin)
+!! direct_int = three_body_ints_bi_ort(mm,p2,p1,mm,h2,h1)
+!! exchange_int = three_body_ints_bi_ort(mm,p2,p1,mm,h1,h2)
+ direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
+ exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
+! print*,direct_int,exchange_int
+ hthree += direct_int - exchange_int
+ enddo
+ do m = 1, Ne(s1) ! pure contribution from s1
+ mm = occ(m,s1)
+ hthree += three_e_double_parrallel_spin(mm,p2,h2,p1,h1)
+ enddo
+ else ! different spin excitation
+ do m = 1, Ne(s1)
+ mm = occ(m,s1) !
+ direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
+ exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1)
+ hthree += direct_int - exchange_int
+ enddo
+ do m = 1, Ne(s2)
+ mm = occ(m,s2) !
+ direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
+ exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1)
+ hthree += direct_int - exchange_int
+ enddo
+ endif
+ endif
+ hthree *= phase
+ end
diff --git a/src/tc_bi_ortho/symmetrized_3_e_int.irp.f b/src/tc_bi_ortho/symmetrized_3_e_int.irp.f
new file mode 100644
index 00000000..e4f7ca93
--- /dev/null
+++ b/src/tc_bi_ortho/symmetrized_3_e_int.irp.f
@@ -0,0 +1,111 @@
+subroutine give_all_perm_for_three_e(n,l,k,m,j,i,idx_list,phase)
+ implicit none
+ BEGIN_DOC
+ ! returns all the list of permutting indices for the antimmetrization of
+ !
+ ! (k^dagger l^dagger n^dagger m j i) when all indices have the same spins
+ !
+ ! idx_list(:,i) == list of the 6 indices corresponding the permutation "i"
+ !
+ ! phase(i) == phase of the permutation "i"
+ !
+ ! there are in total 6 permutations with different indices
+ END_DOC
+ integer, intent(in) :: n,l,k,m,j,i
+ integer, intent(out) :: idx_list(6,6)
+ double precision :: phase(6)
+ integer :: list(6)
+ !!! CYCLIC PERMUTATIONS
+ phase(1:3) = 1.d0
+ !!! IDENTITY PERMUTATION
+ list = (/n,l,k,m,j,i/)
+ idx_list(:,1) = list(:)
+ !!! FIRST CYCLIC PERMUTATION
+ list = (/n,l,k,j,i,m/)
+ idx_list(:,2) = list(:)
+ !!! FIRST CYCLIC PERMUTATION
+ list = (/n,l,k,i,m,j/)
+ idx_list(:,3) = list(:)
+
+ !!! NON CYCLIC PERMUTATIONS
+ phase(1:3) = -1.d0
+ !!! PARTICLE 1 is FIXED
+ list = (/n,l,k,j,m,i/)
+ idx_list(:,4) = list(:)
+ !!! PARTICLE 2 is FIXED
+ list = (/n,l,k,i,j,m/)
+ idx_list(:,5) = list(:)
+ !!! PARTICLE 3 is FIXED
+ list = (/n,l,k,m,i,j/)
+ idx_list(:,6) = list(:)
+
+end
+
+double precision function sym_3_e_int_from_6_idx_tensor(n,l,k,m,j,i)
+ implicit none
+ BEGIN_DOC
+ ! returns all good combinations of permutations of integrals with the good signs
+ !
+ ! for a given (k^dagger l^dagger n^dagger m j i) when all indices have the same spins
+ END_DOC
+ integer, intent(in) :: n,l,k,m,j,i
+ sym_3_e_int_from_6_idx_tensor = three_body_ints_bi_ort(n,l,k,m,j,i) & ! direct
+ + three_body_ints_bi_ort(n,l,k,j,i,m) & ! 1st cyclic permutation
+ + three_body_ints_bi_ort(n,l,k,i,m,j) & ! 2nd cyclic permutation
+ - three_body_ints_bi_ort(n,l,k,j,m,i) & ! elec 1 is kept fixed
+ - three_body_ints_bi_ort(n,l,k,i,j,m) & ! elec 2 is kept fixed
+ - three_body_ints_bi_ort(n,l,k,m,i,j) ! elec 3 is kept fixed
+
+end
+
+double precision function direct_sym_3_e_int(n,l,k,m,j,i)
+ implicit none
+ BEGIN_DOC
+ ! returns all good combinations of permutations of integrals with the good signs
+ !
+ ! for a given (k^dagger l^dagger n^dagger m j i) when all indices have the same spins
+ END_DOC
+ integer, intent(in) :: n,l,k,m,j,i
+ double precision :: integral
+ direct_sym_3_e_int = 0.d0
+ call give_integrals_3_body_bi_ort(n,l,k,m,j,i,integral) ! direct
+ direct_sym_3_e_int += integral
+ call give_integrals_3_body_bi_ort(n,l,k,j,i,m,integral) ! 1st cyclic permutation
+ direct_sym_3_e_int += integral
+ call give_integrals_3_body_bi_ort(n,l,k,i,m,j,integral) ! 2nd cyclic permutation
+ direct_sym_3_e_int += integral
+ call give_integrals_3_body_bi_ort(n,l,k,j,m,i,integral) ! elec 1 is kept fixed
+ direct_sym_3_e_int += -integral
+ call give_integrals_3_body_bi_ort(n,l,k,i,j,m,integral) ! elec 2 is kept fixed
+ direct_sym_3_e_int += -integral
+ call give_integrals_3_body_bi_ort(n,l,k,m,i,j,integral) ! elec 3 is kept fixed
+ direct_sym_3_e_int += -integral
+
+end
+
+double precision function three_e_diag_parrallel_spin(m,j,i)
+ implicit none
+ integer, intent(in) :: i,j,m
+ three_e_diag_parrallel_spin = three_e_3_idx_direct_bi_ort(m,j,i) ! direct
+ three_e_diag_parrallel_spin += three_e_3_idx_cycle_1_bi_ort(m,j,i) + three_e_3_idx_cycle_2_bi_ort(m,j,i) & ! two cyclic permutations
+ - three_e_3_idx_exch23_bi_ort(m,j,i) - three_e_3_idx_exch13_bi_ort(m,j,i) & ! two first exchange
+ - three_e_3_idx_exch12_bi_ort(m,j,i) ! last exchange
+end
+
+double precision function three_e_single_parrallel_spin(m,j,k,i)
+ implicit none
+ integer, intent(in) :: i,k,j,m
+ three_e_single_parrallel_spin = three_e_4_idx_direct_bi_ort(m,j,k,i) ! direct
+ three_e_single_parrallel_spin += three_e_4_idx_cycle_1_bi_ort(m,j,k,i) + three_e_4_idx_cycle_2_bi_ort(m,j,k,i) & ! two cyclic permutations
+ - three_e_4_idx_exch23_bi_ort(m,j,k,i) - three_e_4_idx_exch13_bi_ort(m,j,k,i) & ! two first exchange
+ - three_e_4_idx_exch12_bi_ort(m,j,k,i) ! last exchange
+end
+
+double precision function three_e_double_parrallel_spin(m,l,j,k,i)
+ implicit none
+ integer, intent(in) :: i,k,j,m,l
+ three_e_double_parrallel_spin = three_e_5_idx_direct_bi_ort(m,l,j,k,i) ! direct
+ three_e_double_parrallel_spin += three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) + three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) & ! two cyclic permutations
+ - three_e_5_idx_exch23_bi_ort(m,l,j,k,i) - three_e_5_idx_exch13_bi_ort(m,l,j,k,i) & ! two first exchange
+ - three_e_5_idx_exch12_bi_ort(m,l,j,k,i) ! last exchange
+end
diff --git a/src/tc_bi_ortho/tc_bi_ortho.irp.f b/src/tc_bi_ortho/tc_bi_ortho.irp.f
new file mode 100644
index 00000000..cfa24f3b
--- /dev/null
+++ b/src/tc_bi_ortho/tc_bi_ortho.irp.f
@@ -0,0 +1,61 @@
+program tc_bi_ortho
+ implicit none
+ BEGIN_DOC
+! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end.
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+ call routine_diag
+! call test
+end
+
+subroutine test
+ implicit none
+ integer :: i,j
+ double precision :: hmono,htwoe,hthree,htot
+ use bitmasks
+
+ print*,'test'
+! call htilde_mu_mat_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot)
+ call double_htilde_mu_mat_bi_ortho(N_int,psi_det(1,1,1), psi_det(1,1,2), hmono, htwoe, htot)
+ print*,hmono, htwoe, htot
+
+end
+
+subroutine routine_diag
+ implicit none
+! provide eigval_right_tc_bi_orth
+ provide overlap_bi_ortho
+! provide htilde_matrix_elmt_bi_ortho
+ integer ::i,j
+ print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1)
+ print*,'e_tc_left_right = ',e_tc_left_right
+ print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00
+ print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth
+ print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single
+ print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double
+ print*,'***'
+ print*,'e_corr_bi_orth = ',e_corr_bi_orth
+ print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj
+ print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth
+ print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth
+ print*,'Left/right eigenvectors'
+ do i = 1,N_det
+ write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1)
+ enddo
+ do j=1,N_states
+ do i=1,N_det
+ psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j)
+ psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j)
+ enddo
+ enddo
+ SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho
+ call save_tc_bi_ortho_wavefunction
+! call routine_save_left_right_bi_ortho
+end
+
diff --git a/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f b/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f
new file mode 100644
index 00000000..28f122ee
--- /dev/null
+++ b/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f
@@ -0,0 +1,24 @@
+program tc_bi_ortho_prop
+ implicit none
+ BEGIN_DOC
+! TODO : Put the documentation of the program here
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+! call routine_diag
+ call test
+end
+
+subroutine test
+ implicit none
+ integer :: i
+ print*,'TC Dipole components'
+ do i= 1, 3
+ print*,tc_bi_ortho_dipole(i,1)
+ enddo
+end
diff --git a/src/tc_bi_ortho/tc_cisd_sc2.irp.f b/src/tc_bi_ortho/tc_cisd_sc2.irp.f
new file mode 100644
index 00000000..0fb9f524
--- /dev/null
+++ b/src/tc_bi_ortho/tc_cisd_sc2.irp.f
@@ -0,0 +1,24 @@
+program tc_bi_ortho
+ implicit none
+ BEGIN_DOC
+! TODO : Put the documentation of the program here
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+ call test
+end
+
+subroutine test
+ implicit none
+! double precision, allocatable :: dressing_dets(:),e_corr_dets(:)
+! allocate(dressing_dets(N_det),e_corr_dets(N_det))
+! e_corr_dets = 0.d0
+! call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
+ provide eigval_tc_cisd_sc2_bi_ortho
+end
diff --git a/src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f b/src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f
new file mode 100644
index 00000000..406ee9e3
--- /dev/null
+++ b/src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f
@@ -0,0 +1,110 @@
+ BEGIN_PROVIDER [ double precision, reigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
+&BEGIN_PROVIDER [ double precision, leigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
+&BEGIN_PROVIDER [ double precision, eigval_tc_cisd_sc2_bi_ortho, (N_states)]
+ implicit none
+ integer :: it,n_real,degree,i
+ double precision :: e_before, e_current,thr, hmono,htwoe,hthree
+ double precision, allocatable :: e_corr_dets(:),h0j(:), h_sc2(:,:), dressing_dets(:)
+ double precision, allocatable :: leigvec_tc_bi_orth_tmp(:,:),reigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
+ allocate(leigvec_tc_bi_orth_tmp(N_det,N_det),reigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det))
+ allocate(e_corr_dets(N_det),h0j(N_det),h_sc2(N_det,N_det),dressing_dets(N_det))
+ do i = 1, N_det
+ call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
+ if(degree == 1 .or. degree == 2)then
+ call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,h0j(i))
+ endif
+ enddo
+ do i = 1, N_det
+ e_corr_dets(i) = reigvec_tc_bi_orth(i,1) * h0j(i)/reigvec_tc_bi_orth(1,1)
+ enddo
+ print*,'Starting from ',eigval_right_tc_bi_orth(1)
+
+ e_before = 0.d0
+ e_current = 10.d0
+ thr = 1.d-5
+ it = 0
+ dressing_dets = 0.d0
+ do while (dabs(E_before-E_current).gt.thr)
+ it += 1
+ E_before = E_current
+ h_sc2 = htilde_matrix_elmt_bi_ortho
+ call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
+ do i = 1, N_det
+ print*,'dressing_dets(i) = ',dressing_dets(i)
+ h_sc2(i,i) += dressing_dets(i)
+ enddo
+ call non_hrmt_real_diag(N_det,h_sc2,&
+ leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,&
+ n_real,eigval_right_tmp)
+ do i = 1, N_det
+ e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1)
+ enddo
+ E_current = eigval_right_tmp(1)
+ print*,'it, E(SC)^2 = ',it,E_current
+ enddo
+ eigval_tc_cisd_sc2_bi_ortho(1:N_states) = eigval_right_tmp(1:N_states)
+ reigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = reigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
+ leigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = leigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
+
+END_PROVIDER
+
+subroutine get_cisd_sc2_dressing(dets,e_corr_dets,ndet,dressing_dets)
+ implicit none
+ use bitmasks
+ integer, intent(in) :: ndet
+ integer(bit_kind), intent(in) :: dets(N_int,2,ndet)
+ double precision, intent(in) :: e_corr_dets(ndet)
+ double precision, intent(out) :: dressing_dets(ndet)
+ integer, allocatable :: degree(:),hole(:,:),part(:,:),spin(:,:)
+ integer(bit_kind), allocatable :: hole_part(:,:,:)
+ integer :: i,j,k, exc(0:2,2,2),h1,p1,h2,p2,s1,s2
+ integer(bit_kind) :: xorvec(2,N_int)
+
+ double precision :: phase
+ dressing_dets = 0.d0
+ allocate(degree(ndet),hole(2,ndet),part(2,ndet), spin(2,ndet),hole_part(N_int,2,ndet))
+ do i = 2, ndet
+ call get_excitation_degree(HF_bitmask,dets(1,1,i),degree(i),N_int)
+ do j = 1, N_int
+ hole_part(j,1,i) = xor( HF_bitmask(j,1), dets(j,1,i))
+ hole_part(j,2,i) = xor( HF_bitmask(j,2), dets(j,2,i))
+ enddo
+ if(degree(i) == 1)then
+ call get_single_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
+ else if(degree(i) == 2)then
+ call get_double_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
+ endif
+ call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+ hole(1,i) = h1
+ hole(2,i) = h2
+ part(1,i) = p1
+ part(2,i) = p2
+ spin(1,i) = s1
+ spin(2,i) = s2
+ enddo
+
+ integer :: same
+ if(elec_alpha_num+elec_beta_num<3)return
+ do i = 2, ndet
+ do j = i+1, ndet
+ same = 0
+ if(degree(i) == degree(j) .and. degree(i)==1)cycle
+ do k = 1, N_int
+ xorvec(k,1) = iand(hole_part(k,1,i),hole_part(k,1,j))
+ xorvec(k,2) = iand(hole_part(k,2,i),hole_part(k,2,j))
+ same += popcnt(xorvec(k,1)) + popcnt(xorvec(k,2))
+ enddo
+! print*,'i,j',i,j
+! call debug_det(dets(1,1,i),N_int)
+! call debug_det(hole_part(1,1,i),N_int)
+! call debug_det(dets(1,1,j),N_int)
+! call debug_det(hole_part(1,1,j),N_int)
+! print*,'same = ',same
+ if(same.eq.0)then
+ dressing_dets(i) += e_corr_dets(j)
+ dressing_dets(j) += e_corr_dets(i)
+ endif
+ enddo
+ enddo
+
+end
diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
new file mode 100644
index 00000000..651088be
--- /dev/null
+++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f
@@ -0,0 +1,179 @@
+ use bitmasks
+
+ BEGIN_PROVIDER [ integer, index_HF_psi_det]
+ implicit none
+ integer :: i,degree
+ do i = 1, N_det
+ call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
+ if(degree == 0)then
+ index_HF_psi_det = i
+ exit
+ endif
+ enddo
+ END_PROVIDER
+
+
+
+ BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth, (N_states)]
+&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth, (N_states)]
+&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth, (N_det,N_states)]
+&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth, (N_det,N_states)]
+&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ]
+
+ BEGIN_DOC
+ ! eigenvalues, right and left eigenvectors of the transcorrelated Hamiltonian on the BI-ORTHO basis
+ END_DOC
+
+ implicit none
+ integer :: i, idx_dress, j, istate
+ logical :: converged, dagger
+ integer :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l
+ double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:),leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
+
+ PROVIDE N_det N_int
+
+ if(n_det.le.N_det_max_full)then
+ allocate(reigvec_tc_bi_orth_tmp(N_det,N_det),leigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det))
+ call non_hrmt_real_diag(N_det,htilde_matrix_elmt_bi_ortho,&
+ leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,&
+ n_real_tc_bi_orth_eigval_right,eigval_right_tmp)
+ double precision, allocatable :: coef_hf_r(:),coef_hf_l(:)
+ integer, allocatable :: iorder(:)
+ allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(N_det))
+ do i = 1,N_det
+ iorder(i) = i
+ coef_hf_r(i) = -dabs(reigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
+ enddo
+ call dsort(coef_hf_r,iorder,N_det)
+ igood_r = iorder(1)
+ print*,'igood_r, coef_hf_r = ',igood_r,coef_hf_r(1)
+ do i = 1,N_det
+ iorder(i) = i
+ coef_hf_l(i) = -dabs(leigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
+ enddo
+ call dsort(coef_hf_l,iorder,N_det)
+ igood_l = iorder(1)
+ print*,'igood_l, coef_hf_l = ',igood_l,coef_hf_l(1)
+
+ if(igood_r.ne.igood_l.and.igood_r.ne.1)then
+ print *,''
+ print *,'Warning, the left and right eigenvectors are "not the same" '
+ print *,'Warning, the ground state is not dominated by HF...'
+ print *,'State with largest RIGHT coefficient of HF ',igood_r
+ print *,'coef of HF in RIGHT eigenvector = ',reigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_r)
+ print *,'State with largest LEFT coefficient of HF ',igood_l
+ print *,'coef of HF in LEFT eigenvector = ',leigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_l)
+ endif
+ if(state_following_tc)then
+ print *,'Following the states with the largest coef on HF'
+ print *,'igood_r,igood_l',igood_r,igood_l
+ i= igood_r
+ eigval_right_tc_bi_orth(1) = eigval_right_tmp(i)
+ do j = 1, N_det
+ reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i)
+! print*,reigvec_tc_bi_orth(j,1)
+ enddo
+ i= igood_l
+ eigval_left_tc_bi_orth(1) = eigval_right_tmp(i)
+ do j = 1, N_det
+ leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i)
+ enddo
+ else
+ do i = 1, N_states
+ eigval_right_tc_bi_orth(i) = eigval_right_tmp(i)
+ eigval_left_tc_bi_orth(i) = eigval_right_tmp(i)
+ do j = 1, N_det
+ reigvec_tc_bi_orth(j,i) = reigvec_tc_bi_orth_tmp(j,i)
+ leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i)
+ enddo
+ enddo
+ endif
+ else
+ double precision, allocatable :: H_jj(:),vec_tmp(:,:)
+ external htc_bi_ortho_calc_tdav
+ external htcdag_bi_ortho_calc_tdav
+ allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag))
+ do i = 1, N_det
+ call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
+ enddo
+ !!!! Preparing the left-eigenvector
+ print*,'Computing the left-eigenvector '
+ vec_tmp = 0.d0
+ do istate = 1, N_states
+ vec_tmp(:,istate) = psi_l_coef_bi_ortho(:,istate)
+ enddo
+ do istate = N_states+1, n_states_diag
+ vec_tmp(istate,istate) = 1.d0
+ enddo
+ call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav)
+ do istate = 1, N_states
+ leigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate)
+ enddo
+
+ print*,'Computing the right-eigenvector '
+ !!!! Preparing the right-eigenvector
+ vec_tmp = 0.d0
+ do istate = 1, N_states
+ vec_tmp(:,istate) = psi_r_coef_bi_ortho(:,istate)
+ enddo
+ do istate = N_states+1, n_states_diag
+ vec_tmp(istate,istate) = 1.d0
+ enddo
+ call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav)
+ do istate = 1, N_states
+ reigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate)
+ enddo
+
+ deallocate(H_jj)
+ endif
+ call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,N_det,N_det,N_states)
+ print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1)
+ norm_ground_left_right_bi_orth = 0.d0
+ do j = 1, N_det
+ norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1)
+ enddo
+ print*,'norm l/r = ',norm_ground_left_right_bi_orth
+
+END_PROVIDER
+
+
+
+subroutine bi_normalize(u_l,u_r,n,ld,nstates)
+ !!!! Normalization of the scalar product of the left/right eigenvectors
+ double precision, intent(inout) :: u_l(ld,nstates), u_r(ld,nstates)
+ integer, intent(in) :: n,ld,nstates
+ integer :: i
+ double precision :: accu, tmp
+ do i = 1, nstates
+ !!!! Normalization of right eigenvectors |Phi>
+ accu = 0.d0
+ do j = 1, n
+ accu += u_r(j,i) * u_r(j,i)
+ enddo
+ accu = 1.d0/dsqrt(accu)
+ print*,'accu_r = ',accu
+ do j = 1, n
+ u_r(j,i) *= accu
+ enddo
+ tmp = u_r(1,i) / dabs(u_r(1,i))
+ do j = 1, n
+ u_r(j,i) *= tmp
+ enddo
+ !!!! Adaptation of the norm of the left eigenvector such that = 1
+ accu = 0.d0
+ do j = 1, n
+ accu += u_l(j,i) * u_r(j,i)
+! print*,j, u_l(j,i) , u_r(j,i)
+ enddo
+ if(accu.gt.0.d0)then
+ accu = 1.d0/dsqrt(accu)
+ else
+ accu = 1.d0/dsqrt(-accu)
+ endif
+ tmp = (u_l(1,i) * u_r(1,i) )/dabs(u_l(1,i) * u_r(1,i))
+ do j = 1, n
+ u_l(j,i) *= accu * tmp
+ u_r(j,i) *= accu
+ enddo
+ enddo
+end
diff --git a/src/tc_bi_ortho/tc_hmat.irp.f b/src/tc_bi_ortho/tc_hmat.irp.f
new file mode 100644
index 00000000..bf1388e5
--- /dev/null
+++ b/src/tc_bi_ortho/tc_hmat.irp.f
@@ -0,0 +1,41 @@
+
+ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)]
+
+ BEGIN_DOC
+ ! htilde_matrix_elmt_bi_ortho(j,i) =
+ !
+ ! WARNING !!!!!!!!! IT IS NOT HERMITIAN !!!!!!!!!
+ END_DOC
+
+ implicit none
+ integer :: i, j
+ double precision :: hmono,htwoe,hthree,htot
+
+ PROVIDE N_int
+ !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hmono, htwoe, hthree, htot) &
+ !$OMP SHARED (N_det, psi_det, N_int,htilde_matrix_elmt_bi_ortho)
+ do i = 1, N_det
+ do j = 1, N_det
+ ! < J | Htilde | I >
+ call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
+ htilde_matrix_elmt_bi_ortho(j,i) = htot
+ enddo
+ enddo
+ !$OMP END PARALLEL DO
+! print*,'htilde_matrix_elmt_bi_ortho = '
+! do i = 1, min(100,N_det)
+! write(*,'(100(F16.10,X))')htilde_matrix_elmt_bi_ortho(1:min(100,N_det),i)
+! enddo
+
+
+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
diff --git a/src/tc_bi_ortho/tc_prop.irp.f b/src/tc_bi_ortho/tc_prop.irp.f
new file mode 100644
index 00000000..6db8e80e
--- /dev/null
+++ b/src/tc_bi_ortho/tc_prop.irp.f
@@ -0,0 +1,268 @@
+
+BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_states,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! tc_transition_matrix(p,h,istate,jstate) =
+ !
+ ! where are the left/right eigenvectors on a bi-ortho basis
+ END_DOC
+ integer :: i,j,istate,jstate,m,n,p,h
+ double precision :: phase
+ integer, allocatable :: occ(:,:)
+ integer :: n_occ_ab(2),degree,exc(0:2,2,2)
+ allocate(occ(N_int*bit_kind_size,2))
+ tc_transition_matrix = 0.d0
+ do istate = 1, N_states
+ do jstate = 1, N_states
+ do i = 1, N_det
+ do j = 1, N_det
+ call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+ if(degree.gt.1)then
+ cycle
+ else if (degree == 0)then
+ call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
+ do p = 1, n_occ_ab(1) ! browsing the alpha electrons
+ m = occ(p,1)
+ tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
+ enddo
+ do p = 1, n_occ_ab(2) ! browsing the beta electrons
+ m = occ(p,1)
+ tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
+ enddo
+ else
+ call get_single_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,phase,N_int)
+ if (exc(0,1,1) == 1) then
+ ! Single alpha
+ h = exc(1,1,1) ! hole in psi_det(1,1,j)
+ p = exc(1,2,1) ! particle in psi_det(1,1,j)
+ else
+ ! Single beta
+ h = exc(1,1,2) ! hole in psi_det(1,1,j)
+ p = exc(1,2,2) ! particle in psi_det(1,1,j)
+ endif
+ tc_transition_matrix(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+ END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_mo, (mo_num, mo_num)]
+ &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_mo, (mo_num, mo_num)]
+ &BEGIN_PROVIDER [ double precision, natorb_tc_eigval, (mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! natorb_tc_reigvec_mo : RIGHT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
+ ! natorb_tc_leigvec_mo : LEFT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
+ ! natorb_tc_eigval : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !!
+ END_DOC
+ double precision, allocatable :: dm_tmp(:,:)
+ integer :: i,j,k,n_real
+ allocate( dm_tmp(mo_num,mo_num))
+ dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1)
+ print*,'dm_tmp'
+ do i = 1, mo_num
+ write(*,'(100(F16.10,X))')-dm_tmp(:,i)
+ enddo
+! call non_hrmt_diag_split_degen( mo_num, dm_tmp&
+ call non_hrmt_fock_mat( mo_num, dm_tmp&
+! call non_hrmt_bieig( mo_num, dm_tmp&
+ , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo&
+ , n_real, natorb_tc_eigval )
+ double precision :: accu
+ accu = 0.d0
+ do i = 1, n_real
+ print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)
+ accu += -natorb_tc_eigval(i)
+ enddo
+ print*,'accu = ',accu
+ dm_tmp = 0.d0
+ do i = 1, n_real
+ accu = 0.d0
+ do k = 1, mo_num
+ accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i)
+ enddo
+ accu = 1.d0/dsqrt(dabs(accu))
+ natorb_tc_reigvec_mo(:,i) *= accu
+ natorb_tc_leigvec_mo(:,i) *= accu
+ do j = 1, n_real
+ do k = 1, mo_num
+ dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j)
+ enddo
+ enddo
+ enddo
+ double precision :: accu_d, accu_nd
+ accu_d = 0.d0
+ accu_nd = 0.d0
+ do i = 1, mo_num
+ accu_d += dm_tmp(i,i)
+ ! write(*,'(100(F16.10,X))')dm_tmp(:,i)
+ do j = 1, mo_num
+ if(i==j)cycle
+ accu_nd += dabs(dm_tmp(j,i))
+ enddo
+ enddo
+ print*,'Trace of the overlap between TC natural orbitals ',accu_d
+ print*,'L1 norm of extra diagonal elements of overlap matrix ',accu_nd
+
+
+ END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, fock_diag_sorted_r_natorb, (mo_num, mo_num)]
+ &BEGIN_PROVIDER [ double precision, fock_diag_sorted_l_natorb, (mo_num, mo_num)]
+ &BEGIN_PROVIDER [ double precision, fock_diag_sorted_v_natorb, (mo_num)]
+ implicit none
+ integer ::i,j,k
+ print*,'Diagonal elements of the Fock matrix before '
+ do i = 1, mo_num
+ write(*,*)i,Fock_matrix_tc_mo_tot(i,i)
+ enddo
+ double precision, allocatable :: fock_diag(:)
+ allocate(fock_diag(mo_num))
+ fock_diag = 0.d0
+ do i = 1, mo_num
+ fock_diag(i) = 0.d0
+ do j = 1, mo_num
+ do k = 1, mo_num
+ fock_diag(i) += natorb_tc_leigvec_mo(k,i) * Fock_matrix_tc_mo_tot(k,j) * natorb_tc_reigvec_mo(j,i)
+ enddo
+ enddo
+ enddo
+ integer, allocatable :: iorder(:)
+ allocate(iorder(mo_num))
+ do i = 1, mo_num
+ iorder(i) = i
+ enddo
+ call dsort(fock_diag,iorder,mo_num)
+ print*,'Diagonal elements of the Fock matrix after '
+ do i = 1, mo_num
+ write(*,*)i,fock_diag(i)
+ enddo
+ do i = 1, mo_num
+ fock_diag_sorted_v_natorb(i) = natorb_tc_eigval(iorder(i))
+ do j = 1, mo_num
+ fock_diag_sorted_r_natorb(j,i) = natorb_tc_reigvec_mo(j,iorder(i))
+ fock_diag_sorted_l_natorb(j,i) = natorb_tc_leigvec_mo(j,iorder(i))
+ enddo
+ enddo
+
+ END_PROVIDER
+
+
+
+ BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_ao, (ao_num, mo_num)]
+ &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_ao, (ao_num, mo_num)]
+ &BEGIN_PROVIDER [ double precision, overlap_natorb_tc_eigvec_ao, (mo_num, mo_num) ]
+
+ BEGIN_DOC
+ ! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP
+ !
+ ! THE OVERLAP SHOULD BE THE SAME AS overlap_natorb_tc_eigvec_mo
+ END_DOC
+
+ implicit none
+ integer :: i, j, k, q, p
+ double precision :: accu, accu_d
+ double precision, allocatable :: tmp(:,:)
+
+
+ ! ! MO_R x R
+ call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 &
+ , mo_r_coef, size(mo_r_coef, 1) &
+ , fock_diag_sorted_r_natorb, size(fock_diag_sorted_r_natorb, 1) &
+ , 0.d0, natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) )
+ !
+ ! MO_L x L
+ call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 &
+ , mo_l_coef, size(mo_l_coef, 1) &
+ , fock_diag_sorted_l_natorb, size(fock_diag_sorted_l_natorb, 1) &
+ , 0.d0, natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1) )
+
+
+ allocate( tmp(mo_num,ao_num) )
+
+ ! tmp <-- L.T x S_ao
+ call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 &
+ , natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) &
+ , 0.d0, tmp, size(tmp, 1) )
+
+ ! S <-- tmp x R
+ call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 &
+ , tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) &
+ , 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) )
+
+ deallocate( tmp )
+
+ ! ---
+ double precision :: norm
+ do i = 1, mo_num
+ norm = 1.d0/dsqrt(dabs(overlap_natorb_tc_eigvec_ao(i,i)))
+ do j = 1, mo_num
+ natorb_tc_reigvec_ao(j,i) *= norm
+ natorb_tc_leigvec_ao(j,i) *= norm
+ enddo
+ enddo
+
+ allocate( tmp(mo_num,ao_num) )
+
+ ! tmp <-- L.T x S_ao
+ call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 &
+ , natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) &
+ , 0.d0, tmp, size(tmp, 1) )
+
+ ! S <-- tmp x R
+ call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 &
+ , tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) &
+ , 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) )
+
+
+
+ deallocate( tmp )
+
+ accu_d = 0.d0
+ accu = 0.d0
+ do i = 1, mo_num
+ accu_d += overlap_natorb_tc_eigvec_ao(i,i)
+ do j = 1, mo_num
+ if(i==j)cycle
+ accu += dabs(overlap_natorb_tc_eigvec_ao(j,i))
+ enddo
+ enddo
+ print*,'Trace of the overlap_natorb_tc_eigvec_ao = ',accu_d
+ print*,'mo_num = ',mo_num
+ print*,'L1 norm of extra diagonal elements of overlap matrix ',accu
+ accu = accu / dble(mo_num**2)
+
+ END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, tc_bi_ortho_dipole, (3,N_states)]
+ implicit none
+ integer :: i,j,istate,m
+ double precision :: nuclei_part(3)
+ tc_bi_ortho_dipole = 0.d0
+ do istate = 1, N_states
+ do i = 1, mo_num
+ do j = 1, mo_num
+ tc_bi_ortho_dipole(1,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_x(j,i)
+ tc_bi_ortho_dipole(2,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_y(j,i)
+ tc_bi_ortho_dipole(3,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_z(j,i)
+ enddo
+ enddo
+ enddo
+
+ nuclei_part = 0.d0
+ do m = 1, 3
+ do i = 1,nucl_num
+ nuclei_part(m) += nucl_charge(i) * nucl_coord(i,m)
+ enddo
+ enddo
+!
+ do istate = 1, N_states
+ do m = 1, 3
+ tc_bi_ortho_dipole(m,istate) += nuclei_part(m)
+ enddo
+ enddo
+ END_PROVIDER
+
diff --git a/src/tc_bi_ortho/test_normal_order.irp.f b/src/tc_bi_ortho/test_normal_order.irp.f
new file mode 100644
index 00000000..8bdc57ee
--- /dev/null
+++ b/src/tc_bi_ortho/test_normal_order.irp.f
@@ -0,0 +1,73 @@
+program test_normal_order
+ implicit none
+ BEGIN_DOC
+! TODO : Put the documentation of the program here
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+ call test
+end
+
+subroutine test
+ implicit none
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+ integer :: h1,h2,p1,p2,s1,s2,i_ok,degree,Ne(2)
+ integer :: exc(0:2,2,2)
+ integer(bit_kind), allocatable :: det_i(:,:)
+ double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal
+ integer, allocatable :: occ(:,:)
+ allocate( occ(N_int*bit_kind_size,2) )
+ call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
+ allocate(det_i(N_int,2))
+ s1 = 1
+ s2 = 2
+ accu = 0.d0
+ do h1 = 1, elec_beta_num
+ do p1 = elec_beta_num+1, mo_num
+ do h2 = 1, elec_beta_num
+ do p2 = elec_beta_num+1, mo_num
+ det_i = ref_bitmask
+ call do_single_excitation(det_i,h1,p1,s1,i_ok)
+ call do_single_excitation(det_i,h2,p2,s2,i_ok)
+ call htilde_mu_mat_bi_ortho(det_i,HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+ call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
+ call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
+ hthree *= phase
+ normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1)
+ accu += dabs(hthree-normal)
+ enddo
+ enddo
+ enddo
+ enddo
+ print*,'accu opposite spin = ',accu
+
+ s1 = 2
+ s2 = 2
+ accu = 0.d0
+ do h1 = 1, elec_beta_num
+ do p1 = elec_beta_num+1, mo_num
+ do h2 = h1+1, elec_beta_num
+ do p2 = elec_beta_num+1, mo_num
+ det_i = ref_bitmask
+ call do_single_excitation(det_i,h1,p1,s1,i_ok)
+ call do_single_excitation(det_i,h2,p2,s2,i_ok)
+ if(i_ok.ne.1)cycle
+ call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+ call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
+ call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
+ hthree *= phase
+ normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
+ accu += dabs(hthree-normal)
+ enddo
+ enddo
+ enddo
+ enddo
+ print*,'accu same spin = ',accu
+end
+
+
diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f
new file mode 100644
index 00000000..2d71b6b2
--- /dev/null
+++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f
@@ -0,0 +1,131 @@
+program tc_bi_ortho
+ implicit none
+ BEGIN_DOC
+! TODO : Put the documentation of the program here
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+ ! call routine_2
+ call test_rout
+end
+
+subroutine test_rout
+ implicit none
+ integer :: i,j,ii,jj
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+ integer(bit_kind), allocatable :: det_i(:,:)
+ allocate(det_i(N_int,2))
+ det_i(:,:)= psi_det(:,:,1)
+ call debug_det(det_i,N_int)
+ integer, allocatable :: occ(:,:)
+ integer :: n_occ_ab(2)
+ allocate(occ(N_int*bit_kind_size,2))
+ call bitstring_to_list_ab(det_i, occ, n_occ_ab, N_int)
+ double precision :: hmono, htwoe, htot
+ call diag_htilde_mu_mat_bi_ortho(N_int, det_i, hmono, htwoe, htot)
+ print*,'hmono, htwoe, htot'
+ print*, hmono, htwoe, htot
+ print*,'alpha electrons orbital occupancy'
+ do i = 1, n_occ_ab(1) ! browsing the alpha electrons
+ j = occ(i,1)
+ print*,j,mo_bi_ortho_tc_one_e(j,j)
+ enddo
+ print*,'beta electrons orbital occupancy'
+ do i = 1, n_occ_ab(2) ! browsing the beta electrons
+ j = occ(i,2)
+ print*,j,mo_bi_ortho_tc_one_e(j,j)
+ enddo
+ print*,'alpha beta'
+ do i = 1, n_occ_ab(1)
+ ii = occ(i,1)
+ do j = 1, n_occ_ab(2)
+ jj = occ(j,2)
+ print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii)
+ enddo
+ enddo
+ print*,'alpha alpha'
+ do i = 1, n_occ_ab(1)
+ ii = occ(i,1)
+ do j = 1, n_occ_ab(1)
+ jj = occ(j,1)
+ print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii), mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
+ enddo
+ enddo
+
+ print*,'beta beta'
+ do i = 1, n_occ_ab(2)
+ ii = occ(i,2)
+ do j = 1, n_occ_ab(2)
+ jj = occ(j,2)
+ print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii), mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
+ enddo
+ enddo
+
+
+end
+
+subroutine routine_2
+ implicit none
+ integer :: i
+ double precision :: bi_ortho_mo_ints
+ print*,'H matrix'
+ do i = 1, N_det
+ write(*,'(1000(F16.5,X))')htilde_matrix_elmt_bi_ortho(:,i)
+ enddo
+ i = 1
+ double precision :: phase
+ integer :: degree,h1, p1, h2, p2, s1, s2, exc(0:2,2,2)
+ call get_excitation_degree(ref_bitmask, psi_det(1,1,i), degree, N_int)
+ if(degree==2)then
+ call get_double_excitation(ref_bitmask, psi_det(1,1,i), exc, phase, N_int)
+ call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
+ print*,'h1,h2,p1,p2'
+ print*, h1,h2,p1,p2
+ print*,mo_bi_ortho_tc_two_e(p1,p2,h1,h2),mo_bi_ortho_tc_two_e(h1,h2,p1,p2)
+ endif
+
+
+ print*,'coef'
+ do i = 1, ao_num
+ print*,i,mo_l_coef(i,8),mo_r_coef(i,8)
+ enddo
+! print*,'mdlqfmlqgmqglj'
+! print*,'mo_bi_ortho_tc_two_e()',mo_bi_ortho_tc_two_e(2,2,3,3)
+! print*,'bi_ortho_mo_ints ',bi_ortho_mo_ints(2,2,3,3)
+ print*,'Overlap'
+ do i = 1, mo_num
+ write(*,'(100(F16.10,X))')overlap_bi_ortho(:,i)
+ enddo
+
+end
+
+subroutine routine
+ implicit none
+ double precision :: hmono,htwoe,hthree,htot
+ integer(bit_kind), allocatable :: key1(:,:)
+ integer(bit_kind), allocatable :: key2(:,:)
+ allocate(key1(N_int,2),key2(N_int,2))
+ use bitmasks
+ key1 = ref_bitmask
+ call htilde_mu_mat_bi_ortho(key1,key1, N_int, hmono,htwoe,hthree,htot)
+ key2 = key1
+ integer :: h,p,i_ok
+ h = 1
+ p = 8
+ call do_single_excitation(key2,h,p,1,i_ok)
+ call debug_det(key2,N_int)
+ call htilde_mu_mat_bi_ortho(key2,key1, N_int, hmono,htwoe,hthree,htot)
+! print*,'fock_matrix_tc_mo_alpha(p,h) = ',fock_matrix_tc_mo_alpha(p,h)
+ print*,'htot = ',htot
+ print*,'hmono = ',hmono
+ print*,'htwoe = ',htwoe
+ double precision :: bi_ortho_mo_ints
+ print*,'bi_ortho_mo_ints(1,p,1,h)',bi_ortho_mo_ints(1,p,1,h)
+
+end
diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f
new file mode 100644
index 00000000..badc40b6
--- /dev/null
+++ b/src/tc_bi_ortho/test_tc_fock.irp.f
@@ -0,0 +1,169 @@
+program test_tc_fock
+ implicit none
+ BEGIN_DOC
+! TODO : Put the documentation of the program here
+ END_DOC
+ print *, 'Hello world'
+ my_grid_becke = .True.
+ my_n_pt_r_grid = 30
+ my_n_pt_a_grid = 50
+ read_wf = .True.
+ touch read_wf
+ touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+ !call routine_1
+ !call routine_2
+ call routine_3()
+
+end
+
+! ---
+
+subroutine routine_0
+ implicit none
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+ integer :: i,a,j,m,i_ok
+ integer :: exc(0:2,2,2),h1,p1,s1,h2,p2,s2,degree
+
+ integer(bit_kind), allocatable :: det_i(:,:)
+ double precision :: hmono,htwoe,hthree,htilde_ij,phase
+ double precision :: same, op, tot, accu
+ allocate(det_i(N_int,2))
+ s1 = 1
+ accu = 0.d0
+ do i = 1, elec_alpha_num ! occupied
+ do a = elec_alpha_num+1, mo_num ! virtual
+ det_i = ref_bitmask
+ call do_single_excitation(det_i,i,a,s1,i_ok)
+ if(i_ok == -1)then
+ print*,'PB !!'
+ print*,i,a
+ stop
+ endif
+! call debug_det(det_i,N_int)
+ call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
+ call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+ op = fock_3_mat_a_op_sh_bi_orth(a,i)
+ same = fock_3_mat_a_sa_sh_bi_orth(a,i)
+! same = 0.d0
+ tot = same + op
+ if(dabs(tot - phase*hthree).gt.1.d-10)then
+ print*,'------'
+ print*,i,a,phase
+ print*,'hthree = ',phase*hthree
+ print*,'fock = ',tot
+ print*,'same,op= ',same,op
+ print*,dabs(tot - phase*hthree)
+ stop
+ endif
+ accu += dabs(tot - phase*hthree)
+ enddo
+ enddo
+ print*,'accu = ',accu
+
+end subroutine routine_0
+
+! ---
+
+subroutine routine_1
+
+ implicit none
+ integer :: i, a
+ double precision :: accu
+
+ accu = 0.d0
+ do i = 1, mo_num
+ do a = 1, mo_num
+ accu += dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) )
+ !if(dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) .gt. 1.d-10)then
+ print*, i, a
+ print*, dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) &
+ , fock_3_mat_a_op_sh_bi_orth_old(a,i), fock_3_mat_a_op_sh_bi_orth(a,i)
+ !endif
+ enddo
+ enddo
+
+ print *, 'accu = ', accu
+
+end subroutine routine_1
+
+! ---
+
+subroutine routine_2
+
+ implicit none
+ integer :: i, a
+ double precision :: accu
+
+ accu = 0.d0
+ do i = 1, mo_num
+ do a = 1, mo_num
+ accu += dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) )
+ !if(dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) .gt. 1.d-10)then
+ print*, i, a
+ print*, dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) &
+ , fock_3_mat_a_sa_sh_bi_orth_old(a,i), fock_3_mat_a_sa_sh_bi_orth(a,i)
+ !endif
+ enddo
+ enddo
+
+ print *, 'accu = ', accu
+
+end subroutine routine_2
+
+! ---
+
+subroutine routine_3()
+
+ use bitmasks ! you need to include the bitmasks_module.f90 features
+
+ implicit none
+ integer :: i, a, i_ok, s1
+ double precision :: hmono, htwoe, hthree, htilde_ij
+ double precision :: err_ai, err_tot
+ integer(bit_kind), allocatable :: det_i(:,:)
+
+ allocate(det_i(N_int,2))
+
+ err_tot = 0.d0
+
+ s1 = 1
+
+ det_i = ref_bitmask
+ call debug_det(det_i, N_int)
+ print*, ' HF det'
+ call debug_det(det_i, N_int)
+
+ do i = 1, elec_alpha_num ! occupied
+ do a = elec_alpha_num+1, mo_num ! virtual
+
+
+ det_i = ref_bitmask
+ call do_single_excitation(det_i, i, a, s1, i_ok)
+ if(i_ok == -1) then
+ print*, 'PB !!'
+ print*, i, a
+ stop
+ endif
+ !print*, ' excited det'
+ !call debug_det(det_i, N_int)
+
+ call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
+ err_ai = dabs(htilde_ij)
+ if(err_ai .gt. 1d-7) then
+ print*, ' warning on', i, a
+ print*, hmono, htwoe, htilde_ij
+ endif
+ err_tot += err_ai
+
+ write(22, *) htilde_ij
+ enddo
+ enddo
+
+ print *, ' err_tot = ', err_tot
+
+ deallocate(det_i)
+
+end subroutine routine_3
+
+! ---