quantum_package/plugins/analyze_wf/phi_s.irp.f

129 lines
3.9 KiB
Fortran

BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvalues, (mo_tot_num, 2:N_states) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvectors, (mo_tot_num, mo_tot_num, 2:N_states) ]
implicit none
BEGIN_DOC
! Eigenvalues and eigenvectors of one_body_dm_mo_diff
END_DOC
integer :: i,j,istate
integer :: liwork, lwork, n, info
integer, allocatable :: iwork(:)
double precision, allocatable :: work(:)
one_body_dm_mo_diff_eigvectors(1:mo_tot_num, 1:mo_tot_num, 2:N_states) =&
one_body_dm_mo_diff(1:mo_tot_num, 1:mo_tot_num, 2:N_states)
n = mo_tot_num
lwork = 1+6*n + 2*n*n
liwork = 3 + 5*n
allocate(work(lwork), iwork(liwork))
lwork=-1
liwork=-1
istate=2
call dsyevd( 'V', 'U', mo_tot_num, &
one_body_dm_mo_diff_eigvectors(1,1,istate), &
size(one_body_dm_mo_diff_eigvectors,1), &
one_body_dm_mo_diff_eigvalues(1,istate), &
work, lwork, iwork, liwork, info)
if (info /= 0) then
print *, irp_here//' DSYEVD failed : ', info
stop 1
endif
lwork = int(work(1))
liwork = iwork(1)
deallocate(iwork,work)
allocate(work(lwork), iwork(liwork))
do istate=2,N_states
call dsyevd( 'V', 'U', mo_tot_num, &
one_body_dm_mo_diff_eigvectors(1,1,istate), &
size(one_body_dm_mo_diff_eigvectors,1), &
one_body_dm_mo_diff_eigvalues(1,istate), &
work, lwork, iwork, liwork, info)
if (info /= 0) then
print *, irp_here//' DSYEVD failed : ', info
stop 1
endif
enddo
deallocate(iwork,work)
END_PROVIDER
BEGIN_PROVIDER [ double precision, transition_natorb, (ao_num,mo_tot_num,2:N_states) ]
implicit none
BEGIN_DOC
! Natural transition molecular orbitals
END_DOC
integer :: istate
do istate=2,N_states
call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, &
mo_coef, size(mo_coef,1), &
one_body_dm_mo_diff_eigvectors(1,1,istate), &
size(one_body_dm_mo_diff_eigvectors,1), 0.d0, &
transition_natorb(1,1,istate), size(transition_natorb,1))
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, phi_s, (2:N_states) ]
implicit none
BEGIN_DOC
!
END_DOC
integer :: i,istate
double precision, allocatable :: T(:,:), A(:,:), D(:,:)
double precision :: trace, norm
allocate(T(ao_num,ao_num), A(ao_num,ao_num), D(ao_num,ao_num))
do istate=2,N_states
call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, &
S_half, size(S_half,1), &
one_body_dm_ao_attachment(1,1,istate), size(one_body_dm_ao_attachment,1), 0.d0,&
T, size(T,1))
call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, &
T, size(T,1), &
S_half, size(S_half,1), 0.d0, &
A, size(A,1))
!
call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, &
S_half, size(S_half,1), &
one_body_dm_ao_detachment(1,1,istate), size(one_body_dm_ao_detachment,1), 0.d0,&
T, size(T,1))
call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, &
T, size(T,1), &
S_half, size(S_half,1), 0.d0, &
D, size(D,1))
trace = 0.d0
do i=1,ao_num
trace = trace + A(i,i)
enddo
norm = 0.d0
do i=1,ao_num
norm = norm + D(i,i)
enddo
norm = 0.5d0*(norm + trace)
trace = 0.d0
do i=1,mo_tot_num
trace = trace + dsqrt(A(i,i)*D(i,i))
enddo
phi_s(istate) = trace/norm
enddo
END_PROVIDER