10
1
mirror of https://github.com/pfloos/quack synced 2024-11-05 05:33:50 +01:00

graphical GW

This commit is contained in:
Pierre-Francois Loos 2020-03-12 15:04:16 +01:00
parent c8855d8aed
commit fc73914faa
10 changed files with 286 additions and 87 deletions

View File

@ -1,46 +1,98 @@
1 10 1 15
S 3 S 9
1 82.6400000 0.0020060 1 33980.0000000 0.0000910
2 12.4100000 0.0153430 2 5089.0000000 0.0007040
3 2.8240000 0.0755790 3 1157.0000000 0.0036930
4 326.6000000 0.0153600
5 106.1000000 0.0529290
6 38.1100000 0.1470430
7 14.7500000 0.3056310
8 6.0350000 0.3993450
9 2.5300000 0.2170510
S 9
1 33980.0000000 -0.0000190
2 5089.0000000 -0.0001510
3 1157.0000000 -0.0007850
4 326.6000000 -0.0033240
5 106.1000000 -0.0115120
6 38.1100000 -0.0341600
7 14.7500000 -0.0771730
8 6.0350000 -0.1414930
9 2.5300000 -0.1180190
S 1 S 1
1 0.7977000 1.0000000 1 0.7355000 1.0000000
S 1 S 1
1 0.2581000 1.0000000 1 0.2905000 1.0000000
S 1 S 1
1 0.0898900 1.0000000 1 0.1111000 1.0000000
P 3
1 34.5100000 0.0053780
2 7.9150000 0.0361320
3 2.3680000 0.1424930
P 1 P 1
1 2.2920000 1.0000000 1 0.8132000 1.0000000
P 1 P 1
1 0.8380000 1.0000000 1 0.2890000 1.0000000
P 1 P 1
1 0.2920000 1.0000000 1 0.1007000 1.0000000
D 1 D 1
1 2.0620000 1.0000000 1 1.8480000 1.0000000
D 1 D 1
1 0.6620000 1.0000000 1 0.6490000 1.0000000
D 1
1 0.2280000 1.0000000
F 1 F 1
1 1.3970000 1.0000000 1 1.4190000 1.0000000
2 10
S 3
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.7977000 1.0000000
S 1
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
P 1
1 2.2920000 1.0000000
P 1
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
F 1 F 1
1 1.3970000 1.0000000 1 0.4850000 1.0000000
G 1
1 1.0110000 1.0000000
2 15
S 9
1 61420.0000000 0.0000900
2 9199.0000000 0.0006980
3 2091.0000000 0.0036640
4 590.9000000 0.0152180
5 192.3000000 0.0524230
6 69.3200000 0.1459210
7 26.9700000 0.3052580
8 11.1000000 0.3985080
9 4.6820000 0.2169800
S 9
1 61420.0000000 -0.0000200
2 9199.0000000 -0.0001590
3 2091.0000000 -0.0008290
4 590.9000000 -0.0035080
5 192.3000000 -0.0121560
6 69.3200000 -0.0362610
7 26.9700000 -0.0829920
8 11.1000000 -0.1520900
9 4.6820000 -0.1153310
S 1
1 1.4280000 1.0000000
S 1
1 0.5547000 1.0000000
S 1
1 0.2067000 1.0000000
P 3
1 63.4200000 0.0060440
2 14.6600000 0.0417990
3 4.4590000 0.1611430
P 1
1 1.5310000 1.0000000
P 1
1 0.5302000 1.0000000
P 1
1 0.1750000 1.0000000
D 1
1 3.7750000 1.0000000
D 1
1 1.3000000 1.0000000
D 1
1 0.4440000 1.0000000
F 1
1 2.6660000 1.0000000
F 1
1 0.8590000 1.0000000
G 1
1 1.8460000 1.0000000

View File

@ -5,7 +5,7 @@
# CCD CCSD CCSD(T) ringCCD ladderCCD # CCD CCSD CCSD(T) ringCCD ladderCCD
F F F F F F F F F F
# CIS RPA RPAx ppRPA ADC # CIS RPA RPAx ppRPA ADC
F T T F F F F F F F
# GF2 GF3 # GF2 GF3
F F F F
# G0W0 evGW qsGW # G0W0 evGW qsGW

View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
2 1 1 0 0 2 7 7 0 0
# Znuc x y z # Znuc x y z
H 0. 0. 0. C 0. 0. 0.
H 0. 0. 1.399 O 0. 0. 2.134

View File

@ -1,4 +1,4 @@
2 2
H 0.0000000000 0.0000000000 0.0000000000 C 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 0.0000000000 0.7403189714 O 0.0000000000 0.0000000000 1.1292642494

View File

@ -11,6 +11,6 @@
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 5 F F T F F F F 0.000 256 0.00001 T 5 F F T F F F F 0.000
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
T F F T F T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T 1000000 100000 10 0.3 10000 1234 T

View File

