tc_scf compiles and gives good energy for Ne. Added a test in test_Ne.sh

src/tc_scf/EZFIO.cfg

@@ -0,0 +1,4 @@
+[bitc_energy]
+type: Threshold
+doc: Energy bi-tc HF
+interface: ezfio

src/tc_scf/NEED

@@ -0,0 +1,6 @@
+hartree_fock
+bi_ortho_mos
+three_body_ints
+bi_ort_ints
+tc_keywords
+non_hermit_dav

src/tc_scf/combine_lr_tcscf.irp.f

@@ -0,0 +1,74 @@
+
+! ---
+
+program combine_lr_tcscf
+
+  BEGIN_DOC
+  ! TODO : Put the documentation of the program here
+  END_DOC
+
+  implicit none
+
+  my_grid_becke  = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  bi_ortho = .True.
+  touch bi_ortho
+
+  call comb_orbitals()
+
+end
+
+! ---
+
+subroutine comb_orbitals()
+
+  implicit none
+  integer                       :: i, m, n, nn, mm
+  double precision              :: accu_d, accu_nd
+  double precision, allocatable :: R(:,:), L(:,:), Rnew(:,:), tmp(:,:), S(:,:)
+
+  n  = ao_num
+  m  = mo_num
+  nn = elec_alpha_num
+  mm = m - nn
+
+  allocate(L(n,m), R(n,m), Rnew(n,m), S(m,m))
+  L = mo_l_coef
+  R = mo_r_coef
+
+  call check_weighted_biorthog(n, m, ao_overlap, L, R, accu_d, accu_nd, S, .true.)
+
+  allocate(tmp(n,nn))
+  do i = 1, nn 
+    tmp(1:n,i) = R(1:n,i)
+  enddo
+  call impose_weighted_orthog_svd(n, nn, ao_overlap, tmp)
+  do i = 1, nn
+    Rnew(1:n,i) = tmp(1:n,i)
+  enddo
+  deallocate(tmp)
+
+  allocate(tmp(n,mm))
+  do i = 1, mm
+    tmp(1:n,i) = L(1:n,i+nn)
+  enddo
+  call impose_weighted_orthog_svd(n, mm, ao_overlap, tmp)
+  do i = 1, mm
+    Rnew(1:n,i+nn) = tmp(1:n,i)
+  enddo
+  deallocate(tmp)
+
+  call check_weighted_biorthog(n, m, ao_overlap, Rnew, Rnew, accu_d, accu_nd, S, .true.)
+
+  mo_r_coef = Rnew
+  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+
+  deallocate(L, R, Rnew, S)
+
+end subroutine comb_orbitals
+
+! ---
+

src/tc_scf/diago_bi_ort_tcfock.irp.f

@@ -0,0 +1,229 @@
+! ---
+
+ BEGIN_PROVIDER [ double precision, fock_tc_reigvec_mo, (mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, fock_tc_leigvec_mo, (mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, eigval_fock_tc_mo, (mo_num)]
+&BEGIN_PROVIDER [ double precision, overlap_fock_tc_eigvec_mo, (mo_num, mo_num)]
+
+  BEGIN_DOC
+  ! EIGENVECTORS OF FOCK MATRIX ON THE MO BASIS and their OVERLAP
+  END_DOC
+
+  implicit none
+  integer                       :: n_real_tc 
+  integer                       :: i, j, k, l
+  double precision              :: accu_d, accu_nd, accu_tmp
+  double precision              :: norm
+  double precision, allocatable :: eigval_right_tmp(:)
+  double precision, allocatable :: F_tmp(:,:)
+
+  allocate( eigval_right_tmp(mo_num), F_tmp(mo_num,mo_num) )
+
+  PROVIDE Fock_matrix_tc_mo_tot
+
+  do i = 1, mo_num
+    do j = 1, mo_num
+      F_tmp(j,i) = Fock_matrix_tc_mo_tot(j,i)
+    enddo
+  enddo
+  ! insert level shift here
+  do i = elec_beta_num+1, elec_alpha_num
+    F_tmp(i,i) += 0.5d0 * level_shift_tcscf
+  enddo
+  do i = elec_alpha_num+1, mo_num
+    F_tmp(i,i) += level_shift_tcscf
+  enddo
+
+  call non_hrmt_bieig( mo_num, F_tmp, thresh_biorthog_diag, thresh_biorthog_nondiag &
+                     , fock_tc_leigvec_mo, fock_tc_reigvec_mo                       & 
+                     , n_real_tc, eigval_right_tmp )
+
+  !if(max_ov_tc_scf)then
+  ! call non_hrmt_fock_mat( mo_num, F_tmp, thresh_biorthog_diag, thresh_biorthog_nondiag &
+  !                    , fock_tc_leigvec_mo, fock_tc_reigvec_mo                          & 
+  !                    , n_real_tc, eigval_right_tmp )
+  !else 
+  ! call non_hrmt_diag_split_degen_bi_orthog( mo_num, F_tmp     &
+  !                    , fock_tc_leigvec_mo, fock_tc_reigvec_mo & 
+  !                    , n_real_tc, eigval_right_tmp )
+  !endif
+
+  deallocate(F_tmp)
+
+
+!  if(n_real_tc .ne. mo_num)then
+!   print*,'n_real_tc ne mo_num ! ',n_real_tc
+!   stop
+!  endif
+
+  eigval_fock_tc_mo = eigval_right_tmp
+!  print*,'Eigenvalues of Fock_matrix_tc_mo_tot'
+!  do i = 1, elec_alpha_num
+!    print*, i, eigval_fock_tc_mo(i)
+!  enddo
+!  do i = elec_alpha_num+1, mo_num 
+!    print*, i, eigval_fock_tc_mo(i) - level_shift_tcscf
+!  enddo
+!  deallocate( eigval_right_tmp )
+
+  ! L.T x R 
+  call dgemm( "T", "N", mo_num, mo_num, mo_num, 1.d0          &
+            , fock_tc_leigvec_mo, size(fock_tc_leigvec_mo, 1) &
+            , fock_tc_reigvec_mo, size(fock_tc_reigvec_mo, 1) &
+            , 0.d0, overlap_fock_tc_eigvec_mo, size(overlap_fock_tc_eigvec_mo, 1) )
+
+  ! ---
+
+  accu_d  = 0.d0
+  accu_nd = 0.d0
+  do i = 1, mo_num
+    do k = 1, mo_num
+      if(i==k) then
+        accu_tmp = overlap_fock_tc_eigvec_mo(k,i)
+        accu_d += dabs(accu_tmp )
+      else
+        accu_tmp = overlap_fock_tc_eigvec_mo(k,i)
+        accu_nd += accu_tmp * accu_tmp
+        if(dabs(overlap_fock_tc_eigvec_mo(k,i)) .gt. thresh_biorthog_nondiag)then
+         print *, 'k,i', k, i, overlap_fock_tc_eigvec_mo(k,i)
+        endif
+      endif
+    enddo 
+  enddo
+  accu_nd = dsqrt(accu_nd) / accu_d
+  if(accu_nd .gt. thresh_biorthog_nondiag) then
+    print *, ' bi-orthog failed'
+    print *, ' accu_nd MO = ', accu_nd, thresh_biorthog_nondiag
+    print *, ' overlap_fock_tc_eigvec_mo = '
+    do i = 1, mo_num
+      write(*,'(100(F16.10,X))') overlap_fock_tc_eigvec_mo(i,:)
+    enddo
+    stop
+  endif
+
+  ! ---
+
+  if(dabs(accu_d - dble(mo_num))/dble(mo_num) .gt. thresh_biorthog_diag) then
+
+    print *, ' mo_num     = ', mo_num 
+    print *, ' accu_d  MO = ', accu_d, thresh_biorthog_diag
+    print *, ' normalizing vectors ...'
+    do i = 1, mo_num
+      norm = dsqrt(dabs(overlap_fock_tc_eigvec_mo(i,i)))
+      if(norm .gt. thresh_biorthog_diag) then
+        do k = 1, mo_num
+          fock_tc_reigvec_mo(k,i) *= 1.d0/norm
+          fock_tc_leigvec_mo(k,i) *= 1.d0/norm
+        enddo
+      endif
+    enddo
+
+    call dgemm( "T", "N", mo_num, mo_num, mo_num, 1.d0          &
+              , fock_tc_leigvec_mo, size(fock_tc_leigvec_mo, 1) &
+              , fock_tc_reigvec_mo, size(fock_tc_reigvec_mo, 1) &
+              , 0.d0, overlap_fock_tc_eigvec_mo, size(overlap_fock_tc_eigvec_mo, 1) )
+
+    accu_d  = 0.d0
+    accu_nd = 0.d0
+    do i = 1, mo_num
+      do k = 1, mo_num
+        if(i==k) then
+          accu_tmp = overlap_fock_tc_eigvec_mo(k,i)
+          accu_d  += dabs(accu_tmp)
+        else
+          accu_tmp = overlap_fock_tc_eigvec_mo(k,i)
+          accu_nd += accu_tmp * accu_tmp
+          if(dabs(overlap_fock_tc_eigvec_mo(k,i)) .gt. thresh_biorthog_nondiag)then
+           print *, 'k,i', k, i, overlap_fock_tc_eigvec_mo(k,i)
+          endif
+        endif
+      enddo 
+    enddo
+    accu_nd = dsqrt(accu_nd) / accu_d
+    if(accu_nd .gt. thresh_biorthog_diag) then
+      print *, ' bi-orthog failed'
+      print *, ' accu_nd MO = ', accu_nd, thresh_biorthog_nondiag
+      print *, ' overlap_fock_tc_eigvec_mo = '
+      do i = 1, mo_num
+        write(*,'(100(F16.10,X))') overlap_fock_tc_eigvec_mo(i,:)
+      enddo
+      stop
+    endif
+
+  endif
+ 
+  ! ---
+
+END_PROVIDER 
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, fock_tc_reigvec_ao, (ao_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, fock_tc_leigvec_ao, (ao_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, overlap_fock_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_fock_tc_eigvec_mo
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, k, q, p
+  double precision              :: accu, accu_d
+  double precision, allocatable :: tmp(:,:)
+
+  PROVIDE mo_l_coef mo_r_coef
+
+!  ! 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_tc_reigvec_mo, size(fock_tc_reigvec_mo, 1) &
+             , 0.d0, fock_tc_reigvec_ao, size(fock_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_tc_leigvec_mo, size(fock_tc_leigvec_mo, 1) &
+             , 0.d0, fock_tc_leigvec_ao, size(fock_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                                           &
+            , fock_tc_leigvec_ao, size(fock_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), fock_tc_reigvec_ao, size(fock_tc_reigvec_ao, 1) &
+            , 0.d0, overlap_fock_tc_eigvec_ao, size(overlap_fock_tc_eigvec_ao, 1) )
+
+  deallocate( tmp )
+
+  ! ---
+  double precision :: norm
+  do i = 1, mo_num
+   norm = 1.d0/dsqrt(dabs(overlap_fock_tc_eigvec_ao(i,i)))
+   do j = 1, mo_num
+    fock_tc_reigvec_ao(j,i) *= norm
+    fock_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                                           &
+            , fock_tc_leigvec_ao, size(fock_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), fock_tc_reigvec_ao, size(fock_tc_reigvec_ao, 1) &
+            , 0.d0, overlap_fock_tc_eigvec_ao, size(overlap_fock_tc_eigvec_ao, 1) )
+
+  deallocate( tmp )
+
+END_PROVIDER
+

src/tc_scf/diis_tcscf.irp.f

@@ -0,0 +1,186 @@
+! ---
+
+BEGIN_PROVIDER [ double precision, threshold_DIIS_nonzero_TCSCF ]
+ 
+  implicit none
+ 
+  if(threshold_DIIS_TCSCF == 0.d0) then
+    threshold_DIIS_nonzero_TCSCF = dsqrt(thresh_tcscf)
+  else
+    threshold_DIIS_nonzero_TCSCF = threshold_DIIS_TCSCF
+  endif
+  ASSERT(threshold_DIIS_nonzero_TCSCF >= 0.d0)
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [double precision, Q_alpha, (ao_num, ao_num) ]
+
+  BEGIN_DOC
+  !
+  ! Q_alpha = mo_r_coef x eta_occ_alpha x mo_l_coef.T
+  !
+  ! [Q_alpha]_ij = \sum_{k=1}^{elec_alpha_num} [mo_r_coef]_ik [mo_l_coef]_jk
+  !
+  END_DOC
+
+  implicit none
+
+  Q_alpha = 0.d0
+  call dgemm( 'N', 'T', ao_num, ao_num, elec_alpha_num, 1.d0               &
+            , mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
+            , 0.d0, Q_alpha, size(Q_alpha, 1) )
+
+END_PROVIDER
+  
+! ---
+    
+BEGIN_PROVIDER [ double precision, Q_beta, (ao_num, ao_num) ]
+
+  BEGIN_DOC
+  !
+  ! Q_beta = mo_r_coef x eta_occ_beta x mo_l_coef.T
+  !
+  ! [Q_beta]_ij = \sum_{k=1}^{elec_beta_num} [mo_r_coef]_ik [mo_l_coef]_jk
+  !
+  END_DOC
+
+  implicit none
+
+  Q_beta = 0.d0
+  call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0                &
+            , mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
+            , 0.d0, Q_beta, size(Q_beta, 1) )
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, Q_matrix, (ao_num, ao_num) ]
+
+  BEGIN_DOC
+  !
+  ! Q_matrix = 2 mo_r_coef x eta_occ x mo_l_coef.T
+  ! 
+  ! with: 
+  !                        | 1   if i = j = 1, ..., nb of occ orbitals
+  !        [eta_occ]_ij =  |     
+  !                        | 0   otherwise
+  !
+  ! the diis error is defines as:
+  !                         e = F_ao x Q x ao_overlap - ao_overlap x Q x F_ao
+  ! with: 
+  !       mo_l_coef.T x ao_overlap x mo_r_coef = I
+  !       F_mo = mo_l_coef.T x F_ao x mo_r_coef
+  !       F_ao = (ao_overlap x mo_r_coef) x F_mo x (ao_overlap x mo_l_coef).T
+  !
+  ! ==> e = 2 ao_overlap x mo_r_coef x [ F_mo x eta_occ - eta_occ x F_mo ] x (ao_overlap x mo_l_coef).T
+  !
+  !      at convergence:
+  !                                      F_mo x eta_occ - eta_occ x F_mo = 0
+  !                                  ==> [F_mo]_ij ([eta_occ]_ii - [eta_occ]_jj) = 0  
+  !                                  ==> [F_mo]_ia = [F_mo]_ai = 0 where: i = occ and a = vir
+  !                                  ==> Brillouin conditions
+  !
+  END_DOC
+
+  implicit none
+
+  if(elec_alpha_num == elec_beta_num) then
+    Q_matrix = Q_alpha + Q_alpha
+  else
+    Q_matrix = Q_alpha + Q_beta
+  endif
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)]
+
+  implicit none
+  double precision, allocatable :: tmp(:,:)
+
+  allocate(tmp(ao_num,ao_num))
+
+  ! F x Q
+  call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0                                             &
+            , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), Q_matrix, size(Q_matrix, 1) &
+            , 0.d0, tmp, size(tmp, 1) )
+
+  ! F x Q x S
+  call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0             &
+            , tmp, size(tmp, 1), ao_overlap, size(ao_overlap, 1) &
+            , 0.d0, FQS_SQF_ao, size(FQS_SQF_ao, 1) )
+
+  ! S x Q
+  tmp = 0.d0
+  call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0                       &
+            , ao_overlap, size(ao_overlap, 1), Q_matrix, size(Q_matrix, 1) &
+            , 0.d0, tmp, size(tmp, 1) )
+
+  ! F x Q x S - S x Q x F
+  call dgemm( 'N', 'N', ao_num, ao_num, ao_num, -1.d0                                  &
+            , tmp, size(tmp, 1), Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) &
+            , 1.d0, FQS_SQF_ao, size(FQS_SQF_ao, 1) )
+
+  deallocate(tmp)
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [double precision, FQS_SQF_mo, (mo_num, mo_num)]
+
+  implicit none
+
+  call ao_to_mo_bi_ortho( FQS_SQF_ao, size(FQS_SQF_ao, 1) &
+                        , FQS_SQF_mo, size(FQS_SQF_mo, 1) )
+
+END_PROVIDER
+
+! ---
+
+! BEGIN_PROVIDER [ double precision, eigenval_Fock_tc_ao, (ao_num) ]
+!&BEGIN_PROVIDER [ double precision, eigenvec_Fock_tc_ao, (ao_num,ao_num) ]
+!
+!  BEGIN_DOC
+!  !
+!  ! Eigenvalues and eigenvectors of the Fock matrix over the ao basis
+!  !
+!  ! F' = X.T x F x X   where X = ao_overlap^(-1/2)
+!  ! 
+!  ! F'   x Cr' = Cr' x E ==> F   Cr = Cr x E with Cr = X x Cr'
+!  ! F'.T x Cl' = Cl' x E ==> F.T Cl = Cl x E with Cl = X x Cl'
+!  !
+!  END_DOC
+!
+!  implicit none
+!  double precision, allocatable :: tmp1(:,:), tmp2(:,:)
+!
+!  ! ---
+!  ! Fock matrix in orthogonal basis: F' = X.T x F x X
+!
+!  allocate(tmp1(ao_num,ao_num))
+!  call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0                                                 &
+!            , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), S_half_inv, size(S_half_inv, 1) &
+!            , 0.d0, tmp1, size(tmp1, 1) )
+!
+!  allocate(tmp2(ao_num,ao_num))
+!  call dgemm( 'T', 'N', ao_num, ao_num, ao_num, 1.d0               &
+!            , S_half_inv, size(S_half_inv, 1), tmp1, size(tmp1, 1) &
+!            , 0.d0, tmp2, size(tmp2, 1) )
+!
+!  ! ---
+!
+!  ! Diagonalize F' to obtain eigenvectors in orthogonal basis C' and eigenvalues
+!  ! TODO
+!
+!  ! Back-transform eigenvectors: C =X.C'
+!
+!END_PROVIDER
+
+! ---
+
+~                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    

src/tc_scf/fock_3e_bi_ortho_uhf.irp.f

