Initialization

HaToeV=27.21138602;

HaTokcal=627.509474;

Needs["MaTeX`"]
SetOptions[MaTeX,"Preamble"->{"\\usepackage{amssymb,amsmath,latexsym,amsfonts,amsthm,mathpazo,xcolor,bm,mhchem}"}];

Flush[file_]:=Module[{DataFile,nState,nData,nDet,var,pt2,t3,res,extrap,fit,FormFit,nFit},
(* Import data *)
Print[file];
DataFile=Flatten[Import[StringJoin[file,".dat"],"Data"]];
{nState,nData} ={DataFile[[1]],DataFile[[2]]};
nDet=DataFile[[3;;nData+2]];
var=ArrayReshape[DataFile[[nData+3;;nData+2+nData nState]],{nData,nState}];
pt2=ArrayReshape[DataFile[[nData+3+nData nState;;nData+2+2nData nState]],{nData,nState}];

(* Compute excitation energies *)
nFit=2;
Print["Data to be extrapolated"];Print[TableForm[Table[{pt2[[k,n]]-var[[k,n]],var[[k,n]]},{n,nState},{k,nData}]]];
FormFit=a+b e;
fit=Table[FindFit[Take[Table[{pt2[[k,n]]-var[[k,n]],var[[k,n]]},{k,nData}],-nFit],FormFit,{a,b},e],{n,nState}];
Print[TableForm[Table[{n,fit[[n]]},{n,nState}],TableAlignments->Center,TableHeadings->{None,{"State", "Linear fit"}}]];
t3=Table[FormFit/.fit[[n]]/.e->0,{n,nState}];
res=Table[Flatten[{nDet[[k]],SetAccuracy[Table[pt2[[k,n]],{n,nState}],9]}],{k,nData}];
extrap=Flatten[{∞,SetAccuracy[Table[t3[[n]],{n,nState}],9]}];
Print[TableForm[Append[res,Flatten[extrap]],TableAlignments->Center,TableHeadings->{None,Flatten[{"nDet",Table[ToString[k],{k,nState}]}]}]];
(* Print stuff *)

Print[Show[{ListLogLinearPlot[res,Joined->True,PlotMarkers->Automatic,PlotTheme->"Scientific",FrameLabel->{"Number of determinants","Total energy (Ha)"},BaseStyle->16],Plot[extrap[[2]],{x,0,(10^7)},PlotTheme->"Scientific"]}]];
(*Return[extrap[[2]]];*)
Return[pt2[[nData,1]]];
];

FlushAll[Basis_]:=Module[{Mol,DataFiles,Energy},
(*Mol={"Be","BeH","C","C2H2","C2H4","C2H6","CH","CH2_1A1","CH2_3B1","CH3","CH3Cl","CH4","CN","CO","CO2","CS","Cl","Cl2","ClF","ClO","F","F2","H","H2CO","H2O","H2O2","H2S","H3COH","H3CSH","HCN","HCO","HCl","HF","HOCl","Li","Li2","LiF","LiH","N","N2","N2H4","NH","NH2","NH3","NO","Na","Na2","NaCl","O","O2","OH","P","P2","PH2","PH3","S","S2","SO","SO2","Si","Si2","Si2H6","SiH2_1A1","SiH2_3B1","SiH3","SiH4","SiO"};*)
SetDirectory[StringJoin["~/Dropbox/Manuscripts/srDFT_G2/Data/",Basis]];
Run["for i in *.dat; do echo ${i%_*.*}; done > mol.list"];
Mol=Flatten[Import["mol.list","Table"]];
DataFiles=Table[StringJoin[Mol[[k]],"_",Basis],{k,Length[Mol]}];
Energy=Table[{Mol[[k]],Flush[DataFiles[[k]]]},{k,Length[Mol]}];
Export[StringJoin["G2","_",Basis,".dat"],Energy,"Data"];
Return[Energy]
];

