trying to work on natorb

bin/qp_reset

@@ -97,6 +97,8 @@ if [[ $dets -eq 1 ]] ; then
   rm --force -- ${ezfio}/determinants/psi_{det,coef}.gz
   rm --force -- ${ezfio}/determinants/n_det_qp_edit
   rm --force -- ${ezfio}/determinants/psi_{det,coef}_qp_edit.gz
+  rm --force -- ${ezfio}/tc_bi_ortho/psi_{l,r}_coef_bi_ortho.gz 
+
 fi
 
 if [[ $mos -eq 1 ]] ; then

external/ezfio

@@ -1 +1 @@
-Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93
+Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3

src/non_hermit_dav/biorthog.irp.f

@@ -331,7 +331,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   !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
@@ -405,7 +405,7 @@ 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
 
@@ -422,7 +422,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
 
   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

@@ -1857,7 +1857,7 @@ 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'
 
   ! ---
 
@@ -1865,10 +1865,10 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
             , Vl, size(Vl, 1), Vr, size(Vr, 1) &
             , 0.d0, S, size(S, 1) )
 
-  !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
@@ -1883,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)
 
   ! ---
 
@@ -1987,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
 
   ! ---
 
@@ -2010,7 +2010,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
 
       ! ---
 
-      call impose_orthog_svd(n, m, L)
+!      call impose_orthog_svd(n, m, L)
       call impose_orthog_svd(n, m, R)
       !call impose_orthog_GramSchmidt(n, m, L)
       !call impose_orthog_GramSchmidt(n, m, R)
@@ -2030,7 +2030,8 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
       !call bi_ortho_s_inv_half(m, L, R, S_inv_half)
       !deallocate(S, S_inv_half)
 
-      call impose_biorthog_svd(n, m, L, R)
+!      call impose_biorthog_svd(n, m, L, R)
+!     call impose_biorthog_inverse(n, m, L, R)
 
       !call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
 
@@ -2045,6 +2046,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
 
     endif
   enddo
+  call impose_biorthog_inverse(n, n, L0, R0)
 
 end subroutine impose_biorthog_degen_eigvec 
 
@@ -2420,10 +2422,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
 
   ! ---
  
@@ -2495,10 +2497,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 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)
 
@@ -2506,6 +2508,50 @@ subroutine impose_biorthog_svd(n, m, L, R)
 
 end subroutine impose_biorthog_svd
 
+subroutine impose_biorthog_inverse(n, m, L, R)
+
+  implicit none
+
+  integer,          intent(in)    :: n, m
+  double precision, intent(inout) :: L(n,m)
+  double precision, intent(in)    :: R(n,m)
+  double precision, allocatable :: Lt(:,:),S(:,:)
+  integer :: i,j
+  allocate(Lt(m,n))
+  allocate(S(m,m))
+
+  call dgemm( 'T', 'N', m, m, n, 1.d0      &
+            , L, size(L, 1), R, size(R, 1) &
+            , 0.d0, S, size(S, 1) )
+
+  print *, ' overlap bef SVD: '
+  do i = 1, m
+    write(*, '(1000(F16.10,X))') S(i,:)
+  enddo
+
+  call get_pseudo_inverse(R,n,n,m,Lt,m,1.d-6)
+  do i = 1, m
+   do j = 1, n
+    L(j,i) = Lt(i,j)
+   enddo
+  enddo
+  ! ---
+
+  call dgemm( 'T', 'N', m, m, n, 1.d0      &
+            , L, size(L, 1), R, size(R, 1) &
+            , 0.d0, S, size(S, 1) )
+
+  print *, ' overlap aft SVD: '
+  do i = 1, m
+    write(*, '(1000(F16.10,X))') S(i,:)
+  enddo
+
+  deallocate(S,Lt)
+
+
+end subroutine impose_biorthog_svd
+
+
 ! ---
 
 subroutine impose_weighted_biorthog_qr(m, n, thr_d, thr_nd, Vl, W, Vr)

src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f

@@ -22,6 +22,7 @@ program tc_natorb_bi_ortho
 
   call print_energy_and_mos()
   call save_tc_natorb()
+  call print_angles_tc()
   !call minimize_tc_orb_angles()
 
 end
@@ -35,9 +36,12 @@ subroutine save_tc_natorb()
   print*,'Saving the natorbs '
 
   provide natorb_tc_leigvec_ao natorb_tc_reigvec_ao
+  mo_l_coef = natorb_tc_leigvec_ao 
+  mo_r_coef = natorb_tc_reigvec_ao 
+  touch mo_l_coef mo_r_coef 
 
-  call ezfio_set_bi_ortho_mos_mo_l_coef(natorb_tc_leigvec_ao)
+  call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
-  call ezfio_set_bi_ortho_mos_mo_r_coef(natorb_tc_reigvec_ao)
+  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
   call save_ref_determinant_nstates_1()
   call ezfio_set_determinants_read_wf(.False.)
 

src/tc_scf/routines_rotates.irp.f

@@ -402,6 +402,7 @@ subroutine print_energy_and_mos(good_angles)
   print *, ' TC energy              = ', TC_HF_energy
   print *, ' TC SCF energy gradient = ', grad_non_hermit
   print *, ' Max angle Left/right   = ', max_angle_left_right
+  call print_angles_tc()
 
   if(max_angle_left_right .lt. thresh_lr_angle) then
     print *, ' Maximum angle BELOW 45 degrees, everthing is OK !'