From c379b80eacbbddffd79807d03e969f395930fb4c Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 18 Nov 2021 13:34:08 +0100 Subject: [PATCH] ACFDT for ppRPA --- input/methods | 2 +- input/options | 2 +- src/AOtoMO/AOtoMO_integral_transform.f90 | 2 +- src/AOtoMO/AOtoMO_transform.f90 | 9 ++ src/LR/linear_response_B_pp.f90 | 4 +- src/LR/linear_response_C_pp_od.f90 | 91 ++++++++++++++ src/LR/linear_response_D_pp_od.f90 | 91 ++++++++++++++ src/LR/linear_response_pp.f90 | 2 +- src/MBPT/G0T0.f90 | 2 +- src/MBPT/dynamic_Tmatrix_A.f90 | 4 +- src/MBPT/evGT.f90 | 2 +- src/MBPT/static_Tmatrix_TA.f90 | 2 +- src/MBPT/static_Tmatrix_TB.f90 | 4 +- src/QuAcK/QuAcK.f90 | 4 +- src/RPA/ACFDT_Tmatrix.f90 | 76 ++++++------ src/RPA/ACFDT_pp.f90 | 145 +++++++++++++++++++++++ src/RPA/ACFDT_pp_correlation_energy.f90 | 51 ++++++++ src/RPA/ppRPA.f90 | 14 +-- 18 files changed, 446 insertions(+), 61 deletions(-) create mode 100644 src/LR/linear_response_C_pp_od.f90 create mode 100644 src/LR/linear_response_D_pp_od.f90 create mode 100644 src/RPA/ACFDT_pp.f90 create mode 100644 src/RPA/ACFDT_pp_correlation_energy.f90 diff --git a/input/methods b/input/methods index fa793c5..1e67d49 100644 --- a/input/methods +++ b/input/methods @@ -9,7 +9,7 @@ # CIS* CIS(D) CID CISD FCI F F F F F # RPA* RPAx* crRPA ppRPA - F T T F + F F F T # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW diff --git a/input/options b/input/options index 8dad38d..f54eee2 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.00367493 F F F F F # ACFDT: AC Kx XBS - T T T + T T F # BSE: BSE dBSE dTDA evDyn T F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/src/AOtoMO/AOtoMO_integral_transform.f90 b/src/AOtoMO/AOtoMO_integral_transform.f90 index dabeee8..fbb0d4b 100644 --- a/src/AOtoMO/AOtoMO_integral_transform.f90 +++ b/src/AOtoMO/AOtoMO_integral_transform.f90 @@ -79,7 +79,7 @@ subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI do nu=1,nBas ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket1)*scr(i,nu,k,l) enddo -! print*,i,k,j,l,ERI_MO_basis(i,j,k,l) +! write(11,'(I5,I5,I5,I5,F16.10)') i,j,k,l,ERI_MO_basis(i,j,k,l) enddo enddo enddo diff --git a/src/AOtoMO/AOtoMO_transform.f90 b/src/AOtoMO/AOtoMO_transform.f90 index 7919084..a4cd5a3 100644 --- a/src/AOtoMO/AOtoMO_transform.f90 +++ b/src/AOtoMO/AOtoMO_transform.f90 @@ -9,10 +9,19 @@ subroutine AOtoMO_transform(nBas,c,A) integer,intent(in) :: nBas double precision,intent(in) :: c(nBas,nBas) + integer :: i,j + ! Output variables double precision,intent(inout):: A(nBas,nBas) A = matmul(transpose(c),matmul(A,c)) +! do j=1,nBas +! do i=1,nBas +! write(10,'(I5,I5,F16.10)') i,j,A(i,j) +! enddo +! enddo + + end subroutine AOtoMO_transform diff --git a/src/LR/linear_response_B_pp.f90 b/src/LR/linear_response_B_pp.f90 index 3760373..46c151e 100644 --- a/src/LR/linear_response_B_pp.f90 +++ b/src/LR/linear_response_B_pp.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B_pp) +subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp) ! Compute the B matrix of the pp channel @@ -10,7 +10,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B_pp integer,intent(in) :: ispin integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV double precision,intent(in) :: lambda - double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables diff --git a/src/LR/linear_response_C_pp_od.f90 b/src/LR/linear_response_C_pp_od.f90 new file mode 100644 index 0000000..82b22eb --- /dev/null +++ b/src/LR/linear_response_C_pp_od.f90 @@ -0,0 +1,91 @@ +subroutine linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,C_pp) + +! Compute the C matrix of the pp channel (without diagonal term) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision,external :: Kronecker_delta + + integer :: a,b,c,d,ab,cd + +! Output variables + + double precision,intent(out) :: C_pp(nVV,nVV) + +! Build C matrix for the singlet manifold + + if(ispin == 1) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a,nBas-nR + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=c,nBas-nR + cd = cd + 1 + + C_pp(ab,cd) = lambda*(ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + + end do + end do + end do + end do + + end if + +! Build C matrix for the triplet manifold, or alpha-alpha block, or in the spin-orbital basis + + if(ispin == 2 .or. ispin == 4) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + + C_pp(ab,cd) = lambda*(ERI(a,b,c,d) - ERI(a,b,d,c)) + + end do + end do + end do + end do + + end if + +! Build the alpha-beta block of the C matrix + + if(ispin == 3) then + + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + cd = cd + 1 + + C_pp(ab,cd) = lambda*ERI(a,b,c,d) + + end do + end do + end do + end do + + end if + +end subroutine linear_response_C_pp_od diff --git a/src/LR/linear_response_D_pp_od.f90 b/src/LR/linear_response_D_pp_od.f90 new file mode 100644 index 0000000..d9c902d --- /dev/null +++ b/src/LR/linear_response_D_pp_od.f90 @@ -0,0 +1,91 @@ +subroutine linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,D_pp) + +! Compute the D matrix of the pp channel (without the diagonal term) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision,external :: Kronecker_delta + + integer :: i,j,k,l,ij,kl + +! Output variables + + double precision,intent(out) :: D_pp(nOO,nOO) + +! Build the D matrix for the singlet manifold + + if(ispin == 1) then + + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + + D_pp(ij,kl) = lambda*(ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + + end do + end do + end do + end do + + end if + +! Build the D matrix for the triplet manifold, the alpha-alpha block, or in the spin-orbital basis + + if(ispin == 2 .or. ispin == 4) then + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + + D_pp(ij,kl) = lambda*(ERI(i,j,k,l) - ERI(i,j,l,k)) + + end do + end do + end do + end do + + end if + +! Build the alpha-beta block of the D matrix + + if(ispin == 3) then + + ij = 0 + do i=nC+1,nO + do j=nC+1,nO + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=nC+1,nO + kl = kl + 1 + + D_pp(ij,kl) = lambda*ERI(i,j,k,l) + + end do + end do + end do + end do + + end if + +end subroutine linear_response_D_pp_od diff --git a/src/LR/linear_response_pp.f90 b/src/LR/linear_response_pp.f90 index 2d9873b..a1decea 100644 --- a/src/LR/linear_response_pp.f90 +++ b/src/LR/linear_response_pp.f90 @@ -75,7 +75,7 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om else - call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B) + call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B) ! Diagonal blocks diff --git a/src/MBPT/G0T0.f90 b/src/MBPT/G0T0.f90 index 52a080d..a838c40 100644 --- a/src/MBPT/G0T0.f90 +++ b/src/MBPT/G0T0.f90 @@ -245,7 +245,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing if(exchange_kernel) then EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(2) end if diff --git a/src/MBPT/dynamic_Tmatrix_A.f90 b/src/MBPT/dynamic_Tmatrix_A.f90 index 31e75e2..55d4c44 100644 --- a/src/MBPT/dynamic_Tmatrix_A.f90 +++ b/src/MBPT/dynamic_Tmatrix_A.f90 @@ -76,7 +76,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O do kl=1,nOO eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) - chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) end do A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi @@ -89,7 +89,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O end do do kl=1,nOO - eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + eps = + OmBSE - Omega2(kl) - (eGT(a) + eGT(b)) chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do diff --git a/src/MBPT/evGT.f90 b/src/MBPT/evGT.f90 index 54405e6..12bf2a5 100644 --- a/src/MBPT/evGT.f90 +++ b/src/MBPT/evGT.f90 @@ -287,7 +287,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & if(exchange_kernel) then EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(2) end if diff --git a/src/MBPT/static_Tmatrix_TA.f90 b/src/MBPT/static_Tmatrix_TA.f90 index ac4ed6e..c3f2583 100644 --- a/src/MBPT/static_Tmatrix_TA.f90 +++ b/src/MBPT/static_Tmatrix_TA.f90 @@ -50,7 +50,7 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r do kl=1,nOO ! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) - chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) + chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) enddo TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi diff --git a/src/MBPT/static_Tmatrix_TB.f90 b/src/MBPT/static_Tmatrix_TB.f90 index 3613980..d502c89 100644 --- a/src/MBPT/static_Tmatrix_TB.f90 +++ b/src/MBPT/static_Tmatrix_TB.f90 @@ -49,8 +49,8 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r enddo do kl=1,nOO -! chi = chi - lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 - chi = chi - rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 +! chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 + chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 enddo TB(ia,jb) = TB(ia,jb) + 1d0*lambda*chi diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 7af842b..edaa4ab 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -442,7 +442,7 @@ program QuAcK ket1 = 1 ket2 = 1 call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO) - +! call AOtoMO_transform(nBas,cHF,T+V) end if end if @@ -824,7 +824,7 @@ program QuAcK if(doppRPA) then call cpu_time(start_RPA) - call ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) + call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_RPA) t_RPA = end_RPA - start_RPA diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 index 256abb9..3245e34 100644 --- a/src/RPA/ACFDT_Tmatrix.f90 +++ b/src/RPA/ACFDT_Tmatrix.f90 @@ -64,8 +64,11 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple ! Useful quantities - nOOs = nO*nO - nVVs = nV*nV +! nOOs = nO*nO +! nVVs = nV*nV + + nOOs = nO*(nO+1)/2 + nVVs = nV*(nV+1)/2 nOOt = nO*(nO-1)/2 nVVt = nV*(nV-1)/2 @@ -118,31 +121,31 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple TA(:,:) = 0d0 TB(:,:) = 0d0 - if(doXBS) then +! if(doXBS) then - isp_T = 1 - iblock = 3 +! 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) +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) - if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) +! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) +! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) - isp_T = 2 - iblock = 4 +! 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) +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) - if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) +! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) +! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) - end if +! end if 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)) @@ -183,31 +186,36 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple lambda = rAC(iAC) - if(doXBS) then + ! Initialize T matrix - isp_T = 1 - iblock = 3 + TA(:,:) = 0d0 + TB(:,:) = 0d0 - 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)) +! if(doXBS) then - call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) +! isp_T = 1 +! iblock = 3 - call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) - if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) +! 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)) - isp_T = 2 - iblock = 4 +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - 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 static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) +! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) - call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) +! isp_T = 2 +! iblock = 4 - call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) - if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) +! 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)) - end if +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + +! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) +! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) + +! end if 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)) diff --git a/src/RPA/ACFDT_pp.f90 b/src/RPA/ACFDT_pp.f90 new file mode 100644 index 0000000..540c77c --- /dev/null +++ b/src/RPA/ACFDT_pp.f90 @@ -0,0 +1,145 @@ +subroutine ACFDT_pp(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) + +! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for pp sector + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: singlet + logical,intent(in) :: triplet + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: ispin + integer :: iAC + double precision :: lambda + double precision,allocatable :: Ec(:,:) + + integer :: nOOs,nOOt + integer :: nVVs,nVVt + + double precision,allocatable :: Omega1s(:),Omega1t(:) + double precision,allocatable :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:) + double precision,allocatable :: Omega2s(:),Omega2t(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) + double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) + +! Output variables + + double precision,intent(out) :: EcAC(nspin) + +! Useful quantities + + nOOs = nO*(nO+1)/2 + nVVs = nV*(nV+1)/2 + + nOOt = nO*(nO-1)/2 + nVVt = nV*(nV-1)/2 + +! Memory allocation + + allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & + Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) + allocate(Ec(nAC,nspin)) + +! Antisymmetrized kernel version + + EcAC(:) = 0d0 + Ec(:,:) = 0d0 + +! Singlet manifold + + if(singlet) then + + ispin = 1 + + write(*,*) '--------------' + write(*,*) 'Singlet states' + write(*,*) '--------------' + write(*,*) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' + write(*,*) '-----------------------------------------------------------------------------------' + + do iAC=1,nAC + + lambda = rAC(iAC) + + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,e,ERI,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcAC(ispin)) + + call ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOOs,nVVs,X1s,Y1s,X2s,Y2s,Ec(iAC,ispin)) + + write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) + + end do + + EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) + + end if + +! Triplet manifold + + if(triplet) then + + ispin = 2 + + write(*,*) '--------------' + write(*,*) 'Triplet states' + write(*,*) '--------------' + write(*,*) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' + write(*,*) '-----------------------------------------------------------------------------------' + + do iAC=1,nAC + + lambda = rAC(iAC) + + ! Initialize T matrix + + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,e,ERI,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcAC(ispin)) + + call ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOOt,nVVt,X1t,Y1t,X2t,Y2t,Ec(iAC,ispin)) + + write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) + + end do + + EcAC(ispin) = 1.5d0*dot_product(wAC,Ec(:,ispin)) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) + + end if + +end subroutine ACFDT_pp diff --git a/src/RPA/ACFDT_pp_correlation_energy.f90 b/src/RPA/ACFDT_pp_correlation_energy.f90 new file mode 100644 index 0000000..28991f8 --- /dev/null +++ b/src/RPA/ACFDT_pp_correlation_energy.f90 @@ -0,0 +1,51 @@ +subroutine ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y1,X2,Y2,EcAC) + +! Compute the correlation energy via the adiabatic connection formula for the pp sector + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + +! Local variables + + double precision,allocatable :: B(:,:) + double precision,allocatable :: C(:,:) + double precision,allocatable :: D(:,:) + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: EcAC + +! Memory allocation + + allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO)) + +! Build pp matrices + + call linear_response_B_pp (ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,B) + call linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,C) + call linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,D) + +! Compute Tr(K x P_lambda) + + EcAC = trace_matrix(nVV,matmul(transpose(X1),matmul(B,Y1)) + matmul(transpose(Y1),matmul(transpose(B),X1))) & + + trace_matrix(nVV,matmul(transpose(X1),matmul(C,X1)) + matmul(transpose(Y1),matmul(D,Y1))) & + + trace_matrix(nOO,matmul(transpose(X2),matmul(B,Y2)) + matmul(transpose(Y2),matmul(transpose(B),X2))) & + + trace_matrix(nOO,matmul(transpose(X2),matmul(C,X2)) + matmul(transpose(Y2),matmul(D,Y2))) & + - trace_matrix(nVV,C) - trace_matrix(nOO,D) + +end subroutine ACFDT_pp_correlation_energy + diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index 9066b10..3e0b18a 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -1,4 +1,4 @@ -subroutine ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) +subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Perform pp-RPA calculation @@ -9,10 +9,8 @@ subroutine ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n logical,intent(in) :: TDA logical,intent(in) :: doACFDT - logical,intent(in) :: exchange_kernel logical,intent(in) :: singlet logical,intent(in) :: triplet - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -116,15 +114,7 @@ subroutine ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n write(*,*) '---------------------------------------------------------' write(*,*) - call ACFDT_Tmatrix(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - ERI,e,e,EcAC) - - if(exchange_kernel) then - - EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(1) - - end if + call ACFDT_pp(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------'