Merge branch 'dev-stable-tc-scf' into dev-stable-tc-scf

src/bi_ort_ints/semi_num_ints_mo.irp.f

@@ -138,10 +138,13 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
       enddo
     enddo
 
+    FREE int2_grad1_u12_ao
+
   endif
 
   call wall_time(wall1)
   print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -200,6 +203,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,
     enddo
   enddo
 
+  FREE int2_grad1_u12_bimo_transp
+
 END_PROVIDER 
 
 ! ---

src/bi_ort_ints/three_body_ijm.irp.f

@@ -49,6 +49,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_direct_bi_ort, (mo_num, mo_num,
 
   call wall_time(wall1)
   print *, ' wall time for three_e_3_idx_direct_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -102,6 +103,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_1_bi_ort, (mo_num, mo_num
 
   call wall_time(wall1)
   print *, ' wall time for three_e_3_idx_cycle_1_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -155,6 +157,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_2_bi_ort, (mo_num, mo_num
 
   call wall_time(wall1)
   print *, ' wall time for three_e_3_idx_cycle_2_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -208,6 +211,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch23_bi_ort, (mo_num, mo_num,
 
   call wall_time(wall1)
   print *, ' wall time for three_e_3_idx_exch23_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -261,6 +265,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch13_bi_ort, (mo_num, mo_num,
 
   call wall_time(wall1)
   print *, ' wall time for three_e_3_idx_exch13_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -306,6 +311,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort, (mo_num, mo_num,
 
   call wall_time(wall1)
   print *, ' wall time for three_e_3_idx_exch12_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -359,6 +365,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort_new, (mo_num, mo_
 
   call wall_time(wall1)
   print *, ' wall time for three_e_3_idx_exch12_bi_ort_new', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 

src/bi_ort_ints/three_body_ijmk.irp.f

@@ -43,6 +43,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_direct_bi_ort, (mo_num, mo_num,
 
   call wall_time(wall1)
   print *, ' wall time for three_e_4_idx_direct_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -90,6 +91,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_1_bi_ort, (mo_num, mo_num
 
   call wall_time(wall1)
   print *, ' wall time for three_e_4_idx_cycle_1_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -137,6 +139,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_2_bi_ort, (mo_num, mo_num
 
   call wall_time(wall1)
   print *, ' wall time for three_e_4_idx_cycle_2_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -184,6 +187,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch23_bi_ort, (mo_num, mo_num,
 
   call wall_time(wall1)
   print *, ' wall time for three_e_4_idx_exch23_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER
 
@@ -230,6 +234,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort, (mo_num, mo_num,
 
   call wall_time(wall1)
   print *, ' wall time for three_e_4_idx_exch13_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 
@@ -277,6 +282,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort, (mo_num, mo_num,
 
   call wall_time(wall1)
   print *, ' wall time for three_e_4_idx_exch12_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 

src/bi_ort_ints/three_body_ijmkl.irp.f

@@ -58,9 +58,11 @@
        enddo
     enddo
   enddo
+
   !$OMP END DO
   !$OMP END PARALLEL
 
+
   call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, n_points_final_grid, 1.d0, &
       orb_mat, n_points_final_grid,  &
       grad_mli, n_points_final_grid,  0.d0, &
@@ -81,6 +83,7 @@
 
   deallocate(orb_mat,grad_mli)
 
+
   allocate(lm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
   allocate(rm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
   allocate(rk_grad_im(n_points_final_grid,3,mo_num,mo_num))
@@ -132,6 +135,7 @@
   enddo
   !$OMP END PARALLEL DO
 
+
   call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
       lm_grad_ik, 3*n_points_final_grid,  &
       rk_grad_im, 3*n_points_final_grid,  0.d0, &
@@ -152,6 +156,7 @@
   enddo
   !$OMP END PARALLEL DO
 
+
   deallocate(lm_grad_ik)
 
   allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num))
@@ -198,6 +203,7 @@
   enddo
   !$OMP END PARALLEL DO
 
+
   call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
       lk_grad_mi, 3*n_points_final_grid,  &
       rk_grad_im, 3*n_points_final_grid,  0.d0, &
@@ -223,8 +229,10 @@
   deallocate(rk_grad_im)
   enddo
 
+
   call wall_time(wall1)
   print *, ' wall time for three_e_5_idx_bi_ort', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER
 

src/bi_ort_ints/three_body_ints_bi_ort.irp.f

@@ -57,6 +57,7 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n
 
   call wall_time(wall1)
   print *, ' wall time for three_body_ints_bi_ort', wall1 - wall0
+  call print_memory_usage()
 ! if(write_three_body_ints_bi_ort)then
 !  print*,'Writing three_body_ints_bi_ort on disk ...'
 !  call write_array_6_index_tensor(mo_num,three_body_ints_bi_ort,name_file)

src/non_h_ints_mu/grad_squared.irp.f

@@ -231,6 +231,7 @@ BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_g
   call wall_time(time0)
 
   PROVIDE j1b_type
+  PROVIDE int2_grad1u2_grad2u2_j1b2
 
   do ipoint = 1, n_points_final_grid
     tmp1 = v_1b(ipoint)
@@ -242,6 +243,8 @@ BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_g
     enddo
   enddo
 
+  FREE int2_grad1u2_grad2u2_j1b2
+
   !if(j1b_type .eq. 0) then 
   !  grad12_j12 = 0.d0
   !  do ipoint = 1, n_points_final_grid
@@ -262,6 +265,7 @@ BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_g
 
   call wall_time(time1)
   print*, ' Wall time for grad12_j12 = ', time1 - time0
+  call print_memory_usage()
 
 END_PROVIDER
 
@@ -278,6 +282,9 @@ BEGIN_PROVIDER [double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_g
   print*, ' providing u12sq_j1bsq ...'
   call wall_time(time0)
 
+  ! do not free here
+  PROVIDE int2_u2_j1b2
+
   do ipoint = 1, n_points_final_grid
     tmp_x = v_1b_grad(1,ipoint)
     tmp_y = v_1b_grad(2,ipoint)
@@ -292,6 +299,7 @@ BEGIN_PROVIDER [double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_g
 
   call wall_time(time1)
   print*, ' Wall time for u12sq_j1bsq = ', time1 - time0
+  call print_memory_usage()
 
 END_PROVIDER
 
@@ -310,6 +318,9 @@ BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num,
   print*, ' providing u12_grad1_u12_j1b_grad1_j1b ...'
   call wall_time(time0)
 
+  PROVIDE int2_u_grad1u_j1b2
+  PROVIDE int2_u_grad1u_x_j1b2
+
   do ipoint = 1, n_points_final_grid
 
     x     = final_grid_points(1,ipoint)
@@ -340,14 +351,17 @@ BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num,
     enddo
   enddo
 
+  FREE int2_u_grad1u_j1b2
+  FREE int2_u_grad1u_x_j1b2
+
   call wall_time(time1)
   print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b = ', time1 - time0
+  call print_memory_usage()
 
 END_PROVIDER
 
 ! ---
 
-
 BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao_num)]
 
   BEGIN_DOC
@@ -401,6 +415,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
               , int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
               , 0.d0, tc_grad_square_ao, ao_num*ao_num)
 
+    FREE int2_grad1_u12_square_ao
+
     ! ---
 
     if(((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) .and. use_ipp) then
@@ -442,6 +458,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
       call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0     &
                 , int2_u2_j1b2(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
                 , 1.d0, tc_grad_square_ao, ao_num*ao_num)
+
+      FREE int2_u2_j1b2
     endif
 
     ! ---
@@ -478,6 +496,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 
   call wall_time(time1)
   print*, ' Wall time for tc_grad_square_ao = ', time1 - time0
+  call print_memory_usage()
 
 END_PROVIDER
 

src/non_h_ints_mu/new_grad_tc.irp.f

@@ -284,6 +284,7 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
 
   call wall_time(time1)
   print*, ' Wall time for tc_grad_and_lapl_ao = ', time1 - time0
+  call print_memory_usage()
 
 END_PROVIDER 
 

src/non_h_ints_mu/tc_integ.irp.f

@@ -100,6 +100,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
       !$OMP END DO
       !$OMP END PARALLEL
 
+      FREE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b
+
     elseif(j1b_type .ge. 100) then
 
       PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
@@ -176,6 +178,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
 
   call wall_time(time1)
   print*, ' wall time for int2_grad1_u12_ao =', time1-time0 
+  call print_memory_usage()
 
 END_PROVIDER
 
@@ -242,6 +245,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
       !$OMP END DO
       !$OMP END PARALLEL
 
+      FREE u12sq_j1bsq grad12_j12
+
     else
 
       PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
@@ -262,6 +267,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
       !$OMP END DO
       !$OMP END PARALLEL
 
+      FREE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
+
     endif
 
   elseif(j1b_type .ge. 100) then
@@ -324,6 +331,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
 
   call wall_time(time1)
   print*, ' wall time for int2_grad1_u12_square_ao =', time1-time0 
+  call print_memory_usage()
 
 END_PROVIDER
 

src/non_h_ints_mu/total_tc_int.irp.f

@@ -84,8 +84,11 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
     enddo
   endif
 
+  FREE tc_grad_square_ao tc_grad_and_lapl_ao ao_two_e_coul
+
   call wall_time(wall1)
   print *, ' wall time for ao_tc_int_chemist ', wall1 - wall0
+  call print_memory_usage()
 
 END_PROVIDER 
 

src/tc_bi_ortho/normal_ordered.irp.f

@@ -1,3 +1,6 @@
+
+! ---
+
 BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_num, mo_num)]
 
   BEGIN_DOC 
@@ -11,10 +14,26 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
   integer                        :: i, h1, p1, h2, p2
   integer                        :: hh1, hh2, pp1, pp2
   integer                        :: Ne(2)
-  integer,           allocatable :: occ(:,:)
-  integer(bit_kind), allocatable :: key_i_core(:,:)
   double precision               :: hthree_aba, hthree_aaa, hthree_aab
   double precision               :: wall0, wall1
+  integer,           allocatable :: occ(:,:)
+  integer(bit_kind), allocatable :: key_i_core(:,:)
+
+  print*,' Providing normal_two_body_bi_orth ...'
+  call wall_time(wall0)
+ 
+  PROVIDE N_int
+
+  print*,' Providing normal_two_body_bi_orth ...'
+  call wall_time(wall0)
+
+  if(read_tc_norm_ord) then
+
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="read")
+      read(11) normal_two_body_bi_orth
+    close(11)
+
+  else
 
     PROVIDE N_int
 
@@ -32,8 +51,6 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
     endif
 
     normal_two_body_bi_orth = 0.d0
-  print*,'Providing normal_two_body_bi_orth ...'
-  call wall_time(wall0)
 
     !$OMP PARALLEL                                                                         &
     !$OMP DEFAULT (NONE)                                                                   &
@@ -83,15 +100,27 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
     !$OMP END DO
     !$OMP END PARALLEL
 
+    deallocate( occ )
+    deallocate( key_i_core )
+  endif
+
+  if(write_tc_norm_ord.and.mpi_master) then
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="write")
+      call ezfio_set_work_empty(.False.)
+      write(11) normal_two_body_bi_orth
+      close(11)
+      call ezfio_set_tc_keywords_io_tc_integ('Read')
+  endif
+
   call wall_time(wall1)
   print*,' Wall time for normal_two_body_bi_orth ', wall1-wall0
 