@@ -0,0 +1,405 @@
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
+
+  implicit none
+  integer          :: a, b, i, j
+  double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia
+  double precision :: ti, tf
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  !print *, ' PROVIDING fock_3e_uhf_mo_cs ...'
+  call wall_time(ti)
+
+  fock_3e_uhf_mo_cs = 0.d0
+
+  do a = 1, mo_num
+    do b = 1, mo_num
+   
+      do j = 1, elec_beta_num
+        do i = 1, elec_beta_num
+
+          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
+          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
+          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
+          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
+          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
+          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
+
+          fock_3e_uhf_mo_cs(b,a) -= 0.5d0 * ( 4.d0 * I_bij_aij &
+                                            +        I_bij_ija &
+                                            +        I_bij_jai &
+                                            - 2.d0 * I_bij_aji &
+                                            - 2.d0 * I_bij_iaj &
+                                            - 2.d0 * I_bij_jia )
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  call wall_time(tf)
+  !print *, ' total Wall time for fock_3e_uhf_mo_cs =', tf - ti
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
+
+  implicit none
+  integer          :: a, b, i, j, o
+  double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia
+  double precision :: ti, tf
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  !print *, ' PROVIDING fock_3e_uhf_mo_a ...'
+  call wall_time(ti)
+
+  o = elec_beta_num + 1
+
+  fock_3e_uhf_mo_a = fock_3e_uhf_mo_cs
+
+  do a = 1, mo_num
+    do b = 1, mo_num
+
+      ! ---
+
+      do j = o, elec_alpha_num
+        do i = 1, elec_beta_num
+
+          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
+          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
+          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
+          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
+          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
+          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
+
+          fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
+                                           +        I_bij_ija &
+                                           +        I_bij_jai &
+                                           -        I_bij_aji &
+                                           -        I_bij_iaj &
+                                           - 2.d0 * I_bij_jia )
+
+        enddo
+      enddo
+
+      ! ---
+
+      do j = 1, elec_beta_num
+        do i = o, elec_alpha_num
+
+          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
+          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
+          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
+          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
+          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
+          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
+
+          fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
+                                           +        I_bij_ija &
+                                           +        I_bij_jai &
+                                           -        I_bij_aji &
+                                           - 2.d0 * I_bij_iaj &
+                                           -        I_bij_jia )
+
+        enddo
+      enddo
+
+      ! ---
+
+      do j = o, elec_alpha_num
+        do i = o, elec_alpha_num
+
+          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
+          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
+          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
+          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
+          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
+          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
+
+          fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( I_bij_aij &
+                                           + I_bij_ija &
+                                           + I_bij_jai &
+                                           - I_bij_aji &
+                                           - I_bij_iaj &
+                                           - I_bij_jia )
+
+        enddo
+      enddo
+
+      ! ---
+
+    enddo
+  enddo
+
+  call wall_time(tf)
+  !print *, ' total Wall time for fock_3e_uhf_mo_a =', tf - ti
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)]
+
+  implicit none
+  integer          :: a, b, i, j, o
+  double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia
+  double precision :: ti, tf
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  !print *, ' PROVIDING fock_3e_uhf_mo_b ...'
+  call wall_time(ti)
+
+  o = elec_beta_num + 1
+
+  fock_3e_uhf_mo_b = fock_3e_uhf_mo_cs
+
+  do a = 1, mo_num
+    do b = 1, mo_num
+
+      ! ---
+
+      do j = o, elec_alpha_num
+        do i = 1, elec_beta_num
+
+          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
+          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
+          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
+          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
+          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
+          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
+
+          fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
+                                           -        I_bij_aji &
+                                           -        I_bij_iaj )
+
+        enddo
+      enddo
+
+      ! ---
+
+      do j = 1, elec_beta_num
+        do i = o, elec_alpha_num
+
+          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
+          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
+          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
+          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
+          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
+          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
+
+          fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
+                                           -        I_bij_aji &
+                                           -        I_bij_jia )
+
+        enddo
+      enddo
+
+      ! ---
+
+      do j = o, elec_alpha_num
+        do i = o, elec_alpha_num
+
+          call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij)
+          call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija)
+          call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai)
+          call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji)
+          call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj)
+          call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia)
+
+          fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( I_bij_aij &
+                                           - I_bij_aji )
+
+        enddo
+      enddo
+
+      ! ---
+
+    enddo
+  enddo
+
+  call wall_time(tf)
+  !print *, ' total Wall time for fock_3e_uhf_mo_b =', tf - ti
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_a, (ao_num, ao_num)]
+
+  BEGIN_DOC
+  !
+  ! Equations (B6) and (B7)
+  !
+  ! g <--> gamma
+  ! d <--> delta
+  ! e <--> eta
+  ! k <--> kappa
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: g, d, e, k, mu, nu
+  double precision              :: dm_ge_a, dm_ge_b, dm_ge
+  double precision              :: dm_dk_a, dm_dk_b, dm_dk
+  double precision              :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu
+  double precision              :: ti, tf
+  double precision, allocatable :: f_tmp(:,:)
+
+  print *, ' PROVIDING fock_3e_uhf_ao_a ...'
+  call wall_time(ti)
+
+  fock_3e_uhf_ao_a = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                                                &
+ !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp,  &
+ !$OMP          i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) &
+ !$OMP SHARED  (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_a)
+
+  allocate(f_tmp(ao_num,ao_num))
+  f_tmp = 0.d0
+
+ !$OMP DO
+  do g = 1, ao_num
+    do e = 1, ao_num
+      dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e)
+      dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e)
+      dm_ge   = dm_ge_a + dm_ge_b
+      do d = 1, ao_num
+        do k = 1, ao_num
+          dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k)
+          dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k)
+          dm_dk   = dm_dk_a + dm_dk_b
+          do mu = 1, ao_num
+            do nu = 1, ao_num
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu)
+              f_tmp(mu,nu) -= 0.5d0 * ( dm_ge   * dm_dk   * i_mugd_nuek &
+                                      + dm_ge_a * dm_dk_a * i_mugd_eknu &
+                                      + dm_ge_a * dm_dk_a * i_mugd_knue &
+                                      - dm_ge_a * dm_dk   * i_mugd_enuk &
+                                      - dm_ge   * dm_dk_a * i_mugd_kenu &
+                                      - dm_ge_a * dm_dk_a * i_mugd_nuke &
+                                      - dm_ge_b * dm_dk_b * i_mugd_nuke )
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+ !$OMP END DO NOWAIT
+
+ !$OMP CRITICAL
+  do mu = 1, ao_num
+    do nu = 1, ao_num
+      fock_3e_uhf_ao_a(mu,nu) += f_tmp(mu,nu)
+    enddo
+  enddo
+ !$OMP END CRITICAL
+
+  deallocate(f_tmp)
+ !$OMP END PARALLEL
+
+  call wall_time(tf)
+  print *, ' total Wall time for fock_3e_uhf_ao_a =', tf - ti
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_b, (ao_num, ao_num)]
+
+  BEGIN_DOC
+  !
+  ! Equations (B6) and (B7)
+  !
+  ! g <--> gamma
+  ! d <--> delta
+  ! e <--> eta
+  ! k <--> kappa
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: g, d, e, k, mu, nu
+  double precision              :: dm_ge_a, dm_ge_b, dm_ge
+  double precision              :: dm_dk_a, dm_dk_b, dm_dk
+  double precision              :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu
+  double precision              :: ti, tf
+  double precision, allocatable :: f_tmp(:,:)
+
+  print *, ' PROVIDING fock_3e_uhf_ao_b ...'
+  call wall_time(ti)
+
+  fock_3e_uhf_ao_b = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                                                &
+ !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp,  &
+ !$OMP          i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) &
+ !$OMP SHARED  (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_b)
+
+  allocate(f_tmp(ao_num,ao_num))
+  f_tmp = 0.d0
+
+ !$OMP DO
+  do g = 1, ao_num
+    do e = 1, ao_num
+      dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e)
+      dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e)
+      dm_ge   = dm_ge_a + dm_ge_b
+      do d = 1, ao_num
+        do k = 1, ao_num
+          dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k)
+          dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k)
+          dm_dk   = dm_dk_a + dm_dk_b
+          do mu = 1, ao_num
+            do nu = 1, ao_num
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk)
+              call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu)
+              f_tmp(mu,nu) -= 0.5d0 * ( dm_ge   * dm_dk   * i_mugd_nuek &
+                                      + dm_ge_b * dm_dk_b * i_mugd_eknu &
+                                      + dm_ge_b * dm_dk_b * i_mugd_knue &
+                                      - dm_ge_b * dm_dk   * i_mugd_enuk &
+                                      - dm_ge   * dm_dk_b * i_mugd_kenu &
+                                      - dm_ge_b * dm_dk_b * i_mugd_nuke &
+                                      - dm_ge_a * dm_dk_a * i_mugd_nuke )
+            enddo
+          enddo
+        enddo
+      enddo
+    enddo
+  enddo
+ !$OMP END DO NOWAIT
+
+ !$OMP CRITICAL
+  do mu = 1, ao_num
+    do nu = 1, ao_num
+      fock_3e_uhf_ao_b(mu,nu) += f_tmp(mu,nu)
+    enddo
+  enddo
+ !$OMP END CRITICAL
+
+  deallocate(f_tmp)
+ !$OMP END PARALLEL
+
+  call wall_time(tf)
+  print *, ' total Wall time for fock_3e_uhf_ao_b =', tf - ti
+
+END_PROVIDER 
+
+! ---
+

src/tc_scf/fock_for_right.irp.f

@@ -0,0 +1,107 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, good_hermit_tc_fock_mat, (mo_num, mo_num)]
+
+  BEGIN_DOC
+! good_hermit_tc_fock_mat = Hermitian Upper triangular Fock matrix 
+!
+! The converged eigenvectors of such matrix yield to orthonormal vectors satisfying the left Brillouin theorem
+  END_DOC
+  implicit none
+  integer :: i, j
+
+  good_hermit_tc_fock_mat = Fock_matrix_tc_mo_tot
+  do j = 1, mo_num
+    do i = 1, j-1
+      good_hermit_tc_fock_mat(i,j) = Fock_matrix_tc_mo_tot(j,i) 
+    enddo
+  enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, hermit_average_tc_fock_mat, (mo_num, mo_num)]
+
+  BEGIN_DOC
+! hermit_average_tc_fock_mat = (F + F^\dagger)/2
+  END_DOC
+  implicit none
+  integer :: i, j
+
+  hermit_average_tc_fock_mat = Fock_matrix_tc_mo_tot
+  do j = 1, mo_num
+    do i = 1, mo_num
+      hermit_average_tc_fock_mat(i,j) = 0.5d0 * (Fock_matrix_tc_mo_tot(j,i) + Fock_matrix_tc_mo_tot(i,j))
+    enddo
+  enddo
+
+END_PROVIDER 
+
+
+! ---
+BEGIN_PROVIDER [ double precision, grad_hermit]
+ implicit none
+ BEGIN_DOC
+ ! square of gradient of the energy
+ END_DOC
+ if(symetric_fock_tc)then
+  grad_hermit = grad_hermit_average_tc_fock_mat
+ else
+  grad_hermit = grad_good_hermit_tc_fock_mat
+ endif
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, grad_good_hermit_tc_fock_mat]
+  implicit none
+  BEGIN_DOC
+  ! grad_good_hermit_tc_fock_mat = norm of gradients of the upper triangular TC fock
+  END_DOC
+  integer :: i, j
+  grad_good_hermit_tc_fock_mat = 0.d0
+  do i = 1, elec_alpha_num
+    do j = elec_alpha_num+1, mo_num
+      grad_good_hermit_tc_fock_mat += dabs(good_hermit_tc_fock_mat(i,j))
+    enddo
+  enddo
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, grad_hermit_average_tc_fock_mat]
+  implicit none
+  BEGIN_DOC
+  ! grad_hermit_average_tc_fock_mat = norm of gradients of the upper triangular TC fock
+  END_DOC
+  integer :: i, j
+  grad_hermit_average_tc_fock_mat = 0.d0
+  do i = 1, elec_alpha_num
+    do j = elec_alpha_num+1, mo_num
+      grad_hermit_average_tc_fock_mat += dabs(hermit_average_tc_fock_mat(i,j))
+    enddo
+  enddo
+END_PROVIDER 
+
+
+! ---
+
+subroutine save_good_hermit_tc_eigvectors()
+
+  implicit none
+  integer        :: sign
+  character*(64) :: label
+  logical        :: output
+
+  sign = 1
+  label = "Canonical"
+  output = .False.
+  
+  if(symetric_fock_tc)then
+   call mo_as_eigvectors_of_mo_matrix(hermit_average_tc_fock_mat, mo_num, mo_num, label, sign, output)
+  else
+   call mo_as_eigvectors_of_mo_matrix(good_hermit_tc_fock_mat, mo_num, mo_num, label, sign, output)
+  endif
+end subroutine save_good_hermit_tc_eigvectors
+
+! ---
+

src/tc_scf/fock_tc.irp.f

@@ -0,0 +1,307 @@
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, two_e_tc_non_hermit_integral_seq_alpha, (ao_num, ao_num)]
+&BEGIN_PROVIDER [ double precision, two_e_tc_non_hermit_integral_seq_beta , (ao_num, ao_num)]
+
+  BEGIN_DOC
+  !
+  ! two_e_tc_non_hermit_integral_seq_alpha(k,i) = <k| F^tc_alpha |i> 
+  !
+  ! where F^tc is the two-body part of the TC Fock matrix and k,i are AO basis functions
+  !
+  END_DOC
+
+  implicit none
+  integer          :: i, j, k, l
+  double precision :: density, density_a, density_b
+  double precision :: t0, t1
+
+  !print*, ' providing two_e_tc_non_hermit_integral_seq ...'
+  !call wall_time(t0)
+
+  two_e_tc_non_hermit_integral_seq_alpha = 0.d0
+  two_e_tc_non_hermit_integral_seq_beta  = 0.d0
+
+  do i = 1, ao_num
+    do k = 1, ao_num
+      do j = 1, ao_num
+        do l = 1, ao_num
+
+          density_a = TCSCF_density_matrix_ao_alpha(l,j)
+          density_b = TCSCF_density_matrix_ao_beta (l,j)
+          density   = density_a + density_b
+
+          !!                                         rho(l,j)   *      < k l| T | i j>
+          !two_e_tc_non_hermit_integral_seq_alpha(k,i) += density   * ao_two_e_tc_tot(l,j,k,i)
+          !!                                         rho(l,j)   *      < k l| T | i j>
+          !two_e_tc_non_hermit_integral_seq_beta (k,i) += density   * ao_two_e_tc_tot(l,j,k,i)
+          !!                                         rho_a(l,j) *      < l k| T | i j>
+          !two_e_tc_non_hermit_integral_seq_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i)
+          !!                                         rho_b(l,j) *      < l k| T | i j>
+          !two_e_tc_non_hermit_integral_seq_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i)
+
+          !!                                         rho(l,j)   *      < k l| T | i j>
+          !two_e_tc_non_hermit_integral_alpha(k,i) += density   * ao_two_e_tc_tot(l,j,k,i)
+          !!                                         rho(l,j)   *      < k l| T | i j>
+          !two_e_tc_non_hermit_integral_beta (k,i) += density   * ao_two_e_tc_tot(l,j,k,i)
+          !!                                         rho_a(l,j) *      < l k| T | i j>
+          !two_e_tc_non_hermit_integral_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i)
+          !!                                         rho_b(l,j) *      < l k| T | i j>
+          !two_e_tc_non_hermit_integral_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i)
+
+          !                                         rho(l,j)   *      < k l| T | i j>
+          two_e_tc_non_hermit_integral_seq_alpha(k,i) += density   * ao_two_e_tc_tot(k,i,l,j)
+          !                                         rho(l,j)   *      < k l| T | i j>
+          two_e_tc_non_hermit_integral_seq_beta (k,i) += density   * ao_two_e_tc_tot(k,i,l,j)
+          !                                         rho_a(l,j) *      < k l| T | j i>
+          two_e_tc_non_hermit_integral_seq_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i)
+          !                                         rho_b(l,j) *      < k l| T | j i>
+          two_e_tc_non_hermit_integral_seq_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i)
+
+        enddo
+      enddo
+    enddo
+  enddo
+
+  !call wall_time(t1)
+  !print*, ' wall time for two_e_tc_non_hermit_integral_seq after = ', t1 - t0
+
+END_PROVIDER 
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, two_e_tc_non_hermit_integral_alpha, (ao_num, ao_num)]
+&BEGIN_PROVIDER [ double precision, two_e_tc_non_hermit_integral_beta , (ao_num, ao_num)]
+
+  BEGIN_DOC
+  !
+  ! two_e_tc_non_hermit_integral_alpha(k,i) = <k| F^tc_alpha |i> 
+  !
+  ! where F^tc is the two-body part of the TC Fock matrix and k,i are AO basis functions
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, k, l
+  double precision              :: density, density_a, density_b, I_coul, I_kjli
+  double precision              :: t0, t1
+  double precision, allocatable :: tmp_a(:,:), tmp_b(:,:)
+
+  !print*, ' providing two_e_tc_non_hermit_integral ...'
+  !call wall_time(t0)
+
+  two_e_tc_non_hermit_integral_alpha = 0.d0
+  two_e_tc_non_hermit_integral_beta  = 0.d0
+
+ !$OMP PARALLEL DEFAULT (NONE)                                                                        &
+ !$OMP PRIVATE (i, j, k, l, density_a, density_b, density, tmp_a, tmp_b, I_coul, I_kjli)              &
+ !$OMP SHARED  (ao_num, TCSCF_density_matrix_ao_alpha, TCSCF_density_matrix_ao_beta, ao_two_e_tc_tot, &
+ !$OMP         two_e_tc_non_hermit_integral_alpha, two_e_tc_non_hermit_integral_beta)
+
+  allocate(tmp_a(ao_num,ao_num), tmp_b(ao_num,ao_num))
+  tmp_a = 0.d0
+  tmp_b = 0.d0
+
+ !$OMP DO
+  do j = 1, ao_num
+    do l = 1, ao_num
+      density_a = TCSCF_density_matrix_ao_alpha(l,j)
+      density_b = TCSCF_density_matrix_ao_beta (l,j)
+      density   = density_a + density_b                      
+      do i = 1, ao_num
+        do k = 1, ao_num
+
+          I_coul = density * ao_two_e_tc_tot(k,i,l,j)
+          I_kjli = ao_two_e_tc_tot(k,j,l,i)
+
+          tmp_a(k,i) += I_coul - density_a * I_kjli
+          tmp_b(k,i) += I_coul - density_b * I_kjli
+        enddo
+      enddo
+    enddo
+  enddo
+ !$OMP END DO NOWAIT
+
+ !$OMP CRITICAL
+  do i = 1, ao_num
+    do j = 1, ao_num
+      two_e_tc_non_hermit_integral_alpha(j,i) += tmp_a(j,i)
+      two_e_tc_non_hermit_integral_beta (j,i) += tmp_b(j,i)
+    enddo
+  enddo
+ !$OMP END CRITICAL
+
+  deallocate(tmp_a, tmp_b)
+ !$OMP END PARALLEL
+
+  !call wall_time(t1)
+  !print*, ' wall time for two_e_tc_non_hermit_integral after = ', t1 - t0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_alpha, (ao_num, ao_num)]
+
+  BEGIN_DOC
+  ! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the AO basis
+  END_DOC
+
+  implicit none
+
+  Fock_matrix_tc_ao_alpha =  ao_one_e_integrals_tc_tot + two_e_tc_non_hermit_integral_alpha 
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_beta, (ao_num, ao_num)]
+
+  BEGIN_DOC
+  ! Total beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis
+  END_DOC
+
+  implicit none
+
+  Fock_matrix_tc_ao_beta = ao_one_e_integrals_tc_tot + two_e_tc_non_hermit_integral_beta 
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
+
+  BEGIN_DOC
+  ! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis
+  END_DOC
+
+  implicit none
+  double precision, allocatable :: tmp(:,:)
+
+  if(bi_ortho) then
+
+    !allocate(tmp(ao_num,ao_num))
+    !tmp = Fock_matrix_tc_ao_alpha
+    !if(three_body_h_tc) then
+    !  tmp += fock_3e_uhf_ao_a
+    !endif
+    !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1))
+    !deallocate(tmp)
+
+    call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
+                          , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
+    if(three_body_h_tc) then
+      !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth
+      Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a
+    endif
+
+  else
+    call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
+                 , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
+
+  endif
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
+
+  BEGIN_DOC
+  ! Total beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis
+  END_DOC
+
+  implicit none
+  double precision, allocatable :: tmp(:,:)
+
+  if(bi_ortho) then
+
+    !allocate(tmp(ao_num,ao_num))
+    !tmp = Fock_matrix_tc_ao_beta
+    !if(three_body_h_tc) then
+    !  tmp += fock_3e_uhf_ao_b
+    !endif
+    !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1))
+    !deallocate(tmp)
+
+    call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
+                          , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
+    if(three_body_h_tc) then
+      !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth
+      Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b
+    endif
+
+  else
+
+    call ao_to_mo( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
+                 , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
+
+  endif
+
+END_PROVIDER
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, grad_non_hermit_left]
+&BEGIN_PROVIDER [ double precision, grad_non_hermit_right]
+&BEGIN_PROVIDER [ double precision, grad_non_hermit]
+
+  implicit none
+  integer :: i, k
+
+  grad_non_hermit_left  = 0.d0
+  grad_non_hermit_right = 0.d0
+
+  do i = 1, elec_beta_num ! doc --> SOMO
+    do k = elec_beta_num+1, elec_alpha_num
+      grad_non_hermit_left  = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i)))
+      grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k)))
+      !grad_non_hermit_left  += dabs(Fock_matrix_tc_mo_tot(k,i))
+      !grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k))
+      !grad_non_hermit_left  += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i)
+      !grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k)
+    enddo
+  enddo
+
+  do i = 1, elec_beta_num ! doc --> virt 
+    do k = elec_alpha_num+1, mo_num
+      grad_non_hermit_left  = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i)))
+      grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k)))
+      !grad_non_hermit_left  += dabs(Fock_matrix_tc_mo_tot(k,i))
+      !grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k))
+      grad_non_hermit_left  += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i)
+      grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k)
+    enddo
+  enddo
+
+  do i = elec_beta_num+1, elec_alpha_num ! SOMO --> virt 
+    do k = elec_alpha_num+1, mo_num
+      grad_non_hermit_left  = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i)))
+      grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k)))
+      !grad_non_hermit_left  += dabs(Fock_matrix_tc_mo_tot(k,i))
+      !grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k))
+      grad_non_hermit_left  += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i)
+      grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k)
+    enddo
+  enddo
+
+  !grad_non_hermit = dsqrt(grad_non_hermit_left) + dsqrt(grad_non_hermit_right)
+  grad_non_hermit = grad_non_hermit_left + grad_non_hermit_right
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
+
+  implicit none
+
+  call mo_to_ao_bi_ortho( Fock_matrix_tc_mo_tot, size(Fock_matrix_tc_mo_tot, 1) &
+                        , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) )
+
+END_PROVIDER
+
+! ---
+
+

src/tc_scf/fock_tc_mo_tot.irp.f

@@ -0,0 +1,144 @@
+
+ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num,mo_num) ]
+&BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (mo_num)]
+   implicit none
+   BEGIN_DOC
+   ! Fock matrix on the MO basis.
+   ! For open shells, the ROHF Fock Matrix is ::
+   !
+   !       |   F-K    |  F + K/2  |    F     |
+   !       |---------------------------------|
+   !       | F + K/2  |     F     |  F - K/2 |
+   !       |---------------------------------|
+   !       |    F     |  F - K/2  |  F + K   |
+   !
+   !
+   ! F = 1/2 (Fa + Fb)
+   !
+   ! K = Fb - Fa
+   !
+   END_DOC
+   integer                        :: i,j,n
+   if (elec_alpha_num == elec_beta_num) then
+     Fock_matrix_tc_mo_tot = Fock_matrix_tc_mo_alpha
+   else
+
+     do j=1,elec_beta_num
+       ! F-K
+       do i=1,elec_beta_num !CC
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             - (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+       ! F+K/2
+       do i=elec_beta_num+1,elec_alpha_num  !CA
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+       ! F
+       do i=elec_alpha_num+1, mo_num !CV
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
+       enddo
+     enddo
+
+     do j=elec_beta_num+1,elec_alpha_num
+       ! F+K/2
+       do i=1,elec_beta_num !AC
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+       ! F
+       do i=elec_beta_num+1,elec_alpha_num !AA
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
+       enddo
+       ! F-K/2
+       do i=elec_alpha_num+1, mo_num !AV
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+     enddo
+
+     do j=elec_alpha_num+1, mo_num
+       ! F
+       do i=1,elec_beta_num !VC
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
+       enddo
+       ! F-K/2
+       do i=elec_beta_num+1,elec_alpha_num !VA
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+       ! F+K
+       do i=elec_alpha_num+1,mo_num !VV
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) &
+             + (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+     enddo
+     if(three_body_h_tc)then
+      ! C-O
+      do j = 1, elec_beta_num
+       do i = elec_beta_num+1, elec_alpha_num
+        Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
+        Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
+       enddo
+      enddo
+      ! C-V
+      do j = 1, elec_beta_num
+       do i = elec_alpha_num+1, mo_num
+        Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
+        Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
+       enddo
+      enddo
+      ! O-V
+      do j = elec_beta_num+1, elec_alpha_num
+       do i = elec_alpha_num+1, mo_num
+        Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
+        Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
+       enddo
+      enddo
+     endif
+
+   endif
+
+   do i = 1, mo_num
+     Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i)
+   enddo
+
+
+   if(frozen_orb_scf)then
+     integer                        :: iorb,jorb
+     do i = 1, n_core_orb
+      iorb = list_core(i)
+      do j = 1, n_act_orb
+       jorb = list_act(j)
+       Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
+       Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
+      enddo
+     enddo
+   endif
+
+   if(no_oa_or_av_opt)then
+     do i = 1, n_act_orb
+       iorb = list_act(i)
+       do j = 1, n_inact_orb
+         jorb = list_inact(j)
+         Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
+         Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
+       enddo
+       do j = 1, n_virt_orb
+         jorb = list_virt(j)
+         Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
+         Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
+       enddo
+       do j = 1, n_core_orb
+         jorb = list_core(j)
+         Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
+         Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0                                                                                                                 
+       enddo
+     enddo
+   endif
+  if(.not.bi_ortho .and. three_body_h_tc)then
+   Fock_matrix_tc_mo_tot += fock_3_mat
+  endif
+
+END_PROVIDER
+

