To understand binding energy and mass defects in nuclei, it helps to understand where the mass of the proton comes from.
The news about the recent Higgs discovery emphasizes that the Higgs mechanism gives mass to elementary particles. This is true for electrons and for quarks which are elementary particles (as far as we now know), but it is not true for protons or neutrons or for nuclei. For example, a proton has a mass of approximately $938 \frac{\mathrm{MeV}}{c^2}$, of which the rest mass of its three valence quarks only contributes about $11\frac{\mathrm{MeV}}{c^2}$; much of the remainder can be attributed to the gluons' quantum chromodynamics binding energy. (The gluons themselves have zero rest mass.) So most of the "energy" from the rest mass energy of the universe is actually binding energy of the quarks inside nucleons.
When nucleons bind together to create nuclei it is the "leakage" of this quark/gluon binding energy between the nucleons that determines the overall binding energy of the nucleus. As you state, the electrical repulsion between the protons will tend to decrease this binding energy.
So, I don't think that it is possible to come up with a simple geometrical model to explain the binding energy of nuclei the way you are attempting with your $\left(1\right)$ through $\left(15\right)$ rules. For example, your rules do not account for the varying ratios of neutrons to protons in atomic nuclei. It is possible to have the same total number of nucleons as $\sideset{^{56}}{}{\text{Fe}}$ and the binding energies will be quite different the further you move away from $\sideset{^{56}}{}{\text{Fe}}$ and the more unstable the isotope will be.
To really understand the binding energy of nuclei it would be necessary to fully solve the many body quantum mechanical nucleus problem. This cannot be done exactly but it can be approached through many approximate and numerical calculations. In the 1930's, Bohr did come up with the Liquid Drop model that can give approximations to the binding energy of nuclei, but it does fail to account for the binding energies at the magic numbers where quantum mechanical filled shells make a significant difference. However, the simple model you are talking about will be incapable of making meaningful predictions.
EDIT: The original poster clarified that the sign of the binding energy seems to be confusing. Hopefully this picture will help:
$\hspace{75px}$
.
This graph shows how the potential energy of the neutron and proton that makes up a deuterium nucleus varies as the distance between the neutron and proton changes. The zero value on the vertical axis represents the potential energy when the neutron and proton are far from each other. So when the neutron and proton are bound in a deuteron, the average potential energy will be negative which is why the binding energy per nucleon is a negative number - that is we can get fusion energy by taking the separate neutron and proton and combining them into a deuteron. Note that the binding energy per nucleon of deuterium is $-1.1 \, \mathrm{MeV}$ and how that fits comfortably in the dip of this potential energy curve.
The statement that $\sideset{^{56}}{}{\text{Fe}}$ has the highest binding energy per nucleon means that lighter nuclei fusing towards $\text{Fe}$ will generate energy and heavier elements fissioning towards $\text{Fe}$ will generate energy because the $\text{Fe}$ ground state has the most negative binding energy per nucleon. Hope that makes it clear(er).
By the way, this image is from a very helpful article which should also be helpful for understanding this issue.
Best Answer
It isn't possible to measure potential energy because it has a (global) gauge symmetry. It's like trying to measure the height of a mountain - this could be the height above sea level, the height relative to the deepest sea trench, the height relative to the centre of the earth and so on. Any measurement can only measure the change in potential energy, and of course a change can be positive or negative i.e. the potential energy can increase or decrease.
This is my preferred way to understand the mass defect: start with the separate particles of combined mass $M$ at large distance. As the particles come together they accelerate due to the attractive force between them, so at the moment all the particles meet they have a high kinetic energy, $T$. To get the particles to form a bound state we need to remove this kinetic energy $T$, and that means taking out a mass of $m = T/c^2$ - hence the mass defect.
The total kinetic energy $T$ is just the negative of the potential energy change $\Delta U$ as the particles are brought from infinity into the bound state. It makes no difference what value you assign to $U(\infty)$ because all that matters is the change $\Delta U$.
Response to user52153's comment:
Suppose we start with two particles that attract each other. This could be an electron and proton, or two nucleons, though for convenience I'll assume they have same mass $m_0$. Start with the particles stationary and far enough apart that any interaction is negligable. Then the total energy is just:
$$ E_0 = 2m_0c^2 $$
Now let the particles move together under their mutual attraction. The attractive force will accelerate them, so some time later they will have a kinetic energy $T$. Since we haven't put any energy in or taken any energy out, the total energy must be unchanged so there must be a (negative) potential energy $U$ such that:
$$ 2m_0c^2 + T + U = 2m_0c^2 = E_0 $$
In other words $T + U = 0$ or $T = -U$, and remember that $U$ is a negative number.
At this point there is no mass deficit because the total energy is unchanged. Now user52153 proposes we use some mechanism to take the kinetic energy out of the system. It doesn't matter exactly how we do this - we use the kinetic energy to do work $W$ external to our system of two particles and we continue doing this work until $W = T$ so at this point the particles have been brought to rest. Because we have taken work $W = T$ out of the system the total energy is now:
$$ E_1 = 2m_0c^2 + U $$
The mass deficit $\Delta m = (E_0 - E_1)/c^2$ so:
$$ \Delta_m = \frac{E_0 - E_1}{c^2} = \frac{(2m_0c^2) - (2m_0c^2 + U)}{c^2} = -\frac{U}{c^2} $$
and because $U$ is a negative number that means $\Delta m$ is a positive number i.e. the total mass of the system has decreased.