function gke = gaussian_kinetic_energy(g1,g2)

sum  = g1.exp + g2.exp;
prod = g1.exp * g2.exp;
rat  = prod/sum;
deltaR = norm(g1.R-g2.R);

gke = rat*(3-2*rat*deltaR^2)*2^(3/2)*(rat/sum)^(3/4)*exp(-rat*deltaR^2);