function gnp = gaussian_nuclear_potential(g1,g2,Z_nucl,R_nucl)

sum  = g1.exp + g2.exp;
prod = g1.exp * g2.exp;
rat  = prod/sum;

% Location of product gaussian
Rp = (g1.exp*g1.R + g2.exp*g2.R)/(g1.exp+g2.exp);

deltaR = norm(Rp-R_nucl);

gnp = -2 * (2/pi)^(3/2) * pi/sum * Z_nucl * exp(-rat*deltaR^2)*...
                               gaussian_f0(sum*deltaR^2);