LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:50 +01:00
Code Releases Activity

clean up GF

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-06-03 12:06:16 +02:00
parent ea78befb4f
commit 4b3bb88b84
14 changed files with 223 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.H2
View File

@ -2,4 +2,4 @@
2 1 1 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 1.435
H 0. 0. 10.1

4
examples/molecule.H2O
View File

@ -1,10 +1,6 @@
# nAt nEla nElb nCore nRyd
3 5 5 0 0
# Znuc x y z
O 0.0 0.0 0.226200217
H 0.0 1.43921542 -0.904422923
H 0.0 -1.43921542 -0.904422923
O 0. 0. -0.13209669
H 0. 1.43152878 0.97970006
H 0. -1.43152878 0.97970006

121
input/basis
View File

@ -1,73 +1,64 @@
1 15
1 6
S 9
1 33980.0000000 0.0000910
2 5089.0000000 0.0007040
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
1 9.046000E+03 7.000000E-04
2 1.357000E+03 5.389000E-03
3 3.093000E+02 2.740600E-02
4 8.773000E+01 1.032070E-01
5 2.856000E+01 2.787230E-01
6 1.021000E+01 4.485400E-01
7 3.838000E+00 2.782380E-01
8 7.466000E-01 1.544000E-02
9 2.248000E-01 -2.864000E-03
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
1 9.046000E+03 -1.530000E-04
2 1.357000E+03 -1.208000E-03
3 3.093000E+02 -5.992000E-03
4 8.773000E+01 -2.454400E-02
5 2.856000E+01 -6.745900E-02
6 1.021000E+01 -1.580780E-01
7 3.838000E+00 -1.218310E-01
8 7.466000E-01 5.490030E-01
9 2.248000E-01 5.788150E-01
S 1
1 0.7355000 1.0000000
S 1
1 0.2905000 1.0000000
S 1
1 0.1111000 1.0000000
P 3
1 34.5100000 0.0053780
2 7.9150000 0.0361320
3 2.3680000 0.1424930
1 2.248000E-01 1.000000E+00
P 4
1 1.355000E+01 3.991900E-02
2 2.917000E+00 2.171690E-01
3 7.973000E-01 5.103190E-01
4 2.185000E-01 4.622140E-01
P 1
1 0.8132000 1.0000000
P 1
1 0.2890000 1.0000000
P 1
1 0.1007000 1.0000000
1 2.185000E-01 1.000000E+00
D 1
1 1.8480000 1.0000000
D 1
1 0.6490000 1.0000000
D 1
1 0.2280000 1.0000000
F 1
1 1.4190000 1.0000000
F 1
1 0.4850000 1.0000000
G 1
1 1.0110000 1.0000000
2 10
S 3
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
1 8.170000E-01 1.0000000
2 6
S 9
1 9.046000E+03 7.000000E-04
2 1.357000E+03 5.389000E-03
3 3.093000E+02 2.740600E-02
4 8.773000E+01 1.032070E-01
5 2.856000E+01 2.787230E-01
6 1.021000E+01 4.485400E-01
7 3.838000E+00 2.782380E-01
8 7.466000E-01 1.544000E-02
9 2.248000E-01 -2.864000E-03
S 9
1 9.046000E+03 -1.530000E-04
2 1.357000E+03 -1.208000E-03
3 3.093000E+02 -5.992000E-03
4 8.773000E+01 -2.454400E-02
5 2.856000E+01 -6.745900E-02
6 1.021000E+01 -1.580780E-01
7 3.838000E+00 -1.218310E-01
8 7.466000E-01 5.490030E-01
9 2.248000E-01 5.788150E-01
S 1
1 0.7977000 1.0000000
S 1
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
1 2.248000E-01 1.000000E+00
P 4
1 1.355000E+01 3.991900E-02
2 2.917000E+00 2.171690E-01
3 7.973000E-01 5.103190E-01
4 2.185000E-01 4.622140E-01
P 1
1 2.2920000 1.0000000
P 1
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
1 2.185000E-01 1.000000E+00
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
F 1
1 1.3970000 1.0000000
1 8.170000E-01 1.0000000

4
input/methods
View File

@ -7,9 +7,9 @@
# drCCD rCCD lCCD pCCD
F F F F
# CIS CID CISD
T F F
F F F
# RPA RPAx ppRPA
F T F
F F F
# G0F2 evGF2 G0F3 evGF3
T F F F
# G0W0 evGW qsGW

6
input/molecule
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd
2 3 3 0 0
2 7 7 0 0
# Znuc x y z
C 0. 0. 0.
H 0. 0. 2.13709127
N 0. 0. -1.04008632
N 0. 0. +1.04008632

