Initialization

HaToeV=27.21138602;

SizeLegend=1.5;

SetDirectory["~/Dropbox/Manuscripts/QUESTDB/Data"];

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

Flush[file_,shift_:0]:=Module[{DataFile,nState,nData,nDet,var,pt2,rpt2,t1,t2,rt2,t3,t4,res,res2,extrap,extrapQuad,fit,fitQuad,FormFit,FormFitQuad,nFit,MaxPts,Quality},
(* Import data *)
DataFile=Import[StringJoin[file,".dat"],"Data"];
{nState,nData} ={DataFile[[1,1]],DataFile[[2,1]]};
nDet=DataFile[[3;;nData+2,1]];
var=Table[DataFile[[nData+2+k]],{k,nData}];
pt2=Table[DataFile[[2nData+2+k,l]],{k,nData},{l,1,Length[DataFile[[2nData+2+k]]],2}];
rpt2=Table[DataFile[[3nData+2+k,l]],{k,nData},{l,1,Length[DataFile[[3nData+2+k]]],2}];

(* Compute excitation energies *)
nFit=2; MaxPts=5;
t1=Table[var[[k]][[Ordering[pt2[[k]]]]],{k,nData}];
t2=Table[Sort[pt2[[k]]],{k,nData}];
(*rt2=Table[Sort[rpt2[[k]]],{k,nData}];*)
(*Print["Sorted PT2 energies"];Print[TableForm[t2]];*)
(*Print["Data to be extrapolated"];Print[TableForm[Table[{SetAccuracy[t1[[k,n]],5],SetAccuracy[t2[[k,n]],5]},{n,nState},{k,nData-10,nData}]]];*)
FormFit=a+b e;
FormFitQuad=a2+b2 e+c2 (e^2);
fit=Table[FindFit[Take[Table[{t2[[k,n]]-t1[[k,n]],t2[[k,n]]},{k,nData}],-nPts],FormFit,{a,b},e],{n,nState},{nPts,nFit,nFit+MaxPts-1}];
fitQuad=Table[FindFit[Take[Table[{t2[[k,n]]-t1[[k,n]],t2[[k,n]]},{k,nData}],-(nPts+1)],FormFitQuad,{a2,b2,c2},e],{n,nState},{nPts,nFit,nFit+MaxPts-1}];
Print[TableForm[Table[{n,SetAccuracy[fit[[n,nPts]],5]},{n,nState},{nPts,MaxPts}],TableAlignments->Center,TableHeadings->{None,{"State", "Linear fit"}}]];
Print[TableForm[Table[{n,SetAccuracy[fitQuad[[n,nPts]],5]},{n,nState},{nPts,MaxPts}],TableAlignments->Center,TableHeadings->{None,{"State", "Quadratic fit"}}]];
t3=Table[FormFit/.fit[[n,nPts]]/.e->0,{n,nState},{nPts,MaxPts}];
t4=Table[FormFitQuad/.fitQuad[[n,nPts]]/.e->0,{n,nState},{nPts,MaxPts}];
res=Table[Flatten[{nDet[[k]],SetAccuracy[Table[shift+(t2[[k,n+1]]-t2[[k,1]])HaToeV,{n,nState-1}],3]}],{k,nData-8,nData}];
(*res2=Table[Flatten[{nDet[[k]],SetAccuracy[Table[shift+(rt2[[k,n+1]]-rt2[[k,1]])HaToeV,{n,nState-1}],3]}],{k,nData-8,nData}];*)extrap=SetAccuracy[Table[shift+(t3[[n+1,nPts]]-t3[[1,nPts]])HaToeV,{n,nState-1},{nPts,MaxPts}],3];
extrapQuad=SetAccuracy[Table[shift+(t4[[n+1,nPts]]-t4[[1,nPts]])HaToeV,{n,nState-1},{nPts,MaxPts}],3];
Print["Excitation energies with PT2"];
Print[TableForm[Take[res,-8],TableAlignments->Center,TableHeadings->{None,Flatten[{"nDet",Table[ToString[k],{k,nState-1}]}]}]];
(*Print["Excitation energies with rPT2"];
Print[TableForm[Take[res2,-8],TableAlignments->Center,TableHeadings->{None,Flatten[{"nDet",Table[ToString[k],{k,nState-1}]}]}]];
*)
Print["Linear fits"];Print[TableForm[extrap\[Transpose],TableAlignments->Center,TableHeadings->{Table[StringJoin[ToString[k+1]," pts"],{k,MaxPts}],Table[StringJoin["state ",ToString[k]],{k,nState-1}]}]];
Print["Quadratic fits"];
Print[TableForm[extrapQuad\[Transpose],TableAlignments->Center,TableHeadings->{Table[StringJoin[ToString[k+2]," pts"],{k,MaxPts}],Table[StringJoin["state ",ToString[k]],{k,nState-1}]}]];
(* Print stuff *)
Do[
Print[Show[{ListLogLinearPlot[Take[res[[;;,{1,k+1}]],-8],Joined->True,PlotMarkers->Automatic,PlotTheme->"Scientific",FrameLabel->{"Number of determinants","Excitation energy (eV)"},BaseStyle->16],Plot[{extrap[[k,1]],extrapQuad[[k,1]]},{x,0,(10^7)},PlotTheme->"Scientific"]}]];
,{k,nState-1}];

Return[":-)"]
];

Three non-H molecules

Summary

