PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:42 +01:00
Code Releases Activity

create BSE standalone routine

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-01-14 21:27:34 +01:00
parent d84cff37b2
commit c2ac210e4c
12 changed files with 179 additions and 211 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
input/basis
View File

@ -1,16 +1,9 @@
1 6
S 4 1.00
234.0000000 0.0025870
35.1600000 0.0195330
7.9890000 0.0909980
2.2120000 0.2720500
1 3
S 3 1.00
38.3600000 0.0238090
5.7700000 0.1548910
1.2400000 0.4699870
S 1 1.00
0.6669000 1.0000000
S 1 1.00
0.2089000 1.0000000
0.2976000 1.0000000
P 1 1.00
3.0440000 1.0000000
P 1 1.00
0.7580000 1.0000000
D 1 1.00
1.9650000 1.0000000
1.2750000 1.0000000

4
input/methods
View File

@ -4,12 +4,12 @@
F F F
# CCD CCSD CCSD(T) ringCCD ladderCCD
F F F F F
# CIS RPA TDHF ppRPA ADC
# CIS RPA RPAx ppRPA ADC
F T T F F
# GF2 GF3
F F
# G0W0 evGW qsGW
T F F
T T T
# G0T0 evGT qsGT
F F F
# MCMP2

6
input/options
View File

@ -8,7 +8,9 @@
T T
# GF: maxSCF thresh DIIS n_diis renormalization
1 0.00001 T 10 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta
256 0.00001 T 5 F F T F F F F 0.000
# 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
# ACFDT: AC XBS
T T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

21
input/weight
View File

@ -1,16 +1,9 @@
1 6
S 4 1.00
234.0000000 0.0025870
35.1600000 0.0195330
7.9890000 0.0909980
2.2120000 0.2720500
1 3
S 3 1.00
38.3600000 0.0238090
5.7700000 0.1548910
1.2400000 0.4699870
S 1 1.00
0.6669000 1.0000000
S 1 1.00
0.2089000 1.0000000
0.2976000 1.0000000
P 1 1.00
3.0440000 1.0000000
P 1 1.00
0.7580000 1.0000000
D 1 1.00
1.9650000 1.0000000
1.2750000 1.0000000

68
src/QuAcK/Bethe_Salpeter.f90 Normal file
View File

@ -0,0 +1,68 @@
subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
! Compute the Bethe-Salpeter excitation energies
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: TDA
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eW(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision :: Omega(nS,nspin)
double precision :: XpY(nS,nS,nspin)
double precision :: XmY(nS,nS,nspin)
double precision :: rho(nBas,nBas,nS,nspin)
! Local variables
integer :: ispin
! Output variables
double precision,intent(out) :: EcRPA(nspin)
double precision,intent(out) :: EcBSE(nspin)
! Singlet manifold
if(singlet_manifold) then
ispin = 1
EcBSE(ispin) = 0d0
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,.true.,TDA,.true.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
end if
! Triplet manifold
if(triplet_manifold) then
ispin = 2
EcBSE(ispin) = 0d0
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,.true.,TDA,.true.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
end if
end subroutine Bethe_Salpeter

67
src/QuAcK/G0W0.f90
View File

@ -1,5 +1,5 @@
subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0)
subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW)
! Perform G0W0 calculation
@ -9,6 +9,8 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
! Input variables
logical,intent(in) :: doACFDT
logical,intent(in) :: doXBS
logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
logical,intent(in) :: BSE
@ -41,12 +43,9 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
logical :: adiabatic_connection
logical :: scaled_screening
! Output variables
double precision :: eG0W0(nBas)
double precision :: eGW(nBas)
! Hello world
@ -96,54 +95,23 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
! Solve the quasi-particle equation
eG0W0(:) = eHF(:) + Z(:)*SigC(:)
eGW(:) = eHF(:) + Z(:)*SigC(:)
! Dump results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eG0W0,EcRPA(ispin),EcGM)
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
! Plot stuff
call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eG0W0,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
! Perform BSE calculation
if(BSE) then
! Singlet manifold
if(singlet_manifold) then
ispin = 1
EcBSE(ispin) = 0d0
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,.true.,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
end if
! Triplet manifold
if(triplet_manifold) then
ispin = 2
EcBSE(ispin) = 0d0
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,.true.,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
end if
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -156,25 +124,22 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
! Compute the BSE correlation energy via the adiabatic connection
adiabatic_connection = .true.
scaled_screening = .true.
if(adiabatic_connection) then
if(doACFDT) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of BSE correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(scaled_screening) then
if(doXBS) then
write(*,*) '*** scaled screening version (extended BSE) ***'
write(*,*) '*** scaled screening version (XBS) ***'
write(*,*)
end if
call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eG0W0,Omega,XpY,XmY,rho)
call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho)
end if