FlushBk[Basis_]:=Module[{Mol,DataFiles,Energy},
SetDirectory[StringJoin["~/Dropbox/Manuscripts/G2/Data/",Basis]];
Run["for i in *.dat; do echo ${i%_*.*}; done > mol.list"];
Mol=Flatten[Import["mol.list","Table"]];
Energy=Table[{Mol[[k]],Import[StringJoin[Mol[[k]],"_",Basis,".dat"],"Data"][[1,1]]},{k,Length[Mol]}];
Export[StringJoin["G2","_",Basis,".dat"],Energy,"Data"];
Return[Energy]
];

Carbon dimer

TableForm[
{{"VDZ",130.0,1},{"VTZ",139.9,3},{"VQZ",143.3,2},{"VDZ+F12",142.3,3},{"VTZ+F12",145.3,2},{"Exp",146.9,5}}\[Transpose]
,TableHeadings->{{"FCIQMC","(Subscript[D, e]) (kcal/mol)","Error"},None}]

TableForm[SetAccuracy[{
{"VDZ",HaTokcal(75.73131243-2(37.76050112))},
{"VTZ",HaTokcal(75.78558176-2(37.78101699))},
{"VQZ",HaTokcal(75.80245258-2(37.78853024))},
{"V5Z",HaTokcal(75.80788445-2(37.79540094))},
{"VDZ+(Subscript[srLDA, HF])",HaTokcal(75.73131243+0.134458-2(37.76050112+0.059282))},
{"VDZ+(Subscript[srLDA, CIPSI])",HaTokcal(75.73131243+0.13247-2(37.76050112+0.058016))},
{"VTZ+(Subscript[srLDA, HF])",HaTokcal(75.78558176+0.08395-2(37.78101699+0.039938))},
{"VTZ+(Subscript[srLDA, CIPSI])",HaTokcal(75.78558176+0.08220-2(37.78101699+0.03910))}
}\[Transpose],2]
,TableHeadings->{{"CIPSI"},None}]

TableForm[SetAccuracy[{
{"VDZ",HaTokcal(75.73131243-2(37.76050112))},
{"VTZ",HaTokcal(75.78558176-2(37.78101699))},
{"VQZ",HaTokcal(75.80245258-2(37.78853024))},
{"V5Z",HaTokcal(75.80788445-2(37.79540094))},
{"VDZ+(Subscript[srPBE, HF])",HaTokcal(75.73131243+0.119176-2(37.76050112+0.051074))},
{"VDZ+(Subscript[srPBE, CIPSI])",HaTokcal(75.73131243+0.11798-2(37.76050112+0.050225))},
{"VTZ+(Subscript[srPBE, HF])",HaTokcal(75.78558176+0.08155-2(37.78101699+0.03882))},
{"VTZ+(Subscript[srPBE, CIPSI])",HaTokcal(75.78558176+0.08049-2(37.78101699+0.03835))}
}\[Transpose],2]
,TableHeadings->{{"CIPSI"},None}]

Dunning’s basis

G2 with exCI/pVDZ — DONE —

FlushAll["cc-pVDZ"]

FlushAll["vdz-bfd"]

G2 with exCI/pVTZ

FlushAll["G2_VTZ"]

G2 with exCI/pVQZ

FlushAll["cc-pVQZ"]

G2 with exCI/aug-cc-pVDZ-BFD

FlushAll["avdz-bfd"]

G2 with exCI/aug-cc-pVTZ-BFD

running: 312108

FlushAll["avtz-bfd"]

Shifted-Bk

CIPSI with cc-pVDZ with (Subscript[N, det])=(10^3)

FlushAll["CIPSI3-cc-pVDZ"]

CIPSI with cc-pVDZ with (Subscript[N, det])=(10^4)

FlushAll["CIPSI4-cc-pVDZ"]

shifted-Bk with cc-pVDZ with (Subscript[N, det])=(10^3)

FlushBk["Bk3-cc-pVDZ"]

shifted-Bk with cc-pVDZ with (Subscript[N, det])=(10^4)

FlushBk["Bk4-cc-pVDZ"]

CIPSI with cc-pVTZ with (Subscript[N, det])=(10^3)

