fixed bug in tc_bi_ortho

src/tc_bi_ortho/psi_r_l_prov.irp.f

@@ -192,33 +192,51 @@ subroutine save_tc_wavefunction_general(ndet, nstates, psidet, sze, dim_psicoef,
   endif
 end
 
+! ---
+
 subroutine save_tc_bi_ortho_wavefunction()
+
   implicit none
+
   if(save_sorted_tc_wf) then
+
     call save_tc_wavefunction_general( N_det, N_states, psi_det_sorted_tc, size(psi_det_sorted_tc, 3) &
                                      , size(psi_l_coef_sorted_bi_ortho, 1), psi_l_coef_sorted_bi_ortho, psi_r_coef_sorted_bi_ortho)
+    call routine_save_right_sorted_bi_ortho()
+
   else
+
     call save_tc_wavefunction_general( N_det, N_states, psi_det, size(psi_det, 3) &
                                      , size(psi_l_coef_bi_ortho, 1), psi_l_coef_bi_ortho, psi_r_coef_bi_ortho )
-  endif
     call routine_save_right_bi_ortho()
+
+  endif
+
 end
 
-subroutine routine_save_right_bi_ortho
+! ---
+
+subroutine routine_save_right_sorted_bi_ortho()
+
   implicit none
-  double precision, allocatable :: coef_tmp(:,:)
   integer                       :: i
+  double precision, allocatable :: coef_tmp(:,:)
+
   allocate(coef_tmp(N_det, N_states))
   do i = 1, N_det
     coef_tmp(i,1:N_states) = psi_r_coef_sorted_bi_ortho(i,1:N_states)
   enddo
   call save_wavefunction_general_unormalized(N_det, N_states, psi_det_sorted_tc, size(coef_tmp, 1), coef_tmp(1,1))
+  deallocate(coef_tmp)
+
 end
 
-subroutine routine_save_left_right_bi_ortho
+subroutine routine_save_left_right_sorted_bi_ortho()
+
   implicit none
-  double precision, allocatable :: coef_tmp(:,:)
   integer                       :: i, n_states_tmp
+  double precision, allocatable :: coef_tmp(:,:)
+
   n_states_tmp = 2
   allocate(coef_tmp(N_det, n_states_tmp))
   do i = 1, N_det
@@ -226,5 +244,26 @@ subroutine routine_save_left_right_bi_ortho
     coef_tmp(i,2) = psi_l_coef_bi_ortho(i,1)
   enddo
   call save_wavefunction_general_unormalized(N_det, n_states_tmp, psi_det, size(coef_tmp, 1), coef_tmp(1,1))
+  deallocate(coef_tmp)
 end
 
+! ---
+
+subroutine routine_save_right_bi_ortho()
+
+  implicit none
+  integer                       :: i
+  double precision, allocatable :: coef_tmp(:,:)
+
+  allocate(coef_tmp(N_det, N_states))
+  do i = 1, N_det
+    coef_tmp(i,1:N_states) = psi_r_coef_bi_ortho(i,1:N_states)
+  enddo
+  call save_wavefunction_general_unormalized(N_det, N_states, psi_det, size(coef_tmp, 1), coef_tmp(1,1))
+  deallocate(coef_tmp)
+
+end
+
+! ---
+
+

src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f

@@ -1,15 +0,0 @@
-program tc_bi_ortho
-  implicit none
-  BEGIN_DOC
-! TODO : Put the documentation of the program here
-  END_DOC
-  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 routine_save_left_right_bi_ortho
-! call test
-end

src/tc_bi_ortho/tc_h_eigvectors.irp.f

@@ -266,7 +266,7 @@ end
     converged = .False.
     i_it = 0
     do while (.not. converged)
-      call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
+      call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt)
       i_it += 1
       if(i_it .gt. 5) exit
     enddo
@@ -325,12 +325,19 @@ END_PROVIDER
 
 
 subroutine bi_normalize(u_l, u_r, n, ld, nstates)
+
+  BEGIN_DOC
   !!!! Normalization of the scalar product of the left/right eigenvectors
-  double precision, intent(inout) :: u_l(ld,nstates), u_r(ld,nstates)
+  END_DOC
+
+  implicit none
   integer,          intent(in)    :: n, ld, nstates
-  integer :: i
+  double precision, intent(inout) :: u_l(ld,nstates), u_r(ld,nstates)
+  integer                         :: i, j
   double precision                :: accu, tmp
+
   do i = 1, nstates
+
     !!!! Normalization of right eigenvectors |Phi>
     accu = 0.d0
     do j = 1, n
@@ -345,12 +352,14 @@ subroutine bi_normalize(u_l,u_r,n,ld,nstates)
     do j = 1, n
       u_r(j,i) *= tmp
     enddo
+
     !!!! Adaptation of the norm of the left eigenvector such that <chi|Phi> = 1
     accu = 0.d0
     do j = 1, n
       accu += u_l(j,i) * u_r(j,i)
       !print*,j, u_l(j,i) , u_r(j,i)
     enddo
+    print*,'accu_lr = ', accu
     if(accu.gt.0.d0)then
       accu = 1.d0/dsqrt(accu)
     else
@@ -361,5 +370,8 @@ subroutine bi_normalize(u_l,u_r,n,ld,nstates)
       u_l(j,i) *= accu * tmp
       u_r(j,i) *= accu
     enddo
+
   enddo
+
 end
+