From 435d44391d398ec5ff62da0870ec32dcfc1cde2b Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 25 Sep 2020 15:05:06 +0200 Subject: [PATCH] clean up unrestricted code for one electron --- input/basis | 31 +------------------------------ input/methods | 4 ++-- input/molecule | 5 ++--- input/molecule.xyz | 5 ++--- src/QuAcK/QuAcK.f90 | 5 ++--- src/QuAcK/UG0W0.f90 | 2 +- src/QuAcK/evUGW.f90 | 5 ++--- src/QuAcK/print_UG0W0.f90 | 29 +++++++++++++++++++---------- src/QuAcK/print_UHF.f90 | 30 ++++++++++++++++++------------ src/QuAcK/print_evUGW.f90 | 29 +++++++++++++++++++---------- src/QuAcK/print_excitation.f90 | 2 +- 11 files changed, 69 insertions(+), 78 deletions(-) diff --git a/input/basis b/input/basis index c7ec793..79a747a 100644 --- a/input/basis +++ b/input/basis @@ -1,33 +1,4 @@ -1 6 -S 8 - 1 6665.0000000 0.0006920 - 2 1000.0000000 0.0053290 - 3 228.0000000 0.0270770 - 4 64.7100000 0.1017180 - 5 21.0600000 0.2747400 - 6 7.4950000 0.4485640 - 7 2.7970000 0.2850740 - 8 0.5215000 0.0152040 -S 8 - 1 6665.0000000 -0.0001460 - 2 1000.0000000 -0.0011540 - 3 228.0000000 -0.0057250 - 4 64.7100000 -0.0233120 - 5 21.0600000 -0.0639550 - 6 7.4950000 -0.1499810 - 7 2.7970000 -0.1272620 - 8 0.5215000 0.5445290 -S 1 - 1 0.1596000 1.0000000 -P 3 - 1 9.4390000 0.0381090 - 2 2.0020000 0.2094800 - 3 0.5456000 0.5085570 -P 1 - 1 0.1517000 1.0000000 -D 1 - 1 0.5500000 1.0000000 -2 3 +1 3 S 3 1 13.0100000 0.0196850 2 1.9620000 0.1379770 diff --git a/input/methods b/input/methods index 81c1b6d..58c5d33 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF MOM F T F # MP2* MP3 MP2-F12 - F F F + T F F # CCD CCSD CCSD(T) F F F # drCCD rCCD lCCD pCCD @@ -13,7 +13,7 @@ # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW - T T F + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/molecule b/input/molecule index e2a4fd3..fd4bfbe 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 4 3 0 0 + 1 1 0 0 0 # Znuc x y z - C 0. 0. -0.16245872 - H 0. 0. 1.93436816 + H 0. 0. 0. diff --git a/input/molecule.xyz b/input/molecule.xyz index 7a4f218..3dca9a4 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - C 0.0000000000 0.0000000000 -0.0859694585 - H 0.0000000000 0.0000000000 1.0236236215 + H 0.0000000000 0.0000000000 0.0000000000 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index f7310a1..2efdaa2 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -802,9 +802,8 @@ program QuAcK ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,PHF,cHF,eHF,eG0W0) else - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0) + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & + linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0) end if diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index c13a2a9..5e1b6fc 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -118,7 +118,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - if(print_W) call print_excitation('RPA@UHF',5,nS_sc,OmRPA) + if(print_W) call print_excitation('RPA@UHF ',5,nS_sc,OmRPA) !----------------------! ! Excitation densities ! diff --git a/src/QuAcK/evUGW.f90 b/src/QuAcK/evUGW.f90 index 3c8b2d8..1ba36ac 100644 --- a/src/QuAcK/evUGW.f90 +++ b/src/QuAcK/evUGW.f90 @@ -127,9 +127,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin),OmRPA(nS_sc), & - XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin), & - error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin)) + allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc), & + XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin),error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin)) ! Initialization diff --git a/src/QuAcK/print_UG0W0.f90 b/src/QuAcK/print_UG0W0.f90 index f678a1e..0e8092e 100644 --- a/src/QuAcK/print_UG0W0.f90 +++ b/src/QuAcK/print_UG0W0.f90 @@ -16,15 +16,24 @@ subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA) double precision,intent(in) :: eGW(nBas,nspin) integer :: p - double precision :: HOMO - double precision :: LUMO - double precision :: Gap + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) ! HOMO and LUMO - HOMO = max(eGW(nO(1),1),eGW(nO(2),2)) - LUMO = min(eGW(nO(1)+1,1),eGW(nO(2)+1,2)) - Gap = LUMO - HOMO + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = eGW(nO(ispin),ispin) + LUMO(ispin) = eGW(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = e(1,ispin) + Gap(ispin) = 0d0 + end if + end do ! Dump results @@ -34,7 +43,7 @@ subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA) write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') & - '|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|' + '|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|' write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') & '|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|' write(*,*)'-------------------------------------------------------------------------------& @@ -48,9 +57,9 @@ subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA) write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'G0W0 HOMO energy (eV):',HOMO*HaToeV - write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',LUMO*HaToeV - write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV + write(*,'(2X,A30,F15.6)') 'G0W0 HOMO energy (eV):',maxval(HOMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',minval(LUMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA diff --git a/src/QuAcK/print_UHF.f90 b/src/QuAcK/print_UHF.f90 index 6c52e4d..0ad8c89 100644 --- a/src/QuAcK/print_UHF.f90 +++ b/src/QuAcK/print_UHF.f90 @@ -16,18 +16,24 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) double precision,intent(in) :: Ex(nspin) double precision,intent(in) :: EUHF - integer :: HOMO(nspin) - integer :: LUMO(nspin) + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) double precision :: Gap(nspin) ! HOMO and LUMO - HOMO(:) = nO(:) - - LUMO(:) = HOMO(:) + 1 - - Gap(1) = e(LUMO(1),1) - e(HOMO(1),1) - Gap(2) = e(LUMO(2),2) - e(HOMO(2),2) + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = e(nO(ispin),ispin) + LUMO(ispin) = e(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = e(1,ispin) + Gap(ispin) = 0d0 + end if + end do ! Dump results @@ -62,12 +68,12 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' write(*,'(A40,1X,F16.10,A3)') ' UHF energy: ',EUHF + ENuc,' au' write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,F13.6,A3)') ' UHF HOMO a energy:',e(HOMO(1),1)*HatoeV,' eV' - write(*,'(A40,F13.6,A3)') ' UHF LUMO a energy:',e(LUMO(1),1)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF HOMO a energy:',HOMO(1)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF LUMO a energy:',LUMO(1)*HatoeV,' eV' write(*,'(A40,F13.6,A3)') ' UHF HOMOa-LUMOa gap:',Gap(1)*HatoeV,' eV' write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,F13.6,A3)') ' UHF HOMO b energy:',e(HOMO(2),2)*HatoeV,' eV' - write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',e(LUMO(2),2)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF HOMO b energy:',HOMO(2)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',LUMO(2)*HatoeV,' eV' write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV' write(*,'(A60)') '-------------------------------------------------' write(*,*) diff --git a/src/QuAcK/print_evUGW.f90 b/src/QuAcK/print_evUGW.f90 index 7c32002..984d88c 100644 --- a/src/QuAcK/print_evUGW.f90 +++ b/src/QuAcK/print_evUGW.f90 @@ -18,15 +18,24 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) double precision,intent(in) :: eGW(nBas,nspin) integer :: p - double precision :: HOMO - double precision :: LUMO - double precision :: Gap + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) ! HOMO and LUMO - HOMO = max(eGW(nO(1),1),eGW(nO(2),2)) - LUMO = min(eGW(nO(1)+1,1),eGW(nO(2)+1,2)) - Gap = LUMO - HOMO + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = eGW(nO(ispin),ispin) + LUMO(ispin) = eGW(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = e(1,ispin) + Gap(ispin) = 0d0 + end if + end do ! Dump results @@ -40,7 +49,7 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') & - '|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|' + '|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|' write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') & '|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|' write(*,*)'-------------------------------------------------------------------------------& @@ -58,9 +67,9 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) write(*,'(2X,A14,F15.5)')'Convergence = ',Conv write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'evGW HOMO energy (eV):',HOMO*HaToeV - write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',LUMO*HaToeV - write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW HOMO energy (eV):',maxval(HOMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',minval(LUMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90 index 7aec451..c2bbdb4 100644 --- a/src/QuAcK/print_excitation.f90 +++ b/src/QuAcK/print_excitation.f90 @@ -7,7 +7,7 @@ subroutine print_excitation(method,ispin,nS,Omega) ! Input variables - character*12,intent(in) :: method + character*12,intent(in) :: method integer,intent(in) :: ispin,nS double precision,intent(in) :: Omega(nS)