DNA molecules in the familiar double helical B form are treated here as though they have rod-like structures obtained by stacking the nearly planar base pairs comprising them one on top of another with each rotated by approximately one-tenth of a full turn with respect to its immediate predecessor in the stack. As each base in a base pair is attached to the sugar-phosphate backbone chain of one of the two DNA strands that have come together to form the Watson-Crick structure, and each phosphate group in a backbone chain bears one electronic charge, two such charges are associated with each base pair. Thus, each base pair is subject to not only the elastic forces and moments exerted on it by its neighboring base pairs but also to remote electrostatic forces that, because they are only partially screened out by positively charged counter ions, can render the molecule's equilibrium configurations sensitive to changes in the concentration c of salt in the medium. The observation that the step from one base pair to the next can be one of several distinct types, each having its own mechanical properties that depend on the nucleotide composition of the step, and the assumption that a base pair is rigid, led to the development of a theory of sequence dependent DNA elasticity [Coleman, Olson, and Swigon, J. Chem. Phys. 118, 7127-7140, (2003)]. The theory of DNA molecules in aqueous solution developed here is based on but goes beyond that theory. It takes into account the intramolecular electrostatic interactions of the negatively charged phosphate groups in the molecule and the impenetrability of the DNA molecule for cases in which the electrostatic repulsive forces do not suffice to avoid self penetration. The theory permits one to calculate equilibrium configurations, to determine their stability, and to study the dependence of them on salt concentration and on all kinds of end conditions. When the intramolecular electrostatic forces are taken into account, the equations of mechanical equilibrium for a DNA molecule with N+1 base pairs are a system of mu*N non-linear equations, where mu, the number of kinematical variables describing the relative displacement and orientation of adjacent base pairs is in general 6; it reduces to 3 when base-pair steps are assumed to be inextensible and non-shearable. An efficient numerically stable computational scheme is here presented for solving those equations and determining the mechanical stability of the calculated equilibrium configurations. That scheme is employed to compute and analyze bifurcation diagrams in which c is the bifurcation parameter and to show that, for an intrinsically curved molecule, small changes in c can have a strong effect on stable equilibrium configurations. Cases are presented in which self-contact must be taken into account even though the intramolecular electrostatic forces of repulsion are strong.