LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00
Code Releases Activity

merged version: ok

Browse Source
This commit is contained in:
AbdAmmar 2022-10-31 19:14:27 +01:00
parent 0bea180fc6
commit 8bed36e93e
11 changed files with 1442 additions and 562 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

114
src/bi_ortho_mos/overlap.irp.f
View File

@ -76,67 +76,85 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, overlap_mo_r, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, overlap_mo_l, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! overlap_mo_r_mo(j,i) = <MO_i|MO_R_j>
END_DOC
integer :: i,j,p,q
overlap_mo_r= 0.d0
overlap_mo_l= 0.d0
do i = 1, mo_num
do j = 1, mo_num
do p = 1, ao_num
do q = 1, ao_num
overlap_mo_r(j,i) += mo_r_coef(q,i) * mo_r_coef(p,j) * ao_overlap(q,p)
overlap_mo_l(j,i) += mo_l_coef(q,i) * mo_l_coef(p,j) * ao_overlap(q,p)
enddo
BEGIN_DOC
! overlap_mo_r_mo(j,i) = <MO_i|MO_R_j>
END_DOC
implicit none
integer :: i, j, p, q
overlap_mo_r = 0.d0
overlap_mo_l = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do p = 1, ao_num
do q = 1, ao_num
overlap_mo_r(j,i) += mo_r_coef(q,i) * mo_r_coef(p,j) * ao_overlap(q,p)
overlap_mo_l(j,i) += mo_l_coef(q,i) * mo_l_coef(p,j) * ao_overlap(q,p)
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, overlap_mo_r_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, overlap_mo_l_mo, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! overlap_mo_r_mo(j,i) = <MO_j|MO_R_i>
END_DOC
integer :: i,j,p,q
overlap_mo_r_mo = 0.d0
overlap_mo_l_mo = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do p = 1, ao_num
do q = 1, ao_num
overlap_mo_r_mo(j,i) += mo_coef(p,j) * mo_r_coef(q,i) * ao_overlap(q,p)
overlap_mo_l_mo(j,i) += mo_coef(p,j) * mo_l_coef(q,i) * ao_overlap(q,p)
enddo
BEGIN_DOC
! overlap_mo_r_mo(j,i) = <MO_j|MO_R_i>
END_DOC
implicit none
integer :: i, j, p, q
overlap_mo_r_mo = 0.d0
overlap_mo_l_mo = 0.d0
do i = 1, mo_num
do j = 1, mo_num
do p = 1, ao_num
do q = 1, ao_num
overlap_mo_r_mo(j,i) += mo_coef(p,j) * mo_r_coef(q,i) * ao_overlap(q,p)
overlap_mo_l_mo(j,i) += mo_coef(p,j) * mo_l_coef(q,i) * ao_overlap(q,p)
enddo
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, angle_left_right, (mo_num)]
&BEGIN_PROVIDER [ double precision, max_angle_left_right]
implicit none
BEGIN_DOC
BEGIN_DOC
! angle_left_right(i) = angle between the left-eigenvector chi_i and the right-eigenvector phi_i
END_DOC
integer :: i,j
double precision :: left,right,arg
do i = 1, mo_num
left = overlap_mo_l(i,i)
right = overlap_mo_r(i,i)
arg = min(overlap_bi_ortho(i,i)/(left*right),1.d0)
arg = max(arg,-1.d0)
angle_left_right(i) = dacos(arg) * 180.d0/dacos(-1.d0)
enddo
double precision :: angle(mo_num)
angle(1:mo_num) = dabs(angle_left_right(1:mo_num))
max_angle_left_right = maxval(angle)
END_DOC
implicit none
integer :: i, j
double precision :: left, right, arg
double precision :: angle(mo_num)
do i = 1, mo_num
left = overlap_mo_l(i,i)
right = overlap_mo_r(i,i)
arg = min(overlap_bi_ortho(i,i)/(left*right),1.d0)
arg = max(arg, -1.d0)
angle_left_right(i) = dacos(arg) * 180.d0/dacos(-1.d0)
enddo
angle(1:mo_num) = dabs(angle_left_right(1:mo_num))
max_angle_left_right = maxval(angle)
END_PROVIDER
! ---

3
src/cipsi/selection.irp.f
View File

@ -732,7 +732,10 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
if (delta_E < 0.d0) then
tmp = -tmp
endif
!e_pert(istate) = alpha_h_psi * alpha_h_psi / (E0(istate) - Hii)
e_pert(istate) = 0.5d0 * (tmp - delta_E)
if (dabs(alpha_h_psi) > 1.d-4) then
coef(istate) = e_pert(istate) / alpha_h_psi
else

27
src/non_hermit_dav/biorthog.irp.f
View File

