From 3ecf94bbc9dd888d39bc7acee021baa69daa8d15 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 5 Jan 2022 12:46:59 +0100 Subject: [PATCH] debugging BSE@GT --- input/methods | 8 ++++---- input/options | 2 +- src/CC/CCSD_Ec_nc.f90 | 2 +- src/GT/Bethe_Salpeter_Tmatrix.f90 | 24 ++++++++++++------------ src/GT/dynamic_Tmatrix_A.f90 | 22 +++++++++++----------- src/GT/static_Tmatrix_A.f90 | 15 ++++++++------- src/GT/static_Tmatrix_B.f90 | 11 ++++++----- src/LR/linear_response_Tmatrix.f90 | 22 ++++++++++++---------- src/RPA/ppURPA.f90 | 9 ++++----- 9 files changed, 59 insertions(+), 56 deletions(-) diff --git a/input/methods b/input/methods index 9f505d6..dcb210a 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/input/options b/input/options index d6417cb..2d0e99d 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/src/CC/CCSD_Ec_nc.f90 b/src/CC/CCSD_Ec_nc.f90 index 033ae32..d6214b0 100644 --- a/src/CC/CCSD_Ec_nc.f90 +++ b/src/CC/CCSD_Ec_nc.f90 @@ -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 diff --git a/src/GT/Bethe_Salpeter_Tmatrix.f90 b/src/GT/Bethe_Salpeter_Tmatrix.f90 index 30e76e1..05e019e 100644 --- a/src/GT/Bethe_Salpeter_Tmatrix.f90 +++ b/src/GT/Bethe_Salpeter_Tmatrix.f90 @@ -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 diff --git a/src/GT/dynamic_Tmatrix_A.f90 b/src/GT/dynamic_Tmatrix_A.f90 index 9adef18..e6770a0 100644 --- a/src/GT/dynamic_Tmatrix_A.f90 +++ b/src/GT/dynamic_Tmatrix_A.f90 @@ -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 diff --git a/src/GT/static_Tmatrix_A.f90 b/src/GT/static_Tmatrix_A.f90 index 0710ca2..9e3505e 100644 --- a/src/GT/static_Tmatrix_A.f90 +++ b/src/GT/static_Tmatrix_A.f90 @@ -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 diff --git a/src/GT/static_Tmatrix_B.f90 b/src/GT/static_Tmatrix_B.f90 index 53561be..7cf9bb0 100644 --- a/src/GT/static_Tmatrix_B.f90 +++ b/src/GT/static_Tmatrix_B.f90 @@ -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 diff --git a/src/LR/linear_response_Tmatrix.f90 b/src/LR/linear_response_Tmatrix.f90 index ec087dd..c0a4acf 100644 --- a/src/LR/linear_response_Tmatrix.f90 +++ b/src/LR/linear_response_Tmatrix.f90 @@ -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 diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90 index 395e078..73f960f 100644 --- a/src/RPA/ppURPA.f90 +++ b/src/RPA/ppURPA.f90 @@ -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)