diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90 index 8569885..12a1f90 100644 --- a/src/GT/RG0T0pp.f90 +++ b/src/GT/RG0T0pp.f90 @@ -40,11 +40,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Local variables - logical :: print_T = .false. + logical :: print_T = .true. logical :: plot_self = .false. integer :: isp_T - integer :: iblock integer :: nOOs,nOOt integer :: nVVs,nVVt integer :: n_states, n_states_diag @@ -90,11 +89,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Dimensions of the pp-RPA linear reponse matrices - !nOOs = nO*(nO + 1)/2 - !nVVs = nV*(nV + 1)/2 - - 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 @@ -114,60 +110,54 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d !---------------------------------------------- isp_T = 1 - !iblock = 1 - iblock = 3 ! Compute linear response allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp) + call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) deallocate(Bpp,Cpp,Dpp) - if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s) - if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s) + if(print_T) call print_excitation_energies('ppRPA@RHF','2p (singlets)',nVVs,Om1s) + if(print_T) call print_excitation_energies('ppRPA@RHF','2h (singlets)',nOOs,Om2s) !---------------------------------------------- ! alpha-alpha block !---------------------------------------------- isp_T = 2 -! iblock = 2 - iblock = 4 ! Compute linear response allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp) + call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) deallocate(Bpp,Cpp,Dpp) - if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t) - if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-alpha)',nOOt,Om2t) + if(print_T) call print_excitation_energies('ppRPA@RHF','2p (triplets)',nVVt,Om1t) + if(print_T) call print_excitation_energies('ppRPA@RHF','2h (triplets)',nOOt,Om2t) !---------------------------------------------- ! Compute excitation densities !---------------------------------------------- -! iblock = 1 - iblock = 3 + isp_T = 1 - call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) + call RGTpp_excitation_density(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) -! iblock = 2 - iblock = 4 + isp_T = 2 - call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + call RGTpp_excitation_density(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) !---------------------------------------------- ! Compute T-matrix version of the self-energy @@ -214,34 +204,30 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Compute the ppRPA correlation energy isp_T = 1 -! iblock = 1 - iblock = 3 allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) + call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) deallocate(Bpp,Cpp,Dpp) isp_T = 2 -! iblock = 2 - iblock = 4 allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) - call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) - call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) + call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) deallocate(Bpp,Cpp,Dpp) - EcRPA(1) = EcRPA(1) - EcRPA(2) + EcRPA(1) = 1d0*EcRPA(1) EcRPA(2) = 3d0*EcRPA(2) call print_RG0T0pp(nOrb,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA) diff --git a/src/GT/RGTpp_excitation_density.f90 b/src/GT/RGTpp_excitation_density.f90 index e894d02..8a853f1 100644 --- a/src/GT/RGTpp_excitation_density.f90 +++ b/src/GT/RGTpp_excitation_density.f90 @@ -2,12 +2,6 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho ! Compute excitation densities for T-matrix self-energy - ! TODO - ! debug DGEMM for nC != 0 - ! and nR != 0 - - - implicit none ! Input variables @@ -72,8 +66,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho cd = cd + 1 rho1(p,q,ab) = rho1(p,q,ab) & + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/ & - (1d0 + Kronecker_delta(c,d)) -! sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d))) + sqrt(1d0 + Kronecker_delta(c,d))/sqrt(2d0) end do end do @@ -83,8 +76,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho kl = kl + 1 rho1(p,q,ab) = rho1(p,q,ab) & + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/ & - (1d0 + Kronecker_delta(k,l)) -! sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l))) + sqrt(1d0 + Kronecker_delta(k,l))/sqrt(2d0) end do end do @@ -101,8 +93,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho cd = cd + 1 rho2(p,q,ij) = rho2(p,q,ij) & + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/ & - (1d0 + Kronecker_delta(c,d)) -! sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d))) + sqrt(1d0 + Kronecker_delta(c,d))/sqrt(2d0) end do end do @@ -112,8 +103,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho kl = kl + 1 rho2(p,q,ij) = rho2(p,q,ij) & + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/ & - (1d0 + Kronecker_delta(k,l)) -! sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l))) + sqrt(1d0 + Kronecker_delta(k,l))/sqrt(2d0) end do end do @@ -158,7 +148,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho cd = cd + 1 rho1(p,q,ab) = rho1(p,q,ab) & - + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) + + sqrt(3d0/2d0)*(ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) end do ! d end do ! c @@ -169,7 +159,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho kl = kl + 1 rho1(p,q,ab) = rho1(p,q,ab) & - + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) + + sqrt(3d0/2d0)*(ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) end do ! l end do ! k end do ! b @@ -189,7 +179,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho cd = cd + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) + + sqrt(3d0/2d0)*(ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) end do ! d end do ! c @@ -200,7 +190,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho kl = kl + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) + + sqrt(3d0/2d0)*(ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) end do ! l end do ! k end do ! j @@ -240,7 +230,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho enddo !$OMP END DO !$OMP END PARALLEL - + call dgemm("N", "N", nBas*nBas, dim_1, dim_1, 1.d0, & ERI_1(1,1,1), nBas*nBas, X1(1,1), dim_1, & 0.d0, rho1(1,1,1), nBas*nBas) @@ -258,6 +248,9 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho 1.d0, rho2(1,1,1), nBas*nBas) deallocate(ERI_1, ERI_2) + + rho1 = rho1*sqrt(3d0/2d0) + rho2 = rho2*sqrt(3d0/2d0) endif endif