PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

experimenting with QP_roots

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-02-12 21:44:35 +01:00
parent b09cff2ff7
commit c71ccda4df
5 changed files with 94 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.N2
View File

@ -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

36
input/basis
View File

@ -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

36
input/weight
View File

@ -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

8
src/QuAcK/G0W0.f90
View File

@ -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))

56
src/QuAcK/QP_roots.f90
View File

@ -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