10
1
mirror of https://github.com/pfloos/quack synced 2025-04-01 06:21:37 +02:00

Merge branch 'master' of github.com:pfloos/QuAcK

This commit is contained in:
Antoine Marie 2025-03-26 22:48:54 +01:00
commit 9c923dfcb3
16 changed files with 14 additions and 1357 deletions

13
GoHu
View File

@ -1,13 +0,0 @@
#! /bin/bash
cp input/molecule.Hu input/molecule
cp input/basis.Hu input/basis
cp int/nBas.Hu.dat int/nBas.dat
cp int/ERI.Hu.dat int/ERI.dat
cp int/Kin.Hu.dat int/Kin.dat
cp int/Nuc.Hu.dat int/Nuc.dat
cp int/Ov.Hu.dat int/Ov.dat
cp int/x.Hu.dat int/x.dat
cp int/y.Hu.dat int/y.dat
cp int/z.Hu.dat int/z.dat
./bin/QuAcK

View File

@ -1,6 +0,0 @@
1 1
S 1
1 1.0000000000e+00 1.0000000000e+00
2 1
S 1
1 1.0000000000e+00 1.0000000000e+00

View File

@ -1,5 +0,0 @@
# nAt nEla nElb nCore nRyd
2 1 1 0 0
# Znuc x y z
X 0.0000000000 0.0000000000 0.0000000000
X 0.0000000000 0.0000000000 1.0000000000

View File

@ -1,2 +0,0 @@
# rs
1.0

View File

