PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-01-03 01:55:57 +01:00
Code Releases Activity

clean up

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-10-14 23:31:14 +02:00
parent 5cb83a6bbb
commit bf7cc8dada
6 changed files with 5185 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4431
basis/aug-cc-pvtz.native Normal file
View File

File diff suppressed because it is too large Load Diff

321
scripts/PyDuck Executable file
View File

@ -0,0 +1,321 @@
#!/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)

145
scripts/PyOptions.json Normal file
View File

@ -0,0 +1,145 @@
{
"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
}
}
}
}

145
scripts/PyOptions.template.json Normal file
View File

@ -0,0 +1,145 @@
{
"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
}
}
}
}

119
scripts/scan_w.sh Executable file
View File

@ -0,0 +1,119 @@
#! /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

45
src/eDFT/US51_lda_exchange_energy.f90
View File

@ -1,4 +1,4 @@
subroutine US51_lda_exchange_energy(nGrid,weight,rho,ExLDA)
subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex)
use xc_f90_lib_m
@ -15,40 +15,41 @@ subroutine US51_lda_exchange_energy(nGrid,weight,rho,ExLDA)
! Local variables
integer(8) :: nGri8
integer :: iG
double precision :: r
double precision,allocatable :: Ex(:)
TYPE(xc_f90_func_t) :: xc_func
TYPE(xc_f90_func_info_t) :: xc_info
integer :: func_id = 1
double precision :: alpha,r,alphaw,a2,b2,c2,a1,b1,c1
! Output variables
double precision :: ExLDA
double precision :: Ex
! Memory allocation
! Cxw2 parameters for He N->N+1
! a2 = 0.135068d0
! b2 = -0.00774769d0
! c2 = -0.0278205d0
nGri8 = int(nGrid,8)
print*,nGri8
allocate(Ex(nGrid))
! Cxw1 parameters for He N->N-1
! a1 = 0.420243d0
! b1 = 0.0700561d0
! c1 = -0.288301d0
call xc_f90_func_init(xc_func, func_id, XC_POLARIZED)
xc_info = xc_f90_func_get_info(xc_func)
call xc_f90_lda_exc(xc_func, nGri8, rho(1), Ex(1))
! Cx coefficient for Slater LDA exchange
ExLDA = 0d0
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
! do iG=1,nGrid
! alphaw = alpha*(1d0 - wEns(2)*(1d0 - wEns(2))*(a1 + b1*(wEns(2) - 0.5d0) + c1*(wEns(2) - 0.5d0)**2))
! Compute LDA exchange energy
! write(*,"(F8.6,1X,F9.6)") rho(iG), Ex(iG)
Ex = 0d0
do iG=1,nGrid
! ExLDA = ExLDA + weight(iG)*Ex(iG)
r = max(0d0,rho(iG))
! enddo
if(r > threshold) then
call xc_f90_func_end(xc_func)
Ex = Ex + weight(iG)*alpha*r**(4d0/3d0)
endif
enddo
end subroutine US51_lda_exchange_energy