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");
}
}}