31
src/QuAcK/QuAcK.f90
View File

@ -24,9 +24,13 @@ program QuAcK
double precision,allocatable :: ZNuc(:),rNuc(:,:)
double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:)
double precision,allocatable :: eG0W0(:)
double precision,allocatable :: eG0T0(:)
logical :: doACFDT
logical :: doXBS
integer :: nShell
integer,allocatable :: TotAngMomShell(:),KShell(:)
double precision,allocatable :: CenterShell(:,:),DShell(:,:),ExpShell(:,:)
@ -122,11 +126,13 @@ program QuAcK
! Read options for methods
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,renormalization, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, &
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,renormalization, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, &
COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, &
doACFDT,doXBS, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
! Weird stuff
@ -434,7 +440,7 @@ program QuAcK
if(doRPA) then
call cpu_time(start_RPA)
call RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
call RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_RPA)
t_RPA = end_RPA - start_RPA
@ -450,7 +456,7 @@ program QuAcK
if(doRPAx) then
call cpu_time(start_RPAx)
call RPAx(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
call RPAx(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_RPAx)
t_RPAx = end_RPAx - start_RPAx
@ -466,7 +472,7 @@ program QuAcK
if(doppRPA) then
call cpu_time(start_ppRPA)
call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF)
call ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_ppRPA)
t_ppRPA = end_ppRPA - start_ppRPA
@ -532,7 +538,7 @@ program QuAcK
if(doG0W0) then
call cpu_time(start_G0W0)
call G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
call G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,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)
@ -549,7 +555,8 @@ program QuAcK
if(doevGW) then
call cpu_time(start_evGW)
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize,eta, &
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, &
singlet_manifold,triplet_manifold,linearize,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_evGW)
@ -566,8 +573,8 @@ program QuAcK
if(doqsGW) then
call cpu_time(start_qsGW)
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW, &
COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, &
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, &
singlet_manifold,triplet_manifold,eta, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF)
call cpu_time(end_qsGW)

13
src/QuAcK/RPA.f90
View File

@ -1,4 +1,4 @@
subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e)
subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e)
! Perform a direct random phase approximation calculation
@ -8,6 +8,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
! Input variables
logical,intent(in) :: doACFDT
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBas
@ -31,8 +32,6 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
double precision :: rho
double precision :: EcRPA(nspin)
logical :: adiabatic_connection
! Hello world
write(*,*)
@ -57,7 +56,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
call linear_response(ispin,.true.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
endif
@ -69,7 +68,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
call linear_response(ispin,.true.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
endif
@ -84,9 +83,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
! Compute the correlation energy via the adiabatic connection
adiabatic_connection = .true.
if(adiabatic_connection) then
if(doACFDT) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of RPA correlation energy'

66
src/QuAcK/evGW.f90
View File

@ -1,4 +1,5 @@
subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize,eta, &
subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, &
singlet_manifold,triplet_manifold,linearize,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0)
! Perform self-consistent eigenvalue-only GW calculation
@ -13,6 +14,8 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
double precision,intent(in) :: thresh
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
logical,intent(in) :: doACFDT
logical,intent(in) :: doXBS
logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
logical,intent(in) :: BSE
@ -34,7 +37,6 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Local variables
logical :: dRPA
logical :: linear_mixing
integer :: ispin
integer :: nSCF
@ -57,9 +59,6 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
logical :: adiabatic_connection
logical :: scaled_screening
! Hello world
write(*,*)
@ -80,8 +79,6 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Switch off exchange for G0W0
dRPA = .true.
! Linear mixing
linear_mixing = .false.
@ -89,8 +86,9 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Memory allocation
allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), &
XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), &
XpY(nS,nS,nspin),XmY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization
@ -115,7 +113,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI, &
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif
@ -162,7 +160,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Print results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
! Linear mixing or DIIS extrapolation
@ -217,39 +215,8 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(BSE) then
! Singlet manifold
if(singlet_manifold) then
ispin = 1
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
! Triplet manifold
if(triplet_manifold) then
ispin = 2
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -262,24 +229,21 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Compute the BSE correlation energy via the adiabatic connection
adiabatic_connection = .true.
scaled_screening = .true.
if(adiabatic_connection) then
if(doACFDT) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of BSE correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(scaled_screening) then
if(doXBS) then
write(*,*) '*** scaled screening version (extended BSE) ***'
write(*,*) '*** scaled screening version (XBS) ***'
write(*,*)
end if
call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho)
end if

