diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 81c624a..779d849 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -2,4 +2,4 @@ 2 1 1 0 0 # Znuc x y z H 0. 0. 0. - H 0. 0. 1.4 + H 0. 0. 1.399 diff --git a/input/basis b/input/basis index b2b2293..6f3d2a9 100644 --- a/input/basis +++ b/input/basis @@ -1,30 +1,39 @@ -1 6 +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.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.3300000 1.0000000 D 1 - 1 0.1239000 1.0000000 + 1 4.0140000 1.0000000 +D 1 + 1 1.0960000 1.0000000 +F 1 + 1 2.5440000 1.0000000 + diff --git a/input/methods b/input/methods index ba17592..478a02f 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF MOM - F T F + T F F # MP2 MP3 MP2-F12 F F F # CCD CCSD CCSD(T) @@ -13,7 +13,7 @@ # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0 evGW qsGW - T F F + F F T # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/molecule b/input/molecule index 058d6dd..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 2 1 0 0 + 1 5 5 0 0 # Znuc x y z - Li 0.0 0.0 0.0 + Ne 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index c9a5a65..1c70680 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,3 @@ 1 - Li 0.0000000000 0.0000000000 0.0000000000 + Ne 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/options b/input/options index 337b5e1..6af6fdb 100644 --- a/input/options +++ b/input/options @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T F F F + T F T T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/ACFDT.f90 b/src/QuAcK/ACFDT.f90 index 3e5b6a5..cd8843b 100644 --- a/src/QuAcK/ACFDT.f90 +++ b/src/QuAcK/ACFDT.f90 @@ -1,5 +1,4 @@ -subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eW,e,Omega,XpY,XmY,rho,EcAC) +subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem @@ -15,8 +14,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: BSE - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS @@ -24,11 +23,6 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision :: Omega(nS,nspin) - double precision :: XpY(nS,nS,nspin) - double precision :: XmY(nS,nS,nspin) - double precision :: rho(nBas,nBas,nS,nspin) - ! Local variables integer :: ispin @@ -37,6 +31,16 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl double precision :: lambda double precision,allocatable :: Ec(:,:) + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + + double precision,allocatable :: Omega(:,:) + double precision,allocatable :: XpY(:,:,:) + double precision,allocatable :: XmY(:,:,:) + ! Output variables double precision,intent(out) :: EcAC(nspin) @@ -44,6 +48,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl ! Memory allocation allocate(Ec(nAC,nspin)) + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) + allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) ! Antisymmetrized kernel version @@ -58,12 +64,20 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl EcAC(:) = 0d0 Ec(:,:) = 0d0 +! Compute (singlet) RPA screening + + isp_W = 1 + EcRPA = 0d0 + + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + ! Singlet manifold - if(singlet_manifold) then + if(singlet) then ispin = 1 - isp_W = 1 write(*,*) '--------------' write(*,*) 'Singlet states' @@ -80,17 +94,17 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl if(doXBS) then - call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & - rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & + rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -98,6 +112,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) + if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) write(*,*) '-----------------------------------------------------------------------------------' @@ -107,7 +123,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl ! Triplet manifold - if(triplet_manifold) then + if(triplet) then ispin = 2 isp_W = 1 @@ -127,17 +143,16 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl if(doXBS) then - call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & - rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & + rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -145,6 +160,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) + if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) write(*,*) '-----------------------------------------------------------------------------------' diff --git a/src/QuAcK/BSE2.f90 b/src/QuAcK/BSE2.f90 index cff9aef..64a9e8c 100644 --- a/src/QuAcK/BSE2.f90 +++ b/src/QuAcK/BSE2.f90 @@ -54,7 +54,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & ! Compute BSE2 excitation energies call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), & - rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) ! Compute dynamic correction for BSE via perturbation theory @@ -88,7 +88,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & ! Compute BSE2 excitation energies call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), & - rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) ! Compute dynamic correction for BSE via perturbation theory diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index d01f978..d0763ee 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -1,5 +1,4 @@ -subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE) +subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -13,8 +12,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS @@ -22,51 +21,55 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision :: OmRPA(nS,nspin) - double precision :: XpY_RPA(nS,nS,nspin) - double precision :: XmY_RPA(nS,nS,nspin) - double precision :: rho_RPA(nBas,nBas,nS,nspin) - ! Local variables integer :: ispin integer :: isp_W + + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: OmBSE(:,:) double precision,allocatable :: XpY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:) ! Output variables - double precision,intent(out) :: EcRPA(nspin) double precision,intent(out) :: EcBSE(nspin) ! Memory allocation + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) +!--------------------------------- +! Compute (singlet) RPA screening +!--------------------------------- + + isp_W = 1 + EcRPA = 0d0 + + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + !------------------- ! Singlet manifold !------------------- - if(singlet_manifold) then + if(singlet) then ispin = 1 - isp_W = 1 EcBSE(ispin) = 0d0 - ! Compute (singlet) RPA screening - - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & - rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin)) - ! Compute BSE excitation energies - OmBSE(:,ispin) = OmRPA(:,ispin) - - call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) + call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) !------------------------------------------------- ! Compute the dynamical screening at the BSE level @@ -78,12 +81,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) end if end if @@ -94,25 +97,16 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man ! Triplet manifold !------------------- - if(triplet_manifold) then + if(triplet) then ispin = 2 - isp_W = 1 EcBSE(ispin) = 0d0 - ! Compute (singlet) RPA screening - - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & - rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin)) - ! Compute BSE excitation energies - OmBSE(:,ispin) = OmRPA(:,ispin) - - call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) + call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) !------------------------------------------------- ! Compute the dynamical screening at the BSE level @@ -124,12 +118,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) end if end if diff --git a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 index 25e7dea..2d82750 100644 --- a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 @@ -1,5 +1,4 @@ -subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, & - Ap,Am,Bp,Bm) +subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rhO_RPA,OmBSE,Ap,Am,Bp,Bm) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -13,8 +12,8 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -60,8 +59,8 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O do kc=1,maxS - chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) - chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) + chi_A = chi_A + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) + chi_B = chi_B + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) enddo @@ -80,28 +79,28 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O do kc=1,maxS eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) + chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) + chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) + chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) + chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) + chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) + chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) + chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) + chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) enddo diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 index 61dfa01..634323f 100644 --- a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,A_dyn) +subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -12,8 +12,8 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -48,7 +48,7 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om chi = 0d0 do kc=1,maxS - chi = chi + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) enddo @@ -58,10 +58,10 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om do kc=1,maxS eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2) + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2) eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2) + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2) enddo diff --git a/src/QuAcK/Bethe_Salpeter_B_matrix.f90 b/src/QuAcK/Bethe_Salpeter_B_matrix.f90 index c323753..7e000ee 100644 --- a/src/QuAcK/Bethe_Salpeter_B_matrix.f90 +++ b/src/QuAcK/Bethe_Salpeter_B_matrix.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr) +subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,B_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response @@ -11,8 +11,8 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho, double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) ! Local variables @@ -35,8 +35,8 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho, chi = 0d0 do kc=1,nS - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps + eps = OmRPA(kc)**2 + eta**2 + chi = chi + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/eps enddo B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,j,b,a) + 4d0*lambda*chi diff --git a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 index 14cd120..c7550da 100644 --- a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, & +subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE, & ZAp,ZAm,ZBp,ZBm) ! Compute the dynamic part of the renormalization for the Bethe-Salpeter equation matrices @@ -13,8 +13,8 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -63,28 +63,28 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, do kc=1,maxS eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 + chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 + chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 + chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 + chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 + chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 + chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 + chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 + chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 enddo diff --git a/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 index eb71cb4..e15677d 100644 --- a/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,ZA_dyn) +subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZA_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -12,8 +12,8 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -49,10 +49,10 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O do kc=1,maxS eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 enddo diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index d01e611..ba0d2dd 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) +subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -18,10 +18,10 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW, double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -84,11 +84,11 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW, ! Resonant part of the BSE correction for dynamical TDA - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,Ap_dyn) + call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn) ! Renormalization factor of the resonant parts for dynamical TDA - call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,ZAp_dyn) + call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),ZAp_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X)) @@ -97,12 +97,12 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW, ! Resonant and anti-resonant part of the BSE correction - call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, & + call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), & Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) ! Renormalization factor of the resonant and anti-resonant parts - call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, & + call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), & ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 index 9435c09..72b5e88 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) +subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,OmBSE,XpY,XmY) ! Compute self-consistently the dynamical effects via perturbation theory for BSE @@ -18,10 +18,10 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -97,8 +97,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, ! Resonant part of the BSE correction - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), & - Ap_dyn(:,:)) + call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn) OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) @@ -106,8 +105,8 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, ! Anti-resonant part of the BSE correction - call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), & - Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) + call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia), & + Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index b6ffcbc..a01c21b 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -1,5 +1,5 @@ -subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linearize,eta, & +subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & + dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,PHF,cHF,eHF,eGW) ! Perform G0W0 calculation @@ -21,8 +21,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta @@ -39,17 +39,16 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, logical :: print_W = .true. integer :: ispin - double precision :: EcRPA(nspin) + double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - double precision,allocatable :: rho(:,:,:,:) - double precision,allocatable :: rhox(:,:,:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) double precision,allocatable :: eGWlin(:) @@ -67,13 +66,13 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, ! Initialization - EcRPA(:) = 0d0 + EcRPA = 0d0 ! SOSEX correction if(SOSEX) then - write(*,*) 'SOSEX correction activated!' - write(*,*) + write(*,*) 'SOSEX correction activated but BUG!' + stop end if ! COHSEX approximation @@ -103,34 +102,43 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, ! Memory allocation - allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), & - rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin),eGWlin(nBas)) + allocate(SigC(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGWlin(nBas)) -! Compute linear response +!-------------------! +! Compute screening ! +!-------------------! - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, & + eHF,ERI,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - if(print_W) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) + if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) -! Compute correlation part of the self-energy +!--------------------------! +! Compute spectral weights ! +!--------------------------! - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin)) +!------------------------! +! Compute GW self-energy ! +!------------------------! - call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) + call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) -! Compute renormalization factor +!--------------------------------! +! Compute renormalization factor ! +!--------------------------------! - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) -! Solve the quasi-particle equation +!-----------------------------------! +! Solve the quasi-particle equation ! +!-----------------------------------! eGWlin(:) = eHF(:) + Z(:)*SigC(:) + ! Linearized or graphical solution? + if(linearize) then write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' @@ -138,51 +146,48 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, eGW(:) = eGWlin(:) -! Find all the roots of the QP equation if necessary - -! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin) - else - - ! Find graphical solution of the QP equation - call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + + call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,OmRPA,rho_RPA,eGWlin,eGW) + + ! Find all the roots of the QP equation if necessary + + ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin) end if -! Dump results - - call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM) - ! Compute the RPA correlation energy - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - 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(*,*) +!--------------! +! Dump results ! +!--------------! + + call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) + +! Deallocate memory + + deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) ! Plot stuff -! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin)) +! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA) ! Perform BSE calculation if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcBSE) if(exchange_kernel) then - EcRPA(1) = 0.5d0*EcRPA(1) - EcRPA(2) = 1.5d0*EcRPA(1) + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) end if @@ -211,15 +216,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC) - - if(exchange_kernel) then - - EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(1) - - end if + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/RPAx.f90 b/src/QuAcK/RPAx.f90 index 06c298a..ca9f709 100644 --- a/src/QuAcK/RPAx.f90 +++ b/src/QuAcK/RPAx.f90 @@ -1,4 +1,4 @@ -subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) +subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Perform random phase approximation calculation with exchange (aka TDHF) @@ -21,7 +21,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: e(nBas) + double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables @@ -58,12 +58,11 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, ispin = 1 - call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - endif ! Triplet manifold @@ -72,9 +71,9 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, ispin = 2 - call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif @@ -105,14 +104,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, write(*,*) call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) - - if(exchange_kernel) then - - EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(2) - - end if + nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/dRPA.f90 b/src/QuAcK/dRPA.f90 index b807497..d17ed51 100644 --- a/src/QuAcK/dRPA.f90 +++ b/src/QuAcK/dRPA.f90 @@ -1,5 +1,5 @@ -subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) +subroutine dRPA(doACFDT,exchange_kernel,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Perform a direct random phase approximation calculation @@ -11,8 +11,8 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -22,7 +22,7 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: e(nBas) + double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables @@ -40,7 +40,7 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & write(*,*) write(*,*)'***********************************************' - write(*,*)'| random-phase approximation calculation |' + write(*,*)'| Random-phase approximation calculation |' write(*,*)'***********************************************' write(*,*) @@ -55,32 +55,34 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & ! Singlet manifold - if(singlet_manifold) then + if(singlet) then ispin = 1 - call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif ! Triplet manifold - if(triplet_manifold) then + if(triplet) then ispin = 2 - call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif if(exchange_kernel) then EcRPA(1) = 0.5d0*EcRPA(1) - EcRPA(2) = 1.5d0*EcRPA(1) + EcRPA(2) = 1.5d0*EcRPA(2) end if @@ -103,13 +105,13 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & write(*,*) '------------------------------------------------------' write(*,*) - call ACFDT(exchange_kernel,.false.,.true.,.false.,.false.,.false.,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) + call ACFDT(exchange_kernel,.false.,.true.,.false.,.false.,.false.,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) if(exchange_kernel) then EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(2) end if diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 625c1fa..79f2dc2 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,5 +1,6 @@ -subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0, & - singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) +subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & + G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H, & + ERI,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -26,8 +27,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE logical,intent(in) :: evDyn logical,intent(in) :: G0W logical,intent(in) :: GW0 - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: eHF(nBas) @@ -46,7 +47,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE integer :: n_diis double precision :: rcond double precision :: Conv - double precision :: EcRPA(nspin) + double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM @@ -57,11 +58,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,allocatable :: eOld(:) double precision,allocatable :: Z(:) double precision,allocatable :: SigC(:) - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - double precision,allocatable :: rho(:,:,:,:) - double precision,allocatable :: rhox(:,:,:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) ! Hello world @@ -73,23 +73,45 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! SOSEX correction - if(SOSEX) write(*,*) 'SOSEX correction activated!' - write(*,*) + if(SOSEX) then + write(*,*) 'SOSEX correction activated but BUG!' + stop + end if ! COHSEX approximation - if(COHSEX) write(*,*) 'COHSEX approximation activated!' - write(*,*) + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if ! TDA for W - if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if ! TDA - if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! GW0 + + if(GW0) then + write(*,*) 'GW0 scheme activated!' + write(*,*) + end if + +! G0W + + if(G0W) then + write(*,*) 'G0W scheme activated!' + write(*,*) + end if ! Linear mixing @@ -98,10 +120,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Memory allocation - allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), & - XpY(nS,nS,nspin),XmY(nS,nS,nspin), & - rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & - error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & + rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) ! Initialization @@ -121,36 +141,30 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE do while(Conv > thresh .and. nSCF <= maxSCF) - ! Compute linear response + ! Compute screening if(.not. GW0 .or. nSCF == 0) then - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) endif +! Compute spectral weights + + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + ! Compute correlation part of the self-energy - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - - if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin)) - - ! Correlation self-energy - if(G0W) then - call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) + call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) else - call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) + call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) endif @@ -164,8 +178,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) + call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) ! Linear mixing or DIIS extrapolation @@ -215,23 +228,22 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE endif -! Dump the RPA correlation energy +! Deallocate memory - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) + deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis) ! Perform BSE calculation if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE) + + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) + + end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -258,8 +270,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcAC) + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/linear_response.f90 b/src/QuAcK/linear_response.f90 index c3ea157..e443d52 100644 --- a/src/QuAcK/linear_response.f90 +++ b/src/QuAcK/linear_response.f90 @@ -1,4 +1,4 @@ -subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY,XmY) +subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Omega_RPA,rho_RPA,EcRPA,Omega,XpY,XmY) ! Compute linear response @@ -12,8 +12,10 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS double precision,intent(in) :: lambda double precision,intent(in) :: e(nBas) - double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + + double precision,intent(in) :: Omega_RPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) ! Local variables @@ -42,7 +44,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) - if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A) + if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,A) ! Tamm-Dancoff approximation @@ -51,7 +53,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) - if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B) + if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,B) end if diff --git a/src/QuAcK/plot_GW.f90 b/src/QuAcK/plot_GW.f90 index 11c7f0c..fedc40c 100644 --- a/src/QuAcK/plot_GW.f90 +++ b/src/QuAcK/plot_GW.f90 @@ -1,4 +1,4 @@ -subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho,rhox) +subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho) ! Dump several GW quantities for external plotting @@ -8,7 +8,7 @@ subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho,rhox) ! Input variables integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: eHF(nBas),eGW(nBas),Omega(nS),rho(nBas,nBas,nS),rhox(nBas,nBas,nS) + double precision,intent(in) :: eHF(nBas),eGW(nBas),Omega(nS),rho(nBas,nBas,nS) ! Local variables diff --git a/src/QuAcK/print_G0W0.f90 b/src/QuAcK/print_G0W0.f90 index d1ee3f4..782643b 100644 --- a/src/QuAcK/print_G0W0.f90 +++ b/src/QuAcK/print_G0W0.f90 @@ -40,8 +40,8 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'RPA@HF total energy =',ENuc + EHF + EcRPA - write(*,'(2X,A30,F15.6)') 'RPA@HF correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/print_evGW.f90 b/src/QuAcK/print_evGW.f90 index 0bf95b5..8db2956 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) +subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA,EcGM) ! Print one-electron energies and other stuff for evGW @@ -8,9 +8,15 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW) integer,intent(in) :: nBas,nO,nSCF double precision,intent(in) :: ENuc double precision,intent(in) :: EHF - double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) + double precision,intent(in) :: Conv + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: SigC(nBas) + double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: EcRPA + double precision,intent(in) :: EcGM - integer :: x,HOMO,LUMO + integer :: p,HOMO,LUMO double precision :: Gap ! HOMO and LUMO @@ -32,9 +38,9 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW) '|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do x=1,nBas + do p=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,'|' + '|',p,'|',e(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' enddo write(*,*)'-------------------------------------------------------------------------------' @@ -45,11 +51,11 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW) write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',eGW(LUMO)*HaToeV write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA -! write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA -! write(*,'(2X,A30,F15.6)') 'GM@evGW total energy =',ENuc + EHF + EcGM -! write(*,'(2X,A30,F15.6)') 'GM@evGW correlation energy =',EcGM -! write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'GM@evGW total energy =',ENuc + EHF + EcGM + write(*,'(2X,A30,F15.6)') 'GM@evGW correlation energy =',EcGM + write(*,*)'-------------------------------------------------------------------------------' write(*,*) end subroutine print_evGW diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90 index e5ce17b..7aec451 100644 --- a/src/QuAcK/print_excitation.f90 +++ b/src/QuAcK/print_excitation.f90 @@ -26,7 +26,7 @@ subroutine print_excitation(method,ispin,nS,Omega) write(*,*) write(*,*)'-------------------------------------------------------------' - write(*,'(1X,A1,1X,A14,A14,A14,A9,A13)')'|',method,' calculation: ',spin_manifold,' manifold',' |' + write(*,'(1X,A14,A14,A14,A9)') method,' calculation: ',spin_manifold,' manifold' write(*,*)'-------------------------------------------------------------' write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' diff --git a/src/QuAcK/print_qsGW.f90 b/src/QuAcK/print_qsGW.f90 index b61bdd3..9f3cf62 100644 --- a/src/QuAcK/print_qsGW.f90 +++ b/src/QuAcK/print_qsGW.f90 @@ -70,8 +70,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig write(*,'(2X,A30,F15.6)') 'qsGW GM total energy =',EqsGW + ENuc + EcGM write(*,'(2X,A30,F15.6)') 'qsGW exchange energy =',Ex write(*,'(2X,A30,F15.6)') 'qsGW correlation energy =',Ec -! write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA -! write(*,'(2X,A30,F15.6)') 'GM@qsGW correlation energy =',EcGM + write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'GM@qsGW correlation energy =',EcGM write(*,*)'-------------------------------------------' write(*,*) @@ -95,7 +95,6 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig write(*,'(A32,1X,F16.10)') ' Electronic energy ',EqsGW write(*,'(A32,1X,F16.10)') ' Nuclear repulsion ',ENuc write(*,'(A32,1X,F16.10)') ' qsGW energy ',ENuc + EqsGW -! write(*,'(A32,1X,F16.10)') ' RPA corr. energy ',EcRPA write(*,'(A50)') '---------------------------------------' write(*,*) diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index dad60a5..b1fbfe2 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -1,6 +1,6 @@ -subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & - COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF) +subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & + G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V, & + Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -25,8 +25,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, logical,intent(in) :: evDyn logical,intent(in) :: G0W logical,intent(in) :: GW0 - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: ENuc @@ -49,7 +49,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, integer :: ispin integer :: n_diis double precision :: EqsGW - double precision :: EcRPA(nspin) + double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM @@ -58,11 +58,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, double precision,external :: trace_matrix double precision,allocatable :: error_diis(:,:) double precision,allocatable :: F_diis(:,:) - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - double precision,allocatable :: rho(:,:,:,:) - double precision,allocatable :: rhox(:,:,:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) double precision,allocatable :: c(:,:) double precision,allocatable :: cp(:,:) double precision,allocatable :: eGW(:) @@ -96,29 +95,37 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! SOSEX correction - if(SOSEX) write(*,*) 'SOSEX correction activated!' - write(*,*) + if(SOSEX) then + write(*,*) 'SOSEX correction activated but BUG!' + stop + end if ! COHSEX approximation - if(COHSEX) write(*,*) 'COHSEX approximation activated!' - write(*,*) + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if ! TDA for W - if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if ! TDA - if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if ! Memory allocation allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & - Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & + OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -156,29 +163,23 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, if(.not. GW0 .or. nSCF == 0) then call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) endif ! Compute correlation part of the self-energy - call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) - - if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY_RPA,rho_RPA) if(G0W) then - call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) + call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) else - call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) + call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) endif @@ -221,7 +222,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Print results ! call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW) + call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA,EcGM,EqsGW) ! Increment @@ -258,23 +259,22 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, endif -! Dump RPA correlation energy +! Deallocate memory - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW total energy =',ENuc + EqsGW + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) + deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis) ! Perform BSE calculation if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,EcBSE) + + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) + + end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -301,8 +301,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcAC) + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/renormalization_factor.f90 b/src/QuAcK/renormalization_factor.f90 index 15e3bf2..10f2849 100644 --- a/src/QuAcK/renormalization_factor.f90 +++ b/src/QuAcK/renormalization_factor.f90 @@ -1,4 +1,4 @@ -subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) +subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) ! Compute renormalization factor for GW @@ -8,13 +8,11 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r ! Input variables logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) - double precision,intent(in) :: rhox(nBas,nBas,nS) ! Local variables @@ -36,40 +34,6 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r Z(:) = 1d0 return -! SOSEX correction - - elseif(SOSEX) then - - ! Occupied part of the correlation self-energy - - do x=nC+1,nBas-nR - do i=nC+1,nO - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = e(x) - e(i) + Omega(jb) - Z(x) = Z(x) - (rho(x,i,jb)/eps)*(rhox(x,i,jb)/eps) - end do - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do x=nC+1,nBas-nR - do a=nO+1,nBas-nR - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = e(x) - e(a) - Omega(jb) - Z(x) = Z(x) - (rho(x,a,jb)/eps)*(rhox(x,a,jb)/eps) - end do - end do - end do - end do - else ! Occupied part of the correlation self-energy diff --git a/src/QuAcK/self_energy_correlation.f90 b/src/QuAcK/self_energy_correlation.f90 index 5baccea..e47a9b4 100644 --- a/src/QuAcK/self_energy_correlation.f90 +++ b/src/QuAcK/self_energy_correlation.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) +subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) ! Compute correlation part of the self-energy @@ -8,13 +8,11 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega, ! Input variables logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) - double precision,intent(in) :: rhox(nBas,nBas,nS) ! Local variables @@ -30,7 +28,9 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega, SigC = 0d0 -! COHSEX static approximation +!-----------------------------! +! COHSEX static approximation ! +!-----------------------------! if(COHSEX) then @@ -65,6 +65,10 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega, else +!----------------! +! GW self-energy ! +!----------------! + ! Occupied part of the correlation self-energy do x=nC+1,nBas-nR @@ -91,36 +95,6 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega, enddo enddo - if(SOSEX) then - - ! SOSEX: occupied part of the correlation self-energy - - do x=nC+1,nBas-nR - do y=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - eps = e(x) - e(i) + Omega(jb) - SigC(x,y) = SigC(x,y) - rho(x,i,jb)*rhox(y,i,jb)*eps/(eps**2 + eta**2) - enddo - enddo - enddo - enddo - - ! SOSEX: virtual part of the correlation self-energy - - do x=nC+1,nBas-nR - do y=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(x) - e(a) - Omega(jb) - SigC(x,y) = SigC(x,y) - rho(x,a,jb)*rhox(y,a,jb)*eps/(eps**2 + eta**2) - enddo - enddo - enddo - enddo - - endif - endif end subroutine self_energy_correlation diff --git a/src/QuAcK/self_energy_correlation_diag.f90 b/src/QuAcK/self_energy_correlation_diag.f90 index 441806d..59aab3d 100644 --- a/src/QuAcK/self_energy_correlation_diag.f90 +++ b/src/QuAcK/self_energy_correlation_diag.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) +subroutine self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) ! Compute diagonal of the correlation part of the self-energy @@ -8,7 +8,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O ! Input variables logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -19,7 +18,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) - double precision,intent(in) :: rhox(nBas,nBas,nS) ! Local variables @@ -69,45 +67,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O EcGM = EcGM - SigC(i) end do -!----------------------------- -! SOSEX self-energy *BUG* -!----------------------------- - - elseif(SOSEX) then - - ! SOSEX: occupied part of the correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - eps = e(p) - e(i) + Omega(jb) - SigC(p) = SigC(p) - rho(p,i,jb)*rhox(p,i,jb)*eps/(eps**2 + eta**2) - end do - end do - end do - - ! SOSEX: virtual part of the correlation self-energy - - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(p) - e(a) - Omega(jb) - SigC(p) = SigC(p) - rho(p,a,jb)*rhox(p,a,jb)*eps/(eps**2 + eta**2) - end do - end do - end do - - ! GM correlation energy - - do i=nC+1,nO - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)*eps/(eps**2 + eta**2) - end do - end do - end do - !----------------------------- ! GW self-energy !-----------------------------