@ -252,7 +252,7 @@ end subroutine non_hrmt_real_diag_new
! ---
subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, eigval)
BEGIN_DOC
!
@ -266,13 +266,14 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: A(n,n)
double precision, intent(in) :: thr_d, thr_nd
integer, intent(out) :: n_real_eigv
double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
integer :: i, j
integer :: n_good
double precision :: thr, thr_cut, thr_diag, thr_norm
double precision :: accu_d, accu_nd, thr_d, thr_nd
double precision :: accu_d, accu_nd
integer, allocatable :: list_good(:), iorder(:)
double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
@ -329,7 +330,7 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
n_good = 0
thr = 1.d-5
do i = 1, n
print*, 'Re(i) + Im(i)', WR(i), WI(i)
print*, 'Re(i) + Im(i)', WR(i), WI(i)
if(dabs(WI(i)) .lt. thr) then
n_good += 1
else
@ -393,22 +394,22 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
! -------------------------------------------------------------------------------------
! check bi-orthogonality
thr_d = 1d-10 ! -7
thr_nd = 1d-10 ! -7
thr_diag = 10.d0
thr_norm = 1d+10
allocate( S(n_real_eigv,n_real_eigv) )
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, .false.)
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .lt. thr_d) ) then
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d) ) then
print *, ' lapack vectors are normalized and bi-orthogonalized'
deallocate(S)
return
elseif( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .gt. thr_d) ) then
elseif( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .gt. thr_d) ) then
print *, ' lapack vectors are not normalized but bi-orthogonalized'
call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, .true.)
call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.)
call check_EIGVEC(n, n, A, eigval, leigvec, reigvec, thr_diag, thr_norm, .true.)
@ -427,17 +428,17 @@ subroutine non_hrmt_bieig(n, A, leigvec, reigvec, n_real_eigv, eigval)
call impose_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)
!call impose_orthog_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)
!call impose_orthog_biorthog_degen_eigvec(n, thr_d, thr_nd, eigval, leigvec, reigvec)
!call impose_unique_biorthog_degen_eigvec(n, eigval, mo_coef, ao_overlap, leigvec, reigvec)
! ---
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, .false.)
call check_biorthog(n, n_real_eigv, leigvec, reigvec, 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, reigvec, .true.)
call check_biorthog_binormalize(n, n_real_eigv, leigvec, reigvec, thr_d, thr_nd, .true.)
endif
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, .true.)
call check_biorthog(n, n_real_eigv, leigvec, reigvec, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
!call impose_biorthog_qr(n, n_real_eigv, leigvec, reigvec)
!call impose_biorthog_lu(n, n_real_eigv, leigvec, reigvec)

361
src/non_hermit_dav/lapack_diag_non_hermit.irp.f
View File

@ -573,21 +573,22 @@ end subroutine non_hrmt_general_real_diag
! ---
subroutine impose_biorthog_qr(m, n, Vl, Vr)
subroutine impose_biorthog_qr(m, n, thr_d, thr_nd, Vl, Vr)
implicit none
integer, intent(in) :: m, n
integer, intent(in) :: m, n
double precision, intent(in) :: thr_d, thr_nd
double precision, intent(inout) :: Vl(m,n), Vr(m,n)
integer :: i, j
integer :: LWORK, INFO
double precision :: accu_nd, accu_d, thr_nd, thr_d
double precision :: accu_nd, accu_d
double precision, allocatable :: TAU(:), WORK(:)
double precision, allocatable :: S(:,:), R(:,:), tmp(:,:)
! ---
call check_biorthog_binormalize(m, n, Vl, Vr, .false.)
call check_biorthog_binormalize(m, n, Vl, Vr, thr_d, thr_nd, .false.)
! ---
@ -609,9 +610,7 @@ subroutine impose_biorthog_qr(m, n, Vl, Vr)
enddo
accu_nd = dsqrt(accu_nd)
thr_d = 1d-10
thr_nd = 1d-08
if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n)) .lt. thr_d)) then
if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n))/dble(n) .lt. thr_d)) then
print *, ' bi-orthogonal vectors without QR !'
deallocate(S)
return
@ -1011,7 +1010,7 @@ subroutine check_degen(n, m, eigval, leigvec, reigvec)
double precision :: ei, ej, de, de_thr, accu_nd
double precision, allocatable :: S(:,:)
de_thr = 1d-7
de_thr = 1d-6
do i = 1, m-1
ei = eigval(i)
@ -1082,7 +1081,7 @@ subroutine impose_weighted_orthog_svd(n, m, W, C)
double precision, allocatable :: S(:,:), tmp(:,:)
double precision, allocatable :: U(:,:), Vt(:,:), D(:)
print *, ' apply SVD to orthogonalize vectors'
print *, ' apply SVD to orthogonalize & normalize weighted vectors'
! ---
@ -1097,7 +1096,7 @@ subroutine impose_weighted_orthog_svd(n, m, W, C)
, 0.d0, S, size(S, 1) )
deallocate(tmp)
print *, ' eigenvec overlap bef SVD: '
print *, ' overlap bef SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
@ -1160,7 +1159,7 @@ subroutine impose_weighted_orthog_svd(n, m, W, C)
, 0.d0, S, size(S, 1) )
deallocate(tmp)
print *, ' eigenvec overlap aft SVD: '
print *, ' overlap aft SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
@ -1185,7 +1184,7 @@ subroutine impose_orthog_svd(n, m, C)
double precision, allocatable :: S(:,:), tmp(:,:)
double precision, allocatable :: U(:,:), Vt(:,:), D(:)
print *, ' apply SVD to orthogonalize vectors'
print *, ' apply SVD to orthogonalize & normalize vectors'
! ---
@ -1379,7 +1378,7 @@ subroutine impose_orthog_GramSchmidt(n, m, C)
double precision, allocatable :: S(:,:)
print *, ''
print *, ' apply Gram-Schmidt to orthogonalize vectors'
print *, ' apply Gram-Schmidt to orthogonalize & normalize vectors'
print *, ''
! ---
@ -1663,22 +1662,19 @@ end subroutine get_halfinv_svd
! ---
subroutine check_biorthog_binormalize(n, m, Vl, Vr, stop_ifnot)
subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot)
implicit none
integer, intent(in) :: n, m
logical, intent(in) :: stop_ifnot
double precision, intent(in) :: thr_d, thr_nd
double precision, intent(inout) :: Vl(n,m), Vr(n,m)
integer :: i, j
double precision :: thr_d, thr_nd
double precision :: accu_d, accu_nd, s_tmp
double precision, allocatable :: S(:,:)
thr_d = 1d-6
thr_nd = 1d-7
print *, ' check bi-orthonormality'
! ---
@ -1713,7 +1709,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, stop_ifnot)
endif
enddo
enddo
accu_nd = dsqrt(accu_nd)
accu_nd = dsqrt(accu_nd) / dble(m)
print*, ' diag acc: ', accu_d
print*, ' nondiag acc: ', accu_nd
@ -1755,7 +1751,7 @@ subroutine check_biorthog_binormalize(n, m, Vl, Vr, stop_ifnot)
endif
enddo
enddo
accu_nd = dsqrt(accu_nd)
accu_nd = dsqrt(accu_nd) / dble(m)
print *, ' diag acc: ', accu_d
print *, ' nondiag acc: ', accu_nd
@ -1774,22 +1770,19 @@ end subroutine check_biorthog_binormalize
! ---
subroutine check_weighted_biorthog(n, m, W, Vl, Vr, accu_d, accu_nd, S, stop_ifnot)
subroutine check_weighted_biorthog(n, m, W, Vl, Vr, thr_d, thr_nd, accu_d, accu_nd, S, stop_ifnot)
implicit none
integer, intent(in) :: n, m
double precision, intent(in) :: Vl(n,m), Vr(n,m), W(n,n)
double precision, intent(in) :: thr_d, thr_nd
logical, intent(in) :: stop_ifnot
double precision, intent(out) :: accu_d, accu_nd, S(m,m)
integer :: i, j
double precision :: thr_d, thr_nd
double precision, allocatable :: SS(:,:), tmp(:,:)
thr_d = 1d-6
thr_nd = 1d-08
print *, ' check weighted bi-orthogonality'
! ---
@ -1841,22 +1834,19 @@ end subroutine check_weighted_biorthog
! ---
subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, stop_ifnot)
subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ifnot)
implicit none
integer, intent(in) :: n, m
double precision, intent(in) :: Vl(n,m), Vr(n,m)
logical, intent(in) :: stop_ifnot
double precision, intent(in) :: thr_d, thr_nd
double precision, intent(out) :: accu_d, accu_nd, S(m,m)
integer :: i, j
double precision :: thr_d, thr_nd
double precision, allocatable :: SS(:,:)
thr_d = 1d-6
thr_nd = 1d-08
print *, ' check bi-orthogonality'
! ---
@ -1880,7 +1870,7 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, stop_ifnot)
endif
enddo
enddo
accu_nd = dsqrt(accu_nd)
accu_nd = dsqrt(accu_nd) / dble(m)
print *, ' accu_nd = ', accu_nd
print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m)
@ -2029,7 +2019,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
call impose_biorthog_svd(n, m, L, R)
!call impose_biorthog_qr(n, m, L, R)
!call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
! ---
@ -2047,11 +2037,12 @@ end subroutine impose_biorthog_degen_eigvec
! ---
subroutine impose_orthog_biorthog_degen_eigvec(n, e0, L0, R0)
subroutine impose_orthog_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, L0, R0)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: thr_d, thr_nd
double precision, intent(in) :: e0(n)
double precision, intent(inout) :: L0(n,n), R0(n,n)
@ -2116,12 +2107,12 @@ subroutine impose_orthog_biorthog_degen_eigvec(n, e0, L0, R0)
! ---
call impose_biorthog_qr(n, m, L, R)
call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
allocate(S(m,m))
call check_biorthog(n, m, L, R, accu_d, accu_nd, S, .true.)
!call check_biorthog(n, m, L, L, accu_d, accu_nd, S, .true.)
!call check_biorthog(n, m, R, R, accu_d, accu_nd, S, .false.)
call check_biorthog(n, m, L, R, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
!call check_biorthog(n, m, L, L, accu_d, accu_nd, S, thr_d, thr_nd, .true.)
!call check_biorthog(n, m, R, R, accu_d, accu_nd, S, thr_d, thr_nd, .false.)
deallocate(S)
! ---
@ -2140,11 +2131,12 @@ end subroutine impose_orthog_biorthog_degen_eigvec
! ---
subroutine impose_unique_biorthog_degen_eigvec(n, e0, C0, W0, L0, R0)
subroutine impose_unique_biorthog_degen_eigvec(n, thr_d, thr_nd, e0, C0, W0, L0, R0)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: thr_d, thr_nd
double precision, intent(in) :: e0(n), W0(n,n), C0(n,n)
double precision, intent(inout) :: L0(n,n), R0(n,n)
@ -2255,7 +2247,7 @@ subroutine impose_unique_biorthog_degen_eigvec(n, e0, C0, W0, L0, R0)
call get_inv_half_nonsymmat_diago(S, m, S_inv_half, complex_root)
if(complex_root)then
call impose_biorthog_svd(n, m, L, R)
!call impose_biorthog_qr(n, m, L, R)
!call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
else
call bi_ortho_s_inv_half(m, L, R, S_inv_half)
endif
@ -2502,8 +2494,288 @@ end subroutine impose_biorthog_svd
! ---
subroutine impose_weighted_biorthog_qr(m, n, thr_d, thr_nd, Vl, W, Vr)
subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
implicit none
integer, intent(in) :: m, n
double precision, intent(in) :: thr_d, thr_nd
double precision, intent(inout) :: Vl(m,n), W(m,m), Vr(m,n)
integer :: i, j
integer :: LWORK, INFO
double precision :: accu_nd, accu_d
double precision, allocatable :: TAU(:), WORK(:)
double precision, allocatable :: S(:,:), R(:,:), tmp(:,:), Stmp(:,:)
call check_weighted_biorthog_binormalize(m, n, Vl, W, Vr, thr_d, thr_nd, .false.)
! ---
allocate(Stmp(n,m), S(n,n))
call dgemm( 'T', 'N', n, m, m, 1.d0 &
, Vl, size(Vl, 1), W, size(W, 1) &
, 0.d0, Stmp, size(Stmp, 1) )
call dgemm( 'N', 'N', n, n, m, 1.d0 &
, Stmp, size(Stmp, 1), Vr, size(Vr, 1) &
, 0.d0, S, size(S, 1) )
deallocate(Stmp)
accu_nd = 0.d0
accu_d = 0.d0
do i = 1, n
do j = 1, n
if(i==j) then
accu_d += S(j,i)
else
accu_nd = accu_nd + S(j,i) * S(j,i)
endif
enddo
enddo
accu_nd = dsqrt(accu_nd)
if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n))/dble(n) .lt. thr_d)) then
print *, ' bi-orthogonal vectors without QR !'
deallocate(S)
return
endif
! -------------------------------------------------------------------------------------
! QR factorization of S: S = Q x R
print *, ' apply QR decomposition ...'
allocate( TAU(n), WORK(1) )
LWORK = -1
call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO)
if(INFO .ne. 0) then
print*,'dgeqrf failed !!', INFO
stop
endif
LWORK = max(n, int(WORK(1)))
deallocate(WORK)
allocate( WORK(LWORK) )
call dgeqrf(n, n, S, n, TAU, WORK, LWORK, INFO)
if(INFO .ne. 0) then
print*,'dgeqrf failed !!', INFO
stop
endif
! save the upper triangular R
allocate( R(n,n) )
R(:,:) = S(:,:)
! get Q
LWORK = -1
call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO)
if(INFO .ne. 0) then
print*,'dorgqr failed !!', INFO
stop
endif
LWORK = max(n, int(WORK(1)))
deallocate(WORK)
allocate( WORK(LWORK) )
call dorgqr(n, n, n, S, n, TAU, WORK, LWORK, INFO)
if(INFO .ne. 0) then
print*,'dorgqr failed !!', INFO
stop
endif
deallocate( WORK, TAU )
!
! -------------------------------------------------------------------------------------
! ---
! -------------------------------------------------------------------------------------
! get bi-orhtog left & right vectors:
! Vr' = Vr x inv(R)
! Vl' = inv(Q) x Vl = Q.T x Vl
! Q.T x Vl, where Q = S
allocate( tmp(n,m) )
call dgemm( 'T', 'T', n, m, n, 1.d0 &
, S, size(S, 1), Vl, size(Vl, 1) &
, 0.d0, tmp, size(tmp, 1) )
do i = 1, n
do j = 1, m
Vl(j,i) = tmp(i,j)
enddo
enddo
deallocate(tmp)
! ---
! inv(R)
!print *, ' inversing upper triangular matrix ...'
call dtrtri("U", "N", n, R, n, INFO)
if(INFO .ne. 0) then
print*,'dtrtri failed !!', INFO
stop
endif
!print *, ' inversing upper triangular matrix OK'
do i = 1, n-1
do j = i+1, n
R(j,i) = 0.d0
enddo
enddo
!print *, ' inv(R):'
!do i = 1, n
! write(*, '(1000(F16.10,X))') R(i,:)
!enddo
! Vr x inv(R)
allocate( tmp(m,n) )
call dgemm( 'N', 'N', m, n, n, 1.d0 &
, Vr, size(Vr, 1), R, size(R, 1) &
, 0.d0, tmp, size(tmp, 1) )
deallocate( R )
do i = 1, n
do j = 1, m
Vr(j,i) = tmp(j,i)
enddo
enddo
deallocate(tmp)
call check_weighted_biorthog_binormalize(m, n, Vl, W, Vr, thr_d, thr_nd, .false.)
return
end subroutine impose_weighted_biorthog_qr
! ---
subroutine check_weighted_biorthog_binormalize(n, m, Vl, W, Vr, thr_d, thr_nd, stop_ifnot)
implicit none
integer, intent(in) :: n, m
logical, intent(in) :: stop_ifnot
double precision, intent(in) :: thr_d, thr_nd
double precision, intent(inout) :: Vl(n,m), W(n,n), Vr(n,m)
integer :: i, j
double precision :: accu_d, accu_nd, s_tmp
double precision, allocatable :: S(:,:), Stmp(:,:)
print *, ' check weighted bi-orthonormality'
! ---
allocate(Stmp(m,n), S(m,m))
call dgemm( 'T', 'N', m, n, n, 1.d0 &
, Vl, size(Vl, 1), W, size(W, 1) &
, 0.d0, Stmp, size(Stmp, 1) )
call dgemm( 'N', 'N', m, m, n, 1.d0 &
, Stmp, size(Stmp, 1), Vr, size(Vr, 1) &
, 0.d0, S, size(S, 1) )
deallocate(Stmp)
!print *, ' overlap matrix before:'
!do i = 1, m
! write(*,'(1000(F16.10,X))') S(i,:)
!enddo
! S(i,i) = -1
do i = 1, m
if( (S(i,i) + 1.d0) .lt. thr_d ) then
do j = 1, n
Vl(j,i) = -1.d0 * Vl(j,i)
enddo
S(i,i) = 1.d0
endif
enddo
accu_d = 0.d0
accu_nd = 0.d0
do i = 1, m
do j = 1, m
if(i==j) then
accu_d = accu_d + S(i,i)
else
accu_nd = accu_nd + S(j,i) * S(j,i)
endif
enddo
enddo
accu_nd = dsqrt(accu_nd) / dble(m)
print*, ' diag acc: ', accu_d
print*, ' nondiag acc: ', accu_nd
! ---
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then
do i = 1, m
print *, i, S(i,i)
if(dabs(S(i,i) - 1.d0) .gt. thr_d) then
s_tmp = 1.d0 / dsqrt(S(i,i))
do j = 1, n
Vl(j,i) = Vl(j,i) * s_tmp
Vr(j,i) = Vr(j,i) * s_tmp
enddo
endif
enddo
endif
! ---
allocate(Stmp(m,n))
call dgemm( 'T', 'N', m, n, n, 1.d0 &
, Vl, size(Vl, 1), W, size(W, 1) &
, 0.d0, Stmp, size(Stmp, 1) )
call dgemm( 'N', 'N', m, m, n, 1.d0 &
, Stmp, size(Stmp, 1), Vr, size(Vr, 1) &
, 0.d0, S, size(S, 1) )
deallocate(Stmp)
!print *, ' overlap matrix after:'
!do i = 1, m
! write(*,'(1000(F16.10,X))') S(i,:)
!enddo
accu_d = 0.d0
accu_nd = 0.d0
do i = 1, m
do j = 1, m
if(i==j) then
accu_d = accu_d + S(i,i)
else
accu_nd = accu_nd + S(j,i) * S(j,i)
endif
enddo
enddo
accu_nd = dsqrt(accu_nd) / dble(m)
print *, ' diag acc: ', accu_d
print *, ' nondiag acc: ', accu_nd
deallocate(S)
! ---
if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) ) then
print *, accu_nd, thr_nd
print *, dabs(accu_d-dble(m))/dble(m), thr_d
print *, ' weighted biorthog_binormalize failed !'
stop
endif
end subroutine check_weighted_biorthog_binormalize
! ---
subroutine impose_weighted_biorthog_svd(n, m, overlap, L, R)
implicit none
@ -2527,6 +2799,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
call dgemm( 'T', 'N', m, m, n, 1.d0 &
, L, size(L, 1), Stmp, size(Stmp, 1) &
, 0.d0, S, size(S, 1) )
deallocate(Stmp)
print *, ' overlap bef SVD: '
do i = 1, m
@ -2598,10 +2871,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
! ---
allocate(S(m,m))
! call dgemm( 'T', 'N', m, m, n, 1.d0 &
! , L, size(L, 1), R, size(R, 1) &
! , 0.d0, S, size(S, 1) )
allocate(S(m,m),Stmp(n,m))
! S = C.T x overlap x C
call dgemm( 'N', 'N', n, m, n, 1.d0 &
, overlap, size(overlap, 1), R, size(R, 1) &
@ -2609,6 +2879,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
call dgemm( 'T', 'N', m, m, n, 1.d0 &
, L, size(L, 1), Stmp, size(Stmp, 1) &
, 0.d0, S, size(S, 1) )
deallocate(Stmp)
print *, ' overlap aft SVD with overlap: '
do i = 1, m
@ -2618,7 +2889,7 @@ subroutine impose_biorthog_svd_overlap(n, m, overlap, L, R)
! ---
end subroutine impose_biorthog_svd
end subroutine impose_weighted_biorthog_svd
! ---

185
src/tc_bi_ortho/tc_natorb.irp.f
View File

@ -1,125 +1,150 @@
BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, natorb_tc_eigval, (mo_num)]
! ---
BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, natorb_tc_eigval, (mo_num)]
BEGIN_DOC
!
! natorb_tc_reigvec_mo : RIGHT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
! natorb_tc_leigvec_mo : LEFT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
! natorb_tc_eigval : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !!
!
END_DOC
implicit none
BEGIN_DOC
! natorb_tc_reigvec_mo : RIGHT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
! natorb_tc_leigvec_mo : LEFT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
! natorb_tc_eigval : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !!
END_DOC
double precision, allocatable :: dm_tmp(:,:),fock_diag(:)
double precision :: thr_deg
integer :: i,j,k,n_real
allocate( dm_tmp(mo_num,mo_num),fock_diag(mo_num))
integer :: i, j, k, n_real
double precision :: thr_d, thr_nd, thr_deg, accu
double precision :: accu_d, accu_nd
double precision, allocatable :: dm_tmp(:,:), fock_diag(:)
allocate(dm_tmp(mo_num,mo_num), fock_diag(mo_num))
dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1)
print*,'dm_tmp'
print *, ' dm_tmp'
do i = 1, mo_num
fock_diag(i) = fock_matrix_tc_mo_tot(i,i)
write(*,'(100(F16.10,X))')-dm_tmp(:,i)
fock_diag(i) = fock_matrix_tc_mo_tot(i,i)
write(*, '(100(F16.10,X))') -dm_tmp(:,i)
enddo
thr_d = 1.d-6
thr_nd = 1.d-6
thr_deg = 1.d-3
call diag_mat_per_fock_degen(fock_diag,dm_tmp,mo_num,thr_deg,&
natorb_tc_leigvec_mo,natorb_tc_reigvec_mo,&
natorb_tc_eigval)
call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg &
, natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
! call non_hrmt_bieig( mo_num, dm_tmp&
! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo&
! , n_real, natorb_tc_eigval )
double precision :: accu
accu = 0.d0
do i = 1, n_real
print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)
accu += -natorb_tc_eigval(i)
print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)
accu += -natorb_tc_eigval(i)
enddo
print*,'accu = ',accu
print *, ' accu = ', accu
dm_tmp = 0.d0
do i = 1, n_real
accu = 0.d0
do k = 1, mo_num
accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i)
enddo
accu = 1.d0/dsqrt(dabs(accu))
natorb_tc_reigvec_mo(:,i) *= accu
natorb_tc_leigvec_mo(:,i) *= accu
do j = 1, n_real
accu = 0.d0
do k = 1, mo_num
dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j)
accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i)
enddo
accu = 1.d0/dsqrt(dabs(accu))
natorb_tc_reigvec_mo(:,i) *= accu
natorb_tc_leigvec_mo(:,i) *= accu
do j = 1, n_real
do k = 1, mo_num
dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j)
enddo
enddo
enddo
enddo
double precision :: accu_d, accu_nd
accu_d = 0.d0
accu_d = 0.d0
accu_nd = 0.d0
do i = 1, mo_num
accu_d += dm_tmp(i,i)
! write(*,'(100(F16.10,X))')dm_tmp(:,i)
do j = 1, mo_num
if(i==j)cycle
accu_nd += dabs(dm_tmp(j,i))
enddo
accu_d += dm_tmp(i,i)
!write(*,'(100(F16.10,X))')dm_tmp(:,i)
do j = 1, mo_num
if(i==j)cycle
accu_nd += dabs(dm_tmp(j,i))
enddo
enddo
print*,'Trace of the overlap between TC natural orbitals ',accu_d
print*,'L1 norm of extra diagonal elements of overlap matrix ',accu_nd
print *, ' Trace of the overlap between TC natural orbitals ', accu_d
print *, ' L1 norm of extra diagonal elements of overlap matrix ', accu_nd
deallocate(dm_tmp, fock_diag)
END_PROVIDER
! ---
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_diag_sorted_r_natorb, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, fock_diag_sorted_l_natorb, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, fock_diag_sorted_v_natorb, (mo_num)]
BEGIN_PROVIDER [ double precision, fock_diag_sorted_r_natorb, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, fock_diag_sorted_l_natorb, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, fock_diag_sorted_v_natorb, (mo_num)]
implicit none
integer ::i,j,k
print*,'Diagonal elements of the Fock matrix before '
do i = 1, mo_num
write(*,*)i,Fock_matrix_tc_mo_tot(i,i)
enddo
integer :: i,j,k
integer, allocatable :: iorder(:)
double precision, allocatable :: fock_diag(:)
print *, ' Diagonal elements of the Fock matrix before '
do i = 1, mo_num
write(*,*) i, Fock_matrix_tc_mo_tot(i,i)
enddo
allocate(fock_diag(mo_num))
fock_diag = 0.d0
do i = 1, mo_num
fock_diag(i) = 0.d0
do j = 1, mo_num
do k = 1, mo_num
fock_diag(i) += natorb_tc_leigvec_mo(k,i) * Fock_matrix_tc_mo_tot(k,j) * natorb_tc_reigvec_mo(j,i)
fock_diag(i) = 0.d0
do j = 1, mo_num
do k = 1, mo_num
fock_diag(i) += natorb_tc_leigvec_mo(k,i) * Fock_matrix_tc_mo_tot(k,j) * natorb_tc_reigvec_mo(j,i)
enddo
enddo
enddo
enddo
integer, allocatable :: iorder(:)
allocate(iorder(mo_num))
do i = 1, mo_num
iorder(i) = i
enddo
call dsort(fock_diag,iorder,mo_num)
print*,'Diagonal elements of the Fock matrix after '
call dsort(fock_diag, iorder, mo_num)
print *, ' Diagonal elements of the Fock matrix after '
do i = 1, mo_num
write(*,*)i,fock_diag(i)
write(*,*) i, fock_diag(i)
enddo
deallocate(fock_diag)
do i = 1, mo_num
fock_diag_sorted_v_natorb(i) = natorb_tc_eigval(iorder(i))
do j = 1, mo_num
fock_diag_sorted_r_natorb(j,i) = natorb_tc_reigvec_mo(j,iorder(i))
fock_diag_sorted_l_natorb(j,i) = natorb_tc_leigvec_mo(j,iorder(i))
enddo
fock_diag_sorted_v_natorb(i) = natorb_tc_eigval(iorder(i))
do j = 1, mo_num
fock_diag_sorted_r_natorb(j,i) = natorb_tc_reigvec_mo(j,iorder(i))
fock_diag_sorted_l_natorb(j,i) = natorb_tc_leigvec_mo(j,iorder(i))
enddo
enddo
deallocate(iorder)
END_PROVIDER
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_ao, (ao_num, mo_num)]
&BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_ao, (ao_num, mo_num)]
&BEGIN_PROVIDER [ double precision, overlap_natorb_tc_eigvec_ao, (mo_num, mo_num) ]
BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_ao, (ao_num, mo_num)]
&BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_ao, (ao_num, mo_num)]
&BEGIN_PROVIDER [ double precision, overlap_natorb_tc_eigvec_ao, (mo_num, mo_num) ]
BEGIN_DOC
! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP
!
! THE OVERLAP SHOULD BE THE SAME AS overlap_natorb_tc_eigvec_mo
END_DOC
BEGIN_DOC
! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP
!
! THE OVERLAP SHOULD BE THE SAME AS overlap_natorb_tc_eigvec_mo
END_DOC
implicit none
integer :: i, j, k, q, p
double precision :: accu, accu_d
double precision, allocatable :: tmp(:,:)
implicit none
integer :: i, j, k, q, p
double precision :: accu, accu_d
double precision, allocatable :: tmp(:,:)
! ! MO_R x R

