Cholesky for orthogonalization instead GS

crhf_test/cRHF.py

@@ -98,7 +98,7 @@ def sort_eigenpairs(eigenvalues, eigenvectors):
     return sorted_eigenvalues, sorted_eigenvectors
 
 
-def diagonalize(M):
+def diagonalize_gram_schmidt(M):
     # Diagonalize the matrix
     vals, vecs = np.linalg.eig(M)
     # Sort the eigenvalues and eigenvectors
@@ -108,6 +108,25 @@ def diagonalize(M):
     return vals, vecs
 
 
+def diagonalize(M):
+    # Diagonalize the matrix
+    vals, vecs = np.linalg.eig(M)
+    # Sort the eigenvalues and eigenvectors
+    vals, vecs = sort_eigenpairs(vals, vecs)
+    # Orthonormalize them wrt cTc inner product
+    vecs = orthonormalize(vecs)
+    return vals, vecs
+
+
+def orthonormalize(vecs):
+    # Orthonormalize them wrt cTc inner product
+    R = vecs.T@vecs
+    L = cholesky_decomposition(R)
+    Linv = np.linalg.inv(L)
+    vecs = vecs@Linv.T
+    return vecs
+
+
 def Hartree_matrix_AO_basis(P, ERI):
     # Initialize Hartree matrix with zeros (complex type)
     J = np.zeros((nBas, nBas), dtype=np.complex128)
@@ -151,13 +170,71 @@ def gram_schmidt(vectors):
     return np.column_stack(orthonormal_basis)
 
 
+def DIIS_extrapolation(rcond, n_diis, error, e, error_in, e_inout):
+    """
+    Perform DIIS extrapolation.
+
+    """
+
+    # Update DIIS history by prepending new error and solution vectors
+    error = np.column_stack((error_in, error[:, :-1]))  # Shift history
+    e = np.column_stack((e_inout, e[:, :-1]))          # Shift history
+
+    # Build A matrix
+    A = np.zeros((n_diis + 1, n_diis + 1), dtype=np.complex128)
+    print(np.shape(error))
+    A[:n_diis, :n_diis] = error@error.T
+    A[:n_diis, n_diis] = -1.0
+    A[n_diis, :n_diis] = -1.0
+    A[n_diis, n_diis] = 0.0
+
+    # Build b vector
+    b = np.zeros(n_diis + 1, dtype=np.complex128)
+    b[n_diis] = -1.0
+
+    # Solve the linear system A * w = b
+    try:
+        w = np.linalg.solve(A, b)
+        rcond = np.linalg.cond(A)
+    except np.linalg.LinAlgError:
+        raise ValueError("DIIS linear system is singular or ill-conditioned.")
+
+    # Extrapolate new solution
+    e_inout[:] = w[:n_diis]@e[:, :n_diis].T
+
+    return rcond, n_diis, e_inout
+
+
+def cholesky_decomposition(A):
+    """
+    Performs Cholesky-Decomposition wrt the c product. Returns L such that A = LTL
+    """
+
+    L = np.zeros_like(A)
+    n = np.shape(L)[0]
+    for i in range(n):
+        for j in range(i + 1):
+            s = A[i, j]
+
+            for k in range(j):
+                s -= L[i, k] * L[j, k]
+
+            if i > j:  # Off-diagonal elements
+                L[i, j] = s / L[j, j]
+            else:  # Diagonal elements
+                L[i, i] = s**0.5
+
+    return L
+
+
 if __name__ == "__main__":
-    # Constants
+    # Inputs
     workdir = "../"
     eta = 0.01
     thresh = 0.00001
     maxSCF = 256
     nO = read_nO(workdir)
+    maxDiis = 5
 
     # Read integrals
     T = read_matrix("../int/Kin.dat")
@@ -165,12 +242,18 @@ if __name__ == "__main__":
     V = read_matrix("../int/Nuc.dat")
     ENuc = read_ENuc("../int/ENuc.dat")
     nBas = np.shape(T)[0]
+    nBasSq = nBas*nBas
     W = read_CAP_integrals("../int/CAP.dat", nBas)
     ERI = read_2e_integrals("../int/ERI.bin", nBas)
     X = get_X(S)
     W = -eta*W
     Hc = T + V + W*1j
 
