Sample C language program to calculate hydrogen molecule (H2).

Top page (correct Bohr model including helium).
Back to hydrogen molecule 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 moving in the direction of force.
The initial x-coodinate is automatically increased per calculation until +100.


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

int main(void) 
 {
   int i;
   double nnuc,rr,aa,ke,r,nuc,a, 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 */ 
     
                                      /* 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;
    rra=sqrt(a*a+r*r);       /* rra =distance between electron 1 and n1 */
    rrb=sqrt(r*r+(nuc-a)*(nuc-a));  /* rrb=  between electron 1 and n2 */
    rrc=sqrt(4*r*r+(nuc-2*a)*(nuc-2*a));  /* between two electrons */
 
    poten=ke*(-2.0/rra-2.0/rrb+1.0/nuc+1.0/rrc);   /* poten=potential energy (J) */
    ppot=poten*6.241509*1.0e18;       /* 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= -a/(rra*rra*rra) + (nuc-a)/(rrb*rrb*rrb) - (nuc-2*a)/(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= a/(rra*rra*rra)-1.0/(nuc*nuc)+(nuc-a)/(rrb*rrb*rrb);
   nucforce=nucforce * boh;            /* total force acting on nucleus 1 */

   anotherx=a+(2.0*r*elefx)/elefy; 
   anotherx = anotherx * 1.0e14;      /* another x coordinate of electron 1 */

   elefx=elefx * boh; elefy = elefy * boh;
 
   printf("a:%.1f ", aa);
   printf(" binding-energy: %.4f ", binding); 
   printf(" elefx: %.3f ", elefx);
   printf(" elefy: %.3f \n", elefy);
 
   printf("nucforce: %.3f ", nucforce);
   printf(" another x : %.2f ", anotherx);
   printf(" de Broglie wave: %.3f \n", wave);
   printf("                        \n");

   aa=aa+10;}                             
   
     return 0;
   }