Simple C program to calculate three-electron atoms (ions).

top page
Back to chemistry 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, we first input atomic number Z = 3, 4, 5, 6, 7, 8, 9, or 10.
Next we input the initial x-coordinate r1 (in MM) of electron 1, and the absolute value of the total energy E (in eV) of three-electron atoms (ions).
From the inputted values, this program outputs the y component of electron 1 velocity after a quarter of its orbit, 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.


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

int main(void) 
 {
   int i;
   double Z,r,E,ab,ac,Eb,Rb,Rbb,Rbc ; 
   double vya,vyb,poten,VX,VY,prexx,preyy,WN,ra,rb,rc;
   double xx,yy,vk,preVY,preWN,midWN,leng,wav,rr,disb;
   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 Z,  r1 and |E| */

   printf("Atomic number Z ? (Li=3, Be+=4, B2+=5, C3+=6, N4+=7...) ");  
   scanf("%lf",&Z);

   printf("r1 between nucleus and electron 1 (MM)? ");  
   scanf("%lf",&r);

   printf("total energy |E| of three-electron-atom (eV) ? ");  
   scanf("%lf", &E);
 
  
   ab=(ele*ele)/(4.0*pai*epsi); ac=ab/me;
 
   Eb=-((Z-2)*(Z-2)*ele*ele*ele*ele*me)/(32*epsi*epsi*h*h); /* Eb=2S energy */
   Rb=(4.0*epsi*h*h)/((Z-2.0)*pai*ele*ele*me);           /* Rb=2S radius */
   Rbb=Rb*1.0e14;                              /* change meter to MM */
   Rbc=Rbb/(sqrt(2.0));
   printf("Rbb= %.2f\n", Rbb);

   for (i=1; i < 100 ;i++) {      /* repeat until r1=initial r1+100 */
    rr=r*1.0e-14;   /* change r(MM) to rr(meter) */

                    /* disb=distance (m) between initial e1S and e2S */ 
     disb=sqrt(rr*rr+Rb*Rb);
 
     poten=-(2.0*Z*ab)/rr+ab/(2.0*rr)+(2.0*ab)/disb-(Z*ab)/(Rb);
                             
                           /* vya= total E - potential energy - 2S kinetic energy(-Eb) */  
     vya=-(E*1.60217646e-19)-poten+Eb; 

  if (vya > 0) {
                               /* vyb=electron initial velocity (m/sec) */ 
   vyb=sqrt(vya/me); 
   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 */
    preVY=VY;preWN=WN ;
    vk=VX*VX+VY*VY;                  
    leng=sqrt(vk)*1.0e-14;      /* moving length (m) for 1 SS */
     wav=h/(me*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 1S electron  */   
    rb=sqrt(4.0*prexx*prexx+2.0*preyy*preyy); /* between two 1S electrons  */
                                   /* rc=distance between 1S ele and 2S ele */
    rc=sqrt(prexx*prexx+(Rbb*Rbb)/2.0+(preyy+Rbc)*(preyy+Rbc));

                                   
    ra=ra*1.0e-14; rb=rb*1.0e-14;  rc=rc*1.0e-14;  /* change MM to meter  */
    prexx=prexx*1.0e-14; preyy=preyy*1.0e-14;
    
                                    /* acceleration (MM/SS^2) */
    VX=VX+1.0e-32*ac*prexx*(-Z/(ra*ra*ra)+2.0/(rb*rb*rb)+1.0/(rc*rc*rc));   
    VY=VY+1.0e-32*ac*((-Z*preyy)/(ra*ra*ra)+preyy/(rb*rb*rb)+(preyy+Rbc*1.0e-14)/(rc*rc*rc));

    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("preVY= %.6f ", preVY);
  midWN=(preWN+WN)/2.0; printf("midWN= %.6f\n", midWN);
    }
   }  r=r+1;
   }  return 0;
   }