+    # Initialization
+    F_diis = np.zeros((nBasSq, maxDiis))
+    error_diis = np.zeros((nBasSq, maxDiis))
+    rcond = 0
+
     # core guess
     _, c = diagonalize(X.T @ Hc @ X)
     c = X @ c
@@ -183,6 +266,8 @@ if __name__ == "__main__":
 
     nSCF = 0
     Conv = 1
+    n_diis = 0
+
     while(Conv > thresh and nSCF < maxSCF):
         nSCF += 1
         J = Hartree_matrix_AO_basis(P, ERI)
@@ -197,6 +282,11 @@ if __name__ == "__main__":
         EK = 0.25*np.trace(P@K)
         ERHF = ET + EV + EJ + EK
 
+      #  # DIIS
+      #  n_diis = np.min([n_diis+1, maxDiis])
+      #  rcond, n_diis, F = DIIS_extrapolation(
+      #      rcond, n_diis, error_diis, F_diis, err, F)
+
         Fp = X.T @ F @ X
         eHF, c = diagonalize(Fp)
         c = X @ c

src/HF/cRHF.f90

@@ -167,6 +167,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
     Fp = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
     cp(:,:) = Fp(:,:)
     call complex_diagonalize_matrix(nBas,cp,eHF)
+    call complex_orthogonalize_matrix(nBas,cp)
     c = matmul(X,cp)
     ! Density matrix
     P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
@@ -196,9 +197,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
   end if
 
 ! Compute dipole moments
-
+  call print_cRHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF)
-  call complex_dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
-  call print_cRHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole)
 ! Testing zone
 
   if(dotest) then

src/HF/complex_core_guess.f90

@@ -27,6 +27,7 @@ subroutine complex_core_guess(nBas, nOrb, Hc, X, c)
 
   cp(:,:) = matmul(transpose(X(:,:)), matmul(Hc(:,:), X(:,:)))
   call complex_diagonalize_matrix(nOrb, cp, e)
+  call complex_orthogonalize_matrix(nOrb,cp)
   c(:,:) = matmul(X(:,:), cp(:,:))
   deallocate(cp, e)
 

src/HF/print_cRHF.f90

@@ -1,7 +1,7 @@
 
 ! ---
 
-subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipole)
+subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF)
 
 ! Print one-electron energies and other stuff for G0W0
 
@@ -20,7 +20,6 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipo
   complex*16,intent(in)              :: EK
   complex*16,intent(in)              :: EV
   complex*16,intent(in)              :: ERHF
-  double precision,intent(in)        :: dipole(ncart)
 
 ! Local variables
 
@@ -65,12 +64,6 @@ subroutine print_cRHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipo
   write(*,'(A50)')           '------------------------------------------------------------'
   write(*,'(A33,1X,F16.6)')     ' <Sz>                = ',S
   write(*,'(A33,1X,F16.6)')     ' <S^2>               = ',S2
-  write(*,'(A50)')           '------------------------------------------------------------'
-  write(*,'(A36)')           ' Dipole moment (Debye)    '
-  write(*,'(10X,4A10)')      'X','Y','Z','Tot.'
-  write(*,'(10X,4F10.4)')    (real(dipole(ixyz))*auToD,ixyz=1,ncart),norm2(real(dipole))*auToD
-  write(*,'(A50)')           '---------------------------------------------'
-  write(*,*)
 
 ! Print results
 

src/utils/complex_DIIS_extrapolation.f90

@@ -18,7 +18,6 @@ subroutine complex_DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_
   complex*16,allocatable        :: A(:,:)
   complex*16,allocatable        :: b(:)
   complex*16,allocatable        :: w(:)
-  double precision,allocatable  :: b2(:)
 
 ! Output variables
 
@@ -27,8 +26,6 @@ subroutine complex_DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_
   complex*16,intent(inout)      :: e_inout(n_e)
 
 ! Memory allocation
-  allocate(b2(n_diis+1))
-  print *, "Allocated b2"
   allocate(A(n_diis+1,n_diis+1),b(n_diis+1),w(n_diis+1))
 
 ! Update DIIS "history"
@@ -50,7 +47,8 @@ subroutine complex_DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_
   b(n_diis+1) = cmplx(-1d0,0d0,kind=8)
 
 ! Solve linear system
