Sample C language program to calculate hydrogen-like atoms.

top 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.)
In this program, first we input the atomic number Z ( 1 = hydrogen atom, 2 = helium ion ) and reduced mass condition ( 0 = reduced mass, 1 = usual electron mass ).
From these values, it outputs theoretical values of total energy and (Bohr) radius of the hydrogen-like atoms.
Then we input the initial x-coordinate r1 (in MM) of electron, and the absolute value of the total energy E (in eV) of atoms.
From these inputted values, this program outputs the y component of electron's velocity after a quarter of its orbit, and last y coordinate, and WN (the number of de Broglie's waves included in one quarter of the orbital).
Here 1 SS = 1 × 10-23 second and 1 MM = 1 × 10-14 meter.
The initial x-coodinate is automatically increased per calculation until +100.
This method is the same as neutral helium program.

E = -13.60569 eV, and r1 = 5291.77 MM ( Bohr radius ) ------ Hydrogen atom (H).
E = -13.59829 eV, and r1 = 5294.65 MM -----------------------Reduced mass (H).
E = -54.42277 eV, and r1 = 2645.89 MM -----------------------Helium ion (He+).
E = -54.41531 eV, and r1 = 2646.25 MM -----------------------Reduced mass (He+).
In these values ( |E| = |total energy| and r1 = radius ), last VY is just zero. Try them.


#include <stdio.h>
#include <math.h>
                            /* hydrogen like atoms (ions) */
int main(void) 
 {
   int i;
   double Z,whi,toenergy,radius,rm,r,E;  
   double vya,vyb,poten,VX,VY,prexx,preyy,WN,ra;
   double xx,yy,vk,preVY,preWN,midWN,leng,wav,ac;
   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 mp=1.67262171e-27;       
   double alph = 6.64465650e-27;    
                  
                            /* input atomic number Z and reduced mass condition */
 
   printf("Atomic number Z ( H = 1, He+ = 2 )? ");  
   scanf("%lf",&Z);

   printf("You use reduced mass (= 0 ) or usual electron mass (= 1 ) ?");  
   scanf("%lf", &whi);

   rm=me;                 
    if (Z==1) { rm=(me*mp)/(me+mp); }       /* hydrogen reduced mass */
    if (Z==2) { rm=(me*alph)/(me+alph);}    /* helium ion reduced mass  */
    if (whi == 1 ) { rm=me; }

                           /* theoretical values of total energy and (Bohr) radius */

    toenergy = -(Z*Z*rm*ele*ele*ele*ele)/(8.0*epsi*epsi*h*h); 
    radius = (epsi*h*h)/(pai*rm*Z*ele*ele);
    toenergy = toenergy * 6.241509e18;       /* J to eV */
    radius = radius * 1.0e14;            /* meter to MM  */
    printf("Total energy: %.5f ", toenergy);
    printf("Radius: %.2f \n", radius);

                                         /* input  r1 and |E| */
    printf("r1 between nucleus and electron 1 (MM)? ");  
    scanf("%lf",&r);

    printf("total energy |E| of hydrogen like atoms (eV) ? ");  
    scanf("%lf", &E);
 
   for (i=1; i < 100 ;i++) {      /* repeat until r1=initial r1+100 */
                                
                                /* poten = potential energy  */
    poten=-(Z*ele*ele)/(4.0*pai*epsi*r);
                             
                             /* vya= total E-potential energy */  
     vya=-(E*1.60217646e-19)-poten*1.0e14; 
    if (vya > 0) {
                               /* vyb=electron initial velocity (m/sec) */ 
     vyb=sqrt((2*vya)/rm); 
     VY=vyb*1.0e-9;           /* change m/sec to MM/SS */
     prexx=r;  VX=0.0; WN=0.0; preyy=0.0;
   
  
  do {
    xx=prexx+VX; yy=preyy+VY;        /* electron 1 position after 1SS */
    preWN=WN ;
    vk=VX*VX+VY*VY;                  
    leng=sqrt(vk)*1.0e-14;      /* moving length (m) for 1 SS */
     wav=h/(rm*sqrt(vk)*1.0e9);  /* de Broglie wavelength (m)  */ 
    WN=WN+leng/wav;                  /* add de Broglie wavelength */      
                                 /* calculation of VX,VY from Coulomb force  */
    ra=sqrt(prexx*prexx+preyy*preyy);  /* between nucleus and electron  */   
                                   
    ra=ra*1.0e-14;    /* change MM to meter  */
    prexx=prexx*1.0e-14; preyy=preyy*1.0e-14;
    ac=(ele*ele)/(4.0*pai*epsi*rm);
                                    /* acceleration (MM/SS^2) */
    VX=VX+1.0e-32*ac*prexx*(-Z/(ra*ra*ra));   
    VY=VY+1.0e-32*ac*preyy*(-Z/(ra*ra*ra));
    prexx=xx;preyy=yy;
  
   } while (xx >= 0);              /* electron has moved one quater of an orbit? */ 
   if (VY > -0.0001 && VY < 0.0001) {    /* last VY condition */           
  
  printf("r1= %.2f ", r );
  printf("VX= %.6f ", VX);
  printf("VY= %.6f ", VY);
  printf("last y = %.2f ", yy);
  midWN=(preWN+WN)/2.0; printf("midWN= %.6f\n", midWN);
    }
   }  r=r+1;
   }  return 0;
   }