4
input/molecule.xyz
View File

@ -1,4 +1,4 @@
2
C 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 0.0000000000 1.1309000791
N 0.0000000000 0.0000000000 -0.5503900175
N 0.0000000000 0.0000000000 0.5503900175

6
input/options
View File

@ -6,10 +6,10 @@
64 0.0000001 T 5
# spin: singlet triplet
T T
# GF: maxSCF thresh DIIS n_diis lin renorm
256 0.00001 T 5 T 3
# GF: maxSCF thresh DIIS n_diis lin renorm BSE TDA eta
256 0.00001 T 5 T 3 T T 0.00367493
# GW/GT: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 10 F F T T F F T 0.0
256 0.00001 T 5 F F T F F F T 0.00367493
# ACFDT: AC Kx XBS
F F T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift

8
src/QuAcK/BSE2_dynamic_perturbation.f90
View File

@ -47,11 +47,17 @@ subroutine BSE2_dynamic_perturbation(ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,O
if(.not.dTDA) allocate(B_dyn(nS,nS),ZB_dyn(nS,nS))
! Print main components of transition vectors
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY)
gapGF = eGF(nO+1) - eGF(nO)
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,*) ' First-order dynamical correction to static 2nd-order Bethe-Salpeter excitation energies '
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,'(A58,F10.6,A3)') ' BSE neutral excitation must be lower than the GF2 gap = ',gapGF*HaToeV,' eV'
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
write(*,*) '---------------------------------------------------------------------------------------------------'
@ -95,8 +101,6 @@ subroutine BSE2_dynamic_perturbation(ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,O
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') &
ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,ZDyn(ia)
if(OmBSE(ia) > gapGF) write(*,*) ' !!! BSE2 neutral excitation larger than the GF2 gap !!! '
end do
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,*)

11
src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90
View File

@ -25,8 +25,9 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O
! Local variables
logical :: dTDA = .true.
integer :: ia
logical :: dTDA = .true.
integer,parameter :: maxS = 10
double precision :: gapGW
@ -47,11 +48,17 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O
if(.not.dTDA) allocate(B_dyn(nS,nS),ZB_dyn(nS,nS))
! Print main components of transition vectors
call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY)
gapGW = eGW(nO+1) - eGW(nO)
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies '
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV'
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
write(*,*) '---------------------------------------------------------------------------------------------------'
@ -107,8 +114,6 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') &
ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,ZDyn(ia)
if(OmBSE(ia) > gapGW) write(*,*) ' !!! BSE neutral excitation larger than the GW gap !!! '
end do
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,*)

10
src/QuAcK/G0F2.f90
View File

@ -66,8 +66,8 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,
eps = eHF(p) + eHF(a) - eHF(i) - eHF(j)
V = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
Sig(p) = Sig(p) + V/eps
Z(p) = Z(p) + V/eps**2
Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2)
Z(p) = Z(p) - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
@ -81,15 +81,15 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,
eps = eHF(p) + eHF(i) - eHF(a) - eHF(b)
V = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
Sig(p) = Sig(p) + V/eps
Z(p) = Z(p) + V/eps**2
Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2)
Z(p) = Z(p) - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
end do
end do
Z(:) = 1d0/(1d0 + Z(:))
Z(:) = 1d0/(1d0 - Z(:))
if(linearize) then

51
src/QuAcK/QuAcK.f90
View File

