From 247e887a5134d1e7e90c971865e0616fbe9ce7aa Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 16 Nov 2023 17:56:50 +0100 Subject: [PATCH 1/5] Update README.md --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index eec134f..dc33520 100644 --- a/README.md +++ b/README.md @@ -6,6 +6,7 @@ - [Pierre-Francois Loos](https://pfloos.github.io/WEB_LOOS) - [Enzo Monino](https://enzomonino.github.io) - [Antoine Marie](https://antoine-marie.github.io) +- [Abdallah Amar](https://scholar.google.com/citations?user=y437T5sAAAAJ&hl=en) - [Anthony Scemama](https://scemama.github.io) # What is it? From 784ccf9632eaff9bc6f931f774a24d995c924499 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 16 Nov 2023 17:57:05 +0100 Subject: [PATCH 2/5] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index dc33520..de5cb74 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ - [Pierre-Francois Loos](https://pfloos.github.io/WEB_LOOS) - [Enzo Monino](https://enzomonino.github.io) - [Antoine Marie](https://antoine-marie.github.io) -- [Abdallah Amar](https://scholar.google.com/citations?user=y437T5sAAAAJ&hl=en) +- [Abdallah Ammar](https://scholar.google.com/citations?user=y437T5sAAAAJ&hl=en) - [Anthony Scemama](https://scemama.github.io) # What is it? From 0aabc1a5f4dd5710aca51ae9be49c4a7e31bd146 Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 16 Nov 2023 19:38:32 +0100 Subject: [PATCH 3/5] clean up S^2 GHF --- src/HF/print_GHF.f90 | 57 +++++++------------------------------------- 1 file changed, 9 insertions(+), 48 deletions(-) diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index 4bfeb1a..d9d1e69 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -80,7 +80,6 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) Pba = matmul(transpose(Cb),matmul(S,Ca)) Pbb = matmul(transpose(Cb),matmul(S,Cb)) - ! Compute components of S = (Sx,Sy,Sz) Sx = 0.5d0*(trace_matrix(nO,Pab) + trace_matrix(nO,Pba)) @@ -89,17 +88,14 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) ! Compute = + + - SpSm = 0d0 do i=1,nO do j=1,nO SpSm = SpSm + Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i) end do end do - SpSm = 0.5d0*(trace_matrix(nO,Paa) + SpSm) - -! Sx2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Pab+Pba)**2 & -! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) + matmul(Pab,Pab)) + SpSm = trace_matrix(nO,Paa) + SpSm + print*,' = ',SpSm SmSp = 0d0 do i=1,nO @@ -107,10 +103,8 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) SmSp = SmSp + Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i) end do end do - SmSp = 0.5d0*(trace_matrix(nO,Pbb) + SmSp) - -! Sy2 = 0.25d0*trace_matrix(nO,Paa+Pbb) - 0.25d0*trace_matrix(nO,Pab-Pba)**2 & -! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) - matmul(Pab,Pab)) + SmSp = trace_matrix(nO,Pbb) + SmSp + print*,' = ',SmSp Sz2 = 0d0 do i=1,nO @@ -121,45 +115,17 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) Sz2 = 0.25d0*(dble(nO) + Sz2) print*,' = ',Sz2 - Sz2 = 0.25d0*trace_matrix(nO,Paa+Pbb) & - - 0.25d0*trace_matrix(nO,Paa-Pbb)**2 & - - 0.25d0*trace_matrix(nO,matmul(Paa-Pbb,Paa-Pbb)) +! Compute from Sz^2, S^+S^- and S^-S^+ - S2 = Sz*(Sz+1d0) + trace_matrix(nO,Pbb) + 0.25d0*trace_matrix(nO,Paa+Pbb) - - do i=1,nO - do j=1,nO - S2 = S2 - 0.25d0*(Paa(i,j) - Pbb(i,j))**2 & - + (Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i)) - end do - end do + S2 = Sz2 + 0.5d0*(SpSm + SmSp) print*,' = ',S2 +! Compute and from , , and + Sx2 = 0.5d0*(S2 - Sz2 + 0.5d0*(SmSp + SpSm)) print*,' = ',Sx2 Sy2 = 0.5d0*(S2 - Sz2 - 0.5d0*(SmSp + SpSm)) print*,' = ',Sy2 -! Sx2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Pab+Pba)**2 & -! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) + matmul(Pab,Pab)) - -! Sy2 = 0.25d0*trace_matrix(nO,Paa+Pbb) - 0.25d0*trace_matrix(nO,Pab-Pba)**2 & -! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) - matmul(Pab,Pab)) - -! Sz2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Paa-Pbb)**2 & -! - 0.25d0*trace_matrix(nO,matmul(Paa,Paa) + matmul(Pbb,Pbb)) & -! + 0.25d0*trace_matrix(nO,matmul(Pab,Pba) + matmul(Pba,Pab)) - -! S2 = Sz*(Sz+1d0) + trace_matrix(nO,Pbb) + 0.25d0*trace_matrix(nO,Paa+Pbb) - -! do i=1,nO -! do j=1,nO -! S2 = S2 - 0.25d0*(Paa(i,j) - Pbb(i,j))**2 & -! + (Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i)) -! end do -! end do -! print*,' = ',S2 - - ! na = 0.d0 ! nb = 0.d0 @@ -187,9 +153,6 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) ! enddo ! S2 = Sz * (Sz + 1.d0) + nonco_z + contam_ghf - - - ! deallocate(Paa,Pab,Pba,Pbb) ! Check collinearity and coplanarity @@ -255,10 +218,8 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) write(*,'(A33,1X,F16.6)') ' = ',Sy write(*,'(A33,1X,F16.6)') ' = ',Sz write(*,'(A50)') '---------------------------------------' - write(*,'(A33,1X,F16.6)') ' = ',Sx2 - write(*,'(A33,1X,F16.6)') ' = ',Sy2 + write(*,'(A33,1X,F16.6)') ' = ',S2 - Sz2 write(*,'(A33,1X,F16.6)') ' = ',Sz2 - write(*,'(A33,1X,F16.6)') ' = ',Sx2+Sy2+Sz2 write(*,'(A33,1X,F16.6)') ' = ',S2 write(*,'(A50)') '---------------------------------------' write(*,'(A36)') ' Dipole moment (Debye) ' From 3e8b3ab6da82f596da1b0ea68c13072cbe8a81fb Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 16 Nov 2023 19:45:02 +0100 Subject: [PATCH 4/5] more clean up in GHF print --- src/HF/print_GHF.f90 | 20 ++++---------------- 1 file changed, 4 insertions(+), 16 deletions(-) diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index d9d1e69..fe53865 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -33,11 +33,10 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) integer :: HOMO integer :: LUMO double precision :: Gap - double precision :: Sx ,Sy ,Sz - double precision :: Sx2,Sy2,Sz2 - double precision :: SmSp,SpSm,S2 - double precision :: na, nb - double precision :: nonco_z, contam_uhf, xy_perp, contam_ghf + double precision :: Sx,Sy,Sz + double precision :: SmSp,SpSm,Sz2,S2 +! double precision :: na, nb +! double precision :: nonco_z, contam_uhf, xy_perp, contam_ghf double precision,allocatable :: Ca(:,:) double precision,allocatable :: Cb(:,:) @@ -95,7 +94,6 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) end do end do SpSm = trace_matrix(nO,Paa) + SpSm - print*,' = ',SpSm SmSp = 0d0 do i=1,nO @@ -104,7 +102,6 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) end do end do SmSp = trace_matrix(nO,Pbb) + SmSp - print*,' = ',SmSp Sz2 = 0d0 do i=1,nO @@ -113,19 +110,10 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole) end do end do Sz2 = 0.25d0*(dble(nO) + Sz2) - print*,' = ',Sz2 ! Compute from Sz^2, S^+S^- and S^-S^+ S2 = Sz2 + 0.5d0*(SpSm + SmSp) - print*,' = ',S2 - -! Compute and from , , and - - Sx2 = 0.5d0*(S2 - Sz2 + 0.5d0*(SmSp + SpSm)) - print*,' = ',Sx2 - Sy2 = 0.5d0*(S2 - Sz2 - 0.5d0*(SmSp + SpSm)) - print*,' = ',Sy2 ! na = 0.d0 ! nb = 0.d0 From d5e7be1500e976485111b8bdfaad2f9b1c2ddb4f Mon Sep 17 00:00:00 2001 From: pfloos Date: Fri, 17 Nov 2023 00:26:44 +0100 Subject: [PATCH 5/5] fix scf convergence issues in RHF, UHF, ROHF and GHF --- src/HF/GHF.f90 | 88 +++++++++++++--------------- src/HF/GHF_search.f90 | 2 + src/HF/RHF.f90 | 72 +++++++++++------------ src/HF/ROHF.f90 | 99 ++++++++++++++++---------------- src/HF/UHF.f90 | 88 ++++++++++++++-------------- src/HF/print_ROHF.f90 | 14 +++-- src/utils/DIIS_extrapolation.f90 | 9 ++- src/utils/wrap_lapack.f90 | 11 +++- 8 files changed, 190 insertions(+), 193 deletions(-) diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90 index c3ffd00..56b8058 100644 --- a/src/HF/GHF.f90 +++ b/src/HF/GHF.f90 @@ -182,52 +182,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu if(nSCF > 1) Conv = maxval(abs(err)) -! DIIS extrapolation - - if(max_diis > 1) then - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBas2Sq,nBas2Sq,n_diis,err_diis(:,1:n_diis),F_diis(:,1:n_diis),err,F) - end if - -! Level-shifting - - if(level_shift > 0d0 .and. Conv > thresh) then - call level_shifting(level_shift,nBas,nO,S,C,F) - end if - -! Transform Fock matrix in orthogonal basis - - Fp(:,:) = matmul(transpose(X),matmul(F,X)) - -! Diagonalize Fock matrix to get eigenvectors and eigenvalues - - Cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas2,Cp,eHF) - -! Back-transform eigenvectors in non-orthogonal basis - - C(:,:) = matmul(X,Cp) - -! Mix guess for UHF solution in singlet states - -! if(nSCF == 1) call mix_guess(nBas,nO,mix,c) - -! Form super density matrix - - P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO))) - -! Compute individual density matrices - - Paa(:,:) = P( 1:nBas , 1:nBas ) - Pab(:,:) = P( 1:nBas ,nBas+1:nBas2) - Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) - Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) - -!------------------------------------------------------------------------ -! Compute UHF energy -!------------------------------------------------------------------------ - -! Kinetic energy +! Kinetic energy ETaa = trace_matrix(nBas,matmul(Paa,T)) ETbb = trace_matrix(nBas,matmul(Pbb,T)) @@ -263,9 +218,48 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu EGHF = ET + EV + EJ + EK +! DIIS extrapolation + + if(max_diis > 1) then + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nBas2Sq,nBas2Sq,n_diis,err_diis,F_diis,err,F) + end if + +! Level-shifting + + if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,C,F) + +! Transform Fock matrix in orthogonal basis + + Fp(:,:) = matmul(transpose(X),matmul(F,X)) + +! Diagonalize Fock matrix to get eigenvectors and eigenvalues + + Cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas2,Cp,eHF) + +! Back-transform eigenvectors in non-orthogonal basis + + C(:,:) = matmul(X,Cp) + +! Mix guess for UHF solution in singlet states + +! if(nSCF == 1) call mix_guess(nBas,nO,mix,c) + +! Form super density matrix + + P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO))) + +! Compute individual density matrices + + Paa(:,:) = P( 1:nBas , 1:nBas ) + Pab(:,:) = P( 1:nBas ,nBas+1:nBas2) + Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) + Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) + ! Dump results - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') & '|',nSCF,'|',EGHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|' end do diff --git a/src/HF/GHF_search.f90 b/src/HF/GHF_search.f90 index 5d43822..992d6e8 100644 --- a/src/HF/GHF_search.f90 +++ b/src/HF/GHF_search.f90 @@ -211,6 +211,8 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu call matrix_exponential(nBas2,R,ExpR) c = matmul(c,ExpR) +! call matout(nBas2,nBas2,matmul(transpose(ExpR),ExpR)) + else write(*,'(1X,A40,1X)') 'Well done, GHF solution is stable!' diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 796f6c0..59b1ac0 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -42,7 +42,6 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN double precision :: dipole(ncart) double precision :: Conv - double precision :: Gap double precision :: rcond double precision,external :: trace_matrix double precision,allocatable :: error(:,:) @@ -82,26 +81,25 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN call mo_guess(nBas,guess_type,S,Hc,X,c) P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - + ! Initialization F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - Conv = 1d0 - n_diis = 0 - nSCF = 0 rcond = 0d0 + Conv = 1d0 + nSCF = 0 + !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ + write(*,*) - write(*,*)'----------------------------------------------------' - write(*,*)'| RHF calculation |' - write(*,*)'----------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','HF energy','|','Conv','|','HL Gap','|' - write(*,*)'----------------------------------------------------' + write(*,*)'-----------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(RHF)','|','EJ(RHF)','|','EK(RHF)','|','Conv','|' + write(*,*)'-----------------------------------------------------------------------------' do while(Conv > thresh .and. nSCF < maxSCF) @@ -121,6 +119,26 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) Conv = maxval(abs(error)) +! Kinetic energy + + ET = trace_matrix(nBas,matmul(P,T)) + +! Potential energy + + EV = trace_matrix(nBas,matmul(P,V)) + +! Hartree energy + + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + +! Exchange energy + + EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) + +! Total energy + + ERHF = ET + EV + EJ + EK + ! DIIS extrapolation if(max_diis > 1) then @@ -144,32 +162,14 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! Density matrix P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - -! Compute HF energy - - ERHF = trace_matrix(nBas,matmul(P,Hc)) & - + 0.5d0*trace_matrix(nBas,matmul(P,J)) & - + 0.25d0*trace_matrix(nBas,matmul(P,K)) - -! Compute HOMO-LUMO gap - - if(nBas > nO) then - - Gap = eHF(nO+1) - eHF(nO) - - else - - Gap = 0d0 - - end if ! Dump results - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',ERHF+ENuc,'|',Conv,'|',Gap,'|' + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') & + '|',nSCF,'|',ERHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|' end do - write(*,*)'----------------------------------------------------' + write(*,*)'-----------------------------------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop !------------------------------------------------------------------------ @@ -188,14 +188,6 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN end if -! Compute HF energy - - ET = trace_matrix(nBas,matmul(P,T)) - EV = trace_matrix(nBas,matmul(P,V)) - EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) - EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) - ERHF = ET + EV + EJ + EK - ! Compute dipole moments call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index a84ef20..7d86359 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -1,5 +1,5 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,Ptot) + nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EROHF,eHF,c,Ptot) ! Perform restricted open-shell Hartree-Fock calculation @@ -42,7 +42,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN double precision :: ET(nspin) double precision :: EV(nspin) double precision :: EJ(nsp) - double precision :: Ex(nspin) + double precision :: EK(nspin) double precision :: dipole(ncart) double precision,allocatable :: cp(:,:) @@ -61,9 +61,9 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Output variables - double precision,intent(out) :: EHF - double precision,intent(out) :: e(nBas) - double precision,intent(out) :: c(nBas,nBas) + double precision,intent(out) :: EROHF + double precision,intent(out) :: eHF(nBas) + double precision,intent(inout):: c(nBas,nBas) double precision,intent(out) :: Ptot(nBas,nBas) ! Hello world @@ -80,8 +80,8 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Memory allocation - allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), & - P(nBas,nBas,nspin),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), & + allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), & + P(nBas,nBas,nspin),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), & err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Guess coefficients and demsity matrices @@ -90,6 +90,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN do ispin=1,nspin P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin)))) end do + Ptot(:,:) = P(:,:,1) + P(:,:,2) ! Initialization @@ -106,10 +107,10 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN !------------------------------------------------------------------------ write(*,*) - write(*,*)'----------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','E(ROHF)','|','Ex(ROHF)','|','Conv','|' - write(*,*)'----------------------------------------------------------' + write(*,*)'-----------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(ROHF)','|','EJ(ROHF)','|','EK(ROHF)','|','Conv','|' + write(*,*)'-----------------------------------------------------------------------------' do while(Conv > thresh .and. nSCF < maxSCF) @@ -139,15 +140,43 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Check convergence - err(:,:) = matmul(Ftot(:,:),matmul(Ptot(:,:),S(:,:))) - matmul(matmul(S(:,:),Ptot(:,:)),Ftot(:,:)) + err(:,:) = matmul(Ftot,matmul(Ptot,S)) - matmul(matmul(S,Ptot),Ftot) if(nSCF > 1) Conv = maxval(abs(err(:,:))) +! Kinetic energy + + do ispin=1,nspin + ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) + end do + +! Potential energy + + do ispin=1,nspin + EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) + end do + +! Hartree energy + + EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) + EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) + EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) + +! Exchange energy + + do ispin=1,nspin + EK(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin))) + end do + +! Total energy + + EROHF = sum(ET) + sum(EV) + sum(EJ) + sum(EK) + ! DIIS extrapolation if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis),F_diis(:,1:n_diis),err,Ftot) + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,Ftot) end if @@ -168,7 +197,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Diagonalize Fock matrix to get eigenvectors and eigenvalues cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,e) + call diagonalize_matrix(nBas,cp,eHF) ! Back-transform eigenvectors in non-orthogonal basis @@ -181,45 +210,13 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN end do Ptot(:,:) = P(:,:,1) + P(:,:,2) -!------------------------------------------------------------------------ -! Compute ROHF energy -!------------------------------------------------------------------------ - -! Kinetic energy - - do ispin=1,nspin - ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) - end do - -! Potential energy - - do ispin=1,nspin - EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) - end do - -! Hartree energy - - EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) - EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) - EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) - -! Exchange energy - - do ispin=1,nspin - Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin))) - end do - -! Total energy - - EHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) - ! Dump results - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',Conv,'|' + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') & + '|',nSCF,'|',EROHF + ENuc,'|',sum(EJ),'|',sum(EK),'|',Conv,'|' end do - write(*,*)'----------------------------------------------------------' + write(*,*)'-----------------------------------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop !------------------------------------------------------------------------ @@ -241,13 +238,13 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Compute final UHF energy call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) + call print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,EK,EROHF,dipole) ! Print test values if(dotest) then - call dump_test_value('R','ROHF energy',EHF) + call dump_test_value('R','ROHF energy',EROHF) end if diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 1e30d3f..1d1360d 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -42,7 +42,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu double precision :: ET(nspin) double precision :: EV(nspin) double precision :: EJ(nsp) - double precision :: Ex(nspin) + double precision :: EK(nspin) double precision :: dipole(ncart) double precision,allocatable :: cp(:,:,:) @@ -91,22 +91,22 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Initialization - nSCF = 0 - Conv = 1d0 - n_diis = 0 F_diis(:,:,:) = 0d0 err_diis(:,:,:) = 0d0 + nSCF = 0 + Conv = 1d0 + !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ write(*,*) - write(*,*)'----------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','E(UHF)','|','Ex(UHF)','|','Conv','|' - write(*,*)'----------------------------------------------------------' + write(*,*)'-----------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(UHF)','|','EJ(UHF)','|','EK(UHF)','|','Conv','|' + write(*,*)'-----------------------------------------------------------------------------' do while(Conv > thresh .and. nSCF < maxSCF) @@ -139,7 +139,35 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu end do if(nSCF > 1) Conv = maxval(abs(err(:,:,:))) - + +! Kinetic energy + + do ispin=1,nspin + ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) + end do + +! Potential energy + + do ispin=1,nspin + EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) + end do + +! Hartree energy + + EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) + EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) + EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) + +! Exchange energy + + do ispin=1,nspin + EK(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin))) + end do + +! Total energy + + EUHF = sum(ET) + sum(EV) + sum(EJ) + sum(EK) + ! DIIS extrapolation if(max_diis > 1) then @@ -174,7 +202,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu do ispin=1,nspin call diagonalize_matrix(nBas,cp(:,:,ispin),eHF(:,ispin)) end do - + ! Back-transform eigenvectors in non-orthogonal basis do ispin=1,nspin @@ -191,45 +219,13 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) end do -!------------------------------------------------------------------------ -! Compute UHF energy -!------------------------------------------------------------------------ - -! Kinetic energy - - do ispin=1,nspin - ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) - end do - -! Potential energy - - do ispin=1,nspin - EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) - end do - -! Hartree energy - - EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) - EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) - EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) - -! Exchange energy - - do ispin=1,nspin - Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin))) - end do - -! Total energy - - EUHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) - ! Dump results - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',EUHF + ENuc,'|',sum(Ex(:)),'|',Conv,'|' + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') & + '|',nSCF,'|',EUHF + ENuc,'|',sum(EJ),'|',sum(EK),'|',Conv,'|' end do - write(*,*)'----------------------------------------------------------' + write(*,*)'-----------------------------------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop !------------------------------------------------------------------------ @@ -251,7 +247,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Compute final UHF energy call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_UHF(nBas,nO,S,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) + call print_UHF(nBas,nO,S,eHF,c,ENuc,ET,EV,EJ,EK,EUHF,dipole) ! Print test values diff --git a/src/HF/print_ROHF.f90 b/src/HF/print_ROHF.f90 index db970ad..140bfa0 100644 --- a/src/HF/print_ROHF.f90 +++ b/src/HF/print_ROHF.f90 @@ -24,6 +24,8 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole) double precision :: Gap(nspin) double precision :: S,S2 + logical :: dump_orb = .false. + ! HOMO and LUMO do ispin=1,nspin @@ -96,11 +98,13 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole) ! Print results - write(*,'(A50)') '-----------------------------------------' - write(*,'(A50)') 'ROHF orbital coefficients ' - write(*,'(A50)') '-----------------------------------------' - call matout(nBas,nBas,c) - write(*,*) + if(dump_orb) then + write(*,'(A50)') '-----------------------------------------' + write(*,'(A50)') 'ROHF orbital coefficients ' + write(*,'(A50)') '-----------------------------------------' + call matout(nBas,nBas,c) + write(*,*) + end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' ROHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' diff --git a/src/utils/DIIS_extrapolation.f90 b/src/utils/DIIS_extrapolation.f90 index fdf8019..98b2dfc 100644 --- a/src/utils/DIIS_extrapolation.f90 +++ b/src/utils/DIIS_extrapolation.f90 @@ -9,11 +9,15 @@ subroutine DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_inout) ! Input variables integer,intent(in) :: n_err,n_e - double precision,intent(in) :: error_in(n_err),error(n_err,n_diis),e(n_e,n_diis) + double precision,intent(in) :: error_in(n_err) + double precision,intent(in) :: error(n_err,n_diis) + double precision,intent(in) :: e(n_e,n_diis) ! Local variables - double precision,allocatable :: A(:,:),b(:),w(:) + double precision,allocatable :: A(:,:) + double precision,allocatable :: b(:) + double precision,allocatable :: w(:) ! Output variables @@ -45,7 +49,6 @@ subroutine DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_inout) ! Solve linear system -! call easy_linear_solve(n_diis+1,A,b,w) call linear_solve(n_diis+1,A,b,w,rcond) ! Extrapolate diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 84b8332..c8f3797 100644 --- a/src/utils/wrap_lapack.f90 +++ b/src/utils/wrap_lapack.f90 @@ -175,7 +175,16 @@ subroutine linear_solve(N,A,b,x,rcond) integer,allocatable :: ipiv(:),iwork(:) double precision,allocatable :: AF(:,:),work(:) - lwork = 3*N + ! Find optimal size for temporary arrays + + allocate(work(1)) + + lwork = -1 + call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info) + lwork = int(work(1)) + + deallocate(work) + allocate(AF(N,N),ipiv(N),work(lwork),iwork(N)) call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info)