The parameters are extracted from G3 and 6-311+G(d) computations from Gaussian03 or Gamess. There are four parameters to calibrate: k0,k1,sig0,sig1. We need four independent data and select: Etot, Vne, Vee and C-H bond length. Because these numerical values differ by more than a factor 100 between the largest and the smallest one has to optimize convenient ratios, all in the neighborhood of ~1, see procedure 'init', pardat[3], pardat[4], below. Then a least squares procedure 'Minsqr' adjusts the results of CH4 Kimball optimizations, procedures 'Calculate', 'Minimize' and function subroutine 'E(u)', until Etot, Vne, Vnn, and d(CH) have reasonably converged to the given values. You can view the progress of this convergence in the result link. Here we needed about 151000 CH4 energy minimizations to converge, duration ~10 sec.
{**** CH4 **** ES 29.10.1998}
program CH4exp;
const tb = ' ';
bohr = 0.529177249; {Angstr./bohr}
hartree = 27.2114; {eV/hartree}
kcal = 23.06035; {kcal/eV}
crit = 2.4e-8; {final prec. was 2.4e-10}
minstep = 1.0e-8; {min. stepsize was 1.0e-10}
z = 6.0; {Nuclear Charge 6 }
redc = 0.3457; {Factor for step reduction }
m = 4; {nbr parameters }
type varrar = array[0..2] of extended;
vary = array[0..3] of extended;
var x,dx,x0,x1,x2 : varrar;
y,dy,y0,y1,y2 : vary;
pardat,par,difference,lenlih : array[1..m] of extended;
nbr,nbry,i,ii,ij,ik,flg1,flg2,n,errc : integer;
count : longint;
k0,k1,etot,eold,vir,stepsize,stepsizey,eexper,ecalb,ecal,ccal,
vnn,vne,vee,ekin,sig0,sig1,kk0,kk1,
sum,sumold,sumnew,r,p,p2,r2,c,wd : extended;
df : text;
{-------------- DEFINE CONSTANTS ------------------------}
procedure init;
begin
ecalb:=-40.457618; {G3} ecal:=-40.2122833; {6-311} ccal:=-37.827717;
nbr:=2; nbry:=3; count:=0;
c:=sqrt(2/3); wd:=sqrt(3);
x0[0]:=0.26677; x0[1]:=1.2426; x0[2]:=1.35585;{start values for vars }
k0 := 1.05518; {start values for parameters }
k1 := 1.21738;
sig0 := 0.24406;
sig1 := 0.3045;
y0[0]:=k0; y0[1]:=k1;
y0[2]:=sig0; y0[3]:=sig1;
pardat[1]:=ecalb-ccal+4.0*0.501003; {atomization in G3}
pardat[2]:=1.0894; {re value; r0 1.094}
{ratios from g3//rhf/6-311+g}
pardat[3]:=-119.9142822/ecal;
pardat[4]:=26.0898758/ecal;
assign(df,'ch4_opt.res');
rewrite(df)
end;
{-------------- INPUT OF VARIABLES ----------------------}
procedure inputy;
var yy : extended;
i : integer;
begin
yy:=0.0;
for i:=0 to nbry do begin
y[i]:=y0[i];
if abs(y0[i])>yy then yy:=abs(y0[i])
end;
stepsizey:=0.001*yy;
end;
procedure input;
var xx : extended;
i : integer;
begin
xx:=0.0;
for i:=0 to nbr do begin
x[i]:=x0[i];
if x0[i]>xx then xx:=x0[i]
end;
stepsize:=0.4*xx; { --- 20% of largest parm. as stepsize}
end;
procedure summing;
var i: integer;
begin
sum:=0.0;
par[1]:=etot-ccal+4.0*0.501003; par[2]:=bohr*p*r;
par[3]:=vne/etot; par[4]:=vee/etot;
for i:=1 to m do begin
difference[i]:=par[i]-pardat[i];
sum:=sum+sqr(difference[i])
end;
{sum:=sum+4*sqr(difference[4])} {e.g. put weight 5 on vee}
end;
{----------- Printout of final results ------------------}
procedure results;
var sum1 : extended;
begin
writeln(df,' g3/6-311 Calc. Difference Values');
for i:=1 to m do begin
write(df,' ',pardat[i]:12:6,par[i]:12:6);
write(df,' ',difference[i]:9:6);
if i=1 then write(df,-par[i]*627.50956:14:1,' kcal/mol');
if i=2 then write(df,par[i]:14:6);
if i in [3,4] then write(df,etot*par[i]:14:6);
writeln(df);
end;
writeln(df,' Sum Squares: ',sum:12:8);
writeln(df,' Parameters determined:');
for i:=0 to nbry do
write(df,y[i]:12:8);
writeln(df);
writeln(df);
writeln(df,' etot ekin vnn vee vne vir');
writeln(df,etot:13:7,ekin:13:7,vnn:13:7,vee:13:7,vne:14:7,vir:12:7);
writeln(df,' P Q R');
writeln(df,x[0]:12:7,x[1]:12:7,x[2]:12:7);
writeln(df)
end;
{----------- Adjust stepsize for every parameter --------}
procedure stepy;
var au,yy : extended;
i,ii : integer;
begin
au:=0.0;
for ii:=0 to nbry do begin
yy:=abs(y[ii]);
if au < yy then au:=yy;
end;
au:=stepsizey/au; { --- AU largest Param.}
for i:=0 to nbry do begin
dy[i]:=au*y[i];
if y[i]=0 then dy[i]:=0.000001;
end
end;
procedure step;
var au,xx : extended;
i,ii : integer;
begin
au:=0;
for ii:=0 to nbr do
begin
xx:=abs(x[ii]);
if au < xx then au:=xx;
end;
au:=stepsize/au; { --- AU largest Param.}
for i:=0 to nbr do
begin
dx[i]:=au*x[i];
if x[i]=0.0 then dx[i]:=0.01;
end;
end;
{----------- FUNCTION SUBROUTINE for Total Energy -------}
function E(u:varrar):extended;
begin
ekin:=2.25*k0/sqr(u[0])+9.0*k1/sqr(u[1]);{ total Ekin }
vne:=-3*z/u[0]; vee:=3*sig0/u[0]; { Z-w0, w0-w0 }
p:=u[2]; p2:=p*p; r:=u[0]+u[1]; r2:=r*r;
vne:=vne-4*(3-sqr((p-1)*(1+u[0]/u[1])))/u[1];
vee:=vee+4*3*sig1/u[1];
vee:=vee+8*2/r;
vne:=vne-8*z/r;
vne:=vne-8/(r*p);
vnn:=4*z/(r*p);
vee:=vee+12/(c*r);
vne:=vne-24*wd/sqrt(2*r2*p + 3*r2*p2 + 3*r2);
vnn:=vnn+3/(c*r*p);
vir :=-(vne+vnn+vee)/ekin; { Virial }
etot:=ekin+vne+vnn+vee;
E:=etot { Total Energy}
end;
{----------- MINIMZATION OF TOTAL ENERGY ----------------}
Procedure Minimize;
Var
Eold,Enew : extended;
ii : integer;
Begin
for ii:=0 to nbr do begin
x1:=x;
x2:=x;
x1[ii]:=x[ii]+dx[ii];
x2[ii]:=x[ii]-dx[ii];
if E(x1) > E(x2) then dx[ii]:=-dx[ii]
end;
Enew:=E(x);
Repeat
Eold:=Enew;
for ii:=0 to nbr do begin
x[ii]:=x[ii]+dx[ii];
Enew:=E(x);
if Enew > Eold then x[ii]:=x[ii]-dx[ii];
end;
Until Enew > Eold
End;
{----------- COMPUTATION OF Emin ------------------------}
Procedure Calculate;
Var
Ebegin,Eend :extended;
i :Integer;
Begin
input;
step;
Repeat
Ebegin:=E(x);
Minimize;
Eend :=E(x);
stepsize:=stepsize*Redc;
step;
Until stepsize<1.0e-10;
End;
function S(v : vary):extended;
begin
k0:=v[0];
k1:=v[1];
sig0:=v[2];
sig1:=v[3];
calculate;
summing;
inc(count);
if count>100000 then begin
close(df);
writeln;
writeln('No convergence in >100000 Iterations');
halt(0)
end;
S:=sum
end;
{----------- MINIMZATION OF TOTAL SUMSQUARES ----------------}
Procedure Minsqr;
Var
Sumold,Sumnew : extended;
ii : integer;
Begin
for ii:=0 to nbry do begin
y1:=y;
y2:=y;
y1[ii]:=y[ii]+dy[ii];
y2[ii]:=y[ii]-dy[ii];
if S(y1) > S(y2) then dy[ii]:=-dy[ii]
end;
Sumnew:=S(y);
Repeat
Sumold:=Sumnew;
for ii:=0 to nbry do begin
y[ii]:=y[ii]+dy[ii];
Sumnew:=S(y);
if Sumnew > Sumold then y[ii]:=y[ii]-dy[ii];
end;
Until Sumnew > Sumold
End;
{----------- COMPUTATION OF Summin ------------------------}
Procedure Calcsum;
Var
Sumbegin,Sumend :extended;
i :Integer;
Begin
y:=y0;
ik:=0;
stepy;
Repeat
Sumbegin:=S(y);
Minsqr;
Sumend :=S(y);
stepsizey:=stepsizey*Redc;
stepy;
inc(ik);
writeln(df);
writeln(df,'Nbr. of function evaluations: ',count:5);
writeln(df,ik:3); results;
count:=0;
writeln(df);
write(ik:6);
until stepsizey<1.0e-7;
writeln(df);
end;
{----------- MAIN PROGRAM -------------------------------}
begin
init;
inputy;
calcsum;
writeln(df,' Final Results of Optimization');
results; writeln;
close(df)
end.
|