move orthogonalization matrix

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(X_tmp(nBas,nBas))
-
+  call orthogonalization_matrix(nBas,nOrb,S,X_tmp)
-  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(nBas,nOrb))
-  do j = j0+1, nBas
+  X(1:nBas,1:nOrb) = X_tmp(1:nBas,1:nOrb)
-    do i = 1, nBas
+  deallocate(X_tmp)
-      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)
 
 !---------------------!
 ! Choose QuAcK branch !

src/utils/orthogonalization_matrix.f90

@@ -1,7 +1,4 @@
-
+subroutine orthogonalization_matrix(nBas,nOrb,S,X)
-! ---
-
-subroutine orthogonalization_matrix(nBas, 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
+  Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas)
-
+  call diagonalize_matrix(nBas,Uvec,Uval)
-    !
-    ! 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 
 
+  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
+  deallocate(Uvec,Uval)
-    
-    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
 
 ! Print results