-  deallocate( occ )
+  call wall_time(wall1)
-  deallocate( key_i_core )
+  print*,' Wall time for normal_two_body_bi_orth ', wall1-wall0
 
 END_PROVIDER 
 
-
+! ---
 
 subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
 
@@ -106,30 +135,41 @@ subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
 
   !!!! double alpha/beta
   hthree = 0.d0
+
   do ii = 1, Ne(2) ! purely closed shell part 
     i = occ(ii,2)
+
     call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
     int_direct = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
     int_exc_13 = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
     int_exc_12 = -1.d0 * integral
+
     hthree += 2.d0 * int_direct - 1.d0 * (int_exc_13 + int_exc_12)
   enddo
+
   do ii = Ne(2) + 1, Ne(1) ! purely open-shell part 
     i = occ(ii,1)
+
     call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
     int_direct = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
     int_exc_13 = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
     int_exc_12 = -1.d0 * integral
+
     hthree += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12)
   enddo
 
-end subroutine give_aba_contraction
+  return
-
+end
 
+! ---
 
 BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num, mo_num, mo_num)]
 
@@ -161,6 +201,7 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num,
   else
     call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
   endif
+
   normal_two_body_bi_orth_ab = 0.d0
   do hh1 = 1, n_act_orb
     h1 = list_act(hh1) 
