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 hhmol, so save this text editor as "hhmol.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 (= H2 two orbits in the same plane ), 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 hhmol {
public static void main(String[] args) {
// H2 molecule orbit two same-plane orbit
Scanner stdIn=new Scanner(System.in);
System.out.println("electron 1 (=ea) coordinate (MM) from nucleus-a ？ ");
double ea=stdIn.nextDouble();
System.out.println("electronn 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,avy,bvx,bvy;
double ax=ea; // ax = electron-a x-coodinate
double ay=0.0; // ay=electron-a y-coordinate
double bx=nuc+eb; // bx= electron-b x-coordinate
double by=0; // by = electron-b y coodinate
double fax=0.0; double fay=0.0;
double fbx=0.0; double fby=0.0;
double favx=0.0; double favy=0.0;
ab=Math.sqrt((ax-bx)*(ax-bx)+(ay-by)*(ay-by));
// 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 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);
avy=-vyb*1.0e-9; // change m/s to MM/SS
// avy = initial electron-a velocity in y 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 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 kay=0.0;
double axx,ayy,bxx,byy, leng, wav, vk,nut;
do {
axx=ax+avx; ayy=ay+avy; // 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 )
ara=Math.sqrt(ax*ax+ay*ay)*1.0e-14;
// distance between electron-a-nucleus-a (meter)
arb=Math.sqrt((nuc-ax)*(nuc-ax)+ay*ay)*1.0e-14;
// distance between electron-a-nucleus-b
bra=Math.sqrt(bx*bx+by*by)*1.0e-14;
// distance between electron-b-nucleus-a
brb=Math.sqrt((nuc-bx)*(nuc-bx)+by*by)*1.0e-14;
// distance between electron-b-nucleus-b
ab=Math.sqrt((ax-bx)*(ax-bx)+(ay-by)*(ay-by))*1.0e-14;
// distance between 2 electrons
ax=ax*1.0e-14; ay=ay*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
avx=avx+1.0e-32*acc*hehox;
// add acceleration (MM/SS^2) to velocity
double ahoy=(-ay)/(ara*ara*ara)+(-ay)/(arb*arb*arb)+(ay-by)/(ab*ab*ab);
// ahoy = force acting on electron-a in y direction
double bhoy=(-by)/(bra*bra*bra)+(-by)/(brb*brb*brb)+(by-ay)/(ab*ab*ab);
// bhoy = force acting on electron-b in y direction
double hehoy=(ahoy-bhoy)*0.5;
// average force acting on two electrons in y direction
avy=avy+1.0e-32*acc*hehoy;
// add acceleration in y direction
bvx=avx; bvy=-avy;
nut=(bx)/(bra*bra*bra);
nutba=nutba+nut;
// force acting on nucleus-a from electron-b
nut=(bx-nucc)/(brb*brb*brb);
nutbb=nutbb+nut;
// force acting on nucleus-b from electron-b
nut=(ax)/(ara*ara*ara);
nutaa=nutaa+nut;
// force acting on nucleus-a from electron-a
nut=(ax-nucc)/(arb*arb*arb);
nutab=nutab+nut;
// force acting on nucleus-b from electron-a
vk=Math.sqrt(avx*avx+avy*avy);
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;
// number of de Broglie wave contained in each segment
if (shua < 1.0 ) {
if ( ay > 0 ) { // electron-a moved half orbit
shua = 2.0;
// average force acting on each nucleus
tnutaa=nutaa/nutna; tnutab=nutab/nutna;
tnutba=nutba/nutna; tnutbb=nutbb/nutna;
tnutna=nutna;
fax=axx; fay=ayy; // final coordinate of electron-a
fbx=bxx-nuc; fby=byy; // final coordinate of electron-b
favx=avx; favy=avy; // final velocity of electron-a
twna=wna; // a half orbit de Broglie wavelength
}
}
ax=axx; ay=ayy; bx=bxx; by=byy;
} while (shua < 1.0 ); // moved half orbit
System.out.print("\n");
System.out.printf("final-electron-a-x:%.6f\n", fax);
System.out.printf("final-electron-a-y:%.6f\n", fay);
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);
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");
}
}}