{ LiH, BeH2, BH3, CH4, NH3, H2O, HF  Latest Version 18.12.98
 Calculation of Energy and R1, R2, R3, p, ALPHA, BETA (where applicable);
 Model of G.E.Kimball, parametrized with G3 by E. Schumacher }

uses crt,dos;
const tb='      ';
      Pi  = 3.1415926536;
      bohr=0.529177246;
      hartree=27.2114;
      hatom = 13.633093;
      kcal_ev=23.06053713;
      qs='     様様様様様様様様様様様様様様様様様様様様様様様様様様';
   { Kimball's Parameters }
      c1 : array[0..5,0..6] of extended =
      ((1,1,1,1,1,1,1),
        (1,1,1,1,1,1,1),
        (1,1,1,1,1,1,1),
        (0.4,0.4,0.4,0.4,0.4,0.4,0.4),
        (0.4,0.4,0.4,0.4,0.4,0.4,0.4),
        (0.4,0.4,0.4,0.4,0.4,0.4,0.4));
   {Schumacher's Parameters for G3//6-311g, atomization energy good, structure less}
      c2 : array[0..5,0..6] of extended =
       ((1.01164690,1.05363476,1.06087610,1.04605075,0.95894892,0.92438558,0.91058613),
        (1.28071981,1.13416888,1.15642547,1.20889711,1.34927913,1.42313501,1.41656890),
        (0.0,       0.0,       0.0,       0.0,       1.23272312,1.32139847,1.21497466),
        (0.28512363,0.24543304,0.23296166,0.25273353,0.42320792,0.43732345,0.51017912),
        (0.28991060,0.32013284,0.32783318,0.32960771,0.29390809,0.31779336,0.44166435),
        (0.0,       0.0,       0.0,       0.0,       0.39204779,0.44741336,0.36888374));
        {LiHg3      BeH2g3     BH3g3      CH4g3      NH3g3      H2Og3      HFg3}
type  varrar = array[0..6] of extended;
var   x,dx,x0,x1,x2                    : varrar;
      results                          : array[0..6,1..7] of extended;
      k                                : array[0..5,0..6] of extended;
      nbr,i,ii,ij,j,flg1,flg2,m,n,zz,count,wy,code: integer;
      k0,k1,k2,z,sig0,sig1,sig2,Redc,red,etot,eold,
      vir,stepsize,eexper,c,wd,ekin,e,q2,r2,
      p,p2,qr,acos,bcos,gcos,asin,bsin,asin2,
      sin2,vnn,vne,vee,pr,alfa,beta,gamma,
      q,r,qrp,minstep,be               : extended;
      ch,ch1                           : char;

procedure readdat;
begin
  case zz of
  3: begin
       k0:=k[0,0]*9/4;
       k1:=k[1,0]*9/4;
       sig0:=k[3,0];
       sig1:=k[4,0];
       nbr:=2
     end;
  4: begin
       k0:=k[0,1]*9/4;
       k1:=k[1,1]*9/2;
       sig0:=k[3,1];
       sig1:=k[4,1];
       nbr:=2
     end;
  5: begin
       k0:=k[0,2]*9/4;
       k1:=k[1,2]*3*9/4;
       sig0:=k[3,2];
       sig1:=k[4,2]; wd:=sqrt(3.0);
       nbr:=2
     end;
  6: begin
       k0:=k[0,3]*9/4;           {2*9/8 inner cloud  kin.E.     }
       k1:=k[1,3]*9;             {8*9/8 outer clouds with Proton}
       sig0:=k[3,3];             {Screening for 1s  0.306       }
       sig1:=k[4,3];             {Screening for 2sp(C-H) 0.3830 }
       c:= sqrt(2/3); wd:=sqrt(3.0);
       nbr:=2;
     end;
  7: begin
       k0:=k[0,4]*9/4;           {2*9/8     N(1s) cloud  kin.E. }
       k1:=k[1,4]*3*9/4;         {3*2*9/8   3 N-H clouds        }
       k2:=k[2,4]*9/4;           {2*9/8     N-ep cloud          }
       sig0:=k[3,4];             {Screening N(1s)               }
       sig1:=k[4,4];             {Screening for N-H             }
       sig2:=k[5,4];             {Screening for N-ep            }
       nbr:=4;
     end;
  8: begin
       k0:=k[0,5]*9/4;           {2*9/8     O(1s) cloud  kin.E. }
       k1:=k[1,5]*9/2;           {2*2*9/8   O-H clouds          }
       k2:=k[2,5]*9/2;           {2*2*9/8   O-ep cloud          }
       sig0:=k[3,5];             {Screening O(1s)               }
       sig1:=k[4,5];             {Screening for O-H             }
       sig2:=k[5,5];             {Screening for O-ep            }
       nbr:=5;
     end;
  9: begin
       k0:=k[0,6]*9/4;           {2*9/8     F(1s) cloud  kin.E. }
       k1:=k[1,6]*9/4;           {2*9/8     F-H cloud           }
       k2:=k[2,6]*3*9/4;         {3*2*9/8   F-ep clouds         }
       sig0:=k[3,6];             {Screening F(1s)               }
       sig1:=k[4,6];             {Screening for F-H             }
       sig2:=k[5,6];             {Screening for F-ep            }
       nbr:=4;
     end;
  end; {case}
  Redc :=0.42;                    {Factor for step reduction     }
  minstep:=1.0e-10;
end;

{-------------- INPUT OF VARIABLES ----------------------}
procedure input;
var xx : extended;
    i  : integer;
begin
  clrscr;
  readdat;
  x0[0]:=0.2; x0[1]:=1.0; x0[2]:=1.2;
  x0[3]:=1.5; x0[4]:=1.4; x0[5]:=1.0;
  xx:=0;
  for i:=0 to nbr do
    begin
      x[i]:=x0[i];
      if x0[i]>xx then xx:=x0[i]
    end;
  stepsize:=0.14*xx; { --- 14% of largest parm. as stepsize}
end;

{-------------- Printout of Results ---------------------}
procedure printout;
var i  : integer;
begin
  writeln('   ',etot:16,'      ',2*ekin+e:16,'       ',count:3);
  writeln; writeln(tb,tb,tb,'Values of Variables at Emin');
  for i:=0 to nbr do
    if i mod 2 =1 then begin
      gotoxy(35,wherey);
      writeln(tb,'x[',I,'] =',x[i]:10:7)
    end
    else write(tb,'x[',I,'] =',x[i]:10:7);
  if i mod 2 =0 then writeln;
  writeln;write(tb,'Etot   ',etot:14:7);
  writeln('     -V/T   ',-e/ekin:14:7);
  writeln(tb,'Vnn    ',vnn:14:7,'     Eel    ',vne+ekin+vee:14:7);
  writeln(tb,'T      ',ekin:14:7,'     V      ',vnn+vne+vee:14:7);
  writeln(tb,'Vel    ',vne+vee:14:7,'     Vee    ',vee:14:7);
  writeln(tb,'Vne    ',vne:14:7);
  writeln;
  case zz of
  3: begin
       writeln(tb,'                 LiH - Molecule');
       writeln(qs);
       write(tb,'Li-H distance       ',bohr*(x[0]+x[1]+x[2]):12:4);
       writeln(' Angstrom (calc.)');
       write(tb,'Total Energy        ',hartree*etot:12:4);
       writeln(' eV       (calc.)');
       write(tb,'Energy of Li-Atom   ',-203.138050:12:4);
       writeln(' eV       (G3)');
       write(tb,'Energy of 1 H-Atom  ',-hatom:12:4);
       writeln(' eV       (G3)');
       eexper:=-203.138050-hatom;
       write(tb,'LiH-Bond Energy     ');
       be:=-(hartree*etot-eexper)*kcal_ev;
       write(be:9:1);
       writeln(' kcal/mol    (calc.)');
       writeln(qs);
       wy:=wherey;
       for i:=1 to 5 do begin
         gotoxy(5,wy-7+i);
         write('�');
         gotoxy(58,wy-7+i);
         write('�');
       end;
       gotoxy(5,wy-7);  write('�');
       gotoxy(58,wy-7); write('�');
       gotoxy(5,wy-1);  write('�');
       gotoxy(58,wy-1); write('�');
       results[0,1]:=bohr*x[0]; results[1,1]:=bohr*x[1]; results[2,1]:=0.0;
       results[3,1]:=bohr*(x[0]+x[1]+x[2]); results[4,1]:=0.0;
       results[5,1]:=0.0;  results[6,1]:=be;
     end;
  4: begin
       writeln(tb,'                BeH2 - Molecule');
       writeln(qs);
       write(tb,'Be-H distance       ',bohr*p:12:4);
       writeln(' Angstrom (calc.)');
       write(tb,'Total Energy        ',hartree*etot:12:4);
       writeln(' eV       (calc.)');
       write(tb,'Energy of Be-Atom   ',-398.9145456:12:4);
       writeln(' eV       (G3)');
       write(tb,'Energy of 2 H-Atoms ',-2*hatom:12:4);
       writeln(' eV       (G3)');
       eexper:=-398.9145456-2*hatom;
       write(tb,'BeH2-Atomiz. Energy ');
       be:=-(hartree*etot-eexper)*kcal_ev;
       write(be:9:1);
       writeln(' kcal/mol    (calc.)');
       writeln(qs);
       wy:=wherey;
       for i:=1 to 5 do begin
         gotoxy(5,wy-7+i);
         write('�');
         gotoxy(58,wy-7+i);
         write('�');
       end;
       gotoxy(5,wy-7);  write('�');
       gotoxy(58,wy-7); write('�');
       gotoxy(5,wy-1);  write('�');
       gotoxy(58,wy-1); write('�');
       results[0,2]:=bohr*x[0]; results[1,2]:=bohr*x[1]; results[2,2]:=0.0;
       results[3,2]:=bohr*p; results[4,2]:=180.;
       results[5,2]:=0.0;  results[6,2]:=be/2;
     end;
  5: begin
       writeln(tb,'                 BH3 - Molecule');
       writeln(qs);
       write(tb,'Angle HBH           ',120.0000:12:4);
       writeln(' degrees  (calc.)');
       write(tb,'B-H distance        ',bohr*p:12:4);
       writeln(' Angstrom (calc.)');
       write(tb,'Total Energy        ',hartree*etot:12:4);
       writeln(' eV       (calc.)');
       write(tb,'Energy of B-Atom    ',-670.5637034:12:4);
       writeln(' eV       (G3)');
       write(tb,'Energy of 3 H-Atoms ',-3*hatom:12:4);
       writeln(' eV       (G3)');
       eexper:=-670.5637034-3*hatom;
       write(tb,'BH3-Atomiz. Energy  ');
       be:=-(hartree*etot-eexper)*kcal_ev;
       write(be:9:1);
       writeln(' kcal/mol    (calc.)');
       writeln(qs);
       wy:=wherey;
       for i:=1 to 6 do begin
         gotoxy(5,wy-8+i);
         write('�');
         gotoxy(58,wy-8+i);
         write('�');
       end;
       gotoxy(5,wy-8);  write('�');
       gotoxy(58,wy-8); write('�');
       gotoxy(5,wy-1);  write('�');
       gotoxy(58,wy-1); write('�');
       results[0,3]:=bohr*x[0]; results[1,3]:=bohr*x[1]; results[2,3]:=0.0;
       results[3,3]:=bohr*p; results[4,3]:=120.;
       results[5,3]:=0.0;  results[6,3]:=be/3;
     end;
  6: begin
       writeln(tb,'                 CH4 - Molecule');
       writeln(qs);
       write(tb,'Angle HCH           ',109.4712:12:4);
       writeln(' degrees   (given)');
       write(tb,'C-H distance        ',bohr*p*r:12:4);
       writeln(' Angstr�m  (calc.)');
       write(tb,'Total Energy        ',hartree*etot:12:4);
       writeln(' eV        (calc.)');
       write(tb,'Energy of C-Atom    ',-1029.352704:12:4);
       writeln(' eV        (G3)');
       write(tb,'Energy of 4 H-Atoms ',-4*hatom:12:4);
       writeln(' eV        (G3)');
       eexper:=-1029.352704-4*hatom;
       write(tb,'CH4-Atomiz.  Energy ');
       be:=-(hartree*etot-eexper)*kcal_ev;
       write(be:9:1);
       writeln(' kcal/mol     (calc.)');
       writeln(qs);
       wy:=wherey;
       for i:=1 to 6 do begin
         gotoxy(5,wy-8+i);
         write('�');
         gotoxy(58,wy-8+i);
         write('�');
       end;
       gotoxy(5,wy-8);  write('�');
       gotoxy(58,wy-8); write('�');
       gotoxy(5,wy-1);  write('�');
       gotoxy(58,wy-1); write('�');
       results[0,4]:=bohr*x[0]; results[1,4]:=bohr*x[1]; results[2,4]:=0.0;
       results[3,4]:=bohr*p*r; results[4,4]:=109.4712;
       results[5,4]:=0.0;  results[6,4]:=be/4;
     end;
  7: begin
       writeln(tb,'                 NH3 - Molecule');
       writeln(qs);
       {alfa in the calculations is 180-angle HNL, gamma = angle HNH here}
       gcos:=1-1.5*asin*asin;          {cos gamma in relation to sin alfa}
       gamma:=180.0+180.0/pi*arctan(sqrt(1-gcos*gcos)/gcos); {gamma > 90�!}
       write(tb,'Angle  HNH          ',gamma:12:4);
       writeln(' degrees   (calc.)');
       write(tb,'N-H distance        ',bohr*p*r:12:4);
       writeln(' Angstr�m  (calc.)');
       write(tb,'Total Energy        ',hartree*etot:12:4);
       writeln(' eV        (calc.)');
       write(tb,'Energy of    N-Atom  ',-1484.783076:11:4);
       writeln(' eV        (G3)');
       write(tb,'Energy of  3 H-Atoms ',-3*hatom:11:4);
       writeln(' eV        (G3)');
       eexper:=-1484.783076-3*hatom;
       write(tb,'NH3-Atomiz. Energy   ');
       be:=-(hartree*etot-eexper)*kcal_ev;
       write(be:9:1);
       writeln(' kcal/mol    (calc.)');
       writeln(qs);
       wy:=wherey;
       for i:=1 to 6 do begin
         gotoxy(5,wy-8+i);
         write('�');
         gotoxy(58,wy-8+i);
         write('�');
       end;
       gotoxy(5,wy-8);  write('�');
       gotoxy(58,wy-8); write('�');
       gotoxy(5,wy-1);  write('�');
       gotoxy(58,wy-1); write('�');
       results[0,5]:=bohr*x[0]; results[1,5]:=bohr*x[1];
       results[2,5]:=bohr*x[2];
       results[3,5]:=bohr*p*r; results[4,5]:=gamma;
       results[5,5]:=0.0;  results[6,5]:=be/3;
     end;
  8: begin
       writeln(tb,'                 H2O - Molecule');
       writeln(qs);
       alfa:=180*x[4]/Pi;         {1/2 angle HOH; gamma= angle HOL =}
       beta:=180*x[5]/Pi;         {1/2 angle LOL;  gcos= -acos*bcos }
       write(tb,'Angle  HOH          ',2*alfa:12:4);
       writeln(' degrees   (calc.)');
       write(tb,'Angle  LOL          ',2*beta:12:4);
       writeln(' degrees   (calc.)');
       write(tb,'O-H distance        ',bohr*p*r:12:4);
       writeln(' Angstr�m  (calc.)');
       write(tb,'Total Energy        ',hartree*etot:12:4);
       writeln(' eV        (calc.)');
       write(tb,'Energy of    O-Atom  ',-2041.713315:11:4);
       writeln(' eV        (G3)');
       write(tb,'Energy of  2 H-Atoms ',-2*hatom:11:4);
       writeln(' eV        (G3)');
       eexper:=-2041.713315-2*hatom;
       write(tb,'H2O-Atomiz.  Energy  ');
       be:=-(hartree*etot-eexper)*kcal_ev;
       write(be:9:1);
       writeln(' kcal/mol    (calc.)');
       writeln(qs);
       wy:=wherey;
       for i:=1 to 7 do begin
         gotoxy(5,wy-9+i);
         write('�');
         gotoxy(58,wy-9+i);
         write('�');
       end;
       gotoxy(5,wy-9);  write('�');
       gotoxy(58,wy-9); write('�');
       gotoxy(5,wy-1);  write('�');
       gotoxy(58,wy-1); write('�');
       results[0,6]:=bohr*x[0]; results[1,6]:=bohr*x[1];
       results[2,6]:=bohr*x[2];
       results[3,6]:=bohr*p*r; results[4,6]:=2*alfa;
       results[5,6]:=2*beta;  results[6,6]:=be/2;
     end;
  9: begin
       writeln(tb,'                 HF - Molecule');
       writeln(qs);
       {alfa in the calculations is 180-angle HFL, gamma = angle LFL here}
       gcos:=1-1.5*asin*asin;          {cos gamma in relation to sin alfa}
       gamma:=180.0+180.0/pi*arctan(sqrt(1-gcos*gcos)/gcos); {gamma > 90�!}
       write(tb,'Angle  LFL          ',gamma:12:4);
       writeln(' degrees   (calc.)');
       write(tb,'F-H distance        ',bohr*p*r:12:4);
       writeln(' Angstr�m  (calc.)');
       write(tb,'Total Energy        ',hartree*etot:12:4);
       writeln(' eV        (calc.)');
       write(tb,'Energy of    F-Atom  ',-2712.566713:11:4);
       writeln(' eV        (G3)');
       write(tb,'Energy of  1 H-Atom  ',-hatom:11:4);
       writeln(' eV        (G3)');
       eexper:=-2712.566713-hatom;
       write(tb,'HF-Atomiz. Energy    ');
       be:=-(hartree*etot-eexper)*kcal_ev;
       write(be:9:1);
       writeln(' kcal/mol    (calc.)');
       writeln(qs);
       wy:=wherey;
       for i:=1 to 6 do begin
         gotoxy(5,wy-8+i);
         write('�');
         gotoxy(58,wy-8+i);
         write('�');
       end;
       gotoxy(5,wy-8);  write('�');
       gotoxy(58,wy-8); write('�');
       gotoxy(5,wy-1);  write('�');
       gotoxy(58,wy-1); write('�');
       results[0,7]:=bohr*x[0]; results[1,7]:=bohr*x[1];
       results[2,7]:=bohr*x[2];
       results[3,7]:=bohr*p*r; results[4,7]:=0.0;
       results[5,7]:=gamma; results[6,7]:=be;
     end;
  end; {case}
  gotoxy(5,25);
  textcolor(lightred);
  write('Press any key to continue');
  ch:=readkey;
  textcolor(yellow)
end;

{-------------- Adjust stepsizes ------------------------}
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 then dx[i]:=0.01;
    end;
end;

function U(x: varrar):extended;
begin
  case zz of
  3: begin
       ekin:=k0/sqr(x[0])+k1/sqr(x[1]);      { total Ekin }
       vne:=-3.0*z/x[0]-(3.0-sqr(x[2]/x[1]))/x[1];
       r:=x[0]+x[1]; p:=r+x[2];
       vne:=vne-2.0*z/r-2.0/p;
       vee:=3.0*sig0/x[0]+3.0*sig1/x[1]+4.0/r;
       vnn:=z/p;
     end;
  4: begin
       ekin:=k0/sqr(x[0])+k1/sqr(x[1]);
       vne:=-3.0*z/x[0]-2.0*(3.0-sqr(x[2]/x[1]))/x[1];
       r:=x[0]+x[1]; p:=r+x[2];
       vne:=vne-4.0*z/r-4.0/p-4.0/(r+p);
       vee:=3.0*sig0/x[0]+6.0*sig1/x[1]+8.0/r+4.0/(2.0*r);
       vnn:=2.0*z/p+0.5/p;
     end;
 5: begin
       ekin:=k0/sqr(x[0])+k1/sqr(x[1]);      { total Ekin }
       vne:=-3.0*z/x[0]-3.0*(3.0-sqr(x[2]/x[1]))/x[1];
       r:=x[0]+x[1]; p:=r+x[2];
       vne:=vne-6.0*z/r-6.0/p-12.0/sqrt(r*r+p*p+r*p);
       vee:=3.0*sig0/x[0]+9.0*sig1/x[1]+12.0/r+4.0*wd/r;
       vnn:=3.0*z/p+wd/p;
     end;
  6: begin
       ekin:=k0/sqr(x[0])+k1/sqr(x[1]);        { total Ekin }
       vne:=-3*z/x[0];    vee:=3*sig0/x[0];    { Z-w0, w0-w0 }
       p:=x[2]; p2:=p*p; r:=x[0]+x[1]; r2:=r*r;
       vne:=vne-4*(3-sqr((p-1)*(1+x[0]/x[1])))/x[1];
       vee:=vee+4*3*sig1/x[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);
     end;
  7: begin
       ekin:=k0/sqr(x[0])+k1/sqr(x[1])+k2/sqr(x[2]);
       vne:=-3*z/x[0];  vee:=3*sig0/x[0];
       p:=x[3];  p2:=p*p;  asin:=sin(x[4]); sin2:=sqr(asin);
       acos:=sqrt(1-sin2);
       vne:=vne-3*(3-sqr((p-1)*(1+x[0]/x[1])))/x[1];
       vee:=vee+9*sig1/x[1];
       r:=x[0]+x[1];  q:=x[0]+x[2];  q2:=q*q;  r2:=r*r;
       qr:=q*r; qrp:=qr*p;  pr:=p*r;
       vee:=vee+3*sig2/x[2];
       vee:=vee+4/q;
       vne:=vne-2*z/q;
       vee:=vee+12/r;
       vne:=vne-6*z/r;
       vne:=vne-6/pr;
       vnn:=3*z/pr;
       vee:=vee+12/sqrt(2*qr*acos+q2+r2);
       vne:=vne-6/sqrt(2*qrp*acos+r2*p2+q2);
       vee:=vee+12/(sqrt(3)*r*asin);
       vne:=vne-12/sqrt(3*p*sin2+sqr(p-1))/r;
       vnn:=vnn+3/(sqrt(3)*pr*asin);
     end;
  8: begin
       ekin:=k0/sqr(x[0])+k1/sqr(x[1])+k2/sqr(x[2]);
       vne:=-3*z/x[0];    vee:=3*sig0/x[0];
       p:=x[3];  p2:=p*p;
       vne:=vne-2*(3-sqr((p-1)*(1+x[0]/x[1])))/x[1];
       vee:=vee+6*sig1/x[1];
       r:=x[0]+x[1];  q:=x[0]+x[2];  q2:=q*q;  r2:=r*r;  qr:=q*r;
       vee:=vee+6*sig2/x[2];
       asin:=sin(x[4]);   acos:=cos(x[4]); asin2:=asin*asin;
       bsin:=sin(x[5]);   bcos:=cos(x[5]);
       vee:=vee+8/r;
       vne:=vne-4*z/r;
       vne:=vne-4/(r*p);
       vnn:=2*z/(r*p);
       vee:=vee+8/q;
       vne:=vne-4*z/q;
       vee:=vee+2/(r*asin);
       vne:=vne-4/sqrt(-2*r2*p+4*r2*p*asin2+r2*p2+r2);
       vee:=vee+16/sqrt(2*qr*acos*bcos+q2+r2);
       vnn:=vnn+1/sqrt(4*r2*p2*asin2);
       vne:=vne-8/sqrt(2*qr*p*acos*bcos+r2*p2+q2);
       vee:=vee+2/(q*bsin);
     end;
  9: begin
       ekin:=k0/sqr(x[0])+k1/sqr(x[1])+k2/sqr(x[2]);
       vne:=-3*z/x[0];    vee:=3*sig0/x[0];
       p:=x[3];  p2:=p*p; asin:=sin(x[4]); sin2:=sqr(asin);
       acos:=sqrt(1-sin2);
       vne:=vne-(3-sqr((p-1)*(1+x[0]/x[1])))/x[1];
       vee:=vee+3*sig1/x[1];
       r:=x[0]+x[1];  q:=x[0]+x[2];
       q2:=q*q;  r2:=r*r;  qr:=q*r; pr:=p*r;
       vee:=vee+9*sig2/x[2];
       vee:=vee+4/r;
       vne:=vne-2*z/r;
       vne:=vne-2/pr;
       vnn:=z/pr;
       vee:=vee+12/q;
       vne:=vne-6*z/q;
       vee:=vee+12/sqrt(2*qr*acos+q2+r2);
       vne:=vne-6/sqrt(2*qr*p*acos+r2*p2+q2);
       vee:=vee+12/(sqrt(3)*q*asin);
     end;
  end;
  e:=vnn+vne+vee;                {total Coulomb energy epot}
  vir :=ekin+ekin+e;             {Virial = 2*ekin+epot     }
  etot:=ekin+e;                  {total energy = ekin+epot }
  inc(count);
  U:=etot
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 U(x1) > U(x2) then dx[ii]:=-dx[ii]
  end;
  Enew:=U(x);
  Repeat
    Eold:=Enew;
    for ii:=0 to nbr do begin
      x[ii]:=x[ii]+dx[ii];
      Enew:=U(x);
      if Enew > Eold then x[ii]:=x[ii]-dx[ii];
    end;
  Until Enew > Eold
End; {End MINIMZATION}

{----------- COMPUTATION OF Emin ------------------------}
procedure mindi;
var
  Ebegin,Eend :extended;
begin
  input;
  step;
  count:=0;
  repeat
    Ebegin:=U(x);
    Minimize;
    Eend  :=U(x);
    stepsize:=stepsize*Redc;
    step;
  until stepsize<minstep;
  writeln;
  printout
End;

procedure question;
begin
  repeat
    write('         Give number Z or ''A''(ll) ? ');
    ch:=readkey;
    z:=0.0;
    if ord(ch)<58 then begin
      val(ch,z,code);
      zz:=round(z+0.0001); z:=0.00001*(100000*zz); 
      if (z<3.0) or (z>9.0) then begin
        gotoxy(10,wy); write('             ');
        gotoxy(1,wy)
      end
    end;
  until ((z>=3.0) and (z<=9.0)) or (upcase(ch)='A');
    writeln;
    write('         Kimballs or Schumachers parameters, K or S ?');
    ch1:=readkey;
    if upcase(ch1)='K' then begin
      for i:=0 to 5 do
        for j:=0 to 6 do
          k[i,j]:=c1[i,j]
    end;
    if upcase(ch1)='S' then begin
      for i:=0 to 5 do
        for j:=0 to 6 do
          k[i,j]:=c2[i,j]
    end;
end;

procedure summary;
var i,j : integer;
    Lst : Text;
begin
  clrscr;
  assign(Lst,'hydrk_S.txt');
  rewrite(Lst);
  writeln(Lst);
  writeln(Lst,'                  Summary 1 of computed properties:');
  writeln(Lst,'  陳陳陳陳賃陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳�');
  write(Lst,'           �      LiH       BeH2      BH3      CH4');
  writeln(Lst,'     Property');
  writeln(Lst,'  陳陳陳陳津陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳�');
  wy:=wherey-1;
  gotoxy(3,wy+1);  write(Lst,' R(Z)   ��');
  gotoxy(56,wy+1); write(Lst,'Radius of He core');
  gotoxy(3,wy+2);  write(Lst,' R(Z-H) ��');
  gotoxy(56,wy+2); write(Lst,'Radius of Z-H cloud');
  gotoxy(3,wy+3);  write(Lst,' R(Z-L) ��');
  gotoxy(56,wy+3); write(Lst,'Radius of lone pair');
  gotoxy(3,wy+4);  write(Lst,' d(Z-H) ��');
  gotoxy(56,wy+4); write(Lst,'distance Z-H');
  gotoxy(3,wy+5);  write(Lst,' < HZH  ��');
  gotoxy(56,wy+5); write(Lst,'angle H-Z-H');
  gotoxy(3,wy+6);  write(Lst,' < LZL  ��');
  gotoxy(56,wy+6); write(Lst,'angle LoneP-Z-LoneP');
  gotoxy(3,wy+7);  write(Lst,'BE kc/mol�');
  gotoxy(56,wy+7); write(Lst,'Mean Thermod.BondEn.');
  gotoxy(1,wy+8);
  write(Lst,'  陳陳陳陳珍陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳�');
  for i:=0 to 6 do
    for j:=1 to 4 do begin
      gotoxy(13+10*(j-1),wy+i+1);
      if results[i,j]=0.0 then write(Lst,'      -     ')
      else
        if i<6 then write(Lst,results[i,j]:10:4)
        else write(Lst,results[i,j]:7:1,'   ')
    end;
  writeln(Lst);
  gotoxy(5,25);
  textcolor(lightred);
  write('Press any key to continue');
  ch:=readkey;
  textcolor(yellow);
  clrscr;
  writeln(Lst);
  writeln(Lst,'                  Summary 2 of computed properties:');
  writeln(Lst,'  陳陳陳陳賃陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳�');
  write(Lst,'           �      CH4       NH3       H2O       HF');
  writeln(Lst,'     Property');
  writeln(Lst,'  陳陳陳陳津陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳�');
  wy:=wherey-1;
  gotoxy(3,wy+1);  write(Lst,' R(Z)   ��');
  gotoxy(56,wy+1); write(Lst,'Radius of He core');
  gotoxy(3,wy+2);  write(Lst,' R(Z-H) ��');
  gotoxy(56,wy+2); write(Lst,'Radius of Z-H cloud');
  gotoxy(3,wy+3);  write(Lst,' R(Z-L) ��');
  gotoxy(56,wy+3); write(Lst,'Radius of lone pair');
  gotoxy(3,wy+4);  write(Lst,' d(Z-H) ��');
  gotoxy(56,wy+4); write(Lst,'distance Z-H');
  gotoxy(3,wy+5);  write(Lst,' < HZH  ��');
  gotoxy(56,wy+5); write(Lst,'angle H-Z-H');
  gotoxy(3,wy+6);  write(Lst,' < LZL  ��');
  gotoxy(56,wy+6); write(Lst,'angle LoneP-Z-LoneP');
  gotoxy(3,wy+7);  write(Lst,'BE kc/mol�');
  gotoxy(56,wy+7); write(Lst,'Mean Thermod.BondEn.');
  gotoxy(1,wy+8);
  write(Lst,'  陳陳陳陳珍陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳陳�');
  for i:=0 to 6 do
    for j:=4 to 7 do begin
      gotoxy(13+10*(j-4),wy+i+1);
      if results[i,j]=0.0 then write(Lst,'      -     ')
      else
        if i<6 then write(Lst,results[i,j]:10:4)
        else write(Lst,results[i,j]:7:1,'   ')
    end;
  writeln(Lst);
  close(Lst);
  gotoxy(5,25);
  textcolor(lightred);
  write('Press any key to continue');
  ch:=readkey
end;

{main}
begin
  textbackground(blue);
  textcolor(yellow);
  repeat
    clrscr;
    writeln;
    writeln(tb,'浜様様様様様様様様様様様様様様様様様様様様様様様様様様様様様様融');
    writeln(tb,'�  This Programm computes LiH, BeH2, BH3, CH4, NH3, H2O, or HF �');
    writeln(tb,'�  by the method of KIMBALL, parametrized                      �');
    writeln(tb,'�     by Ernst Schumacher, University of Bern, CH-3000 Bern 9  �');
    writeln(tb,'�  You choose among the seven molecules by giving the          �');
    writeln(tb,'�  nuclear charge Z of the central atom (3,4,5,6,7,8,or 9)     �');
    writeln(tb,'�  or type ''A''(ll) to compute the whole set                    �');
    writeln(tb,'藩様様様様様様様様様様様様様様様様様様様様様様様様様様様様様様夕');
    wy:=wherey;
    question;
    if upcase(ch)='A' then begin
      z:=2.0;
      for m:=1 to 7 do begin
        z:=z+1.0;
        zz:=round(z+0.0001);
        mindi
      end;
      summary
    end
    else begin
      mindi;
      gotoxy(1,25);
      write('  Another one of the seven molecules (Y/..) ? ');
      ch:=readkey;
      write(ch)
    end;
  until upcase(ch)<>'Y';
end.