4
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00
This commit is contained in:
Pierre-Francois Loos 2020-01-15 22:29:43 +01:00
parent 4f37593b5e
commit 4607acc649
22 changed files with 287 additions and 262 deletions

View File

@ -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. 3

View File

@ -2,4 +2,4 @@
2 7 7 0 0
# Znuc x y z
N 0. 0. 0.
N 0. 0. 2.
N 0. 0. 3

View File

@ -10,123 +10,164 @@ INPUT=$1
echo
echo '*** WFT information ***'
echo
grep "Hartree-Fock energy" $INPUT
EHF=`grep "Hartree-Fock energy" $INPUT | cut -f2 -d":"`
EHF=`grep "Hartree-Fock energy" $INPUT | cut -f2 -d"="`
grep "MP2 correlation energy" $INPUT
EcMP2=`grep "MP2 correlation energy" $INPUT | cut -f2 -d"="`
grep "Ec(MP2) =" $INPUT
grep "Ec(CCD) =" $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
# 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 -e "\t" $HF "\t" $G0F2 "\t" $GF2 "\t" $G0W0 "\t" $evGW
echo
echo '*** RPA@TDHF information: RPA, RPA1 & RPA3 ***'
RPA_TDHF=`grep "RPA@TDHF correlation energy =" $INPUT| cut -f2 -d"="`
RPA1_TDHF=`grep "RPA@TDHF correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
RPA3_TDHF=`grep "RPA@TDHF correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
echo '*** RPA information: Tr@RPA (singlet), Tr@RPA (triplet), AC@RPA (singlet), AC@RPA (triplet) ***'
echo
Tr_RPA_1=`grep "Tr@RPA correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
Tr_RPA_3=`grep "Tr@RPA correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
AC_RPA_1=`grep "AC@RPA correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
AC_RPA_3=`grep "AC@RPA correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
echo -e "\t" $RPA_TDHF "\t" $RPA1_TDHF "\t" $RPA3_TDHF
echo -e "\t" $Tr_RPA_1 "\t" $Tr_RPA_3 "\t" $AC_RPA_1 "\t" $AC_RPA_3
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 '*** RPAx information: Tr@RPAx (singlet), Tr@RPAx (triplet), AC@RPAx (singlet), AC@RPAx (triplet) ***'
echo
Tr_RPAx_1=`grep "Tr@RPAx correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
Tr_RPAx_3=`grep "Tr@RPAx correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
AC_RPAx_1=`grep "AC@RPAx correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
AC_RPAx_3=`grep "AC@RPAx correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
echo -e "\t" $RPA_G0W0 "\t" $GM_G0W0 "\t" $BSE1_G0W0 "\t" $BSE3_G0W0
echo -e "\t" $Tr_RPAx_1 "\t" $Tr_RPAx_3 "\t" $AC_RPAx_1 "\t" $AC_RPAx_3
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 '*** G0W0 information: Tr@RPA (singlet), Tr@RPA (triplet), Tr@BSE (singlet), Tr@BSE (triplet), AC@BSE (singlet), AC@BSE (triplet) ***'
echo
Tr_RPA_G0W0_1=`grep "Tr@RPA@G0W0 correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
Tr_RPA_G0W0_3=`grep "Tr@RPA@G0W0 correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
Tr_BSE_G0W0_1=`grep "Tr@BSE@G0W0 correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
Tr_BSE_G0W0_3=`grep "Tr@BSE@G0W0 correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
AC_BSE_G0W0_1=`grep "AC@BSE@G0W0 correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
AC_BSE_G0W0_3=`grep "AC@BSE@G0W0 correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
echo -e "\t" $RPA_evGW "\t" $GM_evGW "\t" $BSE1_evGW "\t" $BSE3_evGW
echo -e "\t" $Tr_RPA_G0W0_1 "\t" $Tr_RPA_G0W0_3 "\t" $Tr_BSE_G0W0_1 "\t" $Tr_BSE_G0W0_3 "\t" $AC_BSE_G0W0_1 "\t" $AC_BSE_G0W0_3
echo
echo '*** Ec@qsGW information: RPA, GM, BSE1 & BSE3 ***'
RPA_qsGW=`grep "RPA@qsGW correlation energy =" $INPUT | tail -1| cut -f2 -d"="`
GM_qsGW=`grep "GM@qsGW correlation energy =" $INPUT | tail -1 | cut -f2 -d"="`
BSE1_qsGW=`grep "BSE@qsGW correlation energy (singlet)" $INPUT | cut -f2 -d"="`
BSE3_qsGW=`grep "BSE@qsGW correlation energy (triplet)" $INPUT | cut -f2 -d"="`
echo '*** evGW information: Tr@RPA (singlet), Tr@RPA (triplet), Tr@BSE (singlet), Tr@BSE (triplet), AC@BSE (singlet), AC@BSE (triplet) ***'
echo
Tr_RPA_evGW_1=`grep "Tr@RPA@evGW correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
Tr_RPA_evGW_3=`grep "Tr@RPA@evGW correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
Tr_BSE_evGW_1=`grep "Tr@BSE@evGW correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
Tr_BSE_evGW_3=`grep "Tr@BSE@evGW correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
AC_BSE_evGW_1=`grep "AC@BSE@evGW correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
AC_BSE_evGW_3=`grep "AC@BSE@evGW correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
echo -e "\t" $RPA_qsGW "\t" $GM_qsGW "\t" $BSE1_qsGW "\t" $BSE3_qsGW
echo -e "\t" $Tr_RPA_evGW_1 "\t" $Tr_RPA_evGW_3 "\t" $Tr_BSE_evGW_1 "\t" $Tr_BSE_evGW_3 "\t" $AC_BSE_evGW_1 "\t" $AC_BSE_evGW_3
echo
echo '*** CIS and TDHF excitation energy (singlet & triplet) ***'
CIS1_1=`grep "| 1 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS3_1=`grep "| 1 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
TDHF1_1=`grep "| 1 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
TDHF3_1=`grep "| 1 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
echo -e "\t" $CIS1_1 "\t" $CIS3_1 "\t" $TDHF1_1 "\t" $TDHF3_1
echo '*** CIS excitation energy (singlet & triplet) ***'
echo
CIS1_2=`grep "| 2 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS3_2=`grep "| 2 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
TDHF1_2=`grep "| 2 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
TDHF3_2=`grep "| 2 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
echo -e "\t" $CIS1_2 "\t" $CIS3_2 "\t" $TDHF1_2 "\t" $TDHF3_2
CIS_1_1=`grep "| 1 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS_1_2=`grep "| 2 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS_1_3=`grep "| 3 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS_1_4=`grep "| 4 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS_1_5=`grep "| 5 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS1_3=`grep "| 3 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS3_3=`grep "| 3 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
TDHF1_3=`grep "| 3 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
TDHF3_3=`grep "| 3 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
echo -e "\t" $CIS1_3 "\t" $CIS3_3 "\t" $TDHF1_3 "\t" $TDHF3_3
CIS_3_1=`grep "| 1 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
CIS_3_2=`grep "| 2 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
CIS_3_3=`grep "| 3 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
CIS_3_4=`grep "| 4 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
CIS_3_5=`grep "| 5 |" $INPUT | head -2 | tail -1 | cut -f3 -d"|"`
CIS1_4=`grep "| 4 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS3_4=`grep "| 4 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
TDHF1_4=`grep "| 4 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
TDHF3_4=`grep "| 4 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
echo -e "\t" $CIS1_4 "\t" $CIS3_4 "\t" $TDHF1_4 "\t" $TDHF3_4
CIS1_5=`grep "| 5 |" $INPUT | head -1 | cut -f3 -d"|"`
CIS3_5=`grep "| 5 |" $INPUT | head -2 | cut -f3 -d"|" | tail -1`
TDHF1_5=`grep "| 5 |" $INPUT | head -3 | cut -f3 -d"|" | tail -1`
TDHF3_5=`grep "| 5 |" $INPUT | head -4 | cut -f3 -d"|" | tail -1`
echo -e "\t" $CIS1_5 "\t" $CIS3_5 "\t" $TDHF1_5 "\t" $TDHF3_5
echo -e "\t" $CIS_1_1 "\t" $CIS_3_1
echo -e "\t" $CIS_1_2 "\t" $CIS_3_2
echo -e "\t" $CIS_1_3 "\t" $CIS_3_3
echo -e "\t" $CIS_1_4 "\t" $CIS_3_4
echo -e "\t" $CIS_1_5 "\t" $CIS_3_5
echo
echo '*** BSE@G0W0 and BSE@evGW excitation energy (singlet & triplet) ***'
G0W01_1=`grep "| 1 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
G0W03_1=`grep "| 1 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
evGW1_1=`grep "| 1 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
evGW3_1=`grep "| 1 |" $INPUT | tail -1 | cut -f3 -d"|"`
echo -e "\t" $G0W01_1 "\t" $G0W03_1 "\t" $evGW1_1 "\t" $evGW3_1
echo '*** RPA excitation energy (singlet & triplet) ***'
echo
G0W01_2=`grep "| 2 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
G0W03_2=`grep "| 2 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
evGW1_2=`grep "| 2 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
evGW3_2=`grep "| 2 |" $INPUT | tail -1 | cut -f3 -d"|"`
echo -e "\t" $G0W01_2 "\t" $G0W03_2 "\t" $evGW1_2 "\t" $evGW3_2
RPA_1_1=`grep "| 1 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
RPA_1_2=`grep "| 2 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
RPA_1_3=`grep "| 3 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
RPA_1_4=`grep "| 4 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
RPA_1_5=`grep "| 5 |" $INPUT | head -3 | tail -1 | cut -f3 -d"|"`
G0W01_3=`grep "| 3 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
G0W03_3=`grep "| 3 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
evGW1_3=`grep "| 3 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
evGW3_3=`grep "| 3 |" $INPUT | tail -1 | cut -f3 -d"|"`
echo -e "\t" $G0W01_3 "\t" $G0W03_3 "\t" $evGW1_3 "\t" $evGW3_3
RPA_3_1=`grep "| 1 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
RPA_3_2=`grep "| 2 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
RPA_3_3=`grep "| 3 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
RPA_3_4=`grep "| 4 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
RPA_3_5=`grep "| 5 |" $INPUT | head -4 | tail -1 | cut -f3 -d"|"`
G0W01_4=`grep "| 4 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
G0W03_4=`grep "| 4 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
evGW1_4=`grep "| 4 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
evGW3_4=`grep "| 4 |" $INPUT | tail -1 | cut -f3 -d"|"`
echo -e "\t" $G0W01_4 "\t" $G0W03_4 "\t" $evGW1_4 "\t" $evGW3_4
echo -e "\t" $RPA_1_1 "\t" $RPA_3_1
echo -e "\t" $RPA_1_2 "\t" $RPA_3_2
echo -e "\t" $RPA_1_3 "\t" $RPA_3_3
echo -e "\t" $RPA_1_4 "\t" $RPA_3_4
echo -e "\t" $RPA_1_5 "\t" $RPA_3_5
G0W01_5=`grep "| 5 |" $INPUT | head -6 | cut -f3 -d"|" | tail -1`
G0W03_5=`grep "| 5 |" $INPUT | head -7 | cut -f3 -d"|" | tail -1`
evGW1_5=`grep "| 5 |" $INPUT | tail -2 | head -1 | cut -f3 -d"|"`
evGW3_5=`grep "| 5 |" $INPUT | tail -1 | cut -f3 -d"|"`
echo -e "\t" $G0W01_5 "\t" $G0W03_5 "\t" $evGW1_5 "\t" $evGW3_5
echo
echo '*** RPAx excitation energy (singlet & triplet) ***'
echo
RPAx_1_1=`grep "| 1 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
RPAx_1_2=`grep "| 2 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
RPAx_1_3=`grep "| 3 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
RPAx_1_4=`grep "| 4 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
RPAx_1_5=`grep "| 5 |" $INPUT | head -5 | tail -1 | cut -f3 -d"|"`
RPAx_3_1=`grep "| 1 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
RPAx_3_2=`grep "| 2 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
RPAx_3_3=`grep "| 3 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
RPAx_3_4=`grep "| 4 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
RPAx_3_5=`grep "| 5 |" $INPUT | head -6 | tail -1 | cut -f3 -d"|"`
echo -e "\t" $RPAx_1_1 "\t" $RPAx_3_1
echo -e "\t" $RPAx_1_2 "\t" $RPAx_3_2
echo -e "\t" $RPAx_1_3 "\t" $RPAx_3_3
echo -e "\t" $RPAx_1_4 "\t" $RPAx_3_4
echo -e "\t" $RPAx_1_5 "\t" $RPAx_3_5
echo
echo '*** BSE@G0W0 excitation energy (singlet & triplet) ***'
echo
G0W0_1_1=`grep "| 1 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
G0W0_1_2=`grep "| 2 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
G0W0_1_3=`grep "| 3 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
G0W0_1_4=`grep "| 4 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
G0W0_1_5=`grep "| 5 |" $INPUT | head -7 | tail -1 | cut -f3 -d"|"`
G0W0_3_1=`grep "| 1 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
G0W0_3_2=`grep "| 2 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
G0W0_3_3=`grep "| 3 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
G0W0_3_4=`grep "| 4 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
G0W0_3_5=`grep "| 5 |" $INPUT | head -8 | tail -1 | cut -f3 -d"|"`
echo -e "\t" $G0W0_1_1 "\t" $G0W0_3_1
echo -e "\t" $G0W0_1_2 "\t" $G0W0_3_2
echo -e "\t" $G0W0_1_3 "\t" $G0W0_3_3
echo -e "\t" $G0W0_1_4 "\t" $G0W0_3_4
echo -e "\t" $G0W0_1_5 "\t" $G0W0_3_5
echo
echo '*** MATHEMATICA OUTPUT ***'
echo -e "\t" $EHF "\t" $CIS1_1 "\t" $CIS1_2 "\t" $CIS1_3 "\t" $CIS1_4 "\t" $CIS1_5 "\t" $CIS3_1 "\t" $CIS3_2 "\t" $CIS3_3 "\t" $CIS3_4 "\t" $CIS3_5 "\t" $RPA_TDHF "\t" $RPA1_TDHF "\t" $RPA3_TDHF "\t" $TDHF1_1 "\t" $TDHF1_2 "\t" $TDHF1_3 "\t" $TDHF1_4 "\t" $TDHF1_5 "\t" $TDHF3_1 "\t" $TDHF3_2 "\t" $TDHF3_3 "\t" $TDHF3_4 "\t" $TDHF3_5 "\t" $RPA_G0W0 "\t" $GM_G0W0 "\t" $BSE1_G0W0 "\t" $BSE3_G0W0 "\t" $G0W01_1 "\t" $G0W01_2 "\t" $G0W01_3 "\t" $G0W01_4 "\t" $G0W01_5 "\t" $G0W03_1 "\t" $G0W03_2 "\t" $G0W03_3 "\t" $G0W03_4 "\t" $G0W03_5 "\t" $RPA_evGW "\t" $GM_evGW "\t" $BSE1_evGW "\t" $BSE3_evGW "\t" $evGW1_1 "\t" $evGW1_2 "\t" $evGW1_3 "\t" $evGW1_4 "\t" $evGW1_5 "\t" $evGW3_1 "\t" $evGW3_2 "\t" $evGW3_3 "\t" $evGW3_4 "\t" $evGW3_5
echo '*** DONE ***'
echo
echo -e "\t" $EHF "\t" $EcMP2 "\t" $Tr_RPA_1 "\t" $Tr_RPA_3 "\t" $AC_RPA_1 "\t" $AC_RPA_3 "\t" $Tr_RPAx_1 "\t" $Tr_RPAx_3 "\t" $AC_RPAx_1 "\t" $AC_RPAx_3 "\t" $Tr_RPA_G0W0_1 "\t" $Tr_RPA_G0W0_3 "\t" $Tr_BSE_G0W0_1 "\t" $Tr_BSE_G0W0_3 "\t" $AC_BSE_G0W0_1 "\t" $AC_BSE_G0W0_3 "\t" $CIS_1_1 "\t" $CIS_1_2 "\t" $CIS_1_3 "\t" $CIS_1_4 "\t" $CIS_1_5 "\t" $CIS_3_1 "\t" $CIS_3_2 "\t" $CIS_3_3 "\t" $CIS_3_4 "\t" $CIS_3_5 "\t" $RPA_1_1 "\t" $RPA_1_2 "\t" $RPA_1_3 "\t" $RPA_1_4 "\t" $RPA_1_5 "\t" $RPA_3_1 "\t" $RPA_3_2 "\t" $RPA_3_3 "\t" $RPA_3_4 "\t" $RPA_3_5 "\t" $RPAx_1_1 "\t" $RPAx_1_2 "\t" $RPAx_1_3 "\t" $RPAx_1_4 "\t" $RPAx_1_5 "\t" $RPAx_3_1 "\t" $RPAx_3_2 "\t" $RPAx_3_3 "\t" $RPAx_3_4 "\t" $RPAx_3_5 "\t" $G0W0_1_1 "\t" $G0W0_1_2 "\t" $G0W0_1_3 "\t" $G0W0_1_4 "\t" $G0W0_1_5 "\t" $G0W0_3_1 "\t" $G0W0_3_2 "\t" $G0W0_3_3 "\t" $G0W0_3_4 "\t" $G0W0_3_5
echo
echo '*** DONE ***'
echo

