fixed bugs + cleaning in TC-SCF DIIS

src/bi_ort_ints/semi_num_ints_mo.irp.f

@@ -122,13 +122,16 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
 
 END_PROVIDER 
 
+! ---
+
 BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)]
 
   implicit none
   integer :: ipoint
-
- print*,'providing int2_grad1_u12_bimo_transp'
   double precision :: wall0, wall1
+
+  !print *, ' providing int2_grad1_u12_bimo_transp'
+
   call wall_time(wall0)
   !$OMP PARALLEL         &
   !$OMP DEFAULT (NONE)   &
@@ -145,12 +148,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num,
    enddo
   !$OMP END DO
   !$OMP END PARALLEL
+
   call wall_time(wall1)
- print*,'Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0
+  !print *, ' Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0
 
 END_PROVIDER 
 
 ! ---
+
 BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,3, mo_num, mo_num )]
  implicit none
  integer          :: i, j, ipoint

src/bi_ortho_mos/mos_rl.irp.f

@@ -61,7 +61,7 @@ subroutine mo_to_ao_bi_ortho(A_mo, LDA_mo, A_ao, LDA_ao)
             , 0.d0, tmp_1, size(tmp_1, 1) )
 
   ! (ao_overlap x mo_r_coef) x A_mo
-  allocate( tmp_1(ao_num,mo_num) )
+  allocate( tmp_2(ao_num,mo_num) )
   call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 &
             , tmp_1, size(tmp_1, 1), A_mo, LDA_mo    &
             , 0.d0, tmp_2, size(tmp_2, 1) )
@@ -73,7 +73,7 @@ subroutine mo_to_ao_bi_ortho(A_mo, LDA_mo, A_ao, LDA_ao)
             , 0.d0, tmp_1, size(tmp_1, 1) )
 
   ! (ao_overlap x mo_r_coef) x A_mo x (ao_overlap x mo_l_coef).T
-  call dgemm( 'N', 'T', ao_num, mo_num, mo_num, 1.d0       &
+  call dgemm( 'N', 'T', ao_num, ao_num, mo_num, 1.d0       &
             , tmp_2, size(tmp_2, 1), tmp_1, size(tmp_1, 1) &
             , 0.d0, A_ao, LDA_ao )
   

src/non_hermit_dav/biorthog.irp.f

@@ -283,16 +283,16 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   ! -------------------------------------------------------------------------------------
   !
 
-  print *, ' '
+  !print *, ' '
-  print *, ' Computing the left/right eigenvectors ...'
+  !print *, ' Computing the left/right eigenvectors ...'
-  print *, ' '
+  !print *, ' '
 
   allocate(WR(n), WI(n), VL(n,n), VR(n,n)) 
   
-  print *, ' fock matrix'
+  !print *, ' fock matrix'
-  do i = 1, n
+  !do i = 1, n
-    write(*, '(1000(F16.10,X))') A(i,:)
+  !  write(*, '(1000(F16.10,X))') A(i,:)
-  enddo
+  !enddo
 
   !thr_cut = 1.d-15
   !call cancel_small_elmts(A, n, thr_cut)
@@ -301,11 +301,11 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   call lapack_diag_non_sym(n, A, WR, WI, VL, VR)
   !call lapack_diag_non_sym_new(n, A, WR, WI, VL, VR)
 
-  print *, ' '
+  !print *, ' '
-  print *, ' eigenvalues'
+  !print *, ' eigenvalues'
-  do i = 1, n
+  !do i = 1, n
-    write(*, '(1000(F16.10,X))') WR(i), WI(i)
+  !  write(*, '(1000(F16.10,X))') WR(i), WI(i)
-  enddo
+  !enddo
   !print *, ' right eigenvect bef' 
   !do i = 1, n
   !  write(*, '(1000(F16.10,X))') VR(:,i)
@@ -328,9 +328,10 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   !                  track & sort the real eigenvalues 
 
   n_good = 0
-  thr    = 1.d-3
+  !thr    = 100d0
+  thr    = Im_thresh_tcscf
   do i = 1, n
-    print*, 'Re(i) + Im(i)', WR(i), WI(i)
+    !print*, 'Re(i) + Im(i)', WR(i), WI(i)
     if(dabs(WI(i)) .lt. thr) then
       n_good += 1
     else
@@ -402,23 +403,24 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
 
   if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d) ) then
 