@ -1,46 +1,98 @@
1 10 1 15
S 3 S 9
1 82.6400000 0.0020060 1 33980.0000000 0.0000910
2 12.4100000 0.0153430 2 5089.0000000 0.0007040
3 2.8240000 0.0755790 3 1157.0000000 0.0036930
4 326.6000000 0.0153600
5 106.1000000 0.0529290
6 38.1100000 0.1470430
7 14.7500000 0.3056310
8 6.0350000 0.3993450
9 2.5300000 0.2170510
S 9
1 33980.0000000 -0.0000190
2 5089.0000000 -0.0001510
3 1157.0000000 -0.0007850
4 326.6000000 -0.0033240
5 106.1000000 -0.0115120
6 38.1100000 -0.0341600
7 14.7500000 -0.0771730
8 6.0350000 -0.1414930
9 2.5300000 -0.1180190
S 1 S 1
1 0.7977000 1.0000000 1 0.7355000 1.0000000
S 1 S 1
1 0.2581000 1.0000000 1 0.2905000 1.0000000
S 1 S 1
1 0.0898900 1.0000000 1 0.1111000 1.0000000
P 3
1 34.5100000 0.0053780
2 7.9150000 0.0361320
3 2.3680000 0.1424930
P 1 P 1
1 2.2920000 1.0000000 1 0.8132000 1.0000000
P 1 P 1
1 0.8380000 1.0000000 1 0.2890000 1.0000000
P 1 P 1
1 0.2920000 1.0000000 1 0.1007000 1.0000000
D 1 D 1
1 2.0620000 1.0000000 1 1.8480000 1.0000000
D 1 D 1
1 0.6620000 1.0000000 1 0.6490000 1.0000000
D 1
1 0.2280000 1.0000000
F 1 F 1
1 1.3970000 1.0000000 1 1.4190000 1.0000000
2 10
S 3
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.7977000 1.0000000
S 1
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
P 1
1 2.2920000 1.0000000
P 1
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
F 1 F 1
1 1.3970000 1.0000000 1 0.4850000 1.0000000
G 1
1 1.0110000 1.0000000
2 15
S 9
1 61420.0000000 0.0000900
2 9199.0000000 0.0006980
3 2091.0000000 0.0036640
4 590.9000000 0.0152180
5 192.3000000 0.0524230
6 69.3200000 0.1459210
7 26.9700000 0.3052580
8 11.1000000 0.3985080
9 4.6820000 0.2169800
S 9
1 61420.0000000 -0.0000200
2 9199.0000000 -0.0001590
3 2091.0000000 -0.0008290
4 590.9000000 -0.0035080
5 192.3000000 -0.0121560
6 69.3200000 -0.0362610
7 26.9700000 -0.0829920
8 11.1000000 -0.1520900
9 4.6820000 -0.1153310
S 1
1 1.4280000 1.0000000
S 1
1 0.5547000 1.0000000
S 1
1 0.2067000 1.0000000
P 3
1 63.4200000 0.0060440
2 14.6600000 0.0417990
3 4.4590000 0.1611430
P 1
1 1.5310000 1.0000000
P 1
1 0.5302000 1.0000000
P 1
1 0.1750000 1.0000000
D 1
1 3.7750000 1.0000000
D 1
1 1.3000000 1.0000000
D 1
1 0.4440000 1.0000000
F 1
1 2.6660000 1.0000000
F 1
1 0.8590000 1.0000000
G 1
1 1.8460000 1.0000000

View File

@ -39,9 +39,9 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
ispin = 1 ispin = 1
EcBSE(ispin) = 0d0 EcBSE(ispin) = 0d0
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & ! call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) ! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) ! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))

View File

@ -45,6 +45,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:) double precision,allocatable :: rhox(:,:,:,:)
double precision,allocatable :: eGWlin(:)
! Output variables ! Output variables
double precision :: eGW(nBas) double precision :: eGW(nBas)
@ -74,7 +76,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
! Memory allocation ! Memory allocation
allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), & allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin)) rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin),eGWlin(nBas))
! Compute linear response ! Compute linear response
@ -97,11 +99,17 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
! Solve the quasi-particle equation ! Solve the quasi-particle equation
eGW(:) = eHF(:) + Z(:)*SigC(:) eGWlin(:) = eHF(:) + Z(:)*SigC(:)
! Find all the roots of the QP equation if necessary ! Find all the roots of the QP equation if necessary
! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGW) ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin)
! Find graphical solution of the QP equation
! call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW)
eGW(:) = eGWlin(:)
! Dump results ! Dump results

82
src/QuAcK/QP_graph.f90 Normal file
View File

@ -0,0 +1,82 @@
subroutine 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)
double precision,intent(in) :: eGWlin(nBas)
! Local variables
integer :: p
integer :: nIt
integer,parameter :: maxIt = 64
double precision,parameter :: thresh = 1d-6
double precision,external :: SigmaC,dSigmaC
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 = SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho)
dsig = dSigmaC(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!'
else
eGW(p) = w
write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV
write(*,*)
end if
end do
end subroutine QP_graph

View File

@ -280,6 +280,11 @@ program QuAcK
call cpu_time(start_AOtoMO) call cpu_time(start_AOtoMO)
write(*,*)
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
if(doSph) then if(doSph) then
ERI_MO_basis = ERI_AO_basis ERI_MO_basis = ERI_AO_basis
@ -542,8 +547,8 @@ program QuAcK
if(doG0W0) then if(doG0W0) then
call cpu_time(start_G0W0) call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0) eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0) call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0 t_G0W0 = end_G0W0 - start_G0W0