-
+  call complex_vecout(b)
+  call complex_matout(A)
   call complex_linear_solve(n_diis+1,A,b,w,rcond)
 
 ! Extrapolate

src/utils/complex_gram_schmidt.f90

@@ -1,55 +0,0 @@
-subroutine complex_gram_schmidt(N, vectors)
-  
-  ! Input variables
-  implicit none
-  integer, intent(in) :: N
-  complex*16, intent(inout) :: vectors(N, N)
-
-  ! Local variables
-  integer :: i, j
-  complex*16 :: Mtmp(N,N)
-  complex*16 :: proj
-  complex*16 :: norm
-
-  ! Copy the input matrix to a temporary matrix
-  Mtmp(:, :) = vectors(:,:)
-  
-  ! Orthonormalize the vectors
-  do i = 1, N
-
-      ! Orthonormalize with respect to the previous vectors in Mtmp
-      do j = 1, i-1
-          ! Compute the dot product (scalar product) of vectors j and i
-          call complex_dot_product(N,Mtmp(:, j), Mtmp(:, i),proj)
-
-          ! Subtract the projection onto the j-th vector
-          Mtmp(:, i) = Mtmp(:, i) - proj * Mtmp(:, j)
-      end do
-
-      ! Normalize the vector
-      call complex_dot_product(N,Mtmp(:, i), Mtmp(:, i),proj)
-      norm = sqrt(proj)
-
-      if (abs(norm) > 1.0d-10) then
-          ! Normalize the i-th vector and store it back in vectors
-          vectors(:, i) = Mtmp(:, i) / norm
-      else
-          print*, "Error: Norm of eigenvector < 1e-10 !!!"
-          stop
-      end if
-  end do
-end subroutine
-
-subroutine complex_dot_product(N,vec1,vec2,proj)
-  ! Input
-  complex*16,intent(in) :: vec1(N),vec2(N)
-  !Output
-  complex*16, intent(out) :: proj
-  ! Local variables
-  integer :: i
-  
-  proj = 0d0
-  do i=1,N
-     proj = proj +vec1(i)*vec2(i)
-  end do
-end subroutine

src/utils/complex_orthogonalization.f90

@@ -0,0 +1,109 @@
+subroutine complex_orthogonalize_matrix(N, vectors)
+  
+  ! Input variables
+  implicit none
+  integer, intent(in) :: N
+  complex*16, intent(inout) :: vectors(N, N)
+
+  ! Local variables
+  integer :: i, j
+  complex*16,allocatable :: L(:,:),Linv(:,:)
+  complex*16 :: proj
+  complex*16 :: norm
+
+  ! Copy the input matrix to a temporary matrix
+  allocate(L(N,N),Linv(N,N))
+  L = matmul(transpose(vectors) ,vectors)
+  call complex_cholesky_decomp(N,L)
+  call complex_inverse_matrix(N,L,Linv)
+  vectors = matmul(vectors,transpose(Linv))
+  deallocate(L,Linv)
+end subroutine
+subroutine complex_gram_schmidt(N, vectors)
+  
+  ! Input variables
+  implicit none
+  integer, intent(in) :: N
+  complex*16, intent(inout) :: vectors(N, N)
+
+  ! Local variables
+  integer :: i, j
+  complex*16,allocatable :: Mtmp(:,:)
+  complex*16 :: proj
+  complex*16 :: norm
+
+  ! Copy the input matrix to a temporary matrix
+  allocate(Mtmp(N,N))
+  Mtmp(:, :) = vectors(:,:)
+  
+  ! Orthonormalize the vectors
+  do i = 1, N
+
+      ! Orthonormalize with respect to the previous vectors in Mtmp
+      do j = 1, i-1
+          ! Compute the dot product (scalar product) of vectors j and i
+          call complex_dot_product(N,Mtmp(:, j), Mtmp(:, i),proj)
+
+          ! Subtract the projection onto the j-th vector
+          Mtmp(:, i) = Mtmp(:, i) - proj * Mtmp(:, j)
+      end do
+
+      ! Normalize the vector
+      call complex_dot_product(N,Mtmp(:, i), Mtmp(:, i),proj)
+      norm = sqrt(proj)
+
+      if (abs(norm) > 1.0d-10) then
+          ! Normalize the i-th vector and store it back in vectors
+          vectors(:, i) = Mtmp(:, i) / norm
+      else
+          print*, "Error: Norm of eigenvector < 1e-10 !!!"
+          stop
+      end if
+  end do
+  deallocate(Mtmp)
+end subroutine
+
+subroutine complex_dot_product(N,vec1,vec2,proj)
+  ! Input
+  complex*16,intent(in) :: vec1(N),vec2(N)
+  !Output
+  complex*16, intent(out) :: proj
+  ! Local variables
+  integer :: i
+  
+  proj = 0d0
+  do i=1,N
+     proj = proj +vec1(i)*vec2(i)
+  end do
+end subroutine
+
+subroutine complex_cholesky_decomp(n,A)
+    implicit none
+    integer, intent(in) :: n               ! Matrix size
+    complex*16, intent(inout) :: A(n, n)   ! Output: Lower triangular Cholesky factor L
+    integer :: i, j, k
+    complex*16 :: s
+
+    ! Perform Cholesky decomposition
+    do i = 1, n
+        do j = 1, i
+            s = A(i, j)
+            do k = 1, j - 1
+                s = s - A(i, k) * A(j, k)
+            end do
+
+            if (i > j) then
+                A(i, j) = s / A(j, j)  ! Compute lower triangular elements
+            else
+                A(i, i) = sqrt(s)
+            end if
+        end do
+    end do
+
+    ! Zero out upper triangular part
+    do i = 1, n
+        do j = i + 1, n
+            A(i, j) = cmplx(0.0d0, 0.0d0, kind=8)
+        end do
+    end do
+end subroutine

