Merge pull request #280 from AbdAmmar/dev-stable-tc-scf

Dev stable tc scf

src/non_h_ints_mu/jast_deriv.irp.f

@@ -272,9 +272,9 @@ subroutine grad1_j1b_nucl(r, grad)
       fact_z += e * z
     enddo
 
-    grad(1) = -2.d0 * fact_x
+    grad(1) = 2.d0 * fact_x
-    grad(2) = -2.d0 * fact_y
+    grad(2) = 2.d0 * fact_y
-    grad(3) = -2.d0 * fact_z
+    grad(3) = 2.d0 * fact_z
 
   else if(j1b_type .eq. 105) then
 
@@ -294,9 +294,9 @@ subroutine grad1_j1b_nucl(r, grad)
       fact_z += e * z
     enddo
 
-    grad(1) = -4.d0 * fact_x
+    grad(1) = 4.d0 * fact_x
-    grad(2) = -4.d0 * fact_y
+    grad(2) = 4.d0 * fact_y
-    grad(3) = -4.d0 * fact_z
+    grad(3) = 4.d0 * fact_z
 
   else
 

src/non_h_ints_mu/tc_integ.irp.f

@@ -46,7 +46,12 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
       PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b
 
       int2_grad1_u12_ao = 0.d0
-      ! TODO OPENMP
+      !$OMP PARALLEL                                                                 &
+      !$OMP DEFAULT (NONE)                                                           &
+      !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp0, tmp1, tmp2, tmp_x, tmp_y, tmp_z)   &
+      !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points, v_1b, v_1b_grad &
+      !$OMP        , v_ij_erf_rk_cst_mu_j1b, v_ij_u_cst_mu_j1b, x_v_ij_erf_rk_cst_mu_j1b, int2_grad1_u12_ao)
+      !$OMP DO SCHEDULE (static)
       do ipoint = 1, n_points_final_grid
         x     = final_grid_points(1,ipoint)
         y     = final_grid_points(2,ipoint)
@@ -65,6 +70,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
           enddo
         enddo
       enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
 
     endif
 

src/tc_bi_ortho/print_tc_var.irp.f

@@ -0,0 +1,20 @@
+program print_tc_var
+
+  BEGIN_DOC
+  ! TODO : Put the documentation of the program here
+  END_DOC
+
+  implicit none
+
+  print *, 'Hello world'
+  my_grid_becke = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  read_wf = .True.
+  touch read_wf
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  call write_tc_var()
+
+end
+

src/tc_bi_ortho/tc_utils.irp.f

@@ -32,3 +32,29 @@ subroutine write_tc_energy()
 
 end
 
+! ---
+
+subroutine write_tc_var()
+
+  implicit none
+  integer          :: i, j, k
+  double precision :: hmono, htwoe, hthree, htot
+  double precision :: SIGMA_TC
+
+  do k = 1, n_states
+
+    SIGMA_TC = 0.d0
+    do j = 2, N_det
+      call htilde_mu_mat_bi_ortho(psi_det(1,1,1), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
+      SIGMA_TC = SIGMA_TC + htot * htot
+    enddo
+
+    print *, " state    : ", k
+    print *, " SIGMA_TC = ", SIGMA_TC
+
+  enddo
+
+end
+
+! ---
+

src/tc_scf/diis_tcscf.irp.f

@@ -87,9 +87,13 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)]
 
   implicit none
   integer                       :: i, j
+  double precision              :: t0, t1
   double precision, allocatable :: tmp(:,:)
   double precision, allocatable :: F(:,:)
 
+  !print *, ' Providing FQS_SQF_ao ...'
+  !call wall_time(t0)
+
   allocate(F(ao_num,ao_num))
   if(var_tc) then
 
@@ -136,6 +140,9 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)]
   deallocate(tmp)
   deallocate(F)
 
+  !call wall_time(t1)
+  !print *, ' Wall time for FQS_SQF_ao =', t1-t0
+
 END_PROVIDER
 
 ! ---
@@ -143,6 +150,10 @@ END_PROVIDER
 BEGIN_PROVIDER [double precision, FQS_SQF_mo, (mo_num, mo_num)]
 
   implicit none