-    print *, ' lapack vectors are normalized and bi-orthogonalized'
+    !print *, ' lapack vectors are normalized and bi-orthogonalized'
     deallocate(S)
     return
 
-  elseif( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .gt. thr_d) ) then
+  ! accu_nd is modified after adding the normalization
+  !elseif( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .gt. thr_d) ) then
 
-    print *, ' lapack vectors are not normalized but bi-orthogonalized'
+  !  print *, ' lapack vectors are not normalized but bi-orthogonalized'
-    call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.)
+  !  call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.)
 
-    call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, thr_diag, thr_norm, .true.)
+  !  call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, thr_diag, thr_norm, .true.)
 
-    deallocate(S)
+  !  deallocate(S)
-    return
+  !  return
 
   else
 
-    print *, ' lapack vectors are not normalized neither bi-orthogonalized'
+    !print *, ' lapack vectors are not normalized neither bi-orthogonalized'
 
     ! ---
 

src/non_hermit_dav/lapack_diag_non_hermit.irp.f

@@ -930,7 +930,7 @@ subroutine check_EIGVEC(n, m, A, eigval, leigvec, reigvec, thr_diag, thr_norm, s
 
     tmp_abs = tmp_abs + tmp
     V_nrm   = V_nrm   + U_nrm 
-    write(*,'(I4,X,(100(F25.16,X)))')j,eigval(j), tmp, U_nrm
+    !write(*,'(I4,X,(100(F25.16,X)))') j,eigval(j), tmp, U_nrm
 
   enddo
 
@@ -973,7 +973,7 @@ subroutine check_EIGVEC(n, m, A, eigval, leigvec, reigvec, thr_diag, thr_norm, s
 
     tmp_abs = tmp_abs + tmp
     V_nrm   = V_nrm   + U_nrm 
-    write(*,'(I4,X,(100(F25.16,X)))')j,eigval(j), tmp, U_nrm
+    !write(*,'(I4,X,(100(F25.16,X)))') j,eigval(j), tmp, U_nrm
 
   enddo
 
@@ -1082,7 +1082,7 @@ subroutine impose_weighted_orthog_svd(n, m, W, C)
   double precision, allocatable   :: S(:,:), tmp(:,:)
   double precision, allocatable   :: U(:,:), Vt(:,:), D(:)
 
-  print *, ' apply SVD to orthogonalize & normalize weighted vectors'
+  !print *, ' apply SVD to orthogonalize & normalize weighted vectors'
 
   ! ---
 
@@ -1097,10 +1097,10 @@ subroutine impose_weighted_orthog_svd(n, m, W, C)
             , 0.d0, S, size(S, 1) )
   deallocate(tmp)
 
-  print *, ' overlap bef SVD: '
+  !print *, ' overlap bef SVD: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+  !  write(*, '(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
 
   ! ---
  
@@ -1160,10 +1160,10 @@ subroutine impose_weighted_orthog_svd(n, m, W, C)
             , 0.d0, S, size(S, 1) )
   deallocate(tmp)
 
-  print *, ' overlap aft SVD: '
+  !print *, ' overlap aft SVD: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+  !  write(*, '(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
 
   deallocate(S)
 
@@ -1185,7 +1185,7 @@ subroutine impose_orthog_svd(n, m, C)
   double precision, allocatable   :: S(:,:), tmp(:,:)
   double precision, allocatable   :: U(:,:), Vt(:,:), D(:)
 
-  print *, ' apply SVD to orthogonalize & normalize vectors'
+  !print *, ' apply SVD to orthogonalize & normalize vectors'
 
   ! ---
 
@@ -1196,10 +1196,10 @@ subroutine impose_orthog_svd(n, m, C)
             , C, size(C, 1), C, size(C, 1) &
             , 0.d0, S, size(S, 1) )
 
-  print *, ' eigenvec overlap bef SVD: '
+  !print *, ' eigenvec overlap bef SVD: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+  !  write(*, '(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
 
   ! ---
  
@@ -1224,6 +1224,7 @@ subroutine impose_orthog_svd(n, m, C)
   if(num_linear_dependencies > 0) then
     write(*,*) ' linear dependencies = ', num_linear_dependencies
     write(*,*) ' m                   = ', m
+    write(*,*) ' try with Graham-Schmidt'
     stop
   endif
 
@@ -1256,10 +1257,10 @@ subroutine impose_orthog_svd(n, m, C)
             , C, size(C, 1), C, size(C, 1) &
             , 0.d0, S, size(S, 1) )
 
-  print *, ' eigenvec overlap aft SVD: '
+  !print *, ' eigenvec overlap aft SVD: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+  !  write(*, '(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
 
   deallocate(S)
 
@@ -1296,10 +1297,10 @@ subroutine impose_orthog_svd_overlap(n, m, C, overlap)
             , 0.d0, S, size(S, 1) )
   deallocate(Stmp)
 
-  print *, ' eigenvec overlap bef SVD: '
+  !print *, ' eigenvec overlap bef SVD: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+  !  write(*, '(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
 
   ! ---
  
@@ -1358,10 +1359,10 @@ subroutine impose_orthog_svd_overlap(n, m, C, overlap)
             , 0.d0, S, size(S, 1) )
   deallocate(Stmp)
 
-  print *, ' eigenvec overlap aft SVD: '
+  !print *, ' eigenvec overlap aft SVD: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+  !  write(*, '(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
   deallocate(S)
 
 end subroutine impose_orthog_svd_overlap
@@ -1528,11 +1529,11 @@ subroutine impose_orthog_degen_eigvec(n, e0, C0)
   enddo
 
   
-  do i = 1, n
+  !do i = 1, n
-    if(deg_num(i).gt.1) then
+  !  if(deg_num(i) .gt. 1) then
-      print *, ' degen on', i, deg_num(i)
+  !    print *, ' degen on', i, deg_num(i)
-    endif
+  !  endif
-  enddo
+  !enddo
 
   ! ---
 
@@ -1677,7 +1678,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
   double precision                :: accu_d, accu_nd, s_tmp
   double precision, allocatable   :: S(:,:)
 
-  print *, ' check bi-orthonormality'
+  !print *, ' check bi-orthonormality'
 
   ! ---
 
@@ -1714,15 +1715,19 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
     enddo
   enddo
   accu_nd = dsqrt(accu_nd) / dble(m)
-  print*, '    diag acc: ', accu_d
+  !print*, '    diag acc bef = ', accu_d
-  print*, ' nondiag acc: ', accu_nd
+  !print*, ' nondiag acc bef = ', accu_nd
 
   ! ---
 
   if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then
 
     do i = 1, m
+      if(S(i,i) <= 0.d0) then
+        print *, ' overap negative'
         print *, i, S(i,i)
+        exit
+      endif
       if(dabs(S(i,i) - 1.d0) .gt. thr_d) then
         s_tmp = 1.d0 / dsqrt(S(i,i))
         do j = 1, n
@@ -1757,8 +1762,8 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
     enddo
   enddo
   accu_nd = dsqrt(accu_nd) / dble(m)
-  print *, '    diag acc: ', accu_d
+  !print *, '    diag acc aft = ', accu_d
-  print *, ' nondiag acc: ', accu_nd
+  !print *, ' nondiag acc aft = ', accu_nd
 
   deallocate(S)
 
@@ -1801,10 +1806,10 @@ subroutine check_weighted_biorthog(n, m, W, Vl, Vr, thr_d, thr_nd, accu_d, accu_
             , 0.d0, S, size(S, 1) )
   deallocate(tmp)
 
-  print *, ' overlap matrix:'
+  !print *, ' overlap matrix:'
-  do i = 1, m
+  !do i = 1, m
-    write(*,'(1000(F16.10,X))') S(i,:)
+  !  write(*,'(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
 
   accu_d  = 0.d0
   accu_nd = 0.d0
@@ -1852,17 +1857,18 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
   integer                       :: i, j
   double precision, allocatable :: SS(:,:)
 
-  print *, ' check bi-orthogonality'
+  !print *, ' check bi-orthogonality'
 
   ! ---
 
   call dgemm( 'T', 'N', m, m, n, 1.d0          &
             , Vl, size(Vl, 1), Vr, size(Vr, 1) &
             , 0.d0, S, size(S, 1) )
-  print *, ' overlap matrix:'
+
-  do i = 1, m
+  !print *, ' overlap matrix:'
-    write(*,'(1000(F16.10,X))') S(i,:)
+  !do i = 1, m
-  enddo
+  !  write(*,'(1000(F16.10,X))') S(i,:)
+  !enddo
 
   accu_d  = 0.d0
   accu_nd = 0.d0
@@ -1877,8 +1883,8 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
   enddo
   accu_nd = dsqrt(accu_nd) / dble(m)
 
-  print *, ' accu_nd = ', accu_nd
+  !print *, ' accu_nd = ', accu_nd
-  print *, ' accu_d  = ', dabs(accu_d-dble(m))/dble(m)
+  !print *, ' accu_d  = ', dabs(accu_d-dble(m))/dble(m)
 
   ! ---
 
@@ -1912,12 +1918,12 @@ subroutine check_orthog(n, m, V, accu_d, accu_nd, S)
             , V, size(V, 1), V, size(V, 1) &
             , 0.d0, S, size(S, 1) )
 
-  print *, ''
+  !print *, ''
-  print *, ' overlap matrix:'
+  !print *, ' overlap matrix:'
-  do i = 1, m
+  !do i = 1, m
-    write(*,'(1000(F16.10,X))') S(i,:)
+  !  write(*,'(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
-  print *, ''
+  !print *, ''
 
   accu_d  = 0.d0
   accu_nd = 0.d0
@@ -1981,11 +1987,11 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
     enddo
   enddo
   
-  do i = 1, n
+  !do i = 1, n
-    if(deg_num(i).gt.1) then
+  !  if(deg_num(i) .gt. 1) then
-      print *, ' degen on', i, deg_num(i), e0(i)
+  !    print *, ' degen on', i, deg_num(i), e0(i)
-    endif
+  !  endif
-  enddo
+  !enddo
 
   ! ---
 
@@ -2181,11 +2187,11 @@ subroutine impose_unique_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, C0, W0, L0,
     enddo
   enddo
   
-  do i = 1, n
+  !do i = 1, n
-    if(deg_num(i).gt.1) then
+  !  if(deg_num(i) .gt. 1) then
-      print *, ' degen on', i, deg_num(i)
+  !    print *, ' degen on', i, deg_num(i)
-    endif
+  !  endif
-  enddo
+  !enddo
 
   ! ---
 
@@ -2414,10 +2420,10 @@ subroutine impose_biorthog_svd(n, m, L, R)
             , L, size(L, 1), R, size(R, 1) &
             , 0.d0, S, size(S, 1) )
 
-  print *, ' overlap bef SVD: '
+  !print *, ' overlap bef SVD: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+  !  write(*, '(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
 
   ! ---
  
@@ -2489,10 +2495,11 @@ subroutine impose_biorthog_svd(n, m, L, R)
             , L, size(L, 1), R, size(R, 1) &
             , 0.d0, S, size(S, 1) )
 
-  print *, ' overlap aft SVD: '
+  !print *, ' overlap aft SVD: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+  !  write(*, '(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
+
   deallocate(S)
 
   ! ---
@@ -2806,10 +2813,10 @@ subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R)
             , 0.d0, S, size(S, 1) )
   deallocate(Stmp)
 
-  print *, ' overlap bef SVD: '
+  !print *, ' overlap bef SVD: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F25.16,X))') S(i,:)
+  !  write(*, '(1000(F25.16,X))') S(i,:)
-  enddo
+  !enddo
 
   ! ---
  
