From 2b1b2096c424f468638ea53f532bb6e9a1c11684 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 23 Sep 2020 14:23:40 +0200 Subject: [PATCH] graph sol for UGW --- input/methods | 4 +- input/options | 4 +- src/QuAcK/SigmaC.f90 | 26 ++---- src/QuAcK/UG0W0.f90 | 13 +-- src/QuAcK/USigmaC.f90 | 48 +++++++++++ src/QuAcK/dUSigmaC.f90 | 50 +++++++++++ src/QuAcK/unrestricted_Bethe_Salpeter.f90 | 2 + .../unrestricted_Bethe_Salpeter_A_matrix.f90 | 19 ++--- .../unrestricted_Bethe_Salpeter_B_matrix.f90 | 8 +- src/QuAcK/unrestricted_QP_graph.f90 | 83 +++++++++++++++++++ .../unrestricted_renormalization_factor.f90 | 5 +- ...estricted_self_energy_correlation_diag.f90 | 4 +- 12 files changed, 218 insertions(+), 48 deletions(-) create mode 100644 src/QuAcK/USigmaC.f90 create mode 100644 src/QuAcK/dUSigmaC.f90 create mode 100644 src/QuAcK/unrestricted_QP_graph.f90 diff --git a/input/methods b/input/methods index 022df90..ba17592 100644 --- a/input/methods +++ b/input/methods @@ -9,11 +9,11 @@ # CIS CID CISD F F F # RPA RPAx ppRPA - T F F + F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0 evGW qsGW - F F F + T F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 23ac893..d55ed53 100644 --- a/input/options +++ b/input/options @@ -5,11 +5,11 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: singlet triplet spin_conserved spin_flip TDA - T T T T F + T T T F F # 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 - 256 0.00001 T 5 T 0.001 F F F F F + 256 0.00001 T 5 F 0.0 F F F F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn diff --git a/src/QuAcK/SigmaC.f90 b/src/QuAcK/SigmaC.f90 index 32a07ca..10eddc2 100644 --- a/src/QuAcK/SigmaC.f90 +++ b/src/QuAcK/SigmaC.f90 @@ -1,4 +1,4 @@ -double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) +double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) ! Compute diagonal of the correlation part of the self-energy @@ -7,7 +7,7 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) ! Input variables - integer,intent(in) :: x + integer,intent(in) :: p double precision,intent(in) :: w double precision,intent(in) :: eta integer,intent(in) :: nBas @@ -22,7 +22,7 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) ! Local variables - integer :: i,j,a,b,p,jb + integer :: i,a,jb double precision :: eps ! Initialize @@ -32,26 +32,18 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) ! Occupied part of the correlation self-energy do i=nC+1,nO - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = w - e(i) + Omega(jb) - SigmaC = SigmaC + 2d0*rho(x,i,jb)**2*eps/(eps**2 + eta**2) - enddo + do jb=1,nS + eps = w - e(i) + Omega(jb) + SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2) enddo enddo ! Virtual part of the correlation self-energy do a=nO+1,nBas-nR - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = w - e(a) - Omega(jb) - SigmaC = SigmaC + 2d0*rho(x,a,jb)**2*eps/(eps**2 + eta**2) - enddo + do jb=1,nS + eps = w - e(a) - Omega(jb) + SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2) enddo enddo diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index 453f326..ac808cd 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -45,6 +45,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Local variables logical :: print_W = .true. + integer :: is integer :: ispin double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) @@ -126,13 +127,13 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Compute self-energy ! !---------------------! - call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,eHF,Omega_sc,rho_sc,SigC) + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Omega_sc,rho_sc,SigC) !--------------------------------! ! Compute renormalization factor ! !--------------------------------! - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,eHF,Omega_sc,rho_sc,Z) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,Omega_sc,rho_sc,Z) !-----------------------------------! ! Solve the quasi-particle equation ! @@ -151,10 +152,10 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Find graphical solution of the QP equation -! do is=1,nspin -! call QP_graph(nBas,nC(:,is),nO(:,is),nV(:,is),nR(:,is),nS(:,is),eta,eHF(:,is),Omega(:,is), & -! rho(:,:,:,ispin),eGWlin(:,is),eGW(:,is)) -! end do + do is=1,nspin + call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),Omega_sc, & + rho_sc,eGWlin(:,is),eGW(:,is)) + end do end if diff --git a/src/QuAcK/USigmaC.f90 b/src/QuAcK/USigmaC.f90 new file mode 100644 index 0000000..d701cd7 --- /dev/null +++ b/src/QuAcK/USigmaC.f90 @@ -0,0 +1,48 @@ +double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + +! Local variables + + integer :: i,a,jb + double precision :: eps + +! Initialize + + USigmaC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do jb=1,nS + eps = w - e(i) + Omega(jb) + USigmaC = uSigmaC + rho(p,i,jb,1)**2*eps/(eps**2 + eta**2) + end do + end do + + do a=nO+1,nBas-nR + do jb=1,nS + eps = w - e(a) - Omega(jb) + USigmaC = USigmaC + rho(p,a,jb,1)**2*eps/(eps**2 + eta**2) + end do + end do + +end function USigmaC diff --git a/src/QuAcK/dUSigmaC.f90 b/src/QuAcK/dUSigmaC.f90 new file mode 100644 index 0000000..cbd0cdb --- /dev/null +++ b/src/QuAcK/dUSigmaC.f90 @@ -0,0 +1,50 @@ +double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) + +! Compute the derivative of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + +! Local variables + + integer :: i,a,jb + double precision :: eps + +! Initialize + + dUSigmaC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do jb=1,nS + eps = w - e(i) + Omega(jb) + dUSigmaC = dUSigmaC + rho(p,i,jb,1)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + +! Virtual part of the correlation self-energy + + do a=nO+1,nBas-nR + do jb=1,nS + eps = w - e(a) - Omega(jb) + dUSigmaC = dUSigmaC + rho(p,a,jb,1)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + +end function dUSigmaC diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 index 38b36ee..42e7d6d 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 @@ -76,6 +76,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcRPA(ispin), & OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc) +! call print_excitation('RPA@UG0W0',5,nS_sc,OmRPA_sc) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb, & XpY_RPA_sc,rho_RPA_sc) diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 index e3f5813..175704a 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 @@ -1,6 +1,5 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & - ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & - Omega,rho,A_lr) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism @@ -40,7 +39,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt ! Build part A of the BSE matrix ! !--------------------------------! - ! alpha-alpha block + ! aaaa block ia = 0 do i=nC(1)+1,nO(1) @@ -55,7 +54,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt do kc=1,nSt eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps & - + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps + + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps enddo A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi @@ -65,7 +64,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt enddo enddo - ! alpha-beta block + ! aabb block ia = 0 do i=nC(1)+1,nO(1) @@ -90,7 +89,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt enddo enddo - ! beta-alpha block + ! bbaa block ia = 0 do i=nC(2)+1,nO(2) @@ -104,8 +103,8 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt chi = 0d0 do kc=1,nSt eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps & - + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps + chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps & + + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps enddo A_lr(nSa+ia,jb) = A_lr(nSa+ia,jb) - lambda*ERI_aabb(b,i,a,j) + 2d0*lambda*chi @@ -115,7 +114,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt enddo enddo - ! beta-beta block + ! bbbb block ia = 0 do i=nC(2)+1,nO(2) @@ -129,7 +128,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt chi = 0d0 do kc=1,nSt eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps & + chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps & + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps enddo diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 index 250978d..20dd14e 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 @@ -51,7 +51,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt do kc=1,nSt eps = Omega(kc)**2 + eta**2 chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps & - + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps + + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps enddo B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi @@ -100,8 +100,8 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt chi = 0d0 do kc=1,nSt eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps & - + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps + chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps & + + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps enddo B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_aabb(j,i,a,b) + 2d0*lambda*chi @@ -125,7 +125,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt chi = 0d0 do kc=1,nSt eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps & + chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps & + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps enddo diff --git a/src/QuAcK/unrestricted_QP_graph.f90 b/src/QuAcK/unrestricted_QP_graph.f90 new file mode 100644 index 0000000..96585a9 --- /dev/null +++ b/src/QuAcK/unrestricted_QP_graph.f90 @@ -0,0 +1,83 @@ +subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) + +! Compute the graphical solution of the QP equation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: eta + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + + double precision,intent(in) :: eGWlin(nBas) + +! Local variables + + integer :: p + integer :: nIt + integer,parameter :: maxIt = 10 + double precision,parameter :: thresh = 1d-6 + double precision,external :: USigmaC,dUSigmaC + double precision :: sig,dsig + double precision :: f,df + double precision :: w + +! Output variables + + double precision,intent(out) :: eGW(nBas) + +! Run Newton's algorithm to find the root + + do p=nC+1,nBas-nR + + write(*,*) '-----------------' + write(*,'(A10,I3)') 'Orbital ',p + write(*,*) '-----------------' + + w = eGWlin(p) + nIt = 0 + f = 1d0 + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f + + do while (abs(f) > thresh .and. nIt < maxIt) + + nIt = nIt + 1 + + sig = USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) + dsig = dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) + f = w - eHF(p) - sig + df = 1d0 - dsig + + w = w - f/df + + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sig + + + end do + + if(nIt == maxIt) then + + write(*,*) 'Newton root search has not converged!' + eGW(p) = eGWlin(p) + + else + + eGW(p) = w + + write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV + write(*,*) + + end if + + end do + +end subroutine unrestricted_QP_graph diff --git a/src/QuAcK/unrestricted_renormalization_factor.f90 b/src/QuAcK/unrestricted_renormalization_factor.f90 index 403060d..f175d93 100644 --- a/src/QuAcK/unrestricted_renormalization_factor.f90 +++ b/src/QuAcK/unrestricted_renormalization_factor.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,e,Omega,rho,Z) +subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,Z) ! Compute the renormalization factor in the unrestricted formalism @@ -13,8 +13,6 @@ subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt, integer,intent(in) :: nO(nspin) integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) - integer,intent(in) :: nSa - integer,intent(in) :: nSb integer,intent(in) :: nSt double precision,intent(in) :: e(nBas,nspin) double precision,intent(in) :: Omega(nSt) @@ -89,5 +87,4 @@ subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt, Z(:,:) = 1d0/(1d0 + Z(:,:)) - end subroutine unrestricted_renormalization_factor diff --git a/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 b/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 index 123ec23..fa4084c 100644 --- a/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 +++ b/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,e,Omega,rho,SigC) +subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC) ! Compute diagonal of the correlation part of the self-energy @@ -13,8 +13,6 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS integer,intent(in) :: nO(nspin) integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) - integer,intent(in) :: nSa - integer,intent(in) :: nSb integer,intent(in) :: nSt double precision,intent(in) :: e(nBas,nspin) double precision,intent(in) :: Omega(nSt)