src/tc_scf/fock_three.irp.f

@@ -0,0 +1,229 @@
+BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)] 
+ implicit none
+  integer :: i,j
+  double precision :: contrib
+  fock_3_mat = 0.d0
+  if(.not.bi_ortho.and.three_body_h_tc)then
+   call give_fock_ia_three_e_total(1,1,contrib)
+!!  !$OMP PARALLEL                  &
+!!  !$OMP DEFAULT (NONE)            &
+!!  !$OMP PRIVATE (i,j,m,integral) & 
+!!  !$OMP SHARED (mo_num,three_body_3_index)
+!!  !$OMP DO SCHEDULE (guided) COLLAPSE(3)
+   do i = 1, mo_num
+    do j = 1, mo_num
+     call give_fock_ia_three_e_total(j,i,contrib)
+     fock_3_mat(j,i) = -contrib
+    enddo
+   enddo
+  else if(bi_ortho.and.three_body_h_tc)then
+!!  !$OMP END DO
+!!  !$OMP END PARALLEL
+!!  do i = 1, mo_num
+!!   do j = 1, i-1
+!!    mat_three(j,i) = mat_three(i,j)
+!!   enddo
+!!  enddo
+ endif
+
+END_PROVIDER 
+
+
+subroutine give_fock_ia_three_e_total(i,a,contrib)
+ implicit none
+ BEGIN_DOC
+! contrib is the TOTAL (same spins / opposite spins) contribution from the three body term to the Fock operator 
+!
+ END_DOC
+ integer, intent(in) :: i,a
+ double precision, intent(out) :: contrib
+ double precision :: int_1, int_2, int_3
+ double precision :: mos_i, mos_a, w_ia
+ double precision :: mos_ia, weight
+
+ integer :: mm, ipoint,k,l
+
+ int_1 = 0.d0
+ int_2 = 0.d0
+ int_3 = 0.d0
+ do mm = 1, 3
+  do ipoint = 1, n_points_final_grid
+   weight = final_weight_at_r_vector(ipoint)                                                                          
+   mos_i  = mos_in_r_array_transp(ipoint,i) 
+   mos_a  = mos_in_r_array_transp(ipoint,a) 
+   mos_ia = mos_a * mos_i
+   w_ia   = x_W_ij_erf_rk(ipoint,mm,i,a) 
+     
+   int_1  += weight * fock_3_w_kk_sum(ipoint,mm) * (4.d0 * fock_3_rho_beta(ipoint) * w_ia               & 
+                                                  + 2.0d0 * mos_ia * fock_3_w_kk_sum(ipoint,mm)         & 
+                                                  - 2.0d0 * fock_3_w_ki_mos_k(ipoint,mm,i) * mos_a      & 
+                                                  - 2.0d0 * fock_3_w_ki_mos_k(ipoint,mm,a) * mos_i      )
+   int_2  += weight * (-1.d0) * ( 2.0d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * w_ia                     & 
+                                + 2.0d0 * fock_3_rho_beta(ipoint) * fock_3_w_ki_wk_a(ipoint,mm,i,a)   & 
+                                + 1.0d0 * mos_ia * fock_3_trace_w_tilde(ipoint,mm)                    )
+
+   int_3  += weight *   1.d0  * (fock_3_w_kl_wla_phi_k(ipoint,mm,i) * mos_a + fock_3_w_kl_wla_phi_k(ipoint,mm,a) * mos_i & 
+                                +fock_3_w_ki_mos_k(ipoint,mm,i)     * fock_3_w_ki_mos_k(ipoint,mm,a)                     )
+  enddo
+ enddo
+ contrib = int_1 + int_2 + int_3
+
+end
+
+! ---
+
+BEGIN_PROVIDER [double precision, diag_three_elem_hf]
+
+  implicit none
+  integer          :: i, j, k, ipoint, mm
+  double precision :: contrib, weight, four_third, one_third, two_third, exchange_int_231
+  double precision :: integral_aaa, hthree, integral_aab, integral_abb, integral_bbb
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  !print *, ' providing diag_three_elem_hf'
+
+  if(.not. three_body_h_tc) then
+
+    diag_three_elem_hf = 0.d0
+
+  else
+
+    if(.not. bi_ortho) then
+
+      ! ---
+
+      one_third  = 1.d0/3.d0
+      two_third  = 2.d0/3.d0
+      four_third = 4.d0/3.d0
+      diag_three_elem_hf = 0.d0
+      do i = 1, elec_beta_num
+        do j = 1, elec_beta_num
+          do k = 1, elec_beta_num
+            call give_integrals_3_body(k, j, i, j, i, k,exchange_int_231)   
+            diag_three_elem_hf += two_third * exchange_int_231
+          enddo
+        enddo
+      enddo
+      do mm = 1, 3
+        do ipoint = 1, n_points_final_grid
+          weight  = final_weight_at_r_vector(ipoint)                                                                          
+          contrib = 3.d0 * fock_3_w_kk_sum(ipoint,mm) * fock_3_rho_beta(ipoint) * fock_3_w_kk_sum(ipoint,mm) & 
+                  - 2.d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * fock_3_w_kk_sum(ipoint,mm)                     & 
+                  - 1.d0 * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm)
+          contrib *= four_third
+          contrib += -two_third  * fock_3_rho_beta(ipoint)    * fock_3_w_kl_w_kl(ipoint,mm) & 
+                     -four_third * fock_3_w_kk_sum(ipoint,mm) * fock_3_w_kl_mo_k_mo_l(ipoint,mm)
+          diag_three_elem_hf += weight * contrib
+       enddo
+      enddo
+
+      diag_three_elem_hf = - diag_three_elem_hf
+
+      ! ---
+
+    else
+
+      provide mo_l_coef mo_r_coef
+      call give_aaa_contrib(integral_aaa)
+      call give_aab_contrib(integral_aab)
+      call give_abb_contrib(integral_abb)
+      call give_bbb_contrib(integral_bbb)
+      diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb
+!      print*,'integral_aaa + integral_aab + integral_abb + integral_bbb'
+!      print*,integral_aaa , integral_aab , integral_abb , integral_bbb
+
+    endif
+
+  endif
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, fock_3_mat_a_op_sh, (mo_num, mo_num)]
+ implicit none 
+ integer :: h,p,i,j
+ double precision :: direct_int, exch_int, exchange_int_231, exchange_int_312
+ double precision :: exchange_int_23, exchange_int_12, exchange_int_13 
+
+ fock_3_mat_a_op_sh = 0.d0
+ do h = 1, mo_num
+  do p = 1, mo_num
+   !F_a^{ab}(h,p) 
+   do i = 1, elec_beta_num ! beta 
+    do j = elec_beta_num+1, elec_alpha_num ! alpha
+     call  give_integrals_3_body(h,j,i,p,j,i,direct_int)    ! <hji|pji>
+     call  give_integrals_3_body(h,j,i,j,p,i,exch_int)   
+     fock_3_mat_a_op_sh(h,p) -= direct_int - exch_int
+    enddo
+   enddo
+   !F_a^{aa}(h,p)
+   do i = 1, elec_beta_num ! alpha 
+    do j = elec_beta_num+1, elec_alpha_num ! alpha
+       direct_int = three_body_4_index(j,i,h,p)                    
+       call  give_integrals_3_body(h,j,i,p,j,i,direct_int) 
+       call  give_integrals_3_body(h,j,i,i,p,j,exchange_int_231)
+       call  give_integrals_3_body(h,j,i,j,i,p,exchange_int_312) 
+       call  give_integrals_3_body(h,j,i,p,i,j,exchange_int_23) 
+       call  give_integrals_3_body(h,j,i,i,j,p,exchange_int_12)
+       call  give_integrals_3_body(h,j,i,j,p,i,exchange_int_13)  
+       fock_3_mat_a_op_sh(h,p) -= ( direct_int + exchange_int_231 + exchange_int_312 & 
+              -  exchange_int_23 & ! i <-> j
+              -  exchange_int_12 & ! p <-> j
+              -  exchange_int_13  )! p <-> i
+    enddo 
+   enddo
+  enddo
+ enddo
+! symmetrized 
+! do p = 1, elec_beta_num
+!  do h = elec_alpha_num +1, mo_num
+!   fock_3_mat_a_op_sh(h,p) = fock_3_mat_a_op_sh(p,h)
+!  enddo
+! enddo
+ 
+! do h = elec_beta_num+1, elec_alpha_num
+!  do p = elec_alpha_num +1, mo_num
+!   !F_a^{bb}(h,p) 
+!   do i = 1, elec_beta_num
+!    do j = i+1, elec_beta_num
+!     call  give_integrals_3_body(h,j,i,p,j,i,direct_int)   
+!     call  give_integrals_3_body(h,j,i,p,i,j,exch_int)   
+!     fock_3_mat_a_op_sh(h,p) -= direct_int - exch_int
+!    enddo
+!   enddo
+!  enddo
+! enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_3_mat_b_op_sh, (mo_num, mo_num)]
+ implicit none 
+ integer :: h,p,i,j
+ double precision :: direct_int, exch_int
+ fock_3_mat_b_op_sh = 0.d0
+ do h = 1, elec_beta_num
+  do p = elec_alpha_num +1, mo_num
+   !F_b^{aa}(h,p) 
+   do i = 1, elec_beta_num
+    do j = elec_beta_num+1, elec_alpha_num
+     call  give_integrals_3_body(h,j,i,p,j,i,direct_int)   
+     call  give_integrals_3_body(h,j,i,p,i,j,exch_int)   
+     fock_3_mat_b_op_sh(h,p) += direct_int - exch_int
+    enddo
+   enddo
+
+   !F_b^{ab}(h,p) 
+   do i = elec_beta_num+1, elec_beta_num
+    do j = 1, elec_beta_num
+     call  give_integrals_3_body(h,j,i,p,j,i,direct_int)   
+     call  give_integrals_3_body(h,j,i,j,p,i,exch_int)   
+     fock_3_mat_b_op_sh(h,p) += direct_int - exch_int
+    enddo
+   enddo
+ 
+  enddo
+ enddo
+
+END_PROVIDER 

src/tc_scf/fock_three_bi_ortho.irp.f

@@ -0,0 +1,178 @@
+BEGIN_PROVIDER [ double precision, fock_a_abb_3e_bi_orth_old, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_a_abb_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,beta,beta contribution
+ END_DOC
+ fock_a_abb_3e_bi_orth_old = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_23_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_beta_num
+    do k = j+1, elec_beta_num
+      ! see contrib_3e_soo
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23
+      fock_a_abb_3e_bi_orth_old(a,i) += direct_int - exch_23_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_a_abb_3e_bi_orth_old = - fock_a_abb_3e_bi_orth_old
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_a_aba_3e_bi_orth_old, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_a_aba_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,beta contribution
+ END_DOC
+ fock_a_aba_3e_bi_orth_old = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_13_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_alpha_num ! a
+    do k = 1, elec_beta_num ! b
+                                                                       !   a b a   a b a
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 
+      fock_a_aba_3e_bi_orth_old(a,i) += direct_int - exch_13_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_a_aba_3e_bi_orth_old = - fock_a_aba_3e_bi_orth_old
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_a_aaa_3e_bi_orth_old, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_a_aaa_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution
+ END_DOC
+ fock_a_aaa_3e_bi_orth_old = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_alpha_num
+    do k = j+1, elec_alpha_num
+      ! positive terms :: cycle contrib 
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int)      ! < a k j | j i k >
+      call  give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i >
+      fock_a_aaa_3e_bi_orth_old(a,i) += direct_int + c_3_int + c_minus_3_int 
+      ! negative terms :: exchange contrib
+      call  give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 
+      call  give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
+      call  give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12
+      fock_a_aaa_3e_bi_orth_old(a,i) += - exch_13_int - exch_23_int  - exch_12_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_a_aaa_3e_bi_orth_old = - fock_a_aaa_3e_bi_orth_old
+END_PROVIDER 
+
+BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth_old, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! fock_a_tot_3e_bi_orth_old = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions 
+ END_DOC
+ fock_a_tot_3e_bi_orth_old = fock_a_abb_3e_bi_orth_old + fock_a_aba_3e_bi_orth_old + fock_a_aaa_3e_bi_orth_old 
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_b_baa_3e_bi_orth_old, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_b_baa_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,alpha contribution
+ END_DOC
+ fock_b_baa_3e_bi_orth_old = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_23_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_alpha_num
+    do k = j+1, elec_alpha_num
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23
+      fock_b_baa_3e_bi_orth_old(a,i) += direct_int - exch_23_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_b_baa_3e_bi_orth_old = - fock_b_baa_3e_bi_orth_old
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_b_bab_3e_bi_orth_old, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_b_bab_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,beta contribution
+ END_DOC
+ fock_b_bab_3e_bi_orth_old = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_13_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_beta_num
+    do k = 1, elec_alpha_num
+      !                                                                    b a b   b a b
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13
+      fock_b_bab_3e_bi_orth_old(a,i) += direct_int - exch_13_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_b_bab_3e_bi_orth_old = - fock_b_bab_3e_bi_orth_old
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_b_bbb_3e_bi_orth_old, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_b_bbb_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution
+ END_DOC
+ fock_b_bbb_3e_bi_orth_old = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_beta_num
+    do k = j+1, elec_beta_num
+      ! positive terms :: cycle contrib 
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int)      ! < a k j | j i k >
+      call  give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i >
+      fock_b_bbb_3e_bi_orth_old(a,i) += direct_int + c_3_int + c_minus_3_int 
+      ! negative terms :: exchange contrib
+      call  give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 
+      call  give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
+      call  give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12
+      fock_b_bbb_3e_bi_orth_old(a,i) += - exch_13_int - exch_23_int  - exch_12_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_b_bbb_3e_bi_orth_old = - fock_b_bbb_3e_bi_orth_old
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_b_tot_3e_bi_orth_old, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! fock_b_tot_3e_bi_orth_old = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions 
+ END_DOC
+ fock_b_tot_3e_bi_orth_old = fock_b_bbb_3e_bi_orth_old + fock_b_bab_3e_bi_orth_old + fock_b_baa_3e_bi_orth_old
+
+END_PROVIDER 

src/tc_scf/fock_three_bi_ortho_new_new.irp.f

@@ -0,0 +1,286 @@
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)]
+
+  implicit none
+  integer :: i, a
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  fock_a_tot_3e_bi_orth = 0.d0
+
+  do i = 1, mo_num
+    do a = 1, mo_num
+      fock_a_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth  (a,i)
+      fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp1_bi_ortho(a,i)
+      fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp2_bi_ortho(a,i)
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_b_tot_3e_bi_orth, (mo_num, mo_num)]
+
+  implicit none
+  integer :: i, a
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  fock_b_tot_3e_bi_orth = 0.d0
+
+  do i = 1, mo_num
+    do a = 1, mo_num
+      fock_b_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth  (a,i)
+      fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp2_bi_ortho(a,i)
+      fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp1_bi_ortho(a,i)
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_cs_3e_bi_orth, (mo_num, mo_num)]
+
+  implicit none
+  integer          :: i, a, j, k
+  double precision :: contrib_sss, contrib_sos, contrib_soo, contrib
+  double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+  double precision :: new
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  fock_cs_3e_bi_orth = 0.d0
+
+  do i = 1, mo_num
+    do a = 1, mo_num
+    
+      do j = 1, elec_beta_num
+        do k = 1, elec_beta_num
+
+          !!call contrib_3e_sss(a,i,j,k,contrib_sss)
+          !!call contrib_3e_soo(a,i,j,k,contrib_soo)
+          !!call contrib_3e_sos(a,i,j,k,contrib_sos)
+          !!contrib = 0.5d0 * (contrib_sss + contrib_soo) + contrib_sos
+ 
+          call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
+          call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int)      ! < a k j | j i k >
+          call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i >
+
+          ! negative terms :: exchange contrib
+          call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 
+          call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
+          call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12
+
+          new = 2.d0 * direct_int + 0.5d0 * (c_3_int + c_minus_3_int - exch_12_int) -1.5d0 * exch_13_int - exch_23_int
+
+          fock_cs_3e_bi_orth(a,i) += new
+        enddo
+      enddo
+    enddo
+  enddo
+ 
+  fock_cs_3e_bi_orth = - fock_cs_3e_bi_orth
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_a_tmp1_bi_ortho, (mo_num, mo_num)]
+
+  implicit none
+  integer          :: i, a, j, k
+  double precision :: contrib_sss, contrib_sos, contrib_soo, contrib
+  double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+  double precision :: new
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  fock_a_tmp1_bi_ortho = 0.d0
+
+  do i = 1, mo_num
+    do a = 1, mo_num
+    
+      do j = elec_beta_num + 1, elec_alpha_num 
+        do k = 1, elec_beta_num
+          call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
+          call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int)      ! < a k j | j i k >
+          call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i >
+          call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 
+          call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
+          call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12
+          
+          fock_a_tmp1_bi_ortho(a,i) += 1.5d0 * (direct_int - exch_13_int) + 0.5d0 * (c_3_int + c_minus_3_int - exch_23_int - exch_12_int)
+        enddo
+      enddo
+    enddo
+  enddo
+
+  fock_a_tmp1_bi_ortho = - fock_a_tmp1_bi_ortho
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_a_tmp2_bi_ortho, (mo_num, mo_num)]
+
+  implicit none
+  integer          :: i, a, j, k
+  double precision :: contrib_sss
+
+  PROVIDE mo_l_coef mo_r_coef
+ 
+  fock_a_tmp2_bi_ortho = 0.d0
+
+  do i = 1, mo_num
+    do a = 1, mo_num
+      do j = 1, elec_alpha_num
+        do k = elec_beta_num+1, elec_alpha_num
+          call contrib_3e_sss(a, i, j, k, contrib_sss)
+
+          fock_a_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_sss
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_b_tmp1_bi_ortho, (mo_num, mo_num)]
+
+  implicit none
+  integer          :: i, a, j, k
+  double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int
+  double precision :: new
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  fock_b_tmp1_bi_ortho = 0.d0
+
+  do i = 1, mo_num
+    do a = 1, mo_num
+      do j = 1, elec_beta_num
+        do k = elec_beta_num+1, elec_alpha_num
+          call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
+          call  give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 
+          call  give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
+
+          fock_b_tmp1_bi_ortho(a,i) += 1.5d0 * direct_int - 0.5d0 * exch_23_int - exch_13_int
+        enddo
+      enddo
+    enddo
+  enddo
+
+  fock_b_tmp1_bi_ortho = - fock_b_tmp1_bi_ortho
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, fock_b_tmp2_bi_ortho, (mo_num, mo_num)]
+
+  implicit none
+  integer          :: i, a, j, k
+  double precision :: contrib_soo
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  fock_b_tmp2_bi_ortho = 0.d0
+
+  do i = 1, mo_num
+    do a = 1, mo_num
+      do j = elec_beta_num + 1, elec_alpha_num 
+        do k = 1, elec_alpha_num
+          call contrib_3e_soo(a, i, j, k, contrib_soo)
+
+          fock_b_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_soo
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+subroutine contrib_3e_sss(a, i, j, k, integral)
+
+  BEGIN_DOC
+  ! returns the pure same spin contribution to F(a,i) from two orbitals j,k
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: a, i, j, k
+  double precision, intent(out) :: integral
+  double precision              :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
+  call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int)      ! < a k j | j i k >
+  call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i >
+  integral = direct_int + c_3_int + c_minus_3_int 
+
+  ! negative terms :: exchange contrib
+  call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 
+  call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
+  call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12
+  integral += - exch_13_int - exch_23_int  - exch_12_int 
+
+  integral = -integral
+
+end
+
+! ---
+
+subroutine contrib_3e_soo(a,i,j,k,integral)
+
+  BEGIN_DOC
+  ! returns the same spin / opposite spin / opposite spin contribution to F(a,i) from two orbitals j,k
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: a, i, j, k
+  double precision, intent(out) :: integral
+  double precision              :: direct_int, exch_23_int
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
+  call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23
+  integral = direct_int - exch_23_int 
+
+  integral = -integral
+
+end
+
+! ---
+
+subroutine contrib_3e_sos(a, i, j, k, integral)
+
+  BEGIN_DOC
+  ! returns the same spin / opposite spin / same spin contribution to F(a,i) from two orbitals j,k
+  END_DOC
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  implicit none
+  integer,          intent(in)  :: a, i, j, k
+  double precision, intent(out) :: integral
+  double precision              :: direct_int, exch_13_int
+
+  call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j >
+  call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 
+  integral = direct_int - exch_13_int 
+
+  integral = -integral
+
+end
+
+! ---
+

src/tc_scf/fock_three_utils.irp.f