@@ -2886,10 +2893,11 @@ subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R)
             , 0.d0, S, size(S, 1) )
   deallocate(Stmp)
 
-  print *, ' overlap aft SVD with overlap: '
+  !print *, ' overlap aft SVD with overlap: '
-  do i = 1, m
+  !do i = 1, m
-    write(*, '(1000(F16.10,X))') S(i,:)
+  !  write(*, '(1000(F16.10,X))') S(i,:)
-  enddo
+  !enddo
+
   deallocate(S)
 
   return

src/scf_utils/roothaan_hall_scf.irp.f

@@ -29,11 +29,11 @@ END_DOC
 
   call write_time(6)
 
-  print*,'Energy of the guess = ',SCF_energy
+  print*,'energy of the guess = ',SCF_energy
   write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)')  &
     '====','================','================','================','================'
   write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)')  &
-    '  N ', 'Energy  ', 'Energy diff  ',  'DIIS error  ', 'Level shift   '
+    '  N ', 'energy  ', 'energy diff  ',  'DIIS error  ', 'Level shift   '
   write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)')  &
     '====','================','================','================','================'
 
@@ -69,9 +69,9 @@ END_DOC
     if ( (scf_algorithm == 'DIIS').and.(dabs(Delta_energy_SCF) > 1.d-6) )  then
 
       ! Store Fock and error matrices at each iteration
