Electronically polarizable atomistic molecular dynamics models have been developed for the 1-alkyl-3-methylimidazolium nitrate ionic liquids with the alkyl-chain length ranging from 2 to 12. The molecular dynamics simulations with these models confirm the spatial heterogeneity previously discovered by the multiscale coarse-grained models (Wang & Voth, J. Am. Chem. Soc. 2005, 127, 12192). The global structures are monitored by a heterogeneity order parameter. The tail groups of cations are found to aggregate and distribute more heterogeneously than the headgroups of cations and the anions, forming discrete tail domains. The anions always stay close to the headgroups. The charged groups retain their local structures relatively unchanged, forming a continuous polar network, and leading to different global structures when varying the side-chain length. The ionic diffusion is slower with increasing alkylchain length. Based on these simulation results, a refined mechanism can be proposed, considering the competition among the Coulombic interactions, the collective short-range interactions, and the geometrical constraint from the intramolecular bonds. This mechanism can explain many experimental and simulation results, and is expected to be general for most kinds of ionic liquids. The finite size effects due to the limited simulation size are also discussed.