LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-18 11:23:38 +01:00
Code Releases Activity

save l/r coef after diag

Browse Source
This commit is contained in:
Abdallah Ammar 2023-04-13 13:03:10 +02:00
parent 31178523fc
commit b1df3d6d03
4 changed files with 186 additions and 181 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
src/tc_bi_ortho/psi_r_l_prov.irp.f
View File

@ -192,21 +192,16 @@ subroutine save_tc_wavefunction_general(ndet, nstates, psidet, sze, dim_psicoef,
endif endif
end end
subroutine save_tc_bi_ortho_wavefunction subroutine save_tc_bi_ortho_wavefunction()
implicit none implicit none
if(save_sorted_tc_wf)then if(save_sorted_tc_wf)then
call save_tc_wavefunction_general( N_det, N_states, psi_det_sorted_tc, N_det & 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 & , size(psi_l_coef_sorted_bi_ortho, 1), psi_l_coef_sorted_bi_ortho, psi_r_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 else
call save_tc_wavefunction_general( N_det, N_states, psi_det, size(psi_det, 3) & 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 & , size(psi_l_coef_bi_ortho, 1), psi_l_coef_bi_ortho, psi_r_coef_bi_ortho )
, psi_r_coef_bi_ortho )
endif endif
call routine_save_right_bi_ortho call routine_save_right_bi_ortho()
end end
subroutine routine_save_right_bi_ortho subroutine routine_save_right_bi_ortho
@ -217,7 +212,7 @@ subroutine routine_save_right_bi_ortho
do i = 1, N_det do i = 1, N_det
coef_tmp(i,1:N_states) = psi_r_coef_sorted_bi_ortho(i,1:N_states) coef_tmp(i,1:N_states) = psi_r_coef_sorted_bi_ortho(i,1:N_states)
enddo enddo
call save_wavefunction_general_unormalized(N_det,N_states,psi_det_sorted_tc,size(coef_tmp,1),coef_tmp(1,1)) call save_wavefunction_general_unormalized(N_det, N_states, psi_det_sorted_tc, size(coef_tmp, 1), coef_tmp(1,1))
end end
subroutine routine_save_left_right_bi_ortho subroutine routine_save_left_right_bi_ortho
@ -230,6 +225,6 @@ subroutine routine_save_left_right_bi_ortho
coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1) coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1)
coef_tmp(i,2) = psi_l_coef_bi_ortho(i,1) coef_tmp(i,2) = psi_l_coef_bi_ortho(i,1)
enddo enddo
call save_wavefunction_general_unormalized(N_det,n_states_tmp,psi_det,size(coef_tmp,1),coef_tmp(1,1)) call save_wavefunction_general_unormalized(N_det, n_states_tmp, psi_det, size(coef_tmp, 1), coef_tmp(1,1))
end end

22
src/tc_bi_ortho/tc_bi_ortho.irp.f
View File

@ -39,7 +39,7 @@ end
subroutine routine_diag() subroutine routine_diag()
implicit none implicit none
integer :: i, j integer :: i, j, k
double precision :: dE double precision :: dE
! provide eigval_right_tc_bi_orth ! provide eigval_right_tc_bi_orth
@ -82,6 +82,26 @@ subroutine routine_diag()
endif 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 end

114
src/tc_bi_ortho/tc_h_eigvectors.irp.f
View File