+      index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
       do j=1,ao_num
         do i=1,ao_num
-          index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
           Fock_matrix_DIIS (i,j,index_dim_DIIS) = Fock_matrix_AO(i,j)
           error_matrix_DIIS(i,j,index_dim_DIIS) = FPS_SPF_matrix_AO(i,j)
         enddo
@@ -106,8 +106,8 @@ END_DOC
 !   SCF energy
 
     energy_SCF = SCF_energy
-    Delta_Energy_SCF = energy_SCF - energy_SCF_previous
+    Delta_energy_SCF = energy_SCF - energy_SCF_previous
-    if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) then
+    if ( (SCF_algorithm == 'DIIS').and.(Delta_energy_SCF > 0.d0) ) then
       Fock_matrix_AO(1:ao_num,1:ao_num) = Fock_matrix_DIIS (1:ao_num,1:ao_num,index_dim_DIIS)
       Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0
       Fock_matrix_AO_beta  = Fock_matrix_AO*0.5d0
@@ -131,15 +131,17 @@ END_DOC
         call initialize_mo_coef_begin_iteration
       endif
       TOUCH mo_coef
-      Delta_Energy_SCF = SCF_energy - energy_SCF_previous
+      Delta_energy_SCF = SCF_energy - energy_SCF_previous
       energy_SCF = SCF_energy
       if (level_shift-level_shift_save > 40.d0) then
         level_shift = level_shift_save * 4.d0
         SOFT_TOUCH level_shift
         exit
       endif