+  double precision :: t0, t1
+
+  !print*, ' Providing FQS_SQF_mo ...'
+  !call wall_time(t0)
 
   PROVIDE mo_r_coef mo_l_coef
   PROVIDE FQS_SQF_ao
@@ -150,6 +161,9 @@ BEGIN_PROVIDER [double precision, FQS_SQF_mo, (mo_num, mo_num)]
   call ao_to_mo_bi_ortho( FQS_SQF_ao, size(FQS_SQF_ao, 1) &
                         , FQS_SQF_mo, size(FQS_SQF_mo, 1) )
 
+  !call wall_time(t1)
+  !print*, ' Wall time for FQS_SQF_mo =', t1-t0 
+
 END_PROVIDER
 
 ! ---

src/tc_scf/fock_3e_bi_ortho_uhf.irp.f

@@ -7,14 +7,24 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
   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
+  double precision, allocatable :: tmp(:,:)
 
   PROVIDE mo_l_coef mo_r_coef
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, I_bij_aij)
 
   !print *, ' PROVIDING fock_3e_uhf_mo_cs ...'
-  call wall_time(ti)
+  !call wall_time(ti)
 
   fock_3e_uhf_mo_cs = 0.d0
 
+  !$OMP PARALLEL DEFAULT (NONE)                                                                     &
+  !$OMP PRIVATE (a, b, i, j, I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia, tmp) &
+  !$OMP SHARED  (mo_num, elec_beta_num, fock_3e_uhf_mo_cs)
+
+  allocate(tmp(mo_num,mo_num))
+  tmp = 0.d0
+
+  !$OMP DO
   do a = 1, mo_num
     do b = 1, mo_num
    
@@ -28,7 +38,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
           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 &
+          tmp(b,a) -= 0.5d0 * ( 4.d0 * I_bij_aij &
                               +        I_bij_ija &
                               +        I_bij_jai &
                               - 2.d0 * I_bij_aji &
@@ -39,8 +49,20 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
       enddo
     enddo
   enddo
+  !$OMP END DO NOWAIT
 
-  call wall_time(tf)
+  !$OMP CRITICAL
+  do a = 1, mo_num
+    do b = 1, mo_num
+      fock_3e_uhf_mo_cs(b,a) += tmp(b,a)
+    enddo
+  enddo
+  !$OMP END CRITICAL
+
+  deallocate(tmp)
+  !$OMP END PARALLEL
+
+  !call wall_time(tf)
   !print *, ' total Wall time for fock_3e_uhf_mo_cs =', tf - ti
 
 END_PROVIDER 
@@ -61,17 +83,27 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
   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
+  double precision, allocatable :: tmp(:,:)
 
   PROVIDE mo_l_coef mo_r_coef
   PROVIDE fock_3e_uhf_mo_cs
 
-  !print *, ' PROVIDING fock_3e_uhf_mo_a ...'
+  !print *, ' Providing fock_3e_uhf_mo_a ...'
   !call wall_time(ti)
 
   o = elec_beta_num + 1
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, I_bij_aij)
 
   fock_3e_uhf_mo_a = fock_3e_uhf_mo_cs
 
+  !$OMP PARALLEL DEFAULT (NONE)                                                                     &
+  !$OMP PRIVATE (a, b, i, j, I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia, tmp) &
+  !$OMP SHARED  (mo_num, o, elec_alpha_num, elec_beta_num, fock_3e_uhf_mo_a)
+
+  allocate(tmp(mo_num,mo_num))
+  tmp = 0.d0
+
+  !$OMP DO
   do a = 1, mo_num
     do b = 1, mo_num
 
@@ -87,7 +119,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
           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 &
