diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index b01002f..df2b8bc 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -30,7 +30,7 @@ program QuAcK double precision,allocatable :: T(:,:) double precision,allocatable :: V(:,:) double precision,allocatable :: Hc(:,:) - double precision,allocatable :: X(:,:) + double precision,allocatable :: X(:,:),X_tmp(:,:) double precision,allocatable :: dipole_int_AO(:,:,:) double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: Uvec(:,:), Uval(:) @@ -71,10 +71,6 @@ program QuAcK logical :: dotest,doRtest,doUtest,doGtest - integer :: i, j, j0 - double precision :: acc_d, acc_nd - double precision, allocatable :: tmp1(:,:), tmp2(:,:) - !-------------! ! Hello World ! !-------------! @@ -180,61 +176,11 @@ program QuAcK ! Compute orthogonalization matrix - !call orthogonalization_matrix(nBas, S, X) - - allocate(Uvec(nBas,nBas), Uval(nBas)) - - Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas) - call diagonalize_matrix(nBas, Uvec, Uval) - - nOrb = 0 - do i = 1, nBas - if(Uval(i) > 1d-6) then - Uval(i) = 1d0 / dsqrt(Uval(i)) - nOrb = nOrb + 1 - else - write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' - end if - end do - - write(*,'(A38)') '--------------------------------------' - write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas - write(*,'(A38,1X,I16)') 'Number of basis functions (MOs)', nOrb - write(*,'(A38,1X,F9.3)') ' % of discarded orbitals = ', 100.d0 * (1.d0 - dble(nOrb)/dble(nBas)) - write(*,'(A38)') '--------------------------------------' - write(*,*) - - j0 = nBas - nOrb + allocate(X_tmp(nBas,nBas)) + call orthogonalization_matrix(nBas,nOrb,S,X_tmp) allocate(X(nBas,nOrb)) - do j = j0+1, nBas - do i = 1, nBas - X(i,j-j0) = Uvec(i,j) * Uval(j) - enddo - enddo - - deallocate(Uvec, Uval) - - !! check if X.T S X = 1_(nOrb,nOrb) - !allocate(tmp1(nOrb,nBas), tmp2(nOrb,nOrb)) - !call dgemm("T", "N", nOrb, nBas, nBas, 1.d0, & - ! X(1,1), nBas, S(1,1), nBas, & - ! 0.d0, tmp1(1,1), nOrb) - !call dgemm("N", "N", nOrb, nOrb, nBas, 1.d0, & - ! tmp1(1,1), nOrb, X(1,1), nBas, & - ! 0.d0, tmp2(1,1), nOrb) - !acc_d = 0.d0 - !acc_nd = 0.d0 - !do i = 1, nOrb - ! !write(*,'(1000(F15.7,2X))') (tmp2(i,j), j = 1, nOrb) - ! acc_d = acc_d + tmp2(i,i) - ! do j = 1, nOrb - ! if(j == i) cycle - ! acc_nd = acc_nd + dabs(tmp2(j,i)) - ! enddo - !enddo - !print*, ' diag part: ', dabs(acc_d - dble(nOrb)) / dble(nOrb) - !print*, ' non-diag part: ', acc_nd - !deallocate(tmp1, tmp2) + X(1:nBas,1:nOrb) = X_tmp(1:nBas,1:nOrb) + deallocate(X_tmp) !---------------------! ! Choose QuAcK branch ! diff --git a/src/utils/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90 index af781e1..fd9e59c 100644 --- a/src/utils/orthogonalization_matrix.f90 +++ b/src/utils/orthogonalization_matrix.f90 @@ -1,7 +1,4 @@ - -! --- - -subroutine orthogonalization_matrix(nBas, S, X) +subroutine orthogonalization_matrix(nBas,nOrb,S,X) ! Compute the orthogonalization matrix X @@ -17,113 +14,47 @@ subroutine orthogonalization_matrix(nBas, S, X) logical :: debug double precision,allocatable :: UVec(:,:),Uval(:) double precision,parameter :: thresh = 1d-6 - integer,parameter :: ortho_type = 1 - integer :: i + integer :: i,j,j0 ! Output variables + integer :: nOrb double precision,intent(out) :: X(nBas,nBas) debug = .false. -! Type of orthogonalization ortho_type -! -! 1 = Lowdin -! 2 = Canonical -! 3 = SVD -! - allocate(Uvec(nBas,nBas),Uval(nBas)) - if(ortho_type == 1) then - - ! - ! S V = V s where - ! - ! V.T V = 1 and s > 0 (S is positive def) - ! - ! S = V s V.T - ! = V s^0.5 s^0.5 V.T - ! = V s^0.5 V.T V s^0.5 V.T - ! = S^0.5 S^0.5 - ! - ! where - ! - ! S^0.5 = V s^0.5 V.T - ! - ! X = S^(-0.5) - ! = V s^(-0.5) V.T - ! - -! write(*,*) -! write(*,*) ' Lowdin orthogonalization' -! write(*,*) - - Uvec = S - call diagonalize_matrix(nBas, Uvec, Uval) - - do i = 1, nBas - - if(Uval(i) < thresh) then - - write(*,*) 'Eigenvalue',i,' is very small in Lowdin orthogonalization = ',Uval(i) - - end if - - Uval(i) = 1d0 / dsqrt(Uval(i)) - - end do - - call ADAt(nBas, Uvec(1,1), Uval(1), X(1,1)) - - elseif(ortho_type == 2) then - -! write(*,*) -! write(*,*) 'Canonical orthogonalization' -! write(*,*) - - Uvec = S - call diagonalize_matrix(nBas, Uvec, Uval) - - do i = 1, nBas - - if(Uval(i) > thresh) then + Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas) + call diagonalize_matrix(nBas,Uvec,Uval) + nOrb = 0 + do i = 1,nBas + if(Uval(i) > thresh) then Uval(i) = 1d0 / dsqrt(Uval(i)) + nOrb = nOrb + 1 + else + write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' + end if + end do - else + write(*,'(A50)') '------------------------------------------------' + write(*,'(A40,1X,I5)') 'Number of basis functions = ',nBas + write(*,'(A40,1X,I5)') 'Number of spatial orbitals = ',nOrb + write(*,'(A40,1X,I5)') 'Number of discarded functions = ',nBas - nOrb + write(*,'(A50)') '------------------------------------------------' + write(*,*) - write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' + j0 = nBas - nOrb - end if + do j = j0+1,nBas + do i = 1,nBas + X(i,j-j0) = Uvec(i,j) * Uval(j) + enddo + enddo - end do - - 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' -! end if -! end do -! -! call AD(nBas,Uvec,Uval) -! X = Uvec - - end if + deallocate(Uvec,Uval) ! Print results