30
src/tc_scf/diago_bi_ort_tcfock.irp.f
View File

@ -11,20 +11,24 @@
integer :: n_real_tc
integer :: i, k, l
double precision :: accu_d, accu_nd, accu_tmp
double precision :: thr_d, thr_nd
double precision :: norm
double precision, allocatable :: eigval_right_tmp(:)
thr_d = 1d-6
thr_nd = 1d-6
allocate( eigval_right_tmp(mo_num) )
PROVIDE Fock_matrix_tc_mo_tot
call non_hrmt_bieig( mo_num, Fock_matrix_tc_mo_tot &
, fock_tc_leigvec_mo, fock_tc_reigvec_mo &
, n_real_tc, eigval_right_tmp )
call non_hrmt_bieig( mo_num, Fock_matrix_tc_mo_tot, thr_d, thr_nd &
, fock_tc_leigvec_mo, fock_tc_reigvec_mo &
, n_real_tc, eigval_right_tmp )
!if(max_ov_tc_scf)then
! call non_hrmt_fock_mat( mo_num, Fock_matrix_tc_mo_tot &
! , fock_tc_leigvec_mo, fock_tc_reigvec_mo &
! , n_real_tc, eigval_right_tmp )
! call non_hrmt_fock_mat( mo_num, Fock_matrix_tc_mo_tot, thr_d, thr_nd &
! , fock_tc_leigvec_mo, fock_tc_reigvec_mo &
! , n_real_tc, eigval_right_tmp )
!else
! call non_hrmt_diag_split_degen_bi_orthog( mo_num, Fock_matrix_tc_mo_tot &
! , fock_tc_leigvec_mo, fock_tc_reigvec_mo &
@ -59,17 +63,17 @@
else
accu_tmp = overlap_fock_tc_eigvec_mo(k,i)
accu_nd += accu_tmp * accu_tmp
if(dabs(overlap_fock_tc_eigvec_mo(k,i)).gt.1.d-10)then
print*,'k,i',k,i,overlap_fock_tc_eigvec_mo(k,i)
if(dabs(overlap_fock_tc_eigvec_mo(k,i)) .gt. thr_nd)then
print *, 'k,i', k, i, overlap_fock_tc_eigvec_mo(k,i)
endif
endif
enddo
enddo
accu_nd = dsqrt(accu_nd)/accu_d
if( accu_nd .gt. 1d-8 ) then
if(accu_nd .gt. thr_nd) then
print *, ' bi-orthog failed'
print*,'accu_nd MO = ', accu_nd
print*,'accu_nd MO = ', accu_nd, thr_nd
print*,'overlap_fock_tc_eigvec_mo = '
do i = 1, mo_num
write(*,'(100(F16.10,X))') overlap_fock_tc_eigvec_mo(i,:)
@ -77,13 +81,13 @@
stop
endif
if( dabs(accu_d - dble(mo_num)) .gt. 1e-7 ) then
if( dabs(accu_d - dble(mo_num))/dble(mo_num) .gt. thr_d ) then
print *, 'mo_num = ', mo_num
print *, 'accu_d MO = ', accu_d
print *, 'accu_d MO = ', accu_d, thr_d
print *, 'normalizing vectors ...'
do i = 1, mo_num
norm = dsqrt(dabs(overlap_fock_tc_eigvec_mo(i,i)))
if( norm.gt.1e-7 ) then
if(norm .gt. thr_d) then
do k = 1, mo_num
fock_tc_reigvec_mo(k,i) *= 1.d0/norm
fock_tc_leigvec_mo(k,i) *= 1.d0/norm

10
src/tc_scf/fock_tc.irp.f
View File

@ -116,17 +116,23 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
endif
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis
! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis
END_DOC
implicit none
Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta)
if(three_body_h_tc) then
Fock_matrix_tc_mo_tot += fock_3_mat
endif
!call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10)
END_PROVIDER
! ---

131
src/tc_scf/rotate_tcscf_orbitals.irp.f
View File

@ -73,7 +73,7 @@ subroutine maximize_overlap()
! ---
call rotate_degen_eigvec(n, m, e, C, W, L, R)
call rotate_degen_eigvec_to_maximize_overlap(n, m, e, C, W, L, R)
! ---
@ -115,7 +115,7 @@ end subroutine maximize_overlap
! ---
subroutine rotate_degen_eigvec(n, m, e0, C0, W0, L0, R0)
subroutine rotate_degen_eigvec_to_maximize_overlap(n, m, e0, C0, W0, L0, R0)
implicit none
@ -124,7 +124,7 @@ subroutine rotate_degen_eigvec(n, m, e0, C0, W0, L0, R0)
double precision, intent(inout) :: L0(n,m), R0(n,m)
integer :: i, j, k, kk, mm, id1
integer :: i, j, k, kk, mm, id1, tot_deg
double precision :: ei, ej, de, de_thr
integer, allocatable :: deg_num(:)
double precision, allocatable :: L(:,:), R(:,:), C(:,:), Lnew(:,:), Rnew(:,:), tmp(:,:)
@ -162,12 +162,19 @@ subroutine rotate_degen_eigvec(n, m, e0, C0, W0, L0, R0)
enddo
enddo
tot_deg = 0
do i = 1, m
if(deg_num(i).gt.1) then
print *, ' degen on', i, deg_num(i)
tot_deg = tot_deg + 1
endif
enddo
if(tot_deg .eq. 0) then
print *, ' no degen'
return
endif
! ---
do i = 1, m
@ -243,6 +250,122 @@ subroutine rotate_degen_eigvec(n, m, e0, C0, W0, L0, R0)
deallocate(S, Snew, T)
end subroutine rotate_degen_eigvec
end subroutine rotate_degen_eigvec_to_maximize_overlap
! ---
subroutine fix_right_to_one()
implicit none
integer :: i, j, m, n, mm, tot_deg
double precision :: accu_d, accu_nd
double precision :: de_thr, ei, ej, de
double precision :: thr_d, thr_nd
integer, allocatable :: deg_num(:)
double precision, allocatable :: R0(:,:), L0(:,:), W(:,:), e0(:)
double precision, allocatable :: R(:,:), L(:,:), S(:,:), Stmp(:,:), tmp(:,:)
thr_d = 1d-7
thr_nd = 1d-7
n = ao_num
m = mo_num
allocate(L0(n,m), R0(n,m), W(n,n), e0(m))
L0 = mo_l_coef
R0 = mo_r_coef
W = ao_overlap
print*, ' fock matrix diag elements'
do i = 1, m
e0(i) = Fock_matrix_tc_mo_tot(i,i)
print*, e0(i)
enddo
! ---
allocate( deg_num(m) )
do i = 1, m
deg_num(i) = 1
enddo
de_thr = 1d-6
do i = 1, m-1
ei = e0(i)
! already considered in degen vectors
if(deg_num(i).eq.0) cycle
do j = i+1, m
ej = e0(j)
de = dabs(ei - ej)
if(de .lt. de_thr) then
deg_num(i) = deg_num(i) + 1
deg_num(j) = 0
endif
enddo
enddo
deallocate(e0)
tot_deg = 0
do i = 1, m
if(deg_num(i).gt.1) then
print *, ' degen on', i, deg_num(i)
tot_deg = tot_deg + 1
endif
enddo
if(tot_deg .eq. 0) then
print *, ' no degen'
return
endif
! ---
do i = 1, m
mm = deg_num(i)
if(mm .gt. 1) then
allocate(L(n,mm), R(n,mm))
do j = 1, mm
L(1:n,j) = L0(1:n,i+j-1)
R(1:n,j) = R0(1:n,i+j-1)
enddo
! ---
call impose_weighted_orthog_svd(n, mm, W, R)
call impose_weighted_biorthog_qr(n, mm, thr_d, thr_nd, R, W, L)
! ---
do j = 1, mm
L0(1:n,i+j-1) = L(1:n,j)
R0(1:n,i+j-1) = R(1:n,j)
enddo
deallocate(L, R)
endif
enddo
call check_weighted_biorthog_binormalize(n, m, L0, W, R0, thr_d, thr_nd, .true.)
deallocate(W, deg_num)
mo_l_coef = L0
mo_r_coef = R0
deallocate(L0, R0)
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
print *, ' orbitals are rotated '
return
end subroutine fix_right_to_one
! ---