@ -1,321 +0,0 @@
#!/usr/bin/env python2
import sys
from termcolor import colored
import shlex
from subprocess import Popen, PIPE
import itertools
import re
import numpy as np
import os
from shutil import copy2
import matplotlib.pyplot as plt
import json
from math import *
from collections import OrderedDict
import csv
import argparse
def GetDuckDir():
return os.path.dirname(os.path.realpath(__file__))
def nNucl(molbaselines):
return float(molbaselines[1].split()[0])
def isMononucle(molbaselines):
return nNucl(molbaselines)==1
def openfileindir(path,readwrite):
mydir=os.path.dirname(path)
if not os.path.exists(mydir) and mydir!="":
os.makedirs(mydir)
return open(path,readwrite)
def outfile(Outdic,item,index=None):
itemdata=Outdic[item]
if itemdata["Enabled"]:
fmt=itemdata["Format"]
if index is not None:
filename=fmt.format(index)
else:
filename=fmt
if "Parent" in Outdic:
path=os.path.join(Outdic["Parent"],filename)
else:
path=filename
return openfileindir(path,'w')
else:
return
def runDuck(mol,basis,x,molbaselines,molbase,basisbase):
#gennerate molecule file
currdir=os.getcwd()
os.chdir(GetDuckDir())
molname='.'.join([mol,str(x)])
lstw=list()
for i,line in enumerate(molbaselines):
if i<3:
lstw.append(line)
else:
if isMononucle(molbaselines):
if i==3:
lstw.append(' '.join([str(x)]+line.split()[1:]))
else:
v=[float(abs(x))/float(2),float(-abs(x)/float(2))]
val=v[i-3]
lstw.append(' '.join([line.split()[0],'0.','0.',str(val)]))
junkfiles=list()
with open(molbase+molname,'w') as n:
junkfiles.append(n.name)
n.write(os.linesep.join(lstw))
#Copy basis
basisfile=basisbase+'.'.join([mol,basis])
newbasisfile=basisbase+'.'.join([molname,basis])
copy2(basisfile,newbasisfile)
junkfiles.append(newbasisfile)
#start child process Goduck
cmd=" ".join(["./GoDuck",molname, basis])
Duck=Popen(shlex.split(cmd),stdout=PIPE)
(DuckOut, DuckErr) = Duck.communicate()
excode=Duck.wait()
for junk in junkfiles:
os.remove(junk)
os.chdir(currdir)
return (excode,DuckOut,DuckErr)
def addvalue(dic,key,x,y):
if key not in dic:
dic[key]=list()
dic[key].append(y)
print(key)
print(x,y)
def main(mol):
#get basepath for files
molbase='examples/molecule.'
basisbase=molbase.replace('molecule','basis')
with open('PyOptions.json','r') as jfile:
options=json.loads(jfile.read())
basis=str(options['Basis'])
#Get mehtod to analyse
methodsdic=options['Methods']
#Get datas to analyse in this method
scandic=options['Scan']
scan=np.arange(scandic['Start'],scandic['Stop']+scandic['Step'],scandic['Step'])
print(scan)
mymethods=dict()
alllabels=list()
for method,methoddatas in methodsdic.iteritems():
if methoddatas['Enabled']:
mymethods[method]=methoddatas
for label,labeldatas in methoddatas['Labels'].iteritems():
if type(labeldatas) is dict:
enabled=labeldatas['Enabled']
else:
enabled=labeldatas
if enabled and label not in alllabels:
alllabels.append(label)
graphdic=dict()
errorconvstring="Convergence failed"
with open(os.path.join(GetDuckDir(),molbase+mol),'r') as b:
molbaselines=b.read().splitlines()
if isMononucle(molbaselines):
print('monoatomic system: variation of the nuclear charge')
else:
print('polyatomic system: variation is on the distance')
for x in scan:
(DuckExit,DuckOut,DuckErr)=runDuck(mol,basis,x,molbaselines,molbase,basisbase)
#print DuckOut on file or not
if "Outputs" in options:
outdat=options["Outputs"]
if 'DuckOutput' in outdat:
outopt=outdat["DuckOutput"]
if outopt['Enabled']:
if outopt['Multiple']:
duckoutf=outfile(outopt,"DuckOutput",x)
else:
if x==scan[0]:
duckoutf=outfile(outdat,"DuckOutput")
duckoutf.write('Z' if isMononucle(molbaselines) else 'Distance'+' '+str(x)+os.linesep+os.linesep)
duckoutf.write(DuckOut)
if outopt['Multiple']:
duckoutf.close()
print("GoDuk exit code " + str(DuckExit))
if DuckExit !=0:
#if GoDuck is not happy
print(DuckErr)
sys.exit(-1)
#get all data for the method
for method,methoddatas in mymethods.iteritems():
isnan=False
if '{0}' in method:
if "index" in methoddatas:
methodheaders=[method.format(str(x)) for x in methoddatas['Index']]
else:
try:
print(method)
reglist=re.findall('(\d+)'.join([re.escape(s) for s in method.split('{0}')]),DuckOut)
print(reglist)
final=max([(int(i[0]) if type(i) is tuple else int(i)) for i in reglist])
print(final)
methodheaders=[method.format(str(final))]
except:
isnan=True
methodheaders=[None]
method=method.replace('{0}','')
else:
methodheaders=list([method])
for methodheader in methodheaders:
if len(methodheaders)!=1:
method=methodheader
lbldic=methoddatas['Labels']
print(methodheader)
if methodheader is None:
methodtxt=''
else:
it=itertools.dropwhile(lambda line: methodheader + ' calculation' not in line , DuckOut.splitlines())
it=itertools.takewhile(lambda line: 'Total CPU time for ' not in line, it)
methodtxt=os.linesep.join(it)
if errorconvstring in methodtxt:
print(colored(' '.join([method, errorconvstring, '!!!!!']),'red'))
isnan=True
if methodtxt=='':
print(colored('No data' +os.linesep+ 'RHF scf not converged or method not enabled','red'))
isnan=True
#find the expected values
for label,labeldatas in lbldic.iteritems():
if type(labeldatas) is dict:
indexed=('Index' in labeldatas)
enabled=labeldatas['Enabled']
graph=labeldatas['Graph'] if 'Graph' in labeldatas else 1
else:
enabled=labeldatas
graph=1
indexed=False
if enabled:
if graph not in graphdic:
graphdic[graph]=OrderedDict()
y=graphdic[graph]
if not indexed:
v=np.nan
print(method)
print(label)
if not isnan:
try:
m=re.search('\s+'.join([re.escape(w) for w in label.split()]) + "\s+(?:"+re.escape("(eV):")+"\s+)?(?:=\s+)?(-?\d+.?\d*)",methodtxt)
v=m.group(1)
except:
v=np.nan
addvalue(y,(method,label),x,v)
else:
startindex=-1
columnindex=-1
linedtxt=methodtxt.split(os.linesep)
for n,line in enumerate(linedtxt):
if all(x in line for x in ['|',' '+label+' ','#']):
startindex=n+2
columnindex=[s.strip() for s in line.split('|')].index(label)
break
with open(os.path.join(GetDuckDir(),'input','molecule'),'r') as molfile:
molfile.readline()
line=molfile.readline()
nel=int(line.split()[1])
print(nel)
HOMO=int(nel/2)
HO=HOMO
LUMO=HOMO+1
BV=LUMO
for i in labeldatas['Index']:
v=np.nan
if type(i) is str or type(i) is unicode:
ival=eval(i)
if type(ival) is not int:
print('Index '+ str(i) + 'must be integer')
sys.exit(-2)
else:
ival=i
v=np.nan
if not isnan:
try:
if startindex!=-1 and columnindex!=-1:
line=linedtxt[startindex+ival-1]
v=float(line.split('|')[columnindex].split()[0])
print(method)
print(label)
print(i)
else:
v=np.nan
except:
v=np.nan
key=(method,label,i)
addvalue(y,key,x,v)
tpl=(x,scan.tolist().index(x)+1,len(y[key]))
print(tpl)
if tpl[1]-tpl[2]:
sys.exit()
#define graph grid
maxgraph=max(graphdic.keys())
maxrow=int(round(sqrt(maxgraph)))
maxcol=int(ceil(float(maxgraph)/float(maxrow)))
#define label ls
for graph,y in graphdic.iteritems():
datas=list()
datas.append(["#x"]+scan.tolist())
if len(y.keys())!=0:
plt.subplot(maxrow,maxcol,graph)
plt.xlabel('Z' if isMononucle(molbaselines) else 'Distance '+mol)
ylbls=list([basis])
for i in range(0,2):
lst=list(set([key[i] for key in y.keys()]))
if len(lst)==1:
ylbls.append(lst[0])
plt.ylabel(' '.join(ylbls))
print('Legend')
print(list(y.keys()))
for key,values in y.iteritems():
legend=list()
for el in key[0:2]:
if el not in ylbls:
legend.append(el)
if len(key)>2:
legend.append(str(key[2]))
#plot curves
lbl=' '.join(legend)
plt.plot(scan,y[key],'-o',label=lbl)
#print("min",x[y.index(min(y))]/2)
#generate legends
plt.legend()
dataout=False
if "Outputs" in options:
outputs=options['Outputs']
if "DataOutput" in outputs:
DataOutput=outputs['DataOutput']
dataout=DataOutput['Enabled']
if dataout:
fmtlegendf='{0}({1})'
datas.append([fmtlegendf.format("y",lbl)]+y[key])
if dataout:
csvdatas=zip(*datas)
with outfile(outputs,"DataOutput",graph) as csvf:
writer = csv.writer(csvf, delimiter=' ')
writer.writerow(['#']+ylbls)
writer.writerows(csvdatas)
#show graph
if "Outputs" in options:
outputs=options['Outputs']
if "FigureOutput" in outputs:
figout=outputs["FigureOutput"]
if figout["Enabled"]:
plt.savefig(figout['Path'])
plt.show()
if __name__ == '__main__':
parser=argparse.ArgumentParser()
parser.add_argument("mol",nargs='?', help="molecule to compute",type=str)
parser.add_argument("-c,--copy", help="Copy sample option file",action="store_true",dest="copy")
args = parser.parse_args()
if len(sys.argv)==1:
parser.print_help()
else:
if args.copy:
copy2(os.path.join(GetDuckDir(),"PyOptions.template.json"),"PyOptions.json")
if args.mol is not None:
os.system("vim PyOptions.json")
if args.mol is not None:
main(args.mol)

