From ff82232fe574d31e786eee8f1da0e0a52b170465 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 23 Sep 2019 13:31:01 +0200 Subject: [PATCH] fix bug in printing qsGW energy after BSE --- examples/molecule.H2 | 2 +- examples/molecule.LiF | 8 ++--- extract_sph.sh | 64 ---------------------------------------- input/methods | 2 +- input/molecule | 2 +- input/options | 2 +- scan_H2.sh | 6 +++- src/QuAcK/G0W0.f90 | 2 -- src/QuAcK/print_RHF.f90 | 18 +++++------ src/QuAcK/print_qsGW.f90 | 7 +++-- src/QuAcK/qsGW.f90 | 5 ++-- 11 files changed, 30 insertions(+), 88 deletions(-) delete mode 100755 extract_sph.sh diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 32aa31b..c248905 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. 0.6 + H 0. 0. 2.4 diff --git a/examples/molecule.LiF b/examples/molecule.LiF index 29e82af..af780ce 100644 --- a/examples/molecule.LiF +++ b/examples/molecule.LiF @@ -1,5 +1,5 @@ -# nAt nEl nCore nRyd - 2 12 4 0 +# nAt nEla nElb nCore nRyd + 2 6 6 0 0 # Znuc x y z - 3. 0. 0. 0.000 - 9. 0. 0. 2.955 + Li 0. 0. 0.000 + F 0. 0. 2.955 diff --git a/extract_sph.sh b/extract_sph.sh deleted file mode 100755 index 1714728..0000000 --- a/extract_sph.sh +++ /dev/null @@ -1,64 +0,0 @@ -#! /bin/bash - -INPUT=$1 - - echo - echo '******************************************' - echo '*** Extracting information of' $INPUT ' ***' - echo '******************************************' - echo - - echo - echo '*** WFT information ***' - grep "MP2 correlation energy" $INPUT - grep "Ec(MP2) =" $INPUT - grep "Ec(CCSD) =" $INPUT - grep "Ec(CCSD(T)) =" $INPUT - - echo - echo '*** Gap information: HF, G0F2, GF2, G0W0 & evGW ***' - HF=`grep "HF HOMO-LUMO gap (eV):" $INPUT | cut -f2 -d":"` - G0F2=`grep "GF2 HOMO-LUMO gap (eV):" $INPUT | head -1 | cut -f2 -d":"` - GF2=`grep "GF2 HOMO-LUMO gap (eV):" $INPUT | tail -1 | cut -f2 -d":"` - G0W0=`grep "G0W0 HOMO-LUMO gap (eV):" $INPUT | cut -f2 -d":"` - evGW=`grep "evGW HOMO-LUMO gap (eV):" $INPUT | tail -1 | cut -f2 -d":"` - - echo -e "\t" $HF "\t" $G0F2 "\t" $GF2 "\t" $G0W0 "\t" $evGW - - echo - echo '*** Ec@G0W0 information: RPA, GM, BSE1 & BSE3 ***' - RPA_G0W0=`grep "RPA@G0W0 correlation energy =" $INPUT| cut -f2 -d"="` - GM_G0W0=`grep "GM@G0W0 correlation energy =" $INPUT| cut -f2 -d"="` - BSE1_G0W0=`grep "BSE@G0W0 correlation energy (singlet)" $INPUT| cut -f2 -d"="` - BSE3_G0W0=`grep "BSE@G0W0 correlation energy (triplet)" $INPUT| cut -f2 -d"="` - - echo -e "\t" $RPA_G0W0 "\t" $GM_G0W0 "\t" $BSE1_G0W0 "\t" $BSE3_G0W0 - - echo - echo '*** Ec@evGW information: RPA, GM, BSE1 & BSE3 ***' - RPA_evGW=`grep "RPA@evGW correlation energy =" $INPUT | tail -1| cut -f2 -d"="` - GM_evGW=`grep "GM@evGW correlation energy =" $INPUT | tail -1 | cut -f2 -d"="` - BSE1_evGW=`grep "BSE@evGW correlation energy (singlet)" $INPUT | cut -f2 -d"="` - BSE3_evGW=`grep "BSE@evGW correlation energy (triplet)" $INPUT | cut -f2 -d"="` - - echo -e "\t" $RPA_evGW "\t" $GM_evGW "\t" $BSE1_evGW "\t" $BSE3_evGW - - echo - echo '*** CIS and TDHF excitation energy (singlet & triplet) ***' - CIS1=`grep "| 1 |" $INPUT | head -1 | cut -f4 -d"|"` - CIS3=`grep "| 1 |" $INPUT | head -2 | cut -f4 -d"|" | tail -1` - TDHF1=`grep "| 1 |" $INPUT | head -3 | cut -f4 -d"|" | tail -1` - TDHF3=`grep "| 1 |" $INPUT | head -4 | cut -f4 -d"|" | tail -1` - echo -e "\t" $CIS1 "\t" $CIS3 "\t" $TDHF1 "\t" $TDHF3 - - echo - echo '*** BSE@G0W0 and BSE@evGW excitation energy (singlet & triplet) ***' - G0W01=`grep "| 1 |" $INPUT | head -6 | cut -f4 -d"|" | tail -1` - G0W03=`grep "| 1 |" $INPUT | head -7 | cut -f4 -d"|" | tail -1` - evGW1=`grep "| 1 |" $INPUT | tail -2 | head -1 | cut -f4 -d"|"` - evGW3=`grep "| 1 |" $INPUT | tail -1 | cut -f4 -d"|"` - echo -e "\t" $G0W01 "\t" $G0W03 "\t" $evGW1 "\t" $evGW3 - - echo - echo '*** DONE ***' - diff --git a/input/methods b/input/methods index 40fd2a7..d324776 100644 --- a/input/methods +++ b/input/methods @@ -9,6 +9,6 @@ # GF2 GF3 F F # G0W0 evGW qsGW - T T F + T F T # MCMP2 F diff --git a/input/molecule b/input/molecule index 32aa31b..c248905 100644 --- a/input/molecule +++ b/input/molecule @@ -2,4 +2,4 @@ 2 1 1 0 0 # Znuc x y z H 0. 0. 0. - H 0. 0. 0.6 + H 0. 0. 2.4 diff --git a/input/options b/input/options index 0fe4f4b..19f6740 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.001 + 256 0.00001 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 622e930..2fbd9ac 100755 --- a/scan_H2.sh +++ b/scan_H2.sh @@ -2,8 +2,11 @@ MOL="H2" BASIS="VDZ" +R_START=0.5 +R_END=2.5 +DR=0.1 -for R in $(seq 0.5 0.1 0.6) +for R in $(seq $R_START $DR $R_END) do echo "# nAt nEla nElb nCore nRyd" > examples/molecule.$MOL echo " 2 1 1 0 0" >> examples/molecule.$MOL @@ -11,5 +14,6 @@ do echo " H 0. 0. 0." >> examples/molecule.$MOL echo " H 0. 0. $R" >> examples/molecule.$MOL ./GoDuck $MOL $BASIS > ${MOL}_${BASIS}_${R}.out + echo $R `./extract.sh ${MOL}_${BASIS}_${R}.out | tail -2 | head -1` done diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index d1bcc76..10474d1 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -46,8 +46,6 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & ! Hello world - print*,eta - write(*,*) write(*,*)'************************************************' write(*,*)'| One-shot G0W0 calculation |' diff --git a/src/QuAcK/print_RHF.f90 b/src/QuAcK/print_RHF.f90 index eef7055..bd0aa51 100644 --- a/src/QuAcK/print_RHF.f90 +++ b/src/QuAcK/print_RHF.f90 @@ -24,17 +24,17 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' Summary ' write(*,'(A50)') '---------------------------------------' - write(*,'(A32,1X,F16.10)') ' One-electron energy ',ET + EV - write(*,'(A32,1X,F16.10)') ' Kinetic energy ',ET - write(*,'(A32,1X,F16.10)') ' Potential energy ',EV + write(*,'(A32,1X,F16.10)') ' One-electron energy: ',ET + EV + write(*,'(A32,1X,F16.10)') ' Kinetic energy: ',ET + write(*,'(A32,1X,F16.10)') ' Potential energy: ',EV write(*,'(A50)') '---------------------------------------' - write(*,'(A32,1X,F16.10)') ' Two-electron energy ',EJ + EK - write(*,'(A32,1X,F16.10)') ' Coulomb energy ',EJ - write(*,'(A32,1X,F16.10)') ' Exchange energy ',EK + write(*,'(A32,1X,F16.10)') ' Two-electron energy: ',EJ + EK + write(*,'(A32,1X,F16.10)') ' Coulomb energy: ',EJ + write(*,'(A32,1X,F16.10)') ' Exchange energy: ',EK write(*,'(A50)') '---------------------------------------' - write(*,'(A32,1X,F16.10)') ' Electronic energy ',ERHF - write(*,'(A32,1X,F16.10)') ' Nuclear repulsion ',ENuc - write(*,'(A32,1X,F16.10)') ' Hartree-Fock energy ',ERHF + ENuc + write(*,'(A32,1X,F16.10)') ' Electronic energy: ',ERHF + write(*,'(A32,1X,F16.10)') ' Nuclear repulsion: ',ENuc + write(*,'(A32,1X,F16.10)') ' Hartree-Fock energy: ',ERHF + ENuc write(*,'(A50)') '---------------------------------------' write(*,'(A36,F13.6)') ' HF HOMO energy (eV):',eHF(HOMO)*HaToeV write(*,'(A36,F13.6)') ' HF LUMO energy (eV):',eHF(LUMO)*HaToeV diff --git a/src/QuAcK/print_qsGW.f90 b/src/QuAcK/print_qsGW.f90 index 81b6bc7..e163765 100644 --- a/src/QuAcK/print_qsGW.f90 +++ b/src/QuAcK/print_qsGW.f90 @@ -1,4 +1,4 @@ -subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigmaC,Z,EcRPA,EcGM) +subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigmaC,Z,EcRPA,EcGM,EqsGW) ! Print one-electron energies and other stuff for qsGW @@ -18,9 +18,12 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig ! Local variables integer :: x,HOMO,LUMO - double precision :: Gap,ET,EV,EJ,Ex,Ec,EqsGW + double precision :: Gap,ET,EV,EJ,Ex,Ec double precision,external :: trace_matrix +! Output variables + + double precision,intent(out) :: EqsGW ! HOMO and LUMO diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 18ac3e9..2e05450 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -40,6 +40,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani integer :: nBasSq integer :: ispin integer :: n_diis + double precision :: EqsGW double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcGM @@ -199,7 +200,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Print results call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,e,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM) + call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,e,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW) ! Increment @@ -273,7 +274,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy (triplet) =',EcBSE(2) write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2) - write(*,'(2X,A40,F15.6)') 'BSE@qsGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*)