444
src/tc_scf/routines_rotates.irp.f
View File

@ -1,203 +1,242 @@
subroutine minimize_tc_orb_angles
implicit none
double precision :: thr_deg
logical :: good_angles
integer :: i
good_angles = .False.
thr_deg = thr_degen_tc
call print_energy_and_mos
i = 1
print*,'Minimizing the angles between the TC orbitals'
do while (.not. good_angles)
print*,'iteration = ',i
call routine_save_rotated_mos(thr_deg,good_angles)
thr_deg *= 10.d0
i+=1
if(i.gt.100)then
print*,'minimize_tc_orb_angles does not seem to converge ..'
print*,'Something is weird in the tc orbitals ...'
print*,'STOPPING'
endif
enddo
print*,'Converged ANGLES MINIMIZATION !!'
call print_angles_tc
call print_energy_and_mos
! ---
subroutine minimize_tc_orb_angles()
implicit none
logical :: good_angles
integer :: i
double precision :: thr_deg
good_angles = .False.
thr_deg = thr_degen_tc
call print_energy_and_mos()
print *, ' Minimizing the angles between the TC orbitals'
i = 1
do while (.not. good_angles)
print *, ' iteration = ', i
call routine_save_rotated_mos(thr_deg, good_angles)
thr_deg *= 10.d0
i += 1
if(i .gt. 100) then
print *, ' minimize_tc_orb_angles does not seem to converge ..'
print *, ' Something is weird in the tc orbitals ...'
print *, ' STOPPING'
stop
endif
enddo
print *, ' Converged ANGLES MINIMIZATION !!'
call print_angles_tc()
call print_energy_and_mos()
end
subroutine routine_save_rotated_mos(thr_deg,good_angles)
implicit none
double precision, intent(in) :: thr_deg
logical, intent(out) :: good_angles
good_angles = .False.
integer :: i,j,k,n_degen_list,m,n,n_degen,ilast,ifirst
double precision, allocatable :: mo_r_coef_good(:,:),mo_l_coef_good(:,:)
allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
double precision, allocatable :: mo_r_coef_new(:,:)
double precision :: norm
print*,'***************************************'
print*,'***************************************'
print*,'THRESHOLD FOR DEGENERACIES ::: ',thr_deg
print*,'***************************************'
print*,'***************************************'
print*,'Starting with the following TC energy gradient :',grad_non_hermit
mo_r_coef_good = mo_r_coef
mo_l_coef_good = mo_l_coef
allocate(mo_r_coef_new(ao_num, mo_num))
mo_r_coef_new = mo_r_coef
do i = 1, mo_num
norm = 1.d0/dsqrt(overlap_mo_r(i,i))
do j = 1, ao_num
mo_r_coef_new(j,i) *= norm
! ---
subroutine routine_save_rotated_mos(thr_deg, good_angles)
implicit none
double precision, intent(in) :: thr_deg
logical, intent(out) :: good_angles
integer :: i, j, k, n_degen_list, m, n, n_degen, ilast, ifirst
double precision :: max_angle, norm
integer, allocatable :: list_degen(:,:)
double precision, allocatable :: new_angles(:)
double precision, allocatable :: mo_r_coef_good(:,:), mo_l_coef_good(:,:)
double precision, allocatable :: mo_r_coef_new(:,:)
double precision, allocatable :: fock_diag(:),s_mat(:,:)
double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:)
double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:)
good_angles = .False.
allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
print *, ' ***************************************'
print *, ' ***************************************'
print *, ' THRESHOLD FOR DEGENERACIES ::: ', thr_deg
print *, ' ***************************************'
print *, ' ***************************************'
print *, ' Starting with the following TC energy gradient :', grad_non_hermit
mo_r_coef_good = mo_r_coef
mo_l_coef_good = mo_l_coef
allocate(mo_r_coef_new(ao_num, mo_num))
mo_r_coef_new = mo_r_coef
do i = 1, mo_num
norm = 1.d0/dsqrt(overlap_mo_r(i,i))
do j = 1, ao_num
mo_r_coef_new(j,i) *= norm
enddo
enddo
enddo
double precision, allocatable :: fock_diag(:),s_mat(:,:)
integer, allocatable :: list_degen(:,:)
allocate(list_degen(mo_num,0:mo_num),s_mat(mo_num,mo_num),fock_diag(mo_num))
do i = 1, mo_num
fock_diag(i) = Fock_matrix_tc_mo_tot(i,i)
enddo
allocate(list_degen(mo_num,0:mo_num), s_mat(mo_num,mo_num), fock_diag(mo_num))
do i = 1, mo_num
fock_diag(i) = Fock_matrix_tc_mo_tot(i,i)
enddo
! compute the overlap between the left and rescaled right
call build_s_matrix(ao_num,mo_num,mo_r_coef_new,mo_r_coef_new,ao_overlap,s_mat)
call build_s_matrix(ao_num, mo_num, mo_r_coef_new, mo_r_coef_new, ao_overlap, s_mat)
! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
call give_degen_full_list(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
print*,'fock_matrix_mo'
do i = 1, mo_num
print*,i,fock_diag(i),angle_left_right(i)
enddo
call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list)
print *, ' fock_matrix_mo'
do i = 1, mo_num
print *, i, fock_diag(i), angle_left_right(i)
enddo
do i = 1, n_degen_list
do i = 1, n_degen_list
! ifirst = list_degen(1,i)
! ilast = list_degen(2,i)
! n_degen = ilast - ifirst +1
n_degen = list_degen(i,0)
double precision, allocatable :: stmp(:,:),T(:,:),Snew(:,:),smat2(:,:)
double precision, allocatable :: mo_l_coef_tmp(:,:),mo_r_coef_tmp(:,:),mo_l_coef_new(:,:)
allocate(stmp(n_degen,n_degen),smat2(n_degen,n_degen))
allocate(mo_r_coef_tmp(ao_num,n_degen),mo_l_coef_tmp(ao_num,n_degen),mo_l_coef_new(ao_num,n_degen))
allocate(T(n_degen,n_degen),Snew(n_degen,n_degen))
do j = 1, n_degen
mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j))
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
enddo
! Orthogonalization of right functions
print*,'Orthogonalization of RIGHT functions'
print*,'------------------------------------'
call orthog_functions(ao_num,n_degen,mo_r_coef_tmp,ao_overlap)
! Orthogonalization of left functions
print*,'Orthogonalization of LEFT functions'
print*,'------------------------------------'
call orthog_functions(ao_num,n_degen,mo_l_coef_tmp,ao_overlap)
print*,'Overlap lef-right '
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_l_coef_tmp,ao_overlap,stmp)
do j = 1, n_degen
write(*,'(100(F8.4,X))')stmp(:,j)
n_degen = list_degen(i,0)
allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen))
allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen))
allocate(T(n_degen,n_degen), Snew(n_degen,n_degen))
do j = 1, n_degen
mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j))
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
enddo
! Orthogonalization of right functions
print *, ' Orthogonalization of RIGHT functions'
print *, ' ------------------------------------'
call orthog_functions(ao_num, n_degen, mo_r_coef_tmp, ao_overlap)
! Orthogonalization of left functions
print *, ' Orthogonalization of LEFT functions'
print *, ' ------------------------------------'
call orthog_functions(ao_num, n_degen, mo_l_coef_tmp, ao_overlap)
print *, ' Overlap lef-right '
call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_l_coef_tmp, ao_overlap, stmp)
do j = 1, n_degen
write(*,'(100(F8.4,X))') stmp(:,j)
enddo
call build_s_matrix(ao_num, n_degen, mo_l_coef_tmp, mo_l_coef_tmp, ao_overlap, stmp)
!print*,'LEFT/LEFT OVERLAP '
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_r_coef_tmp, ao_overlap, stmp)
!print*,'RIGHT/RIGHT OVERLAP '
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
if(maxovl_tc) then
T = 0.d0
Snew = 0.d0
call maxovl(n_degen, n_degen, stmp, T, Snew)
!print*,'overlap after'
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')Snew(:,j)
!enddo
call dgemm( 'N', 'N', ao_num, n_degen, n_degen, 1.d0 &
, mo_l_coef_tmp, size(mo_l_coef_tmp, 1), T(1,1), size(T, 1) &
, 0.d0, mo_l_coef_new, size(mo_l_coef_new, 1) )
call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_r_coef_tmp, ao_overlap, stmp)
!print*,'Overlap test'
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
else
mo_l_coef_new = mo_l_coef_tmp
endif
call impose_weighted_biorthog_svd(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp)
!call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_r_coef_tmp, ao_overlap, stmp)
!print*,'LAST OVERLAP '
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
!call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_l_coef_new, ao_overlap, stmp)
!print*,'LEFT OVERLAP '
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
!call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_r_coef_tmp, ao_overlap, stmp)
!print*,'RIGHT OVERLAP '
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
do j = 1, n_degen
!!! mo_l_coef_good(1:ao_num,j+ifirst-1) = mo_l_coef_new(1:ao_num,j)
!!! mo_r_coef_good(1:ao_num,j+ifirst-1) = mo_r_coef_tmp(1:ao_num,j)
mo_l_coef_good(1:ao_num,list_degen(i,j)) = mo_l_coef_new(1:ao_num,j)
mo_r_coef_good(1:ao_num,list_degen(i,j)) = mo_r_coef_tmp(1:ao_num,j)
enddo
deallocate(stmp, smat2)
deallocate(mo_r_coef_tmp, mo_l_coef_tmp, mo_l_coef_new)
deallocate(T, Snew)
enddo
call build_s_matrix(ao_num,n_degen,mo_l_coef_tmp,mo_l_coef_tmp,ao_overlap,stmp)
!allocate(stmp(mo_num, mo_num))
!call build_s_matrix(ao_num, mo_num, mo_l_coef_good, mo_r_coef_good, ao_overlap, stmp)
!print*,'LEFT/RIGHT OVERLAP '
!do j = 1, mo_num
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
!call build_s_matrix(ao_num, mo_num, mo_l_coef_good, mo_l_coef_good, ao_overlap, stmp)
!print*,'LEFT/LEFT OVERLAP '
!do j = 1, n_degen
!do j = 1, mo_num
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_r_coef_tmp,ao_overlap,stmp)
!call build_s_matrix(ao_num, mo_num, mo_r_coef_good, mo_r_coef_good, ao_overlap, stmp)
!print*,'RIGHT/RIGHT OVERLAP '
!do j = 1, n_degen
!do j = 1, mo_num
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
if(maxovl_tc)then
T = 0.d0
Snew = 0.d0
call maxovl(n_degen, n_degen, stmp, T, Snew)
!print*,'overlap after'
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')Snew(:,j)
!enddo
call dgemm( 'N', 'N', ao_num, n_degen, n_degen, 1.d0 &
, mo_l_coef_tmp, size(mo_l_coef_tmp, 1), T(1,1), size(T, 1) &
, 0.d0, mo_l_coef_new, size(mo_l_coef_new, 1) )
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
!print*,'Overlap test'
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
else
mo_l_coef_new = mo_l_coef_tmp
endif
call impose_biorthog_svd_overlap(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp)
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_r_coef_tmp,ao_overlap,stmp)
!print*,'LAST OVERLAP '
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
call build_s_matrix(ao_num,n_degen,mo_l_coef_new,mo_l_coef_new,ao_overlap,stmp)
!print*,'LEFT OVERLAP '
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
call build_s_matrix(ao_num,n_degen,mo_r_coef_tmp,mo_r_coef_tmp,ao_overlap,stmp)
!print*,'RIGHT OVERLAP '
!do j = 1, n_degen
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
do j = 1, n_degen
! mo_l_coef_good(1:ao_num,j+ifirst-1) = mo_l_coef_new(1:ao_num,j)
! mo_r_coef_good(1:ao_num,j+ifirst-1) = mo_r_coef_tmp(1:ao_num,j)
mo_l_coef_good(1:ao_num,list_degen(i,j)) = mo_l_coef_new(1:ao_num,j)
mo_r_coef_good(1:ao_num,list_degen(i,j)) = mo_r_coef_tmp(1:ao_num,j)
enddo
deallocate(stmp,smat2)
deallocate(mo_r_coef_tmp,mo_l_coef_tmp,mo_l_coef_new)
deallocate(T,Snew)
enddo
allocate(stmp(mo_num, mo_num))
call build_s_matrix(ao_num,mo_num,mo_l_coef_good,mo_r_coef_good,ao_overlap,stmp)
!print*,'LEFT/RIGHT OVERLAP '
!do j = 1, mo_num
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
call build_s_matrix(ao_num,mo_num,mo_l_coef_good,mo_l_coef_good,ao_overlap,stmp)
!print*,'LEFT/LEFT OVERLAP '
!do j = 1, mo_num
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
call build_s_matrix(ao_num,mo_num,mo_r_coef_good,mo_r_coef_good,ao_overlap,stmp)
!print*,'RIGHT/RIGHT OVERLAP '
!do j = 1, mo_num
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
mo_r_coef = mo_r_coef_good
mo_l_coef = mo_l_coef_good
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
TOUCH mo_l_coef mo_r_coef
double precision, allocatable :: new_angles(:)
allocate(new_angles(mo_num))
new_angles(1:mo_num) = dabs(angle_left_right(1:mo_num))
double precision :: max_angle
max_angle = maxval(new_angles)
good_angles = max_angle.lt.45.d0
print*,'max_angle = ',max_angle
print *, ' max_angle = ', max_angle
end
subroutine build_s_matrix(m,n,C1,C2,overlap,smat)
implicit none
integer, intent(in) :: m,n
double precision, intent(in) :: C1(m,n),C2(m,n),overlap(m,m)
double precision, intent(out):: smat(n,n)
integer :: i,j,k,l
smat = 0.D0
! ---
subroutine build_s_matrix(m, n, C1, C2, overlap, smat)
implicit none
integer, intent(in) :: m, n
double precision, intent(in) :: C1(m,n), C2(m,n), overlap(m,m)
double precision, intent(out) :: smat(n,n)
integer :: i, j, k, l
smat = 0.D0
do i = 1, n
do j = 1, n
do k = 1, m
do l = 1, m
smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j)
enddo
do j = 1, n
do k = 1, m
do l = 1, m
smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j)
enddo
enddo
enddo
enddo
enddo
end
! ---
subroutine orthog_functions(m,n,coef,overlap)
implicit none
integer, intent(in) :: m,n
@ -217,45 +256,63 @@ subroutine orthog_functions(m,n,coef,overlap)
coef(1,:m) *= 1.d0/dsqrt(stmp(j,j))
enddo
call build_s_matrix(m,n,coef,coef,overlap,stmp)
!print*,'overlap after'
!do j = 1, n
! write(*,'(100(F16.10,X))')stmp(:,j)
!enddo
deallocate(stmp)
end
subroutine print_angles_tc
implicit none
integer :: i,j
double precision :: left,right
print*,'product of norms, angle between vectors'
do i = 1, mo_num
left = overlap_mo_l(i,i)
right = overlap_mo_r(i,i)
! print*,Fock_matrix_tc_mo_tot(i,i),left*right,angle_left_right(i)
print*,left*right,angle_left_right(i)
enddo
end
! ---
subroutine print_energy_and_mos
implicit none
integer :: i
print*,''
print*,'TC energy = ', TC_HF_energy
print*,'TC SCF energy gradient = ',grad_non_hermit
print*,'Max angle Left/right = ',max_angle_left_right
if(max_angle_left_right.lt.45.d0)then
print*,'Maximum angle BELOW 45 degrees, everthing is OK !'
else if(max_angle_left_right.gt.45.d0.and.max_angle_left_right.lt.75.d0)then
print*,'Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...'
else if(max_angle_left_right.gt.75.d0)then
print*,'Maximum angle between ABOVE 75 degrees, YOU WILL CERTAINLY FIND TROUBLES IN TC-CI calculations ...'
endif
print*,'Diag Fock elem, product of left/right norm, angle left/right '
subroutine print_angles_tc()
implicit none
integer :: i, j
double precision :: left, right
print *, ' product of norms, angle between vectors'
do i = 1, mo_num
write(*,'(I3,X,100(F16.10,X))')i,Fock_matrix_tc_mo_tot(i,i),overlap_mo_l(i,i)*overlap_mo_r(i,i),angle_left_right(i)
left = overlap_mo_l(i,i)
right = overlap_mo_r(i,i)
! print*,Fock_matrix_tc_mo_tot(i,i),left*right,angle_left_right(i)
print *, left*right, angle_left_right(i)
enddo
end
! ---
subroutine print_energy_and_mos()
implicit none
integer :: i
print *, ' '
print *, ' TC energy = ', TC_HF_energy
print *, ' TC SCF energy gradient = ', grad_non_hermit
print *, ' Max angle Left/right = ', max_angle_left_right
if(max_angle_left_right .lt. 45.d0) then
print *, ' Maximum angle BELOW 45 degrees, everthing is OK !'
else if(max_angle_left_right .gt. 45.d0 .and. max_angle_left_right .lt. 75.d0) then
print *, ' Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...'
else if(max_angle_left_right .gt. 75.d0) then
print *, ' Maximum angle between ABOVE 75 degrees, YOU WILL CERTAINLY FIND TROUBLES IN TC-CI calculations ...'
endif
print *, ' Diag Fock elem, product of left/right norm, angle left/right '
do i = 1, mo_num
write(*, '(I3,X,100(F16.10,X))') i, Fock_matrix_tc_mo_tot(i,i), overlap_mo_l(i,i)*overlap_mo_r(i,i), angle_left_right(i)
enddo
end
! ---
subroutine sort_by_tc_fock
implicit none
integer, allocatable :: iorder(:)
@ -276,3 +333,4 @@ subroutine sort_by_tc_fock
touch mo_l_coef mo_r_coef
end