FlushAll["CIPSI3-cc-pVTZ"]

shifted-Bk with cc-pVTZ with (Subscript[N, det])=(10^3)

FlushBk["Bk3-cc-pVTZ"]

CIPSI with cc-pVQZ with (Subscript[N, det])=(10^3)

FlushAll["CIPSI3-cc-pVQZ"]

shifted-Bk with cc-pVQZ with (Subscript[N, det])=(10^3)

FlushBk["Bk3-cc-pVQZ"]

Plots

SetDirectory["/Users/loos/Dropbox/G2_TestSet_CLI"];
Do[
runID=Evaluate[Symbol["run"<>ToString[ID]]];
runID=Import["run"<>ToString[ID]<>".dat","Data"];
Evaluate[Symbol["G2$Mol$"<>ToString[ID]]]=Table[{k,Rotate[MaTeX[StringJoin["\\ce{",runID[[k,1]],"}"]],π/2]},{k,55}];
Evaluate[Symbol["G2$AE$"<>ToString[ID]]]=Table[{k,runID[[k,2]]},{k,55}];
Evaluate[Symbol["G2$Legend$"<>ToString[ID]]]=MaTeX["\\text{"<>ToString[ID]<>"}"];
,{ID,Flatten[{142,144,Range[153,155],Range[166,168],Range[175,192]}]}
];

Show[{
ListPlot[{G2$AE$142,G2$AE$144,G2$AE$183,G2$AE$189},PlotLegends->Placed[{
MaTeX["\\text{exFCI/cc-pVDZ: MAD = $22.86$ kcal.mol$^{-1}$}"],
MaTeX["\\text{exFCI/cc-pVTZ: MAD = $9.38$ kcal.mol$^{-1}$}"],
MaTeX["\\text{exFCI+srLDA($10^{-4}$)/cc-pVDZ: MAD = $12.99$ kcal.mol$^{-1}$}"],
MaTeX["\\text{exFCI+srLDA($10^{-4}$)/cc-pVTZ: MAD = $4.57$ kcal.mol$^{-1}$}"]
},{Left,Bottom}],
Joined->True,PlotMarkers->{Automatic,Medium},PlotStyle->Thick,PlotTheme->"Scientific",GridLines->{Table[k,{k,55}],Automatic},FrameLabel->{None,MaTeX["\\text{Deviation from exp. atomization energy (kcal/mol)}",Magnification->1]},BaseStyle->14,ImageSize->Full,AspectRatio->1/4,FrameTicks->{{Automatic,Automatic},G2$Mol$142},PlotRange->{{0,56},All}],Plot[{-1,1},{x,0,56},Filling->{1->{2}},PlotStyle->{{Red,Thin},{Red,Thin}}]
}]

Show[{
ListPlot[{G2$AE$142,G2$AE$185,G2$AE$184,G2$AE$183},PlotLegends->Placed[{
MaTeX["\\text{exFCI/cc-pVDZ: MAD = $22.86$ kcal.mol$^{-1}$}"],
MaTeX["\\text{exFCI+srLDA(HF)/cc-pVDZ: MAD = $12.81$ kcal.mol$^{-1}$}"],
MaTeX["\\text{exFCI+srLDA($10^{-3}$)/cc-pVDZ: MAD = $12.51$ kcal.mol$^{-1}$}"],
MaTeX["\\text{exFCI+srLDA($10^{-4}$)/cc-pVDZ: MAD = $12.99$ kcal.mol$^{-1}$}"]
},{Left,Bottom}],
Joined->True,PlotMarkers->{Automatic,Medium},PlotStyle->Thick,PlotTheme->"Scientific",GridLines->{Table[k,{k,55}],Automatic},FrameLabel->{None,MaTeX["\\text{Deviation from exp. atomization energy (kcal/mol)}",Magnification->1]},BaseStyle->14,ImageSize->Full,AspectRatio->1/4,FrameTicks->{{Automatic,Automatic},G2$Mol$142},PlotRange->{{0,56},All}],Plot[{-1,1},{x,0,56},Filling->{1->{2}},PlotStyle->{{Red,Thin},{Red,Thin}}]
}]

