4
1
mirror of https://github.com/pfloos/quack synced 2024-06-26 15:12:17 +02:00
This commit is contained in:
Pierre-Francois Loos 2021-03-08 20:09:54 +01:00
parent 72470287f7
commit dd568d8806
4 changed files with 292 additions and 5 deletions

View File

@ -10,8 +10,8 @@
F F F F F
# RPA* RPAx* ppRPA
F F F
# G0F2* evGF2 qsGF2* G0F3 evGF3
T F F F F
# G0F2* evGF2* qsGF2* G0F3 evGF3
F T F F F
# G0W0* evGW* qsGW*
F F F
# G0T0 evGT qsGT

194
src/GF/evUGF2.f90 Normal file
View File

@ -0,0 +1,194 @@
subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_aa,dipole_int_bb,cHF,eHF)
! Perform self-consistent eigenvalue-only GW calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
logical,intent(in) :: BSE
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
! Local variables
logical :: linear_mixing
integer :: is
integer :: ispin
integer :: nSCF
integer :: n_diis
double precision :: rcond(nspin)
double precision :: Conv
double precision :: Ec(nsp)
double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: alpha
double precision,allocatable :: error_diis(:,:,:)
double precision,allocatable :: e_diis(:,:,:)
double precision,allocatable :: eGF2(:,:)
double precision,allocatable :: eOld(:,:)
double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: SigC(:,:)
! Hello world
write(*,*)
write(*,*)'**************************************************'
write(*,*)'| Self-consistent unrestricted evGF2 calculation |'
write(*,*)'**************************************************'
write(*,*)
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Linear mixing
linear_mixing = .false.
alpha = 0.2d0
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(eGF2(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin), &
error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin))
! Initialization
nSCF = 0
ispin = 1
n_diis = 0
Conv = 1d0
e_diis(:,:,:) = 0d0
error_diis(:,:,:) = 0d0
eGF2(:,:) = eHF(:,:)
eOld(:,:) = eHF(:,:)
Z(:,:) = 1d0
!------------------------------------------------------------------------
! Main loop
!------------------------------------------------------------------------
do while(Conv > thresh .and. nSCF <= maxSCF)
!------------------------------------------------!
! Compute self-energy and renormalization factor !
!------------------------------------------------!
call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z)
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
eGF2(:,:) = eHF(:,:) + SigC(:,:)
! Convergence criteria
Conv = maxval(abs(eGF2(:,:) - eOld(:,:)))
! Compute MP2 correlation energy
call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EUHF,eHF,Ec)
! Print results
call print_evUGF2(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec)
! Linear mixing or DIIS extrapolation
if(linear_mixing) then
eGF2(:,:) = alpha*eGF2(:,:) + (1d0 - alpha)*eOld(:,:)
else
n_diis = min(n_diis+1,max_diis)
do is=1,nspin
call DIIS_extrapolation(rcond(ispin),nBas,nBas,n_diis,error_diis(:,1:n_diis,is), &
e_diis(:,1:n_diis,is),eGF2(:,is)-eOld(:,is),eGF2(:,is))
end do
! Reset DIIS if required
if(minval(rcond(:)) < 1d-15) n_diis = 0
endif
! Save quasiparticles energy for next cycle
eOld(:,:) = eGF2(:,:)
! Increment
nSCF = nSCF + 1
enddo
!------------------------------------------------------------------------
! End main loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF+1) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
endif
! Deallocate memory
deallocate(eOld,Z,SigC,error_diis,e_diis)
! Perform BSE calculation
if(BSE) then
print*,'!!! BSE2 NYI for evUGF2 !!!'
endif
end subroutine evUGF2

81
src/GF/print_evUGF2.f90 Normal file
View File

@ -0,0 +1,81 @@
subroutine print_evUGF2(nBas,nO,nSCF,Conv,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec)
! Print one-electron energies and other stuff for evGF2
implicit none
include 'parameters.h'
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: Ec(nsp)
double precision,intent(in) :: Conv
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) :: eGF2(nBas,nspin)
integer :: p
integer :: ispin
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
! HOMO and LUMO
do ispin=1,nspin
if(nO(ispin) > 0) then
HOMO(ispin) = eGF2(nO(ispin),ispin)
LUMO(ispin) = eGF2(nO(ispin)+1,ispin)
Gap(ispin) = LUMO(ispin) - HOMO(ispin)
else
HOMO(ispin) = 0d0
LUMO(ispin) = eGF2(1,ispin)
Gap(ispin) = 0d0
end if
end do
! Dump results
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
if(nSCF < 10) then
write(*,'(1X,A21,I1,A2,A12)')' Self-consistent evG',nSCF,'F2',' calculation'
else
write(*,'(1X,A21,I2,A2,A12)')' Self-consistent evG',nSCF,'F2',' calculation'
endif
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') &
'|',' ','|','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(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
do p=1,nBas
write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') &
'|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', &
Z(p,1),Z(p,2),'|',eGF2(p,1)*HaToeV,eGF2(p,2)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') 'evGF2 HOMO energy:',maxval(HOMO(:))*HaToeV,' eV'
write(*,'(2X,A30,F15.6,A3)') 'evGF2 LUMO energy:',minval(LUMO(:))*HaToeV,' eV'
write(*,'(2X,A30,F15.6,A3)') 'evGF2 HOMO-LUMO gap :',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') ' evGF2 total energy :',ENuc + EUHF + sum(Ec(:)),' au'
write(*,'(2X,A30,F15.6,A3)') ' evGF2 correlation energy:',sum(Ec(:)),' au'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,*)
end subroutine print_evUGF2

View File

@ -841,9 +841,21 @@ program QuAcK
if(doevGF2) then
call cpu_time(start_GF2)
call evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF_GF,thresh_GF,n_diis_GF, &
singlet,triplet,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI_MO,dipole_int_MO,eHF)
if(unrestricted) then
call evUGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, &
dipole_int_aa,dipole_int_bb,cHF,eHF)
else
call evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF_GF,thresh_GF,n_diis_GF, &
singlet,triplet,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, &
ERI_MO,dipole_int_MO,eHF)
end if
call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2