Sample C program of hydrogen molecule - average values.

Top page (correct Bohr model including helium).
Back to molecular bond page.

If you copy and paste the program source code below into a text editor, you can easily compile and run this.
(This program is simple C language, so save this text editor as "filename.c", and compile it.)
Here we use the new unit of 1 MM = 1 × 10-14 meter.

In this program, first we input the internuclear distance (in MM ) of H2.
And then input y-coordinate "r" (in MM) , and x-coordinate "a" ( in MM ) of electron 1.
From the inputted values, this program outputs the binding energy (eV) of H2, total force acting on nucleus 1 (nucforce), x,y components of forces acting on electron 1 (= elefx, elefy ).
These forces are expressed as the ratio to the force between electron and nucleus of H atom's ground state.

From the force acting on electron and its velocity, we calculate the number of de Broglie wave in one orbit.
( Here Virial theorem E = -T = 1/2 V is used. )
"another x" means another x coordinate of electron 1 after moving in the direction of force.
And then we do above calculation again when x-coordinate is another x.
"avebinding" means average binding energies ( when x is "a" and "another x" ), "avewave" means average de Broglie waves ( when x is "a" and "another x" ).
"after-elefx" means the force component in the x direction acting on electron 1 after moving.
The initial x-coodinate is automatically increased per calculation until +100.


#include <stdio.h>
#include <math.h>

void calc(double ca);
void main(void);
  
  double ke,r,nuc, rra,rrb,rrc;
  double poten,ppot, kinetic,velo,binding,elefx,elefy,eleforce;
  double radius, debroglie, wave, nucforce, anotherx;
  double me=9.1093826e-31;
  double pai=3.141592653589793; 
  double epsi=8.85418781787346e-12;
  double h=6.62606896e-34; 
  double ele=1.60217653e-19;  
  double boh=5292.0*5292.0*1.0e-28;  /* boh = ( Bohr radius )^2 */ 
     

void calc(double ca )         /* argument ca = a  */
  {
    rra=sqrt(ca*ca+r*r);       /* rra =distance between electron 1 and n1 */
    rrb=sqrt(r*r+(nuc-ca)*(nuc-ca));  /* rrb=  between electron 1 and n2 */
    rrc=sqrt(4*r*r+(nuc-2*ca)*(nuc-2*ca));  /* between two electrons */
 
    poten=ke*(-2.0/rra-2.0/rrb+1.0/nuc+1.0/rrc);   /* poten=potential energy (J) */
    ppot=poten*6.241509e18;       /* ppot = potential energy (eV) */ 
                                  
    kinetic=-0.5*poten;      /* kinetic = total kinetic energy of two electrons (J) */
    velo=sqrt(kinetic/me);  /* velo = electron's velocity (m/s) */ 
    binding=-ppot*0.5-13.606*2;       /* binding energy (eV) of H2 */

                            /* foeces acting on electron 1 ( x, -y directions ) */
   elefx= -ca/(rra*rra*rra) + (nuc-ca)/(rrb*rrb*rrb) - (nuc-2*ca)/(rrc*rrc*rrc);
   elefy= r/(rra*rra*rra) + r/(rrb*rrb*rrb) - (2*r)/(rrc*rrc*rrc);
   eleforce = ke*sqrt(elefx*elefx + elefy*elefy);  /* total force */

   radius=(me*velo*velo)/eleforce;  /* rotation radius from centrifugal force */
   debroglie = h/(me*velo);       /* de Broglie wavelength of electron */
   wave=(2*pai*radius)/debroglie;  /* wave's number in one orbit */

   nucforce= ca/(rra*rra*rra)-1.0/(nuc*nuc)+(nuc-ca)/(rrb*rrb*rrb);
   nucforce=nucforce * boh;            /* total force acting on nucleus 1 */

   anotherx=ca+(2.0*r*elefx)/elefy;    /* another x coordinate in direction of force  */

  }


void main(void) 
 {
   int i;
   double a,aa,nnuc,rr, wwave, bbind;
     
                                      /* input  nuc, x, y coordinate */
   
    printf(" Internuclear distance (MM) of H2 molecule ? ");  
    scanf("%lf",&nnuc);
 
    printf(" r (MM) ? = y coordinate of electron 1 ");  
    scanf("%lf",&rr);

    printf(" a (MM) ? = x coordinate of electron 1 ");  
    scanf("%lf", &aa);
    printf("                        \n");

    r = rr * 1.0e-14; nuc = nnuc * 1.0e-14;  /*  change MM to meter */ 
    ke = (ele*ele)/(4.0*pai*epsi);
 
   for (i=1; i < 10 ;i++) {      /* repeat until a=initial a+100 */

    a=aa*1.0e-14;
    calc(a);                    /*  to void calc()    */

   elefx=elefx * boh; elefy = elefy * boh; 
   a=anotherx; anotherx=anotherx*1.0e14;     /* a = another x  */
   printf("a:%.1f ", aa);
   printf(" binding: %.4f ", binding); 
   printf(" elefx: %.3f ", elefx);
   printf(" elefy: %.3f ", elefy);
   printf("nucforce: %.3f \n", nucforce);
   printf(" another x : %.2f ", anotherx);
  

    wwave = wave; bbind = binding;
    calc(a);                         /*  repeat calcu, when a = anotherx   */
                                    /*  average waves and binding energy  */
    wave=(wave+wwave)/2.0; binding=(binding+bbind)/2.0;

    elefx=elefx * boh;
    printf("  after-elefx: %.3f ", elefx);   /* elefx at another x */
    printf(" avebinding: %.4f ", binding);
    printf(" avewave: %.4f \n", wave);     
    printf("                        \n");

   aa=aa+10;}                             
     
   }