From 378bfd9356bf630d31a369e1ddedbada2c778624 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 19 May 2021 15:11:44 +0200 Subject: [PATCH] RMOM --- input/methods | 8 +- src/HF/MOM.f90 | 195 ----------------------------------------- src/HF/RMOM.f90 | 209 ++++++++++++++++++++++++++++++++++++++++++++ src/QuAcK/QuAcK.f90 | 14 ++- 4 files changed, 225 insertions(+), 201 deletions(-) delete mode 100644 src/HF/MOM.f90 create mode 100644 src/HF/RMOM.f90 diff --git a/input/methods b/input/methods index 71f9f1c..13ed09a 100644 --- a/input/methods +++ b/input/methods @@ -1,11 +1,11 @@ -# RHF UHF KS MOM - T F F F +# RHF UHF KS MOM + T F F T # MP2* MP3 MP2-F12 F F F # CCD DCD CCSD CCSD(T) - T T T T + F F F F # drCCD rCCD lCCD pCCD - F F F T + F F F F # CIS* CIS(D) CID CISD FCI F F F F F # RPA* RPAx* ppRPA diff --git a/src/HF/MOM.f90 b/src/HF/MOM.f90 deleted file mode 100644 index 3b720e4..0000000 --- a/src/HF/MOM.f90 +++ /dev/null @@ -1,195 +0,0 @@ -subroutine MOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P) - -! Maximum overlap method - - implicit none - -! Input variables - - integer,intent(in) :: maxSCF,max_diis - double precision,intent(in) :: thresh - - integer,intent(in) :: nBas,nO - double precision,intent(in) :: ENuc - double precision,intent(in) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),X(nBas,nBas) - -! Local variables - - integer :: iBas,jBas - integer :: nSCF,nBasSq,n_diis - double precision :: ET,EV,EJ,EK,Conv,Gap - double precision :: rcond - double precision,external :: trace_matrix - double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:) - double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:) - double precision,allocatable :: cG(:,:),ON(:) - -! Output variables - - double precision,intent(inout):: ERHF,c(nBas,nBas),e(nBas),P(nBas,nBas) - -! Hello world - - write(*,*) - write(*,*)'************************************************' - write(*,*)'| Maximum overlap method |' - write(*,*)'************************************************' - write(*,*) - -! Useful quantities - - nBasSq = nBas*nBas - -! Memory allocation - - allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), & - cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas), & - cG(nBas,nBas),ON(nBas), & - error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) - -! Set up guess orbitals - - cG(:,:) = c(:,:) - -! Set up occupation numbers - - ON(1:nO) = 1d0 - ON(nO+1:nBas) = 0d0 - -! HOMO-LUMO transition - - ON(nO) = 0d0 - ON(nO+1) = 1d0 - - write(*,*) - write(*,*) ' --- Initial MO occupations --- ' - write(*,*) - call matout(nBas,1,ON) - write(*,*) - -! Compute density matrix - - call density_matrix(nBas,ON,c,P) - -! Initialization - - n_diis = 0 - F_diis(:,:) = 0d0 - error_diis(:,:) = 0d0 - Conv = 1d0 - nSCF = 0 - -!------------------------------------------------------------------------ -! Main SCF loop -!------------------------------------------------------------------------ - write(*,*) - write(*,*)'----------------------------------------------------' - write(*,*)'| MOM 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(*,*)'----------------------------------------------------' - - do while(Conv > thresh .and. nSCF < maxSCF) - -! Increment - - nSCF = nSCF + 1 - -! Build Fock matrix - - call Coulomb_matrix_AO_basis(nBas,P,ERI,J) - call exchange_matrix_AO_basis(nBas,P,ERI,K) - - F(:,:) = Hc(:,:) + J(:,:) + K(:,:) - -! Check convergence - - error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) - Conv = maxval(abs(error)) - -! DIIS extrapolation - - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) - -! Reset DIIS if required - - if(abs(rcond) < 1d-15) n_diis = 0 - -! Diagonalize Fock matrix - - Fp = matmul(transpose(X),matmul(F,X)) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,e) - c = matmul(X,cp) - -! MOM overlap - - call MOM_overlap(nBas,nO,S,cG,c,ON) - -! Density matrix - - call density_matrix(nBas,ON,c,P) - -! Compute HF energy - - ERHF = trace_matrix(nBas,matmul(P,Hc)) & - + 0.5d0*trace_matrix(nBas,matmul(P,J)) & - + 0.5d0*trace_matrix(nBas,matmul(P,K)) - -! Compute HOMO-LUMO gap - - if(nBas > nO) then - - Gap = e(nO+1) - e(nO) - - else - - Gap = 0d0 - - endif - -! 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,'|' - - enddo - write(*,*)'----------------------------------------------------' -!------------------------------------------------------------------------ -! End of SCF loop -!------------------------------------------------------------------------ - -! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) - - stop - - endif - - write(*,*) - write(*,*) ' --- Final MO occupations --- ' - write(*,*) - call matout(nBas,1,ON) - write(*,*) - -! 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.5d0*trace_matrix(nBas,matmul(P,K)) - ERHF = ET + EV + EJ + EK - - call print_RHF(nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF) - -end subroutine MOM diff --git a/src/HF/RMOM.f90 b/src/HF/RMOM.f90 new file mode 100644 index 0000000..c086b36 --- /dev/null +++ b/src/HF/RMOM.f90 @@ -0,0 +1,209 @@ +subroutine RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P,Vx) + +! Perform restricted Hartree-Fock calculation with MOM algorithm + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF,max_diis,guess_type + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nO + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: nSCF + integer :: nBasSq + integer :: n_diis + double precision :: ET + double precision :: EV + double precision :: EJ + double precision :: EK + double precision :: dipole(ncart) + + double precision :: Conv + double precision :: Gap + double precision :: rcond + double precision,external :: trace_matrix + double precision,allocatable :: error(:,:) + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: F_diis(:,:) + double precision,allocatable :: J(:,:) + double precision,allocatable :: K(:,:) + double precision,allocatable :: cp(:,:) + double precision,allocatable :: F(:,:) + double precision,allocatable :: Fp(:,:) + double precision,allocatable :: ON(:) + +! Output variables + + double precision,intent(out) :: ERHF + double precision,intent(out) :: e(nBas) + double precision,intent(out) :: c(nBas,nBas) + double precision,intent(out) :: P(nBas,nBas) + double precision,intent(out) :: Vx(nBas) + +! Hello world + + write(*,*) + write(*,*)'*********************************************' + write(*,*)'| Restricted Maximum Overlap Method |' + write(*,*)'*********************************************' + write(*,*) + +! Useful quantities + + nBasSq = nBas*nBas + +! Memory allocation + + allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), & + cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas),ON(nBas), & + error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + +! Guess coefficients and eigenvalues + + call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) + +! ON(:) = 0d0 +! do i=1,nO +! ON(i) = 1d0 +! ON(i) = dble(2*i-1) +! end do + +! call density_matrix(nBas,ON,c,P) + +! Initialization + + n_diis = 0 + F_diis(:,:) = 0d0 + error_diis(:,:) = 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(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Build Fock matrix + + call Coulomb_matrix_AO_basis(nBas,P,ERI,J) + call exchange_matrix_AO_basis(nBas,P,ERI,K) + + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + +! Check convergence + + error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + Conv = maxval(abs(error)) + +! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + if(abs(rcond) > 1d-7) then + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + else + n_diis = 0 + end if + +! Diagonalize Fock matrix + + Fp = matmul(transpose(X),matmul(F,X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas,cp,e) + c = matmul(X,cp) + +! Density matrix + + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + +! call density_matrix(nBas,ON,c,P) + +! 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 = e(nO+1) - e(nO) + + else + + Gap = 0d0 + + endif + +! 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,'|' + + enddo + write(*,*)'----------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +! 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) + call print_RHF(nBas,nO,e,C,ENuc,ET,EV,EJ,EK,ERHF,dipole) + +! Compute Vx for post-HF calculations + + call exchange_potential(nBas,c,K,Vx) + +end subroutine RMOM diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 8bbff6e..5ed4a55 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -343,8 +343,18 @@ program QuAcK if(doMOM) then call cpu_time(start_MOM) - call MOM(maxSCF_HF,thresh_HF,n_diis_HF, & - nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,ERHF,cHF,eHF,PHF) + + if(unrestricted) then + +! call UMOM() + + else + + call RMOM(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & + nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc) + + end if + call cpu_time(end_MOM) t_MOM = end_MOM - start_MOM