@@ -171,6 +212,7 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num,
         do pp2 = 1, n_act_orb
           p2 = list_act(pp2)
           call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree)
+
           normal_two_body_bi_orth_ab(p2,h2,p1,h1) = hthree    
         enddo
       enddo
@@ -182,7 +224,7 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num,
 
 END_PROVIDER 
 
-
+! ---
 
 BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_aa_bb, (n_act_orb, n_act_orb, n_act_orb, n_act_orb)]
 
@@ -250,13 +292,14 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_aa_bb, (n_act_orb, n_
 
 END_PROVIDER 
 
-
+! ---
 
 subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
 
   BEGIN_DOC
   ! pure same spin contribution to same spin double excitation s1=h1,p1, s2=h2,p2, with s1==s2
   END_DOC
+
   use bitmasks ! you need to include the bitmasks_module.f90 features
 
   implicit none
@@ -270,46 +313,62 @@ subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
   hthree = 0.d0
   do ii = 1, Ne(2) ! purely closed shell part 
     i = occ(ii,2)
+
     call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
     int_direct = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p2, p1, i, i, h2, h1, integral)
     int_exc_l = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral)
     int_exc_ll= -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
     int_exc_12= -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
     int_exc_13= -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral)
     int_exc_23= -1.d0 * integral
 
     hthree +=  1.d0 * int_direct + int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23)
   enddo
+
   do ii = Ne(2)+1,Ne(1) ! purely open-shell part 
     i = occ(ii,1)
