From dd568d88068f61e130fc91cecd9a16ca510714ed Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 8 Mar 2021 20:09:54 +0100 Subject: [PATCH] evUGF2 --- input/methods | 4 +- src/GF/evUGF2.f90 | 194 ++++++++++++++++++++++++++++++++++++++++ src/GF/print_evUGF2.f90 | 81 +++++++++++++++++ src/QuAcK/QuAcK.f90 | 18 +++- 4 files changed, 292 insertions(+), 5 deletions(-) create mode 100644 src/GF/evUGF2.f90 create mode 100644 src/GF/print_evUGF2.f90 diff --git a/input/methods b/input/methods index 4673c55..7fa0bc0 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/src/GF/evUGF2.f90 b/src/GF/evUGF2.f90 new file mode 100644 index 0000000..5072adb --- /dev/null +++ b/src/GF/evUGF2.f90 @@ -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 diff --git a/src/GF/print_evUGF2.f90 b/src/GF/print_evUGF2.f90 new file mode 100644 index 0000000..f04f6f3 --- /dev/null +++ b/src/GF/print_evUGF2.f90 @@ -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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 1b191d4..8bbff6e 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -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