diff --git a/src/GT/GTpp_excitation_density.f90 b/src/GT/GTpp_excitation_density.f90 index 0cb3a73..a6ed548 100644 --- a/src/GT/GTpp_excitation_density.f90 +++ b/src/GT/GTpp_excitation_density.f90 @@ -53,14 +53,18 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1 do q=nC+1,nBas-nR do p=nC+1,nBas-nR - do ab=1,nVV - + 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 rho1(p,q,ab) = rho1(p,q,ab) & - + ERI(p,q,c,d)*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(c,d))) + + (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))) end do end do @@ -69,20 +73,27 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1 do l=k,nO kl = kl + 1 rho1(p,q,ab) = rho1(p,q,ab) & - + ERI(p,q,k,l)*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(k,l))) + + (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))) end do end do end do + end do - do ij=1,nOO - + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 cd = 0 do c=nO+1,nBas-nR do d=c,nBas-nR cd = cd + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + ERI(p,q,c,d)*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(c,d))) + + (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))) end do end do @@ -91,11 +102,14 @@ subroutine GTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1 do l=k,nO kl = kl + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + ERI(p,q,k,l)*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(k,l))) + + (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))) end do end do end do + end do end do end do diff --git a/src/GT/ufG0T0pp.f90 b/src/GT/ufG0T0pp.f90 index ba498bb..e29c843 100644 --- a/src/GT/ufG0T0pp.f90 +++ b/src/GT/ufG0T0pp.f90 @@ -34,8 +34,8 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) logical :: dRPA integer :: ispin integer :: iblock - integer :: nOO,nOOs,nOOt - integer :: nVV,nVVs,nVVt + integer :: nOOs,nOOt + integer :: nVVs,nVVt double precision :: EcRPA(nspin) integer :: n2h1p,n2p1h,nH double precision,external :: Kronecker_delta @@ -62,6 +62,9 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) double precision :: start_timing,end_timing,timing + double precision :: alpha = sqrt(1d0) + double precision :: beta = sqrt(1d0) + ! Output variables ! Hello world @@ -74,17 +77,17 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Dimensions of the ppRPA linear reponse matrices -! nOOs = nO*(nO + 1)/2 -! nVVs = nV*(nV + 1)/2 + nOOs = nO*(nO + 1)/2 + nVVs = nV*(nV + 1)/2 - nOOs = nO*nO - nVVs = nV*nV +! nOOs = nO*nO +! nVVs = nV*nV nOOt = nO*(nO - 1)/2 nVVt = nV*(nV - 1)/2 - nOO = nO*nO - nVV = nV*nV +! nOO = nO*nO +! nVV = nV*nV ! Dimension of the supermatrix @@ -107,7 +110,7 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Main loop over orbitals ! !-------------------------! - do p=nC+1,nBas-nR + do p=nO,nO+1 H(:,:) = 0d0 @@ -266,8 +269,8 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! alpha-beta block ispin = 1 -! iblock = 1 - iblock = 3 + iblock = 1 +! iblock = 3 ! Compute linear response @@ -287,8 +290,8 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! alpha-alpha block ispin = 2 -! iblock = 2 - iblock = 4 + iblock = 2 +! iblock = 4 ! Compute linear response @@ -309,12 +312,12 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Compute excitation densities !---------------------------------------------- -! iblock = 1 - iblock = 3 + iblock = 1 +! iblock = 3 call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) -! iblock = 2 - iblock = 4 + iblock = 2 +! iblock = 4 call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) call wall_time(start_timing) @@ -357,8 +360,8 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) do a=nO+1,nBas-nR ija = ija + 1 - H(1 ,1+ija) = rho2s(p,a,ij) - H(1+ija,1 ) = rho2s(p,a,ij) + H(1 ,1+ija) = alpha*rho2s(p,a,ij) + H(1+ija,1 ) = alpha*rho2s(p,a,ij) end do end do @@ -367,8 +370,8 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) do a=nO+1,nBas-nR ija = ija + 1 - H(1 ,1+ija) = rho2t(p,a,ij) - H(1+ija,1 ) = rho2t(p,a,ij) + H(1 ,1+ija) = beta*rho2t(p,a,ij) + H(1+ija,1 ) = beta*rho2t(p,a,ij) end do end do @@ -405,8 +408,8 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) do i=nC+1,nO iab = iab + 1 - H(1 ,1+n2h1p+iab) = rho1s(p,i,ab) - H(1+n2h1p+iab,1 ) = rho1s(p,i,ab) + H(1 ,1+n2h1p+iab) = alpha*rho1s(p,i,ab) + H(1+n2h1p+iab,1 ) = alpha*rho1s(p,i,ab) end do end do @@ -415,8 +418,8 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) do i=nC+1,nO iab = iab + 1 - H(1 ,1+n2h1p+iab) = rho1t(p,i,ab) - H(1+n2h1p+iab,1 ) = rho1t(p,i,ab) + H(1 ,1+n2h1p+iab) = beta*rho1t(p,i,ab) + H(1+n2h1p+iab,1 ) = beta*rho1t(p,i,ab) end do end do @@ -564,7 +567,7 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ' (',p,') ',H(1,s),H(1,s)**2 ija = 0 - do ij=1,nOO + do ij=1,nOOs+nOOt do a=nO+1,nBas-nR ija = ija + 1 @@ -576,7 +579,7 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) end do iab = 0 - do ab=1,nVV + do ab=1,nVVs+nVVt do i=nC+1,nO iab = iab + 1