@ -96,12 +96,13 @@ program QuAcK
integer :: maxSCF_GF,n_diis_GF,renormGF
double precision :: thresh_GF
logical :: DIIS_GF,linGF
logical :: DIIS_GF,linGF,BSE_GF,TDA_GF
double precision :: eta_GF
integer :: maxSCF_GW,n_diis_GW
double precision :: thresh_GW
logical :: DIIS_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linGW
double precision :: eta
logical :: DIIS_GW,COHSEX,SOSEX,BSE_GW,TDA_GW,G0W,GW0,linGW
double precision :: eta_GW
integer :: nMC,nEq,nWalk,nPrint,iSeed
double precision :: dt
@ -144,8 +145,9 @@ program QuAcK
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,renormGF, &
BSE_GF,TDA_GF,eta_GF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, &
COHSEX,SOSEX,BSE,TDA,G0W,GW0,linGW,eta, &
COHSEX,SOSEX,BSE_GW,TDA_GW,G0W,GW0,linGW,eta_GW, &
doACFDT,exchange_kernel,doXBS, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
@ -533,7 +535,7 @@ program QuAcK
if(doRPA) then
call cpu_time(start_RPA)
call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_RPA)
@ -550,7 +552,7 @@ program QuAcK
if(doRPAx) then
call cpu_time(start_RPAx)
call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_RPAx)
@ -601,8 +603,8 @@ program QuAcK
if(doG0F2) then
call cpu_time(start_GF2)
call G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linGF, &
eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF)
call G0F2(BSE_GF,TDA_GF,singlet_manifold,triplet_manifold,linGF, &
eta_GF,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2
@ -618,8 +620,8 @@ program QuAcK
if(doevGF2) then
call cpu_time(start_GF2)
call evGF2(BSE,TDA,maxSCF_GF,thresh_GF,n_diis_GF,singlet_manifold,triplet_manifold,linGF, &
eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF)
call evGF2(BSE_GF,TDA_GF,maxSCF_GF,thresh_GF,n_diis_GF,singlet_manifold,triplet_manifold,linGF, &
eta_GF,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2
@ -669,8 +671,8 @@ program QuAcK
if(doG0W0) then
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, &
singlet_manifold,triplet_manifold,linGW,eta, &
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, &
singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0)
@ -687,8 +689,8 @@ program QuAcK
if(doevGW) then
call cpu_time(start_evGW)
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, &
singlet_manifold,triplet_manifold,eta, &
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_evGW)
@ -705,8 +707,8 @@ program QuAcK
if(doqsGW) then
call cpu_time(start_qsGW)
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, &
singlet_manifold,triplet_manifold,eta, &
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, &
BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF)
call cpu_time(end_qsGW)
@ -725,8 +727,8 @@ program QuAcK
if(doG0T0) then
call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, &
singlet_manifold,triplet_manifold,linGW,eta, &
call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, &
singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0)
call cpu_time(end_G0T0)
@ -743,8 +745,9 @@ program QuAcK
if(doevGT) then
call cpu_time(start_evGT)
call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_manifold, &
eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0)
call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, &
BSE_GW,TDA_GW,singlet_manifold,triplet_manifold,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0)
call cpu_time(end_evGT)
t_evGT = end_evGT - start_evGT
@ -859,8 +862,8 @@ program QuAcK
! Long-range G0W0 calculation
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, &
singlet_manifold,triplet_manifold,linGW,eta, &
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, &
singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0)
@ -873,8 +876,8 @@ program QuAcK
ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:)
call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, &
singlet_manifold,triplet_manifold,linGW,eta, &
call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, &
singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0)
call cpu_time(end_G0T0)

32
src/QuAcK/evGF2.f90
View File

@ -1,5 +1,5 @@
subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold,linearize, &
eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,V,e0)
eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,V,eHF)
! Perform eigenvalue self-consistent second-order Green function calculation
@ -25,7 +25,7 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: e0(nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: V(nBas,nBas,nBas,nBas)
! Local variables
@ -63,8 +63,8 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
n_diis = 0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
eGF2(:) = e0(:)
eOld(:) = e0(:)
eGF2(:) = eHF(:)
eOld(:) = eHF(:)
!------------------------------------------------------------------------
! Main SCF loop
@ -81,10 +81,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + e0(a) - e0(i) - e0(j)
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
Sig(p) = Sig(p) &
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)*eps/(eps**2 + eta**2)
end do
end do
@ -96,10 +96,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + e0(i) - e0(a) - e0(b)
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
Sig(p) = Sig(p) &
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)*eps/(eps**2 + eta**2)
end do
end do
@ -115,9 +115,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + e0(a) - e0(i) - e0(j)
eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j)
Z(p) = Z(p) - (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps**2
Z(p) = Z(p) &
- (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
@ -129,9 +130,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + e0(i) - e0(a) - e0(b)
eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b)
Z(p) = Z(p) - (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps**2
Z(p) = Z(p) &
- (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
end do
@ -142,11 +144,11 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
if(linearize) then
eGF2(:) = e0(:) + Z(:)*Sig(:)
eGF2(:) = eHF(:) + Z(:)*Sig(:)
else
eGF2(:) = e0(:) + Sig(:)
eGF2(:) = eHF(:) + Sig(:)
end if
@ -154,7 +156,7 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol
! Print results
call print_evGF2(nBas,nO,nSCF,Conv,e0,eGF2)
call print_evGF2(nBas,nO,nSCF,Conv,eHF,eGF2)
! DIIS extrapolation

66
src/QuAcK/print_transition_vectors.f90 Normal file
View File

@ -0,0 +1,66 @@
subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY)
! Print transition vectors for linear response calculation
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) :: Omega(nS)
double precision,intent(in) :: XpY(nS,nS)
double precision,intent(in) :: XmY(nS,nS)
! Local variables
integer :: ia,jb,i,j,a,b
integer,parameter :: maxS = 10
double precision,parameter :: thres_vec = 0.1d0
double precision,allocatable :: X(:)
double precision,allocatable :: Y(:)
! Memory allocation
allocate(X(nS),Y(nS))
write(*,*)
do ia=1,min(nS,maxS)
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
print*,'--------------------------------'
write(*,'(A15,I3,A2,F10.6,A3)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV'
print*,'--------------------------------'
print*,'Resonant vectors'
print*,'--------------------------------'
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A3,F10.6)') j,' --> ',b,' = ',X(jb)
end do
end do
print*,'--------------------------------'
print*,'Anti-resonant vectors'
print*,'--------------------------------'
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-nR
jb = jb + 1
if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A3,F10.6)') j,' --> ',b,' = ',Y(jb)
end do
end do
write(*,*)
end do
end subroutine print_transition_vectors

