Sample JAVA program to calculate H2 molecule.

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 "hmol.java", 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, 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.in);
 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 = experiment ) ? ");  
 double E=stdIn.nextDouble();

   // input distance between two nuclei of H2 molecule
 System.out.println("distance (= MM ) between two nuclei of H2 (= 7414 MM = experiment ) ? ");  
 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=Math.sqrt((ax-bx)*(ax-bx)+by*by+az*az);

   // 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

// forcei = force attracting electron-a towards nucleus-a (ratio/1/bohr^2)

double forcei=bohr*bohr*1.0e28*(1.0/(ea*ea)-1.0/((nuc-ea)*(nuc-ea))+1.0/((eb+nuc-ea)*(eb+nuc-ea)));

System.out.print("\n");
System.out.printf("electron1-attraction toward nucleus-a:%.6f\n", forcei);

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 potential energy

                             
                    //vya= kinetic energy = total E-potential energy  

double vya=-((E+13.6*2.0)*1.60217646e-19)-poten*1.0e14; 
 if (vya > 0) {

                               
 double vyb=Math.sqrt(vya/me); 
  avz=-vyb*1.0e-9;  // avz = initial electron-a velocity in z direction
  avx=0.0;      // avx = initial electron-a velocity in x direction
 
  bvy=vyb*1.0e-9; // bvy = initial electron-b velocity in y direction
  bvx=0.0;   //  bvx = initial electron-b velocity in x direction
  
 double wna=0.0;    double ora=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

  if ( shua < 0.05 && avz > 0.0 ) { // electron-a moves 1/4-orbit
  shua = 0.1;  kax=ax; kaz=az;
 
  }

   nutna=nutna+1.0;  
           
     // calculate distance ( meter )
                                   
  ara=Math.sqrt(ax*ax+az*az)*1.0e-14; // between electron-a-nucleus-a
  arb=Math.sqrt((nuc-ax)*(nuc-ax)+az*az)*1.0e-14;  // electron-a-nucleus-b
  bra=Math.sqrt(bx*bx+by*by)*1.0e-14; // between electron-b-nucleusa
  brb=Math.sqrt((nuc-bx)*(nuc-bx)+by*by)*1.0e-14; // electron-b-nucleus-b
 ab=Math.sqrt((ax-bx)*(ax-bx)+by*by+az*az)*1.0e-14; // 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);

double bhox=(-bx)/(bra*bra*bra)+(nucc-bx)/(brb*brb*brb)+(bx-ax)/(ab*ab*ab);

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

avx=avx+1.0e-32*acc*hehox; // add acceleration (MM/SS^2) to velocity

double ahoz=(-az)/(ara*ara*ara)+(-az)/(arb*arb*arb)+(az)/(ab*ab*ab);

double bhoy=(-by)/(bra*bra*bra)+(-by)/(brb*brb*brb)+(by)/(ab*ab*ab);

double hehoz=(ahoz-bhoy)*0.5;

avz=avz+1.0e-32*acc*hehoz; // add acceleration in z direction


  bvx=avx;  bvy=-avz;


nut=(bx)/(bra*bra*bra);
  nutba=nutba+nut;  // force between electron-b and nucleus-a
  nut=(bx-nucc)/(brb*brb*brb);
  nutbb=nutbb+nut;

nut=(ax)/(ara*ara*ara);
  nutaa=nutaa+nut;  // force between electron-a and nucleus-a
  nut=(ax-nucc)/(arb*arb*arb);
  nutab=nutab+nut;

vk=Math.sqrt(avx*avx+avz*avz);   ora=ora+vk;               
    leng=vk*1.0e-14;      // moving length (m) for 1 SS
     wav=h/(me*vk*1.0e9);  // de Broglie wavelength (m) 
    wna=wna+leng/wav;
 
  if (shua < 1.0 ) {
                  

  if ( az > 0 ) { // electron-a moves half orbit
  shua = 2.0; tnutaa=nutaa/nutna; tnutab=nutab/nutna;
  tnutna=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

  tnutba=nutba/nutna; tnutbb=nutbb/nutna;
 
  
  }
  }

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

  
   } while (shua < 1.0 );             

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

 System.out.print("\n");
 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.print("\n");
  System.out.printf("one orbit de Broglie wavelength:%.6f\n", twna);


System.out.print("\n");
System.out.printf("1/4-electron-a-x:%.6f\n", kax);
 System.out.printf("1/4-electron-a-z:%.6f\n", kaz);


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

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

System.out.print("\n");
 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);

 System.out.print("\n");
 
   }  

}}