diff --git a/src/tc_bi_ortho/tc_bi_ortho.irp.f b/src/tc_bi_ortho/tc_bi_ortho.irp.f index 2d51f6f0..b0bd6be8 100644 --- a/src/tc_bi_ortho/tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/tc_bi_ortho.irp.f @@ -9,9 +9,13 @@ program tc_bi_ortho 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_diag - call save_tc_bi_ortho_wavefunction + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + 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,26 +32,53 @@ subroutine test end -subroutine routine_diag - implicit none -! provide eigval_right_tc_bi_orth -! provide overlap_bi_ortho -! provide htilde_matrix_elmt_bi_ortho - integer ::i,j - 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 - print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth - print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single - print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double - print*,'***' - print*,'e_corr_bi_orth = ',e_corr_bi_orth - print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj - print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth - print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth - print*,'Left/right eigenvectors' - 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 +subroutine routine_diag() + + implicit none + integer :: i, j + double precision :: dE + + ! 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 + print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth + print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single + print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double + print*,'***' + print*,'e_corr_bi_orth = ',e_corr_bi_orth + print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj + print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth + print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth + print*,'Left/right eigenvectors' + 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 + + diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f index f2cbb637..dc010701 100644 --- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f +++ b/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 - 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,& - leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,& - n_real_tc_bi_orth_eigval_right,eigval_right_tmp) - double precision, allocatable :: coef_hf_r(:),coef_hf_l(:) - integer, allocatable :: iorder(:) - allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(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)) + 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 & + , leigvec_tc_bi_orth_tmp, reigvec_tc_bi_orth_tmp, n_real_tc_bi_orth_eigval_right, eigval_right_tmp) + + allocate(coef_hf_r(N_det), coef_hf_l(N_det), iorder(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) - do i = 1,N_det + print*, 'igood_r, coef_hf_r = ', igood_r, coef_hf_r(1) + 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,31 +85,48 @@ 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) - do j = 1, N_det - reigvec_tc_bi_orth(j,i) = reigvec_tc_bi_orth_tmp(j,i) - leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i) - enddo - enddo + + 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) + do j = 1, N_det + reigvec_tc_bi_orth(j,i) = reigvec_tc_bi_orth_tmp(j,i) + 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 + + else + double precision, allocatable :: H_jj(:),vec_tmp(:,:) external htc_bi_ortho_calc_tdav external htcdag_bi_ortho_calc_tdav