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

Top page (correct Bohr model including helium).
Back to hydrogen molecule ion 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.
From the inputted values, this program outputs the binding energy (eV) of H2+, the ratio of force F1 to centrifugal force, the ratio of forces ( from electron and n1 ) acting on n2.

F1 is the force acting on electron toward the nucleus 1 (n1).
From the force acting on electron and its velocity, we calculate the number of de Broglie wave in one orbit around n1.
( Here Virial theorem E = -T = 1/2 V is used. )
Furthermore, considering precession around x axis by the force of nucleus 2, this de Broglie wave is computed.
Total wave means the sum of them.
We also calculate the de Broglie waves influenced by the force in the x direction by nucleus 2 (= new total waves ).


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

int main(void) 
 {
   int i,la,lb;
   double nnuc,rr,aa,ke,r,nuc,a, rra,rrb,fynb,fxnb,internuc;
   double ntwoforce,debrogli,wave,tim,vp,debropre,wavepre,cf,fgen,cosb;
   double fnb,fone,goux,gouy,gouxx,gouxy,fonea,diele,fxnbb,gouxxx,foneb,rah;
   double newgoux,newx,newradius,newwav,newvelo,newtim,newprewav,twaves;
   double poten,ppot, kinetic,velo,binding;
   double me=9.1093826e-31;
   double pai=3.141592653589793; 
   double epsi=8.85418781787346e-12;
   double h=6.62606896e-34; 
   double ele=1.60217653e-19;  
     
                                      /* 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  and n1 */
    rrb=sqrt(r*r+(nuc-a)*(nuc-a));  /* rrb=  between electron  and n2 */
 
    poten=ke*(-1.0/rra-1.0/rrb+1.0/nuc);   /* 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((2.0*kinetic)/me);  /* velo = electron's velocity (m/s) */ 
    binding=-ppot*0.5-13.606;       /* binding energy (eV) of H2+ */

    fynb=(ke*r)/(rrb*rrb*rrb);  /* n2's force acting on electron in -y direction */
    fxnb=(ke*(nuc-a))/(rrb*rrb*rrb); /* n2's force acting on electron in x direction */
    internuc=ke/(nuc*nuc);   /* force between two nuclei */ 

    ntwoforce = fxnb/internuc; /* ratio of forces acting on n2 */

   debrogli = h/(me*velo);  /* de Broglie wavelength of electron */
   wave=(2.0*pai*rra)/debrogli;    /* de Broglie's waves around nucleus 1 */

   tim=(2.0*pai*rra)/velo;       /* tim = time period of electron  around n1 */
   vp=sqrt((fynb*r)/me);  /* vp = electron precession veleocity around x axis */ 
  
   debropre = h/(me * vp);   /* de Broglie wavelength of precession */
   wavepre=(vp*tim)/debropre;  /*  de Broglie's wave number of precession */ 

   cf=(me*velo*velo)/rra;    /* cf= centrifugal force around n1 */
   fgen=ke/(rra*rra);  /* fgen = n1's force acting on the electron */
   cosb=(rra*rra+nuc*nuc-rrb*rrb)/(2*rra*nuc);  /* cosine b */

   fnb=ke/(rrb*rrb); 
   fnb=(fnb*(rra-nuc*cosb))/rrb;  /* fnb= n2's force * cos (theta) */
 
   fone=fnb+fgen;  /* total force acting on electron toward n1 (=F1) */ 
   goux=fone*cosb;           /* x component of F1 */ 
   gouy=fone*sqrt(1-cosb*cosb);  /* y component of F1 */

   gouxx=goux-fxnb;         /* subtract n2 force from goux  */ 
   fonea=sqrt(gouxx*gouxx+gouy*gouy); /* changed F1  */

   diele=sqrt((nuc+a)*(nuc+a)+r*r);  /* between n2 and electron at position 2 */
   fxnbb=(ke*(nuc+a))/(diele*diele*diele); /* x component of force from n2 */
   gouxxx=goux+fxnbb;        /* add n2 force to goux */
   foneb=sqrt(gouxxx*gouxxx+gouy*gouy);  /* changed F1 at position 2 */

   fonea=(fonea+foneb)/2.0;  /* changed average F1 */
 
   rah=fone/cf;   /* rah = ratio of F1 to centrifugal force around n1 */
 
   newgoux = sqrt(fonea*fonea-gouy*gouy); 
   newx=(r*newgoux)/gouy;
   newradius=sqrt(r*r+newx*newx);   /* radius influened by n2 force in x direction */

   newvelo= sqrt((fonea*newradius)/me);
   newtim=(2.0*pai*newradius)/newvelo;         /* new time period  */
   newprewav=(vp*newtim*me*vp)/h;  /* new de Broglie's wave of precession */ 
   newwav=(2.0*pai*newradius*me*newvelo)/h;    /* new de Broglie's waves number */

   printf("a:%.1f", aa);
   printf(" binding energy (eV): %.4f ", binding); 
   printf("ratio of F1 to centrifugal: %.4f ", rah);
   printf("ratio of forces acting on n2: %.3f \n", ntwoforce);

   printf("waves around n1: %.3f ", wave);
   printf("+ waves around x axis: %.3f ", wavepre);
   la=(int)(1000.0*wave); lb=(int)(1000.0*wavepre);la=la+lb;
   twaves=la/1000.0;
   printf(" = total waves: %.3f \n ", twaves);

   printf("new waves around n1: %.3f ", newwav);
    printf("+ new precession wave: %.3f ", newprewav);
   la=(int)(1000.0*newwav); lb=(int)(1000.0*newprewav);la=la+lb;
   twaves=la/1000.0;
   printf(" = new total waves: %.3f \n ", twaves);
   printf("                        \n");

   aa=aa+10;}                             
   
     return 0;
   }