3
src/QuAcK/ppRPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
! Perform pp-RPA calculation
@ -7,6 +7,7 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
! Input variables
logical,intent(in) :: doACFDT
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBas

62
src/QuAcK/qsGW.f90
View File

@ -1,7 +1,7 @@
subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, &
subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF)
! Compute linear response
! Perform a quasiparticle self-consistent GW calculation
implicit none
include 'parameters.h'
@ -11,6 +11,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
logical,intent(in) :: doACFDT
logical,intent(in) :: doXBS
logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
logical,intent(in) :: BSE
@ -36,7 +38,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Local variables
logical :: dRPA
integer :: nSCF
integer :: nBasSq
integer :: ispin
@ -69,9 +70,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:)
logical :: adiabatic_connection
logical :: scaled_screening
! Hello world
write(*,*)
@ -99,10 +97,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(COHSEX) write(*,*) 'COHSEX approximation activated!'
write(*,*)
! Switch off exchange for qsGW
dRPA = .true.
! Memory allocation
allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
@ -144,7 +138,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif
@ -209,7 +203,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Print results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW)
! Increment
@ -251,37 +245,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(BSE) then
! Singlet manifold
if(singlet_manifold) then
ispin = 1
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
! Triplet manifold
if(triplet_manifold) then
ispin = 2
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -294,24 +259,21 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Compute the BSE correlation energy via the adiabatic connection
adiabatic_connection = .true.
scaled_screening = .true.
if(adiabatic_connection) then
if(doACFDT) then
write(*,*) '------------------------------------------------------'
write(*,*) 'Adiabatic connection version of BSE correlation energy'
write(*,*) '------------------------------------------------------'
write(*,*)
if(scaled_screening) then
if(doXBS) then
write(*,*) '*** scaled screening version (extended BSE) ***'
write(*,*) '*** scaled screening version (XBS) ***'
write(*,*)
end if
call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho)
end if

28
src/QuAcK/read_options.f90
View File

@ -1,8 +1,10 @@
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,renormalization, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, &
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
singlet_manifold,triplet_manifold, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,renormalization, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, &
COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, &
doACFDT,doXBS, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
! Read desired methods
@ -45,6 +47,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
logical,intent(out) :: linearize
double precision,intent(out) :: eta
logical,intent(out) :: doACFDT
logical,intent(out) :: doXBS
integer,intent(out) :: nMC
integer,intent(out) :: nEq
integer,intent(out) :: nWalk
@ -55,7 +60,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
! Local variables
character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7,answer8
character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7,answer8,answer9
! Open file with method specification
@ -150,6 +155,17 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
if(answer8 == 'T') linearize = .true.
if(.not.DIIS_GW) n_diis_GW = 1
! Options for adiabatic connection
doACFDT = .false.
doXBS = .false.
read(1,*)
read(1,*) answer1,answer2
if(answer1 == 'T') doACFDT = .true.
if(answer2 == 'T') doXBS = .true.
! Read options for MC-MP2: Monte Carlo steps, number of equilibration steps, number of walkers,
! Monte Carlo time step, frequency of output results, and seed for random number generator