LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-26 15:12:17 +02:00
Code Releases Activity

BSE Ec

Browse Source
This commit is contained in:
Pierre-Francois Loos 2019-05-04 17:13:50 +02:00
parent 634e0ecf09
commit f345b8e2be
15 changed files with 263 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

170
input/basis
View File

@ -1,4 +1,172 @@
1 16 1 100
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00 S 1 1.00
1.0000000 1.0000000 1.0000000 1.0000000
S 1 1.00 S 1 1.00

4
input/methods
View File

@ -1,13 +1,13 @@
# RHF UHF MOM # RHF UHF MOM
T F F T F F
# MP2 MP3 MP2-F12 # MP2 MP3 MP2-F12
T F F F F F
# CCD CCSD CCSD(T) # CCD CCSD CCSD(T)
F F F F F F
# CIS TDHF ADC # CIS TDHF ADC
F F F F F F
# GF2 GF3 # GF2 GF3
T T F F
# G0W0 evGW qsGW # G0W0 evGW qsGW
T F F T F F
# MCMP2 # MCMP2

2
input/molecule
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd # nAt nEla nElb nCore nRyd
1 4 4 0 0 1 9 9 0 0
# Znuc x y z # Znuc x y z
X 0.0 0.0 0.0 X 0.0 0.0 0.0

4
input/options
View File

@ -5,10 +5,10 @@
# CC: maxSCF thresh DIIS n_diis # CC: maxSCF thresh DIIS n_diis
64 0.00001 F 1 64 0.00001 F 1
# CIS/TDHF: singlet triplet # CIS/TDHF: singlet triplet
T F T T
# GF: maxSCF thresh DIIS n_diis renormalization # GF: maxSCF thresh DIIS n_diis renormalization
64 0.00001 T 5 3 64 0.00001 T 5 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize
64 0.00001 T 5 F F T F F F F 64 0.00001 T 3 F F T F F F F
# 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

4
run_sph.sh
View File