+
     call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
     int_direct = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p2, p1, i , i, h2, h1, integral)
     int_exc_l = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral)
     int_exc_ll = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
     int_exc_12 = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
     int_exc_13 = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral)
     int_exc_23 = -1.d0 * integral
 
     hthree +=  1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23))
   enddo
 
-end subroutine give_aaa_contraction
+  return
-
+end
 
+! ---
 
 subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
-  implicit none
+
   use bitmasks ! you need to include the bitmasks_module.f90 features
+
+  implicit none
   integer,          intent(in)  :: Nint, h1, h2, p1, p2
   integer,          intent(in)  :: Ne(2), occ(Nint*bit_kind_size,2)
   double precision, intent(out) :: hthree
@@ -320,11 +379,18 @@ subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
   hthree = 0.d0
   do ii = 1, Ne(2) ! purely closed shell part 
     i = occ(ii,2)
+
     call give_integrals_3_body_bi_ort(p2, p1, i, h2, h1, i, integral)
     int_direct = -1.d0 * integral
+
     call give_integrals_3_body_bi_ort(p1, p2, i, h2, h1, i, integral)
     int_exc_23= -1.d0 * integral
+
     hthree +=  1.d0 * int_direct - int_exc_23
   enddo
 
-end subroutine give_aab_contraction
+  return
+end
+
+! ---
+

src/tc_keywords/EZFIO.cfg

@@ -226,6 +226,12 @@ doc: Read/Write integrals int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_l
 interface: ezfio,provider,ocaml
 default: None
 
+[io_tc_norm_ord]
+type: Disk_access
+doc: Read/Write normal_two_body_bi_orth from/to disk [ Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
+
 [debug_tc_pt2]
 type: integer
 doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body

src/tc_scf/rh_tcscf_diis.irp.f

@@ -11,6 +11,7 @@ subroutine rh_tcscf_diis()
 
   integer                       :: i, j, it
   integer                       :: dim_DIIS, index_dim_DIIS
+  logical                       :: converged
   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
@@ -92,8 +93,9 @@ subroutine rh_tcscf_diis()
 
   PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot
 
+  converged = .false.
   !do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. dsqrt(thresh_tcscf)))
-  do while(er_DIIS .gt. dsqrt(thresh_tcscf))
+  do while(.not. converged)
 
     call wall_time(t0)
 
@@ -218,21 +220,56 @@ subroutine rh_tcscf_diis()
     !g_delta_th  = dabs(tc_grad) ! g_delta)
     er_delta_th = dabs(er_DIIS) !er_delta)
 
+    converged = er_DIIS .lt. dsqrt(thresh_tcscf)
+
     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
     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, I4,1X, F8.2)')  &
       it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
 
+
+!   Write data in JSON file
+
+    call lock_io
+    if (it == 1) then
+      write(json_unit, json_dict_uopen_fmt)
+    else
+      write(json_unit, json_dict_close_uopen_fmt)
+    endif
+    write(json_unit, json_int_fmt)  ' iteration      ', it
+    write(json_unit, json_real_fmt) ' SCF TC Energy  ', etc_tot
+    write(json_unit, json_real_fmt) ' E(1e)          ', etc_1e
+    write(json_unit, json_real_fmt) ' E(2e)          ', etc_2e
+    write(json_unit, json_real_fmt) ' E(3e)          ', etc_3e
+    write(json_unit, json_real_fmt) ' delta Energy   ', e_delta
+    write(json_unit, json_real_fmt) ' DIIS error     ', er_DIIS
+    write(json_unit, json_real_fmt) ' level_shift    ', level_shift_tcscf
+    write(json_unit, json_real_fmt) ' DIIS           ', dim_DIIS
+    write(json_unit, json_real_fmt) ' Wall time (min)', (t1-t0)/60.d0
+    call unlock_io
+
     if(er_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)
+      write(json_unit, json_true_fmt) 'saved'
+    else
+      write(json_unit, json_false_fmt) 'saved'
     endif
+    call lock_io
 
+    if (converged) then
+      write(json_unit, json_true_fmtx) 'converged'
+    else
+      write(json_unit, json_false_fmtx) 'converged'
+    endif
+    call unlock_io
     if(qp_stop()) exit
   enddo
 
+  write(json_unit, json_dict_close_fmtx)
+
   ! ---
 
   print *, ' TCSCF DIIS converged !'

src/tc_scf/tc_scf.irp.f

@@ -8,6 +8,8 @@ program tc_scf
 
   implicit none
 
+  write(json_unit,json_array_open_fmt) 'tc-scf'
+
   print *, ' starting ...'
 
   my_grid_becke  = .True.
@@ -57,6 +59,8 @@ program tc_scf
 
   endif
 
+  write(json_unit,json_array_close_fmtx)
+  call json_close
 
 end