From b1df3d6d037efe59a3e9272cbf5b7795941e9969 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 13 Apr 2023 13:03:10 +0200 Subject: [PATCH] save l/r coef after diag --- src/tc_bi_ortho/psi_r_l_prov.irp.f | 61 +++--- src/tc_bi_ortho/tc_bi_ortho.irp.f | 22 +- src/tc_bi_ortho/tc_h_eigvectors.irp.f | 282 ++++++++++++-------------- src/tc_bi_ortho/tc_utils.irp.f | 2 + 4 files changed, 186 insertions(+), 181 deletions(-) diff --git a/src/tc_bi_ortho/psi_r_l_prov.irp.f b/src/tc_bi_ortho/psi_r_l_prov.irp.f index 05d132d5..ee8abcec 100644 --- a/src/tc_bi_ortho/psi_r_l_prov.irp.f +++ b/src/tc_bi_ortho/psi_r_l_prov.irp.f @@ -192,44 +192,39 @@ 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, N_det & - , size(psi_l_coef_sorted_bi_ortho, 1), psi_l_coef_sorted_bi_ortho & - , psi_r_coef_sorted_bi_ortho ) - !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 ) - 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 +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) + 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() end subroutine routine_save_right_bi_ortho - implicit none - double precision, allocatable :: coef_tmp(:,:) - integer :: i - 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)) + implicit none + double precision, allocatable :: coef_tmp(:,:) + integer :: i + 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)) end subroutine routine_save_left_right_bi_ortho - implicit none - double precision, allocatable :: coef_tmp(:,:) - integer :: i,n_states_tmp - n_states_tmp = 2 - allocate(coef_tmp(N_det, n_states_tmp)) - do i = 1, N_det - coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1) - 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)) + implicit none + double precision, allocatable :: coef_tmp(:,:) + integer :: i,n_states_tmp + n_states_tmp = 2 + allocate(coef_tmp(N_det, n_states_tmp)) + do i = 1, N_det + coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1) + 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)) end diff --git a/src/tc_bi_ortho/tc_bi_ortho.irp.f b/src/tc_bi_ortho/tc_bi_ortho.irp.f index 9109edc4..2a1919f4 100644 --- a/src/tc_bi_ortho/tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/tc_bi_ortho.irp.f @@ -39,7 +39,7 @@ end subroutine routine_diag() implicit none - integer :: i, j + integer :: i, j, k double precision :: dE ! provide eigval_right_tc_bi_orth @@ -82,6 +82,26 @@ subroutine routine_diag() endif + double precision, allocatable :: buffer(:,:) + allocate(buffer(N_det,N_states)) + do k = 1, N_states + do i = 1, N_det + psi_l_coef_bi_ortho(i,k) = leigvec_tc_bi_orth(i,k) + buffer(i,k) = leigvec_tc_bi_orth(i,k) + enddo + enddo + TOUCH psi_l_coef_bi_ortho + call ezfio_set_tc_bi_ortho_psi_l_coef_bi_ortho(buffer) + do k = 1, N_states + do i = 1, N_det + psi_r_coef_bi_ortho(i,k) = reigvec_tc_bi_orth(i,k) + buffer(i,k) = reigvec_tc_bi_orth(i,k) + enddo + enddo + TOUCH psi_r_coef_bi_ortho + call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(buffer) + deallocate(buffer) + end diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f index 3b346efe..cc236359 100644 --- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f +++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f @@ -58,111 +58,117 @@ end 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),expect_e(N_det)) allocate (H_prime(N_det,N_det),s2_values_tmp(N_det)) + H_prime(1:N_det,1:N_det) = htilde_matrix_elmt_bi_ortho(1:N_det,1:N_det) - if(s2_eig)then - H_prime(1:N_det,1:N_det) += alpha * S2_matrix_all_dets(1:N_det,1:N_det) - do j=1,N_det - H_prime(j,j) = H_prime(j,j) - alpha*expected_s2 - enddo + if(s2_eig) then + H_prime(1:N_det,1:N_det) += alpha * S2_matrix_all_dets(1:N_det,1:N_det) + do j=1,N_det + H_prime(j,j) = H_prime(j,j) - alpha*expected_s2 + enddo endif - call non_hrmt_real_diag(N_det,H_prime,& - 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, H_prime, leigvec_tc_bi_orth_tmp, reigvec_tc_bi_orth_tmp, n_real_tc_bi_orth_eigval_right, eigval_right_tmp) ! do i = 1, N_det ! call get_H_tc_s2_l0_r0(leigvec_tc_bi_orth_tmp(1,i),reigvec_tc_bi_orth_tmp(1,i),1,N_det,expect_e(i), s2_values_tmp(i)) ! enddo call get_H_tc_s2_l0_r0(leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,N_det,N_det,expect_e, s2_values_tmp) + allocate(index_good_state_array(N_det),good_state_array(N_det)) i_state = 0 good_state_array = .False. - if(s2_eig)then - if (only_expected_s2) then - do j=1,N_det + + if(s2_eig) then + + if(only_expected_s2) then + do j = 1, N_det ! Select at least n_states states with S^2 values closed to "expected_s2" ! print*,'s2_values_tmp(j) = ',s2_values_tmp(j),eigval_right_tmp(j),expect_e(j) - if(dabs(s2_values_tmp(j)-expected_s2).le.0.5d0)then - i_state +=1 - index_good_state_array(i_state) = j - good_state_array(j) = .True. - endif - if(i_state.eq.N_states) then - exit - endif - enddo - else - do j=1,N_det - index_good_state_array(j) = j - good_state_array(j) = .True. - enddo - endif - if(i_state .ne.0)then - ! Fill the first "i_state" states that have a correct S^2 value - do j = 1, i_state - do i=1,N_det - reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,index_good_state_array(j)) - leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,index_good_state_array(j)) - enddo - eigval_right_tc_bi_orth(j) = expect_e(index_good_state_array(j)) - eigval_left_tc_bi_orth(j) = expect_e(index_good_state_array(j)) - s2_eigvec_tc_bi_orth(j) = s2_values_tmp(index_good_state_array(j)) - enddo - i_other_state = 0 - do j = 1, N_det - if(good_state_array(j))cycle - i_other_state +=1 - if(i_state+i_other_state.gt.n_states)then - exit - endif - do i=1,N_det - reigvec_tc_bi_orth(i,i_state+i_other_state) = reigvec_tc_bi_orth_tmp(i,j) - leigvec_tc_bi_orth(i,i_state+i_other_state) = leigvec_tc_bi_orth_tmp(i,j) - enddo - eigval_right_tc_bi_orth(i_state+i_other_state) = eigval_right_tmp(j) - eigval_left_tc_bi_orth (i_state+i_other_state) = eigval_right_tmp(j) - s2_eigvec_tc_bi_orth(i_state+i_other_state) = s2_values_tmp(i_state+i_other_state) - enddo - else ! istate == 0 - print*,'' - print*,'!!!!!!!! WARNING !!!!!!!!!' - print*,' Within the ',N_det,'determinants selected' - print*,' and the ',N_states_diag,'states requested' - print*,' We did not find only states with S^2 values close to ',expected_s2 - print*,' We will then set the first N_states eigenvectors of the H matrix' - print*,' as the CI_eigenvectors' - print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' - print*,'' - do j=1,min(N_states_diag,N_det) - do i=1,N_det - leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,j) - reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,j) - enddo - eigval_right_tc_bi_orth(j) = eigval_right_tmp(j) - eigval_left_tc_bi_orth (j) = eigval_right_tmp(j) - s2_eigvec_tc_bi_orth(j) = s2_values_tmp(j) - enddo - endif ! istate .ne. 0 + if(dabs(s2_values_tmp(j) - expected_s2).le.0.5d0)then + i_state +=1 + index_good_state_array(i_state) = j + good_state_array(j) = .True. + endif + if(i_state.eq.N_states) then + exit + endif + enddo + else + do j = 1, N_det + index_good_state_array(j) = j + good_state_array(j) = .True. + enddo + endif + + if(i_state .ne. 0) then + ! Fill the first "i_state" states that have a correct S^2 value + do j = 1, i_state + do i = 1, N_det + reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,index_good_state_array(j)) + leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,index_good_state_array(j)) + enddo + eigval_right_tc_bi_orth(j) = expect_e(index_good_state_array(j)) + eigval_left_tc_bi_orth(j) = expect_e(index_good_state_array(j)) + s2_eigvec_tc_bi_orth(j) = s2_values_tmp(index_good_state_array(j)) + enddo + i_other_state = 0 + do j = 1, N_det + if(good_state_array(j))cycle + i_other_state +=1 + if(i_state+i_other_state.gt.n_states)then + exit + endif + do i = 1, N_det + reigvec_tc_bi_orth(i,i_state+i_other_state) = reigvec_tc_bi_orth_tmp(i,j) + leigvec_tc_bi_orth(i,i_state+i_other_state) = leigvec_tc_bi_orth_tmp(i,j) + enddo + eigval_right_tc_bi_orth(i_state+i_other_state) = eigval_right_tmp(j) + eigval_left_tc_bi_orth (i_state+i_other_state) = eigval_right_tmp(j) + s2_eigvec_tc_bi_orth(i_state+i_other_state) = s2_values_tmp(i_state+i_other_state) + enddo + else ! istate == 0 + print*,'' + print*,'!!!!!!!! WARNING !!!!!!!!!' + print*,' Within the ',N_det,'determinants selected' + print*,' and the ',N_states_diag,'states requested' + print*,' We did not find only states with S^2 values close to ',expected_s2 + print*,' We will then set the first N_states eigenvectors of the H matrix' + print*,' as the CI_eigenvectors' + print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' + print*,'' + do j = 1, min(N_states_diag, N_det) + do i = 1, N_det + leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,j) + reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,j) + enddo + eigval_right_tc_bi_orth(j) = eigval_right_tmp(j) + eigval_left_tc_bi_orth (j) = eigval_right_tmp(j) + s2_eigvec_tc_bi_orth(j) = s2_values_tmp(j) + enddo + endif ! istate .ne. 0 else ! s2_eig - allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(N_det)) - do i = 1,N_det + + 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) - igood_r = iorder(1) - print*,'igood_r, coef_hf_r = ',igood_r,coef_hf_r(1) - do i = 1,N_det + enddo + 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 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) - igood_l = iorder(1) - print*,'igood_l, coef_hf_l = ',igood_l,coef_hf_l(1) + enddo + call dsort(coef_hf_l,iorder,N_det) + igood_l = iorder(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...' @@ -170,22 +176,22 @@ end print *,'coef of HF in RIGHT eigenvector = ',reigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_r) 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 + endif + + 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 + 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 + 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 + 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) @@ -194,22 +200,11 @@ end leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i) enddo enddo - endif - - ! 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(1,1), size(Stmp, 1) ) - print *, ' overlap matrix between states:' - do i = 1, N_states - write(*,'(1000(F16.10,X))') Stmp(i,:) - enddo - deallocate(Stmp) + endif endif - else + else ! n_det > N_det_max_full double precision, allocatable :: H_jj(:),vec_tmp(:,:) external htc_bi_ortho_calc_tdav @@ -218,36 +213,39 @@ end external H_tc_dagger_u_0_opt external H_tc_s2_dagger_u_0_opt external H_tc_s2_u_0_opt + allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag)) + do i = 1, N_det call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i)) enddo - !!!! Preparing the left-eigenvector + print*,'---------------------------------' print*,'---------------------------------' print*,'Computing the left-eigenvector ' print*,'---------------------------------' print*,'---------------------------------' + !!!! Preparing the left-eigenvector vec_tmp = 0.d0 do istate = 1, N_states - vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate) + vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate) enddo do istate = N_states+1, n_states_diag - vec_tmp(istate,istate) = 1.d0 + vec_tmp(istate,istate) = 1.d0 enddo -! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav) -! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_dagger_u_0_opt) + !call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav) + !call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_dagger_u_0_opt) integer :: n_it_max,i_it n_it_max = 1 converged = .False. i_it = 0 do while (.not.converged) - call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt) - i_it += 1 - if(i_it .gt. 5)exit + call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt) + i_it += 1 + if(i_it .gt. 5) exit enddo do istate = 1, N_states - leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) + leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) enddo print*,'---------------------------------' @@ -255,32 +253,43 @@ end print*,'Computing the right-eigenvector ' print*,'---------------------------------' print*,'---------------------------------' - !!!! Preparing the right-eigenvector + !!!! Preparing the right-eigenvector vec_tmp = 0.d0 do istate = 1, N_states - vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate) + vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate) enddo do istate = N_states+1, n_states_diag - vec_tmp(istate,istate) = 1.d0 + vec_tmp(istate,istate) = 1.d0 enddo -! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav) -! call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_u_0_opt) + !call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav) + !call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_u_0_opt) 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) - i_it += 1 - if(i_it .gt. 5)exit + 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) + i_it += 1 + if(i_it .gt. 5) exit enddo do istate = 1, N_states - reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) + reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) enddo deallocate(H_jj) - endif + endif - call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,size(reigvec_tc_bi_orth,1),N_det,N_states) - print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) + call bi_normalize(leigvec_tc_bi_orth, reigvec_tc_bi_orth, size(reigvec_tc_bi_orth, 1), N_det, N_states) + ! 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(1,1), size(Stmp, 1) ) + print *, ' overlap matrix between states:' + do i = 1, N_states + write(*,'(1000(F16.10,X))') Stmp(i,:) + enddo + deallocate(Stmp) + + print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ', leigvec_tc_bi_orth(1,1), reigvec_tc_bi_orth(1,1) do i = 1, N_states norm_ground_left_right_bi_orth = 0.d0 do j = 1, N_det @@ -291,27 +300,6 @@ end print*,' = ', s2_eigvec_tc_bi_orth(i) enddo - ! --- - - double precision, allocatable :: buffer(:,:) - allocate(buffer(N_det,N_states)) - - do k = 1, N_states - do i = 1, N_det - buffer(i,k) = leigvec_tc_bi_orth(i,k) - enddo - enddo - call ezfio_set_tc_bi_ortho_psi_l_coef_bi_ortho(buffer) - - do k = 1, N_states - do i = 1, N_det - buffer(i,k) = reigvec_tc_bi_orth(i,k) - enddo - enddo - call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(buffer) - - deallocate(buffer) - END_PROVIDER diff --git a/src/tc_bi_ortho/tc_utils.irp.f b/src/tc_bi_ortho/tc_utils.irp.f index 92e8639d..594b466c 100644 --- a/src/tc_bi_ortho/tc_utils.irp.f +++ b/src/tc_bi_ortho/tc_utils.irp.f @@ -14,11 +14,13 @@ subroutine write_tc_energy() !htot = htilde_matrix_elmt_bi_ortho(i,j) call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot) E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot + !E_TC = E_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(j,k) * htot enddo enddo O_TC = 0.d0 do i = 1, N_det + !O_TC = O_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(i,k) O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k) enddo