@@ -0,0 +1,140 @@
+
+BEGIN_PROVIDER [ double precision, fock_3_w_kk_sum, (n_points_final_grid,3)]
+ implicit none
+ integer :: mm, ipoint,k
+ double precision :: w_kk
+ fock_3_w_kk_sum = 0.d0
+ do k = 1, elec_beta_num
+  do mm = 1, 3
+   do ipoint = 1, n_points_final_grid
+    w_kk   = x_W_ij_erf_rk(ipoint,mm,k,k) 
+    fock_3_w_kk_sum(ipoint,mm) += w_kk
+   enddo
+  enddo
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_3_w_ki_mos_k, (n_points_final_grid,3,mo_num)]
+ implicit none
+ integer :: mm, ipoint,k,i
+ double precision :: w_ki, mo_k
+ fock_3_w_ki_mos_k = 0.d0
+ do i = 1, mo_num
+  do k = 1, elec_beta_num
+   do mm = 1, 3
+    do ipoint = 1, n_points_final_grid
+     w_ki   = x_W_ij_erf_rk(ipoint,mm,k,i) 
+     mo_k = mos_in_r_array(k,ipoint)
+     fock_3_w_ki_mos_k(ipoint,mm,i) += w_ki * mo_k
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_3_w_kl_w_kl, (n_points_final_grid,3)]
+ implicit none
+ integer :: k,j,ipoint,mm
+ double precision :: w_kj
+ fock_3_w_kl_w_kl = 0.d0
+ do j = 1, elec_beta_num
+  do k = 1, elec_beta_num
+   do mm = 1, 3
+    do ipoint = 1, n_points_final_grid
+     w_kj   = x_W_ij_erf_rk(ipoint,mm,k,j) 
+     fock_3_w_kl_w_kl(ipoint,mm) += w_kj * w_kj
+    enddo
+   enddo
+  enddo
+ enddo
+
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_3_rho_beta, (n_points_final_grid)]
+ implicit none
+ integer :: ipoint,k
+ fock_3_rho_beta = 0.d0
+ do ipoint = 1, n_points_final_grid
+  do k = 1, elec_beta_num
+   fock_3_rho_beta(ipoint) += mos_in_r_array(k,ipoint) * mos_in_r_array(k,ipoint)
+  enddo
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_3_w_kl_mo_k_mo_l, (n_points_final_grid,3)]
+ implicit none
+ integer :: ipoint,k,l,mm
+ double precision :: mos_k, mos_l, w_kl
+ fock_3_w_kl_mo_k_mo_l = 0.d0
+ do k = 1, elec_beta_num
+  do l = 1, elec_beta_num
+   do mm = 1, 3
+    do ipoint = 1, n_points_final_grid
+     mos_k  = mos_in_r_array_transp(ipoint,k) 
+     mos_l  = mos_in_r_array_transp(ipoint,l) 
+     w_kl   = x_W_ij_erf_rk(ipoint,mm,l,k)
+     fock_3_w_kl_mo_k_mo_l(ipoint,mm) += w_kl * mos_k * mos_l 
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_3_w_ki_wk_a, (n_points_final_grid,3,mo_num, mo_num)]
+ implicit none
+ integer :: ipoint,i,a,k,mm
+ double precision :: w_ki,w_ka
+ fock_3_w_ki_wk_a = 0.d0
+ do i = 1, mo_num
+  do a = 1, mo_num
+   do mm = 1, 3
+    do ipoint = 1, n_points_final_grid
+     do k = 1, elec_beta_num
+      w_ki   = x_W_ij_erf_rk(ipoint,mm,k,i)
+      w_ka   = x_W_ij_erf_rk(ipoint,mm,k,a)
+      fock_3_w_ki_wk_a(ipoint,mm,a,i) += w_ki * w_ka
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_3_trace_w_tilde, (n_points_final_grid,3)]
+ implicit none
+ integer :: ipoint,k,mm
+ fock_3_trace_w_tilde = 0.d0
+ do k = 1, elec_beta_num
+   do mm = 1, 3
+    do ipoint = 1, n_points_final_grid
+     fock_3_trace_w_tilde(ipoint,mm) += fock_3_w_ki_wk_a(ipoint,mm,k,k)
+    enddo
+   enddo
+ enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_3_w_kl_wla_phi_k, (n_points_final_grid,3,mo_num)]
+ implicit none
+ integer :: ipoint,a,k,mm,l
+ double precision :: w_kl,w_la, mo_k
+ fock_3_w_kl_wla_phi_k = 0.d0
+ do a = 1, mo_num
+  do k = 1, elec_beta_num 
+   do l = 1, elec_beta_num
+    do mm = 1, 3
+     do ipoint = 1, n_points_final_grid
+      w_kl   = x_W_ij_erf_rk(ipoint,mm,l,k)
+      w_la   = x_W_ij_erf_rk(ipoint,mm,l,a)
+      mo_k  = mos_in_r_array_transp(ipoint,k) 
+      fock_3_w_kl_wla_phi_k(ipoint,mm,a) += w_kl * w_la * mo_k
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+END_PROVIDER 
+

src/tc_scf/integrals_in_r_stuff.irp.f

@@ -0,0 +1,391 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, tc_scf_dm_in_r, (n_points_final_grid) ]
+
+  implicit none
+  integer :: i, j
+
+  tc_scf_dm_in_r = 0.d0
+  do i = 1, n_points_final_grid
+    do j = 1, elec_beta_num
+      tc_scf_dm_in_r(i) += mos_r_in_r_array(j,i) * mos_l_in_r_array(j,i)
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, w_sum_in_r, (n_points_final_grid, 3)]
+
+  implicit none
+  integer :: ipoint, j, xi
+
+  w_sum_in_r = 0.d0
+  do j = 1, elec_beta_num
+    do xi = 1, 3
+      do ipoint = 1, n_points_final_grid
+        !w_sum_in_r(ipoint,xi) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,j)
+        w_sum_in_r(ipoint,xi) += x_W_ki_bi_ortho_erf_rk_diag(ipoint,xi,j)
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, ww_sum_in_r, (n_points_final_grid, 3)]
+
+  implicit none
+  integer          :: ipoint, j, xi
+  double precision :: tmp
+
+  ww_sum_in_r = 0.d0
+  do j = 1, elec_beta_num
+    do xi = 1, 3
+      do ipoint = 1, n_points_final_grid
+        tmp = x_W_ki_bi_ortho_erf_rk_diag(ipoint,xi,j)
+        ww_sum_in_r(ipoint,xi) += tmp * tmp
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, W1_r_in_r, (n_points_final_grid, 3, mo_num)]
+
+  implicit none
+  integer :: i, j, xi, ipoint
+
+  ! TODO: call lapack
+
+  W1_r_in_r = 0.d0
+  do i = 1, mo_num
+    do j = 1, elec_beta_num
+      do xi = 1, 3
+        do ipoint = 1, n_points_final_grid
+          W1_r_in_r(ipoint,xi,i) += mos_r_in_r_array_transp(ipoint,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,i)
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, W1_l_in_r, (n_points_final_grid, 3, mo_num)]
+
+  implicit none
+  integer :: i, j, xi, ipoint
+
+  ! TODO: call lapack
+
+  W1_l_in_r = 0.d0
+  do i = 1, mo_num
+    do j = 1, elec_beta_num
+      do xi = 1, 3
+        do ipoint = 1, n_points_final_grid
+          W1_l_in_r(ipoint,xi,i) += mos_l_in_r_array_transp(ipoint,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,i,j)
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, W1_in_r, (n_points_final_grid, 3)]
+
+  implicit none
+  integer :: j, xi, ipoint
+
+  ! TODO: call lapack
+
+  W1_in_r = 0.d0
+  do j = 1, elec_beta_num
+    do xi = 1, 3
+      do ipoint = 1, n_points_final_grid
+        W1_in_r(ipoint,xi) += W1_l_in_r(ipoint,xi,j) * mos_r_in_r_array_transp(ipoint,j)
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, W1_diag_in_r, (n_points_final_grid, 3)]
+
+  implicit none
+  integer :: j, xi, ipoint
+
+  ! TODO: call lapack
+
+  W1_diag_in_r = 0.d0
+  do j = 1, elec_beta_num
+    do xi = 1, 3
+      do ipoint = 1, n_points_final_grid
+        W1_diag_in_r(ipoint,xi) += mos_r_in_r_array_transp(ipoint,j) * mos_l_in_r_array_transp(ipoint,j) * x_W_ki_bi_ortho_erf_rk_diag(ipoint,xi,j)
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, v_sum_in_r, (n_points_final_grid, 3)]
+
+  implicit none
+  integer :: i, j, xi, ipoint
+
+  ! TODO: call lapack
+  v_sum_in_r = 0.d0
+  do i = 1, elec_beta_num
+    do j = 1, elec_beta_num
+      do xi = 1, 3
+        do ipoint = 1, n_points_final_grid
+          v_sum_in_r(ipoint,xi) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,i,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,i)
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, W1_W1_r_in_r, (n_points_final_grid, 3, mo_num)]
+
+  implicit none
+  integer :: i, m, xi, ipoint
+
+  ! TODO: call lapack
+
+  W1_W1_r_in_r = 0.d0
+  do i = 1, mo_num
+    do m = 1, elec_beta_num
+      do xi = 1, 3
+        do ipoint = 1, n_points_final_grid
+          W1_W1_r_in_r(ipoint,xi,i) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,m,i) * W1_r_in_r(ipoint,xi,m)
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, W1_W1_l_in_r, (n_points_final_grid, 3, mo_num)]
+
+  implicit none
+  integer :: i, j, xi, ipoint
+
+  ! TODO: call lapack
+
+  W1_W1_l_in_r = 0.d0
+  do i = 1, mo_num
+    do j = 1, elec_beta_num
+      do xi = 1, 3
+        do ipoint = 1, n_points_final_grid
+          W1_W1_l_in_r(ipoint,xi,i) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,i,j) * W1_l_in_r(ipoint,xi,j)
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+! ---
+
+subroutine direct_term_imj_bi_ortho(a, i, integral)
+
+  BEGIN_DOC
+  ! computes sum_(j,m = 1, elec_beta_num) < a m j | i m j > with bi ortho mos
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: i, a
+  double precision, intent(out) :: integral
+
+  integer                       :: ipoint, xi
+  double precision              :: weight, tmp
+
+  integral = 0.d0
+  do xi = 1, 3
+    do ipoint = 1, n_points_final_grid
+      weight = final_weight_at_r_vector(ipoint)
+      !integral += ( mos_l_in_r_array(a,ipoint) * mos_r_in_r_array(i,ipoint) * w_sum_in_r(ipoint,xi) * w_sum_in_r(ipoint,xi) & 
+      !            + 2.d0 * tc_scf_dm_in_r(ipoint) * w_sum_in_r(ipoint,xi) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) ) * weight
+
+      tmp = w_sum_in_r(ipoint,xi)
+
+      integral += ( mos_l_in_r_array_transp(ipoint,a) * mos_r_in_r_array_transp(ipoint,i) * tmp * tmp & 
+                  + 2.d0 * tc_scf_dm_in_r(ipoint) * tmp * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i)       &
+                  ) * weight
+    enddo
+  enddo
+
+end 
+
+! ---
+
+subroutine exch_term_jmi_bi_ortho(a, i, integral)
+
+  BEGIN_DOC
+  ! computes sum_(j,m = 1, elec_beta_num) < a m j | j m i > with bi ortho mos
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: i, a
+  double precision, intent(out) :: integral
+
+  integer                       :: ipoint, xi, j
+  double precision              :: weight, tmp
+
+  integral = 0.d0
+  do xi = 1, 3
+    do ipoint = 1, n_points_final_grid
+      weight = final_weight_at_r_vector(ipoint)
+        
+      tmp = 0.d0
+      do j = 1, elec_beta_num
+        tmp = tmp + x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,i) 
+      enddo
+
+      integral += ( mos_l_in_r_array_transp(ipoint,a) * W1_r_in_r(ipoint,xi,i) * w_sum_in_r(ipoint,xi) & 
+                  + tc_scf_dm_in_r(ipoint) * tmp                                                       &
+                  + mos_r_in_r_array_transp(ipoint,i) * W1_l_in_r(ipoint,xi,a) * w_sum_in_r(ipoint,xi) & 
+                  ) * weight
+
+    enddo
+  enddo
+
+end
+
+! ---
+
+subroutine exch_term_ijm_bi_ortho(a, i, integral)
+
+  BEGIN_DOC
+  ! computes sum_(j,m = 1, elec_beta_num) < a m j | i j m > with bi ortho mos
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: i, a
+  double precision, intent(out) :: integral
+
+  integer                       :: ipoint, xi
+  double precision              :: weight
+
+  integral = 0.d0
+  do xi = 1, 3
+    do ipoint = 1, n_points_final_grid
+      weight = final_weight_at_r_vector(ipoint)
+        
+      integral += ( mos_l_in_r_array_transp(ipoint,a) * mos_r_in_r_array_transp(ipoint,i) * v_sum_in_r(ipoint,xi) & 
+                  + 2.d0 * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) * W1_in_r(ipoint,xi)                             &
+                  ) * weight
+
+    enddo
+  enddo
+
+end
+
+! ---
+
+subroutine direct_term_ijj_bi_ortho(a, i, integral)
+
+  BEGIN_DOC
+  ! computes sum_(j = 1, elec_beta_num) < a j j | i j j > with bi ortho mos
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: i, a
+  double precision, intent(out) :: integral
+
+  integer                       :: ipoint, xi
+  double precision              :: weight
+
+  integral = 0.d0
+  do xi = 1, 3
+    do ipoint = 1, n_points_final_grid
+      weight = final_weight_at_r_vector(ipoint)
+
+      integral += ( mos_l_in_r_array_transp(ipoint,a) * mos_r_in_r_array_transp(ipoint,i) * ww_sum_in_r(ipoint,xi) & 
+                  + 2.d0 * W1_diag_in_r(ipoint, xi) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i)                        &
+                  ) * weight
+    enddo
+  enddo
+
+end 
+
+! ---
+
+subroutine cyclic_term_jim_bi_ortho(a, i, integral)
+
+  BEGIN_DOC
+  ! computes sum_(j,m = 1, elec_beta_num) < a m j | j i m > with bi ortho mos
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: i, a
+  double precision, intent(out) :: integral
+
+  integer                       :: ipoint, xi
+  double precision              :: weight
+
+  integral = 0.d0
+  do xi = 1, 3
+    do ipoint = 1, n_points_final_grid
+      weight = final_weight_at_r_vector(ipoint)
+        
+      integral += ( mos_l_in_r_array_transp(ipoint,a) * W1_W1_r_in_r(ipoint,xi,i)  & 
+                  + W1_W1_l_in_r(ipoint,xi,a) * mos_r_in_r_array_transp(ipoint,i)  &
+                  + W1_l_in_r(ipoint,xi,a) * W1_r_in_r(ipoint,xi,i)                &
+                  ) * weight
+
+    enddo
+  enddo
+
+end
+
+! ---
+
+subroutine cyclic_term_mji_bi_ortho(a, i, integral)
+
+  BEGIN_DOC
+  ! computes sum_(j,m = 1, elec_beta_num) < a m j | m j i > with bi ortho mos
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: i, a
+  double precision, intent(out) :: integral
+
+  integer                       :: ipoint, xi
+  double precision              :: weight
+
+  integral = 0.d0
+  do xi = 1, 3
+    do ipoint = 1, n_points_final_grid
+      weight = final_weight_at_r_vector(ipoint)
+        
+      integral += ( mos_l_in_r_array_transp(ipoint,a) * W1_W1_r_in_r(ipoint,xi,i)  & 
+                  + W1_l_in_r(ipoint,xi,a) * W1_r_in_r(ipoint,xi,i)                &
+                  + W1_W1_l_in_r(ipoint,xi,a) * mos_r_in_r_array_transp(ipoint,i)  &
+                  ) * weight
+
+    enddo
+  enddo
+
+end
+
+! ---
+

src/tc_scf/minimize_tc_angles.irp.f

@@ -0,0 +1,12 @@
+program print_angles
+ implicit none
+  my_grid_becke  = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+!  my_n_pt_r_grid = 10 ! small grid for quick debug
+!  my_n_pt_a_grid = 14 ! small grid for quick debug
+  touch my_n_pt_r_grid my_n_pt_a_grid
+!  call sort_by_tc_fock
+  call minimize_tc_orb_angles
+end
+

src/tc_scf/molden_lr_mos.irp.f

@@ -0,0 +1,176 @@
+program molden
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+
+  implicit none
+
+  print *, 'starting ...'
+
+  my_grid_becke  = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+!  my_n_pt_r_grid = 10 ! small grid for quick debug
+!  my_n_pt_a_grid = 26 ! small grid for quick debug
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  call molden_lr
+end
+subroutine molden_lr
+  implicit none
+  BEGIN_DOC
+  ! Produces a Molden file
+  END_DOC
+  character*(128)                :: output
+  integer                        :: i_unit_output,getUnitAndOpen
+  integer                        :: i,j,k,l
+  double precision, parameter :: a0 = 0.529177249d0
+
+  PROVIDE ezfio_filename
+
+  output=trim(ezfio_filename)//'.mol'
+  print*,'output = ',trim(output)
+
+  i_unit_output = getUnitAndOpen(output,'w')
+
+  write(i_unit_output,'(A)') '[Molden Format]'
+
+  write(i_unit_output,'(A)') '[Atoms] Angs'
+  do i = 1, nucl_num
+    write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))')                   &
+        trim(element_name(int(nucl_charge(i)))),                     &
+        i,                                                           &
+        int(nucl_charge(i)),                                         &
+        nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0
+  enddo
+
+  write(i_unit_output,'(A)') '[GTO]'
+
+  character*(1)                  :: character_shell
+  integer                        :: i_shell,i_prim,i_ao
+  integer                        :: iorder(ao_num)
+  integer                        :: nsort(ao_num)
+
+  i_shell = 0
+  i_prim = 0
+  do i=1,nucl_num
+    write(i_unit_output,*) i, 0
+    do j=1,nucl_num_shell_aos(i)
+      i_shell +=1
+      i_ao = nucl_list_shell_aos(i,j)
+      character_shell = trim(ao_l_char(i_ao))
+      write(i_unit_output,*) character_shell, ao_prim_num(i_ao), '1.00'
+      do k = 1, ao_prim_num(i_ao)
+        i_prim +=1
+        write(i_unit_output,'(E20.10,2X,E20.10)') ao_expo(i_ao,k), ao_coef(i_ao,k)
+      enddo
+      l = i_ao
+      do while ( ao_l(l) == ao_l(i_ao) )
+        nsort(l) = i*10000 + j*100
+        l += 1
+        if (l > ao_num) exit
+      enddo
+    enddo
+    write(i_unit_output,*)''
+  enddo
+
+
+  do i=1,ao_num
+    iorder(i) = i
+    ! p
+    if      ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 3
+    ! d
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 6
+    ! f
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 6
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 7
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 8
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 9
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 10
+    ! g
+    else if ((ao_power(i,1) == 4 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 4 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 4 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 6
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 7
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 8
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 9
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 10
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 11
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 12
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 13
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 14
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 15
+    endif
+  enddo
+
+  call isort(nsort,iorder,ao_num)
+  write(i_unit_output,'(A)') '[MO]'
+  do i=1,mo_num
+    write (i_unit_output,*) 'Sym= 1'
+    write (i_unit_output,*) 'Ene=', Fock_matrix_tc_mo_tot(i,i)
+    write (i_unit_output,*) 'Spin= Alpha'
+    write (i_unit_output,*) 'Occup=', mo_occ(i)
+    do j=1,ao_num
+      write(i_unit_output, '(I6,2X,E20.10)') j, mo_r_coef(iorder(j),i)
+    enddo
+
+    write (i_unit_output,*) 'Sym= 1'
+    write (i_unit_output,*) 'Ene=', Fock_matrix_tc_mo_tot(i,i)
+    write (i_unit_output,*) 'Spin= Alpha'
+    write (i_unit_output,*) 'Occup=', mo_occ(i)
+    do j=1,ao_num
+      write(i_unit_output, '(I6,2X,E20.10)') j, mo_l_coef(iorder(j),i)
+    enddo
+  enddo
+  close(i_unit_output)
+end
+

src/tc_scf/print_angle_tc_orb.irp.f

@@ -0,0 +1,9 @@
+program print_angles
+ implicit none
+  my_grid_becke  = .True.
+!  my_n_pt_r_grid = 30
+!  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = 10 ! small grid for quick debug
+  my_n_pt_a_grid = 14 ! small grid for quick debug
+  call print_angles_tc
+end

src/tc_scf/print_fit_param.irp.f

@@ -0,0 +1,60 @@
+program print_fit_param
+
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+
+  implicit none
+
+  my_grid_becke  = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+!  my_n_pt_r_grid = 10 ! small grid for quick debug
+!  my_n_pt_a_grid = 26 ! small grid for quick debug
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  !call create_guess
+  !call orthonormalize_mos
+
+  call main()
+
+end
+
+! ---
+
+subroutine main()
+
+  implicit none
+  integer :: i
+
+  mu_erf = 1.d0
+  touch mu_erf
+
+  print *, ' fit for (1 - erf(x))^2'
+  do i = 1, n_max_fit_slat
+    print*, expo_gauss_1_erf_x_2(i), coef_gauss_1_erf_x_2(i)
+  enddo
+
+  print *, ''
+  print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)]'
+  do i = 1, n_max_fit_slat
+    print *, expo_gauss_j_mu_x(i), 2.d0 * coef_gauss_j_mu_x(i)
+  enddo
+
+  print *, ''
+  print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)]^2'
+  do i = 1, n_max_fit_slat
+    print *, expo_gauss_j_mu_x_2(i), 4.d0 * coef_gauss_j_mu_x_2(i)
+  enddo
+
+  print *, ''
+  print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)] x [1 - erf(mu * r12)]'
+  do i = 1, n_max_fit_slat
+    print *, expo_gauss_j_mu_1_erf(i), 4.d0 * coef_gauss_j_mu_1_erf(i)
+  enddo
+
+  return
+end subroutine main
+
+! ---
+