Show[{
ListPlot[{G2$AE$144,G2$AE$186,G2$AE$190,G2$AE$189},PlotLegends->Placed[{
MaTeX["\\text{exFCI/cc-pVTZ: MAD = $9.38$ kcal.mol$^{-1}$}"],
MaTeX["\\text{exFCI+srLDA(HF)/cc-pVTZ: MAD = $6.56$ kcal.mol$^{-1}$}"],
MaTeX["\\text{exFCI+srLDA($10^{-3}$)/cc-pVTZ: MAD = $6.30$ kcal.mol$^{-1}$}"],
MaTeX["\\text{exFCI+srLDA($10^{-4}$)/cc-pVTZ: MAD = $5.00$ kcal.mol$^{-1}$}"]
},{Left,Bottom}],
Joined->True,PlotMarkers->{Automatic,Medium},PlotStyle->Thick,PlotTheme->"Scientific",GridLines->{Table[k,{k,55}],Automatic},FrameLabel->{None,MaTeX["\\text{Deviation from exp. atomization energy (kcal/mol)}",Magnification->1]},BaseStyle->14,ImageSize->Full,AspectRatio->1/4,FrameTicks->{{Automatic,Automatic},G2$Mol$142},PlotRange->{{0,56},All}],Plot[{-1,1},{x,0,56},Filling->{1->{2}},PlotStyle->{{Red,Thin},{Red,Thin}}]
}]

QMC

SetDirectory["/Users/loos/Dropbox/G2_TestSet_CLI"];
data=Import["G2_QMC_vdz-bfd.dat","Data"];

data2=Table[{{data[[l,3]]-data[[l,2]],data[[l,6]]},{data[[l,4]]-data[[l,2]],data[[l,8]]},{data[[l,5]]-data[[l,2]],data[[l,10]]}},{l,Length[data]}];
data3=Table[Fit[data2[[l]],{1,e},e]/.e->0,{l,Length[data]}];
Table[{data[[l,1]],data3[[l]]},{l,Length[data]}]
Export["toto.dat",%]

data2=Table[{{data[[l,4]]-data[[l,2]],data[[l,8]]},{data[[l,5]]-data[[l,2]],data[[l,10]]}},{l,Length[data]}];
data3=Table[Fit[data2[[l]],{1,e},e]/.e->0,{l,Length[data]}];
Table[{data[[l,1]],data3[[l]]},{l,Length[data]}]
Export["toto.dat",%]

srDFT

VDZ

SetDirectory["/Users/loos/Dropbox/Manuscripts/srDFT_G2/Data"];
exFCI$VDZ=Import["G2_cc-pVDZ.dat","Data"];
srDFT$VDZ=Import["G2_cc-pVDZ_srDFT.dat","Data"];
srLDA$HF$VDZ=Table[{srDFT$VDZ[[k,1]],srDFT$VDZ[[k,2]]},{k,Length[srDFT$VDZ]}];
srPBE$HF$VDZ=Table[{srDFT$VDZ[[k,1]],srDFT$VDZ[[k,3]]},{k,Length[srDFT$VDZ]}];
srLDA$FCI$VDZ=Table[{srDFT$VDZ[[k,1]],srDFT$VDZ[[k,4]]},{k,Length[srDFT$VDZ]}];
srPBE$FCI$VDZ=Table[{srDFT$VDZ[[k,1]],srDFT$VDZ[[k,5]]},{k,Length[srDFT$VDZ]}];

srLDA$HF$VDZ=Table[{exFCI$VDZ[[k,1]],exFCI$VDZ[[k,2]]+srLDA$HF$VDZ[[k,2]]},{k,Length[exFCI$VDZ]}];
srPBE$HF$VDZ=Table[{exFCI$VDZ[[k,1]],exFCI$VDZ[[k,2]]+srPBE$HF$VDZ[[k,2]]},{k,Length[exFCI$VDZ]}];
srLDA$FCI$VDZ=Table[{exFCI$VDZ[[k,1]],exFCI$VDZ[[k,2]]+srLDA$FCI$VDZ[[k,2]]},{k,Length[exFCI$VDZ]}];
srPBE$FCI$VDZ=Table[{exFCI$VDZ[[k,1]],exFCI$VDZ[[k,2]]+srPBE$FCI$VDZ[[k,2]]},{k,Length[exFCI$VDZ]}];

