debug: complex orthogonalization

src/HF/RHF.f90

@@ -89,6 +89,10 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
 ! Guess coefficients and density matrix
 
   call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c)
+  write(*,*)'Initial guess c'
+  call matout(nBas,nBas,c)
+  write(*,*) "verify orthonormalization"
+  call matout(nBas,nBas,matmul(transpose(c),matmul(S,c)))
 
   P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO)))
 ! call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
@@ -127,7 +131,8 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
     call exchange_matrix_AO_basis(nBas,P,ERI,K)
     
     F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
-
+    write(*,*)'Fock matrix'
+    call matout(nBas,nBas,F)
     ! Check convergence 
 
     err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)

src/HF/cRHF.f90

@@ -55,6 +55,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
   complex*16,allocatable        :: err(:,:)
   complex*16,allocatable        :: err_diis(:,:)
   complex*16,allocatable        :: F_diis(:,:)
+  complex*16,allocatable        :: Z(:,:)
 
 ! Output variables
 
@@ -80,7 +81,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
 
   allocate(J(nBas,nBas),K(nBas,nBas),err(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), &
            Fp(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis),     & 
-           Hc(nBas,nBas),W(nBas,nBas))
+           Hc(nBas,nBas),W(nBas,nBas),Z(nBas,nBas))
 
 ! Read CAP integrals from file
   call read_CAP_integrals(nBas,W)
@@ -94,7 +95,17 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
   
 
   call complex_mo_guess(nBas,nBas,guess_type,S,Hc,X,c)
+  write(*,*) "Guess coefficients c"
+  call complex_matout(nBas,nBas,c)
+  write(*,*) "Check if c orthonormal"
+  call complex_matout(nBas,nBas,matmul(transpose(c),matmul(S,c)))
+  c = matmul(c,Z)
+  call complex_orthogonalize_matrix(nBas,matmul(transpose(c),matmul(S,c)),Z)
+  write(*,*) "Check if c tilde orthonormal"
+
+  call complex_matout(nBas,nBas,matmul(transpose(Z),matmul(S,Z)))
   P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
+
 ! Initialization
 
   n_diis          = 0
@@ -108,7 +119,6 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
 !------------------------------------------------------------------------
 ! Main SCF loop
 !------------------------------------------------------------------------
-
   write(*,*)
   write(*,*)'--------------------------------------------------------------------------------------------------'
   write(*,'(1X,A1,1X,A3,1X,A1,1X,A36,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') &
@@ -125,9 +135,10 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
     
     call complex_Hartree_matrix_AO_basis(nBas,P,ERI,J)
     call complex_exchange_matrix_AO_basis(nBas,P,ERI,K)
-
-    F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
     
+    F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
+    write(*,*) "Fock matrix"
+    call complex_matout(nBas,nBas,F)
     ! Check convergence 
 
     err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
@@ -171,6 +182,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
 
     Fp = matmul(transpose(X),matmul(F,X))
     cp(:,:) = Fp(:,:)
+    write(*,*)'t1'
     call complex_diagonalize_matrix(nBas,Fp,eHF)
     c = matmul(X,cp)
     ! Density matrix

src/utils/complex_orthogonalize_matrix.f90

@@ -0,0 +1,72 @@
+subroutine complex_orthogonalize_matrix(nBas,S,X)
+
+! Compute the orthogonalization matrix X
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  complex*16,intent(in)         :: S(nBas,nBas)
+
+! Local variables
+
+  logical                       :: debug
+  complex*16,allocatable        :: UVec(:,:),Uval(:)
+  double precision,parameter    :: thresh = 1d-6
+  complex*16,allocatable        :: D(:,:)
+
+  integer                       :: i
+
+! Output variables
+
+  complex*16,intent(out)        :: X(nBas,nBas)
+
+  debug = .false.
+
+  allocate(Uvec(nBas,nBas),Uval(nBas),D(nBas,nBas))
+
+
+!   write(*,*)
+!   write(*,*) ' Lowdin orthogonalization'
+!   write(*,*)
+
+    Uvec = S
+    call complex_diagonalize_matrix(nBas,Uvec,Uval)
+
+    do i=1,nBas
+
+      if(abs(Uval(i)) < thresh) then 
+
+        write(*,*) 'Eigenvalue',i,' is very small in Lowdin orthogonalization = ',Uval(i)
+
+      endif
+      Uval(i) = 1d0/sqrt(Uval(i))
+
+    enddo
+    call diag(nBas,Uval, D)
+    X = matmul(Uvec,matmul(D,conjg(transpose(Uvec))))
+    deallocate(Uvec,Uval,D)
+end subroutine 
+
+subroutine diag(n,vec,M)
+! Create diag matrix from vector
+
+implicit none
+
+! input variables
+
+integer,intent(in) :: n
+complex*16,intent(in) :: vec(n)
+
+! local variables
+integer :: i
+
+! Output variables
+complex*16   :: M(n,n)
+
+M(:,:) = 0d0
+do i=1,n
+  M(i,i) = vec(i)
+enddo
+end subroutine

src/utils/wrap_lapack.f90

@@ -143,15 +143,15 @@ subroutine complex_diagonalize_matrix(N,A,e)
 ! Local variables
 
   integer                       :: lwork,info
-  double precision              :: rwork(2*N)
+  double precision,allocatable  :: rwork(:)
   complex*16,allocatable        :: work(:)
-  complex*16                    :: VL(N,N)
-  complex*16                    :: VR(N,N)
+  complex*16,allocatable        :: VL(:,:)
+  complex*16,allocatable        :: VR(:,:)
 
 ! Memory allocation
-  allocate(work(1))
+  allocate(work(1),rwork(2*N),VL(1,1),VR(N,N))
   lwork = -1
-  call zgeev('N','V',N,A,N,e,VL,N,VR,N,work,lwork,rwork,info)
+  call zgeev('N','V',N,A,N,e,VL,1,VR,N,work,lwork,rwork,info)
   lwork = max(1,int(real(work(1))))
   
   deallocate(work)