+          tmp(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
                               +        I_bij_ija &
                               +        I_bij_jai &
                               -        I_bij_aji &
@@ -109,7 +141,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
           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 &
+          tmp(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
                               +        I_bij_ija &
                               +        I_bij_jai &
                               -        I_bij_aji &
@@ -131,7 +163,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
           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 &
+          tmp(b,a) -= 0.5d0 * ( I_bij_aij &
                               + I_bij_ija &
                               + I_bij_jai &
                               - I_bij_aji &
@@ -145,35 +177,58 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
 
     enddo
   enddo
+  !$OMP END DO NOWAIT
+
+  !$OMP CRITICAL
+  do a = 1, mo_num
+    do b = 1, mo_num
+      fock_3e_uhf_mo_a(b,a) += tmp(b,a)
+    enddo
+  enddo
+  !$OMP END CRITICAL
+
+  deallocate(tmp)
+  !$OMP END PARALLEL
 
   !call wall_time(tf)
-  !print *, ' total Wall time for fock_3e_uhf_mo_a =', tf - ti
+  !print *, ' 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)]
+
   BEGIN_DOC
-! BETA part of the Fock matrix from three-electron terms 
+  ! BETA part of the Fock matrix from three-electron terms 
-!
+  !
-! WARNING :: non hermitian if bi-ortho MOS used 
+  ! WARNING :: non hermitian if bi-ortho MOS used 
   END_DOC
 
   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
+  double precision, allocatable :: tmp(:,:)
 
   PROVIDE mo_l_coef mo_r_coef
 
   !print *, ' PROVIDING fock_3e_uhf_mo_b ...'
-  call wall_time(ti)
+  !call wall_time(ti)
 
   o = elec_beta_num + 1
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, I_bij_aij)
 
   fock_3e_uhf_mo_b = fock_3e_uhf_mo_cs
 
+  !$OMP PARALLEL DEFAULT (NONE)                                                                     &
+  !$OMP PRIVATE (a, b, i, j, I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia, tmp) &
+  !$OMP SHARED  (mo_num, o, elec_alpha_num, elec_beta_num, fock_3e_uhf_mo_b)
+
+  allocate(tmp(mo_num,mo_num))
+  tmp = 0.d0
+
+  !$OMP DO
   do a = 1, mo_num
     do b = 1, mo_num
 
@@ -189,7 +244,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)]
           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 &
+          tmp(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
                               -        I_bij_aji &
                               -        I_bij_iaj )
 
@@ -208,7 +263,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)]
           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 &
+          tmp(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij &
                               -        I_bij_aji &
                               -        I_bij_jia )
 
@@ -227,7 +282,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)]
           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 &
+          tmp(b,a) -= 0.5d0 * ( I_bij_aij &
                               - I_bij_aji )
 
         enddo
@@ -237,8 +292,20 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)]
 
     enddo
   enddo
+  !$OMP END DO NOWAIT
 
-  call wall_time(tf)
+  !$OMP CRITICAL
+  do a = 1, mo_num
+    do b = 1, mo_num
+      fock_3e_uhf_mo_b(b,a) += tmp(b,a)
+    enddo
+  enddo
+  !$OMP END CRITICAL
+
+  deallocate(tmp)
+  !$OMP END PARALLEL
+
+  !call wall_time(tf)
   !print *, ' total Wall time for fock_3e_uhf_mo_b =', tf - ti
 
 END_PROVIDER 

src/tc_scf/fock_tc.irp.f

@@ -86,7 +86,7 @@ END_PROVIDER
   PROVIDE mo_l_coef mo_r_coef
   PROVIDE TCSCF_density_matrix_ao_alpha TCSCF_density_matrix_ao_beta
 
-  !print*, ' providing two_e_tc_non_hermit_integral ...'
+  !print*, ' Providing two_e_tc_non_hermit_integral ...'
   !call wall_time(t0)
 
   two_e_tc_non_hermit_integral_alpha = 0.d0
@@ -134,7 +134,7 @@ END_PROVIDER
   !$OMP END PARALLEL
 
   !call wall_time(t1)
-  !print*, ' wall time for two_e_tc_non_hermit_integral after = ', t1 - t0
+  !print*, ' Wall time for two_e_tc_non_hermit_integral = ', t1 - t0
 
 END_PROVIDER 
 
@@ -147,9 +147,16 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_alpha, (ao_num, ao_num)]
   END_DOC
 
   implicit none
+  double precision :: t0, t1
+
+  !print*, ' Providing Fock_matrix_tc_ao_alpha ...'
+  !call wall_time(t0)
 
   Fock_matrix_tc_ao_alpha = ao_one_e_integrals_tc_tot + two_e_tc_non_hermit_integral_alpha
 