330
src/tools/print_he_energy.irp.f Normal file
View File

@ -0,0 +1,330 @@
! ---
program print_he_energy
implicit none
!call print_energy1()
call print_energy2()
end
! ---
subroutine print_energy1()
implicit none
integer :: i, j, k, l
double precision :: e, n, e_tmp, n_tmp
double precision, external :: ao_two_e_integral
e = 0.d0
n = 0.d0
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
! < phi_1 phi_1 | h1 | phi_1 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,1)
n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
enddo
enddo
e += e_tmp * n_tmp
! ---
! < phi_1 phi_1 | h2 | phi_1 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,1)
enddo
enddo
e += e_tmp * n_tmp
! ---
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
! ---
do i = 1, ao_num
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
! ( phi_1 phi_1 | phi_1 phi_1 )
e += mo_coef(i,1) * mo_coef(j,1) * ao_two_e_integral(i,j,k,l) * mo_coef(k,1) * mo_coef(l,1)
enddo
enddo
enddo
enddo
! ---
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
! ---
! < phi_1 phi_1 | phi_1 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
enddo
enddo
n += e_tmp * n_tmp
! ---
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
e = e / n
print *, ' energy = ', e
end subroutine print_energy1
! ---
subroutine print_energy2()
implicit none
integer :: i, j, k, l
double precision :: e, n, e_tmp, n_tmp
double precision, external :: ao_two_e_integral
e = 0.d0
n = 0.d0
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
! < phi_1 phi_2 | h1 | phi_1 phi_2 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,1)
n_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,2)
enddo
enddo
e += e_tmp * n_tmp
! ---
! < phi_1 phi_2 | h1 | phi_2 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,2)
n_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,1)
enddo
enddo
e -= e_tmp * n_tmp
! ---
! < phi_2 phi_1 | h1 | phi_1 phi_2 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,2) * ao_one_e_integrals(i,j) * mo_coef(j,1)
n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,2)
enddo
enddo
e -= e_tmp * n_tmp
! ---
! < phi_2 phi_1 | h1 | phi_2 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,2) * ao_one_e_integrals(i,j) * mo_coef(j,2)
n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
enddo
enddo
e += e_tmp * n_tmp
! ---
! < phi_1 phi_2 | h2 | phi_1 phi_2 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
e_tmp += mo_coef(i,2) * ao_one_e_integrals(i,j) * mo_coef(j,2)
enddo
enddo
e += e_tmp * n_tmp
! ---
! < phi_1 phi_2 | h2 | phi_2 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,2)
e_tmp += mo_coef(i,2) * ao_one_e_integrals(i,j) * mo_coef(j,1)
enddo
enddo
e -= e_tmp * n_tmp
! ---
! < phi_2 phi_1 | h2 | phi_1 phi_2 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
n_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,1)
e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,2)
enddo
enddo
e -= e_tmp * n_tmp
! ---
! < phi_2 phi_1 | h2 | phi_2 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
n_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,2)
e_tmp += mo_coef(i,1) * ao_one_e_integrals(i,j) * mo_coef(j,1)
enddo
enddo
e += e_tmp * n_tmp
! ---
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
! ---
do i = 1, ao_num
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
! ( phi_1 phi_1 | phi_2 phi_2 )
e += mo_coef(i,1) * mo_coef(j,1) * ao_two_e_integral(i,j,k,l) * mo_coef(k,2) * mo_coef(l,2)
! ( phi_1 phi_2 | phi_2 phi_1 )
e -= mo_coef(i,1) * mo_coef(j,2) * ao_two_e_integral(i,j,k,l) * mo_coef(k,2) * mo_coef(l,1)
! ( phi_2 phi_1 | phi_1 phi_2 )
e -= mo_coef(i,2) * mo_coef(j,1) * ao_two_e_integral(i,j,k,l) * mo_coef(k,1) * mo_coef(l,2)
! ( phi_2 phi_2 | phi_1 phi_1 )
e += mo_coef(i,2) * mo_coef(j,2) * ao_two_e_integral(i,j,k,l) * mo_coef(k,1) * mo_coef(l,1)
enddo
enddo
enddo
enddo
! ---
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
! ---
! < phi_1 phi_2 | phi_1 phi_2 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
n_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,2)
enddo
enddo
n += e_tmp * n_tmp
! ---
! < phi_1 phi_2 | phi_2 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,2)
n_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,1)
enddo
enddo
n -= e_tmp * n_tmp
! ---
! < phi_2 phi_1 | phi_1 phi_2 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,1)
n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,2)
enddo
enddo
n -= e_tmp * n_tmp
! ---
! < phi_2 phi_1 | phi_2 phi_1 >
e_tmp = 0.d0
n_tmp = 0.d0
do i = 1, ao_num
do j = 1, ao_num
e_tmp += mo_coef(i,2) * ao_overlap(i,j) * mo_coef(j,2)
n_tmp += mo_coef(i,1) * ao_overlap(i,j) * mo_coef(j,1)
enddo
enddo
n += e_tmp * n_tmp
! ---
! --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
e = e / n
print *, ' energy = ', e
end subroutine print_energy2
! ---

