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

fix option mismatch

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-06-17 22:10:08 +02:00
parent 6da5afae6d
commit a61fd351f9
10 changed files with 44 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/methods
View File

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

8
input/options
View File

@ -6,10 +6,10 @@
64 0.0000001 T 5
# spin: singlet triplet TDA
T T F
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.00001 T 5 T 0.0 F F T F F
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.00367493 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.00001 T 5 T 0.00367493 F F T F F
# ACFDT: AC Kx XBS
F F T
# BSE: BSE dBSE dTDA evDyn

12
src/QuAcK/BSE2_B_matrix_dynamic.f90
View File

@ -1,4 +1,4 @@
subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_dyn,ZB_dyn)
subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_dyn)
! Compute the anti-resonant part of the dynamic BSE2 matrix
@ -24,12 +24,10 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_
! Output variables
double precision,intent(out) :: B_dyn(nS,nS)
double precision,intent(out) :: ZB_dyn(nS,nS)
! Initialization
B_dyn(:,:) = 0d0
ZB_dyn(:,:) = 0d0
! Second-order correlation kernel for the block A of the singlet manifold
@ -53,14 +51,12 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_
- ERI(b,k,c,i)*ERI(a,c,j,k) + 2d0*ERI(b,k,c,i)*ERI(a,c,k,j)
B_dyn(ia,jb) = B_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
dem = + eGF(i) - eGF(c) + eGF(k) - eGF(b)
num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) &
- ERI(b,c,k,i)*ERI(a,k,j,c) + 2d0*ERI(b,c,k,i)*ERI(a,k,c,j)
B_dyn(ia,jb) = B_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
end do
end do
@ -73,7 +69,6 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_
- ERI(a,b,d,c)*ERI(c,d,i,j) + 2d0*ERI(a,b,d,c)*ERI(c,d,j,i)
B_dyn(ia,jb) = B_dyn(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2)
ZB_dyn(ia,jb) = ZB_dyn(ia,jb) - 0.5d0*num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
end do
end do
@ -86,7 +81,6 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_
- ERI(a,b,l,k)*ERI(k,l,i,j) + 2d0*ERI(a,b,l,k)*ERI(k,l,j,i)
B_dyn(ia,jb) = B_dyn(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2)
ZB_dyn(ia,jb) = ZB_dyn(ia,jb) - 0.5d0*num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
end do
end do
@ -120,13 +114,11 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_
num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) - ERI(b,k,c,i)*ERI(a,c,j,k)
B_dyn(ia,jb) = B_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
dem = + eGF(i) - eGF(c) + eGF(k) - eGF(b)
num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) - ERI(b,c,k,i)*ERI(a,k,j,c)
B_dyn(ia,jb) = B_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
end do
end do
@ -138,7 +130,6 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_
num = ERI(a,b,c,d)*ERI(c,d,j,i) + ERI(a,b,d,c)*ERI(c,d,i,j)
B_dyn(ia,jb) = B_dyn(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2)
ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + 0.5d0*num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
end do
end do
@ -150,7 +141,6 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_
num = ERI(a,b,k,l)*ERI(k,l,j,i) + ERI(a,b,l,k)*ERI(k,l,i,j)
B_dyn(ia,jb) = B_dyn(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2)
ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + 0.5d0*num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
end do
end do

10
src/QuAcK/BSE2_dynamic_perturbation.f90
View File

@ -41,13 +41,12 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,
double precision,allocatable :: ZAm_dyn(:,:)
double precision,allocatable :: B_dyn(:,:)
double precision,allocatable :: ZB_dyn(:,:)
! Memory allocation
allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS))
if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),B_dyn(nS,nS),ZB_dyn(nS,nS))
if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),B_dyn(nS,nS))
! Print main components of transition vectors
@ -82,14 +81,11 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,
! Second part of the resonant and anti-resonant part of the BSE correction (frequency independent)
call BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),Am_dyn,ZAm_dyn)
ZAm_dyn(:,:) = - ZAm_dyn(:,:)
call BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn,ZB_dyn)
call BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn)
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) &
- dot_product(Y,matmul(ZAm_dyn,Y)) &
+ dot_product(X,matmul(ZB_dyn,Y)) &
- dot_product(Y,matmul(ZB_dyn,X))
+ dot_product(Y,matmul(ZAm_dyn,Y))
OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) &
- dot_product(Y,matmul(Am_dyn,Y)) &

