2019-07-09 16:17:10 +02:00
|
|
|
subroutine orthogonalization_matrix(ortho_type,nBas,S,X)
|
2019-03-13 11:07:31 +01:00
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
! Compute the orthogonalization matrix X
|
2019-03-13 11:07:31 +01:00
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
|
|
! Input variables
|
|
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
integer,intent(in) :: nBas,ortho_type
|
2019-03-13 11:07:31 +01:00
|
|
|
double precision,intent(in) :: S(nBas,nBas)
|
|
|
|
|
|
|
|
|
|
! Local variables
|
|
|
|
|
|
|
|
|
|
logical :: debug
|
|
|
|
|
double precision,allocatable :: UVec(:,:),Uval(:)
|
|
|
|
|
double precision,parameter :: thresh = 1d-6
|
|
|
|
|
|
|
|
|
|
integer :: i
|
|
|
|
|
|
|
|
|
|
! Output variables
|
|
|
|
|
|
|
|
|
|
double precision,intent(out) :: X(nBas,nBas)
|
|
|
|
|
|
|
|
|
|
debug = .false.
|
|
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
! Type of orthogonalization ortho_type
|
|
|
|
|
!
|
|
|
|
|
! 1 = Lowdin
|
|
|
|
|
! 2 = Canonical
|
|
|
|
|
! 3 = SVD
|
|
|
|
|
!
|
|
|
|
|
|
2019-03-13 11:07:31 +01:00
|
|
|
allocate(Uvec(nBas,nBas),Uval(nBas))
|
|
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
if(ortho_type == 1) then
|
|
|
|
|
|
|
|
|
|
write(*,*)
|
|
|
|
|
write(*,*) ' Lowdin orthogonalization'
|
|
|
|
|
write(*,*)
|
|
|
|
|
|
|
|
|
|
Uvec = S
|
|
|
|
|
call diagonalize_matrix(nBas,Uvec,Uval)
|
|
|
|
|
|
|
|
|
|
do i=1,nBas
|
|
|
|
|
|
|
|
|
|
if(Uval(i) > thresh) then
|
|
|
|
|
|
|
|
|
|
Uval(i) = 1d0/sqrt(Uval(i))
|
|
|
|
|
|
|
|
|
|
else
|
|
|
|
|
|
|
|
|
|
write(*,*) 'Eigenvalue',i,'too small for Lowdin orthogonalization'
|
2019-03-13 11:07:31 +01:00
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
endif
|
2019-03-13 11:07:31 +01:00
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
call ADAt(nBas,Uvec,Uval,X)
|
2019-03-13 11:07:31 +01:00
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
elseif(ortho_type == 2) then
|
2019-03-13 11:07:31 +01:00
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
write(*,*)
|
|
|
|
|
write(*,*) 'Canonical orthogonalization'
|
|
|
|
|
write(*,*)
|
|
|
|
|
|
|
|
|
|
Uvec = S
|
|
|
|
|
call diagonalize_matrix(nBas,Uvec,Uval)
|
|
|
|
|
|
|
|
|
|
do i=1,nBas
|
|
|
|
|
|
|
|
|
|
if(Uval(i) > thresh) then
|
2019-03-13 11:07:31 +01:00
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
Uval(i) = 1d0/sqrt(Uval(i))
|
2019-03-13 11:07:31 +01:00
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
else
|
2019-03-13 11:07:31 +01:00
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization'
|
2019-03-13 11:07:31 +01:00
|
|
|
|
2019-07-09 16:17:10 +02:00
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
call AD(nBas,Uvec,Uval)
|
|
|
|
|
X = Uvec
|
|
|
|
|
|
|
|
|
|
elseif(ortho_type == 3) then
|
|
|
|
|
|
|
|
|
|
write(*,*)
|
|
|
|
|
write(*,*) ' SVD-based orthogonalization NYI'
|
|
|
|
|
write(*,*)
|
|
|
|
|
|
|
|
|
|
! Uvec = S
|
|
|
|
|
! call diagonalize_matrix(nBas,Uvec,Uval)
|
|
|
|
|
|
|
|
|
|
! do i=1,nBas
|
|
|
|
|
! if(Uval(i) > thresh) then
|
|
|
|
|
! Uval(i) = 1d0/sqrt(Uval(i))
|
|
|
|
|
! else
|
|
|
|
|
! write(*,*) 'Eigenvalue',i,'too small for canonical orthogonalization'
|
|
|
|
|
! endif
|
|
|
|
|
! enddo
|
|
|
|
|
!
|
|
|
|
|
! call AD(nBas,Uvec,Uval)
|
|
|
|
|
! X = Uvec
|
|
|
|
|
|
|
|
|
|
endif
|
2019-03-13 11:07:31 +01:00
|
|
|
|
|
|
|
|
! Print results
|
|
|
|
|
|
|
|
|
|
if(debug) then
|
|
|
|
|
|
|
|
|
|
write(*,'(A28)') '----------------------'
|
|
|
|
|
write(*,'(A28)') 'Orthogonalization matrix'
|
|
|
|
|
write(*,'(A28)') '----------------------'
|
|
|
|
|
call matout(nBas,nBas,X)
|
|
|
|
|
write(*,*)
|
|
|
|
|
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
end subroutine orthogonalization_matrix
|