@ -1,8 +1,8 @@
#! /bin/bash #! /bin/bash
Lmin=0 Lmin=0
Lmax=2 Lmax=0
Mmax=3 Mmax=9
rs=$1 rs=$1
if [ $# != 1 ] if [ $# != 1 ]

55
src/QuAcK/ADC2.f90
View File

@ -12,7 +12,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
integer,intent(in) :: nBas,nC,nO,nV,nR integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
@ -57,8 +58,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
! Memory allocation ! Memory allocation
allocate(db_ERI(nBas,nBas,nBas,nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eOld(nADC), & allocate(db_ERI(nBas,nBas,nBas,nBas),error_diis(nADC,max_diis),e_diis(nADC,max_diis),eOld(nADC), &
B_ADC(nADC,nADC),X_ADC(nADC,nADC),e_ADC(nADC),G_ADC(nADC,nADC),SigInf(nADC,nADC)) B_ADC(nADC,nADC),X_ADC(nADC,nADC),e_ADC(nADC),G_ADC(nBas,nBas),SigInf(nBas,nBas))
! Create double-bar MO integrals ! Create double-bar MO integrals
@ -99,11 +100,11 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
! !
B_ADC(:,:) = 0d0 B_ADC(:,:) = 0d0
B_ADC(nC+1:nV,nC+1:nV) = SigInf(nC+1:nV,nC+1:nV) B_ADC(nC+1:nBas-nR,nC+1:nBas-nR) = SigInf(nC+1:nBas-nR,nC+1:nBas-nR)
jADC = 0 jADC = 0
do p=nC+1,nV do p=nC+1,nBas-nR
jADC = jADC + 1 jADC = jADC + 1
B_ADC(jADC,jADC) = e(p) B_ADC(jADC,jADC) = e(p)
@ -114,7 +115,7 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
! U matrices -- Eq. (40a) -- ! U matrices -- Eq. (40a) --
! !
do p=nC+1,nV do p=nC+1,nBas-nR
iADC = p - nC iADC = p - nC
jADC = nH + nP jADC = nH + nP
@ -122,8 +123,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
! U^I: 2p-1h -- Eqs. (40a) and (41a) -- ! U^I: 2p-1h -- Eqs. (40a) and (41a) --
do i=nC+1,nO do i=nC+1,nO
do a=nO+1,nV-nR do a=nO+1,nBas-nR
do b=a,nV-nR do b=a,nBas-nR
jADC = jADC + 1 jADC = jADC + 1
B_ADC(iADC,jADC) = db_ERI(p,i,a,b) B_ADC(iADC,jADC) = db_ERI(p,i,a,b)
@ -136,7 +137,7 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
do i=nC+1,nO do i=nC+1,nO
do j=i,nO do j=i,nO
do a=nO+1,nV-nR do a=nO+1,nBas-nR
jADC = jADC + 1 jADC = jADC + 1
B_ADC(iADC,jADC) = db_ERI(p,a,i,j) B_ADC(iADC,jADC) = db_ERI(p,a,i,j)
@ -156,8 +157,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
jADC = nH + nP jADC = nH + nP
do i=nC+1,nO do i=nC+1,nO
do a=nO+1,nV-nR do a=nO+1,nBas-nR
do b=a,nV-nR do b=a,nBas-nR
jADC = jADC + 1 jADC = jADC + 1
B_ADC(jADC,jADC) = e(a) + e(b) - e(i) B_ADC(jADC,jADC) = e(a) + e(b) - e(i)
@ -170,7 +171,7 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
do i=nC+1,nO do i=nC+1,nO
do j=i,nO do j=i,nO
do a=nO+1,nV do a=nO+1,nBas-nR
jADC = jADC + 1 jADC = jADC + 1
B_ADC(jADC,jADC) = e(i) + e(j) - e(a) B_ADC(jADC,jADC) = e(i) + e(j) - e(a)
@ -188,15 +189,15 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
iADC = nH + nP iADC = nH + nP
do i=nC+1,nO do i=nC+1,nO
do a=nO+1,nV-nR do a=nO+1,nBas-nR
do b=a,nV-nR do b=a,nBas-nR
iADC = iADC + 1 iADC = iADC + 1
jADC = nH + nP jADC = nH + nP
do j=nC+1,nO do j=nC+1,nO
do c=nO+1,nV do c=nO+1,nBas-nR
do d=c,nV-nR do d=c,nBas-nR
jADC = jADC + 1 jADC = jADC + 1
B_ADC(iADC,jADC) = B_ADC(iADC,jADC) & B_ADC(iADC,jADC) = B_ADC(iADC,jADC) &
@ -221,14 +222,14 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
do i=nC+1,nO do i=nC+1,nO
do j=i,nO do j=i,nO
do a=nO+1,nV-nR do a=nO+1,nBas-nR
iADC = iADC + 1 iADC = iADC + 1
jADC = nH + nP + nH*nPP jADC = nH + nP + nH*nPP
do k=nC+1,nO do k=nC+1,nO
do l=k,nO do l=k,nO
do b=nO+1,nV-nR do b=nO+1,nBas-nR
jADC = jADC + 1 jADC = jADC + 1
B_ADC(iADC,jADC) = B_ADC(iADC,jADC) & B_ADC(iADC,jADC) = B_ADC(iADC,jADC) &
@ -270,7 +271,7 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
do iADC=1,nADC do iADC=1,nADC
if(NORM2(X_ADC(1:nH+nP,iADC)) > 0.1d0 ) & ! if(NORM2(X_ADC(1:nH+nP,iADC)) > 0.1d0 ) &
write(*,'(2(2X,F12.6))') e_ADC(iADC)*HaToeV,NORM2(X_ADC(1:nH+nP,iADC)) write(*,'(2(2X,F12.6))') e_ADC(iADC)*HaToeV,NORM2(X_ADC(1:nH+nP,iADC))
enddo enddo
@ -290,8 +291,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
SigInf(:,:) = 0d0 SigInf(:,:) = 0d0
do i=nC+1,nO do i=nC+1,nO
do p=nC+1,nV-nR do p=nC+1,nBas-nR
do q=nC+1,nV-nR do q=nC+1,nBas-nR
SigInf(p,q) = SigInf(p,q) - db_ERI(p,i,q,i) SigInf(p,q) = SigInf(p,q) - db_ERI(p,i,q,i)
@ -307,8 +308,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
if(e_ADC(iADC) > 0d0 ) cycle if(e_ADC(iADC) > 0d0 ) cycle
do p=nC+1,nV-nR do p=nC+1,nBas-nR
do q=nC+1,nV-nR do q=nC+1,nBas-nR
G_ADC(p,q) = G_ADC(p,q) + X_ADC(p,iADC)*X_ADC(q,iADC) G_ADC(p,q) = G_ADC(p,q) + X_ADC(p,iADC)*X_ADC(q,iADC)
@ -319,10 +320,10 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI)
! Compute static self-energy for next iteration -- Eq. (25) -- ! Compute static self-energy for next iteration -- Eq. (25) --
do p=nC+1,nV-nR do r=nC+1,nBas-nR
do q=nC+1,nV-nR do s=nC+1,nBas-nR
do r=nC+1,nV-nR do p=nC+1,nBas-nR
do s=nC+1,nV-nR do q=nC+1,nBas-nR
SigInf(p,q) = SigInf(p,q) + db_ERI(p,r,q,s)*G_ADC(r,s) SigInf(p,q) = SigInf(p,q) + db_ERI(p,r,q,s)*G_ADC(r,s)

4
src/QuAcK/CCSD.f90
View File

@ -228,7 +228,7 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
write(*,*)' CCSD energy ' write(*,*)' CCSD energy '
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,'(1X,A20,1X,F15.10)')' E(CCSD) = ',ECCSD write(*,'(1X,A20,1X,F15.10)')' E(CCSD) = ',ECCSD
write(*,'(1X,A20,1X,F10.6)') ' Ec(CCSD) = ',EcCCSD write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD) = ',EcCCSD
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,*) write(*,*)
@ -259,7 +259,7 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
write(*,*)' CCSD(T) energy ' write(*,*)' CCSD(T) energy '
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,'(1X,A20,1X,F15.10)')' E(CCSD(T)) = ',ECCSD + EcCCT write(*,'(1X,A20,1X,F15.10)')' E(CCSD(T)) = ',ECCSD + EcCCT
write(*,'(1X,A20,1X,F10.6)') ' Ec(CCSD(T)) = ',EcCCSD + EcCCT write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD(T)) = ',EcCCSD + EcCCT
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,*) write(*,*)

