From 563aff8286183d0952a2a0a3dc99fbf08ed17f47 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 4 Oct 2020 12:40:44 +0200 Subject: [PATCH] dipole UHF OK --- input/methods | 2 +- src/HF/RHF.f90 | 10 +++++++--- src/HF/UHF.f90 | 13 ++++++++++--- src/HF/dipole_moment.f90 | 2 +- src/HF/print_RHF.f90 | 6 +++--- src/HF/print_UHF.f90 | 8 +++++++- src/QuAcK/QuAcK.f90 | 3 ++- 7 files changed, 31 insertions(+), 13 deletions(-) diff --git a/input/methods b/input/methods index 9a5f0de..3a1c5a4 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF MOM - T F F + F T F # MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 7770314..a907154 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,4 +1,4 @@ -subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,Enuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P) +subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P) ! Perform restricted Hartree-Fock calculation @@ -29,7 +29,12 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,Enuc,nBas,nO,S,T integer :: nSCF integer :: nBasSq integer :: n_diis - double precision :: ET,EV,EJ,EK + 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 @@ -43,7 +48,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,Enuc,nBas,nO,S,T double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) double precision,allocatable :: ON(:) - double precision :: dipole(ncart) ! Output variables diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index ac0e9c5..53678e2 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -1,4 +1,4 @@ -subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUHF,e,c,P) +subroutine UHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P) ! Perform unrestricted Hartree-Fock calculation @@ -13,6 +13,11 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH double precision,intent(in) :: thresh integer,intent(in) :: nBas + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + integer,intent(in) :: nO(nspin) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) @@ -20,7 +25,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH 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) :: ENuc + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -33,6 +38,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH double precision :: EV(nspin) double precision :: EJ(nsp) double precision :: Ex(nspin) + double precision :: dipole(ncart) double precision,allocatable :: cp(:,:,:) double precision,allocatable :: J(:,:,:) @@ -233,6 +239,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH ! Compute final UHF energy - call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF) + call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) + call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) end subroutine UHF diff --git a/src/HF/dipole_moment.f90 b/src/HF/dipole_moment.f90 index 49b11f8..e0a1596 100644 --- a/src/HF/dipole_moment.f90 +++ b/src/HF/dipole_moment.f90 @@ -44,7 +44,7 @@ subroutine dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) do mu=1,nBas do nu=1,nBas - dipole(ixyz) = dipole(ixyz) - P(mu,nu)*dipole_int(mu,nu,ixyz) + dipole(ixyz) = dipole(ixyz) - P(mu,nu)*dipole_int(mu,nu,ixyz) enddo enddo diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 56bd1af..985333c 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -55,9 +55,9 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) write(*,'(A36,F13.6)') ' HF LUMO energy (eV) = ',eHF(LUMO)*HaToeV write(*,'(A36,F13.6)') ' HF HOMO-LUMO gap (eV) = ',Gap*HaToeV write(*,'(A50)') '-----------------------------------------' - write(*,'(A36)') ' Dipole moment (Debye) ' - write(*,'(4A10)') 'X','Y','Z','Tot.' - write(*,'(4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD + 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(*,*) diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index f5e9eb5..064d13c 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -1,4 +1,4 @@ -subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF) +subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) ! Print one- and two-electron energies and other stuff for UHF calculation @@ -16,7 +16,9 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF) double precision,intent(in) :: EJ(nsp) double precision,intent(in) :: Ex(nspin) double precision,intent(in) :: EUHF + double precision,intent(in) :: dipole(ncart) + integer :: ixyz integer :: i,j integer :: ispin double precision :: HOMO(nspin) @@ -93,6 +95,10 @@ subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF) write(*,'(A40,F13.6)') ' (exact) :',S2_exact write(*,'(A40,F13.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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index c84671d..29e771e 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -292,7 +292,8 @@ program QuAcK unrestricted = .true. call cpu_time(start_HF) - call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,EUHF,eHF,cHF,PHF) + call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & + nBas,nO,S,T,V,Hc,ERI_AO,dipole_int,X,EUHF,eHF,cHF,PHF) call cpu_time(end_HF) t_HF = end_HF - start_HF