src/utils/complex_orthogonalize_matrix.f90

@@ -1,60 +0,0 @@
-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))
-
-    Uvec = S
-    call complex_diagonalize_matrix(nBas,Uvec,Uval)
-    call complex_matout(nBas,nBas,matmul(Uvec,transpose(Uvec)))
-    do i=1,nBas
-      Uval(i) = 1d0/sqrt(Uval(i))
-    enddo
-    call diag(nBas,Uval, D)
-    X = matmul(Uvec,matmul(D,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

@@ -159,7 +159,8 @@ subroutine complex_diagonalize_matrix(N,A,e)
 
   call zgeev('N','V',N,A,N,e,VL,N,VR,N,work,lwork,rwork,info)
   call complex_sort_eigenvalues(N,e,VR)
-  call complex_gram_schmidt(N,VR)
+  
+
   deallocate(work)
   A = VR
 
@@ -256,6 +257,46 @@ subroutine inverse_matrix(N,A,B)
 
 end subroutine 
 
+subroutine complex_inverse_matrix(N,A,B)
+
+! Returns the inverse of the complex square matrix A in B
+
+  implicit none
+
+  integer,intent(in)             :: N
+  complex*16, intent(in)         :: A(N,N)
+  complex*16, intent(out)        :: B(N,N)
+
+  integer                        :: info,lwork
+  integer, allocatable           :: ipiv(:)
+  complex*16,allocatable         :: work(:)
+  
+  allocate (ipiv(N),work(N*N))
+  lwork = size(work)
+
+  B(1:N,1:N) = A(1:N,1:N)
+
+  call zgetrf(N,N,B,N,ipiv,info)
+
+  if (info /= 0) then
+
+    print*,info
+    stop 'error in inverse (zgetri)!!'
+
+  end if
+
+  call zgetri(N,B,N,ipiv,work,lwork,info)
+
+  if (info /= 0) then
+
+    print *,  info
+    stop 'error in inverse (zgetri)!!'
+
+  end if
+
+  deallocate(ipiv,work)
+
+end subroutine 
 subroutine linear_solve(N,A,b,x,rcond)
 
 ! Solve the linear system A.x = b where A is a NxN matrix
@@ -315,8 +356,8 @@ subroutine complex_linear_solve(N,A,b,x,rcond)
   ! Find optimal size for temporary arrays
 
   allocate(ipiv(N))
-
   call zgesv(N,1,A,N,ipiv,b,N,info)
+
 end subroutine
 
 subroutine easy_linear_solve(N,A,b,x)