From 87d4c8ac9aa2b75ff64338ea81d9d11c918ab3b5 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 19 Oct 2021 14:09:33 +0200 Subject: [PATCH] fixing bug in BSE@GT --- input/dft | 6 +++--- input/options | 2 +- src/MBPT/Bethe_Salpeter_Tmatrix.f90 | 4 ++-- src/MBPT/static_Tmatrix_TA.f90 | 5 +++-- src/MBPT/static_Tmatrix_TB.f90 | 4 +++- src/RPA/ACFDT_Tmatrix.f90 | 17 +++++++++++------ 6 files changed, 23 insertions(+), 15 deletions(-) diff --git a/input/dft b/input/dft index 1904956..c95d3f4 100644 --- a/input/dft +++ b/input/dft @@ -6,7 +6,7 @@ # GGA = 2: B88,G96,PBE # MGGA = 3: # Hybrid = 4: HF,B3LYP,PBE - 1 S51 + 4 HF # correlation rung: # Hartree = 0: H # LDA = 1: PW92,VWN3,VWN5,eVWN5 @@ -20,10 +20,10 @@ 1 # occupation numbers 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/input/options b/input/options index 681a821..b2f3577 100644 --- a/input/options +++ b/input/options @@ -11,7 +11,7 @@ # 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 F F F # ACFDT: AC Kx XBS - T F T + F F T # BSE: BSE dBSE dTDA evDyn T F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/src/MBPT/Bethe_Salpeter_Tmatrix.f90 b/src/MBPT/Bethe_Salpeter_Tmatrix.f90 index 6303e25..c0a79ff 100644 --- a/src/MBPT/Bethe_Salpeter_Tmatrix.f90 +++ b/src/MBPT/Bethe_Salpeter_Tmatrix.f90 @@ -117,7 +117,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Compute BSE singlet excitation energies - call linear_response_Tmatrix(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & + call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & @@ -157,7 +157,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Compute BSE triplet excitation energies - call linear_response_Tmatrix(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & + call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & diff --git a/src/MBPT/static_Tmatrix_TA.f90 b/src/MBPT/static_Tmatrix_TA.f90 index 38a3c88..a316f56 100644 --- a/src/MBPT/static_Tmatrix_TA.f90 +++ b/src/MBPT/static_Tmatrix_TA.f90 @@ -46,16 +46,17 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r do cd=1,nVV eps = Omega1(cd)**2 + eta**2 +! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/eps chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*Omega1(cd)/eps enddo do kl=1,nOO eps = Omega2(kl)**2 + eta**2 +! chi = chi + lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/eps chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/eps enddo - - TA(ia,jb) = TA(ia,jb) + lambda*ERI(i,b,j,a) - 2d0*lambda*chi + TA(ia,jb) = TA(ia,jb) + 2d0*lambda*chi enddo enddo diff --git a/src/MBPT/static_Tmatrix_TB.f90 b/src/MBPT/static_Tmatrix_TB.f90 index ffa6b0a..d7c8cc6 100644 --- a/src/MBPT/static_Tmatrix_TB.f90 +++ b/src/MBPT/static_Tmatrix_TB.f90 @@ -46,15 +46,17 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r do cd=1,nVV eps = Omega1(cd)**2 + eta**2 +! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/eps chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/eps enddo do kl=1,nOO eps = Omega2(kl)**2 + eta**2 +! chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/eps chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/eps enddo - TB(ia,jb) = TB(ia,jb) + lambda*ERI(i,j,b,a) - 2d0*lambda*chi + TB(ia,jb) = TB(ia,jb) + 2d0*lambda*chi enddo enddo diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 index 711a944..c41d3d7 100644 --- a/src/RPA/ACFDT_Tmatrix.f90 +++ b/src/RPA/ACFDT_Tmatrix.f90 @@ -110,13 +110,18 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple lambda = rAC(iAC) + ! Initialize T matrix + + TA(:,:) = 0d0 + TB(:,:) = 0d0 + if(doXBS) then isp_T = 1 iblock = 3 -! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, & -! Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T)) + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, & + Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T)) call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) @@ -126,8 +131,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple isp_T = 2 iblock = 4 -! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, & -! Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T)) + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, & + Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T)) call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) @@ -136,7 +141,7 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple end if - call linear_response_Tmatrix(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & + call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) @@ -201,7 +206,7 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple end if - call linear_response_Tmatrix(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & + call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))