View File

@ -1,145 +0,0 @@
{
"Scan": {
"Start":1.8,
"Stop":1.9,
"Step":0.1
},
"Basis":"VDZ",
"Outputs": {
"DataOutput": {
"Enabled":true,
"Format":"Duck{0}.dat"
},
"DuckOutput": {
"Enabled":true,
"Multiple":false,
"Format":"DuckOut.out"
},
"FigureOutput":{
"Enabled":false,
"Path":"Figure.png"
}
},
"Methods": {
"RHF":{
"Enabled": true,
"Labels": {
"One-electron energy":false,
"Kinetic energy":false,
"Potential energy":false,
"Two-electron energy":false,
"Coulomb energy":false,
"Exchange energy":false,
"Electronic energy":false,
"Nuclear repulsion":false,
"Hartree-Fock energy":true,
"HF HOMO energy":false,
"HF LUMO energy":false,
"HF HOMO-LUMO gap":false
}
},
"One-shot G0W0": {
"Enabled": true,
"Labels": {
"G0W0 HOMO energy":true,
"G0W0 LUMO energy":true,
"G0W0 HOMO-LUMO gap":false,
"G0W0 total energy":false,
"RPA correlation energy" :false,
"Z": {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":1
},
"Sigma_c (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":2
},
"e_QP (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO+1","LUMO+2"],
"Graph":3
},
"e_HF (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":4
}
}
},
"Self-consistent evG{0}W{0}": {
"Enabled":false,
"Labels": {
"evGW HOMO energy":false,
"evGW LUMO energy":false,
"evGW HOMO-LUMO gap":false,
"evGW total energy":false,
"RPA correlation energy" :false,
"Z": {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":1
},
"Sigma_c (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":2
},
"e_QP (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":3
},
"e_HF (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":4
}
}
},
"Self-consistent qsG{0}W{0}": {
"Enabled": false,
"Labels": {
"qsGW HOMO energy":false,
"qsGW LUMO energy":false,
"qsGW HOMO-LUMO gap":false,
"qsGW total energy":false,
"qsGW exchange energy":false,
"qsGW correlation energy":false,
"RPA correlation energy":{
"Enabled":false,
"Graph":2
},
"Z": {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":4
},
"e_QP-e_HF (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":5
},
"e_QP (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":6
}
}
},
"MP2": {
"Enabled": false,
"Labels": {
"MP2 correlation energy": {
"Enabled":true,
"Graph":4
},
"Direct part":false,
"Exchange part":false,
"MP2 total energy":true,
"MP2 energy":false
}
}
}
}

