Sources and Compilation of Kimball problems with C++

In preparation!

As example, the following optimization of the parameters for CH4.
Compiled with gcc-4.6.3 on Raspberry Pi, runnable under Debian Linux ch4t7.x

/* Parameterization of CH4 with four results of ab initio
   computation of G3 (Gaussian) and HF/6-311G(dp) */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/**** CH4 ****  after pascal source  ES 29.10.1998; with     
gcc ch4t6.cpp -o ch4t6.x -lm   */

const double bohr = 0.529177249;         //Angstr./bohr
const double crit = 2.4e-15;             //final prec. was 2.4e-10
const double minstep = 1.0e-15;          //min. stepsize was 1.0e-10
const double redc = 0.3457;              //Factor for step reduction

/* constants of the specific example, energies in hartree */
const int m = 4;
const int nbr = 3;                       // # of CH4 vars: R1,R2,p
const int nbry = 4;                      // # of parameters
const double z = 6.0;                    // Nuclear Charge 6
const double ecal =-40.2122833;          // Etot 6-311G(dp)
const double ccal =-37.827717;           // C-Atom exp energy ground state 
const double ecalb=-40.457618;           // Etot G3

const double c =0.816496580927726;       // Sqrt(2/3)
const double wd=1.732050807568877;       // Sqrt(3)

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];
double  par[m];
double  difference[m];


int count,i,ii,ij,ik;
double k0,k1,etot,Eold,vir,stepsize,stepsizey,
vnn,vne,vee,ekin,sig0,sig1,sum,r,p,zw;

FILE * df;

/*-------------- DEFINE CONSTANTS ------------------------*/
void init()
{
  count=0;
          
  x0[0] = 0.2647631;
  x0[1] = 1.2527386;
  x0[2] = 0.5424833;
 
  k0  = 1.04310769;
  k1  = 1.21655927;
  sig0 = 0.26308009;
  sig1 = 0.32263502;

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

  pardat[0] = ecalb-ccal+4.0*0.501003;  // atomization energy
  pardat[1] = 1.0894;                   // re value; r0 1.094
  pardat[2] = -119.9142822/ecal;        // ratios g3//rhf/6-311+g
  pardat[3] = 26.0898758/ecal;
  df = fopen("ch4t7opt.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.3*xx;      // --- 30% 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 ++) {
    difference[i] = par[i]-pardat[i];
    sum = sum+pow(difference[i],2);
  }
  sum = sum+7.0*pow(difference[2],2);  // put weight 8 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", difference[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;
  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);
  vne = -3.0*z/u[0];
  vee = (3.0*sig0)/u[0];
  p = u[2]; r = u[0]+u[1];
  zw = p/u[1];
  vne = vne-4.0*(3.0-pow(zw,2))/u[1];
  vee = vee+4.0*3.0*sig1/u[1];
  vee = vee+8.0*2.0/r;
  vne = vne-8.0*z/r;
  vne = vne-8.0/(r+p);
  vnn =     4.0*z/(r+p);
  vee = vee+12.0/(c*r);
  vne = vne-24.0/sqrt(r*r + (r+p)*(r+p) + (2.0/3.0)*r*(r+p));
  vnn = vnn+3.0/(c*(r+p));
  vir = -(vne+vnn+vee)/ekin;
  etot = ekin+vne+vnn+vee;

  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);
    count=0;
    fprintf(df, "#%3d\n",ik);
    results();
    printf ("%6d\n",ik);
    if (ik==7) stepsizey=24*stepsizey;
    if (ik==12) stepsizey=24*stepsizey;
    if (ik==17) stepsizey=24*stepsizey;
    if (ik==22) stepsizey=24*stepsizey;
    if (ik==27) stepsizey=24*stepsizey;
    if (ik==32) stepsizey=24*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;
}