StatMech v2.03 \302\251 D.A. Harrington 2021use DocumentTools in
T:= sscanf(GetProperty(txtTemp,'value'),"%f")[];
if T <= 0 then
Do(%txtError="Non-positive temperature reset to 298.15 K");
T:=298.15;
SetProperty(txtTemp,'value',"298.15");
end if;
end use;
use DocumentTools in
if GetProperty(txtName,'value')=NULL or Trim(GetProperty(txtName,'value'))="" then
Do(%txtError="Molecule must be named")
else
N:=Do(%txtN);
M:=Do(%txtMolarMass);
Info := cat("StatMech Calculation v2.02 ",StringTools:-FormatTime("%d %b %Y %T")," \134nGas phase ");
if N = 1 then
Info := cat(Info,GetProperty(txtName,'value')," atom")
else
if Do(%optLinear) then
Info := cat(Info,"linear molecule ",GetProperty(txtName,'value'))
else
Info := cat(Info,"nonlinear molecule ",GetProperty(txtName,'value'))
end if
end if;
Info := cat(Info,sprintf(".\134nAll data for 1 mole at 1 bar.\134nMolar Mass = %6.4f g/mol.\134n",M));
if N > 1 then
sigma:=Do(%txtSigma);
Info := cat(Info,"Symmetry Number = ",sigma,".\134n");
B:=Do(%txtRotationalConst1);
if Do(%optLinear) then
Info := cat(Info,sprintf("Rotational Constant = %8.6f/cm.\134n",B));
else
A:=Do(%txtRotationalConst2);
C:=Do(%txtRotationalConst3);
Info := cat(Info,sprintf("Rotational Constants: %8.6f/cm, %8.6f/cm, %8.6f/cm.\134n",B,A,C));
end if;
Info := cat(Info,"\134nVibrational data\134n # Energy*cm (degen)\134n");
for i to nfreq do
Info := cat(Info,sprintf("%2d.%9.2f (%1d)\134n",i,Freq[i,1],Freq[i,2]))
end do
end if;
Info := cat(Info,"\134nElectronic Data\134nLevel Energy*cm (degen)\134n");
for i to maxstate do
Info := cat(Info,sprintf("%3d.%12.2f (%1d)\134n",i-1,Energy[i,1],Energy[i,2]))
end do;
Info := cat(Info,"\134n");
#
# Do the partition functions first. Everything we need without knowing the
# temperature comes first
# Collect translational parameters
qtransdivN := (2 * Pi * M * k / (1000 * Na * h * h)) ^ 1.5 * k / P0;
# Collect rotational parameters
if N > 1 then
if Do(%optNonlinear) then
qr := (k / hc) ^ 1.5 * sqrt(Pi / (A * B * C)) / sigma;
Trexp := 1.5;
else
qr := (k / (hc * B)) / sigma;
Trexp := 1;
end if
end if;
ans := "";
#
# Now calculate the data for this temperature
T:= sscanf(GetProperty(txtTemp,'value'),"%f")[];
kT := k * T;
RT := R * T;
# Calc vibrational quantities
qvib := 1;
Uvib := 0;
Cpvib := 0;
if N > 1 then
for i to nfreq do
gi := Freq[i,2];
hcwkT := hc * Freq[i,1] / kT;
qvi := 1 / (1 - exp(-hcwkT));
qvib := qvib * qvi ^ gi;
exphcwkT := exp(hcwkT);
Ui := gi * hcwkT / (exphcwkT - 1);
Cpi := Ui ^ 2 / gi * exphcwkT;
Uvib := Uvib + Ui;
Cpvib := Cpvib + Cpi;
end do;
end if;
Uvib := Uvib * RT;
Cpvib := Cpvib * R;
Svib := Uvib / T + R * ln(qvib);
# calc electronic quantities
qel := 0;
Uel := 0;
Cpel := 0;
for i to maxstate do
gi := Energy[i,2];
q := gi * exp(-hc * Energy[i,1] / kT);
qel := qel + q;
Uel := Uel + q * Energy[i,1];
Cpel := Cpel + q * Energy[i,1] ^ 2;
end do;
Uel := Uel * hc * Na / qel;
Cpel := -Uel ^ 2 / (RT * T) + R * (hc / kT) ^ 2 / qel * Cpel;
Sel := Uel / T + R * ln(qel);
# calc rotational quantities
if N = 1 then
qrot := 1;
Urot := 0;
Srot := 0;
Cprot := 0;
else
qrot := qr * T ^ Trexp;
Urot := Trexp * RT;
Cprot := Trexp * R;
Srot := Urot / T + R * ln(qrot);
end if;
# calc translational quantities
qtdivN := qtransdivN * T ^ 2.5;
Utrans := 1.5 * RT;
Cptrans := 2.5 * R;
Strans := 2.5 * R + R * ln(qtdivN);
# calc totals and other functions
Utot := Utrans + Urot + Uvib + Uel;
Stot := Strans + Srot + Svib + Sel;
Cptot := Cptrans + Cprot + Cpvib + Cpel;
qdivN := qtdivN * qrot * qvib * qel;
Atrans := Utrans - T * Strans;
Arot := Urot - T * Srot;
Avib := Uvib - T * Svib;
Ael := Uel - T * Sel;
Atot := Utot - T * Stot;
Htrans := Utrans + RT;
Hrot := Urot;
Hvib := Uvib;
Hel := Uel;
Htot := Utot + RT;
Gtrans := Htrans - T * Strans;
Grot := Hrot - T * Srot;
Gvib := Hvib - T * Svib;
Gel := Hel - T * Sel;
Gtot := Htot - T * Stot;
# Now format the output
Tstr:=StringTools:-PadRight(sprintf("%4.2f K",T),15);
ans := cat(Tstr,"translational rotational vibrational electronic total\134n");
ans := cat(ans,sprintf(" q : %10.4e%10.3f %12.8f %8.3f %12.4e\134n",qtdivN*Na,qrot,qvib,qel,qdivN*Na));
ans := cat(ans,sprintf(" q/N : %10.4e %12.4e\134n",qtdivN,qdivN));
ans := cat(ans,sprintf(" S/(J/(K mol)): %8.3f %8.3f %8.3f %8.3f %8.3f\134n",Strans,Srot,Svib,Sel,Stot));
ans := cat(ans,sprintf(" Cp/(J/(K mol)): %8.3f %8.3f %8.3f %8.3f %8.3f\134n",Cptrans,Cprot,Cpvib,Cpel,Cptot));
ans := cat(ans,sprintf(" U/(kJ/mol) : %8.3f %8.3f %8.3f %8.3f %8.3f\134n",Utrans/1000,Urot/1000,Uvib/1000,Uel/1000,Utot/1000));
ans := cat(ans,sprintf(" H/(kJ/mol) : %8.3f %8.3f %8.3f %8.3f %8.3f\134n",Htrans/1000,Hrot/1000,Hvib/1000,Hel/1000,Htot/1000));
ans := cat(ans,sprintf(" A/(kJ/mol) : %8.3f %8.3f %8.3f %8.3f %8.3f\134n",Atrans/1000,Arot/1000,Avib/1000,Ael/1000,Atot/1000));
ans := cat(ans,sprintf(" G/(kJ/mol) : %8.3f %8.3f %8.3f %8.3f %8.3f\134n",Gtrans/1000,Grot/1000,Gvib/1000,Gel/1000,Gtot/1000));
Info:=cat(Info,ans);
lines:=StringTools:-CountCharacterOccurrences(Info,"\134n")+1;
Do(%txtOutput(VisibleRows)=lines);
Do(%txtOutput(visible)=true);
Do(%txtOutput=Info);
end if;
end use;use DocumentTools,StringTools in
if GetProperty(txtName,'value')=NULL or Trim(GetProperty(txtName,'value'))="" then Do(%txtError="Molecule must be named") else Do(%txtError="") end if;
end use;
use DocumentTools in
Do(%optLinear=true);
Do(%lblSigma(Enabled)=true);
Do(%txtSigma(Enabled)=true);
Do(%lblRotationalConst(Enabled)=true);
Do(%txtRotationalConst1(Enabled)=true);
Do(%txtRotationalConst2(Visible)=false);
Do(%txtRotationalConst3(Visible)=false);
NFreqChanged();
end use;
use DocumentTools in
M:= parse(GetProperty(txtMolarMass,'value'));
if type(M,positive) then Do(%txtError="") else Do(%txtError="Molar mass must be a positive number") end if;
end use;
use DocumentTools in
Do(%optNonlinear=true);
Do(%lblSigma(Enabled)=true);
Do(%txtSigma(Enabled)=true);
Do(%lblRotationalConst(Enabled)=true);
Do(%txtRotationalConst1(Enabled)=true);
Do(%txtRotationalConst2(Visible)=true);
Do(%txtRotationalConst3(Visible)=true);
Do(%txtRotationalConst2(Enabled)=true);
Do(%txtRotationalConst3(Enabled)=true);
NFreqChanged();
end use;
use DocumentTools in
N:= parse(GetProperty(txtN,'value'));
if not type(N,posint) then Do(%txtError="Number of atoms must be a positive integer")
else Do(%txtError="");
if N=1 then
Do(%optLinear=true);
Do(%optNonlinear(Enabled)=false);
Do(%optLinear(Enabled)=false);
Do(%lblSigma(Enabled)=false);
Do(%txtSigma(Enabled)=false);
Do(%lblRotationalConst(Enabled)=false);
Do(%txtRotationalConst1(Enabled)=false);
Do(%txtRotationalConst2(Visible)=false);
Do(%txtRotationalConst3(Visible)=false);
NFreqChanged();
elif N=2 then
Do(%optLinear=true);
Do(%optNonlinear(Enabled)=false);
Do(%lblSigma(Enabled)=true);
Do(%txtSigma(Enabled)=true);
Do(%lblRotationalConst(Enabled)=true);
Do(%txtRotationalConst1(Enabled)=true);
Do(%txtRotationalConst2(Visible)=false);
Do(%txtRotationalConst3(Visible)=false);
NFreqChanged();
else
Do(%optLinear(Enabled)=true);
Do(%optNonlinear(Enabled)=true);
Do(%optNonlinear=true);
Do(%lblSigma(Enabled)=true);
Do(%txtSigma(Enabled)=true);
Do(%lblRotationalConst(Enabled)=true);
Do(%txtRotationalConst1(Enabled)=true);
Do(%txtRotationalConst2(Visible)=true);
Do(%txtRotationalConst3(Visible)=true);
Do(%txtRotationalConst2(Enabled)=true);
Do(%txtRotationalConst3(Enabled)=true);
NFreqChanged();
end if
end if
end use;
RotationalVibrationalElectronicuse DocumentTools in
sigma:= parse(GetProperty(txtSigma,'value'));
if type(sigma,posint) then Do(%txtError="") else Do(%txtError="Symmetry number must be a positive integer") end if;
end use;
use DocumentTools in
B:= parse(GetProperty(txtRotationalConst1,'value'));
if type(B,positive) then Do(%txtError="") else Do(%txtError="Rotational constants must be positive numbers") end if;
end use;
use DocumentTools in
A:= parse(GetProperty(txtRotationalConst2,'value'));
if type(A,positive) then Do(%txtError="") else Do(%txtError="Rotational constants must be positive numbers") end if;
end use;
use DocumentTools in
C:= parse(GetProperty(txtRotationalConst3,'value'));
if type(C,positive) then Do(%txtError="") else Do(%txtError="Rotational constants must be positive numbers") end if;
end use;
use DocumentTools in
NFreqChanged();
end use;
use DocumentTools in
maxstate:=maxstate+1;
Energy:=Matrix(maxstate,2,Energy);
Energy[maxstate,1]:=0;
Energy[maxstate,2]:=1;
Do(%tblEnergy(VisibleRows)=maxstate);
end use;
use DocumentTools in
maxstate:=1;
Energy:=Matrix(1,2,[0,1]);
Do(%tblEnergy(VisibleRows)=1);
end use;
use DocumentTools in
Energy[1,1]:=0;
SetProperty(tblEnergy,update);
end use;