28
src/QuAcK/G0W0.f90
View File

@ -26,8 +26,9 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
! Local variables ! Local variables
logical :: dRPA logical :: dRPA
integer :: ispin integer :: ispin,jspin
double precision :: EcRPA double precision :: EcRPA(nspin)
double precision :: EcBSE(nspin)
double precision :: EcGM double precision :: EcGM
double precision,allocatable :: H(:,:) double precision,allocatable :: H(:,:)
double precision,allocatable :: SigC(:) double precision,allocatable :: SigC(:)
@ -74,7 +75,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
! Compute linear response ! Compute linear response
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, & call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
@ -104,7 +105,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
! Dump results ! Dump results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eG0W0,EcRPA,EcGM) call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eG0W0,EcRPA(ispin),EcGM)
! Plot stuff ! Plot stuff
@ -119,8 +120,10 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
if(singlet_manifold) then if(singlet_manifold) then
ispin = 1 ispin = 1
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif endif
@ -130,16 +133,27 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
if(triplet_manifold) then if(triplet_manifold) then
ispin = 2 ispin = 2
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, & call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,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,eG0W0,ERI_MO_basis, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, &
rho(:,:,:,1),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,1),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif endif
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
endif endif

26
src/QuAcK/evGW.f90
View File

@ -40,7 +40,8 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
integer :: n_diis integer :: n_diis
double precision :: rcond double precision :: rcond
double precision :: Conv double precision :: Conv
double precision :: EcRPA double precision :: EcRPA(nspin)
double precision :: EcBSE(nspin)
double precision :: EcGM double precision :: EcGM
double precision :: lambda double precision :: lambda
double precision,allocatable :: error_diis(:,:) double precision,allocatable :: error_diis(:,:)
@ -110,7 +111,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(.not. GW0 .or. nSCF == 0) then if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, & call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
endif endif
@ -155,7 +156,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Print results ! 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) call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
! Linear mixing or DIIS extrapolation ! Linear mixing or DIIS extrapolation
@ -213,8 +214,10 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(singlet_manifold) then if(singlet_manifold) then
ispin = 1 ispin = 1
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif endif
@ -223,16 +226,27 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(triplet_manifold) then if(triplet_manifold) then
ispin = 2 ispin = 2
EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, & call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,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,eGW,ERI_MO_basis, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
endif endif
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
endif endif
end subroutine evGW end subroutine evGW

