diff --git a/input/options b/input/options index 107b2f4..681a821 100644 --- a/input/options +++ b/input/options @@ -11,8 +11,8 @@ # 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 - F F T + T F T # BSE: BSE dBSE dTDA evDyn - F T T F + T F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/LR/linear_response_B_pp.f90 b/src/LR/linear_response_B_pp.f90 index 7644897..3760373 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,e,ERI,B_pp) +subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B_pp) ! Compute the B matrix of the pp channel @@ -9,6 +9,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,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) ! Local variables @@ -33,7 +34,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) do j=i,nO ij = ij + 1 - B_pp(ab,ij) = (ERI(a,b,i,j) + ERI(a,b,j,i))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) + B_pp(ab,ij) = lambda*(ERI(a,b,i,j) + ERI(a,b,j,i))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) end do end do @@ -55,7 +56,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) do j=nC+1,nO ij = ij + 1 - B_pp(ab,ij) = ERI(a,b,i,j) + B_pp(ab,ij) = lambda*ERI(a,b,i,j) end do end do @@ -77,7 +78,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) do j=i+1,nO ij = ij + 1 - B_pp(ab,ij) = ERI(a,b,i,j) - ERI(a,b,j,i) + B_pp(ab,ij) = lambda*(ERI(a,b,i,j) - ERI(a,b,j,i)) end do end do @@ -86,7 +87,4 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) end if -! print*,'B pp-matrix' -! call matout(nVV,nOO,B_pp) - end subroutine linear_response_B_pp diff --git a/src/LR/linear_response_C_pp.f90 b/src/LR/linear_response_C_pp.f90 index a8b391d..724dc59 100644 --- a/src/LR/linear_response_C_pp.f90 +++ b/src/LR/linear_response_C_pp.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) +subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp) ! Compute the C matrix of the pp channel @@ -9,6 +9,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_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) ! Local variables @@ -41,7 +42,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) cd = cd + 1 C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & - + (ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + + 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 @@ -64,7 +65,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) cd = cd + 1 C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & - + ERI(a,b,c,d) - ERI(a,b,d,c) + + lambda*(ERI(a,b,c,d) - ERI(a,b,d,c)) end do end do @@ -87,7 +88,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) cd = cd + 1 C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & - + ERI(a,b,c,d) + + lambda*ERI(a,b,c,d) end do end do @@ -96,7 +97,4 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) end if -! print*,'C pp-matrix' -! call matout(nVV,nVV,C_pp) - end subroutine linear_response_C_pp diff --git a/src/LR/linear_response_D_pp.f90 b/src/LR/linear_response_D_pp.f90 index 04a2960..398b96a 100644 --- a/src/LR/linear_response_D_pp.f90 +++ b/src/LR/linear_response_D_pp.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) +subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp) ! Compute the D matrix of the pp channel @@ -9,6 +9,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_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) ! Local variables @@ -41,7 +42,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) kl = kl + 1 D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - + (ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + + 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 @@ -64,7 +65,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) kl = kl + 1 D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - + ERI(i,j,k,l) - ERI(i,j,l,k) + + lambda*(ERI(i,j,k,l) - ERI(i,j,l,k)) end do end do @@ -87,7 +88,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) kl = kl + 1 D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - + ERI(i,j,k,l) + + lambda*ERI(i,j,k,l) end do end do @@ -96,7 +97,4 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) end if -! print*,'D pp-matrix' -! call matout(nOO,nOO,D_pp) - end subroutine linear_response_D_pp diff --git a/src/LR/linear_response_pp.f90 b/src/LR/linear_response_pp.f90 index 4b35ec0..7a14be8 100644 --- a/src/LR/linear_response_pp.f90 +++ b/src/LR/linear_response_pp.f90 @@ -1,5 +1,4 @@ -subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV, & - e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcRPA) +subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcRPA) ! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013) @@ -17,6 +16,7 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV, & integer,intent(in) :: nR integer,intent(in) :: nOO integer,intent(in) :: nVV + double precision,intent(in) :: lambda double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -60,8 +60,8 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV, & ! Build B, C and D matrices for the pp channel - call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C) - call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D) + call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C) + call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D) if(TDA) then @@ -75,7 +75,7 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV, & else - call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B) + call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B) ! Diagonal blocks diff --git a/src/MBPT/Bethe_Salpeter_Tmatrix.f90 b/src/MBPT/Bethe_Salpeter_Tmatrix.f90 index f16a703..6303e25 100644 --- a/src/MBPT/Bethe_Salpeter_Tmatrix.f90 +++ b/src/MBPT/Bethe_Salpeter_Tmatrix.f90 @@ -83,7 +83,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ispin = 1 iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eT,ERI, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eT,ERI, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) ! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) @@ -98,7 +98,7 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ispin = 2 iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eT,ERI, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eT,ERI, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) ! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) diff --git a/src/MBPT/G0T0.f90 b/src/MBPT/G0T0.f90 index 125cde2..ae946c0 100644 --- a/src/MBPT/G0T0.f90 +++ b/src/MBPT/G0T0.f90 @@ -103,7 +103,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ! Compute linear response - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eHF,ERI_MO, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) ! EcRPA(ispin) = 1d0*EcRPA(ispin) @@ -120,7 +120,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ! Compute linear response - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eHF,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) ! EcRPA(ispin) = 2d0*EcRPA(ispin) @@ -183,11 +183,11 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ispin = 1 iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eG0T0,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eG0T0,ERI_MO, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) ispin = 2 iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eG0T0,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eG0T0,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) diff --git a/src/MBPT/evGT.f90 b/src/MBPT/evGT.f90 index e657bf6..2309e94 100644 --- a/src/MBPT/evGT.f90 +++ b/src/MBPT/evGT.f90 @@ -132,7 +132,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! Compute linear response - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) !---------------------------------------------- @@ -144,7 +144,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! Compute linear response - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) !---------------------------------------------- @@ -222,11 +222,11 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ispin = 1 iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) ispin = 2 iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) diff --git a/src/MBPT/qsGT.f90 b/src/MBPT/qsGT.f90 index 59c43c7..1eeff46 100644 --- a/src/MBPT/qsGT.f90 +++ b/src/MBPT/qsGT.f90 @@ -188,13 +188,13 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T ispin = 1 iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) ispin = 2 iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) ! Compute correlation part of the self-energy @@ -305,11 +305,11 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T ispin = 1 iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) ispin = 2 iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,ERI_MO, & + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 new file mode 100644 index 0000000..711a944 --- /dev/null +++ b/src/RPA/ACFDT_Tmatrix.f90 @@ -0,0 +1,224 @@ +subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t, & + ERI,eT,eGT,EcAC) + +! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: doXBS + logical,intent(in) :: exchange_kernel + logical,intent(in) :: dRPA + logical,intent(in) :: TDA_T + logical,intent(in) :: TDA + logical,intent(in) :: BSE + 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 + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + + integer,intent(in) :: nOOs + integer,intent(in) :: nOOt + integer,intent(in) :: nVVs + integer,intent(in) :: nVVt + + double precision,intent(in) :: eT(nBas) + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + + double precision,intent(in) :: Omega1s(nVVs) + double precision,intent(in) :: X1s(nVVs,nVVs) + double precision,intent(in) :: Y1s(nOOs,nVVs) + double precision,intent(in) :: Omega2s(nOOs) + double precision,intent(in) :: X2s(nVVs,nOOs) + double precision,intent(in) :: Y2s(nOOs,nOOs) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs) + double precision,intent(in) :: Omega1t(nVVt) + double precision,intent(in) :: X1t(nVVt,nVVt) + double precision,intent(in) :: Y1t(nOOt,nVVt) + double precision,intent(in) :: Omega2t(nOOt) + double precision,intent(in) :: X2t(nVVt,nOOt) + double precision,intent(in) :: Y2t(nOOt,nOOt) + double precision,intent(in) :: rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: rho2t(nBas,nBas,nOOt) + +! Local variables + + integer :: ispin + integer :: isp_T + integer :: iblock + integer :: iAC + double precision :: lambda + double precision,allocatable :: Ec(:,:) + + double precision :: EcRPA(nspin) + double precision,allocatable :: TA(:,:) + double precision,allocatable :: TB(:,:) + double precision,allocatable :: Omega(:,:) + double precision,allocatable :: XpY(:,:,:) + double precision,allocatable :: XmY(:,:,:) + +! Output variables + + double precision,intent(out) :: EcAC(nspin) + +! Memory allocation + + allocate(Ec(nAC,nspin)) + allocate(TA(nS,nS),TB(nS,nS),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + +! Antisymmetrized kernel version + + if(exchange_kernel) then + + write(*,*) + write(*,*) '*** Exchange kernel version ***' + write(*,*) + + end if + + 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) + + 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 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) + + 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 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,.true.,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)) + + 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)) + + if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(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) + + 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 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) + + 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 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,.true.,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)) + + 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)) + + if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) + + end if + +end subroutine ACFDT_Tmatrix diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index e5f1ea6..a40e297 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -62,9 +62,9 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) - call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & + Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & + Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & Ec_ppRPA(ispin)) call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) @@ -91,9 +91,9 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) - call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & + Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & + Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & Ec_ppRPA(ispin)) call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))