369
src/utils/block_diag_degen.irp.f
View File

@ -1,177 +1,218 @@
subroutine diag_mat_per_fock_degen(fock_diag,mat_ref,n,thr_deg,leigvec,reigvec,eigval)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: fock_diag(n),mat_ref(n,n),thr_deg
double precision, intent(out):: leigvec(n,n),reigvec(n,n),eigval(n)
BEGIN_DOC
! subroutine that diagonalizes a matrix mat_ref BY BLOCK
!
! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag"
!
! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together
!
! : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together
!
! : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together
!
! etc... the advantage is to guarentee no spurious mixing because of numerical problems.
END_DOC
double precision, allocatable :: leigvec_unsrtd(:,:),reigvec_unsrtd(:,:),eigval_unsrtd(:)
integer, allocatable :: list_degen(:,:),list_same_degen(:)
integer, allocatable :: iorder(:),list_degen_sorted(:)
integer :: n_degen_list,n_degen,size_mat,i,j,k,icount,m,index_degen
integer :: ii,jj,i_good,j_good,n_real
double precision, allocatable :: mat_tmp(:,:),eigval_tmp(:),leigvec_tmp(:,:),reigvec_tmp(:,:)
subroutine diag_mat_per_fock_degen(fock_diag, mat_ref, n, thr_d, thr_nd, thr_deg, leigvec, reigvec, eigval)
allocate(leigvec_unsrtd(n,n),reigvec_unsrtd(n,n),eigval_unsrtd(n))
leigvec_unsrtd = 0.d0
reigvec_unsrtd = 0.d0
eigval_unsrtd = 0.d0
allocate(list_degen(n,0:n))
BEGIN_DOC
!
! subroutine that diagonalizes a matrix mat_ref BY BLOCK
!
! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag"
!
! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together
!
! : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together
!
! : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together
!
! etc... the advantage is to guarentee no spurious mixing because of numerical problems.
!
END_DOC
! obtain degeneracies
call give_degen_full_list(fock_diag,n,thr_deg,list_degen,n_degen_list)
allocate(iorder(n_degen_list),list_degen_sorted(n_degen_list))
do i =1, n_degen_list
n_degen = list_degen(i,0)
list_degen_sorted(i) = n_degen
iorder(i)=i
enddo
! sort by number of degeneracies
call isort(list_degen_sorted,iorder,n_degen_list)
integer :: icount_eigval
logical, allocatable :: is_ok(:)
allocate(is_ok(n_degen_list))
is_ok = .True.
icount_eigval = 0
! loop over degeneracies
do i = 1, n_degen_list
if(.not.is_ok(i))cycle
is_ok(i) = .False.
n_degen = list_degen_sorted(i)
print*,'diagonalizing for n_degen = ',n_degen
k = 1
! group all the entries having the same degeneracies
! do while (list_degen_sorted(i+k)==n_degen)
do m = i+1,n_degen_list
if(list_degen_sorted(m)==n_degen)then
is_ok(i+k)=.False.
k += 1
endif
implicit none
integer, intent(in) :: n
double precision, intent(in) :: fock_diag(n), mat_ref(n,n), thr_d, thr_nd, thr_deg
double precision, intent(out) :: leigvec(n,n), reigvec(n,n), eigval(n)
integer :: n_degen_list, n_degen,size_mat, i, j, k, icount, m, index_degen
integer :: ii, jj, i_good, j_good, n_real
integer :: icount_eigval
logical, allocatable :: is_ok(:)
integer, allocatable :: list_degen(:,:), list_same_degen(:)
integer, allocatable :: iorder(:), list_degen_sorted(:)
double precision, allocatable :: leigvec_unsrtd(:,:), reigvec_unsrtd(:,:), eigval_unsrtd(:)
double precision, allocatable :: mat_tmp(:,:), eigval_tmp(:), leigvec_tmp(:,:), reigvec_tmp(:,:)
allocate(leigvec_unsrtd(n,n), reigvec_unsrtd(n,n), eigval_unsrtd(n))
leigvec_unsrtd = 0.d0
reigvec_unsrtd = 0.d0
eigval_unsrtd = 0.d0
! obtain degeneracies
allocate(list_degen(n,0:n))
call give_degen_full_list(fock_diag, n, thr_deg, list_degen, n_degen_list)
allocate(iorder(n_degen_list), list_degen_sorted(n_degen_list))
do i = 1, n_degen_list
n_degen = list_degen(i,0)
list_degen_sorted(i) = n_degen
iorder(i) = i
enddo
print*,'number of identical degeneracies = ',k
size_mat = k*n_degen
print*,'size_mat = ',size_mat
allocate(mat_tmp(size_mat,size_mat),list_same_degen(size_mat))
allocate(eigval_tmp(size_mat),leigvec_tmp(size_mat,size_mat),reigvec_tmp(size_mat,size_mat))
! group all the elements sharing the same degeneracy
icount = 0
do j = 1, k ! jth set of degeneracy
index_degen = iorder(i+j-1)
do m = 1, n_degen
icount += 1
list_same_degen(icount) = list_degen(index_degen,m)
enddo
! sort by number of degeneracies
call isort(list_degen_sorted, iorder, n_degen_list)
allocate(is_ok(n_degen_list))
is_ok = .True.
icount_eigval = 0
! loop over degeneracies
do i = 1, n_degen_list
if(.not.is_ok(i)) cycle
is_ok(i) = .False.
n_degen = list_degen_sorted(i)
print *, ' diagonalizing for n_degen = ', n_degen
k = 1
! group all the entries having the same degeneracies
!! do while (list_degen_sorted(i+k)==n_degen)
do m = i+1, n_degen_list
if(list_degen_sorted(m)==n_degen) then
is_ok(i+k) = .False.
k += 1
endif
enddo
print *, ' number of identical degeneracies = ', k
size_mat = k*n_degen
print *, ' size_mat = ', size_mat
allocate(mat_tmp(size_mat,size_mat), list_same_degen(size_mat))
allocate(eigval_tmp(size_mat), leigvec_tmp(size_mat,size_mat), reigvec_tmp(size_mat,size_mat))
! group all the elements sharing the same degeneracy
icount = 0
do j = 1, k ! jth set of degeneracy
index_degen = iorder(i+j-1)
do m = 1, n_degen
icount += 1
list_same_degen(icount) = list_degen(index_degen,m)
enddo
enddo
print *, ' list of elements '
do icount = 1, size_mat
print *, icount, list_same_degen(icount)
enddo
! you copy subset of matrix elements having all the same degeneracy in mat_tmp
do ii = 1, size_mat
i_good = list_same_degen(ii)
do jj = 1, size_mat
j_good = list_same_degen(jj)
mat_tmp(jj,ii) = mat_ref(j_good,i_good)
enddo
enddo
call non_hrmt_bieig( size_mat, mat_tmp, thr_d, thr_nd &
, leigvec_tmp, reigvec_tmp &
, n_real, eigval_tmp )
do ii = 1, size_mat
icount_eigval += 1
eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues
do jj = 1, size_mat ! copy the eigenvectors
j_good = list_same_degen(jj)
leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii)
reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii)
enddo
enddo
deallocate(mat_tmp, list_same_degen)
deallocate(eigval_tmp, leigvec_tmp, reigvec_tmp)
enddo
print*,'list of elements '
do icount = 1, size_mat
print*,icount,list_same_degen(icount)
enddo
! you copy subset of matrix elements having all the same degeneracy in mat_tmp
do ii = 1, size_mat
i_good = list_same_degen(ii)
do jj = 1, size_mat
j_good = list_same_degen(jj)
mat_tmp(jj,ii) = mat_ref(j_good,i_good)
enddo
enddo
call non_hrmt_bieig( size_mat, mat_tmp&
, leigvec_tmp, reigvec_tmp&
, n_real, eigval_tmp)
do ii = 1, size_mat
icount_eigval += 1
eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues
do jj = 1, size_mat ! copy the eigenvectors
j_good = list_same_degen(jj)
leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii)
reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii)
enddo
enddo
deallocate(mat_tmp,list_same_degen)
deallocate(eigval_tmp,leigvec_tmp,reigvec_tmp)
enddo
if(icount_eigval.ne.n)then
print*,'pb !! (icount_eigval.ne.n)'
print*,'icount_eigval,n',icount_eigval,n
stop
endif
if(icount_eigval .ne. n) then
print *, ' pb !! (icount_eigval.ne.n)'
print *, ' icount_eigval,n', icount_eigval, n
stop
endif
deallocate(iorder)
allocate(iorder(n))
do i = 1, n
iorder(i) = i
enddo
call dsort(eigval_unsrtd,iorder,n)
do i = 1, n
print*,'sorted eigenvalues '
i_good = iorder(i)
eigval(i) = eigval_unsrtd(i)
print*,'i,eigval(i) = ',i,eigval(i)
do j = 1, n
leigvec(j,i) = leigvec_unsrtd(j,i_good)
reigvec(j,i) = reigvec_unsrtd(j,i_good)
deallocate(iorder)
allocate(iorder(n))
do i = 1, n
iorder(i) = i
enddo
enddo
call dsort(eigval_unsrtd, iorder, n)
do i = 1, n
print*,'sorted eigenvalues '
i_good = iorder(i)
eigval(i) = eigval_unsrtd(i)
print*,'i,eigval(i) = ',i,eigval(i)
do j = 1, n
leigvec(j,i) = leigvec_unsrtd(j,i_good)
reigvec(j,i) = reigvec_unsrtd(j,i_good)
enddo
enddo
deallocate(leigvec_unsrtd, reigvec_unsrtd, eigval_unsrtd)
deallocate(list_degen)
deallocate(iorder, list_degen_sorted)
deallocate(is_ok)
end
! ---
subroutine give_degen_full_list(A,n,thr,list_degen,n_degen_list)
implicit none
BEGIN_DOC
! you enter with an array A(n) and spits out all the elements degenerated up to thr
!
! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL
!
! list_degen(i,0) = number of degenerate entries
!
! list_degen(i,1) = index of the first degenerate entry
!
! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries
!
! if list_degen(i,0) == 1 it means that there is no degeneracy for that element
END_DOC
double precision,intent(in) :: A(n)
double precision,intent(in) :: thr
integer, intent(in) :: n
integer, intent(out) :: list_degen(n,0:n),n_degen_list
logical, allocatable :: is_ok(:)
allocate(is_ok(n))
integer :: i,j,icount,icheck
n_degen_list = 0
is_ok = .True.
do i = 1, n
if(.not.is_ok(i))cycle
n_degen_list +=1
is_ok(i) = .False.
list_degen(n_degen_list,1) = i
icount = 1
do j = i+1, n
if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j))then
is_ok(j) = .False.
icount += 1
list_degen(n_degen_list,icount) = j
endif
BEGIN_DOC
! you enter with an array A(n) and spits out all the elements degenerated up to thr
!
! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL
!
! list_degen(i,0) = number of degenerate entries
!
! list_degen(i,1) = index of the first degenerate entry
!
! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries
!
! if list_degen(i,0) == 1 it means that there is no degeneracy for that element
END_DOC
implicit none
double precision, intent(in) :: A(n)
double precision, intent(in) :: thr
integer, intent(in) :: n
integer, intent(out) :: list_degen(n,0:n), n_degen_list
integer :: i, j, icount, icheck
logical, allocatable :: is_ok(:)
allocate(is_ok(n))
n_degen_list = 0
is_ok = .True.
do i = 1, n
if(.not.is_ok(i)) cycle
n_degen_list +=1
is_ok(i) = .False.
list_degen(n_degen_list,1) = i
icount = 1
do j = i+1, n
if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then
is_ok(j) = .False.
icount += 1
list_degen(n_degen_list,icount) = j
endif
enddo
list_degen(n_degen_list,0) = icount
enddo
list_degen(n_degen_list,0) = icount
enddo
icheck = 0
do i = 1, n_degen_list
icheck += list_degen(i,0)
enddo
if(icheck.ne.n)then
print*,'pb ! :: icheck.ne.n'
print*,icheck,n
stop
endif
icheck = 0
do i = 1, n_degen_list
icheck += list_degen(i,0)
enddo
if(icheck.ne.n)then
print *, ' pb ! :: icheck.ne.n'
print *, icheck, n
stop
endif
end
! ---