From 2f835304f70d741559118e0397d7fe25f2f19111 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 13 Apr 2020 14:19:14 +0200 Subject: [PATCH] evGT --- input/basis | 68 ++++----- input/methods | 2 +- input/molecule | 5 +- input/molecule.xyz | 5 +- input/weight | 68 ++++----- src/QuAcK/G0T0.f90 | 22 +++ src/QuAcK/QuAcK.f90 | 19 ++- src/QuAcK/evGT.f90 | 318 +++++++++++++++++++++++++++++++++++++++ src/QuAcK/evGW.f90 | 4 +- src/QuAcK/print_G0T0.f90 | 10 +- src/QuAcK/print_evGT.f90 | 54 +++++++ src/QuAcK/print_evGW.f90 | 4 +- 12 files changed, 482 insertions(+), 97 deletions(-) create mode 100644 src/QuAcK/evGT.f90 create mode 100644 src/QuAcK/print_evGT.f90 diff --git a/input/basis b/input/basis index 20bc6a7..6f3d2a9 100644 --- a/input/basis +++ b/input/basis @@ -1,49 +1,39 @@ -1 9 +1 10 S 8 - 1 1469.0000000 0.0007660 - 2 220.5000000 0.0058920 - 3 50.2600000 0.0296710 - 4 14.2400000 0.1091800 - 5 4.5810000 0.2827890 - 6 1.5800000 0.4531230 - 7 0.5640000 0.2747740 - 8 0.0734500 0.0097510 + 1 24350.0000000 0.0005020 + 2 3650.0000000 0.0038810 + 3 829.6000000 0.0199970 + 4 234.0000000 0.0784180 + 5 75.6100000 0.2296760 + 6 26.7300000 0.4327220 + 7 9.9270000 0.3506420 + 8 1.1020000 -0.0076450 S 8 - 1 1469.0000000 -0.0001200 - 2 220.5000000 -0.0009230 - 3 50.2600000 -0.0046890 - 4 14.2400000 -0.0176820 - 5 4.5810000 -0.0489020 - 6 1.5800000 -0.0960090 - 7 0.5640000 -0.1363800 - 8 0.0734500 0.5751020 + 1 24350.0000000 -0.0001180 + 2 3650.0000000 -0.0009150 + 3 829.6000000 -0.0047370 + 4 234.0000000 -0.0192330 + 5 75.6100000 -0.0603690 + 6 26.7300000 -0.1425080 + 7 9.9270000 -0.1777100 + 8 1.1020000 0.6058360 S 1 - 1 0.0280500 1.0000000 + 1 2.8360000 1.0000000 S 1 - 1 0.0086400 1.0000000 + 1 0.3782000 1.0000000 P 3 - 1 1.5340000 0.0227840 - 2 0.2749000 0.1391070 - 3 0.0736200 0.5003750 + 1 54.7000000 0.0171510 + 2 12.4300000 0.1076560 + 3 3.6790000 0.3216810 P 1 - 1 0.0240300 1.0000000 + 1 1.1430000 1.0000000 P 1 - 1 0.0057900 1.0000000 + 1 0.3300000 1.0000000 D 1 - 1 0.1239000 1.0000000 + 1 4.0140000 1.0000000 D 1 - 1 0.0725000 1.0000000 -2 5 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 -S 1 - 1 0.1220000 1.0000000 -S 1 - 1 0.0297400 1.0000000 -P 1 - 1 0.7270000 1.0000000 -P 1 - 1 0.1410000 1.0000000 + 1 1.0960000 1.0000000 +F 1 + 1 2.5440000 1.0000000 + diff --git a/input/methods b/input/methods index c1d2dda..f5a20e0 100644 --- a/input/methods +++ b/input/methods @@ -13,6 +13,6 @@ # G0W0 evGW qsGW F F F # G0T0 evGT qsGT - T F F + T T F # MCMP2 F diff --git a/input/molecule b/input/molecule index e76247c..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 2 2 0 0 + 1 5 5 0 0 # Znuc x y z - Li 0. 0. 0. - H 0. 0. 3.099 + Ne 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index fb94244..1c70680 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - Li 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 1.6399202947 + Ne 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/weight b/input/weight index 20bc6a7..6f3d2a9 100644 --- a/input/weight +++ b/input/weight @@ -1,49 +1,39 @@ -1 9 +1 10 S 8 - 1 1469.0000000 0.0007660 - 2 220.5000000 0.0058920 - 3 50.2600000 0.0296710 - 4 14.2400000 0.1091800 - 5 4.5810000 0.2827890 - 6 1.5800000 0.4531230 - 7 0.5640000 0.2747740 - 8 0.0734500 0.0097510 + 1 24350.0000000 0.0005020 + 2 3650.0000000 0.0038810 + 3 829.6000000 0.0199970 + 4 234.0000000 0.0784180 + 5 75.6100000 0.2296760 + 6 26.7300000 0.4327220 + 7 9.9270000 0.3506420 + 8 1.1020000 -0.0076450 S 8 - 1 1469.0000000 -0.0001200 - 2 220.5000000 -0.0009230 - 3 50.2600000 -0.0046890 - 4 14.2400000 -0.0176820 - 5 4.5810000 -0.0489020 - 6 1.5800000 -0.0960090 - 7 0.5640000 -0.1363800 - 8 0.0734500 0.5751020 + 1 24350.0000000 -0.0001180 + 2 3650.0000000 -0.0009150 + 3 829.6000000 -0.0047370 + 4 234.0000000 -0.0192330 + 5 75.6100000 -0.0603690 + 6 26.7300000 -0.1425080 + 7 9.9270000 -0.1777100 + 8 1.1020000 0.6058360 S 1 - 1 0.0280500 1.0000000 + 1 2.8360000 1.0000000 S 1 - 1 0.0086400 1.0000000 + 1 0.3782000 1.0000000 P 3 - 1 1.5340000 0.0227840 - 2 0.2749000 0.1391070 - 3 0.0736200 0.5003750 + 1 54.7000000 0.0171510 + 2 12.4300000 0.1076560 + 3 3.6790000 0.3216810 P 1 - 1 0.0240300 1.0000000 + 1 1.1430000 1.0000000 P 1 - 1 0.0057900 1.0000000 + 1 0.3300000 1.0000000 D 1 - 1 0.1239000 1.0000000 + 1 4.0140000 1.0000000 D 1 - 1 0.0725000 1.0000000 -2 5 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 -S 1 - 1 0.1220000 1.0000000 -S 1 - 1 0.0297400 1.0000000 -P 1 - 1 0.7270000 1.0000000 -P 1 - 1 0.1410000 1.0000000 + 1 1.0960000 1.0000000 +F 1 + 1 2.5440000 1.0000000 + diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index d3c471c..d6648d9 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -182,6 +182,28 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m call print_G0T0(nBas,nO,eHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA(:)) +! Compute the ppRPA correlation energy + + ispin = 1 + iblock = 3 + call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eG0T0(:),ERI(:,:,:,:), & + Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin)) + ispin = 2 + iblock = 4 + call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eG0T0(:),ERI(:,:,:,:), & + Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin)) + EcRPA(1) = EcRPA(1) - EcRPA(2) + EcRPA(2) = 3d0*EcRPA(2) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + ! Perform BSE calculation if(BSE) then diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index aa5b86b..ade538d 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -685,7 +685,6 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) -! call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF) call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, & singlet_manifold,triplet_manifold,linGW,eta, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF) @@ -698,6 +697,24 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Perform evGT calculatiom +!------------------------------------------------------------------------ + + if(doevGT) then + + call cpu_time(start_evGT) + call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_manifold, & + eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0) + call cpu_time(end_evGT) + + t_evGT = end_evGT - start_evGT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_evGT,' seconds' + write(*,*) + + end if + + !------------------------------------------------------------------------ ! Information for Monte Carlo calculations !------------------------------------------------------------------------ diff --git a/src/QuAcK/evGT.f90 b/src/QuAcK/evGT.f90 new file mode 100644 index 0000000..acaf92b --- /dev/null +++ b/src/QuAcK/evGT.f90 @@ -0,0 +1,318 @@ +subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_manifold, & + eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0) + +! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: BSE + logical,intent(in) :: TDA + logical,intent(in) :: singlet_manifold + logical,intent(in) :: triplet_manifold + 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 + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eG0T0(nBas) + + +! Local variables + + logical :: linear_mixing + integer :: nSCF + integer :: n_diis + double precision :: rcond + double precision :: Conv + integer :: ispin + integer :: iblock + integer :: nOOs,nOOt + integer :: nVVs,nVVt + double precision :: dERI + double precision :: xERI + double precision :: alpha + double precision :: EcRPA(nspin) + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: e_diis(:,:) + double precision,allocatable :: eGT(:) + double precision,allocatable :: eOld(:) + 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(:) + + double precision,allocatable :: Omega(:,:) + double precision,allocatable :: XpY(:,:,:) + double precision,allocatable :: XmY(:,:,:) + double precision,allocatable :: rho(:,:,:,:) + +! Output variables + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Self-consistent evGT calculation |' + write(*,*)'************************************************' + write(*,*) + +! Dimensions of the pp-RPA linear reponse matrices + + nOOs = nO*nO + nVVs = nV*nV + + nOOt = nO*(nO - 1)/2 + nVVt = nV*(nV - 1)/2 + +! Memory allocation + + allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + rho1s(nBas,nO,nVVs),rho2s(nBas,nV,nOOs), & + Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + rho1t(nBas,nO,nVVt),rho2t(nBas,nV,nOOt), & + eGT(nBas),eOld(nBas),Z(nBas),SigT(nBas), & + error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + +! Initialization + + nSCF = 0 + n_diis = 0 + Conv = 1d0 + e_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + eGT(:) = eG0T0(:) + eOld(:) = eGT(:) + Z(:) = 1d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF <= maxSCF) + + !---------------------------------------------- + ! alpha-beta block + !---------------------------------------------- + + ispin = 1 + iblock = 3 + + ! Compute linear response + + call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF(:),ERI(:,:,:,:), & + Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin)) + +! EcRPA(ispin) = 1d0*EcRPA(ispin) + +! call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:)) +! call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:)) + + !---------------------------------------------- + ! alpha-alpha block + !---------------------------------------------- + + ispin = 2 + iblock = 4 + + ! Compute linear response + + call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:),ERI(:,:,:,:), & + Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin)) + +! EcRPA(ispin) = 2d0*EcRPA(ispin) +! EcRPA(ispin) = 3d0*EcRPA(ispin) + +! call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:)) +! call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:)) + + !---------------------------------------------- + ! Compute T-matrix version of the self-energy + !---------------------------------------------- + + SigT(:) = 0d0 + Z(:) = 0d0 + + iblock = 3 + dERI = +1d0 + xERI = +0d0 + alpha = +1d0 + + call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & + X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) + + call self_energy_Tmatrix_diag(alpha,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:), & + Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:),SigT(:)) + + call renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:), & + Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:),Z(:)) + + iblock = 4 + dERI = +1d0 + xERI = -1d0 + alpha = +1d0 + + call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & + X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) + + call self_energy_Tmatrix_diag(alpha,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:), & + Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),SigT(:)) + + call renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:), & + Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),Z(:)) + + Z(:) = 1d0/(1d0 - Z(:)) + + ! Solve the quasi-particle equation + + !---------------------------------------------- + ! Solve the quasi-particle equation + !---------------------------------------------- + + eGT(:) = eHF(:) + SigT(:) + + ! Convergence criteria + + Conv = maxval(abs(eGT(:) - eOld(:))) + + !---------------------------------------------- + ! Dump results + !---------------------------------------------- + + call print_evGT(nBas,nO,nSCF,Conv,eHF(:),SigT(:),Z(:),eGT(:)) + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGT(:)-eOld(:),eGT(:)) + + ! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 + + ! Save quasiparticles energy for next cycle + + eOld(:) = eGT(:) + + ! Increment + + nSCF = nSCF + 1 + + enddo +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Compute the ppRPA correlation energy + + ispin = 1 + iblock = 3 + call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT(:),ERI(:,:,:,:), & + Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin)) + ispin = 2 + iblock = 4 + call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT(:),ERI(:,:,:,:), & + Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin)) + EcRPA(1) = EcRPA(1) - EcRPA(2) + EcRPA(2) = 3d0*EcRPA(2) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + +! Perform BSE calculation + + if(BSE) then + + allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin)) + + call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcRPA,EcBSE) + + if(exchange_kernel) then + + EcRPA(1) = 0.5d0*EcRPA(1) + EcRPA(2) = 1.5d0*EcRPA(1) + + end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGT correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGT correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGT correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGT total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + + if(doACFDT) then + + write(*,*) '------------------------------------------------------' + write(*,*) 'Adiabatic connection version of BSE correlation energy' + write(*,*) '------------------------------------------------------' + write(*,*) + + if(doXBS) then + + write(*,*) '*** scaled screening version (XBS) ***' + write(*,*) + + end if + + call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcAC) + + if(exchange_kernel) then + + EcAC(1) = 0.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(1) + + end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@BSE@evGT correlation energy (singlet) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@BSE@evGT correlation energy (triplet) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@evGT correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@evGT total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + + end if + +end subroutine evGT diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index c18144c..3db7bd2 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -78,8 +78,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE if(COHSEX) write(*,*) 'COHSEX approximation activated!' write(*,*) -! Switch off exchange for G0W0 - ! Linear mixing linear_mixing = .false. @@ -154,7 +152,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Print results ! call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM) + call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW) ! Linear mixing or DIIS extrapolation diff --git a/src/QuAcK/print_G0T0.f90 b/src/QuAcK/print_G0T0.f90 index 165fa9b..ff33862 100644 --- a/src/QuAcK/print_G0T0.f90 +++ b/src/QuAcK/print_G0T0.f90 @@ -51,11 +51,11 @@ subroutine print_G0T0(nBas,nO,e,ENuc,ERHF,SigT,Z,eGW,EcRPA) 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,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(*,'(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(*,*) end subroutine print_G0T0 diff --git a/src/QuAcK/print_evGT.f90 b/src/QuAcK/print_evGT.f90 new file mode 100644 index 0000000..9b15356 --- /dev/null +++ b/src/QuAcK/print_evGT.f90 @@ -0,0 +1,54 @@ +subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT) + +! Print one-electron energies and other stuff for evGT + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas + integer,intent(in) :: nO + integer,intent(in) :: nSCF + double precision,intent(in) :: Conv + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: SigT(nBas) + double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eGT(nBas) + + integer :: x,HOMO,LUMO + double precision :: Gap + +! HOMO and LUMO + + HOMO = nO + LUMO = HOMO + 1 + Gap = eGT(LUMO)-eGT(HOMO) + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent evG',nSCF,'T',nSCF,' calculation' + else + write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent evG',nSCF,'T',nSCF,' calculation' + endif + 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_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,'|',eHF(x)*HaToeV,'|',SigT(x)*HaToeV,'|',Z(x),'|',eGT(x)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'evGT HOMO energy (eV):',eGT(HOMO)*HaToeV + write(*,'(2X,A30,F15.6)') 'evGT LUMO energy (eV):',eGT(LUMO)*HaToeV + write(*,'(2X,A30,F15.6)') 'evGT HOMO-LUMO gap (eV):',Gap*HaToeV + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +end subroutine print_evGT diff --git a/src/QuAcK/print_evGW.f90 b/src/QuAcK/print_evGW.f90 index 7596655..0bf95b5 100644 --- a/src/QuAcK/print_evGW.f90 +++ b/src/QuAcK/print_evGW.f90 @@ -1,4 +1,4 @@ -subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) +subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW) ! Print one-electron energies and other stuff for evGW @@ -8,8 +8,6 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) integer,intent(in) :: nBas,nO,nSCF double precision,intent(in) :: ENuc double precision,intent(in) :: EHF - double precision,intent(in) :: EcRPA - double precision,intent(in) :: EcGM double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) integer :: x,HOMO,LUMO