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 hymol, so save this text editor as "hymol.java", and compile it.)
In this program, we first input the distance between two H atoms (= angstrom ) and exponent parameter Z, and the program outputs H2 bond energy and antibond energy change from the original 2 × independent H atoms separated by infinite distance.
import java.util.Scanner;
class hymol {
public static void main(String[] args) {
// H2 molecular bond energy, (anti)symmetric
Scanner stdIn=new Scanner(System.in);
System.out.println("distance between two H (= x angstrom )？ ");
double r=stdIn.nextDouble();
System.out.println("exponential parameter Z ？ ");
double z=stdIn.nextDouble();
System.out.println("Euler exponential integral E(-2r) ？ ");
double rha=stdIn.nextDouble();
System.out.println("Euler exponential integral E(-4r) ？ ");
double rhb=stdIn.nextDouble();
r=r*z/0.5291772;
double me=9.1093826e-31; // electron mass
double pai=3.141592653589793;
double epsi=8.85418781787346e-12; // permittivity
double h=6.62606896e-34; // Planck constant
double ele=1.60217653e-19; // electron charge
double aa=5.291772e-11/z; // Bohr radius/z
double ab=5.291772e-11;
double ea = Math.exp(-r); double eb = Math.exp(-2.0*r);
double ec=Math.exp(r);
double ss = ea * ( 1.0 + r + r*r/3.0 ); // overlap integral
System.out.println(" ");
System.out.printf("overlap integral (= S ):%.6f\n", ss);
// nucl= Coulomb between two nuclei (eV)
double nucl=(ele*ele)/(4.0*pai*epsi*r*aa)*6.241509e18;
System.out.printf("Coulomb potential between two nuclei (eV):%.6f\n", nucl);
// hc = Coulomb integral between electron and ther other nucleus
double hc = -(ele*ele)/(4.0*pai*epsi*aa) * ( 1.0/r - eb*( 1.0+ 1.0/r) );
// hcc = electron-nucleus correction when z is not 1.0
double hcc=((z-1.0)*ele*ele)/(aa*4.0*pai*epsi);
double hccd=hc*6.241509e18; // J to eV
System.out.printf("Coulomb between electron and the other nucleus (eV):%.6f\n", hccd);
hc=hc+hcc; // hc = total Coulomb integral
double hee = 1.0/r + 11.0/8.0 + (3.0*r)/4.0+r*r/6.0;
double he= (ele*ele)/(4.0*pai*epsi*aa) * (1.0/r -eb*hee );
// he = Coulomb integral between two electrons
double heea = he*6.241509e18; // J to eV
System.out.printf("Coulomb integral between two electrons (eV):%.6f\n", heea);
double sc = (-2.0+z)*(ele*ele)/(4.0*pai*epsi*aa) * ea *(1.0+r);
// sc = exchange integral between electron and nucleus
hccd = sc*6.241509e18/(2.0-z); // J to eV
System.out.printf("exchange integral between electron and nucleus (eV):%.6f\n", hccd);
// se = exchange integral between two electrons approx-1
double see = 1.0+r+r*r/3.0;
double se = (ele*ele)/(8.0*pai*epsi*aa) *see*see * eb *(5.0/8.0+1.0/r -eb *hee);
double seee=se*6.241509e18;
System.out.printf("exchange integral between two electrons-approx-1 (eV):%.6f\n", seee);
// exchange integral between two electrons (Euler)
double gamma=0.57721566490153286;
double ka=eb*(25.0/8.0-(23.0*r)/4.0-3.0*r*r-r*r*r/3.0);
double ksa=(1.0+r+r*r/3.0)*ea;
double ksb=(1.0-r+r*r/3.0)*ec;
double kb=6.0*(ksa*ksa*(gamma+Math.log(r))+ksb*ksb*rhb-2.0*ksa*ksb*rha)/r;
double kaaa=ele*ele*(ka+kb)/(5.0*aa*4.0*pai*epsi);
double kaaab=kaaa*6.241509e18;
System.out.printf("exchange integral between two electrons (euler):%.6f\n", kaaab);
double hgroun=ele*ele*ele*ele*me/(8.0*epsi*epsi*h*h);
// hgroun = hydrogen ground state ionization energy
double incre = (hgroun*z*z-2*hgroun*z+hgroun)*6.241509e18;
System.out.printf("H atom energy increase by Z contraction (eV):%.6f\n", incre);
System.out.println(" ");
double jj=2.0*hc+he; double kk=2.0*sc*ss+se;
double tene=(ele*ele)/(4.0*pai*epsi*r*aa)-2.0*hgroun*(z*z-1.0);
double enes=(tene+ (jj+kk)/(1+ss*ss))*6.241509e18;
double enea=(tene+ (jj-kk)/(1-ss*ss))*6.241509e18;
System.out.printf("H2bond-appro-1(eV):%.6f\n", enes);
System.out.printf("H2anti-appro-1(eV):%.6f\n", enea);
kk = 2.0*sc*ss+kaaa;
enes=(tene+ (jj+kk)/(1+ss*ss))*6.241509e18;
enea=(tene+ (jj-kk)/(1-ss*ss))*6.241509e18;
System.out.printf("H2bond-Euler(eV):%.6f\n", enes);
System.out.printf("H2anti-Euler(eV):%.6f\n", enea);
}}