+  !call wall_time(t1)
+  !print*, ' Wall time for Fock_matrix_tc_ao_alpha =', t1-t0
+
 END_PROVIDER 
 
 ! ---
@@ -175,8 +182,12 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
   END_DOC
 
   implicit none
+  double precision              :: t0, t1, tt0, tt1
   double precision, allocatable :: tmp(:,:)
 
+  !print*, ' Providing Fock_matrix_tc_mo_alpha ...'
+  !call wall_time(t0)
+
   if(bi_ortho) then
 
     !allocate(tmp(ao_num,ao_num))
@@ -188,12 +199,21 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
     !deallocate(tmp)
 
     PROVIDE mo_l_coef mo_r_coef
+
+    !call wall_time(tt0)
     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 wall_time(tt1)
+    !print*, ' 2-e term:', tt1-tt0
+
     if(three_body_h_tc) then
+      !call wall_time(tt0)
+      !PROVIDE fock_a_tot_3e_bi_orth
       !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth
       PROVIDE fock_3e_uhf_mo_a
       Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a
+      !call wall_time(tt1)
+      !print*, ' 3-e term:', tt1-tt0
     endif
 
   else
@@ -203,6 +223,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
 
   endif
 
+  !call wall_time(t1)
+  !print*, ' Wall time for Fock_matrix_tc_mo_alpha =', t1-t0
+
 END_PROVIDER
 
 ! ---
@@ -229,7 +252,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
     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
+      !PROVIDE fock_b_tot_3e_bi_orth
       !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth
+      PROVIDE fock_3e_uhf_mo_b
       Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b
     endif
 
@@ -284,6 +309,10 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
 
   implicit none
+  double precision :: t0, t1
+
+  !print*, ' Providing Fock_matrix_tc_ao_tot ...'
+  !call wall_time(t0)
 
   PROVIDE mo_l_coef mo_r_coef
   PROVIDE Fock_matrix_tc_mo_tot
@@ -291,6 +320,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
   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) )
 
+  !call wall_time(t1)
+  !print*, ' Wall time for Fock_matrix_tc_ao_tot =', t1-t0
+
 END_PROVIDER
 
 ! ---

src/tc_scf/fock_tc_mo_tot.irp.f

@@ -21,6 +21,10 @@
 
   implicit none
   integer          :: i, j, n
+  double precision :: t0, t1
+
+  !print*, ' Providing Fock_matrix_tc_mo_tot ...'
+  !call wall_time(t0)
 
   if(elec_alpha_num == elec_beta_num) then
 
@@ -154,5 +158,8 @@
     Fock_matrix_tc_mo_tot += fock_3_mat
   endif
 
+  !call wall_time(t1)
+  !print*, ' Wall time for Fock_matrix_tc_mo_tot =', t1-t0
+
 END_PROVIDER
 

src/tc_scf/fock_three_bi_ortho.irp.f

@@ -4,13 +4,21 @@
 BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)]
 
   BEGIN_DOC
+  !
   ! Alpha part of the Fock matrix from three-electron terms 
   !
   ! WARNING :: non hermitian if bi-ortho MOS used 
+  !
+  ! This calculation becomes the dominant part one the integrals are provided
+  ! 
   END_DOC
 
   implicit none
   integer          :: i, a
+  double precision :: t0, t1
+
+  !print*, ' Providing fock_a_tot_3e_bi_orth ...'
+  !call wall_time(t0)
 
   PROVIDE mo_l_coef mo_r_coef
   PROVIDE fock_cs_3e_bi_orth fock_a_tmp1_bi_ortho fock_a_tmp2_bi_ortho
@@ -25,6 +33,9 @@ BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)]
     enddo
   enddo
 
+  !call wall_time(t1)
+  !print*, ' Wall time for fock_a_tot_3e_bi_orth =', t1-t0 
+
 END_PROVIDER 
 
 ! ---
@@ -32,10 +43,15 @@ END_PROVIDER
 BEGIN_PROVIDER [double precision, fock_b_tot_3e_bi_orth, (mo_num, mo_num)]
 
   BEGIN_DOC
