From 114ac2d1d1594a7e193e0079a9c45c8bec3c452c Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sat, 27 Mar 2021 15:03:54 +0100 Subject: [PATCH] debug qsGW and qsGF --- input/methods | 6 +- input/options | 8 +- mol/h2.xyz | 2 +- src/GF/G0W0.f90 | 238 --------------------------------------- src/GF/qsGF2.f90 | 23 ++-- src/GF/qsUGF2.f90 | 35 ++++-- src/HF/RHF.f90 | 10 +- src/HF/UHF.f90 | 14 +-- src/MBPT/evGW.f90 | 10 +- src/MBPT/print_qsGW.f90 | 11 +- src/MBPT/print_qsUGW.f90 | 2 +- src/MBPT/qsGW.f90 | 31 +++-- src/MBPT/qsUGW.f90 | 38 +++++-- 13 files changed, 114 insertions(+), 314 deletions(-) delete mode 100644 src/GF/G0W0.f90 diff --git a/input/methods b/input/methods index c965dd2..098169e 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - T F F F + F T F F # MP2* MP3 MP2-F12 F F F # CCD DCD CCSD CCSD(T) @@ -11,9 +11,9 @@ # RPA* RPAx* ppRPA F F F # G0F2* evGF2* qsGF2* G0F3 evGF3 - F F F F F + F F T F F # G0W0* evGW* qsGW* - T F F + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 367898e..6c25527 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability - 128 0.0000001 T 5 1 1 T F + 1024 0.0000001 T 5 2 1 T F # MP: # CC: maxSCF thresh DIIS n_diis @@ -7,12 +7,12 @@ # spin: TDA singlet triplet spin_conserved spin_flip F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm - 256 0.00001 T 5 T 0.001 3 + 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 - 256 0.00001 T 5 T 0.0 F F F F F + 256 0.0000001 T 5 T 0.0 F F F F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T T T F + F T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/mol/h2.xyz b/mol/h2.xyz index f8e2dab..fe42126 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0.0 0.0 0.0 -H 0.0 0.0 0.71 +H 0.0 0.0 1.0 diff --git a/src/GF/G0W0.f90 b/src/GF/G0W0.f90 deleted file mode 100644 index d2c7fa7..0000000 --- a/src/GF/G0W0.f90 +++ /dev/null @@ -1,238 +0,0 @@ -subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0) - -! Perform G0W0 calculation - - implicit none - include 'parameters.h' - include 'quadrature.h' - -! Input variables - - logical,intent(in) :: doACFDT - logical,intent(in) :: exchange_kernel - logical,intent(in) :: doXBS - logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX - logical,intent(in) :: BSE - logical,intent(in) :: TDA_W - logical,intent(in) :: TDA - logical,intent(in) :: dBSE - logical,intent(in) :: dTDA - logical,intent(in) :: evDyn - logical,intent(in) :: singlet - logical,intent(in) :: triplet - logical,intent(in) :: linearize - double precision,intent(in) :: eta - - integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF - double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: Vxc(nBas) - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) - -! Local variables - - logical :: print_W = .true. - integer :: ispin - double precision :: EcRPA - double precision :: EcBSE(nspin) - double precision :: EcAC(nspin) - double precision :: EcGM - double precision,allocatable :: SigX(:) - double precision,allocatable :: SigC(:) - double precision,allocatable :: Z(:) - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) - - double precision,allocatable :: eG0W0lin(:) - -! Output variables - - double precision :: eG0W0(nBas) - -! Hello world - - write(*,*) - write(*,*)'************************************************' - write(*,*)'| One-shot G0W0 calculation |' - write(*,*)'************************************************' - write(*,*) - -! Initialization - - EcRPA = 0d0 - -! SOSEX correction - - if(SOSEX) then - write(*,*) 'SOSEX correction activated but BUG!' - stop - end if - -! COHSEX approximation - - if(COHSEX) then - write(*,*) 'COHSEX approximation activated!' - write(*,*) - end if - -! TDA for W - - if(TDA_W) then - write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) - end if - -! TDA - - if(TDA) then - write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) - end if - -! Spin manifold - - ispin = 1 - -! Memory allocation - - allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eG0W0lin(nBas)) - -!-------------------! -! Compute screening ! -!-------------------! - - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, & - eHF,ERI_MO,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - - if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) - -!--------------------------! -! Compute spectral weights ! -!--------------------------! - - call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) - -!------------------------! -! Compute GW self-energy ! -!------------------------! - - call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) - call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) - -!--------------------------------! -! Compute renormalization factor ! -!--------------------------------! - - call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) - -!-----------------------------------! -! Solve the quasi-particle equation ! -!-----------------------------------! - - eG0W0lin(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:)) - - ! Linearized or graphical solution? - - if(linearize) then - - write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' - write(*,*) - - eG0W0(:) = eG0W0lin(:) - - else - - write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' - write(*,*) - - call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eG0W0lin,eG0W0) - - ! Find all the roots of the QP equation if necessary - - ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin) - - end if - -! Compute the RPA correlation energy - - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI_MO,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - -!--------------! -! Dump results ! -!--------------! - - call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eG0W0,EcRPA,EcGM) - -! Deallocate memory - - deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eG0W0lin) - -! Plot stuff - -! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA) - -! Perform BSE calculation - - if(BSE) then - - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eG0W0,EcBSE) - - if(exchange_kernel) then - - EcBSE(1) = 0.5d0*EcBSE(1) - EcBSE(2) = 1.5d0*EcBSE(2) - - end if - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1) - write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - -! Compute the BSE correlation energy via the adiabatic connection - - if(doACFDT) then - - write(*,*) '--------------------------------------------------------------' - write(*,*) ' Adiabatic connection version of BSE@UG0W0 correlation energy ' - write(*,*) '--------------------------------------------------------------' - write(*,*) - - if(doXBS) then - - write(*,*) '*** scaled screening version (XBS) ***' - write(*,*) - - end if - - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eG0W0,EcAC) - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (singlet) =',EcAC(1) - write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy (triplet) =',EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 correlation energy =',EcAC(1) + EcAC(2) - write(*,'(2X,A50,F20.10)') 'AC@BSE@UG0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - - end if - - end if - -end subroutine G0W0 diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index 6e8988b..b4b385b 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -64,6 +64,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, double precision,allocatable :: c(:,:) double precision,allocatable :: cp(:,:) double precision,allocatable :: eGF2(:) + double precision,allocatable :: eOld(:) double precision,allocatable :: P(:,:) double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) @@ -101,8 +102,8 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, ! Memory allocation - allocate(eGF2(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & + allocate(eGF2(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & + J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -112,6 +113,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, ispin = 1 Conv = 1d0 P(:,:) = PHF(:,:) + eOld(:) = eHF(:) eGF2(:) = eHF(:) c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 @@ -157,16 +159,15 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, ! Compute commutator and convergence criteria error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) - Conv = maxval(abs(error)) ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) - -! Reset DIIS if required - - if(abs(rcond) < 1d-15) n_diis = 0 + if(abs(rcond) > 1d-7) then + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + else + n_diis = 0 + end if ! Diagonalize Hamiltonian in AO basis @@ -174,6 +175,12 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, cp(:,:) = Fp(:,:) call diagonalize_matrix(nBas,cp,eGF2) c = matmul(X,cp) + SigCp = matmul(transpose(c),matmul(SigCp,c)) + + ! Save quasiparticles energy for next cycle + + Conv = maxval(abs(eGF2 - eOld)) + eOld(:) = eGF2(:) ! Compute new density matrix in the AO basis diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 5c93e01..9c05f6a 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -75,6 +75,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, double precision,allocatable :: F_diis(:,:,:) double precision,allocatable :: c(:,:,:) double precision,allocatable :: cp(:,:,:) + double precision,allocatable :: eOld(:,:) double precision,allocatable :: eGF2(:,:) double precision,allocatable :: P(:,:,:) double precision,allocatable :: F(:,:,:) @@ -117,9 +118,10 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(eGF2(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), & - J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),SigCm(nBas,nBas,nspin), & - Z(nBas,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) + allocate(eGF2(nBas,nspin),eOld(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), & + Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), & + SigCm(nBas,nBas,nspin),Z(nBas,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), & + F_diis(nBasSq,max_diis,nspin)) ! Initialization @@ -205,14 +207,14 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - do is=1,nspin - if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,error_diis(:,1:n_diis,is), & - F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is)) - end do - - ! Reset DIIS if required - - if(minval(rcond(:)) < 1d-15) n_diis = 0 + if(minval(rcond(:)) > 1d-7) then + do is=1,nspin + if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,error_diis(:,1:n_diis,is), & + F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is)) + end do + else + n_diis = 0 + end if ! Transform Fock matrix in orthogonal basis @@ -233,12 +235,23 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, c(:,:,is) = matmul(X(:,:),cp(:,:,is)) end do + ! Back-transform self-energy + + do is=1,nspin + SigCp(:,:,is) = matmul(transpose(c(:,:,is)),matmul(SigCp(:,:,is),c(:,:,is))) + end do + ! Compute density matrix do is=1,nspin P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is))) end do + ! Save quasiparticles energy for next cycle + + Conv = maxval(abs(eGF2(:,:) - eOld(:,:))) + eOld(:,:) = eGF2(:,:) + !------------------------------------------------------------------------ ! Compute total energy !------------------------------------------------------------------------ diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 9fd3d22..05499a5 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -127,11 +127,11 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) - -! Reset DIIS if required - - if(abs(rcond) < 1d-15) n_diis = 0 + if(abs(rcond) > 1d-7) then + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + else + n_diis = 0 + end if ! Diagonalize Fock matrix diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index a55e8bb..9443924 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -179,14 +179,14 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - do ispin=1,nspin - if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,ispin), & + if(minval(rcond(:)) > 1d-7) then + do ispin=1,nspin + if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,ispin), & F_diis(:,1:n_diis,ispin),err(:,:,ispin),F(:,:,ispin)) - end do - -! Reset DIIS if required - - if(minval(rcond(:)) < 1d-15) n_diis = 0 + end do + else + n_diis = 0 + end if !------------------------------------------------------------------------ ! Compute UHF energy diff --git a/src/MBPT/evGW.f90 b/src/MBPT/evGW.f90 index 7d15ef9..85998eb 100644 --- a/src/MBPT/evGW.f90 +++ b/src/MBPT/evGW.f90 @@ -195,11 +195,11 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE else n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW) - -! Reset DIIS if required - - if(abs(rcond) < 1d-15) n_diis = 0 + if(abs(rcond) > 1d-7) then + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW) + else + n_diis = 0 + endif endif diff --git a/src/MBPT/print_qsGW.f90 b/src/MBPT/print_qsGW.f90 index c3995b0..52702be 100644 --- a/src/MBPT/print_qsGW.f90 +++ b/src/MBPT/print_qsGW.f90 @@ -1,5 +1,4 @@ -subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigC,Z, & - ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole) +subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole) ! Print one-electron energies and other stuff for qsGW @@ -23,10 +22,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigC,Z, & double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: c(nBas) - double precision,intent(in) :: P(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas) - double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas) - double precision,intent(in) :: Z(nBas),SigC(nBas,nBas) + double precision,intent(in) :: SigC(nBas,nBas) + double precision,intent(in) :: Z(nBas) double precision,intent(in) :: dipole(ncart) ! Local variables @@ -67,7 +64,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigC,Z, & write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF - write(*,'(2X,A19,F15.5)')'max(|FPS - SPF|) = ',Conv + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv write(*,*)'-------------------------------------------' write(*,'(2X,A30,F15.6,A3)') 'qsGW HOMO energy:',eGW(HOMO)*HaToeV,' eV' write(*,'(2X,A30,F15.6,A3)') 'qsGW LUMO energy:',eGW(LUMO)*HaToeV,' eV' diff --git a/src/MBPT/print_qsUGW.f90 b/src/MBPT/print_qsUGW.f90 index 18594dc..983e148 100644 --- a/src/MBPT/print_qsUGW.f90 +++ b/src/MBPT/print_qsUGW.f90 @@ -92,7 +92,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, & write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' write(*,'(2X,A10,I3)') 'Iteration ',nSCF - write(*,'(2X,A19,F15.5)')'max(|FPS - SPF|) = ',Conv + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv write(*,*)'-------------------------------------------------------------------------------& -------------------------------------------------' write(*,'(2X,A30,F15.6,A3)') 'qsUGW HOMO energy:',maxval(HOMO(:))*HaToeV,' eV' diff --git a/src/MBPT/qsGW.f90 b/src/MBPT/qsGW.f90 index e143ebe..5d54f93 100644 --- a/src/MBPT/qsGW.f90 +++ b/src/MBPT/qsGW.f90 @@ -69,6 +69,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,external :: trace_matrix double precision :: dipole(ncart) + logical :: print_W = .true. double precision,allocatable :: error_diis(:,:) double precision,allocatable :: F_diis(:,:) double precision,allocatable :: OmRPA(:) @@ -78,6 +79,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,allocatable :: c(:,:) double precision,allocatable :: cp(:,:) double precision,allocatable :: eGW(:) + double precision,allocatable :: eOld(:) double precision,allocatable :: P(:,:) double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) @@ -136,9 +138,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Memory allocation - allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & - OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & + J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & + OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -149,6 +151,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE Conv = 1d0 P(:,:) = PHF(:,:) eGW(:) = eHF(:) + eOld(:) = eHF(:) c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 @@ -181,6 +184,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, & OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,OmRPA) endif @@ -211,21 +215,18 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:) - call matout(nBas,nBAs,SigCp) - ! Compute commutator and convergence criteria error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) - Conv = maxval(abs(error)) ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) - -! Reset DIIS if required - - if(abs(rcond) < 1d-15) n_diis = 0 + if(abs(rcond) > 1d-7) then + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + else + n_diis = 0 + end if ! Diagonalize Hamiltonian in AO basis @@ -233,6 +234,12 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE cp(:,:) = Fp(:,:) call diagonalize_matrix(nBas,cp,eGW) c = matmul(X,cp) + SigCp = matmul(transpose(c),matmul(SigCp,c)) + + ! Save quasiparticles energy for next cycle + + Conv = maxval(abs(eGW - eOld)) + eOld(:) = eGW(:) ! Compute new density matrix in the AO basis @@ -265,7 +272,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Print results call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,T,V,J,K,F,SigCp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole) + call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole) enddo !------------------------------------------------------------------------ diff --git a/src/MBPT/qsUGW.f90 b/src/MBPT/qsUGW.f90 index de27201..7e7b051 100644 --- a/src/MBPT/qsUGW.f90 +++ b/src/MBPT/qsUGW.f90 @@ -92,6 +92,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS double precision,allocatable :: c(:,:,:) double precision,allocatable :: cp(:,:,:) double precision,allocatable :: eGW(:,:) + double precision,allocatable :: eOld(:,:) double precision,allocatable :: P(:,:,:) double precision,allocatable :: F(:,:,:) double precision,allocatable :: Fp(:,:,:) @@ -154,10 +155,11 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS nS_bb = nS(2) nS_sc = nS_aa + nS_bb - allocate(eGW(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), & - J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin),SigCm(nBas,nBas,nspin), & - Z(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(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) + allocate(eGW(nBas,nspin),eOld(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), & + Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), & + SigCm(nBas,nBas,nspin),Z(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(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), & + F_diis(nBasSq,max_diis,nspin)) ! Initialization @@ -167,6 +169,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS Conv = 1d0 P(:,:,:) = PHF(:,:,:) eGW(:,:) = eHF(:,:) + eOld(:,:) = eHF(:,:) c(:,:,:) = cHF(:,:,:) F_diis(:,:,:) = 0d0 error_diis(:,:,:) = 0d0 @@ -268,14 +271,14 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - do is=1,nspin - if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,error_diis(:,1:n_diis,is), & - F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is)) - end do - - ! Reset DIIS if required - - if(minval(rcond(:)) < 1d-15) n_diis = 0 + if(minval(rcond(:)) > 1d-7) then + do is=1,nspin + if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,error_diis(:,1:n_diis,is), & + F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is)) + end do + else + n_diis = 0 + end if ! Transform Fock matrix in orthogonal basis @@ -296,12 +299,23 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS c(:,:,is) = matmul(X(:,:),cp(:,:,is)) end do + ! Back-transform self-energy + + do is=1,nspin + SigCp(:,:,is) = matmul(transpose(c(:,:,is)),matmul(SigCp(:,:,is),c(:,:,is))) + end do + ! Compute density matrix do is=1,nspin P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is))) end do + ! Save quasiparticles energy for next cycle + + Conv = maxval(abs(eGW(:,:) - eOld(:,:))) + eOld(:,:) = eGW(:,:) + !------------------------------------------------------------------------ ! Compute total energy !------------------------------------------------------------------------