The properties of short chains of poly(styrene)-co-(styrene sulfonate) are studied using atomistic molecular dynamics simulations with explicit solvent. We study single 8-mers and 16-mers with two species of counterions, Na+ and Mg2+, and for various degrees of sulfonation, f. We find that single trajectories do not efficiently sample configurational space, even for fairly long 100-ns simulations, because of rotational barriers caused by nonbonded interactions. Hamiltonian replica exchange molecular dynamics (HREMD) simulations or averages over multiple trajectories are required in order to obtain equilibrium properties. A polystyrene sulfonate chain adopts collapsed conformations at low f, in which the sulfonate groups are located outside the globule and benzene rings form the inner region, and adopts extended conformations as f is increased. Interestingly, the pair correlation functions between side groups of polystyrene chains are not sensitive to f and species of counterion, i.e., the balance of electrostatic repulsion between charged groups and hydrophobic attraction between benzene rings is achieved by conformational change in a way preserving pair correlations between side groups in a polymer chain. For Na+ counterions, no localization is observed in the simulations. For Mg2+ counterions, there is a large free energy barrier to contact pair formation between the sulfonate groups and the Mg2+ counterions. As a consequence we do not observe the formation or breaking of contact pairs during the course of a simulation. The simulations provide insight into the important interactions and correlations in polyelectrolyte solutions.