View File

@ -1,11 +1,11 @@
# RHF UHF MOM
T F F
# MP2 MP3 MP2-F12
F F F
T F F
# CCD CCSD CCSD(T) ringCCD ladderCCD
F F F F F
# CIS RPA RPAx ppRPA ADC
F T F F F
T T T F F
# GF2 GF3
F F
# G0W0 evGW qsGW

View File

@ -2,4 +2,4 @@
2 7 7 0 0
# Znuc x y z
N 0. 0. 0.
N 0. 0. 2.
N 0. 0. 3

View File

@ -1,5 +1,5 @@
# RHF: maxSCF thresh DIIS n_diis guess_type ortho_type
64 0.0000001 T 5 2 1
64 0.0000001 T 5 1 1
# MP:
# CC: maxSCF thresh DIIS n_diis
@ -7,10 +7,10 @@
# CIS/TDHF/BSE: singlet triplet
T T
# GF: maxSCF thresh DIIS n_diis renormalization
1 0.00001 T 10 3
64 0.00001 T 10 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 5 F F T F F F F 0.000
# ACFDT: AC XBS
T F
T T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

View File

@ -1,10 +1,10 @@
#! /bin/bash
MOL="H2"
BASIS="VDZ"
BASIS="VTZ"
R_START=0.5
R_END=3.0
DR=0.01
DR=0.05
for R in $(seq $R_START $DR $R_END)
do
@ -14,6 +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`
echo $R `./extract.sh ${MOL}_${BASIS}_${R}.out | tail -4 | head -1`
done

View File

@ -3,8 +3,8 @@
MOL="N2"
BASIS="VDZ"
R_START=1.5
R_END=3.5
DR=0.1
R_END=3.0
DR=0.05
for R in $(seq $R_START $DR $R_END)
do
@ -14,6 +14,6 @@ do
echo " N 0. 0. 0." >> examples/molecule.$MOL
echo " N 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`
echo $R `./extract.sh ${MOL}_${BASIS}_${R}.out | tail -4 | head -1`
done

