Chemical bonding/code

Matlab® code to plot figure 2 of chemical bonding (Heitler-London curves). See J. C. Slater, Quantum Theory of Molecules and Solids, vol. 1, McGraw-Hill, New York (1963), Table 3.1 and L. Pauling and E. Bright Wilson, Introduction to Quantum Mechanics, McGraw-Hill, New York (1935), pp. 342, 343.

function HL
close all;
R = [0.5  :0.1  :4.5];
[Eplus, Emin, E0] =  VB(R);

plot(R, [Eplus; Emin; E0], 'linew', 2)
line([R(1) R(end)], [0, 0], 'color' , 'k', 'linew', 1.5, 'linest', '--')

xlabel('R/a_0', 'fonts', 16)
ylabel('E/E_H', 'fonts', 16)
ht = text(0.92,    2, 'E_-', 'fonts', 14 );
ht = text(0.85,  0.5, 'E_0', 'fonts', 14 );
hs = text(2.0, -0.36, 'E_+', 'fonts', 14 );

function [Eplus, Emin, H0] = VB(R)
% Heitler-London equation for H2
% Slater, Molecules and Solids, vol. 1, p. 50
% Slater's equation for K' has a typo: w^2/3  --> w^3/3
a  = 1.0;
g  = -psi(1);     % Euler's constant
w  = a*R;

S  = exp(-w).*(1+w+w.^2/3);
Sp = exp( w).*(1-w+w.^2/3);

J  = 2*(-1./w + exp(-2*w).*(1+1./w));
K  = -2*(exp(-w).*(1+w));
Jp = 2*(1./w - exp(-2*w).*(1./w + 11/8 + 3*w/4 + w.^2/6) );
Kp = (2/5)*(-exp(-2*w).*(-(25/8)+ (23*w)/4 +3*(w.^2) + (w.^3)/3) ...
     +6./w.*(S.^2.*(g+log(w))+Sp.^2.*Ei(-4*w) -2*S.*Sp.*Ei(-2*w) ));

H0 = 2*J+Jp+2./R;
H1 = 2*K.*S + Kp + 2*(S.^2)./R;

Eplus =  (H0+H1)./(1+S.^2);
Emin  =  (H0-H1)./(1-S.^2);
return
function E = Ei(x)
E = real(-expint(-x));