diff --git a/input/options b/input/options index 05cee9f..93da02b 100644 --- a/input/options +++ b/input/options @@ -9,7 +9,7 @@ # 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 - 128 0.00001 T 5 T 0.00367493 F F F F F + 128 0.00001 T 5 T 0.00367493 F F T F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn diff --git a/src/MBPT/print_qsUGW.f90 b/src/MBPT/print_qsUGW.f90 index 38aa146..78fd5e7 100644 --- a/src/MBPT/print_qsUGW.f90 +++ b/src/MBPT/print_qsUGW.f90 @@ -1,4 +1,5 @@ -subroutine print_qsUGW(nBas,nO,Ov,nSCF,Conv,thresh,eGW,cGW,PGW,T,V,J,K,ENuc,EHF,SigC,Z,EcRPA,dipole) +subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,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 @@ -9,16 +10,21 @@ subroutine print_qsUGW(nBas,nO,Ov,nSCF,Conv,thresh,eGW,cGW,PGW,T,V,J,K,ENuc,EHF, integer,intent(in) :: nBas integer,intent(in) :: nO(nspin) - double precision,intent(in) :: Ov(nBas,nBas) integer,intent(in) :: nSCF double precision,intent(in) :: ENuc - double precision,intent(in) :: EHF + double precision,intent(in) :: ET(nspin) + double precision,intent(in) :: EV(nspin) + double precision,intent(in) :: EJ(nsp) + double precision,intent(in) :: Ex(nspin) + double precision,intent(in) :: Ec(nsp) double precision,intent(in) :: EcRPA + double precision,intent(in) :: EqsGW double precision,intent(in) :: Conv double precision,intent(in) :: thresh double precision,intent(in) :: eGW(nBas,nspin) double precision,intent(in) :: cGW(nBas,nBas,nspin) double precision,intent(in) :: PGW(nBas,nBas,nspin) + double precision,intent(in) :: Ov(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: J(nBas,nBas,nspin) @@ -34,12 +40,6 @@ subroutine print_qsUGW(nBas,nO,Ov,nSCF,Conv,thresh,eGW,cGW,PGW,T,V,J,K,ENuc,EHF, double precision :: HOMO(nspin) double precision :: LUMO(nspin) double precision :: Gap(nspin) - double precision :: ET(nspin) - double precision :: EV(nspin) - double precision :: EJ(nsp) - double precision :: Ex(nspin) - double precision :: Ec(nsp) - double precision :: EqsGW double precision :: S_exact,S2_exact double precision :: S,S2 double precision,external :: trace_matrix @@ -64,44 +64,6 @@ subroutine print_qsUGW(nBas,nO,Ov,nSCF,Conv,thresh,eGW,cGW,PGW,T,V,J,K,ENuc,EHF, S_exact = 0.5d0*dble(nO(1) - nO(2)) S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2) -!------------------------------------------------------------------------ -! Compute total energy -!------------------------------------------------------------------------ - -! Kinetic energy - - do ispin=1,nspin - ET(ispin) = trace_matrix(nBas,matmul(PGW(:,:,ispin),T(:,:))) - end do - -! Potential energy - - do ispin=1,nspin - EV(ispin) = trace_matrix(nBas,matmul(PGW(:,:,ispin),V(:,:))) - end do - -! Coulomb energy - - EJ(1) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,1),J(:,:,1))) - EJ(2) = trace_matrix(nBas,matmul(PGW(:,:,1),J(:,:,2))) - EJ(3) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,2),J(:,:,2))) - -! Exchange energy - - do ispin=1,nspin - Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,ispin),K(:,:,ispin))) - end do - -! Correlation energy - - Ec(1) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,1),SigC(:,:,1))) - Ec(2) = trace_matrix(nBas,matmul(PGW(:,:,1),SigC(:,:,2))) - Ec(3) = 0.5d0*trace_matrix(nBas,matmul(PGW(:,:,2),SigC(:,:,2))) - -! Total energy - - EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)) - ! Dump results write(*,*)'-------------------------------------------------------------------------------& diff --git a/src/MBPT/qsUGW.f90 b/src/MBPT/qsUGW.f90 index 3be8732..1852e4e 100644 --- a/src/MBPT/qsUGW.f90 +++ b/src/MBPT/qsUGW.f90 @@ -71,7 +71,13 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS integer :: nS_aa,nS_bb,nS_sc double precision :: dipole(ncart) + double precision :: ET(nspin) + double precision :: EV(nspin) + double precision :: EJ(nsp) + double precision :: Ex(nspin) + double precision :: Ec(nsp) double precision :: EcRPA + double precision :: EqsGW double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: Conv @@ -296,10 +302,50 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is))) end do + !------------------------------------------------------------------------ + ! Compute total energy + !------------------------------------------------------------------------ + + ! Kinetic energy + + do ispin=1,nspin + ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) + end do + + ! Potential energy + + do ispin=1,nspin + EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),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(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))) + 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))) + + ! Total energy + + EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)) + + !------------------------------------------------------------------------ ! Print results + !------------------------------------------------------------------------ call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsUGW(nBas,nO,S,nSCF,Conv,thresh,eGW,c,P,T,V,J,K,ENuc,EHF,SigC,Z,EcRPA,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) enddo !------------------------------------------------------------------------