Python program to calculate two-electron 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 Python ( 2.7 version ), so save this text editor as "filename.py" and compile it.)
In this program, first we input the initial x-coordinate r1 (in MM) of electron 1, and the absolute value of the total energy E (in eV) of helium atom.
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-22 second and 1 MM = 1 × 10-14 meter. Velocity 1 MM/SS = 1 × 108 m/s.
Python is a little slower than C and JAVA, so here we change the definition of 1 SS.
Orbits are a little rough, and calculation value becomes a little different from C and JAVA. ( Though almost same. )
So it's better to run C and JAVA programs, if you can.
( In this program, inputting energy |E| = 79.0054 and r1 = 3074.0 give WN = 0.250000. )
The initial x-coodinate is automatically increased per calculation until +50.


import math
                         # input initial r1 and |E|
r=float(raw_input("r1 between nucleus and electron 1 (MM) ? ")) 
E=float(raw_input("total energy |E| of helium atom (eV) ? "))
me=9.1093826e-31;
pai=3.141592653589793; epsi=8.85418781787346e-12
h=6.62606896e-34; ele=1.60217653e-19
                           
Z = 2.0;                  # helium atomic number
nucle = 6.64465650e-27;   # alpha particle
                          # rm = reduced mass
rm=(2.0*me*nucle)/(2.0*me+nucle); rm=rm*0.5
                        
                          # repeat until r1 = r1 + 50      
for num in range(1,50):
                          # initial potential energy
    poten=-(2.0*Z*ele*ele)/(4.0*pai*epsi*r)+(ele*ele)/(4.0*pai*epsi*2.0*r)
                          # vya = total kinetic energy
    vya=-(E*1.60217646e-19)-poten*1.0e14 
    if vya > 0:
        vyb=math.sqrt(vya/me);  # vyb = initial velocity (m/s) 
        VY=vyb*1.0e-8;          # change m/s to MM/SS = initial VY
        prexx=r; xx = r         # initial x coordinate = r
        preyy=0.0               # initial y coordinate = 0
                                # WN = de Broglie wave number                             
        VX=0.0; WN=0.0          # initial VX ( x velocity ) = 0

        while xx > 0:           # repeat until  move 1/4 orbit
            xx=prexx+VX; yy=preyy+VY   # coordinate after 1 SS      
            preVY=VY
            vk=VX*VX+VY*VY                  
            leng=math.sqrt(vk)*1.0e-14      # moving distance (m) for 1 SS
            wav=h/(rm*math.sqrt(vk)*1.0e8)   # de Broglie wavelength (m)  
            WN=WN+leng/wav                 #  add de Broglie wave number for 1 SS
 
            ra=math.sqrt(prexx*prexx+preyy*preyy) # distance nucleus and electron     
            rb=math.sqrt(4.0*prexx*prexx+2.0*preyy*preyy)  # between two electrons
            ra=ra*1.0e-14; rb=rb*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)
                                             # add acceleration ( MM/SS^2 )
            VX=VX+1.0e-30*ac*prexx*(-Z/(ra*ra*ra)+2.0/(rb*rb*rb))   
            VY=VY+1.0e-30*ac*preyy*(-Z/(ra*ra*ra)+1.0/(rb*rb*rb))
            prexx=xx;preyy=yy
        if VY > -0.0001 and VY < 0.0001:    # last VY condition       
            print "r1: "+str(r)+"   ",
            print("VX:%.6f " % VX),
            print("VY:%.6f " % VY),
            print("preVY:%.6f " % preVY),
            print("WN:%.6f " % WN)
    r=r+1