From ddb2ccfb54b047013407b7bad05b4665b2d46c8b Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 9 Nov 2023 23:23:00 +0100 Subject: [PATCH] fix ROHF --- input/methods | 4 +-- src/HF/ROHF.f90 | 18 ++++++------- src/HF/ROHF_fock_matrix.f90 | 6 ++--- src/HF/print_GHF.f90 | 50 ++++++++++++++++++------------------- src/HF/print_RHF.f90 | 4 +-- src/HF/print_ROHF.f90 | 4 +-- src/QuAcK/QuAcK.f90 | 2 +- 7 files changed, 44 insertions(+), 44 deletions(-) diff --git a/input/methods b/input/methods index 97c42d5..24d535e 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF GHF ROHF - F T F F + F F F T # MP2 MP3 F F # CCD pCCD DCD CCSD CCSD(T) @@ -13,6 +13,6 @@ # G0F2 evGF2 qsGF2 G0F3 evGF3 F F F F F # G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW - F F T F F F + F F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh F F F F F F diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index 03d6c1a..bee532a 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -35,7 +35,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc integer :: nSCF integer :: nBasSq integer :: n_diis - double precision :: conv + double precision :: Conv double precision :: rcond double precision :: ET(nspin) double precision :: EV(nspin) @@ -67,9 +67,9 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc ! Hello world write(*,*) - write(*,*)'************************************************' - write(*,*)'* Restricted Open-Shell Hartree-Fock *' - write(*,*)'************************************************' + write(*,*)'****************************************' + write(*,*)'* Restricted Open-Shell HF Calculation *' + write(*,*)'****************************************' write(*,*) ! Useful stuff @@ -91,13 +91,13 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc ! Initialization - n_diis = 0 + n_diis = 0 F_diis(:,:) = 0d0 err_diis(:,:) = 0d0 rcond = 0d0 nSCF = 0 - conv = 1d0 + Conv = 1d0 !------------------------------------------------------------------------ ! Main SCF loop @@ -109,7 +109,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc '|','#','|','E(ROHF)','|','Ex(ROHF)','|','Conv','|' write(*,*)'----------------------------------------------------------' - do while(conv > thresh .and. nSCF < maxSCF) + do while(Conv > thresh .and. nSCF < maxSCF) ! Increment @@ -138,7 +138,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc ! Check convergence err(:,:) = matmul(Ftot(:,:),matmul(Ptot(:,:),S(:,:))) - matmul(matmul(S(:,:),Ptot(:,:)),Ftot(:,:)) - if(nSCF > 1) conv = maxval(abs(err(:,:))) + if(nSCF > 1) Conv = maxval(abs(err(:,:))) ! DIIS extrapolation @@ -214,7 +214,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc ! 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,'|' + '|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',Conv,'|' end do write(*,*)'----------------------------------------------------------' diff --git a/src/HF/ROHF_fock_matrix.f90 b/src/HF/ROHF_fock_matrix.f90 index ef2c4c7..9272ef3 100644 --- a/src/HF/ROHF_fock_matrix.f90 +++ b/src/HF/ROHF_fock_matrix.f90 @@ -74,8 +74,8 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,Fa,Fb,F) F(nO+nC+1:nC+nO+nV, nC+1:nC+nO ) = Fa(nO+nC+1:nC+nO+nV, nC+1:nC+nO ) F(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) = aV*Fa(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) + bV*Fb(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) - call MOtoAO_transform(nBas,S,c,F) - call MOtoAO_transform(nBas,S,c,Fa) - call MOtoAO_transform(nBas,S,c,Fb) + call MOtoAO_transform(nBas,S,c,F,F) + call MOtoAO_transform(nBas,S,c,Fa,Fa) + call MOtoAO_transform(nBas,S,c,Fb,Fb) end subroutine diff --git a/src/HF/print_GHF.f90 b/src/HF/print_GHF.f90 index d296f39..8fde002 100644 --- a/src/HF/print_GHF.f90 +++ b/src/HF/print_GHF.f90 @@ -90,44 +90,44 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole) ! 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: ',EHF,' au' - write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' - write(*,'(A32,1X,F16.10,A3)') ' GHF energy: ',EHF + ENuc,' au' - write(*,'(A50)') '-----------------------------------------' - write(*,'(A32,1X,F16.6,A3)') ' GHF HOMO energy: ',e(HOMO)*HaToeV,' eV' - write(*,'(A32,1X,F16.6,A3)') ' GHF LUMO energy: ',e(LUMO)*HaToeV,' eV' - write(*,'(A32,1X,F16.6,A3)') ' GHF HOMO-LUMO gap : ',Gap*HaToeV,' eV' - write(*,'(A50)') '-----------------------------------------' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33)') ' Summary ' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',ET + EV,' au' + write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',ET,' au' + write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',EV,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',EJ + EK,' au' + write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',EJ,' au' + write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',EK,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',EHF,' au' + write(*,'(A33,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au' + write(*,'(A33,1X,F16.10,A3)') ' GHF energy = ',EHF + ENuc,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.6,A3)') ' GHF HOMO energy = ',e(HOMO)*HaToeV,' eV' + write(*,'(A33,1X,F16.6,A3)') ' GHF LUMO energy = ',e(LUMO)*HaToeV,' eV' + write(*,'(A33,1X,F16.6,A3)') ' GHF HOMO-LUMO gap = ',Gap*HaToeV,' eV' + write(*,'(A50)') '---------------------------------------' ! write(*,'(A32,1X,F16.6)') ' :',S2 -! write(*,'(A50)') '-----------------------------------------' - write(*,'(A35)') ' Dipole moment (Debye) ' +! write(*,'(A50)') '---------------------------------------' + write(*,'(A36)') ' 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(*,'(10X,4F10.4)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD + write(*,'(A50)') '---------------------------------------' write(*,*) ! Print results if(dump_orb) then write(*,'(A50)') '---------------------------------------' - write(*,'(A32)') ' GHF orbital coefficients' + write(*,'(A50)') ' GHF orbital coefficients ' write(*,'(A50)') '---------------------------------------' call matout(nBas2,nBas2,c) write(*,*) end if write(*,'(A50)') '---------------------------------------' - write(*,'(A32)') ' GHF orbital energies' + write(*,'(A50)') ' GHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' call matout(nBas2,1,e) write(*,*) diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 51ff774..ce14458 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -73,7 +73,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) if(dump_orb) then write(*,'(A50)') '---------------------------------------' - write(*,'(A50)') ' RHF orbital coefficients' + write(*,'(A50)') ' RHF orbital coefficients ' write(*,'(A50)') '---------------------------------------' call matout(nBas,nBas,cHF) write(*,*) @@ -81,7 +81,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' RHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call matout(nBas,1,eHF) + call vecout(nBas,eHF) write(*,*) end subroutine diff --git a/src/HF/print_ROHF.f90 b/src/HF/print_ROHF.f90 index acfc694..e818c3d 100644 --- a/src/HF/print_ROHF.f90 +++ b/src/HF/print_ROHF.f90 @@ -90,7 +90,7 @@ subroutine print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) 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(*,'(19X,4F10.4)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD write(*,'(A60)') '-------------------------------------------------' write(*,*) @@ -102,7 +102,7 @@ subroutine print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) call matout(nBas,nBas,c) write(*,*) write(*,'(A50)') '---------------------------------------' - write(*,'(A50)') ' ROHF orbital energies ' + write(*,'(A50)') ' ROHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' call matout(nBas,1,e) write(*,*) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 81b179d..543d206 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -172,7 +172,7 @@ program QuAcK !---------------------! doRQuAcK = .false. - if(doRHF) doRQuAcK = .true. + if(doRHF .or. doROHF) doRQuAcK = .true. doUQuAcK = .false. if(doUHF) doUQuAcK = .true.