From f345b8e2be957d52cf4060e8adff6e1880446330 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sat, 4 May 2019 17:13:50 +0200 Subject: [PATCH] BSE Ec --- input/basis | 170 ++++++++++++++++++++- input/methods | 4 +- input/molecule | 2 +- input/options | 4 +- run_sph.sh | 4 +- src/QuAcK/ADC2.f90 | 55 +++---- src/QuAcK/CCSD.f90 | 4 +- src/QuAcK/G0W0.f90 | 28 +++- src/QuAcK/evGW.f90 | 26 +++- src/QuAcK/linear_response_B_matrix.f90 | 2 +- src/QuAcK/print_G0W0.f90 | 5 +- src/QuAcK/print_RHF.f90 | 2 +- src/QuAcK/print_evGW.f90 | 5 +- src/QuAcK/self_energy_correlation.f90 | 5 +- src/QuAcK/self_energy_correlation_diag.f90 | 6 +- 15 files changed, 263 insertions(+), 59 deletions(-) diff --git a/input/basis b/input/basis index 9c24319..18db43a 100644 --- a/input/basis +++ b/input/basis @@ -1,4 +1,172 @@ -1 16 +1 100 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 +S 1 1.00 + 1.0000000 1.0000000 S 1 1.00 1.0000000 1.0000000 S 1 1.00 diff --git a/input/methods b/input/methods index fd7a88f..341e057 100644 --- a/input/methods +++ b/input/methods @@ -1,13 +1,13 @@ # RHF UHF MOM T F F # MP2 MP3 MP2-F12 - T F F + F F F # CCD CCSD CCSD(T) F F F # CIS TDHF ADC F F F # GF2 GF3 - T T + F F # G0W0 evGW qsGW T F F # MCMP2 diff --git a/input/molecule b/input/molecule index 39a84f4..321523e 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 4 4 0 0 + 1 9 9 0 0 # Znuc x y z X 0.0 0.0 0.0 diff --git a/input/options b/input/options index 3699ef2..5095ea0 100644 --- a/input/options +++ b/input/options @@ -5,10 +5,10 @@ # CC: maxSCF thresh DIIS n_diis 64 0.00001 F 1 # CIS/TDHF: singlet triplet - T F + T T # GF: maxSCF thresh DIIS n_diis renormalization 64 0.00001 T 5 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize - 64 0.00001 T 5 F F T F F F F + 64 0.00001 T 3 F F T F F F F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/run_sph.sh b/run_sph.sh index afaa813..bfe1c48 100755 --- a/run_sph.sh +++ b/run_sph.sh @@ -1,8 +1,8 @@ #! /bin/bash Lmin=0 -Lmax=2 -Mmax=3 +Lmax=0 +Mmax=9 rs=$1 if [ $# != 1 ] diff --git a/src/QuAcK/ADC2.f90 b/src/QuAcK/ADC2.f90 index 85d5469..1cb2dbe 100644 --- a/src/QuAcK/ADC2.f90 +++ b/src/QuAcK/ADC2.f90 @@ -12,7 +12,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) double precision,intent(in) :: thresh integer,intent(in) :: max_diis integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables @@ -57,8 +58,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) ! Memory allocation - allocate(db_ERI(nBas,nBas,nBas,nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eOld(nADC), & - B_ADC(nADC,nADC),X_ADC(nADC,nADC),e_ADC(nADC),G_ADC(nADC,nADC),SigInf(nADC,nADC)) + allocate(db_ERI(nBas,nBas,nBas,nBas),error_diis(nADC,max_diis),e_diis(nADC,max_diis),eOld(nADC), & + B_ADC(nADC,nADC),X_ADC(nADC,nADC),e_ADC(nADC),G_ADC(nBas,nBas),SigInf(nBas,nBas)) ! Create double-bar MO integrals @@ -99,11 +100,11 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) ! B_ADC(:,:) = 0d0 - B_ADC(nC+1:nV,nC+1:nV) = SigInf(nC+1:nV,nC+1:nV) + B_ADC(nC+1:nBas-nR,nC+1:nBas-nR) = SigInf(nC+1:nBas-nR,nC+1:nBas-nR) jADC = 0 - do p=nC+1,nV + do p=nC+1,nBas-nR jADC = jADC + 1 B_ADC(jADC,jADC) = e(p) @@ -114,7 +115,7 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) ! U matrices -- Eq. (40a) -- ! - do p=nC+1,nV + do p=nC+1,nBas-nR iADC = p - nC jADC = nH + nP @@ -122,8 +123,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) ! U^I: 2p-1h -- Eqs. (40a) and (41a) -- do i=nC+1,nO - do a=nO+1,nV-nR - do b=a,nV-nR + do a=nO+1,nBas-nR + do b=a,nBas-nR jADC = jADC + 1 B_ADC(iADC,jADC) = db_ERI(p,i,a,b) @@ -136,7 +137,7 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) do i=nC+1,nO do j=i,nO - do a=nO+1,nV-nR + do a=nO+1,nBas-nR jADC = jADC + 1 B_ADC(iADC,jADC) = db_ERI(p,a,i,j) @@ -156,8 +157,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) jADC = nH + nP do i=nC+1,nO - do a=nO+1,nV-nR - do b=a,nV-nR + do a=nO+1,nBas-nR + do b=a,nBas-nR jADC = jADC + 1 B_ADC(jADC,jADC) = e(a) + e(b) - e(i) @@ -170,7 +171,7 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) do i=nC+1,nO do j=i,nO - do a=nO+1,nV + do a=nO+1,nBas-nR jADC = jADC + 1 B_ADC(jADC,jADC) = e(i) + e(j) - e(a) @@ -188,15 +189,15 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) iADC = nH + nP do i=nC+1,nO - do a=nO+1,nV-nR - do b=a,nV-nR + do a=nO+1,nBas-nR + do b=a,nBas-nR iADC = iADC + 1 jADC = nH + nP do j=nC+1,nO - do c=nO+1,nV - do d=c,nV-nR + do c=nO+1,nBas-nR + do d=c,nBas-nR jADC = jADC + 1 B_ADC(iADC,jADC) = B_ADC(iADC,jADC) & @@ -221,14 +222,14 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) do i=nC+1,nO do j=i,nO - do a=nO+1,nV-nR + do a=nO+1,nBas-nR iADC = iADC + 1 jADC = nH + nP + nH*nPP do k=nC+1,nO do l=k,nO - do b=nO+1,nV-nR + do b=nO+1,nBas-nR jADC = jADC + 1 B_ADC(iADC,jADC) = B_ADC(iADC,jADC) & @@ -270,7 +271,7 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) do iADC=1,nADC - if(NORM2(X_ADC(1:nH+nP,iADC)) > 0.1d0 ) & +! if(NORM2(X_ADC(1:nH+nP,iADC)) > 0.1d0 ) & write(*,'(2(2X,F12.6))') e_ADC(iADC)*HaToeV,NORM2(X_ADC(1:nH+nP,iADC)) enddo @@ -290,8 +291,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) SigInf(:,:) = 0d0 do i=nC+1,nO - do p=nC+1,nV-nR - do q=nC+1,nV-nR + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR SigInf(p,q) = SigInf(p,q) - db_ERI(p,i,q,i) @@ -307,8 +308,8 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) if(e_ADC(iADC) > 0d0 ) cycle - do p=nC+1,nV-nR - do q=nC+1,nV-nR + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR G_ADC(p,q) = G_ADC(p,q) + X_ADC(p,iADC)*X_ADC(q,iADC) @@ -319,10 +320,10 @@ subroutine ADC2(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) ! Compute static self-energy for next iteration -- Eq. (25) -- - do p=nC+1,nV-nR - do q=nC+1,nV-nR - do r=nC+1,nV-nR - do s=nC+1,nV-nR + do r=nC+1,nBas-nR + do s=nC+1,nBas-nR + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR SigInf(p,q) = SigInf(p,q) + db_ERI(p,r,q,s)*G_ADC(r,s) diff --git a/src/QuAcK/CCSD.f90 b/src/QuAcK/CCSD.f90 index 039d368..bbe6d46 100644 --- a/src/QuAcK/CCSD.f90 +++ b/src/QuAcK/CCSD.f90 @@ -228,7 +228,7 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF) write(*,*)' CCSD energy ' write(*,*)'----------------------------------------------------' write(*,'(1X,A20,1X,F15.10)')' E(CCSD) = ',ECCSD - write(*,'(1X,A20,1X,F10.6)') ' Ec(CCSD) = ',EcCCSD + write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD) = ',EcCCSD write(*,*)'----------------------------------------------------' write(*,*) @@ -259,7 +259,7 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF) write(*,*)' CCSD(T) energy ' write(*,*)'----------------------------------------------------' write(*,'(1X,A20,1X,F15.10)')' E(CCSD(T)) = ',ECCSD + EcCCT - write(*,'(1X,A20,1X,F10.6)') ' Ec(CCSD(T)) = ',EcCCSD + EcCCT + write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD(T)) = ',EcCCSD + EcCCT write(*,*)'----------------------------------------------------' write(*,*) diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 2b22280..25770ed 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -26,8 +26,9 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & ! Local variables logical :: dRPA - integer :: ispin - double precision :: EcRPA + integer :: ispin,jspin + double precision :: EcRPA(nspin) + double precision :: EcBSE(nspin) double precision :: EcGM double precision,allocatable :: H(:,:) double precision,allocatable :: SigC(:) @@ -74,7 +75,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & ! Compute linear response call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) ! Compute correlation part of the self-energy @@ -104,7 +105,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & ! Dump results call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eG0W0,EcRPA,EcGM) + call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eG0W0,EcRPA(ispin),EcGM) ! Plot stuff @@ -119,8 +120,10 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & if(singlet_manifold) then ispin = 1 + EcBSE(ispin) = 0d0 + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif @@ -130,16 +133,27 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & if(triplet_manifold) then ispin = 2 + EcBSE(ispin) = 0d0 + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, & - rho(:,:,:,1),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,1),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + endif diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 7aa7204..e7443c0 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -40,7 +40,8 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani integer :: n_diis double precision :: rcond double precision :: Conv - double precision :: EcRPA + double precision :: EcRPA(nspin) + double precision :: EcBSE(nspin) double precision :: EcGM double precision :: lambda double precision,allocatable :: error_diis(:,:) @@ -110,7 +111,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(.not. GW0 .or. nSCF == 0) then call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) endif @@ -155,7 +156,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Print results call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA) + call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM) ! Linear mixing or DIIS extrapolation @@ -213,8 +214,10 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(singlet_manifold) then ispin = 1 + EcBSE(ispin) = 0d0 + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif @@ -223,16 +226,27 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(triplet_manifold) then ispin = 2 + EcBSE(ispin) = 0d0 + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + endif end subroutine evGW diff --git a/src/QuAcK/linear_response_B_matrix.f90 b/src/QuAcK/linear_response_B_matrix.f90 index 17e5a85..c91dc0a 100644 --- a/src/QuAcK/linear_response_B_matrix.f90 +++ b/src/QuAcK/linear_response_B_matrix.f90 @@ -32,7 +32,7 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,B_lr) delta_dRPA = 0d0 if(dRPA) delta_dRPA = 1d0 -! Build A matrix +! Build B matrix ia = 0 do i=nC+1,nO diff --git a/src/QuAcK/print_G0W0.f90 b/src/QuAcK/print_G0W0.f90 index e069f4c..d16769a 100644 --- a/src/QuAcK/print_G0W0.f90 +++ b/src/QuAcK/print_G0W0.f90 @@ -6,7 +6,10 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) include 'parameters.h' integer,intent(in) :: nBas,nO - double precision,intent(in) :: ENuc,EHF,EcRPA,EcGM + double precision,intent(in) :: ENuc + double precision,intent(in) :: EHF + double precision,intent(in) :: EcRPA + double precision,intent(in) :: EcGM double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) integer :: x,HOMO,LUMO diff --git a/src/QuAcK/print_RHF.f90 b/src/QuAcK/print_RHF.f90 index eef7055..1cda0ec 100644 --- a/src/QuAcK/print_RHF.f90 +++ b/src/QuAcK/print_RHF.f90 @@ -47,7 +47,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') 'MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,cHF) +! call matout(nBas,nBas,cHF) write(*,*) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') 'MO energies' diff --git a/src/QuAcK/print_evGW.f90 b/src/QuAcK/print_evGW.f90 index 54acbcb..eaed01f 100644 --- a/src/QuAcK/print_evGW.f90 +++ b/src/QuAcK/print_evGW.f90 @@ -6,7 +6,10 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) include 'parameters.h' integer,intent(in) :: nBas,nO,nSCF - double precision,intent(in) :: ENuc,EHF,EcRPA,EcGM + 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 diff --git a/src/QuAcK/self_energy_correlation.f90 b/src/QuAcK/self_energy_correlation.f90 index ccc3524..8f78d68 100644 --- a/src/QuAcK/self_energy_correlation.f90 +++ b/src/QuAcK/self_energy_correlation.f90 @@ -27,7 +27,8 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, ! Infinitesimal - eta = 0.001d0 + eta = 0.0d0 +! eta = 0.001d0 ! COHSEX static approximation @@ -68,7 +69,7 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, EcGM = 0d0 do i=nC+1,nO - EcGM = EcGM + SigC(i,i) + EcGM = EcGM + 0.5d0*SigC(i,i) enddo else diff --git a/src/QuAcK/self_energy_correlation_diag.f90 b/src/QuAcK/self_energy_correlation_diag.f90 index 8b4679b..ad3f6ef 100644 --- a/src/QuAcK/self_energy_correlation_diag.f90 +++ b/src/QuAcK/self_energy_correlation_diag.f90 @@ -68,7 +68,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega EcGM = 0d0 do i=nC+1,nO - EcGM = EcGM + SigC(i) + EcGM = EcGM + 0.5d0*SigC(i) enddo else @@ -117,7 +117,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega do b=nO+1,nBas-nR jb = jb + 1 eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)/eps + EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)/eps enddo enddo enddo @@ -164,7 +164,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega do b=nO+1,nBas-nR jb = jb + 1 eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM + 2d0*rho(a,i,jb)*rhox(a,i,jb)/eps + EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)/eps enddo enddo enddo