View File

@ -1,145 +0,0 @@
{
"Scan": {
"Start":1.8,
"Stop":1.9,
"Step":0.1
},
"Basis":"VDZ",
"Outputs": {
"DataOutput": {
"Enabled":true,
"Format":"Duck{0}.dat"
},
"DuckOutput": {
"Enabled":true,
"Multiple":false,
"Format":"DuckOut.out"
},
"FigureOutput":{
"Enabled":false,
"Path":"Figure.png"
}
},
"Methods": {
"RHF":{
"Enabled": true,
"Labels": {
"One-electron energy":false,
"Kinetic energy":false,
"Potential energy":false,
"Two-electron energy":false,
"Coulomb energy":false,
"Exchange energy":false,
"Electronic energy":false,
"Nuclear repulsion":false,
"Hartree-Fock energy":true,
"HF HOMO energy":false,
"HF LUMO energy":false,
"HF HOMO-LUMO gap":false
}
},
"One-shot G0W0": {
"Enabled": true,
"Labels": {
"G0W0 HOMO energy":true,
"G0W0 LUMO energy":true,
"G0W0 HOMO-LUMO gap":false,
"G0W0 total energy":false,
"RPA correlation energy" :false,
"Z": {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":1
},
"Sigma_c (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":2
},
"e_QP (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO+1","LUMO+2"],
"Graph":3
},
"e_HF (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":4
}
}
},
"Self-consistent evG{0}W{0}": {
"Enabled":false,
"Labels": {
"evGW HOMO energy":false,
"evGW LUMO energy":false,
"evGW HOMO-LUMO gap":false,
"evGW total energy":false,
"RPA correlation energy" :false,
"Z": {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":1
},
"Sigma_c (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":2
},
"e_QP (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":3
},
"e_HF (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":4
}
}
},
"Self-consistent qsG{0}W{0}": {
"Enabled": false,
"Labels": {
"qsGW HOMO energy":false,
"qsGW LUMO energy":false,
"qsGW HOMO-LUMO gap":false,
"qsGW total energy":false,
"qsGW exchange energy":false,
"qsGW correlation energy":false,
"RPA correlation energy":{
"Enabled":false,
"Graph":2
},
"Z": {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":4
},
"e_QP-e_HF (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":5
},
"e_QP (eV)" : {
"Enabled":true,
"Index":["HOMO","LUMO","LUMO+1","LUMO+2"],
"Graph":6
}
}
},
"MP2": {
"Enabled": false,
"Labels": {
"MP2 correlation energy": {
"Enabled":true,
"Graph":4
},
"Direct part":false,
"Exchange part":false,
"MP2 total energy":true,
"MP2 energy":false
}
}
}
}

View File

