From c22278f497c0b1347d2e97ad2953c47aaf88a0f0 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 20 Oct 2019 08:18:58 +0200 Subject: [PATCH] singlet and triplet G0T0 --- input/basis | 100 +++++------ input/molecule | 8 +- input/weight | 100 +++++------ src/QuAcK/G0T0.f90 | 106 ++++++++---- src/QuAcK/excitation_density_Tmatrix.f90 | 167 +++++++++++++------ src/QuAcK/print_G0T0.f90 | 18 +- src/QuAcK/renormalization_factor_Tmatrix.f90 | 100 +++++++++-- src/QuAcK/self_energy_Tmatrix_diag.f90 | 64 +++++-- 8 files changed, 423 insertions(+), 240 deletions(-) diff --git a/input/basis b/input/basis index 11cb61a..b18e129 100644 --- a/input/basis +++ b/input/basis @@ -1,76 +1,52 @@ 1 9 S 9 1.00 - 1.471000D+04 7.210000D-04 - 2.207000D+03 5.553000D-03 - 5.028000D+02 2.826700D-02 - 1.426000D+02 1.064440D-01 - 4.647000D+01 2.868140D-01 - 1.670000D+01 4.486410D-01 - 6.356000D+00 2.647610D-01 - 1.316000D+00 1.533300D-02 - 3.897000D-01 -2.332000D-03 + 1.469000D+03 7.660000D-04 + 2.205000D+02 5.892000D-03 + 5.026000D+01 2.967100D-02 + 1.424000D+01 1.091800D-01 + 4.581000D+00 2.827890D-01 + 1.580000D+00 4.531230D-01 + 5.640000D-01 2.747740D-01 + 7.345000D-02 9.751000D-03 + 2.805000D-02 -3.180000D-03 S 9 1.00 - 1.471000D+04 -1.650000D-04 - 2.207000D+03 -1.308000D-03 - 5.028000D+02 -6.495000D-03 - 1.426000D+02 -2.669100D-02 - 4.647000D+01 -7.369000D-02 - 1.670000D+01 -1.707760D-01 - 6.356000D+00 -1.123270D-01 - 1.316000D+00 5.628140D-01 - 3.897000D-01 5.687780D-01 + 1.469000D+03 -1.200000D-04 + 2.205000D+02 -9.230000D-04 + 5.026000D+01 -4.689000D-03 + 1.424000D+01 -1.768200D-02 + 4.581000D+00 -4.890200D-02 + 1.580000D+00 -9.600900D-02 + 5.640000D-01 -1.363800D-01 + 7.345000D-02 5.751020D-01 + 2.805000D-02 5.176610D-01 S 1 1.00 - 3.897000D-01 1.000000D+00 + 2.805000D-02 1.000000D+00 S 1 1.00 - 0.0986300 1.0000000 + 8.600000D-03 1.000000D+00 P 4 1.00 - 2.267000D+01 4.487800D-02 - 4.977000D+00 2.357180D-01 - 1.347000D+00 5.085210D-01 - 3.471000D-01 4.581200D-01 + 1.534000D+00 2.278400D-02 + 2.749000D-01 1.391070D-01 + 7.362000D-02 5.003750D-01 + 2.403000D-02 5.084740D-01 P 1 1.00 - 3.471000D-01 1.000000D+00 + 2.403000D-02 1.000000D+00 P 1 1.00 - 0.0850200 1.0000000 + 5.800000D-03 1.000000D+00 D 1 1.00 - 1.640000D+00 1.0000000 + 1.144000D-01 1.0000000 D 1 1.00 - 0.4640000 1.0000000 -2 9 -S 9 1.00 - 1.471000D+04 7.210000D-04 - 2.207000D+03 5.553000D-03 - 5.028000D+02 2.826700D-02 - 1.426000D+02 1.064440D-01 - 4.647000D+01 2.868140D-01 - 1.670000D+01 4.486410D-01 - 6.356000D+00 2.647610D-01 - 1.316000D+00 1.533300D-02 - 3.897000D-01 -2.332000D-03 -S 9 1.00 - 1.471000D+04 -1.650000D-04 - 2.207000D+03 -1.308000D-03 - 5.028000D+02 -6.495000D-03 - 1.426000D+02 -2.669100D-02 - 4.647000D+01 -7.369000D-02 - 1.670000D+01 -1.707760D-01 - 6.356000D+00 -1.123270D-01 - 1.316000D+00 5.628140D-01 - 3.897000D-01 5.687780D-01 + 7.330000D-02 1.000000D+00 +2 5 +S 4 1.00 + 1.301000D+01 1.968500D-02 + 1.962000D+00 1.379770D-01 + 4.446000D-01 4.781480D-01 + 1.220000D-01 5.012400D-01 S 1 1.00 - 3.897000D-01 1.000000D+00 + 1.220000D-01 1.000000D+00 S 1 1.00 - 0.0986300 1.0000000 -P 4 1.00 - 2.267000D+01 4.487800D-02 - 4.977000D+00 2.357180D-01 - 1.347000D+00 5.085210D-01 - 3.471000D-01 4.581200D-01 + 0.0297400 1.0000000 P 1 1.00 - 3.471000D-01 1.000000D+00 + 7.270000D-01 1.0000000 P 1 1.00 - 0.0850200 1.0000000 -D 1 1.00 - 1.640000D+00 1.0000000 -D 1 1.00 - 0.4640000 1.0000000 + 0.1410000 1.0000000 diff --git a/input/molecule b/input/molecule index 7ef2def..c38d7a8 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,5 @@ -# nAt nEl nCore nRyd - 2 9 9 0 0 +# nAt nEla nElb nCore nRyd + 2 2 2 0 0 # Znuc x y z - F 0. 0. 0.000 - F 0. 0. 2.68454493 + Li 0. 0. 0. + H 0. 0. 3.09839495 diff --git a/input/weight b/input/weight index 11cb61a..b18e129 100644 --- a/input/weight +++ b/input/weight @@ -1,76 +1,52 @@ 1 9 S 9 1.00 - 1.471000D+04 7.210000D-04 - 2.207000D+03 5.553000D-03 - 5.028000D+02 2.826700D-02 - 1.426000D+02 1.064440D-01 - 4.647000D+01 2.868140D-01 - 1.670000D+01 4.486410D-01 - 6.356000D+00 2.647610D-01 - 1.316000D+00 1.533300D-02 - 3.897000D-01 -2.332000D-03 + 1.469000D+03 7.660000D-04 + 2.205000D+02 5.892000D-03 + 5.026000D+01 2.967100D-02 + 1.424000D+01 1.091800D-01 + 4.581000D+00 2.827890D-01 + 1.580000D+00 4.531230D-01 + 5.640000D-01 2.747740D-01 + 7.345000D-02 9.751000D-03 + 2.805000D-02 -3.180000D-03 S 9 1.00 - 1.471000D+04 -1.650000D-04 - 2.207000D+03 -1.308000D-03 - 5.028000D+02 -6.495000D-03 - 1.426000D+02 -2.669100D-02 - 4.647000D+01 -7.369000D-02 - 1.670000D+01 -1.707760D-01 - 6.356000D+00 -1.123270D-01 - 1.316000D+00 5.628140D-01 - 3.897000D-01 5.687780D-01 + 1.469000D+03 -1.200000D-04 + 2.205000D+02 -9.230000D-04 + 5.026000D+01 -4.689000D-03 + 1.424000D+01 -1.768200D-02 + 4.581000D+00 -4.890200D-02 + 1.580000D+00 -9.600900D-02 + 5.640000D-01 -1.363800D-01 + 7.345000D-02 5.751020D-01 + 2.805000D-02 5.176610D-01 S 1 1.00 - 3.897000D-01 1.000000D+00 + 2.805000D-02 1.000000D+00 S 1 1.00 - 0.0986300 1.0000000 + 8.600000D-03 1.000000D+00 P 4 1.00 - 2.267000D+01 4.487800D-02 - 4.977000D+00 2.357180D-01 - 1.347000D+00 5.085210D-01 - 3.471000D-01 4.581200D-01 + 1.534000D+00 2.278400D-02 + 2.749000D-01 1.391070D-01 + 7.362000D-02 5.003750D-01 + 2.403000D-02 5.084740D-01 P 1 1.00 - 3.471000D-01 1.000000D+00 + 2.403000D-02 1.000000D+00 P 1 1.00 - 0.0850200 1.0000000 + 5.800000D-03 1.000000D+00 D 1 1.00 - 1.640000D+00 1.0000000 + 1.144000D-01 1.0000000 D 1 1.00 - 0.4640000 1.0000000 -2 9 -S 9 1.00 - 1.471000D+04 7.210000D-04 - 2.207000D+03 5.553000D-03 - 5.028000D+02 2.826700D-02 - 1.426000D+02 1.064440D-01 - 4.647000D+01 2.868140D-01 - 1.670000D+01 4.486410D-01 - 6.356000D+00 2.647610D-01 - 1.316000D+00 1.533300D-02 - 3.897000D-01 -2.332000D-03 -S 9 1.00 - 1.471000D+04 -1.650000D-04 - 2.207000D+03 -1.308000D-03 - 5.028000D+02 -6.495000D-03 - 1.426000D+02 -2.669100D-02 - 4.647000D+01 -7.369000D-02 - 1.670000D+01 -1.707760D-01 - 6.356000D+00 -1.123270D-01 - 1.316000D+00 5.628140D-01 - 3.897000D-01 5.687780D-01 + 7.330000D-02 1.000000D+00 +2 5 +S 4 1.00 + 1.301000D+01 1.968500D-02 + 1.962000D+00 1.379770D-01 + 4.446000D-01 4.781480D-01 + 1.220000D-01 5.012400D-01 S 1 1.00 - 3.897000D-01 1.000000D+00 + 1.220000D-01 1.000000D+00 S 1 1.00 - 0.0986300 1.0000000 -P 4 1.00 - 2.267000D+01 4.487800D-02 - 4.977000D+00 2.357180D-01 - 1.347000D+00 5.085210D-01 - 3.471000D-01 4.581200D-01 + 0.0297400 1.0000000 P 1 1.00 - 3.471000D-01 1.000000D+00 + 7.270000D-01 1.0000000 P 1 1.00 - 0.0850200 1.0000000 -D 1 1.00 - 1.640000D+00 1.0000000 -D 1 1.00 - 0.4640000 1.0000000 + 0.1410000 1.0000000 diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index 301b31d..2c784d9 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -21,18 +21,18 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc, ! Local variables integer :: ispin - integer :: nOO - integer :: nVV + integer :: nOOs,nOOt + integer :: nVVs,nVVt double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) - double precision,allocatable :: Omega1(:,:) - double precision,allocatable :: X1(:,:,:) - double precision,allocatable :: Y1(:,:,:) - double precision,allocatable :: rho1(:,:,:,:) - double precision,allocatable :: Omega2(:,:) - double precision,allocatable :: X2(:,:,:) - double precision,allocatable :: Y2(:,:,:) - double precision,allocatable :: rho2(:,:,:,:) + 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(:,:,:) double precision,allocatable :: SigT(:) double precision,allocatable :: Z(:) @@ -48,57 +48,97 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc, write(*,*)'************************************************' write(*,*) -! Spin manifold +! Dimensions of the rr-RPA linear reponse matrices - ispin = 1 + nOOs = nO*(nO+1)/2 + nVVs = nV*(nV+1)/2 - nOO = nO*(nO+1)/2 - nVV = nV*(nV+1)/2 + nOOt = nO*(nO-1)/2 + nVVt = nV*(nV-1)/2 ! Memory allocation - allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & - Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin), & - rho1(nBas,nBas,nVV,nspin),rho2(nBas,nBas,nOO,nspin), & + 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(nOOs),X2t(nVVs,nOOs),Y2t(nOOs,nOOs), & + rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), & SigT(nBas),Z(nBas)) +!---------------------------------------------- +! Singlet manifold +!---------------------------------------------- + + ispin = 1 + ! Compute linear response - call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,eHF(:),ERI(:,:,:,:), & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR, & + nOOs,nVVs,eHF(:),ERI(:,:,:,:), & + Omega1s(:),X1s(:,:),Y1s(:,:), & + Omega2s(:),X2s(:,:),Y2s(:,:), & EcRPA(ispin)) + call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s(:)) + call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s(:)) + ! Compute excitation densities for the T-matrix - call excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI(:,:,:,:), & - X1(:,:,ispin),Y1(:,:,ispin),rho1(:,:,:,ispin), & - X2(:,:,ispin),Y2(:,:,ispin),rho2(:,:,:,ispin)) + call excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOOs,nVVs,ERI(:,:,:,:), & + X1s(:,:),Y1s(:,:),rho1s(:,:,:), & + X2s(:,:),Y2s(:,:),rho2s(:,:,:)) +!---------------------------------------------- +! Triplet manifold +!---------------------------------------------- + + ispin = 2 + + ! Compute linear response + + call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR, & + nOOt,nVVt,eHF(:),ERI(:,:,:,:), & + Omega1t(:),X1t(:,:),Y1t(:,:), & + Omega2t(:),X2t(:,:),Y2t(:,:), & + EcRPA(ispin)) + + call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t(:)) + call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t(:)) + + ! Compute excitation densities for the T-matrix + + call excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOOt,nVVt,ERI(:,:,:,:), & + X1t(:,:),Y1t(:,:),rho1t(:,:,:), & + X2t(:,:),Y2t(:,:),rho2t(:,:,:)) + +!---------------------------------------------- ! Compute T-matrix version of the self-energy +!---------------------------------------------- - call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF(:), & - Omega1(:,ispin),rho1(:,:,:,ispin), & - Omega2(:,ispin),rho2(:,:,:,ispin), & + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & + Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & + Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & SigT(:)) ! Compute renormalization factor for T-matrix self-energy - call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,eHF(:), & - Omega1(:,ispin),rho1(:,:,:,ispin), & - Omega2(:,ispin),rho2(:,:,:,ispin), & + call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & + Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & + Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & Z(:)) +!---------------------------------------------- ! Solve the quasi-particle equation +!---------------------------------------------- eG0T0(:) = eHF(:) + Z(:)*SigT(:) +!---------------------------------------------- ! Dump results +!---------------------------------------------- - call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) - call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) - - call print_G0T0(nBas,nO,eHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA(ispin)) + call print_G0T0(nBas,nO,eHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA(:)) ! Perform BSE calculation diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index 9b594d6..751690d 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -1,4 +1,4 @@ -subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) +subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) ! Compute excitation densities for T-matrix self-energy @@ -6,6 +6,7 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2 ! Input variables + integer,intent(in) :: ispin integer,intent(in) :: nBas,nC,nO,nR,nOO,nVV double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(out) :: X1(nVV,nVV) @@ -29,60 +30,132 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2 rho1(:,:,:) = 0d0 rho2(:,:,:) = 0d0 - do p=nC+1,nBas-nR +!---------------------------------------------- +! Singlet manifold +!---------------------------------------------- - do i=nC+1,nO - do ab=1,nVV + if(ispin == 1) then - cd = 0 - do c=nO+1,nBas-nR - do d=nO+1,c - cd = cd + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + 1.0d0*(ERI(p,i,c,d) - 1.0d0*ERI(p,i,d,c))*X1(cd,ab)!& -! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) - enddo - enddo + do p=nC+1,nBas-nR - kl = 0 - do k=nC+1,nO - do l=nC+1,k - kl = kl + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + 1.0d0*(ERI(p,i,k,l) - 1.0d0*ERI(p,i,l,k))*Y1(kl,ab)!& -! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) - enddo - enddo + do i=nC+1,nO + do ab=1,nVV - enddo - enddo + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,c + cd = cd + 1 + rho1(p,i,ab) = rho1(p,i,ab) & + + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) +! + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) & +! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) + end do + end do - do a=nO+1,nBas-nR - do ij=1,nOO + kl = 0 + do k=nC+1,nO + do l=nC+1,k + kl = kl + 1 + rho1(p,i,ab) = rho1(p,i,ab) & + + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) +! + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) & +! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) + end do + end do - cd = 0 - do c=nO+1,nBas-nR - do d=nO+1,c - cd = cd + 1 - rho2(p,a,ij) = rho2(p,a,ij) & - + 1.0d0*(ERI(p,a,c,d) - 1.0d0*ERI(p,a,d,c))*X2(cd,ij)!& -! /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(c,d))) - enddo - enddo + end do + end do - kl = 0 - do k=nC+1,nO - do l=nC+1,k - kl = kl + 1 - rho2(p,a,ij) = rho2(p,a,ij) & - + 1.0d0*(ERI(p,a,k,l) - 1.0d0*ERI(p,a,l,k))*Y2(kl,ij)!& -! /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(k,l))) - enddo - enddo + do a=nO+1,nBas-nR + do ij=1,nOO - enddo - enddo - enddo + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,c + cd = cd + 1 + rho2(p,a,ij) = rho2(p,a,ij) & + + (ERI(p,a,c,d) - ERI(p,a,d,c))*X2(cd,ij) +! + (ERI(p,a,c,d) + ERI(p,a,d,c))*X2(cd,ij) & +! /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(c,d))) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=nC+1,k + kl = kl + 1 + rho2(p,a,ij) = rho2(p,a,ij) & + + (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij) +! + (ERI(p,a,k,l) + ERI(p,a,l,k))*Y2(kl,ij) & +! /sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(k,l))) + + end do + end do + + end do + end do + end do + + end if + +!---------------------------------------------- +! Triplet manifold +!---------------------------------------------- + + if(ispin == 2) then + + do p=nC+1,nBas-nR + + do i=nC+1,nO + do ab=1,nVV + + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,c-1 + cd = cd + 1 + rho1(p,i,ab) = rho1(p,i,ab) & + + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=nC+1,k-1 + kl = kl + 1 + rho1(p,i,ab) = rho1(p,i,ab) & + + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) + end do + end do + + end do + end do + + do a=nO+1,nBas-nR + do ij=1,nOO + + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,c-1 + cd = cd + 1 + rho2(p,a,ij) = rho2(p,a,ij) & + + (ERI(p,a,c,d) - ERI(p,a,d,c))*X2(cd,ij) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=nC+1,k-1 + kl = kl + 1 + rho2(p,a,ij) = rho2(p,a,ij) & + + (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij) + end do + end do + + end do + end do + end do + + end if end subroutine excitation_density_Tmatrix diff --git a/src/QuAcK/print_G0T0.f90 b/src/QuAcK/print_G0T0.f90 index 4cea81e..0d4f025 100644 --- a/src/QuAcK/print_G0T0.f90 +++ b/src/QuAcK/print_G0T0.f90 @@ -1,4 +1,4 @@ -subroutine print_G0T0(nBas,nO,e,ENuc,EHF,SigT,Z,eGW,EcRPA) +subroutine print_G0T0(nBas,nO,e,ENuc,ERHF,SigT,Z,eGW,EcRPA) ! Print one-electron energies and other stuff for G0T0 @@ -7,8 +7,8 @@ subroutine print_G0T0(nBas,nO,e,ENuc,EHF,SigT,Z,eGW,EcRPA) integer,intent(in) :: nBas,nO double precision,intent(in) :: ENuc - double precision,intent(in) :: EHF - double precision,intent(in) :: EcRPA + double precision,intent(in) :: ERHF + double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: e(nBas),SigT(nBas),Z(nBas),eGW(nBas) integer :: x,HOMO,LUMO @@ -35,12 +35,14 @@ subroutine print_G0T0(nBas,nO,e,ENuc,EHF,SigT,Z,eGW,EcRPA) enddo write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'G0T0 HOMO energy (eV):',eGW(HOMO)*HaToeV - write(*,'(2X,A30,F15.6)') 'G0T0 LUMO energy (eV):',eGW(LUMO)*HaToeV - write(*,'(2X,A30,F15.6)') 'G0T0 HOMO-LUMO gap (eV):',Gap*HaToeV + write(*,'(2X,A40,F15.6)') 'G0T0 HOMO energy (eV) :',eGW(HOMO)*HaToeV + write(*,'(2X,A40,F15.6)') 'G0T0 LUMO energy (eV) :',eGW(LUMO)*HaToeV + write(*,'(2X,A40,F15.6)') 'G0T0 HOMO-LUMO gap (eV) :',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'RPA@G0T0 total energy =',ENuc + EHF + EcRPA - write(*,'(2X,A30,F15.6)') 'RPA@G0T0 correlation energy =',EcRPA + write(*,'(2X,A40,F15.6)') 'RPA@G0T0 correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A40,F15.6)') 'RPA@G0T0 correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A40,F15.6)') 'RPA@G0T0 correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A40,F15.6)') 'RPA@G0T0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90 index 71466cc..118ac64 100644 --- a/src/QuAcK/renormalization_factor_Tmatrix.f90 +++ b/src/QuAcK/renormalization_factor_Tmatrix.f90 @@ -1,4 +1,6 @@ -subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,Z) +subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, & + Omega1s,rho1s,Omega2s,rho2s,Omega1t,rho1t,Omega2t,rho2t, & + Z) ! Compute renormalization factor of the T-matrix self-energy @@ -8,12 +10,14 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1, ! Input variables double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + integer,intent(in) :: nBas,nC,nO,nV,nR + integer,intent(in) :: nOOs,nVVs + integer,intent(in) :: nOOt,nVVt double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Omega2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: Omega1s(nVVs),Omega1t(nVVt) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: Omega2s(nOOs),Omega2t(nOOt) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt) ! Local variables @@ -28,7 +32,11 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1, Z(:) = 0d0 - ! Occupied part of the T-matrix self-energy +!---------------------------------------------- +! Singlet part of the T-matrix self-energy +!---------------------------------------------- + +! Occupied part of the T-matrix self-energy do p=nC+1,nBas-nR do i=nC+1,nO @@ -36,14 +44,14 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1, do c=nO+1,nBas-nR do d=nO+1,c cd = cd + 1 - eps = e(p) + e(i) - Omega1(cd) - Z(p) = Z(p) - 2d0*rho1(p,i,cd)**2/eps**2 + eps = e(p) + e(i) - Omega1s(cd) + Z(p) = Z(p) - 2d0*rho1s(p,i,cd)**2/eps**2 enddo enddo enddo enddo - ! Virtual part of the T-matrix self-energy +! Virtual part of the T-matrix self-energy do p=nC+1,nBas-nR do a=nO+1,nBas-nR @@ -51,8 +59,76 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1, do k=nC+1,nO do l=nC+1,k kl = kl + 1 - eps = e(p) + e(a) - Omega2(kl) - Z(p) = Z(p) - 2d0*rho2(p,a,kl)**2/eps**2 + eps = e(p) + e(a) - Omega2s(kl) + Z(p) = Z(p) - 2d0*rho2s(p,a,kl)**2/eps**2 + enddo + enddo + enddo + enddo + +!---------------------------------------------- +! Triplet part of the T-matrix self-energy +!---------------------------------------------- + +! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,c-1 + cd = cd + 1 + eps = e(p) + e(i) - Omega1t(cd) + Z(p) = Z(p) - 2d0*rho1t(p,i,cd)**2/eps**2 + enddo + enddo + enddo + enddo + +! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + kl = 0 + do k=nC+1,nO + do l=nC+1,k-1 + kl = kl + 1 + eps = e(p) + e(a) - Omega2t(kl) + Z(p) = Z(p) - 2d0*rho2t(p,a,kl)**2/eps**2 + enddo + enddo + enddo + enddo + +!---------------------------------------------- +! Singlet part of the T-matrix self-energy +!---------------------------------------------- + +! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,c + cd = cd + 1 + eps = e(p) + e(i) - Omega1s(cd) + Z(p) = Z(p) - 2d0*rho1s(p,i,cd)**2/eps**2 + enddo + enddo + enddo + enddo + +! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + kl = 0 + do k=nC+1,nO + do l=nC+1,k + kl = kl + 1 + eps = e(p) + e(a) - Omega2s(kl) + Z(p) = Z(p) - 2d0*rho2s(p,a,kl)**2/eps**2 enddo enddo enddo diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90 index 73eb8dc..c01784b 100644 --- a/src/QuAcK/self_energy_Tmatrix_diag.f90 +++ b/src/QuAcK/self_energy_Tmatrix_diag.f90 @@ -1,4 +1,6 @@ -subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT) +subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, & + Omega1s,rho1s,Omega2s,rho2s,Omega1t,rho1t,Omega2t,rho2t, & + SigT) ! Compute diagonal of the correlation part of the T-matrix self-energy @@ -13,13 +15,13 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR - integer,intent(in) :: nOO - integer,intent(in) :: nVV + integer,intent(in) :: nOOs,nOOt + integer,intent(in) :: nVVs,nVVt double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega1(nVV) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: Omega2(nOO) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: Omega1s(nVVs),Omega1t(nVVt) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs),rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: Omega2s(nOOs),Omega2t(nOOt) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs),rho2t(nBas,nBas,nOOt) ! Local variables @@ -34,7 +36,11 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O SigT(:) = 0d0 - ! Occupied part of the T-matrix self-energy +!---------------------------------------------- +! Singlet part of the T-matrix self-energy +!---------------------------------------------- + +! Occupied part of the T-matrix self-energy do p=nC+1,nBas-nR do i=nC+1,nO @@ -42,8 +48,8 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O do c=nO+1,nBas-nR do d=nO+1,c cd = cd + 1 - eps = e(p) + e(i) - Omega1(cd) - SigT(p) = SigT(p) + 2d0*rho1(p,i,cd)**2/eps + eps = e(p) + e(i) - Omega1s(cd) + SigT(p) = SigT(p) + 2d0*rho1s(p,i,cd)**2/eps enddo enddo enddo @@ -57,8 +63,42 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O do k=nC+1,nO do l=nC+1,k kl = kl + 1 - eps = e(p) + e(a) - Omega2(kl) - SigT(p) = SigT(p) + 2d0*rho2(p,a,kl)**2/eps + eps = e(p) + e(a) - Omega2s(kl) + SigT(p) = SigT(p) + 2d0*rho2s(p,a,kl)**2/eps + enddo + enddo + enddo + enddo + +!---------------------------------------------- +! Triplet part of the T-matrix self-energy +!---------------------------------------------- + +! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,c-1 + cd = cd + 1 + eps = e(p) + e(i) - Omega1t(cd) + SigT(p) = SigT(p) + 2d0*rho1t(p,i,cd)**2/eps + enddo + enddo + enddo + enddo + + ! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + kl = 0 + do k=nC+1,nO + do l=nC+1,k-1 + kl = kl + 1 + eps = e(p) + e(a) - Omega2t(kl) + SigT(p) = SigT(p) + 2d0*rho2t(p,a,kl)**2/eps enddo enddo enddo