9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-13 07:03:29 +01:00
qp2/plugins/local/tc_bi_ortho/tc_h_eigvectors.irp.f

403 lines
14 KiB
Fortran
Raw Normal View History

2023-07-02 21:49:25 +02:00
! ---
2023-02-07 17:07:49 +01:00
use bitmasks
2023-07-02 21:49:25 +02:00
! ---
BEGIN_PROVIDER [integer, index_HF_psi_det]
implicit none
2023-07-02 21:49:25 +02:00
integer :: i, degree
do i = 1, N_det
call get_excitation_degree(HF_bitmask, psi_det(1,1,i), degree, N_int)
if(degree == 0) then
index_HF_psi_det = i
exit
endif
enddo
END_PROVIDER
! ---
subroutine diagonalize_CI_tc()
BEGIN_DOC
2023-07-02 21:49:25 +02:00
! Replace the coefficients of the |CI| states by the coefficients of the
! eigenstates of the |CI| matrix.
END_DOC
2023-07-02 21:49:25 +02:00
implicit none
integer :: i, j
do j = 1, N_states
do i = 1, N_det
psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j)
psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j)
enddo
enddo
SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho
2023-07-02 21:49:25 +02:00
end
2023-02-07 17:07:49 +01:00
2023-07-02 21:49:25 +02:00
! ---
2023-02-07 17:07:49 +01:00
BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth , (N_states) ]
&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth , (N_states) ]
&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth , (N_det,N_states)]
&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth , (N_det,N_states)]
&BEGIN_PROVIDER [double precision, s2_eigvec_tc_bi_orth , (N_states) ]
&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth , (N_states) ]
2023-02-07 17:07:49 +01:00
BEGIN_DOC
! eigenvalues, right and left eigenvectors of the transcorrelated Hamiltonian on the BI-ORTHO basis
END_DOC
implicit none
2023-04-03 14:55:02 +02:00
integer :: i, idx_dress, j, istate, k
2023-07-02 21:49:25 +02:00
integer :: i_good_state, i_other_state, i_state
integer :: n_real_tc_bi_orth_eigval_right, igood_r, igood_l
2023-02-07 17:07:49 +01:00
logical :: converged, dagger
2023-07-02 21:49:25 +02:00
double precision, parameter :: alpha = 0.1d0
integer, allocatable :: index_good_state_array(:)
integer, allocatable :: iorder(:)
logical, allocatable :: good_state_array(:)
double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:), leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
2023-04-10 20:41:16 +02:00
double precision, allocatable :: s2_values_tmp(:), H_prime(:,:), expect_e(:)
2023-07-02 21:49:25 +02:00
double precision, allocatable :: coef_hf_r(:),coef_hf_l(:)
double precision, allocatable :: Stmp(:,:)
2023-02-07 17:07:49 +01:00
PROVIDE N_det N_int
2023-07-02 21:49:25 +02:00
if(N_det .le. N_det_max_full) then
2023-04-13 13:03:10 +02:00
2023-07-02 21:49:25 +02:00
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))
2023-04-13 13:03:10 +02:00
2023-04-10 20:41:16 +02:00
H_prime(1:N_det,1:N_det) = htilde_matrix_elmt_bi_ortho(1:N_det,1:N_det)
2023-04-13 13:03:10 +02:00
if(s2_eig) then
H_prime(1:N_det,1:N_det) += alpha * S2_matrix_all_dets(1:N_det,1:N_det)
2023-07-02 21:49:25 +02:00
do j = 1, N_det
2023-04-13 13:03:10 +02:00
H_prime(j,j) = H_prime(j,j) - alpha*expected_s2
enddo
2023-02-07 17:07:49 +01:00
endif
2023-04-13 13:03:10 +02:00
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)
2023-04-10 20:41:16 +02:00
! 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)
2023-04-13 13:03:10 +02:00
2023-04-10 20:41:16 +02:00
allocate(index_good_state_array(N_det),good_state_array(N_det))
i_state = 0
good_state_array = .False.
2023-04-13 13:03:10 +02:00
if(s2_eig) then
if(only_expected_s2) then
do j = 1, N_det
2023-04-10 20:41:16 +02:00
! 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)
2023-04-13 13:03:10 +02:00
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
2023-04-10 20:41:16 +02:00
else ! s2_eig
2023-04-13 13:03:10 +02:00
allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(N_det))
do i = 1,N_det
2023-04-10 20:41:16 +02:00
iorder(i) = i
coef_hf_r(i) = -dabs(reigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
2023-04-13 13:03:10 +02:00
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
2023-04-10 20:41:16 +02:00
iorder(i) = i
coef_hf_l(i) = -dabs(leigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
2023-04-13 13:03:10 +02:00
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)
2023-04-10 20:41:16 +02:00
2023-04-13 13:03:10 +02:00
if(igood_r.ne.igood_l .and. igood_r.ne.1) then
2023-04-10 20:41:16 +02:00
print *,''
print *,'Warning, the left and right eigenvectors are "not the same" '
print *,'Warning, the ground state is not dominated by HF...'
print *,'State with largest RIGHT coefficient of HF ',igood_r
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)
2023-04-13 13:03:10 +02:00
endif
if(state_following_tc) then
2023-04-10 20:41:16 +02:00
print *,'Following the states with the largest coef on HF'
print *,'igood_r,igood_l',igood_r,igood_l
2023-04-13 13:03:10 +02:00
i = igood_r
2023-04-10 20:41:16 +02:00
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)
enddo
2023-04-13 13:03:10 +02:00
i = igood_l
2023-04-10 20:41:16 +02:00
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
2023-04-13 13:03:10 +02:00
else
2023-04-10 20:41:16 +02:00
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
2023-03-28 12:02:28 +02:00
enddo
2023-04-13 13:03:10 +02:00
endif
2023-03-28 12:02:28 +02:00
2023-02-07 17:07:49 +01:00
endif
2023-03-28 12:02:28 +02:00
2023-04-13 13:03:10 +02:00
else ! n_det > N_det_max_full
2023-03-28 12:02:28 +02:00
2023-02-07 17:07:49 +01:00
double precision, allocatable :: H_jj(:),vec_tmp(:,:)
external H_tc_u_0_opt
external H_tc_dagger_u_0_opt
2023-04-10 20:41:16 +02:00
external H_tc_s2_dagger_u_0_opt
external H_tc_s2_u_0_opt
2023-09-15 01:06:32 +02:00
external H_tc_s2_dagger_u_0_with_pure_three_omp
external H_tc_s2_u_0_with_pure_three_omp
2023-04-13 13:03:10 +02:00
2023-02-07 17:07:49 +01:00
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag))
2023-04-13 13:03:10 +02:00
2024-03-05 17:24:29 +01:00
! TODO : OPEN-MP
2023-02-07 17:07:49 +01:00
do i = 1, N_det
call htilde_mu_mat_opt_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
2023-02-07 17:07:49 +01:00
enddo
2023-04-13 13:03:10 +02:00
2023-04-12 17:10:06 +02:00
print*,'---------------------------------'
print*,'---------------------------------'
2023-02-07 17:07:49 +01:00
print*,'Computing the left-eigenvector '
2023-04-12 17:10:06 +02:00
print*,'---------------------------------'
print*,'---------------------------------'
2023-04-13 13:03:10 +02:00
!!!! Preparing the left-eigenvector
2023-02-07 17:07:49 +01:00
vec_tmp = 0.d0
do istate = 1, N_states
2023-04-13 13:03:10 +02:00
vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate)
2023-02-07 17:07:49 +01:00
enddo
do istate = N_states+1, n_states_diag
2023-04-13 13:03:10 +02:00
vec_tmp(istate,istate) = 1.d0
2023-02-07 17:07:49 +01:00
enddo
2023-04-12 22:44:14 +02:00
integer :: n_it_max,i_it
2023-04-12 17:10:06 +02:00
n_it_max = 1
converged = .False.
2023-04-12 22:44:14 +02:00
i_it = 0
2023-04-12 17:10:06 +02:00
do while (.not.converged)
2023-09-13 18:28:52 +02:00
if(.not.pure_three_body_h_tc)then
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)
else
2023-09-15 01:06:32 +02:00
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_with_pure_three_omp)
2023-09-13 18:28:52 +02:00
endif
2023-04-13 13:03:10 +02:00
i_it += 1
if(i_it .gt. 5) exit
2023-04-12 17:10:06 +02:00
enddo
2023-02-07 17:07:49 +01:00
do istate = 1, N_states
2023-04-13 13:03:10 +02:00
leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
2023-02-07 17:07:49 +01:00
enddo
2023-04-12 17:10:06 +02:00
print*,'---------------------------------'
print*,'---------------------------------'
2023-02-07 17:07:49 +01:00
print*,'Computing the right-eigenvector '
2023-04-12 17:10:06 +02:00
print*,'---------------------------------'
print*,'---------------------------------'
2023-04-13 13:03:10 +02:00
!!!! Preparing the right-eigenvector
2023-02-07 17:07:49 +01:00
vec_tmp = 0.d0
do istate = 1, N_states
2023-04-13 13:03:10 +02:00
vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate)
2023-02-07 17:07:49 +01:00
enddo
do istate = N_states+1, n_states_diag
2023-04-13 13:03:10 +02:00
vec_tmp(istate,istate) = 1.d0
2023-02-07 17:07:49 +01:00
enddo
2023-04-12 17:10:06 +02:00
converged = .False.
2023-04-12 22:44:14 +02:00
i_it = 0
2023-04-13 13:03:10 +02:00
do while (.not. converged)
2023-09-13 18:28:52 +02:00
if(.not.pure_three_body_h_tc)then
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)
else
2023-09-15 01:06:32 +02:00
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_with_pure_three_omp)
2023-09-13 18:28:52 +02:00
endif
2023-04-13 13:03:10 +02:00
i_it += 1
if(i_it .gt. 5) exit
2023-04-12 17:10:06 +02:00
enddo
2023-02-07 17:07:49 +01:00
do istate = 1, N_states
2023-04-13 13:03:10 +02:00
reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
2023-02-07 17:07:49 +01:00
enddo
deallocate(H_jj)
2023-04-13 13:03:10 +02:00
endif
2023-04-03 14:55:02 +02:00
2023-04-13 13:03:10 +02:00
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)
norm_ground_left_right_bi_orth = 0.d0
2023-04-03 14:55:02 +02:00
do i = 1, N_states
do j = 1, N_det
norm_ground_left_right_bi_orth(i) += leigvec_tc_bi_orth(j,i) * reigvec_tc_bi_orth(j,i)
2023-04-03 14:55:02 +02:00
enddo
2023-04-10 20:41:16 +02:00
print*,' state ', i
print*,' norm l/r = ', norm_ground_left_right_bi_orth(i)
2023-04-10 20:41:16 +02:00
print*,' <S2> = ', s2_eigvec_tc_bi_orth(i)
2023-04-03 14:55:02 +02:00
enddo
2023-05-13 12:18:36 +02:00
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)
! print*,'After diag'
! do i = 1, N_det! old version
! print*,'i',i,psi_l_coef_bi_ortho(i,1),psi_r_coef_bi_ortho(i,1)
! call debug_det(psi_det(1,1,i),N_int)
! enddo
2023-05-13 12:18:36 +02:00
2023-02-07 17:07:49 +01:00
END_PROVIDER
2023-04-15 00:59:22 +02:00
subroutine bi_normalize(u_l, u_r, n, ld, nstates)
BEGIN_DOC
2023-02-07 17:07:49 +01:00
!!!! Normalization of the scalar product of the left/right eigenvectors
2023-04-15 00:59:22 +02:00
END_DOC
implicit none
integer, intent(in) :: n, ld, nstates
2023-02-07 17:07:49 +01:00
double precision, intent(inout) :: u_l(ld,nstates), u_r(ld,nstates)
2023-04-15 00:59:22 +02:00
integer :: i, j
double precision :: accu, tmp
2023-02-07 17:07:49 +01:00
do i = 1, nstates
2023-04-15 00:59:22 +02:00
!!!! Normalization of right eigenvectors |Phi>
accu = 0.d0
2024-03-05 17:24:29 +01:00
! TODO: dot product lapack
2023-04-15 00:59:22 +02:00
do j = 1, n
accu += u_r(j,i) * u_r(j,i)
enddo
2023-02-07 17:07:49 +01:00
accu = 1.d0/dsqrt(accu)
2023-04-15 00:59:22 +02:00
print*,'accu_r = ',accu
do j = 1, n
u_r(j,i) *= accu
enddo
tmp = u_r(1,i) / dabs(u_r(1,i))
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
accu = 1.d0/dsqrt(-accu)
endif
tmp = (u_l(1,i) * u_r(1,i) )/dabs(u_l(1,i) * u_r(1,i))
do j = 1, n
u_l(j,i) *= accu * tmp
u_r(j,i) *= accu
enddo
2023-02-07 17:07:49 +01:00
enddo
2023-04-15 00:59:22 +02:00
2023-02-07 17:07:49 +01:00
end
2023-04-15 00:59:22 +02:00