29
src/QuAcK/Bethe_Salpeter.f90
View File

@ -29,13 +29,11 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
! Local variables
logical :: W_BSE = .false.
integer :: ispin
integer :: isp_W
double precision,allocatable :: OmBSE(:,:)
double precision,allocatable :: XpY_BSE(:,:,:)
double precision,allocatable :: XmY_BSE(:,:,:)
double precision,allocatable :: rho_BSE(:,:,:,:)
! Output variables
@ -45,7 +43,6 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
! Memory allocation
allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin))
if(W_BSE) allocate(rho_BSE(nBas,nBas,nS,nspin))
!-------------------
! Singlet manifold
@ -77,8 +74,6 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
if(dBSE) then
if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin))
! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized)
if(evDyn) then
@ -87,14 +82,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
else
if(W_BSE) then
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin))
else
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
end if
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
end if
end if
@ -131,24 +120,16 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man
if(dBSE) then
if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin))
! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized)
if(evDyn) then
call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
else
if(W_BSE) then
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin))
else
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
end if
call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), &
XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin))
end if
end if

4
src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90
View File

@ -57,10 +57,10 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om
chi = 0d0
do kc=1,maxS
eps = OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2)
eps = OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2)
enddo

4
src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90
View File

@ -48,10 +48,10 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O
chi = 0d0
do kc=1,maxS
eps = OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
eps = OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i))
chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
enddo

36
src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90
View File

@ -84,42 +84,40 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,
! Resonant part of the BSE correction for dynamical TDA
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), &
Ap_dyn(:,:))
call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,Ap_dyn)
! Renormalization factor of the resonant parts for dynamical TDA
call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), &
ZAp_dyn(:,:))
call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,ZAp_dyn)
ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:)))
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:)))
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X))
OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X))
else
! Resonant and anti-resonant part of the BSE correction
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), &
Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:))
call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, &
Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn)
! Renormalization factor of the resonant and anti-resonant parts
call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), &
ZAp_dyn(:,:),ZAm_dyn(:,:),ZBp_dyn(:,:),ZBm_dyn(:,:))
call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, &
ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn)
ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) &
- dot_product(Y(:),matmul(ZAm_dyn(:,:),Y(:))) &
+ dot_product(X(:),matmul(ZBp_dyn(:,:),Y(:))) &
- dot_product(Y(:),matmul(ZBm_dyn(:,:),X(:)))
ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) &
- dot_product(Y,matmul(ZAm_dyn,Y)) &
+ dot_product(X,matmul(ZBp_dyn,Y)) &
- dot_product(Y,matmul(ZBm_dyn,X))
OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) &
- dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) &
+ dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) &
- dot_product(Y(:),matmul(Bm_dyn(:,:),X(:)))
OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) &
- dot_product(Y,matmul(Am_dyn,Y)) &
+ dot_product(X,matmul(Bp_dyn,Y)) &
- dot_product(Y,matmul(Bm_dyn,X))
end if
ZDyn(ia) = 1d0/(1d0 - ZDyn(ia))
ZDyn(ia) = 1d0/(1d0 - ZDyn(ia))
OmDyn(ia) = ZDyn(ia)*OmDyn(ia)
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') &

5
src/QuAcK/G0W0.f90
View File

@ -66,6 +66,10 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
write(*,*)'************************************************'
write(*,*)
! Initialization
EcRPA(:) = 0d0
! SOSEX correction
if(SOSEX) write(*,*) 'SOSEX correction activated!'
@ -148,7 +152,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)

8
src/QuAcK/read_options.f90
View File

@ -147,13 +147,13 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
thresh_GW = 1d-5
DIIS_GW = .false.
n_diis_GW = 5
linGW = .false.
eta_GW = 0d0
COHSEX = .false.
SOSEX = .false.
TDA_W = .false.
G0W = .false.
GW0 = .false.
linGW = .false.
eta_GW = 0d0
read(1,*)
read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2,eta_GW, &
@ -163,9 +163,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
if(answer2 == 'T') linGW = .true.
if(answer3 == 'T') COHSEX = .true.
if(answer4 == 'T') SOSEX = .true.
if(answer5 == 'T') TDA_W = .true.
if(answer6 == 'T') G0W = .true.
if(answer5 == 'T') G0W = .true.
if(answer7 == 'T') GW0 = .true.
if(answer7 == 'T') TDA_W = .true.
if(.not.DIIS_GW) n_diis_GW = 1
! Options for adiabatic connection