diff --git a/examples/molecule.F2 b/examples/molecule.F2 index e2925d6..6f3bd1d 100644 --- a/examples/molecule.F2 +++ b/examples/molecule.F2 @@ -1,5 +1,5 @@ # nAt nEl nCore nRyd - 2 18 4 0 + 2 9 9 0 0 # Znuc x y z - 9. 0. 0. 0.000 - 9. 0. 0. 2.668 + F 0. 0. 0.000 + F 0. 0. 2.668 diff --git a/input/basis b/input/basis index 4c61da7..27724b4 100644 --- a/input/basis +++ b/input/basis @@ -1,16 +1,37 @@ -1 6 -S 4 1.00 - 234.0000000 0.0025870 - 35.1600000 0.0195330 - 7.9890000 0.0909980 - 2.2120000 0.2720500 +1 10 +S 8 1.00 + 24350.0000000 0.0005020 + 3650.0000000 0.0038810 + 829.6000000 0.0199970 + 234.0000000 0.0784180 + 75.6100000 0.2296760 + 26.7300000 0.4327220 + 9.9270000 0.3506420 + 1.1020000 -0.0076450 +S 8 1.00 + 24350.0000000 -0.0001180 + 3650.0000000 -0.0009150 + 829.6000000 -0.0047370 + 234.0000000 -0.0192330 + 75.6100000 -0.0603690 + 26.7300000 -0.1425080 + 9.9270000 -0.1777100 + 1.1020000 0.6058360 S 1 1.00 - 0.6669000 1.0000000 + 2.8360000 1.0000000 S 1 1.00 - 0.2089000 1.0000000 + 0.3782000 1.0000000 +P 3 1.00 + 54.7000000 0.0171510 + 12.4300000 0.1076560 + 3.6790000 0.3216810 P 1 1.00 - 3.0440000 1.0000000 + 1.1430000 1.0000000 P 1 1.00 - 0.7580000 1.0000000 + 0.3300000 1.0000000 D 1 1.00 - 1.9650000 1.0000000 + 4.0140000 1.0000000 +D 1 1.00 + 1.0960000 1.0000000 +F 1 1.00 + 2.5440000 1.0000000 diff --git a/input/methods b/input/methods index 9ac4f33..513d614 100644 --- a/input/methods +++ b/input/methods @@ -7,7 +7,7 @@ # CIS TDHF ppRPA ADC F F F F # GF2 GF3 - F F + T F # G0W0 evGW qsGW F F F # G0T0 evGT qsGT diff --git a/input/molecule b/input/molecule index c78e87e..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 1 5 5 0 0 # Znuc x y z - He 0.0 0.0 0.0 + Ne 0.0 0.0 0.0 diff --git a/input/weight b/input/weight index 4c61da7..27724b4 100644 --- a/input/weight +++ b/input/weight @@ -1,16 +1,37 @@ -1 6 -S 4 1.00 - 234.0000000 0.0025870 - 35.1600000 0.0195330 - 7.9890000 0.0909980 - 2.2120000 0.2720500 +1 10 +S 8 1.00 + 24350.0000000 0.0005020 + 3650.0000000 0.0038810 + 829.6000000 0.0199970 + 234.0000000 0.0784180 + 75.6100000 0.2296760 + 26.7300000 0.4327220 + 9.9270000 0.3506420 + 1.1020000 -0.0076450 +S 8 1.00 + 24350.0000000 -0.0001180 + 3650.0000000 -0.0009150 + 829.6000000 -0.0047370 + 234.0000000 -0.0192330 + 75.6100000 -0.0603690 + 26.7300000 -0.1425080 + 9.9270000 -0.1777100 + 1.1020000 0.6058360 S 1 1.00 - 0.6669000 1.0000000 + 2.8360000 1.0000000 S 1 1.00 - 0.2089000 1.0000000 + 0.3782000 1.0000000 +P 3 1.00 + 54.7000000 0.0171510 + 12.4300000 0.1076560 + 3.6790000 0.3216810 P 1 1.00 - 3.0440000 1.0000000 + 1.1430000 1.0000000 P 1 1.00 - 0.7580000 1.0000000 + 0.3300000 1.0000000 D 1 1.00 - 1.9650000 1.0000000 + 4.0140000 1.0000000 +D 1 1.00 + 1.0960000 1.0000000 +F 1 1.00 + 2.5440000 1.0000000 diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index dbabd5c..530623a 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -19,6 +19,7 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2 integer :: a,b,c,d integer :: p integer :: ab,cd,ij,kl + double precision,external :: Kronecker_delta ! Output variables @@ -29,22 +30,25 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2 rho2(:,:,:) = 0d0 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=c+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) + rho1(p,i,ab) = rho1(p,i,ab) & + + (ERI(p,i,c,d) - 0.5d0*ERI(p,i,d,c))*X1(cd,ab)!/sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) enddo enddo kl = 0 do k=nC+1,nO - do l=k+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) + rho1(p,i,ab) = rho1(p,i,ab) & + + (ERI(p,i,k,l) - 0.5d0*ERI(p,i,l,k))*Y1(kl,ab)!/sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) enddo enddo @@ -56,17 +60,19 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2 cd = 0 do c=nO+1,nBas-nR - do d=c+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) + rho2(p,a,ij) = rho2(p,a,ij) & + + (ERI(p,a,c,d) - 0.5d0*ERI(p,a,d,c))*X2(cd,ij)!/sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(c,d))) enddo enddo kl = 0 do k=nC+1,nO - do l=k+1,nO + do l=nC+1,k kl = kl + 1 - rho1(p,a,ij) = rho1(p,a,ij) + (ERI(p,a,k,l) - ERI(p,a,l,k))*Y2(kl,ij) + rho2(p,a,ij) = rho2(p,a,ij) & + + (ERI(p,a,k,l) - 0.5d0*ERI(p,a,l,k))*Y2(kl,ij)!/sqrt((1d0 + Kronecker_delta(p,a))*(1d0 + Kronecker_delta(k,l))) enddo enddo @@ -74,4 +80,5 @@ subroutine excitation_density_Tmatrix(nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2 enddo enddo + end subroutine excitation_density_Tmatrix diff --git a/src/QuAcK/linear_response_B_pp.f90 b/src/QuAcK/linear_response_B_pp.f90 index 0dc6cbd..696b490 100644 --- a/src/QuAcK/linear_response_B_pp.f90 +++ b/src/QuAcK/linear_response_B_pp.f90 @@ -14,7 +14,6 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) ! Local variables double precision,external :: Kronecker_delta - integer :: a,b,i,j,ab,ij ! Output variables diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index 6bef720..6671e15 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -70,7 +70,6 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1 Z(:,:) = M(:,:) call diagonalize_general_matrix(nOO+nVV,M(:,:),Omega(:),Z(:,:)) -! call diagonalize_matrix(nOO+nVV,Z(:,:),Omega(:)) ! write(*,*) 'pp-RPA excitation energies' ! call matout(nOO+nVV,1,Omega(:)) diff --git a/src/QuAcK/print_G0T0.f90 b/src/QuAcK/print_G0T0.f90 index 53272b7..4cea81e 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,SigmaC,Z,eGW,EcRPA) +subroutine print_G0T0(nBas,nO,e,ENuc,EHF,SigT,Z,eGW,EcRPA) ! Print one-electron energies and other stuff for G0T0 @@ -9,7 +9,7 @@ subroutine print_G0T0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) double precision,intent(in) :: ENuc double precision,intent(in) :: EHF double precision,intent(in) :: EcRPA - double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) + double precision,intent(in) :: e(nBas),SigT(nBas),Z(nBas),eGW(nBas) integer :: x,HOMO,LUMO double precision :: Gap @@ -26,12 +26,12 @@ subroutine print_G0T0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) write(*,*)' One-shot G0T0 calculation (T-matrix self-energy) ' write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & - '|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|' + '|','#','|','e_HF (eV)','|','Sigma_T (eV)','|','Z','|','e_QP (eV)','|' write(*,*)'-------------------------------------------------------------------------------' do x=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' + '|',x,'|',e(x)*HaToeV,'|',SigT(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' enddo write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90 index d95f8c5..71466cc 100644 --- a/src/QuAcK/renormalization_factor_Tmatrix.f90 +++ b/src/QuAcK/renormalization_factor_Tmatrix.f90 @@ -34,10 +34,10 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1, do i=nC+1,nO cd = 0 do c=nO+1,nBas-nR - do d=c+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/(eps**2 + eta**2) + Z(p) = Z(p) - 2d0*rho1(p,i,cd)**2/eps**2 enddo enddo enddo @@ -49,10 +49,10 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1, do a=nO+1,nBas-nR kl = 0 do k=nC+1,nO - do l=k+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/(eps**2 + eta**2))**2 + Z(p) = Z(p) - 2d0*rho2(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 0c9e8db..73eb8dc 100644 --- a/src/QuAcK/self_energy_Tmatrix_diag.f90 +++ b/src/QuAcK/self_energy_Tmatrix_diag.f90 @@ -32,7 +32,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O ! Initialize - SigT = 0d0 + SigT(:) = 0d0 ! Occupied part of the T-matrix self-energy @@ -40,10 +40,10 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O do i=nC+1,nO cd = 0 do c=nO+1,nBas-nR - do d=c+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**2 + eta**2) + SigT(p) = SigT(p) + 2d0*rho1(p,i,cd)**2/eps enddo enddo enddo @@ -55,10 +55,10 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O do a=nO+1,nBas-nR kl = 0 do k=nC+1,nO - do l=k+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**2 + eta**2) + SigT(p) = SigT(p) + 2d0*rho2(p,a,kl)**2/eps enddo enddo enddo