@ -1,229 +0,0 @@
#! /bin/bash
INPUT=$1
echo
echo '******************************************'
echo '*** Extracting information of' $INPUT ' ***'
echo '******************************************'
echo
echo
echo '*** WFT information ***'
echo
grep "Hartree-Fock energy" $INPUT
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 -e "\t" $HF "\t" $G0F2 "\t" $GF2 "\t" $G0W0 "\t" $evGW
echo
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" $Tr_RPA_1 "\t" $Tr_RPA_3 "\t" $AC_RPA_1 "\t" $AC_RPA_3
echo
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" $Tr_RPAx_1 "\t" $Tr_RPAx_3 "\t" $AC_RPAx_1 "\t" $AC_RPAx_3
echo
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" $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 '*** 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" $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 '*** qsGW information: Tr@RPA (singlet), Tr@RPA (triplet), Tr@BSE (singlet), Tr@BSE (triplet), AC@BSE (singlet), AC@BSE (triplet) ***'
echo
Tr_RPA_qsGW_1=`grep "Tr@RPA@qsGW correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
Tr_RPA_qsGW_3=`grep "Tr@RPA@qsGW correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
Tr_BSE_qsGW_1=`grep "Tr@BSE@qsGW correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
Tr_BSE_qsGW_3=`grep "Tr@BSE@qsGW correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
AC_BSE_qsGW_1=`grep "AC@BSE@qsGW correlation energy (singlet) =" $INPUT| cut -f2 -d"="`
AC_BSE_qsGW_3=`grep "AC@BSE@qsGW correlation energy (triplet) =" $INPUT| cut -f2 -d"="`
echo -e "\t" $Tr_RPA_qsGW_1 "\t" $Tr_RPA_qsGW_3 "\t" $Tr_BSE_qsGW_1 "\t" $Tr_BSE_qsGW_3 "\t" $AC_BSE_qsGW_1 "\t" $AC_BSE_qsGW_3
echo
echo '*** CIS excitation energy (singlet & triplet) ***'
echo
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"|"`
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"|"`
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 '*** RPA excitation energy (singlet & triplet) ***'
echo
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"|"`
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"|"`
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
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 '*** BSE@evGW excitation energy (singlet & triplet) ***'
echo
evGW_1_1=`grep "| 1 |" $INPUT | head -9 | tail -1 | cut -f3 -d"|"`
evGW_1_2=`grep "| 2 |" $INPUT | head -9 | tail -1 | cut -f3 -d"|"`
evGW_1_3=`grep "| 3 |" $INPUT | head -9 | tail -1 | cut -f3 -d"|"`
evGW_1_4=`grep "| 4 |" $INPUT | head -9 | tail -1 | cut -f3 -d"|"`
evGW_1_5=`grep "| 5 |" $INPUT | head -9 | tail -1 | cut -f3 -d"|"`
evGW_3_1=`grep "| 1 |" $INPUT | head -10 | tail -1 | cut -f3 -d"|"`
evGW_3_2=`grep "| 2 |" $INPUT | head -10 | tail -1 | cut -f3 -d"|"`
evGW_3_3=`grep "| 3 |" $INPUT | head -10 | tail -1 | cut -f3 -d"|"`
evGW_3_4=`grep "| 4 |" $INPUT | head -10 | tail -1 | cut -f3 -d"|"`
evGW_3_5=`grep "| 5 |" $INPUT | head -10 | tail -1 | cut -f3 -d"|"`
echo -e "\t" $evGW_1_1 "\t" $evGW_3_1
echo -e "\t" $evGW_1_2 "\t" $evGW_3_2
echo -e "\t" $evGW_1_3 "\t" $evGW_3_3
echo -e "\t" $evGW_1_4 "\t" $evGW_3_4
echo -e "\t" $evGW_1_5 "\t" $evGW_3_5
echo
echo '*** BSE@qsGW excitation energy (singlet & triplet) ***'
echo
qsGW_1_1=`grep "| 1 |" $INPUT | head -11 | tail -1 | cut -f3 -d"|"`
qsGW_1_2=`grep "| 2 |" $INPUT | head -11 | tail -1 | cut -f3 -d"|"`
qsGW_1_3=`grep "| 3 |" $INPUT | head -11 | tail -1 | cut -f3 -d"|"`
qsGW_1_4=`grep "| 4 |" $INPUT | head -11 | tail -1 | cut -f3 -d"|"`
qsGW_1_5=`grep "| 5 |" $INPUT | head -11 | tail -1 | cut -f3 -d"|"`
qsGW_3_1=`grep "| 1 |" $INPUT | head -12 | tail -1 | cut -f3 -d"|"`
qsGW_3_2=`grep "| 2 |" $INPUT | head -12 | tail -1 | cut -f3 -d"|"`
qsGW_3_3=`grep "| 3 |" $INPUT | head -12 | tail -1 | cut -f3 -d"|"`
qsGW_3_4=`grep "| 4 |" $INPUT | head -12 | tail -1 | cut -f3 -d"|"`
qsGW_3_5=`grep "| 5 |" $INPUT | head -12 | tail -1 | cut -f3 -d"|"`
echo -e "\t" $qsGW_1_1 "\t" $qsGW_3_1
echo -e "\t" $qsGW_1_2 "\t" $qsGW_3_2
echo -e "\t" $qsGW_1_3 "\t" $qsGW_3_3
echo -e "\t" $qsGW_1_4 "\t" $qsGW_3_4
echo -e "\t" $qsGW_1_5 "\t" $qsGW_3_5
echo
echo '*** MATHEMATICA OUTPUT ***'
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 "\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 "\t" $evGW_1_1 "\t" $evGW_1_2 "\t" $evGW_1_3 "\t" $evGW_1_4 "\t" $evGW_1_5 "\t" $evGW_3_1 "\t" $evGW_3_2 "\t" $evGW_3_3 "\t" $evGW_3_4 "\t" $evGW_3_5 "\t" $Tr_RPA_qsGW_1 "\t" $Tr_RPA_qsGW_3 "\t" $Tr_BSE_qsGW_1 "\t" $Tr_BSE_qsGW_3 "\t" $AC_BSE_qsGW_1 "\t" $AC_BSE_qsGW_3 "\t" $qsGW_1_1 "\t" $qsGW_1_2 "\t" $qsGW_1_3 "\t" $qsGW_1_4 "\t" $qsGW_1_5 "\t" $qsGW_3_1 "\t" $qsGW_3_2 "\t" $qsGW_3_3 "\t" $qsGW_3_4 "\t" $qsGW_3_5
echo
echo '*** DONE ***'
echo

View File