-! Beta part of the Fock matrix from three-electron terms 
+  !
-!
+  ! Beta part of the Fock matrix from three-electron terms 
-! WARNING :: non hermitian if bi-ortho MOS used 
+  !
+  ! WARNING :: non hermitian if bi-ortho MOS used 
+  !
+  ! This calculation becomes the dominant part one the integrals are provided
+  !
   END_DOC
+
   implicit none
   integer :: i, a
 
@@ -61,12 +77,27 @@ BEGIN_PROVIDER [double precision, fock_cs_3e_bi_orth, (mo_num, mo_num)]
   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
+  double precision              :: t0, t1
+  double precision, allocatable :: tmp(:,:)
+
+  !print*, ' Providing fock_cs_3e_bi_orth ...'
+  !call wall_time(t0)
 
   PROVIDE mo_l_coef mo_r_coef
 
+  ! to PROVIDE stuffs
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, contrib) 
+
   fock_cs_3e_bi_orth = 0.d0
 
+  !$OMP PARALLEL DEFAULT (NONE)                                                                              &
+  !$OMP PRIVATE (i, a, j, k, direct_int, c_3_int, c_minus_3_int, exch_13_int, exch_23_int, exch_12_int, tmp) &
+  !$OMP SHARED  (mo_num, elec_beta_num, fock_cs_3e_bi_orth)
+
+  allocate(tmp(mo_num,mo_num))
+  tmp = 0.d0
+
+  !$OMP DO
   do i = 1, mo_num
     do a = 1, mo_num
     
@@ -87,16 +118,29 @@ BEGIN_PROVIDER [double precision, fock_cs_3e_bi_orth, (mo_num, mo_num)]
           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
+          tmp(a,i) += 2.d0 * direct_int + 0.5d0 * (c_3_int + c_minus_3_int - exch_12_int) -1.5d0 * exch_13_int - exch_23_int
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO NOWAIT
 
-          fock_cs_3e_bi_orth(a,i) += new
+  !$OMP CRITICAL
-        enddo
+  do i = 1, mo_num
-      enddo
+    do a = 1, mo_num
+      fock_cs_3e_bi_orth(a,i) += tmp(a,i)
     enddo
   enddo
+  !$OMP END CRITICAL
+
+  deallocate(tmp)
+  !$OMP END PARALLEL
  
   fock_cs_3e_bi_orth = - fock_cs_3e_bi_orth
 
+  !call wall_time(t1)
+  !print*, ' Wall time for fock_cs_3e_bi_orth =', t1-t0
+
 END_PROVIDER 
 
 ! ---
@@ -104,20 +148,37 @@ END_PROVIDER
 BEGIN_PROVIDER [double precision, fock_a_tmp1_bi_ortho, (mo_num, mo_num)]
 
   implicit none
-  integer          :: i, a, j, k
+  integer                       :: i, a, j, k, ee
   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
+  double precision              :: t0, t1
+  double precision, allocatable :: tmp(:,:)
+
+  !print*, ' Providing fock_a_tmp1_bi_ortho ...'
+  !call wall_time(t0)
 
   PROVIDE mo_l_coef mo_r_coef
 
+  ! to PROVIDE stuffs
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, contrib) 
+
+  ee = elec_beta_num + 1
   fock_a_tmp1_bi_ortho = 0.d0
 
+  !$OMP PARALLEL DEFAULT (NONE)                                                                              &
+  !$OMP PRIVATE (i, a, j, k, direct_int, c_3_int, c_minus_3_int, exch_13_int, exch_23_int, exch_12_int, tmp) &
+  !$OMP SHARED  (mo_num, elec_alpha_num, elec_beta_num, ee, fock_a_tmp1_bi_ortho)
+
+  allocate(tmp(mo_num,mo_num))
+  tmp = 0.d0
+
+  !$OMP DO
   do i = 1, mo_num
     do a = 1, mo_num
 
-      do j = elec_beta_num + 1, elec_alpha_num 
+      do j = ee, 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 >
@@ -125,14 +186,29 @@ BEGIN_PROVIDER [double precision, fock_a_tmp1_bi_ortho, (mo_num, mo_num)]
           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)
