From 55f7bffd762be9be3525079a1f8bfeba604f3cf3 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 1 Oct 2019 14:53:36 +0200 Subject: [PATCH] fix bug in evGW for BSE Ec --- examples/molecule.H2 | 2 +- examples/molecule.N2 | 2 +- input/methods | 2 +- input/molecule | 2 +- input/options | 2 +- scan_H2.sh | 4 ++-- scan_N2.sh | 2 +- src/QuAcK/G0W0.f90 | 7 +++++-- src/QuAcK/evGW.f90 | 24 +++++++++++++++--------- src/QuAcK/qsGW.f90 | 4 ++++ 10 files changed, 32 insertions(+), 19 deletions(-) diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 53ae340..32aa31b 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. 2.9 + H 0. 0. 0.6 diff --git a/examples/molecule.N2 b/examples/molecule.N2 index 35e31a1..93b448f 100644 --- a/examples/molecule.N2 +++ b/examples/molecule.N2 @@ -2,4 +2,4 @@ 2 7 7 0 0 # Znuc x y z N 0. 0. 0. - N 0. 0. 2.0 + N 0. 0. 3.4 diff --git a/input/methods b/input/methods index 116afdf..46c8861 100644 --- a/input/methods +++ b/input/methods @@ -9,6 +9,6 @@ # GF2 GF3 F F # G0W0 evGW qsGW - T F T + T T F # MCMP2 F diff --git a/input/molecule b/input/molecule index 35e31a1..93b448f 100644 --- a/input/molecule +++ b/input/molecule @@ -2,4 +2,4 @@ 2 7 7 0 0 # Znuc x y z N 0. 0. 0. - N 0. 0. 2.0 + N 0. 0. 3.4 diff --git a/input/options b/input/options index 2d124bb..23aacb6 100644 --- a/input/options +++ b/input/options @@ -9,6 +9,6 @@ # GF: maxSCF thresh DIIS n_diis renormalization 64 0.00001 T 10 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta - 256 0.00001 T 5 F F T F F F F 0.000 + 256 0.0000001 T 5 F F T F F F F 0.000 # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/scan_H2.sh b/scan_H2.sh index 6cb632a..57b12ae 100755 --- a/scan_H2.sh +++ b/scan_H2.sh @@ -1,10 +1,10 @@ #! /bin/bash MOL="H2" -BASIS="VDZ" +BASIS="AVTZ" R_START=0.5 R_END=3.0 -DR=0.1 +DR=0.01 for R in $(seq $R_START $DR $R_END) do diff --git a/scan_N2.sh b/scan_N2.sh index 9fbb688..86a445e 100755 --- a/scan_N2.sh +++ b/scan_N2.sh @@ -2,7 +2,7 @@ MOL="N2" BASIS="VDZ" -R_START=3.1 +R_START=1.5 R_END=3.5 DR=0.1 diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 10474d1..9812534 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -126,6 +126,10 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & ispin = 1 EcBSE(ispin) = 0d0 + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, & + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI, & rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) @@ -144,7 +148,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI, & - rho(:,:,:,1),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif @@ -160,5 +164,4 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & endif - end subroutine G0W0 diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 16ca47f..651d398 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,5 +1,5 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -30,7 +30,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani double precision,intent(in) :: eG0W0(nBas) double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: H(nBas,nBas) - double precision,intent(in) :: ERI_MO_basis(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables @@ -111,16 +111,16 @@ 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, & + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI, & rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) endif ! Compute correlation part of the self-energy - call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) + 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_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin)) + if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin)) ! Correlation self-energy @@ -212,28 +212,34 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(BSE) then ! Singlet manifold + 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, & + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI, & + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI, & rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif ! Triplet manifold + 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, & + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI, & rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, & + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI, & rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 4f3421b..dbd32bf 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -247,6 +247,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ispin = 1 EcBSE(ispin) = 0d0 + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, & + 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,e,ERI_MO_basis, & rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin))