mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-07 14:03:37 +01:00
177 lines
5.1 KiB
Fortran
177 lines
5.1 KiB
Fortran
|
|
|
|
BEGIN_PROVIDER [ double precision, overlap_bi_ortho, (mo_num, mo_num)]
|
|
&BEGIN_PROVIDER [ double precision, overlap_diag_bi_ortho, (mo_num)]
|
|
|
|
BEGIN_DOC
|
|
! Overlap matrix between the RIGHT and LEFT MOs. Should be the identity matrix
|
|
END_DOC
|
|
|
|
implicit none
|
|
integer :: i, k, m, n
|
|
double precision :: accu_d, accu_nd
|
|
double precision, allocatable :: tmp(:,:)
|
|
|
|
! overlap_bi_ortho = 0.d0
|
|
! do i = 1, mo_num
|
|
! do k = 1, mo_num
|
|
! do m = 1, ao_num
|
|
! do n = 1, ao_num
|
|
! overlap_bi_ortho(k,i) += ao_overlap(n,m) * mo_l_coef(n,k) * mo_r_coef(m,i)
|
|
! enddo
|
|
! enddo
|
|
! enddo
|
|
! enddo
|
|
|
|
allocate( tmp(mo_num,ao_num) )
|
|
! tmp <-- L.T x S_ao
|
|
call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 &
|
|
, mo_l_coef(1,1), size(mo_l_coef, 1), ao_overlap(1,1), size(ao_overlap, 1) &
|
|
, 0.d0, tmp(1,1), size(tmp, 1) )
|
|
! S <-- tmp x R
|
|
call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 &
|
|
, tmp(1,1), size(tmp, 1), mo_r_coef(1,1), size(mo_r_coef, 1) &
|
|
, 0.d0, overlap_bi_ortho(1,1), size(overlap_bi_ortho, 1) )
|
|
deallocate(tmp)
|
|
|
|
do i = 1, mo_num
|
|
overlap_diag_bi_ortho(i) = overlap_bi_ortho(i,i)
|
|
enddo
|
|
|
|
accu_d = 0.d0
|
|
accu_nd = 0.d0
|
|
do i = 1, mo_num
|
|
do k = 1, mo_num
|
|
if(i==k) then
|
|
accu_d += dabs(overlap_bi_ortho(k,i))
|
|
else
|
|
accu_nd += dabs(overlap_bi_ortho(k,i))
|
|
endif
|
|
enddo
|
|
enddo
|
|
accu_d = accu_d/dble(mo_num)
|
|
accu_nd = accu_nd/dble(mo_num**2-mo_num)
|
|
if(dabs(accu_d-1.d0).gt.1.d-10.or.dabs(accu_nd).gt.1.d-10)then
|
|
print*,'Warning !!!'
|
|
print*,'Average trace of overlap_bi_ortho is different from 1 by ', dabs(accu_d-1.d0)
|
|
print*,'And bi orthogonality is off by an average of ',accu_nd
|
|
print*,'****************'
|
|
!print*,'Overlap matrix betwee mo_l_coef and mo_r_coef '
|
|
!do i = 1, mo_num
|
|
! write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:)
|
|
!enddo
|
|
endif
|
|
print*,'Average trace of overlap_bi_ortho (should be 1.)'
|
|
print*,'accu_d = ',accu_d
|
|
print*,'Sum of off diagonal terms of overlap_bi_ortho (should be zero)'
|
|
print*,'accu_nd = ',accu_nd
|
|
print*,'****************'
|
|
|
|
END_PROVIDER
|
|
|
|
! ---
|
|
|
|
BEGIN_PROVIDER [ double precision, overlap_mo_r, (mo_num, mo_num)]
|
|
&BEGIN_PROVIDER [ double precision, overlap_mo_l, (mo_num, mo_num)]
|
|
|
|
BEGIN_DOC
|
|
! overlap_mo_r_mo(j,i) = <MO_i|MO_R_j>
|
|
END_DOC
|
|
|
|
implicit none
|
|
integer :: i, j, p, q
|
|
double precision, allocatable :: tmp(:,:)
|
|
|
|
!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
|
|
|
|
allocate( tmp(mo_num,ao_num) )
|
|
|
|
tmp = 0.d0
|
|
call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 &
|
|
, mo_r_coef(1,1), size(mo_r_coef, 1), ao_overlap(1,1), size(ao_overlap, 1) &
|
|
, 0.d0, tmp(1,1), size(tmp, 1) )
|
|
call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 &
|
|
, tmp(1,1), size(tmp, 1), mo_r_coef(1,1), size(mo_r_coef, 1) &
|
|
, 0.d0, overlap_mo_r(1,1), size(overlap_mo_r, 1) )
|
|
|
|
tmp = 0.d0
|
|
call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 &
|
|
, mo_l_coef(1,1), size(mo_l_coef, 1), ao_overlap(1,1), size(ao_overlap, 1) &
|
|
, 0.d0, tmp(1,1), size(tmp, 1) )
|
|
call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 &
|
|
, tmp(1,1), size(tmp, 1), mo_l_coef(1,1), size(mo_l_coef, 1) &
|
|
, 0.d0, overlap_mo_l(1,1), size(overlap_mo_l, 1) )
|
|
|
|
deallocate(tmp)
|
|
|
|
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)]
|
|
|
|
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
|
|
|
|
END_PROVIDER
|
|
|
|
! ---
|
|
|
|
BEGIN_PROVIDER [ double precision, angle_left_right, (mo_num)]
|
|
&BEGIN_PROVIDER [ double precision, max_angle_left_right]
|
|
|
|
BEGIN_DOC
|
|
! angle_left_right(i) = angle between the left-eigenvector chi_i and the right-eigenvector phi_i
|
|
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
|
|
|
|
! ---
|
|
|
|
|