From c71ccda4df808774b9a754400b2cdb1ed8f1236f Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 12 Feb 2020 21:44:35 +0100 Subject: [PATCH] experimenting with QP_roots --- examples/molecule.N2 | 2 +- input/basis | 36 ++++++++++++++++++--------- input/weight | 36 ++++++++++++++++++--------- src/QuAcK/G0W0.f90 | 8 +++--- src/QuAcK/QP_roots.f90 | 56 +++++++++++++++++++++++++++++++----------- 5 files changed, 94 insertions(+), 44 deletions(-) diff --git a/examples/molecule.N2 b/examples/molecule.N2 index 85fe0df..35e31a1 100644 --- a/examples/molecule.N2 +++ b/examples/molecule.N2 @@ -2,4 +2,4 @@ 2 7 7 0 0 # Znuc x y z N 0. 0. 0. - N 0. 0. 2.072 + N 0. 0. 2.0 diff --git a/input/basis b/input/basis index fb05e68..69b24c8 100644 --- a/input/basis +++ b/input/basis @@ -1,18 +1,30 @@ -1 3 +1 6 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 +S 1 + 1 0.1027000 1.0000000 P 1 - 1 0.7270000 1.0000000 -2 3 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +2 6 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 +S 1 + 1 0.1027000 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 +D 1 + 1 1.0570000 1.0000000 diff --git a/input/weight b/input/weight index fb05e68..69b24c8 100644 --- a/input/weight +++ b/input/weight @@ -1,18 +1,30 @@ -1 3 +1 6 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 +S 1 + 1 0.1027000 1.0000000 P 1 - 1 0.7270000 1.0000000 -2 3 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +2 6 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 +S 1 + 1 0.1027000 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 +D 1 + 1 1.0570000 1.0000000 diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index cd20734..d2d55c1 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -95,14 +95,14 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) -! Find all the roots of the QP equation if necessary - - call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho) - ! Solve the quasi-particle equation eGW(:) = eHF(:) + Z(:)*SigC(:) +! Find all the roots of the QP equation if necessary + + call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGW) + ! Dump results ! call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) diff --git a/src/QuAcK/QP_roots.f90 b/src/QuAcK/QP_roots.f90 index 2dca6fe..fc1cde5 100644 --- a/src/QuAcK/QP_roots.f90 +++ b/src/QuAcK/QP_roots.f90 @@ -1,4 +1,4 @@ -subroutine QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho) +subroutine QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGW) ! Compute all the roots of the QP equation for each orbital @@ -21,10 +21,9 @@ subroutine QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho) ! Local variables integer :: i,j,a,b,x,jb,g - integer :: nRoot - integer,parameter :: nGrid = 10000 - double precision,parameter :: wmin = -10d0 - double precision,parameter :: wmax = +10d0 + integer,parameter :: nGrid = 100000 + double precision,parameter :: wmin = -50d0 + double precision,parameter :: wmax = +50d0 double precision :: dw double precision :: eps double precision :: left_sign,right_sign @@ -37,9 +36,17 @@ subroutine QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho) integer,parameter :: maxIt = 64 double precision,parameter :: thresh = 1d-6 double precision,external :: SigmaC,dSigmaC - double precision :: f,df - double precision :: s,ds + double precision :: sig,dsig + double precision :: sat,dsat double precision :: om + integer :: nRoot + double precision :: Z + double precision :: Ztot + double precision :: QP + +! Output variables + + double precision,intent(out) :: eGW(nBas) ! Construct grid @@ -101,6 +108,8 @@ subroutine QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho) ! Find the zeros of the QP equation nRoot = 0 + QP = 0d0 + Ztot = 0d0 write(*,*) '-----------------' write(*,'(A10,I3)') 'Orbital ',x @@ -124,31 +133,48 @@ subroutine QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho) om = w(g-1) nIt = 0 - f = 1d0 + sat = 1d0 - do while (abs(f) > thresh .and. nIt < maxIt) + do while (abs(sat) > thresh .and. nIt < maxIt) nIt = nIt + 1 - s = SigmaC(x,om,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) - ds = dSigmaC(x,om,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) - f = om - eHF(x) - s - df = 1d0 - ds + sig = SigmaC(x,om,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) + dsig = dSigmaC(x,om,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) + sat = om - eHF(x) - sig + dsat = 1d0 - dsig + Z = 1d0/dsat - write(*,'(A3,I3,A1,1X,4F15.9)') 'It.',nIt,':',om,f,s,1d0/(1d0-ds) + write(*,'(A3,I3,A1,1X,4F15.9)') 'It.',nIt,':',om,sat,sig,Z - om = om - f/df + om = om - sat/dsat end do + + if(nIt == maxIt) then + write(*,*) 'Newton root search has not converged!' + + else + + Ztot = Ztot + Z + QP = QP + Z*om + + end if + end if left_sign = right_sign left_QP = right_QP end do + + eGW(x) = QP/Ztot + write(*,*) write(*,'(A32,I3,A1,I3)') 'Number of roots for orbital ',x,':',nRoot + write(*,'(A32,F16.10)') 'Total weight ',Ztot + write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(x)*HaToeV write(*,*) end do