TableForm[{
{"2.18","2.14(0)","2.15","2.11(1)"},
{"3.81","3.79(1)","3.78","3.78(1)"},
{"3.08","3.06(0)","3.07","3.05(0)"},
{"2.47","2.46(0)","2.46","2.47(0)"},
{"1.60","1.60(0)","1.59","1.59(0)"},
{"5.06","5.04(0)","4.85","4.84(0)"},
{"5.33","5.32(0)","5.15","5.15(0)"},
{"3.47","3.49(0)","3.45","3.47(0)"},
{"4.35","4.34(0)","4.21","4.22(0)"},
{"5.12","5.12(0)","5.07","5.09(1)"},
{"2.71","2.71(1)","2.76","2.77(1)"},
{"2.69","2.68(2)","2.61","NA"},
{"4.40","4.46(1)","4.35","NA"},
{"1.20","1.22(3)","1.20","1.22(5)"},
{"4.34","4.36(1)","4.28","4.35(1)"},
{"6.03","6.00(3)","5.99","NA"},
{"5.65","5.65(1)","5.62","NA"},
{"3.63","3.62(1)","3.56","3.63(6)"},
{"3.93","3.95(2)","3.90","3.88(3)"},
{"2.45","2.47(5)","2.48","2.49(4)"},
{"","2.54(0)","","2.49(2)"},
{"","7.30(4)","","7.32(5)"},
{"","7.08(1)","","7.04(10)"},
{"","2.04(1)","","1.98(0)"}
},TableHeadings->{{"Silylidene","","HSiF","HPO","HPS","HPC","","","","(Subscript[CF, 2])","","(Subscript[CCl, 2])","","","","Formylfluoride","","CClF","(Subscript[SiCl, 2])","","HCF","F2CO","","HCCl"},{"CC3/AVDZ","exFCI/AVDZ","CC3/AVTZ","exFCI/AVTZ"}}]

Acetaldehyde

AVDZ ✓

Singlets

Flush["acetaldehyde_1_AVDZ"];

Triplets

Flush["acetaldehyde_3_AVDZ"];

AVTZ ✓

Singlets

Flush["acetaldehyde_1_AVTZ"];

Triplets

Flush["acetaldehyde_3_AVTZ"];

Cyanine

AVDZ ✓

Singlets

Flush["cyanine_1_AVDZ"];

Triplets

Flush["cyanine_3_AVDZ"];

AVTZ ✓

Singlets

Flush["cyanine_1_AVTZ"];

Triplets

Flush["cyanine_3_AVTZ"];

Cyclopropene

AVDZ ✓

Singlets

Flush["cyclopropene_1_AVDZ"];

Triplets

Flush["cyclopropene_3_AVDZ"];

AVTZ ✓

Singlets

Flush["cyclopropene_1_AVTZ"];

Triplets

Flush["cyclopropene_3_AVTZ"];

Diazomethane

AVDZ ✓

Singlets

Flush["diazomethane_1_AVDZ"];

Triplets

Flush["diazomethane_3_AVDZ"];

Fluo

Flush["diazomethane_Fl_AVDZ"];

AVTZ ✓

Singlets

Flush["diazomethane_1_AVTZ"];

Triplets

Flush["diazomethane_3_AVTZ"];

Fluo

Flush["diazomethane_Fl_AVTZ"];

Formamide

AVDZ ✓

Singlets

Flush["formamide_1_AVDZ"];

Triplets

Flush["formamide_3_AVDZ"];

AVTZ ✓

Singlets

Flush["formamide_1_AVTZ"];

Triplets

Flush["formamide_3_AVTZ"];

Ketene

AVDZ ✓

Singlets

Flush["ketene_1_AVDZ"];

Triplets

Flush["ketene_3_AVDZ"];

Fluo

Flush["ketene_Fl_AVDZ"];

AVTZ ✓

Singlets

Flush["ketene_1_AVTZ"];

Triplets

Flush["ketene_3_AVTZ"];

Fluo

Flush["ketene_Fl_AVTZ"];

Nitrosomethane

AVDZ ✓

Singlets

Flush["nitrosomethane_1_AVDZ"];

Triplets

Flush["nitrosomethane_3_AVDZ"];

Fluo

Flush["nitrosomethane_Fl_AVDZ"];

AVTZ ✓

Singlets

Flush["nitrosomethane_1_AVTZ"];

Triplets

Flush["nitrosomethane_3_AVTZ"];

Fluo

Flush["nitrosomethane_Fl_AVTZ"];

Figure

Table[{k,Rotate[MaTeX[ring[[k,1]]],π/2]},{k,Length[ring],2}]

ring=Import["rings.dat"];
ListPlot[{Table[{k,Around[ring[[k,5]]-ring[[k,4]],ring[[k,6]]]},{k,Length[ring]}],Table[{k,Around[ring[[k,7]]-ring[[k,4]],ring[[k,8]]]},{k,Length[ring]}]},PlotTheme->"Scientific",ImageSize->Large,BaseStyle->18,FrameTicks->{{Automatic,None},{Table[{k,Rotate[MaTeX[ring[[k,2]]],π/2]},{k,1,Length[ring]}],Table[{k,Rotate[MaTeX[ring[[k,1]],FontSize->14],π/2]},{k,1,Length[ring],2}]}},FrameLabel->{None,MaTeX["\\Delta E^\\text{FCI} - \\Delta E^\\text{CCSDT} \\text{ (eV)}",FontSize->18]}]
Export["errors.pdf",%]