+          tmp(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
+  !$OMP END DO NOWAIT
+
+  !$OMP CRITICAL
+  do i = 1, mo_num
+    do a = 1, mo_num
+      fock_a_tmp1_bi_ortho(a,i) += tmp(a,i)
+    enddo
+  enddo
+  !$OMP END CRITICAL
+
+  deallocate(tmp)
+  !$OMP END PARALLEL
 
   fock_a_tmp1_bi_ortho = - fock_a_tmp1_bi_ortho
 
+  !call wall_time(t1)
+  !print*, ' Wall time for fock_a_tmp1_bi_ortho =', t1-t0
+
 END_PROVIDER 
 
 ! ---
@@ -140,24 +216,56 @@ END_PROVIDER
 BEGIN_PROVIDER [double precision, fock_a_tmp2_bi_ortho, (mo_num, mo_num)]
 
   implicit none
-  integer          :: i, a, j, k
+  integer                       :: i, a, j, k, ee
   double precision              :: contrib_sss
+  double precision              :: t0, t1
+  double precision, allocatable :: tmp(:,:)
+
+  !print*, ' Providing fock_a_tmp2_bi_ortho ...'
+  !call wall_time(t0)
 
   PROVIDE mo_l_coef mo_r_coef
  
+  ! to PROVIDE stuffs
+  call contrib_3e_sss(1, 1, 1, 1, contrib_sss)
+
+  ee = elec_beta_num + 1
   fock_a_tmp2_bi_ortho = 0.d0
 
+  !$OMP PARALLEL DEFAULT (NONE)                &
+  !$OMP PRIVATE (i, a, j, k, contrib_sss, tmp) &
+  !$OMP SHARED  (mo_num, elec_alpha_num, ee, fock_a_tmp2_bi_ortho)
+
+  allocate(tmp(mo_num,mo_num))
+  tmp = 0.d0
+
+  !$OMP DO
   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
+        do k = ee, elec_alpha_num
           call contrib_3e_sss(a, i, j, k, contrib_sss)
 
-          fock_a_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_sss
+          tmp(a,i) += 0.5d0 * contrib_sss
         enddo
       enddo
     enddo
   enddo
+  !$OMP END DO NOWAIT
+
+  !$OMP CRITICAL
+  do i = 1, mo_num
+    do a = 1, mo_num
+      fock_a_tmp2_bi_ortho(a,i) += tmp(a,i)
+    enddo
+  enddo
+  !$OMP END CRITICAL
+
+  deallocate(tmp)
+  !$OMP END PARALLEL
+
+  !call wall_time(t1)
+  !print*, ' Wall time for fock_a_tmp2_bi_ortho =', t1-t0
 
 END_PROVIDER 
 
@@ -166,30 +274,61 @@ END_PROVIDER
 BEGIN_PROVIDER [double precision, fock_b_tmp1_bi_ortho, (mo_num, mo_num)]
 
   implicit none
-  integer          :: i, a, j, k
+  integer                       :: i, a, j, k, ee
   double precision              :: direct_int, exch_13_int, exch_23_int, exch_12_int
-  double precision :: new
+  double precision              :: t0, t1
+  double precision, allocatable :: tmp(:,:)
+
+  !print*, ' Providing fock_b_tmp1_bi_ortho ...'
+  !call wall_time(t0)
 
   PROVIDE mo_l_coef mo_r_coef
 
+  ! to PROVIDE stuffs
+  call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, direct_int) 
+
+  ee = elec_beta_num + 1
   fock_b_tmp1_bi_ortho = 0.d0
 
+  !$OMP PARALLEL DEFAULT (NONE)                                         &
+  !$OMP PRIVATE (i, a, j, k, direct_int, exch_13_int, exch_23_int, tmp) &
+  !$OMP SHARED  (mo_num, elec_beta_num, elec_alpha_num, ee, fock_b_tmp1_bi_ortho)
+
+  allocate(tmp(mo_num,mo_num))
+  tmp = 0.d0
+
+  !$OMP DO
   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
+        do k = ee, 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
+          tmp(a,i) += 1.5d0 * direct_int - 0.5d0 * exch_23_int - exch_13_int
         enddo
       enddo
     enddo
   enddo
