diff --git a/crhf_test/cRHF.py b/crhf_test/cRHF.py index f4fc729..17b9c2e 100755 --- a/crhf_test/cRHF.py +++ b/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 diff --git a/src/HF/cRHF.f90 b/src/HF/cRHF.f90 index 0f42ba9..ad6bb04 100644 --- a/src/HF/cRHF.f90 +++ b/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 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) + call print_cRHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF) ! Testing zone if(dotest) then diff --git a/src/HF/complex_core_guess.f90 b/src/HF/complex_core_guess.f90 index 622cb36..2754fee 100644 --- a/src/HF/complex_core_guess.f90 +++ b/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) diff --git a/src/HF/print_cRHF.f90 b/src/HF/print_cRHF.f90 index 04808b5..d8de170 100644 --- a/src/HF/print_cRHF.f90 +++ b/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)') ' = ',S write(*,'(A33,1X,F16.6)') ' = ',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 diff --git a/src/utils/complex_DIIS_extrapolation.f90 b/src/utils/complex_DIIS_extrapolation.f90 index 33436b3..81192b0 100644 --- a/src/utils/complex_DIIS_extrapolation.f90 +++ b/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 diff --git a/src/utils/complex_gram_schmidt.f90 b/src/utils/complex_gram_schmidt.f90 deleted file mode 100644 index 7ce152e..0000000 --- a/src/utils/complex_gram_schmidt.f90 +++ /dev/null @@ -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 diff --git a/src/utils/complex_orthogonalization.f90 b/src/utils/complex_orthogonalization.f90 new file mode 100644 index 0000000..7112063 --- /dev/null +++ b/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 diff --git a/src/utils/complex_orthogonalize_matrix.f90 b/src/utils/complex_orthogonalize_matrix.f90 deleted file mode 100644 index 3fbc431..0000000 --- a/src/utils/complex_orthogonalize_matrix.f90 +++ /dev/null @@ -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 diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 4036eb4..f6704ca 100644 --- a/src/utils/wrap_lapack.f90 +++ b/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)