mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-06 21:43:39 +01:00
added non-sym diag for tc-rpa
This commit is contained in:
Abdallah Ammar
parent
9fc4b6d63b
commit
6235c2015d
@ -346,7 +346,7 @@ subroutine davidson_general_ext_rout_nonsym_b1space(u_in, H_jj, energies, sze, N
|
|||||||
endif
|
endif
|
||||||
|
|
||||||
if(i_omax(l) .ne. l) then
|
if(i_omax(l) .ne. l) then
|
||||||
print *, ' !!! WARNONG !!!'
|
print *, ' !!! WARNING !!!'
|
||||||
print *, ' index of state', l, i_omax(l)
|
print *, ' index of state', l, i_omax(l)
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
|||||||
114
src/hartree_fock/print_scf_int.irp.f
Normal file
114
src/hartree_fock/print_scf_int.irp.f
Normal file
@ -0,0 +1,114 @@
|
|||||||
|
|
||||||
|
program print_scf_int
|
||||||
|
|
||||||
|
call main()
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine main()
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer :: i, j, k, l
|
||||||
|
|
||||||
|
print *, " Hcore:"
|
||||||
|
do j = 1, ao_num
|
||||||
|
do i = 1, ao_num
|
||||||
|
print *, i, j, ao_one_e_integrals(i,j)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
print *, " P:"
|
||||||
|
do j = 1, ao_num
|
||||||
|
do i = 1, ao_num
|
||||||
|
print *, i, j, SCF_density_matrix_ao_alpha(i,j)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
double precision :: integ, density_a, density_b, density
|
||||||
|
double precision :: J_scf(ao_num, ao_num)
|
||||||
|
double precision :: K_scf(ao_num, ao_num)
|
||||||
|
|
||||||
|
|
||||||
|
double precision, external :: get_ao_two_e_integral
|
||||||
|
PROVIDE ao_integrals_map
|
||||||
|
|
||||||
|
print *, " J:"
|
||||||
|
!do j = 1, ao_num
|
||||||
|
! do l = 1, ao_num
|
||||||
|
! do i = 1, ao_num
|
||||||
|
! do k = 1, ao_num
|
||||||
|
! ! < 1:k, 2:l | 1:i, 2:j >
|
||||||
|
! print *, '< k l | i j >', k, l, i, j
|
||||||
|
! print *, get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
!enddo
|
||||||
|
|
||||||
|
!do k = 1, ao_num
|
||||||
|
! do i = 1, ao_num
|
||||||
|
! do j = 1, ao_num
|
||||||
|
! do l = 1, ao_num
|
||||||
|
! ! ( 1:k, 1:i | 2:l, 2:j )
|
||||||
|
! print *, '(k i | l j)', k, i, l, j
|
||||||
|
! print *, get_ao_two_e_integral(l, j, k, i, ao_integrals_map)
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! print *, ''
|
||||||
|
! enddo
|
||||||
|
!enddo
|
||||||
|
|
||||||
|
J_scf = 0.d0
|
||||||
|
K_scf = 0.d0
|
||||||
|
do i = 1, ao_num
|
||||||
|
do k = 1, ao_num
|
||||||
|
do j = 1, ao_num
|
||||||
|
do l = 1, ao_num
|
||||||
|
|
||||||
|
density_a = SCF_density_matrix_ao_alpha(l,j)
|
||||||
|
density_b = SCF_density_matrix_ao_beta (l,j)
|
||||||
|
density = density_a + density_b
|
||||||
|
|
||||||
|
integ = get_ao_two_e_integral(l, j, k, i, ao_integrals_map)
|
||||||
|
J_scf(k,i) += density * integ
|
||||||
|
integ = get_ao_two_e_integral(l, i, k, j, ao_integrals_map)
|
||||||
|
K_scf(k,i) -= density_a * integ
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
print *, 'J x P'
|
||||||
|
do i = 1, ao_num
|
||||||
|
do k = 1, ao_num
|
||||||
|
print *, k, i, J_scf(k,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
print *, ''
|
||||||
|
print *, 'K x P'
|
||||||
|
do i = 1, ao_num
|
||||||
|
do k = 1, ao_num
|
||||||
|
print *, k, i, K_scf(k,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
print *, ''
|
||||||
|
print *, 'F in AO'
|
||||||
|
do i = 1, ao_num
|
||||||
|
do k = 1, ao_num
|
||||||
|
print *, k, i, Fock_matrix_ao(k,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
print *, ''
|
||||||
|
print *, 'F in MO'
|
||||||
|
do i = 1, ao_num
|
||||||
|
do k = 1, ao_num
|
||||||
|
print *, k, i, 2.d0 * Fock_matrix_mo_alpha(k,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
@ -1883,8 +1883,13 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
|
|||||||
enddo
|
enddo
|
||||||
accu_nd = dsqrt(accu_nd) / dble(m)
|
accu_nd = dsqrt(accu_nd) / dble(m)
|
||||||
|
|
||||||
|
if((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d) then
|
||||||
|
print *, ' non bi-orthogonal vectors !'
|
||||||
print *, ' accu_nd = ', accu_nd
|
print *, ' accu_nd = ', accu_nd
|
||||||
print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m)
|
print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m)
|
||||||
|
else
|
||||||
|
print *, ' vectors are bi-orthogonaly'
|
||||||
|
endif
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
@ -1994,10 +1999,13 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0)
|
|||||||
ii = ii + 1
|
ii = ii + 1
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
if(ii .eq. 0) then
|
if(ii .eq. 0) then
|
||||||
print*, ' WARNING: bi-orthogonality is lost but there is no degeneracies'
|
print*, ' WARNING: bi-orthogonality is lost but there is no degeneracies'
|
||||||
print*, ' rotations may change energy'
|
print*, ' rotations may change energy'
|
||||||
|
stop
|
||||||
endif
|
endif
|
||||||
|
|
||||||
print *, ii, ' type of degeneracies'
|
print *, ii, ' type of degeneracies'
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
@ -2018,17 +2026,18 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0)
|
|||||||
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
call dgemm( 'T', 'N', m, m, n, 1.d0 &
|
||||||
, L, size(L, 1), R, size(R, 1) &
|
, L, size(L, 1), R, size(R, 1) &
|
||||||
, 0.d0, S, size(S, 1) )
|
, 0.d0, S, size(S, 1) )
|
||||||
print*,'Overlap matrix '
|
|
||||||
accu_nd = 0.D0
|
print*, 'Overlap matrix '
|
||||||
|
accu_nd = 0.d0
|
||||||
do j = 1, m
|
do j = 1, m
|
||||||
write(*,'(100(F16.10,X))') S(1:m,j)
|
write(*,'(100(F16.10,X))') S(1:m,j)
|
||||||
do k = 1, m
|
do k = 1, m
|
||||||
if(j==k)cycle
|
if(j==k) cycle
|
||||||
accu_nd += dabs(S(j,k))
|
accu_nd += dabs(S(j,k))
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
print*,'accu_nd = ',accu_nd
|
print*,'accu_nd = ',accu_nd
|
||||||
! if(accu_nd .gt.1.d-10)then
|
! if(accu_nd .gt.1.d-10) then
|
||||||
! stop
|
! stop
|
||||||
! endif
|
! endif
|
||||||
do j = 1, m
|
do j = 1, m
|
||||||
@ -2036,13 +2045,15 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0)
|
|||||||
R0(1:n,i+j-1) = R(1:n,j)
|
R0(1:n,i+j-1) = R(1:n,j)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
deallocate(L, R,S)
|
deallocate(L, R, S)
|
||||||
|
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
end subroutine reorder_degen_eigvec
|
end subroutine reorder_degen_eigvec
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
|
subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
@ -2108,8 +2119,10 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
|
|||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
! call impose_orthog_svd(n, m, L)
|
|
||||||
call impose_orthog_svd(n, m, R)
|
call impose_orthog_svd(n, m, R)
|
||||||
|
L(:,:) = R(:,:)
|
||||||
|
|
||||||
|
!call impose_orthog_svd(n, m, L)
|
||||||
!call impose_orthog_GramSchmidt(n, m, L)
|
!call impose_orthog_GramSchmidt(n, m, L)
|
||||||
!call impose_orthog_GramSchmidt(n, m, R)
|
!call impose_orthog_GramSchmidt(n, m, R)
|
||||||
|
|
||||||
@ -2128,8 +2141,8 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
|
|||||||
!call bi_ortho_s_inv_half(m, L, R, S_inv_half)
|
!call bi_ortho_s_inv_half(m, L, R, S_inv_half)
|
||||||
!deallocate(S, S_inv_half)
|
!deallocate(S, S_inv_half)
|
||||||
|
|
||||||
call impose_biorthog_svd(n, m, L, R)
|
!call impose_biorthog_svd(n, m, L, R)
|
||||||
! call impose_biorthog_inverse(n, m, L, R)
|
!call impose_biorthog_inverse(n, m, L, R)
|
||||||
|
|
||||||
!call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
|
!call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
|
||||||
|
|
||||||
|
|||||||
@ -1,116 +0,0 @@
|
|||||||
|
|
||||||
BEGIN_PROVIDER [double precision, M_RPA, (2*nS_exc, 2*nS_exc)]
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
!
|
|
||||||
! full matrix for direct RPA calculation
|
|
||||||
! with the TC-Hamiltonian
|
|
||||||
!
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: ia, i, a, jb, j, b
|
|
||||||
double precision :: e(mo_num)
|
|
||||||
double precision, external :: Kronecker_delta
|
|
||||||
|
|
||||||
PROVIDE mo_tc_effec2e_int
|
|
||||||
PROVIDE Fock_matrix_tc_diag_mo_tot
|
|
||||||
|
|
||||||
e(1:mo_num) = Fock_matrix_tc_diag_mo_tot(1:mo_num)
|
|
||||||
|
|
||||||
|
|
||||||
! --- --- ---
|
|
||||||
! block A
|
|
||||||
|
|
||||||
ia = 0
|
|
||||||
do i = nC_orb+1, nO_orb
|
|
||||||
do a = nO_orb+1, mo_num-nR_orb
|
|
||||||
ia = ia + 1
|
|
||||||
|
|
||||||
jb = 0
|
|
||||||
do j = nC_orb+1, nO_orb
|
|
||||||
do b = nO_orb+1, mo_num-nR_orb
|
|
||||||
jb = jb + 1
|
|
||||||
|
|
||||||
M_RPA(ia,jb) = (e(a) - e(i)) * Kronecker_delta(i,j) * Kronecker_delta(a,b) + 2.d0 * mo_tc_effec2e_int(a,j,i,b)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
!
|
|
||||||
! --- --- ---
|
|
||||||
|
|
||||||
|
|
||||||
! --- --- ---
|
|
||||||
! block B
|
|
||||||
|
|
||||||
ia = 0
|
|
||||||
do i = nC_orb+1, nO_orb
|
|
||||||
do a = nO_orb+1, mo_num-nR_orb
|
|
||||||
ia = ia + 1
|
|
||||||
|
|
||||||
jb = nS_exc
|
|
||||||
do j = nC_orb+1, nO_orb
|
|
||||||
do b = nO_orb+1, mo_num-nR_orb
|
|
||||||
jb = jb + 1
|
|
||||||
|
|
||||||
M_RPA(ia,jb) = 2.d0 * mo_tc_effec2e_int(a,b,i,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
!
|
|
||||||
! --- --- ---
|
|
||||||
|
|
||||||
|
|
||||||
! --- --- ---
|
|
||||||
! block C
|
|
||||||
|
|
||||||
ia = nS_exc
|
|
||||||
do i = nC_orb+1, nO_orb
|
|
||||||
do a = nO_orb+1, mo_num-nR_orb
|
|
||||||
ia = ia + 1
|
|
||||||
|
|
||||||
jb = 0
|
|
||||||
do j = nC_orb+1, nO_orb
|
|
||||||
do b = nO_orb+1, mo_num-nR_orb
|
|
||||||
jb = jb + 1
|
|
||||||
|
|
||||||
M_RPA(ia,jb) = 2.d0 * mo_tc_effec2e_int(i,j,a,b)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
!
|
|
||||||
! --- --- ---
|
|
||||||
|
|
||||||
|
|
||||||
! --- --- ---
|
|
||||||
! block D
|
|
||||||
|
|
||||||
ia = nS_exc
|
|
||||||
do i = nC_orb+1, nO_orb
|
|
||||||
do a = nO_orb+1, mo_num-nR_orb
|
|
||||||
ia = ia + 1
|
|
||||||
|
|
||||||
jb = nS_exc
|
|
||||||
do j = nC_orb+1, nO_orb
|
|
||||||
do b = nO_orb+1, mo_num-nR_orb
|
|
||||||
jb = jb + 1
|
|
||||||
|
|
||||||
M_RPA(ia,jb) = (e(a) - e(i)) * Kronecker_delta(i,j) * Kronecker_delta(a,b) + 2.d0 * mo_tc_effec2e_int(i,b,a,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
!
|
|
||||||
! --- --- ---
|
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
@ -1,39 +0,0 @@
|
|||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [double precision, mo_tc_effec2e_int, (mo_num, mo_num, mo_num, mo_num)]
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
!
|
|
||||||
! mo_tc_effec2e_int(p,q,s,t) = < p q| V(12) | s t > + \sum_i < p q i | L(123)| s t i >
|
|
||||||
!
|
|
||||||
! the potential V(12) contains ALL TWO-E CONTRIBUTION OF THE TC-HAMILTONIAN
|
|
||||||
!
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, k, l, ii
|
|
||||||
double precision :: integral
|
|
||||||
|
|
||||||
PROVIDE mo_bi_ortho_tc_two_e_chemist
|
|
||||||
|
|
||||||
do j = 1, mo_num
|
|
||||||
do i = 1, mo_num
|
|
||||||
do l = 1, mo_num
|
|
||||||
do k = 1, mo_num
|
|
||||||
mo_tc_effec2e_int(k,l,i,j) = mo_bi_ortho_tc_two_e_chemist(k,i,l,j)
|
|
||||||
|
|
||||||
do ii = 1, elec_alpha_num
|
|
||||||
call give_integrals_3_body_bi_ort(k, l, ii, i, j, ii, integral)
|
|
||||||
mo_tc_effec2e_int(k,l,i,j) -= 2.d0 * integral
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
FREE mo_bi_ortho_tc_two_e_chemist
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
@ -1,181 +0,0 @@
|
|||||||
program tc_rpa
|
|
||||||
|
|
||||||
BEGIN_DOC
|
|
||||||
!
|
|
||||||
!
|
|
||||||
!
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
my_grid_becke = .True.
|
|
||||||
PROVIDE tc_grid1_a tc_grid1_r
|
|
||||||
my_n_pt_r_grid = tc_grid1_r
|
|
||||||
my_n_pt_a_grid = tc_grid1_a
|
|
||||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
|
||||||
|
|
||||||
if(j1b_type .ge. 100) then
|
|
||||||
my_extra_grid_becke = .True.
|
|
||||||
PROVIDE tc_grid2_a tc_grid2_r
|
|
||||||
my_n_pt_r_extra_grid = tc_grid2_r
|
|
||||||
my_n_pt_a_extra_grid = tc_grid2_a
|
|
||||||
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
|
|
||||||
|
|
||||||
call write_int(6, my_n_pt_r_extra_grid, 'radial internal grid over')
|
|
||||||
call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over')
|
|
||||||
endif
|
|
||||||
|
|
||||||
|
|
||||||
call main()
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
subroutine main()
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, j, n
|
|
||||||
integer :: n_good, n_real_eigv
|
|
||||||
double precision :: thr_cpx, thr_d, thr_nd
|
|
||||||
double precision :: accu_d, accu_nd
|
|
||||||
integer, allocatable :: list_good(:), iorder(:)
|
|
||||||
double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
|
|
||||||
double precision, allocatable :: Omega_p(:), Reigvec_p(:,:), Leigvec_p(:,:)
|
|
||||||
double precision, allocatable :: Omega_m(:), Reigvec_m(:,:), Leigvec_m(:,:)
|
|
||||||
double precision, allocatable :: S(:,:)
|
|
||||||
|
|
||||||
PROVIDE M_RPA
|
|
||||||
|
|
||||||
print *, ' '
|
|
||||||
print *, ' Computing left/right eigenvectors for TC-RPA ...'
|
|
||||||
print *, ' '
|
|
||||||
|
|
||||||
|
|
||||||
n = 2 * nS_exc
|
|
||||||
|
|
||||||
thr_cpx = 1d-7
|
|
||||||
thr_d = 1d-07
|
|
||||||
thr_nd = 1d-07
|
|
||||||
|
|
||||||
|
|
||||||
allocate(WR(n), WI(n), VL(n,n), VR(n,n))
|
|
||||||
call lapack_diag_non_sym(n, M_RPA, WR, WI, VL, VR)
|
|
||||||
FREE M_RPA
|
|
||||||
|
|
||||||
print *, ' excitation energies:'
|
|
||||||
do i = 1, nS_exc
|
|
||||||
write(*, '(I3, X, 1000(F16.10,X))') i, WR(i), WI(i)
|
|
||||||
if(dabs(WI(i)) .gt. thr_cpx) then
|
|
||||||
print *, ' WARNING ! IMAGINARY EIGENVALUES !!!'
|
|
||||||
write(*, '(1000(F16.10,X))') WR(i), WI(i+1)
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
print *, ' '
|
|
||||||
print *, ' desexcitation energies:'
|
|
||||||
do i = nS_exc+1, n
|
|
||||||
write(*, '(I3, X, 1000(F16.10,X))') i, WR(i), WI(i)
|
|
||||||
if(dabs(WI(i)) .gt. thr_cpx) then
|
|
||||||
print *, ' WARNING ! IMAGINARY EIGENVALUES !!!'
|
|
||||||
write(*, '(1000(F16.10,X))') WR(i), WI(i+1)
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
|
|
||||||
! track & sort the real eigenvalues
|
|
||||||
|
|
||||||
n_good = 0
|
|
||||||
do i = 1, nS_exc
|
|
||||||
if(dabs(WI(i)) .lt. thr_cpx) then
|
|
||||||
if(dabs(WI(nS_exc+i)) .lt. thr_cpx) then
|
|
||||||
n_good += 1
|
|
||||||
endif
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
n_real_eigv = n_good
|
|
||||||
|
|
||||||
print *, ' '
|
|
||||||
print *, ' nb of real eigenvalues = ', n_real_eigv
|
|
||||||
print *, ' total nb of eigenvalues = ', nS_exc
|
|
||||||
|
|
||||||
allocate(Omega_p(n_real_eigv), Reigvec_p(n,n_real_eigv), Leigvec_p(n,n_real_eigv))
|
|
||||||
allocate(Omega_m(n_real_eigv), Reigvec_m(n,n_real_eigv), Leigvec_m(n,n_real_eigv))
|
|
||||||
|
|
||||||
n_good = 0
|
|
||||||
do i = 1, nS_exc
|
|
||||||
if(dabs(WI(i)) .lt. thr_cpx) then
|
|
||||||
if(dabs(WI(nS_exc+i)) .lt. thr_cpx) then
|
|
||||||
n_good += 1
|
|
||||||
|
|
||||||
Omega_p(n_good) = WR(i)
|
|
||||||
do j = 1, n
|
|
||||||
Reigvec_p(j,n_good) = VR(j,n_good)
|
|
||||||
Leigvec_p(j,n_good) = VL(j,n_good)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
Omega_m(n_good) = WR(nS_exc+i)
|
|
||||||
do j = 1, n
|
|
||||||
Reigvec_m(j,n_good) = VR(j,nS_exc+n_good)
|
|
||||||
Leigvec_m(j,n_good) = VL(j,nS_exc+n_good)
|
|
||||||
enddo
|
|
||||||
endif
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
|
|
||||||
deallocate(WR, WI, VL, VR)
|
|
||||||
|
|
||||||
|
|
||||||
! check bi-orthogonality
|
|
||||||
|
|
||||||
! first block
|
|
||||||
|
|
||||||
allocate(S(n_real_eigv,n_real_eigv))
|
|
||||||
|
|
||||||
call check_biorthog(n, n_real_eigv, Leigvec_p, Reigvec_p, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
|
|
||||||
print *, ' accu_d = ', accu_d
|
|
||||||
print *, ' accu_nd = ', accu_nd
|
|
||||||
|
|
||||||
if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d)) then
|
|
||||||
print *, ' RPA first-block eigenvectors are normalized and bi-orthogonalized'
|
|
||||||
else
|
|
||||||
print *, ' RPA first-block eigenvectors are neither normalized nor bi-orthogonalized'
|
|
||||||
|
|
||||||
call reorder_degen_eigvec(n, Omega_p, Leigvec_p, Reigvec_p)
|
|
||||||
call impose_biorthog_degen_eigvec(n, Omega_p, Leigvec_p, Reigvec_p)
|
|
||||||
|
|
||||||
call check_biorthog(n, n_real_eigv, Leigvec_p, Reigvec_p, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
|
|
||||||
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .gt. thr_d) ) then
|
|
||||||
call check_biorthog_binormalize(n, n_real_eigv, Leigvec_p, Reigvec_p, thr_d, thr_nd, .true.)
|
|
||||||
endif
|
|
||||||
call check_biorthog(n, n_real_eigv, Leigvec_p, Reigvec_p, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
|
|
||||||
endif
|
|
||||||
|
|
||||||
|
|
||||||
! second block
|
|
||||||
|
|
||||||
call check_biorthog(n, n_real_eigv, Leigvec_m, Reigvec_m, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
|
|
||||||
print *, ' accu_d = ', accu_d
|
|
||||||
print *, ' accu_nd = ', accu_nd
|
|
||||||
|
|
||||||
if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d)) then
|
|
||||||
print *, ' RPA first-block eigenvectors are normalized and bi-orthogonalized'
|
|
||||||
else
|
|
||||||
print *, ' RPA first-block eigenvectors are neither normalized nor bi-orthogonalized'
|
|
||||||
|
|
||||||
call reorder_degen_eigvec(n, Omega_m, Leigvec_m, Reigvec_m)
|
|
||||||
call impose_biorthog_degen_eigvec(n, Omega_m, Leigvec_m, Reigvec_m)
|
|
||||||
|
|
||||||
call check_biorthog(n, n_real_eigv, Leigvec_m, Reigvec_m, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
|
|
||||||
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .gt. thr_d) ) then
|
|
||||||
call check_biorthog_binormalize(n, n_real_eigv, Leigvec_m, Reigvec_m, thr_d, thr_nd, .true.)
|
|
||||||
endif
|
|
||||||
call check_biorthog(n, n_real_eigv, Leigvec_m, Reigvec_m, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
|
|
||||||
endif
|
|
||||||
|
|
||||||
deallocate(S)
|
|
||||||
|
|
||||||
return
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
@ -600,3 +600,43 @@ end function Kronecker_delta
|
|||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
|
subroutine diagonalize_sym_matrix(N, A, e)
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
!
|
||||||
|
! Diagonalize a symmetric matrix
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
integer, intent(in) :: N
|
||||||
|
double precision, intent(inout) :: A(N,N)
|
||||||
|
double precision, intent(out) :: e(N)
|
||||||
|
|
||||||
|
integer :: lwork, info
|
||||||
|
double precision, allocatable :: work(:)
|
||||||
|
|
||||||
|
allocate(work(1))
|
||||||
|
|
||||||
|
lwork = -1
|
||||||
|
call dsyev('V', 'U', N, A, N, e, work, lwork, info)
|
||||||
|
lwork = int(work(1))
|
||||||
|
|
||||||
|
deallocate(work)
|
||||||
|
|
||||||
|
allocate(work(lwork))
|
||||||
|
|
||||||
|
call dsyev('V', 'U', N, A, N, e, work, lwork, info)
|
||||||
|
deallocate(work)
|
||||||
|
|
||||||
|
if(info /= 0) then
|
||||||
|
print*,'Problem in diagonalize_sym_matrix (dsyev)!!'
|
||||||
|
endif
|
||||||
|
|
||||||
|
end subroutine diagonalize_sym_matrix
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user