+  !$OMP END DO NOWAIT
+
+  !$OMP CRITICAL
+  do i = 1, mo_num
+    do a = 1, mo_num
+      fock_b_tmp1_bi_ortho(a,i) += tmp(a,i)
+    enddo
+  enddo
+  !$OMP END CRITICAL
+
+  deallocate(tmp)
+  !$OMP END PARALLEL
 
   fock_b_tmp1_bi_ortho = - fock_b_tmp1_bi_ortho
 
+  !call wall_time(t1)
+  !print*, ' Wall time for fock_b_tmp1_bi_ortho =', t1-t0
+
 END_PROVIDER 
 
 ! ---
@@ -197,24 +336,56 @@ END_PROVIDER
 BEGIN_PROVIDER [double precision, fock_b_tmp2_bi_ortho, (mo_num, mo_num)]
 
   implicit none
-  integer          :: i, a, j, k
+  integer                       :: i, a, j, k, ee
   double precision              :: contrib_soo
+  double precision              :: t0, t1
+  double precision, allocatable :: tmp(:,:)
+
+  !print*, ' Providing fock_b_tmp2_bi_ortho ...'
+  !call wall_time(t0)
 
   PROVIDE mo_l_coef mo_r_coef
 
+  ! to PROVIDE stuffs
+  call contrib_3e_soo(1, 1, 1, 1, contrib_soo)
+
+  ee = elec_beta_num + 1
   fock_b_tmp2_bi_ortho = 0.d0
 
+  !$OMP PARALLEL DEFAULT (NONE)                &
+  !$OMP PRIVATE (i, a, j, k, contrib_soo, tmp) &
+  !$OMP SHARED  (mo_num, elec_alpha_num, ee, fock_b_tmp2_bi_ortho)
+
+  allocate(tmp(mo_num,mo_num))
+  tmp = 0.d0
+
+  !$OMP DO
   do i = 1, mo_num
     do a = 1, mo_num
-      do j = elec_beta_num + 1, elec_alpha_num 
+      do j = ee, 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
+          tmp(a,i) += 0.5d0 * contrib_soo
         enddo
       enddo
     enddo
   enddo
+  !$OMP END DO NOWAIT
+
+  !$OMP CRITICAL
+  do i = 1, mo_num
+    do a = 1, mo_num
+      fock_b_tmp2_bi_ortho(a,i) += tmp(a,i)
+    enddo
+  enddo
+  !$OMP END CRITICAL
+
+  deallocate(tmp)
+  !$OMP END PARALLEL
+
+  !call wall_time(t1)
+  !print*, ' Wall time for fock_b_tmp2_bi_ortho =', t1-t0
 
 END_PROVIDER 
 

src/tc_scf/rh_tcscf_diis.irp.f

@@ -231,8 +231,7 @@ subroutine rh_tcscf_diis()
   ! ---
 
   print *, ' TCSCF DIIS converged !'
-  call print_energy_and_mos()
+  !call print_energy_and_mos()
-
   call write_time(6)
 
   deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, E_DIIS)

src/tc_scf/routines_rotates.irp.f

@@ -46,6 +46,7 @@ end subroutine LTxR
 ! ---
 
 subroutine minimize_tc_orb_angles()
+
   BEGIN_DOC
   ! routine that minimizes the angle between left- and right-orbitals when degeneracies are found
   END_DOC

src/tc_scf/tc_scf.irp.f

@@ -54,7 +54,6 @@ program tc_scf
     endif
 
     call minimize_tc_orb_angles()
-    call print_energy_and_mos()
 
   endif
 

src/tc_scf/tc_scf_energy.irp.f

@@ -9,6 +9,10 @@
 
   implicit none
   integer          :: i, j
+  double precision :: t0, t1
+
+  !print*, ' Providing TC energy ...'
+  !call wall_time(t0)
 
   PROVIDE mo_l_coef mo_r_coef
   PROVIDE two_e_tc_non_hermit_integral_alpha two_e_tc_non_hermit_integral_beta
@@ -29,6 +33,9 @@
   TC_HF_energy += TC_HF_one_e_energy + TC_HF_two_e_energy
   TC_HF_energy += diag_three_elem_hf
 
+  !call wall_time(t1)
+  !print*, ' Wall time for TC energy=', t1-t0
+
 END_PROVIDER
 
 ! ---