diff --git a/input/methods b/input/methods index 86c7c20..af2bd81 100644 --- a/input/methods +++ b/input/methods @@ -13,7 +13,7 @@ # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW* - F F F + F F T # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 4c47176..02bfdc4 100644 --- a/input/options +++ b/input/options @@ -1,11 +1,11 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess - 128 0.0000001 T 10 1 1 T + 128 0.0000001 T 5 1 1 F # MP: # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - T T T T T + F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T T T F + F F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/mol/furan.xyz b/mol/furan.xyz index 665b07f..9de0f7d 100644 --- a/mol/furan.xyz +++ b/mol/furan.xyz @@ -1,5 +1,5 @@ 9 -Furan,^1A_1,CC3,aug-cc-pVTZ + C 0.00000000 1.09204092 -0.31777753 C 0.00000000 -1.09204092 -0.31777753 C 0.00000000 0.71623383 0.98604985 @@ -8,4 +8,4 @@ O 0.00000000 0.00000000 -1.13214994 H 0.00000000 2.04440888 -0.81369302 H 0.00000000 -2.04440888 -0.81369302 H 0.00000000 1.37146217 1.83713421 -H 0.00000000 -1.37146217 1.83713421 \ No newline at end of file +H 0.00000000 -1.37146217 1.83713421 diff --git a/src/MBPT/print_evUGW.f90 b/src/MBPT/print_evUGW.f90 index 984d88c..7a8f23b 100644 --- a/src/MBPT/print_evUGW.f90 +++ b/src/MBPT/print_evUGW.f90 @@ -1,4 +1,4 @@ -subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) +subroutine print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA) ! Print one-electron energies and other stuff for evGW @@ -9,10 +9,10 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) integer,intent(in) :: nO(nspin) integer,intent(in) :: nSCF double precision,intent(in) :: ENuc - double precision,intent(in) :: EHF + double precision,intent(in) :: ERHF double precision,intent(in) :: EcRPA double precision,intent(in) :: Conv - double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: SigC(nBas,nspin) double precision,intent(in) :: Z(nBas,nspin) double precision,intent(in) :: eGW(nBas,nspin) @@ -32,7 +32,7 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) Gap(ispin) = LUMO(ispin) - HOMO(ispin) else HOMO(ispin) = 0d0 - LUMO(ispin) = e(1,ispin) + LUMO(ispin) = eGW(1,ispin) Gap(ispin) = 0d0 end if end do @@ -57,7 +57,7 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) do p=1,nBas write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & - '|',p,'|',e(p,1)*HaToeV,e(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & + '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|' enddo @@ -72,7 +72,7 @@ subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) 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 + write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + ERHF + EcRPA write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' diff --git a/src/MBPT/print_qsGW.f90 b/src/MBPT/print_qsGW.f90 index 74ecc10..ecf6224 100644 --- a/src/MBPT/print_qsGW.f90 +++ b/src/MBPT/print_qsGW.f90 @@ -1,6 +1,5 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcRPA,EqsGW,dipole) - ! Print one-electron energies and other stuff for qsGW implicit none @@ -32,11 +31,13 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z LUMO = HOMO + 1 Gap = eGW(LUMO)-eGW(HOMO) +! Compute energies + ET = trace_matrix(nBas,matmul(P,T)) EV = trace_matrix(nBas,matmul(P,V)) EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) - Ec = 0.5d0*trace_matrix(nBas,matmul(P,SigC)) + Ec = 0.50d0*trace_matrix(nBas,matmul(P,SigC)) EqsGW = ET + EV + EJ + Ex + Ec ! Dump results @@ -49,12 +50,12 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z endif write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & - '|','#','|','e_HF (eV)','|','e_QP-e_HF (eV)','|','Z','|','e_QP (eV)','|' + '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' write(*,*)'-------------------------------------------------------------------------------' do x=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',x,'|',eHF(x)*HaToeV,'|',(eGW(x)-eHF(x))*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' + '|',x,'|',eHF(x)*HaToeV,'|',SigC(x,x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' enddo write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/MBPT/print_qsUGW.f90 b/src/MBPT/print_qsUGW.f90 index 78fd5e7..a779af9 100644 --- a/src/MBPT/print_qsUGW.f90 +++ b/src/MBPT/print_qsUGW.f90 @@ -1,4 +1,4 @@ -subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eGW,cGW,PGW,Ov,T,V,J,K, & +subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, & ENuc,ET,EV,EJ,Ex,Ec,EcRPA,EqsGW,SigC,Z,dipole) ! Print one-electron energies and other stuff for qsUGW @@ -21,6 +21,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eGW,cGW,PGW,Ov,T,V,J,K, & double precision,intent(in) :: EqsGW double precision,intent(in) :: Conv double precision,intent(in) :: thresh + double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: eGW(nBas,nspin) double precision,intent(in) :: cGW(nBas,nBas,nspin) double precision,intent(in) :: PGW(nBas,nBas,nspin) @@ -84,7 +85,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eGW,cGW,PGW,Ov,T,V,J,K, & do p=1,nBas write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & - '|',p,'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|',SigC(p,p,1)*HaToeV,SigC(p,p,2)*HaToeV,'|', & + '|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,p,1)*HaToeV,SigC(p,p,2)*HaToeV,'|', & Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|' enddo @@ -94,15 +95,15 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eGW,cGW,PGW,Ov,T,V,J,K, & write(*,'(2X,A19,F15.5)')'max(|FPS - SPF|) = ',Conv write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'qsGW HOMO energy (eV):',maxval(HOMO(:))*HaToeV - write(*,'(2X,A30,F15.6)') 'qsGW LUMO energy (eV):',minval(LUMO(:))*HaToeV - write(*,'(2X,A30,F15.6)') 'qsGW HOMO-LUMO gap (eV):',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV + write(*,'(2X,A30,F15.6)') 'qsUGW HOMO energy (eV):',maxval(HOMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'qsUGW LUMO energy (eV):',minval(LUMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'qsUGW HOMO-LUMO gap (eV):',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'qsGW total energy =',EqsGW + ENuc - write(*,'(2X,A30,F15.6)') 'qsGW exchange energy =',sum(Ex(:)) - write(*,'(2X,A30,F15.6)') 'qsGW correlation energy =',sum(Ec(:)) - write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') ' qsUGW total energy =',EqsGW + ENuc + write(*,'(2X,A30,F15.6)') ' qsUGW exchange energy =',sum(Ex(:)) + write(*,'(2X,A30,F15.6)') ' qsUGW correlation energy =',sum(Ec(:)) + write(*,'(2X,A30,F15.6)') 'RPA@qsUGW correlation energy =',EcRPA write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' write(*,*) diff --git a/src/MBPT/qsUGW.f90 b/src/MBPT/qsUGW.f90 index 1852e4e..0da7d9c 100644 --- a/src/MBPT/qsUGW.f90 +++ b/src/MBPT/qsUGW.f90 @@ -308,33 +308,34 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS ! Kinetic energy - do ispin=1,nspin - ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) + do is=1,nspin + ET(is) = trace_matrix(nBas,matmul(P(:,:,is),T(:,:))) end do ! Potential energy - do ispin=1,nspin - EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) + do is=1,nspin + EV(is) = trace_matrix(nBas,matmul(P(:,:,is),V(:,:))) end do ! Coulomb energy EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) - EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) + EJ(2) = 1.0d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) ! Exchange energy - do ispin=1,nspin - Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin))) + do is=1,nspin + Ex(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is))) end do ! Correlation energy - Ec(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),SigC(:,:,1))) - Ec(2) = trace_matrix(nBas,matmul(P(:,:,1),SigC(:,:,2))) - Ec(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),SigC(:,:,2))) + Ec(1) = 0.25d0*trace_matrix(nBas,matmul(P(:,:,1),SigCp(:,:,1))) + Ec(2) = 0.25d0*trace_matrix(nBas,matmul(P(:,:,1),SigCp(:,:,2))) & + + 0.25d0*trace_matrix(nBas,matmul(P(:,:,2),SigCp(:,:,1))) + Ec(3) = 0.25d0*trace_matrix(nBas,matmul(P(:,:,2),SigCp(:,:,2))) ! Total energy @@ -345,7 +346,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS !------------------------------------------------------------------------ call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EcRPA,EqsGW,SigC,Z,dipole) + call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,Ec,EcRPA,EqsGW,SigCp,Z,dipole) enddo !------------------------------------------------------------------------