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

problem minus sign solved

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-06-01 17:14:42 +02:00
parent 503cd430ae
commit 96207d9c58
7 changed files with 113 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
input/methods
View File

@ -9,11 +9,11 @@
# CIS CID CISD # CIS CID CISD
F F F F F F
# RPA RPAx ppRPA # RPA RPAx ppRPA
F F F F T F
# G0F2 evGF2 G0F3 evGF3 # G0F2 evGF2 G0F3 evGF3
F F F F T F F F
# G0W0 evGW qsGW # G0W0 evGW qsGW
T F F F F F
# G0T0 evGT qsGT # G0T0 evGT qsGT
F F F F F F
# MCMP2 # MCMP2

2
src/QuAcK/Bethe_Salpeter.f90
View File

@ -1,5 +1,5 @@
subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE) nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE)
! Compute the Bethe-Salpeter excitation energies ! Compute the Bethe-Salpeter excitation energies

16
src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90
View File

@ -76,24 +76,24 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
do kc=1,maxS do kc=1,maxS
eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2
eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 eps_Am = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am
eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2
eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 eps_Am = (- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am
eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2
eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bp chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bp
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm
eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2
eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm
enddo enddo

20
src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90
View File

@ -59,32 +59,32 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,
do kc=1,maxS do kc=1,maxS
eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2
eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2 eps_Am = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))**2 + eta**2
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap)**2 chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Ap)**2
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am)**2 chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((- OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)))/eps_Am)**2
eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 eps_Ap = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2
eps_Am = (+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2 eps_Am = (- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))**2 + eta**2
chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap)**2 chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Ap)**2
chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((+ OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am)**2 chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*((- OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)))/eps_Am)**2
eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2
eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2 eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))**2 + eta**2
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bp)**2 chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bp)**2
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm)**2 chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((- OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)))/eps_Bm)**2
eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 eps_Bp = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2
eps_Bm = (+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2 eps_Bm = (- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))**2 + eta**2
chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp)**2 chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bp)**2
chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((+ OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm)**2 chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*((- OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)))/eps_Bm)**2
enddo enddo
ZAp(ia,jb) = ZAp(ia,jb) + 2d0*lambda*chi_Ap ZAp(ia,jb) = ZAp(ia,jb) + 2d0*lambda*chi_Ap
ZAm(ia,jb) = ZAm(ia,jb) + 2d0*lambda*chi_Am ZAm(ia,jb) = ZAm(ia,jb) - 2d0*lambda*chi_Am
ZBp(ia,jb) = ZBp(ia,jb) + 2d0*lambda*chi_Bp ZBp(ia,jb) = ZBp(ia,jb) + 2d0*lambda*chi_Bp
ZBm(ia,jb) = ZBm(ia,jb) + 2d0*lambda*chi_Bm ZBm(ia,jb) = ZBm(ia,jb) - 2d0*lambda*chi_Bm
enddo enddo
enddo enddo

92
src/QuAcK/G0F2.f90
View File

@ -1,4 +1,4 @@
subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0) subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! Perform a one-shot second-order Green function calculation ! Perform a one-shot second-order Green function calculation
@ -7,19 +7,28 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
! Input variables ! Input variables
logical,intent(in) :: BSE
logical,intent(in) :: TDA
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: linearize logical,intent(in) :: linearize
double precision,intent(in) :: eta
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nV integer,intent(in) :: nV
integer,intent(in) :: nR integer,intent(in) :: nR
double precision,intent(in) :: e0(nBas) integer,intent(in) :: nS
double precision,intent(in) :: V(nBas,nBas,nBas,nBas) double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
double precision :: eps double precision :: eps
double precision :: VV double precision :: V
double precision :: EcBSE(nspin)
double precision,allocatable :: eGF2(:) double precision,allocatable :: eGF2(:)
double precision,allocatable :: Sig(:) double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
@ -55,10 +64,10 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
do j=nC+1,nO do j=nC+1,nO
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
eps = e0(p) + e0(a) - e0(i) - e0(j) eps = eHF(p) + eHF(a) - eHF(i) - eHF(j)
VV = (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j) V = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
Sig(p) = Sig(p) + VV/eps Sig(p) = Sig(p) + V/eps
Z(p) = Z(p) + VV/eps**2 Z(p) = Z(p) + V/eps**2
end do end do
end do end do
@ -70,10 +79,10 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
do a=nO+1,nBas-nR do a=nO+1,nBas-nR
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
eps = e0(p) + e0(i) - e0(a) - e0(b) eps = eHF(p) + eHF(i) - eHF(a) - eHF(b)
VV = (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b) V = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
Sig(p) = Sig(p) + VV/eps Sig(p) = Sig(p) + V/eps
Z(p) = Z(p) + VV/eps**2 Z(p) = Z(p) + V/eps**2
end do end do
end do end do
@ -84,16 +93,69 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
if(linearize) then if(linearize) then
eGF2(:) = e0(:) + Z(:)*Sig(:) eGF2(:) = eHF(:) + Z(:)*Sig(:)
else else
eGF2(:) = e0(:) + Sig(:) eGF2(:) = eHF(:) + Sig(:)
end if end if
! Print results ! Print results
call print_G0F2(nBas,nO,e0,Sig,eGF2,Z) call print_G0F2(nBas,nO,eHF,Sig,eGF2,Z)
! Perform BSE2 calculation
if(BSE) then
call BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE2@G0F total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Compute the BSE correlation energy via the adiabatic connection
! if(doACFDT) then
! write(*,*) '------------------------------------------------------'
! write(*,*) 'Adiabatic connection version of BSE correlation energy'
! write(*,*) '------------------------------------------------------'
! write(*,*)
! if(doXBS) then
! write(*,*) '*** scaled screening version (XBS) ***'
! write(*,*)
! end if
! call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
! nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC)
! if(exchange_kernel) then
!
! EcAC(1) = 0.5d0*EcAC(1)
! EcAC(2) = 1.5d0*EcAC(1)
!
! end if
! write(*,*)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2)
! write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
! write(*,*)'-------------------------------------------------------------------------------'
! write(*,*)
! end if
end if
end subroutine G0F2 end subroutine G0F2

6
src/QuAcK/QuAcK.f90
View File

@ -601,7 +601,8 @@ program QuAcK
if(doG0F2) then if(doG0F2) then
call cpu_time(start_GF2) call cpu_time(start_GF2)
call G0F2(linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) 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 cpu_time(end_GF2) call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2 t_GF2 = end_GF2 - start_GF2
@ -617,7 +618,8 @@ program QuAcK
if(doevGF2) then if(doevGF2) then
call cpu_time(start_GF2) call cpu_time(start_GF2)
call evGF2(maxSCF_GF,thresh_GF,n_diis_GF,linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) 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 cpu_time(end_GF2) call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2 t_GF2 = end_GF2 - start_GF2

11
src/QuAcK/evGF2.f90
View File

@ -1,4 +1,5 @@
subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0) subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold,linearize, &
eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,V,e0)
! Perform eigenvalue self-consistent second-order Green function calculation ! Perform eigenvalue self-consistent second-order Green function calculation
@ -7,15 +8,23 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
! Input variables ! Input variables
logical,intent(in) :: BSE
logical,intent(in) :: TDA
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: linearize logical,intent(in) :: linearize
double precision,intent(in) :: eta
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nV integer,intent(in) :: nV
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: e0(nBas) double precision,intent(in) :: e0(nBas)
double precision,intent(in) :: V(nBas,nBas,nBas,nBas) double precision,intent(in) :: V(nBas,nBas,nBas,nBas)