2
src/QuAcK/linear_response_B_matrix.f90
View File

@ -32,7 +32,7 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B_lr)
delta_dRPA = 0d0 delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0 if(dRPA) delta_dRPA = 1d0
! Build A matrix ! Build B matrix
ia = 0 ia = 0
do i=nC+1,nO do i=nC+1,nO

5
src/QuAcK/print_G0W0.f90
View File

@ -6,7 +6,10 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
include 'parameters.h' include 'parameters.h'
integer,intent(in) :: nBas,nO integer,intent(in) :: nBas,nO
double precision,intent(in) :: ENuc,EHF,EcRPA,EcGM double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: EcRPA
double precision,intent(in) :: EcGM
double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
integer :: x,HOMO,LUMO integer :: x,HOMO,LUMO

2
src/QuAcK/print_RHF.f90
View File

@ -47,7 +47,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF)
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') 'MO coefficients' write(*,'(A32)') 'MO coefficients'
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,cHF) ! call matout(nBas,nBas,cHF)
write(*,*) write(*,*)
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') 'MO energies' write(*,'(A32)') 'MO energies'

5
src/QuAcK/print_evGW.f90
View File

@ -6,7 +6,10 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
include 'parameters.h' include 'parameters.h'
integer,intent(in) :: nBas,nO,nSCF integer,intent(in) :: nBas,nO,nSCF
double precision,intent(in) :: ENuc,EHF,EcRPA,EcGM double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: EcRPA
double precision,intent(in) :: EcGM
double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
integer :: x,HOMO,LUMO integer :: x,HOMO,LUMO

5
src/QuAcK/self_energy_correlation.f90
View File

@ -27,7 +27,8 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,
! Infinitesimal ! Infinitesimal
eta = 0.001d0 eta = 0.0d0
! eta = 0.001d0
! COHSEX static approximation ! COHSEX static approximation
@ -68,7 +69,7 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,
EcGM = 0d0 EcGM = 0d0
do i=nC+1,nO do i=nC+1,nO
EcGM = EcGM + SigC(i,i) EcGM = EcGM + 0.5d0*SigC(i,i)
enddo enddo
else else

6
src/QuAcK/self_energy_correlation_diag.f90
View File

@ -68,7 +68,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega
EcGM = 0d0 EcGM = 0d0
do i=nC+1,nO do i=nC+1,nO
EcGM = EcGM + SigC(i) EcGM = EcGM + 0.5d0*SigC(i)
enddo enddo
else else
@ -117,7 +117,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
jb = jb + 1 jb = jb + 1
eps = e(a) - e(i) + Omega(jb) eps = e(a) - e(i) + Omega(jb)
EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)/eps EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)/eps
enddo enddo
enddo enddo
enddo enddo
@ -164,7 +164,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
jb = jb + 1 jb = jb + 1
eps = e(a) - e(i) + Omega(jb) eps = e(a) - e(i) + Omega(jb)
EcGM = EcGM + 2d0*rho(a,i,jb)*rhox(a,i,jb)/eps EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)/eps
enddo enddo
enddo enddo
enddo enddo