@ -58,32 +58,36 @@ end
PROVIDE N_det N_int 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(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)) 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) H_prime(1:N_det,1:N_det) = htilde_matrix_elmt_bi_ortho(1:N_det,1:N_det)
if(s2_eig)then if(s2_eig) then
H_prime(1:N_det,1:N_det) += alpha * S2_matrix_all_dets(1:N_det,1:N_det) H_prime(1:N_det,1:N_det) += alpha * S2_matrix_all_dets(1:N_det,1:N_det)
do j=1,N_det do j=1,N_det
H_prime(j,j) = H_prime(j,j) - alpha*expected_s2 H_prime(j,j) = H_prime(j,j) - alpha*expected_s2
enddo enddo
endif endif
call non_hrmt_real_diag(N_det,H_prime,&
leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_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)
n_real_tc_bi_orth_eigval_right,eigval_right_tmp)
! do i = 1, N_det ! 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)) ! 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 ! 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) 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)) allocate(index_good_state_array(N_det),good_state_array(N_det))
i_state = 0 i_state = 0
good_state_array = .False. good_state_array = .False.
if(s2_eig)then
if (only_expected_s2) then if(s2_eig) then
do j=1,N_det
if(only_expected_s2) then
do j = 1, N_det
! Select at least n_states states with S^2 values closed to "expected_s2" ! 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) ! 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 if(dabs(s2_values_tmp(j) - expected_s2).le.0.5d0)then
i_state +=1 i_state +=1
index_good_state_array(i_state) = j index_good_state_array(i_state) = j
good_state_array(j) = .True. good_state_array(j) = .True.
@ -93,15 +97,16 @@ end
endif endif
enddo enddo
else else
do j=1,N_det do j = 1, N_det
index_good_state_array(j) = j index_good_state_array(j) = j
good_state_array(j) = .True. good_state_array(j) = .True.
enddo enddo
endif endif
if(i_state .ne.0)then
if(i_state .ne. 0) then
! Fill the first "i_state" states that have a correct S^2 value ! Fill the first "i_state" states that have a correct S^2 value
do j = 1, i_state do j = 1, i_state
do i=1,N_det do i = 1, N_det
reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,index_good_state_array(j)) 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)) leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,index_good_state_array(j))
enddo enddo
@ -116,7 +121,7 @@ end
if(i_state+i_other_state.gt.n_states)then if(i_state+i_other_state.gt.n_states)then
exit exit
endif endif
do i=1,N_det do i = 1, N_det
reigvec_tc_bi_orth(i,i_state+i_other_state) = reigvec_tc_bi_orth_tmp(i,j) 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) leigvec_tc_bi_orth(i,i_state+i_other_state) = leigvec_tc_bi_orth_tmp(i,j)
enddo enddo
@ -134,8 +139,8 @@ end
print*,' as the CI_eigenvectors' 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*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
print*,'' print*,''
do j=1,min(N_states_diag,N_det) do j = 1, min(N_states_diag, N_det)
do i=1,N_det do i = 1, N_det
leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,j) 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) reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,j)
enddo enddo
@ -146,6 +151,7 @@ end
endif ! istate .ne. 0 endif ! istate .ne. 0
else ! s2_eig else ! s2_eig
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 iorder(i) = i
@ -162,7 +168,7 @@ end
igood_l = iorder(1) 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 *,''
print *,'Warning, the left and right eigenvectors are "not the same" ' print *,'Warning, the left and right eigenvectors are "not the same" '
print *,'Warning, the ground state is not dominated by HF...' print *,'Warning, the ground state is not dominated by HF...'
@ -171,16 +177,16 @@ end
print *,'State with largest LEFT coefficient of HF ',igood_l 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) print *,'coef of HF in LEFT eigenvector = ',leigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_l)
endif endif
if(state_following_tc)then
if(state_following_tc) then
print *,'Following the states with the largest coef on HF' print *,'Following the states with the largest coef on HF'
print *,'igood_r,igood_l',igood_r,igood_l 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) eigval_right_tc_bi_orth(1) = eigval_right_tmp(i)
do j = 1, N_det do j = 1, N_det
reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i) reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i)
! print*,reigvec_tc_bi_orth(j,1)
enddo enddo
i= igood_l i = igood_l
eigval_left_tc_bi_orth(1) = eigval_right_tmp(i) eigval_left_tc_bi_orth(1) = eigval_right_tmp(i)
do j = 1, N_det do j = 1, N_det
leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i) leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i)
@ -196,20 +202,9 @@ end
enddo enddo
endif 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(:,:) double precision, allocatable :: H_jj(:),vec_tmp(:,:)
external htc_bi_ortho_calc_tdav external htc_bi_ortho_calc_tdav
@ -218,16 +213,19 @@ end
external H_tc_dagger_u_0_opt external H_tc_dagger_u_0_opt
external H_tc_s2_dagger_u_0_opt external H_tc_s2_dagger_u_0_opt
external H_tc_s2_u_0_opt external H_tc_s2_u_0_opt
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag)) allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag))
do i = 1, N_det 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)) call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
enddo enddo
!!!! Preparing the left-eigenvector
print*,'---------------------------------' print*,'---------------------------------'
print*,'---------------------------------' print*,'---------------------------------'
print*,'Computing the left-eigenvector ' print*,'Computing the left-eigenvector '
print*,'---------------------------------' print*,'---------------------------------'
print*,'---------------------------------' print*,'---------------------------------'
!!!! Preparing the left-eigenvector
vec_tmp = 0.d0 vec_tmp = 0.d0
do istate = 1, N_states 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)
@ -235,8 +233,8 @@ end
do istate = N_states+1, n_states_diag do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0 vec_tmp(istate,istate) = 1.d0
enddo 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, 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, H_tc_dagger_u_0_opt)
integer :: n_it_max,i_it integer :: n_it_max,i_it
n_it_max = 1 n_it_max = 1
converged = .False. converged = .False.
@ -244,7 +242,7 @@ end
do while (.not.converged) 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) 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 i_it += 1
if(i_it .gt. 5)exit if(i_it .gt. 5) exit
enddo enddo
do istate = 1, N_states 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)
@ -263,14 +261,14 @@ end
do istate = N_states+1, n_states_diag do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0 vec_tmp(istate,istate) = 1.d0
enddo 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, 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, H_tc_u_0_opt)
converged = .False. converged = .False.
i_it = 0 i_it = 0
do while (.not.converged) 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_dagger_u_0_opt)
i_it += 1 i_it += 1
if(i_it .gt. 5)exit if(i_it .gt. 5) exit
enddo enddo
do istate = 1, N_states 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)
@ -279,8 +277,19 @@ end
deallocate(H_jj) 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) 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) ! 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 do i = 1, N_states
norm_ground_left_right_bi_orth = 0.d0 norm_ground_left_right_bi_orth = 0.d0
do j = 1, N_det do j = 1, N_det
@ -291,27 +300,6 @@ end
print*,' <S2> = ', s2_eigvec_tc_bi_orth(i) print*,' <S2> = ', s2_eigvec_tc_bi_orth(i)
enddo 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 END_PROVIDER

2
src/tc_bi_ortho/tc_utils.irp.f
View File

@ -14,11 +14,13 @@ subroutine write_tc_energy()
!htot = htilde_matrix_elmt_bi_ortho(i,j) !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) 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 + 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
enddo enddo
O_TC = 0.d0 O_TC = 0.d0
do i = 1, N_det 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) O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k)
enddo enddo