src/tc_scf/rh_tcscf.irp.f

@@ -0,0 +1,336 @@
+! ---
+
+subroutine rh_tcscf()
+
+  BEGIN_DOC
+  !
+  ! Roothaan-Hall algorithm for TC-SCF calculation
+  !
+  END_DOC
+
+  implicit none
+
+  integer                       :: i, j
+  integer                       :: iteration_TCSCF, dim_DIIS, index_dim_DIIS
+  double precision              :: energy_TCSCF, energy_TCSCF_1e, energy_TCSCF_2e, energy_TCSCF_3e, gradie_TCSCF
+  double precision              :: energy_TCSCF_previous, delta_energy_TCSCF
+  double precision              :: gradie_TCSCF_previous, delta_gradie_TCSCF
+  double precision              :: max_error_DIIS_TCSCF
+  double precision              :: level_shift_save
+  double precision              :: delta_energy_tmp, delta_gradie_tmp
+  double precision, allocatable :: F_DIIS(:,:,:), e_DIIS(:,:,:)
+  double precision, allocatable :: mo_r_coef_save(:,:), mo_l_coef_save(:,:)
+
+  logical, external             :: qp_stop
+
+
+  !PROVIDE ao_md5 mo_occ
+  PROVIDE level_shift_TCSCF
+
+  allocate( mo_r_coef_save(ao_num,mo_num), mo_l_coef_save(ao_num,mo_num) &
+          , F_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF), e_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF) )
+
+  F_DIIS         = 0.d0
+  e_DIIS         = 0.d0
+  mo_l_coef_save = 0.d0
+  mo_r_coef_save = 0.d0
+
+  call write_time(6)
+
+  ! ---
+  ! Initialize energies and density matrices
+
+  energy_TCSCF_previous = TC_HF_energy
+  energy_TCSCF_1e       = TC_HF_one_e_energy
+  energy_TCSCF_2e       = TC_HF_two_e_energy
+  energy_TCSCF_3e       = 0.d0
+  if(three_body_h_tc) then
+    energy_TCSCF_3e     = diag_three_elem_hf
+  endif
+  gradie_TCSCF_previous = grad_non_hermit
+  delta_energy_TCSCF    = 1.d0
+  delta_gradie_TCSCF    = 1.d0
+  iteration_TCSCF       = 0
+  dim_DIIS              = 0
+  max_error_DIIS_TCSCF  = 1.d0
+
+  ! ---
+
+  ! Start of main SCF loop
+
+  PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot
+
+  do while( (max_error_DIIS_TCSCF > threshold_DIIS_nonzero_TCSCF) .or. &
+            !(dabs(delta_energy_TCSCF) > thresh_TCSCF)             .or. &
+            (dabs(gradie_TCSCF_previous) > dsqrt(thresh_TCSCF))        )
+
+    iteration_TCSCF += 1
+    if(iteration_TCSCF > n_it_TCSCF_max) then
+      print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max
+      stop
+    endif
+
+    dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF)
+
+    ! ---
+
+    if((tcscf_algorithm == 'DIIS') .and. (dabs(delta_energy_TCSCF) > 1.d-6))  then
+
+      ! store Fock and error matrices at each iteration
+      index_dim_DIIS = mod(dim_DIIS-1, max_dim_DIIS_TCSCF) + 1
+      do j = 1, ao_num
+        do i = 1, ao_num
+          F_DIIS(i,j,index_dim_DIIS) = Fock_matrix_tc_ao_tot(i,j)
+          e_DIIS(i,j,index_dim_DIIS) = FQS_SQF_ao(i,j)
+        enddo
+      enddo
+
+      call extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), iteration_TCSCF, dim_DIIS)
+
+      Fock_matrix_tc_ao_alpha = 0.5d0 * Fock_matrix_tc_ao_tot
+      Fock_matrix_tc_ao_beta  = 0.5d0 * Fock_matrix_tc_ao_tot
+      !TOUCH Fock_matrix_tc_ao_alpha Fock_matrix_tc_ao_beta
+
+      call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
+                            , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
+      call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta , size(Fock_matrix_tc_ao_beta , 1) &
+                            , Fock_matrix_tc_mo_beta , size(Fock_matrix_tc_mo_beta , 1) )
+      TOUCH Fock_matrix_tc_mo_alpha Fock_matrix_tc_mo_beta
+    endif
+
+    ! ---
+
+    mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num)
+    mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num)
+    TOUCH mo_l_coef mo_r_coef
+
+    ! ---
+
+    ! calculate error vectors
+    max_error_DIIS_TCSCF = maxval(abs(FQS_SQF_mo))
+
+    ! ---
+
+    delta_energy_tmp = TC_HF_energy    - energy_TCSCF_previous
+    delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous
+
+    ! ---
+
+    do while((delta_gradie_tmp > 1.d-7) .and. (iteration_TCSCF > 1))
+    !do while((dabs(delta_energy_tmp) > 0.5d0) .and. (iteration_TCSCF > 1))
+      print *, ' very big or bad step  : ', delta_energy_tmp, delta_gradie_tmp
+      print *, ' TC level shift = ', level_shift_TCSCF
+
+      mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num) 
+      mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num) 
+
+      if(level_shift_TCSCF <= .1d0) then
+        level_shift_TCSCF = 1.d0
+      else
+        level_shift_TCSCF = level_shift_TCSCF * 3.0d0
+      endif
+      TOUCH mo_l_coef mo_r_coef level_shift_TCSCF
+
+      mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num)
+      mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num)
+      TOUCH mo_l_coef mo_r_coef
+
+      delta_energy_tmp = TC_HF_energy    - energy_TCSCF_previous
+      delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous
+
+      if(level_shift_TCSCF - level_shift_save > 40.d0) then
+        level_shift_TCSCF = level_shift_save * 4.d0
+        SOFT_TOUCH level_shift_TCSCF
+        exit
+      endif
+
+      dim_DIIS = 0
+    enddo
+!    print *, ' very big step  : ', delta_energy_tmp
+!    print *, ' TC level shift = ', level_shift_TCSCF
+
+    ! ---
+
+    level_shift_TCSCF = 0.d0
+    !level_shift_TCSCF = level_shift_TCSCF * 0.5d0
+    SOFT_TOUCH level_shift_TCSCF
+
+    gradie_TCSCF       = grad_non_hermit
+    energy_TCSCF       = TC_HF_energy
+    energy_TCSCF_1e    = TC_HF_one_e_energy
+    energy_TCSCF_2e    = TC_HF_two_e_energy
+    energy_TCSCF_3e    = 0.d0
+    if(three_body_h_tc) then
+      energy_TCSCF_3e  = diag_three_elem_hf
+    endif
+    delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous
+    delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous
+
+    energy_TCSCF_previous = energy_TCSCF
+    gradie_TCSCF_previous = gradie_TCSCF
+
+
+    level_shift_save = level_shift_TCSCF
+    mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num)
+    mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num)
+
+
+    print *, ' iteration         = ', iteration_TCSCF
+    print *, ' total TC energy   = ', energy_TCSCF 
+    print *, ' 1-e   TC energy   = ', energy_TCSCF_1e
+    print *, ' 2-e   TC energy   = ', energy_TCSCF_2e
+    print *, ' 3-e   TC energy   = ', energy_TCSCF_3e
+    print *, ' |delta TC energy| = ', dabs(delta_energy_TCSCF)
+    print *, ' TC gradient       = ', gradie_TCSCF
+    print *, ' delta TC gradient = ', delta_gradie_TCSCF
+    print *, ' max TC DIIS error = ', max_error_DIIS_TCSCF 
+    print *, ' TC DIIS dim       = ', dim_DIIS
+    print *, ' TC level shift    = ', level_shift_TCSCF
+    print *, ' '
+
+    call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+    call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+
+    if(qp_stop()) exit
+  enddo
+
+  ! ---
+
+  print *, ' TCSCF DIIS converged !'
+  call print_energy_and_mos()
+
+  call write_time(6)
+
+  deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, e_DIIS)
+
+end
+
+! ---
+
+subroutine extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, F_ao, size_F_ao, iteration_TCSCF, dim_DIIS)
+
+  BEGIN_DOC
+  !
+  ! Compute the extrapolated Fock matrix using the DIIS procedure
+  !
+  ! e = \sum_i c_i e_i and \sum_i c_i = 1 
+  ! ==> lagrange multiplier with L = |e|^2 - \lambda (\sum_i c_i = 1)
+  !
+  END_DOC
+
+  implicit none
+
+  integer,          intent(in)    :: iteration_TCSCF, size_F_ao
+  integer,          intent(inout) :: dim_DIIS
+  double precision, intent(in)    :: F_DIIS(ao_num,ao_num,dim_DIIS)
+  double precision, intent(in)    :: e_DIIS(ao_num,ao_num,dim_DIIS)
+  double precision, intent(inout) :: F_ao(size_F_ao,ao_num)
+
+  double precision, allocatable   :: B_matrix_DIIS(:,:), X_vector_DIIS(:), C_vector_DIIS(:)
+
+  integer                         :: i, j, k, l, i_DIIS, j_DIIS
+  integer                         :: lwork
+  double precision                :: rcond, ferr, berr
+  integer,          allocatable   :: iwork(:)
+  double precision, allocatable   :: scratch(:,:)
+
+  if(dim_DIIS < 1) then
+    return
+  endif
+
+  allocate( B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), X_vector_DIIS(dim_DIIS+1) &
+          , C_vector_DIIS(dim_DIIS+1), scratch(ao_num,ao_num) )
+
+  ! Compute the matrices B and X
+  B_matrix_DIIS(:,:) = 0.d0
+  do j = 1, dim_DIIS
+    j_DIIS = min(dim_DIIS, mod(iteration_TCSCF-j, max_dim_DIIS_TCSCF)+1)
+
+    do i = 1, dim_DIIS
+      i_DIIS = min(dim_DIIS, mod(iteration_TCSCF-i, max_dim_DIIS_TCSCF)+1)
+
+      ! Compute product of two errors vectors
+      do l = 1, ao_num
+        do k = 1, ao_num
+          B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + e_DIIS(k,l,i_DIIS) * e_DIIS(k,l,j_DIIS)
+        enddo
+      enddo
+
+    enddo
+  enddo
+
+  ! Pad B matrix and build the X matrix
+
+  C_vector_DIIS(:) = 0.d0
+  do i = 1, dim_DIIS
+    B_matrix_DIIS(i,dim_DIIS+1) = -1.d0
+    B_matrix_DIIS(dim_DIIS+1,i) = -1.d0
+  enddo
+  C_vector_DIIS(dim_DIIS+1) = -1.d0
+
+  deallocate(scratch)
+
+  ! Estimate condition number of B
+  integer                       :: info
+  double precision              :: anorm
+  integer,          allocatable :: ipiv(:)
+  double precision, allocatable :: AF(:,:)
+  double precision, external :: dlange
+
+  lwork = max((dim_DIIS+1)**2, (dim_DIIS+1)*5)
+  allocate(AF(dim_DIIS+1,dim_DIIS+1))
+  allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+1)) )
+  allocate(scratch(lwork,1))
+  scratch(:,1) = 0.d0
+
+  anorm = dlange('1', dim_DIIS+1, dim_DIIS+1, B_matrix_DIIS, size(B_matrix_DIIS, 1), scratch(1,1))
+
+  AF(:,:) = B_matrix_DIIS(:,:)
+  call dgetrf(dim_DIIS+1, dim_DIIS+1, AF, size(AF, 1), ipiv, info)
+  if(info /= 0) then
+    dim_DIIS = 0
+    return
+  endif
+
+  call dgecon('1', dim_DIIS+1, AF, size(AF, 1), anorm, rcond, scratch, iwork, info)
+  if(info /= 0) then
+    dim_DIIS = 0
+    return
+  endif
+
+  if(rcond < 1.d-14) then
+    dim_DIIS = 0
+    return
+  endif
+
+  ! solve the linear system C = B x X
+
+  X_vector_DIIS = C_vector_DIIS
+  call dgesv(dim_DIIS+1, 1, B_matrix_DIIS, size(B_matrix_DIIS, 1), ipiv , X_vector_DIIS, size(X_vector_DIIS, 1), info)
+
+  deallocate(scratch, AF, iwork)
+  if(info < 0) then
+    stop ' bug in TC-DIIS'
+  endif
+
+  ! Compute extrapolated Fock matrix
+
+  !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200)
+  do j = 1, ao_num
+    do i = 1, ao_num
+      F_ao(i,j) = 0.d0
+    enddo
+    do k = 1, dim_DIIS
+      if(dabs(X_vector_DIIS(k)) < 1.d-10) cycle
+      do i = 1,ao_num
+        ! FPE here
+        F_ao(i,j) = F_ao(i,j) + X_vector_DIIS(k) * F_DIIS(i,j,dim_DIIS-k+1)
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+end
+
+! ---
+

src/tc_scf/rh_tcscf_diis.irp.f

@@ -0,0 +1,362 @@
+! ---
+
+subroutine rh_tcscf_diis()
+
+  implicit none
+
+  integer                       :: i, j, it
+  integer                       :: dim_DIIS, index_dim_DIIS
+  double precision              :: etc_tot, etc_1e, etc_2e, etc_3e, e_save, e_delta
+  double precision              :: tc_grad, g_save, g_delta, g_delta_th
+  double precision              :: level_shift_save, rate_th
+  double precision              :: t0, t1
+  double precision              :: er_DIIS, er_delta, er_save, er_delta_th
+  double precision, allocatable :: F_DIIS(:,:,:), E_DIIS(:,:,:)
+  double precision, allocatable :: mo_r_coef_save(:,:), mo_l_coef_save(:,:)
+
+  logical, external             :: qp_stop
+
+  it          = 0
+  e_save      = 0.d0
+  dim_DIIS    = 0
+  g_delta_th  = 1d0
+  er_delta_th = 1d0
+  rate_th     = 100.d0 !0.01d0 !0.2d0
+
+  allocate(mo_r_coef_save(ao_num,mo_num), mo_l_coef_save(ao_num,mo_num))
+  mo_l_coef_save = 0.d0
+  mo_r_coef_save = 0.d0
+
+  allocate(F_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF), E_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF))
+  F_DIIS = 0.d0
+  E_DIIS = 0.d0
+
+  call write_time(6)
+
+  ! ---
+
+  PROVIDE level_shift_TCSCF
+  PROVIDE mo_l_coef mo_r_coef
+
+  write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')              &
+    '====', '================', '================', '================', '================', '================' &
+          , '================', '================', '================', '====', '========'
+
+  write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')              &
+    ' it ', '  SCF TC Energy ', '      E(1e)     ', '      E(2e)     ', '      E(3e)     ', '   energy diff  ' &
+          , '    gradient    ', '    DIIS error  ', '  level shift   ', 'DIIS', '  WT (m)'
+
+  write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')              &
+    '====', '================', '================', '================', '================', '================' &
+          , '================', '================', '================', '====', '========'
+
+
+  ! first iteration (HF orbitals)
+  call wall_time(t0)
+
+  etc_tot = TC_HF_energy
+  etc_1e  = TC_HF_one_e_energy
+  etc_2e  = TC_HF_two_e_energy
+  etc_3e  = 0.d0
+  if(three_body_h_tc) then
+    etc_3e = diag_three_elem_hf
+  endif
+  tc_grad = grad_non_hermit
+  er_DIIS = maxval(abs(FQS_SQF_mo))
+  e_delta = dabs(etc_tot - e_save)
+
+  e_save  = etc_tot
+  g_save  = tc_grad
+  er_save = er_DIIS
+
+  call wall_time(t1)
+  write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)')  &
+    it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
+
+  ! ---
+
+  PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot
+
+  do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. threshold_DIIS_nonzero_TCSCF))
+
+    call wall_time(t0)
+
+    it += 1
+    if(it > n_it_TCSCF_max) then
+      print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max
+      stop
+    endif
+
+    dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF)
+
+    ! ---
+
+    if(dabs(e_delta) > 1.d-12) then
+
+      index_dim_DIIS = mod(dim_DIIS-1, max_dim_DIIS_TCSCF) + 1
+      do j = 1, ao_num
+        do i = 1, ao_num
+          F_DIIS(i,j,index_dim_DIIS) = Fock_matrix_tc_ao_tot(i,j)
+          E_DIIS(i,j,index_dim_DIIS) = FQS_SQF_ao           (i,j)
+        enddo
+      enddo
+
+      call extrapolate_TC_Fock_matrix(E_DIIS, F_DIIS, Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), it, dim_DIIS)
+
+      call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) &
+                            , Fock_matrix_tc_mo_tot, size(Fock_matrix_tc_mo_tot, 1) )
+      TOUCH Fock_matrix_tc_mo_tot fock_matrix_tc_diag_mo_tot
+    endif
+
+    ! ---
+
+    mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num)
+    mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num)
+    !call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+    !call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+    TOUCH mo_l_coef mo_r_coef
+
+    ! ---
+
+    g_delta  = grad_non_hermit         -  g_save
+    er_delta = maxval(abs(FQS_SQF_mo)) - er_save
+
+    !if((g_delta > rate_th * g_delta_th) .and. (er_delta > rate_th * er_delta_th) .and. (it > 1)) then
+    if((g_delta > rate_th * g_delta_th) .and. (it > 1)) then
+    !if((g_delta > 0.d0) .and. (it > 1)) then
+
+      Fock_matrix_tc_ao_tot(1:ao_num,1:ao_num) = F_DIIS(1:ao_num,1:ao_num,index_dim_DIIS)
+      call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) &
+                            , Fock_matrix_tc_mo_tot, size(Fock_matrix_tc_mo_tot, 1) )
+      TOUCH Fock_matrix_tc_mo_tot fock_matrix_tc_diag_mo_tot
+
+      mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num)
+      mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num)
+      !call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+      !call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+      TOUCH mo_l_coef mo_r_coef
+
+    endif
+
+    ! ---
+
+    g_delta  = grad_non_hermit         -  g_save
+    er_delta = maxval(abs(FQS_SQF_mo)) - er_save
+
+    mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num)
+    mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num)
+
+    !do while((g_delta > rate_th * g_delta_th) .and. (er_delta > rate_th * er_delta_th) .and. (it > 1))
+    do while((g_delta > rate_th * g_delta_th) .and. (it > 1))
+      print *, ' big or bad step : ', g_delta, rate_th * g_delta_th
+
+      mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num) 
+      mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num) 
+      if(level_shift_TCSCF <= .1d0) then
+        level_shift_TCSCF = 1.d0
+      else
+        level_shift_TCSCF = level_shift_TCSCF * 3.0d0
+      endif
+      TOUCH mo_l_coef mo_r_coef level_shift_TCSCF
+
+      mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num)
+      mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num)
+      !call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+      !call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+      TOUCH mo_l_coef mo_r_coef
+
+      g_delta  = grad_non_hermit         -  g_save
+      er_delta = maxval(abs(FQS_SQF_mo)) - er_save
+
+      if(level_shift_TCSCF - level_shift_save > 40.d0) then
+        level_shift_TCSCF = level_shift_save * 4.d0
+        SOFT_TOUCH level_shift_TCSCF
+        exit
+      endif
+
+      dim_DIIS = 0
+    enddo
+
+    ! ---
+
+    level_shift_TCSCF = level_shift_TCSCF * 0.5d0
+    SOFT_TOUCH level_shift_TCSCF
+
+    etc_tot = TC_HF_energy
+    etc_1e  = TC_HF_one_e_energy
+    etc_2e  = TC_HF_two_e_energy
+    etc_3e  = 0.d0
+    if(three_body_h_tc) then
+      etc_3e = diag_three_elem_hf
+    endif
+    tc_grad  = grad_non_hermit
+    er_DIIS  = maxval(abs(FQS_SQF_mo))
+    e_delta  = dabs(etc_tot - e_save)
+    g_delta  = tc_grad - g_save
+    er_delta = er_DIIS - er_save
+    
+    e_save           = etc_tot
+    g_save           = tc_grad
+    level_shift_save = level_shift_TCSCF
+    er_save          = er_DIIS
+
+    g_delta_th  = dabs(tc_grad) ! g_delta)
+    er_delta_th = dabs(er_DIIS) !er_delta)
+
+    call wall_time(t1)
+    write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)')  &
+      it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
+
+    if(g_delta .lt. 0.d0) then
+      call ezfio_set_tc_scf_bitc_energy(etc_tot)
+      call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+      call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+    endif
+
+    if(qp_stop()) exit
+  enddo
+
+  ! ---
+
+  print *, ' TCSCF DIIS converged !'
+  call print_energy_and_mos()
+
+  call write_time(6)
+
+  deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, E_DIIS)
+
+  call ezfio_set_tc_scf_bitc_energy(TC_HF_energy)
+  call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+
+end
+
+! ---
+
+subroutine extrapolate_TC_Fock_matrix(E_DIIS, F_DIIS, F_ao, size_F_ao, it, dim_DIIS)
+
+  BEGIN_DOC
+  !
+  ! Compute the extrapolated Fock matrix using the DIIS procedure
+  !
+  ! e = \sum_i c_i e_i and \sum_i c_i = 1 
+  ! ==> lagrange multiplier with L = |e|^2 - \lambda (\sum_i c_i = 1)
+  !
+  END_DOC
+
+  implicit none
+
+  integer,          intent(in)    :: it, size_F_ao
+  integer,          intent(inout) :: dim_DIIS
+  double precision, intent(in)    :: F_DIIS(ao_num,ao_num,dim_DIIS)
+  double precision, intent(in)    :: E_DIIS(ao_num,ao_num,dim_DIIS)
+  double precision, intent(inout) :: F_ao(size_F_ao,ao_num)
+
+  double precision, allocatable   :: B_matrix_DIIS(:,:), X_vector_DIIS(:), C_vector_DIIS(:)
+
+  integer                         :: i, j, k, l, i_DIIS, j_DIIS
+  integer                         :: lwork
+  double precision                :: rcond, ferr, berr
+  integer,          allocatable   :: iwork(:)
+  double precision, allocatable   :: scratch(:,:)
+
+  if(dim_DIIS < 1) then
+    return
+  endif
+
+  allocate( B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), X_vector_DIIS(dim_DIIS+1) &
+          , C_vector_DIIS(dim_DIIS+1), scratch(ao_num,ao_num) )
+
+  ! Compute the matrices B and X
+  B_matrix_DIIS(:,:) = 0.d0
+  do j = 1, dim_DIIS
+    j_DIIS = min(dim_DIIS, mod(it-j, max_dim_DIIS_TCSCF)+1)
+
+    do i = 1, dim_DIIS
+      i_DIIS = min(dim_DIIS, mod(it-i, max_dim_DIIS_TCSCF)+1)
+
+      ! Compute product of two errors vectors
+      do l = 1, ao_num
+        do k = 1, ao_num
+          B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + E_DIIS(k,l,i_DIIS) * E_DIIS(k,l,j_DIIS)
+        enddo
+      enddo
+
+    enddo
+  enddo
+
+  ! Pad B matrix and build the X matrix
+
+  C_vector_DIIS(:) = 0.d0
+  do i = 1, dim_DIIS
+    B_matrix_DIIS(i,dim_DIIS+1) = -1.d0
+    B_matrix_DIIS(dim_DIIS+1,i) = -1.d0
+  enddo
+  C_vector_DIIS(dim_DIIS+1) = -1.d0
+
+  deallocate(scratch)
+
+  ! Estimate condition number of B
+  integer                       :: info
+  double precision              :: anorm
+  integer,          allocatable :: ipiv(:)
+  double precision, allocatable :: AF(:,:)
+  double precision, external :: dlange
+
+  lwork = max((dim_DIIS+1)**2, (dim_DIIS+1)*5)
+  allocate(AF(dim_DIIS+1,dim_DIIS+1))
+  allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+1)) )
+  allocate(scratch(lwork,1))
+  scratch(:,1) = 0.d0
+
+  anorm = dlange('1', dim_DIIS+1, dim_DIIS+1, B_matrix_DIIS, size(B_matrix_DIIS, 1), scratch(1,1))
+
+  AF(:,:) = B_matrix_DIIS(:,:)
+  call dgetrf(dim_DIIS+1, dim_DIIS+1, AF, size(AF, 1), ipiv, info)
+  if(info /= 0) then
+    dim_DIIS = 0
+    return
+  endif
+
+  call dgecon('1', dim_DIIS+1, AF, size(AF, 1), anorm, rcond, scratch, iwork, info)
+  if(info /= 0) then
+    dim_DIIS = 0
+    return
+  endif
+
+  if(rcond < 1.d-14) then
+    dim_DIIS = 0
+    return
+  endif
+
+  ! solve the linear system C = B x X
+
+  X_vector_DIIS = C_vector_DIIS
+  call dgesv(dim_DIIS+1, 1, B_matrix_DIIS, size(B_matrix_DIIS, 1), ipiv , X_vector_DIIS, size(X_vector_DIIS, 1), info)
+
+  deallocate(scratch, AF, iwork)
+  if(info < 0) then
+    stop ' bug in TC-DIIS'
+  endif
+
+  ! Compute extrapolated Fock matrix
+
+  !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200)
+  do j = 1, ao_num
+    do i = 1, ao_num
+      F_ao(i,j) = 0.d0
+    enddo
+    do k = 1, dim_DIIS
+      if(dabs(X_vector_DIIS(k)) < 1.d-10) cycle
+      do i = 1,ao_num
+        ! FPE here
+        F_ao(i,j) = F_ao(i,j) + X_vector_DIIS(k) * F_DIIS(i,j,dim_DIIS-k+1)
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+end
+
+! ---
+