Export["G2_cc-pVDZ_srLDA_HF.dat",srLDA$HF$VDZ,"Data"];
Export["G2_cc-pVDZ_srPBE_HF.dat",srPBE$HF$VDZ,"Data"];
Export["G2_cc-pVDZ_srLDA_FCI.dat",srLDA$FCI$VDZ,"Data"];
Export["G2_cc-pVDZ_srPBE_FCI.dat",srPBE$FCI$VDZ,"Data"];

SetDirectory["/Users/loos/Dropbox/G2_TestSet_CLI"];
Do[
runID=Evaluate[Symbol["run"<>ToString[ID]]];
runID=Import["run"<>ToString[ID]<>".dat","Data"];
Evaluate[Symbol["G2$Mol$"<>ToString[ID]]]=Table[{k,Rotate[MaTeX[StringJoin["\\ce{",runID[[k,1]],"}"]],π/2]},{k,55}];
Evaluate[Symbol["G2$AE$"<>ToString[ID]]]=Table[{k,runID[[k,2]]},{k,55}];
Evaluate[Symbol["G2$Legend$"<>ToString[ID]]]=MaTeX["\\text{"<>ToString[ID]<>"}"];
,{ID,141,144}
];

Show[{
ListPlot[{G2$AE$141,G2$AE$142,G2$AE$143,G2$AE$144,G2$AE$150},PlotLegends->Placed[{
MaTeX["\\text{srLDA(HF)/cc-pVDZ: MAD = $10.83$ kcal.mol$^{-1}$}",Magnification->1.5],
MaTeX["\\text{srPBE(HF)/cc-pVDZ: MAD = $12.82$ kcal.mol$^{-1}$}",Magnification->1.5],
MaTeX["\\text{srLDA(FCI)/cc-pVDZ: MAD = $10.85$ kcal.mol$^{-1}$}",Magnification->1.5],
MaTeX["\\text{srPBE(FCI)/cc-pVDZ: MAD = $12.94$ kcal.mol$^{-1}$}",Magnification->1.5],
MaTeX["\\text{srPBE(FCI)/cc-pVDZ: MAD = $22.86$ kcal.mol$^{-1}$}",Magnification->1.5]
},{Left,Bottom}],
Joined->True,PlotMarkers->{Automatic,Medium},PlotStyle->Thick,PlotTheme->"Scientific",GridLines->{Table[k,{k,55}],Automatic},FrameLabel->{None,MaTeX["\\text{Deviation from exp. atomization energy (kcal.mol$^{-1}$)}",Magnification->1.25]},BaseStyle->14,ImageSize->Full,AspectRatio->1/4,FrameTicks->{{Automatic,Automatic},G2$Mol$141},PlotRange->{{0,56},All}],Plot[{-1,1},{x,0,56},Filling->{1->{2}},PlotStyle->{{Red,Thin},{Red,Thin}}]
}]

VTZ

SetDirectory["/Users/loos/Dropbox/Manuscripts/srDFT_G2/Data"];
exFCI$VTZ=Import["G2_cc-pVTZ.dat","Data"];
srDFT$VTZ=Import["G2_cc-pVTZ_srDFT.dat","Data"];
srLDA$HF$VTZ=Table[{srDFT$VTZ[[k,1]],srDFT$VTZ[[k,2]]},{k,Length[srDFT$VTZ]}];
srPBE$HF$VTZ=Table[{srDFT$VTZ[[k,1]],srDFT$VTZ[[k,3]]},{k,Length[srDFT$VTZ]}];
srLDA$FCI$VTZ=Table[{srDFT$VTZ[[k,1]],srDFT$VTZ[[k,4]]},{k,Length[srDFT$VTZ]}];
srPBE$FCI$VTZ=Table[{srDFT$VTZ[[k,1]],srDFT$VTZ[[k,5]]},{k,Length[srDFT$VTZ]}];