+
       dim_DIIS=0
     enddo
+
     level_shift = level_shift * 0.5d0
     SOFT_TOUCH level_shift
     energy_SCF_previous = energy_SCF
@@ -175,7 +177,7 @@ END_DOC
    call save_mos
   endif
 
-  call write_double(6, Energy_SCF, 'SCF energy')
+  call write_double(6, energy_SCF, 'SCF energy')
 
   call write_time(6)
 

src/tc_keywords/EZFIO.cfg

@@ -86,13 +86,13 @@ default: False
 type: Threshold
 doc: Threshold on the convergence of the Hartree Fock energy.
 interface: ezfio,provider,ocaml
-default: 1.e-10
+default: 1.e-12
 
 [n_it_tcscf_max]
 type: Strictly_positive_int
 doc: Maximum number of SCF iterations
 interface: ezfio,provider,ocaml
-default: 500
+default: 100
 
 [selection_tc]
 type: integer
@@ -160,3 +160,9 @@ doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS]
 interface: ezfio,provider,ocaml
 default: DIIS
 
+[im_thresh_tcscf]
+type: Threshold
+doc: Thresholds on the Imag part of energy
+interface: ezfio,provider,ocaml
+default: 1.e-7
+

src/tc_scf/diago_bi_ort_tcfock.irp.f

@@ -76,6 +76,8 @@
             , 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
@@ -93,7 +95,6 @@
     enddo 
   enddo
   accu_nd = dsqrt(accu_nd) / accu_d
-
   if(accu_nd .gt. thr_nd) then
     print *, ' bi-orthog failed'
     print *, ' accu_nd MO = ', accu_nd, thr_nd
@@ -104,7 +105,10 @@
     stop
   endif
 