src/tc_scf/rh_tcscf_simple.irp.f

@@ -0,0 +1,129 @@
+! ---
+
+subroutine rh_tcscf_simple()
+
+  implicit none
+  integer                       :: i, j, it, dim_DIIS
+  double precision              :: t0, t1
+  double precision              :: e_save, e_delta, rho_delta
+  double precision              :: etc_tot, etc_1e, etc_2e, etc_3e, tc_grad
+  double precision              :: er_DIIS
+  double precision, allocatable :: rho_old(:,:), rho_new(:,:)
+
+  allocate(rho_old(ao_num,ao_num), rho_new(ao_num,ao_num))
+
+  it       = 0
+  e_save   = 0.d0
+  dim_DIIS = 0
+
+  ! ---
+
+  if(.not. bi_ortho) then
+   print *, ' grad_hermit = ', grad_hermit
+   call save_good_hermit_tc_eigvectors
+   TOUCH mo_coef 
+   call save_mos
+  endif
+
+  ! ---
+
+  if(bi_ortho) then
+
+    PROVIDE level_shift_tcscf
+    PROVIDE mo_l_coef mo_r_coef
+
+    write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')              &
+      '====', '================', '================', '================', '================', '================' &
+            , '================', '================', '================', '====', '========'
+
+    write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')              &
+      ' it ', '  SCF TC Energy ', '      E(1e)     ', '      E(2e)     ', '      E(3e)     ', '   energy diff  ' &
+            , '    gradient    ', '    DIIS error  ', '  level shift   ', 'DIIS', '  WT (m)'
+
+    write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')              &
+      '====', '================', '================', '================', '================', '================' &
+            , '================', '================', '================', '====', '========'
+
+
+    ! first iteration (HF orbitals)
+    call wall_time(t0)
+
+    etc_tot = TC_HF_energy
+    etc_1e  = TC_HF_one_e_energy
+    etc_2e  = TC_HF_two_e_energy
+    etc_3e  = 0.d0
+    if(three_body_h_tc) then
+      etc_3e = diag_three_elem_hf
+    endif
+    tc_grad = grad_non_hermit
+    er_DIIS = maxval(abs(FQS_SQF_mo))
+    e_delta = dabs(etc_tot - e_save)
+    e_save  = etc_tot
+
+    call wall_time(t1)
+    write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)')  &
+      it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
+
+    do while(tc_grad .gt. dsqrt(thresh_tcscf))
+      call wall_time(t0)
+
+      it += 1
+      if(it > n_it_tcscf_max) then
+        print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max
+        stop
+      endif
+
+      mo_l_coef = fock_tc_leigvec_ao
+      mo_r_coef = fock_tc_reigvec_ao
+      call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+      call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+      TOUCH mo_l_coef mo_r_coef
+
+      etc_tot = TC_HF_energy
+      etc_1e  = TC_HF_one_e_energy
+      etc_2e  = TC_HF_two_e_energy
+      etc_3e  = 0.d0
+      if(three_body_h_tc) then
+        etc_3e = diag_three_elem_hf
+      endif
+      tc_grad = grad_non_hermit
+      er_DIIS = maxval(abs(FQS_SQF_mo))
+      e_delta = dabs(etc_tot - e_save)
+      e_save  = etc_tot
+
+      call ezfio_set_tc_scf_bitc_energy(etc_tot)
+
+      call wall_time(t1)
+      write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)')  &
+        it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
+    enddo
+
+  else
+
+   do while( (grad_hermit.gt.dsqrt(thresh_tcscf)) .and. (it.lt.n_it_tcscf_max) )
+      print*,'grad_hermit = ',grad_hermit
+      it += 1
+      print *, 'iteration = ', it
+      print *, '***'
+      print *, 'TC HF total energy = ', TC_HF_energy
+      print *, 'TC HF 1 e   energy = ', TC_HF_one_e_energy
+      print *, 'TC HF 2 e   energy = ', TC_HF_two_e_energy
+      print *, 'TC HF 3 body       = ', diag_three_elem_hf
+      print *, '***'
+      print *, ''
+      call save_good_hermit_tc_eigvectors
+      TOUCH mo_coef 
+      call save_mos
+    enddo
+
+  endif
+
+  print *, ' TCSCF Simple converged !'
+  call print_energy_and_mos()
+
+  deallocate(rho_old, rho_new)
+
+end
+
+! ---
+

src/tc_scf/rotate_tcscf_orbitals.irp.f

@@ -0,0 +1,367 @@
+
+! ---
+
+program rotate_tcscf_orbitals
+
+  BEGIN_DOC
+  ! TODO : Put the documentation of the program here
+  END_DOC
+
+  implicit none
+
+  my_grid_becke  = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  bi_ortho = .True.
+  touch bi_ortho
+
+  call maximize_overlap()
+
+end
+
+! ---
+
+subroutine maximize_overlap()
+
+  implicit none
+  integer                       :: i, m, n
+  double precision              :: accu_d, accu_nd
+  double precision, allocatable :: C(:,:), R(:,:), L(:,:), W(:,:), e(:)
+  double precision, allocatable :: S(:,:)
+
+  n = ao_num
+  m = mo_num
+
+  allocate(L(n,m), R(n,m), C(n,m), W(n,n), e(m))
+  L = mo_l_coef
+  R = mo_r_coef
+  C = mo_coef
+  W = ao_overlap
+
+  print*, ' fock matrix diag elements'
+  do i = 1, m
+    e(i) = Fock_matrix_tc_mo_tot(i,i)
+    print*, e(i)
+  enddo
+
+  ! ---
+   
+  print *, ' overlap before :'
+  print *, ' '
+
+  allocate(S(m,m)) 
+
+  call LTxSxR(n, m, L, W, R, S)
+  !print*, " L.T x R"
+  !do i = 1, m
+  !  write(*, '(100(F16.10,X))') S(i,i)
+  !enddo
+  call LTxSxR(n, m, L, W, C, S)
+  print*, " L.T x C"
+  do i = 1, m
+    write(*, '(100(F16.10,X))') S(i,:)
+  enddo
+  call LTxSxR(n, m, C, W, R, S)
+  print*, " C.T x R"
+  do i = 1, m
+    write(*, '(100(F16.10,X))') S(i,:)
+  enddo
+
+  deallocate(S)
+
+  ! ---
+
+  call rotate_degen_eigvec_to_maximize_overlap(n, m, e, C, W, L, R)
+
+  ! ---
+   
+  print *, ' overlap after :'
+  print *, ' '
+
+  allocate(S(m,m)) 
+
+  call LTxSxR(n, m, L, W, R, S)
+  !print*, " L.T x R"
+  !do i = 1, m
+  !  write(*, '(100(F16.10,X))') S(i,i)
+  !enddo
+  call LTxSxR(n, m, L, W, C, S)
+  print*, " L.T x C"
+  do i = 1, m
+    write(*, '(100(F16.10,X))') S(i,:)
+  enddo
+  call LTxSxR(n, m, C, W, R, S)
+  print*, " C.T x R"
+  do i = 1, m
+    write(*, '(100(F16.10,X))') S(i,:)
+  enddo
+
+  deallocate(S)
+
+  ! ---
+
+  mo_l_coef = L
+  mo_r_coef = R
+  call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+
+  ! ---
+
+  deallocate(L, R, C, W, e)
+
+end subroutine maximize_overlap
+
+! ---
+
+subroutine rotate_degen_eigvec_to_maximize_overlap(n, m, e0, C0, W0, L0, R0)
+
+  implicit none
+
+  integer,          intent(in)    :: n, m
+  double precision, intent(in)    :: e0(m), W0(n,n), C0(n,m)
+  double precision, intent(inout) :: L0(n,m), R0(n,m)
+
+
+  integer                         :: i, j, k, kk, mm, id1, tot_deg
+  double precision                :: ei, ej, de, de_thr
+  integer,          allocatable   :: deg_num(:)
+  double precision, allocatable   :: L(:,:), R(:,:), C(:,:), Lnew(:,:), Rnew(:,:), tmp(:,:)
+  !double precision, allocatable   :: S(:,:), Snew(:,:), T(:,:), Ttmp(:,:), Stmp(:,:)
+  double precision, allocatable   :: S(:,:), Snew(:,:), T(:,:), Ttmp(:,:), Stmp(:,:)
+  !real*8                          :: S(m,m), Snew(m,m), T(m,m)
+
+  id1 = 700
+  allocate(S(id1,id1), Snew(id1,id1), T(id1,id1))
+
+  ! ---
+
+  allocate( deg_num(m) )
+  do i = 1, m
+    deg_num(i) = 1
+  enddo
+
+  de_thr = thr_degen_tc
+
+  do i = 1, m-1
+    ei = e0(i)
+
+    ! already considered in degen vectors
+    if(deg_num(i).eq.0) cycle
+
+    do j = i+1, m
+      ej = e0(j)
+      de = dabs(ei - ej)
+
+      if(de .lt. de_thr) then
+        deg_num(i) = deg_num(i) + 1
+        deg_num(j) = 0
+      endif
+
+    enddo
+  enddo
+
+  tot_deg = 0
+  do i = 1, m
+    if(deg_num(i).gt.1) then
+      print *, ' degen on', i, deg_num(i)
+      tot_deg = tot_deg + 1
+    endif
+  enddo
+
+  if(tot_deg .eq. 0) then
+    print *, ' no degen'
+    return
+  endif
+
+  ! ---
+
+  do i = 1, m
+    mm = deg_num(i)
+
+    if(mm .gt. 1) then
+
+      allocate(L(n,mm), R(n,mm), C(n,mm))
+      do j = 1, mm
+        L(1:n,j) = L0(1:n,i+j-1)
+        R(1:n,j) = R0(1:n,i+j-1)
+        C(1:n,j) = C0(1:n,i+j-1)
+      enddo
+
+      ! ---
+
+      ! C.T x W0 x R
+      allocate(tmp(mm,n), Stmp(mm,mm))
+      call dgemm( 'T', 'N', mm, n, n, 1.d0       &
+                , C, size(C, 1), W0, size(W0, 1) &
+                , 0.d0, tmp, size(tmp, 1) )
+      call dgemm( 'N', 'N', mm, mm, n, 1.d0        &
+                , tmp, size(tmp, 1), R, size(R, 1) &
+                , 0.d0, Stmp, size(Stmp, 1) )
+      deallocate(C, tmp)
+
+      S = 0.d0
+      do k = 1, mm
+        do kk = 1, mm
+          S(kk,k) = Stmp(kk,k)
+        enddo
+      enddo
+      deallocate(Stmp)
+
+      !print*, " overlap bef"
+      !do k = 1, mm
+      !  write(*, '(100(F16.10,X))') (S(k,kk), kk=1, mm)
+      !enddo
+    
+      T    = 0.d0
+      Snew = 0.d0
+      call maxovl(mm, mm, S, T, Snew)
+
+      !print*, " overlap aft"
+      !do k = 1, mm
+      !  write(*, '(100(F16.10,X))') (Snew(k,kk), kk=1, mm)
+      !enddo
+
+      allocate(Ttmp(mm,mm))
+      Ttmp(1:mm,1:mm) = T(1:mm,1:mm)
+
+      allocate(Lnew(n,mm), Rnew(n,mm))
+      call dgemm( 'N', 'N', n, mm, mm, 1.d0               &
+                , R, size(R, 1), Ttmp(1,1), size(Ttmp, 1) &
+                , 0.d0, Rnew, size(Rnew, 1) )
+      call dgemm( 'N', 'N', n, mm, mm, 1.d0               &
+                , L, size(L, 1), Ttmp(1,1), size(Ttmp, 1) &
+                , 0.d0, Lnew, size(Lnew, 1) )
+
+      deallocate(L, R)
+      deallocate(Ttmp)
+
+      ! ---
+
+      do j = 1, mm
+        L0(1:n,i+j-1) = Lnew(1:n,j)
+        R0(1:n,i+j-1) = Rnew(1:n,j)
+      enddo
+      deallocate(Lnew, Rnew)
+
+    endif
+  enddo
+
+  deallocate(S, Snew, T)
+
+end subroutine rotate_degen_eigvec_to_maximize_overlap
+
+! ---
+
+subroutine fix_right_to_one()
+
+  implicit none
+  integer                       :: i, j, m, n, mm, tot_deg
+  double precision              :: accu_d, accu_nd
+  double precision              :: de_thr, ei, ej, de
+  integer,          allocatable :: deg_num(:)
+  double precision, allocatable :: R0(:,:), L0(:,:), W(:,:), e0(:)
+  double precision, allocatable :: R(:,:), L(:,:), S(:,:), Stmp(:,:), tmp(:,:)
+
+  n = ao_num
+  m = mo_num
+
+  allocate(L0(n,m), R0(n,m), W(n,n), e0(m))
+  L0 = mo_l_coef
+  R0 = mo_r_coef
+  W  = ao_overlap
+
+  print*, ' fock matrix diag elements'
+  do i = 1, m
+    e0(i) = Fock_matrix_tc_mo_tot(i,i)
+    print*, e0(i)
+  enddo
+
+  ! ---
+
+  allocate( deg_num(m) )
+  do i = 1, m
+    deg_num(i) = 1
+  enddo
+
+  de_thr = 1d-6
+
+  do i = 1, m-1
+    ei = e0(i)
+
+    ! already considered in degen vectors
+    if(deg_num(i).eq.0) cycle
+
+    do j = i+1, m
+      ej = e0(j)
+      de = dabs(ei - ej)
+
+      if(de .lt. de_thr) then
+        deg_num(i) = deg_num(i) + 1
+        deg_num(j) = 0
+      endif
+
+    enddo
+  enddo
+
+  deallocate(e0)
+
+  tot_deg = 0
+  do i = 1, m
+    if(deg_num(i).gt.1) then
+      print *, ' degen on', i, deg_num(i)
+      tot_deg = tot_deg + 1
+    endif
+  enddo
+
+  if(tot_deg .eq. 0) then
+    print *, ' no degen'
+    return
+  endif
+
+  ! ---
+
+  do i = 1, m
+    mm = deg_num(i)
+
+    if(mm .gt. 1) then
+
+      allocate(L(n,mm), R(n,mm))
+      do j = 1, mm
+        L(1:n,j) = L0(1:n,i+j-1)
+        R(1:n,j) = R0(1:n,i+j-1)
+      enddo
+
+      ! ---
+
+      call impose_weighted_orthog_svd(n, mm, W, R)
+      call impose_weighted_biorthog_qr(n, mm, thresh_biorthog_diag, thresh_biorthog_nondiag, R, W, L)
+
+      ! ---
+
+      do j = 1, mm
+        L0(1:n,i+j-1) = L(1:n,j)
+        R0(1:n,i+j-1) = R(1:n,j)
+      enddo
+      deallocate(L, R)
+
+    endif
+  enddo
+
+  call check_weighted_biorthog_binormalize(n, m, L0, W, R0, thresh_biorthog_diag, thresh_biorthog_nondiag, .true.)
+
+  deallocate(W, deg_num)
+
+  mo_l_coef = L0
+  mo_r_coef = R0
+  deallocate(L0, R0)
+
+  call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+  print *, ' orbitals are rotated '
+
+  return
+end subroutine fix_right_to_one
+
+! ---

src/tc_scf/routines_rotates.irp.f

