PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-03 20:53:53 +01:00
Code Releases Activity

fix ROHF

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-11-09 23:23:00 +01:00
parent de2ace91fc
commit ddb2ccfb54
7 changed files with 44 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF GHF ROHF # RHF UHF GHF ROHF
F T F F F F F T
# MP2 MP3 # MP2 MP3
F F F F
# CCD pCCD DCD CCSD CCSD(T) # CCD pCCD DCD CCSD CCSD(T)
@ -13,6 +13,6 @@
# G0F2 evGF2 qsGF2 G0F3 evGF3 # G0F2 evGF2 qsGF2 G0F3 evGF3
F F F F F F F F F F
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW # 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 # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh
F F F F F F F F F F F F

18
src/HF/ROHF.f90
View File

@ -35,7 +35,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
integer :: nSCF integer :: nSCF
integer :: nBasSq integer :: nBasSq
integer :: n_diis integer :: n_diis
double precision :: conv double precision :: Conv
double precision :: rcond double precision :: rcond
double precision :: ET(nspin) double precision :: ET(nspin)
double precision :: EV(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 ! Hello world
write(*,*) write(*,*)
write(*,*)'************************************************' write(*,*)'****************************************'
write(*,*)'* Restricted Open-Shell Hartree-Fock *' write(*,*)'* Restricted Open-Shell HF Calculation *'
write(*,*)'************************************************' write(*,*)'****************************************'
write(*,*) write(*,*)
! Useful stuff ! Useful stuff
@ -91,13 +91,13 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
! Initialization ! Initialization
n_diis = 0 n_diis = 0
F_diis(:,:) = 0d0 F_diis(:,:) = 0d0
err_diis(:,:) = 0d0 err_diis(:,:) = 0d0
rcond = 0d0 rcond = 0d0
nSCF = 0 nSCF = 0
conv = 1d0 Conv = 1d0
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Main SCF loop ! 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','|' '|','#','|','E(ROHF)','|','Ex(ROHF)','|','Conv','|'
write(*,*)'----------------------------------------------------------' write(*,*)'----------------------------------------------------------'
do while(conv > thresh .and. nSCF < maxSCF) do while(Conv > thresh .and. nSCF < maxSCF)
! Increment ! Increment
@ -138,7 +138,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
! Check convergence ! 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(:,:))) if(nSCF > 1) Conv = maxval(abs(err(:,:)))
! DIIS extrapolation ! DIIS extrapolation
@ -214,7 +214,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
! Dump results ! 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)') & 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 end do
write(*,*)'----------------------------------------------------------' write(*,*)'----------------------------------------------------------'

6
src/HF/ROHF_fock_matrix.f90
View File

@ -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, 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) 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,F,F)
call MOtoAO_transform(nBas,S,c,Fa) call MOtoAO_transform(nBas,S,c,Fa,Fa)
call MOtoAO_transform(nBas,S,c,Fb) call MOtoAO_transform(nBas,S,c,Fb,Fb)
end subroutine end subroutine

50
src/HF/print_GHF.f90
View File

@ -90,44 +90,44 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
! Dump results ! Dump results
write(*,*) write(*,*)
write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' Summary ' write(*,'(A33)') ' Summary '
write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' One-electron energy: ',ET + EV,' au' write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',ET + EV,' au'
write(*,'(A32,1X,F16.10,A3)') ' Kinetic energy: ',ET,' au' write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',ET,' au'
write(*,'(A32,1X,F16.10,A3)') ' Potential energy: ',EV,' au' write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',EV,' au'
write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + EK,' au' write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',EJ + EK,' au'
write(*,'(A32,1X,F16.10,A3)') ' Hartree energy: ',EJ,' au' write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',EJ,' au'
write(*,'(A32,1X,F16.10,A3)') ' Exchange energy: ',EK,' au' write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',EK,' au'
write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EHF,' au' write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',EHF,' au'
write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' write(*,'(A33,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au'
write(*,'(A32,1X,F16.10,A3)') ' GHF energy: ',EHF + ENuc,' au' write(*,'(A33,1X,F16.10,A3)') ' GHF energy = ',EHF + ENuc,' au'
write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.6,A3)') ' GHF HOMO energy: ',e(HOMO)*HaToeV,' eV' write(*,'(A33,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(*,'(A33,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(*,'(A33,1X,F16.6,A3)') ' GHF HOMO-LUMO gap = ',Gap*HaToeV,' eV'
write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') '---------------------------------------'
! write(*,'(A32,1X,F16.6)') ' <S**2> :',S2 ! write(*,'(A32,1X,F16.6)') ' <S**2> :',S2
! write(*,'(A50)') '-----------------------------------------' ! write(*,'(A50)') '---------------------------------------'
write(*,'(A35)') ' Dipole moment (Debye) ' write(*,'(A36)') ' Dipole moment (Debye) '
write(*,'(10X,4A10)') 'X','Y','Z','Tot.' write(*,'(10X,4A10)') 'X','Y','Z','Tot.'
write(*,'(10X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD write(*,'(10X,4F10.4)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,*) write(*,*)
! Print results ! Print results
if(dump_orb) then if(dump_orb) then
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' GHF orbital coefficients' write(*,'(A50)') ' GHF orbital coefficients '
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas2,nBas2,c) call matout(nBas2,nBas2,c)
write(*,*) write(*,*)
end if end if
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' GHF orbital energies' write(*,'(A50)') ' GHF orbital energies (au) '
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas2,1,e) call matout(nBas2,1,e)
write(*,*) write(*,*)

4
src/HF/print_RHF.f90
View File

@ -73,7 +73,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
if(dump_orb) then if(dump_orb) then
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' RHF orbital coefficients' write(*,'(A50)') ' RHF orbital coefficients '
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,cHF) call matout(nBas,nBas,cHF)
write(*,*) write(*,*)
@ -81,7 +81,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' RHF orbital energies (au) ' write(*,'(A50)') ' RHF orbital energies (au) '
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eHF) call vecout(nBas,eHF)
write(*,*) write(*,*)
end subroutine end subroutine

4
src/HF/print_ROHF.f90
View File

@ -90,7 +90,7 @@ subroutine print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') '-------------------------------------------------'
write(*,'(A45)') ' Dipole moment (Debye) ' write(*,'(A45)') ' Dipole moment (Debye) '
write(*,'(19X,4A10)') 'X','Y','Z','Tot.' 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(*,'(A60)') '-------------------------------------------------'
write(*,*) write(*,*)
@ -102,7 +102,7 @@ subroutine print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
call matout(nBas,nBas,c) call matout(nBas,nBas,c)
write(*,*) write(*,*)
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' ROHF orbital energies ' write(*,'(A50)') ' ROHF orbital energies (au) '
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,e) call matout(nBas,1,e)
write(*,*) write(*,*)

2
src/QuAcK/QuAcK.f90
View File

@ -172,7 +172,7 @@ program QuAcK
!---------------------! !---------------------!
doRQuAcK = .false. doRQuAcK = .false.
if(doRHF) doRQuAcK = .true. if(doRHF .or. doROHF) doRQuAcK = .true.
doUQuAcK = .false. doUQuAcK = .false.
if(doUHF) doUQuAcK = .true. if(doUHF) doUQuAcK = .true.