% <g1 g2|1/r_12|g3 g4>

function gte = gaussian_two_electron(g1,g2,g3,g4)

% r1
sum1  = g1.exp + g3.exp;
prod1 = g1.exp * g3.exp;
rat1  = prod1/sum1;

deltaR1 = norm(g1.R - g3.R);
Rp1 = (g1.exp*g1.R + g3.exp*g3.R)/(g1.exp+g3.exp);

% r2
sum2  = g2.exp + g4.exp;
prod2 = g2.exp * g4.exp;
rat2  = prod2/sum2;

deltaR2 = norm(g2.R - g4.R);
Rp2 = (g2.exp*g2.R + g4.exp*g4.R)/(g2.exp+g4.exp);

% r1/r2 combined
sum_tot    = sum1 + sum2;
prod_tot   = prod1 * prod2;
deltaR_tot = norm(Rp1-Rp2); 

gte = 16/sqrt(pi) * prod_tot^(3/4) / (sum1*sum2*sqrt(sum_tot)) * ...
        exp(-rat1*deltaR1^2 - rat2*deltaR2^2) * ...
        gaussian_f0(sum1*sum2/sum_tot*deltaR_tot^2);