srLDA$HF$VTZ=Table[{exFCI$VTZ[[k,1]],exFCI$VTZ[[k,2]]+srLDA$HF$VTZ[[k,2]]},{k,Length[exFCI$VTZ]}];
srPBE$HF$VTZ=Table[{exFCI$VTZ[[k,1]],exFCI$VTZ[[k,2]]+srPBE$HF$VTZ[[k,2]]},{k,Length[exFCI$VTZ]}];
srLDA$FCI$VTZ=Table[{exFCI$VTZ[[k,1]],exFCI$VTZ[[k,2]]+srLDA$FCI$VTZ[[k,2]]},{k,Length[exFCI$VTZ]}];
srPBE$FCI$VTZ=Table[{exFCI$VTZ[[k,1]],exFCI$VTZ[[k,2]]+srPBE$FCI$VTZ[[k,2]]},{k,Length[exFCI$VTZ]}];

Export["G2_cc-pVTZ_srLDA_HF.dat",srLDA$HF$VTZ,"Data"];
Export["G2_cc-pVTZ_srPBE_HF.dat",srPBE$HF$VTZ,"Data"];
Export["G2_cc-pVTZ_srLDA_FCI.dat",srLDA$FCI$VTZ,"Data"];
Export["G2_cc-pVTZ_srPBE_FCI.dat",srPBE$FCI$VTZ,"Data"];

SetDirectory["/Users/loos/Dropbox/G2_TestSet_CLI"];
Do[
runID=Evaluate[Symbol["run"<>ToString[ID]]];
runID=Import["run"<>ToString[ID]<>".dat","Data"];
Evaluate[Symbol["G2$Mol$"<>ToString[ID]]]=Table[{k,Rotate[MaTeX[StringJoin["\\ce{",runID[[k,1]],"}"]],π/2]},{k,55}];
Evaluate[Symbol["G2$AE$"<>ToString[ID]]]=Table[{k,runID[[k,2]]},{k,55}];
Evaluate[Symbol["G2$Legend$"<>ToString[ID]]]=MaTeX["\\text{"<>ToString[ID]<>"}"];
,{ID,145,150}
];

SizeLegend=1.5;
Show[{
ListPlot[{G2$AE$145,G2$AE$146,G2$AE$147,G2$AE$148,G2$AE$149},PlotLegends->Placed[{
MaTeX["\\text{srLDA(HF)/cc-pVTZ: MAD = $5.20$ kcal.mol$^{-1}$}",Magnification->SizeLegend],
MaTeX["\\text{srPBE(HF)/cc-pVTZ: MAD = $4.11$ kcal.mol$^{-1}$}",Magnification->SizeLegend],
MaTeX["\\text{srLDA(FCI)/cc-pVTZ: MAD = $5.49$ kcal.mol$^{-1}$}",Magnification->SizeLegend],
MaTeX["\\text{srPBE(FCI)/cc-pVTZ: MAD = $4.29$ kcal.mol$^{-1}$}",Magnification->SizeLegend],
MaTeX["\\text{exFCI/cc-pVTZ: MAD = $8.03$ kcal.mol$^{-1}$}",Magnification->SizeLegend]
},{Left,Bottom}],
Joined->True,PlotMarkers->{Automatic,Medium},PlotStyle->Thick,PlotTheme->"Scientific",GridLines->{Table[k,{k,55}],Automatic},FrameLabel->{None,MaTeX["\\text{Deviation from exp. atomization energy (kcal.mol$^{-1}$)}",Magnification->1.25]},BaseStyle->14,ImageSize->Full,AspectRatio->1/4,FrameTicks->{{Automatic,Automatic},G2$Mol$145},PlotRange->{{0,56},All}],Plot[{-1,1},{x,0,56},Filling->{1->{2}},PlotStyle->{{Red,Thin},{Red,Thin}}]
}]