+  ! ---
+
   if(dabs(accu_d - dble(mo_num))/dble(mo_num) .gt. thr_d) then
+
     print *, ' mo_num     = ', mo_num 
     print *, ' accu_d  MO = ', accu_d, thr_d
     print *, ' normalizing vectors ...'
@@ -117,11 +121,42 @@
         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. thr_nd)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. thr_nd) then
+      print *, ' bi-orthog failed'
+      print *, ' accu_nd MO = ', accu_nd, thr_nd
+      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 
 

src/tc_scf/diis_tcscf.irp.f

@@ -27,6 +27,7 @@ BEGIN_PROVIDER [double precision, Q_alpha, (ao_num, ao_num) ]
 
   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) )
@@ -47,6 +48,7 @@ BEGIN_PROVIDER [ double precision, Q_beta, (ao_num, ao_num) ]
 
   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) )
@@ -113,15 +115,18 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)]
             , 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 P x S - S x P x F
+  ! 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
 
 ! ---

src/tc_scf/fock_tc.irp.f

@@ -74,68 +74,61 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_beta, (ao_num, ao_num)]
                          + two_e_tc_non_hermit_integral_beta 
 
 END_PROVIDER 
-! ---
-
-!BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
-!  implicit none
-!  BEGIN_DOC
-! ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis
-!  END_DOC
-!  Fock_matrix_tc_ao_tot = 0.5d0 * (Fock_matrix_tc_ao_alpha + Fock_matrix_tc_ao_beta)
-!END_PROVIDER
 
 ! ---
 
 BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
-  implicit none
+
   BEGIN_DOC
   ! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis
   END_DOC
+
+  implicit none
+
   if(bi_ortho) then
+
     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
     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) ]
-  implicit none
+
   BEGIN_DOC
   ! Total beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis
   END_DOC
+
+  implicit none
+
   if(bi_ortho) then
+
    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
     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, Fock_matrix_tc_mo_tot, (mo_num, mo_num)]
-!  implicit none
-!  BEGIN_DOC
-! ! Total alpha+beta  TC Fock matrix : h_c + Two-e^TC terms on the MO basis
-!  END_DOC
-!  Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta)
-!  if(three_body_h_tc) then
-!    Fock_matrix_tc_mo_tot += fock_3_mat
-!  endif
-!  !call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10)
-!END_PROVIDER 
-
 ! ---
 
  BEGIN_PROVIDER [ double precision, grad_non_hermit_left]

src/tc_scf/fock_three.irp.f

@@ -70,15 +70,27 @@ subroutine give_fock_ia_three_e_total(i,a,contrib)
 
 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
- print*,'providing diag_three_elem_hf'
+  double precision :: integral_aaa, hthree, integral_aab, integral_abb, integral_bbb
+
+  !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
@@ -103,19 +115,27 @@ BEGIN_PROVIDER [double precision, diag_three_elem_hf]
           diag_three_elem_hf += weight * contrib
        enddo
       enddo
+
       diag_three_elem_hf = - diag_three_elem_hf
+
+      ! ---
+
     else
-   double precision :: integral_aaa,hthree, integral_aab,integral_abb,integral_bbb
+
       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
+
     endif
+
   endif
+
 END_PROVIDER 
 
+! ---
 
 BEGIN_PROVIDER [ double precision, fock_3_mat_a_op_sh, (mo_num, mo_num)]
  implicit none 

src/tc_scf/rh_tcscf.irp.f

@@ -16,7 +16,8 @@ subroutine rh_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_TCSCF_save
+  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(:,:)
 
@@ -60,8 +61,8 @@ subroutine rh_tcscf()
   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(delta_energy_TCSCF) > thresh_TCSCF)             .or. &
-            (dabs(delta_gradie_TCSCF) > dsqrt(thresh_TCSCF))           )
+            (dabs(gradie_TCSCF_previous) > dsqrt(thresh_TCSCF))        )
 
     iteration_TCSCF += 1
     if(iteration_TCSCF > n_it_TCSCF_max) then
@@ -69,11 +70,6 @@ subroutine rh_tcscf()
       exit
     endif
 
-    ! TODO
-    !if(frozen_orb_scf) then
-    !  call initialize_mo_coef_begin_iteration
-    !endif
-
     ! current size of the DIIS space
     dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF)
 
