Sample JAVA program to calculate H2 molecule (= perpendicular two orbits ).

top page

If you copy and paste the program source code below into a text editor, you can easily compile and run this.
(This class file name is hmol, so save this text editor as "", and compile it.)
In this program, first we input the initial x-coordinates of electron-a (= ea in MM) and electron-b (= eb in MM ), the absolute value of H2 binding energy (in eV = experimental value = 4.746 eV ), distance between two nuclei ( in MM, experimental value 7414 MM ).

From the inputted values, after two electrons move a half orbit (= H2 two orbits perpendicular to each other ), this program outputs final coordinates of two electrons, de Broglie wavelength of one orbit, and average forces acting on each nucleus.
Here 1 SS = 1 × 10-23 second and 1 MM = 1 × 10-14 meter.

import java.util.Scanner;
 class hmol {
 public static void main(String[] args) {

           // H2 molecule two perpendicular orbits

            // input initial x-coordinates of two electrons
 Scanner stdIn=new Scanner(;
 System.out.println("electron 1 (=ea) coordinate (MM) from nucleus-a? ");  
 double ea=stdIn.nextDouble();

System.out.println("electron 2 (=eb) coordinate (MM) from nucleus-b? ");  
 double eb=stdIn.nextDouble();

     // input absolute value of bond energy of H2 molecule
System.out.println("binding energy |E| of H2 (= 4.746 eV = experimental value ) ? ");  
 double E=stdIn.nextDouble();

   // input distance between two nuclei of H2 molecule
 System.out.println("distance (= MM ) between two nuclei of H2 (= 7414 MM = experimental value ) ? ");  
 double nuc=stdIn.nextDouble();

           // physical constants
           // me = electron mass, h = Planck constant
           // bohr = Bohr radius ( meter )
           //  epsi = permittivity, ele = electron charge

 double me=9.1093826e-31; double bohr=5291.0*1.0e-14;
 double pai=3.141592653589793; 
 double epsi=8.85418781787346e-12;
 double h=6.62606896e-34; double ele=1.60217653e-19;  
 double nucc=nuc*1.0e-14;  
     // change length unit from MM to meter
                             // poten = potential energy 
double ac=(ele*ele)/(4.0*pai*epsi); double acc=ac/me;
double ara=ea; double brb=eb; 
double arb,bra,ab,avx,avz,bvx,bvy;

double ax=ea;  // ax = electron-a x-coordinate
double az=0.0; // az = electron-a z-coordinate

double bx=nuc+eb;  // bx = electron-b x-coordinate 
double by=0.0;     // by = electron-b y-coordinate

double fax=0.0; double faz=0.0; 
double fbx=0.0; double fby=0.0; 


   // ab = distance between two electrons

double potena=-ac*(1.0/ea+1.0/(nuc-ea)); 

// potena = initial potential energy between electron-a and two nuclei

double potenb=-ac*(1.0/eb+1.0/(nuc+eb)); 

// potenb = initial potential energy between electron-b and two nuclei

double potenc=ac*(1.0/nuc + 1.0/ab);

// potenc = initial potential energies between two nuclei, two electrons

double pnuc=ac/nuc;  double pele=ac/ab;  

double poten=potena+potenb+potenc; 
 // poten = initial total potential energy
 //vya= kinetic energy(J) = total E-potential energy  

double vya=-((E+13.6057*2.0)*1.60217653e-19)-poten*1.0e14; 
 if (vya > 0) {
 double vyb=Math.sqrt(vya/me); 
  avz=-vyb*1.0e-9; // change m/s to MM/SS
 // avz = initial electron-a velocity in z direction
 // avx = initial electron-a velocity in x direction
// bvy = initial electron-b velocity in y direction

 //  bvx = initial electron-b velocity in x direction
 double wna=0.0;    
 double twna=0.0; 
 double shua=0.0;  
 double nutaa=0.0; double nutab=0.0; 
 double nutba=0.0; double nutbb=0.0;
 double nutna=0.0; 

 double tnutaa=0.0; double tnutab=0.0; 
 double tnutba=0.0; double tnutbb=0.0;
 double tnutna=1.0;

 double kax=0.0; double kaz=0.0; 

 double axx,azz,bxx,byy, leng, wav, vk,nut;
 do {
    axx=ax+avx; azz=az+avz; // electron-a moves for 1SS
    bxx=bx+bvx; byy=by+bvy; // electron-b moves for 1SS

   nutna=nutna+1.0;  // calculation cycle number
     // calculate distance ( meter )
// distance between electron-a-nucleus-a (meter)

// distance between electron-a-nucleus-b

// distance between electron-b-nucleus-a

// distance between electron-b-nucleus-b

// between 2 electrons

  ax=ax*1.0e-14; az=az*1.0e-14; bx=bx*1.0e-14; by=by*1.0e-14;

double ahox=(-ax)/(ara*ara*ara)+(nucc-ax)/(arb*arb*arb)-(bx-ax)/(ab*ab*ab);
// ahox = force acting on electron-a in x direction

double bhox=(-bx)/(bra*bra*bra)+(nucc-bx)/(brb*brb*brb)+(bx-ax)/(ab*ab*ab);
// bhox = force acting on electron-b in x direction

double hehox=(ahox+bhox)*0.5;
 // average x-direction force of 2 electrons

// add acceleration (MM/SS^2) to velocity

double ahoz=(-az)/(ara*ara*ara)+(-az)/(arb*arb*arb)+(az)/(ab*ab*ab);
// ahoz = force acting on electron-a in z direction 

double bhoy=(-by)/(bra*bra*bra)+(-by)/(brb*brb*brb)+(by)/(ab*ab*ab);
// bhoy = force acting on electron-b in y direction 

double hehoz=(ahoz-bhoy)*0.5;
// average force acting on two electrons in z,y directions

// add acceleration in z direction

  bvx=avx;  bvy=-avz;

// force acting on nucleus-a from electron-b 

// force acting on nucleus-b from electron-a

// force acting on nucleus-a from electron-a

// force acting on nucleus-b from electron-a

    leng=vk*1.0e-14;      // moving length (m) for 1 SS
     wav=h/(me*vk*1.0e9);  // de Broglie wavelength (m) 
 // number of de Broglie wave contained in each segment
  if (shua < 1.0 ) {

  if ( az > 0 ) { // electron-a moves half orbit
  shua = 2.0; 

// average force acting on each nucleus
 tnutaa=nutaa/nutna; tnutab=nutab/nutna;
 tnutba=nutba/nutna; tnutbb=nutbb/nutna;

  fax=axx; faz=azz; // final coordinate of electron-a
  fbx=bxx-nuc; fby=byy; // final coordinate of electron-b

  twna=wna;   // a half orbit de Broglie wavelength

 ax=axx; az=azz; bx=bxx; by=byy;

   } while (shua < 1.0 );  // moved half orbit           

 System.out.printf("final-electron-a-x:%.6f\n", fax);
 System.out.printf("final-electron-a-z:%.6f\n", faz);

 System.out.printf("final-electron-b-x:%.6f\n", fbx);
 System.out.printf("final-electron-b-y:%.6f\n", fby);

  twna = twna *2.0;

  System.out.printf("one orbit de Broglie wavelength:%.6f\n", twna);

 double force=1.0/(bohr*bohr); double fnuc=1.0/(nucc*nucc*force);

double forcer = (tnutaa+tnutba)/(force)-fnuc;

 System.out.printf("force on nucleus-a:%.6f\n", forcer);

forcer = -(tnutab+tnutbb)/(force)-fnuc;

 System.out.printf("force on nucleus-b:%.6f\n", forcer);