@ -1,41 +0,0 @@
#! /bin/bash
Lmin=1
Lmax=1
Mmax=10
rs=$1
if [ $# != 1 ]
then
echo "Please, specify rs value"
else
echo "------------------------"
echo "Maxmium L value = " $Lmax
echo "Maxmium M value = " $Mmax
echo "------------------------"
echo
for (( L=$Lmin ; L<=$Lmax ; L++ )) ; do
ne=$(bc -l <<< "(2*($L+1)*($L+1))")
echo
echo "------------------------"
echo "Number of electrons = " $ne
echo "------------------------"
echo
for (( M=$L+1 ; M<=$Mmax ; M++ )) ; do
nb=$(bc -l <<< "(($M+1)*($M+1))")
echo "Number of basis functions = " $nb
echo -e "# rs \n" $rs > input/sph
./GoSph $ne $M > out/Sph_${ne}_${nb}.out
grep "Total CPU time for QuAcK =" out/Sph_${ne}_${nb}.out
done
done
fi

View File

@ -1,119 +0,0 @@
#! /bin/bash
MOL=$1
BASIS=$2
w_start=0.00
w_end=1.05
dw=0.05
w1=0.00
XF=$3
CF=$4
# for H
#aw1="1.49852 7.79815 25.1445"
#aw2="0.424545 -0.0382349 -0.32472"
# for He
#aw1="0.429447 0.053506 -0.339391"
#aw2="0.254939 -0.0893396 0.00765453"
# for H2
aw1="0.445525 0.0901503 -0.286898"
aw2="0.191734 -0.0364788 -0.017035"
# for Li
#aw1="0.055105 -0.00943825 -0.0267771"
#aw2="0.0359827 0.0096623 -0.0173542"
# for Li+
#aw1="0.503566, 0.137076, -0.348529"
#aw2="0.0553828, 0.00830375, -0.0234602"
# for B
#aw1="0.052676 -0.00624118 -0.000368825"
#aw2="0.0385558 -0.0015764 -0.000894297"
# for O
#aw1="-0.0187067 -0.0141017 -0.0100849"
#aw2="0.00544868 -0.0000118236 -0.000163245"
# for Al
#aw1="-0.00201219 -0.00371002 -0.00212719"
#aw2="-0.00117715 0.00188738 -0.000414761"
# for Be
#aw1="0.0663282, -0.0117682, -0.0335909"
#aw2="0.0479262, 0.00966351, -0.0208712"
DATA=${MOL}_${BASIS}_${XF}_${CF}_${w2}.dat
rm $DATA
touch $DATA
for w2 in $(seq $w_start $dw $w_end)
do
## w2=${w1}
echo "# Restricted or unrestricted KS calculation" > input/dft
echo " eDFT-UKS" >> input/dft
echo "# exchange rung:" >> input/dft
echo "# Hartree = 0" >> input/dft
echo "# LDA = 1: RS51,RMFL20" >> input/dft
echo "# GGA = 2: RB88" >> input/dft
echo "# Hybrid = 4" >> input/dft
echo "# Hartree-Fock = 666" >> input/dft
echo " 1 $XF " >> input/dft
echo "# correlation rung: " >> input/dft
echo "# Hartree = 0" >> input/dft
echo "# LDA = 1: RVWN5,RMFL20" >> input/dft
echo "# GGA = 2: " >> input/dft
echo "# Hybrid = 4: " >> input/dft
echo "# Hartree-Fock = 666" >> input/dft
echo " 0 $CF " >> input/dft
echo "# quadrature grid SG-n" >> input/dft
echo " 1" >> input/dft
echo "# Number of states in ensemble (nEns)" >> input/dft
echo " 3" >> input/dft
echo "# occupation numbers of orbitals nO and nO+1" >> input/dft
echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " " >> input/dft
echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " " >> input/dft
echo " 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo " 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 " >> input/dft
echo "# Ensemble weights: wEns(1),...,wEns(nEns-1)" >> input/dft
echo " ${w1} ${w2} " >> input/dft
echo "# Ncentered ? 0 for NO " >> input/dft
echo " 0 " >> input/dft
echo "# Parameters for CC weight-dependent exchange functional" >> input/dft
echo ${aw1} >> input/dft
echo ${aw2} >> input/dft
echo "# choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2" >> input/dft
echo "2" >> input/dft
echo "# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type" >> input/dft
echo " 1000 0.00001 T 5 1 1" >> input/dft
OUTPUT=${MOL}_${BASIS}_${XF}_${CF}_${w2}.out
./GoXC $MOL $BASIS > ${OUTPUT}
Ew=`grep "Ensemble energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
E0=`grep "Individual energy state 1:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
E1=`grep "Individual energy state 2:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
E2=`grep "Individual energy state 3:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
IP=`grep "Ionization Potential" ${OUTPUT} | grep " au" | tail -1 | cut -d":" -f 2 | sed 's/au//'`
EA=`grep "Electronic Affinity" ${OUTPUT} | grep " au" | tail -1 | cut -d":" -f 2 | sed 's/au//'`
FG=`grep "Fundamental Gap" ${OUTPUT} | grep " au" | tail -1 | cut -d":" -f 2 | sed 's/au//'`
Ex=`grep "Exchange energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/au//'`
HOMOa=`grep "HOMO a energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
LUMOa=`grep "LUMO a energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
HOMOb=`grep "HOMO a energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
LUMOb=`grep "LUMO b energy:" ${OUTPUT} | cut -d":" -f 2 | sed 's/eV//'`
echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex $HOMOa $LUMOa $HOMOb $LUMOb
echo $w1 $w2 $Ew $E0 $E1 $E2 $IP $EA $FG $Ex $HOMOa $LUMOa $HOMOb $LUMOb >> ${DATA}
done

View File

@ -134,9 +134,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
if(doufG0W0) then
call wall_time(start_GW)
! TODO
call ufRG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
! call eomRG0W0(dotest,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
@ -152,7 +150,6 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
if(doufGW) then
call wall_time(start_GW)
! TODO
call ufRGW(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)

View File

@ -1,315 +0,0 @@
subroutine eomRG0W0(dotest,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
! EOM version of G0W0
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: eHF(nOrb)
! Local variables
integer :: p
integer :: s
integer :: i,j,k,l
integer :: a,b,c,d
integer :: jb,kc,ia,ja
integer :: klc,kcd,ija,ijb,iab,jab
logical :: print_W = .false.
logical :: dRPA
integer :: isp_W
double precision :: EcRPA
integer :: n2h1p,n2p1h,nH
double precision,external :: Kronecker_delta
double precision,allocatable :: H(:,:)
double precision,allocatable :: cGW(:,:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: Z(:)
integer,allocatable :: order(:)
logical :: verbose = .false.
double precision,parameter :: cutoff1 = 0.01d0
double precision,parameter :: cutoff2 = 0.01d0
double precision :: eF
double precision,parameter :: window = 2.5d0
double precision :: start_timing,end_timing,timing
! Output variables
! Hello world
write(*,*)
write(*,*)'***********************************'
write(*,*)'* Restricted EOM-G0W0 Calculation *'
write(*,*)'***********************************'
write(*,*)
! Dimension of the supermatrix
n2h1p = nO*nO*nV
n2p1h = nV*nV*nO
nH = 1 + n2h1p + n2p1h
! Memory allocation
allocate(H(nH,nH),eGW(nH),cGW(nH,nH),Z(nH),order(nH))
! Initialization
dRPA = .true.
EcRPA = 0d0
eF = 0.5d0*(eHF(nO+1) + eHF(nO))
!-------------------------!
! Main loop over orbitals !
!-------------------------!
do p=nO,nO+1
H(:,:) = 0d0
!-----------------------------------------!
! Compute BSE supermatrix !
!-----------------------------------------!
! !
! | A V2h1p V2p1h 0 0 | !
! | | !
! | V2h1p A2h2p 0 B2h1p 0 | !
! | | !
! H = | V2p1h 0 A2p2h 0 B2p1h | !
! | | !
! | 0 0 0 0 0 | !
! | | !
! | 0 0 0 0 0 | !
! !
!-----------------------------------------!
call wall_time(start_timing)
!---------!
! Block F !
!---------!
H(1,1) = eHF(p)
!-------------!
! Block V2h1p !
!-------------!
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nOrb-nR
ija = ija + 1
H(1 ,1+ija) = sqrt(2d0)*ERI(p,a,i,j)
H(1+ija,1 ) = sqrt(2d0)*ERI(p,a,i,j)
! H(1+n2h1p+n2p1h+ija,1 ) = sqrt(2d0)*ERI(p,a,i,j)
! H(1+ija,1+n2h1p+n2p1h ) = sqrt(2d0)*ERI(p,a,i,j)
end do
end do
end do
!-------------!
! Block V2p1h !
!-------------!
iab = 0
do i=nC+1,nO
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
H(1 ,1+n2h1p+iab) = sqrt(2d0)*ERI(p,i,b,a)
H(1+n2h1p+iab,1 ) = sqrt(2d0)*ERI(p,i,b,a)
! H(1 ,1+2*n2h1p+n2p1h+iab) = sqrt(2d0)*ERI(p,i,b,a)
! H(1+2*n2h1p+n2p1h+iab,1 ) = sqrt(2d0)*ERI(p,i,b,a)
end do
end do
end do
!-------------!
! Block A2h1p !
!-------------!
ija = 0
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nOrb-nR
ija = ija + 1
klc = 0
do k=nC+1,nO
do l=nC+1,nO
do c=nO+1,nOrb-nR
klc = klc + 1
H(1+ija,1+klc) &
= ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) &
- 2d0*ERI(j,c,a,l) - 2d0*ERI(j,l,a,c))*Kronecker_delta(i,k)
! H(1+n2h1p+n2p1h+ija,1+n2h1p+n2p1h+klc) &
! = ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) &
! - 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k)
end do
end do
end do
end do
end do
end do
!-------------!
! Block A2p1h !
!-------------!
iab = 0
do i=nC+1,nO
do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR
iab = iab + 1
kcd = 0
do k=nC+1,nO
do c=nO+1,nOrb-nR
do d=nO+1,nOrb-nR
kcd = kcd + 1
H(1+n2h1p+iab,1+n2h1p+kcd) &
= ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
+ 2d0*ERI(a,k,i,c) + 2d0*ERI(a,c,i,k))*Kronecker_delta(b,d)
! H(1+2*n2h1p+n2p1h+iab,1+2*n2h1p+n2p1h+kcd) &
! = ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
! + 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d)
end do
end do
end do
end do
end do
end do
!-------------!
! Block B2h1p !
!-------------!
! ija = 0
! do i=nC+1,nO
! do j=nC+1,nO
! do a=nO+1,nOrb-nR
! ija = ija + 1
! kcd = 0
! do k=nC+1,nO
! do c=nO+1,nOrb-nR
! do d=nO+1,nOrb-nR
! kcd = kcd + 1
!
! H(1+ija,1+n2h1p+kcd) = - 2d0*ERI(j,k,a,c)
!
! end do
! end do
! end do
!
! end do
! end do
! end do
!-------------!
! Block B2p1h !
!-------------!
! iab = 0
! do i=nC+1,nO
! do a=nO+1,nOrb-nR
! do b=nO+1,nOrb-nR
! iab = iab + 1
! klc = 0
! do k=nC+1,nO
! do l=nC+1,nO
! do c=nO+1,nOrb-nR
! klc = klc + 1
! H(1+n2h1p+iab,1+klc) = - 2d0*ERI(a,c,i,l)
!
! end do
! end do
! end do
!
! end do
! end do
! end do
!-------------------------!
! Diagonalize supermatrix !
!-------------------------!
call wall_time(start_timing)
call diagonalize_general_matrix(nH,H,eGW,cGW)
do s=1,nH
order(s) = s
end do
call quick_sort(eGW,order,nH)
call set_order(cGW,order,nH,nH)
call wall_time(end_timing)
timing = end_timing - start_timing
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for construction of supermatrix = ',timing,' seconds'
write(*,*)
!-----------------!
! Compute weights !
!-----------------!
do s=1,nH
Z(s) = cGW(1,s)**2
end do
write(*,*)'-------------------------------------------'
write(*,'(1X,A32,I3,A8)')'| G0W0 energies (eV) for orbital',p,' |'
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
'|','#','|','e_QP','|','Z','|'
write(*,*)'-------------------------------------------'
do s=1,nH
! if(eGW(s) < eF .and. eGW(s) > eF - window) then
if(Z(s) > cutoff1) then
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|'
end if
end do
write(*,*)'-------------------------------------------'
write(*,*)
end do ! Loop on the orbital in the e block
end subroutine

View File

@ -1,6 +1,3 @@
! ---
subroutine diagonalize_nonsym_matrix(N, A, L, e_re, thr_d, thr_nd, thr_deg, imp_bio, verbose)
! Diagonalize a non-symmetric matrix A
@ -626,5 +623,3 @@ subroutine svd_local(A, LDA, U, LDU, D, Vt, LDVt, m, n)
enddo
end

View File

@ -1,5 +1,5 @@
# RHF UHF GHF ROHF
T F F F
# RHF UHF GHF ROHF HFB
T F F F F
# MP2 MP3
T T
# CCD pCCD DCD CCSD CCSD(T)
@ -12,11 +12,13 @@
T T T T
# G0F2 evGF2 qsGF2 ufGF2 G0F3 evGF3
T F F F F F
# G0W0 evGW qsGW SRG-qsGW ufG0W0 ufGW
T T F F F F
# G0W0 evGW qsGW ufG0W0 ufGW
T T F F F
# G0T0pp evGTpp qsGTpp ufG0T0pp
T F F F
# G0T0eh evGTeh qsGTeh
F F F
# Parquet
F
# Rtest Utest Gtest
T F F

View File

@ -1,11 +1,11 @@
# HF: maxSCF thresh DIIS guess mix shift stab search
10000 0.0000001 5 1 0.0 0.0 F F
# HF: maxSCF thresh DIIS guess mix shift stab search
10000 0.0000001 5 1 0.0 0.0 F F
# MP: reg
F
# CC: maxSCF thresh DIIS
64 0.0000001 5
# spin: TDA singlet triplet
F T T
# LR: TDA singlet triplet
F T T
# GF: maxSCF thresh DIIS lin eta renorm reg
256 0.00001 5 F 0.0 0 F
# GW: maxSCF thresh DIIS lin eta TDA_W reg
@ -16,3 +16,7 @@
F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA
F F F F T
# HFB: temperature sigma chem_pot_HF restart_HFB
0.05 1.00 T F
# Parquet: max_it_macro conv_one_body max_it_micro conv_two_body
1 0.00001 1 0 0.00001

Binary file not shown.