  We construct general relativistic models of stationary, strongly magnetized
neutron stars. The magnetic field configuration, obtained by solving the
relativistic Grad-Shafranov equation, is a generalization of the twisted torus
model recently proposed in the literature; the stellar deformations induced by
the magnetic field are computed by solving the perturbed Einstein's equations;
stellar matter is modeled using realistic equations of state. We find that in
these configurations the poloidal field dominates over the toroidal field and
that, if the magnetic field is sufficiently strong during the first phases of
the stellar life, it can produce large deformations.
