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 (= 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.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 = 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=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
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=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 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 )
ara=Math.sqrt(ax*ax+az*az)*1.0e-14;
// distance between electron-a-nucleus-a (meter)
arb=Math.sqrt((nuc-ax)*(nuc-ax)+az*az)*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)+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);
// 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 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
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 acting on nucleus-a from electron-b
nut=(bx-nucc)/(brb*brb*brb);
nutbb=nutbb+nut;
// force acting on nucleus-b from electron-a
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+avz*avz);
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 ( 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;
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
}
}
ax=axx; az=azz; 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-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);
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");
}
}}