Does the negative potential energy in the gravitational field have to be considered in calculating the total mass of the system in question (because of $E=mc^2$)?
If so it seems to me that the adding this negative energy to the calculation would slightly diminish the strength of the gravitational force over long distances.
You are talking about relativity and gravity together so the question can only be answered in the context of general relativity, but concepts like gravitational potential energy and gravitational force acting over a distance are Newtonian and do not really carry over to general relativity.
However, the gravitational field does contribute to total energy and for gravitationally bound systems it is a negative contribution. If you bring two masses together from a distance energy will be released and the total mass of the combined system will be slightly less. Here mass means the ADM mass of an isolated system measured at a distance.
The statement about the force getting less with distance does not really make sense but the non-linear self interaction of gravity with itself certainly modifies the interaction leading to effects such as the procession of perihelion.
Depends on what you're doing. General relativity handles it for you, in the sense that the Einstein field equation links geometry to the non-gravitational stress-energy tensor. That general relativity is non-linear can be interpreted in part as gravity itself contributing to gravity, but it's generally not even possible to localize gravitational energy in a coordinate-independent manner. The equivalence principle forbids it, since freefall is locally inertial and so any gravitational field can be locally 'transformed away'.
That said, for weak fields the parametrized post-Newtonian formalism includes nonlinearity. For the static case, the first few terms are: $$\mathrm{d}s^2 = -\left(1+2\Phi + 2\beta\Phi^2\right)\mathrm{d}t^2 + \left(1-2\gamma\Phi\right)\left(\mathrm{d}x^2+\mathrm{d}y^2+\mathrm{d}z^2\right)\text{,}$$ where general relativity predicts $\beta = \gamma = 1$.
This is treated in a lot more detail in, e.g., the textbook of Misner, Thorne, and Wheeler. In their notation, the parameter describing the contribution of the gravitational potential itself is called $\beta_2$, cf. Box 39.2 and Section 39.8 to see how it works quantitatively. The wikipedia link provides an overview.
PBS Space Time has a really clear video on this subject: The Real Meaning of E=mc² In the video, Gabe points out that the mass of a Hydrogen atom is less than the sum of the masses of a free proton and electron, and that the mass of an oxygen molecule is less than the sum of the masses of two free oxygen atoms, in both cases because of negative potential energy. Gravitational attraction to the overall system would therefore be less.
Note that this assumes some perfect way of changing the position of the objects within the system.
For a counter-example, if a rock falls off a cliff on Earth, its potential energy is reduced, but other forms of energy result from that loss, in equal measure, and the Earth's gravitational attraction for the Moon is not affected.
