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

debugging BSE@GT

Browse Source
This commit is contained in:
Pierre-Francois Loos 2022-01-05 12:46:59 +01:00
parent f6270a0ba5
commit 3ecf94bbc9
9 changed files with 59 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
input/methods
View File

@ -7,15 +7,15 @@
# drCCD rCCD crCCD lCCD
F F F F
# CIS* CIS(D) CID CISD FCI
F F T T F
F F F F F
# RPA* RPAx* crRPA ppRPA
F F F F
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
T F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW
F F F F F
T F F F F
# G0T0 evGT qsGT
F F F
T F F
# MCMP2
F
# * unrestricted version available

2
input/options
View File

@ -5,7 +5,7 @@
# CC: maxSCF thresh DIIS n_diis
64 0.00001 T 5
# spin: TDA singlet triplet spin_conserved spin_flip
F T T T T
F T F T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm reg
256 0.00001 T 5 T 0.0 3 F
# GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg

2
src/CC/CCSD_Ec_nc.f90
View File

@ -1,6 +1,6 @@
subroutine CCSD_Ec_nc(nO,nV,t1,t2,Fov,OOVV,EcCCSD)
! Compute the CCSD correlatio energy in non-conanical form
! Compute the CCSD correlatio energy in non-canonical form
implicit none

24
src/GT/Bethe_Salpeter_Tmatrix.f90
View File

@ -88,13 +88,13 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
print*,'aa block of TA'
call matout(nS,nS,TA)
print*,'aa block of TB'
call matout(nS,nS,TB)
! print*,'aa block of TA'
! call matout(nS,nS,TA)
! print*,'aa block of TB'
! call matout(nS,nS,TB)
!----------------------------------------------
! Compute T-matrix for alpha-alpha block
@ -108,13 +108,13 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,
! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
print*,'aa+ab block of TA'
call matout(nS,nS,TA)
print*,'aa+ab block of TB'
call matout(nS,nS,TB)
! print*,'aa+ab block of TA'
! call matout(nS,nS,TA)
! print*,'aa+ab block of TB'
! call matout(nS,nS,TB)
!-------------------
! Singlet manifold

22
src/GT/dynamic_Tmatrix_A.f90
View File

@ -58,44 +58,44 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
chi = 0d0
do cd=1,nVV
eps = + Omega1(cd)
chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*eps/(eps**2 + eta**2)
eps = - Omega1(cd)
chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*eps/(eps**2 + eta**2)
end do
do kl=1,nOO
eps = - Omega2(kl)
chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*eps/(eps**2 + eta**2)
eps = + Omega2(kl)
chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*eps/(eps**2 + eta**2)
end do
A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi
A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi
chi = 0d0
do cd=1,nVV
eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*eps/(eps**2 + eta**2)
chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*eps/(eps**2 + eta**2)
end do
do kl=1,nOO
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*eps/(eps**2 + eta**2)
chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*eps/(eps**2 + eta**2)
end do
A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi
A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi
chi = 0d0
do cd=1,nVV
eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
do kl=1,nOO
eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + 1d0*lambda*chi
ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - 1d0*lambda*chi
end do
end do

15
src/GT/static_Tmatrix_A.f90
View File

@ -1,4 +1,4 @@
subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TA)
subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TA)
! Compute the OOVV block of the static T-matrix for the resonant block
@ -7,6 +7,8 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh
! Input variables
integer,intent(in) :: ispin
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -15,7 +17,6 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh
integer,intent(in) :: nS
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega1(nVV)
@ -45,18 +46,18 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh
chi = 0d0
do cd=1,nVV
eps = + Omega1(cd)
eps = - Omega1(cd)
! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*eps/(eps**2 + eta**2)
chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2)
enddo
do kl=1,nOO
eps = - Omega2(kl)
eps = + Omega2(kl)
! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*eps/(eps**2 + eta**2)
chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2)
enddo
TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi
TA(ia,jb) = TA(ia,jb) + lambda*chi
enddo
enddo

11
src/GT/static_Tmatrix_B.f90
View File

@ -1,4 +1,4 @@
subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TB)
subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TB)
! Compute the OVVO block of the static T-matrix for the coupling block
@ -7,6 +7,8 @@ subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh
! Input variables
integer,intent(in) :: ispin
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -15,7 +17,6 @@ subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh
integer,intent(in) :: nS
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega1(nVV)
@ -45,18 +46,18 @@ subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh
chi = 0d0
do cd=1,nVV
eps = + Omega1(cd)
eps = - Omega1(cd)
! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
enddo
do kl=1,nOO
eps = - Omega2(kl)
eps = + Omega2(kl)
! chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
enddo
TB(ia,jb) = TB(ia,jb) + 1d0*lambda*chi
TB(ia,jb) = TB(ia,jb) + lambda*chi
enddo
enddo

22
src/LR/linear_response_Tmatrix.f90
View File

@ -42,12 +42,13 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
A(:,:) = A(:,:) + A_BSE(:,:)
if(ispin == 1) A(:,:) = A(:,:) + A_BSE(:,:)
if(ispin == 2) A(:,:) = A(:,:) - A_BSE(:,:)
print*,'A'
call matout(nS,nS,A)
print*,'TA'
call matout(nS,nS,A_BSE)
! print*,'A'
! call matout(nS,nS,A)
! print*,'TA'
! call matout(nS,nS,A_BSE)
! Tamm-Dancoff approximation
@ -63,12 +64,13 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
B(:,:) = B(:,:) + B_BSE(:,:)
if(ispin == 1) B(:,:) = B(:,:) + B_BSE(:,:)
if(ispin == 2) B(:,:) = B(:,:) - B_BSE(:,:)
print*,'B'
call matout(nS,nS,B)
print*,'TB'
call matout(nS,nS,B_BSE)
! print*,'B'
! call matout(nS,nS,B)
! print*,'TB'
! call matout(nS,nS,B_BSE)
! Build A + B and A - B matrices

9
src/RPA/ppURPA.f90
View File

@ -67,12 +67,11 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH
! Memory allocation
allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), &
Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc))
allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), &
Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc))
call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sc, &
nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc,Omega2sc,X2sc,Y2sc,&
Ec_ppURPA(ispin))
call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sc,nHaa,nHab,nHbb,nH_sc,1d0, &
e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc,Omega2sc,X2sc,Y2sc,Ec_ppURPA(ispin))
call print_excitation('pp-RPA (N+2)',5,nP_sc,Omega1sc)
call print_excitation('pp-RPA (N-2)',5,nH_sc,Omega2sc)