@@ -97,7 +93,13 @@ subroutine rh_tcscf()
 
       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
+      !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
 
@@ -109,61 +111,54 @@ subroutine rh_tcscf()
 
     ! ---
 
-    ! TODO
-    !if(frozen_orb_scf) then
-    ! call reorder_core_orb
-    ! call initialize_mo_coef_begin_iteration
-    !endif
-
     ! calculate error vectors
     max_error_DIIS_TCSCF = maxval(abs(FQS_SQF_mo))
 
-    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
-    gradie_TCSCF       = grad_non_hermit
-    delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous
-    delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous
 
-    if((TCSCF_algorithm == 'DIIS') .and. (delta_gradie_TCSCF > 0.d0)) then
+    delta_energy_tmp = TC_HF_energy    - energy_TCSCF_previous
-      Fock_matrix_tc_ao_tot(1:ao_num,1:ao_num) = F_DIIS(1:ao_num,1:ao_num,index_dim_DIIS)
+    delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous
-      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
-    endif
 
     ! ---
 
-    level_shift_TCSCF_save = level_shift_TCSCF
+    do while((dabs(delta_energy_tmp) > 0.1d0) .and. (iteration_TCSCF > 1))
+!      print *, ' very big step  : ', delta_energy_tmp
+!      print *, ' TC level shift = ', level_shift_TCSCF
 
-    mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num)
-    mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num)
-
-    do while(delta_gradie_TCSCF > 0.d0)
-
-      mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num)
       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_r_coef mo_l_coef level_shift_TCSCF
+      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)
-
-      !if(frozen_orb_scf) then
-      !  call reorder_core_orb
-      !  call initialize_mo_coef_begin_iteration
-      !endif
       TOUCH mo_l_coef mo_r_coef
 
+      delta_energy_tmp = TC_HF_energy - energy_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
@@ -171,32 +166,17 @@ subroutine rh_tcscf()
     if(three_body_h_tc) then
       energy_TCSCF_3e  = diag_three_elem_hf
     endif
-      gradie_TCSCF       = grad_non_hermit
     delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous
     delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous
 
-      if(level_shift_TCSCF - level_shift_TCSCF_save > 40.d0) then
-        level_shift_TCSCF = level_shift_TCSCF_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
-
     energy_TCSCF_previous = energy_TCSCF
-    energy_TCSCF_1e       = TC_HF_one_e_energy
+    gradie_TCSCF_previous = gradie_TCSCF
-    energy_TCSCF_2e       = TC_HF_two_e_energy
+
-    energy_TCSCF_3e       = 0.d0
+
-    if(three_body_h_tc) then
+    level_shift_save = level_shift_TCSCF
-      energy_TCSCF_3e     = diag_three_elem_hf
+    mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num)
-    endif
+    mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num)
-    gradie_TCSCF_previous = grad_non_hermit
+
 
     print *, ' iteration         = ', iteration_TCSCF
     print *, ' total TC energy   = ', energy_TCSCF 
@@ -204,36 +184,25 @@ subroutine rh_tcscf()
     print *, ' 2-e   TC energy   = ', energy_TCSCF_2e
     print *, ' 3-e   TC energy   = ', energy_TCSCF_3e
     print *, ' |delta TC energy| = ', 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 *, ' '
 
-    if(delta_gradie_TCSCF < 0.d0) then
     call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
     call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
-      call ezfio_set_tc_scf_bitc_energy(energy_TCSCF)
-    endif
 
     if(qp_stop()) exit
   enddo
 
   ! ---
 
-  !if(iteration_TCSCF < n_it_TCSCF_max) then
-  !  mo_label = 'Canonical'
-  !endif
-
-  !if(.not.frozen_orb_scf) then
-  ! call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo, size(Fock_matrix_mo,1), size(Fock_matrix_mo, 2), mo_label, 1, .true.)
-  ! call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef, 1), 1.d-10)
-  ! call orthonormalize_mos
-  ! call save_mos
-  !endif
-  !call write_double(6, energy_TCSCF, 'TCSCF energy')
-
   call write_time(6)
 