View File

@ -8,13 +8,14 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
! Input variables
logical,intent(in) :: singlet_manifold,triplet_manifold
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),eHF(nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
logical :: dRPA
logical :: dump_matrix = .false.
logical :: dump_trans = .false.
integer :: ispin
@ -33,10 +34,6 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
lambda = 1d0
! Switch on exchange for CIS
dRPA = .false.
! Memory allocation
allocate(A(nS,nS),Omega(nS))
@ -46,7 +43,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
if(singlet_manifold) then
ispin = 1
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
if(dump_matrix) then
print*,'CIS matrix (singlet state)'
@ -55,7 +52,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
endif
call diagonalize_matrix(nS,A,Omega)
call print_excitation('CIS ',ispin,nS,Omega)
call print_excitation('CIS ',ispin,nS,Omega)
if(dump_trans) then
print*,'Singlet CIS transition vectors'
@ -68,7 +65,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
if(triplet_manifold) then
ispin = 2
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
if(dump_matrix) then
print*,'CIS matrix (triplet state)'
@ -77,7 +74,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
endif
call diagonalize_matrix(nS,A,Omega)
call print_excitation('CIS ',ispin,nS,Omega)
call print_excitation('CIS ',ispin,nS,Omega)
if(dump_trans) then
print*,'Triplet CIS transition vectors'

View File

@ -97,16 +97,30 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani
! Solve the quasi-particle equation
eGW(:) = eHF(:) + Z(:)*SigC(:)
! Dump results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
! call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
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,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! 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,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
! Perform BSE calculation
if(BSE) then
@ -116,19 +130,10 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -153,10 +158,10 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

View File

@ -56,16 +56,16 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
EcMP2(1) = EcMP2(2) + EcMP2(3)
write(*,*)
write(*,'(A32)') '-----------------------'
write(*,'(A32)') ' MP2 calculation '
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy',EcMP2(1)
write(*,'(A32,1X,F16.10)') ' Direct part ',EcMP2(2)
write(*,'(A32,1X,F16.10)') ' Exchange part ',EcMP2(3)
write(*,'(A32)') '-----------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy',EHF + EcMP2(1)
write(*,'(A32,1X,F16.10)') ' MP2 total energy',ENuc + EHF + EcMP2(1)
write(*,'(A32)') '-----------------------'
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcMP2(1)
write(*,'(A32,1X,F16.10)') ' Direct part = ',EcMP2(2)
write(*,'(A32,1X,F16.10)') ' Exchange part = ',EcMP2(3)
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcMP2(1)
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcMP2(1)
write(*,'(A32)') '--------------------------'
write(*,*)
end subroutine MP2

View File

@ -61,7 +61,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERH
! Guess coefficients and eigenvalues
call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
call mo_guess(nBas,nO,guess_type,S,Hc,ERI,X,cp,F,Fp,e,c,P)
! ON(:) = 0d0
! do i=1,nO

View File

@ -76,10 +76,10 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@RPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@RPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F15.6)') 'Tr@RPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -98,10 +98,10 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'AC@RPA correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A40,F15.6)') 'AC@RPA correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@RPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@RPA correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F15.6)') 'AC@RPA correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@RPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

