The rate for the fusion process p+p -->d+e(+)+nu (e) is calculated using no
nrelativistic 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. C
oulomb interactions between the incoming protons are included nonperturbati
vely in a systematic way. The resulting fusion rate is independent of speci
fic models and wave functions for the interacting nucleons. At this order i
n the effective Lagrangian there is an unknown counterterm which limits the
numerical accuracy. Assuming the counterterm to have a natural magnitude,
the result is consistent with previous nuclear physics calculations.