  The rate for the fusion process $p + p \ra d + e^+ + \nu_e$ is calculated
using non-relativistic effective field theory. Including the four-nucleon
derivative interaction, results are obtained in next-to-leading order in the
momentum expansion. This reproduces the effects of the effective range
parameter. Coulomb interactions between the incoming protons are included
non-perturbatively in a systematic way. The resulting fusion rate is
independent of specific models and wavefunctions for the interacting nucleons.
At this order in the effective Lagrangian there is an unknown counterterm which
limits the numerical accuracy of the calculated rate given by the squared
reduced matrix element $\Lambda^2(0) = 7.37$. Assuming the counterterm to have
a natural magnitude, we estimate the accuracy of this result to be 6\% - 8\%.
This is consistent with previous nuclear physics calculations based on
effective range theory and inclusion of axial two-body weak currents. The true
magnitude of the counterterm can be determined from a precise measurement of
the cross-section for low-energy neutrino scattering on deuterons.