45
src/QuAcK/read_options.f90
View File

@ -2,8 +2,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,renormGF, &
BSE_GF,TDA_GF,eta_GF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, &
COHSEX,SOSEX,BSE,TDA,G0W,GW0,linGW,eta, &
COHSEX,SOSEX,BSE_GW,TDA_GW,G0W,GW0,linGW,eta_GW, &
doACFDT,exchange_kernel,doXBS, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
@ -34,6 +35,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
integer,intent(out) :: n_diis_GF
logical,intent(out) :: linGF
integer,intent(out) :: renormGF
logical,intent(out) :: BSE_GF
logical,intent(out) :: TDA_GF
double precision,intent(out) :: eta_GF
integer,intent(out) :: maxSCF_GW
double precision,intent(out) :: thresh_GW
@ -41,12 +45,12 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
integer,intent(out) :: n_diis_GW
logical,intent(out) :: COHSEX
logical,intent(out) :: SOSEX
logical,intent(out) :: BSE
logical,intent(out) :: TDA
logical,intent(out) :: BSE_GW
logical,intent(out) :: TDA_GW
logical,intent(out) :: G0W
logical,intent(out) :: GW0
logical,intent(out) :: linGW
double precision,intent(out) :: eta
double precision,intent(out) :: eta_GW
logical,intent(out) :: doACFDT
logical,intent(out) :: exchange_kernel
@ -120,14 +124,19 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
thresh_GF = 1d-5
DIIS_GF = .false.
n_diis_GF = 5
linGF = .false.
renormGF = 0
linGF = .false.
renormGF = 0
BSE_GF = .false.
TDA_GF = .false.
eta_GF = 0d0
read(1,*)
read(1,*) maxSCF_GF,thresh_GW,answer1,n_diis_GF,answer2,renormGF
read(1,*) maxSCF_GF,thresh_GW,answer1,n_diis_GF,answer2,renormGF,answer3,answer4,eta_GF
if(answer1 == 'T') DIIS_GF = .true.
if(answer2 == 'T') linGF = .true.
if(answer3 == 'T') BSE_GF = .true.
if(answer4 == 'T') TDA_GF = .true.
if(.not.DIIS_GF) n_diis_GF = 1
! Read GW options
@ -138,25 +147,25 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
n_diis_GW = 5
COHSEX = .false.
SOSEX = .false.
BSE = .false.
TDA = .false.
BSE_GW = .false.
TDA_GW = .false.
G0W = .false.
GW0 = .false.
linGW = .false.
eta = 0d0
eta_GW = 0d0
read(1,*)
read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2, &
answer3,answer4,answer5,answer6,answer7,answer8,eta
answer3,answer4,answer5,answer6,answer7,answer8,eta_GW
if(answer1 == 'T') DIIS_GW = .true.
if(answer2 == 'T') COHSEX = .true.
if(answer3 == 'T') SOSEX = .true.
if(answer4 == 'T') BSE = .true.
if(answer5 == 'T') TDA = .true.
if(answer6 == 'T') G0W = .true.
if(answer7 == 'T') GW0 = .true.
if(answer8 == 'T') linGW = .true.
if(answer2 == 'T') COHSEX = .true.
if(answer3 == 'T') SOSEX = .true.
if(answer4 == 'T') BSE_GW = .true.
if(answer5 == 'T') TDA_GW = .true.
if(answer6 == 'T') G0W = .true.
if(answer7 == 'T') GW0 = .true.
if(answer8 == 'T') linGW = .true.
if(.not.DIIS_GW) n_diis_GW = 1
! Options for adiabatic connection