diff --git a/input/methods b/input/methods index b71abbf..f582e3d 100644 --- a/input/methods +++ b/input/methods @@ -1,19 +1,18 @@ # RHF UHF GHF ROHF - F T F F -# MP2* MP3 + F F T F +# MP2 MP3 F F # CCD pCCD DCD CCSD CCSD(T) F F F F F # drCCD rCCD crCCD lCCD F F F F -# CIS* CIS(D) CID CISD FCI - F F F F F -# phRPA* phRPAx* crRPA ppRPA - F F F F -# G0F2* evGF2* qsGF2* G0F3 evGF3 - F F F F F -# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW - T F F F F F -# G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh - F F F F F F -# * unrestricted version available +# CIS CIS(D) CID CISD FCI + F F F F F +# phRPA phRPAx crRPA ppRPA + F F F F +# G0F2 evGF2 qsGF2 G0F3 evGF3 + F F F F F +# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW + F F F F F F +# G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh + F F F F F F diff --git a/input/options b/input/options index e86f9b7..8a5de02 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS guess mix_guess level_shift stability - 10000 0.000001 5 2 0.0 0.0 T + 1000 0.0000001 5 3 0.0 0.0 T # MP: reg F # CC: maxSCF thresh DIIS diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90 index 8d8f38f..93ac208 100644 --- a/src/HF/GHF.f90 +++ b/src/HF/GHF.f90 @@ -42,7 +42,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, double precision :: ET,ETaa,ETbb double precision :: EV,EVaa,EVbb double precision :: EJ,EJaaaa,EJaabb,EJbbaa,EJbbbb - double precision :: Ex,Exaaaa,Exabba,Exbaab,Exbbbb + double precision :: EK,EKaaaa,EKabba,EKbaab,EKbbbb double precision :: dipole(ncart) double precision,allocatable :: Jaa(:,:),Jbb(:,:) @@ -139,7 +139,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, 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(GHF)','|','EJ(GHF)','|','Ex(GHF)','|','Conv','|' + '|','#','|','E(GHF)','|','EJ(GHF)','|','EK(GHF)','|','Conv','|' write(*,*)'-----------------------------------------------------------------------------' do while(Conv > thresh .and. nSCF < maxSCF) @@ -239,7 +239,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, EV = EVaa + EVbb -! Hartree energy: 16 terms? +! Hartree energy EJaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Jaa)) EJaabb = 0.5d0*trace_matrix(nBas,matmul(Paa,Jbb)) @@ -250,21 +250,21 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Exchange energy - Exaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Kaa)) - Exabba = 0.5d0*trace_matrix(nBas,matmul(Pab,Kba)) - Exbaab = 0.5d0*trace_matrix(nBas,matmul(Pba,Kab)) - Exbbbb = 0.5d0*trace_matrix(nBas,matmul(Pbb,Kbb)) + EKaaaa = 0.5d0*trace_matrix(nBas,matmul(Paa,Kaa)) + EKabba = 0.5d0*trace_matrix(nBas,matmul(Pab,Kba)) + EKbaab = 0.5d0*trace_matrix(nBas,matmul(Pba,Kab)) + EKbbbb = 0.5d0*trace_matrix(nBas,matmul(Pbb,Kbb)) - Ex = Exaaaa + Exabba + Exbaab + Exbbbb + EK = EKaaaa + EKabba + EKbaab + EKbbbb ! Total energy - EHF = ET + EV + EJ + Ex + EHF = ET + EV + EJ + EK ! 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)') & - '|',nSCF,'|',EHF + ENuc,'|',EJ,'|',Ex,'|',Conv,'|' + '|',nSCF,'|',EHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|' end do write(*,*)'-----------------------------------------------------------------------------' @@ -292,6 +292,6 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Compute final GHF energy -! call print_GHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) + call print_GHF(nBas2,nO,e,C,ENuc,ET,EV,EJ,EK,EHF,dipole) end subroutine diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 new file mode 100644 index 0000000..8ea8fa7 --- /dev/null +++ b/src/HF/print_GHF.f90 @@ -0,0 +1,77 @@ +subroutine print_GHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,EGHF,dipole) + +! Print one-electron energies and other stuff for GHF + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nO + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: ENuc + double precision,intent(in) :: ET + double precision,intent(in) :: EV + double precision,intent(in) :: EJ + double precision,intent(in) :: EK + double precision,intent(in) :: EGHF + double precision,intent(in) :: dipole(ncart) + +! Local variables + + integer :: ixyz + integer :: HOMO + integer :: LUMO + double precision :: Gap + +! HOMO and LUMO + + HOMO = nO + LUMO = HOMO + 1 + Gap = eHF(LUMO)-eHF(HOMO) + +! Dump results + + + write(*,*) + write(*,'(A50)') '-----------------------------------------' + write(*,'(A32)') ' Summary ' + write(*,'(A50)') '-----------------------------------------' + write(*,'(A32,1X,F16.10,A3)') ' One-electron energy: ',ET + EV,' au' + write(*,'(A32,1X,F16.10,A3)') ' Kinetic energy: ',ET,' au' + write(*,'(A32,1X,F16.10,A3)') ' Potential energy: ',EV,' au' + write(*,'(A50)') '-----------------------------------------' + write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + EK,' au' + write(*,'(A32,1X,F16.10,A3)') ' Hartree energy: ',EJ,' au' + write(*,'(A32,1X,F16.10,A3)') ' Exchange energy: ',EK,' au' + write(*,'(A50)') '-----------------------------------------' + write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EGHF,' au' + write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' + write(*,'(A32,1X,F16.10,A3)') ' GHF energy: ',EGHF + ENuc,' au' + write(*,'(A50)') '-----------------------------------------' + write(*,'(A32,1X,F16.6,A3)') ' HF HOMO energy: ',eHF(HOMO)*HaToeV,' eV' + write(*,'(A32,1X,F16.6,A3)') ' HF LUMO energy: ',eHF(LUMO)*HaToeV,' eV' + write(*,'(A32,1X,F16.6,A3)') ' HF HOMO-LUMO gap : ',Gap*HaToeV,' eV' + write(*,'(A50)') '-----------------------------------------' + write(*,'(A35)') ' Dipole moment (Debye) ' + write(*,'(10X,4A10)') 'X','Y','Z','Tot.' + write(*,'(10X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD + write(*,'(A50)') '-----------------------------------------' + write(*,*) + +! Print results + + write(*,'(A50)') '---------------------------------------' + write(*,'(A32)') 'MO coefficients' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,nBas,cHF) + write(*,*) + write(*,'(A50)') '---------------------------------------' + write(*,'(A32)') 'MO energies' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,1,eHF) + write(*,*) + +end subroutine