LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-07 06:33:55 +01:00
Code Releases Activity

debugging BSE@GT

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-10-15 22:32:22 +02:00
parent afdf87dc58
commit ec4ff2cc9f
5 changed files with 66 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
GoHu
View File

@ -6,4 +6,7 @@ cp int/ERI.Hu.dat int/ERI.dat
cp int/Kin.Hu.dat int/Kin.dat cp int/Kin.Hu.dat int/Kin.dat
cp int/Nuc.Hu.dat int/Nuc.dat cp int/Nuc.Hu.dat int/Nuc.dat
cp int/Ov.Hu.dat int/Ov.dat cp int/Ov.Hu.dat int/Ov.dat
cp int/x.Hu.dat int/x.dat
cp int/y.Hu.dat int/y.dat
cp int/z.Hu.dat int/z.dat
./bin/QuAcK ./bin/QuAcK

4
src/LR/linear_response_Tmatrix.f90
View File

@ -13,14 +13,14 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
double precision,intent(in) :: lambda double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: A_BSE(nS,nS)
double precision,intent(in) :: B_BSE(nS,nS)
! Local variables ! Local variables
double precision :: trace_matrix double precision :: trace_matrix
double precision,allocatable :: A(:,:) double precision,allocatable :: A(:,:)
double precision,allocatable :: B(:,:) double precision,allocatable :: B(:,:)
double precision,allocatable :: A_BSE(:,:)
double precision,allocatable :: B_BSE(:,:)
double precision,allocatable :: ApB(:,:) double precision,allocatable :: ApB(:,:)
double precision,allocatable :: AmB(:,:) double precision,allocatable :: AmB(:,:)
double precision,allocatable :: AmBSq(:,:) double precision,allocatable :: AmBSq(:,:)

65
src/MBPT/Bethe_Salpeter_Tmatrix.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eT,eGT,EcBSE) subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eT,eGT,EcBSE)
! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel ! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel
@ -9,6 +9,9 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR
logical,intent(in) :: TDA_T logical,intent(in) :: TDA_T
logical,intent(in) :: TDA logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: singlet logical,intent(in) :: singlet
logical,intent(in) :: triplet logical,intent(in) :: triplet
@ -88,11 +91,11 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR
dERI = +1d0 dERI = +1d0
xERI = +0d0 xERI = +0d0
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eT(:),ERI(:,:,:,:), & call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eT,ERI, &
Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin)) Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin))
call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI, &
X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
@ -106,11 +109,11 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR
dERI = +1d0 dERI = +1d0
xERI = -1d0 xERI = -1d0
call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eT(:),ERI(:,:,:,:), & call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eT,ERI, &
Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin)) Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin))
call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI, &
X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
@ -126,12 +129,32 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR
! Compute BSE singlet excitation energies ! Compute BSE singlet excitation energies
call linear_response_Tmatrix(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & call linear_response_Tmatrix(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, &
EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin))
call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, &
OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
!-------------------------------------------------
! Compute the dynamical screening at the BSE level
!-------------------------------------------------
if(dBSE) then
! 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,dipole_int,OmRPA,rho_RPA, &
! OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
else
! call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, &
! OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
end if
end if
end if end if
!------------------- !-------------------
@ -145,12 +168,32 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR
! Compute BSE triplet excitation energies ! Compute BSE triplet excitation energies
call linear_response_Tmatrix(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & call linear_response_Tmatrix(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, &
EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin))
call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, &
OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
!-------------------------------------------------
! Compute the dynamical screening at the BSE level
!-------------------------------------------------
if(dBSE) then
! 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,dipole_int,OmRPA,rho_RPA, &
! OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
else
! call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, &
! OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin))
end if
end if
end if end if
end subroutine Bethe_Salpeter_Tmatrix end subroutine Bethe_Salpeter_Tmatrix

8
src/MBPT/G0T0.f90
View File

@ -1,4 +1,4 @@
subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0) linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0)
! Perform one-shot calculation with a T-matrix self-energy (G0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -12,7 +12,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing
logical,intent(in) :: exchange_kernel logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS logical,intent(in) :: doXBS
logical,intent(in) :: BSE logical,intent(in) :: BSE
logical,intent(in) :: TDA_W logical,intent(in) :: TDA_T
logical,intent(in) :: TDA logical,intent(in) :: TDA
logical,intent(in) :: dBSE logical,intent(in) :: dBSE
logical,intent(in) :: dTDA logical,intent(in) :: dTDA
@ -222,7 +222,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing
if(BSE) then if(BSE) then
call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, & call Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eG0T0,EcBSE) nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eG0T0,EcBSE)
if(exchange_kernel) then if(exchange_kernel) then
@ -257,7 +257,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing
end if end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta, & call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta, &
nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eG0T0,EcAC) nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eG0T0,EcAC)
if(exchange_kernel) then if(exchange_kernel) then

3
src/MBPT/static_Tmatrix_TA.f90
View File

@ -47,13 +47,16 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r
do cd=1,nVV do cd=1,nVV
eps = Omega1(cd)**2 + eta**2 eps = Omega1(cd)**2 + eta**2
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/eps chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/eps
print*,rho1(i,j,cd),rho1(a,b,cd)
enddo enddo
do kl=1,nOO do kl=1,nOO
eps = Omega2(kl)**2 + eta**2 eps = Omega2(kl)**2 + eta**2
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/eps chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/eps
print*,rho2(i,j,kl),rho2(a,b,kl)
enddo enddo
TA(ia,jb) = TA(ia,jb) + lambda*ERI(i,b,j,a) - 2d0*lambda*chi TA(ia,jb) = TA(ia,jb) + lambda*ERI(i,b,j,a) - 2d0*lambda*chi
enddo enddo