@@ -0,0 +1,359 @@
+
+! ---
+
+subroutine minimize_tc_orb_angles()
+
+  implicit none
+  logical          :: good_angles
+  integer          :: i
+  double precision :: thr_deg
+
+  good_angles = .False.
+  thr_deg = thr_degen_tc
+
+  call print_energy_and_mos()
+
+  print *, ' Minimizing the angles between the TC orbitals'
+  i = 1
+  do while (.not. good_angles)
+    print *, ' iteration = ', i
+    call routine_save_rotated_mos(thr_deg, good_angles)
+    thr_deg *= 10.d0
+    i += 1
+    if(i .gt. 100) then
+      print *, ' minimize_tc_orb_angles does not seem to converge ..'
+      print *, ' Something is weird in the tc orbitals ...'
+      print *, ' STOPPING'
+      stop
+    endif
+  enddo
+  print *, ' Converged ANGLES MINIMIZATION !!'
+
+  call print_angles_tc()
+  call print_energy_and_mos()
+
+end
+
+! ---
+
+subroutine routine_save_rotated_mos(thr_deg, good_angles)
+
+  implicit none
+
+  double precision, intent(in)  :: thr_deg
+  logical,          intent(out) :: good_angles
+
+  integer                       :: i, j, k, n_degen_list, m, n, n_degen, ilast, ifirst
+  double precision              :: max_angle, norm
+  integer,          allocatable :: list_degen(:,:)
+  double precision, allocatable :: new_angles(:)
+  double precision, allocatable :: mo_r_coef_good(:,:), mo_l_coef_good(:,:)
+  double precision, allocatable :: mo_r_coef_new(:,:)
+  double precision, allocatable :: fock_diag(:),s_mat(:,:)
+  double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:)
+  double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:)
+
+  good_angles = .False.
+
+  allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
+
+  print *, ' ***************************************'
+  print *, ' ***************************************'
+  print *, ' THRESHOLD FOR DEGENERACIES ::: ', thr_deg
+  print *, ' ***************************************'
+  print *, ' ***************************************'
+  print *, ' Starting with the following TC energy gradient :', grad_non_hermit
+
+  mo_r_coef_good = mo_r_coef
+  mo_l_coef_good = mo_l_coef
+
+  allocate(mo_r_coef_new(ao_num, mo_num))
+  mo_r_coef_new = mo_r_coef
+  do i = 1, mo_num
+    norm = 1.d0/dsqrt(overlap_mo_r(i,i))
+    do j = 1, ao_num
+      mo_r_coef_new(j,i) *= norm
+    enddo
+  enddo
+
+  allocate(list_degen(mo_num,0:mo_num), s_mat(mo_num,mo_num), fock_diag(mo_num))
+  do i = 1, mo_num
+    fock_diag(i) = Fock_matrix_tc_mo_tot(i,i)
+  enddo
+
+ ! compute the overlap between the left and rescaled right
+  call build_s_matrix(ao_num, mo_num, mo_r_coef_new, mo_r_coef_new, ao_overlap, s_mat)
+! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
+  call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list)
+  print *, ' fock_matrix_mo'
+  do i = 1, mo_num
+    print *, i, fock_diag(i), angle_left_right(i)
+  enddo
+   
+  do i = 1, n_degen_list
+!  ifirst = list_degen(1,i)
+!  ilast  = list_degen(2,i)
+!  n_degen = ilast - ifirst +1
+
+    n_degen = list_degen(i,0)
+    if(n_degen .eq. 1) cycle
+
+    allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen))
+    allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen))
+    allocate(T(n_degen,n_degen), Snew(n_degen,n_degen))
+
+    do j = 1, n_degen
+      mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j))
+      mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
+    enddo
+    ! Orthogonalization of right functions
+    print *, ' Orthogonalization of RIGHT functions'
+    print *, ' ------------------------------------'
+    call orthog_functions(ao_num, n_degen, mo_r_coef_tmp, ao_overlap)
+  
+    ! Orthogonalization of left functions
+    print *, ' Orthogonalization of LEFT functions'
+    print *, ' ------------------------------------'
+    call orthog_functions(ao_num, n_degen, mo_l_coef_tmp, ao_overlap)
+
+    print *, ' Overlap left-right '
+    call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_l_coef_tmp, ao_overlap, stmp)
+    do j = 1, n_degen
+     write(*,'(100(F8.4,X))') stmp(:,j)
+    enddo
+    call build_s_matrix(ao_num, n_degen, mo_l_coef_tmp, mo_l_coef_tmp, ao_overlap, stmp)
+
+    !print*,'LEFT/LEFT OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_r_coef_tmp, ao_overlap, stmp)
+    !print*,'RIGHT/RIGHT OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+
+    if(maxovl_tc) then
+      T    = 0.d0
+      Snew = 0.d0
+      call maxovl(n_degen, n_degen, stmp, T, Snew)
+    !print*,'overlap after'
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')Snew(:,j)
+    !enddo
+      call dgemm( 'N', 'N', ao_num, n_degen, n_degen, 1.d0                  &
+                , mo_l_coef_tmp, size(mo_l_coef_tmp, 1), T(1,1), size(T, 1) &
+                , 0.d0, mo_l_coef_new, size(mo_l_coef_new, 1) )
+     call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_r_coef_tmp, ao_overlap, stmp)
+    !print*,'Overlap test'
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    else 
+      mo_l_coef_new = mo_l_coef_tmp
+    endif
+
+    call impose_weighted_biorthog_svd(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp)
+
+    !call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_r_coef_tmp, ao_overlap, stmp)
+    !print*,'LAST OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    !call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_l_coef_new, ao_overlap, stmp)
+    !print*,'LEFT OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    !call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_r_coef_tmp, ao_overlap, stmp)
+    !print*,'RIGHT OVERLAP '
+    !do j = 1, n_degen
+    ! write(*,'(100(F16.10,X))')stmp(:,j)
+    !enddo
+    do j = 1, n_degen
+!!!   mo_l_coef_good(1:ao_num,j+ifirst-1) = mo_l_coef_new(1:ao_num,j)
+!!!   mo_r_coef_good(1:ao_num,j+ifirst-1) = mo_r_coef_tmp(1:ao_num,j)
+      mo_l_coef_good(1:ao_num,list_degen(i,j)) = mo_l_coef_new(1:ao_num,j)
+      mo_r_coef_good(1:ao_num,list_degen(i,j)) = mo_r_coef_tmp(1:ao_num,j)
+    enddo
+
+    deallocate(stmp, smat2)
+    deallocate(mo_r_coef_tmp, mo_l_coef_tmp, mo_l_coef_new)
+    deallocate(T, Snew)
+  enddo
+
+  !allocate(stmp(mo_num, mo_num))
+  !call build_s_matrix(ao_num, mo_num, mo_l_coef_good, mo_r_coef_good, ao_overlap, stmp)
+  !print*,'LEFT/RIGHT OVERLAP '
+  !do j = 1, mo_num
+  ! write(*,'(100(F16.10,X))')stmp(:,j)
+  !enddo
+  !call build_s_matrix(ao_num, mo_num, mo_l_coef_good, mo_l_coef_good, ao_overlap, stmp)
+  !print*,'LEFT/LEFT OVERLAP '
+  !do j = 1, mo_num
+  ! write(*,'(100(F16.10,X))')stmp(:,j)
+  !enddo
+  !call build_s_matrix(ao_num, mo_num, mo_r_coef_good, mo_r_coef_good, ao_overlap, stmp)
+  !print*,'RIGHT/RIGHT OVERLAP '
+  !do j = 1, mo_num
+  ! write(*,'(100(F16.10,X))')stmp(:,j)
+  !enddo
+
+  mo_r_coef = mo_r_coef_good
+  mo_l_coef = mo_l_coef_good
+  call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
+  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
+  TOUCH mo_l_coef mo_r_coef
+
+  allocate(new_angles(mo_num))
+  new_angles(1:mo_num) = dabs(angle_left_right(1:mo_num))
+  max_angle = maxval(new_angles)
+  good_angles = max_angle.lt.45.d0
+  print *, ' max_angle = ', max_angle
+  
+end
+
+! ---
+
+subroutine build_s_matrix(m, n, C1, C2, overlap, smat)
+
+  implicit none
+  integer,          intent(in)  :: m, n
+  double precision, intent(in)  :: C1(m,n), C2(m,n), overlap(m,m)
+  double precision, intent(out) :: smat(n,n)
+  integer                       :: i, j, k, l
+  double precision, allocatable :: S_tmp(:,:)
+
+  smat = 0.d0
+
+  !do i = 1, n
+  !  do j = 1, n
+  !    do k = 1, m
+  !      do l = 1, m
+  !        smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j) 
+  !      enddo
+  !    enddo
+  !  enddo
+  !enddo
+
+  ! C1.T x overlap
+  allocate(S_tmp(n,m))
+  call dgemm( 'T', 'N', n, m, m, 1.d0                    &
+            , C1, size(C1, 1), overlap, size(overlap, 1) &
+            , 0.d0, S_tmp, size(S_tmp, 1) )
+  ! C1.T x overlap x C2
+  call dgemm( 'N', 'N', n, n, m, 1.d0                     &
+            , S_tmp, size(S_tmp, 1), C2(1,1), size(C2, 1) &
+            , 0.d0, smat, size(smat, 1) )
+  deallocate(S_tmp)
+
+end
+
+! ---
+
+subroutine orthog_functions(m, n, coef, overlap)
+
+  implicit none
+
+  integer,          intent(in)    :: m, n
+  double precision, intent(in)    :: overlap(m,m)
+  double precision, intent(inout) :: coef(m,n)
+  double precision, allocatable   :: stmp(:,:)
+  integer                         :: j, k
+
+  allocate(stmp(n,n))
+  call build_s_matrix(m, n, coef, coef, overlap, stmp)
+! print*,'overlap before'
+! do j = 1, n
+!  write(*,'(100(F16.10,X))')stmp(:,j)
+! enddo
+  call impose_orthog_svd_overlap(m, n, coef, overlap)
+  call build_s_matrix(m, n, coef, coef, overlap, stmp)
+  do j = 1, n
+    ! ---
+    ! TODO: MANU check ici
+    !coef(1,:m) *= 1.d0/dsqrt(stmp(j,j))
+    do k = 1, m
+      coef(k,j) *= 1.d0/dsqrt(stmp(j,j))
+    enddo
+    ! ---
+  enddo
+  call build_s_matrix(m, n, coef, coef, overlap, stmp)
+
+ !print*,'overlap after'
+ !do j = 1, n
+ ! write(*,'(100(F16.10,X))')stmp(:,j)
+ !enddo
+
+ deallocate(stmp)
+
+end
+
+! ---
+
+subroutine print_angles_tc()
+
+  implicit none
+  integer          :: i, j
+  double precision :: left, right
+
+  print *, ' product of norms, angle between vectors'                                                                  
+  do i = 1, mo_num
+    left  = overlap_mo_l(i,i)
+    right = overlap_mo_r(i,i)
+!  print*,Fock_matrix_tc_mo_tot(i,i),left*right,angle_left_right(i)
+    print *, left*right, angle_left_right(i)
+  enddo
+
+end
+
+! ---
+
+subroutine print_energy_and_mos()
+
+  implicit none
+  integer :: i
+
+  print *, ' '
+  print *, ' TC energy = ', TC_HF_energy
+  print *, ' TC SCF energy gradient = ', grad_non_hermit
+  print *, ' Max angle Left/right   = ', max_angle_left_right
+
+  if(max_angle_left_right .lt. 45.d0) then
+    print *, ' Maximum angle BELOW 45 degrees, everthing is OK !'
+  else if(max_angle_left_right .gt. 45.d0 .and. max_angle_left_right .lt. 75.d0) then
+    print *, ' Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...'
+  else if(max_angle_left_right .gt. 75.d0) then
+    print *, ' Maximum angle between ABOVE 75 degrees, YOU WILL CERTAINLY FIND TROUBLES IN TC-CI calculations ...'
+  endif
+
+  print *, ' Diag Fock elem, product of left/right norm, angle left/right '
+  do i = 1, mo_num
+    write(*, '(I3,X,100(F16.10,X))') i, Fock_matrix_tc_mo_tot(i,i), overlap_mo_l(i,i)*overlap_mo_r(i,i), angle_left_right(i)
+  enddo
+
+end
+
+! ---
+
+subroutine sort_by_tc_fock
+ implicit none 
+ integer, allocatable :: iorder(:)
+ double precision, allocatable :: mo_l_tmp(:,:), mo_r_tmp(:,:),fock(:)
+ allocate(iorder(mo_num),fock(mo_num),mo_l_tmp(ao_num, mo_num),mo_r_tmp(ao_num,mo_num))
+ integer :: i
+ mo_l_tmp = mo_l_coef
+ mo_r_tmp = mo_r_coef
+ do i = 1, mo_num
+  iorder(i) = i
+  fock(i) = Fock_matrix_tc_mo_tot(i,i)
+ enddo
+ call dsort(fock,iorder,mo_num)
+ do i = 1, mo_num
+  mo_l_coef(1:ao_num,i) = mo_l_tmp(1:ao_num,iorder(i))
+  mo_r_coef(1:ao_num,i) = mo_r_tmp(1:ao_num,iorder(i))
+ enddo
+ touch mo_l_coef mo_r_coef
+
+end
+

src/tc_scf/tc_petermann_factor.irp.f

@@ -0,0 +1,78 @@
+
+! ---
+
+program tc_petermann_factor
+
+  BEGIN_DOC
+  ! TODO : Put the documentation of the program here
+  END_DOC
+
+  implicit none
+
+  my_grid_becke  = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+!  my_n_pt_r_grid = 10 ! small grid for quick debug
+!  my_n_pt_a_grid = 26 ! small grid for quick debug
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  call main()
+
+end
+
+! ---
+
+subroutine main()
+
+  implicit none
+  integer                       :: i, j
+  double precision              :: Pf_diag_av
+  double precision, allocatable :: Sl(:,:), Sr(:,:), Pf(:,:)
+
+  allocate(Sl(mo_num,mo_num), Sr(mo_num,mo_num), Pf(mo_num,mo_num))
+
+  call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0                       &
+            , mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
+            , 0.d0, Sl, size(Sl, 1) )
+
+  print *, ''
+  print *, ' left-orthog matrix:'
+  do i = 1, mo_num
+    write(*,'(100(F8.4,X))') Sl(:,i)
+  enddo
+
+  call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0                       &
+            , mo_r_coef, size(mo_r_coef, 1), mo_r_coef, size(mo_r_coef, 1) &
+            , 0.d0, Sr, size(Sr, 1) )
+
+  print *, ''
+  print *, ' right-orthog matrix:'
+  do i = 1, mo_num
+    write(*,'(100(F8.4,X))') Sr(:,i)
+  enddo
+
+  print *, ''
+  print *, ' Petermann matrix:'
+  do i = 1, mo_num
+    do j = 1, mo_num
+      Pf(j,i) = Sl(j,i) * Sr(j,i)
+    enddo
+    write(*,'(100(F8.4,X))') Pf(:,i)
+  enddo
+
+  Pf_diag_av = 0.d0
+  do i = 1, mo_num
+    Pf_diag_av = Pf_diag_av + Pf(i,i)
+  enddo
+  Pf_diag_av = Pf_diag_av / dble(mo_num)
+
+  print *, ''
+  print *, ' mean of the diagonal Petermann factor = ', Pf_diag_av
+
+  deallocate(Sl, Sr, Pf)
+
+  return
+end subroutine
+
+! ---
+

src/tc_scf/tc_scf.irp.f

@@ -0,0 +1,75 @@
+! ---
+
+program tc_scf
+
+  BEGIN_DOC
+  ! TODO : Put the documentation of the program here
+  END_DOC
+
+  implicit none
+
+  print *, 'starting ...'
+
+  my_grid_becke  = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+!  my_n_pt_r_grid = 10 ! small grid for quick debug
+!  my_n_pt_a_grid = 26 ! small grid for quick debug
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  PROVIDE mu_erf 
+  print *, ' mu = ', mu_erf
+  PROVIDE j1b_type
+  print *, ' j1b_type = ', j1b_type
+  print *, j1b_pen
+
+  !call create_guess()
+  !call orthonormalize_mos()
+
+  PROVIDE tcscf_algorithm
+  if(tcscf_algorithm == 'DIIS') then
+    call rh_tcscf_diis()
+  elseif(tcscf_algorithm == 'Simple') then
+    call rh_tcscf_simple()
+  else
+    print *, ' not implemented yet', tcscf_algorithm
+    stop
+  endif
+
+  call minimize_tc_orb_angles()
+  call print_energy_and_mos()
+
+end
+
+! ---
+
+subroutine create_guess()
+
+  implicit none
+  logical :: exists
+
+  PROVIDE ezfio_filename
+  !call ezfio_has_mo_basis_mo_coef(exists)
+  exists = .false.
+
+  if(.not.exists) then
+    mo_label = 'Guess'
+    if(mo_guess_type == "HCore") then
+      mo_coef = ao_ortho_lowdin_coef
+      call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef, 1), 1.d-10)
+      TOUCH mo_coef
+      call mo_as_eigvectors_of_mo_matrix(mo_one_e_integrals, size(mo_one_e_integrals, 1), size(mo_one_e_integrals, 2), mo_label, 1, .false.)
+      call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef, 1), 1.d-10)
+      SOFT_TOUCH mo_coef
+    elseif (mo_guess_type == "Huckel") then
+      call huckel_guess
+    else
+      print *,  'Unrecognized MO guess type : '//mo_guess_type
+      stop 1
+    endif
+    SOFT_TOUCH mo_label
+  endif
+
+end subroutine create_guess
+
+! ---

src/tc_scf/tc_scf_dm.irp.f

@@ -0,0 +1,37 @@
+! ---
+
+BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_beta, (ao_num, ao_num) ]
+
+  implicit none
+
+  if(bi_ortho) then
+    PROVIDE mo_l_coef mo_r_coef
+    TCSCF_density_matrix_ao_beta = TCSCF_bi_ort_dm_ao_beta
+  else
+    TCSCF_density_matrix_ao_beta = SCF_density_matrix_ao_beta
+  endif
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_alpha, (ao_num, ao_num) ]
+
+  implicit none
+
+  if(bi_ortho) then
+    PROVIDE mo_l_coef mo_r_coef
+    TCSCF_density_matrix_ao_alpha = TCSCF_bi_ort_dm_ao_alpha
+  else
+    TCSCF_density_matrix_ao_alpha = SCF_density_matrix_ao_alpha
+  endif
+END_PROVIDER 
+
+
+! ---
+
+BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_tot, (ao_num, ao_num) ]
+  implicit none
+  TCSCF_density_matrix_ao_tot = TCSCF_density_matrix_ao_beta + TCSCF_density_matrix_ao_alpha
+END_PROVIDER 
+
+

src/tc_scf/tc_scf_energy.irp.f

@@ -0,0 +1,34 @@
+
+ BEGIN_PROVIDER [ double precision, TC_HF_energy]
+&BEGIN_PROVIDER [ double precision, TC_HF_one_e_energy]
+&BEGIN_PROVIDER [ double precision, TC_HF_two_e_energy]
+
+  BEGIN_DOC
+  ! TC Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components.
+  END_DOC
+
+  implicit none
+  integer :: i, j
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  TC_HF_energy = nuclear_repulsion
+  TC_HF_one_e_energy = 0.d0
+  TC_HF_two_e_energy = 0.d0
+
+  do j = 1, ao_num
+    do i = 1, ao_num
+      TC_HF_two_e_energy += 0.5d0 * ( two_e_tc_non_hermit_integral_alpha(i,j) * TCSCF_density_matrix_ao_alpha(i,j) &
+                                    + two_e_tc_non_hermit_integral_beta (i,j) * TCSCF_density_matrix_ao_beta (i,j) )
+      TC_HF_one_e_energy += ao_one_e_integrals_tc_tot(i,j) &
+                          * (TCSCF_density_matrix_ao_alpha(i,j) + TCSCF_density_matrix_ao_beta (i,j) )
+    enddo
+  enddo
+
+  TC_HF_energy += TC_HF_one_e_energy + TC_HF_two_e_energy
+  TC_HF_energy += diag_three_elem_hf
+
+END_PROVIDER
+
+! ---
+

src/tc_scf/tc_scf_utils.irp.f

@@ -0,0 +1,43 @@
+
+! ---
+
+subroutine LTxSxR(n, m, L, S, R, C)
+
+  implicit none
+  integer,          intent(in)  :: n, m
+  double precision, intent(in)  :: L(n,m), S(n,n), R(n,m)
+  double precision, intent(out) :: C(m,m)
+  integer                       :: i, j
+  double precision              :: accu_d, accu_nd
+  double precision, allocatable :: tmp(:,:)
+
+  ! L.T x S x R
+  allocate(tmp(m,n))
+  call dgemm( 'T', 'N', m, n, n, 1.d0      &
+            , L, size(L, 1), S, size(S, 1) &
+            , 0.d0, tmp, size(tmp, 1) )
+  call dgemm( 'N', 'N', m, m, n, 1.d0          &
+            , tmp, size(tmp, 1), R, size(R, 1) &
+            , 0.d0, C, size(C, 1) )
+  deallocate(tmp)
+
+  accu_d  = 0.d0
+  accu_nd = 0.d0
+  do i = 1, m
+    do j = 1, m
+      if(j.eq.i) then
+        accu_d += dabs(C(j,i))
+      else
+        accu_nd += C(j,i) * C(j,i)
+      endif
+    enddo
+  enddo
+  accu_nd = dsqrt(accu_nd)
+
+  print*, ' accu_d  = ', accu_d
+  print*, ' accu_nd = ', accu_nd
+
+end subroutine LTxR
+
+! ---
+

src/tc_scf/test_Ne.sh

@@ -0,0 +1,13 @@
+QP_ROOT=/home/eginer/new_qp2/qp2
+source ${QP_ROOT}/quantum_package.rc
+  echo Ne > Ne.xyz
+  echo $QP_ROOT
+  qp create_ezfio -b cc-pcvdz Ne.xyz 
+  qp run scf 
+  qp set tc_keywords bi_ortho True 
+  qp set ao_two_e_erf_ints mu_erf 0.87 
+  qp set tc_keywords j1b_pen [1.5]
+  qp set tc_keywords j1b_type 3 
+  qp run tc_scf | tee Ne.ezfio.tc_scf.out 
+  grep "TC energy =" Ne.ezfio.tc_scf.out | tail -1 
+  eref=-128.552134

src/tc_scf/three_e_energy_bi_ortho.irp.f

