mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-21 11:03:29 +01:00
v0 of tc-dRPA
This commit is contained in:
Abdallah Ammar
parent
f437a84adb
commit
9fc4b6d63b
@ -1944,6 +1944,7 @@ subroutine check_orthog(n, m, V, accu_d, accu_nd, S)
|
||||
end subroutine check_orthog
|
||||
|
||||
! ---
|
||||
|
||||
subroutine reorder_degen_eigvec(n, e0, L0, R0)
|
||||
|
||||
implicit none
|
||||
@ -1953,7 +1954,7 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0)
|
||||
double precision, intent(inout) :: L0(n,n), R0(n,n)
|
||||
|
||||
logical :: complex_root
|
||||
integer :: i, j, k, m
|
||||
integer :: i, j, k, m, ii
|
||||
double precision :: ei, ej, de, de_thr
|
||||
double precision :: accu_d, accu_nd
|
||||
integer, allocatable :: deg_num(:)
|
||||
@ -1986,11 +1987,18 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
ii = 0
|
||||
do i = 1, n
|
||||
if(deg_num(i) .gt. 1) then
|
||||
print *, ' degen on', i, deg_num(i), e0(i)
|
||||
ii = ii + 1
|
||||
endif
|
||||
enddo
|
||||
if(ii .eq. 0) then
|
||||
print*, ' WARNING: bi-orthogonality is lost but there is no degeneracies'
|
||||
print*, ' rotations may change energy'
|
||||
endif
|
||||
print *, ii, ' type of degeneracies'
|
||||
|
||||
! ---
|
||||
|
||||
@ -2013,7 +2021,7 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0)
|
||||
print*,'Overlap matrix '
|
||||
accu_nd = 0.D0
|
||||
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
|
||||
if(j==k)cycle
|
||||
accu_nd += dabs(S(j,k))
|
||||
|
||||
38
src/tc_bi_ortho/ORBITALS.irp.f
Normal file
38
src/tc_bi_ortho/ORBITALS.irp.f
Normal file
@ -0,0 +1,38 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [integer, nC_orb]
|
||||
&BEGIN_PROVIDER [integer, nO_orb]
|
||||
&BEGIN_PROVIDER [integer, nV_orb]
|
||||
&BEGIN_PROVIDER [integer, nR_orb]
|
||||
&BEGIN_PROVIDER [integer, nS_exc]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! nC_orb = number of core orbitals
|
||||
! nO_orb = number of occupied orbitals
|
||||
! nV_orb = number of virtual orbitals
|
||||
! nR_orb = number of Rydberg orbitals
|
||||
! nS_exc = number of single excitation
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
nC_orb = 0
|
||||
nO_orb = elec_beta_num - nC_orb
|
||||
nV_orb = mo_num - (nC_orb + nO_orb)
|
||||
nR_orb = 0
|
||||
nS_exc = (nO_orb-nC_orb) * (nV_orb-nR_orb)
|
||||
|
||||
print *, ' nC_orb = ', nC_orb
|
||||
print *, ' nO_orb = ', nO_orb
|
||||
print *, ' nV_orb = ', nV_orb
|
||||
print *, ' nR_orb = ', nR_orb
|
||||
print *, ' nS_exc = ', nS_exc
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
|
||||
116
src/tc_bi_ortho/drpa_matrix.irp.f
Normal file
116
src/tc_bi_ortho/drpa_matrix.irp.f
Normal file
@ -0,0 +1,116 @@
|
||||
|
||||
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
|
||||
|
||||
|
||||
39
src/tc_bi_ortho/tc_effect_int.irp.f
Normal file
39
src/tc_bi_ortho/tc_effect_int.irp.f
Normal file
@ -0,0 +1,39 @@
|
||||
|
||||
|
||||
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
|
||||
|
||||
! ---
|
||||
|
||||
181
src/tc_bi_ortho/tc_rpa.irp.f
Normal file
181
src/tc_bi_ortho/tc_rpa.irp.f
Normal file
@ -0,0 +1,181 @@
|
||||
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
|
||||
|
||||
! ---
|
||||
|
||||
@ -580,4 +580,23 @@ end function is_same_spin
|
||||
|
||||
! ---
|
||||
|
||||
function Kronecker_delta(i, j) result(delta)
|
||||
|
||||
BEGIN_DOC
|
||||
! Kronecker Delta
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i, j
|
||||
double precision :: delta
|
||||
|
||||
if(i == j) then
|
||||
delta = 1.d0
|
||||
else
|
||||
delta = 0.d0
|
||||
endif
|
||||
|
||||
end function Kronecker_delta
|
||||
|
||||
! ---
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user