Merge branch 'dev-stable-tc-scf' of https://github.com/AbdAmmar/qp2 into dev-stable-tc-scf

src/tc_bi_ortho/tc_bi_ortho.irp.f

@@ -10,8 +10,12 @@ program tc_bi_ortho
   read_wf = .True.
   touch read_wf
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-  call routine_diag
+
-  call save_tc_bi_ortho_wavefunction
+  print*, ' nb of states = ', N_states
+  print*, ' nb of det    = ', N_det
+
+  call routine_diag()
+  call save_tc_bi_ortho_wavefunction()
 end
 
 subroutine test
@@ -28,12 +32,18 @@ subroutine test
 
 end
 
-subroutine routine_diag
+subroutine routine_diag()
+
   implicit none
-! provide eigval_right_tc_bi_orth
+  integer          :: i, j
-!  provide overlap_bi_ortho
+  double precision :: dE
-!  provide htilde_matrix_elmt_bi_ortho
+
- integer ::i,j
+  ! provide eigval_right_tc_bi_orth
+  ! provide overlap_bi_ortho
+  ! provide htilde_matrix_elmt_bi_ortho
+
+  if(N_states .eq. 1) then
+
     print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1)
     print*,'e_tc_left_right         = ',e_tc_left_right
     print*,'e_tilde_bi_orth_00      = ',e_tilde_bi_orth_00
@@ -49,5 +59,26 @@ subroutine routine_diag
     do i = 1,N_det
       write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1)
     enddo
+
+  else
+
+    print*,'eigval_right_tc_bi_orth : '
+    do i = 1, N_states
+      print*, i, eigval_right_tc_bi_orth(i)
+    enddo
+
+    print*,''
+    print*,'******************************************************'
+    print*,'Excitation energies (au)                     (eV)'
+    do i = 2, N_states
+      dE = eigval_right_tc_bi_orth(i) - eigval_right_tc_bi_orth(1)
+      print*, i, dE, dE/0.0367502d0
+    enddo
+    print*,''
+
+  endif
+
 end
 
+
+

src/tc_bi_ortho/tc_h_eigvectors.irp.f

@@ -47,34 +47,36 @@ end
   integer                       :: i, idx_dress, j, istate
   logical                       :: converged, dagger
   integer                       :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l
-  double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:),leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
+  integer,          allocatable :: iorder(:)
+  double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:), leigvec_tc_bi_orth_tmp(:,:), eigval_right_tmp(:)
+  double precision, allocatable :: coef_hf_r(:),coef_hf_l(:), Stmp(:,:)
 
   PROVIDE N_det N_int
 
-   if(n_det.le.N_det_max_full)then
+  if(n_det .le. N_det_max_full) then
-    allocate(reigvec_tc_bi_orth_tmp(N_det,N_det),leigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det))
+
-    call non_hrmt_real_diag(N_det,htilde_matrix_elmt_bi_ortho,& 
+    allocate(reigvec_tc_bi_orth_tmp(N_det,N_det), leigvec_tc_bi_orth_tmp(N_det,N_det), eigval_right_tmp(N_det))
-         leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,& 
+
-         n_real_tc_bi_orth_eigval_right,eigval_right_tmp)
+    call non_hrmt_real_diag( N_det, htilde_matrix_elmt_bi_ortho &
-    double precision, allocatable :: coef_hf_r(:),coef_hf_l(:)
+                           , leigvec_tc_bi_orth_tmp, reigvec_tc_bi_orth_tmp, n_real_tc_bi_orth_eigval_right, eigval_right_tmp)
-    integer, allocatable :: iorder(:)
+
-    allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(N_det))
+    allocate(coef_hf_r(N_det), coef_hf_l(N_det), iorder(N_det))
-    do i = 1,N_det
+    do i = 1, N_det
       iorder(i) = i
       coef_hf_r(i) = -dabs(reigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
     enddo
-    call dsort(coef_hf_r,iorder,N_det)
+    call dsort(coef_hf_r, iorder, N_det)
     igood_r = iorder(1)
-    print*,'igood_r, coef_hf_r = ',igood_r,coef_hf_r(1)
+    print*, 'igood_r, coef_hf_r = ', igood_r, coef_hf_r(1)
-    do i = 1,N_det
+    do i = 1, N_det
      iorder(i) = i
      coef_hf_l(i) = -dabs(leigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
     enddo
-    call dsort(coef_hf_l,iorder,N_det)
+    call dsort(coef_hf_l, iorder, N_det)
     igood_l = iorder(1)
-    print*,'igood_l, coef_hf_l = ',igood_l,coef_hf_l(1)
+    print*, 'igood_l, coef_hf_l = ', igood_l, coef_hf_l(1)
 
-    if(igood_r.ne.igood_l.and.igood_r.ne.1)then
+    if(igood_r .ne. igood_l .and. igood_r .ne. 1)then
      print *,''
      print *,'Warning, the left and right eigenvectors are "not the same" '
      print *,'Warning, the ground state is not dominated by HF...'
@@ -83,21 +85,24 @@ end
      print *,'State with largest LEFT  coefficient of HF ',igood_l
      print *,'coef of HF in LEFT  eigenvector = ',leigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_l)
     endif
-    if(state_following_tc)then
+
+    if(state_following_tc) then
+
      print *,'Following the states with the largest coef on HF'
      print *,'igood_r,igood_l',igood_r,igood_l
      i= igood_r
      eigval_right_tc_bi_orth(1) = eigval_right_tmp(i)
      do j = 1, N_det
        reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i)
-!       print*,reigvec_tc_bi_orth(j,1)
      enddo
      i= igood_l
      eigval_left_tc_bi_orth(1)  = eigval_right_tmp(i)
      do j = 1, N_det
        leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i)
      enddo
+
     else 
+
       do i = 1, N_states
         eigval_right_tc_bi_orth(i) = eigval_right_tmp(i)
         eigval_left_tc_bi_orth(i)  = eigval_right_tmp(i)
@@ -106,8 +111,22 @@ end
           leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i)
         enddo
       enddo
+
+      ! check bi-orthogonality
+      allocate(Stmp(N_states,N_states))
+      call dgemm( 'T', 'N', N_states, N_states, N_det, 1.d0                                                              &
+            , leigvec_tc_bi_orth(1,1), size(leigvec_tc_bi_orth, 1), reigvec_tc_bi_orth(1,1), size(reigvec_tc_bi_orth, 1) &
+            , 0.d0, Stmp, size(Stmp, 1) )
+      print *, ' overlap matrix between states:'
+      do i = 1, N_states
+        write(*,'(1000(F16.10,X))') Stmp(i,:)
+      enddo
+      deallocate(Stmp)
+
     endif
+
   else 
+
     double precision, allocatable :: H_jj(:),vec_tmp(:,:)
     external                         htc_bi_ortho_calc_tdav
     external                         htcdag_bi_ortho_calc_tdav