Static and dynamic properties of low density outer envelopes of neutron stars are calculated within the nonlinear magnetic Thomas-Fermi model, assuming degenerate electrons. A novel domain decomposition enables proper description of lattice symmetry and may be seen as a prototype for the general class of problems involving nonlinear charge screening of periodic, quasi-low-dimensionality structures, e.g. liquid crystals. We describe a scalable implementation of the method using Hypre. Phase velocity of long wavelength transverse phonons is found to be a factor of 5-7 larger than in the corresponding Coulomb crystal model, which could have implications for low temperature phonon-mediated thermal conductivity. Other findings include $c'<0$ elastic instabilities for both bcc and fcc lattices, reminiscent of the situation in some light actinides, and suggestive of a symmetry-lowering transition to a tetragonal or orthorhombic lattice.