From 4e057b71f5b6ef82fb4e5b81cbc46f1e2f2e4975 Mon Sep 17 00:00:00 2001 From: pfloos Date: Wed, 6 Sep 2023 12:22:53 +0200 Subject: [PATCH] print ROHF --- input/methods | 2 +- src/HF/ROHF.f90 | 4 +- src/HF/print_ROHF.f90 | 117 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 120 insertions(+), 3 deletions(-) create mode 100644 src/HF/print_ROHF.f90 diff --git a/input/methods b/input/methods index 9843b39..4a9d502 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF ROHF RMOM UMOM KS - T F T F F F + F F T F F F # MP2* MP3 F F # CCD pCCD DCD CCSD CCSD(T) diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index b77308f..7c22a26 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -238,7 +238,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc ! Compute final UHF energy -! call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) -! call print_ROHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) + call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole) + call print_ROHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) end subroutine diff --git a/src/HF/print_ROHF.f90 b/src/HF/print_ROHF.f90 new file mode 100644 index 0000000..e6d165b --- /dev/null +++ b/src/HF/print_ROHF.f90 @@ -0,0 +1,117 @@ +subroutine print_ROHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) + +! Print one- and two-electron energies and other stuff for RoHF calculation + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + double precision,intent(in) :: Ov(nBas,nBas) + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: ENuc + double precision,intent(in) :: ET(nspin) + double precision,intent(in) :: EV(nspin) + double precision,intent(in) :: EJ(nsp) + double precision,intent(in) :: Ex(nspin) + double precision,intent(in) :: EHF + double precision,intent(in) :: dipole(ncart) + + integer :: ixyz + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) + double precision :: S_exact,S2_exact + double precision :: S,S2 + +! HOMO and LUMO + + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = e(nO(ispin)) + if(nO(ispin) < nBas) then + LUMO(ispin) = e(nO(ispin)+1) + else + LUMO(ispin) = 0d0 + end if + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = e(1) + Gap(ispin) = 0d0 + end if + end do + + S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0) + S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1))),matmul(Ov,c(:,1:nO(2))))**2) + + S_exact = 0.5d0*dble(nO(1) - nO(2)) + S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2) + +! Dump results + + write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40)') ' Summary ' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' One-electron energy: ',sum(ET(:)) + sum(EV(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' One-electron a energy: ',ET(1) + EV(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' One-electron b energy: ',ET(2) + EV(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic energy: ',sum(ET(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic a energy: ',ET(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic b energy: ',ET(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential energy: ',sum(EV(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential a energy: ',EV(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential b energy: ',EV(2),' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron energy: ',sum(EJ(:)) + sum(Ex(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree energy: ',sum(EJ(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree aa energy: ',EJ(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree ab energy: ',EJ(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Hartree bb energy: ',EJ(3),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',EHF,' au' + write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' + write(*,'(A40,1X,F16.10,A3)') ' UHF energy: ',EHF + ENuc,' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.6,A3)') ' UHF HOMO a energy:',HOMO(1)*HatoeV,' eV' + write(*,'(A40,1X,F16.6,A3)') ' UHF LUMO a energy:',LUMO(1)*HatoeV,' eV' + write(*,'(A40,1X,F16.6,A3)') ' UHF HOMOa-LUMOa gap:',Gap(1)*HatoeV,' eV' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.6,A3)') ' UHF HOMO b energy:',HOMO(2)*HatoeV,' eV' + write(*,'(A40,1X,F16.6,A3)') ' UHF LUMO b energy:',LUMO(2)*HatoeV,' eV' + write(*,'(A40,1X,F16.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.6)') ' S (exact) :',2d0*S_exact + 1d0 + write(*,'(A40,1X,F16.6)') ' S :',2d0*S + 1d0 + write(*,'(A40,1X,F16.6)') ' (exact) :',S2_exact + write(*,'(A40,1X,F16.6)') ' :',S2 + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A45)') ' Dipole moment (Debye) ' + write(*,'(19X,4A10)') 'X','Y','Z','Tot.' + write(*,'(19X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +! Print results + + write(*,'(A50)') '-----------------------------------------' + write(*,'(A50)') 'ROHF orbital coefficients ' + write(*,'(A50)') '-----------------------------------------' + call matout(nBas,nBas,c) + write(*,*) + write(*,'(A50)') '---------------------------------------' + write(*,'(A50)') ' ROHF orbital energies ' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,1,e) + write(*,*) + +end subroutine