@@ -0,0 +1,174 @@
+
+subroutine contrib_3e_diag_sss(i,j,k,integral)
+ implicit none
+ integer, intent(in) :: i,j,k
+ BEGIN_DOC
+ ! returns the pure same spin contribution to diagonal matrix element of 3e term
+ END_DOC
+ double precision, intent(out) :: integral
+ double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+ call  give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int )!!! < i k j | i k j >
+ call  give_integrals_3_body_bi_ort(i, k, j, j, i, k, c_3_int)      ! < i k j | j i k >
+ call  give_integrals_3_body_bi_ort(i, k, j, k, j, i, c_minus_3_int)! < i k j | k j i >
+ integral = direct_int + c_3_int + c_minus_3_int 
+ ! negative terms :: exchange contrib
+ call  give_integrals_3_body_bi_ort(i, k, j, j, k, i, exch_13_int)!!! < i k j | j k i > : E_13 
+ call  give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)!!! < i k j | i j k > : E_23
+ call  give_integrals_3_body_bi_ort(i, k, j, k, i, j, exch_12_int)!!! < i k j | k i j > : E_12
+ integral += - exch_13_int - exch_23_int  - exch_12_int 
+ integral = -integral
+end
+
+subroutine contrib_3e_diag_soo(i,j,k,integral)
+ implicit none
+ integer, intent(in) :: i,j,k
+ BEGIN_DOC
+ ! returns the pure same spin contribution to diagonal matrix element of 3e term
+ END_DOC
+ double precision, intent(out) :: integral
+ double precision :: direct_int, exch_23_int
+ call  give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int) ! < i k j | i k j >
+ call  give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)! < i k j | i j k > : E_23
+ integral = direct_int - exch_23_int 
+ integral = -integral
+end
+
+
+subroutine give_aaa_contrib_bis(integral_aaa)
+ implicit none
+ double precision, intent(out) :: integral_aaa
+ double precision :: integral
+ integer :: i,j,k
+ integral_aaa = 0.d0
+ do i = 1, elec_alpha_num
+  do j = i+1, elec_alpha_num
+   do k = j+1, elec_alpha_num
+    call contrib_3e_diag_sss(i,j,k,integral)
+    integral_aaa += integral
+   enddo
+  enddo
+ enddo
+
+end
+
+subroutine give_aaa_contrib(integral_aaa)
+ implicit none
+ double precision, intent(out) :: integral_aaa
+ double precision :: integral
+ integer :: i,j,k
+ integral_aaa = 0.d0
+ do i = 1, elec_alpha_num
+  do j = 1, elec_alpha_num
+   do k = 1, elec_alpha_num
+    call contrib_3e_diag_sss(i,j,k,integral)
+    integral_aaa += integral
+   enddo
+  enddo
+ enddo
+ integral_aaa *= 1.d0/6.d0 
+end
+
+
+subroutine give_aab_contrib(integral_aab)
+ implicit none
+ double precision, intent(out) :: integral_aab
+ double precision :: integral
+ integer :: i,j,k
+ integral_aab = 0.d0
+ do i = 1, elec_beta_num
+  do j = 1, elec_alpha_num
+   do k = 1, elec_alpha_num
+    call contrib_3e_diag_soo(i,j,k,integral)
+    integral_aab += integral
+   enddo
+  enddo
+ enddo
+ integral_aab *= 0.5d0
+end
+
+
+subroutine give_aab_contrib_bis(integral_aab)
+ implicit none
+ double precision, intent(out) :: integral_aab
+ double precision :: integral
+ integer :: i,j,k
+ integral_aab = 0.d0
+ do i = 1, elec_beta_num
+  do j = 1, elec_alpha_num
+   do k = j+1, elec_alpha_num
+    call contrib_3e_diag_soo(i,j,k,integral)
+    integral_aab += integral
+   enddo
+  enddo
+ enddo
+end
+
+
+subroutine give_abb_contrib(integral_abb)
+ implicit none
+ double precision, intent(out) :: integral_abb
+ double precision :: integral
+ integer :: i,j,k
+ integral_abb = 0.d0
+ do i = 1, elec_alpha_num
+  do j = 1, elec_beta_num
+   do k = 1, elec_beta_num
+    call contrib_3e_diag_soo(i,j,k,integral)
+    integral_abb += integral
+   enddo
+  enddo
+ enddo
+ integral_abb *= 0.5d0
+end
+
+subroutine give_abb_contrib_bis(integral_abb)
+ implicit none
+ double precision, intent(out) :: integral_abb
+ double precision :: integral
+ integer :: i,j,k
+ integral_abb = 0.d0
+ do i = 1, elec_alpha_num
+  do j = 1, elec_beta_num
+   do k = j+1, elec_beta_num
+    call contrib_3e_diag_soo(i,j,k,integral)
+    integral_abb += integral
+   enddo
+  enddo
+ enddo
+end
+
+subroutine give_bbb_contrib_bis(integral_bbb)
+ implicit none
+ double precision, intent(out) :: integral_bbb
+ double precision :: integral
+ integer :: i,j,k
+ integral_bbb = 0.d0
+ do i = 1, elec_beta_num
+  do j = i+1, elec_beta_num
+   do k = j+1, elec_beta_num
+    call contrib_3e_diag_sss(i,j,k,integral)
+    integral_bbb += integral
+   enddo
+  enddo
+ enddo
+
+end
+
+subroutine give_bbb_contrib(integral_bbb)
+ implicit none
+ double precision, intent(out) :: integral_bbb
+ double precision :: integral
+ integer :: i,j,k
+ integral_bbb = 0.d0
+ do i = 1, elec_beta_num
+  do j = 1, elec_beta_num
+   do k = 1, elec_beta_num
+    call contrib_3e_diag_sss(i,j,k,integral)
+    integral_bbb += integral
+   enddo
+  enddo
+ enddo
+ integral_bbb *= 1.d0/6.d0 
+end
+
+

src/utils/block_diag_degen.irp.f

@@ -0,0 +1,218 @@
+
+subroutine diag_mat_per_fock_degen(fock_diag, mat_ref, n, thr_d, thr_nd, thr_deg, leigvec, reigvec, eigval)
+
+
+  BEGIN_DOC
+  !
+  ! subroutine that diagonalizes a matrix mat_ref BY BLOCK
+  !
+  ! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag"
+  !
+  ! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together
+  !
+  !          : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together
+  !
+  !          : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together
+  !
+  ! etc... the advantage is to guarentee no spurious mixing because of numerical problems. 
+  !
+  END_DOC
+
+  implicit none
+  integer,          intent(in)  :: n
+  double precision, intent(in)  :: fock_diag(n), mat_ref(n,n), thr_d, thr_nd, thr_deg
+  double precision, intent(out) :: leigvec(n,n), reigvec(n,n), eigval(n)
+
+  integer                       :: n_degen_list, n_degen,size_mat, i, j, k, icount, m, index_degen
+  integer                       :: ii, jj, i_good, j_good, n_real
+  integer                       :: icount_eigval
+  logical,          allocatable :: is_ok(:)
+  integer,          allocatable :: list_degen(:,:), list_same_degen(:)
+  integer,          allocatable :: iorder(:), list_degen_sorted(:)
+  double precision, allocatable :: leigvec_unsrtd(:,:), reigvec_unsrtd(:,:), eigval_unsrtd(:)
+  double precision, allocatable :: mat_tmp(:,:), eigval_tmp(:), leigvec_tmp(:,:), reigvec_tmp(:,:)
+
+  allocate(leigvec_unsrtd(n,n), reigvec_unsrtd(n,n), eigval_unsrtd(n))
+  leigvec_unsrtd = 0.d0
+  reigvec_unsrtd = 0.d0
+  eigval_unsrtd  = 0.d0
+
+  ! obtain degeneracies 
+  allocate(list_degen(n,0:n))
+  call give_degen_full_list(fock_diag, n, thr_deg, list_degen, n_degen_list)
+
+  allocate(iorder(n_degen_list), list_degen_sorted(n_degen_list))
+  do i = 1, n_degen_list
+    n_degen = list_degen(i,0) 
+    list_degen_sorted(i) = n_degen
+    iorder(i) = i
+  enddo
+
+  ! sort by number of degeneracies 
+  call isort(list_degen_sorted, iorder, n_degen_list)
+
+  allocate(is_ok(n_degen_list))
+  is_ok = .True.
+  icount_eigval = 0
+
+  ! loop over degeneracies 
+  do i = 1, n_degen_list
+    if(.not.is_ok(i)) cycle
+
+    is_ok(i) = .False.
+    n_degen  = list_degen_sorted(i)
+
+    print *, ' diagonalizing for n_degen = ', n_degen
+
+    k = 1
+
+   ! group all the entries having the same degeneracies 
+!!  do while (list_degen_sorted(i+k)==n_degen)
+    do m = i+1, n_degen_list
+      if(list_degen_sorted(m)==n_degen) then
+        is_ok(i+k) = .False.
+        k += 1
+      endif
+    enddo
+
+    print *, ' number of identical degeneracies = ', k
+    size_mat = k*n_degen  
+    print *, ' size_mat = ', size_mat
+    allocate(mat_tmp(size_mat,size_mat), list_same_degen(size_mat))
+    allocate(eigval_tmp(size_mat), leigvec_tmp(size_mat,size_mat), reigvec_tmp(size_mat,size_mat))
+    ! group all the elements sharing the same degeneracy
+    icount = 0
+    do j = 1, k ! jth set of degeneracy
+      index_degen = iorder(i+j-1)
+      do m = 1, n_degen
+        icount += 1
+        list_same_degen(icount) = list_degen(index_degen,m)
+      enddo
+    enddo
+
+    print *, ' list of elements '
+    do icount = 1, size_mat
+      print *, icount, list_same_degen(icount)
+    enddo
+
+    ! you copy subset of matrix elements having all the same degeneracy in mat_tmp
+    do ii = 1, size_mat
+      i_good = list_same_degen(ii)
+      do jj = 1, size_mat
+        j_good = list_same_degen(jj)
+        mat_tmp(jj,ii) = mat_ref(j_good,i_good)
+      enddo
+    enddo
+
+    call non_hrmt_bieig( size_mat, mat_tmp, thr_d, thr_nd &
+                       , leigvec_tmp, reigvec_tmp         & 
+                       , n_real, eigval_tmp )
+
+    do ii = 1, size_mat
+      icount_eigval += 1
+      eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues 
+      do jj = 1, size_mat ! copy the eigenvectors 
+        j_good = list_same_degen(jj)
+        leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii)
+        reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii)
+      enddo
+    enddo
+
+    deallocate(mat_tmp, list_same_degen)
+    deallocate(eigval_tmp, leigvec_tmp, reigvec_tmp)
+  enddo
+
+  if(icount_eigval .ne. n) then
+    print *, ' pb !! (icount_eigval.ne.n)'
+    print *, ' icount_eigval,n', icount_eigval, n
+    stop
+  endif
+ 
+  deallocate(iorder)
+  allocate(iorder(n))
+  do i = 1, n
+    iorder(i) = i
+  enddo
+  call dsort(eigval_unsrtd, iorder, n)
+
+  do i = 1, n
+    print*,'sorted eigenvalues '
+    i_good = iorder(i)
+    eigval(i) = eigval_unsrtd(i)
+    print*,'i,eigval(i) = ',i,eigval(i)
+    do j = 1, n
+      leigvec(j,i) = leigvec_unsrtd(j,i_good)
+      reigvec(j,i) = reigvec_unsrtd(j,i_good)
+    enddo
+  enddo
+
+  deallocate(leigvec_unsrtd, reigvec_unsrtd, eigval_unsrtd)
+  deallocate(list_degen)
+  deallocate(iorder, list_degen_sorted)
+  deallocate(is_ok)
+
+end
+
+! ---
+
+subroutine give_degen_full_list(A, n, thr, list_degen, n_degen_list)
+
+  BEGIN_DOC
+  ! you enter with an array A(n) and spits out all the elements degenerated up to thr
+  !
+  ! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL 
+  !
+  ! list_degen(i,0) = number of degenerate entries 
+  !
+  ! list_degen(i,1) = index of the first degenerate entry
+  !
+  ! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries 
+  !
+  ! if list_degen(i,0) == 1 it means that there is no degeneracy for that element
+  END_DOC
+
+  implicit none
+
+  double precision, intent(in)  :: A(n)
+  double precision, intent(in)  :: thr
+  integer,          intent(in)  :: n
+  integer,          intent(out) :: list_degen(n,0:n), n_degen_list
+  integer                       :: i, j, icount, icheck
+  logical, allocatable          :: is_ok(:)
+
+
+  allocate(is_ok(n))
+  n_degen_list = 0
+  is_ok = .True.
+  do i = 1, n
+    if(.not.is_ok(i)) cycle
+    n_degen_list +=1
+    is_ok(i) = .False.
+    list_degen(n_degen_list,1) = i
+    icount = 1
+    do j = i+1, n
+      if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then
+        is_ok(j) = .False.
+        icount += 1
+        list_degen(n_degen_list,icount) = j
+      endif
+    enddo
+
+    list_degen(n_degen_list,0) = icount
+  enddo
+
+  icheck = 0
+  do i = 1, n_degen_list
+    icheck += list_degen(i,0)
+  enddo
+
+  if(icheck.ne.n)then
+    print *, ' pb ! :: icheck.ne.n'
+    print *, icheck, n
+    stop
+  endif
+
+end
+
+! ---
+

src/utils/loc.f

@@ -0,0 +1,327 @@
+c************************************************************************
+      subroutine maxovl(n,m,s,t,w)
+C
+C     This subprogram contains an iterative procedure to find the
+C     unitary transformation of a set of n vectors which maximizes
+C     the sum of their square overlaps with a set of m reference
+C     vectors (m.le.n)
+C
+C     S: overlap matrix <ref|vec>
+C     T: rotation matrix
+C     W: new overlap matrix
+C
+C
+      implicit real*8(a-h,o-y),logical*1(z)
+!      parameter (id1=700)
+!      dimension s(id1,id1),t(id1,id1),w(id1,id1)
+      double precision, intent(inout) :: s(n,n)
+      double precision, intent(out) :: t(n,n),w(n,n)
+      data small/1.d-6/
+
+      zprt=.true.
+      niter=1000000
+      conv=1.d-12
+
+C      niter=1000000
+C      conv=1.d-6
+      write (6,5) n,m,conv
+    5 format (//5x,'Unitary transformation of',i3,'  vectors'/
+     * 5x,'following the principle of maximum overlap with a set of',
+     * i3,' reference vectors'/5x,'required convergence on rotation ',
+     * 'angle =',f13.10///5x,'Starting overlap matrix'/)
+      do i=1,m
+        write (6,145) i
+        write (6,150) (s(i,j),j=1,n)
+      end do
+    8 mm=m-1
+      if (m.lt.n) mm=m
+      iter=0
+      do j=1,n
+        do i=1,n
+          t(i,j)=0.d0
+        end do
+        do i=1,m
+          w(i,j)=s(i,j)
+        enddo
+        t(j,j)=1.d0
+      enddo
+      sum=0.d0
+      do i=1,m
+        sum=sum+s(i,i)*s(i,i)
+      end do
+      sum=sum/m
+      if (zprt) write (6,12) sum
+   12 format (//5x,'Average square overlap =',f10.6)
+      if (n.eq.1) goto 100
+      last=n
+      j=1
+   21 if (j.ge.last) goto 30
+      sum=0.d0
+      do i=1,n
+        sum=sum+s(i,j)*s(i,j)
+      enddo
+      if (sum.gt.small) goto 28
+      do i=1,n
+        sij=s(i,j)
+        s(i,j)=-s(i,last)
+        s(i,last)=sij
+        tij=t(i,j)
+        t(i,j)=-t(i,last)
+        t(i,last)=tij
+      end do
+      last=last-1
+      goto 21
+   28 j=j+1
+      goto 21
+   30 iter=iter+1
+      imax=0
+      jmax=0
+      dmax=0.d0
+      amax=0.d0
+      do 60 i=1,mm
+      ip=i+1
+      do 50 j=ip,n
+      a=s(i,j)*s(i,j)-s(i,i)*s(i,i)
+      b=-s(i,i)*s(i,j)
+      if (j.gt.m) goto 31
+      a=a+s(j,i)*s(j,i)-s(j,j)*s(j,j)
+      b=b+s(j,i)*s(j,j)
+   31 b=b+b
+      if (a.eq.0.d0) goto 32
+      ba=b/a
+      if (dabs(ba).gt.small) goto 32
+      if (a.gt.0.d0) goto 33
+      tang=-0.5d0*ba
+      cosine=1.d0/dsqrt(1.d0+tang*tang)
+      sine=tang*cosine
+      goto 34
+   32 tang=0.d0
+      if (b.ne.0.d0) tang=(a+dsqrt(a*a+b*b))/b
+      cosine=1.d0/dsqrt(1.d0+tang*tang)
+      sine=tang*cosine
+      goto 34
+   33 cosine=0.d0
+      sine=1.d0
+   34 delta=sine*(a*sine+b*cosine)
+      if (zprt.and.delta.lt.0.d0) write (6,71) i,j,a,b,sine,cosine,delta
+      do k=1,m
+        p=s(k,i)*cosine-s(k,j)*sine
+        q=s(k,i)*sine+s(k,j)*cosine
+        s(k,i)=p
+        s(k,j)=q
+      enddo
+      do k=1,n
+        p=t(k,i)*cosine-t(k,j)*sine
+        q=t(k,i)*sine+t(k,j)*cosine
+        t(k,i)=p
+        t(k,j)=q
+      enddo
+   45 d=dabs(sine)
+      if (d.le.amax) goto 50
+      imax=i
+      jmax=j
+      amax=d
+      dmax=delta
+   50 continue
+   60 continue
+      if (zprt) write (6,70) iter,amax,imax,jmax,dmax
+   70 format (' iter=',i4,' largest rotation=',f12.8,
+     * ', vectors',i3,' and',i3,', incr. of diag. squares=',g12.5)
+   71 format (' i,j,a,b,sin,cos,delta =',2i3,5f10.5)
+      if (amax.lt.conv) goto 100
+      if (iter.lt.niter) goto 30
+      write (6,80)
+      write (6,*) 'niter=',niter
+   80 format (//5x,'*** maximum number of cycles exceeded ',
+     * 'in subroutine maxovl ***'//)
+      stop
+  100 continue
+      do j=1,n
+        if (s(j,j).gt.0.d0) cycle
+        do i=1,m
+          s(i,j)=-s(i,j)
+        enddo
+        do i=1,n
+          t(i,j)=-t(i,j)
+        enddo
+      enddo
+      sum=0.d0
+      do i=1,m
+        sum=sum+s(i,i)*s(i,i)
+      enddo
+      sum=sum/m
+      do i=1,m
+        do j=1,n
+          sw=s(i,j)
+          s(i,j)=w(i,j)
+          w(i,j)=sw
+        enddo
+      enddo
+      if (.not.zprt) return
+      write (6,12) sum
+      write (6,130)
+  130 format (//5x,'transformation matrix')
+      do i=1,n
+        write (6,145) i
+        write (6,150) (t(i,j),j=1,n)
+      enddo
+  145 format (i8)
+  150 format (2x,10f12.8)
+      write (6,160)
+  160 format (//5x,'new overlap matrix'/)
+      do i=1,m
+        write (6,145) i
+        write (6,150) (w(i,j),j=1,n)
+      enddo
+      return
+      end
+
+
+c************************************************************************
+      subroutine maxovl_no_print(n,m,s,t,w)
+C
+C     This subprogram contains an iterative procedure to find the
+C     unitary transformation of a set of n vectors which maximizes
+C     the sum of their square overlaps with a set of m reference
+C     vectors (m.le.n)
+C
+C     S: overlap matrix <ref|vec>
+C     T: rotation matrix
+C     W: new overlap matrix
+C
+C
+      implicit real*8(a-h,o-y),logical*1(z)
+      parameter (id1=300)
+      dimension s(id1,id1),t(id1,id1),w(id1,id1)
+      data small/1.d-6/
+
+      zprt=.false.
+      niter=1000000
+      conv=1.d-8
+
+C      niter=1000000
+C      conv=1.d-6
+    8 mm=m-1
+      if (m.lt.n) mm=m
+      iter=0
+      do j=1,n
+        do i=1,n
+          t(i,j)=0.d0
+        enddo
+        do i=1,m
+          w(i,j)=s(i,j)
+        enddo
+        t(j,j)=1.d0
+      enddo
+      sum=0.d0
+      do i=1,m
+        sum=sum+s(i,i)*s(i,i)
+      enddo
+      sum=sum/m
+   12 format (//5x,'Average square overlap =',f10.6)
+      if (n.eq.1) goto 100
+      last=n
+      j=1
+   21 if (j.ge.last) goto 30
+      sum=0.d0
+
+      do i=1,n
+        sum=sum+s(i,j)*s(i,j)
+      enddo
+      if (sum.gt.small) goto 28
+      do i=1,n
+        sij=s(i,j)
+        s(i,j)=-s(i,last)
+        s(i,last)=sij
+        tij=t(i,j)
+        t(i,j)=-t(i,last)
+        t(i,last)=tij
+      end do
+      last=last-1
+      goto 21
+   28 j=j+1
+      goto 21
+   30 iter=iter+1
+      imax=0
+      jmax=0
+      dmax=0.d0
+      amax=0.d0
+      do i=1,mm
+        ip=i+1
+        do j=ip,n
+          a=s(i,j)*s(i,j)-s(i,i)*s(i,i)
+          b=-s(i,i)*s(i,j)
+          if (j.gt.m) goto 31
+          a=a+s(j,i)*s(j,i)-s(j,j)*s(j,j)
+          b=b+s(j,i)*s(j,j)
+  31      b=b+b
+          if (a.eq.0.d0) goto 32
+          ba=b/a
+          if (dabs(ba).gt.small) goto 32
+          if (a.gt.0.d0) goto 33
+          tang=-0.5d0*ba
+          cosine=1.d0/dsqrt(1.d0+tang*tang)
+          sine=tang*cosine
+          goto 34
+  32      tang=0.d0
+          if (b.ne.0.d0) tang=(a+dsqrt(a*a+b*b))/b
+          cosine=1.d0/dsqrt(1.d0+tang*tang)
+          sine=tang*cosine
+          goto 34
+  33      cosine=0.d0
+          sine=1.d0
+  34      delta=sine*(a*sine+b*cosine)
+          do k=1,m
+            p=s(k,i)*cosine-s(k,j)*sine
+            q=s(k,i)*sine+s(k,j)*cosine
+            s(k,i)=p
+            s(k,j)=q
+          enddo
+          do k=1,n
+            p=t(k,i)*cosine-t(k,j)*sine
+            q=t(k,i)*sine+t(k,j)*cosine
+            t(k,i)=p
+            t(k,j)=q
+          enddo
+  45      d=dabs(sine)
+          if (d.le.amax) goto 50
+          imax=i
+          jmax=j
+          amax=d
+          dmax=delta
+  50      continue
+        end do
+      end do
+   70 format (' iter=',i4,' largest rotation=',f12.8,
+     * ', vectors',i3,' and',i3,', incr. of diag. squares=',g12.5)
+   71 format (' i,j,a,b,sin,cos,delta =',2i3,5f10.5)
+      if (amax.lt.conv) goto 100
+      if (iter.lt.niter) goto 30
+   80 format (//5x,'*** maximum number of cycles exceeded ',
+     * 'in subroutine maxovl ***'//)
+      stop
+  100 continue
+      do j=1,n
+        if (s(j,j).gt.0.d0) cycle
+        do i=1,m
+          s(i,j)=-s(i,j)
+        enddo
+        do i=1,n
+          t(i,j)=-t(i,j)
+        enddo
+      enddo
+      sum=0.d0
+      do i=1,m
+        sum=sum+s(i,i)*s(i,i)
+      enddo
+      sum=sum/m
+      do i=1,m
+        do j=1,n
+          sw=s(i,j)
+          s(i,j)=w(i,j)
+          w(i,j)=sw
+        enddo
+      enddo
+      return
+      end
+