+  deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, e_DIIS)
+
 end
 
 ! ---

src/tc_scf/routines_rotates.irp.f

@@ -116,7 +116,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
     print *, ' ------------------------------------'
     call orthog_functions(ao_num, n_degen, mo_l_coef_tmp, ao_overlap)
 
-    print *, ' Overlap lef-right '
+    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)

src/tc_scf/tc_scf.irp.f

@@ -15,8 +15,8 @@ program tc_scf
 !  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 create_guess()
-  !call orthonormalize_mos
+  call orthonormalize_mos()
 
   PROVIDE tcscf_algorithm
   if(tcscf_algorithm == 'DIIS') then
@@ -42,7 +42,8 @@ subroutine create_guess
   logical :: exists
 
   PROVIDE ezfio_filename
-  call ezfio_has_mo_basis_mo_coef(exists)
+  !call ezfio_has_mo_basis_mo_coef(exists)
+  exists = .false.
 
   if (.not.exists) then
     mo_label = 'Guess'
@@ -117,45 +118,56 @@ subroutine simple_tcscf()
 
   if(bi_ortho) then
 
-    !do while( it .lt. n_it_tcscf_max .and. (e_delta .gt. dsqrt(thresh_tcscf)) )
+    !do while(e_delta .gt. dsqrt(thresh_tcscf)) )
-    !do while( it .lt. n_it_tcscf_max .and. (e_delta .gt. thresh_tcscf) )
+    !do while(e_delta .gt. thresh_tcscf) )
-    !do while( it .lt. n_it_tcscf_max .and. (rho_delta .gt. thresh_tcscf) )
+    !do while(rho_delta .gt. thresh_tcscf) )
-    do while( it .lt. n_it_tcscf_max .and. (grad_non_hermit_right.gt. dsqrt(thresh_tcscf)) )
+    !do while(grad_non_hermit_right .gt. dsqrt(thresh_tcscf))
+    do while(grad_non_hermit .gt. dsqrt(thresh_tcscf))
 
       it += 1
-      print*,'iteration = ', it
+      if(it > n_it_tcscf_max) then
+        print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max
+        exit
+      endif
+
+
       print *, ' ***'
+      print *, ' iteration = ', it
+
       print *, ' TC HF total energy = ', TC_HF_energy
       print *, ' TC HF 1 e   energy = ', TC_HF_one_e_energy
       print *, ' TC HF 2 non hermit = ', TC_HF_two_e_energy
       if(three_body_h_tc) then
         print *, ' TC HF 3 body       = ', diag_three_elem_hf
       endif
-      print*,'***'
       e_delta = dabs(TC_HF_energy - e_save)
-      print*, 'it, delta E = ', it, e_delta
+
-      print*, 'it, gradient= ',grad_non_hermit_right
+      print *, ' delta E  = ', e_delta
+      print *, ' gradient = ', grad_non_hermit
+      !print *, ' gradient= ', grad_non_hermit_right
+
+      !rho_new   = TCSCF_bi_ort_dm_ao
+      !!print*, rho_new
+      !rho_delta = 0.d0
+      !do i = 1, ao_num 
+      !  do j = 1, ao_num 
+      !    rho_delta += dabs(rho_new(j,i) - rho_old(j,i))
+      !  enddo
+      !enddo
+      !print *, ' rho_delta =', rho_delta
+      !rho_old = rho_new
+
       e_save    = TC_HF_energy
       mo_l_coef = fock_tc_leigvec_ao
       mo_r_coef = fock_tc_reigvec_ao
-
-      rho_new   = TCSCF_bi_ort_dm_ao
-      !print*, rho_new
-      rho_delta = 0.d0
-      do i = 1, ao_num 
-        do j = 1, ao_num 
-          rho_delta += dabs(rho_new(j,i) - rho_old(j,i))
-        enddo
-      enddo
-      print*, ' rho_delta =', rho_delta
-      rho_old = rho_new
-
       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
-
       call ezfio_set_tc_scf_bitc_energy(TC_HF_energy)
 
+      print *, ' ***'
+      print *, ''
+
     enddo
 
   else
@@ -169,6 +181,7 @@ subroutine simple_tcscf()
       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