View File

@ -161,7 +161,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
! Print results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
! call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
! Linear mixing or DIIS extrapolation
@ -196,7 +196,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
! 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,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin))
! Did it actually converge?
@ -212,6 +212,17 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
endif
! Dump the RPA correlation energy
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Perform BSE calculation
if(BSE) then
@ -221,19 +232,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -258,10 +260,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

View File

@ -1,4 +1,4 @@
subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
subroutine huckel_guess(nBas,nO,S,Hc,ERI,X,cp,F,Fp,e,c)
! Hickel guess of the molecular orbitals for HF calculation
@ -11,13 +11,11 @@ subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: J(nBas,nBas)
double precision,intent(inout):: K(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(inout):: cp(nBas,nBas)
double precision,intent(inout):: F(nBas,nBas)
double precision,intent(inout):: Fp(nBas,nBas)
double precision,intent(inout):: e(nBas)
double precision,intent(inout):: P(nBas,nBas)
! Local variables
@ -30,43 +28,20 @@ subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
a = 1.75d0
Fp(:,:) = Hc(:,:)
F(:,:) = Hc(:,:)
do mu=1,nBas
do nu=mu+1,nBas
Fp(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
Fp(nu,mu) = Fp(mu,nu)
F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
F(nu,mu) = F(mu,nu)
enddo
enddo
Fp(:,:) = matmul(transpose(X(:,:)),matmul(Fp(:,:),X(:,:)))
Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! call Coulomb_matrix_AO_basis(nBas,P,ERI,J)
! call exchange_matrix_AO_basis(nBas,P,ERI,K)
! do mu=1,nBas
! Fp(mu,mu) = Hc(mu,mu) + J(mu,mu) + 0.5d0*K(mu,mu)
! do nu=mu+1,nBas
! Fp(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
! Fp(nu,mu) = Fp(mu,nu)
! enddo
! enddo
! Fp = matmul(transpose(X),matmul(Fp,X))
! cp(:,:) = Fp(:,:)
! call diagonalize_matrix(nBas,cp,e)
! c = matmul(X,cp)
end subroutine

View File

@ -59,6 +59,10 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
if(minval(Omega) < 0d0) &
call print_warning('You may have instabilities in linear response!!')
do ia=1,nS
if(Omega(ia) < 0d0) Omega(ia) = 0d0
end do
call ADAt(nS,AmB,sqrt(Omega),AmBSq)
Z = matmul(AmBSq,matmul(ApB,AmBSq))
@ -74,11 +78,11 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
Omega = sqrt(Omega)
XpY = matmul(transpose(Z),AmBSq)
call DA(nS,1d0/sqrt(abs(Omega)),XpY)
call DA(nS,1d0/sqrt(Omega),XpY)
call ADAt(nS,AmB,1d0/sqrt(Omega),AmBSq)
XmY = matmul(transpose(Z),AmBSq)
call DA(nS,sqrt(abs(Omega)),XmY)
call DA(nS,sqrt(Omega),XmY)
! Compute the RPA correlation energy

View File

@ -1,4 +1,4 @@
subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,X,cp,F,Fp,e,c,P)
! Guess of the molecular orbitals for HF calculation
@ -12,10 +12,9 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: J(nBas,nBas)
double precision,intent(inout):: K(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(inout):: cp(nBas,nBas)
double precision,intent(inout):: F(nBas,nBas)
double precision,intent(inout):: Fp(nBas,nBas)
double precision,intent(inout):: e(nBas)
double precision,intent(inout):: P(nBas,nBas)
@ -37,7 +36,7 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
elseif(guess_type == 2) then
call huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
call huckel_guess(nBas,nO,S,Hc,ERI,X,cp,F,Fp,e,c)
elseif(guess_type == 3) then
@ -52,5 +51,4 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P)
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
end subroutine

View File

@ -40,11 +40,11 @@ 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@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(*,*)'-------------------------------------------------------------------------------'
! 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(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end subroutine print_G0W0

View File

@ -21,25 +21,25 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF)
write(*,*)
write(*,'(A50)') '---------------------------------------'
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(*,'(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(*,'(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(*,'(A50)') '---------------------------------------'
write(*,'(A36,F13.6)') ' HF HOMO energy (eV):',eHF(HOMO)*HaToeV
write(*,'(A36,F13.6)') ' HF LUMO energy (eV):',eHF(LUMO)*HaToeV
write(*,'(A36,F13.6)') ' HF HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,'(A50)') '---------------------------------------'
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(*,'(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(*,'(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(*,'(A50)') '-----------------------------------------'
write(*,'(A36,F13.6)') ' HF HOMO energy (eV) = ',eHF(HOMO)*HaToeV
write(*,'(A36,F13.6)') ' HF LUMO energy (eV) = ',eHF(LUMO)*HaToeV
write(*,'(A36,F13.6)') ' HF HOMO-LUMO gap (eV) = ',Gap*HaToeV
write(*,'(A50)') '-----------------------------------------'
write(*,*)
! Print results

View File

@ -47,11 +47,11 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM)
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

View File

@ -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,EqsGW)
subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigC,Z,EcRPA,EcGM,EqsGW)
! Print one-electron energies and other stuff for qsGW
@ -13,7 +13,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
double precision,intent(in) :: eHF(nBas),eGW(nBas),c(nBas),P(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas)
double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas)
double precision,intent(in) :: Z(nBas),SigmaC(nBas,nBas)
double precision,intent(in) :: Z(nBas),SigC(nBas,nBas)
! Local variables
@ -35,8 +35,9 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
EV = trace_matrix(nBas,matmul(P,V))
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
EqsGW = ET + EV + EJ + Ex
! Ec = 0.5d0*trace_matrix(nBas,matmul(P,SigC))
Ec = 0d0
EqsGW = ET + EV + EJ + Ex + Ec
! Dump results
@ -69,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(*,*)
@ -94,7 +95,7 @@ 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(*,'(A32,1X,F16.10)') ' RPA corr. energy ',EcRPA
write(*,'(A50)') '---------------------------------------'
write(*,*)

View File

@ -204,7 +204,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
! Print results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
! 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)
! Increment
@ -218,7 +218,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
! Compute second-order correction of the Hermitization error
call qsGW_PT(nBas,nC,nO,nV,nR,nS,eGW,SigCm)
! call qsGW_PT(nBas,nC,nO,nV,nR,nS,eGW,SigCm)
! Compute the overlap between HF and GW orbitals
@ -242,6 +242,17 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
endif
! Dump RPA correlation energy
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW total energy =',ENuc + EqsGW + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Perform BSE calculation
if(BSE) then
@ -251,19 +262,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW total energy =',ENuc + EqsGW + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -288,10 +290,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)