Parametrization of CH4 with Gaussian G3

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. Below the Pascal source you can find a C++ source. Compiled with gcc this runs much faster than the Pascal program.
{**** 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.


C++ source for the same code:

I thank John Beveridge, Calgary, for a draft of the Pascal conversion.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/**** CH4 ****    ES 29.10.1998     
gcc ch4.cpp -o ch4 -lm   */

const char tb[] = "             ";
const double bohr = 0.529177249;            /*Angstr./bohr*/
const double hartree = 27.2114;             /*eV/hartree*/
const double kcal = 23.06037;               /*kcal/eV*/
const double crit = 2.4e-10;                /*final prec. was 2.4e-10*/
const double minstep = 1.0e-10;             /*min. stepsize was 1.0e-10*/
const double z = 6.0;                       /*Nuclear Charge 6          */
const double redc = 0.3457;                 /*Factor for step reduction */
const int m = 4;
const int nbr = 3;
const int nbry = 4;
const double ecal =-40.2122833;
const double ccal =-37.827717;  
const double ecalb=-40.457618;              /*G3*/
const double c=0.816496580927726;
const double wd=1.73205080756887729;

/*typedef double varrar[3];  /* pascal[0..2] == c array [3] 
typedef double vary[4];  /* pascal[0..3] == c array [4] */

double x[3];
double dx[3];
double x0[3];
double x1[3];
double x2[3];
double u[3];


double y[4];
double dy[4];
double y00[4];
double y01[4];
double y02[4];
double v[4];

double  pardat[m];       /* pascal array [1..m] == c [m] */
double  par[m];
double  difference1[m];


int count,i,ii,ij,ik,flg1,flg2,n,errc;
double k0,k1,e_result,s_result,etot,Eold,vir,stepsize,stepsizey,eexper,
vnn,vne,vee,ekin,sig0,sig1,kk0,kk1,
sum,sumold,sumnew,r,p,p2,r2,zw;

FILE * df;

/*-------------- DEFINE CONSTANTS ------------------------*/
void init()
{
  count=0;

  x0[0] = 0.26;
  x0[1] = 1.24;
  x0[2] = 1.35;

  k0  = 1.0;
  k1  = 1.0;
  sig0 = 0.40;
  sig1 = 0.40;

  y00[0] = k0;
  y00[1] = k1;
  y00[2] = sig0;
  y00[3] = sig1;

  pardat[0] = ecalb-ccal+4.0*0.501003;
  pardat[1] = 1.0894;    /*re value; r0 1.094*/
  pardat[2] = -119.9142822/ecal;  /*ratios from g3//rhf/6-311+g*/
  pardat[3] = 26.0898758/ecal;
  df = fopen("ch4_opt6.res","w");
  }

/*-------------- INPUT OF VARIABLES ----------------------*/
void inputy()
{
  double yy;
  int i;

  yy = 0.0;
  for( i = 0; i < nbry; i ++) {
    y[i] = y00[i];
    if (abs(y[i]) > yy)  yy = abs(y[i]);
  }
  stepsizey = 0.001*yy;
}

void input1()
{
  double xx;
  int i;

  xx=0.0;
  for( i = 0; i < nbr; i ++) {
    x[i] = x0[i];
    if (x0[i] > xx)  xx = x0[i];
  }
  stepsize=0.4*xx;  /* --- 20% of largest parm. as stepsize*/
}

void summing()
{
  int i;

  sum=0.0;
  par[0] = etot-ccal+4.0*0.501003;
  par[1] = bohr*p*r;
  par[2] = (vne)/etot;
  par[3] = (vee)/etot;
  for( i = 0; i < m; i ++) {
    difference1[i] = par[i]-pardat[i];
    sum = sum+pow(difference1[i],2);
  }
  /*sum:=sum+4*sqr(difference[4])*/  /*put weight 5 on vee*/
}

/*----------- Printout of final results ------------------*/
void results()
{
  double sum1;

  fprintf(df, "%s", "     g3/6-311         Calc.       Difference     Values \n");
  for( i = 0; i < m; i ++) {
    fprintf(df, "%14.6f %14.6f",pardat[i],par[i]);
    fprintf(df, "%14.6f", difference1[i]);
    if (i == 0) { fprintf(df, "%9.1f",-par[i]*627.5096);
                  fprintf(df,"     kcal/mol"); }
    if (i == 1)  fprintf(df, "%14.6f",par[i]);
    if (i == 2 || i == 3)  fprintf(df, "%14.6f",etot*par[i]);
    fprintf(df, "\n");
  }
  fprintf(df, "%s%12.8f\n", "  Sum Squares:       ",sum);
  fprintf(df, "%s\n","  Parameters determined:");
  for( i=0; i < nbry; i ++)
    fprintf(df, "%12.8f",y[i]);
  fprintf(df, "\n");
  fprintf(df, "%s\n","   ekin          vnn           vee           vne            vir");
  fprintf(df, "%13.7f %13.7f %13.7f %14.7f %12.7f\n",ekin,vnn,vee,vne,vir);
  fprintf(df, "%s","    P            Q            R\n");
  fprintf(df, "%12.7f %12.7f %12.7f\n",x[0],x[1],x[2]);
}

/*----------- Adjust stepsize for every parameter --------*/
void stepy()
{
    double au,yy;
    int i,ii;

  au=0.0;
  for( ii = 0; ii < nbry; ii ++) {
    yy = abs(y[ii]);
    if (au < yy)  au = yy;
  }
  au = stepsizey/au;         /* --- AU largest Param.*/
  for( i = 0; i < nbry; i ++) {
    dy[i] = au*y[i];
    if (y[i] == 0)  dy[i] = 0.000001;
  }
}

void step()
{
    double au,xx;
    int i,ii;

  au = 0.0;
  for( ii = 0; ii < nbr; ii ++)
    {
      xx = abs(x[ii]);
      if (au < xx)  au = xx;
    }
    au = stepsize/au;         /* --- AU largest Param.*/
  for( i = 0; i < nbr; i ++)
    {
      dx[i] = au*x[i];
      if (x[i] == 0.0)  dx[i] = 0.01;
    }
}


/*----------- FUNCTION SUBROUTINE for Total Energy -------*/
double E(double u[])
 
{
  ekin = 2.25*k0/pow(u[0],2)+9.0*k1/pow(u[1],2);         /* 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;
  zw = (p-1)*(1+(u[0]/u[1]));
  vne = vne-4*(3-pow(zw,2))/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;                        /* Total Energy*/

  return etot;
}

/*----------- MINIMZATION OF TOTAL ENERGY ----------------*/
void minimize()
{
  double Eold,Enew,E1,E2;
  int ii,i;

  for( ii = 0; ii < nbr; ii ++) {
     for( i = 0; i < nbr; i ++) {
       x1[i] = x[i];
       x2[i] = x[i];
     }
    x1[ii] = x[ii]+dx[ii];
    x2[ii] = x[ii]-dx[ii];
    memcpy(u,x1,sizeof x1);
    E1=E(u);
    memcpy(u,x2,sizeof x2);
    E2=E(u);
    if (E1 > E2)  dx[ii] = -dx[ii];
  }
  memcpy(u,x,sizeof x);
  Enew=E(u);
  do {
    Eold = Enew;
    for( ii = 0; ii < nbr; ii ++) {
      x[ii] = x[ii]+dx[ii];
      memcpy(u,x,sizeof x);
      Enew = E(u);
      if (Enew > Eold)  x[ii] = x[ii]-dx[ii];
    }
  } while (Enew < Eold);
}

/*----------- COMPUTATION OF Emin ------------------------*/
void calculate()
{
  double ebegin,eend;
  int i;

  input1();
  step();
  do {
    memcpy(u,x,sizeof x);
    ebegin = E(u);
    minimize();
    memcpy(u,x,sizeof x);
    eend = E(u);
    stepsize = stepsize*redc;
    step();
  } while (stepsize > minstep);
}

double S(double v[])
{
  k0 = v[0];
  k1 = v[1];
  sig0 = v[2];
  sig1 = v[3];
  calculate();
  summing();
  count += 1;
  if (count>1000000)  {
    fclose(df);
    printf ("\nNo convergence in >1000000 Iterations \n");
    exit(0);
  }
  return sum;
}

/*----------- MINIMZATION OF TOTAL SUMSQARES ----------------*/
void minsqr()
{
  double sumold,sumnew,S1,S2;
  int ii,i;

  for( ii = 0; ii < nbry; ii ++) {
    for( i = 0; i < nbry; i ++) {
      y01[i] = y[i];
      y02[i] = y[i];
    }
    y01[ii] = y[ii]+dy[ii];
    y02[ii] = y[ii]-dy[ii];
    memcpy(v,y01,sizeof y01);
    S1 = S(v);
    memcpy(v,y02,sizeof y02);
    S2 = S(v);
    if (S1 > S2)  dy[ii]=-dy[ii];
  }
    memcpy(v,y,sizeof y);
    sumnew = S(v);

  do {
    sumold = sumnew;
    for( ii = 0; ii < nbry; ii ++) {
      y[ii] = y[ii]+dy[ii];
      memcpy(v,y,sizeof y);
      sumnew = S(v);
      if (sumnew > sumold)  y[ii] = y[ii]-dy[ii];
    }
  } while (sumnew < sumold);
}

/*----------- COMPUTATION OF Summin ------------------------*/
void calcsum()
{
  double sumbegin,sumend;
  int i;

  for( i = 0; i < nbry; i ++ ) {
    y[i] = y00[i];
  }

  ik=0;
  stepy();
  do {
    memcpy(v,y,sizeof y);
    sumbegin = S(v); 
    minsqr(); 
    memcpy(v,y,sizeof y);
    sumend  = S(v);
    stepsizey = stepsizey*redc;  
    stepy();
    ik += 1;
    fprintf(df, "\n%s%5d\n", "Nbr. of function evaluations: " ,count);
    fprintf(df, "#%3d\n",ik);
    results();
    printf ("%6d",ik);
    if (ik==10) stepsizey=1000*stepsizey;
  } while (stepsizey > crit);
  fprintf(df, "\n");
}

/*----------- MAIN PROGRAM -------------------------------*/
int main()
{
  init();
  inputy();
  calcsum();
  fprintf(df,"%s\n", "  Final Results of Optimization ");
  results();
  fclose(df);
  return EXIT_SUCCESS;
}